Automatic needle detection using improved random sample consensus in CT image‐guided lung interstitial brachytherapy

Abstract Purpose To develop a method for automatically detecting needles from CT images, which can be used in image‐guided lung interstitial brachytherapy to assist needle placement assessment and dose distribution optimization. Material and Methods Based on the preview model parameters evaluation, local optimization combining local random sample consensus, and principal component analysis, the needle shaft was detected quickly, accurately, and robustly through the modified random sample consensus algorithm. By tracing intensities along the axis, the needle tip was determined. Furthermore, multineedles in a single slice were segmented at once using successive inliers deletion. Results The simulation data show that the segmentation efficiency is much higher than the original random sample consensus and yet maintains a stable submillimeter accuracy. Experiments with physical phantom demonstrate that the segmentation accuracy of described algorithm depends on the needle insertion depth into the CT image. Application to permanent lung brachytherapy image is also validated, where manual segmentation is the counterparts of the estimated needle shape. Conclusions From the results, the mean errors in determining needle orientation and endpoint are regulated within 2° and 1 mm, respectively. The average segmentation time is 0.238 s per needle.

deformation and needle deviation. 8 Suboptimal dose distribution leads to decreased biochemical control and high risk of damage to surrounding tissues. An image-based automatic needle detection algorithm makes it possible to perform efficient intraoperative dose correction based on the actual needle paths. In the treatment planning system, the planned seeds can be imported into the corresponding detected needle paths for dose calculation. If any dose deficiency and hot spot are identified within the lung target, position adjustment of seeds and additional needles will be done. However, needle artifacts, image attenuation, signal dropouts, and low contrast between the needle and the ribs have brought many obstacles to automatic needle detection in CT data, as shown in Fig. 1. Therefore, accurately and quickly identifying the needle from a large amount complex of CT data has become the nodus of this technique.
A wide variety of research studies have been performed on automatic long straight surgical instruments segmentation from two-dimensional (2D) and three-dimensional (3D) ultrasound (US) images.
Radon, Hough transform (HT), and parallel projections are three common detection algorithms. An example that using a form of the Radon transform, Barva et al. 9 proposed that the electrode axis can be found through maximizing parallel integral projection (PIP).
Although the accuracy of this method is high with tip error within 0.2-0.3 mm, it is not a pragmatic approach in a clinical environment due to a long processing time in MATLAB (tens of minutes). Additional PIP speedup achieved through harnessing the computational power of the GPU and using a fast search schema to track the surgical instruments real time. 10 Uhercik et al. 11 presented a multiresolution PIP method, achieving faster electrode axis localization in 3D US data than generalized PIP and yet maintain the same accuracy.
Qiu et al. 12 used the 3D HT-based method for locating the prostate therapy needle, which achieved endpoint localization accuracy of 1 mm in vivo patient images and the mean time cost is 2 s. By modifying the HT algorithm with coarse-fine search strategy and a four parameter representation of lines in 3D space, Zhou et al. 13 segmented needles quickly, accurately, and robustly. Quick randomized Hough transform (RHT) reduced the computational effort to 1 s in a 35 MB data by doing the RHT only on coarse resolution volumes. 14 Okazawa et al. 15 generalized the HT to address needle curvature in 2D ultrasound images. Ding et al. 16 used two orthogonal image projections to segment needles from 3D US images and a similar segment algorithm can be used to detect curved needle found by Aboofazeli et al. 17 For higher preciseness, these methods of projection need images with low noise and good contrast. Regarding localizing the needle position and estimate the endpoint in real time, Yan et al. 18 presented the shape information and level set technique, which had been validated in prostate phantoms.
The general disadvantages of aforementioned methods are their high sensitivity to interference, computational complexity, and focusing on single needle separation. Hence, a more robust and computational efficient approach that functions well within a complicated environment is still required. In this paper, a modified needle detection algorithm has been proposed and evaluated, which estimates the entire needle shape utilizing image preprocessing and an improved random sample consensus (RANSAC) algorithm. RANSAC has been successfully used in surgical tool localization from 3D US images 19 and multiple transverse US images, 20 but so far has not been applied to CT transverse data. In this algorithm, the iteration process has been accelerated through pretesting technology to improve the overall needle segmentation efficiency. Meanwhile, the accuracy is guaranteed by combining locally optimized RANSAC 21 with principal component analysis (PCA). 22 Furthermore, by cycling out of the current best model inliers, all needles in a single transverse image are detected. Finally, the algorithm is comprehensively evaluated in simulation data, physical phantom, and in vivo CTguided lung brachytherapy images.

2.A | CT image preprocessing
The distance between adjacent needles is usually 5 or 10 mm, which has been proved with minimal damage to the human body. 5 The 1. An instance of the needle appearance during lung interstitial brachytherapy performed on CT images. The ROI is displayed on the right. thickness of CT slices used for lung brachytherapy is normally 5 mm to ensure that all needles and seeds can be detected while obtaining distinct tissue imaging. Using a locating device and lasers to adjust and position the coplanar template before scanning can ensure that each needle is exactly on the CT slice. It is assumed that the acquired images contain all the needles, and the count of needles on each slice is known. The planned contours are manually adjusted by the physicist to account for the anatomic deformation, and then imported into intraoperative CT images. According to literature, 23 we found that the deformation of the 18-gauge brachyneedle occurs mainly in the direction perpendicular to the bevel of the needle and no bending occurs during insertion.

2.A.1 | Region of interest (ROI) extraction
Comparing with the manual extraction, a semiautomatic ROI detection method is used based on the modified target contour. First, a smallest enclosing rectangle of the target area is calculated through the pixels on the contour. Then, the rectangle is expanded by 3 mm to ensure all the needle information within the ROI. If any needle point candidates are located outside the enlarged rectangle, the expansion distance will be increased manually. In the end, the image surrounded by the updated rectangle is cropped for subsequent processing. A ROI for one of the clinical CT images is shown in Fig. 1.

2.A.2 | Contrast enhancement
In order to improve the contrast between the needle and the exterior background noise, an intensity mapping function is applied to the cropped image after localizing the search area. The function is In the transformation process, the original pixels with intensities smaller than i low_in are given the value of intensity i low_out . Likewise, pixels having the intensities higher than i high_in are assigned to i high_out . As for the parameter γ, which defines the shape of the convert curve. When γ <1, the pixel intensities of the low-intensity are expanded and compressing those of the high-intensity pixels. Contrarily, when γ >1, the low-intensity pixels are compressed and the high-intensity pixels are expanded.
In our algorithm, the values i low_in and i high_in are chosen to 15% and 100% of the maximum intensity value 1, and the desired spectrum is set to range 0,1 ½ . By setting γ >1, additional expansion toward higher intensity pixels, allowing for enhancing contrast between the needle and external background noise. Fig. 2b,c present the comparison between the original cropped image and the contrast enhanced image.

2.A.3 | Thresholding process
After process of contrast enhancement, the histogram of the gray distribution in the ROI is obtained. Based on thresholding α, intensities IðxÞ∈ 0,1 ½ is separated into two disjoint sets: X n (needle pixels) and X b (background pixels).
where X is a series of pixel coordinate. Threshold α is identified based on empirical results, because the needle grayscale distribution of each patient case is different. Through Eq. (2), pixels with intensity larger than α are identified as candidate needle points. Note that X n just provides an approximate estimation of the candidate needle points and contains some false points. A modified RANSAC method for needle axis identification is described in the following section. | 123 among them m is the minimum data collection for determining the parameterized model and its value is respond to the complexity of the geometry model. Once a model is presumed by calculation from the minimum dataset, entire data points in dataset X are used to examine the model and then determine the inliers of it. The sampling process will be repeated until meeting the final evaluation criteria. It means that the minimum samples required for a smallest subset without outliers can be screened at least on the given confidence probability P. Finally, the model with the highest support rate (most inliers) is the desired model. Supposing that the incorrect outlier occupancy rate is ɛ, and the probability of choosing an un-contaminated m inliers sample point is 1 À ɛ ð Þ m . Similarly, k samples are chosen and the probability of contamination of these sample data (at least one outlier) is ð1 À ð1 À ɛÞ m Þ k . Therefore, to make sure that the probability of confidence is above P, the least number of samples k must satisfy the following relationship: Through aforementioned step, the time complexity of RANSAC algorithm is considered to T: where T X is the time of each random sampling, T C is the average time of validating the model for each data point; N is total number of data points in X. In general, T X and T C are considered invariant for a specific issue, so the time complexity of RANSAC algorithm is decided by k and N.
From Eq. (3), we can see that k is exponential with ɛ and m.
Under a certain confidence probability, the value of k and N will heighten with the increasing complexity of the data set and the model. Meanwhile, the time intricacy T will be increased accordingly.
Hence, the technology of combination of preview model parameters evaluation and local optimization is proposed to improve the speed of the algorithm with a high accuracy.

2.B.2 | Preview model parameters evaluation technology
Assuming that the sample data for pretesting is n and the probability of existence at least n f inliers is P t . In other words, at the confidence probability P t , the inliers are at least n f for correct model under n tests. If the number of inliers is less than n f , this model's parameters are contaminated by outliers, which no need to participate in the subsequent verification of all the data. For the correct model parameters, the probability of passing the pretesting is: where C i n is a combination of i data selected from sample data n. In the actual calculation, the maximum n f and the corresponding P t , that satisfy the condition are calculated by testing different n f under the minimal limit of P t . According to the standard deviation of the error, the point whether belonging to inliers is determined 25 and as shown in Eq. (6): when n ≥ 2m, the error standard deviation σ can be estimated as follows: where d i is the data error shown in Fig. 3, med i d i j j is the median of the absolute value of error.
After pretesting, the correct model might be regarded as wrong and removed. The probability of choosing an un-contaminated m inliers sample point is changed to ð1 À ɛÞ m P t , and then the Eq. (3) becomes: The impact of the pretesting process on the calculation results is transferred to the sampling process through the Eq. (8). When setting P t < 1, the corresponding samples k is more than that used in general RANSAC algorithm. However, the added model parameters are very few compared with the model that can be filtered out during the pretesting, since most of the model parameters were affected by outliers.

2.B.3 | Local optimization technique
As the criteria for stopping operation of algorithm are too idealized, a contaminated model may also get high support (more inliers) and lead to an incorrect result. The combination of local RANSAC and PCA is proposed to tackle the accuracy problem. Through pretesting technology, an initial solution and a consensus set that match can be F I G . 3. Diagram illustrating needle and d i for the ith point in the dataset X.
acquired. The local RANSAC is sampling on this suboptimal inlier dataset and the new model is also tested on these data points.
Unlike the ordinary RANSAC algorithms, there is no need to minimize the number of samples because sampling is done at inliers. In general, to reduce the overall computation time, the number of samples k L in local RANSAC is set between 10 and 20 times.
After using local RANSAC, an optimized model parameters and dataset of inliers can be obtained. Nevertheless, the parameters are not the best solution because it is calculated from the smallest sample subset. Therefore, the PCA is applied to obtaining the desired model fit by minimizing the error estimation of the consensus dataset. The optimized approximation of the needle shape will be performed in next experiment section.

2.C | Needle tip localization
Once the needle axis has been identified, the intensities of all pixels along the axis are obtained and the endpoint of the needle can be found as a mutation point by the method of Barva et al. 9 An instance of intensity profile along the estimated axis is shown in Fig. 4. Given the intensity I, the priori estimated probability dis-

2.D | Multineedles detection
Note that the RANSAC program can only take one needle at a time.
To extract all the needles in a slice, we have adopted a method of successive deletion. It refers to labeling all inliers corresponding to the current optimal model after a modified RANSAC routine and removing these data from the overall sample, then using the updated sample dataset for next detection. As the sample data are updated, the outlier occupancy rate ϵ for current needle must also be changed: where n needles is the total needle number in current slice, ɛ noise is the noise occupancy rate after image processing, n en refers to the maximum inliers number of each needles, N 0 is the initial count of the sample dataset and N j is the points number after j times updating.
The successive deletion method is repeated until there is no subset in the total data that satisfies the best model conditions. A flowchart showing the complete needle shape estimation steps is shown in

2.E | Experimental methods
Three different experiments including simulation data, physical phantom images, and case images of lung interstitial brachytherapy were applied to test the performance of proposed approach.
In this study, the confidence probability P was set as 99%, the sample data n for pre-testing was chosen as 15 and the probability P t was set to 80%. According to empirical results, ɛ noise was selected as 15%. For each experimental data, the detection algorithm was

2.E.1 | Evaluation on simulation data
Datasets with a size of 1500 points were generated by the computer to test our method. All the inliers came from a given straight line in each dataset. Moreover, uniformly distributed different proportions of Gaussian noise varying between 0.1 and 0.8 was added with MATLAB function "randn" in a certain area. Note that the synthetic data generated were applied only to verify the feasibility and improvement of proposed method in the data with multiplicative Gaussian noise.

2.E.3 | Evaluation on in vivo experiment
Fourteen patient CT images of lung cancer brachytherapy were also used in our experiments. All the patients were recruited from the Second Hospital of Tianjin Medical University. Moreover, details of the needle provided by medical physicist are presented in Table 1.

3.A | Results on simulation data
The results of simulation data between the proposed RANSAC (P-RANSAC) and the general RANSAC (G-RANSAC) are shown in Fig. 7.

3.B | Results on physical phantom
For each needle insertion distance h investigated, the observed variation of β and ξ with different insertion angulation θ is plotted in     Table 1, it can be found that the maximum angular deviation is 1.243 AE 0.562°and maximum tip location error is 0.757 AE 0.218 mm in all identified needles. Comparing with the previous experimental data, the error of needle configuration in patient data becomes larger. This is mainly because the mistake is not only derived from the segmentation algorithm but also from the medical physicist.

3.C | Results on in vivo experiments
Additionally, the segmentation efficiency of the algorithm is also model in the first few needles extraction. In other words, the outlier ratio is higher in algorithm early period.
In our physical phantom experiment, the influence of different insertion angles and insertion distances on the performance of detection method are analyzed. It can be seen that the difference in accuracy of needle axis and endpoint segmentation is mainly caused by target amount of information, which depends on the insertion depth. Despite aforementioned factor, the algorithm we proposed still has higher accuracy and better robustness, since the error in determining the needle tip is less than 0.78 mm and the needle orientation is less than 1.35°. Comparing with previous related study in needle orientation, the experimental result illustrates more competitive than 2.3°(Cool's method). 26 Moreover, for the error of needle tip, the evaluation value is smaller than 0.8 mm (Ding's algorithm) 16 and 1.43 mm (Qiu's algorithm). 27 For the actual brachytherapy images, the needle tip location and orientation error are controlled within 1 mm and 2°for most cases.
The submillimeter accuracy is sufficient for intended brachytherapy clinical practice. Meanwhile, the whole failure rate is 6.75%, which is completely reasonable for a random iterative algorithm (RANSAC), and this value is more competitive than 33% (Uhercík's result) 19 and 16% (Qiu's result). 27 As for the segmentation efficient, the average processing cost by our approach is 0.238 s per needle. It is much shorter than other algorithms, the mean time in reference 11 is 7.3 s, 26 is 3-5 s, and 19 is 0.64 s. This time consumption is perfectly acceptable compared to a surgery with few hours. In addition, based on our knowledge, the application of needle detection based on CT images for brachytherapy has not yet been reported.
It is important for the best therapeutic effect in lung brachytherapy to realize intraoperative dose replanning before seeds implanted into pathological tissue. In this paper, a CT image-based quick and accurate needle reconstruction technology has been described, which provides a new idea for dose correction. Once identified all needles, the original plan (all the seeds) can be imported based on the correspondence between preoperative and intraoperative needle, and then ameliorate dose distribution by adjusting primary seeds position as well as adding extra needles and seeds. On the other hand, an image measurement tool in auto-determining the needle tip position and needle orientation is available for medical physicist. If obvious positioning mistakes are detected, the remedy like withdraw or reinsert the needle will be considered. In addition, the proposed method provides the flexibility to segment other parametric objects through redefining the model function and corresponding error term.
Although our algorithm can detect multiple needles precisely at once, there are still some failures. Thus, to guarantee a robust segmentation, an advanced initialization (preimage processing) is desired in our application. It is important to note that the needle deformation in Z direction was not considered in this article. There are a few reasons for current lung brachytherapy: from one perspective, the needle insertion depth will be decreased as much as possible by adjusting the patient position in surgery to increase the control of needle. From another perspective, changing the bevel orientation of the needle tip during needle insertion contributes to diminishing the deformation problem. Through these intraoperative methods, the deviation of needle position in Z direction is much less than the deviation in X and Y direction. In our future work, the possibility of combining with the coronal or sagittal slice will be considered to further tackle this problem and achieve a 3D parametric representation of the needle shape.

| CONCLUSION
An improved RANSAC algorithm has been presented in this paper for CT image-based needle segmentation. Preview model parameters evaluation, local optimization combining local RANSAC and PCA, and successive inlier deletion are the key technologies applied in this method, which were tested in simulation data, physical phantom, and brachytherapy case. Meanwhile, the promising experimental results with short calculation consumption and higher accuracy meet the requirement of clinical treatment. We expect that, with the image-based automatic needle placement, a more reliable dose optimization module will be integrated into the brachytherapy treatment planning system.

CONFLI CT OF INTEREST
The authors declare that they have no conflict of interest.

ETH ICAL APPROVAL
All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards.