Backward Planning for a Multi-Stage Steerable Needle Lung Robot

I. Introduction

LUNG cancer is responsible for the most cancer-related deaths in the United States [1]. Physicians can identify suspicious nodules in a patient’s lung using medical imaging, but it is currently not possible to definitively diagnose those nodules with sufficient accuracy via imaging alone. While biopsies are the gold standard for diagnosis, existing approaches are not always sufficient: (i) transoral approaches (i.e., through the mouth) cannot reach many regions in the lung as these approaches are restricted to the airways accessible by a bronchoscope; and (ii) percutaneous approaches that insert a needle through the chest wall carry high risks such as lung collapse [2]. As early diagnosis increases a patient’s chance of survival, a reliable procedure bearing lower risk could be transformative for lung cancer diagnostics and patient care.

To enable a minimally invasive lung biopsy procedure that is both effective and safe, we are developing a new steerable needle lung robot [3]. The robot combines the advantages of both the transoral and percutaneous approaches by combining low-risk transoral deployment with the ability to reach targets in the peripheral lung using a steered needle. Our robot includes three stages that must be sequentially deployed during a procedure. First, a physician deploys a bronchoscope into the patient’s airways. Then, the physician utilizes a sharpened piercing tube to pierce through the airway wall and into the lung parenchyma. Subsequently, the system extends a steerable needle from the piercing tube and into the lung parenchyma, where it automatically steers to the target location.

In this paper, we introduce a new motion planner for the steerable needle lung robot that integrates planning for all three robot stages to enable fast robot deployment to a target location. In this process, it is crucial to minimize patient risk. Therefore, the planner has to avoid piercing critical anatomic structures–such as large blood vessels or airways–which could lead to severe complications. Fig. 1 shows planned trajectories in lung anatomy extracted from a preoperative 3D CT scan. Motion planning for the steerable needle lung robot requires creating plans for each of the three stages, where each stage’s start pose depends on the previous stage’s end pose. When computing plans, prior planning algorithms considered the stages in the order of their deployment. Based on the observation that planning from a point (the target) to a region (the bronchial tubes) is faster than planning from a region to a point, our new planner considers the stages in reverse order to achieve significant computation time reduction.

Specifically, prior work [4] used a forward planning approach that first randomly samples a bronchoscope and piercing tube pose. From the sampled piercing tube’s tip pose the needle planner then starts to iteratively build a search tree toward the target. In case no plan is found from this randomly sampled start pose, the pose is rejected and the full process is repeated. This includes repeated construction of new needle search trees, which is very computationally expensive.

In this work, we propose the inversion of the problem by constructing a needle search tree from the target location and by integrating planning for all three stages into one search tree. Once this tree approaches a region near the airways, the algorithm attempts to locate a straight piercing trajectory that leads into the center of an airway reachable by the bronchoscope. If no successful piercing pose is found, the search tree is grown further, extending existing paths until the planner is successful. This is much more efficient than restarting the entire planning process and relinquishing all previous exploration. This process relies on pre-computing a model of regions reachable by the bronchoscope inside the airway structure. Furthermore, we introduce two speedup strategies for the planner: one is based on combined geometric constraints of all three robot stages, and the other biases search tree growth toward target regions.

We evaluate our new motion planner for a three-stage robot in simulation using hardware specifications of the robot shown in Fig. 2. We use anatomy from ex-vivo and in-vivo porcine lungs, which provide a good model for human lungs [5]. On average, our approach finds a plan to a random target more than five times faster than prior work.

II. Related Work

Chest CT scans are the standard for lung cancer imaging [6]. To plan safe deployment of our steerable needle lung robot, we require an accurate geometric representation of the planning environment, which we obtain by segmenting a patient-specific CT scan. A large number of algorithms for segmentation of different lung structures are available as summarized by Van Rikxoort et al. [7]. The segmentation approach used in this work is described in [8].

Graph and tree data structures are common ways for describing airway branching structures [9]. Applications that use such tree representations include motion planning for bronchoscopies [10], branch point matching for lung registration [11], and trachial stenosis diagnostics through imaging [12]. Popular methods for extracting such a tree structure are thinning or skeletonization algorithms. We apply the parallelized thinning algorithm by Lee et al. [13] in an implementation as developed by Homann et al. [14] to our airway segmentation. Skeleton representations can be parsed voxel by voxel to extract a directed graph representation [15]. The links in such a graph are often modeled as cylinders representing the varying radii of airway branches [16], [17], [18]. We use such a representation to model the region in the airways that can be reached by the bronchoscope.

Robotic steerable needles are a promising approach for minimally invasive procedures to operate in parts of the body not safely reachable with conventional instruments [19]. Bevel tip steerable needles apply asymmetric force on tissue when they are inserted, which causes them to curve in the direction of the bevel [20]. By axially rotating the needle during insertion, different curving directions can be achieved. Using duty cycling with axial rotation results in varying insertion curvatures [21].

Different motion planning approaches for steerable needles have been previously introduced, including fractal-based planning [22] and optimization methods [23]. Sampling-based motion planning is another planning approach that has been proven to be effective in many different high-dimensional applications [24]. The Rapidly-exploring Random Tree algorithm (RRT) [25] provides an efficient way to plan a path from a start to a specific target by constructing a search tree structure. One source of difficulty in sampling-based motion planning for steerable needles is the needle’s non-holonomic constraints (i.e., its minimum radius of curvature). Connecting two samples in SE(3) under these constraints requires solving a difficult two-point boundary value problem [26]. Many sampling-based algorithms such as RRT* [27] that require connecting two arbitrary configurations can therefore not be directly applied to needle planning. One solution to this problem is to first find a plan and then smooth it such that it complies with the constraints [28]. Our work builds on a sampling-based method by Patil et al. [29], which creates a search tree fully compliant with non-holonomic constraints. This needle planner was extended to a 3-stage motion planner for robot lung procedures by Kuntz et al. [4]. To our knowledge, this is the only existing work integrating a needle planner into a multi-stage planning process. In addition, needle planning methods exist that explicitly consider deformations [30], [31] or consider uncertainty [32], [33]. Berenson et al. [34] present an RRT approach for planning toward a goal region instead of a single target which is based on inverse kinematics and bidirectional planning from the goal region and the start location simultaneously. This method is not directly applicable to our planning problem as our goal regions are large and a bidirectional search would therefore be very inefficient. Nevertheless, it inspired our strategy for target region biasing, which encourages RRT growth toward the airways and results in a significant speedup.

III. Problem Definition

Each of the three robotic system stages has physical constraints that need to be considered in the planning process. We visualize these constant hardware limitations in Fig. 3. A bronchoscope has a diameter d_bronch that determines how far it can be extended into the airways. A piercing tube is represented by a maximum insertion length l_tube beyond the bronchoscope tip, a diameter d_tube, and a maximum piercing angle θ_tube that signifies how much the piercing tube’s direction can deviate from the bronchoscope tip’s medial axis. A steerable needle’s hardware parameters are defined by its minimum steering radius of curvature r_curve, cross-sectional diameter d_needle, and maximum insertion length l_needle.

The planning algorithm produces a deployment path ρ to the target location p_target based on the consecutive system stages: ρ = [ρ_bronch, ρ_tube, ρ_needle]. We define a path of stage s as ρ_s = [Q_s,1, …, Q_s,I], an ordered list of $I\in \mathbb{N}$ configurations. The bronchoscope, piercing tube, and needle paths consist of $L,M,N\in \mathbb{N}$ configurations, respectively. Configurations are described by poses, denoted as 4 × 4 homogeneous transformation matrices

where ${\mathbf{p}}_{s,i}\in {ℝ}^3$ is the position and R_s,i ∈ SO(3) is the orientation relative to a world coordinate frame.

The first piercing tube position p_tube,1 is equal to the final bronchoscope position p_bronch,L and deviates from the bronchoscope tip’s orientation by up to θ_tube degrees. Q_tube,pierce denotes the pose in which the piercing tube pierces out of the airways, and the last piercing tube pose Q_tube,M is equal to the first needle pose Q_needle,1. The needle plan ends at the target: p_needle,N = p_target.

Each point along a plan has an associated cost $C:{ℝ}^3\to [0,\infty )$. The cost for a plan ρ is defined as

where s ∈ [0, 1] and $p:[0,1]\to {ℝ}^3$ is a mapping of configurations along ρ into 3D space. A plan is only valid if its total cost is C(ρ) < ∞, which signifies that it pierces no critical anatomical structures.

Finding an optimal motion plan in this setting can be expressed as an optimization problem:

Subject to:

where ρ* is an optimal plan. The general inequalities g_bronch(ρ_bronch), g_tube(ρ_tube), and g_needle(ρ_needle) represent the stages’ kinematic constraints, in particular the bronchoscope’s ability to reach areas in the airways limited by its diameter, the maximum piercing angle, and the needle’s minimum radius of curvature.

IV. Methods

Planning a safe robot path toward a target includes sampling configurations for the bronchoscope, piercing tube, and steerable needle while avoiding critical obstacles and aiming to find the lowest cost plan. When planning backwards from the target, the planner attempts to find a path to the reachable bronchoscope region. Moreover, the planner has to consider hardware constraints such as a limitation for piercing angles and the steerable needle’s minimum radius of curvature. In the following, we will describe each element of the planning process in detail.

A. Anatomical Environment Representation

To ensure a safe procedure, collisions with critical structures need to be avoided. Using the method in [8], we create a cost map I_C where each voxel is assigned a cost value:

$$C(p)=\{\begin{array}{ll}\infty \hfill & \text{if}\text{voxel}p\in I_A\cup I_B\cup I_R\hfill \\ [0,1]\hfill & \text{otherwise}\hfill \end{array}$$

where I_A, I_B, and I_R are sets of voxels representing the airway, large blood vessels and the region outside the lung, respectively. Collisions with smaller blood vessels are common during clinical procedures, but they should still be avoided when possible. We choose a cost metric that reduces the risk of piercing vessels during robot deployment. To compute this cost, we extract smaller blood vessels I_V, which are assigned values [0,1] in I_C. Higher values indicate a higher risk of piercing through small vessels at each voxel [8]. Our framework also supports other cost metrics to evaluate the path quality such as clearance from obstacles and path length. We determine all locations inside the airways that are reachable by the bronchoscope tip. We model this region by first extracting the airways’ medial axes by applying a 3D thinning algorithm [13] to binary image I_A. A resulting skeleton binary representation is shown in Fig. 4a. This representation is parsed to extract a directed tree representation. We construct such a tree T(n, e), where each branch point in the skeleton is a node n and each branch is an edge e. We parse all voxels in the skeleton image starting from a manually selected voxel in the trachea. Bronchoscope deployment is limited by its diameter in comparison to the diameter of corresponding airway branches. To determine bronchoscope reachability, we create a simplified airway model by approximating each branch by a cylinder. For each edge e in the tree structure we connect the start and end node locations to form a cylinder’s medial axis and iteratively test positions along this axis for radii inside the airway. For each position we determine the largest circle around it in the plane perpendicular to the medial axis that is fully contained inside the airways. To avoid overestimation of the airway diameter, the smallest radius found along the axis is chosen to be the representing cylinder’s radius.

Fig. 4.

Open in a new tab

Positions in the airway that can be reached by the bronchoscope are determined in a recursive way utilizing tree T. Starting in root node n₁ located in the trachea, we conduct a depth-first search following each edge until its associated diameter becomes smaller than the bronchoscope’s diameter. The result of this search is a set of cylinders B = {b₁, …, b_n} that represents all reachable airway branches. A cylinder b ∈ B representing the reachable part of edge e is defined as $b=\left\{{\mathbf{p}}_b,{\mathbf{p}}_b^{\prime },{\mathbf{v}}_b,l_b,r_b\right\}$. The start position ${\mathbf{p}}_b\in {ℝ}^3$ is equivalent to the coordinates of n_i, the node the edge starts from. The end position ${\mathbf{p}}_b^{\prime }\in {ℝ}^3$ is the last reachable position along e. The cylinder’s direction ${\mathbf{v}}_b\in {ℝ}^3$ is defined as the vector connecting the edge’s start node and end node v_b = n_j −n_i. The cylinder’s length is $l_b={\Vert {\mathbf{p}}_b^{\prime }-{\mathbf{p}}_b\Vert }_2$. The radius r_b is the smallest radius determined along the reachable part of e subtracted by 0.5d_bronch to represent regions inside the branch that the center of the bronchoscope tip can reach. In addition, we define the piercing region P as a cylinder set representing the region the piercing tube tip can reach. Here, we extend the length of each cylinder in B by the piercing tube’s maximum length l_tube and the radius of each cylinder by l_tube sin(θ_tube) to define an outer boundary of all points the piercing tube can reach if extended to its maximum. The cylinder sets are visualized in Fig. 4b. It is highly discouraged to puncture the pleura between lung lobes due to potential lung collapse. Therefore, we always restrict the planning process to the lobe the target is located in.

Algorithm 1.

Backward Planning

Open in a new tab

Algorithm 2.

Sampling a Configuration

Open in a new tab

V. Evaluation

We evaluated our algorithm based on CT scans of porcine lungs, which closely match the anatomy of human lungs [5] and will be used in pre-human trials of the steerable needle lung robot. We acquired CT scans of porcine lungs using a Siemens Biograph mCT scanner (Siemens Healthineers, Erlangen, Germany). Scans were reconstructed with a voxel size of (0.6 × 0.6 × 1.0) applying a pulmonary B70F kernel. We used scans of two inflated ex-vivo porcine lungs and one in-vivo porcine lung with 150 randomly selected targets in each. In the two ex-vivo lungs we sampled targets in 5 lobes (right caudal, left caudal, middle, accessory, and left cranial lobes) that correspond to the 5 lobes in a human lung [5]. The in-vivo lung is significantly smaller (approximate lung volume 1.7 liters) than an adult human lung. Therefore, only the two caudal lobes are reachable with our system and we sampled targets in these lobes only. The in-vivo scan was taken during a single breath hold.

We only selected targets with a safety margin of at least 1 mm in addition to the minimum distance to the closest obstacle voxel as defined in Equation 1. We assumed the bronchoscope’s diameter is d_bronch = 6mm. The piercing tube had a diameter of d_tube = 2mm and a maximum insertion length of l_tube = 50mm. We chose a maximum piercing angle of θ_tube = 45 degrees based on experimental capabilities. The needle’s diameter was d_needle = 1.0mm with a maximum insertion length of l_needle = 120mm. A conservative measure of its minimum radius of curvature is r_curve = 100mm [36], [3], [37].

We compare our backward planning approach for a 3-stage robot against prior work, which used a forward planning method [4]. The prior work forward planner randomly samples bronchoscope and piercing tube poses that serve as a start pose for building an RRT for the needle stage. Needle configuration samples for the RRT are geometrically constrained to the needle’s reachable workspace. The forward planner employs goal biasing by adding p_target to the RRT with probability Pr_bias. This is different from our biasing strategy, which consists of adding random samples from within the reachable airway region to the search tree. For fair comparison we set Pr_bias = 0.1 for both algorithms and we set Pr_start = 0.05.

Additionally, we set a timeout t = 3 sec for the needle planning stage in the forward planner to avoid wasting compute time if the needle planning problem is infeasible for a specific combination of start pose and target position. The backward planner does not require timeouts for any of its subroutines as it keeps extending the RRT from the target until a solution is found.

For each target, we ran the planning process 10 times for a duration of 60 seconds for each algorithm, respectively. All simulations were run on a 3.2GHz 16-core Intel Xeon E5–2667 CPU with 64GB of RAM. We implemented all motion planner variants discussed in this paper based on the Motion Planning Template (MPT) library [38], which provides functionality for parallelized search tree construction. For efficient collision detection we used the nigh library [39] for nearest neighbor search. We evaluated the algorithms with respect to three different performance metrics: the ability to find a plan to a variety of targets, computational speed, and path costs.

First, we analyzed planning success for each target location. A successful planning process finds at least one plan within 60 seconds. The ex-vivo scenarios are simpler planning problems since blood vessel segmentation is limited as vessels are no longer filled with blood. Furthermore, the inflated ex-vivo lungs are significantly larger than the in-vivo lung providing more space for the needle to curve towards a target. In the two ex-vivo scenarios our backward algorithm found at least one plan for 290 out of 300 target location, whereas the forward planner only found solutions for 234 targets. In the in-vivo scenario our backward planner found plans for 128 out of 150 targets, while the forward planner only found plans to 61 targets. There is no target the forward algorithm found a plan for that our backward algorithm did not. Fig. 7 depicts that most targets for which only the backward algorithm was successful are located close to the airways. This illustrates a fundamental disadvantage of the forward planner: finding a plan to a target close to the airways is highly dependent on the highly random poses sampled in the first two stages.

Second, we were interested in the duration of the planning process as it is critical to quickly find a safe plan during a procedure. We define success rate as the percentage of targets for which a plan has been found during a single run evaluated for all 450 targets. We recorded how the success rate changes throughout the time interval of 60 seconds and display an average over 10 such runs in Fig. 8a. It is evident that even the basic backward algorithm is significantly faster at finding plans than the prior forward approach. In addition, we see that both the geometric constraints and the biasing speedup strategies contribute to faster planning. Combining both strategies resulted in finding plans for 92.53% (±0.37% standard deviation) of targets in under 60 seconds. The backward approach with speedup takes on average 2.48 seconds to find the first plan per target, whereas the forward approach takes 13.90 seconds taking into account only targets with a plan found at least once by both algorithms.

Third, we were interested in minimizing path cost, which is equivalent to minimizing the piercing of smaller blood vessels. Thus, plan cost is a measure of plan safety. During each planning process, the cost to reach the target is updated whenever a less expensive plan is found. We computed the average cost progression across all targets for each of the 10 runs. For each time step we determined all targets for which at least one plan was found and we computed the average of their current costs, resulting in the average cost progression over time for one run. We computed the average and standard deviation across all 10 runs, as shown in Fig. 8b. For fair comparison we considered only targets for which all algorithm variants found at least one plan. Fig. 8b shows the average cost progression and its standard deviation across all runs. It can be seen that our backward algorithm found less expensive plans significantly faster than the forward algorithm and that both speedup strategies contributed to better results.

VI. Conclusion

In this work, we introduce a new planning approach for a multi-stage steerable needle lung robot designed for lung interventions. Our approach is based on backward planning, which allows for more efficiency. In addition, we introduce two speedup strategies. One strategy is based on imposing geometric constraints on the search space and the other one biases tree extension toward regions of interest. We compared our approach to an existing forward planning method. For three different anatomical scenarios we demonstrated (i) that our approach is more likely to find a path to a target; (ii) it does so faster; and (iii) it produces lower-risk paths in shorter time. The main advantage of our new backward approach over the forward approach is that it performs the most complex part of the planning process - the creation of a needle search tree - only once instead of repeating it many times.

While our backward planning approach shows promising results, there are still some features to be optimized. For instance, to accelerate the process of finding lower cost plans, locally optimizing existing plans for lower cost is an option to be explored.

Similar to previous experiments [40], we are planning on evaluating the algorithm by performing a biopsy procedure in an ex-vivo porcine lung. Clinical translation introduces new challenges including but not limited to deformable image registration of a pre-operative scan to the lung and modeling tissue motion from breathing and heartbeats. From a planning perspective, rapid re-planning for some or all of the robot stages may be required to adjust for uncertainty. We hope that the increased efficiency of our 3-stage planner will be an enabler to successful clinical experiments.