PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
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

Feedback Control for Steering Needles Through 3D Deformable Tissue Using Helical Paths

Abstract

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.

I. INTRODUCTION

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%.

Fig. 1
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/ ...

II. RELATED WORK

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.

III. PROBLEM STATEMENT AND MODELING ASSUMPTIONS

A. Problem Statement

We wish to place the needle tip at a certain goal location pg 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 pg 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 W [subset, dbl equals] R3 in which the needle is permitted to travel. In this work, W is taken to be either all of R3 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 [var phi]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.

B. Rigid Tissue Motion Model

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 [var phi](t). In screw notation, the instantaneous velocity of the frame is given by:

V^(t)=[0φ(t)00φ(t)0υ(t)/r00υ(t)/r0υ(t)0001]

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 [var phi](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(dV) [12].

Fig. 2
The coordinate frame attached to the needle tip. Figure adapted with permission from [13].

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

IV. FEEDBACK CONTROLLER WITH HELICAL PATHS

A. Controller Framework

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

  1. Propose. Generate a set P of proposal trajectories that start from the current estimated needle state.
  2. Select. Find the trajectory in P with control [var phi](d) that achieves the minimal distance to the target.
  3. Execute. The needle is then inserted with twists according to [var phi](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).

Fig. 3
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 P. Convergence is difficult to prove, and depends on the selection of P. Empirical results in Section VI suggest that setting P 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 P (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 [var phi]. We let the set of proposal trajectories P contain all such helices for [var phi] [set membership] [−[var phi]max, [var phi]max]. Using a visualization of P, we argue that the controller converges for sufficiently distant targets.

B. Constant-Twist-Rate Helical Paths

Let t[var phi](d) be the needle tip trajectory followed under a constant twist rate [var phi]. 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

h(d)=[acos((d/a)cosθ)asin((d/a)cosθ)dsinθ].
(1)

Here we derive the helix parameters a, θ and a rigid transformation A such that Ah(d) = t[var phi](d). Denote hx,hy,hz, as the coordinate axes of A, and denote ho as its origin. We proceed by matching derivatives:

Ah(0)=tφ(0)ahx+ho=[0,0,0]T
(2)

Ah(0)=tφ(0)cosθhy+sinθhz=[0,0,1]T
(3)

Ah(0)=tφ(0)cos2θ/ahx=[0,1/r,0]T
(4)

Let ω = [1/r, 0, [var phi]]T be the angular velocity of the needle tip. The axis of rotation must be aligned with the helix axis, so hz = ω/‖ω‖. Let q=ω=1/r2+φ2. Taking the dot product of both sides of (3) with hz, we get

sinθ=φ/q   and   cosθ=1/(rq).

From (4) we see

a=rcos2θ=1/(rq2)   and   hx=[0,1,0]T.

Setting hy orthogonal to hx and hz, we get

hy=[φ,0,1/r]T/q.

Finally, from (2) we have ho = [0, −a, 0]T.

In summary, letting tan θ = r[var phi], we have

tφ(d)=[asinθsin((d/r)cscθ)+dsinθcosθacos((d/r)cscθ)aacosθsin((d/r)cscθ)+dsin2θ]
(5)

Where a = r cos2 θ = r/(1/r2 + [var phi]2) is the helix radius, and [cos θ, 0, sin θ]T is the helix axis. The tip moves at the speed sinθ=φ/1/r2+φ2 with respect to the helix axis.

C. Constant-Twist-Rate Reachable Set

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

Fig. 4
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].

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

This visualization helps understand why the controller should converge. The projection of R(∞) 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 R(∞) is rotated by a full turn about any axis, it will sweep out a large portion of R3. 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.

D. Boosting Accuracy with Alternating-Twist Maneuvers

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

Fig. 5
Alternating-twist maneuvers for various twist rates. Units in rad/r.

Intuitively, the maneuver allows the needle to travel approximately straight forward. After each full turn, the needle tip is displaced by (2πr sin θ cos2 θ, 0, 2πr sin2 θ cos θ) where θ = arctan r[var phi]max. Orientation is unchanged. Executing both turns places the needle tip at (0, 0, 4πr sin2 θ 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 P is augmented with this maneuver.

V. FAST TRAJECTORY SELECTION

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.

A. Branch-and-Bound Closest Point Computation

As can be seen in Figure 4, the reachable set R([var phi]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([var phi], d) = ‖t[var phi](d) − p‖ over the space S = {([var phi], d) | |[var phi]| ≤ [var phi]max, d > 0}. Here p denotes the coordinates of pg relative to the needle tip frame. We use an auxiliary function fL(R), described in Section V-C, that computes a lower bound on the function value attained in a region R [subset, dbl equals] S. 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 [var phi][open star] and insertion distance d[open star] that give the minimum value of seen thus far. If fL(R) is found to be larger than , then subregion R does not contain the optimum and can be safely pruned. The process continues until 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 [var phi][open star] and insertion distance d[open star] from the previous iteration, and initialize = f([var phi][open star], d[open star]). 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 S by selecting n = 100 possible values for [var phi]. For a more uniform coverage of the reachable set given a fixed n, we discretize evenly in the θ = arctan r[var phi] coordinate. To include other proposal trajectories in the search, we divide them into helix segments, and add the segments to S. Then, for each helix, we narrow S 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.

B. Initial intervals for helix closest 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 hc(u) = [a, u, ub]T. Let pc = [rp, θp, zp]T.

First, the helix attains z-coordinate zp at u = zp/b. The closest point to p must be within a half turn of this parameter; that is, |u[open star]zp/b| ≤ π (Figure 6a). The proof is by contradiction: if u[open star] were actually outside this interval, then a full turn around the helix could move the z coordinate closer to zp while keeping x and y fixed. Therefore u[open star] is not optimal.

Fig. 6
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 = zp plane. (b) The interval containing the z-range of the intersection between a sphere of radius dR centered ...

Second, we examine the cylinder bounding the helix. If dR is the distance from p to any point on the reachable set (say, h(zp/b)), then the closest point must lie within a ball of radius dR around p. Projecting the intersection of the bounding cylinder with the sphere onto the z axis, we see that the z-coordinate z[open star] of the closest point must satisfy |zpz|dR2(rpa)2 (Figure 6b). This bound performs well when dR is close to the actual helix-point distance.

C. Lower bounding the helix-point distance

Now we describe a lower bound fL for u [set membership] [, ū] and [var phi] given. We compute the distance from the point to a patch (in cylindrical coordinates) that contains the portion of the helix [a, a] × [, ū] × [u̠b, ūb]. Let C(x, [a, b]) denote the clamping function max(a, min(b, x)), and implicitly assume operations in the θ coordinate are modulo 2π. Then the closest point to this patch is given by

qc=[aC(θp,[u¯,u¯])C(zp,[u¯b,u¯b])]
(6)

The Cartesian distance ‖qp‖ is a lower bound.

D. Workspace Bounds

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 W. 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, d], where d 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 [set membership] I using interval computations. If the box lies within W, 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.

E. Running Time

Since no interval will have size less than [set membership], the number of intervals examined by the branch and bound algorithm is no more than [left ceiling]2L[right ceiling], 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 [left ceiling]rn[right ceiling]. 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.

VI. SIMULATIONS IN RIGID TISSUE

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.

A. Accuracy Without Perturbations

For comparison, we instantiate a reference controller that uses the following parameters: a refresh occurs every 2%r of insertion distance, maximum twist rate [var phi]max = 10π rad/r, and the workspace W is all of R3.

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 pg 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.

Fig. 7
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:

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

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] × [−3r, 3r] × [1.5r, 6r], which is the rightmost 3/4 of the x = 0r 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.

TABLE I
Accuracy under varying controller parameters.

B. Under Perturbations and Modeling Errors

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.

TABLE II
Accuracy under simulated perturbations.

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.

  1. 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.
  2. 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°.
  3. 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.44r and σo = 11.1° per each r distance the needle is inserted.
  4. Curvature Estimation Errors: We set the actual radius of curvature to [r with tilde] = 0.8r and [r with tilde] = 1.2r, while the estimated radius remains at r.
  5. 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 [theta w/ dot above]b = [var phi], and [theta w/ dot above]t = 0 if |θb − θt| < ηL, or [var phi] otherwise. We set η = 6°/r, which is consistent with the experimental data in [21] assuming r=5 cm.

VII. SIMULATIONS IN DEFORMABLE TISSUE

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 W.

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).

Fig. 8
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. ...
Fig. 9
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.

VIII. CONCLUSION

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.

ACKNOWLEDGMENT

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.

REFERENCES

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.