Tip Estimation Method in Phantoms for Curved Needle Using 2D Transverse Ultrasound Images

Abstract

Flexible needles have been widely used in minimally invasive surgeries, especially in percutaneous interventions. Among the interventions, tip position of the curved needle is very important, since it directly affects the success of the surgeries. In this paper, we present a method to estimate the tip position of a long-curved needle by using 2D transverse ultrasound images from a robotic ultrasound system. Ultrasound is first used to detect the cross section of long-flexible needle. A new imaging approach is proposed based on the selection of numbers of pixels with a higher gray level, which can directly remove the lower gray level to highlight the needle. After that, the needle shape tracking method is proposed by combining the image processing with the Kalman filter by using 3D needle positions, which develop a robust needle tracking procedure from 1 mm to 8 mm scan intervals. Shape reconstruction is then achieved using the curve fitting method. Finally, the needle tip position is estimated based on the curve fitting result. Experimental results showed that the estimation error of tip position is less than 1 mm within 4 mm scan intervals. The advantage of the proposed method is that the shape and tip position can be estimated through scanning the needle’s cross sections at intervals along the direction of needle insertion without detecting the tip.

1. Introduction

With the help of the beveled-tip needle, percutaneous interventions and therapies have been widely involved in current clinical procedures, such as brachytherapy [1,2], tissue biopsy [3,4], and drug delivery [5,6]. In intervention procedures, less needle misplacement will lead to a more reliable treatment and a more accurate medical practice. According to the clinical studies [7,8], the needle is easy to be deflected, which will cause needle tip misplacement and may lead to unsafe procedures. Due to the needle-tissue interaction, improper insertion force or physiological motions, such as breathing may cause targets or obstacles to be unstable, which will lead to an unexpected error. To address the challenge, real-time feedback is highly required. Usually, medical imaging devices are used, such as ultrasound (US) [9], computerized tomography (CT) [10,11], or magnetic resonance imaging (MRI) [2,12]. Generally, the image-guided percutaneous interventions are conducted with the use of CT or MRI. However, ultrasound-guided procedures are more attractive due to their advantages such as none ionizing radiations and real-time detection.

Many studies for the guidance of the needle during the insertion operation have been conducted with the help of ultrasound devices, and 2D ultrasound images are quite general to use, especially for the sagittal one (shown in Figure 1). Elif et al. proposed to use circular Hough transform to locate the needle tip accurately, even when the imaging is out-of-plane [13]. Kaya et al. localized the needle axis and estimated the needle tip by using a Gabor Filter in sagittal US images [14]. To execute in real-time, they improved the processing time by applying the bin packing method [15]. Recently, a template-based tracking method with the efficient second-order minimization optimization method has been used to track the needle [16]. In recent studies, more and more novel ideas have been used to locate the needle and evaluate its tip to sagittal US images, such as the use of signal attenuation maps [17], convolution neural networks (CNN) [18], and maximum likelihood estimation sample consensus (MLESAC) method [19]. However, a demerit of using sagittal US images is that out-of-plane bending of the needle cannot be detected. Therefore, the methods applied on sagittal US images are not suitable for the needle which may be deflected by the inevitable factors, especially for long needles.

An alternate solution for this problem is to use a 3D US image, which has been widely studied in recent researches. Yue et al. used a RANSAC method to detect the needle in a 3D US situation and Kalman filter has been used to reduce the error [20]. Chatelain et al. used the particle filtering to track a robot-guided flexible needle by using 3D US [21]. In addition, a convolutional neural network with conventional image processing techniques has also been used to track and detect the needle [22] and a naive Bayesian classifier was used to localize the needle among 3D US volume voxels [23]. However, the large 3D US volumetric dataset would make it difficult to obtain and process the online data.

Due to the above disadvantages, sagittal US images and 3D US volume are not suitable for a long flexible needle. To locate the needle accurately, methods that use transverse US images (shown in Figure 1) have been used successfully in some studies. For example, Vrooijink et al. [24] present a method to track the flexible needle during the insertion into a gelatin tissue by using 2D US images perpendicular to its needle tip. However, the background is pure without noise, which makes it impractical. Waine et al. [25,26,27] focus on the research about the needle insertion, in permanent prostate brachytherapy (PPB), where needles are typically 200 mm and easily to be deflected, indicating the fact that the rectum limits the movement of US probe. As a result, it is hard to acquire the sagittal images of the curved needle to observe the deflection during needle insertion. This is because the sagittal method has a strong relationship to the movement of the US probe when the needle is out of the view-field of the US images, and this movement maybe deforms the prostate as well as affects the needle target. For a deflected needle, the transverse US image is a better choice for its detection. Compared to the sagittal images, the transverse US images are easy to be acquired when the US probe scanning along the needle, no matter how much the needle is curved.

In this paper, we present a method to track a long-curved needle from the 2D transverse US images and estimate its tip for the guidance of needle insertion. Ultrasound is first used to detect the cross-sections of the long-flexible needle (shown in Figure 2 STEP 1). The needle shape tracking method combined needle detection with Kalman filter develops an accurate location and a robust tracking procedure with scan intervals from 1 mm to 8 mm (shown in Figure 2 STEP 2). Unlike the previous study [26], the 3D needle positions obtained from 2D US images and optical tracking systems have been used in KF for the precise location. The curve fitting method is then used to achieve the shape reconstruction and the needle tip position is estimated based on its length and the curve fitting result (shown in Figure 2 STEP 3). The advantage of the proposed method is that the shape and tip position can be estimated through scanning the needle’s cross-sections at intervals along the direction of needle insertion without detecting the tip. Besides, a novel histogram method is introduced to detect the needle in image processing, which can improve the needle localization under the effect of needle comet tail and the poor reflection, despite of the abrupt intensity changes. In addition, a robotic ultrasound system (RUS) [28] is built to evaluate the proposed needle tip estimation method. Results showed that the estimation of the tip position is less than 1 mm with 4 mm scan intervals.

The rest of this paper is organized as follows. The proposed methods will be introduced in detail in Section 2. Section 3 intends to represent the experimental setup and the results. Finally, the discussion and conclusions are detailed in Section 4.

2. Materials and Methods

The proposed needle tip estimation method in successive transverse US images can be divided into three stages: needle detection, needle shape tracking, and tip estimation. The processing diagram are shown in the Figure 3. The needle location will be manually selected as an initial region of interest (ROI) by the binary method in the first US image. After that, the prediction of needle position from a Kalman filter can be transformed in the transverse US images. At the same time, the next needle position in this ROI can be found through the histogram method. The KF is then updated for the current precise needle position and prediction of the next position. After all the cross-sections of the needle have been collected, the needle shape can be fitted by the curve fitting method (part C in Figure 3). Finally, the position of the needle tip can be estimated from the curve fitting result based on the length of the needle. In this study, the ROI is set as a square window with a length of three times larger than the needle diameter and its center represents the needle position. Needle detection (part B in Figure 3) is mainly about image processing, while the Kalman filter is used for needle shape tracking (part A in Figure 3).

2.1. Needle Detection

Needle detection is mainly for identifying the cross-section of the needle in the US images by using the binary method and the histogram method. As the ultrasound is quite sensitive to the metal, the needle-inside area can be brighter than others. A binary method [26] is first used to select ROI in the first image and then a histogram method is used to locate the needle despite the US shadows and poor reflection.

Binary method intends to strengthen the contrast of brightness to highlight the brighter area to select them. This method is used for the initialization which supposes to find the candidates in the first image. It includes intensity normalization, background reinforcement, and brightness enhancement. The center of the area is the location of the needle. During the experiment, a histogram method is proposed to find the needle accurately. The histogram method contains intensity normalization and background reinforcement. The histogram method tends to find an area of high-intensity pixels, which intends to find the upper face of the needle and then locate the needle based on the diameter. The background reinforcement part can be described as follows:

(1)$\begin{array}{cc}\hfill \min & I_t\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \sum _{I_t}^{255}{n}_{I_j}\le M\delta \hfill \\ \hfill & I_t\in [0,255]\hfill \end{array}$

where $I_j$ is the gray level of the pixel and $n_{I_j}$ is the number of pixels which have the gray level of $I_j$. M is the size of ROI, and $\delta$ is the manually selected parameter to limit the bright pixels. In this work, M is a square of $45\times 45$, and $\delta$ is set to 0.08 based on empirical tests.

There are two unexpected situations that may affect the position accuracy, namely the comet tail and the poor reflection. The comet tail will affect the size of the needle area and usually lead to a larger area than the actual size (Figure 4a). On the contrary, the poor reflection makes the needle area look much smaller in the image (Figure 4b). Therefore, the accurate location should be intended to eliminate the effect of shadows and poor reflection. In Figure 4, there are two examples which are used the two methods relatively. As the example shown in Figure 4a, the noise (yellow circle in the histogram of ROI) may probably be concerned as the needle while the needle is just concerned about a few pixels (red circle in histogram of ROI). The two methods can both filter the noise and locate the needle correctly. However, in ROI configuration (shown in Figure 3), the histogram method intends to find more candidates than the binary method, since the former focuses on the higher intensity pixels while the latter focuses on the area and intensity. Therefore, the binary method is more feasible in ROI configuration.

However, in the case of poor reflection, the needle displays a little in the image, and the area of the needle is smaller than the expected. Because the needle would reflect as long as the image gain is high enough or the sound power is big enough, it would reveal apparently compared to its surroundings. An example of the histogram of ROI is shown in Figure 4b, the red circle represents the upper surface of the needle. Moreover, the ROI square is darker than the one in Figure 4a, while the settings of the ultrasound are the same in Figure 4. From the figure, the histogram method seems to be more accurate than the binary method. In fact, the error of two methods, in this case, is 0.34 mm with 1.2 mm diameter of the needle. It is not that obvious to judge the accuracy. In this study, we use the histogram method during the experiment.

2.2. Needle Shape Tracking

As indicated in previous researches [20,29,30], Kalman filter has been successfully used for tracking needle in the successive ultrasound images. In this study, the Kalman filter is used to improve the estimation of the needle location in successive frames. As shown in Figure 5, the applied Kalman filter has two processes, prediction and update. The prediction stage intends to locate the needle position previously and set the ROI (red and yellow square in Figure 5) to find the needle precisely with a small window, which is supposed to reduce the computation. The update stage is the result of the needle position after the measurement from the histogram method.

The state prediction ${\widehat{t}}_i$ intends to represent the prediction state of the transverse needle center position $(x,y,z)$ in the reference frame with the change of the needle position $(▵x,▵y,▵z)$ at sample i according to the state t. $(▵x,▵y,▵z)$ are the difference between the previous needle position and the current needle position. $t_i$ is the result from the previous iteration, which is as follows:

(2)$t_i=\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ ▵x_i\\ ▵y_i\\ ▵z_i\end{array}\right]$

where $▵x_1$ and $▵y_1$ are set to be 0, while $▵z_1$ is equal to the scan interval. The prediction equations are as follows:

(3)${\widehat{t}}_i=A{t}_{i-1},$

(4)${\widehat{P}}_i=A{P}_{i-1}{A}^T+Q.$

The measurement update equations are as follows:

(5)$K_i={\widehat{P}}_i{H}^T{(H{\widehat{P}}_i{H}^T+R)}^{-1},$

(6)$t_i={\widehat{t}}_i+K_i(m_i-H{\widehat{t}}_i),$

(7)$P_i=(I-K_{i}H){\widehat{P}}_i,$

where A, H, R, and Q are as follows:

(8)$A=\left[\begin{array}{c}100100\\ 010010\\ 001001\\ 000100\\ 000010\\ 000001\end{array}\right],$

(9)$H=\left[\begin{array}{c}{I}_{6\times 6}\end{array}\right],R=Q={10}^{-6}\times \left[\begin{array}{c}{I}_{6\times 6}\end{array}\right].$

A is the state transition matrix, H is the measurement matrix. ${\widehat{P}}_i$ and $P_i$ are the priori and posteriori estimate error covariance, and R and Q are the measurement error covariance and processing error covariance, respectively. $K_i$ is the Kalman gain at sample i. $m_i$ is the measurement state from the needle detection.

The 3D prediction position $(x_{prediction},y_{prediction},z_{prediction})$ is obtained from the previous state. Before needle detection in the current US image, the 3D prediction position should be transformed on the image plane frame as 2D position $(I_{x\_prediction},I_{y\_prediction})$. After the update, the needle position $(I_{x\_update},I_{y\_update})$ in the image will be transformed into the 3D position $(x_m,y_m,z_m)$. Meanwhile, we can get $(▵x_m,▵y_m,▵z_m)$ from:

(10)$\begin{array}{c}\hfill (▵x_m,▵y_m,▵z_m)=(x_m,y_m,z_m)-(x_{previous},y_{previous},z_{previous}).\end{array}$

As a result, the measurement state $m_i$ in this frame is acquired as $(x_m,y_m,z_m,▵x_m,▵y_m,▵z_m)$. Through the measurement update, the current state can be obtained as $(x_c,y_c,z_c,▵x_c,▵y_c,▵z_c)$ from Equation (6). The relationship between these transformations will be described in the next subsection. In the previous study [26], the data from the image has only two dimensions, lacking the data from the direction along the movement of the probe, which leads to an incomplete location. Moreover, the space information is more capable to locate the needle accurately than the plane information. Therefore, 3D positions have been used in the KF for the precise location. The KF in this work is not only used for filtering, but also for predicting the next needle position in the US image. The ROI for the next iteration is centered around the needle position of the Kalman Filtering prediction, which can help to remove the outliers from the ROI.

2.3. Tip Estimation

Before tip estimation, 2D points should be transformed into 3D points based on the relationship of each frame. The relationship among the reference frame, probe frame, marker frame, and image frame are specified in Figure 6.

As shown in Figure 6, the image has one plane with 2 axes (axis x and axis y) and axis z is vertical to the image. Moreover, the probe has the same frame with image, except that the probe frame is designed in millimeters and the image frame is set in pixels. Equation (9) implies the transformation from image to reference:

(11)${Point}_1=T_{marker}^{ref}\times T_{probe}^{marker}\times T_{image}^{probe}\times {Point}_2$

where ${Point}_1$ and ${Point}_2$ are the points on the reference frame and image frame, respectively, $T_{marker}^{ref}$ is the transformation from the reference frame to the marker frame, $T_{probe}^{marker}$ is the transformation from the marker frame to the probe frame, $T_{image}^{probe}$ is the transformation from the probe frame to the image frame. Through this transformation, the needle position in the image can be directly re-defined in the reference frame for the needle tracking and curve fitting.

The tip estimation has two steps: curve fitting and tip estimation. In this study, the third-order curve line is used for the shape construction where the equations can be written as:

(12)$f(x)=\sum _{k=0}^3{a}_k{x}^k,$

(13)$g(x)=\sum _{k=0}^3{b}_k{x}^k,$

where $f(x)$ and $g(x)$ are the equations to fit the line along the x, which is the axis with the same direction of the insertion. $(a_0,a_1,a_2,a_3)$ and $(b_0,b_1,b_2,b_3)$ are the free parameters of the needle shape model.

After sample points of the inflected needle have been obtained, the least-square curve fitting method will be used to fit these points as a cubic line. The target functions to fit the cubic line can be defined as follows:

(14)$F(a_0,a_1,a_2,a_3)=min\sum _{i=1}^n{(f(x_i)-y_i)}^2=min\sum _{i=1}^n(\sum _{k=0}^3{a}_k{x}_i^k-y_i){}^2,$

(15)$G(b_0,b_1,b_2,b_3)=min\sum _{i=1}^n{(g(x_i)-z_i)}^2=min\sum _{i=1}^n(\sum _{k=0}^3{b}_k{x}_i^k-z_i){}^2,$

where n is the number of the points $(n\ge 4)$ and $(x_i,y_i,z_i)$ is the position of point i. By applying the $l_2$ norm minimization in the two-dimensional Euclidean space, it can be formulated as:

(16)$\underset{t\in R}{argmin}{\parallel X\mathit{a}-Y\parallel }_2,$

where $\mathit{a}={[a_0,a_1,a_2,a_3]}^{\prime }$, $Y={[y_1,y_2,\dots ,y_n]}^{\prime }$ and X can be written as:

(17)$X=\left[\begin{array}{cccc}1& x_1& {x_1}^2& {x_1}^3\\ 1& x_2& {x_2}^2& {x_2}^3\\ ⋮& ⋮& ⋮& ⋮\\ 1& x_n& {x_n}^2& {x_n}^3\end{array}\right].$

The solution can be estimated as follows:

(18)$\mathit{a}={(X^{T}X)}^{-1}{X}^{T}Y.$

From this solution, $f(x)$ and $g(x)$ can be obtained to construct needle shape. The tip position can then be estimated by the following optimum solutions based on the length of the needle:

(19)$\begin{array}{cc}\hfill \min & \parallel {\int }_{{tail}_x}^{{tip}_x}\sqrt{1+{f^{\prime }(x)}^2+{g^{\prime }(x)}^2}{dx-L\parallel }_2\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {tip}_x>{tail}_x\hfill \end{array},$

where ${tail}_x$ is the measured value of the tail position from the optical tracker, ${tip}_x$ is the expected value of the tip position in axis x, L is the length of the needle.

3. Results

3.1. Experimental Platform Setup

To verify the proposed tip tracking and shape sensing method, a robotic ultrasound system has been built, which includes a KUKA IIWA robot arm, a Wisonic ultrasound scanner, an NDI optical tracker, an NDI electromagnetic (EM) tracker, and a computer. As shown in Figure 7, the US probe is mounted on the effector of the robot arm by the gripper attached to the passive marker. The phantom or ex-vivo (like chicken breast in Figure 7) is punctured by an 18G beveled-tip needle with 200 mm long, while the needle tip is completely exposed for validation. The diameter of the needle is 15 pixels in the image. The NDI optical tracker is used to localize the marker bound with the probe, while NDI electromagnetic tracker is used to validate the tip position.

Experiments have been taken in a water tank, which provides a liquid environment for the ultrasound. And the needle is placed in water or inserted in the silica gel phantom (shown in Figure 8a), pork and chicken breast. The depth of the ultrasound is set to 4 cm. In this study, the needle is usually detected in 1 to 3 cm from the US probe. During the experiment, the robot arm automatically moves with the ultrasound probe along the direction of the needle insertion without any contact with the tissue (shown in Figure 8b). The whole scan length is at most 160 mm which depends on the scan intervals (shown in Table 1). The scan interval decreases with the increasing collected points.

Before data collection, the US image and the marker need to be calibrated. After that, the experiment starts after the needle finished puncturing manually. The robot arm is used to move the US probe scanning along the needle. Meanwhile, pose data are collected from the optical tracker and US images from the ultrasound scanner. Finally, the tail of the needle and its tip are measured by the optical and electromagnetic sensors, respectively, for the curve fitting and the tip validation. The needle is inserted manually, imitating the real situation of percutaneous intervention.

3.2. Tip Estimation

Four kinds of platforms have been used in the experiments: water, phantom, chicken, and pork. Each platform was tested several times. Figure 9 shows one test in chicken breast. In this case, the US probe moved along the needle in the chicken breast every 4 mm. The black square point on the left is the needle tail position and the blue line is the estimated needle shape. The green points are the detected needle and considered as the center of the needle, the blue point is the estimated tip position and the red point is measured tip position from EM. The estimation error is 0.69 mm in this test.

The error of the algorithm is shown in Table 2, which suggests that the errors increase with the increase of scan intervals. Figure 10 shows the results of the experiment on four platforms. The mean errors are all under 0.4 mm with a 1 mm scan interval in the four experiments, while the error is around 1 mm with an 8 mm interval.

4. Discussion and Conclusions

Needle insertion guided by ultrasound images is widely used for percutaneous interventions. However, the needle detection due to its deflection by the inevitable factors is challenging during the needle insertion. Such factors include needle-tissue interaction, improper insertion force, physiological motions, and so on. Automatic needle detection with needle tracking in 2D transverse US images could overcome these limitations and estimate needle tip through the curve-fitting method. The target of this study is to develop a robust needle detection and tracking method with the help of ultrasound images to estimate the needle tip precisely and accurately. We used a histogram method to detect the needle in transverse US images to decrease the effects of comet tail and poor reflection. In subsequent post-processing, the needle was tracked by the Kalman filter tracking method in consecutive US images with the help of the displacement of the probe. A third-order curve fitting method has been used to estimate the needle tip. When the probe is moved by the robot arm, the scanning time is different. We assume that the time when the probe stops to collect the data is the same. The less scan interval we choose, the more collection points we can obtain and the more scanning time it takes. Therefore, the scanning time mainly depends on the number of scanning points while the accuracy lies in how short the scan interval is. In other words, the accuracy of the tip estimation can be improved by reducing the scan interval to collect more needle positions. However, this will consume more scanning time and reduce the efficiency of tracking. Inversely, fewer collecting positions would cost less time but may lead to a more possibility of the failure of shape construction and a more possibility of large error of the tip estimation. As a result, how to balance precision and scanning time is very important to make the proposed method more efficient. In our experiment, it is found that a 4 mm scan interval has an error less than 1 mm, which is a better choice to satisfy both requirements.

In the proposed method, needle shape tracking has a great contribution to the accurate localization, since the needle can be tracked precisely by Kalman filter through its prediction and update. However, needle shape tracking is heavily dependent on the scan interval, especially for a large curved needle, since Kalman filter is generally well functioned in the lineal system. As a result, if the needle is deflected during insertion, the Kalman filter would make mistakes and wrongly predict the needle position when the scan interval is large. In this study, it is found that Kalman filter would fail if the scan interval is more than 8 mm. This may due to the impact on the prediction of KF with a large deviation. Moreover, the deviation will not be eliminated even with the change of the ROI size. However, the Histogram method showed the accurate and effective detection of the needle, but it relies on the brightness of the image as the needle could not be easily detected where there are plenty of pixels with the highest intensity (which has a max value of 255). However, this condition can be controlled by the setting of the ultrasound scanner to expand the gray level of the image properly.

The proposed method still has its limitations. During the experiment, time is needed for the data collection from the US scanner and optical tracking system, and the movement of the robot arm, which is affected by the scan length and scan interval. However, it is very hard to acquire the whole position of a long needle in one scan for any medical image sensor. Therefore, in the future, it is valuable to find a method to reduce the times of needle detection in order to reduce the time for the tip estimation. Moreover, when the robot arm moves with the US probe precisely, the tissue and needle have a possibility to be deformed by the probe motion. Hence, it is necessary to make the robot arm move smoothly as well as correctly on the surface of tissue. In addition, patient motion is the biggest uncertain problem, which leads to the failure of needle insertion and detection.

In the proposed system, we use the 2D US scanner for the detection of the needle in various kinds of tissue. However, the 3D US scanner can also be used in this system. Although it has volume data and the detection method is different, the tracking method is able to be similar, as we also use the 3D positions for tracking in this study. Moreover, time is also needed for data collection and the movement of the robot arm, especially for the long needle, which is easily out of view-field of US volumes or images. Therefore, we use 2D US images in this study.

In this paper, a method for tracking a long-curved needle from the 2D transverse US images and the tip estimation is represented and demonstrated with RUS. Ultrasound is first used to detect the cross section, with the probe moving along a long-flexible needle. Needle shape tracking method combined needle detection with Kalman filter by using 3D needle positions develops an accurate location and a robust tracking procedure. Needle shape is then constructed by using the curve fitting method and its tip position is estimated based on the former result. A histogram method is introduced to detect the needle in image processing, which can improve the needle localization despite the abrupt intensity changes. This new imaging approach is proposed based on the selection of numbers of pixels with a higher gray level, which can directly remove the lower gray level to highlight the needle. Results of the experiments suggest that the detection of the needle by the histogram method and Kalman Filter has high precision with minimum error 0.13 mm with a 1 mm scan interval in the phantom experiment and maximum error 1.35 mm with a 8 mm scan interval in the pork experiment. With the increase of the scan interval, the mean error would rise. Moreover, results showed that the estimation of the tip position is less than 1 mm within 4 mm scan intervals. We suggest choosing a 4 mm scan interval to balance the precision and scanning time to maximize efficiency. In the future, we would make the experiments of how long the scan length is the best length to estimate the needle tip. The proposed method would be a great assist to surgeons to locate the needle tip when they perform percutaneous insertion procedures with a long flexible needle, such as prostate brachytherapy.

Author Contributions

Conceptualization, M.Q.-H.M.; validation, Z.L.; writing—original draft, Z.L.; writing—review and editing, S.S. and L.L.

Funding

This research was funded in part by the National Natural Science Foundation of China, grant number 61803123, and in part by the National Key R&D Program of China, grant number 2018YFB1307700, and in part by the Natural Science Foundation of Guangdong Province, China, grant number 2018A030310565.

Acknowledgments

For this kind of study, no formal ethics approval is required by the institutional ethics com. Informed consent was obtained from all individual participants included in the study.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

KF	Kalman filter
US	Ultrasound
RUS	Robot ultrasound system
CT	Computed tomography
MRI	Magnetic resonance imaging
RANSAC	Random sample consensus

Figures and Tables

Figure 1. Two methods by using 2D ultrasound for detection.

Figure 2. The proposed tip estimation method. STEP 1: 2D transverse US images with needle cross-sections are collected by using RUS; STEP 2: Needle cross-sections are detected and tracked in the successive US images; STEP 3: Needle shape is constructed and its tip is estimated.

Figure 3. The pipeline of curved needle tip estimation. As the ROI configuration finished, needle tracking with needle detection begins step by step. Part A shows the needle tracking by KF, while Part B shows the needle detection to locate the needle tip. Part C shows the shape construction and tip estimation.

Figure 4. Two cases may generate the errors: (a) the comet tail of needle with binary method process, histogram method process and the histogram of ROI; (b) the poor reflection of the needle with histogram method process and the histogram of ROI.

Figure 5. The two steps of the Kalman filter. As the next US image is acquired, the previous state (xprevious,yprevious,zprevious,▵xprevious,▵yprevious,▵zprevious) is used to predict the needle position (xprediction,yprediction,zprediction, which then will be transformed to (Ix_prediction,Iy_prediction) in the image as the center of ROI. The yellow square is the ROI corresponding to (Ix_prediction,Iy_prediction) and the red one is the update step in KF by using the measurement data from needle detection to locate the needle with its center (Ix_update,Iy_update) as the needle position. Finally, the measurement state (xm,ym,zm,▵xm,▵ym,▵zm) and the current state (xc,yc,zc,▵xc,▵yc,▵zc) can be obtained.

Figure 6. The relationship of the frames.

Figure 7. The devices of the experiment.

Figure 8. The phantom used in the experiment and the movement of the probe.

Figure 9. Experiment in chicken breast with a 4 mm scan interval.

Figure 10. The experiment with different scan intervals on the four platforms.

Table 1

Scan lengths with different scan intervals.

Scan Interval	Scan Length	Points
1 mm	160 mm	160
2 mm	159 mm	80
3 mm	160 mm	54
4 mm	157 mm	40
5 mm	156 mm	32
6 mm	157 mm	27
7 mm	155 mm	23
8 mm	153 mm	20

Table 2

The results of tip estimation (mm).

Accuracy
Intervals	Water	Phantom	Pork	Chicken
1 mm	$0.32\pm 0.10$	$0.33\pm 0.15$	$0.37\pm 0.07$	$0.29\pm 0.09$
2 mm	$0.36\pm 0.21$	$0.41\pm 0.16$	$0.31\pm 0.05$	$0.31\pm 0.12$
3 mm	$0.45\pm 0.16$	$0.46\pm 0.18$	$0.56\pm 0.08$	$0.33\pm 0.09$
4 mm	$0.55\pm 0.16$	$0.59\pm 0.31$	$0.68\pm 0.14$	$0.38\pm 0.21$
5 mm	$0.60\pm 0.13$	$0.59\pm 0.29$	$0.90\pm 0.20$	$0.40\pm 0.14$
6 mm	$0.69\pm 0.13$	$0.48\pm 0.09$	$0.99\pm 0.13$	$0.55\pm 0.11$
7 mm	$0.83\pm 0.17$	$0.62\pm 0.14$	$0.84\pm 0.30$	$0.52\pm 0.17$
8 mm	$0.95\pm 0.11$	$0.73\pm 0.26$	$1.06\pm 0.18$	$0.81\pm 0.19$