Needle Insertion Simulation and Path Planning for Prostate Brachytherapy by Ehsan Dehghan Marvast B.Sc., Sharif University of Technology, 2001 M.Sc., Sharif University of Technology, 2003 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April, 2009 c© Ehsan Dehghan Marvast 2009 Abstract Low dose rate prostate brachytherapy has emerged as a treatment option for localized prostate cancer. During prostate brachytherapy, tiny radioactive capsules - seeds - are implanted inside the tissue using long needles. The quality of the treatment depends on the accuracy of seed delivery to their desired positions. Prostate deformation and displacement during insertion and lack of sufficient visual feedback complicate accurate targeting and ne- cessitates extensive training on the part of the physician. Needle insertion simulators can be useful for physician training. In addition, insertion of the needle with optimized parameters can compensate for prostate deforma- tion, can decrease the targeting errors and, subsequently, can increase the post-treatment quality of life of the patients. Therefore, needle insertion simulation and path planning have gained a lot of attention in the research community in the past decade. Moreover, several robots have been designed for brachytherapy; however, they are yet to be coupled with proper needle insertion path planning algorithms. In this thesis, steps toward a path planning algorithm for needles are taken. An optimization method is proposed that updates the initial inser- tion parameters of a rigid needle, iteratively, based on the simulated posi- tions of the targets, in order to reduce the error between the needle and several targets in a 3D tissue model. The finite element method is used in the needle insertion simulator. The simulator requires a model for the needle-tissue interaction. Therefore, an experimental method has been de- veloped to identify the force model and the tissue elasticity parameters using measurements of insertion force, tissue displacement and needle position. Ultrasound imaging is used to measure tissue displacement. Ultrasound is a common imaging modality in the operating room and does not need beads or markers to track the tissue motion. Therefore, the experimental method can be used in patient studies. The needle-tissue modeling method and insertion parameter optimization were validated in experiments with tissue-mimicking phantoms. In order to facilitate the accommodation of needle flexibility for future applications, three flexible needle models have been compared. Based on these comparisons, it was determined that an efficient and accurate angular ii Abstract spring model is best suited for future studies. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Thesis Objectives and Contributions . . . . . . . . . . . . . . 3 1.2 Needle Insertion Simulators . . . . . . . . . . . . . . . . . . . 4 1.2.1 Deformable Material Models . . . . . . . . . . . . . . 5 1.2.2 Needle-Tissue Interaction Models . . . . . . . . . . . 7 1.2.3 Flexible Needle Model . . . . . . . . . . . . . . . . . . 8 1.3 Needle Path Planning and Steering . . . . . . . . . . . . . . 9 1.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . 11 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2 Needle-Tissue Interaction Modeling using Ultrasound based Motion Estimation: Phantom Study . . . . . . . . . . . . . . 22 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . 26 2.2.1 Experimental Setup . . . . . . . . . . . . . . . . . . . 26 2.2.2 Phantom Design . . . . . . . . . . . . . . . . . . . . . 26 2.2.3 Tissue Displacement Measurement . . . . . . . . . . . 29 2.2.4 TDPE Validation . . . . . . . . . . . . . . . . . . . . 31 2.2.5 Force and Displacement Measurement . . . . . . . . . 34 2.2.6 Needle Shaft Force Distribution . . . . . . . . . . . . 38 iv Table of Contents 2.2.7 Parameter Identification Method . . . . . . . . . . . . 40 2.2.8 Simulation Method . . . . . . . . . . . . . . . . . . . 43 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.5 Conclusions and Future Work . . . . . . . . . . . . . . . . . 51 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3 Needle Insertion Study using Ultrasound-based 2D Motion Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2 Motion Estimation Algorithm . . . . . . . . . . . . . . . . . 61 3.3 Validation Experiment . . . . . . . . . . . . . . . . . . . . . 63 3.4 Needle Insertion Experiment . . . . . . . . . . . . . . . . . . 64 3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.6 Conclusions and Future Work . . . . . . . . . . . . . . . . . 68 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4 Needle Insertion Parameter Optimization for Brachyther- apy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.3 Optimization Method . . . . . . . . . . . . . . . . . . . . . . 80 4.4 Performance in Simulation . . . . . . . . . . . . . . . . . . . 84 4.5 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . 85 4.5.1 Posterior Surface Fixed . . . . . . . . . . . . . . . . . 85 4.5.2 Superior Surface Fixed . . . . . . . . . . . . . . . . . 88 4.5.3 Pubic Arch Boundary Condition Only . . . . . . . . . 88 4.5.4 Lower Young’s Modulus for the Inclusion . . . . . . . 88 4.5.5 Neo-Hookean Material Model for Tissue . . . . . . . . 91 4.6 Experimental Validation . . . . . . . . . . . . . . . . . . . . 94 4.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.8 Conclusions and Future Work . . . . . . . . . . . . . . . . . 105 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5 Modeling and Simulation of Flexible Needles . . . . . . . . 112 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.2.1 Finite Element Method using Tetrahedral Elements . 115 v Table of Contents 5.2.2 Finite Element Method using Non-Linear Beam Ele- ments . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.2.3 Angular Springs Model . . . . . . . . . . . . . . . . . 121 5.2.4 Convergence Criterion . . . . . . . . . . . . . . . . . . 125 5.3 Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.4.1 Parameter Identification . . . . . . . . . . . . . . . . 130 5.4.2 Performance in Modeling Needle Deflection . . . . . . 134 5.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 5.5.1 Approximating the Angular Spring Constants . . . . 139 5.6 Conclusions and Future Work . . . . . . . . . . . . . . . . . 139 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 6 Conclusions and Future Work . . . . . . . . . . . . . . . . . . 145 6.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 146 6.1.1 Needle-Tissue Interaction Modeling . . . . . . . . . . 146 6.1.2 Needle Insertion Initial Parameter Optimization . . . 147 6.1.3 Investigation on the Effects of Boundary Conditions and Tissue Non-linearities . . . . . . . . . . . . . . . 147 6.1.4 Flexible Needle Modeling . . . . . . . . . . . . . . . . 148 6.2 The Road Ahead . . . . . . . . . . . . . . . . . . . . . . . . . 149 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 Appendices A 3D Line Fitting to Displaced Targets . . . . . . . . . . . . . 154 B Neo-Hookean Hyperelastic Material Model . . . . . . . . . 156 C Euler-Bernoulli Beam Element . . . . . . . . . . . . . . . . . 158 D Angle of Deflection for a Bending Moment . . . . . . . . . 161 E Angle of Deflection for a Torsional (Twisting) Moment . . 163 F Code Structures for Simulations and Modeling . . . . . . . 164 F.1 Force Model Parameter Identification . . . . . . . . . . . . . 164 F.2 2D RF-Envelope Block Matching Algorithm . . . . . . . . . 168 vi Table of Contents F.3 Needle Insertion Simulation and Optimization . . . . . . . . 168 F.3.1 Needle Insertion Simulator with 3-Parameter Force Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 F.4 Simulation of Non-linear Euler-Bernoulli Beam Element . . . 177 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 vii List of Tables 2.1 Validation results for TDPE. . . . . . . . . . . . . . . . . . . 37 2.2 Needle shaft force density and elastic parameters. . . . . . . . 47 2.3 Simulation errors. . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.1 Validation results for lateral motion estimation algorithm. . . 66 3.2 Needle shaft force distribution and elastic parameters. . . . . 71 4.1 Simulation parameters for prostate model. . . . . . . . . . . . 87 4.2 Simulation results for case A. Posterior mesh surface is fixed. 87 4.3 Simulation results for case B. Superior mesh surface is fixed. . 89 4.4 Simulation results for case C. The tissue is fixed only anteri- orly near the pubic arch. E=60 kPa. . . . . . . . . . . . . . . 89 4.5 Simulation results for case D. The tissue is fixed only anteri- orly near the pubic arch. E=40kPa. . . . . . . . . . . . . . . 90 4.6 Simulation results for case E. Neo-Hookean material model. . 93 4.7 Optimization steps for insertion into tissue phantom, simu- lated and measured errors. . . . . . . . . . . . . . . . . . . . . 99 5.1 Tip deflection for various lateral tip forces. . . . . . . . . . . 129 F.1 Parameters used in the codes . . . . . . . . . . . . . . . . . . 165 F.2 User defined parameters used in F2DXCorr.mat . . . . . . . . 169 F.3 User defined parameters used in Condensed3D.m . . . . . . . . 171 F.4 Functions called inside Condensed3D.m . . . . . . . . . . . . . 173 F.5 Inputs and outputs of neoHookeanSolve3D.m . . . . . . . . . 176 F.6 User defined parameters in EGB3D.m . . . . . . . . . . . . . . . 179 F.7 Inputs and outputs of StiffnessProduceEB3D.m . . . . . . . 180 viii List of Figures 2.1 Needle insertion during prostate brachytherapy. . . . . . . . . 23 2.2 The experimental setup for needle insertion into prostate phan- tom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3 Schematic diagram of data flow in the experimental system. . 28 2.4 The tissue phantom, showing the inclusion and the hollow cylinder. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5 Cross section of the tissue phantom, showing the position of the probe and US imaging plane. . . . . . . . . . . . . . . . . 32 2.6 A portion of a normalized RF line and three overlapping win- dows. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.7 TDPE validation setup. . . . . . . . . . . . . . . . . . . . . . 35 2.8 TDPE validation results. . . . . . . . . . . . . . . . . . . . . . 36 2.9 Measured experimental data . . . . . . . . . . . . . . . . . . . 39 2.10 The needle shaft force density. . . . . . . . . . . . . . . . . . 41 2.11 Tissue mesh with tetrahedral elements. . . . . . . . . . . . . . 42 2.12 The measured force and the piece-wise linear model. . . . . . 44 2.13 The length L over which the force is integrated in the de- formed configuration. . . . . . . . . . . . . . . . . . . . . . . . 46 2.14 Simulated and measured needle forces for two different Pois- son’s ratios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 2.15 Simulated and measured node positions. . . . . . . . . . . . . 50 3.1 A sample cross-correlation matrix and the sub-sampling process. 62 3.2 The relative position of probe and its movement respect to the tissue phantom. . . . . . . . . . . . . . . . . . . . . . . . 65 3.3 The cross-section of the phantom with the trans-rectal probe and the experimental setup. . . . . . . . . . . . . . . . . . . . 67 3.4 2D displacement of the mesh nodes. . . . . . . . . . . . . . . 69 3.5 Simulated and measured insertion force and node displace- ments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.1 The prostate displacement during needle insertion. . . . . . . 77 ix List of Figures 4.2 Simplified mesh of prostate and surrounding tissue. . . . . . . 82 4.3 Target sets inside the prostate model. . . . . . . . . . . . . . 86 4.4 Optimization iterations for three targets. . . . . . . . . . . . . 92 4.5 Cross-section of the tissue phantoms and the targets. . . . . . 95 4.6 The 3D mesh of the tissue phantom used in the simulations. . 96 4.7 The 6-DoF stage, used for the positioning and orientation of the needle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.8 The target positions and the inserted needle tracks with and without optimization. . . . . . . . . . . . . . . . . . . . . . . 101 4.9 Optimization error in different iterations. . . . . . . . . . . . 104 5.1 The forces on the needle shaft caused by tissue deformation and needle base manipulation. . . . . . . . . . . . . . . . . . . 113 5.2 Euler-Bernoulli beam element under deformation. . . . . . . . 119 5.3 Angles of bending and twisting between two needle segments. 120 5.4 Three joints of our spherical wrist model in (a) initial and (b) bent configurations. . . . . . . . . . . . . . . . . . . . . . . . 123 5.5 The entire needle in equilibrium state constrained at the base and under nodal forces at its joints. . . . . . . . . . . . . . . . 124 5.6 For a sample bending simulation, (a) shows the sum of joint angles γ, which is the deflection of the needle tip, and (b) shows the residual torque at each iteration. . . . . . . . . . . 126 5.7 Needle shaft clamped at its base a vertical force is applied at its tip. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.8 Tip error and area error between the simulated model and the observed needle postures. . . . . . . . . . . . . . . . . . . 131 5.9 Needle discretized by triangular elements on a 11×2 mesh. . . 132 5.10 Distribution of identified parameters for needle models. . . . 133 5.11 The experimental deformed needle posture (circles) and the simulated needle (thicker solid lines) using the mean value of identified parameters. . . . . . . . . . . . . . . . . . . . . . . 135 5.12 Comparison of lateral tip position (left) and area (right) er- rors of different models with 10 and 20 segments using the mean value of parameters. . . . . . . . . . . . . . . . . . . . . 136 5.13 (a) Identified angular spring constants for all experiments (b) the mean spring values shown as a function of segment length in a log-log plot. . . . . . . . . . . . . . . . . . . . . . . . . . 140 D.1 A short section of a bent cantilever. . . . . . . . . . . . . . . 162 x List of Figures F.1 The measured force and the related parameters for model identification. . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 F.2 Flowchart of 2D motion estimation algorithm. . . . . . . . . . 170 F.3 Flowchart for needle insertion simulator (Condensed3D.m). . . 174 F.4 Flowchart for integration length update method. . . . . . . . 178 xi Acknowledgments I would like to thank my thesis advisor, Professor Tim Salcudean, for his support and guidance throughout the course of this research. I would also like to appreciate all my colleagues at the Robotics and Control laboratory, especially, Hani Eskandari, Reza Zahiri-Azar, Sara Khosravi, Sara Badiei, Sara Mahdavi, Ali Baghani, Julian Guerrero, Diego Prananta and Maud Marchal. Special thanks to Orcun Goksel who was always ready to help and shared with me his great knowledge of computer programming. I am thankful to my wonderful friends for the beautiful memories they created for me. I will cherish these memories for the rest of my life. I am eternally grateful to my parents and my sisters whose love and support travels around half of the globe. xii Statement of Co-Authorship This thesis is written in a manuscript-based style. Each manuscript is the result of a collaboration between several researchers. A version of Chapter 2 appeared in the Journal of Computer Aided Surgery in 2008. This paper was co-authored with Xu Wen, Reza Zahiri- Azar, Maud Marchal and Septimiu E. Salcudean. The research program for this chapter was designed based on discussions between the author, Septimiu E. Salcudean and Maud Marchal. The co- authors’ contributions to Chapter 2 are listed as follows: • The manuscript for this chapter was written by the author and revised by Septimiu E. Salcudean and Maud Marchal. • The use of ultrasound for tissue motion measurement during insertion was proposed by Septimiu E. Salcudean. • The force model and the identification algorithm used in this chapter were proposed and implemented by the author. • The numerical data analysis using the measured insertion force and displacement during the experiments were performed by the author. • The simulations for the parameter identification were carried out by the author. • The software developed by Reza Zahiri-Azar during his M.A.Sc. pro- gram was used in this chapter for tissue motion tracking. • Xu Wen and Reza Zahiri-Azar used the aforementioned software to measure tissue displacement during needle insertion experiments and model validation. • The PVC phantom used in the experiment was manufactured by Xu Wen. xiii Statement of Co-Authorship A version of Chapter 3 appeared in the International Conference on Medical Image Computing and Computer Assisted Intervention in 2008 and was co-authored with Septimiu E. Salcudean. The research program was designed based on discussions between the au- thor and Septimiu E. Salcudean. The co-authors’ contributions to Chapter 3 are summarized bellow: • The manuscript for this chapter was written by the author and revised by Septimiu E. Salcudean. • The tissue motion estimation, simulations and model parameter iden- tification were performed by the author. A version of Chapter 4 appeared in the IEEE Transactions on Robotics. This paper was co-authored with Septimiu E. Salcudean. The research program for this chapter was designed based on discussions between the author and Septimiu E. Salcudean. The co-authors contribu- tions to Chapter 4 are detailed as follows: • The manuscript for this chapter was prepared by the author and re- vised by Septimiu E. Salcudean. • The optimization algorithm was proposed and implemented by the author. • The simulations and experimental validations were performed by the author. A version of Chapter 5 has been submitted for publication. This paper was co-authored with Orcun Goksel and Septimiu E. Salcudean. The co-authors’ specific contributions in Chapter 5 are detailed below. • The experiment in Chapter 5 was conducted jointly with Orcun Gok- sel. The experiment included taking several photos of a bent needle and measuring its displacement from the photos. The displacements were measured by Orcun Goksel. • The simulations and parameter identification for the finite element based models reported in Sections 5.2.1 and 5.2.2 were carried out by the author. • The simulations and parameter identification for the angular spring model reported in Section 5.2.3 were performed by Orcun Goksel. xiv Statement of Co-Authorship • The manuscript for Chapter 5 was prepared, jointly, with Orcun Gok- sel and revised by Septimiu E. Salcudean. In detail: – Sections 5.1, 5.2.1, 5.2.2, 5.3, 5.4 and 5.4.1, and Appendix C were written by the author and revised by Orcun Goksel and Septimiu E. Salcudean. – Sections 5.2.3, 5.2.4, 5.4.2 and 5.5.1, and Appendices D and E were written by Orcun Goksel and revised by the author and Septimiu E. Salcudean. – The parts of the Sections 5.5 and 5.6 relating to the finite element based models were written by the author and the parts relating to the angular spring model were written by Orcun Goksel. • The angular spring model is a contribution of Orcun Goksel. There- fore, it is not mentioned as author’s contribution in Chapters 1 (In- troduction) or 6 (Conclusions). However, since this thesis is prepared in a manuscript based style, Chapter 5 is written exactly as the paper which has been submitted. Septimiu E. Salcudean supervised the research reported in this thesis, and assisted with editing and suggestions and improvements throughout this thesis. xv Chapter 1 Introduction Low dose-rate prostate brachytherapy – or simply brachytherapy - emerged as an effective treatment option for patients with localized prostate cancer. At least 45,000 patients were treated with brachytherapy in the US in 2004 and the numbers are increasing [1]. The procedure has been demonstrated to be “curative” in several large studies including ones carried out by our partners at the British Columbia Cancer Agency [2, 3]. This procedure entails permanent placement of radio-active capsules – seeds – of 125I or 103Pd in the prostate and peri-prostatic tissue using needles inserted through the perineum, to kill the cancer by radiation. Generally, the seed locations inside the tissue are pre-operatively planned to provide sufficient radiation on the cancerous prostate tissue but spare the bladder, rectum and urethra. In order to implant the seeds based on the plan, needles are inserted through a guiding grid with a mesh of holes with a 5mm resolution, with the aid of real-time guidance from trans-rectal ultrasound (TRUS). Non-conforming execution of the plan will result in under-radiation of the cancerous tissue or increased dose to the bladder, rectum and urethra. Under-radiation leads to treatment failure which in turn may require repeated treatment, while over-radiation of the urethra, bladder and rectum leads to toxicity after brachytherapy and causes side effects such as rectal ulceration, incontinence and painful urination [4]. Under-radiation and post-operation toxicity are still common and reduce the quality of life of the patients [5]. Generally seed misplacement errors of more that 4mm are common; however, errors of less than 3mm are clinically acceptable. A number of technical difficulties impede accurate seed delivery (con- forming to the plan) which is essential for high quality treatment with min- imal side effects. Prostate displacement, rotation and deformation, lack of sufficient 3D visual feedback, intra-operative edema (swelling caused by ex- cessive fluid in the tissue), organ registration and instrument calibration errors, needle bending and human errors are the major culprits for compli- cations due to the inaccurate execution of a plan. Needle insertion exerts forces on the soft tissue during penetration, caus- ing it to deform, displace and rotate. The rotation of the prostate during 1 Chapter 1. Introduction brachytherapy was measured by Lagerburg et al. and showed to be more than 10 degrees in both the sagittal and coronal planes [6]. Prostate rotation and deformation are major sources of error as the planned target positions move away from their initial position during the insertion of the needle. Different techniques and devices have been proposed to increase the efficacy of brachytherapy and decrease its long term side effects. Gener- ally two approaches can be taken in order to tackle this problem. The first approach is to help the physicians execute the current dosimetry plan more accurately using simulators for training [12, 13] or assisting robots [21, 22, 23, 24, 25, 26]; and the second is to update the dosimetry plan intra-operatively to compensate for the errors, should they occur [7]. In order to perform intra-operative dosimetry planning, first, the seeds should be localized after each insertion [8, 9, 10] and the radiation field should be analyzed; later, the position of the remaining seeds can be re- arranged to modify the radiation dose or new seeds can be introduced to cover the under-dosed regions. Several research groups have been working with different approaches, in order to help the physicians perform a more accurate implant. As the first step, needle insertion simulators received a lot of attention to train physicians for the complicated task of percutaneous needle insertion. Needle insertion simulators have been introduced to simulate insertion into soft tissue [11, 12, 13, 14], especially into the prostate [12, 13]. The aim of these simulators is to provide an easy and repeatable training platform that allows the physician to visualize tissue deformation during needle insertion. These needle insertion simulators were later used for needle insertion path planning [15, 16, 17, 18, 19, 20]. Needle insertion path planners, provide optimal insertion parameters and/or needle movements in order to reach a target and avoid the obstacles (bones, sensitive organs, blood vessels, etc.), in the presence of tissue deformation. It should be noted that in the current brachytherapy procedure, the seed positions are planned to be on straight lines parallel to the long axis of the TRUS probe. The guiding grid limits the orientation of the needles to lines parallel to the probe axis. The spatial insertion location is also limited to a grid of holes 5mm apart. However, the physicians often reach beyond the guiding grid and bend the needle manually to affect its path in the tissue. The optimal insertion path or insertion point will not necessarily comply to the grid constraints. Therefore, replacing the guiding grid with a guiding robot has been investigated and several brachytherapy robots have been designed to orient and localize the needle and help the physician to insert the needle or insert it automatically [21, 22, 23, 24, 25, 26]. Some of these robots are able to rotate the needle during 2 1.1. Thesis Objectives and Contributions automatic insertion in order to reduce the needle bending [24]. In addition, by providing freedom of insertion orientation, they enable the physician to reach underneath the pubic bone, which in some patients with larger prostates obstructs the needle if it is inserted parallel to the probe axis. These robots can be co-registered with the ultrasound probe in order to decrease instrument calibration errors too [22, 26]. Although these robots are mostly capable of inserting the needle with arbitrary orientation, they are yet to be coupled with a 3D prostate brachytherapy path planner. Therefore, such path planners would be extremely useful. A model based needle path planner can be based on simulations. In these path planners needle insertion simulations are used to make decisions on the insertion parameters such as initial location and orientation of the needle or the base movements during insertion. Needle insertion simulators are composed of a deformable tissue model, a flexible needle model and a needle-tissue interaction model. This thesis addresses the development of these models. 1.1 Thesis Objectives and Contributions The goal of this work is to develop a 3D simulation-based needle insertion path planner for prostate brachytherapy. In order to do so, the technical and theoretical requirements of such a planner were investigated and novel techniques were introduced to improve on current limitations. As the plan- ner is a simulation-based one, several intermediate objectives were set to be achieved for each component of the simulator. The intermediate objectives were: 1. To introduce an experimental method to investigate needle-tissue in- teraction which is also feasible for patient studies. 2. Introduce a planning algorithm for needle insertion into a 3D tissue model with multiple targets. 3. Study the effects of boundary conditions and non-linearities on the outcome of the path planner. 4. Study the effects of needle flexibility on the resulting plan. The following contributions were made in this thesis in order to achieve the objectives. 3 1.2. Needle Insertion Simulators • A new ultrasound-based experimental method was introduced to iden- tify the force distribution along the needle inside the tissue as a needle- tissue interaction model. The experimental method and the proposed modeling algorithm are able to identify the force model and tissue elas- ticity parameters using insertion force and tissue displacement data acquired during a single insertion. • The experimental method was modified to be appropriate for patient studies during prostate brachytherapy. For this purpose a 2D radio- frequency envelope block matching algorithm was used to measure 2D tissue displacement. • A new simulation-based insertion planning algorithm was introduced that optimizes the needle insertion heading, insertion point and depth for a rigid needle to be inserted into a 3D tissue model with multiple targets. • The effects of tissue non-linearities and boundary conditions on the outcome of the planning algorithm were investigated. It was shown that the accurate identification of boundary conditions is necessary in order to decide whether to use a non-linear model. • Different flexible needle models were compared in terms of their accu- racy and computational efficiency in the simulation of the flexibility of a brachytherapy needle. It was shown that the non-linear beam model is computationally more efficient than the FEM based model with triangular elements and that the latter is slightly more accurate than the former. It was also shown that the angular spring model is more accurate and more efficient compared to the other two. In the remainder of this chapter, a literature review of needle insertion simulators will be provided first. Different techniques used for each com- ponent of a needle insertion simulator will be reviewed. Later, previously developed needle path planning and steering methods will be discussed. Fi- nally, an overview of the following chapters of this thesis will be presented to clarify the logical flow of the work. 1.2 Needle Insertion Simulators In order to provide a thorough literature review of needle insertion simula- tors, previous research on each component of needle insertion simulators are 4 1.2. Needle Insertion Simulators discussed independently in the following sub-sections. The first sub-section reviews the research on deformable models used in surgery simulation or needle insertion. Next, the needle-tissue interaction models are reviewed, followed by an overview of the previous research on the flexible needle mod- els. 1.2.1 Deformable Material Models Although deformable material models have a long history in mechanical engineering and computer graphics [27], during the past few decades they have found a new purpose – visual and haptic interaction simulation for surgical simulation and planning. In this section we limit ourselves to the models used for tissue simulations. However, as a common approach among researchers, tissue mimicking materials ( plastics, gels, etc.) are used as a substitute for real human body tissue during ongoing steps of the research. Since tissue mimicking materials, especially PVC, are used in this thesis, this literature review includes the models devised for these materials as well. Mass-Spring-Damper (MSD) models that simulate the deformable ma- terial as a set of mass points connected through a lattice of Hookean springs and dampers have been used to simulate tissue deformations [16, 28, 29]. This configuration leads to a set of linear algebraic equations between nodal forces and displacements. This set of equations can be easily solved with high speed. MSD models are straightforward, easy to implement and consis- tent with the data structures used in computer graphics. They can be used for dynamic or static modeling and easily accommodate geometric changes such as cutting. However, they are not based on continuum mechanics equa- tions and therefore are prone to unrealistic deformations. Their behavior is also sensitive to the network topology (the position of the nodes). The Finite Element Method (FEM) is a powerful tool for approximately solving the complicated continuum mechanics equations [30]. In this method the global equations over the entire body Ω, are discretized over finite el- ements Ωi. The discrete problems will be approximately solved for each element and then combined to provide an approximation for the global so- lution. In contrast to MSD models, finite element models are based on continuum mechanics equations and hence are physically realistic. In addi- tion, these models are not geometry dependent as the relation between force and displacements are governed by material physical parameters such as stiffness moduli, which are independent of the position of the nodes. Finite element models are also capable of simulating the non-linear visco-elastic behavior of tissue, as it has been shown that real tissue is non-linear and 5 1.2. Needle Insertion Simulators visco-elastic [31]. However, finite element models can be computationally expensive, especially if non-linearities are included in the models. Although real tissue shows non-linear visco-elastic behavior [31], linear finite element models are used in the majority of simulators, since they are significantly faster for computation and can be used for real-time simula- tions. The condensation method, introduced by Bro-Nielsen et al. [32], paved the road for the use of linear finite element models in real-time surgery simulation. In this method, unknown variables are only computed for nodes of interest (i.e. surface nodes for graphic rendering, or nodes on the needle for force feedback), resulting in a much smaller system of linear equations that can be solved in real-time. Real-time simulation for human liver surgery was achieved by Cotin et al. using a linear elastic model in quasi-static mode at the expense of time-consuming pre-computations [33]. Linear models were specifically used in needle insertion simulators by DiMaio and Salcudean [11, 34] who achieved refresh rates suitable for haptics (≥ 500Hz), Alterovitz et al. [13] who used linear models for 2D brachytherapy simulation and Goksel et al. [12] who presented the first 3D brachytherapy simulator. Linear finite element models, based on linear geometry (linear relation between strain and the derivatives of displacements) and mechanical linear- ity (linear relation between stress and strain) are not capable of accurate simulation of tissue deformation for the large rotations and deformations that may occur during surgery or needle insertion. Therefore, models with quadratic strain measures [35, 36], hyper-elastic material models [14, 37, 38] and interpolated non-linear functions [39, 40] were used to model tissue during large deformations and/or high strains. Nienhuys and van der Stappen [14] used neo-Hookean material models in a needle insertion simulator and showed the effect of the material model on targeting accuracy. A great deal of effort has been dedicated to identification of parameters required in linear or non-linear finite element models. Hendriks et al. [38] modeled human skin as an extended Mooney material ( a generalization of neo-Hookean solids) using a numerical-experimental method in-vivo. Elastic moduli of breast and prostate ex-vivo were identified by Krouskop et al. [41]. Miller and Chinzei [42] experimentally constituted a model for brain tissue. Ophir and colleagues have introduced elastography [43, 44, 45] as a method to identify tissue elasticity parameters in-vivo. 6 1.2. Needle Insertion Simulators 1.2.2 Needle-Tissue Interaction Models Upon availability of a tissue deformation model, a needle-tissue interaction model is required to apply the effects of the needle insertion on the de- formable tissue model. Generally, the effects of insertion are applied by force or displacement boundary conditions on the nodes which are located on the needle path or very close to it. The friction force between the nee- dle and the tissue is generally modeled by a force distribution along the needle shaft. Different approaches have been taken to identify such a force distribution model. The first approach is to measure only the insertion force during the insertion of the needle. Okamura et al. measured the insertion force for a bovine liver [46]. Based on the measured force, they divided the force into three components: (1) capsule stiffness, modeled by a non-linear spring to simulate the force before capsule puncture; (2) a cutting force applied at the needle tip, a constant for each type of tissue, and (3) friction simulated by a modified Karnopp model applied along the needle shaft. The measurements took place during an ex-vivo experiment. In order to separate the friction force from the cutting force, the friction force was measured while the needle tip stuck out through the tissue, therefore no cutting occurred. Due to these constraints this method is unlikely to be used in patient studies. Misra et al. measured the insertion force of a highly flexible needle in- serted into several soft tissue mimicking phantoms and chicken tissue to investigate the sensitivity of the needle bevel force to tissue rupture tough- ness, tissue non-linearities and needle tip bevel angle [47]. A patient- and procedure-specific statistical model was devised by Pod- der et al. [48] using the experimental data recorded during brachytherapy. This statistical model predicts the maximum force applied to a brachyther- apy needle, based on some patient-specific parameters such as age, body- mass index, ethnicity, prior treatment, stage of cancer, prostate-specific anti- gen (PSA) and Gleason score, and some procedure-specific parameters such as needle diameter and average of needle insertion speed and acceleration. However, this model is limited to the prediction of maximum applied force. The second approach to needle-tissue interaction modeling is to combine the insertion force data with the tissue deformation information. While the insertion force can be simply measured using a force sensor or load cell, the tissue deformation measurement can be a challenging problem. DiMaio and Salcudean [11] used a video-camera system to track the posi- tion of several markers placed on the top surface of a slab of PVC phantom, while the phantom was perforated by a needle. Based on the force and dis- 7 1.2. Needle Insertion Simulators placement information, they introduced a force distribution model with a peak close to the needle tip followed by a constant shaft force density. In order to apply this force model in high speed simulations, they modified this model to a stick-slip interaction model [34]. In this method, as the needle tip reaches a node, the node sticks to the needle tip. In this condition the node moves with the needle tip, as long as its reaction force is smaller than a threshold. If the reaction force exceeds the threshold, the node’s state will be changed to slip state. In this state the node slips along the needle shaft, while a friction force is applied to it. Crouch et al. used a stereo-camera system and measured 3D displace- ment of a grid of markers, implanted inside a homogeneous and transparent silicon gel tissue [49]. They introduced a force distribution model with a peak at the tip followed by a dip and a constant shaft force density. Hing et al. also implanted beads inside tissue and measured their 3D displacement using a dual C-arm fluoroscope system [50]. The tissue motion measurement techniques used in the aforementioned works, relied on markers or beads [11, 49, 50] or transparency of the tissue [49]. These constraints impede the use of the above mentioned techniques for patient studies. In order to have a technique to measure tissue displacement in a patient study, especially for prostate brachytherapy, ultrasound imaging can be the best modality. Ultrasound imaging is non-ionizing and is readily available for brachytherapy. In addition, human tissue structure provides adequate scattering effects which makes tissue motion tracking possible [43, 51, 52, 53]. 1.2.3 Flexible Needle Model Although most needles are flexible, the assumption that the needle is stiff enough to be considered rigid has been widely used [11, 12, 13, 14, 34]. Depending on the mechanical properties of the needle and insertion circum- stances this assumption may be reasonable or a good approximation. Rigid needles can be used for simulation and planning when the needle bending inside tissue is negligible, for example, when the needle is stiff compared to the tissue, the needle is not manipulated during insertion, and the bevel force is negligible. In order to accommodate needle bending in models and simulation, flex- ible needles are generally divided into two groups; the first group being the group of highly flexible needles. This group includes thin and bevel-tipped needles that bend in the direction of their bevel – usually with constant curvature- without applying considerable lateral force to the tissue. Mod- 8 1.3. Needle Path Planning and Steering eled as a nonholonomic system [54], these needles were used in insertion simulation and path planning [19, 20, 55]. The members of the second group – such as brachytherapy needles – are not rigid; however, a considerable amount of lateral force is necessary to bend them. Different models have been used for these needles. DiMaio and Salcudean [15] simulated this type of needle using a finite element model with triangular elements. They took the effects of non-linear geometry into account in order to simulate large lateral deflection more accurately. This model was later extended to 3D by Goksel et al. using tetrahedral ele- ments [12]. Modeling a needle with triangular or tetrahedral elements leads to a computationally costly system of equations, especially if needle axial symmetry is preserved [12]. Modeling the needle with beam elements is computationally less demanding compared to triangular and tetrahedral el- ements, especially when large deformations are ignored. Therefore, linear beam theory has been used to model needles in simulations [56], parameter estimation [57] and bevel force effect identification [58, 59]. However, lin- ear beam models are not rotation invariant and cannot preserve the needle length during large deformations. 1.3 Needle Path Planning and Steering The needle path planning and steering algorithms can be divided into two major groups. The first group consists of the algorithms for optimization of initial parameters of insertion. These parameters include the needle heading angles, starting insertion point, insertion depth and/or bevel tip direction. In these algorithms, it is assumed that when the initial insertion parameters are set, the needle will be inserted without any manipulation. The second group manipulate the needle during insertion to guide the needle tip through tissue. Needle base movements or rotations are used to guide the needle. In the first group, the work of Alterovitz et al. on path planning for rigid [17] and highly flexible needles [19] should be noted. In [17], Alterovitz et al. introduced a sensorless path planning algorithm that optimized the insertion parameters of a rigid needle inserted inside a 2D tissue to reduce the targeting error for a single target. The insertion parameters were the insertion depth and height of the insertion point. The algorithm, which was a simulation-based one, used a search-based method to optimize the insertion parameters. In another work [19], the insertion point, insertion angle and bevel- direction were optimized for a highly flexible needle – with a nonholonomic 9 1.3. Needle Path Planning and Steering model - in 2D tissue, to reach the target and avoid obstacles. Based on sim- ulations, the insertion point and needle initial angle were optimized using a gradient descent method, with the derivatives being approximated using perturbations. To optimize for the bevel direction, the two other parameters were optimized once for each bevel direction possibility (bevel left or bevel right). The bevel direction with the better overall result was chosen as the optimal one. Manipulation of needle base during insertion was proposed by DiMaio and Salcudean [15]. They used attractive and repulsive potential fields to calculate appropriate needle tip motion inside the tissue [15]. They also introduced a needle Jacobian with which they calculated the necessary base movements that resulted in the planned tip motions. Glozman and Shoham introduced a needle steering algorithm based on the inverse kinematics of the needle [16]. They validated their algorithm using a robotic system in a closed loop setup using fluoroscope imaging as feedback [56]. Alterovitz et al. steered a highly flexible needle with a beveled tip to reach a target and avoid obstacles under Markov uncertainty using the bevel tip direction as input [20]. Kallem and Cowan developed a non-linear observer-based controller and stabilized a highly flexible bevel-tip needle in a 2D plane [60]. Reed et al. [61] combined the methods in [20] and [60] and presented a functional needle steering system that integrates the path plan- ning algorithm in [20] and the low level controller in [60] to reach a target and avoid obstacles and keep the needle in a plane during the insertion. In [20, 60, 61] a nonholonomic needle model was used. The previously developed needle path planning and steering methods assumed needles inserted into a planar environment. However, 3D deforma- tion of tissue during needle insertion in brachytherapy cannot be neglected. In addition, the real-time and closed loop control algorithms proposed in literature need real-time visual feedback which is not readily available in the operating room. Mallapragada et al. introduced the concept of tumor manipulation in- stead of needle manipulation [62]. In this method, the tumor is kept in the needle path by applying external forces to the tissue. They implemented this method on a breast phantom with a hard inclusion. This method is less likely to be adopted for brachytherapy since it requires significant changes to the current brachytherapy procedure. In addition, the prostate manip- ulation seems to be more difficult than breast due to the presence of other tissues and bones surrounding the prostate. In the rest of this thesis, we introduce and validate an algorithm to 10 1.4. Chapter Summary optimize the needle insertion parameters for a rigid needle inserted into a 3D tissue model. 1.4 Chapter Summary The goal of this thesis is to introduce a simulation-based needle insertion optimization method to be used in brachytherapy. Since the optimization algorithm is a simulation-based method, choosing a needle insertion simula- tion is the first part of the work. Since fast and efficient 2D [11] and 3D [12] needle insertion simulators were developed previously, these methods were adopted for our work. However, all these methods need a needle-tissue inter- action model, which should be applicable in later patient studies. Therefore, such a needle-tissue interaction model and its identification algorithm were set as intermediate objectives. Each chapter of this thesis presents the re- search conducted to achieve these objectives. This thesis is presented in a manuscript style, each chapter of which presents an independent research. Each chapter is accompanied by its own overview, motivation, literature review and conclusions, so it can be read independently. In Chapter 2 a new experimental method is introduced to model needle and tissue interaction and also to identify tissue elasticity parameters [63]. In this method the tissue is imaged using ultrasound during insertion of a needle while the insertion force is measured using a load cell. Ultrasound imaging does not rely on beads, markers or tissue transparency and is the most common imaging modality during prostate brachytherapy. The tissue displacement is estimated from the radio-frequency (RF) data using the time-delay cross-correlation with prior estimates (TDPE) algorithm [52]. For experiments, a two-layered PVC prostate phantom – a harder inclusion surrounded by softer tissue - was designed. The needle was inserted into the tissue phantom while the tissue was imaged using a linear US probe and the insertion force was measured. The tissue motion was measured using TDPE. A three-parameter force distribution model was introduced to model the needle-tissue interactions. The force distribution model parameters are a peak force density close to the needle tip, a shaft force density and the parameter to define the spacing between these two levels. The parameters of the force model and the elasticity parameters of the PVC phantom were identified through successive simulations in two nested loops using the force and tissue displacement data from a single insertion. In this work the needle was inserted in the axial direction of the ultra- 11 1.4. Chapter Summary sound image. The TDPE algorithm showed high resolution and accuracy in estimation of axial tissue motion. However, a trans-rectal probe is used during brachytherapy which images the tissue in the sagittal plane. There- fore, the needle insertion direction and hence the major tissue displacement are in the lateral direction of the US images. TDPE showed low resolution and accuracy in the lateral direction. In Chapter 3 an RF envelope block matching algorithm was used in order to accurately estimate 2D tissue displacement [64]. The RF envelope block matching algorithm was validated in an experiment, during which a linear stage was used to move a US probe against a piece of PVC phantom. Due to the position of the probe and PVC phantom and the direction of the movement, the PVC phantom seemed to be rigidly moving in the lateral direction of the US probe. The estimated lateral motion of the tissue was compared to the known position of the probe and showed high accuracy. In another experiment which resembled brachytherapy, the needle was inserted into the two-part PVC prostate phantom while the tissue was im- aged in the sagittal plane using a trans-rectal ultrasound (TRUS) probe. The TRUS probe was placed in a cylindrical hole in the phantom mimicking the rectum. In this situation the major tissue displacement occurred in the lateral direction of the US image. 2D tissue displacement was measured us- ing RF envelope block matching. The insertion force and tissue displacement were used to identify the force model and tissue elasticity parameters. In Chapter 4 a new needle insertion parameter optimization algorithm is introduced and experimentally validated [65]. The optimization algorithm is a simulation-based one that iteratively improves the initial insertion param- eters for a rigid needle. Assuming several targets inside a 3D tissue model, the algorithm starts with the best 3D line passing through the targets as the first insertion line. The intersection of this line and any surface in front of the tissue can be assumed as an insertion point. The simulation is per- formed with these given insertion parameters and the displaced positions of the targets are calculated. Then the best 3D line fitted to displaced tar- gets is chosen as the new insertion line. After several iterations the optimal insertion heading angles, the insertion point and depth are calculated. The optimization is performed for a simplified 3D mesh of the prostate and surrounding tissue. The effects of boundary conditions, tissue proper- ties and non-linearities are investigated through extensive simulations. It is shown that boundary conditions play an important role on the outcome of the optimization and the choice of model used – linear versus non-linear. The algorithm was used in an experiment to optimize the insertion pa- rameters to reduce the targeting error for three targets inside the prostate 12 1.4. Chapter Summary phantom, the force model for which was identified in Chapter 3. The opti- mized insertion showed significant increase in the insertion accuracy com- pared to the initial one. Although the needle is assumed to be rigid in the work presented in Chapter 4, brachytherapy needles are flexible. In Chapter 5 different mod- els are compared to model a brachytherapy needle [66]. The first model uses the finite element method with tetrahedral/triangular elements, while the second one uses Euler-Bernoulli beam elements. Nonlinear geometry is accommodated in both models to accurately simulate large deformations. During an experiment, variable lateral forces are applied to the tip of a brachytherapy needle and the consequent displacements are recorded. The model parameters are identified to fit the simulated results to the measured ones. The models are later compared in terms of their accuracy in simulating the 2D displacement of the needle that occurred during the experiment. Chapter 6 presents the collective results of each chapter. The strengths and weaknesses of the thesis and future challenges are discussed. The struc- ture of the codes used in this thesis can be found in the Appendix. 13 References [1] M. R. Cooperberg, J. M. Broering, M. S. Litwin, D. P. Lubeck, S. S. Mehta, J. M. Henning, and P. R. Carroll, “The contemporary man- agement of prostate cancer in the united states: Lessons from the can- cer of the prostate strategic urologic research endeavor (CAPSURE), a national disease registry,” Journal of Urology, vol. 171, no. 4, pp. 1393–1401, 2004. [2] W. J. Morris, M. Keyes, D. Palma, I. Spadinger, M. R. McKenzie, A. Agranovich, T. Pickles, M. Liu, W. Kwan, J. Wu, E. Berthelet, and H. Pai, “Population-based study of biochemical and survival out- comes after permanent 125I brachytherapy for low- and intermediate- risk prostate cancer,” Journal of Urology, 2008, In Press. [3] W. J. Morris, M. Keyes, D. Palma and et al., “Evaluation of dosi- metric parameters and disease response after 125Iodine transperineal brachytherapy for low- and intermediate-risk prostate cancer,” Int. J. Radiation Oncology, Biology, Physics, 2007, In Press. [4] R. Taschereau, J. Pouliot, J. Roy, and D. Tremblay, “Seed misplacement and stabilizing needles in transperineal permanent prostate implants,” Radiotherapy and Oncology, vol. 55, no. 1, pp. 59–63, 2000. [5] M. Zelefsky, Y. Yamada, G. Cohen, E. S. Venkatraman, A. Y. Fung, E. Furhang, and M. Z. D. Silvern, “Postimplantation dosimetric analy- sis of permanent transperineal prostate implantation: improved dose distributions with an intraoperative computer-optimized conformal planning technique,” Int. J. Radiation Oncology, Biology, Physics, vol. 48, no. 2, pp. 601–608, 2000. [6] V. Lagerburg, M. A. Moerland, J. J. Lagendijk, and J. J. Batter- mann, “Measurement of prostate rotation during insertion of needles for brachytherapy,” Radiotherapy and Oncology, vol. 77, no. 3, pp. 318– 323, 2005. 14 Chapter 1. References [7] S. Nag, J. P. Ciezki, R. Cormak, S. Doggett, K. Dewyngaert, G. K. Edmundson, R. G. Stock, N. N. Stone, Y. Yan, and M. J. Zelefsky, “In- traoperative planning and evaluation of permanent prostate brachyther- apy: Report of the American brachytherapy society,” Int. J. Radiation Oncology, Biology, Physics, vol. 51, no. 5, pp. 1422–1430, 2001. [8] Z. Wei, L. Gardi, D. B. Downey, and A. Fenster, “Automated lo- calization of implanted seeds in 3D TRUS images used for prostate brachytherapy,” Medical Physics, vol. 33, no. 7, pp. 24042417, 2006. [9] D. F. French, J. Morris, M. Keyes, O. Goksel, and S. E. Salcudean, “Intraoperative dosimetry for prostate brachytherapy from fused ul- trasound and fluoroscopy images,” Academic Radiology, vol. 12, pp. 1262–1272, 2005. [10] A. Jain, Y. Zhou, T. Mustufa, E. C. Burdette, G. S. Chirikjian, and G. Fichtinger, “Matching and reconstruction of brachytherapy seeds using the Hungarian algorithm (MARSHAL),”Medical Physics, vol. 32, no. 11, pp. 3475–3492, 2005. [11] S. P. DiMaio and S. E. Salcudean, “Needle insertion modeling and simulation,” IEEE Trans. Robotics and Automation: Special Issue on Medical Robotics, vol. 19, no. 5, pp. 864 – 875, 2003. [12] O. Goksel, S. E. Salcudean, and S. P. DiMaio, “3D simulation of needle- tissue interaction with application to prostate brachytherapy,” Com- puter Aided Surgery, vol. 11, no. 6, pp. 279–288, 2006. [13] R. Alterovitz, J. Pouliot, R. Taschereau, I. C. Hsu, and K. Gold- berg, “Simulating needle insertion and radioactive seed implantation for prostate brachytherapy,” in Medicine Meets Virtual Reality 11, J. Woods. et al., Eds. IOS Press, January 2003, pp. 19–25. [14] H. W. Nienhuys and A. F. van der Stappen, “A computational tech- nique for interactive needle insertions in 3D nonlinear material,” in Proc. IEEE Int. Conf. Robotics and Automation, vol. 2, 2004, pp. 2061 – 2067. [15] S. P. DiMaio and S. E. Salcudean, “Needle steering and motion planning in soft tissue,” IEEE Trans. Biomedical Engineering, vol. 52, no. 6, pp. 965 – 974, 2005. 15 Chapter 1. References [16] D. Glozman and M. Shoham, “Flexible needle steering and optimal tra- jectory planning for percutaneous therapies,” in Medical Image Com- puting and Computer Assisted Intervention, C. Barillot, D. Haynor, and P. Hellier, Eds. Springer-Verlag, 2004, pp. 137–144. [17] R. Alterovitz, K. Goldberg, J. Pouliot, R. Taschereau, and I. C. Hsu, “Sensorless planning for medical needle insertion procedures,” in Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems, vol. 3, 2003, pp. 3337–3343. [18] E. Dehghan and S. E. Salcudean, “Needle insertion point and orienta- tion optimization in non-linear tissue with application to brachyther- apy,” in Proc. IEEE Int. Conf. Robotics and Automation, 2007. [19] R. Alterovitz, K. Goldberg, and A. Okamura, “Planning for steerable bevel-tip needle insertion through 2D soft tissue with obstacles,” in Proc. IEEE Int. Conf. Robotics and Automation, 2005, pp. 1640 – 1645. [20] R. Alterovitz, A. Lim, K. Goldberg, G. S. Chirikjian, and A. M. Oka- mura, “Steering flexible needles under Markov motion uncertainty,” in Int. Conf. Intelligent Robots and Systems, 2005, pp. 120 – 125. [21] V. Lagerburg, M. A. Moerland, M. van Vulpen, and J. J. Lagendijk, “A new robotic needle insertion method to minimize attendant prostate motion,” Radiotherapy and Oncology, vol. 80, pp. 73–77, 2006. [22] G. Fichtinger, E. C. Burdette, A. Tanacs, A. Patriciu, D. Mazilu, L. L. Whitcomb, and D. Stoianovici, “Robotically assisted prostate brachytherapy with transrectal ultrasound guidance – phantom exper- iment,” Brachytherapy, vol. 5, no. 1, pp. 14–26, 2006. [23] S. P. DiMaio, S. Pieper, K. Chinzei, N. Hata, S. J. Haker, D. F. Kacher, G. Fichtinger, C. M. Tempany, and R. Kikinis, “Robot-assisted needle placement in open MRI: System architecture, integration and valida- tion,” Computer Aided Surgery, vol. 12, no. 1, pp. 15–24, 2007. [24] H. Bassan, T. Hayes, R. Patel, and M. Moallem, “A novel manipulator for 3D ultrasound guided percutaneous needle insertion,” in Proc. IEEE International Conference on Robotics and Automation, 2007, pp. 617– 622. [25] Y. D. Zhang, T. K. Podder, W. S. Ng, J. Sherman, V. Misic, D. Fuller, E. M. Messing, D. J. Rubens, J. G. Strang, R. Brasacchio, and 16 Chapter 1. References Y. Yu, “Semi-automated needling and seed delivery device for prostate brachytherapy,” in Proc. IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 2006, pp. 1279–1284. [26] S. E. Salcudean, T. D. Prananta, W. J. Morris, and I. Spadinger, “A robotic needle guide for prostate brachytherapy,” in Proc. IEEE Int. Conf. Robotics and Automation, 2008, pp. 2975 – 2981. [27] S. F. Gibson and B. Mirtich, “A survey of deformable models in com- puter graphics,” 1997, technical Report, MERL, TR-97-19. [28] F. B. de Casson and C. Laugier, “Modeling the dynamics of the human liver for minimally invasive surgery simulator,” in Medical Image Com- puting and Computer Assisted Intervention. Springer-Verlag, 1999, pp. 1156–1165. [29] D. d’Aulignac, M. C. Cavusoglu, and C. Laugier, “Modeling the dy- namics of the human thigh for realistic echographic simulator with force feedback,” in Medical Image Computing and Computer Assisted Inter- vention. Springer-Verlag, 1999, pp. 1191–1198. [30] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method: The Basis, 5th ed. Butterworht-Heinemann, 2000, vol. 1. [31] Y. C. Fung, Biomechanics: Mechanical Properties of Living Tissues, 2nd ed. New York: Springer-Verlag, 1993. [32] M. Bro-Nielsen and S. Cotin, “Real-time volumetric deformable models for surgery simulation using finite elements and condensation,” in Proc. Computer Graphics Forum, Eurographics’96, vol. 15, no. 3, 1996, pp. 57–66. [33] S. Cotin and H. Dellingate, “Real-time surgery simulation with haptic feedback using finite elements,” in Proc. IEEE Int. Conf. Robotics and Automation, vol. 4, May 1998, pp. 3739 – 3744. [34] S. P. DiMaio and S. E. Salcudean, “Interactive simulation of needle insertion models,” IEEE Trans. Biomedical Engineering, vol. 52, no. 7, pp. 1167 – 1179, 2005. [35] G. Picinbono, H. Delingette, and N. Ayache, “Nonlinear and anisotropic elastic soft tissue models for medical simulation,” in Proc. IEEE Int. Conf. Robotics and Automation, vol. 2, 2001, pp. 1370–1375. 17 Chapter 1. References [36] Y. Zhuang, “Real-time Simulation of Physically Realistic Global Defor- mations,” Ph.D. dissertation, UC Berkeley, 2000. [37] X. Wu, M. S. Downes, T. Goktekin, and F. Tendick, “Adaptive non- linear finite elements for deformable body simulation using dynamic progressive meshes,” Computer Graphics Forum, vol. 20, no. 3, pp. 349–358, 2001. [38] F. M. Hendriks, D. Brokken, J. T. van Eemeren, C. W. Oomens, F. P. Baaijens, and J. B. Horsten, “A numerical-experimental method to characterize the non-linear mechanical behavior of human skin,” Skin Research Technology, vol. 9, no. 3, pp. 274–283, 2003. [39] F. S. Azar, D. N. Metaxas, and M. D. Schnall, “A deformable finite element model of the breast for predicting mechanical deformations under external perturbations,” Journal of Academic Radiology, vol. 8, no. 10, pp. 965–975, 2001. [40] M. Mahvash and V. Hayward, “Haptic simulation of a tool in contact with a nonlinear deformable body,” in IS4TM: International Sympo- sium on Surgery Simulation and Soft Tissue Modelling, Lecture Notes in Computer Science (LNSC 2673), Springer Verlag, 2003, pp. 311– 320. [41] T. A. Krouskop, T. M. Wheeler, F. Kallel, B. S. Garria, and T. Hall, “Elastic moduli of breast and prostate tissue under compression,” Ul- trasonic Imaging, vol. 20, no. 4, pp. 260–274, 1998. [42] K. Miller, and K. Chinzei, “Constitutive modelling of brain tissue; ex- periment and theory,” Journal of Biomechanics, vol. 30, no. 11, pp. 1115–1121, 1997. [43] J. Ophir, “Elastography, a quantitative method for imaging the elastic- ity of biological tissue,” Ultrasound Imaging, vol. 13, no. 2, p. 111134, 1991. [44] J. Ophir, S. Alam, B. Garra, F. Kallel, E. Konofagou, T. Krouskop, and T. Varghese, “Elastography: ultrasonic estimation and imaging of the elastic properties of tissues,” Journal of Engineering in Medicine, vol. 213, no. 3, pp. 203–233, 1999. [45] J. Ophir, S. K. Alam, B. S. Garra, F. Kallel, E. E. Konofagou, T. Krouskop, C. R. B. Merritt, R. Righetti, R. Souchon, S. Srinivasan, 18 Chapter 1. References and T. Varghese, “Elastography: Imaging the elastic properties of soft tissues with ultrasound,” Medical Ultrasonics, vol. 29, no. 4, pp. 155– 171, 2002. [46] A. Okamura, C. Simone, and M. O’Leary, “Force modeling for needle insertion into soft tissue,” IEEE Trans. Biomedical Engineering, vol. 51, no. 10, pp. 1707 – 1716, 2004. [47] S. Misra, K. B. Reed, A. S. Douglas, K. T. Ramesh, and A. M. Okamura,“Needle-tissue interaction forces for bevel-tip steerable nee- dles,” Proc. IEEE/RAS-EMBS Int. Conf. Biomedical Robotics and Biomechatronics, pp. 224–231, 2008. [48] T. Podder, J. Sherman, E. Messing, D. Rubens, D. Fuller, J. Strang, R. Brasacchio, and Y. Yu, “Needle insertion force estimation model using procedure-specific and patient-specific criteria,” in Proc. IEEE EMBS Int. Conf., 2006, pp. 555–558. [49] J. R. Crouch, C. M. Schneider, J. Wainer, and A. M. Okamura, “A velocity-dependent model for needle insertion in soft tissue,” in Proc. Medical Image Computing and Computer Assisted Intervention (MIC- CAI), 2005, pp. 624–632. [50] J. T. Hing, A. D. Brooks, and J. P. Desai, “A biplanar fluoroscopic approach for the measurement, modeling , and simulation of needle and soft-tissue interaction,” Medical Image Analysis, vol. 11, no. 1, pp. 62–78, 2007. [51] H. M. Langevine, E. E. Konofagou, G. J. Badger, D. L. Churchill, J. R. Fox, J. Ophir, and B. S. Garra, “Tissue displacements during acupuncture using ultrasound elastography techniques,” Ultrasound in Medicine and Biology, vol. 30, no. 9, pp. 1173–1183, 2004. [52] R. Zahiri-Azar and S. E. Salcudean, “Motion estimation in ultrasound images using time domain cross correlation with prior estimates,” IEEE Trans. Biomedical Engineering, vol. 53, no. 10, pp. 1990–2000, 2006. [53] ——, “Real-time estimation of lateral motion using time domain cross correlation with prior estimates,” in Proc. IEEE Ultrasonics Sympo- sium,, 2006, pp. 1209–1212. [54] R. J. Webster, J. S. Kim, N. J. Cowan, G. S. Chirikjian, and A. M. Okamura, “Nonholonomic modeling of needle steering,” Int. J. Robotics Research , vol. 25, no. 5-6, pp. 509 – 525, 2006. 19 Chapter 1. References [55] W. Park, J. S. Kim, Y. Zhou, N. J. Cowan, A. M. Okamura, and G. S. Chirikjian, “Diffusion-based motion planning for a nonholonomic flexi- ble needle model,” in Proc. IEEE Int. Conf. Robotics and Automation, 2005, pp. 4600–4605. [56] D. Glozman and M. Shoham, “Image-guided robotic flexible needle steering,” IEEE Trans. Robotics, vol. 23, no. 3, pp. 459–466, 2007. [57] K. G. Yan, T. Podder, D. Xiao, Y. Yu, T.-I. Liu, C. Cheng, and W. S. Ng, “An improved needle steering model with online parameter estima- tor,” Int. J. Computer Assisted Radiology and Surgery, vol. 1, no. 4, pp. 205–212, 2006. [58] N. Abolhassani and R. V. Patel, “Deflection of a flexible needle during insertion into soft tissue,” in Int. Conf. IEEE Engineering in Medicine and Biology Society, 2006. [59] H. Kataoka, T. Washio, M. Audette, and K. K. Mizuhara, “A model for relation between needle deflection, force, and thickness on needle penetration,” in Medical image Computing and Computer Assisted In- tervention, LNCS 2208, W. Niessen and M. Wiergever, Eds. Springer Verlag, 2001, pp. 966–974. [60] V. Kallem and N. J. Cowan, “Image-guided control of flexible bevel-tip needles,” in Proc. IEEE Int. Conf. Robotics and Automation, 2007. [61] K. B. Reed, V. Kallem, R. Alterovitz, K. Goldberg, A. M. Okamura, and N. J. Cowan, “Integrated planning and image-guided control for planar needle steering,” Proc. IEEE/RAS-EMBS Int. Conf. Biomedical Robotics and Biomechatronics, pp. 819–824, 2008. [62] V. G. Mallapragada, N. Sarkar, and T. K. Podder, “Robot Assisted Real-time Tumor Manipulation for Breast Biopsy.” Proc. IEEE Int. Conf. Robotics and Automation, 2008, pp. 2515–2520. [63] E. Dehghan, X. Wen, R. Zahiri-Azar, M. Marchal, and S. E. Salcudean, “Needle-tissue interaction modeling using ultrasound-based motion es- timation: Phantom study,” Computer Aided Surgery, vol. 13, no. 5, pp. 265–280, 2008. [64] E. Dehghan and S. E. Salcudean, “Needle insertion study using ultrasound-based 2D motion tracking,” in Proc. Medical Image Com- puting and Computer Assisted Intervention (MICCAI), 2008, pp. 660 – 667. 20 Chapter 1. References [65] ——, “Needle insertion parameter optimization for brachytherapy,” IEEE Trans. Robotics, vol. 25, no. 2, pp. 303–315, 2009. [66] O. Goksel, E. Dehghan and S. E. Salcudean, “Modeling and simulation of flexible needles,” Medical Engineering and Physics, 2009, submitted. 21 Chapter 2 Needle-Tissue Interaction Modeling using Ultrasound based Motion Estimation: Phantom Study 1 2.1 Introduction Advances in medical robotics during the last decade have contributed to the widespread application of robots in surgical procedures [1]. The use of robotics for percutaneous needle insertion and guidance has been proposed by several groups (e.g. [2, 3, 4, 5]). There has been extensive research in the utilization of robots for prostate brachytherapy, a treatment for localized prostate cancer, during which many radioactive capsules (seeds) of 125I or 103Pd are delivered to the prostate and peri-prostatic tissue using several long needles. The seed positions are carefully planned in the pre-operative stage to maintain a sufficient level of radiation to kill cancerous cells while having a tolerable dose delivered to the urethra, rectum and bladder. During the operation, the physician inserts the needles with the aid of a guiding grid. The guiding grid allows the needle to move in parallel to the trans- rectal ultrasound (TRUS) probe axis (see Fig. 2.1). Trans-rectal ultrasound and real-time X-ray fluoroscopy are used for image guidance during the operation. The needle insertion forces deform and displace the prostate [6]. These deformations and displacements may lead to misplacement of seeds that can cause treatment errors (over-dosed or under-dosed areas) and complications (impotence or urinary incontinence [7, 8]). 1A version of this chapter has been published as: E. Dehghan, X. Wen, R. Zahiri-Azar, M. Marchal, and S. E. Salcudean, “Needle-tissue interaction modeling using ultrasound- based motion estimation: Phantom study”. Computer Aided Surgery, vol. 13, no. 5, pp. 265–280, 2008 22 2.1. Introduction Prostate Bladder PubicBone Ultrasound Probe Brachytherapy Needle Guiding Grid x y z Figure 2.1: Needle insertion during prostate brachytherapy. 23 2.1. Introduction Considering the tissue deformation during needle insertion requires con- siderable experience and skill on the part of the physician and therefore, extensive training. Needle insertion simulators can be used for physician training, path planning and robotic insertion. A number of needle insertion simulators have been developed recently [9, 10, 11], as well as path planners [12, 13, 14, 15, 16, 17] and robots to insert the needle automatically [2, 3, 18] or assist the physician to orient the needle before manual insertion [4, 5, 19, 20]. A needle insertion simulator has three major components: a deformable model to simulate tissue defor- mation, a model to simulate the needle bending if the needle is flexible, and a needle-tissue interaction model to combine needle and tissue models into one simulator. The Finite Element Method (FEM) [21] is widely used to simulate tis- sue deformation during surgery [22] and needle insertion [11, 23, 24, 25]. Different flexible needle models have also been developed and used in path planning [13, 17, 26, 27]. This chapter introduces a new experimental method to study and iden- tify a model of the tissue and the tissue-needle interaction. Needle-tissue interaction models have been studied extensively [23, 28, 29, 30, 31, 32]. Okamura et al. divided the force during insertion of a needle into a bovine liver into three components: a force due to the capsule and organ stiffness, modeled by a non-linear spring, a friction force, modeled by a modified Karnopp model, and a cutting force occurring at the needle tip, modeled by a constant force dependent on a particular tissue property [28]. The model parameters were identified using the needle insertion force, without measuring tissue displacement. The needle tip and friction forces applied by a needle during insertion into a canine prostate were reported by Kataoka et al. [29]. A patient- and procedure-specific model was devised by Podder et al. to predict the maximum force that is applied to a needle during inser- tion into the prostate and perineum [33]. A Best Linear Unbiased Estimator (BLUE) was used for this purpose. The patient-specific parameters in the model were age, body-mass index, ethnicity, prior treatment, stage of can- cer, prostate-specific antigen (PSA) and Gleason score, while the procedure- specific ones were needle diameter and average of needle insertion speed and acceleration. The insertion force data collected during brachytherapy of 25 patients were used to identify the model coefficients. Another approach to modeling the needle-tissue interaction is to com- bine the measurements of tissue deformation with measurements of needle insertion forces. Tissue deformation can be measured with different modal- ities and is used to identify the model parameters. DiMaio and Salcudean 24 2.1. Introduction used this approach to model the needle-tissue interaction during insertion of a needle into a PVC phantom [23]. They measured the tissue deformation by tracking several tissue surface markers in 2D using a video camera. The force and displacement data were incorporated through a non-linear finite element model for the tissue phantom. The elastic properties of the tissue phantom were identified prior to the experiment. They introduced a force model with a peak force density at the needle tip followed by a constant force density along the rest of the needle shaft. In order to be used in a simulator, the force density was integrated with a stick-slip model [34]. The stick-slip interaction method was later used in several other works [11, 14]. Crouch et al. extended this work to 3D [32]. They implemented a grid of fiducial beads inside a homogeneous silicon gel tissue phantom and measured their displacements in 3D during the needle insertion, using a stereo video camera system. The tissue phantom was transparent to enable the video cameras to track the implanted markers. A linear finite element analysis was carried out to identify the force model parameters from the force and displacement data. They introduced a force model composed of a peak force density near the needle tip followed by a dip and a constant shaft force density. They also studied the velocity-dependency of these parameters and the effects of tissue relaxation. Hing et al. implanted several fiducial beads inside tissue and measured their displacement during insertion using a dual C-arm fluo- roscope system [24]. They identified a local effective modulus for the tissue before puncture and an approximate cutting force for soft tissue samples. In this chapter, a new experimental method is introduced to model needle-tissue interactions. In this approach, ultrasound radio-frequency (RF) signals are used to measure the tissue displacements [35, 36]. The needle position within the tissue, the needle base force and the tissue de- formation measurements are combined to identify the force distribution pa- rameters for a two-layer phantom using a linear finite element model. The elastic moduli of the layers are identified at the same time. In contrast to other works, ultrasound imaging enables us to measure the internal tis- sue motion without implanted markers. Furthermore, ultrasound imaging is a common imaging modality in many image guided interventions includ- ing biopsies and prostate brachytherapy. Therefore, this method is easily adaptable to patient studies. Previously, Langevin et al., used RF signals to measure tissue displacements during acupuncture [37]. However, they did not derive a needle-tissue interaction model. To demonstrate the feasibility of our method, we carried out experiments with a prostate phantom. A brachytherapy needle was inserted with con- trolled speed into a non-homogeneous, non-transparent two-layer prostate 25 2.2. Materials and Methods phantom with unknown elasticity parameters. Using the recorded needle position and base force and the tissue displacement estimates using RF sig- nals, a finite element model was used to identify the force model parameters and elastic moduli jointly, from a single insertion. The experimental setup designed to insert the needle into the tissue phantom and collect data is detailed in the next section, followed by a de- scription of the construction of the phantom, an introduction to the tissue displacement measurement method and its validation. Following this, a three-parameter needle shaft force distribution model is introduced. The parameter identification algorithm and the simulation method used are de- scribed next. Experimental results and a discussion are provided, followed by conclusions and avenues for future work. 2.2 Materials and Methods 2.2.1 Experimental Setup In order to insert the needle into the phantom and measure the insertion forces and tissue displacements, a translational lead-screw stage equipped with a load cell was used to insert an 18 gauge brachytherapy needle (Bard, NJ, USA) into a prostate phantom with controlled speed. The linear stage was powered by a Maxon DC motor with an optical encoder. The linear stage and DC motor has a resolution of approximately 1µm. A proportional controller was used to control the heavily geared DC motor at a rate of 500Hz. The needle insertion force was measured via a load cell (MBD-2.5 Transducer Techniques, CA, USA) mounted at the base of the needle. B-mode images and digitized radio-frequency signals were recorded si- multaneously using a Sonix RP PC-based ultrasound machine and an L12-5 38-mm linear probe (Ultrasonix Medical Corp., BC, Canada) mounted at the back side of the tissue phantom. The experiment setup and the position of the probe can be seen in Fig. 2.2. The ultrasound machine recorded the RF data at 20 frames per second. The computer that controlled the insertion speed also transmitted the encoder and load cell data to the US machine PC through the serial port connection at 20Hz. The schematic diagram of the data flow between different devices is shown in Fig. 2.3. 2.2.2 Phantom Design A non-homogeneous phantom was constructed for the experiment. The phantom consists of an inclusion made of 100% polyvinyl chloride (PVC) 26 2.2. Materials and Methods TissuePhantomLoad Cell Ultrasound ProbeDC Motor and Encoder Needle Figure 2.2: The experimental setup, showing the tissue phantom, the needle, load cell and the linear stage. The US probe is located behind the tissue phantom. 27 2.2. Materials and Methods DC Motor 500Hz Encoder 500Hz Load Cell 20Hz 20fpsSerial connection 20Hz Figure 2.3: Schematic diagram of data flow in the experimental system. 28 2.2. Materials and Methods plasticizer (M-F Manufacturing Co., Inc. Fort Worth, TX, USA), a softer surrounding made of 66.7% PVC plasticizer and 33.3% plastic softner (M-F Manufacturing Co., Inc. Fort Worth, TX, USA) and a hollow cylinder. The phantom is shown in Fig.2.4. The harder inclusion is composed of a cylinder with two hemispheres at the two ends connected to the base by another cylinder. The inclusion is designed to mimic the prostate and its rotation around the pubic bone. The hollow cylinder is designed to simulate the rectum. During the experiment, a stiff cylinder was inserted inside this hole to mimic the effects of the trans-rectal probe on the motion of the prostate. In order to create enough scattering effect for ultrasound imaging, Cellulose (Sigma-Aldrich Inc., St Louis, MO, USA) was added to the inclusion and the surrounding. It is important to choose the phantom material carefully, as needle fric- tion forces and material brittleness vary significantly for materials that have similar Young’s moduli. PVC was chosen over agar and bovine gels, since it creates needle forces closer to the insertion forces perceptible during brachytherapy (recorded forces during brachytherapy are reported in [38]). Commonly used agar (1-3% agar added to water) and bovine gels (up to 20% gelatin added to water) create very small needle insertion forces com- pared to PVC. The needle insertion forces for agar, gel and PVC tissue phantoms with variable compositions are reported in [39]. It is shown that the insertion force for PVC is significantly higher than that of the strongest agar or gelatin phantoms. In addition, since the tissue displacement algo- rithm used in this work is based on ultrasound imaging, the quality of US images directly affects the outcome of the algorithm. PVC materials with added Cellulose, create US images with higher quality compared to the other materials. 2.2.3 Tissue Displacement Measurement The phantom was imaged to a depth of 75mm using a linear-array probe with 128 elements. The line density and sector were set to 256 and 70% to capture RF frames with 178 lines in the lateral direction and a density of 6.7 lines per millimeter. The transmit center frequency was set to 5MHz and digitized at 40MHz. The echo signals were captured at 20 frames per second. The RF signals were up-sampled with a factor of two to increase the accuracy of the estimation algorithm. The relative position of the ultrasound probe and the imaging plane with respect to the tissue phantom are shown in Fig. 2.5. The recorded data were processed off-line using a normalized cross- 29 2.2. Materials and Methods Figure 2.4: The tissue phantom, showing the inclusion and the hollow cylin- der. 30 2.2. Materials and Methods correlation based time delay estimator. To reduce the computation time, the standard algorithm was replaced with the time delay estimator with prior estimates (TDPE) [35, 36]. TDPE works similarly to the standard algorithm but it reduces the size of the search region by predicting the mo- tion using previous estimates. In order to ensure that the prediction will not result in losing track of motion, TDPE also monitors the corresponding correlation coefficients. If the estimation based on the prediction results in poor correlation coefficients, the prediction is ignored and a standard cross- correlation algorithm over a larger search region is used to recover the actual displacement. In order to process the data, each RF line was divided into 120 overlapping windows with the axial length of 1mm each (52 samples). A portion of a typical RF line and three overlapping windows are shown in Fig. 2.6. To reduce the bias and variance introduced by finite sampling in- tervals, sub-sampling using a cosine fit was used in both directions [40]. The axial and lateral components of the displacements (along the y and z axes according to Fig. 2.5) between the reference and displaced frame were esti- mated using TDPE. Initially the first frame was set as the reference frame. In order to compensate for the de-correlation caused by the large displace- ments and deformations, a dynamic reference update method was used. In this method, every time the average cross-correlation value corresponding to the latest simulated displacement dropped below a threshold (e.g. 0.95% in this work), the reference frame was replaced with the latest frame. There- fore, despite the deformation in the tissue, the deformation between each window and its corresponding window in the last reference frame is small. Using this algorithm, the absolute displacement for every spatial location was the sum of the relative displacement of the current location from the corresponding location in the last reference frame plus the displacement of the current reference frame from the initial frame. This can be shown as, d(k) = δ(k) + d0, (2.1) where, d(k) is the absolute displacement at time k, δ(k) is the relative dis- placement of current location from the last reference frame, and d0 is the displacement of the corresponding node from the initial frame to the last reference frame. 2.2.4 TDPE Validation An independent experiment was conducted to verify the accuracy of the TDPE algorithm in measuring phantom material displacement. In order to be able to measure the material displacement with another device and 31 2.2. Materials and Methods 89mm 75 mm 26.6 mm Y Z US Probe Hollow Cylinder Insertion Direction Figure 2.5: Cross section of the tissue phantom, showing the position of the probe and US imaging plane. 32 2.2. Materials and Methods 28 29 30 31 32 33−1 −0.5 0 0.5 1 Depth(mm) N or m al iz ed A m pl itu de i−1 i i+1 Figure 2.6: A portion of a normalized RF line and three overlapping win- dows. 33 2.2. Materials and Methods compare it to the TDPE results without addition of markers, a known rigid motion was applied to the material with respect to the probe. The ultra- sound probe was mounted on the lead screw stage that provides controlled motion for the probe. In addition, the tissue phantom was placed inside a box, the front side of which was closed with a thin layer of wrapping plastic. The box was filled with water, and the probe was in contact with the plastic membrane in front of the box (see Fig.2.7). In this configuration, a layer of water and a thin plastic membrane separated the probe from the tissue phantom. Therefore, the probe was able to move axially without deforming the tissue phantom, causing a relative rigid motion for the tissue phantom in the opposite direction. The probe was moved axially in sinusoidal paths with variable amplitudes and frequencies, while the tissue was imaged at 20 frames per second. The displacement of the phantom material within a region of interest was measured using TDPE. Since the motion was rigid, the displacement within this region of interest was spatially constant. The probe position was measured by the encoder and was considered as the real displacement. Figure 2.8 shows the recorded probe position and the average of the tissue displacement in the region of interest. The maximum average and standard deviation of the error are summarized in Table 2.1, showing the accuracy of the TDPE method. Please note that in one case the frequency was changed linearly in time from 0.25Hz to 2.5Hz as shown in Fig. 2.8(d). 2.2.5 Force and Displacement Measurement During the experiment, the needle was inserted using the linear stage with a controlled position along the y axis as shown in Fig. 2.9(a). In order to avoid the deteriorating effects of the metallic needle on the ultrasound images and in order to increase the accuracy of the tissue motion tracking algorithm, the needle was inserted 5mm out of the ultrasound image plane. The tis- sue phantom was meshed with 4-node tetrahedral elements to be used in a finite element analysis. The axial and lateral displacement of the mesh nodes located in the ultrasound image plane were measured with TDPE and reported in Figs. 2.9(c) and 2.9(d). The lateral motion estimates are insignificant compared to the axial motions. In addition, the resolution of the lateral motion estimation is one order of magnitude lower than the res- olution in the axial direction. Therefore, the accuracy of the lateral motion estimation can be low, especially when the lateral motion is less than a line width. As a result, in the rest of this work, only the axial motion estimation was used for modeling. The axial insertion force applied to the needle dur- 34 2.2. Materials and Methods Tissuephantom Water Soft plastic membrane US probe Figure 2.7: The schematic of the setup used for TDPE validation. Tissue phantom is in a box filled with water. The probe motion is not transmitted to the tissue phantom due to presence of water between them. 35 2.2. Materials and Methods 0 5 10 15 20 −1 −0.5 0 0.5 1 Time(s) D is pl ac em en t(m m) (a) 0 5 10 15 20 −1 −0.5 0 0.5 1 Time(s) D is pl ac em en t(m m) (b) 0 5 10 15 20 −2 −1 0 1 2 Time(s) D is pl ac em en t(m m) (c) 0 5 10 15 20 −0.5 0 0.5 Time(s) D is pl ac em en t(m m) (d) Figure 2.8: The probe motion and estimated tissue displacement for different frequencies and amplitudes. (dash-dot, blue) the probe position; (solid-red) the average of tissue displacement in the region of interest. 36 2.2. Materials and Methods Table 2.1: Validation results for TDPE. Frequency(Hz) 0.25 0.25 1.0 0.25-2.5 Amplitude(mm) 1.0 2.0 1.0 0.5 Max. Ave. Error(mm) 0.08 0.15 0.4 0.15 Max. StD. Error(mm) 5e-3 0.03 0.025 0.019 37 2.2. Materials and Methods ing the insertion is shown in Fig. 2.9(b). A needle base force equal to zero after t = 118.5 s indicates the full retraction of the needle. However, the tissue displacement measurements show an average displacement of 0.5mm after the full retraction in Fig. 2.9(c). This error may be caused by the ac- cumulation of residuals caused by the dynamic frame update method, by the slight rigid movement of tissue during the retraction of the needle, and by the minor topological changes in the tissue caused by the needle. As can be seen in Fig. 2.9(a), the needle motion has two stop periods. Relaxation effects in the tissue can be seen as a decrease in force and displacement during these two periods. Dealing with the viscoelastic properties of tissue is beyond the scope of this thesis and is the subject of further research. The needle was partially retracted and inserted after each stop. Since the needle was inserted again along the same path, tissue cutting did not occur to the same extent. Therefore, the second and third force peaks (at t = 60 s and t = 80 s) are smaller than the main force peak (t = 40 s). 2.2.6 Needle Shaft Force Distribution The needle base force during penetration (t ≤ 40 s) consists of four temporal regions which are marked in Fig.2.9(b). When the needle punctures the substrate the force increases more rapidly compared to the time that the needle penetrates in the substrate after puncture. As it can be seen in Fig. 2.9(b) the slope of force in region 1 is greater than that of region 2. The same pattern can be seen when the needle punctures and then penetrates into the inclusion as shown in regions 3 and 4 of Fig.2.9(b). Knowing that the tissue phantom has two components, these four regions can be reduced to two major components. Regions 1 and 3 belong to a tissue puncturing phase, while regions 2 and 4 occur after puncturing and can contribute to normal penetration. Inspired by the shape of the needle base force and the model introduced in [23], a three-parameter force distribution – two consecutive pulses with variable amplitude and extent – was chosen to model the needle-tissue in- teractions for each layer of the tissue. As shown in Fig. 2.10, the model is composed of a peak force density near the needle tip fp(N/m), a shaft force density fs(N/m) and the width of the the peak force density w(mm). The peak force density contributes to cutting close to the needle tip, while the shaft force density contributes to friction along the needle shaft. This model is identified to simulate the force during the main insertion phase (i.e. 0 ≤ t ≤ 40 s). The needle-tissue interactions during relaxation or retraction of the needle are not modeled. 38 2 .2 . M a teria ls a n d M eth o d s 0 20 40 60 80 100 1200 20 40 60 80 Time(s) N e e d l e P o s i t i o n ( m m ) (a) 0 20 40 60 80 100 120−4 −2 0 2 4 6 Time(s) F o r c e ( N ) 1 2 3 4 (b) 0 20 40 60 80 100 120−4 −2 0 2 4 Time (s) N o d e D i s p l a c e m e n t ( m m ) (c) 0 20 40 60 80 100 120 −0.01 0 0.01 0.02 Time(s) A v e . L a t e r a l D i s p . ( m m ) Nodes below the needle Nodes above the needle (d) Figure 2.9: (a) Needle tip position, (b) Measured insertion force, (c) Axial displacement of the mesh nodes located in the US imaging plane, and (d) Average lateral displacements of the nodes above and below the needle. In parts 1 and 3 in (b) the needle is puncturing the substrate and the inclusion, while in parts 2 and 4 the needle is penetrating into the substrate and the inclusion, respectively. The needle was partially retracted and inserted again after the main insertion. Since in the second and third insertions the needle was inserted along the same path as during the first insertion, less tissue cutting occurred. Therefore, the second and third peak forces (t =60 and 80 s) are smaller than the first one (t =40 s).39 2.2. Materials and Methods 2.2.7 Parameter Identification Method The elastic properties of the tissue phantom were not known prior to the experiment. Therefore, in addition to the force model parameters, the elas- ticity parameters of the phantom material had to be identified. A finite element model was used to identify these parameters jointly, through simu- lation. The tissue phantom was meshed for FEM simulation (see Fig. 2.11). The needle tip position, the needle insertion force and axial displacement of the mesh nodes located in the US plane were fed to the identification algorithm as experimental data. Through successive simulations the force model and elasticity parameters of the tissue phantom were identified to fit the simulated forces and displacements to the measured ones. The data from the main insertion phase were chosen for this purpose. The measured displacement drift is less than 1mm at t = 120 s. From the three sources suggested for this drift, the integration of residuals in- creases over time and the rigid movement of the tissue phantom during the retraction occurred after t = 40 s. Therefore, the drift in the displacement during 0≤ t ≤40 s is assumed to be negligible compared to the extent of tissue material displacement in this period. The unknowns to be identified for each tissue phantom region are three parameters for the force model and two elasticity parameters (Young’s mod- ulus E and Poisson’s ratio ν, assuming a linear elastic behavior for the tissue phantom material). In order to simplify the identification, the Poisson’s ra- tios were pre-assigned to decrease the number of unknowns to four for each layer of phantom. As shown in Fig. 2.9(b), the measured force in the main insertion phase should be divided into four regions. To clearly define the width and slope of the force in each of these regions a piece-wise linear model was fitted to the measured force. This model can be seen in Fig. 2.12. This model has seven parameters to be optimized. Four of these are the slopes of each region and the three remaining ones are the needle positions in which the transition from one region to the next occurs. These parameters were optimized using the Nelder-Mead method [41]. In the rest of this work, the force model parameters are optimized to fit the simulated insertion force to this piece- wise linear model. In the identification algorithm, an initial set of force model parameters and Young’s moduli are given as the inputs to a needle insertion simulator. The needle insertion simulation is achieved using the given parameters and provides a needle base force and a set of nodal displacements. In an inner loop, the force model parameters are updated to match the simulated and 40 2.2. Materials and Methods fs[N/m] fp[N/m] w y Figure 2.10: The needle shaft force density. 41 2.2. Materials and Methods Figure 2.11: Tissue mesh with tetrahedral elements. 42 2.2. Materials and Methods measured insertion forces. In order to do so, the parameter fp for each layer of the tissue phantom is adjusted so that the slopes of the simulated forces in parts 1 and 3 closely match the slopes of the corresponding recorded forces (see Fig. 2.12). The same procedure is applied to the parameter fs for each layer to match the slopes of the recorded and simulated forces in parts 2 and 4. The parameter w for each layer of the tissue phantom is also adjusted so that simulated regions 1 and 3 have similar width to the corresponding regions in the measured force. The inner loop iterations are terminated when the error between simulated and measured forces is smaller than a threshold. The parameters fp and fs are rounded to the closest integer, since higher precisions did not show a significant effect on the results. In an outer loop, the Young’s moduli are updated using the Nelder-Mead search algorithm [41]. The cost function to be minimized is: C = N∑ i=0 ∫ (dsi (l)− d m i (l)) 2W (l)dl, (2.2) where dsi and d m i are the simulated and measured axial displacements of the mesh node i in the US plane, l is the needle tip position and W (l) is a weight function that puts more weight on regions 2 and 4 since in this work, the accuracy of simulation during the normal penetration phase was more important than during the puncture phase. The outer loop iterations were terminated when the size of the final triangle in the Nelder-Mead method [41] was sufficiently small. The final values for the identified Young’s moduli were tested using perturbations to ensure they belong to local minima. 2.2.8 Simulation Method The tissue phantom was meshed using 4453 linear tetrahedral elements with 991 nodes as shown in Fig. 2.11. The linear finite element method with linear strain was used for the simulation. The tissue phantom was fixed at the back and the bottom. The hollow cylinder surface was also fixed using a hard plastic bar. Therefore, mesh nodes located on these surfaces were assumed as fixed boundary conditions. With these boundary conditions, the tissue cannot rotate or deform easily during needle insertion. It has been shown in [42] and other papers (e.g. [43]) that boundary conditions play a major rule in deciding whether non-linear models are necessary. If the tissue is confined and cannot undergo large deformations and rotations the difference between linear and non-linear models is insignificant. Under these conditions, the use of an efficient linear finite element model does not compromise the simulation accuracy [42]. 43 2.2. Materials and Methods 0 10 20 30 40 50 60 0 1 2 3 4 5 6 Needle Position(mm) Fo rc e (N ) 1 2 3 4 Figure 2.12: The measured force (dashed) and the piece-wise linear model (solid). 44 2.3. Results Due to the low speed of insertion, the simulation is performed in a quasi- static mode. The force profile shown in Fig. 2.10 was implemented in the simulation. This profile was integrated on the part of the needle which was inside the tissue phantom. This part is shown as L in Fig. 2.13. As it can be seen, this part is measured in the deformed configuration. Therefore, at each iteration, L is part of the solution. Thus, an iterative method was used to integrate the force over the equilibrium length. To distribute the force on the mesh nodes in contact with the needle the following procedure was used. The resultant force from integration of the force density from the needle tip to the first node was applied to the first node and the integration result from the first node to the second node was applied to the second node. The same procedure was used to apply forces to the other nodes in contact with the needle. In this simulation approach the forces are applied only to the nodes. However, a needle may enter an element from any point on the surface of the element. Therefore, mesh adaptation [11] was used to increase the accuracy of the simulation. With mesh adaptation, every time the needle tip reaches the surface of an element, the closest node on that surface is re-located on the needle tip in the undeformed mesh. The stiffness matrix is updated using the Woodbury matrix identity to accommodate the topological change in the mesh with low computational costs [11]. No external force was applied to re-locate the node, since the modifications are applied to the undeformed mesh. A force boundary condition is applied to every node in contact with the needle, in the axial direction of the needle. Two displacement boundary conditions are applied to the same nodes in the other two orthogonal directions to maintain the nodes on the needle. 2.3 Results The identification method takes the measured needle force and nodal axial displacements as input data. In addition, the Poisson’s ratio has to be set as a known parameter. Two values of Poisson’s ratios were considered - 0.45 and 0.49 - for both the inclusion and the substrate. In the second case, with Poisson’s ratio of 0.49, the tissue is nearly incompressible, which is one of the characteristics of real tissue [44]. The identified force model parameters and Young’s moduli for each case are presented in Table 2.2. Figure 2.14 shows the measured and simulated needle base forces for both cases. Table 2.3 presents the maximum mean and standard deviation of the errors in the prediction of the displacements as functions of time. Figure 2.15 illustrates the simulated and measured position of the nodes 45 2.3. Results L Figure 2.13: The length L over which the force is integrated in the deformed configuration. The dotted lines show the elements in the undeformed con- figuration. 46 2.3. Results Table 2.2: Needle shaft force density and elastic parameters. ν Tissue fs(N/m) fp(N/m) w(mm) E(kPa) 0.45 Inclusion Surrounding 70 61 200 168 12.5 2.5 41.0 8.5 0.49 Inclusion Surrounding 69 63 206 117 12.0 3.0 24.0 7.0 Table 2.3: Simulation errors. ν Max. Mean Disp. Error(mm) Max. Std Error(mm) 0.45 0.5 0.5 0.49 0.4 0.6 47 2.4. Discussion located in the US plane in the deformed and undeformed configurations for ν = 0.49. 2.4 Discussion In this work, the force model parameters and the Young’s moduli were identified jointly for a two-layer phantom using the experimental data from a single needle insertion. In the presence of several layers of tissue, the joint identification can be very complex due to the mutual effects of different layers of tissue on each other. However, for specific procedures such as brachytherapy, it will be useful to identify both the tissue model and the interaction model at the same time. Furthermore, the methodology introduced in this chapter can be used to identify the force model parameters upon the availability of elasticity prop- erties from other sources. The inner loop of the identification algorithm can deliver the force model parameters for many layers of tissue, if the Young’s moduli are known. Therefore, it is useful to identify tissue elasticity param- eters independently, if possible. Research in the field of vibro-elastography is in progress to provide the elasticity parameters of the prostate and sur- rounding tissue [45]. The Young’s modulus of the inclusion material was independently mea- sured using rheometry and was between 30 and 36 kPa. The Young’s mod- ulus varies based on the heating and cooling procedure during PVC fab- rication. It also changes with time after the fabrication. Note that the identified Young’s moduli are not close to the real Young’s moduli of the tissue. As shown in this work, the identified Young’s moduli depend on the assumed Poisson’s ratio. The mesh quality affects the accuracy of the identified parameters too. Large tetrahedral elements show over-stiffness since they impose their constant strain field over a large volume of tissue [46]. A mesh with higher resolution can yield identified elasticity parame- ters with higher accuracy at the expense of higher computational burden. The minimization cost function can also change the identified parameters. As a comparison, the elasticity parameters in [47] were adjusted to fit the average of the simulated axial displacements to the measured ones. The identified parameters in the mentioned work are different from the results provided in this chapter. However, it should be noted that the ultimate goal of this work is to achieve accurate insertion force and tissue displacement prediction. Therefore, mechanically inaccurate values for Young’s moduli are acceptable if they lead to accurate simulation of tissue deformations. 48 2.4. Discussion 0 10 20 30 40 50 60 0 2 4 6 0 10 20 30 40 50 60 0 2 4 6 Needle Position(mm) N ee dl e Fo rc e(N ) Sim. Force ν=0.45 Measured Force Sim. Force ν=0.49 Measured Force Figure 2.14: Simulated and measured needle forces for two different Pois- son’s ratios. 49 2.4. Discussion 20 40 60 80 100 25 35 45 55 Y(mm) Z( mm ) ν = 0.49 Figure 2.15: Position of the nodes in the US field of view. Squares denote the original positions, stars simulated positions and circles the positions measured with TDPE. The dotted line shows the projection of the needle on the US plane. Only axial displacements are considered. The Poisson’s ratio is equal to 0.49. 50 2.5. Conclusions and Future Work In the identification algorithm, the cost function in (2.2) was used as a similarity measure. However, other similarity measurements can be used with the same method. The similarity measure in [47] was the average of nodal displacements. Since positive and negative errors can cancel each other in the average of nodal displacements, the simulated displacements in [47] have higher standard deviations of the errors. A linear FEM model was used to increase the computational speed. However, it is possible to substitute this linear model with a non-linear one without increasing the number of unknowns. As an example, the Lamé parameters (interchangeable with Young’s modulus and Poisson’s ratio) are the only parameters necessary to build the constitutive equations for neo- Hookean material models [48]. Therefore such a model can replace the linear model if higher accuracy is necessary and computational time can be increased. The tissue displacement estimation method was validated when the tis- sue motion was rigid. However, tissue deforms during needle insertion. The reference frame update algorithm was used to compensate for this deforma- tion. The TDPE algorithm can be validated in the presence of deformation, by tracking the position of several implanted beads and comparing the re- sults with the results of other tracking methods. 2.5 Conclusions and Future Work In this work, a new methodology for the investigation of needle-tissue inter- action was introduced. In this method, ultrasound radio-frequency data are recorded during the insertion and are used for tissue displacement measure- ment. The tissue displacement data, together with the needle base position and needle insertion force data, are used to identify the parameters of a three-parameter force distribution along the needle. In this chapter, these parameters were identified for a two-layer non-homogeneous PVC prostate phantom. The tissue phantom was composed of a harder inclusion mimick- ing the prostate surrounded by a softer tissue. The tissue phantom was not transparent and had no fiducial beads or markers. Prior knowledge of the elastic properties of the phantom were not used to determine the needle- tissue interaction model. The Young’s moduli of the tissue phantom were identified jointly with the force model parameters with the aid of a linear finite element simulation. The tissue phantom was meshed for this purpose and the displacements of the mesh nodes located in the US plane were mea- sured using a cross-correlation method. In an inner loop the force model 51 2.5. Conclusions and Future Work parameters were adjusted to match the simulated and the measured forces, and in the outer loop the Young’s moduli were adjusted to minimize the cost function in (2.2). Two different cases with two different Poisson’s ratios were assumed. The identified force parameters and elastic properties can be used to build a needle insertion simulator for physician training, path planning and robotic insertion. In this work a linear probe was used to acquire the images during in- sertion. However, during the brachytherapy procedure, a trans-rectal probe is used to image the prostate in sagittal/parasagittal and transverse planes. Therefore, during brachytherapy, the major displacement of tissue can be seen at best in the lateral direction of the US image, which has significantly lower resolution (approximately one order of magnitude) than the axial mo- tion estimation. The effect of this discrepancy of resolution on the model depends on the tissue properties and will be the subject of further research. The proposed algorithm will be used with finer meshes to increase the accuracy of the identification. The inner loop of the suggested method should be combined with other tissue identification algorithms (e.g. vibro- elastography) to identify the force model parameters for real tissue with several layers. The visco-elastic properties of the tissue and their effect on the identified parameters will be investigated in the future. Effects of needle insertion speed on the force model will be studied. The tissue and force model parameters were identified using a single insertion. Further validation to show the efficacy of the model to simulate tissue deformation and insertion force for other insertion trajectories is part of the future work. 52 References [1] R. H. Taylor and D. Stoianovici, “Medical robotics in computer- integrated surgery,” IEEE Trans. Robotics and Automation, vol. 19, no. 5, pp. 765–781, 2003. [2] H. Bassan, T. Hayes, R. Patel, and M. Moallem, “A novel manipulator for 3D ultrasound guided percutaneous needle insertion,” in Proc. IEEE Int. Conf. Robotics and Automation, 2007, pp. 617–622. [3] V. Lagerburg, M. A. Moerland, M. van Vulpen, and J. J. Lagendijk, “A new robotic needle insertion method to minimise attendant prostate motion,” Radiotherapy and Oncology, vol. 80, pp. 73–77, 2006. [4] G. Fichtinger, J. Fiene, C. W. Kennedy, G. Kronreif, I. Iordachita, D. Y. Song, E. C. Burdette, and P. Kazanzides, “Robotic assistance for ultrasound guided prostate brachytherapy,” in Proc. Medical Image Computation and Computer Assisted Intervention (MICCAI), N. Ay- ache, S. Ourselin, and A. Maeder, Eds., 2007, pp. 119–127. [5] S. E. Salcudean, T. D. Prananta, W. J. Morris, and I. Spadinger, “A robotic needle guide for prostate brachytherapy,” in Proc. IEEE Int. Conf. Robotics and Automation (ICRA), 2008, 2975 – 2981. [6] V. Lagerburg, M. A. Moerland, J. J. Lagendijk, and J. J. Batter- mann, “Measurement of prostate rotation during insertion of needles for brachytherapy,” Radiotherapy and Oncology, vol. 77, no. 3, pp. 318– 323, 2005. [7] J. E. Dawson, T. Wu, T. Roy, J. Y. Gy, and H. Kim, “Dose effects of seeds placement deviation from pre-planned position in ultrasound guided prostate implants,” Radiotherapy and Oncology, vol. 32, no. 2, pp. 268–270, 1994. [8] R. Taschereau, J. Pouliot, J. Roy, and D. Tremblay, “Seed misplacement and stabilizing needles in transperineal permanent prostate implants,” Radiotherapy and Oncology, vol. 55, no. 1, pp. 59–63, 2000. 53 Chapter 2. References [9] R. Alterovitz, J. Pouliot, R. Taschereau, I. C. Hsu, and K. Goldberg, “Needle insertion and radioactive seed implantation in human tissue: Simulation and sensitivity analysis,” in Proc. IEEE Int. Conf. Robotics and Automation (ICRA), vol. 2, 2003, pp. 1793–1799. [10] ——, “Simulating needle insertion and radioactive seed implantation for prostate brachytherapy,” in Medicine Meets Virtual Reality 11, West- wood J.D. et al., Ed. IOS Press, January 2003, pp. 19–25. [11] O. Goksel, S. E. Salcudean, and S. P. DiMaio, “3D simulation of needle- tissue interaction with application to prostate brachytherapy,” Com- puter Aided Surgery, vol. 11, no. 6, pp. 279–288, 2006. [12] R. Alterovitz, K. Goldberg, J. Pouliot, R. Taschereau, and I. C. Hsu, “Sensorless planning for medical needle insertion procedures,” in Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems, vol. 3, 2003, pp. 3337–3343. [13] S. P. DiMaio and S. E. Salcudean, “Needle steering and motion planning in soft tissue,” IEEE Trans. Biomedical Engineering, vol. 52, no. 6, pp. 965 – 974, 2005. [14] E. Dehghan and S. E. Salcudean, “Needle insertion point and orienta- tion optimization in non-linear tissue with application to brachyther- apy,” in Proc. IEEE Int. Conf. Robotics and Automation, 2007, pp. 2267–2272. [15] K. G. Yan, T. Podder, D. Xiao, Y. Yu, T.-I. Liu, C. Cheng, and W. S. Ng, “An improved needle steering model with online parameter estima- tor,” Int. J. Computer Assisted Radiology and Surgery, vol. 1, no. 4, pp. 205–212, 2006. [16] R. Alterovitz, A. Lim, K. Goldberg, G. S. Chirikjian, and A. M. Oka- mura, “Steering flexible needles under Markov motion uncertainty,” in Int. Conf. Intelligent Robots and Systems, 2005, pp. 1570 –1575. [17] D. Glozman and M. Shoham, “Flexible needle steering and optimal trajectory planning for percutaneous therapies,” in Proc: Medical Image Computing and Computer Assisted Intervention, C. Barillot, D. Haynor, and P. Hellier, Eds. Springer-Verlag, 2004, pp. 137–144. [18] Y. D. Zhang, T. K. Podder, W. S. Ng, J. Sherman, V. Misic, D. Fuller, E. M. Messing, D. J. Rubens, J. G. Strang, R. Brasacchio, and 54 Chapter 2. References Y. Yu, “Semi-automated needling and seed delivery device for prostate brachytherapy,” in Proc. IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 2006, pp. 1279–1284. [19] G. Fichtinger, E. C. Burdette, A. Tanacs, A. Patriciu, D. Mazilu, L. L. Whitcomb, and D. Stoianovici, “Robotically assisted prostate brachytherapy with transrectal ultrasound guidance – phantom exper- iment,” Brachytherapy, vol. 5, no. 1, pp. 14–26, 2006. [20] J. Hong, T. Dohi, M. Hashizume, K. Konishi, and N. Hata, “An ultrasound-driven needle-insertion robot for percutaneous cholecys- tostomy,” Physics in Medicine and Biology, vol. 49, pp. 441–455, 2004. [21] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method: The Basis, 5th ed. Butterworht-Heinemann, 2000, vol. 1. [22] M. Bro-Nielsen and S. Cotin, “Real-time volumetric deformable models for surgery simulation using finite elements and condensation,” in Proc. Computer Graphics Forum, Eurographics’96, vol. 15, no. 3, 1996, pp. 57–66. [23] S. P. DiMaio and S. E. Salcudean, “Needle insertion modeling and simulation,” IEEE Trans. Robotics and Automation: Special Issue on Medical Robotics, vol. 19, no. 5, pp. 864 – 875, 2003. [24] J. T. Hing, A. D. Brooks, and J. P. Desai, “A biplanar fluoroscopic approach for the measurement, modeling , and simulation of needle and soft-tissue interaction,” Medical Image Analysis, vol. 11, no. 1, pp. 62–78, 2007. [25] H. W. Nienhuys and A. F. van der Stappen, “A computational tech- nique for interactive needle insertions in 3D nonlinear material,” in Proc. IEEE Int. Conf. Robotics and Automation, vol. 2, 2004, pp. 2061 – 2067. [26] E. Dehghan, O. Goksel, and S. E. Salcudean, “A comparison of needle bending models,” in Proc. Medical Image Computing and Computer As- sisted Intervention (MICCAI), R. Larsen, M. Nielsen, and J. Sporring, Eds., 2006, pp. 305–312. [27] R. J. Webster, J. S. Kim, N. J. Cowan, G. S. Chirikjian, and A. M. Okamura, “Nonholonomic modeling of needle steering,” Int. J. Robotics Research, vol. 25, no. 5-6, pp. 509 – 525, 2006. 55 Chapter 2. References [28] A. Okamura, C. Simone, and M. O’Leary, “Force modeling for needle insertion into soft tissue,” IEEE Trans. Biomedical Engineering, vol. 51, no. 10, pp. 1707 – 1716, 2004. [29] H. Kataoka, T. Washio, K. Chinzei, K. Mizuhara, C. Simone, and A. Okamura, “Measurement of tip and friction force acting on a needle during penetration,” in Proc. Medical Image Computing and Computer Assisted Intervention - MICCAI, T. Dohi and R. Kikinis, Eds., 2002, pp. 216–223. [30] T. Podder, J. Sherman, E. Messing, D. Rubens, D. Fuller, J. Strang, R. Brasacchio, and Y. Yu, “Needle insertion force estimation model us- ing procedure-specific and patient-specific criteria,” in Proc. Int. Conf. IEEE Engineering in Medicine and Biology Society, 2006, pp. 555–558. [31] J. T. Hing, A. D. Brooks, and J. P. Desai, “Reality-based needle inser- tion simulation for haptic feedback in prostate brachytherapy,” in Proc. IEEE Int. Conf. Robotics and Automation, 2006, pp. 619–624. [32] J. R. Crouch, C. M. Schneider, J. Wainer, and A. M. Okamura, “A velocity-dependent model for needle insertion in soft tissue,” in Proc. Medical Image Computing and Computer Assited Intervention (MIC- CAI), J. Duncan and G. Gerig, Eds., 2005, pp. 624–632. [33] T. Podder, J. Sherman, E. Messing, D. Rubens, D. Fuller, J. Strang, R. Brasacchio, and Y. Yu, “Needle insertion force estimation model using procedure-specific and patient-specific criteria,” in Proc. IEEE EMBS Int. Conf., 2006, pp. 555–558. [34] S. P. DiMaio and S. E. Salcudean, “Interactive simulation of needle insertion models,” IEEE Trans. Biomedical Engineering, vol. 52, no. 7, pp. 1167 – 1179, 2005. [35] R. Zahiri-Azar and S. E. Salcudean, “Motion estimation in ultrasound images using time domain cross correlation with prior estimates,” IEEE Trans. Biomedical Engineering, vol. 53, no. 10, pp. 1990–2000, 2006. [36] ——, “Real-time estimation of lateral motion using time domain cross correlation with prior estimates,” in Proc. IEEE Ultrasonics Sympo- sium,, 2006, pp. 1209–1212. [37] H. M. Langevine, E. E. Konofagou, G. J. Badger, D. L. Churchill, J. R. Fox, J. Ophir, and B. S. Garra, “Tissue displacements during 56 Chapter 2. References acupuncture using ultrasound elastography techniques,” Ultrasound in Medicine and Biology, vol. 30, no. 9, pp. 1173–1183, 2004. [38] T. Podder, D. Clark, J. Sherman, D. Fuller, E. Messing, D. Rubens, J. Strang, R. Brasacchio, L. Liao, W.-S. Ng, and Y. Yu, “In vivo motion and force measurement of surgical needle intervention during prostate brachytherapy,” Medical Physics, vol. 33, no. 8, pp. 2915–2922, 2006. [39] T. Podder, D. Clark, J. Sherman, D. Fuller, D. Rubins, W. Ng, E. Mess- ing, W. O’Dell, J. Strang, Y. Zhang, and Y. Yu, “Robotic needle in- sertion in soft material phantoms: An evaluation of properties of com- monly used soft materials,” in IEEE Int. Conf. Biomedical Engineering (ICBME), 2005. [40] I. Cespedes, Y. Huang, J. Ophir, and S. Spratt, “Methods for the esti- mation of subsample time-delays of digitized echo signals,” Ultrasonic Imaging, vol. 17, no. 2, pp. 142–171, 1995. [41] M. Avriel, Nonlinear Programming: Analysis and Methods. Dover Publications, 2003. [42] E. Dehghan and S. E. Salcudean, “Comparison of linear and non-linear models in 2D needle insertion simulation,” in Proc. Workshop on Com- putational Biomechanics for Medicine, 2006, in conjunction with MIC- CAI 2006. [43] M. Chabanas, Y. Payan, C. Marecaux, P. Swider, and F. Boutault, “Comparison of linear and non-linear soft tissue models with post- operative ct scan in maxillofacial surgery,” Lecture Notes in Computer Science, vol. 3078, pp. 19–27, 2004. [44] J. Ophir, S. K. Alam, B. S. Garra, F. Kallel, E. E. Konofagou, T. Krouskop, C. R. B. Merritt, R. Righetti, R. Souchon, S. Srinivasan, and T. Varghese, “Elastography: Imaging the elastic properties of soft tissues with ultrasound,” Medical Ultrasonics, vol. 29, no. 4, pp. 155– 171, 2002. [45] S. E. Salcudean, D. French, S. Bachmann, R. Zahiri-Azar, X. Wen, and W. J. Morris, “Viscoelasticity modeling of the prostate region us- ing vibro-elastography,” in Proc. Medical Image Computing and Com- puter Assisted Intervention (MICCAI), R. Larsen, M. Nielsen, and J. Sporring, Eds., 2006, pp. 389–396. 57 Chapter 2. References [46] R. D. Cook, D. S. Malcus, and M. E. Plesha, Concepts and Applications of Finite Element Analysis. John Wiley & Sons, 1989. [47] E. Dehghan, X. Wen, R. Zahiri-Azar, M. Marchal, and S. E. Salcudean, “Modeling of needle-tissue interaction using ultrasound-based motion estimation,” in Proc: Medical Image Computing and Computer Assisted Intervention (MICCAI), 2007. [48] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method: Solid Mechanics, 5th ed. Butterworht-Heinemann, 2000, vol. 2. 58 Chapter 3 Needle Insertion Study using Ultrasound-based 2D Motion Tracking 2 3.1 Introduction Prostate brachytherapy is a common treatment for prostate cancer during which radioactive capsules are implanted in the prostate and surrounding tissue using needles. Accurate placement of the capsules leads to a proper dose and increases the quality of the treatment [1]. However, due to prostate deformation and rotation [2] and lack of visual feedback, accurate needle placement is still a challenge. Brachytherapy simulators have the potential to train physicians and provide pre-operative planning. Every needle inser- tion simulator requires a needle-tissue interaction model. Extensive research has been conducted to provide such a model [3, 4, 5], including studies on the liver [6, 7] and in-vivo studies on the prostate [8]. In some model- ing approaches, only the needle force measurements were used to derive the needle-tissue interaction parameters [6, 8], while in others, tissue dis- placements were measured and employed in addition to force measurements [3, 4, 5, 7]. In [3], 2D displacements of a PVC phantom were measured using surface markers tracked by a camera. By combining the force and displace- ment data, a needle force model was introduced that combines a peak at the tip of the needle followed by a constant force distribution along the rest of the needle. Crouch et al. measured the 3D displacement of several markers implanted inside a cube of homogeneous, transparent silicon gel phantom using a stereo-camera system [4]. The force and displacement data were used to identify the parameters of a force distribution model with a peak and a dip at the tip of the needle followed by a constant shaft force density. 2A version of this chapter has been published as: E. Dehghan, and S. E. Salcudean, “Needle Insertion Study using Ultrasound-based 2D Motion Tracking”. In: D. Metaxas et al. (Eds.), MICCAI 2008, Part II, LNCS 5242, pp. 660-667, 2008 59 3.1. Introduction Hing et al. [7] used a dual C-arm fluoroscope to measure the displacement of fiducial beads inside the tissue in 3D. They identified a local effective mod- ulus during puncture and a cutting force for several samples of tissue. The tissue displacement measurement techniques used by the aforementioned re- search groups rely on the implanted markers, beads and/or transparency of the tissue. In the previous chapter, we used ultrasound (US) radio-frequency (RF) signals to estimate the tissue displacement in the US field of view. Ultra- sound imaging does not need any beads or markers and is a common imaging modality in many medical procedures including prostate brachytherapy. An experiment was carried out in which a brachytherapy needle was inserted into a PVC prostate phantom. The insertion forces were measured via a load cell and the RF signals were collected using a linear probe mounted at the back of the tissue. Correlation-based methods [9, 10] were used to measure the tissue displacement from the (RF) data. Since the major displacement of the tissue occurred parallel to the axial direction of the US images, the axial tissue displacement was used to identify the unknown parameters of the model. The parameters of a three-parameter force distribution model - tip and shaft force densities and the spacing between them - and the Young’s moduli of the two parts of the tissue phantom were identified in a single in- sertion. During prostate brachytherapy, a trans-rectal ultrasound (TRUS) probe is used to image the tissue in transverse or sagittal/parasagittal planes. If the tissue is imaged in sagittal/parasagittal planes, the major tissue dis- placement (along the needle insertion axis) occurs in the lateral direction of the US image (lateral and axial directions of the US image and the axis of insertion are shown in Fig. 3.3(a)). Therefore, to study needle-tissue inter- action during brachytherapy, one should be able to estimate tissue displace- ment in the lateral direction with high accuracy. The decoupled correlation methods employed in the previous chapter did not have high accuracy in the lateral direction due to image de-correlation over large displacements. In this chapter, an RF envelope block matching technique is used to estimate tissue deformation in 2D. Two dimensional motion estimation us- ing speckle tracking methods performed on B-Scan (RF envelope) images was reported in [11, 12, 13]. RF envelope cross-correlation outperforms the RF data cross-correlation method for large lateral displacements and high strains [14]. An experiment was carried out in which the TRUS probe was used to image a non-homogeneous PVC tissue phantom during insertion of a needle. This experiment mimics the practical case in which a needle inserted during prostate brachytherapy is imaged with the sagittal (axial) array of the TRUS probe. The most significant tissue displacements, which 60 3.2. Motion Estimation Algorithm occur along the needle insertion axis in the lateral direction of the sagittal ultrasound array, the needle position and needle insertion forces are used to model the needle-tissue interaction. In particular, the gathered data is used to identify the three parameters of the force model introduced in the previous chapter and the Young’s moduli of the two-layered tissue phantom. This tissue displacement measurement method can be used in-vivo, during prostate brachytherapy, to measure prostate motion and rotation during in- sertion, which in turn can lead to prostate deformation simulators and path planners for intra-operative treatment planners. 3.2 Motion Estimation Algorithm The motion estimation algorithm uses the envelope of the RF data. The RF envelope is calculated using: ei = √ |h(RFi)|, (3.1) where RFi is the i th RF data line, h(·) denotes the Hilbert transform and ei is the i th line of the RF envelope. The algorithm tracks the position of the points of interest in several consecutive frames. At first, a Block of Interest (BI) is assumed around every point of interest in the reference frame. Later, the algorithm locates the BIs in the subsequent frames using 2D cross-correlation. For this purpose, each BI is cross-correlated with a Search Region (SR) in the subsequent frames. The position of the BI in each frame and its displacement are found using the maximum values of the cross-correlation matrix. Cross-correlation of the BIs with the entire frame is time consuming and results in multiple false peaks. Therefore, each BI is cross-correlated with a corresponding SR, the size of which is larger than the BI but smaller than the entire frame. The positions of the SRs in the target frame are calculated based on the estimated positions of the corresponding BIs up to the target frame. To find the lateral displacement with higher resolution sub-sampling using quadratic functions was employed [15]. In this method, at first, the maximum value of the cross-correlation matrix is chosen. Later, a quadratic function is fitted to the corresponding row of the matrix. The maximum value of this quadratic function is assumed to correspond to the true displacement. The cross-correlation matrix and sub-sampling processes are shown in Fig.3.1. 61 3.2. Motion Estimation Algorithm 1 2 3 5 10 15 20 0.9 0.92 0.94 0.96 0.98 1 Lateral PositionAxial Position Cr os s− Co rre la tio n Re su lts (a) −1 0 1 0.95 0.975 1 Lateral Displacement (Line) Cr os s− Co rre la tio n Re su lts Estimated Displacement (b) Figure 3.1: (a) a sample cross-correlation matrix and (b) the sub-sampling process. 62 3.3. Validation Experiment 3.3 Validation Experiment In order to verify the accuracy of the block matching algorithm, a rigid known movement was applied to the tissue phantom and the motion was measured using block matching. A linear probe was mounted on a linear stage powered by a DC motor. A PVC phantom was placed inside a box, the front side of which was covered by soft wrapping plastic. The box was filled with water. The probe was placed in touch with the plastic membrane (see Fig. 3.2). In this situation, the probe can move without applying any force on the PVC phantom. The probe was moved sinusoidally (x(t) = A sin(ωt)) with variable amplitudes and frequencies. The direction of the probe movement was parallel to the lateral direction of the US image. Due to presence of the plastic membrane and water, it can be assumed that the tissue phantom moves rigidly in the lateral direction of the probe. The displacements of several points in the US imaging plane were mea- sured using the envelope block matching. At each time sample the average of measured displacements was compared with the probe displacement. There- fore, at each time sample, the algorithm showed an error with an average and a standard deviation over the entire US image. Table 3.1 shows the max- imum of these values over the time of the experiment (10 seconds for each frequency and amplitude). In our experiments the maximum error occurred at the peak of the displacement and did not increase over time. The error increases with amplitude and frequency as higher speeds of tissue motion decrease the motion estimation accuracy. The only exception occurred when the frequency was 0.25Hz and the amplitude was increased from 4.0mm to 6.0mm. When the amplitude of the movement was 6.0mm the tracking algorithm lost track of some nodes due to the large displacement. These lost nodes resulted in a decrease in the average of the error, compared to the displacement with amplitude of 4.0mm, as positive and negative errors canceled each other. However, these lost nodes increased the standard de- viation of the error, significantly. This table shows the accuracy of the RF envelope block matching in the estimation of large lateral displacements. The linear probe used for the experiment had an array of 128 elements with a density of 2.56 elements per millimeter. The recorded RF frames were up-sampled by a factor of two before the Hilbert transform was applied in order to increase the accuracy of the algorithm. The BIs were three lines in the lateral direction and 56 samples in the axial direction of the image. The SRs were five lines wide and 72 samples long. Increasing the size of BI can decrease the sensitivity of the algorithm to noise at the expense of lower resolution [16]. In the presence of tissue deformation (for example, 63 3.4. Needle Insertion Experiment caused by needle insertion), de-correlation can occur. To compensate for this de-correlation the reference frame was updated every 50th frame. The motion estimation algorithm under-estimates the lateral displace- ment if the real movement is less than the line width of the US image. This inaccuracy is due to the fact that the cross-correlation of one BI in the ref- erence frame and the same size blocks in the SR is not only a function of the spatial proximity of the two blocks, but also a function of their texture sim- ilarity. This fact also leads to non-zero motion estimation when the tissue is stationary. In this situation, the cross-correlation matrix is asymmetric respect to its middle column due to asymmetric texture of the SR. The asymmetric cross-correlation matrix leads to a non-zero motion displacement after the sub-sampling process is applied. 3.4 Needle Insertion Experiment An experiment was conducted to record the insertion force and tissue dis- placement data during insertion of a needle into a two-layered prostate phan- tom. In this experiment, a brachytherapy needle (Bard, USA) was mounted on a linear stage powered by a DC motor equipped with an optical encoder. The axial insertion forces were measured using a 1 DoF load cell (MBD-2.5 Transducer Techniques, USA) housed at the base of the needle. A two- layered prostate phantom was constructed which was composed of a hard inclusion (50% plastic hardner and 50% PVC plasticizer) and a softer sub- strate (33% plastic softner and 67% PVC plasticizer). Cellulose was added to both parts to produce scatterers. A PC-based Sonix RP (Ultrasonix Med- ical Corp., Canada) ultrasound machine was used to image the tissue in the sagittal plane using a TRUS probe located inside a hollow cylinder in the phantom (see Fig. 3.3(a)). The experimental setup is shown in Fig. 3.3(b). The needle was inserted with a constant speed of 5mm/s. The needle stylet was fully inserted to create a symmetric tip. No significant needle bending was observed during the insertion. The TRUS probe imaged the tissue to the depth of 50mm with 128 elements in the lateral direction (60mm) at a rate of 25 frames per second. The needle position, the insertion force and a time stamp were sent to the US machine computer via a serial port. In a loop, the US machine records the RF data and then reads the values of the load cell, needle position and the accompanying time stamp. Since the tissue phantom was meshed to be used in an FEM analysis, the displacements of mesh nodes located on the US plane were measured in 2D 64 3.4. Needle Insertion Experiment PVCPhantom Water Probe Movement Direction Plastic Membrane Figure 3.2: The relative position of probe and its movement respect to the tissue phantom. 65 3.4. Needle Insertion Experiment Table 3.1: Validation results for lateral motion estimation algorithm. Frequency(Hz) 0.25 0.25 0.25 0.5 1.0 Peak to peak(mm) 2.0 4.0 6.0 2.0 2.0 Max. Ave. Error(mm) 0.11 0.26 0.23 0.14 0.30 Max. StD. Error(mm) 0.07 0.12 0.40 0.08 0.11 66 3.4. Needle Insertion Experiment USProbe Insertion Direction Z Y US Plane US gel 1 2 (a) Tissue Phantom Ultrasound Probe Needle Load Cell DC Motor & Encoder (b) Figure 3.3: (a) The cross-section of the tissue phantom, showing the US field of view. Arrows 1 and 2 show the lateral and axial directions in the US image, respectively. In the shown configuration, the US plane coincides with the sagittal plane of the tissue. (b) The experimental setup. 67 3.5. Results using the block matching algorithm. The lateral and axial displacement of the mesh nodes as functions of time are shown in Fig. 3.4. 3.5 Results Similarly to the method presented in the previous chapter, the experimental data were used in a linear FEM-based needle insertion simulator. The three parameters of the force distribution model (peak force, shaft force, and width of the peak force) were adopted to match the simulated and measured insertion forces. The Young’s moduli of the two parts of the phantom were identified to match the average of the simulated and measured displacement along the y axis (the major displacement). The PVC phantom Poisson’s ratio was taken to be 0.49 to simulate a nearly incompressible material. The identified parameters are presented in Table 3.2. The simulated and measured forces are shown in Fig. 3.5(a). The original positions of the mesh nodes in the US image and the measured and simulated displaced positions are shown in Figs. 3.5(b) and 3.5(c). The accuracy of the identified parameters depends on the accuracy of the simulation and the measured data. The accuracy of the simulation can be increased by using finer meshes and/or non-linear models (i.e. neo- Hookean) at the cost of higher computational time. The identified Young’s moduli may be different from the real Young’s moduli of the tissue. However, for simulation purposes, they are acceptable as long as they yield accurate simulation results. 3.6 Conclusions and Future Work In this chapter a 2D block matching algorithm was used to measure the displacement of a tissue phantom during needle insertion. The algorithm estimates the displacement using 2D correlation maximization of the RF envelopes of the reference and target frames. The accuracy of the algorithm in measuring the lateral motion was tested during an experiment in which a linear probe was used to measure the rigid lateral motion of a PVC phantom. The 2D tracking algorithm showed high accuracy in tracking the large lateral displacements. However, lateral displacements of less than one line width were under-estimated. Without sub-sampling in the lateral direction, the algorithm shows sharp jumps when the displacement is larger than one line width. Sub-sampling creates a continuous output. Due to under-estimation of small displacements, the reference frame update rate should be chosen 68 3.6. Conclusions and Future Work 0 10 20 30 40 −2 0 2 4 6 8 10 12 Needle Position (mm) La te ra l D is pl ac em en t (m m) (a) 0 10 20 30 40 −1 −0.5 0 0.5 1 1.5 Needle Position (mm) Ax ia l D is pa lce m en t (m m) (b) Figure 3.4: Displacements of the mesh nodes (a) along the axis of insertion (lateral direction of the US image) and (b) along the axial direction of the US images. 69 3.6. Conclusions and Future Work 0 10 20 30 40 0 0.5 1 1.5 2 2.5 3 3.5 4 Needle Position (mm) Fo rc e (N ) (a) 40 60 80 100 120 500 1000 1500 2000 2500 (b) 40 60 80 100 120 500 1000 1500 2000 2500 Needle (c) Figure 3.5: (a) Simulated and measured insertion forces vs. needle tip position (needles moves from the left to the right), (b) The original positions of the mesh nodes and (c) the measured (squares) and simulated (circles) positions of the displaced nodes. 70 3.6. Conclusions and Future Work Table 3.2: Needle shaft force distribution and elastic parameters. fs(N/m) fp(N/m) w(mm) E(kPa) Inclusion 90 380 4.5 30 Surrounding tissue 50 100 4.0 7 71 3.6. Conclusions and Future Work carefully. Very high update rates lead to large errors, since at each time sample the displacement of the current frame with respect to the reference frame can be less than one line. Very low update rates lead to de-correlation and poor tracking. The optimal reference update rate requires further study. During a needle insertion experiment a two-layered PVC phantom was imaged using a TRUS probe in the sagittal plane. The tissue displacement data, the needle position and the insertion force measured during the exper- iment were fed to an FEM needle insertion simulator to identify the param- eters of a needle-tissue interaction model. The results show good agreement between the simulated and measured force and displacements. A PVC phantom was used to make the experiment repeatable. However, insertion forces in this phantom do not show the sharp dips observable during insertion of the needle into real tissue. This may be due to the absence of a hard capsule in the phantoms. The extension of this work to real tissue will be undertaken in the future. In the future, this method will be used during prostate brachytherapy to measure the prostate rotation and deformation during needle insertion. The in-vivo measurements can be used later to generate a needle-tissue interaction model to be used in brachytherapy simulators and path planners. 72 References [1] J. E. Dawson, T. Wu, T. Roy, J. Y. Gy, and H. Kim, “Dose effects of seeds placement deviation from pre-planned position in ultrasound guided prostate implants,” Radiotherapy and Oncology, vol. 32, no. 2, 1994. [2] V. Lagerburg, M. A. Moerland, J. J. Lagendijk, and J. J. Batter- mann, “Measurement of prostate rotation during insertion of needles for brachytherapy,” Radiotherapy and Oncology, vol. 77, no. 3, pp. 318– 323, 2005. [3] S. P. DiMaio and S. E. Salcudean, “Needle insertion modeling and simulation,” IEEE Trans. on Robot. Autom.: Special Issue on Medical Robotics, vol. 19, no. 5, pp. 864–875, 2003. [4] J. R. Crouch, C. M. Schneider, J. Wainer, and A. M. Okamura, “A velocity-dependent model for needle insertion in soft tissue,” in Proc: Medical Image Computing and Computer Assisted Intervention (MIC- CAI), 2005, pp. 624–632. [5] E. Dehghan, X. Wen, R. Zahiri-Azar, M. Marchal, and S. E. Salcudean, “Needle-tissue interaction modeling using ultrasound-based motion es- timation: Phantom study,” in Computer Aided Surgery, vol. 13, no. 5, pp. 265-280, 2008 [6] A. Okamura, C. Simone, and M. O’Leary, “Force modeling for needle insertion into soft tissue,” IEEE Trans. Biomed. Eng., vol. 51, no. 10, pp. 1707–1716, 2004. [7] J. T. Hing, A. D. Brooks, and J. P. Desai, “A biplanar fluoroscopic approach for the measurement, modeling , and simulation of needle and soft-tissue interaction,” Medical Image Analysis, vol. 11, no. 1, pp. 62–78, 2007. [8] T. Podder, J. Sherman, E. Messing, D. Rubens, D. Fuller, J. Strang, R. Brasacchio, and Y. Yu, “Needle insertion force estimation model 73 Chapter 3. References using procedure-specific and patient-specific criteria,” in Proc. IEEE EMBS Int. Conf., 2006, pp. 555–558. [9] R. Zahiri-Azar and S. E. Salcudean, “Motion estimation in ultrasound images using time domain cross correlation with prior estimates,” IEEE Trans. Biomedical Engineering, vol. 53, no. 10, pp. 1990–2000, 2006. [10] ——, “Real-time estimation of lateral motion using time domain cross correlation with prior estimates,” in Proc. IEEE Ultrasonics Sympo- sium, 2006, pp. 1209–1212. [11] D. E. Robinson, F. Chen, and L. S. Wilson, “Measurement of velocity of propagation from ultrasound pulse-echo data,” Ultrasound in Medicine and Biology, vol. 8, no. 4, pp. 413–420, 1982. [12] I. Akiyama, A. Hayama, and M. Nakjima, “Movement analysis of soft tissue by speckle patterns’ fluctuation,” JSUM Proc., 1986, pp.615–616. [13] E. J. Chen, R. S. Adler, P. L. Carson, W. K. Jenkins, and W. D. O’Brien, “Ultrasound tissue displacement imaging with application to breast cancer,” Ultrasound in Medicine and Biology, vol. 21, pp. 1135– 1162, 1995. [14] T. Varghese, and J. Ophir, “Characterization of elastographic noise using the envelope of echo signals,” Ultrasound in Medicine and Biology, vol. 24, no. 4, pp. 543–555, 1998. [15] I. Cespedes, Y. Huang, J. Ophir, and S. Spratt, “Methods for the esti- mation of subsample time-delays of digitized echo signals,” Ultrasonic Imaging, vol. 17, no. 2, pp. 142–171, 1995. [16] B. S. Ramamurthy and G. E. Trahey, “Potential and limitations of angle-dependent flow detection algorithms using radio-frequency and detected echo signals,” Ultrasonic Imaging, vol. 13, pp. 252–268, 1991. 74 Chapter 4 Needle Insertion Parameter Optimization for Brachytherapy 3 4.1 Introduction Needle insertion is involved in many medical interventions. In spite of the use of image guidance, accurate needle placement is still a challenge in many procedures. One example illustrating the need for accurate placement is low dose rate prostate brachytherapy. This effective treatment for local- ized prostate cancer consists of implanting many tiny radioactive capsules or “seeds” of 125I or 103Pd inside the prostate and the peri-prostatic tissue using long needles (see Fig.2.1). In order to kill the cancerous tissue while maintaining a tolerable radiation dose on other organs, the seed distribution is carefully planned prior to the intervention using brachytherapy treatment planning software. The planning software requires a 3D prostate model, typ- ically constructed from several parallel transversal trans-rectal ultrasound (TRUS) images which are manually segmented. This 3D model is used to plan for the seed positions which are usually along straight lines parallel to the long axis of the ultrasound probe (the z axis in Fig. 2.1). In the operat- ing room, the physician registers the real-time TRUS images of the prostate to the pre-operatively acquired images used for planning and then implants the seeds, typically using pre-loaded needles. To deliver the seeds according to the plan, the physician inserts the needle through a guiding template - a grid with equidistant 5mm apart parallel guidance holes in the x− y plane. The grid confines the needle movement along the z axis. Real time ultra- sound and occasional X-ray fluoroscopy are used for visual feedback during the operation. 3A version of this chapter has been published as: E. Dehghan, and S. E. Salcudean, “Needle Insertion Parameter Optimization for Brachytherapy ”. IEEE Transactions on Robotics, vol. 20, no. 2, pp. 303–315, 2009 75 4.1. Introduction During the insertion, the prostate moves and deforms due to needle forces. The rotation is caused by the configuration of boundary conditions (pubic arch, bladder) and tissue deformation parameters and is likely caused by a stiffer anterior anchoring of the prostate to the pubic arch. The prostate displacement during needle insertion is illustrated in Fig. 4.1 using a US image of the region in the transverse view. At first, only the bladder can be seen in the US image as in Fig. 4.1(a). When the needle penetrates more into the prostate, it pushes the prostate back into the bladder; hence, in Fig. 4.1(b) the prostate and the needle tip are visible in the US image. Lagerburg et al. measured the rotation of the prostate during brachytherapy and found it to be significant, often in excess of 10 degrees [1]. As a result of prostate motion, seed misplacements are still common in brachytherapy [2]. Misplacements can lead to under-dosed regions that may require re- peated procedures or to higher doses that can cause complications - such as impotence or urinary incontinence. The guiding template allows the needle to move along the z axis only. However, the target positions diverge from this line during insertion due to tissue deformation. Modeling the tissue deformation and allowing for an arbitrary needle insertion direction as opposed to only parallel directions may allow for more accurate delivery of seeds. A number of brachytherapy robots have been developed (e.g., [4, 5, 6, 7, 8, 9]) that are able to insert the needle with desired pre-planned ori- entation. However, the problems of pre-planning for brachytherapy in the presence of tissue deformation and path planning for brachytherapy robots have not been solved. A few methods have been introduced to optimize the needle insertion depth, initial position and/or orientation for 2D deformable tissue models with single targets [10, 11]. However, no method has been reported previously for insertion point and heading optimization for needle path planning in 3D and in the presence of multiple targets. This paper proposes such a method based on 3D simulations that use the Finite Element Method (FEM). The paper is organized as follows: an overview of the previous work on needle insertion is presented in Section 4.2; the details of the optimization method are described in Section 4.3, followed by the description of the simulation method in Section 4.4. Needle insertion planning simulation results for a brachytherapy simulation and experimental validation on a tissue phantom are presented in Sections 4.5 and 4.6, respectively, followed by a discussion in Section 4.7. Conclusions and future work are presented in Section 4.8. 76 4.2. Related Work Prostate Needle Bladder TRUS imaging plane Bladder Prostate Needle TRUS imaging plane Bladder Novisible prostate (a) Bladder Soft Tissue Prostate (b) Figure 4.1: The prostate displacement during needle insertion. The image is taken from [3]. 77 4.2. Related Work 4.2 Related Work Needle insertion into soft tissue has been a topic of many research activities in the past decade (for a survey please see [12]). Among different methods available, the Finite Element Method – a powerful approximate solution to solid mechanics problems – has been widely used to simulate tissue defor- mation during surgery and needle insertion [10, 13, 14, 17, 18]. Alterovitz et al. [17] simulated brachytherapy using a 2D linear FEM based model. Goksel et al. [3, 14] generated a 3D mesh of the prostate – constructed from manually segmented transverse images – and used it in a 3D brachytherapy simulator. Computation acceleration techniques have been presented for 2D and 3D linear FEM based needle insertion models to achieve high speeds [13, 14]. The effects of tissue non-linearities, boundary conditions and organ geometry on the outcome of the insertion simulators and path planners were investigated in [15] and [16]. Alterovitz et al. optimized the needle insertion height and depth to de- crease the seed misplacement error for a single target, in a 2D brachyther- apy simulator [10]. They used a search-based algorithm to optimize the two parameters. For optimization of the insertion depth, they simulated the insertion of a needle with a given insertion height to the maximum depth allowed. The optimal depth for the given height is the depth with minimum error. To optimize the insertion height they used the golden section search method. Extension of this work to 3D will result in a system with three more optimization parameters (width of the insertion point and two an- gles for needle orientation). Optimizing for these parameters using a search based algorithm requires a large number of simulations. In another work, Alterovitz et al. optimized the initial point and angle for insertion of a highly flexible needle into a 2D tissue to reach a target while avoiding ob- stacles [11]. Their technique involves the transformation of the constrained optimization problem into an unconstrained one using a penalty method, then applying a gradient descent algorithm to solve the problem. Due to its complexity, the cost function cannot be directly differentiated with present techniques. Therefore, the authors used a perturbation method to find the approximate derivatives of the cost function at each iteration, each requir- ing multiple simulations. As opposed to these gradient-based approaches, the algorithm presented in this paper does not need the derivatives of the cost function and addresses the 3D problem. An average plan for prostate brachytherapy treatment involves 20-30 needle insertions. The planning al- gorithm can be used pre-operatively to optimize the insertion parameters in order to achieve accurate implants. Furthermore, a number of studies have 78 4.2. Related Work shown that there may be a benefit in terms of improved outcomes from the simultaneous planning and execution of needle insertions intra-operatively, as a function of intra-operatively computed dose computations [19]. Thus the need for high-speed planning is clear. The needle model used in [11] is an nonholonomic model fitted to a highly flexible bevel-tip needle by Webster et al. [20]. Alterovitz et al. [21] used the same needle model to steer the needle tip inside a 2D tissue model under Markov uncertainty. They used the needle bevel tip direction as input and guided the needle to reach a target while avoiding obstacles. However, they did not optimize the initial needle position or orientation. Park et al. for- mulated the problem of steering a highly flexible needle as a nonholonomic kinematics problem and addressed the planning problem by diffusion-based motion planning on the Euclidean group, SE(3) [22]. However, they did not simulate the tissue deformation caused by insertion forces. The highly flexi- ble needle model, which assumes the needle follows its beveled tip trajectory direction without applying considerable lateral force to the tissue is not ap- plicable to brachytherapy needles. In Chapter 5, three different models for brachytherapy needles were compared. DiMaio and Salcudean used a 2D non-linear FEM based model for the needle and manipulated the needle base to steer the tip inside tissue [24]. They related the movements of the needle tip to the movements of the base using a numerically computed needle Jacobian. They steered the flexible needle to hit a target while avoiding obstacles using potential fields. Using a spring mesh for the tissue and a linear beam model for the needle, Glozman and Shoham steered a needle inside the tissue by using the inverse kinematics of the needle [25]. The initial location and orientation of the needle were not optimized in [24] or [25]. Needle rotation and insertion speed were used to reduce tissue defor- mation during insertion in [26] and [27]. Heverly et al. [27] showed that velocity dependent properties of the tissue can be used to reduce tissue de- formations at high needle insertion speed, essentially hammering the needle into prostate. In contrast to other published work, this paper proposes a method to optimize the needle insertion point and heading into a 3D tissue model containing multiple targets. The presented method which, for now, assumes a rigid needle, uses the simulation results to modify the insertion parameters. To simulate the tissue deformation as a result of needle insertion, in addition to a model for the tissue, a model for needle-tissue interaction needs to be developed. Okamura et al. [28] divided the forces applied by the needle into three parts: 1) a force due to capsule stiffness occurring 79 4.3. Optimization Method before puncture of the capsule; 2) a friction force acting on the needle shaft; and 3) a cutting force on the needle tip. They modeled the capsule stiffness force with a non-linear spring, the friction force with a modified Karnopp model, and the cutting force with a constant for a given tissue. DiMaio and Salcudean [13] introduced the stick-slip needle tissue interaction model, based on the force measurements during insertion of a needle into a slab of PVC. They used this model to accurately predict insertion trajectories in PVC phantoms, and also to achieve high computational speed. Kataoka et al. [30] measured the force applied to a needle during insertion into a canine prostate. Needle motion and forces during brachytherapy of several patients were reported by Podder et al. in [31]. 4.3 Optimization Method As done in prior work - see, for example, [32], [10]- the optimization method in this paper is based on simulations. In each iteration of the optimization, the simulation program finds the location of the displaced targets after the deformation caused by needle insertion with given initial orientation, tip location and depth. Therefore, the simulation acts as a function in which the needle insertion point ps, needle heading – denoted as a unit vector v – and needle insertion depth d are the input parameters and the target dis- placed positions are the outputs. Other parameters such as tissue elasticity parameters, target initial positions and needle-tissue interaction parameters are assumed to be known constants in this function. Assuming N targets located in the tissue, the input-output relation can be written as: ui = Fi(ps, v, d) i ∈ {1, 2, · · · , N}, (4.1) where ui are the displaced locations of the targets (uN is the displaced position of the distal target). The optimization goal is to minimize the distance between a rigid needle and the displaced targets in the deformed configuration, by optimizing the needle insertion point, heading and depth, i.e.: (p∗s , v ∗, d∗) = arg min ps,v,d ( N∑ i=1 min α ‖ps + αv − ui‖ 2 ) (4.2) Subject to ui = Fi(ps, v, d), α ≤ d, ‖v‖ = 1 80 4.3. Optimization Method The constraint α ≤ d assures that for the targets deeper than the needle tip, the cost function is calculated as the distance of these targets from the needle tip. Assume a global coordinate system Oxyz and a needle- attached coordinate system O N x N y N z N as shown in Fig. 4.2. Since the needle is axially symmetric, the pitch and yaw angles are sufficient to relate the bases of two coordinate systems. Therefore, the two angles, pitch and yaw (interchangeable with v), the needle insertion point on any given plane ( for example the guiding template in brachytherapy or alternatively the front surface of the tissue) and the insertion depth are the parameters to be optimized. Assuming that the targets are originally located on a straight line, the following iterative algorithm can be used to optimize the insertion parameters. 1. Initialize the parameters by passing a line p1 + αv1 through the tar- gets in the undeformed configuration. The insertion point p1s is the intersection of this line and a desired plane (the front surface of the tissue in this paper). The insertion depth is the distance between the insertion point and the distal target. The parameters ps and v at the first iteration are shown in Fig. 4.4(a). 2. Simulate the needle insertion with the given needle heading parame- ters. The simulation program yields the position of the targets (uki = Fi(p k s , v k, dk), i = 1, · · · , N) in the deformed configuration. The su- perscript k is the index of the iteration. This step is illustrated in Figs. 4.4(b) and 4.4(c). 3. Find the new insertion line parameters (vk+1 and pk+1) by fitting a line through the target positions in the deformed configuration (see Figs. 4.4(b) and 4.4(c)). 4. Find the new insertion point pk+1s as the intersection of this line and the desired plane. 5. Find the new insertion depth dk+1 as the distance between pk+1s and the closest point along the insertion line to the distal target. 6. If the convergence criterion is met, then terminate the algorithm. Oth- erwise go to step 2. In Step 1), the line is trivial to compute assuming no tissue deformation. The new line parameters in Step 3) can be calculated by optimizing the line 81 4.3. Optimization Method z y x Nx Ny Nz O NO Figure 4.2: The simplified mesh of prostate and surrounding tissue. The area connected to the pubic bone is shown in black color. 82 4.3. Optimization Method parameters pk+1 and vk+1 to minimize the sum of distances between the displaced targets and the line. Therefore: (pk+1, vk+1) = arg min p,v ( N∑ i=1 min α ‖p+ αv − uki ‖ 2 ) , (4.3) Subject to ‖v‖ = 1 It is shown in AppendixA that: pk+1 = 1 N N∑ i=1 uki , (4.4) and vk+1 is the right singular vector corresponding to the largest singular value of A, where A = [ uk1 − p k+1 uk2 − p k+1 · · · ukN − p k+1 ]′ , (4.5) and (·)′ denotes the transpose of the operand. The convergence criterion used is the satisfaction of a bound on the sum of distances between the needle and the targets in the deformed tissue. In general, due to tissue deformation and target displacement, the targets will not lie on a straight line after needle insertion. Therefore, the above mentioned convergence criterion usually does not converge to zero. Alternatively, the Euclidean norm of differences in the vector of optimization parameters can be considered as a convergence criterion. With this convergence criterion, if the norm of the difference between two consecutive vectors of parameters is smaller than a threshold, the algorithm is terminated (see Fig. 4.4(d)). The thresholds in the above criteria can be determined by experience and to match actual geometrical parameters. For example, there is no need to drive the targeting error to less than the needle diameter or the parameter convergence to less than the robot resolution. In fact, organ registration errors are more likely to determine the actual thresholds used in practice. Using this convergence criteria, the optimization and simulation program can be seen as a single function in which the insertion parameters at iteration k are the inputs and the insertion parameters at iteration k + 1 are the outputs. The goal is to find a fixed point of this function. 83 4.4. Performance in Simulation 4.4 Performance in Simulation To demonstrate the feasibility of the optimization algorithm, extensive sim- ulations were carried out. The optimization method was implemented in a brachytherapy simulator to reduce the targeting errors caused by defor- mation. A simplified mesh of the prostate and surrounding tissue [3, 14] – shown in Fig. 4.2 – is used in the simulations. For simplicity, it is assumed that the surrounding tissue is homogeneous. Insertion of a rigid needle into the 3D tissue model is simulated using an FEM based model in a quasi-static mode. The needle is modeled as a straight line without thickness. The needle-tissue interaction is described by a stick-slip force model as in [13] due to its simplicity of implementation and computational speed. In this model a node is stuck to the needle as long as its reaction force is smaller than a predefined threshold. In this state, three displacement boundary conditions are applied to the stuck node. If the reaction force is greater than the threshold, the node state changes to slip mode. In this state the node can slide along the needle. A force boundary condition is applied to this node in the needle axial direction. To confine the node movement to be along the needle, two displacement boundary conditions are applied orthogonally to the needle axis. A mesh adaptation algorithm [11, 14] is used to increase the accuracy of the simulation. In this method, when the needle tip arrives at a new element boundary, the closest node to the needle tip is re-located on the tip. To do so, the undeformed mesh is adapted in such a way that at the last de- formation configuration the corresponding node position coincides with the needle tip. No external force is applied to re-locate the node. Mesh adapta- tion does not lead to inverted elements since the modification is applied to the reference mesh. The stiffness matrix is regenerated accordingly. When using linear FEM, the inverse stiffness matrix can be modified directly, with low computational cost. More details about mesh adaptation algorithm and low cost inverse stiffness matrix modification are available in [13] and [14]. During the insertion simulation, the displaced location of the targets are of interest and should be computed. Therefore, it is sufficient to compute the displacement of the nodes in contact with the needle and the nodes at the vertices of elements in which the targets are located. When linear FEM is used, this computational shortcut called condensation can be exploited to increase the simulation speed significantly [13]. In this case, the constant global stiffness matrix can be pre-computed, inverted and saved in the mem- ory. During the simulation, only a small part of the inverse stiffness matrix which corresponds to the nodes of interest – the operating matrix – is en- 84 4.5. Simulation Results gaged in the calculations [13, 14]. However, in non-linear FEM which leads to a non-linear set of algebraic equations, the displacements of all nodes should be computed at each iteration. 4.5 Simulation Results Nine sets of targets were assumed inside the prostate model, each of which consisted of three targets located along the z axis – inspired by the location of the targets in brachytherapy. In each set, the targets were located at depth z = 15, z = 25, and z = 35mm in the undeformed tissue (see for example Fig. 4.4(a)). The x− y coordinates of the target sets are shown in Fig. 4.3 as a cross section of the prostate at z = 25mm. Although the organs surrounding the prostate are known with good de- tail from Magnetic Resonance Imaging or anatomy literature, the boundary condition of the prostate in the tissue is not accurately known. This is due to the fact that the interaction of the prostate with the surrounding or- gans during the deformation is not well known and depends on how firmly or loosely it is attached to the surrounding tissue. Other factors such as bladder filling can affect the movement of the prostate and impose compli- cated boundary conditions. However, the prostate is more stiffly anchored around the pubic arch and therefore rotates during needle insertion due to needle torques about the anchor region. To model this, a number of prostate mesh nodes in the anterior/superior region are fixed and therefore determine displacement boundary conditions as marked in Fig. 4.2. Different bound- ary conditions can lead to different simulation results and consequently to different optimized parameters. The optimization algorithm is tested for different boundary condition assumptions and tissue models. The tissue elastic parameters - the Young’s modulus E and the Poisson’s ratio ν – were typical of soft tissue [33] and are presented, together with other simulation parameters, in Table 4.1. 4.5.1 Posterior Surface Fixed In this case, in addition to nodes in connection with the pubic arch, the nodes located at the back of the mesh are fixed. Therefore, the prostate motion is confined and hence, one can expect small rotation of the prostate. The parameter optimization was performed for each of the target sets. The iteration was continued until the change in the vector of optimization parameters was less than 1% or the error – defined as the sum of distances 85 4.5. Simulation Results 8070605040302010 15 20 25 30 35 40 45 50 X(mm) Y( mm ) 89 7 456 3 2 1 Figure 4.3: The x − y position of the target sets in a cross section of the prostate at z = 25mm. 86 4.5. Simulation Results Table 4.1: Simulation parameters for prostate model. (E, ν) for the prostate tissue (60 kPa, 0.49) (E, ν) for the surrounding tissue (20 kPa, 0.49) Number of Nodes and Elements 570, 2801 Element Type Tetrahedron Simulation type Linear elasticity, quasi-static Table 4.2: Simulation results for case A. Posterior mesh surface is fixed. Set Num. 1st Err. Last Err. Opt. angles Opt. insertion Opt. depth Itrs. (mm) (mm) (yaw, pitch) offset (x, y) (mm) 1 2 1.8 0.8 (0.5o, 0.4o) (0.7, -0.2) 56.4 2 2 1.6 0.4 (0.1o, 0.6o) (0.1, -0.6) 56.3 3 2 1.8 0.5 (0.2o, 1.3o) (-0.2, -1.1) 56.4 4 2 1.8 0.5 (0.1o, 0.4o) (0.1, -0.4) 56.4 5 2 1.9 0.3 (0.1o, 0.1o) (-0.2, -0.1) 56.3 6 2 1.4 0.2 (0.2o, 0.8o) (0.2, -0.7) 56.1 7 2 1.4 0.4 (0.3o, 0.0o) (0.3, 0.0) 56.2 8 2 1.4 0.4 (0.5o, 0.0o) (0.4, 0.0) 56.2 9 2 1.1 0.1 (0.0o, 0.2o) (0.1, -0.1) 55.6 87 4.5. Simulation Results between the needle and the targets – was less than the brachytherapy nee- dle diameter (1.28 mm). The linear FEM was used for the simulation. The errors in the first and last iterations, the number of iterations and the op- timized parameters are reported in Table 4.2 for each target set. It can be seen that the algorithm converged to sub-millimeter errors in two iterations for all the target sets. This is due to the small rotations of the prostate caused by the confined boundary conditions. 4.5.2 Superior Surface Fixed In the second boundary condition scenario, the posterior mesh nodes are free while the nodes located on the top of the mesh are fixed. The simu- lation is performed with the same elastic parameters as in the former case. Simulation results are summarized in Table 4.3. In the two simulated boundary condition scenarios the tissue is confined and it cannot easily rotate or displace. Therefore, the optimal insertion parameters are close to the initial values, especially for the targets close to the pubic bone where the tissue is more confined (e.g. target sets 8 and 9). 4.5.3 Pubic Arch Boundary Condition Only In another case, the posterior and superior surfaces of the mesh are both free and the prostate is only confined in the area of the pubic arch. In this case the prostate is able to rotate and deform more than in the two former cases. Therefore, the number of iterations necessary for the optimization algorithm convergence increased. Also, the optimized insertion parameters are considerably different from the initial insertion parameters. However for target sets 8 and 9 the rotation and consequent displacements are low, since these target sets are close to the center of rotation and therefore, needle forces produce less torque for these targets compared to targets farther from the center of rotation. The simulation results are summarized in Table 4.4. The optimization algorithm takes 2 to 5 iterations to converge to sub- millimeter errors. 4.5.4 Lower Young’s Modulus for the Inclusion To allow the prostate to deform even more, the elastic modulus of the in- clusion is decreased to E=40 kPa to represent a softer tissue. The prostate is only fixed at the connection with the pubic bone. The simulation results in this case are summarized in Table 4.5. 88 4.5. Simulation Results Table 4.3: Simulation results for case B. Superior mesh surface is fixed. Set Num. 1st Err. Last Err. Opt. angles Opt. insertion Opt. depth Itrs. (mm) (mm) (yaw, pitch) offset (x, y) (mm) 1 2 5.0 0.7 (1.8o, 1.9o) (1.2, -1.1) 59.6 2 2 4.2 0.4 (0.3o, 2.4o) (0.3, -1.7) 59.2 3 2 5.0 0.6 (1.9o, 2.7o) (0.5, -2.1) 59.8 4 2 3.5 0.7 (1.5o, 1.3o) (0.9, -0.9) 58.5 5 2 3.1 0.7 (0.9o, 0.9o) (0.4, -0.5) 58.1 6 2 2.8 0.4 (0.3o, 1.3o) (0.1, -0.9) 57.6 7 2 3.4 0.5 (2.0o, 0.6o) (1.2, -0.5) 58.2 8 2 2.9 0.3 (1.3o, 0.3o) (0.8, -0.2) 57.4 9 2 1.7 0.2 (0.0o, 0.3o) (0.1, -0.1) 56.5 Table 4.4: Simulation results for case C. The tissue is fixed only anteriorly near the pubic arch. E=60 kPa. Set Num. 1st Err. Last Err. Opt. angles Opt. insertion Opt. depth Itrs. (mm) (mm) (yaw, pitch) offset (x, y) (mm) 1 3 13.5 0.6 (3.6o, 9.8o) (2.6, -5.8) 64.6 2 5 11.3 0.9 (0.9o, 10.7o) (0.7, -6.8) 65.4 3 3 12.3 0.7 (4.9o, 9.3o) (3.3, -5.6) 64.6 4 3 8.0 0.5 (4.3o, 6.6o) (2.6, -4.1) 61.3 5 3 7.1 0.3 (3.2o, 6.5o) (1.8, -3.9) 60.4 6 3 6.0 0.4 (0.9o, 6.0o) (0.6, -3.7) 59.5 7 3 7.4 0.4 (6.3o, 3.9o) (3.8, -2.0) 59.8 8 2 4.4 0.9 (4.3o, 1.8o) (2.6, -1.1) 58.7 9 2 2.4 0.5 (0.5o, 2.5o) (0.4, -1.5) 56.9 89 4.5. Simulation Results Table 4.5: Simulation results for case D. The tissue is fixed only anteriorly near the pubic arch. E=40kPa. Set Num. 1st Err. Last Err. Opt. angles Opt. insertion Opt. depth Itrs. (mm) (mm) (yaw, pitch) offset (x, y) (mm) 1 3 15.7 0.7 (3.9o, 10.6o) (2.8, -6.3) 66.1 2 7 13.3 1.1 (1.1o, 11.6o) (0.8, -7.4) 65.6 3 4 13.5 0.7 (5.0o, 10.4o) (3.3, -6.4) 65.9 4 4 9.6 0.5 (4.7o, 7.4o) (2.7, -4.4) 60.7 5 3 8.9 0.6 (3.7o, 6.9o) (2.0, -1.5) 61.7 6 3 5.8 0.6 (0.8o, 5.9o) (0.6, -3.6) 58.4 7 3 8.7 0.7 (6.8o, 3.4o) (4.2, -1.8) 61.8 8 3 5.1 0.9 (4.9o, 1.8o) (3.0, -1.0) 59.5 9 2 2.9 0.6 (0.5o, 2.7o) (0.4, -1.6) 57.5 90 4.5. Simulation Results The optimization method shows stability and successful decrease in the targeting error in a maximum of 7 iterations. As an example, the optimiza- tion steps for the 5th target set in this case are shown in Fig. 4.4. Fig- ure 4.4(a) shows the position of the targets in the undeformed configuration inside the prostate. The first candidate for the insertion line is shown as the line passing through the targets. Figures 4.4(b) and 4.4(c) show the de- formed configuration after the first and the second iterations. The displaced position of the targets, the needle and the new insertion line are shown in each figure. Figure 4.4(d) shows the final iteration in which the targets are very close to the needle and the next insertion path is very close to the old one. Comparing Tables 4.5 and 4.4 shows that in spite of a 33% decrease in the Young’s moduli of the prostate in this case, the optimized parameters are not significantly different. Based on this information, one can deduce that optimizing the needle insertion parameters can be helpful even in the presence of uncertainaty in the elastic properties of the tissue. Decreasing the Young’s modulus and increasing the friction forces have similar effects. Therefore, these cases are not presented here. 4.5.5 Neo-Hookean Material Model for Tissue When the prostate is mainly fixed by the pubic bone only, the tissue can un- dergo large rotations and deformations. To simulate these conditions more accurately a non-linear FEM based model is used which can accommodate both the geometric and mechanical non-linearities. Neo-Hookean hyperelas- tic material models [34] are assumed for the prostate and the surrounding tissue. Neo-Hookean material models were previously used to simulate tissue during needle insertion [35]. They have been used also for tissue modeling for liver simulation [36] and brain image registration [37]. Lamé parameters (λ and µ) are usually used in constitutive equations for neo-Hookean ma- terial models (please see AppendixB). The material property parameters for the prostate model are E =60kPa for the inclusion and E = 20 kPa for the surrounding tissue (ν = 0.49) which corresponds to λ = 985.6 kPa and µ = 20.1 kPa for the inclusion and λ = 328.8 kPa and µ = 6.7 kPa for the surroundings. The Newton-Raphson iterative method [34] is used to solve the non-linear set of algebraic equations. Due to the non-linearity of the problem, pre-computation and condensation cannot be used and the displacement of all the nodes should be computed at each iteration. The optimization is performed for all the 9 target sets and the results are summa- rized in Table 4.6. The optimization algorithm converged in 3 to 6 iterations 91 4.5. Simulation Results 100 50 0 −10 0 20 40 60 80 0 20 40 60 70 X(mm) Z(mm) Y( mm ) Insertion line Insertion depth p s 1 v1 (a) The undeformed configuration. 100 50 0 −20 0 20 40 60 80 0 20 40 60 80 X(mm) Z(mm) Y( mm ) New insertion lineNew insertion depth (b) After the first iteration. 100 50 0 −20 0 20 40 60 80 0 20 40 60 80 X(mm) Z(mm) Y( mm ) New insertion line New insertion depth (c) After the second iteration. 100 50 0 −20 0 20 40 60 80 0 20 40 60 80 X(mm) Z(mm) Y( mm ) (d) After the third iteration. Figure 4.4: Optimization iterations for three targets. (b), (c) and (d) show the position of the targets in the deformed tissue after insertion of the needle with the insertion parameters calculated in (a), (b) and (c), respectively. 3D fitted line is shown as a dotted line. The front surface of the prostate mesh is removed to show the position of the targets and the needle inside. 92 4.5. Simulation Results Table 4.6: Simulation results for case E. Neo-Hookean material model. Set Num. 1st Err. Last Err. Opt. angles Opt. insertion Opt. depth Itrs. (mm) (mm) (yaw, pitch) offset (x, y) (mm) 1 3 15.5 0.6 (4.3o, 11.0o) (2.9, -6.3) 64.5 2 4 12.1 0.8 (1.0o, 11.7o) (0.2, -7.0) 64.8 3 6 13.0 0.4 (5.6o, 10.1o) (3.5, -6.0) 64.0 4 3 7.4 0.7 (4.6o, 6.7o) (2.6, -3.9) 61.0 5 3 7.5 0.8 (2.9o, 7.8o) (1.9, -4.1) 59.4 6 3 6.3 0.6 (0.9o, 6.2o) (0.7, -3.7) 59.3 7 3 7.8 0.3 (6.1o, 4.2o) (3.6, -2.3) 60.5 8 3 4.4 0.7 (4.2o, 1.8o) (2.5, -1.1) 58.6 9 3 2.4 0.5 (0.4o, 2.3o) (0.4, -1.4) 56.9 93 4.6. Experimental Validation to sub-millimeter errors. 4.6 Experimental Validation In order to experimentally validate the optimization algorithm, a prostate phantom was constructed. The tissue phantom is composed of a prostate mimicking inclusion made of a mixture of PVC plasticizer and plastic hard- ener surrounded by a substrate made of a mixture of PVC plasticizer and plastic softener. The inclusion consists of a cylinder with two hemispheres at its ends, connected to the base of the phantom with another cylinder of the same material. The combination of the inclusion and its connective cylinder mimics the prostate and its rotation around the pubic bone. A hollow cylinder was provided in the phantom to simulate the rectum. A cross-section of the tissue phantom is shown in Fig. 4.5. The PVC was chosen over bovine gels and agar phantoms because it produces insertion forces that are of the same order of magnitude as the insertion force during prostate brachytherapy, as reported in [31]). In order to identify the parameters required for simulation of needle in- sertion into this tissue phantom a pilot insertion was executed. An 18 gauge brachytherapy needle (Bard, NJ, USA) was inserted into the tissue phantom with a constant speed, while the tissue phantom was imaged using the sagit- tal (axial) array of a trans-rectal ultrasound (TRUS) probe housed in the hollow cylinder provided in the tissue phantom. The sequence of recorded radio-frequency signals was later used to measure tissue displacement dur- ing needle insertion. The insertion forces and needle positions were also measured. Using the method described in Chapters 2 and 3 the parameters of a three-parameter force distribution model (see Fig.2.10) and the elas- ticity modulus were identified for each layer of the tissue by incorporating the force-displacement data gathered during the experiment. The three- parameter force model consists of a peak force density fp contributing to cutting, shaft force density fs contributing to friction and the width of the peak force density w as shown in Fig.2.10. Using the identified force model and the elastic moduli, a linear FEM based simulator was constructed to simulate the needle insertion into the tissue phantom (for more detail see Chapters 2 and 3). The 3D mesh of the phantom is shown in Fig. 4.6 and the simulation parameters are presented in Table 3.2. The simulation was implemented in quasi-static mode. The bottom side of the tissue phantom was fixed using a wire-mesh and was assumed to be the fixed displacement boundary condition in the simulation. The needle was assumed rigid and 94 4.6. Experimental Validation HollowCylinder Insertion Direction z y 76.2 mm 8 9 m m Targets Figure 4.5: Cross-section of the phantom designed for the experiment, show- ing the inclusion, its connective cylinder, the hollow cylinder and the posi- tions of the targets. 95 4.6. Experimental Validation Figure 4.6: The 3D mesh of the tissue phantom used in the simulations. 96 4.6. Experimental Validation without thickness. A six degrees of freedom position stage was used to orient and insert the needle into the tissue phantom. The stage consists of two perpendicular manually adjustable rotational stages to orient the needle, connected to a manually adjustable linear stage with three degrees of freedom to localize the needle tip on the tissue surface. The sixth degree of freedom is a motorized linear stage to insert the needle to the desired depth with constant speed. The stage is shown in Fig. 4.7. Three targets were located inside the tissue phantom in the central sagit- tal plane (x = 31.8mm). To create the targets, a needle was inserted into the phantom from three points on its side. The intersections of needle tracks and the central sagittal plane were assumed to be the targets. The position of the targets are shown in Fig. 4.5. Since the targets were manually created, they were not located on a straight line. The needle insertion parameters were optimized using the proposed algorithm for these three targets. The needle was inserted twice at the same speed of 5mm/s, once with the initial insertion parameters and another time with the optimized ones. Table 4.7 shows the iterations of the optimization algorithm. The insertion parame- ters and simulated errors are shown for all the iterations, while the actual errors are reported for the two real insertions. The targets were located in the central sagittal plane of the tissue phantom. Considering the symmetry of the tissue phantom the targets will not move laterally (in the x direction) if the needle is inserted in the central sagittal plane (yaw = 0). However, due to the asymmetric mesh used for the simulation, after each insertion of the needle, small lateral displacements were observed resulting in a non-zero yaw angle of insertion for the next iteration. This non-zero yaw angle was assumed as simulation error and was set to zero before advancing to the next iteration. The optimization iterations were terminated, when the error was less than the needle diameter (1.28mm). The brachytherapy needle consists of a stylet sliding inside a cannula. The stylet is used to push the seeds out of the needle after the insertion of the needle into tissue. In order to reduce needle bending due to the beveled-tip of the needle, the needle stylet was fully inserted into the cannula to create a symmetric tip for the needle. No significant needle bending was observed. The TRUS was removed before insertions. After the needle was retracted, the tissue phantom was imaged using a hand-held ultrasound probe. Figure 4.8 shows the position of the targets and the track of the needle. These figures show that insertion with optimal insertion values decreased the targeting error significantly. For the simulations used in this part, the three-parameter force model – shown in Fig. 2.10– was implemented. It was shown that such a model can 97 4.6. Experimental Validation Pitch Yaw x y z Figure 4.7: The 6-DoF stage, used for the positioning and orientation of the needle. 98 4.6. Experimental Validation Table 4.7: Optimization steps for insertion into tissue phantom, simulated and measured errors. Itr. Angles Ins. point Ins. depth Sim. error Real error Num. (Yaw,Pitch) (x, y) (mm) (mm) (mm) 1 (0o,2.5o) (31.8,33.5) 47.7 10.0 13.5 2 (0o,3.9o) (31.8,34.3) 55.7 3.3 – 3 (0o,3.9o) (31.8,34.6) 58.1 1.6 – 4 (0o,4.0o) (31.8,34.9) 58.8 1.2 1.7 99 4.7. Discussion accurately predict the needle insertion forces and the tissue deformations during insertion of a needle into a two-layered PVC prostate phantom simi- lar to the one used in this paper [29, 38]. The parameters of this model were identified in [38] for the same tissue phantom. In order to implement this force distribution model, at each sample time of the simulation, this force profile should be integrated over the part of the needle which is inside the tissue in the deformed configuration. The resultant force is distributed on the mesh nodes that are located on the needle, based on their relative dis- tances from each other and their distances from the needle tip. These forces are applied as a force boundary conditions in the simulation and result in nodal displacements. However, due to tissue deformation, especially at the point of entry of the needle, the length of the part of the needle which is inside the tissue is a function of the displaced nodal positions, which are part of the simulation results. Therefore, an iterative method was employed to implement the three-parameter force model in each time sample [29, 38]. In contrast, the stick-slip method can be implemented in a single step for each sample time of the simulation. However, the three-parameter force models delivers higher accuracy in predicting the insertion force and the tis- sue deformation. In general the force recorded during insertion of a needle into a PVC phantom shows a high slope rise at the start of the penetra- tion followed by a more gradual increase as the insertion continues (see for example measured insertion forces, shown in Fig. 2.12). Such a force can be easily modeled by the three-parameter force model. In contrast, upon implementation of the stick-slip model, anytime a node sticks to the needle tip, a displacement boundary condition is applied, thus the force increases with a gradient which is a function of the elastic properties of the tissue. When the node state changes to slip mode, the force gradient is a function of the friction. Therefore, the insertion force that results from stick-slip model cannot be smooth as the one recorded from a PVC phantom. 4.7 Discussion For each iteration of the optimization algorithm, a needle insertion pro- cess was simulated. When the tissue model was assumed to have a linear elastic behavior, pre-computation and accelerated simulation methods were employed to increase the simulation speed. As an example, the original stiffness matrix for the prostate model was of dimension 1710× 1710 which was inverted and saved in the memory. However, only a part of the inverted stiffness matrix (operating matrix) was used when the tissue was assumed to 100 4.7. Discussion Targets Needle Track (a) Targets Needle Track (b) Figure 4.8: Experimental results. The target positions and the needle track inserted (a) without optimization, and (b) with optimized parameters. The images are brightened and cropped for better visibility. 101 4.7. Discussion be linear elastic. In this work, the operating matrix dimension was 30×30 at the maximum penetration of the needle and maximum 50×50 at the final it- eration to find the displaced locations of the targets. Fast update algorithms from [13, 14] were used to implement the mesh adaptation and update the inverted stiffness matrix. However, the adaptation process is time consum- ing due to the large dimension of the original matrix and huge number of multiplications necessary to update the inverse stiffness matrix. The simu- lation part including the finite element code was programed in Matlab R©and the mean time per simulation was 50 seconds. The mesh adaptation algo- rithm, the application of force and displacement boundary conditions, and the change of coordinate system from global to needle coordinate system are included in each simulation. The time per simulation is the time nec- essary to simulate insertion of a needle to a desired depth. Since variable insertion depths were simulated in this work, the mean time per simulation was reported. With implementation of the code in C++ as in [14] higher speeds can be achieved. For more details about computational issues in needle insertion simulation please refer to [13] and [14]. In the simulation with prostate model the stick-slip interaction model was used. However, in the tissue phantom simulation, the three-parameter force model was imple- mented. As discussed in Section 4.6, implementation of this model needs several iterations for each time sample and therefore, requires significantly more computational time. The original stiffness matrix in this case was of dimension 2361 × 2361. With a non-linear system, it is not possible to use the pre-computation methods. A linear set of equations resulted from Newton-Raphson should be solved at each iteration the coefficient matrix of which – named tangent stiffness matrix - is a function of displacements in all the nodes. Therefore a set of algebraic equations with matrix of di- mensions 1710 × 1710 was solved at each iteration for each sample time of the simulation. The Newton-Raphson method never showed instability in our simulations and always converged to the desired errors within less than 10 iterations, while convergence was achieved in 5 iterations in most sim- ulations. The mean time per simulation was 70 minutes using Matlab R©. The tangent stiffness matrix is sparse therefore, sparse matrix calculation methods can be exploited to increase the simulation speed. In addition, a computationally efficient method for simulation of needle insertion into neo- Hookean tissue models was introduced in [35]. However, linear models are always faster. In Sections 4.5.3 and 4.5.5, the tissue parameters and the boundary conditions are similar while the simulation method is different. Since the tissue was only anchored around the pubic bone and the rotation was high, 102 4.7. Discussion one can expect the largest difference between linear and non-linear models. However, a comparison between Table 4.4 and Table 4.6 shows only minor differences between optimal parameters obtained from linear and non-linear models. Therefore, when very high accuracy is not required, linear models can be used instead of non-linear ones due to their significantly higher speed. The optimization method showed convergence in a few iterations in all the simulation cases. When the simulation model was linear, the error was monotonically decreasing in each iteration. The optimization error at each iteration for the linear simulation model is shown in Fig. 4.9(a) and Fig. 4.9(b). In these figures the prostate mesh was anchored around the pubic bone only. Using the non-linear model, the error monotonically de- creased in all the cases except for sample targets in set 3. The error in optimization for the non-linear simulation model is shown in Fig. 4.9(c). It should be noted that mesh adaptation algorithm applies sudden topological changes to the problem that may lead to sudden changes in the computed displacements of nodes. The relative position of target 3 respect to one of the moving nodes can be the reason for the increase in the error for this target at the 4th iteration as shown in Fig. 4.9(c). The boundary conditions applied in the simulations are not realistic. However, they cover a wide range of possible boundary conditions, from confined situation as in Sec.4.5.1 and Sec.4.5.2 to free to rotate situation as in Sec.4.5.3 to Sec.4.5.5. The optimization algorithm was stable in all cases. Note that it is possible to stabilize such a fixed point algorithm by switching between the proposed iteration and a gradient based descent, as done to stabilize Newton’s method. It should be noted that Lagerburg et al. [1] measured prostate rotations of maximum 13.8o in the coronal plane and rotations between −8.5o to 10.2o in the sagittal plane. The rotation in coronal plane showed strong correlation to the needle insertion point and the method of insertion (with or without locking needles), while the rotation in the sagittal plane did not show significant correlation with the insertion point [1]. The extent of rotations in the simulations reported in this paper shows the most resemblance to the ones reported in [1], when the prostate is only anchored around the pubic bone. However, the simulated rotations in the sagittal plane are strongly correlated with the insertion point of the needle and negative values for rotations in this plane are impossible. Therefore, prior to clinical use of the optimization algorithm, accurate boundary conditions of the prostate within the tissue should be determined as the boundary conditions play a significant role in the outcomes of the optimization. In addition, the elastic properties of the prostate and the surrounding tissue should be identified. Magnetic 103 4.7. Discussion 1 2 3 4 5 0 2 4 6 8 10 12 14 Iteration Number Er ro r ( mm ) Set #1 Set #2 Set #3 Set #4 Set #5 Set #6 Set #7 Set #8 Set #9 (a) 1 2 3 4 5 6 7 0 2 4 6 8 10 12 14 16 Iteration Number Er ro r ( mm ) Set #1 Set #2 Set #3 Set #4 Set #5 Set #6 Set #7 Set #8 Set #9 (b) 1 2 3 4 5 6 0 2 4 6 8 10 12 14 16 Iteration Number Er ro r ( mm ) Set #1 Set #2 Set #3 Set #4 Set #5 Set #6 Set #7 Set #89 (c) Figure 4.9: Optimization error in different iterations. (a): Linear model with E=60kPa, (b): Linear model with E=40kPa, and (c): neo-Hookean nonlinear model. The boundary conditions are the fixed points in connection with the pubic bone. 104 4.8. Conclusions and Future Work Resonance Imaging does produce anatomical detail of the prostate region with very good detail and could be used to produce the required boundary conditions for our simulations. It can also be used to generate a mesh of the prostate and surrounding organs with more details. Alternatively, elastography is also being explored for such modeling as it can produce both the boundary conditions and the tissue elasticity parameters effectively [39, 40]. The proposed optimization algorithm was combined with FEM simula- tion methods in this work. Two different needle-tissue interaction models were employed in these simulations. Since the simulation and optimization parts are independent and the assumptions in the simulation part do not put any restriction on the optimization algorithm, the optimization approach can be combined with any simulation method. In addition, similar methodology can be exploited for other medical interventions involving needle insertion. The needle was assumed to be rigid and straight. In the experiment, the needle stylet was fully extended to create a symmetric tip. Needle bending was negligible even though the needle was flexible. Common brachytherapy needles are flexible and beveled-tip, and may bend due to the bevel force. However, a symmetric-tip needle can be used, or alternatively, the needle can be rotated during insertion as it has been shown to significantly decrease the needle bending due to bevel forces [41]. Under these assumptions, a rigid needle model can be used for insertion optimization. In addition, the opti- mized parameters can also be used as a starting insertion point for flexible needle steering methods such as those presented in [24] or [25]. At first, the initial position of the needle tip and needle orientation can be optimized using the presented algorithm assuming rigidity of the needle. Later the needle can be steered during penetration to further decrease the errors. 4.8 Conclusions and Future Work An approach to the optimization of needle insertion point, heading and depth was proposed for the planning of needle insertions into deformable tissue. The approach was demonstrated in a number of simulations that have explored the effect of boundary conditions, elasticity parameters and nonlinearity on optimization performance. It was shown that boundary conditions play a key role in deciding between linear and non-linear models. If the tissue is confined and is not free to rotate, linear models can be used instead of non-linear models due to their computational efficiency. Unlike in previously published work, our optimization approach employs a fixed- 105 4.8. Conclusions and Future Work point type of iteration instead of descent or search methods and therefore works faster. In fact, more simulation iterations are required in a single step of a descent method than are required to attain good convergence in almost all the simulations studied. The simulation was performed using a finite element based model in quasi-static mode. The algorithm converged to sub-millimeter errors in few iterations in all the cases. The optimization method was also validated in an experiment on a two- layered prostate phantom and showed substantial decrease in the targeting error. In the future, the optimization algorithm will be applied to a more re- alistic model of the prostate based on segmented patient images that in- clude boundary conditions such as the pubic bone and the bladder. We will include more accurate needle models accounting for needle flexibility and lateral bevel forces. 106 References [1] V. Lagerburg, M. A. Moerland, J. J. Lagendijk, and J. J. Batter- mann, “Measurement of prostate rotation during insertion of needles for brachytherapy,” Radiotherapy and Oncology, vol. 77, no. 3, pp. 318– 323, 2005. [2] R. Taschereau, J. Pouliot, J. Roy, and D. Tremblay, “Seed misplacement and stabilizing needles in transperineal permanent prostate implants,” Radiotherapy and Oncology, vol. 55, no. 1, pp. 59–63, 2000. [3] O. Goksel, S. E. Salcudean, S. P. DiMaio, R. Rohling, and J. Morris, “3D needle-tissue interaction simulation for prostate brachytherapy,” in Proc. Int. Conf. Medical Image Computing and Computer Assisted Intervention, 2005, pp. 827–834. [4] V. Lagerburg, M. A. Moerland, M. van Vulpen, and J. J. Lagendijk, “A new robotic needle insertion method to minimize attendant prostate motion,” Radiotherapy and Oncology, vol. 80, pp. 73–77, 2006. [5] G. Fichtinger, E. C. Burdette, A. Tanacs, A. Patriciu, D. Mazilu, L. L. Whitcomb, and D. Stoianovici, “Robotically assisted prostate brachytherapy with transrectal ultrasound guidance – phantom exper- iment,” Brachytherapy, vol. 5, no. 1, pp. 14–26, 2006. [6] S. P. DiMaio, S. Pieper, K. Chinzei, N. Hata, S. J. Haker, D. F. Kacher, G. Fichtinger, C. M. Tempany, and R. Kikinis, “Robot-assisted needle placement in open MRI: System architecture, integration and valida- tion,” Computer Aided Surgery, vol. 12, no. 1, pp. 15–24, 2007. [7] H. Bassan, T. Hayes, R. Patel, and M. Moallem, “A novel manipulator for 3D ultrasound guided percutaneous needle insertion,” in Proc. IEEE Int. Conf. Robotics and Automation, 2007, pp. 617–622. [8] Y. D. Zhang, T. K. Podder, W. S. Ng, J. Sherman, V. Misic, D. Fuller, E. M. Messing, D. J. Rubens, J. G. Strang, R. Brasacchio, and Y. Yu, “Semi-automated needling and seed delivery device for prostate 107 Chapter 4. References brachytherapy,” in Proc. IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 2006, pp. 1279–1284. [9] S. E. Salcudean, T. D. Prananta, W. J. Morris, and I. Spadinger, “A robotic needle guide for prostate brachytherapy,” in Proc. IEEE Int. Conf. Robotics and Automation, 2008, pp. 2975 – 2981. [10] R. Alterovitz, K. Goldberg, J. Pouliot, R. Taschereau, and I. C. Hsu, “Sensorless planning for medical needle insertion procedures,” in Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems, vol. 3, 2003, pp. 3337–3343. [11] R. Alterovitz, K. Goldberg, and A. Okamura, “Planning for steerable bevel-tip needle insertion through 2D soft tissue with obstacles,” in Proc. IEEE Int. Conf. Robotics and Automation, 2005, pp. 1640 – 1645. [12] N. Abolhassani, R. Patel, and M. Moallem, “Needle insertion into soft tissue: A survey,” Medical Engineering and Physics, vol. 29, no. 4, pp. 413–431, 2007. [13] S. P. DiMaio and S. E. Salcudean, “Needle insertion modeling and simulation,” IEEE Trans. Robotics and Automation: Special Issue on Medical Robotics, vol. 19, no. 5, pp. 864 – 875, 2003. [14] O. Goksel, S. E. Salcudean, and S. P. DiMaio, “3D simulation of needle- tissue interaction with application to prostate brachytherapy,” Com- puter Aided Surgery, vol. 11, no. 6, pp. 279–288, 2006. [15] E. Dehghan, and S. E. Salcudean, “Comparison of Linear and Non- linear Models in 2D Needle Insertion Simulation,” in Proc. Workshop on Computational Biomechanics for Medicine in Conjunction with MIC- CAI’06, pp. 117–124, 2006. [16] S. Misra, K. J. Macura, K. T. Ramesh, and A. M. Okamura, “The importance of organ geometry and boundary constraints for planning of medical interventions,” Medical Engineering & Physics, vol. 31, no. 2, pp. 195–206, 2009. [17] R. Alterovitz, J. Pouliot, R. Taschereau, I. C. Hsu, and K. Gold- berg, “Simulating needle insertion and radioactive seed implantation for prostate brachytherapy,” inMedicine Meets Virtual Reality 11. IOS Press, 2003, pp. 19–25. 108 Chapter 4. References [18] E. Dehghan and S. E. Salcudean, “Needle insertion point and orienta- tion optimization in non-linear tissue with application to brachyther- apy,” in Proc. IEEE Int. Conf. Robotics and Automation, 2007, pp. 2267–2272. [19] S. Nag, J. P. Ciezki, R. Cormak, S. Doggett, K. Dewyngaert, G. K. Edmundson, R. G. Stock, N. N. Stone, Y. Yan, and M. J. Zelefsky, “In- traoperative planning and evaluation of permanent prostate brachyther- apy: Report of the American brachytherapy society,” Int. J. Radiation Oncology, Biology, Physics, vol. 51, no. 5, pp. 1422–1430, 2001. [20] R. J. Webster, J. S. Kim, N. J. Cowan, G. S. Chirikjian, and A. M. Okamura, “Nonholonomic modeling of needle steering,” Int. J. Robotics Research, vol. 25, no. 5-6, pp. 509 – 525, 2006. [21] R. Alterovitz, A. Lim, K. Goldberg, G. S. Chirikjian, and A. M. Oka- mura, “Steering flexible needles under Markov motion uncertainty,” in Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems, 2005, pp. 120 – 125. [22] W. Park, J. S. Kim, Y. Zhou, N. J. Cowan, A. M. Okamura, and G. S. Chirikjian, “Diffusion-based motion planning for a nonholonomic flexi- ble needle model,” in Proc. IEEE Int.l Conf. Robotics and Automation, 2005, pp. 4600–4605. [23] E. Dehghan, O. Goksel, and S. E. Salcudean, “A comparison of needle bending models,” in Proc. Medical Image Computing and Computer Assisted Intervention, 2006, pp. 305–312. [24] S. P. DiMaio and S. E. Salcudean, “Needle steering and motion planning in soft tissues,” IEEE Trans. Biomedical Engineering, vol. 52, no. 6, pp. 965 – 974, 2005. [25] D. Glozman and M. Shoham, “Flexible needle steering and optimal trajectory planning for percutaneous therapies,” in Proc. Medical Im- age Computing and Computer Assisted Intervention. Springer-Verlag, 2004, pp. 137–144. [26] N. Abolhassani, R. Patel, and M. Moallem, “Trajectory generation for robotic needle insertion in soft tissue,” in Proc. IEEE Engineering in Medicine and Biology Society Int. Conf., 2004, pp. 2730–2733. 109 Chapter 4. References [27] M. Heverly, P. Dupont, and J. Triedman, “Trajectory optimization for dynamic needle insertion,” in Proc. of IEEE Int. Conf. Robotics and Automation, 2005, pp. 1646 – 1651. [28] A. Okamura, C. Simone, and M. O’Leary, “Force modeling for needle insertion into soft tissue,” IEEE Tran. Biomedical Engineering, vol. 51, no. 10, pp. 1707 – 1716, 2004. [29] E. Dehghan, X. Wen, R. Zahiri-Azar, M. Marchal, and S. E. Salcudean, “Needle-tissue interaction modeling using ultrasound-based motion es- timation: Phantom study,” Computer Aided Surgery, vol. 13, no. 5, pp. 265–280, 2008. [30] H. Kataoka, T. Washio, K. Chinzei, K. Mizuhara, C. Simone, and A. Okamura, “Measurement of tip and friction force acting on a needle during penetration,” in Proc. Medical Image Computing and Computer Assisted Intervention, vol. 2488, 2002, pp. 216–223. [31] T. Podder, D. Clark, J. Sherman, D. Fuller, E. Messing, D. Rubens, J. Strang, R. Brasacchio, L. Liao, W.-S. Ng, and Y. Yu, “In vivo motion and force measurement of surgical needle intervention during prostate brachytherapy,” Medical Physics, vol. 33, no. 8, pp. 2915–2922, 2006. [32] S. P. DiMaio and S. E. Salcudean, “Needle steering and model-based trajectory planning,” in Proc. Medical Image Computing and Computer Assisted Intervention, 2003, pp. 33–40. [33] J. Ophir, S. Alam, B. Garra, F. Kallel, E. Konofagou, T. Krouskop, and T. Varghese, “Elastography: ultrasonic estimation and imaging of the elastic properties of tissues,” Engineering in Medicine, vol. 213, no. 3, pp. 203–233, 1999. [34] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method: Solid Mechanics, 5th ed. Butterworht-Heinemann, 2000, vol. 2. [35] H. W. Nienhuys and A. F. van der Stappen, “A computational tech- nique for interactive needle insertions in 3D nonlinear material,” in Proc. IEEE Int. Conf. Robotics and Automation, vol. 2, 2004, pp. 2061 – 2067. [36] X. Wu, M. S. Downes, T. Goktekin, and F. Tendick, “Adaptive non- linear finite elements for deformable body simulation using dynamic progressive meshes,” Computer Graphics Forum, vol. 20, no. 1, pp. 349 – 358, 2001. 110 Chapter 4. References [37] S. K. Kyriacou and C. Davatzikos, “A biomechanical model of soft tissue deformation, with applications to non-rigid registration of brain images with tumor pathology,” in Proc. Medical Image Computing and Computer Assisted Intervention, 1998, pp. 531–538. [38] E. Dehghan and S. E. Salcudean, “Needle insertion study using ultrasound-based 2D motion tracking,” in Proc. Medical Image Com- puting and Computer Assisted Intervention, 2008, pp. 660 – 667. [39] J. Ophir, S. K. Alam, B. S. Garra, F. Kallel, E. E. Konofagou, T. Krouskop, C. R. B. Merritt, R. Righetti, R. Souchon, S. Srinivasan, and T. Varghese, “Elastography: Imaging the elastic properties of soft tissues with ultrasound,” Medical Ultrasounics, vol. 29, no. 4, pp. 155– 171, 2002. [40] S. E. Salcudean, D. French, S. Bachmann, R. Zahiri-Azar, X. Wen, and W. J. Morris, “Viscoelasticity modeling of the prostate region using vibro-elastography,” in Proc. Medical Image Computing and Computer Assisted Intervention, 2006, pp. 389 – 396. [41] N. Abolhassani, R. V. Patel, and F. Ayazi, “Minimization of nee- dle deflection in robot-assisted percutaneous therapy,” Int. J. Medical Robotics and Computer Assisted Surgery, vol. 3, no. 2, pp. 140–148, 2007. 111 Chapter 5 Modeling and Simulation of Flexible Needles 4 5.1 Introduction Percutaneous insertion of long and flexible needles into soft tissue is in- volved in many clinical and therapeutic procedures such as biopsy and prostate brachytherapy. As a result of needle-tissue interaction, the needle base movements and, for needles with a beveled tip, the asymmetric cutting force on the tip, the tissue deforms and the needle bends (see Fig. 5.1). The targeting procedure is complicated by the bending of the needle shaft, tar- get displacement due to tissue deformation, and insufficient visual feedback from medical imaging modalities. Accurate needle insertion requires signif- icant skill and training of the performing physician. Modeling, simulation, and path planning of needle insertion are emerging fields of research aimed at providing the physicians with training devices and accurate pre-surgery plans. Needle insertion simulators usually include a soft tissue model, a flexible needle model and a needle-tissue interaction model. The simulation of needle bending and tissue deformation together in one combined model is generally not feasible for two reasons: (i) a fast solution for a large tissue mesh requires exploiting simplifications in deformation equations, which are not suitable for estimating the needle behavior and (ii) the interaction surface of the needle and tissue changes during insertion, which many techniques like the FEM cannot inherently accommodate. Therefore, generally two separate models for the tissue and the needle are employed [1, 2, 3] with a third needle-tissue contact model governing their interaction. Deformable tissue models have been studied extensively in simulating tis- sue deformation during surgery and needle insertion [1, 2, 3, 4, 5]. Models based on the Finite Element Method (FEM) are the most common tech- 4A version of this chapter is submitted as: O. Goksel, E. Dehghan, and S. E. Salcudean, “Modeling and Simulation of Flexible Needles”, 2009 112 5.1. Introduction F1 F2 F3 Figure 5.1: The forces on the needle shaft caused by tissue deformation and needle base manipulation. 113 5.1. Introduction niques employed for tissue deformation simulation during needle insertion [1, 2, 3, 6]. Mass-spring models [7] have also been used for this purpose. The interaction of needles with such deformable tissue models has also been studied widely [2, 5, 9, 10]. The aim of this paper is to compare different models for simulating the bending of flexible needles. In general, medical needles can be categorized into three major groups: rigid needles, highly flexible needles, and (moder- ately) flexible needles. Rigid needles keep a straight posture regardless of the forces applied on them during insertion. Due to their simplicity, rigid needle models have been used when the needle physical properties and the insertion procedure lead to negligible bending [1, 2, 6]. On the opposite end of the flexibility range are the highly flexible needles. These needles are as- sumed to bend in the direction of their tip bevel with a constant curvature, without applying considerable force to the tissue in the lateral direction. These needles were modeled as a non-holonomic system [4] and were used for needle insertion simulations and planning in [11, 12, 13, 14]. Some needles, such as brachytherapy needles, cannot be categorized as either rigid or highly flexible. They are not rigid, since their deflection dur- ing procedures can be significant. They are not highly flexible either, since a considerable lateral force is necessary to bend them. Several groups have modeled this type of needles. DiMaio and Salcudean simulated the nee- dle as an elastic material using the FEM with geometric non-linearity and 3-node triangular elements [15]. This method was extended to 3D using 4-node tetrahedral elements by Goksel et al. [3]. These simulation methods were validated in phantom studies. Glozman and Shoham used linear 2D beam elements to simulate the needle bending for needle steering [7]. Linear beam theory was also used to introduce a needle steering model with online parameter estimator [17], to estimate the needle tip deflection during inser- tion due to tip bevel [18], and to identify the shaft force profile due to the bevel [19]. Models based on linear beam theory are relatively simple and fast. However, they are not rotationally invariant and cannot preserve the needle length during large deformations/deflections. Many needles, such as brachytherapy needles, consist of a stylet sliding inside a hollow cannula. Physical modeling of this combination without any simplifications is indeed very complicated as it requires separate models for the stylet and the cannula and an interaction model to simulate their interface. In this paper, flexible needles are approximated as solid bars and, accordingly, models applicable to solid bars are examined to simulate their flexion. Three different models were used to simulate needle bending. The first 114 5.2. Methods two models are based on the FEM and were chosen due to their frequent use in the literature, while the third model is an angular springs model. The first model uses 4-node tetrahedral elements, where non-linear geometry is accommodated to simulate large deformations. The second model also ac- commodates non-linear geometry and uses Euler-Bernoulli non-linear beam elements. In this work, non-linear beam elements were preferred over more common linear ones in the literature due to their superiority in modeling large deformations. The third model is novel and utilizes angular springs for the quasi-static simulation of needle bending. In the literature, angular springs have been used to model hair and cloth [20] since their dynamics behavior can be computed easily. They have been also incorporated in 3D mass-spring models to simulate vibration and deformation [21]. The Young’s modulus is the parameter that describes needle bending in the first two models. Similarly, the third model is identified by its spring constant. In this paper, all models are devised for needle bending in 3D. The parameter of each model is identified for a brachytherapy needle through experiments, where several lateral forces were applied to the needle tip and the shaft deflection was recorded. The parameters defining each model were identified and the models were studied for their accuracy in simulating the actual needle deflection observed during experiments. For the single force applied at the needle tip during the experiments, the needle deformation is entirely planar. Therefore, the parameter identi- fication and the model validations were performed in 2D. Nevertheless, the same identified parameters also describe a needle in 3D, since the shafts of most medical needles are built with axial symmetry resulting in symmetric deflection for the same force rotated around their long axis. The following section derives the models in 3D. Next, the experimental method to validate the models is described in section 5.3. The results and a discussion follow in Sections 5.4 and 5.5, respectively. Finally, conclusions are presented in the last section. 5.2 Methods 5.2.1 Finite Element Method using Tetrahedral Elements The Finite Element method is a powerful tool for approximating a solution to the continuum mechanics equations. In this method, an entire body Ω is divided into several discrete elements Ωe. Then, the constitutive equations are approximated over each element and combined to give an approximation to the global solution. Various types of elements can be chosen depending on 115 5.2. Methods the nature of the problem. 4-node tetrahedral (TET4) or 3-node triangular (TRI3) elements are the simplest elements to use in 3D or 2D deformation analysis, respectively. Considering a deformable body, the 3D coordinates of a point in the un- deformed configuration can be written as [x1 x2 x3] in the reference frame. The coordinates of that point in the deformed configuration can then be expressed using the displacements [u1 u2 u3] added onto the undeformed co- ordinates. Using this definition, the Green strain and the second Piola- Kirchhoff stress tensors are given as follows: ij = 1 2 ( ∂ui ∂xj + ∂uj ∂xi + 3∑ k=1 ∂uk ∂xi ∂uk ∂xj ) , (5.1) σij = 3∑ k=1 3∑ l=1 Cijkl kl, (5.2) where i, j ∈ {1, 2, 3} and Cijkl is the material moduli tensor for linear elastic materials which is a function of Young’s modulus and Poisson’s ratio. When the deformation is small, the second order terms in (5.1) can be neglected leading to a linear relation between strain and displacement. This linear relation leads to a set of linear algebraic equations with a con- stant coefficient matrix. Such a system can be solved rapidly using off-line computation and condensation [2, 22]. However, the assumption of small deformation is not valid for needle bending, where the deformation is large with respect to the needle diameter [15]. Thus, non-linear geometry as in (5.1) must be included in the model. This leads to a non-linear relation between the nodal forces and displacements: f = K̂(u)u (5.3) where u and f are the vectors of nodal displacements and forces, respectively, and K̂ is the stiffness matrix. This non-linear equation can accommodate the axial displacements during lateral deformation and, therefore, preserves the needle length. However, the use of condensation is not possible in the non-linear case and this leads to a time-consuming solution. The Newton-Raphson method is an effective iterative technique for solv- ing the non-linear algebraic equations in (5.3) [23]. In this method, a tangent stiffness matrix KT is computed at each iteration and used to find a descent 116 5.2. Methods direction for the solution using the iteration scheme below: KT∆u t = f − K̂(ut)ut (5.4) ut+1 = ut +∆ut (5.5) to update the nodal displacements at each iteration t. In our simulations, the Newton-Raphson method converged to a solution for (5.3) in a few iterations. Unfortunately, computing the tangent stiffness matrix and solving (5.4) at each iteration are both computationally intensive operations. 5.2.2 Finite Element Method using Non-Linear Beam Elements Beams elements were designed and introduced in the FEM literature in order to model the deformation of bars, rods, and other beam-like structures [24]. These elements are expected to be better-suited for needle deformations than triangular/tetrahedral elements. Thus, a beam-element model was also implemented to compare its performance in modeling needles. The Euler-Bernoulli beam element is a common element formulation to model thin bars and therefore is used as a needle model in this paper. Euler-Bernoulli beam theory suggests that each plane perpendicular to the beam axis prior to any deformation remains as a plane that is still perpendicular to the beam axis after the deformation. Under this condition the displacements of each material point can be written as a function of the displacement along the neutral axis of the beam as follows: u(x, y, z) = u0(x)− y ∂v0(x) ∂x − z ∂w0(x) ∂x , (5.6) v(x, y, z) = v0(x)− zα0(x), (5.7) w(x, y, z) = w0(x) + yα0(x), (5.8) where u0, v0, and w0 are the neutral axis displacements along the x, y, and z axes, respectively, and α0 is the twist around the beam axis. It is assumed that the beam axis is lying along the x axis as shown in Fig. 5.2. Please note that, in order to conform to the beam-element literature, the node positions [x1 x2 x3] and the node displacements [u1 u2 u3] in the previous subsection are referred here as [x y z] and [u v w], respectively. Assuming that the beam twist is small while the lateral deformations are 117 5.2. Methods moderately large, strains in the beam element can be approximated as [25]: xx ≈ du0 dx − y d2v0 dx2 − z d2w0 dx2 + 1 2 (( dv0 dx )2 + ( dw0 dx )2) , (5.9) xy ≈ − 1 2 z dα0 dx , (5.10) xz ≈ 1 2 y dα0 dx . (5.11) In the beam element, the displacement of the material nodes on the beam axis are interpolated from the nodal variables of the two end-points (nodes) of the beam. The nodal variables at node p are the axial displacement up, the twist αp, the transverse displacements vp and wp (along the y and z axes, respectively), and their derivatives φp = −dv/dx and ψp = −dw/dx. Euler-Bernoulli beam theory employs cubic interpolation functions for lat- eral displacements, and linear interpolation functions for axial displacements and twist angle. Therefore, within each element between nodes p and p+1, the variables are interpolated as follows: u0(x) = N1(x)up +N2(x)up+1 (5.12) α0(x) = N1(x)αp +N2(x)αp+1 (5.13) v0(x) = M1(x)vp +M2(x)φp +M3(x)vp+1 +M4(x)φp+1 (5.14) w0(x) = M1(x)wp +M2(x)ψp +M3(x)wp+1 +M4(x)ψp+1 (5.15) where Ni(x) and Mi(x) denote the linear and cubic interpolation functions, respectively. An element equilibrium equation relating loads on the nodes to the nodal variables is written as: K(u)u = f, (5.16) where u is the vector of nodal variables described above and f is the vector of corresponding nodal forces and torques. Please see Appendix C for detailed information on the calculation of the stiffness matrix and its integration. To solve (5.16), Picard’s method [24] is used in this paper. This method employs the following iterative steps: K(ut)ut+1 = f (5.17) where t is the iteration number. 118 5.2. Methods z x v0 y u0 w0 dx dw0 Figure 5.2: Euler-Bernoulli beam element under deformation. 119 5.2. Methods y z xèP P á P - + Figure 5.3: Angles of bending and twisting between two needle segments. 120 5.2. Methods 5.2.3 Angular Springs Model In this section, the needle is modeled as a discrete structure composed of n connected rigid rods. A joint P between two consecutive segments of such a model is shown in Fig. 5.3, where the segment PP+ is bent and twisted relative to the segment before it, P−P , under external loads. The magnitude of bending and twisting can be described by two angles: the bending (zenith) angle θ from the x axis to the segment PP+, and the twisting angle α about the segment PP+ due to torsion around its axis. Naturally, a deformed (bent and twisted) needle structure applies inter- nal reaction (“straightening”) torques/forces to un-bend and un-twist itself. These are modeled as torques applied to link PP+ at joint P to restore its rest configuration. The effect of these torques can be visualized as two ro- tational springs —one un-bending and the other un-twisting the segment— shown in Fig. 5.3. It is evident that the magnitudes of such restoring torques are functions of the deviations (angles θ and α) from the rest configuration. Let us define the magnitudes of these reaction torques, τb for bending and τt for twisting, as functions of the angles as follows: τb = hb(θ) , (5.18) τt = ht(α) . (5.19) Although the actual characteristics of hb and ht may be complex and possibly non-linear, for the elastic range of deformations, a linear relation- ship is expected. Indeed, when an infinitesimal section of a cantilever shaft is analyzed for its bending (in Appendix D) and its twisting (in Appendix E), linear angle-torque dependencies are observed for both of them. These insights motivate us to devise our model composed of linearly-acting springs as follows: τb = kb θ , (5.20) τt = kt α . (5.21) Then, the energy stored in these springs can be written as: V (θ, α) = 1 2 kbθ 2 + 1 2 ktα 2 . (5.22) As a model of the connections between each pair of segments, a spherical wrist —the rest configuration of which is shown in Fig. 5.4(a)— is chosen in this paper. These three wrist joints, seen in a bent/twisted configuration 121 5.2. Methods in Fig. 5.4(b), define the orientation of the segment to its right with respect to the segment to its left. Note that, whereas the third joint α defines the twisting angle independently, the bending angle θ is combination of the joint angles γ and β. It can be geometrically shown that: cos θ = cos γ cos β . (5.23) In order to derive compliance models for the joints γ and β, let us rewrite the energy as a function of these joint variables, e.g. V (γ, β, α). For suffi- ciently small angles, θ ≈ sin θ. This greatly simplifies the following deriva- tion and the solution of final system of equations. Thus, replacing θ2 = sin2 θ and (5.23) in (5.22): V (θ, α) = 1 2 kb(1− cos 2 θ) + 1 2 ktα 2 , (5.24) V (γ, β, α) = 1 2 kb(1− cos 2 γ cos2 β) + 1 2 ktα 2 . (5.25) Considering the invariance of the Langrangian at equilibrium, the partial derivatives of joint angles are equal to joint torques, which are the gradients of the potential with respect to the joint angle coordinates, as follows: τγ = ∂V ∂γ = kb sin γ cos γ cos 2 β , (5.26) τβ = ∂V ∂β = kb sin β cos β cos 2 γ , (5.27) τα = ∂V ∂α = ktα . (5.28) The small angle approximation above also applies to γ and β, i.e. sin γ ≈ γ, sin β ≈ β, cos γ ≈ 1, and cos β ≈ 1, thus yielding: τγ ≈ kbγ , (5.29) τβ ≈ kbβ . (5.30) Consider a 3D bent needle posture in equilibrium with force loads f i ap- plied at its joints as in Fig. 5.5. Using our spherical wrist model, a Jacobian matrix J can be written for this 3n-joint robotic arm [26], where the joint torques τ = [τ ′1 . . . τ ′ n] ′ relate to forces f as follows: τ = J ′f . (5.31) 122 5.2. Methods ° á¯ PP - + y z x P (a) ° á¯ P P - +y z x P (b) Figure 5.4: Three joints of our spherical wrist model in (a) initial and (b) bent configurations. 123 5.2. Methods fnfifi-1 fn-1f2f1 ( )° ¯, ,á1 11 ( )° ¯, ,á2 22 ( )° ¯, ,á3 33 ( )° ¯, ,ái-1 i-1i-1 ( )° ¯, ,ái ii ( )° ¯, ,án nn Figure 5.5: The entire needle in equilibrium state constrained at the base and under nodal forces at its joints. 124 5.2. Methods From (5.28–5.30), the 3n vector of joint torques τ can be written as a combination of 3n joint angles Φ = [γ1 β1 α1 · · · γn βn αn ] ′ and the corre- sponding joint stiffnesses as follows: τ =KsΦ (5.32) whereKs is a diagonal matrix formed by joint stiffnesses {kb,1 kb,1 kt,1 · · · kb,n kb,n kt,n}. Note that since the Jacobian J is configuration dependent, (5.31) is indeed a nonlinear system of equations as follows: KsΦ = J ′(Φ)f . (5.33) In this paper, this nonlinear system is solved iteratively using Picard’s method: Φt+1 =Ks −1 J ′(Φt)f (5.34) which consists of two phases: (i) Given the needle configuration, the Ja- cobian is composed, and (ii) the torques found by J ′f are divided by cor- responding joint stiffnesses to find the joint angles and hence the needle posture at the next iteration. 5.2.4 Convergence Criterion In order to have a fair comparison between the three needle models, a com- mon convergence criterion was devised to be used in the simulations as follows: Error = ‖ut − ut−1‖ ‖ut‖ < γ (5.35) where ut is the vector of nodal displacements at the iteration t and γ is the tolerance. In this paper, a tolerance of 10−4 is used. When the convergence criterion in (5.35) is met, the residual of error for all the models were checked to ensure convergence. For the proposed angular-springs model in (5.33), the torque residual is: e = ∥∥J ′(Φt)f −KsΦt∥∥ . (5.36) For all the simulations in this paper, the terminal residuals of the three models were sufficiently small indicating model convergence to equilibrium. Figure 5.6 shows the convergence during a sample planar deflection with a tip force applied in −y direction. The sum of angles γ is plotted in 125 5.2. Methods 1 2 3 4 5 6 0.39 0.395 0.4 0.405 0.41 0.415 0.42 iteration # Sum of bending angles [rad] (a) 1 2 3 4 5 6 10−2 10−1 100 101 102 iteration # Residual torque [mN] (b) Figure 5.6: For a sample bending simulation, (a) shows the sum of joint angles γ, which is the deflection of the needle tip, and (b) shows the residual torque at each iteration. 126 5.3. Experiment Fig. 5.6(a). Note that for this planar bending, this sum is equal to the deflection angle of the needle tip with respect to the needle base frame lying along x axis. The torque residual during this simulation is also plotted in Fig. 5.6(b). 5.3 Experiment In order to show the feasibility of the models and compare their accu- racy, the following experiments were conducted using an 18 gauge 20 cm Bard R©BrachyStar R©needle (C.R.Bard, Inc. , Covington, GA), that is used in prostate brachytherapy seed-implant procedures. In these experiments, the needle was bent under several known forces and its bent shaft form is recorded for evaluating our model simulations. During the experiments, the needle was clamped at its base while its shaft lying horizontal to the ground. Due to the clamping mechanism, the effective bending shaft-length of the needle was reduced to 18.7 cm. Sub- sequently, several vertical forces were applied to its tip using combinations of scaling weights. The weights were placed on a stage hanging on a tiny hook bent at the tip in order to achieve a perpendicular force to the needle base at all times. The stage setup for hanging weighs 1.3 g and the weights applied were 5, 10, 15, 20, 25, 30, and 40 g. For each force applied, the needle was imaged as in Fig. 5.7 over a white background using a digital camera, the shutter of which was controlled via computer. Ordinary brachytherapy needles are marked with black stripes at every 1 cm. These markers were located on each image to determine the bent configuration of the needle. Needle-tip lateral deflections from the nominal axis, measured using this procedure, are shown in Table 5.1. To minimize the effect of lens warping on the data, the needle was cen- tered close to the image center and a far camera focal distance was used. By imaging straight lines, this effect was found to be negligible in the given setup. In a calibration step using a ruler, the image resolution on the needle plane was measured to be better than 9pixels/mm. The observations of shaft bending in these experiments were next used to (i) identify the model parameters for the given brachytherapy needle, and (ii) evaluate the accuracy of each model by comparing the simulated bent needles to the experimental data. 127 5.3. Experiment Figure 5.7: Needle shaft clamped at its base a vertical force is applied at its tip. 128 5.4. Results Table 5.1: Tip deflection for various lateral tip forces. Force [mN] 63 113 163 213 263 313 413 Deflection [mm] 8.3 14.5 20.9 27.0 33.0 39.0 50.2 129 5.4. Results 5.4 Results In the experiments, the needle deformations were all planar. Therefore, the model equations were reduced to their 2D equivalents. The equations for a 2D non-linear beam element are obtained by simply removing the rows and columns of the stiffness matrix in (5.16) that correspond to v, φ and α in the vector of nodal variables (reducing K to its parts K11, K13, K31, and K33 in Appendix C). The angular spring model can be used in 2D simply by neglecting the torsion. The FEM with tetrahedral elements reduces to triangular elements in 2D. A 2D plane-stress analysis was then performed by also taking the out-of-plane thickness of the elements into account. The Poisson’s ratio was taken 0.3, which is the Poisson’s ratio of steel. Therefore, the model parameters to be identified are the Young’s modulus E in the two FEM-based models and the bending spring constant kb in the angular spring model. Equal-length segments were chosen so that the needle bending is described by a single parameter in each model. 5.4.1 Parameter Identification For parameter identification, two independent goals were specified: (i) a best fit to the experimental tip position and (ii) a best fit of the entire nee- dle shaft. Accordingly, two cost functions relating an actual bending to a simulated one were defined as seen in Fig. 5.8: (i) a lateral tip position error, that is the vertical difference between the simulated and the measured tip position, and (ii) the area lying between the simulated and the observed nee- dle shafts. These goals, namely the tip fit and the area fit, were formulated as two separate optimization problems. Subsequently, the model parameters were found for each of them, independently, using the Nelder-Mead simplex search method. In the first phase the model parameters were identified for different num- bers of elements using each pair of force-displacement data. The effect of the number of elements on the identified parameters and the accuracy of the simulation were studied using models with 10 and 20 segments. In the triangular FEM model, 11×2 and 21×2 node meshes were employed corre- sponding to 10 and 20 segments, respectively. The 11×2 triangular mesh employed is seen in Fig. 5.9, where the shaft thickness was set to that of a brachytherapy needle. Figure 5.10 shows the range of the identified Young’s moduli and the identified angular spring constants over the range of applied forces for dif- ferent numbers of segments. In this figure, the identified parameters for both 130 5.4. Results AreaError Tip Error Figure 5.8: Tip error and area error between the simulated model and the observed needle postures. 131 5.4. Results 186.7 mm 1.28 mma segment Figure 5.9: Needle discretized by triangular elements on a 11×2 mesh. 132 5.4. Results Tip fit Area fit 4.75 4.8 4.85 x 108 E (Pa) (a) Tip fit Area fit 1.84 1.86 1.88 1.9 x 109 E (Pa) (b) Tip fit Area fit 5 5.1 5.2 5.3 5.4 x 1011 E (Pa) (c) Tip fit Area fit 1020 1030 1040 K ( Nm)µspr (d) Tip fit Area fit 1900 1910 1920 1930 K ( Nm)µspr (e) Tip fit Area fit 4520 4560 4600 K ( Nm)µspr (f) Figure 5.10: Distribution of identified Young’s moduli (a-c) and spring con- stants (d-f) by minimizing tip and area errors for each bending experiment. Boxes extend from the lower quartile to the upper quartile value with the median marked. Whiskers show the full extent of the data. Models shown are (a) 10 segment triangular, (b) 20 segment triangular, (c) 10/20 segment beam (having equal values), (d) 10 segment spring, (e) 20 segment spring, and (f) 50 segment spring. 133 5.4. Results optimization goals are presented as box-plots displaying the distribution of parameters for varying tip loads. Figures 5.10(a) and 5.10(b) present the identified Young’s moduli for 10 and 20 triangular elements, respectively. The Young’s moduli identified for both 10 and 20 segments are the same and therefore they are presented using a single plot in Fig. 5.10(c). Similarly, Figs. 5.10(d), 5.10(e), and 5.10(f) present the identified bending spring constants for 10, 20, and 50 segments, respectively. As seen in the plots, Young’s modulus for triangular elements and spring constants for angular springs are highly dependent on the number of segments. This is due to the nature of these models where adding segments increases the shaft flexibility since the deflections of each segment add up. As opposed to these two models above, the beam model is not sensitive to the number of segments. Note that, for a given number of segments and optimization goal (e.g., the tip fit), the variation of the identified parameters with the changing force is relatively small in all three models. This slight variation is observed to be mainly due to experimental errors such as the resolution of the imaging system and the localization of the black stripes on the shaft. 5.4.2 Performance in Modeling Needle Deflection During a needle insertion simulation, the forces are not known a priori. Thus, the goal of a needle deformation model is to successfully predict nee- dle deflections for a wide range of possible loads using only a single set of fixed parameters. In this paper, we use a single identified parameter for each model. These parameters, identified in the previous section, vary slightly with the applied force. Consequently, the effect of using a single fixed pa- rameter for all applied loads is studied next. In total, six parameters have been identified for three different models with 10 and 20 segments each. As seen in Fig. 5.10, the parameters identified using the tip fit and the area fit are close, thus either of these error measures could be used to determine a fixed parameter for each model. The mean of the parameter values identified using the tip fit method was chosen as the fixed parameter in simulating each model for evaluation. Each fixed parameter was then used to simulate the needle deflection for all the seven different tip loads and these simulated deflections were compared with the experimental results. Fig. 5.11 shows these simulated deflections for the three models with 20 segments along with the experimen- tally observed needle configurations. The tip and area errors for 10 and 20 segments using fixed parameters are presented in Fig. 5.12. 134 5.4. Results 0 20 40 60 80 100 120 140 160 180 −50 −40 −30 −20 −10 0 X(mm) Y( mm ) (a) The triangular element model with a 21× 2 mesh. 0 20 40 60 80 100 120 140 160 180 −50 −40 −30 −20 −10 0 X(mm) Y( mm ) (b) The non-linear beam model with 20 elements. 0 20 40 60 80 100 120 140 160 180 −50 −40 −30 −20 −10 0 X(mm) Y( mm ) (c) The angular springs model with 20 segments. Figure 5.11: The experimental deformed needle posture (circles) and the simulated needle (thicker solid lines) using the mean value of identified pa- rameters. 135 5.5. Discussion 65 113 163 213 263 313 413 0 0.5 1 1.5 2 2.5 Force(mN) Ti p Er ro r(m m) Non−linear Beam, 10 segments Non−linear Beam, 20 segments Angular Springs, 10 segments Angular Springs, 20 segments Triangular Element, 10 segments Triangular Element, 20 segments 63 113 163 213 263 313 413 0 20 40 60 80 100 120 140 Force(mN) Ar ea E rro r(m m2 ) Non−linear Beam, 10 segments Non−linear Beam, 20 segments Angular Springs, 10 segments Angular Springs, 20 segments Triangular Element, 10 segments Triangular Element, 20 segments Figure 5.12: Comparison of lateral tip position (left) and area (right) er- rors of different models with 10 and 20 segments using the mean value of parameters. 136 5.5. Discussion 5.5 Discussion Thanks to the fast computation of the presented models, and in particular, the angular springs model, they can easily be integrated into real-time med- ical training simulations systems. Due to their high accuracy, they can also be used in simulations for needle steering and path planning. Computationally, the beam element model is more efficient than the triangular/tetrahedral element model. Note that a beam element and a triangular element both have 6 × 6 stiffness matrices in 2D and, similarly in 3D, a tetrahedral element and a 3D beam element both have 12 × 12 stiffness matrices. However, modeling a needle using triangles/tetrahedra requires more elements than the beam element approach. For instance, if the needle is modeled in 2D using a mesh of 2n nodes (resulting in n − 1 segments), the tangent stiffness matrix will be 4n × 4n. In contrast, an equivalent beam model with n− 1 elements will result in a 3n× 3n stiffness matrix. Extension to 3D escalates this computational discrepancy, since a sig- nificantly larger number of elements are required for the tetrahedral FEM to model the same number of segments in 3D, whereas the beam model and, similarly, the angular springs model use the same number of nodes and elements both in 2D and 3D. Furthermore, achieving an accurate and symmetric [3] 3D model requires many tetrahedral elements, whereas the beam and angular spring models are inherently symmetric. Also note that, assuming the needle twist to be negligible facilitates the beam-element and the angular-springs models since the dimensions of the beam-element stiff- ness matrix can be reduced as shown in Appendix C and, in the angular springs case, skipping the computations involving the twist angle acceler- ates the convergence of the system in (5.33) by decreasing the number of unknowns. However, the tetrahedral model does not benefit from this sim- plification. In our simulations, the beam-element and angular-springs models often converged in fewer iterations than the triangular element model. Not only is the angular spring model robust requiring at most 6 iterations in our simulations, but its each iteration is also computed very fast since it consists of forward vector-algebraic equations with few trigonometric functions, that can be implemented using look-up tables, if necessary. Indeed, during a coupled tissue-needle simulation, the forces on the needle change smoothly. Therefore, an initial needle-posture guess from the previous time instant is often close to its final solution and thus generally only a single iteration is required to find the new equilibrium state. 137 5.5. Discussion There are alternative methods to Picard’s method that was used solve (5.34) for the angular springs model. In particular, a physically-based method based on potential energy could be developed using (5.25) as a starting model. The equilibrium configuration of the system would satisfy (5.33), with guaranteed convergence. The identified Young’s modulus did not change with the number of ele- ments in the beam element model. Therefore, the beam element model can be used in an adaptive deformation simulation scheme. In this scheme, the simulation starts with a minimum number of elements and, as the needle perforates the tissue and interacts with more tissue, the number of elements is increased (as done in [7]). Note that the given angular springs model can be simplified to neglect the twist along the shaft by ignoring/zeroing the angles αi. This is a valid assumption for many needles since the twist is often negligible due to the high shear modulus of the needle shaft and/or small twisting torques acting on it. The angular springs model as presented deals with force loads only. Nev- ertheless, torque loads can easily be integrated in the formulation by extend- ing the Jacobian definition in (5.31) as follows: τ = J ′ [ f g ] (5.37) where J i is now a 6×3n Jacobian matrix and g is the vector of torque loads at joints. Although in our formulation of the angular springs model the force loads were considered only at nodes/joints, in practice, forces at any location on the model can be accommodated, since their corresponding moment arms can be calculated regardless of them being applied at a joint or not. In this paper, equal length segments were used in the angular springs model that led to equal spring constants for all the joints. It should be noted that both bending and twisting spring constants of a segment are defined for a specific segment length and they are both inversely proportional to this length as formulated in Appendices D and E. This relation between segment length and the bending spring constant was further studied (Section 5.5.1) to devise a spring-constant approximation scheme for any given length. It was confirmed that the bending spring constant is inversely proportional to the segment length. The experiments were conducted with the needle stylet within the needle cannula. The needle is expected to bend more without the stylet (as also 138 5.6. Conclusions and Future Work observed in our preliminary observations). Experiments can be repeated in order to identify the (expectedly lower) model parameters for the same needle without the stylet. Subsequently, for a bending simulation where the stylet is partially inside the cannula, the corresponding parameter can be used for the sections of the needle with and without the stylet. 5.5.1 Approximating the Angular Spring Constants It is not always possible to determine the best segment length to use in a needle simulation at the time the bending model is developed. For a spe- cific simulation, a finer or coarser model than the one for which parameters were initially identified may be required. For each segment length, running simulations to identify individual parameters may not be feasible. Further- more, segment lengths may need to be adjusted online during a simulation. Therefore, identifying a relation between the segment length and the spring parameter is beneficial. Considering the angular springs model being the most accurate in modeling needle bending and the fastest to simulate, such a relation, that was also expected from Appendix D, is studied using the simulation results of this model. The identified angular spring constants using 10, 20, 30, 40, and 50 seg- ments were plotted in Fig. 5.13(a) for each of the vertical load experiments. In this figure, it is seen that the spring constant has a linear relation with the number of segments, which is inversely proportional to the segment length, that is in turn expected in (D.3) to be inversely proportional to the spring constant. A mean spring constant was next found for each of these five segment lengths from the values identified for all different tip loads. The change of this mean spring constant with segment length is presented in Fig. 5.13(b) in a log-log scale. Recall that the effective bending needle shaft is 186.7mm. The angular spring parameter for a different segment length can be inter- polated easily using this plot. For instance, if the parameters for 20, 30, or 40 segments were to be identified using a line fit (seen in Fig. 5.13(b)) to the values of 10 and 50 segments, these values would be within 2.5% of the parameters that were actually identified through individual simulations. 5.6 Conclusions and Future Work Three different models were presented to simulate the deformations of a needle. The first two models use tetrahedral and beam elements, while the third model uses rigid bars connected through angular springs. All the 139 5.6. Conclusions and Future Work 10 20 30 40 50 1000 2000 3000 4000 5000 Number of segments [#] Identified spring constants [µNm] (a) 3.74 4.67 6.23 9.35 18.7 103.1 103.3 103.5 Length of a segment [mm] Mean spring constant [µNm] (b) Figure 5.13: (a) Identified angular spring constants for all experiments (b) the mean spring values shown as a function of segment length in a log-log plot. 140 5.6. Conclusions and Future Work models can preserve the needle length during moderately large deformations. The efficacy of the models in simulations of needle bending was evaluated through experiments during which several lateral forces were applied to a brachytherapy needle and the resulting deformations were recorded. The model parameters —Young’s modulus for the FEM-based models and spring constants for the angular springs model— were identified to fit each model to these experimental data. The lateral tip error and the area error were independently minimized to find the parameters defining each model. Later, one single parameter value was extracted for each model and the needle de- formations simulated using these values were compared to the experimental deformations. The angular springs model demonstrated the highest accuracy compared to the other two models. Furthermore, this method is computationally the most efficient and it is the simplest to implement. It was shown that the Young’s modulus in the tetrahedral/triangular model is significantly depen- dent on the number of elements. However, the same parameter is indepen- dent of the segment length in the beam element. This property can be used in adaptive simulation to increase the speed. The models are capable of simulating the needle twist. However, twist is often negligible in medical procedures and it is not validated through our experimental setup. A single nodal force was applied to the needle in our experiments. How- ever, during insertion into soft tissue, the tissue applies a distributed force to the needle. To simulate the needle bending as a result of distributed forces using the proposed models, the distributed forces should be first integrated into nodal forces. Such conversion is common in the FEM literature. In ad- dition, a tissue model should be coupled with the needle model to simulate the forces that the tissue applies to the needle and take the insertion history into account. Such a needle and tissue coupled system is presented in [3]. Further validation of the needle models in 3D and investigation on the bevel tip force on needle bending will be carried out in the future. 141 References [1] R. Alterovitz, J. Pouliot, R. Taschereau, I. C. Hsu, and K. Goldberg, “Needle insertion and radioactive seed implantation in human tissue: Simulation and sensitivity analysis,” in Proc. IEEE Int. Conf. Robotics and Automation (ICRA), vol. 2, 2003, pp. 1793–1799. [2] S. P. DiMaio and S. E. Salcudean, “Needle insertion modeling and simulation,” IEEE Transaction on Robotics and Automation: Special Issue on Medical Robotics, vol. 19, no. 5, pp. 864 – 875, 2003. [3] O. Goksel, S. E. Salcudean, and S. P. DiMaio, “3D simulation of needle- tissue interaction with application to prostate brachytherapy,” Com- puter Aided Surgery, vol. 11, no. 6, pp. 279–288, 2006. [4] R. J. Webster, J. S. Kim, N. J. Cowan, G. S. Chirikjian, and A. M. Okamura, “Nonholonomic modeling of needle steering,” Int. J. Robotics Research, vol. 25, no. 5-6, pp. 509 – 525, 2006. [5] E. Dehghan, X. Wen, R. Zahiri-Azar, M. Marchal, and S. E. Salcudean, “Needle-tissue interaction modeling using ultrasound-based motion es- timation: Phantom study,” in Computer Aided Surgery, vol. 13, no. 5, pp. 265 – 280, 2008. [6] E. Dehghan and S. E. Salcudean, “Needle Insertion Parameter Op- timization for Brachytherapy ,” in IEEE Transactions on Robotics, vol. 25, no. 2 ,pp. 303–315, 2009. [7] D. Glozman and M. Shoham, “Image-guided robotic flexible needle steering,” IEEE Transactions on Robotics, vol. 23, no. 3, pp. 459–466, 2007. [8] X. Wang and A. Fenster, “Haptic-enhanced 3D real-time interactive needle insertion simulation system for prostate brachytherapy,” in Proc. SPIE Int. Symposium on Medical Imaging, 2004. 142 Chapter 5. References [9] A. Okamura, C. Simone, and M. O’Leary, “Force modeling for needle insertion into soft tissue,” IEEE Trans. Biomedical Engineering, vol. 51, no. 10, pp. 1707 – 1716, 2004. [10] J. R. Crouch, C. M. Schneider, J. Wainer, and A. M. Okamura, “A velocity-dependent model for needle insertion in soft tissue,” in Proc. Medical Image Computing and Computer Assisted Intervention (MIC- CAI), 2005, pp. 624–632. [11] R. Alterovitz, K. Goldberg, and A. Okamura, “Planning for steerable bevel-tip needle insertion through 2D soft tissue with obstacles,” in Proc. IEEE Int. Conf. Robotics and Automation, 2005, pp. 1640 – 1645. [12] R. Alterovitz, A. Lim, K. Goldberg, G. S. Chirikjian, and A. M. Oka- mura, “Steering flexible needles under Markov motion uncertainty,” in Int. Conf. Intelligent Robots and Systems, 2005. [13] W. Park, J. S. Kim, Y. Zhou, N. J. Cowan, A. M. Okamura, and G. S. Chirikjian, “Diffusion-based motion planning for a nonholonomic flex- ible needle model,” in Proc. IEEE Int. Conf. Robotics and Autmation, 2005, pp. 4600–4605. [14] K. B. Reed, V. Kallem, R. Alterovitz, K. Goldberg, A. M. Okamura, and N. J. Cowan, “Integrated planning and image-guided control for planar needle steering,” Proc. IEEE/RAS-EMBS Int. Conf. Biomedical Robotics and Biomechatronics, pp. 819–824, 2008. [15] S. P. DiMaio, “Modeling, simulation and planning of needle motion in soft tissue,” Ph.D. Thesis, University of British Columbia , 2003. [16] D. Glozman and M. Shoham, “Flexible needle steering and optimal tra- jectory planning for percutaneous therapies,” in Medical Image Com- puting and Computer Assisted Intervention, C. Barillot, D. Haynor, and P. Hellier, Eds. Springer-Verlag, 2004, pp. 137–144. [17] K. G. Yan, T. Podder, D. Xiao, Y. Yu, T.-I. Liu, C. Cheng, and W. S. Ng, “An improved needle steering model with online parameter estima- tor,” Int. J. Computer Assisted Radiology and Surgery, vol. 1, no. 4, pp. 205–212, 2006. [18] N. Abolhassani, R. V. Patel, and A. Farzam “Minimization of needle de- flection in robot-assisted percutaneous therapy,” Medical Robotics and Computer Assisted Surgery, vol. 3, no. 2, pp. 140 – 148, 2007. 143 Chapter 5. References [19] H. Kataoka, T. Washio, M. Audette, and K. K. Mizuhara, “A model for relation between needle deflection, force, and thickness on needle penetration,” in Medical image Computing and Computer Assisted In- tervention, LNCS 2208, W. Niessen and M. Wiergever, Eds. Springer Verlag, 2001, pp. 966–974. [20] K. Anjyo, Y. Usami, and T. Kurihara, “Simple method for extracting the natural beauty of hair,” Computer Graphics (ACM), vol. 26, no. 2, pp. 111–120, 1992. [21] E. Anshelevich, S. Owens, F. Lamiraux, and L. E. Kavraki, “De- formable volumes in path planning applications,” in Proc. IEEE Int. Conf. Robotics and Automation, 2000. [22] M. Bro-Nielsen and S. Cotin, “Real-time volumetric deformable models for surgery simulation using finite elements and condensation,” in Proc. Computer Graphics Forum, Eurographics’96, vol. 15, no. 3, 1996, pp. 57–66. [23] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method: Solid Mechanics, 5th ed. Butterworht-Heinemann, 2000, vol. 2. [24] J. N. Reddy, An Introduction to Nonlinear Finite Element Analysis. Oxford University Press, 2004. [25] I. Sharf, “Geometrically non-linear beam element for dynamics simula- tion of multibody systems,” Int. J. Numerical Methods in Engineering, vol. 39, no. 5, pp. 763–786, 1996. [26] L. Sciavicco, and B. Siciliano, Modeling and Control of Robot Manipu- lators , 2nd ed. Springer, 2000. 144 Chapter 6 Conclusions and Future Work The long term goal of this research is to develop a needle insertion path planner for prostate brachytherapy. In present practice, all needles are in- serted in the prostate along the TRUS axis, assuming that the prostate does not deform or move. A needle path planning system that makes use of in- sertion optimization could be used in a robotic system for brachytherapy. It is envisaged that such a insertion optimization algorithm will be used in path planning for brachytherapy robots to compensate for the prostate deformation and displacement and increase the targeting accuracy. In this research a simulation based approach was chosen for the needle insertion optimization. Therefore, several requirements had to be fulfilled in order to create the simulator. The desired needle insertion optimization algorithm along with the fulfillment of the requirements for the needle in- sertion simulator were set as the objectives of this research. In particular, the objectives were: 1. An experimental method to study needle-tissue interactions, necessary for the needle insertion simulator. 2. A path planning algorithm for optimized needle insertion. 3. Investigation on the effects of boundary conditions and tissue non- linearities on the optimization output. 4. Study the effects of needle flexibility on the path planning. In different chapters of this thesis steps toward this goal were taken. Partic- ularly, the first objective was achieved in Chapters 2 and 3, by introduction of an ultrasound based needle-tissue interaction study method. The inser- tion optimization method was introduced in Chapter 4 to meet the second objective of the research. The effects of the boundary conditions and tissue non-linearities were investigated in the same chapter. The needle was as- sumed to be rigid, in the force modeling and optimization. However, to pave 145 6.1. Contributions the road for future accommodation of needle flexibility in the optimization algorithm, different flexible needle models were compared in Chapter 5 to be used for brachytherapy needles. The algorithms and methods introduced in this thesis for model identi- fication and needle path planning were validated on tissue-mimicking mate- rials. In Section 6.1, the conclusions of previous chapters are discussed in the context of the thesis objectives and the contributions are highlighted. The strengths of the thesis are described in the same section. Ideas for improve- ments and the potential future works are discussed in Section 6.2. 6.1 Contributions 6.1.1 Needle-Tissue Interaction Modeling With the introduction of a novel experimental method to study needle-tissue interactions and the force model in Chapters 2 and 3, the first objective of the thesis was accomplished. Similarly to some previous work [1, 2, 3], the method was based on measurements of insertion forces and tissue displacements. However, the novelty was in the use of ultrasound imaging for the needle-tissue interaction study. Unlike other approaches that use markers [1] or beads [2, 3] or rely on the transparency of the tissue [2], this method can measure internal tissue motions without beads or markers, and hence is more readily adoptable for patient studies. A three-parameter force distribution along the needle shaft was intro- duced and used to model the interaction of the needle and the tissue. The force model has a peak force level close to the tip that contributes to cut- ting, a friction force level along the rest of the needle shaft and a parameter for the spacing between these two levels. The suitability of the model was demonstrated in an experiment during which a two-layer prostate phantom was perforated by a needle. During the insertion, the tissue was imaged using ultrasound and the insertion force was measured using a load cell. By successive simulations in two nested loops the force model and elasticity pa- rameters of the two-layer phantom were identified using data from a single insertion. However, the elasticity parameters had to be known in previous work [1, 2, 3] in order to identify the force model. At first TDPE was used to estimate tissue displacements from RF data. This method was experimentally validated and showed accurate displace- ment estimation in the axial direction of the ultrasound image. However, its 146 6.1. Contributions resolution and accuracy were lower in the lateral direction, especially when large displacements occurred. Therefore, in order to be used for patient studies during brachytherapy, the tissue displacement estimation method was changed to an RF envelope block matching method in Chapter 3. This method, which was also exper- imentally validated, showed higher accuracy – compared to TDPE – in 2D motion tracking, especially in the estimation of large displacements in the lateral direction of the US image. The insertion experiment was conducted using a TRUS to resemble the situation during prostate brachytherapy. Sim- ilarly to brachytherapy, the major tissue displacement in this case occurred in the lateral direction of the ultrasound image. The force model and elas- ticity parameters for the two-layer prostate phantom were identified in the same manner as in Chapter 2. 6.1.2 Needle Insertion Initial Parameter Optimization A simulation-based needle insertion optimization algorithm was introduced in Chapter 4. This method optimizes the insertion heading angles, insertion point and depth for a rigid needle inserted into a 3D tissue model contain- ing multiple targets which is similar to the case in prostate brachytherapy. The novelty of this method is the use of a fixed point approach to opti- mize the insertion parameters. The previously published work used search methods [4] or approximated derivatives of the cost function in a gradient descent method [5] to optimize the parameters. Optimization of the inser- tion parameters using fixed point methods increases the speed of algorithm, significantly. In addition, the proposed optimization algorithm optimizes the parameters for a 3D tissue model with multiple targets. Moreover, the previously published methods for needle path planning assumed 2D tissue models [4, 5]. The optimization algorithm was validated through extensive simulations and also using a two-layer prostate phantom, and showed significant decrease in needle targeting error. The needle-tissue interaction model and elasticity parameters for this phantom were identified in Chapter 3. 6.1.3 Investigation on the Effects of Boundary Conditions and Tissue Non-linearities In order to build a simulator for needle insertion or any other application that includes a deformable tissue model, the boundary conditions of the problem should be accurately defined. However, identification of the bound- 147 6.1. Contributions ary conditions of the prostate inside the body is challenging due to the limited available data. The effects of boundary conditions and tissue non- linearities on the output of the optimization algorithm were investigated in Chapter 4. The optimization was performed on a simplified mesh of the prostate and the surrounding tissue, with variable boundary conditions and elasticity parameters. The effects of tissue non-linearities were also investigated using a neo-Hookean material model that accommodated both mechanical and geometrical non-linearities. It was shown that boundary conditions have a significant effect on the outcome of the optimization algorithm and are the key factors in deciding whether to use a non-linear model. Indeed, when the tissue is confined and is not allowed to rotate, linear models can be used due to their low computational cost. Therefore, prior to the selection between linear and non-linear models the boundary conditions of the tissue should be known with high accuracy. While applied boundary conditions in our simulations were not realistic, a constraint around the pubic bone showed reasonable resemblance to the reported prostate rotations during brachytherapy [6]. 6.1.4 Flexible Needle Modeling The needle was assumed to be rigid in Chapters 2 to 4. However, brachyther- apy needles are flexible and may bend due to the asymmetric force on their beveled tip or physician manipulation. Therefore, modeling and simulation of flexible needles are required for brachytherapy simulators, when needle flexibility, asymmetric forces on the needle tip or physician manipulations should be taken into account. Two different finite element based needle models were compared in this thesis. Tetrahedral elements and their 2D counterparts – triangular elements - were the first models used in Chapter 5 followed by the non-linear Euler-Bernoulli beam element model. While these models were derived in 3D, they were only validated in 2D using variable tip loads. All the needle models in Chapter 5 accommodated the geomet- ric non-linearity and thus were able to preserve the needle length during moderately large deformations. The parameters for each model – the Young’s moduli - were identified in an experiment during which variable vertical forces were applied to the needle tip and resultant deflections were recorded. The models showed high accuracy in the simulation of the bending of a brachytherapy needle. Most of the time the error was less than the needle diameter. The finite ele- ment model with triangular elements showed slightly higher accuracy than 148 6.2. The Road Ahead the non-linear beam element model, but the non-linear beam element is computationally more efficient. In addition, the FEM-based models were compared with an angular springs model proposed by my colleague Orcun Goksel. It was shown that the angular springs model is the most accurate and computationally the most efficient model. 6.2 The Road Ahead The needle insertion parameter optimization is envisaged to work within a path planner for brachytherapy robots. In this scenario, the pre- or intra- operatively calculated dose delivery plan that shows the desired locations of the seeds will be given. Then a mesh of the prostate and the surrounding tissue will be generated from imaging modalities such as MRI or Ultrasound. This mesh will be used in an FEM-based needle insertion simulator. The needle insertion path planner will use this simulator to optimize the inser- tion parameters of the needle in order to compensate for tissue deformation and increase the targeting accuracy. Using these optimized parameters a robot will assist the physician to insert the needles. Therefore, in order to achieve the long term goal of this thesis – a clinically approved needle path planning algorithm for brachytherapy – the proposed methods should be validated in-vivo. The modeling and path planning algorithms of this thesis were validated using PVC phantoms, since they provide a controllable and repeatable experimental platform. Several technical difficulties should be resolved prior to in-vivo validation. In order to construct a prostate brachytherapy simulator that is able to provide sufficiently accurate results, tissue elasticity parameters, needle- tissue interaction models and the boundary conditions should be available. The needle-tissue interaction model and tissue elasticity parameters were identified jointly in Chapters 2 and 3. However, the complexities of hav- ing several layers of tissue in the needle path during the brachytherapy, precludes such joint identifications. Therefore, the elasticity parameters of tissue should be identified using other methods. Elastography methods us- ing different imaging modalities (ultrasound or magnetic resonance imaging) are useful in identification of tissue elasticity parameters, in-vivo. The simulation errors reported in Chapters 2 and 3 were measured by comparing the simulation outputs against the data which was used for pa- rameter identification. Further cross validation is required to measure the simulation error by inserting the needle with variable orientation and entry points, and comparing the simulated and measured data. However, it should 149 6.2. The Road Ahead be noted that the needle insertion optimization in Chapter 4 was performed using the simulator generated in Chapter 3. Therefore, validation results in Chapter 4 can show the accuracy of the simulator in simulating the inser- tion of the needle with two other entry points and orientation. Further cross validation and statistical analysis of the error are part of the future work. The tissue displacement measurement methods in this thesis measure tissue displacement in 2D. However, during brachytherapy, out of plane motion occurs depending on the insertion point. Therefore, 3D motion esti- mation algorithms should be developed for 3D ultrasound or other imaging modalities to measure prostate motion during brachytherapy. A very important unknown in prostate modeling is the interaction of the prostate with its surrounding tissue, in other words, the boundary condi- tions. Prostate surrounding organs and tissues are well known from anatom- ical literature and can be well studied using common imaging modalities such as MRI, CT or Ultrasound. The prostate can be easily segmented in US or MRI images. However, its interaction with surrounding tissue is not fully understood. For example, it is not known how firmly the prostate is attached to the surrounding tissue. It was shown in this thesis that bound- ary conditions have a significant effect on the outcome of the simulator and path planner. Therefore, the boundary condition for the prostate should be accurately identified for brachytherapy simulators. Investigating the move- ments of the prostate during variable mechanical excitation can be helpful for a better understanding of the interaction between the prostate and its surroundings. For example, vibro-elastography images of the prostate re- ported in [11] show a relatively softer layer of tissue around the prostate. This layer can be responsible for rotations of prostate during brachytherapy. 3D imaging of the prostate, especially during brachytherapy can provide the answer to the above question. The force model parameters were identified for constant speed of inser- tion. However, during brachytherapy the physician inserts the needle with arbitrary speed. If robots are used for insertion, the same method can be used for force model identification as insertion speed can be controlled. How- ever, force model identification for temporally variable insertion speeds is a complicated problem, which needs dynamic deformable models. Application of dynamic deformable models for tissue simulation requires identification of parameters such as mass and damping or viscosity [12] in addition to the stiffness. Needle flexibility should be taken into account in the simulators and path planners. The major reasons for needle bending during insertions are the lateral force applied to the beveled tip of the needle and needle manipulation 150 6.2. The Road Ahead by the physician. The physician manipulation effects can be simulated using the flexible needle models introduced in Chapter 5. However, modeling of the asymmetric bevel force remains a subject for further research. It should be noted that symmetric tip needles are available that do not suffer from asymmetric tip forces and therefore show less bending. In addition rotation of the needle during insertion will decrease the needle bending. In order to incorporate needle bending in the path planning algorithm, at first, the initial insertion parameters can be optimized assuming the rigid- ity of the needle. Later a flexible needle steering method such as the ones introduced in [14] or [7] can be used to steer the needle and further reduce the errors. If possible, analytical solutions for needle Jacobian should be developed (similar to [7]) to increase the speed of needle trajectory calcu- lation. Another approach is to insert the needle using the optimal initial parameters to half of the insertion depth. At this point the displaced po- sition of the target can be compared with the simulated one. Using the same optimization algorithm the new heading parameters of the needle can be calculated based on the current position of the target. Later, the needle base can be manipulated to align the needle with the new insertion path. The needle bending models were derived in 3D, although only the 2D counterparts were verified by a single lateral tip force. Verification of these models in 3D and inside tissue requires further research. 151 References [1] S. P. DiMaio and S. E. Salcudean, “Needle insertion modeling and simulation,” IEEE Trans. Robotics and Automation: Special Issue on Medical Robotics, vol. 19, no. 5, pp. 864 – 875, 2003. [2] J. R. Crouch, C. M. Schneider, J. Wainer, and A. M. Okamura, “A velocity-dependent model for needle insertion in soft tissue,” in Proc. Medical Image Computing and Computer Assisted Intervention (MIC- CAI), 2005, pp. 624–632. [3] J. T. Hing, A. D. Brooks, and J. P. Desai, “A biplanar fluoroscopic approach for the measurement, modeling , and simulation of needle and soft-tissue interaction,” Medical Image Analysis, vol. 11, no. 1, pp. 62–78, 2007. [4] R. Alterovitz, K. Goldberg, J. Pouliot, R. Taschereau, and I. C. Hsu, “Sensorless planning for medical needle insertion procedures,” in Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems, vol. 3, 2003, pp. 3337–3343. [5] R. Alterovitz, K. Goldberg, and A. Okamura, “Planning for steerable bevel-tip needle insertion through 2D soft tissue with obstacles,” in Proc. IEEE Int. Conf. Robotics and Automation, 2005, pp. 1640 – 1645. [6] V. Lagerburg, M. A. Moerland, J. J. Lagendijk, and J. J. Batter- mann, “Measurement of prostate rotation during insertion of needles for brachytherapy,” Radiotherapy and Oncology, vol. 77, no. 3, pp. 318– 323, 2005. [7] D. Glozman and M. Shoham, “Flexible needle steering and optimal tra- jectory planning for percutaneous therapies,” in Medical Image Com- puting and Computer Assisted Intervention, C. Barillot, D. Haynor, and P. Hellier, Eds. Springer-Verlag, 2004, pp. 137–144. [8] K. G. Yan, T. Podder, D. Xiao, Y. Yu, T.-I. Liu, C. Cheng, and W. S. Ng, “An improved needle steering model with online parameter estima- 152 Chapter 6. References tor,” Int. J. Computer Assisted Radiology and Surgery, vol. 1, no. 4, pp. 205–212, 2006. [9] N. Abolhassani and R. V. Patel, “Deflection of a flexible needle during insertion into soft tissue,” in Int. Conf. IEEE Engineering in Medicine and Biology Society, 2006. [10] H. Kataoka, T. Washio, M. Audette, and K. K. Mizuhara, “A model for relation between needle deflection, force, and thickness on needle penetration,” in Medical image Computing and Computer Assisted In- tervention, LNCS 2208, W. Niessen and M. Wiergever, Eds. Springer Verlag, 2001, pp. 966–974. [11] S. E. Salcudean, D. French, S. Bachmann, R. Zahiri-Azar, X. Wen, and W. J. Morris, “Viscoelasticity modeling of the prostate region using vibro-elastography,” in Proc. Medical Image Computing and Computer Assisted Intervention (MICCAI), 2006, pp. 389 – 396. [12] H. Eskandari, S. E. Salcudean, and R. Rohling, “Viscoelastic parameter estimation based on spectral analysis,” IEEE Trans. on Ultrasonics, Ferroelectrics and Frequency Control, vol. 55, pp. 1611–1625, 2008. [13] T. Podder, D. Clark, J. Sherman, D. Fuller, E. Messing, D. Rubens, J. Strang, R. Brasacchio, L. Liao, W.-S. Ng, and Y. Yu, “In vivo motion and force measurement of surgical needle intervention during prostate brachytherapy,” Medical Physics, vol. 33, no. 8, pp. 2915–2922, 2006. [14] S. P. DiMaio and S. E. Salcudean, “Needle steering and motion planning in soft tissue,” IEEE Trans. Biomedical Engineering, vol. 52, no. 6, pp. 965 – 974, 2005. 153 Appendix A 3D Line Fitting to Displaced Targets The optimization problem to be solved to fit a 3D line to some targets is pre- sented in (4.3). The optimization problem in (4.3) includes a minimization problem over the parameter α for each displaced target uki , i = 1, · · · , N . Let: qi = min α Ji = min α ‖p + αv − uki ‖ 2 (A.1) Subject to : ‖v‖ = 1, i = 1, · · · , N Taking the derivative of Ji respect to α we have: ∂Ji ∂α = 2αv′v − 2v′(uki − p) = 0⇒ α = v ′(uki − p) (A.2) Substituting for α in (A.1) we have: qi = ‖(I − vv ′)(uki − p)‖ 2, (A.3) where I is a 3× 3 identity matrix. Substituting (A.3) into (4.3), we have: (pk+1, vk+1) = arg min p,v N∑ i=1 ‖(I − vv′)(uki − p)‖ 2 (A.4) Subject to ‖v‖ = 1. Please note that (I − vv′)(uki − p) is the normal projection of the vector (uki −p) on the vector v. By changing the constrained optimization problem to an unconstrained one, taking the derivative of the new cost function 154 Appendix A. 3D Line Fitting to Displaced Targets respect to p and equating it to zero, we have: − 2(I − vv′)′(I − vv′) N∑ i=1 (uki − p k+1) = 0 ⇒ pk+1 = 1 N N∑ i=1 uki . (A.5) Now that pk+1 is known, the optimization problem can be changed to finding the optimal line – parameterized by v – that passes through pk+1 and has the minimized accumulated distances to the displaced targets. Since pk+1 and uki are known and constant, the problem changes to minimization of the accumulated magnitudes of the normal projections of uki − p k+1, i = 1, · · · , N to the prospective unit vector v. This minimization problem can be exchanged with the problem of maximization of the accumulated magnitudes of the parallel projections of uki − p k+1, i = 1, · · · , N on the unit vector v. Therefore: vk+1 = argmax v N∑ i=1 ( (uki − p k+1)′v )2 = argmax v ‖Av‖2 (A.6) Subject to : ‖v‖ = 1 where A is defined in (4.5). The right side of the (A.6) is the singular value problem and therefore, vk+1 is the right singular vector of A corresponding to its largest singular value. 155 Appendix B Neo-Hookean Hyperelastic Material Model The stored energy density for a neo-Hookean material that includes the compressibility effect is expressed as [1]: W (I,J ) = 1 2 µ (I − 3− 2 lnJ ) + 1 2 λ (J − 1)2 , (B.1) where I is the first invariant of the deformation tensor and J is the deter- minant of the deformation gradient tensor. λ and µ are Lamé parameters which are interchangeable with Young’s modulus and Poisson’s ratio. The Green strain, the second Piola-Kirchhoff stress [1] and the material moduli tensors with respect to the reference configuration can be written as: ij = 1 2 (Cij − δij), (B.2) σij = µ ( δij −C −1 ij ) + λJ (J − 1)C−1ij , (B.3) Dijkl = λJ (2J − 1)C −1 ij C −1 kl + 2 [µ− λJ (J − 1)] C −1 ijkl, (B.4) where i, j, k, l ∈ {1, 2, 3}, δij is Kronecker delta, C is the right Cauchy-Green deformation tensor, C−1ij is the element (i, j) of C −1 and: C−1ijkl = 1 2 [ C−1ik C −1 jl + C −1 il C −1 jk ] . (B.5) Discretization and approximation of these equations over several finite ele- ments lead to a set of non-linear algebraic equations of the form [1]: Φ(u)− f = 0, (B.6) where u and f are vectors of nodal displacements and external forces, re- spectively. The Newton-Raphson method is used to iteratively solve the non-linear equations of deformation [1]. Therefore, the following equation 156 Appendix B. Neo-Hookean Hyperelastic Material Model should be solved for δu at each iteration: KT(u t) δut+1 = −(Φ(ut)− f) = −et. (B.7) KT is the tangent stiffness matrix constructible from (B.2)-(B.5), u t is the vector of nodal displacements at iteration t, et is the error vector and ut+1 = ut + δut+1. (B.8) 157 Appendix C Euler-Bernoulli Beam Element The strain energy in a beam can be written as: U = 1 2 ∫ V [ E2xx + 4G ( 2xy + 2 xz )] dV, (C.1) where E and G are the Young’s and shear moduli, respectively. Substituting (5.9)-(5.11) in (C.1), one can write the strain energy as [2]: U = EA 2 ∫ ( ∂u0 ∂x )2 dx+ EIz 2 ∫ ( ∂2v0 ∂x2 )2 dx+ EIy 2 ∫ ( ∂2w0 ∂x2 )2 dx + EA 2 ∫ [( ∂u0 ∂x )( ∂v0 ∂x )2 + ( ∂u0 ∂x )( ∂w0 ∂x )2] dx + EA 8 ∫ [( ∂v0 ∂x )4 + ( ∂w0 ∂x )4 + 2 ( ∂v0 ∂x )2(∂w0 ∂x )2] dx + GJ 2 ∫ ( ∂α0 ∂x )2 dx, (C.2) where A is the cross-sectional area, and J = Iy + Iz where Iy and Iz are the second moments of inertia about the y and z axes, respectively. By substituting (5.12)-(5.15) in (C.2), the beam strain energy can be written as a function of nodal variables. Using the virtual work principle, the re- lation between the nodal variables and nodal forces/torques can be derived by symbolically differentiating the strain energy respect to nodal variables. Based on the structure of the equations for beam element (interpolation 158 Appendix C. Euler-Bernoulli Beam Element functions and energy equation), this can be written in the following form: f = K11(2×2) K 12 (2×2) K 13 (2×4) K 14 (2×4) K21(2×2) K 22 (2×2) K 23 (2×4) K 24 (2×4) K31(4×2) K 32 (4×2) K 33 (4×4) K 34 (4×4) K41(4×2) K 42 (4×2) K 43 (4×4) K 44 (4×4) u , (C.3) u = [up up+1 αp αp+1 vp φp vp+1 φp+1 wp ψp wp+1 ψp+1] ′ , (C.4) f = [fxp f x p+1 τ x p τ x p+1 f y p τ y p f y p+1 τ y p+1 f z p τ z p f z p+1 τ z p+1] ′ (C.5) where u is the vector of nodal variables, fx, fy, and f z are the nodal forces, τx, τy and τ z are the nodal torques. The blocks of the stiffness matrix K can be calculated [3] as follows: K11(i, j) = EA ∫ L Ṅi Ṅjdx , (C.6) K22(i, j) = GJ ∫ L Ṅi Ṅj dx , K 2t = Kt2 = 0, t 6= 2 (C.7) K13(i, j) = 1 2 EA ∫ L v̇0 Ṅi Ṁj dx , K 31(i, j) = 2K13(i, j) , (C.8) K33(i, j) = EIz ∫ L M̈i M̈j dx+ 1 2 EA ∫ L v̇20 Ṁi Ṁj dx , (C.9) K34(i, j) = 1 2 EA ∫ L v̇0 ẇ0 Ṁi Ṁj dx , K 43(i, j) = K34(j, i) (C.10) where L is the element length; ˙( ) = ∂ ∂x and (̈ ) = ∂ 2 ∂x2 . Blocks K14, K41, and K44 are computed similarly to (C.8) and (C.9) by substituting w0 for v0 and Iy for Iz. During the bending of a cantilever beam, such as the needle, no stretch- ing occurs on the beam’s neutral axis. Therefore, the strain component xx in (5.9) should be equal to zero on the beam axis. However, this constraint cannot be met since v0 and w0 are interpolated by cubic functions while a linear interpolation function is employed for u0. As a result, membrane locking occurs resulting in an over-stiff element [3]. One solution for the membrane locking problem is to use one-point Gauss quadrature or reduced integration for the calculation of the sub-matrices of K that include v̇0 or ẇ0 in (C.8)-(C.10) [3]. Other parts can still be integrated using two-point Gauss quadrature. As seen in (C.7), nodal twists and axial torques (the second row/column blocks of the K matrix) form a separate linear set of equations that is inde- 159 Appendix C. Euler-Bernoulli Beam Element pendent of the other variables, i.e. the lateral forces and torques. Therefore, if the twist angles or the twisting torques are negligible, these rows/columns responsible for twist can be removed to decrease the dimension of the stiff- ness matrix. 160 Appendix D Angle of Deflection for a Bending Moment Consider a short section of length L of a bent beam under a constant uniform bending moment M . Figure D.1 shows its neutral axis of length L, that is located along a radius of curvature ρ. Let c be the distance of the neutral axis from the outer fibre. From the similarity of the arcs, we obtain the following: L+ δ L = ρ+ c ρ ⇒ δ L = c ρ (D.1) The fibre strain (x) of a fibre at a radial distance x to the neutral fibre is given as: (x) = Mx EI (D.2) where E is the Young’s modulus of the material and I is the moment of inertia of the section. Substituting (D.1) in the definition of the outer fibre strain (c) = δ/L, which is the amount of elongation divided by the initial length, gives (c) = c/ρ. Substituting this in (D.2) for a fiber distance of c results in: (c) = Mc EI = c ρ ⇒ M EI = 1 ρ . (D.3) An arc of angle ω at radius ρ has a length L = ωρ. Thus, (D.3) can be rewritten as follows: M EI = L ω ⇒ ω = ML EI . (D.4) Note that the bending angle ω is equal to the angle of curvature θ. Thus, the linear relation between the bending moment and the angle of deflection can be written as θ =M/kb where kb = EI/L and, by our definition,M = ‖τ b‖ . 161 Appendix D. Angle of Deflection for a Bending Moment ñ c ¿b ! L L+ä µ ¿b Figure D.1: A short section of a bent cantilever. 162 Appendix E Angle of Deflection for a Torsional (Twisting) Moment The angular deflection α of a torsion shaft under torsional moment M can be expressed as: α = 584LM GD4 = M kα (E.1) where L is the length of the twisting section, G is the shear modulus of the material, and D is the shaft radius. Note that, according to our definition, M = ‖τ t‖ . 163 Appendix F Code Structures for Simulations and Modeling F.1 Force Model Parameter Identification Several codes have been developed to be used in different chapters of this thesis. Although the codes in each chapter can be used independently, there are several parameters that have been used by almost all the codes. These parameters are defined in Table F.1. Other parameters are defined in their related sections. This code that was used to find the force model parameters for the two- layer tissue phantom is located at: C:\MATLAB7\work\Ehsan\Phantom Data\ForceParamFitt.m. This file includes the insertion simulation code too. The user defined parameters in the code are: 1. TissueParams: a 6× 2 matrix, each row of which consists of Young’s modulus, Poisson’s ratio, shaft force density, peak force density and the peak force width of a layer of the tissue phantom. 2. USPNode: the index of nodes located on the ultrasound plane. The outputs of the code are: 1. FitError: The error between the measured (piece-wise linear) and the simulated forces, using the given force model parameters. 2. NewParams: Improved force model parameters to be used in the next iteration. These parameters will be used as part of TissueParams in the next iteration. 3. NodesOnNeedle: The nodes which are in contact with the needle. 4. X: The needle tip positions over simulation time. 164 F.1. Force Model Parameter Identification Table F.1: Parameters used in the codes NodePlaces or Xn3 The 3D position of nodes in the space. NumOfElements Number of elements in a mesh. N Num Number of nodes in a mesh. NodeList or ElementNodes The indexes of the nodes that create a tetrahedron. NodeElements All the elements that share a node. ElemTissueType The tissue type of an element, when the tissue has more than one type. NodeTissueType The tissue type of a node, when the tissue has more than one type. BoundaryNodes List of nodes that define the fixed boundaries. BoundaryCndInd The indexes of fixed boundary conditions. Example: if node i is a fixed node, the boundary index includes: 3i− 2, 3i − 1, 3i. K Inv or InvStiffness The inverted stiffness matrix for the mesh. NeedleTip The position of needle tip over the time of insertion. 165 F.1. Force Model Parameter Identification 5. IntegratedForce: The simulated insertion force with the given pa- rameters. 6. SimUSPlaneDSP: The simulated displacements of nodes located on the US plane. Several other parameters are used in this code as global parameters. The parameters are defined in Table F.1. Another input to the code is the parameter set of the piece-wise linear model fitted to the measured force. These 7 parameters are the slopes of the insertion force in each part and the starting point of each part. The parameters are saved in a file and are loaded in the program. These parameters can be easily identified using MATLAB’s fminsearch. In the first step, the program simulates the insertion of the needle into the PVC phantom using the given parameters. The major outputs are the simulated insertion force and simulated nodal displacements. In the second step, the program finds the parameters of a piece-wise linear model fitted to the simulated force to be compared to the parameters of the piece-wise linear model fitted to the measured data. In an iterative scheme the force model parameters are updated as: f t+1pi = f t pi + α( ms1i − s st1i) (F.1) f t+1si = f t si + α( ms2i − s st2i) (F.2) wt+1i = w t i + α( mpi − s pti), (F.3) where fsi and fpi are the shaft and peak force densities for the i th layer of the tissue, respectively. ms1i and ss1i are the measured and simulated slopes of force in the part contributing to puncture for the ith layer of the tissue, respectively. ms2i and ss2i are the measured and simulated slopes of force in the parts contributing to penetration for each part of the tissue. wi is the width of the peak force density for each part and mpi and spi are the measured and simulated width of the penetration part in the measured and simulated force. t is the index of iteration. Figure F.1 shows the parameters that belong to the measured force. The ForceParamFitt.m runs as a loop inside the function DispFitt.m that measures the error in simulated displacement. The code for this function is located at: C:\MATLAB7\work\Ehsan\Phantom Data\DispFitt.m. The phantom mesh parameters, including the mesh structure, the el- ement connection list, node positions, element volumes and the element interpolation parameters are initially loaded in this function. These values 166 F.1. Force Model Parameter Identification 11s m 21s m 12s m 22s m 1p m 2p m Figure F.1: The measured force and the related parameters for model iden- tification. 167 F.2. 2D RF-Envelope Block Matching Algorithm are shared with the function ForceParamFitt.m. At first, the DispFitt.m calculates the stiffness matrix with the given Young’s moduli. Then it calls ForceParamFitt.m to identify the force model parameters for the given pair of Young’s moduli and to calculate the simulated nodal displacements. At the end the function measures the error between the simulated and nodal displacements to be used for updating the Young’s moduli pair using the Nelder-Mead method. F.2 2D RF-Envelope Block Matching Algorithm This code calculates 2D displacement of several nodes using the RF-envelope block matching. The code is located at: C:\MATLAB7\work\Ehsan\Reza Code\F2DXCorr.m. The user defined parameters used in this function are shown in Table F.2: At the beginning, the function loads the file containing the RF data. The nodes of interest can be loaded into the function from a file, using a grid or by the user by clicking on the nodes of interest on the first frame. Different lines of the function should be commented/un-commented for this purpose. In a loop the function reads the RF frames and then takes the Hilbert transform to create the envelope. The envelope can be shown on the screen as an image. The RIs are taken from the initial frame, while, the SRs are taken from the subsequent frames. The RI and the corresponding SR are cross-correlated in 2D using XCorrImage2D.m. Sub-sampling in the lateral direction is performed using quadratic functions. The initial position of the SRs in each step are calculated using the previous estimated displacement of the corresponding RI. In order to reduce the de-correlation every 50th frame the initial frame is replaced by the current frame. Fig.F.2 shows the flowchart of the code. F.3 Needle Insertion Simulation and Optimization In this part the codes for needle insertion simulation and optimization are described. At first the needle insertion simulator with stick-slip force model is described. The code is located at: C:\MATLAB7\work\Ehsan\FEM\NeedleInsertion3D\Linear\ Condensed3D.m 168 F.3. Needle Insertion Simulation and Optimization Table F.2: User defined parameters used in F2DXCorr.mat LineSize The number of samples in an RF lane (depends on the US machine setting). LineCount The number of lines in an RF frame (depends on the probe and US machine settings). upSamp Upsampling factor in the axial direction (= 2). upSampL Upsampling factor in the lateral direction (= 1). BlockWidthX The length of the region of interest (RI) in the axial direction (number of samples). BlockWidthY The width of the region of interest (RI) in lateral direction (number of lines). MarginX The margin around the RI in the axial direction to create the SR (= 10 samples). MarginY The margin around the RI in the lateral direction to create the SR (= 1 line). scale Scaling factor between the axial sample and mm. Lscale Scaling factor between the lateral sample and mm. NodeNumber The number of nodes of interest. alpha Temporal filtering parameter. alpha=1: No temporal filtering. 169 F.3. Needle Insertion Simulation and Optimization Start Readthe first frame Create RIs around nodes of interest Counter =1 Read the next frame Create SRs around nodes of interest Mode(Counter, 50) =0? Calculate the cross- correlation matrix. Find the maximum of the matrix. Subsample in lateral direction. Save the displacments. End of file? Counter= Counter+1 End Replace the first frame with current frame 1 1 Yes Yes No No Figure F.2: Flowchart of 2D motion estimation algorithm. 170 F.3. Needle Insertion Simulation and Optimization Table F.3: User defined parameters used in Condensed3D.m . TissueParams A 2×4 matrix, each row of which includes the Young’s modulus, the Poisson’s ratio, the stuck node reaction force threshold and the shaft force density for each part of the mesh (first row, the prostate and the second row, the inclusion). BndryCndNodes The nodes connected to the pubic bone. NeedleVelocity Needle velocity during insertion as meter per sample time. NeedleYawRoll The yaw and roll angles of needle heading. NeedleTip The initial position of the needle tip on the surface of the tissue. Target The initial positions of targets inside the tissue. 171 F.3. Needle Insertion Simulation and Optimization This program is able to simulate the insertion of a rigid needle with the desired heading angles into a 3D tissue model to the desired depth. The user defined parameters are shown in Table F.3. In the first step the program loads the finite element model parameters including the parameters in Table F.1 from a file. These parameters are calculated using another program, TissueCreate3D.m, located at: C:\MATLAB7\work\Ehsan\FEM\NeedleInsertion3D\neoHookean\ TissueCreate3D.m. Several other functions are called inside Condensed3D.m, and are de- scribed in Table F.4. For needle insertion optimization, the initial parameters should be set by the user. After the simulation, the improved insertion parameters are shown as outputs of the simulation. For the next iteration, these parameters should be replaced the new ones manually.The flowchart for this function is presented in Fig. F.3. 172 F.3. Needle Insertion Simulation and Optimization Table F.4: Functions called inside Condensed3D.m . BChange3D.m Applies the force and displacement boundary conditions to the nodes connected to the needle in the needle coordinate frame. DSP ForceG.m Calculates the displacements and forces of the nodes connected to the needle in the global coordinate frame. FindElement3DCondensed.m Finds the element in which the needle tip is located. NH NodeJump3D.m The mesh adaptation algorithm. FindNewPlace3DCondensed.m Finds the new place of the moving node to be used NH NodeJump3D.m. 173 F .3 . N eed le In sertio n S im u la tio n a n d O p tim iza tio n Start Readthe mesh parameters Find the closest node to the tip on the surface. Relocate the node and adapt the mesh. Stuck mode? Add the node to the contact list Advance the Needle Calculate the working matrix Apply the boundary conditions in needle coordinate frame. Solve the linear equations Reaction force exceeded the threshold? Change to slip mode. Tip close to an element surface? Find the force and displacement in global coordinate frame Find the closest node on that surface Find the new position for that node Relocate the node and adapt the mesh. Add the node to the contact list Change to stuck mode. Desired depth reached? Find the positions of the targets. Fit a 3D line to the targets. Find the new parameters. End Yes Yes Yes Yes No No No No NH_NodeJump3D.m BChange3D.m DSP_ForceG.m FindElement3DCondensed.m FindNewPlace3DCondensed.m See Appendix A Figure F.3: Flowchart for needle insertion simulator (Condensed3D.m). 174 F.3. Needle Insertion Simulation and Optimization This code uses condensation to decrease the computational time. To simulate the non-linearities a neo-Hookean material model was used. A needle insertion simulator similar to Condensed3D.m was programed for this model. The code is located at: C:\MATLAB7\work\Ehsan\FEM\NeedleInsertion3D\neoHookean\ nHInsertOriented.m The user defined parameters in this code are similar to those of Condensed3D.m. The function neoHookeanSolve3D.m is called inside the program and solves the non-linear set of equations. The inputs and outputs of this file are shown in Table F.5. The function neoHookeanSolve3D.m uses Newton-Raphson iterations us- ing the tangent stiffness matrix constructed by NHTangant3D.m F.3.1 Needle Insertion Simulator with 3-Parameter Force Model This section describes the function that integrates the 3-parameter force model and uses it for needle insertion simulation. This part is not an in- dependent function, but it is part of the code used in ForceParamFitt.m explained in Section 1. The most important part of this code is its iterative nature. Assume that we integrate the force density over the length from the tip of the needle to the first node located on needle and apply it as a force boundary condition to this node. This length can be shown as: l(t) = pn(t)− (p0 + u(t)) (F.4) where l is the integration length, pn is the position of the needle tip, p0 is the initial position of the first node on the needle, u is its displacement and t is the discrete time. As it can be seen at each sample time this length depends on the value of u(t) which is part of the solution. Therefore, an iterative method is used in which, for each time sample t, the value of l(t) is iteratively updated as: lk(t) = αlk−1(t) + (1− α) ( pn(t)− ( p0 + u k−1(t) )) (F.5) l0(t) = l(t− 1), α = 0.8. 175 F.3. Needle Insertion Simulation and Optimization Table F.5: Inputs and outputs of neoHookeanSolve3D.m . Inputs Force The forces on the mesh nodes. FixedNodesIndex Indexes of fixed variables. FixedNodesDSP Values for the above indexes. U0 Initial displacement for Newton- Raphson method. Outputs U Vector of displacements. Error The final residual of the error. EE A vector including the residual of error in each iteration. 176 F.4. Simulation of Non-linear Euler-Bernoulli Beam Element where k is the index of iteration. The iterations are terminated when: lk(t)− lk−1(t) lk(t) ≤ 0.01. (F.6) The function DistributeForce.m integrates the force distribution on the needle and applies it to the nodes. The following flowchart shows the itera- tive scheme used. F.4 Simulation of Non-linear Euler-Bernoulli Beam Element The code for simulation of needle using Euler-Bernoulli beam model is lo- cated at: C:\MATLAB7\work\Ehsan\FEM\New\BeamNeedle\NewBeam\EGB3D.m. The user defined parameters in this function are shown in Table F.6. This program uses the Picard’s method to solve the non-linear equations, iteratively. The stiffness matrix is calculated by the function StiffnessProduceEB3D.m. The inputs and outputs of this function are reported in Table F.7. A re- duced integration method is used in the calculation of some of the elements of the stiffness matrix to avoid the membrane locking. 177 F.4. Simulation of Non-linear Euler-Bernoulli Beam Element Convergence criterionmet? Advance the Needle Set l (t)=l(t-1) k=1 0 Update l (t) k Integrate the force and apply boundary k=k+1 Solve the linear system Yes No DistributeForce.m Figure F.4: Flowchart for integration length update method. 178 F.4. Simulation of Non-linear Euler-Bernoulli Beam Element Table F.6: User defined parameters in EGB3D.m . BeamLength The length of the needle. NodeNumber The number of nodes used in the model. E Young’s modulus of the needle I y, I z Second moments of inertia of the needle. A The area of the cross section of the needle. IndexCnst The indexes of fixed boundary conditions. OldVector Initial value for displacement vector. Force The vector of nodal forces and torques. 179 F.4. Simulation of Non-linear Euler-Bernoulli Beam Element Table F.7: Inputs and outputs of StiffnessProduceEB3D.m . Inputs OldVector The vector of displacements. Force The vector of nodal forces and torques. NumOfElements Number of elements in the model. L Length of elements. A, E, I y, I z Defined in Table F.6. Outputs Stiffness The stiffness matrix of the model. Error The residual error. 180 References [1] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method: Solid Mechanics, 5th ed. Butterworht-Heinemann, 2000, vol. 2. [2] I. Sharf, “Geometrically non-linear beam element for dynamics simula- tion of multibody systems,” Int. J. Numerical Methods in Engineering, vol. 39, no. 5, pp. 763–786, 1996. [3] J. N. Reddy, An Introduction to Nonlinear Finite Element Analysis. Oxford University Press, 2004. 181
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Needle insertion simulation and path planning for prostate...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Needle insertion simulation and path planning for prostate brachytherapy Dehghan Marvast, Ehsan 2009
pdf
Page Metadata
Item Metadata
Title | Needle insertion simulation and path planning for prostate brachytherapy |
Creator |
Dehghan Marvast, Ehsan |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | Low dose rate prostate brachytherapy has emerged as a treatment option for localized prostate cancer. During prostate brachytherapy, tiny radioactive capsules - seeds - are implanted inside the tissue using long needles. The quality of the treatment depends on the accuracy of seed delivery to their desired positions. Prostate deformation and displacement during insertion and lack of sufficient visual feedback complicate accurate targeting and necessitates extensive training on the part of the physician. Needle insertion simulators can be useful for physician training. In addition, insertion of the needle with optimized parameters can compensate for prostate deformation, can decrease the targeting errors and, subsequently, can increase the post-treatment quality of life of the patients. Therefore, needle insertion simulation and path planning have gained a lot of attention in the research community in the past decade. Moreover, several robots have been designed for brachytherapy; however, they are yet to be coupled with proper needle insertion path planning algorithms. In this thesis, steps toward a path planning algorithm for needles are taken. An optimization method is proposed that updates the initial insertion parameters of a rigid needle, iteratively, based on the simulated positions of the targets, in order to reduce the error between the needle and several targets in a 3D tissue model. The finite element method is used in the needle insertion simulator. The simulator requires a model for the needle-tissue interaction. Therefore, an experimental method has been developed to identify the force model and the tissue elasticity parameters using measurements of insertion force, tissue displacement and needle position. Ultrasound imaging is used to measure tissue displacement. Ultrasound is a common imaging modality in the operating room and does not need beads or markers to track the tissue motion. Therefore, the experimental method can be used in patient studies. The needle-tissue modeling method and insertion parameter optimization were validated in experiments with tissue-mimicking phantoms. In order to facilitate the accommodation of needle flexibility for future applications, three flexible needle models have been compared. Based on these comparisons, it was determined that an efficient and accurate angular spring model is best suited for future studies. |
Extent | 2901593 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0065470 |
URI | http://hdl.handle.net/2429/7788 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_fall_dehghan marvast_ehsan.pdf [ 2.77MB ]
- Metadata
- JSON: 24-1.0065470.json
- JSON-LD: 24-1.0065470-ld.json
- RDF/XML (Pretty): 24-1.0065470-rdf.xml
- RDF/JSON: 24-1.0065470-rdf.json
- Turtle: 24-1.0065470-turtle.txt
- N-Triples: 24-1.0065470-rdf-ntriples.txt
- Original Record: 24-1.0065470-source.json
- Full Text
- 24-1.0065470-fulltext.txt
- Citation
- 24-1.0065470.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0065470/manifest