Gabor filter-based localization of straight and curved needles in 2D ultrasound images

2D ultrasound (US) is one of the most commonly used medical imaging devices for needle localization in biopsies. However, the produced images are low-resolution and contain an excessive number of artifacts, which makes the needle localization challenging. Image processing techniques can help resolve this issue. This paper presents a novel Gabor filter-based method for needle localization in 2D US images, which enhances the needle outline in the images while suppressing other structures. The scheme works in two stages: First, the Gabor filter is applied to the image, the needle insertion angle is estimated, and the needle trajectory is found using a RANSAC line fitting; then, the Gabor filter is repeated using the estimated insertion angle and the location of the needle tip is estimated using probability mapping. The proposed scheme works with both straight and curved needles. An automatic parameter tuning method, which is used to optimize the threshold value of the Otsu’s method, is also presented here. The tests on this needle localization scheme were done using two different water mixtures, three different gelatin-based phantoms, and ex vivo experiments. The accuracy was assessed using an infrared-camera-based tracking sensor and computed tomography images. The results showed that the suggested localization scheme can be effectively used in the 2D US image-guided needle procedures.


Introduction
Percutaneous needle procedures, such as drug delivery or biopsy, are frequently performed in medical operations.
Accurate needle placement plays a vital role in the success of percutaneous procedures [1]. Incorrect placement of needle might result in an incorrect sample collection and this might result in multiple insertions that might harm the tissue and cause hemorrhage. Needle placement accuracy could be increased by following a predetermined path by detecting and tracking the needle tip.
Percutaneous procedures could be guided by medical imaging such as ultrasound (US), fluoroscopy, computed tomography (CT), and magnetic resonance imaging (MR or MRI). While MRI is a good option for getting high-quality images, but it requires custom-made needles as ferromagnetic equipment cannot be used in the operations. Another disadvantage is its size because it does not allow the needle localization to be run efficiently. The other options, namely the CT and fluoroscopy imaging, are not desirable due to their risks involving radiation. Moreover, as in the case of MRI, CT provides limited workspace for the operations. Considering these drawbacks, US imaging comes up as the better option, without any known harmful effects, in which a small probe can be used to track the needle. The US images, however, generally contain reverberation effects, undesirable artifacts, speckle noise, and acoustic clutter, which reduce the image quality, making the visualization of the needle, especially its tip, difficult. Furthermore, the appearance of the needle is dependent on its insertion angle. Unless the needle is inserted at a proper angle, it would be difficult to image the needle. To address this problem, image processing algorithms are needed to increase needle visibility and eliminate undesirable artifacts to localize the insertion axis of the needle and distinguish the needle's tip.
This paper presents a novel image processing method for automated needle localization in single 2D US images without any prior information. The proposed method uses a Gabor filter to segment the needles in US images. The scheme is a label-free localization method and requires no training or labeling during segmentation. The Gabor filter works in the wavelet domain, transforming the image frame into a sparse matrix that provides computational advantage over other methods. The basis to use the Gabor filter is its ease of implementation and its performance in real time to segment needle pixels in the desired range. We presented an early version of the novel proposed scheme for detecting axis of the straight needles in [2]. A more comprehensive and detailed detection method which also works for curved needles is provided in this paper, and its effectiveness is shown using different types of tissue phantoms. In the improved method, the Gabor filtered images are binarized with an automatic parameter tuning method and the needle trajectory is obtained using curve fitting. Finally, the location of the needle tip is estimated using the probabilistic mapping method. This paper is organized as follows. Relevant studies are described in the next section. Then, Gabor-based line filtering methods for needle segmentation and the needle axis localization are presented. Next, the needle tip estimation methods are described and the modifications to the algorithm to localize curved needles are given. Finally, experiment results are given and the conclusion is presented.

Related studies in the literature
Summary of the studies about needle localization is provided here starting with works on localization algorithms.
Draper et al. [3] segmented the needle pixels using variance mapping method. In order to determine the needle axis and its tip, they used a thresholding technique and principal component analysis. Ding [4] localized the needles using Hough transform in real time using 2D US images. Cheung and Rohling [5] enhanced the needle visibility in 2D US images by steering the US beam. The needle was localized using Hough transform and specific beam patterns were used to maximize the reflections from the needle. Using 3D US images, Barva et al. [6] localized curvilinear objects. First, they converted the 3D US images into binary images and then applied curve fitting to localize the objects. Okazawa et al. [7] localized curved needles in 2D US images using Hough transform. They segmented the needle using the rough insertion angle. Ayvali et al. [8] detected the needle tip using circular Hough transform and tracked the needle tip with optical flow.
Neshat and Patel [9] modeled curved needles using Bezier polynomials, where the coefficient estimations were done using radon transform. They also implemented their curved needle segmentation in real time using GPU. Adebar and Okamura [10] tracked the needle tip using a 3D color-Doppler US machine. Using a piezoelectric buzzer, they vibrated the needle with a high frequency, and then they used color Doppler filters to detect and track its tip.
Dong et al. [11] proposed a framework for real-time needle localization in highly noisy 2D US images.
They used level set and partial differential equation to localize and track the biopsy needle. Uhercik et al. [12] and Zhao et al. [13] first filtered lines in the 3D US images using Frangi's vesselness filter and then determined the axes of the needles using RANSAC line fitting. Aboofazeli et al. [14] used an anisotropic filter to segment the curve shaped needles in 3D US images and ray casting to project 3D images to 2D. Then the curved needles were localized by Hough transform.
Waine et al. [15] visualized the 3D shape of the needles from 2D US images. First, they transversely scanned the medium along the needle path and obtained the possible 3D needle pixels' positions by localizing the needle pixels in each scan. Then, they applied RANSAC curve fitting to eliminate false detections and obtain the shape of the needle. Wen et al. [16] detected the brachytherapy seeds in 3D US images. They used reflected power images instead of conventional B-mode US images. After the seed insertion was detected using Hough transform, the seeds were segmented and localized in the predetermined local search space. Ayvaci et al. [17] used 2D transrectal US images to identify biopsy needles. Prior position of the needle and a model of the image background were used to segment the needle pixels in the images. Cao et al. [18] proposed an automated method to detect catheters in 3D US images. A probability map was created by the catheter's physical model in 3D US images, and then the maximum intensity approach was used to project the map onto the image to distinguish the catheter quicker.
Vrooijink et al. [19] built an automated framework for tracking the needle tip in 3D with 2D transverse US images. An actuated device was used to place the needles while the imaging probe was servoed simultaneously with the needle. Comet tail artifact, a particular view of the needle seen in transverse images, was used to localize the needle tip. Chatelain et al. [20] tracked the needles in real time using 3D US images. The US probe was moved by a 6-DOF robot. The volume intensity differences of the US images were used to detect needle pixels, and then the needle axis was determined using RANSAC.
Renfrew et al. [21] proposed a framework to identify the needles and targets in image-guided robotic interventions. They maximized the information to localize the needle and targets and used a particle filter to estimate the target and needle conditions. The proposed approach was validated by a simulation study.
Mathiassen et al. [22] used difference information of sequential frames and second-order derivative of a 2D Gaussian filter for needle tip estimation in 2D US images. Needle tip prediction accuracy was measured using an optical tracking system. Beigi [23] located a needle in a 2D US image sequence using motion information.
After the motion was detected, the needle axis was obtained using polynomial fitting and Hough transform on the pixel extracted from the difference image. The last detected point along the needle axis in the difference image was considered to be the needle tip.
Mwikirize et al. proposed a convolutional neural network-based needle detection and localization algorithm in 2D US images. They fine-tuned and trained a network to detect regions of interest in images and used these regions to detect the needle trajectory and the tip. In another study, signal transmission maps and statistical optimization were used to model attenuation information and detect the needle tip. They combined this pipeline with a deep learning model to filter out motion events [24][25]. Daoud et al. [26] used the needleinduced motion and edge detection to localize needles in US images. The initial trajectory region is analyzed with ranklet and Radon transforms. Lee et al. [27] used neural networks and spatial information to segment and track the needle tip. Xu et al. [28] proposed a pipeline that includes thresholding, morphological operations, and maximum likelihood estimation sample consensus algorithm.
Novotny et al. [29] tracked marker-attached surgical instruments in 3D and obtained their orientation angles. Novotny used a modified version of Radon transformation implemented in real time to detect the surgical instrument. Chatelain et al. [30] also tracked needles in real time using 3D US imaging. Volume intensity differences between the US images were used to detect tubular pixels. Then, the insertion axis was determined using RANSAC and its tip was localized. The authors also attached the US probe onto a 6-DOF robotic arm and it was visually servoed with respect to needle tip coordinates. Pourtaherian et al. [31] detected and tracked the needle in 3D US images. The authors used Gabor filter for needle detection. After the needle was detected, the needle was tracked using a Gradient descent based tracking algorithm. Hacihaliloglu et al. [32] localized the needle in 2D US images. First, Log-Gabor filter was used for needle segmentation and then segmented pixels were binarized. Finally, a modified version of RANSAC algorithm was used for needle axis localization. Hatt et al. [33] segmented the needle images using a machine learning method. Log-Gabor filter was used to create the classification features and a statistical approach was used to train the classifier. After the needle pixels were classified, Radon transformation was used to find the needle axis and direction.
A primitive and incomplete variant of the Gabor filter-based needle segmentation method was presented in our earlier work [2]. There, the RANSAC line fitting was used for the detection of the needle axis. In [34], the Gabor filter was automatically binarized using US images and the needle tip position was estimated with a probabilistic method. In [35], the Kalman filter was used to track the tip of the needle lost in noisy US images.
This paper extends the work conducted in [2] by providing the complete Gabor filter-based needle localization method. The effectiveness of the method is shown for both straight and curved needles using different types of tissue phantoms. In the next sections of this paper, contribution to the literature is explained in detail.

Gabor filter
Gabor filter is widely used to copy fingerprint and iris shape in the literature [36][37]. Gabor filter is also used to segment the anatomical structures such as liver and the veins [38][39]. In this research, the Gabor filter formulation is adapted from [38]. In the spatial domain, 2D Gabor filter function is defined with the multiplication of a 2D Gaussian envelope and a complex carrier sinusoid: where θ is the orientation of the Gabor filter, σ is the standard deviation, and λ is the wavelength of the modulating sinusoid. σ /λ ratio relates to the bandwidth of the filter to the Gaussian spatial envelope and hence to the spatial frequency bandwidth of the Gabor filter. If σ / λ is equal to 0.56, the maximum filter response is attained [40].
2D Gabor filter can be rewritten using Euler's equation ( e jx = cos(x) + j sin(x) ): The shape of the needle can be assumed to be straight in the most 2D US-guided biopsy procedures. The line structures in the image, I g (x, y) , is revealed by convolution operation (*) given below.
In Figures 1a and 1b the raw image of a needle in gelatin phantom and the image's Gabor filtered output are shown, respectively. The raw 2D US image in gelatin phantom contains artifacts, such as the reverberation artifact placed at the bottom of the image. The Gabor filter enhances the needle pixels while significantly filtering the artifacts.

Needle localization using gabor filter
There are four commonly used methods for needle detection in US images. The first two methods are based on the Frangi vesselness filter [41], Hough [7], and Radon transforms [29]. The third method is variance mapping [3]. These methods improve tube-like shapes, including undesired structures and artifacts in images. Due to this shortcoming, it might be challenging to differentiate the needles from other structures (see Figures 1c and   1d). The fourth method is based on empirical and adaptive thresholding. The threshold value can vary by image and segmenting the needle pixels using this method can create discontinuities in the needle structure in the US image (see Figure 1e and 1f). Because of the drawbacks of these methods, a new method based on Gabor filter [42] is proposed.
The important difference between the proposed method and the literature is that the Gabor filter cancels the structures that are not aligned with the needle direction while improving tube-shaped structures along the needle direction. As a result, brightness of the needle pixels becomes higher than the brightness of the image background pixels. This makes the needle structures prominent and gives better results compared to other needle segmentation methods.

Needle insertion angle estimation
The Gabor filter is applied in 2D US images at an orientation angle, θ , to localize the needle axis and its tip. The filter orientation should be close to the needle insertion angle for desired filtering. Figure 2 shows filtering results of the same image by changing the orientation angle of the Gabor filter, θ , with increments of 30 • . When the Gabor filter angle, θ , is close to the needle insertion angle the artifacts are minimized; hence, needle pixel's continuity and needle visibility is improved.
If unknown, the estimation of the needle insertion angle is quite complex. During US image collection, the bottom of the US probe is pressed against the tissue to keep continuous contact with the skin. The needle can only be inserted from either side of the probe so that the needle tip can be visible in the imaging plane. In medical practice, the needle insertion angle cannot be too steep or too shallow. Therefore, the needle has to be inserted at an angle close to 135°or 225°. Therefore, a binary estimate of the insertion angle is collected as

Image binarization
Image binarization is required for line fitting once the Gabor filter kernel is convolved with the input image to segment the needle pixels. Binarization is performed in three consecutive steps: median filtering, automatic thresholding, and morphological operations. The details of the binarization are given in [2] but are summarized here for completeness of the method. In the first step, a 7 × 7 sized median filter is applied which is empirically tuned for optimal noise suppression and structure clearance. The median filter was preferred over other softening filters due to its satisfactory performance and suitability for use in real time. Thresholding is then applied to binarize the image. The threshold is determined using an entropy-based automatic parameter setting method [34] based on [43]. The method eliminates faults and reduces the processing time during the RANSAC line fitting. Finally, morphologic erosion and dilation processes are performed using a 3 ×3 square-shaped structuring element to ensure continuity of the structures in the image.

Line detection
In the binarization step, the needle pixels can be classified as background, which forms gaps between the needle pixels. This condition makes RANSAC line fitting a good solution. Briefly, RANSAC algorithm starts by randomly selecting points from a point cloud and a straight line is fitted onto selected points. Then, the model parameters are calculated and the consistency of the parameters are evaluated. The same process is applied to the remaining points in the point cloud for the best fit. Even if the needle pixels are discontinuous in the binarized image, the needle axis can effectively be localized. Figures 4a and 4b show the 2D US images obtained after the needles are inserted to the agar, and agar-gelatine phantoms, respectively. The Gabor filter is applied with the user-supplied rough insertion angle estimates, θ 1 , and the images are binarized by Otsu's thresholding method (Figures 4c and 4d). The needles are detected and the exact insertion angles, θ 2 , are calculated in the US image plane (Figures 4e and 4f).
In the second iteration step, the Gabor filter with updated angle, flowed by image binarization (Figures 4g   and 4h) and RANSAC line fitting processes are applied to localize the needle axis (Figures 4i and 4j). Upon experimentation, 2 cm is determined to be the minimum needle insertion length for successful localization for RANSAC, which is below the target depth in the US image.

Estimation of needle tip location
Needle tip localization can be challenging in 2D US images because the needle tip visibility changes with needle's orientation. The bevel-tip needle is used in most biopsy procedures, and it is imaged laterally. In such cases, the needle and the tip may be viewed as a single structure. However, the needle tip may appear as an infinitely small circle when close to anatomical structures or artifacts, where it cannot be distinguished as a separate structure. For these reasons, the exact location of the needle tip cannot be determined only using information obtained from 2D US images, but is estimated from the US images.
In the first stage of the proposed localization method (Figure 3), the needle axis is localized in the US image. The steps of the proposed needle tip estimation method for a straight needle is shown in Figure 5. A region of interest (ROI) is created around the needle pixels with a 6 pixels wide rotated rectangle to decrease the computation time (ROI is shown in Figure 5c). In the second stage, the needle tip is estimated using the information obtained in the previous stage. The filtering process is repeated using the refined needle insertion angle ( θ 2 ). In the final step, probability mapping is used to estimate the needle tip location. The probability of each location within the needle ROI shows belongingness to the actual needle pixels. The number of white pixels and black pixels before and after the corresponding position is used to calculate the probability of a pixel being on the needle. White and black pixels prior the position are classified as confirmed (1:White) and unconfirmed (0:Black) needle pixels, respectively. White and black pixels beyond the location are classified as falsely identified (1:White) and nonneedle (0:Black) pixels, respectively (see Figure 6). The conditional probability of each location within the needle ROI can be calculated using the counts of (1:White) and (0:Black) values. The probability of being an actual needle pixel is calculated as: where P (N, n x,y ) = (P 1 + P 3 )P 1 + (P 2 + P 4 )P 4 (7) P (+ | N, n x,y ) = P 1 (8) P (+, n x,y ) = (P 1 + P 3 )P 1 + (P 2 + P 4 )P 2 It should be noted that with this scheme nonneedle pixels can also be binarized as needle pixels as the probability of these pixels can be high resulting in false identification. In order to eliminate false identifications, the absolute intensity value of the Gabor filter output, | I g (n x,y ) |, is multiplied with the conditional probability of each location (10). This way, the intensity values of the needle pixels become higher than nonneedle pixels. This multiplication increases the probability of being a needle pixel. P (T | +, n xy ) = P (N | +, n xy )· | I g (n x,y ) | In the final step, the probability of each location within the needle ROI is calculated and the needle tip location ( p tip (x, y)) is estimated from the maximum probability: The probability map is drawn in 3D and overlaid onto the raw 2D US image in Figure 7 (left image).

Gabor filtering algorithm for curved needle localization
Up to this point, the geometrical model of the needle is assumed to be straight and the needles are localized accordingly. However, even straight needles slightly bend in tissue, and also flexible needles are also used in biopsies. If the needle segmentation and needle axis localization method (Figure 3 -1st Stage) is modified, the proposed needle localization algorithm can also be used to localize curved needles.
The Gabor filter is applied with a fixed orientation angle, θ , to segment the needle's pixels. For straight needles, all pixels lie on the same line, and all needle pixels can be enhanced. In the case of curved needles, if the Gabor filter is applied with a fixed angle, only some parts of the curved needle pixels will be enhanced, while the rest will be suppressed. To segment all pixels of the curved needle, the Gabor filter is applied with two different orientation angles ( θ a , θ b ) and then the filtered images are summed. As a result, all pixels of the curved needles are segmented. After the segmentation process, the summed image is smoothened using a 7 ×7 sized median filter and then the smoothened image is binarized. In our experiments, θ a and θ b are selected to be 135 • and 225 • , respectively. The needle pixels lying along cross directions can be segmented using these angles. The steps of the proposed localization method for curved shaped needles are shown in Figure 8.
The curvature of the needle can be modeled using a high-order polynomial if the image only consists of pure needle pixels. However, there can be noise in the binarized image, which makes polynomial fitting less desirable for localizing the needle axis. RANSAC curve fitting can be used here, but it is hard to generalize RANSAC for higher dimensions. However, if the binarized image is divided with grids, RANSAC line fitting can be applied to each section (see Figure 8g). Fitting a line to binarized pixels using RANSAC is a robust fitting method because it is a simple parametric model. The needle axis can then be localized by gathering the detected lines. After the needle axis is localized, needle tip estimation (Figure 3 -2nd stage) can be used. Since the geometrical model is known, the needle tip can be estimated using probability mapping. The 3D probability map overlaid on the 2D US image is given in right image of Figure 7, where the point with the highest probable point is chosen as the needle tip.

Results
This section presents the findings of the accuracy and execution time of the proposed needle localization algorithm on 2D US images.

US machine
In this study, a GE 11L linear 2D US probe and a GE LOGIQ P5 ultrasound machine were used to obtain the images. The size of the obtained images were 640 × 480 pixels. During imaging, despeckle filters of the device were closed to obtain raw US images. Hence, most of the acquired images contained different types of artifacts, especially reverberation and mirror artifacts. Moreover, needles in the acquired images were suppressed by artifacts. In addition, needles did not appear as complete structures and their intensity values were close to the intensity values of the background pixels. 22G, 15-cm-long biopsy needles were used for the needle insertion experiments (Wipak Medical).

Phantoms
Distilled water and four different types of phantoms were used to test the proposed method in different consistencies. Details of the preparation of the phantoms can be obtained from [34,44,45].

Accuracy of linear needle localization
In our experiments, the needle was inserted into different types of phantoms, and 8491 single 2D US images were acquired for evaluation of the algorithm. During the experiments, needles were inserted both manually and using a biopsy robot [46]. The needles were placed at least 2 cm in the phantoms at various angles from 15°to 75°to demonstrate that the proposed method can localize the needle tip at different angles. In all images, the needle was successfully detected and the tip was localized.
An infrared motion capture system (OptiTrack, Natural Point, Inc. OR, USA) was used to evaluate the accuracy of the straight needle localization method (see Figure 9)

Accuracy of the curved needle localization
To evaluate the accuracy of curved needle localization method, a computer tomography (CT) scanner with 10 µm accuracy (YXLON MU2000-D) was used. A mold was fabricated to fixate a curved needle. The needle was inserted into the mold to keep it fixed during the experiments. First, the curved needle was scanned by a CT scanner and a 3D needle shape was obtained. Then, the needle was imaged by US in water and its axis was localized using the proposed localization method. Similar to the needle axis localization in CT image, the needle axis was formulated using a quadratic order polynomial equation. The two equations were discretized using 0.01 mm intervals. Coordinates obtained using the US were transformed to the CT coordinate frame. Finally, the Euclidean distance was calculated to measure the axis localization error. Fifty images obtained from both systems were compared. The RMS error of the axis localization was 0.97 mm and the error in needle tip estimation was 1.08 mm.

Discussion of the results
In medical practice, tumors between 5 and 30 mm in diameter are biopsied. Smaller sized tumors are considered benign and larger than 40 mm are considered malignant tumors [47]. There are no specific guidelines for the accuracy of reaching the tumors. During the insertion of the needle, the acceptable localization error for the transverse directions is 1.5 mm and 2.5 mm for the longitudinal direction [16]. In the proposed method, the tip localization error is around 1 mm for both curved and straight needles and is in the acceptable range for a successful biopsy. In the literature, the needle tip localization errors are in the range of 1-3 mm, but it should be noted that the error may vary depending on the phantom or tissue type tested.

Execution time
The proposed needle localization algorithm shown in Figure 3 was implemented in both C++ using the OpenCV library and MATLAB, and was run on a 64-bit Windows 7 workstation (Intel Xeon E5-2620 CPU running at 2 GHz). Tests were conducted using single image frames or frames obtained from image sequences without a priori information where every frame is treated as new. The execution time of the localization algorithm for a single image frames was 35.01 ± 3.04 ms using C++, and 490.01 ± 9.04 ms using MATLAB.

Conclusion
This study presents a needle localization method for straight and curved needles in 2D US images. The needle insertion angle and the needle trajectory are estimated without using an external sensor. The major difference of the proposed method compared to others in the literature is segmenting the needle pixels with a Gabor filter, which can separate the needle from the background image. Moreover, the proposed localization method binarizes the US image using an entropy-based parameter tuning method. Images from different types of backgrounds are binarized automatically with an optimum threshold value.
The proposed method works effectively and fast in real time. Since the method is label-free and deterministic, it was possible to implement it to run automatically. The method is the first in the literature to propose the segmentation of curved or straight needles using the Gabor filter. This method is designed for use with a robotic system so that the 2D US image is always aligned with the inserted needle.
The method is optimized to find a single needle in the US image if multiple needles or needle-like fine objects appear in the image and they are positioned parallel to the needle, identification of the correct needle would be challenging; moreover, it is unlikely that more than one needle will be simultaneously inserted during a biopsy. Another challenge encountered during the development of this method was to find a suitable thresholding method that would work for different contrast images. In all segmentation methods, applying a solution similar to Otsu-based thresholding, as in this study, is necessary. In addition, the proposed algorithm works with needles of different shapes and sizes and can be used to detect other medical tools imaged under the US.
When the experimental results are evaluated, the accuracy and the speed of the proposed method is sufficient to be used in actual percutaneous procedures in real time. The RMS error of the optical tracking system in position and the angle measurements were 1 mm and 1 • , respectively. The accuracy of the proposed method is similar to the external sensor's accuracy, which shows that the proposed needle localization method works reliably.