Patent Publication Number: US-11638569-B2

Title: Computer vision systems and methods for real-time needle detection, enhancement and localization in ultrasound

Description:
RELATED APPLICATIONS 
     This application is a continuation application of and claims priority to U.S. Provisional Patent Application No. 62/682,456 filed on Jun. 8, 2018, the entire disclosure of which is expressly incorporated herein by reference. 
    
    
     BACKGROUND 
     Technical Field 
     The present disclosure relates generally to the field of computer vision technology. More specifically, the present disclosure relates to computer vision systems and methods for real-time needle detection, enhancement and localization in ultrasound. 
     Related Art 
     Minimally invasive procedures such as regional anesthesia and interventional oncology involve the insertion of a needle toward target anatomy. In practice, image guidance is used to improve targeting accuracy. Of all imaging modules, ultrasound (“US”) is ubiquitously used due to its real-time, low-cost and radiation-free capabilities. However, with the conventional 2D ultrasound, aligning the needle with the ultrasound imaging plane at steep angles and deep insertions is difficult. 
     Therefore, there is a need for computer vision systems and methods for real-time needle detection, enhancement and localization in ultrasound, thereby improving the ability of computer vision systems to detect a needle at steep angles and deep insertions. These and other needs are addressed by the computer vision systems and methods of the present disclosure. 
     SUMMARY 
     The present disclosure relates to computer vision systems and methods for real-time needle detection, enhancement and localization in ultrasound. In a first embodiment, the computer vision system will detect a placement of an object in a digital image. Specifically, the computer vision system will perform a needle detection phase, a needle trajectory localization phase, and a tip localization phase. In a second embodiment, the computer vision system will enhance and localize a steeply inserted needle. Specifically, the computer vision system will perform a tip enhancement phase, a tip localization phase, and a shaft enhancement phase. In a third embodiment, the computer vision system will automatically detect and enhance needles under 3D ultrasound guidance. Specifically, the computer vision system will perform a needle detection phase, a needle enhancement phase, and a tip localization phase. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The foregoing features of the invention will be apparent from the following Detailed Description of the Invention, taken in connection with the accompanying drawings, in which: 
         FIG.  1    is a flowchart illustrating overall process steps carried out by a first embodiment of the system of the present disclosure; 
         FIG.  2    is a flowchart illustrating step  12  of  FIG.  1    in greater detail; 
         FIG.  3    is a diagram illustrating the layers of the first convolution network and a second convolution network; 
         FIGS.  4 A,  4 B,  4 C, and  4 D  are images illustrating a digital image and feature maps from convolution layers; 
         FIG.  5    is a table illustrating an example architecture of the fast R-CNN; 
         FIG.  6    is a flowchart illustrating step  26  of  FIG.  2    in greater detail; 
         FIGS.  7 A,  7 B,  7 C,  7 D,  7 E, and  7 F  are images illustrating needle detection results; 
         FIGS.  8 A- 8 B  are images illustrating needles in an ultrasound image; 
         FIG.  9    is a flowchart illustrating step  14  in greater detail; 
         FIGS.  10 A,  10 B,  10 C, and  10 D  are images illustrating a needle trajectory estimation; 
         FIG.  11    is a set of images illustrating a prediction of a needle trajectory; 
         FIG.  12    is a flowchart illustrating step  16  of  FIG.  1    in greater detail; 
         FIGS.  13 - 14    are images illustrating the needle tip localization process; 
         FIG.  15    is a flowchart illustrating overall process steps being carried out by a second embodiment of the system of the present disclosure; 
         FIGS.  16 - 17    is a diagram illustrating the tip enhancement phase, the tip localization phase and the shaft enhancement phase; 
         FIG.  18    is a flowchart illustrating step  152  of  FIG.  15    in greater detail; 
         FIG.  19    is a set of images illustrating the effect of γ on the patch-wise transmission map; 
         FIG.  20    is a set of images illustrating the derived signal transmission map; 
         FIG.  21    is a set of images illustrating the output of the tip enhancement phase; 
         FIG.  22    is a flowchart illustrating step  154  of  FIG.  15    in greater detail; 
         FIG.  23    is a set of images illustrating the automatic tip detection; 
         FIG.  24    is a set of images illustrating tissue independence of a PS(x,y) image for four different tissue types; 
         FIG.  25    is a flowchart illustrating step  156  of  FIG.  15    in greater detail; 
         FIG.  26    is a set of images illustrating the needle shaft enhancement process; 
         FIG.  27    is a set of images illustrating the shaft enhancement and tip localization at steep insertion angles; 
         FIG.  28    is a flowchart illustrating overall process steps being carried out by a third embodiment of the system of the present disclosure; 
         FIG.  29    is a diagram illustrating the phases of the method  200 ; 
         FIG.  30    is a flowchart illustrating step  202  of  FIG.  28    in greater detail; 
         FIG.  31    is a set of images illustrating the needle phase detection phase; 
         FIG.  32    is a flowchart illustrating step  204  of  FIG.  28    in greater detail; 
         FIG.  33    is a flowchart illustrating step  206  of  FIG.  28    in greater detail; 
         FIG.  34    is a set of images illustrating the tip localization phase; and 
         FIG.  35    is a diagram illustrating hardware and software component capable of implementing the system of the present disclosure. 
     
    
    
     DETAILED DESCRIPTION 
     The present disclosure relates to computer vision systems and methods for real-time needle detection, enhancement and localization in ultrasound as described in detail below in connection with  FIGS.  1 - 35   . 
       FIG.  1    shows a flowchart illustrating overall process steps being carried out by a first embodiment of the system of the present disclosure, indicated generally at  10 . Specifically, method  10  will detect a placement of an object in a digital image. The digital image can be a ultrasound image, a magnetic resonance imaging (“MRI”) image, a functional MRI (“fMRI”) image, a nuclear magnetic resonance (“NMR”) image, a positron emission tomography (“PET”) scan image, a computer tomography (“CT”) image, an X-ray image, or any other related image. The embodiments described below relate to ultrasound images. It should be understood that any reference to the ultrasound image is only by way of example and that the systems, method and embodiments discussed throughout this disclosure may be applied to any image, including, but not limited to, the digital images listed above. Furthermore, the object in the digital image can be a needle, a catheter, a probe, a measurement tool, a medical implant, a medical rod, or any other insertable object. The embodiments described below will be related to a needle and will refer to the object as a needle. It should be understood that any reference to the needle is only by way of example and that the systems, method and embodiments discussed throughout this disclosure may be applied to any object, including but not limited to the objects listed above. 
     In step  12 , the system performs a needle detection phase. The needle detection phase uses one or more convolutional neural networks (“CNN”) to detect the needle. A neural network, such as a CNN, is a multiple layer network with learnable weights and biases that can be used for, among other things, analyzing visual imagery. CNNs are widely used in machine learning and are an effective tool in various image processing tasks, such as classification of objects. In particular, CNNs can be used as feature extractors to extract different details from images to identify objects in the digital images. In step  14 , the system performs a needle trajectory localization phase. In step  16 , the method performs a tip localization phase. 
     CNNs will be discussed and referenced in this disclosure by way of example. It should be understood that those skilled in the art are capable of adapting the systems and methods discussed herein to use other machine learning systems, including, but not limited to, deep neural networks (“DNNs”), recurrent neural networks (“RNNs”), etc. For example, a DNN using a deep learning based method and the methods discussed herein can detect a placement of an object, such as a needle, in a digital image, such as an ultrasound image, while also reducing irrelevant background and/or background noise. 
     It should be understood that  FIG.  1    is only one potential configuration, and the system of the present disclosure can be implemented using a number of different configurations. The process steps of the invention disclosed herein could be embodied as computer-readable software code executed by one or more computer systems, and could be programmed using any suitable programming languages including, but not limited to, C, C++, C#, Java, Python or any other suitable language. Additionally, the computer system(s) on which the present disclosure may be embodied includes, but is not limited to, one or more personal computers, servers, mobile devices, cloud-based computing platforms, etc., each having one or more suitably powerful microprocessors and associated operating system(s) such as Linux, UNIX, Microsoft Windows, MacOS, etc. Still further, the invention could be embodied as a customized hardware component such as a field-programmable gate array (“FPGA”), application-specific integrated circuit (“ASIC”), embedded system, or other customized hardware component without departing from the spirit or scope of the present disclosure. 
       FIG.  2    shows a flowchart illustrating step  12  of  FIG.  1    in greater detail. In particular,  FIG.  2    illustrates process steps performed during the needle detection phase  12 . In step  22 , the system receives an input image. For example, the system receives a digital image, such as an ultrasound image. A size of the input image size can be chosen to be similar to a smallest detectable object. For example, a 32×32×3 input image is chosen. The image may be preprocessed. For example, ultrasound images from minimally invasive procedures such as biopsies and epidural spinal injections may contain high-intensity artifacts which increase the likelihood of false positives, thus reducing accuracy of needle detection. Preprocessing the image can reduce the influence of the artifacts. In an example, the system can subject a B-mode image I (x, y) to a Top-hat filter using a linear structuring element L (x, y), and the filtered image F (x, y)=I (x, y)−[I (x, y) o L (x, y)], where I (x, y) o L (x, y) denotes an erosion operation. 
     In step  24 , the system processes the digital image through a first set of layers of a first convolution network. The process is illustrated in  FIG.  3   , which is an illustration showing the layers of the first convolution network and a second convolution network. As shown in  FIG.  3   , the first set of layers includes convolution layer  1  ( 32 ), max pooling layer  1  ( 34 ), convolution layer  2  ( 36 ), average pooling layer  1  ( 38 ), and convolution layer  3  ( 40 ). In an example, convolution layer  1  ( 32 ) learns the distinct linear features specific to the needle, while convolution layer  2  ( 36 ) and convolution layer  3  ( 40 ) learn semantic features associated with the needle. A linear feature refers to a coarse and correctly orientated line-like structure (such as, for example, a needle) in the digital image (e.g., ultrasound image). A semantic feature refers to information that gives finer grain inference on the needle. This can include needle tip data, additional data relating to needle dimensions, etc. Semantic features can provide can be more relied upon when needle shaft data is scantily available. Those skilled in the art would understand that any number of layers can be included in the first set of layers.  FIG.  4 A  illustrates an example of a digital image  52 ,  FIG.  4 B  illustrates a feature map from convolution layer  1  ( 54 ),  FIG.  4 C  illustrates a feature map from convolution layer  2  ( 56 ), and  FIG.  4 D  illustrates a feature map from convolution layer  3  ( 58 ). The first convolutional network can be a fast region-based convolutional neural network (“R-CNN”). Further,  FIG.  5    illustrates an example architecture of the fast R-CNN. Those skilled in the art would understand that a different architecture can be used, which includes one or more of, a different number of layers, different layer types, different dimension, etc. 
     Returning to  FIG.  2   , in step  26 , the system generates a set of proposal regions in a feature map of an output from the first convolutional network using the second convolutional network. The output is from a last layer of the first set of layers (e.g., convolution layer  3 ). The one or more regions can be, for example, regions of interest. The second convolutional network can be a regional proposal network (“RPN”). The RPN can generate the one or more regions of interest by ranking potential bounding boxes for the needle (e.g., anchors). 
       FIG.  6    shows a flowchart illustrating step  26  of  FIG.  2    in greater detail. In particular,  FIG.  6    illustrates process steps performed to generate the areas of interest. In step  62 , the system applies a sliding 3×3 convolution window over the feature map to generate anchor boxes. At each window, a maximum of 9 anchor boxes are predicted, generated from 3 scales with a scaling stride of 1.5×(d m , 1.5×d m , 2.25×d m ) and three aspect ratios (1:1, 1:2, 2:1), where d m  corresponds to the minimum dimension of the bounding boxes in the labeled training images. For a feature map of size w×h, the maximum number of anchor boxes is equal to w×h×9. In this example, the maximum number of anchor boxes is about 600. 
     In step  64 , the system minimizes the number of anchor boxes. In an example, the system eliminates cross-boundary anchor boxes. This will leave approximately 60 anchor boxes per ultrasound image. In step  66 , the system assigns a label to each anchor box. For example, the system can assign a positive class label (which indicates a needle is present) when the intersection-over-union (“IoU”) overlap with the needle bounding box in a labeled ultrasound image is greater than 0.7. The system can further assign a negative class label (which indicates no needle is present) if the IoU is less than 0.7. By assigning a class label to each anchor box, the system detects whether a region centered at each sliding window location contains needle data, and the sliding window location encodes coarse localization information with reference to an input ultrasound image. 
     In step  68 , the system maps each window location to an intermediate layer and to two fully connected layers. The intermediate layer can be a 256-dimensional intermediate layer. The two fully connected layers can be two sibling 1×1 fully connected layers, shared across all sliding window locations, for box-classification (e.g., needle or no needle) and box-regression (e.g., for finer localization information). This is also illustrated in  FIG.  3   , which shows a intermediate layer  42 , a classification layer  44 , and a regression layer  46 . In an example, the classification layer can output a maximum of 18 scores, and the regression layer can output a maximum of 36 outputs encoding coordinates of the 9 anchor boxes, for each sliding window location. The outputs of the fully connected layers can be determined by minimizing an RPN multi-task loss function, by using the following formula: 
     
       
         
           
             
               L 
               ⁡ 
               ( 
               
                 
                   p 
                   j 
                 
                 , 
                 
                   p 
                   j 
                   * 
                 
                 , 
                 
                   t 
                   j 
                 
                 , 
                 
                   t 
                   j 
                   * 
                 
               
               ) 
             
             = 
             
               
                 
                   1 
                   
                     N 
                     c 
                   
                 
                 ⁢ 
                 
                   
                     ∑ 
                     j 
                   
                   
                     
                       L 
                       c 
                     
                     ( 
                     
                       
                         p 
                         j 
                       
                       , 
                       
                         p 
                         j 
                         * 
                       
                     
                     ) 
                   
                 
               
               + 
               
                 
                   λ 
                   
                     N 
                     r 
                   
                 
                 ⁢ 
                 
                   
                     ∑ 
                     j 
                   
                   
                     
                       p 
                       j 
                       * 
                     
                     ⁢ 
                     
                       
                         L 
                         r 
                       
                       ( 
                       
                         
                           t 
                           j 
                         
                         , 
                         
                           t 
                           j 
                           * 
                         
                       
                       ) 
                     
                   
                 
               
             
           
         
       
     
     The first term of the above formula describers a box-classifier and the second term of the above formula describes the box aggressor. An anchor index is represented by j, p j  denotes an associated predicted probability, and t j  is the predicted location. A ground-truth label arising from IoU scores is denoted by p* j  ϵ [0, 1], and t* j  denotes an associated location. λ denotes a regularization parameter. L c  denotes a log loss over two classes (needle and no needle), and L r  denotes a regression loss. L r  is a smooth L 1  loss, which is denoted by the following formula: 
     
       
         
           
             
               
                 L 
                 r 
               
               ( 
               
                 t 
                 , 
                 
                   t 
                   * 
                 
               
               ) 
             
             = 
             
               { 
               
                 
                   
                     
                       
                         0.5 
                            
                         
                           
                             ( 
                             
                               t 
                               - 
                               
                                 t 
                                 * 
                               
                             
                             ) 
                           
                           2 
                         
                       
                       , 
                           
                       
                         
                           if 
                           ⁢ 
                               
                           
                             
                               ❘ 
                               &#34;\[LeftBracketingBar]&#34; 
                             
                             
                               t 
                               - 
                               
                                 t 
                                 * 
                               
                             
                             
                               ❘ 
                               &#34;\[RightBracketingBar]&#34; 
                             
                           
                         
                         &lt; 
                         1 
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           
                             ❘ 
                             &#34;\[LeftBracketingBar]&#34; 
                           
                           
                             t 
                             - 
                             
                               t 
                               * 
                             
                           
                           
                             ❘ 
                             &#34;\[RightBracketingBar]&#34; 
                           
                         
                         - 
                         0.5 
                       
                       , 
                           
                       otherwise 
                     
                   
                 
               
             
           
         
       
     
     In an example, N c  is set to 256, N r  is set to 600 (the approximate number of total anchor boxes locations, and the regularization parameter λ is set to 10). Since the derived RPN proposals may overlap, the system uses a non-maximum suppression based on p j , using a threshold of 0.8. In step  60 , the system generates a set of proposal regions. In an example, the proposal regions are top-N ranked. 
     Returning to  FIG.  2   , in step  28 , the system processes the proposal regions using the first convolution network (e.g., the fast R-CNN). Specifically, the R-CNN generates an overall classification and tightened bounding boxes. In step  30 , the system produces needle detection results. The needle detection results can include a marked region or interest with a detection score. This is shown in  FIGS.  7 A,  7 B,  7 C,  7 D,  7 E, and  7 F , which illustrate images of the needle detection results. The needle shaft is accurately localized, despite, for example, low needle shaft intensity (as seen in  FIGS.  7 A,  7 B, and  7 C ) or imperceptible needle shaft intensity (as seen in the bottom row of  FIGS.  7 D,  7 E, and  7 F ). The numbers associated with the bounding boxes in  FIGS.  7 A-F  are the detection score(s). The detection scores are a measure of the confidence of detection. When multiple detection exist, the system selects the highest detection score.  FIG.  8 A  shows an image of a needle in an ultrasound image, and  FIG.  8 B  shows an image of two needle detections for a same needle, each with a confidence score. Again, the system selects the highest detection score. 
     To achieve a unified network, the RPN and the fast R-CNN can be trained using a 4-step alternating process. In step 1, the RPN is trained end-to-end. In step 2, the fast R-CNN is trained using region proposals derived from step 1. In step 3, the RPN is retrained using shared weights of the fast R-CNN from step 2. In step 4, while keeping the shared layers fixed, the fast R-CNN is retrained using the updated region proposals from step 3. 
       FIG.  9    shows a flowchart illustrating step  14  of  FIG.  1    in greater detail. In particular,  FIG.  9    illustrates process steps performed during the needle trajectory localization phase. In step  82 , the system constructs a phase-based image descriptor (phase symmetry (PS (x, y))). In an example, the phase-based image descriptor is constructed using an orientation-tuned 2D Log-Gabor filter bank applied to the marked region of interest on an ultrasound image (e.g., the output from the fast R-CNN). A function of the filter bank is denoted as follows: 
     
       
         
           
             
               LG 
               ⁡ 
               ( 
               
                 ω 
                 , 
                 θ 
               
               ) 
             
             = 
             
               
                 exp 
                 ( 
                 
                   
                     - 
                     
                       
                         log 
                         ⁡ 
                         ( 
                         
                           ω 
                           κ 
                         
                         ) 
                       
                       2 
                     
                   
                   
                     2 
                     ⁢ 
                     
                       
                         log 
                         ⁡ 
                         ( 
                         
                           σ 
                           ω 
                         
                         ) 
                       
                       2 
                     
                   
                 
                 ) 
               
               ⁢ 
               
                 exp 
                 ⁡ 
                 ( 
                 
                   
                     - 
                     
                       
                         ( 
                         
                           θ 
                           - 
                           
                             θ 
                             m 
                           
                         
                         ) 
                       
                       2 
                     
                   
                   
                     2 
                     ⁢ 
                     
                       
                         ( 
                         
                           σ 
                           θ 
                         
                         ) 
                       
                       2 
                     
                   
                 
                 ) 
               
             
           
         
       
     
     Frequency/orientation coordinates are denoted as (ω, θ), k is a center frequency, σ ω  is a bandwidth of the frequency spectrum, σ θ  is an angular bandwidth, and θ m  is a specific filter orientation. The filter parameters can be either tuned or based on priori data. For example, the priori data can include an estimate of the needle trajectory or insertion side of the needle (for θ m ) and a fixed region of interest containing the shaft, close to a transducer surface. 
     In step  84 , the system determines a region of interest. In an example, the region of interest corresponds to the needle bounding box from the needle detection phase. In step  86 , the system determines an insertion side and an insertion angle. In an example, the system determines the insertion side and the insertion angle from the bounding box. Specifically, the system defines the bounding box with parameters x i , y i , L, and W. The system then sets (x c , y c ) as a center of the ultrasound image. For right side insertions, x i &lt;x c , and an estimate of the needle trajectory, β is given by β=tan −1  (W/L). For left side insertions, x i &gt;x c , and an estimate of the needle trajectory, β is given by β=90+tan −1  (W/L). The filter bank is applied with 3 scales and 3 orientations. Θ m =[β−10, β, β+10], and yields a PS (x, y) image containing a prominent needle feature. 
     In step  88 , the system eliminates artifacts not belonging to needle from the PS (x, y) image. In an example, the system extracts the longest connected component. This yields a PS L  (x, y) image, which includes a distinct, intensity invariant straight feature. In step  90 , the system applies Hough transform (“HT”) to the ultrasound image to determine the trajectory. In step  92 , the system determines a trajectory error. In an example, the system determines the trajectory error by comparing the determined trajectory with a gold-standard trajectory estimated by an expert sonographer. For example, the system can denote an original ultrasound image as I (x, y) m×n  where m and n are horizontal and vertical dimensions, and the center of the image is estimated as (x c , y c )=(m/2, n/2). The system then calculates an angle subtended by the automatically detected trajectory on the horizontal axis (γ 1 ), an angle subtended by the trajectory labeled by an expert on the horizontal axis (γ 2 ), a shortest distance between the detected trajectory and the center of the image (λ 1 ), and a shortest distance between the expert-labeled trajectory and the center of the image (λ 2 ). The trajectory error is then quantified using γ 1 −γ 2  and λ 1 −λ 2 . 
       FIGS.  10 A,  10 B,  10 C, and  10 D  are images showing a needle trajectory estimation.  FIG.  10 A  is an illustration  92  showing ultrasound image,  FIG.  10 B  is an illustration  94  showing a preprocessed image with a marked region of interest,  FIG.  10 C  is an illustration  96  showing an image after the extraction of local-phase features and a determined longest connected component, and  FIG.  10 D  is an illustration  98  showing an estimated trajectory after applying the Hough transform. 
       FIG.  11    is an illustration showing a prediction of a needle trajectory. Illustration  102  shows an estimation of needle insertion side and trajectory. Illustration  104  shows the parameters use in calculating needle trajectory error. Illustration  106  shows a trajectory from the method described in step  14  in solid line, and an expert-labeled overlaid on an ultrasound image in dashed line. 
       FIG.  12    shows a flowchart illustrating step  16  of  FIG.  1    in greater detail. In particular,  FIG.  12    illustrates process steps performed during the tip localization phase. In step  112 , the system generates a mask of the trajectory region. In step  114 , the system convolves a mask of the trajectory region with the preprocessed ultrasound image. In step  116 , the system performs a line fitting function using a maximum likelihood estimation sample and consensus (“MLESAC”) algorithm. In step  118 , the system localizes the needle tip by filtering the resulting image with a 2D Log-Gabor filter, and performs a statistical search along the trajectory line. 
       FIG.  13    is an illustration showing the needle tip localization process. Illustration  122  shows an ultrasound image. Illustration  124  shows an extended trajectory region computed with the Hough transform, which is used to form a trajectory mask. Illustration  126  shows an output of the MLESAC algorithm. Illustration  128  shows an enhanced needle image. Illustration  130  shows the localized tip  132 .  FIG.  14    is another illustration showing the needle tip localization process. 
       FIG.  15    shows a flowchart illustrating overall process steps being carried out by a second embodiment of the system, indicated generally at  150 . Specifically, method  150  will enhance and localize a steeply inserted needle. In step  152 , the system performs a tip enhancement phase. In step  154 , the system performs a tip localization phase. In step  156 , the system performs a shaft enhancement phase. Each phase will be discussed in greater detail below.  FIGS.  16  and  17    are an illustration showing the tip enhancement phase, the tip localization phase and the shaft enhancement phase in more detail. 
       FIG.  18    shows a flowchart illustrating step  152  of  FIG.  15    in greater detail. In particular,  FIG.  18    illustrates process steps performed during the needle detection phase. In step  162 , the system performs an image restoration. In an example, the system can account for attenuation and scattering by using a linear interpolation model. For example, the following formula can be used: I (x, y)=t (x, y) I e  (x, y)+(1−t (x, y))v, where I (x, y) is an ultrasound image, t (x, y) is a depth-dependent signal transmission map function (representing response of a loss field in a transmission medium), I e  (x, y) is an ultrasound image intensity to be recovered and v is a constant intensity equal to echogenicity of the tissue confining the needle. It should be noted that a bold notation denotes a matrix. An enhanced image is generated by the system using the following formula: 
     
       
         
           
             
               
                 I 
                 e 
               
               ( 
               
                 x 
                 , 
                 y 
               
               ) 
             
             = 
             
               
                 
                   
                     I 
                     ⁡ 
                     ( 
                     
                       x 
                       , 
                       y 
                     
                     ) 
                   
                   - 
                   v 
                 
                 
                   
                     [ 
                     
                       max 
                       ⁡ 
                       ( 
                       
                         
                           t 
                           ⁡ 
                           ( 
                           
                             x 
                             , 
                             y 
                           
                           ) 
                         
                         , 
                         κ 
                       
                       ) 
                     
                     ] 
                   
                   ∝ 
                 
               
               + 
               v 
             
           
         
       
     
     In the enhanced image formula, k=0.001, where k is a small constant that prevents division by zero. 
     In step  164 , the system generates a patch-wise transmission map. In an example, the patch-wise transmission map is generated by the following formula: ψ (x, y)=I c (x, y)*. I c (x, y) denotes a confidence map, which results from a probability density function that assigns to each pixel in I(x, y) a probability that a random walk emanating from that pixel would be able to reach virtual transducer elements at the top of the ultrasound image, given ultrasound specific constraints, and * denotes a complement. The behavior of the confidence map is controlled by three free parameters, α, β, and γ. The free parameter α denote the attenuation coefficient which controls the depth-dependent attenuation rate. The free parameter β denotes an algorithm constant which affects the robustness and accuracy of segmentation. The free parameter γ models the beam width, imposing a penalty on random walks crossing a horizontal/diagonal edge in the graph with increasing corresponding distance from the starting scanline. 
       FIG.  19    is an illustration showing the effect of γ on the patch-wise transmission map, ψ (x, y), while fixing α=2 and β=90. It should be noted that γ achieves a distinct function with minimum horizontal discontinuities. 
     In step  166  ( FIG.  18   ), the system derives a signal transmission map. In an example, the signal transmission map, t(x,y), is derived by using the following formula: 
     
       
         
           
             
               
                 λ 
                 2 
               
               ⁢ 
               
                 
                    
                   
                     
                       t 
                       ⁡ 
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                     - 
                     
                       ψ 
                       ⁡ 
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                   
                    
                 
                 2 
                 2 
               
             
             + 
             
               
                 ∑ 
                 
                   i 
                   ∈ 
                   Ω 
                 
               
               
                 
                    
                   
                     
                       W 
                       i 
                     
                     ∘ 
                     
                       ( 
                       
                         
                           G 
                           i 
                         
                         * 
                         
                           t 
                           ⁡ 
                           ( 
                           
                             x 
                             , 
                             y 
                           
                           ) 
                         
                       
                       ) 
                     
                   
                    
                 
                 1 
               
             
           
         
       
     
     The equation has two components. The first is the data which measures the closeness of t(x,y) to ψ (x, y), while the second introduces additional contextual constraints on t(x, y). A regulation parameter λ is used to balance the two parts. G i  is a bank of high-order differential operators, consisting of 8 Kirsch filters. The 8 Kirsch filters consist of the same kernel mask rotated in 45-degree increments through all 8 compass directions. Combining the first-order derivative Kirsch filters with a second-order derivative Laplacian mask preserves edge features associated with the needle. W is a weighing function to further constrain t(x,y) in a local image patch. Considering two neighboring pixels, the weighing function is such that W(t 2 (x,y)−t 1 (x,y))≈0. If the two pixels are far apart, then W should be small, and vice versa. The system computes the weighing function from W i (q)=exp(−|G i *I(x,y))) Q | 2 , where Q is a given location in the image. When t 2 (x,y)=t 1  (x,y), W=0. When W=0, the constraint on t(x,y) between neighboring pixels is eliminated.  FIG.  20    is an illustration showing the result of deriving the signal transmission map t (x, y) from ψ (x, y) using various λ. 
     In step  168  ( FIG.  18   ), the system applies a top-hat filter to yield a final tip-enhanced image. In an example, the top-hat filter (TF) computes the morphological opening of I e (x,y) and subtracts it from I e (x,y): TF(I e (x,y))=I e (x,y)−D l [E L (I e (x,y))], where L is a linear structuring element, D L  denotes dilation operations, and E L  denotes erosion operations.  FIG.  21    is an illustration showing the output (e.g., the final tip-enhanced image) of the tip enhancement phase. 
       FIG.  22    shows a flowchart illustrating step  154  of  FIG.  15    in greater detail. In particular,  FIG.  22    illustrates process steps performed during the tip localization phase. In step  172 , the system estimates the needle trajectory. In an example, the system estimates the needle trajectory by, first, extracting shaft information from the enhanced image I e  (x,y) tip . For example, for a 500×500 I e  (x,y) tip  enhanced image, a fixed 100×100 region of interest (I e  ROI (x,y)), is defined on the intersection side of the needle.  FIG.  23    is an illustration showing the automatic tip detection, where, in illustration (a), a fixed region of interest is defined. It should be noted that I e  ROI (x,y) must contain part of the needle shaft, although it need not contain the needle tip. To extract the shaft information from I e  ROI (x,y), the system uses a bank of orientation-tuned band-pass 2D Log-Gabor filters. The filters facilitate local phase-based processing, from which the resulting image descriptor is intensity invariant and therefore insensitive to ultrasound imaging variable such as tissue type. The resulting phase symmetry image PS (x, y), contains distinct local phase features for the needle shaft. The filter parameters are tuned to be optimal. Next, from the PS (x, y), the needle trajectory is estimated using the Radon transform with an angular range of 0°-179°. 
     In step  174  ( FIG.  22   ), the system expands the needle trajectory over the entire I e  (x,y) tip  image, as can be seen in illustration (b) of  FIG.  23   . It should be noted that knowledge of the trajectory region can aid the system in extracting only data lying along the trajectory in I e  (x,y) tip  by convolution. The system then trims the data using a maximum likelihood estimation sample consensus (MLESAC) algorithm, which performs inlier detection and geometrical optimization. The resulting image (I MLESAC (x,y)) is shown in illustration (c) of  FIG.  23   . Next, the system distributes the resulting collinear candidate intensities lying along a line L among a set of line segments, each defined by a set of points or knots denotes as μ 1  . . . μ n . The system then extracts the needle tip using the following formula: 
     
       
         
           
             
               
                 
                   I 
                   needle 
                 
                 ( 
                 
                   
                     I 
                     BP 
                   
                   ( 
                   
                     x 
                     , 
                     y 
                   
                   ) 
                 
                 ) 
               
               = 
               
                 
                   
                     ∫ 
                     
                       μ 
                       i 
                     
                     
                       
                         μ 
                         i 
                       
                       + 
                       1 
                     
                   
                   
                     
                       
                         I 
                         BP 
                       
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                     ⁢ 
                     d 
                     ⁢ 
                     μ 
                   
                 
                 
                   
                      
                     
                       
                         L 
                         
                           ( 
                           
                             
                               μ 
                               i 
                             
                             + 
                             1 
                           
                           ) 
                         
                       
                       - 
                       
                         L 
                         μ 
                       
                     
                      
                   
                   2 
                 
               
             
             ; 
             
               μ 
               ∈ 
               
                 [ 
                 
                   
                     μ 
                     i 
                   
                   , 
                   
                     μ 
                     
                       i 
                       + 
                       1 
                     
                   
                 
                 ] 
               
             
           
         
       
     
     The system obtains I BP  (x, y) by applying a Log-Gabor filter without orientation selectivity to I e  (x, y) tip , where μ 1  and μ i+1  are successive knots. The system, using the above needle tip formula, assigns to pixels between knots μ 1  and μ i+1 , a mean intensity value along L. In step  176 , the system generates I needle  (x, y), which can be seen in illustration (d) of  FIG.  23   . The system then, from I needle  (x, y), localizes the needle tip as the farthest maximum intensity pixel at the distal end of the needle trajectory, which is illustrated in illustration (e) of  FIG.  23   . 
       FIG.  24    is an illustration showing tissue independence of PS (x, y) image for four different tissue types: bovine, porcine, kidney and chicken. The top row shows an ultrasound image I (x, y). The middle row shows a tip enhanced image I e  (x, y) tip . The enhanced tip is surrounded by the circle. The bottom row shows the PS (x, y) image. 
       FIG.  25    shows a flowchart illustrating step  156  of  FIG.  15    in greater detail. In particular,  FIG.  25    illustrates process steps performed during the shaft enhancement phase. In step  182 , the system performs a shaft enhancement using the following formula: 
     
       
         
           
             
               
                 ψ 
                 ⁡ 
                 ( 
                 
                   x 
                   , 
                   y 
                 
                 ) 
               
               shaft 
             
             = 
             
               
                 
                   ❘ 
                   &#34;\[LeftBracketingBar]&#34; 
                 
                 
                   
                     σ 
                     ⁡ 
                     ( 
                     
                       x 
                       , 
                       y 
                     
                     ) 
                   
                   - 
                   
                     max 
                     ⁡ 
                     ( 
                     
                       σ 
                       ⁡ 
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                     ) 
                   
                 
                 
                   ❘ 
                   &#34;\[RightBracketingBar]&#34; 
                 
               
               
                 max 
                 ⁢ 
                 
                   ( 
                   
                     σ 
                     ⁡ 
                     ( 
                     
                       x 
                       , 
                       y 
                     
                     ) 
                   
                   ) 
                 
               
             
           
         
       
     
     In step  184 , the system constrains ψ (x, y) shaft  not to exceed the tip position. The system obtains the signal transmission map for the shaft, t(x,y) shaft , from the signal transmission map formula and shaft restoration is performed using enhanced image formula, where v=max(I(x,y)).  FIG.  26    is an illustration showing the needle shaft enhancement process. It should be noted that t(x,y) shaft  has low intensities along the needle axis and higher intensities for image regions away from the axis. The enhanced shaft arises from a local average of pixels belongings to the shaft along the trajectory.  FIG.  27    is an illustration showing shaft enhancement and tip localization at steep insertion angles. 
       FIG.  28    shows a flowchart illustrating overall process steps being carried out by a third embodiment of the system, indicated generally at  200 . Specifically, method  200  will automatically detect and enhance needles under 3D ultrasound guidance. In step  202 , the system performs a needle detection phase. In step  204 , the system performs a needle enhancement phase. In step  206 , the system performs a tip localization phase.  FIG.  29    is an illustration showing the phases of method  200 . 
       FIG.  30    shows a flowchart illustrating step  202  of  FIG.  28    in greater detail. In particular,  FIG.  30    illustrates process steps performed during the needle detection phase. In step  212 , the system receives a set of ultrasound slices. In step  214 , the system uses a needle detector to classify slices that contain needle data. In an example, the system applies orientation tuned intensity invariant local phase filter banks to each ultrasound slice in the set (hereafter denoted as US volume ) to extract a needle phase descriptor, hereafter denoted as NPD(x,y). The filter banks are constructed from 2D Log-Gabor filters, whose parameters are selected automatically. It should be noted that the insertion side of the needle is known a priori, and the calculation is limited to an automatically selected region of interest. On the insertion side it is expected that the ROI contains a visible part of the shaft. The output of the filter operation generates a phase-based descriptor called phase symmetry, PS(x,y), which is used as an input to the MLESAC algorithm. The system can use the MLESAC algorithm to prune false positive pixels and connect inliers to yield NPD(x,y).  FIG.  31    is an illustration showing the needle detection phase and examples of slices with and without NPD(x,y). The top row shows B-mode US slices constituent of US volume . The middle row shows respective NPD(x,y) images. The slices with needle data possess a salient straight feature with minimum bending. The slices without needle data lack such features. The bottom row shows slice classification results. 
     The needle detector uses an L 2 -Hys (Lowe-style clipped L 2 -norm) contrast normalization on overlapping 33 cells blocks of 4×4 pixel cells. The detector is applied to each of the slices in the US volume  after preprocessing to elicit needle phase descriptors similar to those used in training the needle detector. In step  216 , the system compiles the slices with needle data into a further set. The further set consists of only slices that contain needle data. 
       FIG.  32    shows a flowchart illustrating step  204  of  FIG.  28    in greater detail. In particular,  FIG.  32    illustrates process steps performed during the needle enhancement phase. In step  222 , the system models ultrasound signal transmission in each slice. In an example, each slice is modeled as S(x,y)=S t (x,y)S e (x,y)+(1−S t (x,y))k. S(x,y) is a slice in the US volume , S e (x,y) is the signal transmission map, S e (x,y) is the desired enhanced image, and k is the average intensity of the tissue surrounding the needle in attenuated regions. S t (x,y) is obtained by minimizing the following objection function: 
     
       
         
           
             
               
                 λ 
                 2 
               
               ⁢ 
               
                 
                    
                   
                     
                       
                         S 
                         t 
                       
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                     - 
                     
                       
                         S 
                         a 
                       
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                   
                    
                 
                 2 
                 2 
               
             
             + 
             
               
                 ∑ 
                 
                   i 
                   ∈ 
                   Ϛ 
                 
               
               
                 
                    
                   
                     
                       Γ 
                       i 
                     
                     ∘ 
                     
                       ( 
                       
                         
                           R 
                           i 
                         
                         ⁢ 
                         
                           
                             ★S 
                             t 
                           
                           ( 
                           
                             x 
                             , 
                             y 
                           
                           ) 
                         
                       
                       ) 
                     
                   
                    
                 
                 1 
               
             
           
         
       
     
     S a (x,y) is a patch-wise transmission function representing boundary constraints imposed on the image by attenuation and orientation of the needle, ζ is an index set of image pixels, ° is element wise multiplication, and   is a convolution operator. R i  is a bank of high order differential filters consisting of eight Kirsch filters and one Laplacian filter, and Γ i  is a weighing matrix calculated from Γ i  (x,y)=exp(−|R i (x,y) S(x,y)| 2 ). 
     In step  224 , the system extracts S e (x,y). In an example, the system extracts S e (x,y) from S e (x,y)=[(S(x,y)−k)/[max(S t (x,y),ε)] ρ ]+k. ε is a small constant and ρ is related to the attenuation co-efficient of the tissue. In step  226 , the system subjects each enhanced slice to a Top-hat filter operation using a linear structure element. The final enhanced slices constitute the enhanced sub-volume denoted as USE* volume . 
       FIG.  33    shows a flowchart illustrating step  206  of  FIG.  28    in greater detail. In particular,  FIG.  33    illustrates process steps performed during the tip localization phase. In step  232 , the system performs a 2D tip localization. In an example, when the needle insertion is in the y-z plane, then the x-y plane is parallel to the needle insertion direction. The system determines x′ and y′ from a projection P x,y  since x′ and y′ have the same value in all slices. P x,y  is calculated as the maximum intensity projection (MIP) of USE* volume , by extracting maximum intensity values along optical paths in the z direction. From the projection, the needle tip is localized by determining the phase symmetry PS(x,y) of P x,y  in a region limited to the needle trajectory, applying the MLESAC algorithm for inlier detection and geometrical optimization, and feature extraction on the resultant point cloud using a combination of spatially distributed image statistics which enhance the needle tip. The system then yields the projection enhanced needle image denoted as PE(x,y). (x′,y′) is determined from the first maximum intensity pixel at the distal end of the needle trajectory in PE(x,y). 
     In step  234 , the system performs a scan plane determination to determine a 3D tip location. In an example, the system determines the scan plane by calculating elevation (z) direction of the volume. The system can calculate the scan plane by using the following formula: 
     
       
         
           
             
               ∑ 
               
                 i 
                 = 
                 
                   - 
                   γ 
                 
               
               
                 + 
                 γ 
               
             
               
             
               
                 ∑ 
                 
                   j 
                   = 
                   
                     - 
                     γ 
                   
                 
                 
                   + 
                   γ 
                 
               
                 
               
                 I 
                 ⁡ 
                 ( 
                 
                   
                     
                       x 
                       ′ 
                     
                     + 
                     i 
                   
                   , 
                   
                     
                       y 
                       ′ 
                     
                     + 
                     j 
                   
                 
                 ) 
               
             
           
         
       
     
     The sum of pixel intensities in a bounded square patch of length 2γ centered at (x′,y′) in each slice within USE* volume . The system estimates the scan plane as the slice with the maximum intensity sum. The result produces z′.  FIG.  34    is an illustration showing the tip localization phase. The first column shown a P x,y  image. The second column shows a PE(x,y) image. The automatically localized tip is overlaid on x-y in the third column, and on y-z in the fourth column. The fifth column shows 3D imaging coordinates. The top row shows a moderate insertion angle and the needle aligned with the ultrasound beam. The middle row shows a moderate insertion angle and the needle not aligned with the ultrasound beam. The bottom row shows a steep insertion angle and the needle aligned with the ultrasound beam. 
       FIG.  35    is a diagram showing a hardware and software components of a computer system  250  on which the system of the present disclosure can be implemented. The computer system  250  can include a storage device  252 , computer vision software code  254 , a network interface  256 , a communications bus  258 , a central processing unit (CPU) (microprocessor)  260 , a random access memory (RAM)  264 , and one or more input devices  266 , such as a keyboard, mouse, etc. The computer system  250  could also include a display (e.g., liquid crystal display (LCD), cathode ray tube (CRT), etc.). The storage device  252  could comprise any suitable, computer-readable storage medium such as disk, non-volatile memory (e.g., read-only memory (ROM), eraseable programmable ROM (EPROM), electrically-eraseable programmable ROM (EEPROM), flash memory, field-programmable gate array (FPGA), etc.). The computer system  802  could be a networked computer system, a personal computer, a server, a smart phone, tablet computer etc. It is noted that the computer system  250  need not be a networked server, and indeed, could be a stand-alone computer system. 
     Having thus described the system and method in detail, it is to be understood that the foregoing description is not intended to limit the spirit or scope thereof. It will be understood that the embodiments of the present disclosure described herein are merely exemplary and that a person skilled in the art can make any variations and modification without departing from the spirit and scope of the disclosure. All such variations and modifications, including those discussed above, are intended to be included within the scope of the disclosure. What is intended to be protected by Letters Patent is set forth in the following claims.