Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3004482

Formats

Article sections

- Abstract
- I. INTRODUCTION
- II. RELATED WORK
- III. PROBLEM STATEMENT AND MODELING ASSUMPTIONS
- IV. FEEDBACK CONTROLLER WITH HELICAL PATHS
- V. FAST TRAJECTORY SELECTION
- VI. SIMULATIONS IN RIGID TISSUE
- VII. SIMULATIONS IN DEFORMABLE TISSUE
- VIII. CONCLUSION
- REFERENCES

Authors

Related links

Robot Sci Syst. Author manuscript; available in PMC 2010 December 20.

Published in final edited form as:

Robot Sci Syst. 2009 June 28; V: 37.

PMCID: PMC3004482

NIHMSID: NIHMS191732

Kris Hauser,^{*} Ron Alterovitz,^{†} Nuttapong Chentanez,^{*} Allison Okamura,^{‡} and Ken Goldberg^{*}

Kris Hauser: ude.yelekreb@resuah; Ron Alterovitz: ude.cnu.sc@nor; Nuttapong Chentanez: ude.yelekreb@natnehcn; Allison Okamura: ude.uhj@arumakoa; Ken Goldberg: ude.yelekreb@grebdlog

See other articles in PMC that cite the published article.

Bevel-tip steerable needles are a promising new technology for improving accuracy and accessibility in minimally invasive medical procedures. As yet, 3D needle steering has not been demonstrated in the presence of tissue deformation and uncertainty, despite the application of progressively more sophisticated planning algorithms. This paper presents a feedback controller that steers a needle along 3D helical paths, and varies the helix radius to correct for perturbations. It achieves high accuracy for targets sufficiently far from the needle insertion point; this is counterintuitive because the system is highly under-actuated and not locally controllable. The controller uses a model predictive control framework that chooses a needle twist rate such that the predicted helical trajectory minimizes the distance to the target. Fast branch and bound techniques enable execution at kilohertz rates on a 2GHz PC. We evaluate the controller under a variety of simulated perturbations, including imaging noise, needle deflections, and curvature estimation errors. We also test the controller in a 3D finite element simulator that incorporates deformation in the tissue as well as the needle. In deformable tissue examples, the controller reduced targeting error by up to 88% compared to open-loop execution.

Needles are used in medicine for a wide range of diagnostic and therapy delivery procedures. Many applications, such as biopsies and prostate brachytherapy, require the needle tip to be positioned accurately at a target in the tissue. Such procedures require a great deal of clinician skill and/or trial and error, and can be difficult even under image guidance. Errors are introduced when positioning the needle for insertion, and the needle can be deflected from its intended path by tissue inhomogeneities. Furthermore, the target may shift due to tissue deformations caused by the needle or patient motion.

Asymmetric-tip flexible needles are a new class of needles, developed in collaboration between Johns Hopkins and U.C. Berkeley, that can steer through soft tissue [24]. As it is inserted, the needle travels along a curved path due to the asymmetric cutting forces generated by either a bevel cutting tip or a pre-bent kink. The arc direction can be controlled by twisting the base of the needle, which allows the needle to travel in circular arcs [24] (by holding twist constant) and helical trajectories (by simultaneously twisting and inserting the base at constant velocity) [12]. In principle, these inputs can be used to correct for needle placement errors and deflections during insertion. They may also steer around obstacles such as bones or sensitive organs to reach targets that are inaccessible to standard rigid needles.

A variety of techniques have been used for 3D path planning of steerable needle trajectories [12, 13, 20, 26]. But these techniques assume an idealized rigid tissue model; in practice, the needle will deviate from the planned path due to deformation effects (of both tissue and needle), tissue inhomogeneities, and patient variability. Planners that address uncertainty and deformations are currently limited to 2D needle steering, and perform extensive precomputations [1, 2].

Our work breaks with the trend toward increasingly precise modeling and more sophisticated (and computationally expensive) planning techniques. We present a closed-loop, model predictive controller that steers the needle along helical paths to reach a desired target position. Perturbations are sensed using imaging feedback, and the controller reacts by adjusting the radius and heading of the helix to pass as close as possible to the target. The optimization step searches only the trajectory space of constant controls, and does not predict future changes in the control input. A fast branch and bound search enables trajectory corrections to be computed at kilohertz rates.

We evaluate the controller on a wide variety of simulation experiments in both rigid tissue and in a finite element deformable tissue simulator (Figure 1). In light of the fact that steerable needles are underactuated and non-controllable, our results are surprising. Assuming no perturbations, the controller achieves high accuracy (less than 1% of the needle’s radius of curvature *r*) for all targets sufficiently far from the insertion point (approximately twice *r*). It attains reasonable accuracy (less than 5% of *r*) under a variety of simulated Markovian and non-Markovian perturbations, including imaging noise, needle deflection, and curvature estimation errors. In deformable tissue, the controller can compensate for deformation effects to achieve much higher accuracy than open-loop execution. In the example of Figure 1, the controller reduces targeting error by 88%.

Needle steering in a soft block of tissue. Left: the needle pierces the surface. Dotted curve is the predicted helical path of the needle. Right: feedback control corrects for errors due to deformation. More examples can be found at http://automation.berkeley.edu/projects/needlesteering/ **...**

Several mechanisms for needle steering have been proposed. Symmetric-tip flexible needles can be steered by translation and rotation at the needle base [10, 14]. Bevel-tip flexible needles achieve steering with a tip that produces asymmetric cutting forces, causing them to bend when inserted into soft tissue [24, 25]. Recent experiments with pre-bent needle tips have achieved up to a 6.1 cm radius of curvature [22]. New needle designs are expected to further decrease this radius.

Bevel-tip and pre-bent needles are controlled by two degrees of freedom: insertion and twisting at the needle base. When restricted to a plane, the needle moves like a Dubins car that can only turn left or right. In 3D, the needle moves in approximately helical trajectories. The radius of curvature of the path can be adjusted during insertion using duty cycling, enabling the needle to be driven straight forward [18].

Planning motions of 1D flexible objects has been studied in the context of redundant manipulators [6, 7, 8], ropes [23] and symmetric-tip flexible needles [10, 14]. Several researchers have developed motion planners for steerable needles in 3D tissues [12, 13, 20, 26]. These methods have all assumed that tissue is rigid. Duindam et al. derived a closed-form inverse kinematics solution for reaching a desired position orientation in 3D [13]. The solution uses up to three stop-and-turn motions in which the needle, between turns, follows a circular arc. Inverse kinematics solutions for helical paths have not been found. Planning with helices in 3D has been achieved using numerical optimization [12], an error diffusion technique [20], and a sampling-based planner [26].

Needle and tissue deformation will deflect the needle from an idealized trajectory. In such cases a controller may be used to keep a needle on its intended path. Kallem and Cowan developed controllers to stabilize a needle to a given plane [17]. For 2D, Alterovitz et al. developed an optimization-based planner that predicts tissue deformations using finite element simulation [2], although this assumes no uncertainty in the motion model or tissue properties. Alterovitz et al. developed a 2D motion planner for needle steering with obstacles and Markov motion uncertainty [1], but did not model large-scale tissue deformations or curvature estimation errors.

We wish to place the needle tip at a certain goal location **p**_{g} in the tissue. The desired needle orientation is not specified. We assume we have access to an imaging system that periodically provides an estimate of the world-space coordinate of **p**_{g} as it moves due to tissue deformations, as well as the position, heading, and bevel direction of the needle tip. We also assume a given workspace ^{3} in which the needle is permitted to travel. In this work, is taken to be either all of ^{3} or a bounding box.

We assume an insertion velocity is chosen for the duration of the procedure (on the order of 1 cm/s). The controller outputs a twist rate to be applied at the needle base. We impose a maximum twist rate _{max} (specified in rad/m) because fast spinning may damage the tissue. A robot or physician will simultaneously insert and twist for Δ*t* seconds at the twist rate computed by the controller. The controller then repeats.

Though the controller is designed to operate in deformable tissue, it uses a purely kinematic rigid tissue motion model [24] for computational efficiency. This model assumes the dynamics of the needle are fully determined at the needle tip, such that the tip travels with constant curvature 1/*r* in the bevel direction. Thus, if inserted without twisting, the needle trajectory is a circle with radius of curvature *r*. Experiments suggest that this assumption approximately holds in homogeneous materials [24]. The model also assumes that a twist at the base of the needle is transmitted directly to the needle tip, and does not deflect the tip position. This assumption neglects torsional friction along the needle shaft, which causes a “lag” before twists are felt at the tip. Control techniques have been developed to counteract these effects [21].

We place a frame **x, y, z** at the needle tip such that **z** is forward and **y** is the bevel direction (Figure 2). Inserting the needle causes the tip to move along **z** with velocity υ(*t*) and rotate about x with angular velocity υ (*t*)/*r*. A twist rotates the tip about **z** with angular velocity (*t*). In screw notation, the instantaneous velocity of the frame is given by:

$$\widehat{V}(t)=\left[\begin{array}{cccc}\hfill 0\hfill & \hfill -\phi (t)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \phi (t)\hfill & \hfill 0\hfill & \hfill -\upsilon (t)/r\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \upsilon (t)/r\hfill & \hfill 0\hfill & \hfill \upsilon (t)\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right]$$

If υ(*t*) ≠ 0, then we can factor out υ(*t*) so that the needle travels at unit velocity, and the only variable to control is the twist rate (*t*). If the twist rate is held constant, then after inserting the needle a distance *d*, the tip frame is transformed by the quantity *T*(*d*) = exp(*d*) [12].

We model a trajectory as a sequence of screw motions with piecewise constant twist rates (*d*) parameterized by insertion distance. Denote **t**_{(d)} (*d*) as the trajectory of the tip position under twists (*d*). We drop the subscript if the twists are unimportant.

We assume a model predictive control (MPC) framework that iterates the following steps (Figure 3):

*Propose*. Generate a set of proposal trajectories that start from the current estimated needle state.*Select*. Find the trajectory in with control (*d*) that achieves the minimal distance to the target.*Execute*. The needle is then inserted with twists according to (*d*) and constant velocity for time Δ*t*.

The process repeats until no trajectory can improve the distance to the target by more than ε (we set ε=0.2%*r*).

Overview of the control policy. Curves emanating from the needle depict proposal trajectories. Circle depicts the target.

The controller meets basic stability criteria: first, it will terminate, because the predicted distance to the target must decrease by at least ε every iteration; second, the final error will be no larger than the distance from the target to . Convergence is difficult to prove, and depends on the selection of . Empirical results in Section VI suggest that setting to the set of constant-twist-rate helices is sufficient to reach most target locations with high accuracy, even under perturbations.

This is surprising and counterintuitive for two reasons. First, the system is not locally controllable, even if the needle were allowed to move backward [17]. Even though our controller cannot reduce errors asymptotically to zero, it often has low error upon termination. Second, traditional MPC techniques use expensive numerical optimizations at every time step to compute (e.g. receding horizon optimal control). Our controller performs well even without considering future changes in twist rate, so each time step is computed efficiently.

In the remainder of this section, we derive a closed form expression for the helical trajectories followed by the needle tip when inserted with constant velocity and twist rate . We let the set of proposal trajectories contain all such helices for [−_{max}, _{max}]. Using a visualization of , we argue that the controller converges for sufficiently distant targets.

Let **t**_{}(*d*) be the needle tip trajectory followed under a constant twist rate . We wish to find the helix that describes this trajectory relative to the initial frame of the tip. A helix with radius *a*, slope θ, and oriented in the **z** axis has the arc-length parameterization

$$\mathbf{h}(d)=\left[\begin{array}{c}\hfill a\phantom{\rule{thinmathspace}{0ex}}\text{cos}((d/a)\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\theta )\hfill \\ \hfill a\phantom{\rule{thinmathspace}{0ex}}\text{sin}((d/a)\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\theta )\hfill \\ \hfill d\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta \hfill \end{array}\right]\phantom{\rule{thinmathspace}{0ex}}.$$

(1)

Here we derive the helix parameters *a*, θ and a rigid transformation *A* such that *A***h**(*d*) = **t**_{}(*d*). Denote **h _{x}**,

$$A\mathbf{h}(0)={\mathbf{t}}_{\phi}(0)\Rightarrow a{\mathbf{h}}_{\mathbf{x}}+{\mathbf{h}}_{\mathbf{o}}={[0,0,0]}^{T}$$

(2)

$$A\mathbf{h}\prime (0)={\mathbf{t}}_{\phi}^{\prime}(0)\Rightarrow \text{cos}\phantom{\rule{thinmathspace}{0ex}}\theta {\mathbf{h}}_{\mathbf{y}}+\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta {\mathbf{h}}_{\mathbf{z}}={[0,0,1]}^{T}$$

(3)

$$A\mathbf{h}\u2033(0)={\mathbf{t}}_{\phi}^{\u2033}(0)\Rightarrow -{\text{cos}}^{2}\theta /a{\mathbf{h}}_{\mathbf{x}}={[0,-1/r,0]}^{T}$$

(4)

Let ω = [1/*r*, 0, ]^{T} be the angular velocity of the needle tip. The axis of rotation must be aligned with the helix axis, so **h _{z}** = ω/‖ω‖. Let $q=\Vert \omega \Vert =\sqrt{1/{r}^{2}+{\phi}^{2}}$. Taking the dot product of both sides of (3) with

$$\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta =\phi /q\text{andcos}\phantom{\rule{thinmathspace}{0ex}}\theta =1/(\mathit{\text{rq}}).$$

From (4) we see

$$a=r\phantom{\rule{thinmathspace}{0ex}}{\text{cos}}^{2}\theta =1/({\mathit{\text{rq}}}^{2})\text{and}{\mathbf{h}}_{\mathbf{x}}={[0,1,0]}^{T}.$$

Setting **h _{y}** orthogonal to

$${\mathbf{h}}_{\mathbf{y}}={[-\phi ,0,1/r]}^{T}/q.$$

Finally, from (2) we have **h _{o}** = [0, −

In summary, letting tan θ = *r*, we have

$${\mathbf{t}}_{\phi}(d)=\left[\begin{array}{c}\hfill -a\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta \phantom{\rule{thinmathspace}{0ex}}\text{sin}((d/r)\phantom{\rule{thinmathspace}{0ex}}\text{csc}\phantom{\rule{thinmathspace}{0ex}}\theta )+d\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta \phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \hfill \\ \hfill a\phantom{\rule{thinmathspace}{0ex}}\text{cos}((d/r)\phantom{\rule{thinmathspace}{0ex}}\text{csc}\phantom{\rule{thinmathspace}{0ex}}\theta )-a\hfill \\ \hfill a\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\theta \phantom{\rule{thinmathspace}{0ex}}\text{sin}((d/r)\phantom{\rule{thinmathspace}{0ex}}\text{csc}\phantom{\rule{thinmathspace}{0ex}}\theta )+d\phantom{\rule{thinmathspace}{0ex}}{\text{sin}}^{2}\phantom{\rule{thinmathspace}{0ex}}\theta \hfill \end{array}\right]$$

(5)

Where *a* = *r* cos^{2} θ = *r*/(1/*r*^{2} + ^{2}) is the helix radius, and [cos θ, 0, sin θ]^{T} is the helix axis. The tip moves at the speed
$\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta =\phi /\sqrt{1/{r}^{2}+{\phi}^{2}}$ with respect to the helix axis.

Let (_{max}) = {**t**_{}(*d*)| [−_{max}, _{max}], *d* ≥ 0} be the set of all needle tip positions reachable under constant twist rate, given twist rate bounds [−_{max}, _{max}]. Figure 4 plots a portion of (∞).

Four views of the constant-twist-rate reachable set: *z*, −*y* plane, *z*, *x* plane, perspective, and surface shaded. Twist rate in (−∞, ∞), insertion distance in [0, 2π*r*].

(_{max}) is scale-invariant w.r.t. radius of curvature, in the following sense. Let _{α}(*d*) denote the tip trajectory with curvature α/*r* and constant twist rate α. It is easily verified that _{α}(α*d*) = α**t**_{}(*d*). Thus, (α_{max}) = α(_{max}).

This visualization helps understand why the controller should converge. The projection of (∞) on the *x, z* plane covers the entire *z* > 0 halfplane, and all twists produce axes of rotation that are parallel to the *x, z* plane. Thus, if (∞) is rotated by a full turn about any axis, it will sweep out a large portion of ^{3}. Most trajectories make a full turn quickly — in fact, constructing a trajectory that does not ever make a full turn is extremely difficult. For any target point in this swept volume, errors will converge to zero.

For a finite maximum twist rate _{max}, (_{max}) has a gap of coverage along the **z** axis, shaped approximately like a wedge of angle | arctan *r*_{max}|. Consider filling this gap with a maneuver that makes a full turn of the helix with twist rate _{max}, and another with twist rate −_{max} (Figure 5).

Intuitively, the maneuver allows the needle to travel approximately straight forward. After each full turn, the needle tip is displaced by (2π*r* sin θ cos^{2} θ, 0, 2π*r* sin^{2} θ cos θ) where θ = arctan *r*_{max}. Orientation is unchanged. Executing both turns places the needle tip at (0, 0, 4π*r* sin^{2} θ cos θ). The path length of the maneuver is *d* = 4π*r* cos θ.

Because it is composed of helices, the maneuver can be included in the optimization of Section V with only minor changes. Experiments in Section VI-A demonstrate that accuracy is improved when is augmented with this maneuver.

The speed of computing the proposal and selection steps limit the speed at which the needle can be inserted. If more than Δ*t* time is spent, the needle must be halted while the computation finishes. Thus, it is critical for the closest point optimization in the selection step to be fast. In this section we present an efficient branch-and-bound technique for computing the closest point. We also present a bounding volume method for incorporating workspace constraints.

As can be seen in Figure 4, the reachable set (_{max}) has many valleys and ridges, so local optimization techniques easily get stuck in local minima. Instead, we use branch-and-bound optimization to ensure that a global optimum is found.

We seek to minimize the function *f*(, *d*) = ‖**t**_{}(*d*) − **p**‖ over the space = {(, *d*) | || ≤ _{max}, *d* > 0}. Here **p** denotes the coordinates of **p**_{g} relative to the needle tip frame. We use an auxiliary function *f*_{L}(*R*), described in Section V-C, that computes a lower bound on the function value attained in a region *R* . A search tree is then built by recursively splitting the space into subregions which may contain the optimum. When a subregion is split, we also test the value of *f* at its midpoint. We maintain the helix ^{} and insertion distance *d*^{} that give the minimum value of *f̠* seen thus far. If *f*_{L}(*R*) is found to be larger than *f̠*, then subregion *R* does not contain the optimum and can be safely pruned. The process continues until *f̠* achieves an ε tolerance.

In subsequent iterations of the control loop, the closest point does not change very much, so we can exploit this temporal coherence to help prune many intervals quickly. We keep the optimal twist ^{} and insertion distance *d*^{} from the previous iteration, and initialize *f̠* = *f*(^{}, *d*^{}). In our experiments this reduces overall running time by about a factor of four.

The efficiency of branch and bound depends greatly on the tightness of the lower bound. Interval analysis could be used to automatically compute the bounds [16], but the bounds can be extremely poor, depending on how (5) is formulated. In our experiments, the running time of a standard interval analysis implementation ranged from milliseconds to several seconds, and was sensitive to the location query point.

Our implementation is as follows. We partially discretize by selecting *n* = 100 possible values for . For a more uniform coverage of the reachable set given a fixed *n*, we discretize evenly in the θ = arctan *r* coordinate. To include other proposal trajectories in the search, we divide them into helix segments, and add the segments to . Then, for each helix, we narrow further by computing a narrow initial interval of insertion distances that is guaranteed to contain the closest point. We then employ a good lower bound on the distance from a helix segment to a target point.

We derive two intervals that must contain the closest point using geometrical reasoning; their intersection provides a tight initial interval for the search. Reparameterize (1) such that *u* = (*d/a*) cos θ, and let *b* = *a* tan θ. Write the helix and point **p** in cylindrical coordinates with respect to the helix origin. Denoting cylindrical coordinates with the superscript ·^{c}, we have **h**^{c}(*u*) = [*a, u, ub*]^{T}. Let **p**^{c} = [*r _{p}*, θ

First, the helix attains *z*-coordinate *z _{p}* at

Two initial intervals for finding the closest point from **p** to a helix. (a) The interval containing half a turn of the helix, centered about the *z* = *z*_{p} plane. (b) The interval containing the *z*-range of the intersection between a sphere of radius *d*_{R} centered **...**

Second, we examine the cylinder bounding the helix. If *d _{R}* is the distance from

Now we describe a lower bound *f _{L}* for

$${\mathbf{q}}^{c}=\left[\begin{array}{c}\hfill a\hfill \\ \hfill C({\theta}_{p},[\underset{\xaf}{u},\overline{u}])\hfill \\ \hfill C({z}_{p},[\underset{\xaf}{u}b,\overline{u}b])\hfill \end{array}\right]$$

(6)

The Cartesian distance ‖**q** − **p**‖ is a lower bound.

We consider workspace constraints using a simple extension to the selection step. Before computing the helix-point distance, we truncate the trajectory **t**(*d*) if it collides with the workspace boundary. In other words, we find the maximum *D* for which **t**(*d*) remains in . Then, the initial branch-and-bound search interval for the closest point (Section V-B) is restricted to lie in [0, *D*].

We use an recursive bisection to find *D* up to tolerance ε. We start the recursion with interval *I* = [0, ], where is the upper bound on the insertion distance from Section V-B. We then recurse as follows. If the size of *I* is less than ε we set *D* to its lower bound, as we cannot rule out the possibility of collision in *I*. Otherwise, we compute the axis-aligned bounding box containing the portion of **t**(*d*) for all *d* *I* using interval computations. If the box lies within , the trajectory segment lies entirely in the workspace. If not, we bisect *I*, and recurse on the lower half. If no collision is reported in the lower half, then we recurse on the upper half.

Since no interval will have size less than , the number of intervals examined by the branch and bound algorithm is no more than 2*L*/ε, where *L* is the total length of all initial intervals. Each of the *n* initial intervals has length at most 2π*r*, so the number of intervals is bounded by 4π*rn*/ε. In practice, this bound is extremely loose; the number of examined intervals rarely exceeds a few hundred.

We performed timing experiments on a 2 GHz laptop PC. With an unbounded workspace, the average controller iteration was computed in 0.11 ms, with standard deviation of 1.3 ms. On rare occasions, particularly as the target nears the needle tip, it is much slower, at one point taking 63 ms. More computation is required for bounded workspaces with obstacles. In a box domain with 10 spherical obstacles, each iteration averaged 0.41 ms, with standard deviation 2.52 ms.

This section performs a set of simulation experiments to evaluate the accuracy of the controller, as a function of the target position. We simulate the needle using the rigid tissue motion model, first under varying simulation parameters, and then under a variety of perturbations. We define accuracy as the final distance from the needle tip to the target when the controller terminates.

For comparison, we instantiate a *reference controller* that uses the following parameters: a refresh occurs every 2%*r* of insertion distance, maximum twist rate _{max} = 10π rad/*r*, and the workspace is all of ^{3}.

Our first experiment evaluates how accuracy of the reference controller varies over the space of targets, under the rigid-tissue motion model. Figure 7 plots accuracy as the target position varies in 2D slices of space. Grids are defined parallel to the *z*, −*y* plane (see Figure 7). We then exhaustively enumerate each grid cell, set the target **p**_{g} to the cell center, and then run the controller from the initial configuration to completion. The final error from the needle tip to the target is recorded as the accuracy for that cell. The plots show convergence to sub-millimeter distances almost everywhere, except for a narrow strip along the +*y* axis, and a circular spot directly above the needle tip in the −*y* direction. As the *x* value deviates further from 0, the radius of this spot shrinks. Comparing to Figure 4, we see this region corresponds to the large “lobe” of the reachable set.

Controller accuracy in rigid tissue as a function of target’s initial relative location. For a given target, the reported accuracy is the final error *e* achieved by the controller when the needle is inserted at (0, 0, 0). Targets are varied over **...**

The next experiments evaluate accuracy when the controller still has an exact motion model, but with the following control parameters varied:

- Maximum twist rate reduced to 2.5πrad/
*r*. - Without alternating-twist maneuvers (see Section IV-D).
- In a bounded workspace [−3
*r*, 3*r*] × [−3*r*, 3*r*] × [0*r*, 6*r*].

For each setting of parameters, we evaluate the accuracy distribution for targets in a region *R* in space. We chose *R* to be the rectangle [0, 0] × [−3*r*, 3*r*] × [1.5*r*, 6*r*], which is the rightmost 3/4 of the *x* = 0*r* plot of Figure 7. This region was chosen to exclude the circular spot where the reference controller does not converge. For each point in a 50 × 50 grid over *R*, we set the target to the point, and ran the controller from start to finish. Table I reports the average accuracy over all grid points and its standard deviation.

The second set of experiments introduce a variety of perturbations and modeling errors into the simulation. The perturbations are discussed in detail below. Table II summarizes the results. Accuracy is measured exactly as in Table I.

All perturbations degrade accuracy to some extent, but rarely drive the error above 10%*r*. This suggests that the controller is stable. Noise in the target position and needle pose estimation can cause significant reductions in accuracy, but in practice a smoothing filter will help mitigate these effects. Twist lag causes high variance in accuracy, in some cases exceeding 20%*r* error.

*Noise introduced into the target position estimation:*Tracking the target position may require the use of deformable image registration techniques (e.g. [3]), so target position estimates will contain errors due to imaging noise and deformation modeling errors. To simulate, we randomly perturb the target position estimate at each time step with isotropic Gaussian noise with standard deviation σ. We performed trials with σ = 2%*r*and with σ = 4%*r*. The random number generator is seeded with the result of the C`time`procedure at the beginning of each trial.*Noise introduced into the needle pose estimation:*Limited resolution of medical imaging (approximately 0.8 mm for modern ultrasound [11]) will cause errors in estimating the tip position and orientation. The position and direction of the tip may be estimated directly from the image, but more sophisticated state estimation models are needed to estimate the twist about the needle shaft because the bevel tip is too small to be seen [17]. To simulate these errors, we randomly perturb the needle position by isotropic Gaussian noise with standard deviation σ_{p}= 2%*r*and orientation by a rotation about a random axis and a random rotation angle with standard deviation σ_{o}= 2°.*Random deflections of the needle tip:*We simulate random deflections of the needle tip that accumulate over time, and can be mathematically described as a Wiener process. The standard deviation of the position and orientation deflections are respectively σ_{p}= 0.44*r*and σ_{o}= 11.1° per each*r*distance the needle is inserted.*Curvature Estimation Errors:*We set the actual radius of curvature to = 0.8*r*and = 1.2*r*, while the estimated radius remains at*r*.*Twist Lag:*As a needle is twisted inside tissue, needle-tissue friction causes internal torsion to build up along the length of the needle. This causes twist lag, a phenomenon where twists applied at the base are not directly commanded to the tip [21]. We model twist lag by nullifying all base twists if |θ_{b}− θ_{t}| < η*L*, where θ_{b}and θ_{t}are the base and tip angle,*L*is the length of needle inside the tissue, and η is a constant (with units of angle over distance). In other words, we evolve θ_{b}and θ_{t}with the equations_{b}= , and_{t}= 0 if |θ_{b}− θ_{t}| < η*L*, or otherwise. We set η = 6°/*r*, which is consistent with the experimental data in [21] assuming*r*=5 cm.

We evaluate the convergence of our controller in a finite element (FEM) simulation of a flexible needle in deformable tissue. The simulator is presented in detail in [5]. Past work has addressed simulation of rigid needles in 2D and 3D deformable FEM models [9, 15, 19], and symmetric tip flexible needles using a linear beam model [14]. Our simulator couples a discrete elastic rod flexible needle model with a FEM tissue simulation. The needle is modeled as a piecewise linear curve connecting a sequence of nodes. These nodes are coupled to nodes of the tissue, which is modeled as a tetrahedral mesh. Viscous and sticking friction along the needle shaft is simulated, but torsional friction is not.

We performed experiments with an *r*=5 cm needle in a 10 cm cube of tissue constrained on the bottom face. The tissue mesh is discretized at 1.25 cm intervals, and contains 729 nodes and 3072 tetrahedra. Tissue constitutive properties were set to approximately those of prostate tissue (e.g., elastic modulus 60 kPa) [27]. The friction coefficient between the tissue and needle was set to 0.1 and the minimum cutting force was set to 1 N. We tuned needle stiffness to achieve a turning radius of approximately 5 cm (±0.3 cm) on a small number of experiments, where we varied the insertion point and angle. The tissue cube is used as the workspace .

A single trial is shown in Figure 8. The needle tip, initially at (0 cm,5 cm,5 cm) and pointing horizontally, reaches the target (8 cm,2 cm,8 cm) with 0.88 mm accuracy. Open-loop execution of a trajectory planned in rigid tissue achieved 5.9 mm error, so the controller’s improvement is 85%. Figure 9.a plots accuracy for targets in the *x, z* slice at *y*=5 cm, showing good agreement with accuracy in rigid tissue (Figure 9.b).

A simulated trial of the controller in deformable tissue: initial setup, cutting through the surface, tracking the target, final needle configuration. Bottom plane is fixed. Grid used for visualization is at half the resolution of the simulation mesh. **...**

Accuracy images for a deformable 10 cm cube compared to a rigid cube. Needle was inserted at (0 cm, 5 cm). Slice taken down the middle of the cube. Because deformable tissue simulation is computationally expensive, a lower resolution grid was used. Units **...**

We also tested our controller under exaggerated deformation. In Figure 1 we increased the cutting force, friction coefficient, and viscous damping by factor of 10. The tip entered the tissue at a downward angle of 30°. The final error achieved by open loop execution is 13.5 mm, while the controller achieved 1.59 mm, an 88% improvement.

Clinical imaging may impart excessive radiation when used continuously, so we tested controller performance when the imaging rate is reduced. Experiments show less than 2 mm error in both of the above examples even when up to 30 mm of needle travel occurs between imaging refreshes.

We presented a closed-loop control policy for steering bevel-tip flexible needles in deformable tissue. It rapidly iterates a single rule: execute the helical trajectory that minimizes the distance to the target. Fast branch and bound techniques that enable the controller to run in real-time on standard PC hardware. We evaluated the accuracy of the controller in rigid-tissue simulations under simulated disturbances as well as simulations of 3D finite element deformable tissue. The controller reaches relatively distant targets with high accuracy, which is surprising because steerable needles have nonholonomic and non-controllable dynamics.

A primary goal for future work is integration and evaluation on robotic hardware and tissue phantoms. We also plan to address obstacles in the workspace. Preliminary work suggests that the presented controller can already avoid simple obstacles by excluding them from the workspace. Complex workspaces might be addressed by decomposition into convex regions, or more sophisticated planners to generate obstacle-avoiding proposal trajectories. We will also seek to understand the effects of different proposal trajectories, which may enable better workspace coverage, or enable reaching targets with orientation specified. Finally, traditional asymptotic stability criteria (e.g., [4]) do not apply to our system, so we are pursuing alternate methods to further understand why and when the controller converges.

This work is supported in part by NIH grant R01 EB 006435. We thank the reviewers for their helpful comments. We thank our collaborators and colleagues N. Cowan and G. Chirikjian of Johns Hopkins, J. O’Brien, J. Shewchuk, V. Duindam, and R. Jansen of U.C. Berkeley, and J. Pouliot, I-C. Hsu, and A. Cunha of UCSF for their valuable insights and suggestions.

1. Alterovitz R, Simeon T, Goldberg K. The stochastic motion roadmap: A sampling framework for planning with markov motion uncertainty. Robotics: Science and Systems. 2007 June

2. Alterovitz R, Goldberg K, Okamura AM. Planning for steerable bevel-tip needle insertion through 2D soft tissue with obstacles; IEEE Intl. Conf. on Robot. and Autom; 2005. Apr., pp. 1652–1657.

3. Bowden A, Rabbitt R, Weiss J. Anatomical registration and segmentation by warping template finite element models. Proc. SPIE (Laser-Tissue Interaction IX) 1998;vol. 3254:469–476.

4. Chen H, Allgöwer F. Nonlinear model predictive control schemes with guaranteed stability. In: Berber R, Kravaris C, editors. Nonlinear Model Based Process Control. Dodrecht: Kluwer Academic Publishers; 1998. pp. 465–494.

5. Chentanez N, Alterovitz R, Ritchie D, Cho L, Hauser K, Goldberg K, Shewchuk J, O’Brien J. Interactive simulation of surgical needle insertion and steering. ACM SIGGRAPH. 2009 Aug

6. Chirikjian G, Burdick J. A geometric approach to hyper-redundant manipulator obstacle avoidance. ASME Journal of Mechanical Design. 1992 December;vol. 114:580–585.

7. Chirikjian G, Burdick J. A modal approach to hyper-redundant manipulator kinematics. IEEE Trans. Robot. and Autom. 1994;vol. 10:343–354.

8. Chirikjian G, Ebert-Uphoff I. Numerical convolution on the euclidean group with applications to workspace generation. IEEE Trans. Robot. and Autom. 1998 February;vol. 14(no. 1):123–136.

9. DiMaio S, Salcudean S. Needle insertion modelling and simulation; IEEE Intl. Conf. on Robot. and Autom; 2002. pp. 2098–2105.

10. DiMaio S, Salcudean S. Needle steering and motion planning in soft tissues. IEEE Trans. Biomedical Engineering. 2005;vol. 52(no. 6):965–974. [PubMed]

11. Ding M, Cardinal H. Automatic needle segmentation in three-dimensional ultrasound image using two orthogonal two-dimensional image projections. Medical Physics. 2003;vol. 30(no. 2):222–234. [PubMed]

12. Duindam V, Alterovitz R, Sastry S, Goldberg K. Screw-based motion planning for bevel-tip flexible needles in 3d environments with obstacles; IEEE Intl. Conf. on Robot. and Autom; 2008. May, pp. 2483–2488. [PMC free article] [PubMed]

13. Duindam V, Xu J, Alterovitz R, Sastry S, Goldberg K. 3d motion planning algorithms for steerable needles using inverse kinematics. Workshop on the Algorithmic Foundations of Robotics. 2008 [PMC free article] [PubMed]

14. Glozman D, Shoham M. Flexible needle steering for percutaneous therapies. Computer Aided Surgery. 2006;vol. 11(no. 4):194–201. [PubMed]

15. Goksel O, Salcudean SE, DiMaio SP, Rohling R, Morris J. 3d needle-tissue interaction simulation for prostate brachytherapy. Medical Image Computing and Computer-Assisted Intervention MICCAI. 2005;2005:827–834. [PubMed]

16. Hansen ER. Global Optimization Using Interval Analysis. New York: Marcel Dekker; 1992.

17. Kallem V, Cowan NJ. Image-guided control of flexible bevel-tip needles; IEEE Intl. Conf. on Robot. and Autom; 2007. [PMC free article] [PubMed]

18. Minhas D, Engh JA, Fenske MM, Riviere C. Modeling of needle steering via duty-cycled spinning; Proc. 29th Annual Intl. Conf. of the IEEE Eng. in Medicine and Biology Society; 2007. Aug, pp. 2756–2759. [PubMed]

19. Nienhuys H-W, van der Stappen A. A computational technique for interactive needle insertions in 3d nonlinear material; IEEE Intl. Conf. on Robot. and Autom; 2004. May, pp. 2061–2067.

20. Park W, Liu Y, Zhou Y, Moses M, Chirikjian GS. Kinematic state estimation and motion planning for stochastic nonholonomic systems using the exponential map. Robotica. 2008 July;vol. 26:419–434. [PMC free article] [PubMed]

21. Reed KB. Compensating for torsion windup in steerable needles. Proc. 2nd Biennial IEEE/RAS-EMBS Intl. Conf. on Biomedical Rob. and Biomechatronics; Scottsdale, Arizona. 2008. [PMC free article] [PubMed]

22. Reed KB, Kallem V, Alterovitz R, Goldberg K, Okamura AM, Cowan NJ. Integrated planning and image-guided control for planar needle steering. Proc. 2nd Biennial IEEE/RAS-EMBS Intl. Conf. on Biomedical Rob. and Biomechatronics; Scottsdale, Arizona. 2008. [PMC free article] [PubMed]

23. Saha M, Isto P. Manipulation planning for deformable linear objects. IEEE Trans. Robot. 2007 Dec.vol. 23(no. 6):1141–1150.

24. Webster RJ, III, Cowan NJ, Chirikjian G, Okamura AM. Nonholonomic modeling of needle steering; 9th International Symposium on Experimental Robotics; 2004. Jun,

25. Webster RJ, III, Memisevic J, Okamura AM. Design considerations for robotic needle steering; IEEE Intl. Conf. on Robot. and Autom; 2005. Apr, pp. 3588–3594.

26. Xu J, Duindam V, Alterovitz R, Goldberg K. Motion planning for steerable needles in 3d environments with obstacles using rapidly-exploring random trees and backchaining; Proc. IEEE Conf. on Automation Science and Engineering; 2008. Aug, pp. 41–46.

27. Yamada H. Strength of Biological Materials. Baltimore: Williams & Wilkins; 1990.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |