"Trajectory planning method for reduced patient risk in image-guided neurosurgery: concept and

preliminary results" Reuben R. Shamir, Leo Joskowicz, Luca Antiga, Roberto I. Foroni, and Yigal Shoshan. Proc. SPIE, Vol. 7625, 726250I (2010); DOI:10.1117/12.843991 Copyright Medical Imaging (2010) Society of Photo-Optical Instrumentation Engineers. One print or electronic copy may be made for personal use only. Systematic reproduction and distribution, duplication of any material in this paper for a fee or for commercial purposes, or modification of the content of the paper are prohibited. DOI abstract link : http://dx.doi.org/10.1117/12.843991

Trajectory planning method for reduced patient risk in image-guided neurosurgery: concept and preliminary results Reuben R. Shamir1*, Leo Joskowicz1, Luca Antiga2, Roberto I. Foroni3, and Yigal Shoshan4 1 School of Engineering and Computer Science, The Hebrew University of Jerusalem, Israel. 2 Bioengineering Department, Mario Negri Institute for Pharmacological Research, Ranica, Italy. 3 Department of Neurosurgery, Verona University Hospital, Verona, Italy. 4 Department of Neurosurgery, Hadassah University Hospital, Ein-Karem, Jerusalem, Israel. ABSTRACT We present a new preoperative planning method to quantify and help reduce the risk associated with needle and tool insertion trajectories in image-guided keyhole neurosurgery. The goal is to quantify the risk of a proposed straight trajectory, and/or to find the trajectory with the lowest risk to nearby brain structures based on pre-operative CT/MRI images. The method automatically computes the risk associated with a given trajectory, or finds the trajectory with the lowest risk to nearby brain structures based on preoperative image segmentation and on a risk volume map. The surgeon can revise the suggested trajectory, add a new one using interactive 3D visualization, and obtain a quantitative risk measure. The trajectory risk is evaluated based on the tool placement uncertainty, on the proximity of critical brain structures, and on a predefined table of quantitative geometric risk measures. Our preliminary results on a clinical dataset with eight targets show a significant reduction in trajectory risk and a shortening of the preoperative planning time as compared to the conventional method. Keywords: Image-Guided Therapy, Therapy Planning, Neurosurgical Procedure, Neuronavigation

1. INTRODUCTION Precise targeting of tumors, lesions, and anatomical structures with a probe, needle, catheter, or electrode inside the brain based on pre-operative CT/MRI images is the standard of care in many keyhole neurosurgical procedures. The procedures include tumor biopsies, treatment of hydrocephalus, aspiration and evacuation of deep brain hematomas, Ommaya catheter insertion, deep brain stimulation, and minimal access craniotomies, among many others [1, 2]. In all cases, misplacement of the surgical instrument may result in non-diagnostic tissue samples, hemorrhage and/or severe neurological complications [3]. Trajectory planning in current image-guided keyhole neurosurgery is usually performed manually on 2D cross-sections of the patient pre-operative CT/MRI head images. This requires the surgeon to mentally reconstruct complex 3D brain structures and functional zones, and derive the spatial relations between them. The evaluation of the treatment consequences of various candidate trajectories is thus a complex, challenging, and timeconsuming task even for experienced neurosurgeons. While volume visualization and spatial segmentation of critical brain structures are sometimes used to help the neurosurgeon with the spatial perception, the insertion trajectory is selected manually and does not include any quantitative measures or trajectory-specific visualization of nearby critical structures. The resulting trajectory is thus surgeon-dependent and may not be optimal. Several methods have been proposed to better assess and reduce risk in image guided neurosurgery [4-7]. Some of them automatically compute the linear trajectory with the smallest risk, where the trajectory risk is defined as the weighted sum of voxels in the segmented image that are intersected by the defined trajectory [4, 7]. Brunenberg et al. [6] argue that this cumulative trajectory risk computation can be misleading and suggest instead to compute the maximum risk value among voxels along the tool trajectory. The Euclidean distance of the trajectory from critical brain structures is used to compute the risk of each voxel. Their method outputs tens to hundreds of trajectories associated with distances above a predefined given threshold. Although the method significantly reduces the number of possible trajectories, it * [email protected]; phone +972-77-200-5991; http://www.cs.huji.ac.il/~caslab/site/

(a)

conventional and computed trajectories

(b) conventional method

(c) our method

Fig. 1. Illustration of trajectory planning: (a) conventional (manually selected) and automatically computed trajectories for a given target -- the green crosses indicate the trajectories entry points. Detail of (b) the conventional trajectory, and (c) the trajectory computed with our method showing the uncertainty location cylinder and closest blood vessel distance. The uncertainty cylinder diameter represents the estimated tool trajectory uncertainty. Our method yields a trajectory that is further from blood vessels and hence has a lower risk of damaging them. The tool trajectory localization uncertainty visualization supports the examination of nearby critical structures.

still leaves a considerable amount of manual work. Moreover, no quantitative results are reported. Popovic et al. [5] use a distance map for the risk computation and show how to compute the optimal trajectory for nested cannula devices using the A* algorithm. In their application, the number of cannula configurations is large, so the A* algorithm may be appropriate, but is computationally time-consuming. However, in our application, the entry and target points directly define the trajectory, so simpler, faster, and more robust algorithms can be used instead. Additionally, their method provides no convenient method for modifying the computed trajectory or for assessing its risk. We have developed a new preoperative trajectory planning method to quantify and reduce risk in image-guided keyhole neurosurgery that aims at overcoming the limitations described above. The method automatically computes the straightline trajectory that is farthest from critical structures that are segmented on the preoperative head image. The surgeon can then revise the suggested trajectory, and add or edit trajectories with visual feedback and updated risk information (Fig. 1). A novel feature of our method is that it incorporates interactive 3D visualization of critical structures and of surgical tool placement uncertainty. The computed trajectory is presented to the physician along with a „risk card‟ that includes quantitative risk measures such as the length of the trajectory and distance between the trajectory and closest blood vessels (Fig. 2). Our preliminary results show a reduction in trajectory risk and a shortening of planning time compared to the conventional method.

2. METHODS 2.1 Method overview We propose the following preoperative planning workflow for image guided keyhole neurosurgery (Fig. 3). Initially, the neurosurgeon selects the target location on a preoperative Magnetic Resonance Angiography (MRA) image (1). Then, the head outer surface is computed automatically (2) and the neurosurgeon interactively defines the surface region on which the entry point should reside (3). Afterwards, the relevant structures (blood vessels and no-go zones) are segmented (4, 5) and assigned a risk value based on the potential damage associated with them. The input segmentations are automatically combined into a single volume, called the risk volume, in which voxels are associated with values representing their risk (6). A trajectory risk value is then computed for each trajectory automatically (7). The candidate trajectories are then stored in a priority queue so that they can be efficiently retrieved in increasing order of risk (8). The best trajectory is retrieved first, and a „risk card‟ is computed and presented to the user (Figs. 1,2).

Fig. 2. The risk card shows the values of various geometrical risk parameters identified by a senior neurosurgeon. The first column is the trajectory ID; the following columns show the trajectory length, the distance of trajectory to the closest blood vessel, the distance of trajectory to the closest „no-go‟ regions, the distance of the selected target to its closest blood vessel, and the distance of the selected target to its closest no-go zone, respectively (all values are in mm). The table conveniently supports the direct comparison of candidate trajectories risk parameters.

The risk card is a summary of geometrical parameters regarding critical structures and planned trajectories. It provides the neurosurgeon with valuable information for the assessment of an intervention‟s risk and enables the direct quantitative comparison between candidate trajectories. The relevant risk parameters were identified by a senior neurosurgeon (last author) as valuable measures for assessing trajectory‟s risk and include trajectory‟s length, and distances of trajectory, target and entry to closest blood vessels and other manually defined restricted and “no-go” zones. Trajectories can be added and modified based on the 3D visualization of their localization uncertainty and predefined critical structures (Fig. 1). The visualization of tool localization uncertainty eases the identification of cases where a planned line trajectory does not cross a critical structure, including the placement uncertainty, and indicates the possible damage to the critical structure. We model the localization uncertainty with a cylinder as in the Medtronic cranial application (Medtronic, USA). However, other uncertainty geometries can be incorporated as well. In the following section, we describe our method for the computation of risk volume and trajectory risk. Then, we briefly describe the automatic blood vessels segmentation method that we used in the experimental study. 2.2 Computation of risk volume and trajectory’s risk The risk volume associates each voxel with the estimated cost incurred in penetrating that voxel with a surgical tool, thereby possibly damaging the associated tissue or organ. The risk volume is generated based on two key guidelines: 1) the risk value that is assigned for each voxel is directly related with the estimated consequences and severity of the damage to the corresponding brain tissue or organ – severe complications and high morbidity regions are assigned with a higher risk values than tissues with minor and reversible complications, and; 2) voxels near critical structures are assigned with high risk values to reflect the intrinsic localization error of the procedure, be it freehand, frame-based stereotaxy, or with an image-guided surgery system [9]. Therefore, voxels that are closer to a critical structure are associated with a higher risk value than those that reside further from them. One can envisage several ways to model the risk to fulfill the above guidelines. We describe next our first implementation of them. The input is a set of critical structures for which insertion of a surgical tool is forbidden or not desired, S  {S1 , S2 , , S p } , and their associated risk values R  {r1 , r2 , , rp } . A structure Si is a segmented image, a risk value ri is a non-negative scalar. In our implementation, each voxel in the risk volume is defined as:



   dist  x,Sk    k

riskVolume  x   max 

rk

(1)

1 Target definition

2

MRA image

4

5

Head surface extraction

Blood vessels segmentation

3

Entry points surface definition

7

Trajectories risk computation

8 6

Risk volume computation

Trajectory selection

No-Go zones segmentation

Fig. 3. Preoperative planning workflow: initially, the neurosurgeon selects the target location on the preoperative MRA image (1). Then, the head outer surface is computed automatically (2), and the neurosurgeon interactively defines the surface regions on which the entry point should reside (3). Next, the blood vessels are automatically segmented (4), and the restricted regions (No-Go zones) are interactively segmented by the neurosurgeon. The risk volume is then automatically computed on the segmented images (6). The trajectories risks are then computed automatically for all the possible trajectories based on the given risk volume, entry points surfaces, and target location (7). Finally, the neurosurgeon selects the preferred trajectory based on the results of this computation and on additional visual and quantitative feedback (8). Bold frames indicate automatic processes, while light frames indicate interactive processes.

where, x is the voxel center location,  is a non-negative scalar. Eq. 1 assigns to each voxel the maximal expected risk computed with the above cost function and with respect to the input structures and risk values. For   1 , the distance is zero (e.g. voxel is located on the structure), so the voxel value is the same as the input risk value rk; it decreases as the voxel is further from the structure. Note that we chose the maximal value of the cost function over other measures, such as the sum or the average. To see why, consider the following two cases: in the first case, a voxel is located on a blood vessel that is associated with a high risk value (ri) but it is further from all the other structures; in the second case, a voxel resides near many structures associated with a low risk value indicating recoverable damage. Taking the sum or average of the cost function assigns the same risk value to both scenarios. However, the insertion of a surgical tool through the blood vessel can result in severe and unrecoverable complications (case 1) while damaging the other structures will result in less severe consequences (case 2). A given trajectory is assigned with a risk value that is the maximal value between voxels in the risk volume that intersect the given trajectory. The input is a target location, t, a set of candidate entry points {e1 , e2 , , en } , and the risk

  . The trajectory risk is:

volume riskVolume. Each target and candidate entry point pair defines a trajectory tri  ei ; t



xtr 

risk  tri , riskVolume   max riskVolume  x 

i

(2)

Blood vessel

(a) blood vessels segmentation and risk volume

(b) close-up view

Fig. 4. Cross-section image showing the blood vessels segmentation and the computed risk volume. Voxels that are near blood vessels (blue) are associated with a high risk value; the risk values decrease as the voxel is further from the blood vessel (white indicates high risk, gray indicates lower risk).

In practice, we expect the neurosurgeon to define only a handful of risk levels, i.e, risk values R  {r1 , r2 ,

, rp } are

selected from a small group of numbers ri  0,1,2,..., c . Therefore, structures with the same risk level will be assigned a risk value depending on their proximity to the closest structure (Eq. 1). We define the distance map as a volume in which each voxel is assigned with its distance to the closest structure. The distance map covers all the input structures associated with same level of risk and computed once for those structures. This allows the efficient computation of the risk volume based on a few distance maps tha cover all structures associated with various levels of risk. In our implementation, the distance map is computed once on a single volume including all segmentations of the same risk level based on a linear-time approximation method [8]. 2.3 Automatic blood vessels segmentation The trajectory planning strategy proposed in this work is independent from the actual segmentation method used to extract anatomical structures, such as blood vessels, from images. Clearly, the effectiveness of trajectory planning vitally depends on reliable segmentation methods. In the following, we will briefly describe the segmentation method used to extract blood vessels in this work. As a first step, images are segmented using an expectation-maximization (EM) algorithm [14], in which image intensity is modeled as a mixture of Gaussian Probability Density Functions (PDF). The method is initialized using Otsu's thresholding method, and voxels having maximum a posteriori probability associated to the Gaussian PDF having the highest mean are considered to be candidate vessel voxels. Voxels located outside the skull were not considered in the classification process. For typical MRA images with good signal-to-noise ratio, two Gaussian PDFs, one for vessels and one for non-vascluar tissue, are sufficient to obtain the initial segmentation. Successively, multi-scale Frangi's vesselness [15] is computed from MRA data. With this filter, Hessian eigenvalues are evaluated at each neighborhood and at different scales, and values ranging from 0 to 1 are assigned according to how close the local eigenstructure is to that of a tube. Candidate voxels obtained in the previous step corresponding to a vesselness value lower than 0.5 are discarded, so that bright structures not associated to tubelike neighborhoods at any scale are not included in the segmentation. The remaining voxels are employed to initialize a gradient-driven level-set method [10] for the optimization of the location of vessel boundaries at the ridges of gradient magnitude of image intensity.

Target #

Time (min)

1 2 3 4 5 6 7 8 Average

30 6 12 11.5 7 9 16 8 12.4

Conventional method Trajectory Closest Closest length BV NA (mm) (mm) (mm) 24 3.7 21.9 24 2.9 29.3 47 4.8 20.2 37.2 2.2 30.8 64 0.3 3.74 35.5 7.2 13.6 43 1.6 25.4 36.3 3.2 32.8 38.9 3.2 22.2

Time (min) 5 3 2 4 4 3.5 3 2 3.3

Our method Trajectory Closest length BV (mm) (mm) 67 8.6 45 4.5 45 9.2 61.2 4.36 83.9 8.96 70.4 7.2 86.2 6.56 55.6 7.04 57.4 7.0

Closest NA (mm) 24.4 19.9 20.2 21.7 3.74 24.1 25.2 32.1 21.4

Table. 1. Summary of the experimental results. Eight trajectories selected with the conventional method and with our methods are compared with respect to planning time, trajectory length, and distances of the selected trajectory from a Blood Vessel (BV) and from predefined areas for which insertion of a surgical tool is Not Allowed (NA).

3. EXPERIMENTAL RESULTS To validate our method, we conducted the following experiment on a clinical case. We compared the conventional approach for trajectory planning with our method on 8 targets selected at various locations on a clinical Magnetic Resonance Angiography (MRA) head image. The image is 256320166 voxels3 with voxel size of 0.690.690.8 mm3. For each target, the user selected two trajectories: one with the conventional method based on the axial, sagittal, and coronal 2D views of the original MRA image, and the second one with our method. The planning protocol was as in Fig 3. Eight targets were initially defined on the MRA image. Then, the outer surface of the head was automatically segmented and extracted with a method similar to that in Joskowicz et al. [3], and sampled with 40,000 points. For each target, the user defined surface areas on the outer head surface from which the entry point can be chosen. Each candidate entry point defines a candidate trajectory with the predefined target. Then, the blood vessels were automatically segmented with the method described in section 2.3, and their surfaces were reconstructed. Then, for each target, the user provided a rough indication of the no-go zones to be avoided, those for which insertion of a surgical tool is forbidden based on additional clinical knowledge. The no-go zones and the blood vessels were associated with the same risk level. The risk volume was computed using Eq. 1 with ri  1 and   0 (Fig 4). For each optional trajectory, the risk volume voxels that were intersected by the trajectory were identified, and the trajectory‟s risk was computed using Eq. 2. Next, the optimal trajectory was retrieved from the priority queue. Its risk card was automatically computed, and a visualization of the localization uncertainty and blood vessels was generated and presented to the user. The method was implemented with the Visualization ToolKit (VTK) [11] and the Insight segmentation and registration ToolKit (ITK) [12] (both by Kitware Inc., NY, USA) and integrated as a module in Slicer [13] on a standard PC computer running windows XP OS. Table 1 summarizes the results. The mean trajectory planning time using the conventional method was 12.4 min (min=6min, max=30min) compared to a mean of 3.3min (min=2min, max=5min) using our method. Using the conventional method, the mean distance of a planned trajectory to closest Blood Vessel (BV) and closest forbidden areas for which insertion of a surgical tool is Not Allowed (NA) are 3.2mm (min=0.3mm, max=7.2mm) and 22.2mm (min=3.7mm, max=32.8mm), respectively. Our method yielded mean distances of 7.0mm (min=4.36mm, max=9.2mm) from a blood vessel and 21.4mm (min=3.7mm, max=32.1mm) from no-go zones. The mean trajectory length with the conventional method was 38.9mm (min=24mm, max=64mm) compared to 57.4mm using our method. The segmentation of blood vessels, the creation of risk volume and outer head‟s surface were fully automatic, and required no user interaction. Their computation took less than 12mins.

4. DISCUSSION Our preliminary experimental study indicates that our method results in a significant shortening of the planning time and larger distances between trajectories and their closest blood vessel. The preliminary results suggest that our method produces safer trajectories for which a misplaced surgical tool is less likely to damage a critical structure. In addition, our method increases the neurosurgeon control and confidence level, and improves risk assessment. The 3D visualization of blood vessels greatly helped in understanding their complex structure and their spatial relations with respect to the planned trajectory. The uncertainty visualization enabled to check if tool misplacement may harm critical structures in a convenient manner. The risk card presented quantitative measures and was considered before the selection of the trajectory. Thus, it is easier to estimate the possible surgery outcomes and their likelihood to occur. The computed trajectories were usually longer than the manually selected ones. One reason is that the trajectory is computed to maximize the distance to closest blood vessel and no penalty was given to trajectories length. Another reason is the tendency of the user to select entry points from a nearby surface patch and ignore further areas. We conjecture that this is because it is hard to understand the complex blood vessels structure and their relation with the trajectory from the 2D views, a fortiori optimizing trajectory pose and position in this setup. What constitutes a good cost function is still an open question. Indeed, the defined guidelines describe general rules that any cost function should follow, but there are other „gray‟ areas in which the decision may be subjective. It is unclear to us, for example, if a voxel that is located on anatomy associated with a low risk value for which it is allowed but not desired to pass through the anatomy, and it is also located 8mm from another anatomy with a high risk value for which it is not allowed to pass through the anatomy, what should be the cost function? On one hand, damaging the high risk anatomy may have fatal results even if the probability for such large tool misplacement is very small. With this consideration, the risk value should be higher than the low risk value associated with the other anatomy. On the other hand, such a large misplacement is not likely to occur and raising the voxel risk value can cause disqualification of trajectories and result in a sub-optimal solution. We believe that the users should incorporate their own experience and preferences into the cost function. The presented cost function allows some level of adjustment using different risk values and α, but a more intuitive interface should be developed to ease the parameters adjustment for the specific neurosurgeon preferences. One limitation of our method, especially if additional brain structures will be incorporated, is that it requires an efficient method for the segmentation of the head image. Manual segmentation of multiple brain structures is difficult and very time-consuming, and is thus inappropriate for daily clinical use. Although some automatic and semi-automatic segmentation methods produce good results, these are usually on healthy subjects or for a specific type of disease. It remains an open question to determine if the methods can handle the large variability and irregular structures that characterize some actual patient‟s images.

5. CONCLUSIONS We have presented a novel method to enhance the conventional trajectory planning method by automatically suggesting minimal risk trajectories and by providing quantitative risk information and interactive 3D visualization of localization uncertainty and structures associated with a high risk for better assessment of the possible risks in image guided keyhole neurosurgery. Our preliminary experiment shows that trajectories selected with our method were significantly further from critical structures (e.g. blood vessels) with comparison to trajectories selected in the conventional method. While the average distance of a trajectory from a blood vessel was 3.2mm using the conventional method, it was 7.0mm using the suggested method. In addition, the suggested method is associated with a shorter planning time. Compare the average of 12.4min required for the selection of a trajectory using the conventional method and the 3.3min average time using our method. Since the experiment was performed only with one user and one head image, what are the benefits for the experienced neurosurgeon with state-of-the-art planning system is still not well studied. To address this question we are currently developing software that is better compatible with the commercial planning system. Once complete, we will test our

approach with several experts and novices neurosurgeons. We also plan to segment other critical structures and to increase our head-images database. Finally, we investigate further methods to improve risk assessment in image guided keyhole neurosurgery.

Acknowledgements This study is supported by the ROBOCAST project, EU contract FP7-ICT-215190.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

Greenberg, M. S., [Handbook of Neurosurgery], sixth edition, Thieme, New York, 545-546 (2006). Moore A. J., Newell D. W., [Neurosurgery principles and practice], Springer, Berlin, 123-138 (2005). Joskowicz, L., Shamir, R., Freiman, M., Shoham, M., Zehavi, E., Umanski F., and Shoshan, Y., "Imageguided system with miniature robot for precise positioning and targeting in keyhole neurosurgery," Comput Aided Surg, 11(4), 181-93, (2006). Tirelli, P., De Momi, E., Borghese, N. A., and Ferrigno, G., "An intelligent atlas-based planning system for keyhole neurosurgery," Computer Assisted Radiology and Surgery (CARS 2009), supplemental S85-S86, (2009). Popovic, A., and Trovato, K., "Path planning for reducing tissue damage in minimally invasive brain access" in Computer Assisted Radiology and Surgery (CARS 2009), supplemental S132-S133, (2009). Brunenberg, E. J., Vilanova, A., Visser-Vandewalle, V., Temel, Y., Ackermans, L., Platel B., and ter Haar Romeny, B. M., "Automatic trajectory planning for deep brain stimulation: a feasibility study," Int Conf Med Image Comput Comput Assist Interv (MICCAI 2007), 10, 584-92, (2007). Vaillant, M., Davatzikos, C., Taylor, R. H., Bryan, R. N., "A path-planning algorithm for image-guided neurosurgery," in First Joint Conference Computer Vision, Virtual Reality and Robotics in Medicine and Medical Robotics and Computer-Assisted Surgery, 467-476, (1997). Danielsson, P. E., "Euclidean distance mapping," Computer Vision,Graphics,and Image Processing, 14, 227– 248, (1980). Shamir, R. R., Joskowicz, L. Spektor, S. and Shoshan, Y., "Localization and registration accuracy in image guided neurosurgery: a clinical study," International Journal of Computer Assisted Radiology and Surgery, 4(1), 45-52, (2009). Piccinelli, M., Veneziani, A., Remuzzi, A., Antiga, L., "A Framework for Geometric Analysis of Vascular Structures: Application to Cerebral Aneurysms," IEEE Trans Med Imaging, 28(8), 1141-1155, (2009). The Visualization Toolkit (VTK), http://www.vtk.org/. The Insight Segmentation and Registration Toolkit (ITK), http://www.itk.org/. 3D Slicer, http://www.slicer.org/. Van Leemput, K., Maes, F., Vandermeulen, D. and Suetens, P. “Automated model-based tissue classification of MR images of the brain,” IEEE Transactions On Medical Imaging, 18(10), 885–896, (1999). Frangi, A. F., Niessen, W. J., Hoogeveen, R. M., van Walsum, T., and Viergever, M. A. “Model-based quantification of 3-D magnetic resonance angiographic images,” IEEE Transactions On Medical Imaging, 18 (10), 946–956, (1999).