Research Article Force-Sensor-Based Estimation of Needle Tip Deflection in Brachytherapy

A virtual sensor is developed for the online estimation of needle tip deflection during permanent interstitial brachytherapy needle insertion. Permanent interstitial brachytherapy is an effective, minimally invasive, and patient friendly cancer treatment procedure. The deflection of the needles used in the procedure, however, undermines the treatment efficiency and, therefore, needs to be minimized. Any feedback control technique to minimize the needle deflection will require feedback of this quantity, which is not easy to provide. The proposed virtual sensor for needle deflection incorporates a force/torque sensor, mounted at the base of the needle that always remains outside the patient. The measured forces/torques are used by a mathematical model, developed based on mechanical needle properties. The resulting estimation of tip deflection 
in real time during needle insertion is the main contribution of this paper. The proposed approach solely relies on the measured forces and torques without a need for any other invasive/noninvasive sensing devices. A few mechanical models have been introduced previously regarding the way the forces are composed along the needle during insertion; we will compare our model to those approaches in terms of accuracy. In order to conduct experiments to verify the deflection model, a custom-built, 2-DOF robotic system for needle insertion is developed and discussed. This system is a prototype of an intelligent, hand-held surgical assistant tool that incorporates the virtual sensor proposed in this paper.


Introduction
Permanent interstitial brachytherapy is a cancer treatment procedure, in which radioactive seeds are implanted in tissue (e.g., the prostate, see Figure 1) in order to eliminate the cancer from inside.This procedure has emerged as an effective, minimally invasive, patient friendly, and cost-effective treatment option.For maximum treatment efficiency, the seeds have to be placed at exact locations inside and around the tumor that are determined in the preoperative planning stage.The radioactive seeds are initially loaded inside special needles.Intraoperatively, the seed-carrying needles are manually advanced toward planned locations where the seeds are deposited.
Two critical assumptions in the above procedure are that the needles will remain parallel across the entire length of their insertion in tissue and that the tissue will not deform as the needles penetrate it.However, in practice, neither assumption holds well, causing the actual needle trajectories to not pass through the planned locations.In fact, current manual needle insertion techniques for prostate treatment can place seeds with an accuracy of only 5 mm, which is a substantial error given the average prostate size.As a result, due to delivery of a different radiation dose, the radiation might not have the desired effect on the cancerous tissue and could instead undermine healthy tissue.
The predominant causes of inaccuracy in seed placement are needle deflection and tissue deformation during needle insertion/retraction [1,2].In general, needle bending is a function of the needle geometry [3,4].For instance, needles experience more bending with smaller diameters and with beveled tips.Beveled tips (Figure 2) are needed for easily cutting into tissue [5].With a beveled tip, needle bending is larger for smaller bevel angles [6].Needle deflection and tissue deformation are coupled effects and influence each other.For simplicity, past studies have tried to decouple these effects by considering one of two extreme cases: a rigid needle in deformable tissue or a flexible needle in rigid tissue.For the second case, a bicycle-like kinematic model predicting the needle tip position was proposed [7].A result pertinent to this paper is that when a beveled-tip flexible needle is pushed through tissue, the asymmetry of the tip causes the needle to bend (see Figure 2(a)).This is intuitively understood because the asymmetric geometry of the beveled tip leads to a force on the tip in the same direction as the bevel and a displacement of tissue to that side [5].
To improve needle targeting accuracy despite the bevel effect, there is active research on high-resolution imaging, tissue modeling, and robotic insertion preplanning.However, the resulting systems are still too complex and costly for the operating room.There is currently a need for a simple, inexpensive, and effective instrument in the surgeons' hands that can enable them to place the seeds more accurately.
We are working toward a low-cost and effective handheld instrument for use in prostate needle insertion that automatically corrects the needle tip's bevel location during manual insertion in order to improve needle targeting and seed placement accuracy; the bevel location is defined as the angular position of the bevel tip relative to the needle's longitudinal axis (Figure 2(b)).A main idea related to needle steering is to leverage the very tendency of the bevel-tip needle for traveling on a curve to bring the needle tip back on the intended path once the needle has been bent, regardless of what caused the needle deflection [8].To this end, the needle's bevel location can be changed by a feedback control algorithm based on the real time feedback of needle deflection.The feedback control algorithm will steer the needle into the opposite direction when needle deflection exceeds a set threshold by simply turning the needle about its longitudinal axis by 180 ∘ .This causes the forces acting at the needle tip to point in the opposite direction and steer the needle back toward the unbent ideal path.Since the surgeon is fully in charge of inserting the needle and the computerized system only modifies the bevel location, the hand-held instrument will be intrinsically safe to use.
The contribution of this paper lies in real time measurement of needle deflection using our proposed novel virtual sensor; incorporating a physical sensor for measuring this quantity is next to impossible in a clinical setting at least with the current technology.To measure the needle deflection during insertion into soft tissue, one may consider using 3D ultrasound imaging.However, the detection of the typical needle tip deflections in the order of several millimeters will be difficult due to the low resolution of ultrasound images and the associated high computational demands, notwithstanding the high costs of integrating a 3D imager and the required image processors in a computerized needle insertion system.Using 2D ultrasound imaging is not suitable because deflection measurement requires that the needle be always visible in the image plane while in practice the surgeon needs to intermittently switch between various sagittal and axial planes to monitor the insertion progress.Another possibility to measure the needle deflection is to use an electromagnetic tracking system, such as the Aurora Electromagnetic Tracking System by NDI, which involves a sensor embedded in the needle.However, this is only appropriate for in vitro testing where invasiveness is not an issue; sterilization requirements and the preloading of seeds inside the needle make in vivo utilization of a tracker inside the needle prohibitive.Therefore, in this paper, we pursue a noninvasive virtual sensing approach to needle deflection measurement in soft-tissue needle insertion.The foundation of our proposed virtual sensor is the hypothesis that the transverse force and the bending moment at the needle base predict the needle's longitudinal deflection inside the tissue.Therefore, in this paper, as the needle is being inserted into tissue, the forces and moments reflected at the needle base are measured by a physical sensor.Then, a needle tip deflection model and virtual sensor are developed solely based on this force/moment at the needle base.This virtual needle deflection sensor can be applied to a wide range of percutaneous medical procedures that involve needles prone to deflection.

The Needle Deflection Model.
The developed needle model is fundamentally based on a static model for deflection of a beam.The needle can be regarded as a cantilever beam since its only fixation is at the base.The equations are based on the Euler-Bernoulli beam theory, which relates the loads or forces applied to the beam to its deflection.As the needle deflects increasingly during insertion in soft tissue, it acts as a loaded spring that tries to return to its initial unbent state but is kept in place by the tissue.The needle exerts a distributed load perpendicular to the needle axis onto the tissue.To keep the needle in its bent state, that is, to maintain the equilibrium condition, the tissue in return reacts with a distributed load along the needle (Figure 3).In general, distributed loads can be replaced by a resultant point force (concentrated load) that acts at a specific point along the needle.This reduction allows for a simplification of the equations used to calculate the deflection.The location of the resultant force is the geometric center or centroid of the area of the distributed load [9].The magnitude of the force acting on the needle base (  ) is the integral of the load and can be expressed as or for an equally distributed load, where  is the load per unit length,  0 represents length portion outside of tissue, and  is the length portion of the needle, which is inserted into tissue.It should be noted that the indicated force at the needle clamping is the resistance force exerted by the clamp onto the needle.While the measured force (  ) is equal in magnitude to   , the directions are opposed.The relation between   and   can be expressed as Since the only support of the needle is at its base, the full magnitude of the resulting force and moment as well as the direction are reflected at the needle base and can be Tissue Clamping Needle Centroid of the triangle measured.This means that the expected force measurement results should have the same direction as the deflection.Insertion experiments, however, showed a force at the needle base in opposite direction.Figure 8(a) shows forces and torques for a deflection in negative -direction.This means that the forces exerted by the tissue cannot only be one distributed load acting in negative -direction (cf. Figure 3) but must be comprised of multiple loads acting in positive and negative -direction.Tissue below the needle can be considered as acting as a support, and tissue above the needle still exerts a distributed reaction force on the needle.The triangular distribution of  1 in Figure 4(a) can be explained by assuming that tissue is pushed down more towards the side where the needle enters the tissue.As the needle penetrates deeper into tissue, the deflection increases, which leads to pressure being applied on the tissue in negative -direction.Tissue, which is horizontally closer to the point of insertion (here: left side, see Figure 4(a)), is exposed to this pressure for an increased period.Therefore, the distribution of the load applied on the needle by tissue is at its maximum at the insertion point and decreases in -direction.This assumption was also made by Abolhassani et al. [10].They however did not include it in their model, as it had very little impact on their case.For the load above the needle, a reverse effect applies.As the tissue closer to the needle entry point is pushed towards the direction of deflection increasingly, the tissue above the needle in this area cannot apply any resistance force from above.
As before mentioned, the loads  1 and  2 can be reduced to the concentrated loads  1 and  2 , which act at the centroid of the triangularly distributed load.The centroid of a triangle and its location is illustrated in Figure 4(b).The figure shows that the centroid of a triangularly distributed load is at /3, where  is the length of the triangle's leg.Thus, the factors  and  should have the values 1/3 and 2/3, respectively (see Figure 4(a)).It should be noted that the triangular distribution of loads  1 and  2 is an assumption.
A common method to finding deflections in beam structures is twice integrating the bending moment equation: where  and V are the moment and deflection at a distance  from the base and  and  are the Young's modulus of stainless steel (200 GPa) and the area moment of inertia of the needle.
The equation for the moment  can be found by analysis of the free body diagram.In our case there are multiple forces acting simultaneously ( 1 and  2 ) at different locations.Therefore, the method of superposition needs to be applied.This method allows for regarding multiple forces acting separately on the beam and superimposing the deflections resulting from each force [11].The deflections for each force are where  1,2 is the deflection at the needle tip for forces  1 and  2 combined.The up to now unknown forces  1 and  2 , with which a relation to   and   can be established, can be obtained from equilibrium conditions.These equilibrium conditions apply in the force and moment system of the free body diagrams in Figure 5.The cut in the free body diagram in Figure 5(a) was chosen because, in this case, the force  1 will cancel out at both sides, as  1 acts in the same direction at both sides of the cut.Force  1 at point B can also be regarded as support, which means that no shear forces are acting at the cut.Since in every point during the insertion equilibrium conditions are assumed, the moments acting at both the left and right sides of B must sum up to zero: As the bending moments at each side of the cut must be equal, ( 7) can be equalized to and  2 can be obtained.The force  1 can be obtained from the free body diagram in Figure 5(b).The sum of the forces acting perpendicularly along the needle must also be zero: Although here only deflection in the 2D plane is regarded, the model can be easily transferred to 3D, provided that a 4 DOF force/torque sensor is used.Also the deflections can be measured in both directions.If the needle deflects into positive -direction, all forces and moments, including forces/torques measured at the needle base, will simply point into the opposite direction.

Experimental Setup.
The setup for conducting insertion experiments into soft tissue consists of a robotic system with two degrees of freedom (DOF), which is capable of translational and rotational motion (see Figure 6).The translational motion is along the needle axis and rotational about the needle axis.The restriction to translational rotation along the needle axis is granted by a linear stage, which consists of a ball bearing mounted rail-carriage system.The system provides high precision while maintaining very low friction between rail and carriage (see Figure 6(a)).The carriage, which holds the rotational part, can be connected to a motor, which ensures linear insertion velocities (see Figure 6(b)).To carry out manual insertion experiments, the carriage needs to be entirely decoupled from the motor such that motor inertia does not affect the manual insertion procedure.A belt is used to connect motor and carriage as this provides the  most convenient way to couple or completely decouple motor and carriage by simply clamping the belt to the carriage when necessary.Mounted on the carriage is the motor for needle rotation, which holds a force/torque sensor, and the needle holder, which fixates the hand piece of the needle via two opposing set screws.Also attached to the carriage is a handle to perform manual insertions.
The resulting needle bevel location correction system is a marked improvement over current practice, which typically involves nonoverlapping insertion-only and rotation-only phases; simply making rotations at predetermined insertion depths is an example of the strategies used by implanters.The same is the case in robotic needle insertion research.For instance, in [12], the robot maintains a constant insertion velocity for the needle except when the insertion is paused to rotate the needle by a fixed amount.Or, in [13], highspeed rotation reversal at half the insertion depth was tested.In contrast to the above, we adjust the needle bevel location during the insertion and believe that such controlled, smooth, and incremental corrections can maximize needle tip targeting accuracy.Our research should not be confused with continuous, drilling-like rotation of a needle during insertion.Given that a bevel-tip needle's trajectory in soft tissue is curved, constantly spinning the needle at a fast rate helps to keep the needle straight but also increases tissue trauma.As opposed to this, we will only adjust the needle's rotational position on a calculated and intelligent basis.
For sensing forces, a JR3 6 DOF sensor of type 50M31A is used.The maximally admissible forces are 100, 100, and 200 N in , , and  direction, and maximum torques are 5 Nm for all three axes.The sensor's ADC has a resolution of 14 bits.The moments are measured about the center of the sensor.Since the relevant moments for deflection measurement occur at the tip of the needle holder (see Figure 6), the measured moments have to be recalculated in our setup.The moment (  ) acting at the tip of the needle holder can be expressed as where   and   are the moment and force measured by the sensor and  represents the distance from the center of the force sensor to the tip of the needle holder (52.75 mm).Both motors for rotational and translational motion are unipolar stepper motors.Stepper motors are used due to their ease of control and availability.The motors are controlled via L297 stepper controller ICs.The clock, enable, and direction signal are provided by a HILINK data acquisition (DAQ) board (see Figure 6(b)).The HILINK board interfaces a PC via RS232 and is controlled in MATLAB/Simulink.
The force sensor is interfaced via a C API.To read the data provided by the force sensor in Simulink, a C S-Function was written.Simulink S-Functions provide a way to run C/C++code in Simulink models.The S-Function block outputs two 3D vectors for forces and torques.By including the reading of sensor data in Simulink, the setup can entirely be controlled in Simulink.Real time control of the assembly is however only possible in soft real time, as the force sensor's S-Function cannot be used together with MATLAB's Real Time Windows kernel.This limits the maximally admissible sampling rate to roughly 20 Hz, which is a relatively low but sufficient rate for our purposes.
The used phantom tissue for the insertion experiments is liquid plastic, which is made of plastisol and produced by M-F Manufacturing Co.The stiffness of the plastic can be adjusted by the amount of added plastic softener.The thickness of the homogeneous tissue is roughly four centimeters, which provides enough weight to prevent too much shifting along the insertion axis during insertion.A key factor of the phantom tissue is its transparency since the needle needs to be tracked while being inserted in order to measure its deflection.To determine the deflection, images of the needle inside tissue are recorded from above by a Logitech C270 webcam during insertion.Since the needle deflection can only be monitored in the horizontal two-dimensional plane, the deflection also needs to be kept to this plane.This can be achieved by aligning the bevel vertically.This way it can be safely assumed that the needle will deflect in the horizontal plane only.Potentially occurring gravity effects, which could lead to deflection in the vertical plane, can be neglected as the magnitude is below the range of the sensor.Furthermore, the supporting effect of the tissue prevents the needle from being pulled down by gravity.
To track the needle tip during insertion for comparison to estimated deflection, template matching is used.The template is illustrated in Figure 7(a).A MATLAB implementation by Dirk-Jan Kroon is used.The used matching method is normalized cross-correlation (NCC).This approach provides a high robustness against brightness variations and specifically bubbles in the tissue (see Figure 7(b), bright spots on the needle).The images have a resolution of 800 × 448 pixels and the camera observes a width of 80 mm at the height of the needle.This results in a resolution of 0.179 mm/pixel.The deflection is measured in respect to the unbent needle in air.
The used needle is a standard 18 G brachytherapy needle with an outer diameter of 1.27 mm and inner diameter of

Results
In this section we present experimental results for the validation of the proposed deflection model, and for the illustration of the impact that needle turning has on the deflection.

Model Verification.
To verify the developed deflection model, insertion experiments were conducted consisting of six trials.Manual insertions with varying velocities and automated insertion with 10 mm/s and 15 mm/s were performed with two homogeneous tissue samples, which varied in stiffness.For this set of experiments, the needle was not turned at any point during the insertion.Forces and moments at the needle base and images were recorded during the insertion.The insertion depth was 120 mm throughout all the trials.Each trial consisted of 6 runs, in each of which a new point of insertion was used and the needle bevel was adjusted to deflect into the right direction as seen from the side of the robot.The recorded forces and moments were then used to estimate the deflection.The measured tip deflection was finally compared to the estimated tip deflection.For the manual insertion it was tried to keep the velocity at a constant level since the insertion speed was not measurable in real time.Before the start of insertion, the needle was inserted approximately 10 mm into the tissue.
Figure 8 shows an unfiltered sample deflection curve of the needle tip (see Figure 8(b)) and forces and moments at the needle base during insertion (see Figure 8(a)).The sensor data was filtered offline by a zero-phase lowpass filter, using MATLAB's filtfilt(⋅) function.Figure 8(b) also shows the deflections estimated by a model proposed by Abolhassani et al. [10] and the estimation of our proposed model.As the plot shows, the estimations of both models are very similar until later in the insertion process.At this point, our model maintains a relatively high precision whereas Abolhassani's model increasingly overestimates the deflection.
In Figures from 9(a) to 9(f), the mean estimations ( = 6) at insertion depths of 60 mm and 120 mm are illustrated for all 6 trials.The data shows that the aforementioned overestimation later on during the insertion can be observed throughout all the trials.It furthermore shows that our proposed model also slightly overestimates the deflection at a depth of 120 mm.At a depth of 60 mm, the estimations of both models are very close to the measured deflection throughout most of the trials.
To show where the estimations are in fact significantly different to the measured deflection, a paired t-test was carried out using MATLAB's function ttest(⋅).The null hypothesis of a paired t-test is that the mean of the difference of two data sets is not significantly different at a significance level of 5%.The test was carried out between measured and estimated samples ( = 6) for each model, for each   trial, and for the two observed insertion depths separately.
The results of the t-test can be found in Table 1.They show that, at the depth of 60 mm, all trials do not reject the null hypothesis, which means that both model estimations are not significantly different at this point.At the depth of 120 mm however, the test shows that both model estimations are significantly different from the measured deflection for most cases.

Correcting the Deflection.
A second type of experiment was performed to study the corrective effect of turning the needle about 180 ∘ when a certain estimated deflection threshold is reached.The needle is turned with an angular velocity of 180 deg/s.As the needle is being turned, the insertion is not paused but continues at the same speed.
To estimate the threshold, the virtual sensor, which utilizes the developed model, was used.Two trials were executed, also consisting of six runs each.In each trial the threshold was set to a different level, the first being 1 mm and the second 5 mm.The needle also deflected in each run to the right side, or negative -direction, until it was turned.One of the tissue samples, which was used for the first experiment (tissue 1), was used for this set of experiments.The runs involving turning the needle were conducted with a velocity of 5 mm/s, and, for the comparison with deflection without turning, data from the model verification

Figure 1 :Figure 2 :
Figure 1: View of a prostate brachytherapy procedure as it is currently performed.

Figure 3 :
Figure 3: The distributed load  acting along the inserted portion of the needle as a reaction force to the needle, as proposed by Kataoka et al. [3].  is the resulting force at the needle base.The 2D diagram plane is the deflection plane of the needle.The dashed arrow represents the resulting force acting at the geometric center (centroid) of the area of the distributed load .

Figure 4 :
Figure 4: (a) The triangularly distributed loads  1 and  2 acting along the needle while inserting into tissue.  is the resulting force at the needle base. 1 and  2 -shown with dashed lines-are the resulting forces from the distributed loads  1 and  2 .(b) The location of the centroid of a triangle.

Figure 5 :
Figure 5: Free body diagrams of the needle while inserting into tissue.

Figure 6 :
Figure 6: The setup of the intelligent instrument prototype.
(a) The template used for tracking the needle tip (b) Tracked needle tip, indicated by the red cross at the tip

Figure 7 :
Figure 7: Tracking of the needle tip.

ForceFigure 8 :
Figure 8: Sample plots for forces and moments at the needle base and tip deflection during insertion into tissue 1 with a velocity of 10 mm/s.
Figures 9(a) to 9(f) and the t-test confirm the previously made assertion based on Figure 8(b).
([11, p. 1084]) where  1 and  2 are the deflections at the needle tip for  1 and  2 , respectively,  in and  out are the length proportions inside and outside of tissue,  represents the total needle length, and  and  are factors to adjust the points at which forces  1 and  2 act.Superimposing the two deflections results in

Table 1 :
Results of a paired -test performed on estimated and measured data.
a Insertion speed in mm/s or "man" for manual.b Model identifier, "Abl" for Abolhassani.c Insertion depth in mm.d "" if null hypothesis is not rejected.