Modelling, Simulation and Planning of Needle M o t i o n i n Soft Tissues by Simon P . DiMaio B.Sc. (Electrical Engineering), University of Cape Town, 1995 A.Sc. (Electrical Engineering), University of British Columbia, 1998 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September 2003 © Simon P . DiMaio, 2003 UBC Rare Books and Special Collections - Thesis Authorisation Form Page 1 of 1 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of t h e r e q u i r e m e n t s f o r an advanced degree a t t h e U n i v e r s i t y o f B r i t i s h Columbia, I agree t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y purposes may be g r a n t e d by t h e head o f my department o r by h i s o r h e r r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . Department of The U n i v e r s i t y o f B r i t i s h Columbia Vancouver, Canada Date http://www.library.ubc.ca/spcoll/thesauth.html 06/10/2003 Abstract Precise needle placement is required for the success of a wide variety of percutaneous interventions in medicine. Insertions into soft tissues can be difficult to learn and to perform, due to tissue deformation, needle deflection and limited visual feedback. Little quantitative information is known about the interaction between needles and soft tissues during puncture, and no effective physicallybased training, planning and guidance systems exist for such procedures. This work aims to characterise needle-tissue interaction by measuring contact forces and deformations that are applied during insertions into soft tissue phantoms. A new methodology for estimating the forces that occur along the needle shaft during insertion is described. The approach is based on physical experiments, as well as on linear elastic phantom models that are discretised by traditional Finite Element Methods. Shaft force distributions are derived from insertions into homogeneous and simple layered inhomogeneous tissue phantoms at several driving velocities, and are applied as boundary conditions to tissue models for physically-based simulations of needle insertion trajectories. A large-strain elastic needle model is coupled to the tissue models to account for needle deflection and bending during simulated insertion. Since the force-displacement relationship is only of interest along the needle shaft, a condensation technique is shown to reduce the computational complexity of linear simulation models significantly. The boundary conditions that determine the tissue and needle motion change as the needle penetrates, or is withdrawn from the tissue model. Boundary condition and local material coordinate changes are facilitated by fast low-rank matrix updates. Such numerical schemes have been seen in prior work involving point and surface interaction; however, in this work the condensation state, boundary conditions and material coordinates evolve as the needle penetrates the tissue volume, and as internal contact states change. These novel interactive simulation techniques allow users to manipulate a three-degree-of-freedom virtual needle as it penetrates virtual tissue models, while experiencing steering torques and forces through a planar haptic interface. Models and simulations are also used to formulate needle insertion as a trajectory planning and control problem. The concept of needle steering is developed, and a Needle Manipulation Jacobian is defined to express the relationship between the needle base and tip velocities. This concept is used in conjunction with a potential-field-based path planning technique to demonstrate needle tip placement and obstacle avoidance. Results from open loop insertion experiments are also provided. n Contents Abstract ii T a b l e of C o n t e n t s iii L i s t of T a b l e s vi L i s t of F i g u r e s vii N o t a t i o n Index xii Acknowledgements xiii Preface xiv 1 Introduction 1 1.1 1.2 4 5 5 7 7 8 1.3 2 Thesis Objectives Synopsis of Thesis Results and Approach 1.2.1 Physically-based Modelling 1.2.2 Interactive Simulation 1.2.3 Steering and Planning Thesis Organisation Soft T i s s u e M o d e l s and S i m u l a t i o n M e t h o d s 2.1 2.2 Deformable Models in the Literature 2.1.1 Graphics Rendering and Animation 2.1.2 Medical Simulation . 2.1.2.1 Surgical Analysis and Planning 2.1.2.2 Surgical Training 2.1.3 Medical Image Analysis and Registration Models and Simulation Methods for Deformable Materials 2.2.1 Surface Models versus Volumetric Models 2.2.2 Implicit Surface Models 2.2.3 ChainMail 2.2.4 Mass-Spring Models 2.2.5 Continuum Models 2.2.5.1 The Finite Element Method 2.2.5.2 Condensation and the Boundary Element Method iii 10 10 11 12 12 15 18 18 19 20 21 21 24 24 26 CONTENTS 2.3 2.4 iv The Biomechanics of Tissue Needle Insertion Models and Motion Planning 27 27 3 Needle Insertion Experiments and Model Characterisation 3.1 Motion and Force Sensing 3.2 Image-based Deformation Estimation 3.3 Tissue Phantom Marker Tracking 3.4 Compression Testing of the Tissue Phantom 3.5 A Planar Linear Elastostatic Model for Tissue Phantoms 3.6 Model Parameter Identification 32 35 35 38 43 43 45 4 Needle Shaft Force Model 4.1 Force Estimates from a 2-D Model 4.2 Force Estimates from 3-D Analysis 4.3 Force Distribution Versus Driving Velocity 4.4 Force Distributions in Inhomogeneous, Layered Phantoms 4.5 Discussion 49 50 54 59 65 65 5 Physically-Based Numerical Simulation of Needle Insertion Mechanics 5.1 A Simple Insertion Simulation Algorithm 5.2 Needle Contact and Tissue Model Boundary Conditions 5.3 Non-conforming Needle Trajectories 5.3.1 Mesh Subdivision 5.3.2 Mesh Adaptation 5.4 General Non-linear Needle Trajectories 5.5 Needle Flexibility 5.6 Solving the Coupled Needle and Tissue Models 5.7 Discussion 68 68 71 73 74 75 75 78 81 85 6 Interactive Simulation of Needle Insertion Models 6.1 Condensing the Numerical System 6.2 Low Rank Boundary Condition Updates 6.3 Low Rank Material Coordinate Updates 6.4 Solving the Tissue Model 6.5 Node Intercepts and Mesh Adaptation 6.6 Haptic Needle Insertion 6.7 Discussion 88 88 89 90 91 92 94 96 7 Needle Steerability and Model-Based Trajectory Selection 7.1 Motion Planning and Needle Steering 7.2 Needle Motion Planning using Potential Fields 7.2.1 The Attraction Potential 7.2.2 The Repulsion Potential 7.2.3 Minimising Total Potential 7.2.4 The Planning Algorithm 7.2.5 Simulation Results 7.3 Experimental Results . 98 98 101 101 102 105 108 108 109 CONTENTS 7.4 8 Conclusion 8.1 Summary and Conclusions 8.2 Future Work 8.3 The Big Picture A The A.l A.2 A.3 A.4 A. 5 B Discussion Planar Pantograph Mechanism for Manipulation and Haptics Mechanism Actuation and Dynamics Force Observer Closed-Loop Position Control Haptic Control Architecture Redundant Actuation v Ill 112 112 115 117 137 138 139 140 142 145 Motion Imaging and Camera Calibration B . 1 Camera Specifications and Image Acquisition B . 2 Camera Calibration B.2.1 Camera Model B.2.2 Estimating the Intrinsic Camera Parameters B . 2.3 Image Rectification 148 148 151 151 153 156 C Linear Elasticity and the Finite Element Method C. l State of Stress C.2 Strain and Deformation C.3 Constitutive Equations C.4 The Discrete Solution C.5 Large Strain Analysis C.6 Linear Elastic Models in Two Dimensions C. 6.1 Linear Strain C.6.2 Quadratic Strain C.7 Beyond Elastostatic Models C.7.1 Dynamic Elastic Materials C.7.2 Plasticity C.7.3 Stokes Fluids C.7.4 Composite Models 158 158 160 163 164 168 170 170 171 172 172 173 173 174 List o f Tables 3.1 Summary of the 2-D and 3-D phantom model parameters used in ANSYS® 46 4.1 4.2 4.3 Summary of the planar phantom model and model discretisation parameters (2-D model) Summary of phantom model and model discretisation parameters (3-D model). . . . Summary of layered inhomogeneous phantom model parameters 50 54 65 5.1 Simulation model specifications 69 5.2 Flowchart Index 86 6.1 Model and algorithm specifications for haptic simulation 95 B.1 SI-3170 Camera Specifications (Silicon Imaging) B.2 Intrinsic parameters estimated from camera calibration images 148 156 vi List o f Figures 1.1 The Celiac Plexus Blockade (images reprinted from [1]) 1 1.2 Brachytherapy of the prostate: (top right) transverse trans-rectal ultrasound image of the prostate, and (bottom right) a fluoroscope image (images reprinted from [2]). 2 1.3 Needle insertion into a soft body 4 1.4 Needle insertion into a 2-D elastic body: (a) the contact forces exerted by the needle on the body are difficult to measure directly; however, the motion and deformation of the body can easily be measured, as shown in (b) 5 1.5 A point force Fb is applied at the boundary of a deformable body 6 2.1 Representations of a continuous body in R as (a) a surface, or (b) a solid 19 2.2 The deformation of non-rigid structures using: (a) a surface representation, and (b) a solid representation f9 2.3 A curved boundary described by a spline 20 2.4 A one dimensional chainmail structure undergoing deformation 21 2.5 A mass-spring representation of region Q 22 2.6 An under constrained spring network, with multiple equilibrium configurations. ... 23 2.7 The addition of mesh nodes requires the re-computation of adjacent spring constants and masses 24 2.8 The continuous domain O is divided into afinitenumber of discrete elements Q . . . 25 2.9 Typical loading-unloading curves of rabbit mesentery tissue (supportive membrane). 28 2.10 (a) A steered four-wheel vehicle with articulated trailer, and (b) an analogous needle approximation 30 2 e 3.1 Needle insertion and probing - experimental procedure 3.2 The complete experimental setup: a robotic manipulator with instrumented epidural needle mounted to it, a tissue phantom and a CMOS camera 3.3 A tissue phantom with markers painted on its top surface 3.4 (a) Instrumented probe, and (b) Instrumented 17 gauge Tuohy epidural needle. ... 3.5 Image-based marker tracking during (a) boundary probing and (b) needle insertion. Markers are located by a correlation-based template matching scheme (see Section 3.3). The crosses overlayed onto the images indicate the marker coordinates estimated by the algorithm 3.6 (a) Actual and imaged calibration marker locations, (b) Actual and rectified calibration marker locations (mean error of approximately 1 pixel = 0.079mm) 3.7 (a) Marker displacements due to probing, (b) Marker displacements due to needle insertion vii 32 33 34 34 36 36 37 LIST OF 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16 3.17 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 FIGURES vm A profile view of the tissue phantom, illustrating deformation perpendicular to the imaging plane A typical image taken during a needle insertion experiment. Markers are clearly visible, but vary in size and quality throughout the image A binary template is cross-correlated with image regions in order to locate markers. Tracking results for three different markers after: a) raw image capture, b) R G B to Grayscale conversion, c) erosion and dilation, d) cross-correlation with the template, and e) final coordinate estimate The tracking algorithm locates 441 markers within each frame. Estimated marker coordinates are superimposed as crosses on the camera image Compressive loading of a 30mm x 30mm x 10mm material sample. The measured stress-strain relationship is shown on the right The continuous domain O is divided into a finite number of discrete elements J7 . Element vertices constitute nodes in the mesh The tissue phantom model can be represented in two or three dimensions; however, planar material displacement measurements are only available from the top surface of the phantom Measured marker locations and estimated model node locations for a typical homogeneous tissue phantom are superimposed. A 2-D plane strain model (a) and a 3-D model (b) are compared. Error vectors are shown at each model node (scaled by 10). This phantom model has two layers of differing material stiffnesses, and is shown here with a point force applied at its boundary 38 39 40 41 42 43 e External nodal forces computed by a numerical phantom model, based on measured nodal displacements Nodal forces in the vicinity of the needle shaft are of interest. A cross-section of the force plot is taken along the row of markers that is coincident with the needle. Note that forces are also applied along the rear edge of the model; these are reaction forces from the platform constraint Estimated y-axis forces at nodes that are coincident with the needle shaft, taken at various instants during a l m m / s insertion. The needle location is shown below each graph Comparison of the measured base needle driving force and the integrated needle force distribution (2-D model) Nodal force loads are applied to a three-dimensional model of the tissue phantom in order to evaluate the 2-D analysis method A three-dimensional phantom model is used to investigate the effect of out-of-plane deformations A three dimensional model of the phantom is used to compare 2-D and 3-D analysis methods, (a) Known nodal force loads are applied to the 3-D tissue phantom model, (b) Nodal forces are estimated by a 2-D model, using the nodal displacements computed as a result of loads applied to the 3-D model, (c) Estimated nodal forces based on the 3-D analysis. As expected, the applied nodal loads are recovered with very small errors. B y contrast, note the significant errors produced by the 2-D model. . . Illustration of the cross-section of a tissue phantom during needle insertion. The top surface of the phantom does not remain flat and marker motion is not perfectly planar. Estimated force distribution along the needle shaft during an insertion at l m m / s . . . 44 46 47 48 51 52 52 53 54 55 56 57 57 LIST OF FIGURES ix 4.10 Comparison of the measured base needle driving force and the integrated needle force distribution (3-D analysis) 4.11 Estimated nodal displacements shown superimposed on actual measured nodes. The mean displacement error for 441 nodes is 0.18mm 4.12 External nodal forces computed by the planar phantom model, based on node displacements measured during insertions at 1, 3 and 5mm/s 4.13 External nodal forces computed by the planar phantom model, based on node displacements measured during insertions at 7 and 9mm/s 4.14 Estimated y-axis forces at nodes coincident with the needle shaft, taken at various instants during a 3mm/s insertion 4.15 Estimated y-axis forces at nodes coincident with the needle shaft, taken at various instants during a 5mm/s insertion 4.16 Estimated y-axis forces at nodes coincident with the needle shaft, taken at various instants during a 7mm/s insertion 4.17 Estimated y-axis forces at nodes coincident with the needle shaft, taken at various instants during a 9mm/s insertion 4.18 Estimated force distributions along the needle shaft during insertions at 1, 3, 5, 7 and 9mm/s 4.19 Shaft force density versus insertion velocity, where the shaft forces are those uniform forces that occur behind the needle tip. A 3 degree polynomial is fitted to the datapoints 4.20 Nodal forces computed by a planar phantom model, based on nodal displacements measured during insertion into an inhomogeneous tissue phantom 4.21 Estimated y-axis forces at nodes coincident with the needle shaft, taken at various instants during insertion into an inhomogeneous tissue phantom. The needle location is shown below each graph 4.22 Shaft force distribution during insertion into a layered inhomogeneous phantom. . . 58 58 60 61 62 62 63 63 64 r d 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 Comparison between measured and simulated needle insertions, using the estimated force distribution. Left: a sample of the measured and simulated node positions from one half of the phantom Comparison between measured and simulated needle insertions, with altered force distribution. Left: a sample of the measured and simulated node positions from one half of the phantom Mesh nodes lying along the needle are constrained along the ^c-axis, and either slip or stick along the ^-axis A finite state automaton for stick-slip contact between tissue mesh nodes and the needle. Such an automation governs the contact behaviour at each needle node. . . . A boundary-condition change at a contact node involves the exchange of displacement and force variables between the left- and right-hand sides of the tissue model equations. New intercept nodes are identified by searching within a small neighbourhood centred at the most distal needle node Mesh subdivision Simulated needle intercept of a small target embedded within elastic tissue 64 66 67 67 69 70 71 72 72 74 74 75 LIST OF 5.9 5.10 5.11 5.12 5.13 5.14 5.15 5.16 5.17 5.18 5.19 FIGURES x (a) During simulation, as the needle tip moves from one element into another, the intersection point between the needle axis and the common edge is determined, (b) The node closest to this intersection point, node i, is projected onto the needle axis, along the direction of the common edge, by moving the nominal position of the candidate node in the nominal reference mesh such that it lies on the needle axis in the deformed mesh state Trajectory changes (b) and deflection (c) cause local changes in needle orientation. . A 10cm, 22 gauge biopsy needle deflects as it is manipulated within a soft tissue phantom The needle is represented by a linear elastic beam Needle elasticity model with: (a) a comparison of linear and quadratic strain models, and (b) physical validation of the quadratic strain model against measured deflections of a 15cm, 20 gauge biopsy needle A step test illustrates the convergence of coupled needle and tissue models Simulated needle steering with needle flexibility. Lateral needle motion causes needle deflection, steering the needle toward the target (shown as a disk embedded in the tissue model) Contact forces that occur along the needle shaft are referred to the needle base and integrated to obtain the resultant driving force and torque A n example of the base forces and torques required to drive the needle along a simulated insertion trajectory. The trajectory proceeds as follows: (1) the needle tip makes contact with the tissue model, (2) the needle punctures the model surface and insertion continues along the y-axis, (3) insertion continues along the y-axis, with lateral base motion along the —x-axis, (4) insertion continues along the y-axis only, and (5) end of insertion. Small force spikes occur while the system converges after new working nodes are added The bevel deflection effect A flowchart of the needle insertion simulation algorithm. References for numbered boxes are indexed in Table 5.2 6.1 6.2 6.3 6.4 The relocation of node i affects the set of connected elements £ (shaded) The haptic simulation system Interactive virtual needle insertion in a planar environment Interactive virtual lung nodule biopsy 7.1 A flexible needle inserted into a body of soft tissue, toward the target T. The needle is manipulated from its base, while numerical tissue and needle models characterise motion and force at node points (a) Motion of the needle tip [it Vt 8t] with respect to the motion applied to the 7.2 7.3 7.4 76 76 79 79 80 82 83 83 84 85 87 93 94 95 96 99 T base [if, yt, determines needle "steerability". Needle steerability depends on the tissue and needle configuration q^ 100 A simple 2D tissue domain that contains a target location, as well as an obstacle that is to be avoided by the needle. Clearly, multiple starting and goal positions exist. . . 101 (a) Attraction potential over the tissue domain, and (b) the local direction of the field gradient 102 LIST OF FIGURES xi 7.5 Examples of repulsive potential field shapes: (a) circular, (b) elliptical, and (c) hemielliptical. In (d) two distinct points near the needle tip experience different repulsion forces, and a moment is applied to the line segment connecting these two points. This is the repulsive torque at the tip 7.6 (a) Repulsion potential over the tissue domain, and (b) the local direction of the field gradient 7.7 (a) Total summed potential over the tissue domain, and (b) the local direction of the field gradient at time step k 7.8 Potential field gradient directions when three obstacles are present 7.9 Simulated needle insertion along a path selected by the potential function approach. 7.10 The needle and tissue state during a simulated needle insertion along the path automatically prescribed by the potential field gradient 7.11 Simulated needle trajectory plan 7.12 Robot controlled trajectory in a transparent P V C phantom. The arrows indicate the initial insertion point. Note slight divergence due to out-of-plane deformations. . . . A.f A.2 A.3 A.4 A.5 A.6 A. 7 C.2 C.3 C.4 C.5 C.6 105 106 106 109 109 110 lfO The planar pantograph mechanism for 3-DOF manipulation and haptics f37 Pantograph dimensions and parameters f 38 Force observation using only applied actuator torques and measured joint angles. . . f40 Force observer performance f41 P I D controller for needle trajectory control 141 Four-channel architecture with local force feedback f42 The virtual environment system architecture f44 B . 1 The pinhole model projects scene points along a straight line through a projection centre (the pin hole), onto the image plane B.2 Image projection through a lens system B.3 Sample images of the planar calibration pattern, acquired using the SI-3170RGB camera B.4 Corner feature points on the checkerboard pattern are located in the calibration images, and used to evaluate the reprojection error B.5 Radial and tangential distortion fields. Image coordinate displacements due to radial distortions are significantly larger than those due to tangential distortions. Contour annotations are measured in pixels B . 6 A sample image shown before and after rectification. Straight lines and edges are restored in the corrected image C. l 103 Force and moment vectors at a point, with respect to an area element, within a deformable body under load Components of stress A non-rigid body before and after deformation The domain Q of the elastic body is discretised into elements £l Triangular and tetrahedral elements A composite fluid-elastic body under compressive and tensile loads. The grey areas are regions of fluid, while the remainder of the body is simulated using a linear elastic model e 151 152 154 155 156 157 159 159 161 165 166 176 Notation Index Notation Description Page Q a function that represents a tissue deformation model 5 u compound vector of displacements at all model nodes 5 f compound vector of forces at all model nodes 5 Q domain of a solid body 19 u compound vector of displacements at vertices of element e 166 e compound vector of forces at vertices of element e 167 cr stress vector 160 e engineering strain vector f e 160 B transformation between e and u u, x 1-tensors (vectors) 162 Qoi G.I, • • • coordinate frames 162 C material matrix that characterises Hookean elastic materials 163 E, v Young's modulus and Poisson's ratio 164 /j,, A Lame constants 164 Estrain strain energy of a solid elastic body 164 iVf(x) i 166 K stiffness matrix for element e 167 K global stiffness matrix 167 VV the set of working e e th 167 e barycentric coordinate of element e nodes in a condensed tissue model 89 Kyy condensed working (K )o precomputed stiffness matrix inverse 89 J Needle Manipulation Jacobian 99 C the configuration space of a system 29 CB obstacle regions in the configuration space 29 _ 1 stiffness matrix 89 U attractive potential field 101 F attractive field gradient 101 U repulsive potential field 102 F repulsive field gradient 102 xii Acknowledgements I would like to acknowledge my thesis advisor, Professor T i m Salcudean, who has contributed a great deal to this project. T i m has been a mentor and friend throughout my time at U B C , and I greatly appreciate the guidance and support that he has provided. In addition, I'd like to thank the thesis committee members, as well as the department, university and external examiners, whose ideas and suggestions have improved this thesis significantly. For technical assistance, and for the countless mechanical gizmos and components that have been built during this project, I am grateful to Simon Bachmann and the team of E C E workshop technicians. Many thanks to Rudiger Six, Luca Filipozzi and the IT staff for their computing support, as well as to the E C E office staff for their administrative support. Many thanks to the group of wonderful and talented people who shared the laboratory with me, and so enriched my time there. I have learned a great deal from them all, wish them the very best of success in their lives and careers, and hope to be able to stay in touch for many years to come. I am deeply grateful to my parents and family who have extended their love and encouragement over a great distance, despite the frustrations of e-mail communication and telephone calls across multiple time-zones. It is their motivation and support, shared over the many years of study, that have enabled me to challenge my goals. Some good genes may also have something to do with it! I would like to thank my wonderful group of friends, who occasionally managed to coax me away from the thesis work, in order to enjoy more social pursuits. Because of them, I now have a box full of photographs, and a head full of memories of camping trips, hiking, skiing, kayaking, sailing, climbing, horseback riding, dancing and bowling, to name but a few. Very special thanks to Francesca for her understanding and support during this time; the thesis work led to a number of emotional highs and lows, and this support was invaluable to me. Finally, I would like to acknowledge several sponsors for their generous financial support, including the K i l l a m Trust, the J . K . Zee Trust, and Raytheon Systems L t d . Support was also provided by the Canadian IRIS/Precarn Networks of Centres of Excellence, as well as the National Science and Engineering Council of Canada (NSERC). I am eternally grateful to you all. ^ xiii Preface This dissertation presents a body of work that has spanned a handful of years. It represents the final focused state of a project that has evolved from nothing but the desire to find and to develop an idea. It provides little insight into the journey from start to finish; therefore, I'd like to share—in this short preface—a few snapshots, as well as a little of what I learned from the adventure, both to aid my own memory, and to fill you in on the story. What does it mean to be a "Doctor of Philosophy", and how is it that students of literature, biology, history, education and applied science can all aspire to the same title? Their worlds seem so far apart, yet they all come together for this small pursuit. Perhaps it is all about a process of experience, rather than a specific result or discipline (of course, this is a very good definition for those who worry that the final result will not be very good, and it certainly kept me going). A senior doctoral student in our laboratory once told me that he had taken an introductory course on philosophy, so that he could feel better justified to eventually call himself a Doctor of Philosophy. This made sense at the time, although I never got around to doing the same. The Merriam-Webster dictionary defines "philosophy" as follows: f. (a) a l l learning exclusive of technical precepts and practical arts: the sciences and liberal arts exclusive of medicine, law, and theology; (b) a discipline comprising as its core logic, aesthetics, ethics, metaphysics, and epistemology; 2. (a) pursuit of wisdom; (b) a search for a general understanding of values and reality by chiefly speculative rather than observational means; (c) an analysis of the grounds of and concepts expressing fundamental beliefs; 3. (a) a system of philosophical concepts; (b) a theory underlying or regarding a sphere of activity or thought; 4. (a) the most general beliefs, concepts, and attitudes of an individual or group; (b) calmness of temper and judgment befitting a philosopher. To me, definitions f and 2 fit best, since they include words like learning, study, pursuit, search and analysis. They represent the action of learning, rather than the knowledge itself. This is why I am interested in remembering the journey, in addition to the destination: J N the beginning my thesis advisor, Professor T i m Salcudean, gave me a lot of freedom while choosing a thesis topic. I had already completed a Masters degree under his supervision, and we knew each other quite well at that stage, ft was a great opportunity to be given such freedom, but in some respects it was a disaster (at first). I agonised over the search for a suitable topic for several months without much success. Ideas ranged from haptics to computer vision, and each idea came to life with a burst of excitement, followed by a period of intense reading and justification, xiv LIST OF FIGURES XV but every one of them screeched to a halt in a blaze of disappointment after only a week or two (I probably should have paid closer attention to the etymology of the word "research"). It was at about that time that T i m , Riidiger Six (a Research Engineer in our group) and I decided to learn how to snowboard, under the watchful eye of fellow graduate student, Leo Stocco. Needless to say, there was more falling than sliding, but on one of our rides back up the hill on a ski lift, T i m recalled a research project to design a better sport brassiere for women. It had something to do with vibration analysis and resonance modes, but I no longer recall the details. Anyway, the idea of soft tissue modelling was born. Surgical simulation was very new at the time; however, several people had already done some rather interesting work i n the area. The problem of tissue modelling was—and still is—so involved, that it was difficult to imagine how we could make a meaningful contribution beyond the obvious medical training applications, for which it was not clear that particularly complex models were indeed necessary. Nevertheless, the deadline for a thesis proposal was fast approaching, so I pulled together a review of the state of the art in tissue modelling and simulation, implemented one of the more complicated modelling methods, and made a case for: "The Physical Modelling of Deformable Tissues for Surgical Analysis and Simulation." The committee did not seem to be entirely convinced at the time, but they accepted that there was some scope for a good thesis in the area, and approved my rather vague proposal. That was i n December 1999, just before the world did not come to an end. The following year was a difficult and frustrating one. The idea of surgical simulation was far too general, and we went around and around in circles looking for an interesting piece of the problem to work on. Then, after a visit to Baltimore in May 2000, a number of interesting ideas started to take shape. One of these was the modelling and simulation of needle insertions for percutaneous therapies in soft tissues. Very little work had been done in the area, but after some discussions with medical people in Vancouver, we realised that soft tissue deformation, needle deflection and rather poor visual feedback make these kinds of procedures very difficult to learn and to perform. Unfortunately, very few clinical people admitted this in the literature, so we remained uncertain about the scope of specific contributions i n the area for a long time. Nevertheless, work on the needle models continued, and by early 2001 I had a simple simulation scheme worked out, as well as plans for a series of experiments that involved sticking long needles into some sort of synthetic tissue model, while measuring forces and deformations. The search for a suitable material for tissue phantoms began. Tissue phantoms had been used extensively in medical imaging studies, where water-based materials like gelatine and agar were popular. While these materials tend to exhibit good properties for Ultrasound, M R I and X-Ray experiments, they are quite unrealistic for most applications of physical needle insertions. Both gelatine and agar offer little resistance to puncture, and crumble quite easily. O n the other hand, latex and some silicone rubbers proved to be far too tough. A search of the Internet, using keywords like: "synthetic tissue", "soft rubber", "soft plastic", etc., unearthed some rather obscene websites, as well as an abundantly confident fellow who liked to demonstrate the effectiveness of his bulletproof vests by shooting himself i n the stomach with a high-powered handgun, at close range. Eventually, I happened upon the topic of fishing lures. These are generally plastic worms, and such, that are used as bait by fishing enthusiasts. This led me to Berry's Bait and Tackle Shop i n Richmond, where I purchased a small selection of pink plastic worms and grubs, hoping that I'd be able to melt them down and re-mold them. Of course, I rushed home, threw the plastic creatures into a frying pan and proceeded to shock and horrify my roommates. The worms did melt and then re-solidify, so a week later I was back at Berry's LIST OF FIGURES xvi to purchased their entire stock of pink plastic Salmon Roe lures. The cashier enquired about the incredible fishing spot where I'd need so many lures, so I started to explain my idea to melt down the plastic eggs, and then re-mold them into body parts, but quickly abandoned the explanation, paid and left feeling rather embarrassed. The Salmon Roe did a marvelous job, but they were soon finished, and I wasn't ready to return to Berry's Bait and Tackle. I managed to find—I no longer remember how—an inexplicable character who lived on the East side of Vancouver, and who had in his possession two bags full of orange and green lures that were shaped like shrimp. He seemed to have lost interest in these soft plastic crustaceans, and so I managed to get quite a good deal on them. The little shrimp were melted down and re-molded without much trouble, but produced some rather stiff and sickly looking brown tissue phantoms. After further searching, I finally found Zeiner's Bass Shop in Wichita, Kansas. They sell liquid plastic that is easy to manage and to mold, and they sell it by the gallon. Of course the molten plastic produced some rather copious gasses, and my roommates very quickly banned me from the kitchen in our apartment. This led to a rather awkward fume hood incident in a biology laboratory on campus; however, I cannot recount that story here. The tissue phantoms led to a series of experiments and analysis, and the final set of these, the results of which are published in this thesis, were taken using a new high resolution C M O S camera. In order to reduce blurring in the camera images, I had to use a very short frame integration time. To get enough light onto the sensor during the brief frame time, I ended up with 300-500 watts worth of halogen spotlights illuminating the tissue phantom (an area of just 20-by-20 centimetres). These lights had to be driven from a D C source, as I found that under A C power, they generated light fluctuations that showed up in the images as interference. Three hundred watts was quite a lot of power for a few lights, and I needed a power supply that could produce 25-40 amps at f 2 volts. We eventually found a twenty-year-old voltage supply that weighed more than me, and sounded like an jet aircraft during take-off. The light from the halogen bulbs was so bright that I had to wear dark sunglasses while performing the experiments. The bulbs also became very hot, and would heat up the tissue phantom within minutes, so a set of small fans had to be set up to keep the phantom cool during the performance. The idea was to have a robot insert a large epidural needle into several homogeneous and inhomogeneous tissue phantoms, at several velocities. Everything had to be properly synchronised; the robot, the camera, the lights, the force sensor and two computers. The robot was a little temperamental—it did not like to be rushed—and on one occasion retaliated quite violently. 1 had just started an experiment and the robot was about to insert the needle into a sample, when it swung around and plunged the large epidural needle into the side of my hand. The experiment was over, we were into the real medicine. I would do these experiments late at night when the laboratory was quiet, and I could control the lighting. Of course there was often a lot of screaming and cursing, and I didn't want anybody to hear that. It was probably a bizarre sight: the laboratory dark, except for a point of incredibly intense light in one corner of the room, me tending to a block of re-molded plastic shrimp, wearing dark sunglasses, hair tossed about by countless cooling fans, murmuring soothing words to a robot wielding a large spinal needle. In the end, the experiments worked out quite well, and the results were used to create an interactive simulation of needle insertion. The simulator allowed users to manipulate a virtual needle on the computer screen, and feel force feedback from a haptic interface when the virtual needle was inserted into the virtual tissue model—not for the faint of heart. We ended up taking this demonstration to quite a number of conferences and exhibitions, and hundreds of people tried it. This is when I realised that most people are sadists. There was very little mercy for that helpless little virtual tissue model, and the needle quickly turned into a weapon, rather LIST OF FIGURES xvii than a medical tool. The final phase of the project was to investigate needle steering and trajectory planning using the models and simulations that were built using the experimental measurements. The idea was to find a suitable path, within a tissue model, that would take the needle tip to a specified target (e.g., a biopsy site or a tumor in the clinical world), while avoiding obstacles (e.g., blood vessels and nerves). Using the simulations, I found that regular flexible needles could be steered by manipulating them from the base. This is all explained in the thesis, but the method worked quite well, and resulted in some interesting needle trajectories. In one of the last tests, I wanted to see what would happen if there were several obstacles, so I created a tissue model with a deeply embedded target and several obstacles arranged like a slalom course. The motion plan succeeded, but this was when I began to suspect that I may also have sadistic tendencies. Other than some further tweaking of the needle insertion simulation algorithm, and some final experiments to validate the motion planning technique, the thesis work was just about done, and the writing began. Luckily, T i m and I had submitted several papers along the way, so the dissertation came together rather quickly, starting from a compilation of these papers. A n d that brings us to September 2003, the thesis is complete and the journey has been full of l e a r n i n g , s e a r c h i n g , r e s e a r c h i n g , p u r s u i t , s p e c u l a t i o n , and even epistemology, but very little calmness of temper. I cannot call myself a philosopher, but I've earned my degree by taking the journey, along with the biologists, writers, historians, teachers and scientists who were probably just as naive as I was when the adventure began. No doubt they also have some interesting snapshots from along the way... Chapter 1 Introduction Some of the most common procedures employed in modern clinical practice involve subcutaneous insertions of needles and catheters. Such procedures range in complexity from superficial needle sticks to the biopsy of deep-seated tumours, and involve subcutaneous insertions of long, slender surgical tools and needles into soft inhomogeneous tissue, usually with limited visual feedback from below the skin's surface. Needle Trajectory Figure 1.1: The Celiac Plexus Blockade (images reprinted from [1]). Physicians and surgeons often rely only upon kinesthetic feedback from the tool, correlated with 2 Figure 1.2: Brachytherapy of the prostate: (top right) transverse transrectal ultrasound image of the prostate, and (bottom right) a fluoroscope image (images reprinted from [2]). their own mental 3-D visualisation of anatomic structures. Examples of deep needle insertions for anaesthesia and brachytherapy are shown in Figures 1.1 and 1.2. Complications arising from the complexity of such interventions have been studied in biopsy [3,4], brachytherapy [5] and particularly in anaesthesia [6-10], where it is found that complications are due, in large part, to poor technique and needle placement [6]. Percutaneous procedures involving needle insertions include vaccinations, blood/fluid sampling, regional anaesthesia, tissue biopsy, abscess drainage, catheter insertion, cryogenic ablation, electrolytic ablation, brachytherapy, neurosurgery, deep brain stimulation and minimally invasive surgeries, to name but a few. Figure 1.1 illustrates a needle insertion for the Celiac Plexus Blockade, which is performed at the level of the L I vertebral body, with the patient in the prone position. The needle is passed between the vertebra and the kidney until it is anterior to the L I vertebral body. Complications related to needle misplacement can be serious, and can include puncture of the kidney, accidental intra-arterial injection, and paralysis (due to spread of the neurolytic agent into the spinal or epidural space, or due to damage of the main arterial supply to the spinal cord). Figure 1.2 illustrates needle insertion during a prostate brachytherapy procedure. Prostate brachytherapy carries the risk of seed misplacement, 3 incomplete tumor eradication, tumor seeding, and damage to the urethra and rectum, causing incontinence. Furthermore, the risks associated with the procedure must be minimised, since prostate cancer appears to progress slowly in elderly patients. Although the incidence of such complications is on the order of a few percent [6,8,9], their effects can be significant, long-lasting and even fatal [11]. Medical students and residents typically learn such procedures by performing them on real patients, under close supervision. It is estimated that between 45 and 90 attempts are required to become proficient in epidural and spinal procedures, and an even larger number of attempts may be required for peripheral nerve blocks that have less well-defined tactile queues [10]. While mannequins, cadavers and animal models do exist for training, they are often unrealistic or ethically questionable. Virtual-reality-based training systems for catheter insertion [12-14], endoscopy [15,16], epidural lumbar puncture [7,17-24], spine biopsy [25], breast biopsy [3], neurosurgical probe insertion [26], interstitial brachytherapy [27], prostate needle biopsy [28], etc. constitute part of the present trend toward computer-based simulators for medical and surgical training [29-3f ]. Experienced physicians and surgeons are also taking advantage of these developments by using simulators to plan and to rehearse complex procedures, to predict the outcomes of procedures, to design and evaluate new methods and equipment [32,33], and to control complex medical robotic systems [34-38]. The majority of the simulation and manipulation systems for needle- and catheter-based procedures have been built using largely phenomenological and heuristic models that have not been validated, and that are not generalizable. In some cases, lookup tables based on empirical values determine needle force feedback versus depth for f - D O F interaction [f3,19]. Others employ elastic, viscous and viscoelastic models to approximate needle driving impedance along the axial direction [17,18]. While perhaps effective for the simulation of predominantly f - D O F problems, these approaches will not be suitable for problems involving more complex soft tissue anatomy, needle placement optimisation, trajectory planning and automatic control. Indeed, in these cases, more detailed verifiable knowledge of the biomechanical interaction between surgical needles and soft tissues is required. In prior work, needle insertion forces were determined for gelatine [20], ex vivo porcine, human and bovine samples [f 7, 39]. fn each case, only the resultant force acting at the proximal end of the needle was measured, while in fact penetration forces are distributed along the entire length of 1.1 Thesis Objectives 4 the needle axis, resulting from physical phenomena such as cutting, elastic deformation and friction (static, kinetic and viscous) [17]. The needle driving forces measured previously are the integration of this force distribution along the needle shaft. Consider the insertion of a needle into a volume of soft material, as depicted in Figure 1.3. In Figure 1.3: Needle insertion into a soft body. general, as the needle penetrates the tissue, forces are applied both axially and laterally and the tissue deforms, affecting the path of the needle and in turn the force distribution itself. The resultant proximal force does not provide sufficient information for the tissue deformation to be determined. As it will be shown in this thesis, the distribution of forces along the needle plays an important part and should be determined in order to gain insight into the more general biomechanics of needle insertion procedures. Recent work on model-based sensitivity analysis for prostate brachytherapy has already made use of this concept [40]. 1.1 Thesis Objectives This thesis aims to explore modelling and simulation methods for needle insertion applications by physically observing insertions into soft tissues. Methods for measuring the forces applied by the needle shaft, and the resulting tissue deformations, are required in order to quantify the mechanics of puncture in a generalisable and verifiable way. Such quantification has not been available in prior studies. A further objective of this work is to show that a physically-based simulation model can describe the observed insertion mechanics in a satisfactory manner, and that such models and simulations can be applied to training, guidance and planning tasks. 1.2 Synopsis of Thesis Results and Approach 1.2 5 Synopsis of Thesis Results and Approach The scope of this work includes physically-based modelling, interactive simulation, needle steering and motion planning. 1.2.1 Physically-based Modelling In order to better quantify the mechanics of needle insertion, a system has been developed to observe planar insertions during physical experiments. The objective of these experiments is to determine the contact forces that occur along the needle shaft during puncture, as illustrated in Figure 1.4(a). The methodology has been developed in two dimensions, for the sake of exposition and experimental validation, but can easily be extended to problems in three-dimensions. (a) (b) Figure 1.4: Needle insertion into a 2-D elastic body: (a) the contact forces exerted by the needle on the body are difficult to measure directly; however, the motion and deformation of the body can easily be measured, as shown in (b). Consider a homogeneous, elastic, deformable body as forces are applied to its boundary by a blunt probe. This is illustrated in Figure 1.5. Probe forces and body deformations are measured, and used to characterise the following relationship: u = 0(f) f = G~ (n). l (1.1) If the function Q exists, then it represents a model that relates local material displacements u and 1.2 Synopsis of Thesis Results and Approach 6 Figure 1.5: A point force Ft, is applied at the boundary of a deformable body. external forces f over the entire body. Q is not a function of time, since it is assumed that body accelerations and velocities are relatively small during typical needle insertions. Equation (1.1) expresses the quasi-static behaviour of the deformable body, but does not exclude second-order dynamic effects along the needle-body contact. During a needle insertion (see Figure 1.4(b)), the body deformation can be measured, but the contact forces applied by the needle in the interior of the body are not known. Providing that the inverse model Q~ exists, the set of local externally applied forces f can be computed using l the displacements u that are observed. The external forces computed along the needle shaft are estimates of the actual needle forces applied to the body, since Q is an approximation. In this thesis, tissue phantoms are used instead of biological tissues during experiments because they offer a controlled environment for repeatable measurements; whereas, the inhomogeneous nature of biological tissues makes such experiments and validation infeasible. This approach also ensures that an effective, yet tractable model Q exists. It is unclear whether cx vivo tissue samples can provide more relevant results, due to mechanical changes that occur post mortem and during sample preparation. In addition, determining tissue parameters, either ex vivo or in vivo, is not trivial [41-43]. In the phantoms selected for this work, the base driving forces measured are of the same order of magnitude as some of the ex vivo forces presented in the literature, and that we have observed in porcine liver samples. Insertion forces measured in vivo are likely to be lower due to lubrication, and higher order effects may be more significant. Nevertheless, results from phantom experiments will help to refine the methodology for future in vivo measurements. 1.2 Synopsis of Thesis Results and Approach 1.2.2 7 Interactive Simulation Simulations of multi-degree-of-freedom trajectories of flexible needles into elastic tissue models are demonstrated using estimated shaft force distributions that are applied as contact boundary conditions to the tissue model Q. Physically-based material models and contact force distributions are used for this purpose. Others have already adopted this approach for sensitivity analysis and trajectory estimation [40,44]. Interactive training systems require model solution rates that match human sensory perception at visually and kinesthetically convincing rates on the order of 30Hz and 1000Hz respectively. Techniques for achieving such performance for insertion simulation are presented. 1.2.3 Steering and P l a n n i n g The placement of needles for percutaneous intervention is usually performed using straight-line needle trajectories, likely due to the absence of sufficient visual feedback or knowledge of needle steering mechanics required to make trajectory compensation feasible. The insertion point (at the skin surface) and needle orientation are selected in order to set up a straight-line intercept with the target. If a needle fails to reach its goal, due to tissue or needle motion, then it must be retracted and re-inserted. Several attempts may be required before precise placement is achieved. In some cases, attempts are made to steer the needle. For example, there is anecdotal evidence of needle steering using bevel deflection and direct manipulation of the needle between the template and the perineum in prostate brachytherapy. Steering mechanisms have been developed for catheters and needles [45-48]. Needle placement planners have been used to identify optimal seed placement and biopsy sites [4,28], as well as bone interference [49], for prostate brachytherapy and biopsy. However, misplacement due to tissue deformations and needle deflection have been largely ignored. Similar targeting issues exist in other percutaneous procedures, particularly when fine needles are used. The insertion models and simulation algorithms developed in this work are used to investigate needle steering and trajectory manipulation techniques for precise needle placement. The needle insertion problem exhibits a large number of degrees-of-freedom, corresponding to an extremely high-dimensional configuration space; therefore, as discussed in Section 2.4, global search methods are not feasible. A local search method, based on potential fields, is demonstrated in this work. The approach is shown to produce feasible needle trajectories in simulation, and is 1.3 Thesis Organisation 8 demonstrated experimentally. Needle motion planning in the presence of deformable tissues and obstacles has not been considered before. 1.3 Thesis Organisation The thesis chapters are organised as follows: Chapter 2: Prior work in the area of deformation modelling, simulation and motion planning is presented, primarily in the context of medical and surgical applications. Chapter 3: This chapter describes a measurement system that has been developed to model needle insertion mechanics. Forces and planar tissue phantom deformations are measured during experiments in which probing and puncture trajectories are effected by a robotic mechanism, while motion and deformation are measured by a high-resolution imaging system. A n elastic model is shown to adequately describe the tissue phantom deformations observed. The equipment, experimental methods and model characterisation are presented. Chapter 4: This chapter presents a new model-based method for determining forces that occur along the needle shaft during soft tissue puncture. Tissue phantom deformation measurements and an inverse material model are used to estimate material forces that are applied by the needle. Needle shaft force distributions are derived over a range of insertion velocities, as well as in layered inhomogeneous phantoms. Chapter 5 : A numerical simulation of needle insertion is described, based upon the experimentally derived needle force distributions. The algorithm simulates complex multi-degree-of-freedom needle trajectories using an elastic tissue model and a large-strain needle deflection model. Chapter 6: This chapter details an interactive virtual needle insertion environment, in which a user can experience needle forces while manipulating a virtual needle within soft tissue models. The environment is based on the simulation methods described in Chapter 5. Fast numerical techniques are developed to enable the simulation of reasonably complex models at fast "haptic" rates. 1.3 Thesis Organisation Chapter 7: 9 Needle manipulation, steering and motion planning are discussed in the context of needle placement. This chapter outlines a new motion planning approach that is based on potential fields, and demonstrates its effectiveness in simulation and experiments. Chapter 8: Thesis contributions and limitations are summarised, and the scope for future work is discussed. Appendix A : The 3-DOF planar pantograph mechanism is used as both a manipulator and a haptic interface in this work. Details of its kinematics, dynamics and control are provided. Appendix B: This appendix provides detailed specifications for the imaging system that is used to measure tissue motion during experiments. Acquisition and calibration issues are presented. Appendix C: A primer on the basic continuum mechanics of linear elastic materials and the Finite Element Method. Details of both small and large strain techniques are provided. A number of inelastic material models are also discussed. To date, portions of this thesis have been published in [50-54]. Chapter 2 Soft Tissue M o d e l s a n d Simulation Methods Models that describe the behaviour of soft, deformable materials and tissues are of general interest for surgical simulation and planning systems. The needle insertion modelling and simulation methodologies described in this work rely upon the existence of such models. This chapter provides an overview of deformable models and simulation methods that have been presented in the literature, with emphasis on medical applications. 2.1 Deformable Models in the Literature The dynamics of physically deformable bodies are of interest in a number of disciplines including mechanics, biomechanics, computer graphics, image processing, computer vision and medicine. In each of these fields, techniques have been proposed for measuring, modelling, characterising and simulating the behaviour of non-rigid materials. While the mechanics and dynamics of rigid bodies and rigid body contact have been widely studied, knowledge of deformable materials—their behaviour and co-interaction—is still limited. Studies of the constitutive equations for metals, polymers, fluids, gasses, magnetic and electromagnetic fields have led to an understanding of the basic principles of these and other continua, and hence the study of continuum mechanics. Progress in this area has catalyzed several significant advancements in both the applied sciences and physics. New and innovative applications of these principles are currently emerging, particularly within the field of medicine, where severe demands are placed upon existing methods. This section outlines a variety 10 2.1 Deformable Models in the Literature 11 of these applications and the approaches that have been adopted. 2.1.1 Graphics Rendering and Animation In the quest for ultimate realism in graphical modelling and animation, researchers are turning to physically-based models of natural phenomena for image synthesis and rendering. While previous work employed rigid bodies and sometimes several articulated or arbitrarily constrained rigid bodies to simulate physical environments, models based upon continuum mechanics have provided the means to express and to simulate active deformations of non-rigid bodies in a natural manner. Terzopoulos and Fleischer [55, 56] explored "three canonical genres of inelastic behaviour", namely viscoelasticity, plasticity and fracture to visually convey the qualitative dynamic characteristics of rubber, cloth, paper, string, Plasticine and flexible metals. Although the emphasis was on the development of realistic, yet computationally tractable animations, their Finite-Element-based numerical simulation techniques become expensive for the dynamic simulation of complex inelastic models in three-dimensions. Nevertheless, Terzopoulos and Fleischer's techniques do yield visually appealing renderings of a range of natural phenomena and non-rigid deformations. Chen and Zeltzer [57] extended some of these ideas to the modelling of muscles for the realistic animation of humans and animals. A n elastic Finite Element model simulated the passive behaviour of the muscle tissue, while non-linear empirically motivated force generators were connected to the mesh node points. The total number of mesh nodes was limited to 56 in order to achieve reasonable simulation times. The computational complexity associated with models based on Finite Element Analysis has spawned a variety of alternative approaches, including deformation modelling using implicit surfaces [58], spline-based free-form and extended free-form deformations [59], super quadratics, deformation polynomials and modal analysis [60]. In each case the deformable shapes are described by a relatively small set of model parameters. This parametric efficiency is offset by difficulties incorporating physically plausible mechanical properties. Mass-spring network representations of non-rigid bodies are able to exhibit elementary physical properties and have been popular for animation and image synthesis [61] due to their simplicity and scope for efficient computation. James and Pai [62] recently demonstrated the effectiveness of the Boundary Integral Formulation for simulating linear elastic volumetric deformations. In this approach, only the state of the object boundary is solved, thus significantly improving computational efficiency, while maintaining global volumetric behaviour. They recently extended the approach to multi-resolution, multi-body and dynamic systems [63-65]. 2.1 Deformable Models in the Literature 2.1.2 12 Medical Simulation The human body is a system of organs, tissues and fluids that are coupled by interactions that are extremely complex. Medical science is often required to understand the dynamics of such interaction for the improved diagnosis, analysis and treatment of pathologies. In many cases this involves the modelling and simulation of soft tissues. Two categories of soft tissue simulation have been explored in prior work, and are indicated to have distinctly different constraints, in terms of model accuracy and computation time [66]: Surgical Analysis and Planning: For the analysis and evaluation of new procedures, surgical tools, implants etc., as well as for effectively evaluating and comparing the outcomes of interventions and techniques prior to surgery. Model accuracy and fidelity are of paramount importance, while simulation time is secondary. Surgical Training: Real-time performance (i.e., >30Hz update rates) is important for interactive surgical training, while accuracy and detail may be compromised to some degree. The extent of this compromise appears to vary from procedure to procedure. 2.1.2.1 Surgical Analysis and Planning Soft tissue modelling has been of particular interest for computer-assisted planning of craniofacial surgery [32,67-72], where injuries and congenital deformities are repaired by modifying the shape of the skull. In a typical procedure, the surgeon cuts (osteotomises) the skull into several fragments and resets them into a normal configuration. In the case of Maxillofacial or Orthognathic surgery, this remodelling of the skull is used to correct a malaligned jaw, an affliction that often affects the patient's ability to talk or chew properly. Early work in this area involved the manual integration of physical examination results, cephalometric radiographs, anthropometric measurements, and dental models by a team of specialists. Measurements from these images were compared with normal values in order for the surgeon to formulate a surgical plan, and to be able to predict the expected post-operative facial contour. Advances in computing and medical imaging made it possible to integrate, visualise and to manipulate this data more effectively and intuitively. In 1993, Altobelli et al. [67] presented a computer-assisted three-dimensional planning system that enabled the surgical team to simulate osteotomies and skeletal movements in three-dimensions using overlays of normative data (Bolton standards) as a guide. A n important aspect of the planning process is the 2.1 Deformable Models in the Literature 13 ability to estimate functional and aesthetic changes to the skin surface, resulting from motion of the underlying bone. Altobelli et al. based their estimates of skin changes on ratios between bone and skin movements found in the orthognathic literature. Herve Delingette et al. [68] used skull and skin meshes, interconnected by "muscle springs", to simulate the dependence of the skin surface shape on bone movements. Their system employed a "virtual hand" glove interface to facilitate interactive cutting and positioning of the skull. Keeve et al. [70] used a linear elastic tissue model to better predict skin behaviour during the planning process. Although this approach was shown to work quite well, the authors note that better results may be attained by using a non-linear elastic function to describe tissue behaviour. Sarti et al. [71] followed a similar approach, electing to use an isotropic elastic model. Their numerical simulation of a "Sagittal Split Mandibular Remus Osteotomy", with a 2,160,000 node volumetric mesh, was processed on a 128-node Cray T3E supercomputer. Schutyser et al. [72] traded realism for computation efficiency by selecting a voxel-based method for their tissue models, while Teschner et al. [32] used an elastic soft tissue model based upon a spring representation. The latter approach is reported to require 3.3 seconds to compute soft-tissue changes using a relatively standard computing workstation, with a mesh size of approximately 2000 nodes. This system also simulates jaw movement for the evaluation of post-operative physiological functionality. In the field of Neurosurgery, the now standard frameless stereotactic guidance technique allows the surgeon to use pre-operative calculations regarding intracranial targets, such as tumors and lesions, intra-operatively. This strategy assumes that the medical images and the brain remain in perfect rigid correspondence during the surgical procedure, so that there is a precise mapping of anatomic structures between the images and the patient. In practice, significant brain shifts and deformations occur in open microsurgical operations [73], and are shown to be on the order of 4-7mm during tumor resection procedures [74], and greater than 10mm following a craniotomy [75]. In response to this problem, several image-based compensation schemes have been proposed [73, 75-78]. In each case, the pre-operative images and navigation models are updated based upon intra-operative CT, MRI or X-Ray images that can be expensive, impractical and even dangerous. A computational model-guided approach provides an alternative means for coping with the brain shift problem, by attempting to predict brain tissue deformations due to intracranial pressure changes and surgical loads. Auer et al. [33] used a spring-damper approximation of brain tissue to predict 2.1 Deformable Models in the Literature 14 deformations due to surgical load. This model was not physically-based and was not validated against measured brain deformations. Miga et al. [79] proposed a more advanced tissue model described as a homogeneous biphasic medium—a solid elastic matrix saturated with an interstitial fluid. It was reported to account for 70-80% of brain surface motion, based upon interventional imaging in a human clinical case; however, substantial differences in subsurface deformation were noted. These experiments were shown to be particularly sensitive to the incorporation of the Falx Cerebri, a prominent structure in the brain, into the tissue model. This observation motivated the study of other support structures as a means of further improving model accuracy. For example, studies that are based upon experimental loading of brain tissue samples in vitro, provided some material properties (e.g. [80]); however, the measurements from uniaxial loading of homogeneous tissue samples does little to characterise the dynamic behaviour of a heterogeneous brain in vivo. Similar modelling issues emerged in the study of cardiac motion mechanics for diagnostic purposes. The characterisation of heart wall motion is helpful for understanding the process of diseases like ischemia, and can be used to estimate the location and extent of ischemic myocardial injury. The Spatial Modulation of Magnetization ( S P A M M ) method is one of a variety of non-invasive techniques for tracking the motion of the myocardium. The three dimensional locations of several tissue points are tracked simultaneously by selectively altering tissue magnetisation at discrete locations within the heart muscle, and tracking these over a series of magnetic resonance images [81-83]. Park et al. [84] use this motion information in conjunction with a volumetric model to capture the deformation and motion of the left ventricle. A set of locally varying global parametric functions was fitted to the magnetically tracked points, in order to capture motion and thus describe complex volumetric behaviour, such as twisting and ejection. Bardinet et al. [85] pursued a similar approach for the analysis of left ventrical deformation, based on the combination of parametric superquadratic and free-form deformation functions. Their volumetric model was shown to track the deformation of both the endocardium and epicardium walls, from which quantitative wall thickness, volume, ejection and twist could be inferred, in much the same way as described by Park et al. Their reliance on pre-segmented inner- and outer-wall surfaces, rather than the volumetric information acquired by S P A M M , limited the accuracy of motion estimates. Papademetris et al. addressed this problem by using a linear elastic tissue model to regularize cardiac deformation tracking [86]. Both the endo- and epi-cardial surfaces were interactively extracted from 3D ultrasound images 2.1 Deformable Models in the Literature 15 and fitted to a geometric model by means of a snake-like tracking algorithm. In contrast to earlier schemes, the registration of images and models was regulated by a physically-based tissue deformation constraint. For this, a simple anisotropic linear elastic tissue model was adopted. The anisotropic elastic property allows the deformation model to account for the effect of fibers in the ventricle wall, which are known to cause twisting [87]. This strategy was shown to identify areas of myocardial dysfunction—caused by coronary occlusion—by locating low-strain regions in the heart wall (see [83] for a brief discussion of myocardial dysfunction and infarct). Computation times on the order of 3-4 hours were required to process 13-17 ultrasound frames taken over a single cardiac cycle (Silicon Graphics Octane R10000, 195MHz), using a Finite Element model containing approximately 2500 hexahedral elements. Soft tissue models have also been of use in orthopaedics. Calisse et al. [88] were able to predict trunk muscle forces and spinal loads during standing and flexion of the upper body. Their 3000-node model included elastic models of the L I to L5 vertebrae, incompressible fluid-filled intervertebral disks, dorsal and ventral trunk muscles and homogeneous visco-elastic abdominal tissue. The performance of the resulting model was limited by long computation times, numerous assumptions and simplifications regarding material properties, and the crude abdominal model. Nevertheless, predicted loads matched implant measurements fairly well (measurements were taken using instrumented fixators), and the estimated angles between adjacent vertebrae were within the ranges derived from in vivo measurements. 2.1.2.2 Surgical Training Training doctors and surgeons for complex procedures is often difficult and slow. Skills are acquired by using textbooks, mannequins, cadavers and live animals, or by observing and assisting an experienced surgeon. Cadavers and mannequins provide a poor alternative to living tissue since they exhibit unrealistic macroscopic behaviour and there is no bleeding. The anatomies of animals differ from that of humans and there are ethical issues related to their use for training purposes. Training on actual human patients poses increased risks of complications, misdiagnosis, mistreatment, greater pain, morbidity and even death. In addition, greater demands are made upon the instructor and the surgical unit during training. This has prompted the development of a variety of virtual-reality-based training devices ranging from crude mechanical apparatus [89] to hi-tech 2.1 Deformable Models in the Literature 16 interactive multimedia-based simulators [15]. In addition to surgical training, such systems can also help to plan complex surgical procedures, predict the outcomes of surgical procedures, or to design and evaluate new procedures, tools and hardware [29]. Typical surgical simulation systems comprise visual displays and input devices for measuring the surgeon's actions. Some form of kinesthetic or tactile feedback may also be provided [90]. A good example of the importance and difficulty of training is the diagnosis of prostate cancer, which is the second leading cause of cancer death amongst men [91]. The most common method of detection is the digital rectal examination, during which the physician physically feels and palpates the prostate for tissue nodules that characterise tumors. This is uncomfortable for the patient, particularly when being performed by an inexperienced examiner, fn addition, it is difficult for the instructor to evaluate performance. Burdea et al. [9f] used a simple deformable mesh to describe the prostate surface, and created an interactive palpation trainer using a P H A N T o M ® haptic interface [92]. Interaction forces were simulated by means of a variation on Hooke's law, and a simple "region of influence" or propagation law that was also used by Popescu et al. [93] and Schutyser et al. [72]. The model was not physically based, but led to a computationally efficient palpation model. Human factors studies indicated a need for more realistic physical modelling. Similar heuristic models were employed by Dinsmore et al. [94] for the palpation of subsurface tumors in the liver; and by Langrana et al. [95] for virtual knee palpation. Others have studied the mechanics of palpation for the location of tumors in the lungs and breasts, and have reported tissue behaviour far more complex than those described by the abovementioned models [96-98]. A related diagnostic problem is the detection of Deep Vein Thrombosis (DVT), which can lead to fatal pulmonary emboli. The most popular non-invasive means of detection is Doppler Ultrasonography with Compression, where the compliance of the blood vessels of the lower limbs is measured by palpating or compressing the tissue. B y observing vessel wall deformations in the ultrasound images, an experienced examiner can detect thrombi [99-102]. Approximately 1000 echographic examinations are typically required in order to acquire acceptable competence in diagnosis. d'Aulignac et al. [103] proposed a massspring model of the human thigh—based upon empirical data—for the haptic simulation of thigh echography. The non-linear spring parameters that characterised their 2-layer mass-spring network were determined by means of a least-squares fit of several measured force-displacement relationships. The resulting surface model characterised only skin surface deformation and did not include vessels. 2.1 Deformable Models in the Literature 17 There has been significant interest in the development of simulators for laparoscopic surgeries, due to their high level of difficulty. During a laparoscopic procedure, large incisions are avoided by inserting the surgical tools and a small camera into the patient through small skin punctures. This approach reduces patient trauma and recovery time, however significantly increases the difficulty of the operation. The surgeon does not have direct contact with the tissue, and visibility is limited by monoscopic endoscope images. The technique requires a great deal of practice and training before it can be attempted by a novice surgeon [104]. Cotin et al. implemented a deformable liver model for the simulation of laparoscopic hepatic surgery [105-107]. Their volumetric model, based on a simplex mesh of the liver, exhibited quasi-static linear elastic behaviour approximated by the Finite Element Method; a 1400-node liver model could be computed in approximately 7ms ( D E C A l p h a 400MHz workstation) by taking advantage of superposition. The pre-computation step took approximately 8 hours. B y employing a hybrid model, Cotin et al. were subsequently able to incorporate simple mesh topology changes resulting from cutting and fracturing [108]. Small regions of tissue lying within the desired neighbourhood for cutting or fracture were simulated volumetrically without pre-computation. The additional load was minimised by restricting the number of mesh nodes in these regions. Cotin also included quasi-nonlinear deformation behaviour in order to slightly improve model correspondence with rheological experiments conducted by Chinsei and Miller [109]. Boux de Casson and Laugier also modelled the dynamics of the human liver for the simulation of minimally invasive surgery [110]. Their heterogeneous network of spring-damper elements simulated the outer layer of the liver ("Capsule of Glisson") as an elastic surface, while the internal material was treated as a visco-elastic volume. The resulting deformation exhibited somewhat more complex behaviour than that of Cotin et ai, including visco-elasticity and plasticity; however, the ad hoc selection of network parameters was not physically based. The resulting computation, performed by a 300MHz dual-Pentium II system, required 3.3 seconds to process 1224 surface facets and 6291 internal tetrahedrons. Finally, a system for Bronchoscopy training was developed at HT Medical Systems by Bro-Nielsen et al. [15]. Their PC-based system blends multimedia, 3D graphics and passive force feedback in what is essentially a state-machine. A database of visual effects are stored for playback during the simulated procedure, in response to the student's actions. Passive force feedback reflects the effects 2.2 Models and Simulation Methods for Deformable Materials 18 of friction due to the motion of the flexible endoscope. 2.1.3 M e d i c a l Image Analysis and Registration Medical images taken from several different imaging modalities, from different patients, or at different time instants, can contain useful complimentary information, provided that they are satisfactorily aligned. The image registration problem has been solved in a variety of different ways [111-114], but is generally posed as an optimisation problem, the purpose of which is to find a rigid, affine, projective or non-rigid transformation that results in optimal global or local image registration according to some similarity metric [73, 77, 78,115-118]. Non-rigid registration schemes have become important tools for accurate intraoperative guidance [78]; intra-subject and subject-model matching [116,119,120]; progression analysis [121] and tissue motion tracking [82,85,86,122]. There is particular interest in maintaining accurate correspondence between intra-operative images and pre-operative medical images, models and plans. Because of the deformable nature of soft tissues, non-rigid transformation schemes are necessary. Early non-rigid transformation models used stabilising or regularising terms, such as smoothness or continuity constraints to guarantee convergence [118,123]; however, these typically produced physically implausible deformations. Physically-based tissue deformation models were introduced to stabilise non-rigid registration by establishing a physical basis for the deformations. Bajcsy and Kovacic; Gee et al. and Ferrant et al. used simple elastic models [116,120,124]; however, these were found to perform no better than standard Gaussian smoothing or propagation methods, due to their linear homogeneous behaviour. Recent results obtained by Hagemann et al. suggest that more detailed physically-based models do in fact lead to more accurate and robust non-rigid registration [78,125,126]. Their tissue deformation models incorporate coupled rigid, elastic and fluid structures that are based upon local tissue characteristics. Davatzikos et al. used a combination of statistical and biomechanical models for predicting anatomical deformations in medical images [127]. 2.2 Models and Simulation Methods for Deformable Materials Although it is clear that the incorporation of non-rigid tissue behaviour is important for a large number of applications, it is the numerical representation and simulation of these behaviours that 2.2 Models and Simulation Methods for Deformable Materials 19 limit widespread use. This section outlines several different deformation models, their performance and utility. 2.2.1 Surface Models versus Volumetric Models Figure 2.1 depicts a body Q, in R , that is either represented by the set of all points that belong in Q, 2 or by the set of points that lie on its boundary F. In R these are solid and surface representations of 3 a body, respectively. The selection of a surface or a solid representation greatly influences achiev- (a) (b) Figure 2.1: Representations of a continuous body in R (b) a solid. 2 as (a) a surface, or able deformation model behaviour. This is shown in Figure 2.2, which illustrates the deformation of elastic materials using: (a) a surface model, and (b) a solid model. Surface-based models are clearly (a) (b) Figure 2.2: The deformation of non-rigid structures using: (a) a surface representation, and (b) a solid representation. 2.2 Models and Simulation Methods for Deformable Materials 20 unable to characterise volumetric behaviour such as incompressibility or the bending of thin plates, leading to invalid global behaviour. While their representations are usually far more compact, these models are generally not suitable for the physical modelling of non-rigid solids, and are only useful under special circumstances, e.g., the simulation of thin membranes and highly compressible bodies. The continuous nature of physical materials poses a problem for the numerical simulation of deformations by finite precision computing hardware. In practice, a non-rigid body is approximated by a discrete representation Cl' that can be represented by a finite set of parameters. The remainder of this section describes a variety of representations that have been selected for the simulation of non-rigid bodies. 2.2.2 I m p l i c i t Surface M o d e l s Splines and super-quadratics are commonly used to represent geometric curves and surfaces, as shown in Figure 2.3. In essence, these models exploit parametric functions to describe complex Figure 2.3: A curved boundary described by a spline. geometry. Just as an ellipsoid can be used to approximate the shape of an egg, the combination of several parametric curves and surfaces can be used to describe arbitrarily complex shapes with a finite set of parameters, e.g., control points in the case of splines, polynomial coefficients for quadratics and superquadratics, deformation polynomials for free-form deformations [59,85], and modal analysis [60]. The primary drawback of this approach is that it is difficult to express the mechanics of physical bodies in terms of these abstracted parameters. As a result, material defer- 2.2 Models and Simulation Methods for Deformable Materials 21 mations are often "tuned" by making manual adjustments to the model parameters. In addition it is difficult to use this approach for volumetric model representations. Nevertheless, because relatively few parameters are required to describe complex shapes, computations such as collision detection can be performed extremely efficiently. 2.2.3 ChainMail The Chainmail algorithm represents a non-rigid body as a set of linked elements that resemble a chain, as shown in Figure 2.4 [128,129]. The motion of a single mesh node is absorbed by the maximally stretched maximally compressed Figure 2.4: A one dimensional chainmail structure undergoing deformation. link interface between adjacent nodes. If the link between two nodes is stretched or compressed beyond predefined limits, then displacements are transferred to neighbouring nodes. In this way, small displacements of a single node within a relaxed medium results in local deformations, while displacements within a fully stressed medium cause global motion. This approach is computationally efficient and accommodates simple mesh topology changes; however, it is not clear how physically realistic deformation behaviour can be modelled. 2.2.4 Mass-Spring Models Mass-spring models attempt to approximate the elastic behaviour of continuous bodies by means of a finite set of mass nodes connected by a network of springs, as shown in Figures 2.5(a) and 2.5(b) [32, 33, 61, 68,103,110,130]. The network of springs, each obeying Hooke's law, leads to a set of linear equations describing the relationship between the applied force and displacement at each 2.2 Models and Simulation Methods for Deformable Materials 22 r' (a) (b) Figure 2.5: A mass-spring representation of region Q. network node n: Pn Pi Pn ~ Pi n (2.1) M where f„ is the internal elastic force at vertex n, P is the position of vertex n, N(n) is the set of n vertices adjacent to vertex n, and is the rest length of the spring that connects P and Pj. n The advantages of this formulation are that: • the network of springs and masses is easily implemented, • representation as a network of connected nodes is straightforward and consistent with the data structures used in graphics rendering, • it can be used for both static and dynamic simulations, and • geometric changes due to cutting, tearing, fracturing and suturing are accommodated by simply adding or removing connections between mesh vertices. Several disadvantages also exist: • Model behaviour is particularly sensitive to the network topology, which can result in an under constrained system that exhibits local minima (multiple rest positions—see Figure 2.6), or an over constrained system that tends to decrease the range of deformation. 2.2 Models and Simulation Methods for Deformable Materials 23 Figure 2.6: A n under constrained spring network, with multiple equilibrium configurations. The range of achievable mesh stiffnesses is limited by node density and the simulation sample time. Considering a single mesh spring, Shannon-Nyquist constraints limit its maximum stiffness to: k « ^ c > (2-2) beyond which the system response may become unstable, where m is the node mass and At is the simulation time step. In fact, undesirable network behaviour can occur at far lower spring stiffnesses within a mass-spring-damper system [131]. This is an issue primarily when explicit integration methods are used. In general, stability problems can be ameliorated by using semi-implicit integration schemes. Stability issues can also emerge with changes in mesh topology that may occur while cutting or suturing, where vertices are to be added to the system. Figure 2.7 shows that the introduction of a new node to the mesh model requires the re-computation of spring constants and masses within the vicinity of the addition. As vertices are added to the system, total body mass is preserved by reducing all or some of the vertex masses, leading to inconsistent dynamic behaviour, and further constraining spring stiffness by Equation (2.2). A n increase in the number of mesh springs requires an increase i n spring stiffness i n order to maintain constant global elastic behaviour (e.g., max^i,/^) > 2k in Figure 2.7). Therefore, changes to the mesh topology can destabilise its dynamic performance if the simulation time step At remains constant. In addition, the introduction of new mesh nodes increases computation time, which may decrease the simulation sampling rate as a consequence. Therefore significant topological changes are difficult to implement for dynamic mass-spring models. 2.2 Models and Simulation Methods for Deformable Materials 24 m (a) (b) Figure 2.7: The addition of mesh nodes requires the re-computation of adjacent spring constants and masses. • Mass-spring models dp not rely upon continuum mechanics, therefore their behaviour cannot easily be compared with the results of biomechanical studies. For large deformations, spring models do not behave like linear elastic materials; however, with the assumption of small deformations, spring parameters can be identified in order to model the desired behaviour, e.g., Radetzky et al. use a fuzzy system to initialise model parameters from coarse descriptions of deformation behaviour, followed by back-propagation from experimental data for parameter refinement [33,132]. 2.2.5 Continuum Models The constitutive equations that characterise metals, polymers, fluids, gasses, magnetic and electromagnetic fields are well established. They are able to describe the dynamics of many continuous materials, often with extremely high fidelity. The Finite Element and Boundary Element methods of analysis are popular and effective tools for approximating solutions to these equations, and are described in this section. 2.2.5.1 The Finite Element Method The continuous constitutive equations that describe the physical behaviour of a solid body can be analysed by dividing the body Q into a set of discrete sub-domains called elements ( 0 ) , as shown e in Figure 2.8. Constitutive behaviour can be solved analytically in each element, and combined to estimate the global solution. This method of analysis is detailed in Section 3.5 and Appendix C . The advantages of the Finite Element Method, in the context of soft tissue modelling, are: • Finite element models are based upon continuum mechanics and the physics of materials; therefore, it is easy to incorporate the results of biomechanical studies. For example, the 2.2 Models and Simulation Methods for Deformable Materials 25 Figure 2.8: The continuous domain 0, is divided into a finite number of discrete elements fl . e Young's Modulus and Poisson's Ratio characterise linear elastic materials, are used to construct the system matrix directly, regardless of the mesh topology. • Dynamic behaviours are easily accommodated (e.g., by using a Newtonian formulation), and integrated using semi-implicit or explicit methods. The simulation time step At does not restrict the model stiffness when using a semi-implicit integration technique, unlike the massspring model approach. Therefore, re-meshing due to cutting, tearing, fracture or suturing does not affect simulation stability. • The wealth of existing knowledge regarding both linear and non-linear continuum mechanics and Finite Element Analysis methods can be employed for the simulation of complex tissue behaviours. For example, such methods are also adept at describing fluid behaviours, which is difficult for many of the methods discussed earlier. Furthermore, the limits of performance and accuracy are well known. • Representation as a network of connected nodes is consistent with the data structures used in graphics rendering. The disadvantages of the finite element approach are that: • It is extremely difficult to attain real-time (i.e., >30Hz update rates) performance for significant models using current off-the-shelf computing hardware, because of the extreme computational complexity involved in solving a system containing more than a few hundred mesh 2.2 Models and Simulation Methods for Deformable Materials 26 nodes. A few significant optimisation methods exist; however, many of these are only applicable to linear problems (see Section 2.2.5.2). • Its implementation is not as straightforward as for mass-spring models. The finite element approach was adopted for several of the projects described in Section 2.1 [55,57, 66,70, 71, 78,79,86,88,97,98,105-109,124,133]. 2.2.5.2 Condensation and the Boundary Element Method In some cases only the motion of visible surface nodes is of interest, while that of internal nodes is irrelevant. While this does not mean that a volumetric representation can be dispensed with completely, due to the desire to describe realistic global behaviour, it does mean that the simulation need not solve for all the internal states explicitly. Bro-Nielsen [134] showed that the linear system of equations derived from a Finite Element discretisation of linear elasticity can be simplified by means of condensation, or the elimination of internal nodes from the system equations, leading to smaller, denser models. This strategy results in a system with the same dimensionality as a Finite Element surface model, but with precisely the same behaviour as the original volumetric model, with reported decreases in computation times by at least a factor of 100. Condensation can only be applied to linear systems and is not feasible when topological changes are likely. James and Pai achieved real-time deformation performance for graphics rendering by means of a similar approach effected in the continuous domain prior to discretisation [62]. They expressed Navier's equation using a Boundary Integral Formulation, and discretised it using the Boundary Element Method. The set of fundamental solutions were computed offline, prior to simulation, with pre-computation costs similar to those of the condensation technique (both methods essentially compute the Green's Functions). This approach permits the real-time interactive simulation of linear elastic bodies with a fixed mesh topology. De and Srinivasan use an alternative approach for reducing the solid modelling problem to a surface modelling problem for compliant organ models (e.g., stomach, gall bladder). They treat the solid body as a thin membrane enclosing a fluid that is characterised by a Bulk Modulus, hence reducing computational complexity significantly [135]. 2.3 The Biomechanics of Tissue 2.3 27 T h e B i o m e c h a n i c s of T i s s u e Although the behaviour of structures and materials in engineering can be relatively accurately described by abstracted models such as those for Hookean elastic solids, non-viscous fluids and Newtonian viscous fluids, most biological materials exhibit far more complex dynamic behaviour [136]. For example, blood is better characterised as a visco-elastic material that exhibits hysteresis, relaxation and creep under dynamic loading, than as a Newtonian fluid. Mucus, sputum and synovial fluids present quasi-elastic behaviour, along with most bioviscoelastic solids, including abductin, resilin, elastin and collagen. Collagen is the basic structural element for soft and hard tissues, but is present in a large variety of structural forms within different tissues and organs. Skeletal muscle is an anisotropic viscoelastic material. Bone exhibits anisotropic linear elastic behavior, and cartilage behaves in a visco-elastic manner. fn summary, biological tissues exhibit complex behaviours, particularly when combined within heterogeneous bodies. For example Figure 2.9 illustrates the non-linear viscoelastic behaviour of a membrane sample (graph reprinted from [136]). Much ofthe contemporary work on tissue modelling and simulation has focused on the elastic component of deformation, particularly linear elasticity. Work to extend these models is ongoing. 2.4 Needle Insertion M o d e l s and M o t i o n P l a n n i n g To date, there has been very little work on the modelling and simulation of the physical mechanics of needle insertion. Those systems that do exist are primarily aimed at spinal procedures that involve stiff tissues and single axis motion. For example, simulations of needle insertions through the ligaments ofthe spine—for teaching epidural lumbar puncture [7,17-24] and spine biopsy [25]— have been reported. The majority of these systems have been built using ad-hoc models. Many employ lookup tables to determine needle force feedback versus depth for 1-DOF interaction [f3,19], while others use elastic, viscous and viscoelastic models to approximate needle driving impedance along the axial direction [17,18]. Soft tissue interaction, deformation and needle flexibility have not been considered. Needle insertion forces have been measured in gelatine [20], ex vivo porcine, human and bovine samples [17, 39]; however, only the resultant forces acting at the proximal ends of the needles were 2.4 Needle Insertion Models and Motion Planning 0 .02 .04 ' i .06 28 i i .08 i i .10 i i .12 i 1 .14 EXTENSION, INCHES II - 2 RABBIT MES. l_o 0.58' L 1.34" STRAIN RATE ± 0.1 in./min. p n Figure 2.9: Typical loading-unloading curves of rabbit mesentery tissue (supportive membrane). measured. Simone and Okamura attempted to isolate stiffness, friction and cutting forces from the driving forces measured during insertions into ex vivo bovine liver tissue [39]. Their results varied significantly from one insertion to the next, due to the inhomogeneous nature of the liver, and force components could not be isolated from a single insertion. Segmented tissue models were not considered; therefore, it is difficult to generalise or to validate their results. Kataoka et al. investigated needle deflections during linear insertions into soft muscle tissue ex vivo [137], and later measured tip and friction forces by means of an instrumentation approach [138]. They used an instrumented needle comprising an inner needle shaft and an outer sleeve that extends to just behind the inner needle tip. Axial tip forces were measured by the inner shaft, and the remaining forces were measured by the outer sleeve. Two insertions into a canine prostate were reported. While both the tip and shaft forces were measured simultaneously, tissue structures were not segmented; therefore, the result is also difficult to generalise. Tissue deformation was not measured. Alterovitz et al. recently presented a simulation system for analysing the sensitivity of seed placement errors to relative insertion depth, insertion height, needle sharpness, friction, velocity 2.4 Needle Insertion Models and Motion Planning 29 and tissue parameters (e.g., stiffness) during prostate brachytherapy procedures [40]. They used a linear elastic model, discretised by a Finite Element approach similar to that presented i n this work, as well as a simple linear needle insertion model. Tissue and insertion model parameters were varied in order to study placement error sensitivity, and were loosely based on the measurements reported in [39]. Needle deflection was not considered. Nienhuys extended the two dimensional simulation models presented in this thesis to three-dimensional neo-Hookean models, using an element subdivision approach (based on edge-bisection) to ameliorate the effects of added computational complexity. They consider homogeneous models and rigid needles. Needle motion planning i n the presence of deformable tissues and obstacles has not been considered before. A variety of algorithms are available for solving motion planning problems involving rigid bodies and articulated rigid bodies with kinematic and dynamic constraints (see [f 39] for a review of classical techniques). For needle manipulation, the control parameters are defined to be the linear and angular velocities that are applied at the needle base, as a function of time. The system configuration q lies i n the configuration space C , which usually contains a goal configuration, or a family of goal configurations. Obstacles are mapped to regions in C called C-obstacles, denoted CB . The motion planning problem is solved by finding a satisfactory path from the initial configuration qo to a goal configuration q { in C — CB. goa Cell decomposition search methods involve discretisation of the configuration and parameter spaces, followed by a search for the optimal path that connects go and q i i goa n C, subject to a cost function, defined constraints, and CB. Prior motion planning work has focused on efficient discretisation and search methods in C. Because the entire C-space is considered, global solutions can be sought; however, the major disadvantage of such approaches is that the search spaces can become extremely large or unbounded, even for modest problems with multi-dimensional configuration and control parameter spaces. The "tractor-trailer" is a classic example of motion planning with nonholonomic kinematic constraints [139], and is depicted i n Figure 2.f 0(a). The configuration of the tractor is specified using its orientation angle 6\, and the coordinates (x\, y\) of point R, which lies in the middle of the tractor's rear axle. The trailer is coupled to point R by a joint that introduces one additional state variable, 62. The control parameters are the velocity v of the tractor along its axial direction, and steering angle <f>. For constant 0 and v, it is possible to solve for the system state {x\, X2, 61, 62} after some time interval St. The goal of the tractor-trailer problem is to find 2.4 Needle Insertion Models and Motion Planning 30 Figure 2.10: (a) A steered four-wheel vehicle with articulated trailer, and (b) an analogous needle approximation. a feasible free-space path that will take the tractor and its trailer from an initial configuration to a final goal configuration, subject to the non-holonomic constraints that are introduced by the steering mechanism. Barraquand and Latombe [140] solved this problem by decomposing the configuration space into discrete cells, and searching a directed graph with nodes representing cells in C. fn the graph, cells are adjacent if a feasible free-space connecting path exists in C—arcs are constructed by discretising the control parameters (e.g., v 6 {-vo, «o} and <f> € {-(f>max, 0, <j> ax})- They use m an A* search with a suitable cost function (e.g., minimum reversals, minimum path length) to find the best graph traversal, and hence motion plan. An analogous approximation ofthe needle motion planning problem is illustrated in Figure 2.10(b), in which the needle is represented by a series of "articulated trailers". Kinematic constraints include those imposed by needle flexibility, as well as by the tissue. Non-holonomic constraints appear when the Needle Manipulation Jacobian J becomes singular; however, it is difficult to show that the resulting needle insertion constraints are indeed non-holonomic, at least using the notion of integrability [139]. A C-space search—similar to that outlined for the tractor-trailer problem—may be applied to the needle insertion problem; however, even coarse discretisations of the needle and the tissue result in extremely large configuration spaces, since there are many more states than in the tractor-trailer problem. Potential field methods treat the system as a particle under the influence of an artificial potential field in the configuration space. The potential field is typically defined over C as the sum of 2.4 Needle Insertion Models and Motion Planning 31 an attractive potential that pulls the system toward a goal configuration, and a repulsive potential that pushes the system away from constraints and obstacles. Potentialfieldmethods are often implemented as "local methods", where the potential functions are defined such that their value at any point in the configuration space does not depend upon structures in CB outside of the local neighbourhood. This approach has been used for on-line planning [141]. The method resembles a gradient-descent optimisation method, and as such is sensitive to local minima, and cannot guarantee a path to the goal. Nevertheless, such methods are computationally efficient, particularly for problems with high-dimensional configuration spaces. Several methods for dealing with local minima are discussed in [139]. Several soft tissue modelling applications and methods were presented in this chapter, culminating in a review of prior work that involves needle insertion models, simulation and planning. The following chapter details a system for measuring the interaction between needles and soft materials during insertions, as well as the characterisation of tissue phantom models. Chapter 3 Needle Insertion Experiments and Model Characterisation A system for experimentally observing and characterising needle insertions is described in this chapter. The materials and methods required for each block shown in Figure 3.1 are discussed. Experimental measurements are used to identify parameters for a tissue phantom model that forms the basis for subsequent chapters. Instrumented Probe/Needle Robotic Manipulator Tissue Phantom measure force apply trajectory probe/penetrate Image-based Marker Tracking Algorithm track marker motion Figure 3.1: Needle insertion and probing - experimental procedure. A 3-DOF, planar robotic manipulator [142] is used to drive an end-effector (a blunt probe or a needle) into a thin rectangular slab of soft elastic material that acts as a tissue phantom. Both tool force and tissue phantom deformation are sampled during the manipulator trajectory, to be used in model identification and verification, as described in Section 1.2.1. The equipment is shown in Figure 3.2, as it is configured during experiments. 32 The soft tissue phantom is molded from 33 polyvinyl chloride ( P V C ) and liquid plasticizer - a phthalate ester (Super Soft Liquid Plastic, M - F Manufacturing, T X ) . By varying the proportion of plasticizer added to the P V C , one can produce soft elastic phantoms with elasticities (Young's Modulus) ranging from below lOkPa (approximately that of breast tissue [3]) to well over lOOkPa, a range that covers many of the body's soft elastic tissues [136]. The tissue phantom lies upon a horizontal platform and is fixed (attached to the platform) along its rear edge, as shown in Figures 3.2, 3.3 and 3.8. It is mounted directly below a colour C M O S camera, and is marked with a grid of black dots, each approximately 2mm in diameter, that can be tracked by the camera as they move. A silicone lubricant allows the phantom to slide on the platform, parallel to the camera image plane and needle insertion axis, with very little friction. This chapter describes the actuation and sensing systems that are used to measure forces and deformations that occur during probing and needle insertion experiments. characterisation, based on experimental results, are also discussed. Model selection and 3.1 Motion and Force Sensing Figure 3.3: A tissue phantom with markers painted on its top surface. Figure 3.4: (a) Instrumented probe, and (b) Instrumented 17 gauge Tuohy epidural needle. 34 3.1 Motion and Force Sensing 3.1 35 M o t i o n and Force Sensing A blunt probing tool and a needle are instrumented with a 6-axis force/torque sensor (ATI Nano-17 SI-12-0.12), and can be attached to the 3-DOF planar manipulator as shown in Figure 3.4. The sensor has force and torque resolutions of f2.5mN and 0.0625Nmm force/torque resolution, respectively. The manipulator mechanism moves the end-effector along probing or insertion trajectories by means of a redundantly actuated pantograph-based mechanism that has been presented i n the context of haptics (see [142]). It has a translational workspace greater than 100mm x 100mm, with unlimited rotation and a spatial resolution of approximately 40/um and 0.05°. W i t h the addition of motor gearheads, the manipulator is capable of exerting peak forces greater than 30N, and 12N continuously. Mechanism and controller details are provided in Appendix A . 3.2 Image-based Deformation Estimation Images from a single C M O S camera are processed to estimate the 2-D deformation of the top surface of the tissue phantom. The camera is mounted as shown in Figure 3.2, and is interfaced with a P C by means of a Camera Link Interface. Images of up to 2056 x 1544 pixels are captured by the camera and streamed to the P C at a maximum rate of fOO megapixels per second, or 30 frames per second at full resolution. Typical images are shown in Figure 3.5. Camera specifications, calibration details and intrinsic camera parameters are provided in Appendix B . After calibration, the standard deviation in pixel reprojection error is approximately 0.4 pixels along each axis. For the given geometry, one pixel spacing corresponds to a distance of 0.079mm at the top surface of the tissue phantom. Markers are arranged on a regular 2f x 2f grid, spaced at 5mm intervals, and are each tracked from frame to frame throughout a sequence of images taken during probing and insertion experiments. A correlation-based template matching scheme is used to locate and track the material markers over the entire sequence of image frames, as detailed in Section 3.3. This approach provides robust tracking despite non-uniform lighting conditions, specular reflections on surface markers, and non-stationary needle shadows that are visible through the transparent phantom. Tracking results for two image frames are given in Figure 3.5. The camera lens introduces optical aberrations such as radial and tangential distortions that significantly affect marker position estimates. Distortions of more than 80 pixels are present, as 3.2 Image-based Deformation Estimation 36 (a) (b) Figure 3.5: Image-based marker tracking during (a) boundary probing and (b) needle insertion. Markers are located by a correlation-based template matching scheme (see Section 3.3). The crosses overlayed onto the images indicate the marker coordinates estimated by the algorithm. measured node position actual node position S00 rectified node position actual node position 4501 •100 350 | 300 ) 250 200 100 200 300 400 X [pixels] (a) 500 600 700 400 X [pixels] (b) Figure 3.6: (a) Actual and imaged calibration marker locations, (b) Actual and rectified calibration marker locations (mean error of approximately 1 pixel = 0.079mm). 3.2 Image-based Deformation Estimation 37 is evident in Figure 3.5 (notice that straight edges of the platform become curved). Much of this distortion is compensated by inverting the standard radial-decentering distortion model [143], based on the intrinsic camera parameters that are identified as described in Appendix B . This is shown using a calibration grid, as illustrated in Figure 3.6. After rectification, optical aberrations account for systematic tracking errors of less than one pixel (<0.1% of tissue phantom width), which has very little impact on estimated marker displacements from frame to frame. Tissue motion at the markers/nodes is sampled during both boundary probing and needle insertion. Sample node displacements, measured from an initial unperturbed state, are shown in Figure 3.7. In (a), a point force is applied to the tissue boundary by the probe, while in (b) the needle intercepts several internal nodes as it penetrates the material. Such measurements are used to identify material models and contact forces in subsequent chapters. 0.06 0.06 0.04 0.04 0.02 0.02 V • • • V t k * • * ^ •» f •. o x * * " • « • » r 0.04 ' r r f f • t r f f J» -0.04 -0.06 -0.04 -0.02 0 0.02 0.06 s *— *—• <L *• • •< ** .• . * 0.02 -0.02 -0.06 0 •* *» .•• •i. •• - •» f -0.06 -0.04 —• —• -0.02 Y[m] Y[m] (a) (b) 0 0.02 Figure 3.7: (a) Marker displacements due to probing, (b) Marker displacements due to needle insertion. The camera system and marker tracking algorithm measure the motion of the phantom's upper surface; however, the sample has finite thickness so that needle insertions and boundary probes occur some distance below this surface, as illustrated in Figure 3.8. Therefore, the camera images do not measure motion in the needle plane. The slab of phantom material is made as thin as possible in order to minimise out-of-plane deformations that cause the motion in the insertion plane and the motion at the top surface of the phantom to differ. A thickness of one tenth of the phantom width 3.3 Tissue Phantom Marker Tracking 38 produces small measurement errors while preventing buckling. Fiducial markers could be placed marked tissue surface tissue thickness lubricated platform Figure 3.8: A profile view of the tissue phantom, illustrating deformation perpendicular to the imaging plane. or injected into the insertion plane [f44]; however, this complicates imaging and may affect local material properties. 3.3 Tissue Phantom Marker Tracking Markers are applied using a paint mask; however, their size, colour density and shape can vary over the top surface of the tissue phantom. During experiments, the camera's frame integration-time was set to approximately 0.6ms to minimise motion blur and distortion artifacts due to the rolling shutter. Such effects are of concern at high needle insertion rates approaching 9mm/s. Halogen flood lights, powered by a direct current (DC) source, illuminate the phantom from above and below the platform (both the platform and the phantom are semi-transparent). Approximately 300W of lighting is required to sufficiently illuminate the scene for the short integration time. Non-uniform lighting conditions cause intensity variations and specular reflections within the image scene. In addition, the needle is visible through the transparent phantoms during insertion, as is evident in Figure 3.9. A correlation-based template matching technique was used to locate individual markers, and was found to be quite robust under the given experimental conditions. Marker tracking is initialised using the first image frame of an experiment sequence, and proceeds as follows: 1. Capture a sequence of images during a probing or needle insertion trajectory. 3.3 Tissue Phantom Marker Tracking 39 Figure 3.9: A typical image taken during a needle insertion experiment. Markers are clearly visible, but vary in size and quality throughout the image. 2. Ask the user to identify the approximate locations of the four corner markers in the f i r s t image (the tissue phantom i s i n i t s undeformed state). 3. Estimate the coordinates of a l l markers using the known undeformed grid geometry. 4. Convert the RGB pixel values to intensity values by choosing the minimum colour component at each p i x e l . 5. Erode and then dilate the image to f i l l the markers. "holes" due to specular reflections on This is a closing operation. 6. For each marker, cross-correlate a binary template (shown i n Figure 3.10) with a 41 x 41 p i x e l neighbourhood around the estimated coordinates. The template is cross-correlated with the complement of the image in this neighbourhood. 7. Find the peak correlation coefficient that indicates the centre of the marker. Interpolate the correlation coefficients using cubic polynomials along each image axis to locate the peak on a sub-pixel scale. 8. Repeat steps 6 and 7 u n t i l the marker is in the centre of the search neighbourhood. 9. Load the next image in the sequence, and use the estimated marker coordinates from the previous frame to locate i n i t i a l search neighbourhoods. 10. Repeat steps 4 to 9 u n t i l the last image is reached. 3.3 Tissue Phantom Marker Tracking 40 Erosion, dilation and template matching operations are discussed in most standard textbooks on digital image processing and computer vision (e.g., [145] and [146]). Figure 3.11 shows the results Figure 3.10: A binary template is cross-correlated with image regions in order to locate markers. at several stages during the localisation of three different markers taken from the frame shown in Figure 3.9. Note that despite significant differences in colour density, shape and lighting, the centre of each marker is located using the same template matching algorithm - marker 1 has low contrast with the background, marker 2 has a 'hole' due to specular reflection and marker 3 is directly above the needle. The full set of tracked markers is shown in Figure 3.12. Figure 3.11: Tracking results for three different markers after: a) raw image capture, b) RGB to Grayscale conversion, c) erosion and dilation, d) cross-correlation with the template, and e) final coordinate estimate. 3.3 Tissue Phantom Marker Tracking Figure 3.12: The tracking algorithm locates 441 markers within each frame. Estimated marker coordinates are superimposed as crosses on the camera image. 42 3.4 Compression Testing of the Tissue Phantom 3.4 43 Compression Testing of the Tissue P h a n t o m The static stress-strain characteristics of a 30 x 30 x 10mm sample of the tissue phantom were measured by means of a parallel-plate unconfined compression test, in order to facilitate modelling. The result is shown in Figure 3.13. The stress-strain curve* is close to linear up to strains of 15-20%. The relationship is non-linear for strains greater than 20%, due to a combination of material and geometric non-linearities (e.g., local strains are not uniform throughout the sample). The initial slope (Young's Modulus) is of interest, and is approximately 27kPa. Figure 3.13: Compressive loading of a 30mm x 30mm x 10mm material sample. The measured stress-strain relationship is shown on the right. 3.5 A P l a n a r L i n e a r E l a s t o s t a t i c M o d e l for Tissue P h a n t o m s From the review of applications and methods in Sections 2.1 and 2.2, it is clear that tissue deformation is complex, and is still the subject of much research. In general, tissue modelling is complex because of inhomogeneous, non-linear, anisotropic elastic and viscous behaviour. Some models approximate elastic behaviour and are suitable for predicting the deformation of skin and muscle fascia, while other models approximate viscous or viscoelastic behaviour and can estimate the motion of fluids and adipose tissue [17]. *Cauchy Stress is used with the assumption of incompressibility. Area measurements taken during compression confirm that the material is close to incompressible. 3.5 A Planar Linear Elastostatic Model for Tissue Phantoms 44 As a first approximation, this study focuses on homogeneous, linear elastostatic models. The soft plastic material selected for the tissue phantoms does indeed exhibit linear elastic behaviour over a large strain range, as indicated by the stress-strain curve in Figure 3.13. Higher-order dynamic effects (i.e., mass and damping) are neglected, as it is assumed that body accelerations and velocities are relatively small during typical needle insertions. The phantom is discretised in two dimensions, as shown in Figure 3.14. For an elastic continuum, the total strain 0.1 0.08 0.06 & >- 0.04 0.02 0 -0.06 -0.04 0 X[m] -0.02 0.02 0.04 0.06 Figure 3.14: The continuous domain is divided into a finite number of discrete elements Q . Element vertices constitute nodes in the mesh. e energy E ain over a solid body Q, as a function of stress cr and engineering strain e, is given by: stT Estrain = \ I e (x)o-(x) dx , r 1 (3.1) JQ and is minimised at static equilibrium. W i t h the assumption of a linear relationship between cr and e, and after discretising (3.1) using linear shape functions, the static equilibrium condition for a single element under load is expressed as: SE =[ e B CB u eT e e dy.-f e = 0 , (3.2) 3.6 Model Parameter Identification 45 where each element e reaches its equilibrium state when the first variation of the energy functional 5E vanishes. (B e C B ) characterises the elastic behaviour of element £l , u contains the dis- eT e e e placement vectors and f contains the force vectors for those mesh nodes that constitute element e e (the applied forces f are concentrated at mesh nodes). Over the entire set of elements in the body, e this leads to a set of 2n linear equations: K (2«x2n)M = £ • (3-3) If Green-Lagrange strain is used, then the relationship becomes non-linear: (3.4) K(u)u = f. Equations (3.3) and (3.4) constitute the function Q that was introduced in Section 1.2.1. The constitutive equations for linear elastic materials, as well as the Finite Element Method of discretisation, are presented in detail in Appendix C . 3.6 Model Parameter Identification In biomechanics studies, tissue model parameters such as Young's modulus (E) and the Poisson ratio (u) are typically determined by testing small homogeneous tissue samples using rheometers and similar materials testing equipment [136]. Measurements of elasticity and viscoelasticity are usually obtained by loading the tissue sample between two parallel plates, in a fashion similar to that shown in Figure 3.13. In this section, model parameters for tissue phantoms are identified by a boundary probing approach, in which a constant point force is applied at the phantom's boundary. Applied forces and material deformations are measured and analysed in order to identify both the Young's modulus and the Poisson ratio that characterise the elastic behaviour of the phantom. This approach allows for the measurement of non-homogeneous tissue samples. The set of model parameters 0 = {E u} that minimises the squared error between measured and predicted node positions is determined by the following minimisation: min (u-g(E,u,f)) .(u-g(E,v,i)) {E,w} T , (3.5) where u is the vector of node/marker displacements measured at the top surface of the tissue 3.6 Model Parameter Identification 46 phantom, Q is the elastic material model that is a function of model parameters E and v, and f is the vector of external forces applied at model nodes (in this case, the probe is in contact with just one model node on the boundary, as shown in Figure 1.5). Function Q can be a two dimensional plane stress model, or a three dimensional model; however, only the node displacement errors at the top surface of the phantom are available to create the objective function to be minimised in (3.5), as illustrated in Figure 3.15. Figure 3.15: The tissue phantom model can be represented in two or three dimensions; however, planar material displacement measurements are only available from the top surface of the phantom. The ANSYS® Finite Element analysis software was used to simulate the model, as well as to implement an optimisation loop that solves (3.5) by means of a first-order gradient descent optimisation algorithm. Model parameters are provided in Table 3.1. The least squares result Table 3.1: Summary of the 2-D and 3-D phantom model parameters used in ANSYS®. 2-D 3-D linear elastic, isotropic linear elastic, isotropic 100 x 100mm 100 x 100 x 11mm Material Young's Modulus: 25.3kPa 25.3kPa Material Poisson Ratio: 0.48 0.48 Number of Nodes: 441 (21 x 21) 2205 (21 x 21 x 5) Element Type: 4-node quadrilateral, 8-node brick Material Type: plane stress, thickness 11mm Strain: Green-Lagrange, Green-Lagrange, large strain analysis large strain analysis 3.6 Model Parameter 47 Identification produces consistent parameter estimates over a large range of initial values. Parameters identified over several probing forces are more consistent when the three dimensional phantom model is used. Measured marker positions and node positions estimated by the optimal model are superimposed in Figure 3.16 (for a typical homogeneous phantom). The 3-D model produces mean and maximum nodal displacement errors of 0.117mm and 0.547mm respectively, while the 2-D model produces mean and maximum nodal displacement errors of 0.147mm and 3.125mm. For the homogeneous tissue phantoms used in experiments, the optimal material parameters derived by three-dimensional analysis are E = 25.3kPa and v = 0.48. The Young's Modulus estimate is approximately 6% lower than that measured by the compression test in Section 3.4, due to changes in the material stiffness that took place during the period of several days between the boundary probing experiment and the compression test. The plasticizer tends to leak out of the P V C phantoms, so that they become stiffer with time. This has been confirmed in follow-up tests. 0.06 • x 0.06 measured node position model node position (2D) 0.04 0.02 -0.02 x x <i n * » * * X X * X X X -0.06 -0.04 0.02 * X x—px—x ' X X X x x X x X -0.02 Y[m] x x V X « J « < — * » I -0.02 x' x* * J ' X -0.06 0.04 x<^x X -0.04 measured node position model node position (3D) 1 X* -0.04 x f tl 0.02 -0.06 -0.06 -0.04 -0.02 Y[m] 0.02 Figure 3.16: Measured marker locations and estimated model node locations for a typical homogeneous tissue phantom are superimposed. A 2-D plane strain model (a) and a 3-D model (b) are compared. Error vectors are shown at each model node (scaled by 10). A two-layer inhomogeneous tissue phantom was tested in similar fashion, and material parameters identified from boundary probing experiments. Model elements conform to the tissue layers, 3.6 Model Parameter Identification 48 as shown in Figure 3.17. The model parameters for the outer, thinner layer were identified as E = 63kPa and v — 0.49. Parameters for the thicker layer were found to be E = f8kPa and v — 0.49. This parameter identification approach may also be applicable to multi-phase tissue £=6~3kPa v = 0.49 E= 18kPa v = 0.49 Figure 3.17: This phantom model has two layers of differing material stiffnesses, and is shown here with a point force applied at its boundary. samples and tissue phantoms, as well as for dynamic analysis. In this work, the inhomogeneous phantom samples have known segmentation/structure; however, in general this may not be the case. When the structure of the sample is unknown, the number of parameters to be identified is far greater, and issues of model excitation may be relevant. The following chapter presents the needle shaft forces that are estimated from needle insertion experiments, based on the tissue phantom model derived here. Chapter 4 Needle Shaft Force M o d e l The direct measurement of needle forces by an instrumentation technique is a challenging problem. A n alternative, indirect method for estimating needle shaft forces, using measured experiment data, was described in Section 1.2.1. A constitutive model of the tissue phantom was selected and parameterised using elastic deformations measured during probing experiments (see Chapter 2 and Appendix C ) , resulting in the following relationships: u = g(o,i) f = g~\e,-a), (4.1) where Q represents a model that relates local material displacements u and local external forces f, and is parameterised by a set of constant material parameters 6 that have already been identified. G is a non-linear function o f f if the quadratic strain terms are considered (see Appendix C ) . The inverse model Q~ exists; therefore, the set of local external forces f can be computed using the l displacements u that are observed during a needle insertion, where external forces are applied to the tissue phantom by the needle shaft. These are estimates of the actual forces applied to the phantom, since Q is an approximation. This chapter presents needle force estimates based on the linear elastostatic material model described in Section 3.5. Shaft force distributions are derived from displacement measurements taken during needle insertions at several driving velocities, as well as during insertions into a layered inhomogeneous tissue phantom. 49 4.1 Force Estimates from a 2-D Model 4.1 50 Force Estimates from a 2-D Model Clinical needle insertion velocities vary between 0.4mm and 10mm per second [20]. In the experiments reported here, a 17 gauge Tuohy epidural needle was inserted into the tissue phantom at several linear velocities from l m m / s to 9mm/s, while needle position, insertion force (at the needle base) and tissue phantom deformation were sampled as described in Chapter 3. It is assumed that each sample is a snapshot of static phantom deformation, for which the local strains are low enough (<15%) that we may consider linear elastic behaviour, with a Young's modulus of 25.3kPa and a Poisson ratio of 0.48. Needle flexion and base moments are small and are neglected. A planar material model (4.1) computes applied external force estimates—model parameters are summarised in Table 4.1. While the material model is quasi-static, second-order effects may be introduced at the needle-phantom contact. Insertions are observed at several different velocities in order to measure such dynamic effects. Table 4.1: Summary of the planar phantom model and model discretisation parameters (2-D model). Material Type Material Young's Modulus Material Poisson Ratio Number of Nodes Element Type Strain linear elastic, isotropic, 100 x 100mm 25.3kPa 0.48 441 (21 x 21) 4-node quadrilateral, plane stress with thickness of 11mm Green-Lagrange, large strain analysis Figure 4.1 shows the y-axis components of the nodal forces (the y-axis is the direction of insertion) after a needle displacement of approximately 70mm (from boundary intercept), during a l m m / s insertion. Significant external forces are found to occur at nodes lying along the needle shaft. Reaction forces are also evident at the rear edge of the model (y = 0.075m), where the phantom is fixed to the platform. A cross-section of the force plot, at the known location of the needle shaft, is illustrated in Figure 4.2. The material model provides force estimates along the needle shaft for each image frame, which are taken at 2mm needle displacement intervals, amounting to a total of 36 frames during a 70mm insertion. Shaft forces estimated at several instants during a l m m / s insertion are graphed in 4.1 Force Estimates from a 2-D Model 51 Y [ml Figure 4.1: External nodal forces computed by a numerical phantom model, based on measured nodal displacements. Figure 4.3. Each insertion is repeated several times, and four sets of nodal force estimates are averaged in order to reduce the effects of random measurement errors. The needle is inserted into the same phantom, but at different locations for each iteration. The distributions are integrated by summing nodal forces that lie along the needle shaft, and compared with the driving forces measured by the force/torque sensor at the base of the needle. This result is shown in Figure 4.4, where there is clearly a large discrepancy between measured and estimated driving force magnitudes. In order to investigate this, the ANSYS® Finite Element analysis software was used to compute a three dimensional model of the tissue phantom. 4.1 Force Estimates from a 2-D Model 52 Y|m] Figure 4.2: Nodal forces in the vicinity of the needle shaft are of interest. A cross-section of the force plot is taken along the row of markers that is coincident with the needle. Note that forces are also applied along the rear edge of the model; these are reaction forces from the platform constraint. Figure 4.3: Estimated y-axis forces at nodes that are coincident with the needle shaft, taken at various instants during a l m m / s insertion. The needle location is shown below each graph. 4.1 Force Estimates from a 2-D Model Figure 4.4: Comparison of the measured base needle driving force and the integrated needle force distribution (2-D model). 53 4.2 Force Estimates from 3-D Analysis 4.2 54 Force Estimates from 3-D Analysis The effects of out-of-plane deformations are investigated by computing a three dimensional model of the phantom, with parameters given in Table 4.2. The 3-D model is used to simulate phantom Table 4.2: Summary of phantom model and model discretisation parameters (3-D model). Material Type linear elastic, isotropic, 100 x 100 x 11mm Material Young's Modulus Material Poisson Ratio 25.3kPa 0.48 Number of Nodes 2205 (21 x 21 x 5) Element Type 8-node brick Strain Green-Lagrange, large strain analysis deformations due to known force loads, and then to check the force estimation technique that is based on a planar model. Hypothetical nodal force loads are applied along a row of mesh nodes in the medial plane of the 3-D model (i.e., the "insertion plane" shown in Figure 4.5), the appropriate displacement boundary conditions are applied, and the model is solved. The resulting nodal displacements at the top surface of the model are projected onto the plane stress model described in Table 4.1, which is used to solve for the nodal forces in much the same way as the experiment data was processed in Section 4.1. This is illustrated in Figure 4.6. node 1 ... node 10 The applied node 20 Figure 4.5: Nodal force loads are applied to a three-dimensional model of the tissue phantom in order to evaluate the 2-D analysis method. nodal force loads and estimated nodal forces are graphed in Figure 4.7. Two hypothetical load force profiles were applied, one with uniform force, and another that included a peak at the distal node, as shown. B y comparing the force estimates produced by the 2-D model, with the known force loads 4.2 Force Estimates from 3-D Analysis 55 Y[m] Figure 4.6: A three-dimensional phantom model is used to investigate the effect of out-of-plane deformations. applied to the 3-D model, it was determined that the 2-D analysis tends to underestimate actual needle shaft forces by approximately 26%. This is because out-of-plane deformations at the top surface of the phantom are ignored, and because the top surface of the phantom deforms less than the needle plane, as shown in Figure 4.8. Nodal force estimates based on a 3-D analysis technique, also shown in Figure 4.7, match the applied forces far more closely than those estimated using the planar model. The approach uses a 3-D model to correct nodal force estimates by the following least-squares minimisation: min (u-g- (E, l u,&)f,(u- Q~ (E,v,f )) x s , (4.2) where the squared difference between measured and estimated nodal displacements at the top surface of the phantom is minimised with respect to f^, the set of forces that occur at nodes along the needle shaft. Q(E,v,f ) s is a model of the tissue phantom, similar to that used in (3.5), as shown in Figure 4.6. Q is a non-linear function of f , since the quadratic strain terms are present; s therefore, a standard least squares algorithm is not applicable. ANSYS® was used to simulate the model Q, and to solve for load forces f_ using a first-order s gradient descent optimisation algorithm to solve (4.2). Nodal force estimates from the planar model (see Figure 4.7) were used to initialise f , and Equation (4.2) was solved with some manual s 4.2 Force Estimates from 3-D Analysis 5 10 15 20 ' 56 5 10 15 20 5 10 Node Node Node (a) (b) (c) 15 20 Figure 4.7: A three dimensional model of the phantom is used to compare 2-D and 3-D analysis methods, (a) Known nodal force loads are applied to the 3-D tissue phantom model, (b) Nodal forces are estimated by a 2-D model, using the nodal displacements computed as a result of loads applied to the 3-D model, (c) Estimated nodal forces based on the 3-D analysis. As expected, the applied nodal loads are recovered with very small errors. B y contrast, note the significant errors produced by the 2-D model. intervention to avoid local minima. Note that spurious node forces introduced by two dimensional analysis, such as the peak that appears close to the insertion point (node 20 in Figure 4.5), disappear after three dimensional analysis. Such artifacts appear in the two-dimensional analysis because the motion at the top surface of the phantom is not identical to the motion in the plane of the needle. In fact, as illustrated in Figure 4.8, motion at the top surface is not even planar. The corrected needle force distributions, one of which is shown in Figure 4.9, are integrated and once again compared with the measured driving forces, as shown in Figure 4.10, where there is very good agreement between experiment and model results. The shaft force distribution is found to be approximately constant at 3 4 N m , with a peak at the _ 1 most distal node of approximately 6 8 N m . Both the estimated and measured nodal displacements _ 1 at the top surface of the tissue phantom, resulting from a needle penetration of approximately 13 nodes, are plotted in Figure 4.11. The maximum and mean nodal displacement errors are within experimental error (see Chapter 3), and are 0.33mm and 0.18mm respectively. 57 4.2 Force Estimates from 3-D Analysis camera r 4 marked tissue surface tissue phantom lubricated platform Figure 4.8: Illustration of the cross-section of a tissue phantom during needle insertion. The top surface of the phantom does not remain flat and marker motion is not perfectly planar. Force per unit length 68N/m 34N/m Penetration depth Figure 4.9: Estimated force distribution along the needle shaft during an insertion at lmm/s. 4.2 Force Estimates from 3-D Analysis 58 1 — CD O 1 r i I i Measured Insertion Force Integrated Force Distribution (corrected) _ 1.5 o LL D) C > 0.5 1 1 1 0.01 0 0.02 0.03 0.04 0.05 0.06 Needle Displacement [m] Figure 4.10: Comparison of the measured base needle driving force and the integrated needle force distribution (3-D analysis). 0.06 x • measured node position model node position (3D) 0.04 0.02 X -0.02 X K K X K X X K * * * * X X X K X K X K X X * X X X X X K x X X X X X X X X X x X X X * v x x x x x x x x -0.04 -0.06 -0.08 x * X X X x x x x x x x x * x x x x x x x X -0.06 X X -0.04 X X -0.02 X x * x x x x x x x 0.02 Y[m] Figure 4.11: Estimated nodal displacements shown superimposed on actual measured nodes. The mean displacement error for 441 nodes is 0.18mm. 4.3 Force Distribution 4.3 Versus Driving Velocity 59 Force Distribution Versus D r i v i n g Velocity Needle insertions were repeated at several driving velocities from l m m / s to 9mm/s, and shaft force distributions were computed for each case. Y-axis nodal forces estimated for insertions at l m m / s , 3mm/s, 5mm/s, 7mm/s and 9mm/s, computed using the plane stress model, are presented in Figures 4.12 and 4.13. Nodal forces occurring along the needle axis, at several instances during each insertion, are graphed in Figures 4.14-4.17. Shaft force distributions computed by the threedimensional phantom model are illustrated in Figure 4.18, and the relationship between insertion velocity and shaft force "density" (force per unit length behind the force peak at the needle tip) is graphed in Figure 4.19. Shaft forces (those behind the tip force) increase with insertion velocity, as shown in Figure 4.19, while the force peak at the distal end of the needle appears to be somewhat independent of velocity, becoming less prominent at higher velocities. At a velocity of 9mm/s, this force peak is hardly discernible, and contributes less than 0.1% of the total driving force after a needle displacement of 70mm. The force estimated at the needle tip during the 7mm/s insertion (shown in Figure 4.18) is slightly higher than expected, and is considered to be an outlier result. Figure 4.12: External nodal forces computed by the planar phantom model, based on node displacements measured during insertions at 1, 3 and 5mm/s. 4.3 Force Distribution Versus Driving Velocity Figure 4.13: External nodal forces computed by the planar phantom model, based on node displacements measured during insertions at 7 and 9mm/s. 61 4.3 Force Distribution Versus Driving Velocity Figure 4.15: Estimated y-axis forces at nodes coincident with the needle shaft, taken at various instants during a 5mm/s insertion. 62 4.3 Force Distribution Versus Driving Velocity Figure 4.17: Estimated y-axis forces at nodes coincident with the needle shaft, taken at various instants during a 9mm/s insertion. 63 4.3 Force Distribution Versus Driving Velocity 64 . 9mm/s ' 7mm/s • 5mm/s • 3mm/s CO c 01 lmm/s CD a. Figure 4.18: Estimated force distributions along the needle shaft during insertions at 1, 3, 5, 7 and 9mm/s. Figure 4.19: Shaft force density versus insertion velocity, where the shaft forces are those uniform forces that occur behind the needle tip. A 3 degree polynomial is fitted to the datapoints. rd 4.4 Force Distributions in Inhomogeneous, Layered Phantoms 4.4 65 Force Distributions in Inhomogeneous, Layered Phantoms The force estimation method generalises to inhomogeneous materials. A two-part, layered phantom was used in experiments—dimensions and material parameters for the corresponding model are given in Table 4.3. Table 4.3: Summary of layered inhomogeneous phantom model parameters. Material Type linear elastic, isotropic, 100 x 100mm Layer 1 width 35mm Layer 2 width 65mm Layer 1: Material Young's Modulus 63kPa Layer 1: Material Poisson Ratio Layer 2: Material Young's Modulus Layer 2: Material Poisson Ratio 0.49 18kPa 0.49 Similarly to the homogeneous case, insertion experiments and analysis produce material force, shaft force and shaft force distribution estimates, as shown in Figures 4.20, 4.21 and 4.22 respectively. The relatively coarse mesh and planar analysis produce noisy results (Figures 4.20 and 4.21). Three dimensional analysis produces a clearer picture of the shaft forces; however, for narrower layers and finer structures, finer sampling will certainly be necessary. 4.5 Discussion The 2-D measurements and models have been very helpful in this first attempt to model needle insertions into soft tissue phantoms, as well as to establish the shape of the needle shaft force distribution. Errors due to the planar material model, planar deformation measurements and associated assumptions were limited by using a thin tissue phantom and by correcting force estimates using a 3-D model of the phantom. The needle force distribution that is estimated from experimental measurements, and used for simulation, indicates that axial forces between the needle and the tissue phantom are relatively uniform along the needle shaft. The force peak at the needle tip may be attributable to tissue cutting. In other work, researchers have attempted to separate individual physical contact phenomena such as cutting/fracture and friction, usually by inferring models from forces observed at the base of the 4.5 Discussion 66 Y[ml Figure 4.20: Nodal forces computed by a planar phantom model, based on nodal displacements measured during insertion into an inhomogeneous tissue phantom. needle. In contrast, this work aims to identify contact boundary conditions as a shaft force distribution, and to understand how this distribution changes during insertion. This force distribution may incorporate a number of physical contact phenomena such as cutting, elastic deformation and friction (static, kinetic and viscous) [17, 39]; however, it is not clear how to decouple these effects, or whether such separation is necessary. The force distribution itself may be useful for developing such models. Instrumentation-based approaches may also be used to separate force phenomena (e.g., [138]). Needle lubrication by blood and other fluids will affect contact friction, and can also be observed in the estimated force distributions. Further study of the effects of lubrication on needle shaft forces can be carried out using the approach outlined in this chapter. This chapter has presented a model-based approach for estimating the force distributions that occur along the needle shaft during insertions into tissue phantoms. In the following chapter, these force distributions are applied to the material model as boundary conditions in order to simulate needle insertion. While the needle insertion forces were estimated using three dimensional phantom models, the simulation methods presented in subsequent chapters are introduced using 2-D planar models. 67 4.5 Discussion -0.08 -0.04 0 0.04 -0.08 -0.04 0 0.04 -0.08 -0.04 0 0.04 -0.08 -0.04 0 0.04 -0.08 -0.04 0 0.04 -0.08 -0.04 0 0.04 y [m] Figure 4.21: Estimated y-axis forces at nodes coincident with the needle shaft, taken at various instants during insertion into an inhomogeneous tissue phantom. The needle location is shown below each graph. Force per unit length Figure 4.22: Shaft force distribution during insertion into a layered inhomogeneous phantom. Chapter 5 Physically-Based Numerical Simulation of Needle Insertion Mechanics This chapter describes a general, multi-degree-of-freedom simulation algorithm that simulates arbitrary trajectories of flexible needles as they are inserted into linear elastic tissue models. 5.1 A Simple Insertion Simulation A l g o r i t h m The estimated needle shaft force distributions (e.g., Figure 4.9) provide force boundary conditions that can be applied to a numerical soft tissue model for simulation purposes. For example, a simple simulation algorithm proceeds as follows: 1. Choose a tissue geometry represented by a mesh of elements and nodes, 2. Select a function Q to describe the relationship between external forces f applied to the nodes, and t h e i r resulting displacements u, 3. Consider a r i g i d needle shaft represented by a straight line segment, 4. Advance the needle shaft along a straight line into the numerical tissue model, 5. Identify tissue nodes that are close to, or coincident with the needle shaft and apply forces to them according to the shaft force d i s t r i b u t i o n , 6. Compute the r e s u l t i n g tissue node displacements u according to the function Q, 68 5.1 A Simple Insertion Simulation Algorithm 69 7. Update tissue node positions, go to step 4 and iterate. Figure 5.1 shows a simulated needle insertion to a depth of 70mm. The FEM-based linear elastostatic tissue model presented in Section 3.5 is used to describe tissue deformation, constituting the function Q. The model specifications are provided in Table 5.1. Simulated node displacements were found Table 5.1: Simulation model specifications. Material Type linear elastic, isotropic, 100 x 100mm Material Young's Modulus 25.3kPa Material Poisson Ratio 0.48 up to 441 (21 x 21) Number of Nodes Element Type 3-node triangle, plane stress with thickness of 11mm Strain Cauchy, linear strain analysis to closely match the tracked marker motion measured in experiments. Maximum and mean nodal displacement errors of 1.4mm and 0.6mm are observed, respectively. In a second simulation, the 0.34 0.32 0.3 0.28 0.26 0.24 -0.06 -0.04 -0.02 X[m] 0 Figure 5.1: Comparison between measured and simulated needle insertions, using the estimated force distribution. Left: a sample of the measured and simulated node positions from one half of the phantom. 5.1 A Simple Insertion Simulation Algorithm 70 needle force distribution is adjusted by increasing the tip force to approximately five times the shaft force, while keeping the total integrated force constant (the base driving force is the same as in Figure 5.1). Figure 5.2 shows that the simulated deformation differs markedly from the actual tissue phantom deformation, thus indicating the significance of the force distribution's shape. 0.34 0.32 0.3 E, >- 0.28 0.26 0.24 -0.06 -0.04 -0.02 X[m] 0 Figure 5.2: Comparison between measured and simulated needle insertions, with altered force distribution. Left: a sample of the measured and simulated node positions from one half of the phantom. This rather crude simulation algorithm works for the two examples given above, due to symmetry and the fact that the virtual needle is inserted along a set of colinear element edges. The algorithm is not suitable for general trajectories because: 1. the needle shaft might not always intercept control points or nodes on the discretised tissue model, 2. there is no contact model (nodal forces are simply applied according to the estimated shaft force distribution), 3. the needle may not move in a straight line or with uniform velocity, and 4. the needle may bend in practice. Each of these issues is considered in the remainder of this chapter. 5.2 Needle Contact and Tissue Model Boundary 5.2 Conditions 71 Needle Contact and Tissue M o d e l Boundary Conditions Tissue model nodes that are in contact with the needle are constrained as shown in Figure 5.3. If the needle is rigid, then the lateral position of the node is fixed along the ^-axis, constituting a displacement boundary condition. Either the node force or node displacement may be constrained along the needle shaft, depending upon its state of contact with the needle (i.e., sticking to the needle, or slipping) [50]. If the node is free to slide along the needle shaft, then a force boundary condition is applied along the ^/-axis, and a constant force consistent with the force distribution is applied to the slipping node. If it is in the stuck state, then the node is constrained to lie at a fixed point on the needle, namely along the Vaxis. The physical contact behaviour between needle and Figure 5.3: Mesh nodes lying along the needle are constrained along the ^-axis, and either slip or stick along the \/-axis. tissue is not known; therefore, a stick-slip model is adopted in order to be able to apply a needle force profile consistent with that determined experimentally (see Chapter 4), as well as to be able to accommodate direction changes and needle retraction. Afinite-stateautomaton is applied at each contact node (a tissue mesh node that is in contact with the needle) to simulate stick-slip behaviour (see Figure 5.4). Axial force thresholds are set according to the force distribution (e.g., Figure 4.9), such that the experimentally derived needle shaft forces are applied when simulated needle nodes are in the slip state. During simulated insertion, the boundary conditions at contact nodes will change according to their contact state. The system of equations in u = K _ 1 f is rearranged in order to reflect the resulting inhomogeneous collection of boundary conditions: u = K" f 1 -> x = (K- )'y , 1 (5.1) 5.2 Needle Contact and Tissue Model Boundary 72 axial force threshold exceeded force below threshold f j V_yv Conditions axial force opposes motion Node \ Stuck J axial needle motion is reversed Figure 5.4: A finite state automaton for stick-slip contact between tissue mesh nodes and the needle. Such an automation governs the contact behaviour at each needle node. where x and y are formed by exchanging elements between u and f, as illustrated in Figure 5.5. K' u Figure 5.5: A boundary-condition change at a contact node involves the exchange of displacement and force variables between the leftand right-hand sides of the tissue model equations. The x vector contains those variables that are to be solved, while y contains parameters (forces and displacements) that are specified as boundary constraints. For example, if a displacement constraint is applied at the i constraint, and the force vector is specified as a node, then the displacement vector Ji is to be solved. The vectors x and y become: Ji x uf u\ y = ft f! / f u\ fx Ji fV Ji V fx in (5.2) fV J i Similarly, if a displacement constraint is applied only along the ^-axis at the i th node, then the 5.3 Non-conforming Needle Trajectories is known, and vector 73 is to be solved. The vectors x and y become: ft l T 1 fx JI wt "1 uP "2 ••• fV Jl fx J2 ••• fi u i ••• «n „.x ft) h ••• fx Jn i u n u (5.3) fV In A single boundary condition change can be expressed as an inexpensive low-rank update: (K- )' = K - 1 where pi is the i th pivot of K the exception of their i th _ 1 ; Cj and r, are the i 1 th (5.4) Pi column and i th row of K coordinates, which are both set to (pi — f). _ 1 , respectively, with ( K ) ' is the new system -1 matrix. Displacement boundary conditions constrain tissue nodes to the needle geometry; whereas, force boundary conditions and thresholds impose the needle shaft force distribution estimated from physical experiments. The desired force distribution is integrated along the needle geometry, from tip to insertion point, and lumped at tissue nodes that are in contact with the shaft. 5.3 Non-conforming Needle Trajectories Figure 5.f illustrates a special case for which the tissue mesh exactly conforms to the needle geometry; however, in general the needle may not pass along element edges, or be coincident with mesh nodes. A simple remedy is to force the tissue model geometry to conform to the needle, such that mesh nodes "snap" to the needle shaft as it passes. This amounts to setting the lateral displacement constraints along the ^-axis in order to place the closest tissue nodes onto the needle shaft, as illustrated in Figure 5.6. When the needle tip begins to slip past the most distal contact node, the locations of neighbouring nodes are computed, such that subsequent intercepts can be identified. A precomputed, indexed connectivity list is used to determine the set of nodes Af that are connected to the most distal node, as shown. When the needle tip moves outside of the elements defined by the vertices in jV, then the closest neighbouring node is selected, and constrained to lie at the needle 5.3 Non-conforming Needle Trajectories 74 Figure 5.6: New intercept nodes are identified by searching within a small neighbourhood centred at the most distal needle node. tip, starting in the stuck state. This approach can produce large local tissue deformations and force discontinuities, particularly when the tissue mesh is coarse and the tissue stiffness is high. This section discusses two approaches for matching needle and tissue geometries using mesh subdivision and adaptation techniques. 5.3.1 Mesh Subdivision Mesh elements can be progressively subdivided in order to concentrate node density around the needle axis, as shown in Figure 5.7. Mesh elements can either be subdivided in a uniform manner that places new mesh nodes at predetermined locations within divided elements, or they can be subdivided in order to introduce new nodes that are exactly coincident with the needle shaft. Figure 5.7: Mesh subdivision. In the simulation example shown in Figure 5.8, a needle is inserted into the side of a tissue model that is rigidly fixed along one edge, in order to investigate the effect of tissue deformation on needle placement accuracy. The needle axis is initially coincident with a "virtual biopsy target". Elements in the vicinity of the needle are subdivided using edge bisection; therefore, new mesh nodes are introduced closer to the needle shaft, but not necessarily coincident with it. Mesh elements can be 5.4 General Non-linear Needle Trajectories 75 recursively subdivided until the mesh nodes adjacent to the needle shaft are sufficiently close (i.e., on the order of the needle diameter). In the example, elements are subdivided twice and thereafter the closest nodes are simply constrained to the needle shaft. The series of simulation frames shown in Figure 5.8 illustrates the needle's path through the tissue sample, and its failure to intercept the target (shown as a black disk). -0.06 -0.04 -0.02 0 X [m] 0.02 0.04 0.06 -0.06 -0.04 -0.02 0 XH 0.02 0.04 0.06 0.08 -0.06 -0.04 -0.02 O 0.02 0.04 0.06 0.08 X [m] Figure 5.8: Simulated needle intercept of a small target embedded within elastic tissue. 5.3.2 Mesh Adaptation The drawback of the subdivision approach is that it can cause changes in mesh topology that require significant recomputation of the discretised material models. A n alternative approach is to locally adapt existing mesh nodes so that they are coincident with the needle geometry. Instead of pulling a mesh node onto the needle by deforming the tissue model, the nominal reference mesh is adjusted so that the node is coincident with the needle. Thereafter, the node can be constrained to the needle shaft without consequence. This is illustrated in Figure 5.9. A low-rank mesh adaptation algorithm is detailed in Section 6.5. 5.4 General Non-linear Needle Trajectories The coordinate system shown in Figure 5.3 is fixed to the needle; therefore, as its orientation changes it is necessary to effect local coordinate changes in K _ 1 . Needle orientation changes, curved and flexible needles cause local changes in the constraint frames, as illustrated in Figure 5.f0. Consider a material coordinate change at a single tissue model node i. If the boundary conditions are uniform 5.4 General Non-linear Needle Trajectories 76 Reference mesh: Deformed mesh: Figure 5.9: (a) During simulation, as the needle tip moves from one element into another, the intersection point between the needle axis and the common edge is determined, (b) The node closest to this intersection point, node i, is projected onto the needle axis, along the direction of the common edge, by moving the nominal position of the candidate node in the nominal reference mesh such that it lies on the needle axis in the deformed mesh state. Figure 5.10: Trajectory changes (b) and deflection (c) cause local changes in needle orientation. (i.e., either a force or a displacement constraint along both the x- and ^-axes), then K l _ 1 undergoes a transformation of coordinates as follows: °u = K - 1 0 f =» u = A K~ A i 1 T 1 (5.5) l where °u and °f are displacement and force vectors in a nominal system coordinate frame, while u and f are the vectors after rotating the coordinate frame at the i l L th node by an angle 6. Matrix 77 5.4 General Non-linear Needle Trajectories A is composed of (2 x 2) rotation submatrices on its diagonal, and has the following form: A= 1 0 0 1 ... 0 0 ... 0 0 0 0 0 cos(0) - sin(0) 0 0 0 sin(0) cos(0) 0 0 0 0 0 ... 0 (5.6) 1 If node i has different boundary conditions along its two coordinate axes, then such a transformation is not possible, due to the mixed force and displacement variables in x and y . For example, if the i th node is sliding along the needle axis, then the model equations are expressed as shown in Equation (5.1), with x and y defined i n (5.3). A coordinate frame rotation at the i th node is expressed as follows: V cos(0) - sin(0) 0 V sin(0) cos(0) 0 fx Ji cos(0) - sin(0) 0 fV sin(0) cos(0) x_ where °uf °u¥ and o jx o jy Vf (5.7) are the displacement and force vectors of node i before the T coordinate frame change, and and lfx 1 fV are the displacement and force vectors of node i after the frame is rotated by an angle 0. The local coordinate transformation is derived by substituting (5.7) into (5.1): M x + J V y = ( K " ) ' [N x + M y ] 1 1 x 1 x = [M - ( K : _ 1 : )'JV] _ 1 [K~ M l - N] V , (5.8) 5.5 Needle 78 Flexibility where M and ./V are sparse transformation matrices that have the following form: 1 0 ... 0 f M N = = 0 0 0 0 0 0 cos(0) 0 0 0 0 cos(6») 0 0 ... 0 0 ... 0 ... 0 0 0 0 0 0 0 0 0 0 - sin(6>) 0 0 0 sin(0) 0 0 0 0 0 0 0 ... (5.9) Node coordinate frames are updated incrementally from one simulation step to the next, according to the change in needle orientation angle AO. If the needle is curved or flexible, then the local coordinate system transformations will vary along the length of the needle shaft. 5.5 Needle Flexibility Needles used for fine needle aspiration biopsy, certain types of nerve blocks, and brachytherapy are long, thin and flexible. In many cases, needle deflection has a significant impact on placement accuracy. This is demonstrated in Figure 5.11, where a steering motion is applied at the needle base, causing the needle to bend. A n elastic model is used to simulate the bending of stainless steel needles in a manner similar to that used to describe the tissue phantom in Section 3.5. In two dimensions, the needle is modelled as a thin elastic beam, as illustrated in Figure 5.f 2. The constitutive equations are discretised using the Finite Element Method (see Appendix C ) . Since the needle is long and slender, even modest 5.5 Needle Flexibility 79 Figure 5.11; A 10cm, 22 gauge biopsy needle deflects as it is manipulated within a soft tissue phantom. Figure 5.12: The needle is represented by a linear elastic beam deflections result in large global deformations and high strains. The small-strain assumption that is made in order to use the Cauchy strain tensor is not valid for such global deformations; therefore, the Green-Lagrange strain formulation is used. This is also referred to as quadratic strain [147,148] (see Appendix C for a detailed derivation). The practical relevance of these two formulations is illustrated in Figure 5.13(a), which shows simulated needle deformations due to a range of tip forces applied along the X-axis, using both linear (Cauchy) and quadratic (Green-Lagrange) strain. The linear strain model does not preserve element volume and hence needle length for large strains, and is not rotation independent; therefore, it performs poorly. Physical measurements of planar needle deflections were obtained experimentally, and compared with a two-dimensional needle model derived using linear elasticity and quadratic strain. A 20 gauge, f5cm Chiba biopsy needle was fixed at its base, and deflected by tip forces of varying magnitude, 5.5 Needle Flexibility 0.18 80 0.18 Linear Strain Quadratic Strain 0.16 - - Needle Model — Needle Measurement 0.16 0.14 0.14 0.12 0.12 0.1 0.1 0.08 0.08 0.06 0.06 0.04 0.04 0.02 0.02 E, > -0.1 0 -0.08 -0.06 -0.04 -0.02 X[m] -0.1 -0.08 -0.06 -0.04 -0.02 X[m] (a) 0 (b) Figure 5.13: Needle elasticity model with: (a) a comparison of linear and quadratic strain models, and (b) physical validation of the quadratic strain model against measured deflections of a 15cm, 20 gauge biopsy needle. from 0.05N to 0.5N, while planar deformation was measured. These measurements are shown in Figure 5.13(b), along with deflections computed by the numerical model. The needle model uses a Young's modulus of 980MPa, a Poisson ratio of 0.3, and a mesh of triangular elements with 16 x 4 nodes. Measured and computed needle displacements agree within approximately 2% over a large range of needle deflections. For tip deflections greater than 60mm, the steel needle begins to stiffen and the model becomes inaccurate. The set of non-linear model equations that describe the relationship between applied needle node forces and nodal displacements are expressed as: In = K ( u j u n => f n = n ((K )o + ^ S A A C 7 d ^ u T n =» f„ = ( K n K + t , , , n (5.10) 81 5.6 Solving the Coupled Needle and Tissue Models where the "stiffness matrix" K (u„) is a non-linear function of displacement. (K )o is the matrix n n of linear stiffness components, while the integral introduces higher order terms due to quadratic strain, as detailed in Appendix C. Needle node forces f are expressed as the sum of forces due to n linear strain and forces due to quadratic strain, f . nq ff force boundary conditions are applied to the needle model, then (5.10) can be inverted by means of the following iterative numerical method: ( u j ; i - (K„) + 1 0 f n - l ((u )i) nq n , (5.H) where (u )o = 0, and f„ is a vector of external force constraints applied at needle nodes. For the n model shown in Figure 5.13, (u )i converges to within 0.1% of u after 10 iterations of (5.11), when n n a tip force of — 0.3iV" is applied along the X-axis. 5.6 Solving the Coupled Needle and Tissue Models At each simulation time step, after all node boundary conditions and local coordinate frames have been updated, the linear elastic tissue model has the following form: x = (K- )V, , 1 (5.12) where (K )' is the updated system matrix. The specified constraint vector y' and the vector -1 of unknown variables x' contain node displacement and force coordinates for each observed node, expressed in their local material coordinate frames. The vector of unknown node forces and displacements, x', is computed by (5.12), and then reconstituted into displacement and force vectors expressed in a reference coordinate frame, based on the known boundary conditions and coordinate frames at each node. The needle is represented by the midline of the needle model, as shown in Figure 5.12. Forces at contact nodes on the tissue model are projected onto the needle model in order to compute the deflection of the needle nodes u„ according to Equation (5.11). The resulting configuration of the needle line segment in turn applies new displacement constraints on the tissue contact nodes. With appropriate damping, the coupled system reaches equilibrium within several simulation steps. In the examples given here, a velocity gain of 0.033 is applied at each needle model node (i.e., 5.6 Solving the Coupled Needle and Tissue Models (x )fc n the k th = 82 (x )fc_i + 0.033Ax„, where (x )/i; is the compound vector of needle node coordinates at n n simulation step), with a sample time of 2ms. Figure 5.14 illustrates the system convergence S 0.8 100 50 150 200 250 300 350 Simulation Step Number Figure 5.14: A step test illustrates the convergence of coupled needle and tissue models. following a step test; after a straight-line insertion, the needle base is stepped 3mm to the left. A normalised plot of lateral displacement at the needle tip indicates a damped response as the coupled system converges to within 1% of equilibrium after approximately 50 simulation steps (100ms). Figure 5.15 shows coupled needle and tissue models during a steering maneuver (to correct an off-target needle). Note that the response resembles that of a non-minimum phase system. Contact forces that occur along the needle shaft are integrated to determine the resultant driving force at the needle base. The force vector computed at each contact node is decomposed as shown in Figure 5.16, and referred to a coordinate frame located at the base. Base forces and torques resulting from each contact node are summed to produce the resultant driving force ft and torque Th- Examples of base driving forces and torques for a simulated needle insertion trajectory are graphed in Figure 5.17. For rigid needle models, the force summation can be expressed as: 0 -sin(<9j) 0 cos(0j) o ~o 0 'fi , (S.13) 5.6 Solving the Coupled Needle and Tissue Models 83 Figure 5.15: Simulated needle steering with needle flexibility. Lateral needle motion causes needle deflection, steering the needle toward the target (shown as a disk embedded in the tissue model). Figure 5.16: Contact forces that occur along the needle shaft are referred to the needle base and integrated to obtain the resultant driving force and torque. where o is the location of the i ih node and o is the location of the needle base, so that o - o the axial distance between the node and the needle base. The total wrench °[f T b in terms of a fixed environment coordinate frame, while the node force vector \ T^] T IS is expressed lies in the local coordinate frame at node i , which is rotated by an angle of Qi with respect to the fixed reference frame. 5.6 Solving the Coupled Needle and Tissue Models 84 Figure 5.17: A n example of the base forces and torques required to drive the needle along a simulated insertion trajectory. The trajectory proceeds as follows: (1) the needle tip makes contact with the tissue model, (2) the needle punctures the model surface and insertion continues along the y-axis, (3) insertion continues along the y-axis, with lateral base motion along the -x-axis, (4) insertion continues along the y-axis only, and (5) end of insertion. Small force spikes occur while the system converges after new working nodes are added. 5.7 Discussion 5.7 85 Discussion The simulation algorithm is physically-based, but under specific conditions only. A user may command arbitrary needle velocities, and may choose to insert and then retract the needle. A slip-stick force model was selected in order to permit such interaction, and the slip force threshold is set according to the shaft force distribution derived for insertions at a single velocity [50]. Therefore, simulated needle forces closely match measured values for insertions at this velocity. The forcevelocity relationships shown in Figures 4.18 and 4.19 can be used to adjust the force distribution according to insertion velocity. Velocity dependence and the slip-stick contact model introduce dynamic behaviour into the simulation, despite the quasi-static tissue model. A n asymmetrical bevel at the needle tip can cause additional deflections during insertion. The needle bevel introduces a lateral tip force, as shown in Figure 5.18(a); however, the relationship between the tip force, the physical shape of the needle tip, and the lateral force has not been established in experiments. The simulation algorithm assumes that the tip force is applied axially, (a) (b) Figure 5.18: The bevel deflection effect. as is expected to occur for needles with a symmetrical bevel (e.g., a trihedral tip, as illustrated in Figure 5.18(b)). A n algorithm flow chart for the simulation algorithm is given in Figure 5.19. 5.7 Discussion 86 Table 5.2: Flowchart Index. Box Number Reference 1 Section 5.3, Figure 5.6 2 Section 5.3.1, 5.3.2 3 Equation (5.4), Equation (5.5), Equation (5.8) 4 Section 5.2, Figure 5.4 5 Equation (5.4) 6 Equation (5.5), Equation (5.8) 7 Equation (5.12) 8 Section 5.5, Equation (5.11) 9 Section 5.6, Equation (5.13) 10 Section 7.1, Equation (7.2) 5.7 Discussion 87 C START Needle Insertion Simulation L ) 7 INPUT needle base coordinates commanded by user. Subdivide OR adapt mesh so that the needle is coincident with a candidate tissue node. Transform needle model geometry and motion into the tissue model frame of reference. Re-arrange boundary conditions and material coordinates at the new contact node. Compute constraint displacements. Compute forces and force thresholds from the shaft force distribution. Switch node to SLIP state: Rearrange boundary conditions. Switch node to STUCK, state: Rearrange boundary conditions. Compute displacement and force boundary conditions. Update material coordinate frames. Solve tissue model. Compute displacements and forces at contact nodes. Project forces at tissue contact nodes onto the needle model nodes. Add lateral bevel force if required. This is a function of axial force and tip geometry. Solve large strain needle deformation model. DAMP needle displacement. Integrate needle forces to compute resultant driving forces and moments at needle base. Compute Needle Manipulation JACOBIAN. L OUTPUT needle base force. DISPLAY tissue and needle model configurations. C 7 STOP Needle Insertion Simulation ) Figure 5.19: A flowchart of the needle insertion simulation algorithm. References for numbered boxes are indexed in Table 5.2. Chapter 6 Interactive Simulation of Needle Insertion Models Real-time computation of the insertion model described in Chapter 5 is complicated by the "curse of dimensionality" that is established by the large number of linear equations required to describe even small models. For example, a 10 x 10 node 2-D model, discretised as shown in Figure 3.14, results in 200 linear equations, and a system matrix K containing 40,000 elements. For a 30 x 30 node discretisation, the number of matrix elements swells to 3.24 million, while a 10 x 10 x 10 node 3-D body results in 3000 linear equations in as many unknowns (9 million matrix elements). It is clearly infeasible to solve such systems at haptic sample rates of 500Hz and higher. This chapter presents several numerical methods for reducing the computational complexity of the needle simulation algorithm, and describes an interactive haptic environment that has been implemented using these novel methods. 6.1 Condensing the Numerical System The behaviour of a continuum model that is discretised by the Finite Element Method is observed through the behaviour of a finite set of mesh nodes. For large volumes of tissue that are finely discretised, the size of matrix K in Equation (3.3) becomes large due to the large number of material nodes at which force or displacement need to be computed. For needle insertion simulation it is not necessary to consider the motion of nodes that are not visible (e.g., interior nodes), or the forces applied at nodes that are not in direct contact with the needle shaft. In Figure 5.8 it is evident 88 6.2 Low Rank Boundary Condition Updates 89 that the large majority of mesh nodes are neither visible nor palpable; therefore, the system of linear equations can be reduced to u w = K ^ f y y , where only the behaviour at a small subset W of mesh nodes (called working nodes) is explicitly considered. These are the tissue nodes that are in contact with the needle. A t run-time, the subset of working nodes W is selected and the system matrix reduced to K^y . As the needle penetrates the tissue surface, it intercepts hidden 1 nodes that need to be re-introduced into the reduced system by simply adding new matrix rows and columns to K ^ from K 1 _ 1 , which is precomputed. The matrix reduction approach is similar to the condensation techniques discussed in [134], and the Boundary Element Method selected by [62]; however, in this work access to the interior of the tissue volumes is retained for quick inclusion when needle penetration occurs. Most tissue model nodes that are not in W (i.e., that are eliminated during condensation) are traction nodes that require a zero force boundary constraint. Therefore, it is convenient to work with K _ 1 , rather than K . Prior to simulation, the entire stiffness matrix K is inverted and stored as ( K ) o . The condensed working matrix K^y is dynamically constituted from the rows and _ 1 1 columns of ( K ) o , as required during simulation (i.e., when a node is added). _ 1 6.2 Low R a n k Boundary Condition Updates The slip-stick contact model described in Section 5.2 results in frequent boundary condition changes, each expressed as a low-rank update to K^y : 1 ( K ^ ) ' = K ^ - SLEi , Pi 1 where pi is the i th pivot of K ^ ; Cj and the exception of their i 1 th (6.1) are the i 1 th column and i th row of K^y , respectively, with 1 coordinates, which are both set to (pi — 1). ( K ^ ) ' is the new system 1 matrix. This approach to boundary condition changes results in an 0 ( n ) computation rather than 2 the 0(n ) 3 operation required to re-invert stiffness matrix Kyv, where n is the number of mesh nodes in VV. This is similar to the capacitance matrix strategy presented by James and Pai in [62]. 6.3 Low Rank Material Coordinate Updates 6.3 90 Low R a n k M a t e r i a l Coordinate Updates Local material coordinate frames at all contact nodes in W are updated during each simulation step, based on the needle configuration, as described in Section 5.4. From Equation (5.5), local coordinate changes at free traction nodes, or fully displacement constrained nodes are computed as follows: °Mw = K vv °fw X Mw = A K^A T 1 f = (K^y )' f 1 v v . l w (6.2) The matrix A is (2 x 2) block diagonal (see Equation 5.6); therefore, (Kyy )' is inexpensively com1 puted. In Equation (5.8) the local coordinate transformation for a single node that is in the slip state requires that ( M — K ^ T V ) be inverted. This matrix inverse is sparsely populated, and is expressed analytically as follows: (M-K^iV)" 1 0 1 = ... 10 1 a (l,2i) a (l,2i+l) a (2,2i) a (2,2i+l) (6.3) 0 0 0 0 0 0 a a (2i,2i) (2i+l,2i) (2s,2i) a a a (2i,2i+l) (2i+l,2i+l) a (2s,2i+l) ••• 1 6.4 Solving the Tissue Model 91 where, a sin6>[fc (i,2i) = (1)2i+1) (cos6> + k i+i,2i) sin0)] ^ + {2 sin0[fc 2j)fc(2i+i,2i+i) (1) 6] sm E °(l,2i+l) _sin6'[fc j (cos6» ^ (1)2 fy^i+i) ) sin *)] 6 h sin^[fc(i,2i+i)fc(2i,2i) shig] E cos0 +fc(2i+i,2t)sinf9 (2i,2i) a a fc(2i,2i) (2i,2j+l) a a E E ' fc(2i+l,2i+l) (2i+l,2i) s i n 6 > E (2i+l,2i+l) — ^ s m 6 > - cosfl +fc(2i,2i+i)sing > = - cos 0 + cosr?sin6'(A; 2i,2j+i) - fc(2i+i,2i)) 2 ( + Sin 0(fe(2i,2i+l)fc(2i+l,2i) - k(2i,2i)k(2i+l,2i+l)) , 2 and fc( ) is the element at the m th m;n row and n t f t column of K ^ , . If the needle is rigid and straight, 1 then the rotation angle 0 is constant for all tissue nodes in contact with the needle shaft. If the needle is curved or flexible, then 0 varies according to the shaft geometry. 6.4 Solving the Tissue M o d e l Once all boundary condition changes and local coordinate frame changes have been computed at each contact tissue node, using the fast low-rank update methods described in this chapter, the vector of unknown forces and displacements at working nodes are calculated as follows: xw = (K^ )'y' , , (6.5) 1 v vv where (Kyy )' is the updated system matrix after all boundary condition changes and local coordinate 1 frame changes have been made. The specified constraint vector y ' w and resulting solution vector x'yy contain displacements and forces for each working node. The set V of tissue model displacements at nodes that are not in W , but that are required for 6.5 Node Intercepts and Mesh Adaptation 92 visualisation, are computed as follows: u where u v = [ ( K ) ] v w fw , (6-6) _ 1 v 0 is a vector of displacement variables that are required and f w is the vector of forces at working nodes (known from (6.5) and the applied boundary constraints). [ ( K ) o ] ^ ^ is the _1 v w matrix that contains rows and columns corresponding to the variables in V , and the nodes in W respectively. 6.5 Node Intercepts and M e s h A d a p t a t i o n New needle-tissue contact nodes are identified as described in Section 5.3, with the displacement of each neighbouring node i G jV computed as follows: fx (6.7) 0 where [/? /J] is the force vector at the j th fc(2i+l,2j+l) f V node (j € W). k(ij) is the element at row i and column j of ( K ) o - A fast implementation of the mesh adaptation scheme outlined in Section 5.3.2 _ 1 avoids having to recompute ( K ) o from scratch when reference node coordinates are relocated in _ 1 order to be coincident with the needle shaft prior to being added to the list of working nodes W . The relocation of one mesh node in the discrete tissue model affects only the elements that share it as a vertex, as shown in Figure 6.1. This set of elements £ corresponds to a submatrix within the stiffness matrix K that must be updated due to the new discretisation. Each discrete element on the material domain $7 results in a relationship of the following form: f = Kf e u , (6.8) e 6 x 6 ) as described in Appendix C. The element stiffness matrix K e describes the force-displacement relationship for a single element in isolation. K is the composite of element stiffness matrices over the entire body. The relocation of one node, as shown in Figure 6.1, requires changes in the stiffness matrices for all elements in £. For example, consider the element e with vertices i, m and n in Figure 6.1. K e and (K )' e are the (6 x 6) stiffness matrices for element e, computed before and 6.5 Node Intercepts and Mesh Adaptation (K'') 93 0 —< 0 « .r A* — \ (K y e \ J \ Figure 6.1: The relocation of node i affects the set of connected elements £ (shaded). after node i is relocated, respectively, and AK = (K )' — K . A single element update results in e e e a rank-6 update of K : K' = K + VAK V . e (6.9) T V is a (2N x 6) selection matrix that samples K from K , and is defined as follows: e [V]kl (6.10) = where N is the total number of nodes in the model discretisation, C is an ordered list of row indices in K that correspond to the rows of K , and Ci is the I e matrix K - 1 th entry in this list. The updated inverse is required for simulation, and can be computed using the matrix inversion lemma (Sherman-Morrison-Woodbury formula) [149] as follows: (K + V AK V )~ e T 1 = K " - K~ VAK (I 1 1 + V K' e Equation (6.11) requires the inverse of (/ + V K~ T T K(2AT 2A0X ^ o r e l e m e n t e 1 1 V Aif )( e 6x6 V AK )~ V e l T K~ ] (6.11) ^ , but avoids having to re-invert bi the example, £ = (2i - 1, 2i, 2m - 1, 2 m , 2n - 1, 2 n ) . The remaining five elements in £ are updated in the same manner. Note that it is more efficient to update the six element stiffness matrices simultaneously, by means of a single rank-14 update, rather than six individual rank-6 updates. The resulting updated matrix K _ 1 replaces ( K ) o . -1 This procedure can also be implemented as a sequence of boundary condition changes (5.4) and 6.6 Haptic Needle Insertion matrix updates AK e 6.6 94 made to (K Haptic Needle Insertion A real-time implementation of the needle insertion simulation algorithm allows users to experience both visual and kinesthetic feedback while executing a virtual planar needle insertion. The system comprises a host processor with a graphical display, a real-time processor, and a haptic interface with its associated I/O interfaces and power amplifiers, as illustrated in Figure 6.2. Both the host and real-time processors are ordinary Pentium PC's. The haptic interface is similar to the planar Host Processor Current Driver DAC Decoder Ethernet Graphical Display Haptic Interface Figure 6.2: The haptic simulation system. pantograph mechanism that was used to manipulate the needles during insertion experiments, and is controlled by a real-time processor, as described in Appendix A. A Quanser MultiQ-PCI interface facilitates I/O between the haptic interface and the real-time processor. The visual environment is rendered by the host processor, which receives the needle and tissue model states from the real-time processor via an Ethernet connection. The low-rank numerical operations described in this chapter require 0(s ) computations per 2 sample period in order to compute the needle insertion models, where s is the number of contact or working nodes in W. In contrast, a simulation involving all nodes in the discretised tissue volume, with full matrix inversions, would require 0(n ) computations, where n is the total number of mesh 3 nodes (n S> s). Displacements at all tissue model nodes are computed in order to graphically render 6.6 Haptic Needle Insertion 95 the complete deformation state of the model; therefore, a further 0(ns) computations are required. Figure 6.3 presents a series of display snapshots taken during interactive virtual needle insertion in a planar environment, with a rigid needle and a 361-node tissue model computed at 512Hz. The model and algorithm specifications are summarised in Table 6.1. Both the model simula- Table 6.1: Model and algorithm specifications for haptic simulation. Material Type linear elastic, isotropic, 100 x 100mm Material Young's Modulus inhomogeneous Material Poisson Ratio inhomogeneous Number of Nodes Element Type 361 (19 x 19) 3-node triangle, plane stress with thickness of 11mm Strain Cauchy, linear strain analysis Needle Model rigid Mesh Subdivision none Mesh Adaptation none (nodes "snap" to needle) Model Sample Rate 512Hz tion and haptic interface control algorithms are computed by the real-time processor under the VxWorks® operating system, without any particular effort to optimise code. Contact forces and torques are fed back to the operator via the haptic interface, as facilitated by a four-channel control architecture that is described in Appendix A . Figure 6.3: Interactive virtual needle insertion in a planar environment. A n example of a simulated lung nodule biopsy is shown in Figure 6.4. The tissue geometry 6.7 Discussion 96 was taken from a segmented C T image, while tissue parameters were chosen in an ad-hoc manner, guided by values published in the literature. Figure 6.4: Interactive virtual lung nodule biopsy. 6.7 Discussion The discretised linear elastostatic model has worked well for haptic needle insertion. Local phantom stresses due to the observed needle force distribution were fairly low, so that axial strains are low. There is little lateral motion and rotation during typical needle trajectories; however, users are not constrained during haptic simulation and may generate large strains. More complex material models may be substituted in order to account for these large strains, and perhaps non-linear constitutive behaviour, without affecting the force distribution modelling methodology presented here. However, not all models will be amenable to real-time simulation. Feedback from the medical community has been very positive, despite rather simple model geometries. For example, upon trying the interactive haptic simulation, a respected radiation oncologist in Vancouver remarked that he "had not felt a physical phantom to be as realistic". A controlled user survey has not yet been undertaken to quantify such feedback. 6.7 Discussion 97 This chapter has presented several fast numerical techniques that allow the needle insertion simulation algorithms to run at interactive rates. The algorithm was implemented as a haptic simulation that allows users to perform virtual needle simulations with both visual and kinesthetic feedback. The following chapter considers the problem of needle steering and trajectory planning using the insertion models and algorithms developed thus far. Chapter 7 Needle Steerability and Model-Based Trajectory Selection This chapter outlines the problem of motion planning and steering for the placement of flexible needles in soft tissues, and presents a manipulation and planning framework that is based on the needle insertion models presented in earlier chapters. 7.1 M o t i o n Planning and Needle Steering For the sake of exposition, we define the concepts in 2-D space, for which tissue and needle models have already been derived. Consider elastostatic model approximations for both tissue and needle. A suitable discretisation, based upon Finite Element Analysis yields a finite set of tissue nodes n t and a finite set of needle nodes n n [52]. A t time step k, the system is said to be in configuration q^, for which the positions of all nodes nj and n n are known; the external forces applied at all nodes nt and n are known; and the boundary conditions and material coordinate frames at each node in n nt are known. The following equations are satisfied at each step k: it = G ~ (ut) u = Au in^Gn'^n) L = Bi_ 1 t where u and u t ra t n (7.1) t are vectors of tissue and needle node displacements, and f_ and f t n are vectors of tissue and needle node forces, respectively. We use Cauchy strain for the tissue model, resulting in 98 7.1 Motion Planning and Needle Steering 99 a linear function G~j~ (u ); and Green-Lagrange strain for the needle model, resulting in a non-linear t function ^ ( u ) . A and B represent the coupling conditions between needle and tissue models, 1 n based on the contact states at each contact node (see Chapter 5). The needle is manipulated from its base. Given a needle configuration and base velocity [if, y\ 6b] at sample time k, the new T configuration q^+i after time interval St can be computed. The needle insertion simulation evolves as described in Chapter 5. The motion planning problem is to find a sequence of control parameters (needle base velocities) that guides the needle tip node ( n ) i to a target tissue node T = (nj)j, as n shown in Figure 7.1. Such a condition can be detected by observing the system configuration q^ Figure 7.1: A flexible needle inserted into a body of soft tissue, toward the target T. The needle is manipulated from its base, while numerical tissue and needle models characterise motion and force at node points. at time step k. There may be additional constraints due to tissue force limits and displacement constraints that avoid the puncture of obstacles within the tissue. If the needle configuration is specified using only the position and orientation of its tip, then the problem of determining a needle trajectory becomes that of guiding its tip along a feasible path. We say that a needle is "steerable" if base motion can be used to control the position and orientation of the needle tip with respect to the target. This is illustrated in Figure 7.2. A Needle Manipulation Jacobian J embodies the concept of needle steerability, and is defined by the relationship between the needle tip and base velocities, [it yt 6t] T and [if, i)b 9b] T respectively. In the 2-D, three-degree-of-freedom case, J is a 3 x 3 matrix and is a function of the 7.1 Motion Planning and Needle Steering 100 Figure 7.2: (a) Motion ofthe needle tip [x yt 0t\ with respect to the motion applied to the base [x y 6b] determines needle "steerability". Needle steerability depends on the tissue and needle configuration q^. T t T b b needle and tissue configuration q^: dxt dx dxt. dyb dyt dx dyt dyb oe dB _dx dB 9y dBt d9 . b J(Qk) = b t b t b dxt ae b (7.2) dyt b b The velocity relationship can be written in terms of infinitesimal changes in position, as follows: Ax Ax t b ~ J{Qk) _A0 _ t Ay 6 as || (Ax , Ay , b A$ ) b (7.3) b _Ad _ b A closed-form expression for J is not available due to the complexity of needle-tissue interaction. Indeed, the Jacobian depends upon the needle and the tissue configuration at a given instant. However, needle tip displacements caused by perturbations [Ax b Ay A6 ] T b b applied to the needle base can be measured in order to numerically determine the needle manipulation Jacobian matrix. For example, if the needle base is perturbed by [Ax [Axt Ay A9 ] T t t is observed after several simulation steps, then by (7.3) the first column of J is estimated to be approximately -^[Axt in similar fashion. 0 0 ] , and a tip displacement of T b Ayt A9t] '. T The remaining columns of J are computed 7.2 Needle Motion Planning using Potential Fields 101 As long as the needle manipulation Jacobian matrix J is non-singular and well conditioned, the tip velocity can be controlled by the base velocity. This means that J~ exists and Equation l (7.3) can be inverted in order to compute the needle base displacement that will result in a desired displacement at the needle tip, as long as the displacements are "small enough". 7.2 Needle M o t i o n Planning using Potential Fields Consider the simple tissue domain illustrated in Figure 7.3, and the problem of determining a path that will take the needle tip from an initial position to the target, while avoiding the obstacle. The needle tip configuration space C is parameterised by the position and orientation of the needle tip, while the obstacle is mapped to CB. The potentialfieldapproach proceeds by placing an attractive potential at the target location in C, and repulsive potentials at obstacles that correspond to CB. The gradient of the resultant potentialfieldguides the motion planner. \j ^ Possible Goal Configuration \ • Obstacle • Target Tissue Domain Figure 7.3: A simple 2D tissue domain that contains a target location, as well as an obstacle that is to be avoided by the needle. Clearly, multiple starting and goal positions exist. 7.2.1 The Attraction Potential For an approach based on potential functions, an attractive parabolic potential well U is placed at the target [139]: U(q ) k = ?k\\q -q oai\\ 2 k 9 , (7.4) 7.2 Needle Motion Planning using Potential Fields where k is a positive scaling factor and q k £ C. 102 The function has a minimum at U positive elsewhere, monotonically increasing in all directions away from q i- goa G C and is q i goa The gradient of this potential function is a generalised attraction force F: F(q ) k Note that q, q l k goa GR , 3 a Q heading component, and q = -VU{q ) k = -k{q - q i) k . goa (7.5) d include an orientation component. The configuration i goa q k includes a tip includes the heading angle between the needle tip and the target. An example of an attraction potential and the local directions of the field gradient, for the environment shown in in Figure 7.3, are plotted in Figure 7.4 (orientation potentials are not shown). -0.015 \ . \ . \ \ \ - 0 02 \ \ \ \ \ \ \ -0025 Si \ \ \ \ Si / / \ \ / / \ \ I I i / / • . . \ Si S i Si \ \ i / ^ ^ ^ ^ - - 0 03 • • Y[m] X[m] i • • • 0 00!. / 7 M 0 X[m] (a) \ K. K * . 0.005 (b) Figure 7.4: (a) Attraction potential over the tissue domain, and (b) the local direction of the field gradient. 7.2.2 The Repulsion Potential A repulsive potential U , placed around each obstacle region, prevents the needle tip from intercepting an obstacle [139]. This barrier must be carefully selected so that it does not affect needle motion when it is sufficiently far from the obstacle. A simple repulsive potential function may be written as: fcG-A)' 0 (7.6) if p > p . 0 7.2 Needle Motion Planning using Potential Fields 103 where n is a positive constant, p is the Euclidean distance || q — q b ||, q bs is the position of the k 0 s 0 obstacle, and po defines a maximum region of repulsion (orientation states are neglected here, and are considered later). If the obstacle is not a circular/spherical point object, then |<ft; QobsW becomes — min 11(fa — q bs\\- The scaling factors k in (7.4) and n in (7.6) are chosen such that the attraction 0 and repulsion potentials are of the same order of magnitude. For a point obstacle, the resulting repulsion force is isotropic within a circular region of influence. Such a field is not feasible for needle motion planning because of the large lateral deviation that may be required at the needle tip (illustrated in Figure 7.5(a)). Needle motion becomes increasingly constrained as the needle penetrates further into the tissue (i.e., the condition number of J increases with needle depth). Therefore, needle trajectory compensation must be made early, and is therefore facilitated by selecting an elliptical or hemi-elliptical repulsive potential field shape, as shown in Figures 7.5(b) and (c). (a) A n elliptical region of repulsion is defined in Cartesian coordinates by (b) (c) (d) Figure 7.5: Examples of repulsive potential field shapes: (a) circular, (b) elliptical, and (c) hemi-elliptical. In (d) two distinct points near the needle tip experience different repulsion forces, and a moment is applied to the line segment connecting these two points. This is the repulsive torque at the tip. 7.2 Needle Motion Planning using Potential Fields 104 h(x,y), and by h(r,9) i n polar coordinates: M*,v) = (!) + (if) = i 2 2 ab 2 2 2 (7.7) P0 , r = V 6 cos (0-0)+a sin(0-</>) 2 where a and b are the radii at the major and minor axes of the ellipse. The normal n at a point (a;, y) on h is determined by finding the gradient of h at (x, y): n(x,y) "2s 2y = Vh{x,y) = T 2 r b r cos(0 - 4>) 2 ab I 2 a?r sin(<9 - <f>) (7.8) 2 where 4> is the orientation angle of the major axis, chosen to be /.{X b 0 — X s g o a The anisotropic, i). elliptical repulsive force field is described as: ^\P PO ) = I VU{q ) k 7? II cos(e-<p) [b P0 2 •(g-0)J if (p ^ p ) and (d 0 ap 2 ' ( « - 0* ) ]s i IIn if (p > po) or where p is the Euclidean distance || X is 0 ^d ) flt (7.9) > dgt) ||, and the radius po that defines the region of repulsion — Xb t ffo s ab V& cos (fl-</>)+a sin {9-4,)' 2 2 2 2 d g o " 1 S X s o a / ~ X o b s " ^ 1 S " X g o a l ~ X t i p "' X o b s " X g o a l ^* are the Cartesian coordinates of the obstacle, target and needle tip, respectively. In practice, it is necessary to match the field gradient along the major and minor ellipse axes by scaling both p and po by f Q in (7.9): r 2 2 ii>, cos(0-0) o ~\T if (/> ^ Po) and (d ll fl0 <d ) flt (7.10) if (p > po) or (d > (i ) po 9t This defines the translational case, where CB lies in R . The effect of a repulsive potential field 2 on needle tip heading is included by considering the repulsion force at two points P i and Pi on the needle shaft, both close to the tip. This is shown in Figure 7.5(d). The force "exerted" by the field at each end of the line segment that connects P i and P results in an effective moment that is used 2 7.2 Needle Motion Planning using Potential Fields 105 to guide angular motion at the needle tip. where X\ and X are the positions of P i and P , while U\ and f/ are the potentials computed at 2 2 2 P i and P respectively. 2 Examples of the repulsive potential field and the local directions of the field gradient, for the environment shown in Figure 7.3, are plotted in Figure 7.6. Angular potentials and forces are not shown. The total generalised repulsive force, including both translation and orientation components, Figure 7.6: (a) Repulsion potential over the tissue domain, and (b) the local direction of the field gradient. is written as follows: r — 7.2.3 (7.12) Minimising Total Potential The sum of the attraction and repulsion potentials results in a field that attracts the needle tip toward a goal configuration (tip coinciding with target), while repelling it from obstacles. This is shown in Figure 7.7. Similarly, multiple obstacles can be accommodated by adding each of their repulsive fields together (see Figure 7.8). 7.2 Needle Motion Planning using Potential Fields 106 Figure 7.8: Potential field gradient directions when three obstacles are present. In classical problems that consider path planning for rigid or articulated bodies with mass, the potential field gradient force is often applied and acceleration derived to obtain motion. The notion of mass at a needle tip is not meaningful; therefore, the force is used to indicate the direction of 107 7.2 Needle Motion Planning using Potential Fields motion only. The motion of the needle base AQ& can then be expressed as: AQ = -J-' (W f f o l +fo*)) , b where F(q ) k and F(q ) k V \\F(q )+F( )\\) k (7.13) qk are generalised attraction and repulsion forces respectively, and W is a diagonal weighting matrix that accounts for differing units. In practice, the tip motion is expressed in a coordinate frame centred at the target location, and the Jacobian is computed accordingly. 7.2 Needle Motion Planning using Potential 7.2.4 Fields 108 The Planning Algorithm The model simulation and trajectory plan updates are computed according to the following algorithm: 1. Select the i n i t i a l needle configuration. 2. Locate the target and obstacle in the tissue model. 3. Compute the gradient of the t o t a l potential f i e l d at the needle t i p , to (7.5) P i , according and (7.10). 4. Compute the t i p heading repulsion torque according to Equation (7.11). 5. Set the desired t i p displacement in the direction of the f i e l d gradient. 6. If J - 1 exists ( i . e . , cond(J) ^ e) , compute the required base displacement by (7.13), else generate error condition. 7. Advance the needle base by the displacement (7.13) and solve (7.1) to find the new needle/tissue configuration. 8. Go to 2 and iterate u n t i l the needle t i p and the target are coincident. 7.2.5 Simulation Results The motion planning approach is shown in simulation, where an updated potential field is computed at each time step, and used to determine base displacements AQb- Figures 7.9 and 7.10 show the simulated insertion trajectory for an environment that contains three obstacles, with a sample potential field illustrated in Figure 7.8. In practice, the planning algorithm is sensitive to the selection of the free parameters a, 6, po and W. The shapes of the elliptical repulsive fields are governed by a and b, and are chosen so that trajectory adjustments are made early (before reaching the obstacle), and with care not to introduce local minima. The maximum region of repulsion po is selected based on the size of each obstacle. The diagonal weighting matrix W scales the desired tip motion in the direction of the potential field gradient. It accounts for measurement units and the relative importance of motion in each degree of freedom at the needle tip. 7.3 Experimental Results 109 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 Figure 7.9: Simulated needle insertion along a path selected by the potential function approach. Figure 7.10: The needle and tissue state during a simulated needle insertion along the path automatically prescribed by the potential field gradient. 7.3 Experimental Results Figure 7.11 illustrates a planned trajectory for a 10cm, 22G Franseen biopsy needle in a soft tissue phantom (E — 20KPa) with one obstacle and one target as shown in frame 1. The simulated needle trajectory was applied in experiments using the physical measurement system presented in Chapter 3. A trihedral Franseen tip was used in order to minimise deflection due to the bevel 7.3 Experimental Results 110 © CD j f j | ; : - © j • 1—1 :::::::::•::::::::: • -0.06 -0.04 -0.02 0 • 0.02 0.04 0.06 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 Figure 7.11: Simulated needle trajectory plan. (see Figure 5.18). Camera images corresponding to frames 2 and 3 in Figure 7.11 are shown in Figure 7.12. The measured needle trajectory closely resembles the simulated trajectory; however, Figure 7.12: Robot controlled trajectory in a transparent P V C phantom. The arrows indicate the initial insertion point. Note slight divergence due to out-of-plane deformations. this is not immediately apparent from the images due to out-of-plane deformations. The grid of black markers (visible in the camera images) is applied to the top surface of a 12mm thick P V C phantom, whereas the needle moves in a plane that is 6mm below this surface (see Section 3.2). Out-of-plane deformations cause the needle and marker planes to diverge slightly (clearly evident 7.4 Discussion 111 at the insertion point indicated). Based on the magnitude of the divergence at the insertion point, we estimate approximate target and obstacle locations in the needle plane. These are superimposed on the image frames. 7.4 Discussion A feasible path from the current needle configuration to a goal configuration in each static potential field is not guaranteed to exist, e.g., if the free-space in C is disjoint. Both the tissue and the needle deform depending upon the path, thus the evolution of the potential function from one time step to the next is very important. Even if a feasible path is shown to exist at each time step, convergence is not guaranteed, e.g., a target that is embedded in soft tissue might be pushed away by the needle as it approaches, and continue to move away as the needle tip is steered toward it. In this case, the manipulation Jacobian J is close to singular, a condition that is easily detected (cond(J) > e, where e is the threshold at which J becomes numerically ill-conditioned; its value depends upon units and numerical precision). Such a state may be used to indicate that the initial configuration was inappropriate, or that the potential field parameters need adjustment. In some cases, a trajectory may fail to reach the target despite the manipulation Jacobian remaining nonsingular, due to limited achievable needle curvature. The planning technique may also be useful for evaluating placement uncertainty, target "reachability" , as well as to prescribe insertion points and steering strategies. Local minima can occur in the potential function, and may be avoided by selecting an alternative set of potential functions. Randomised planning, or hybrid strategies that include some form of configuration space search [139] may also be required in some cases. In addition, elliptical repulsive fields may not be suitable for all obstacle geometries, and may need to be modified for some tissue structures. In general, the existence of solutions to such motion planning problems is difficult to prove. Indeed, such proofs do not exist for problems with far fewer degrees of freedom, such as the tractor-trailer problem outlined in Section 1.2.3. Additional force constraints may be important in clinical applications, in order to avoid needle trajectories that exert excessive force on local tissue structures. For details on incorporating such constraints into motion planning problems, see [139]. Chapter 8 Conclusion 8.1 Summary and Conclusions Needle insertion models and simulation are of interest for the development of physically-based virtual training, planning and guidance systems that are aimed at reducing the incidence of complications and improving needle placement in clinical practice. This thesis presents a novel approach for estimating needle shaft forces and tissue phantom behaviour using a measurement system and physical deformation models. Needle shaft force distributions estimated using this approach were used to simulate needle and tissue phantom motion, as well as to predict target motion and needle trajectories due to phantom and needle deformations. The integral of the estimated shaft force distribution was found to be in close agreement with driving forces measured at the base of the needle during experiments. In addition, simulated phantom deformations matched measured phantom deformations, indicating that the method works well. A fast, optimised simulation algorithm for interactive needle insertion has been developed by judicious application of matrix condensation and low-rank update techniques that reduce computation and permit the simulation of significant models at high sample rates on an ordinary P C . Unlike existing single axis simulations, steering torques and lateral needle forces can be felt, while tissue deformation is observed—illustrating the potential of physically-based needle insertion simulations for planning and training purposes. The concept of needle steering has been introduced, and a model-based needle trajectory planning technique proposed. It is shown, both in simulation and experiments, that complex steering of the needle tip in soft elastic phantoms can be performed by manipulating the needle base. 112 8.1 Summary and Conclusions 113 The specific contributions of this thesis are summarised as follows: Modelling: A method for estimating needle forces during insertions into soft tissue phantoms has been developed. In prior work, driving forces at the needle base were measured during linear insertions, and soft tissue motion was largely ignored. In this work, the relationship between needle and tissue phantom motion is characterised and used to derive an estimate of the force distribution along the needle shaft, rather than the point driving force at the base of the needle. The approach has been demonstrated using homogeneous and layered inhomogeneous soft elastic tissue phantoms, and validated by means of experimental measurements and simulations. The methodology may be extended to in vivo experiments using emerging imaging and tissue characterisation techniques that are based on Computed Tomography (CT), Magnetic Resonance Imaging (MRI) and Ultrasound Imaging. Simulation: Needle shaft forces, estimated from experimental measurements, are used to constitute a physically-based simulation of needle and soft tissue motion during needle insertion. Soft tissue interaction, multi-degree-of-freedom needle trajectories and flexible needle models have not been demonstrated in previous work on this topic. Interactive Simulation: Fast numerical techniques have been developed in order to solve coupled linear needle and tissue models for insertion simulations at interactive rates. This has led to a number of analytical advances: • It is shown that boundary conditions need only be explicitly considered at the small subset of tissue model nodes that are in contact with the needle. The number of system variables is reduced by means of a matrix condensation approach; however, unlike previous applications of this technique, the set of variables that are explicitly considered in the condensed system evolves as the needle penetrates the tissue volume. Variables are added to, and removed from the condensed representation as the state of penetration changes during simulation. • Boundary condition changes are applied using low-rank system updates, according to a slip-stick contact model. The boundary conditions are more complex than those seen previously for point and surface interactions, due to the volumetric interaction and the continuously evolving condensed system representation. In addition, varying constraint S.J Summary and Conclusions 114 directions—due to needle orientation, curved and flexible needles—require local material coordinate changes. A n analytical inverse has been derived to efficiently compute material coordinate updates at contact nodes. • Fast local mesh adaptations using low-rank matrix updates were developed i n order to match needle and tissue mesh geometries, while avoiding computationally expensive matrix re-inversions. • A large-strain needle model is coupled to the tissue model by matching boundary conditions, and the coupled system is solved iteratively. Needle deflections in soft tissues have been shown to have a significant impact on placement in practice, and are simulated here for the first time. • A compelling haptic environment, implemented based on the interactive needle insertion simulation algorithm, has been developed and demonstrated. This is the first such system to include physically-based needle forces and torques, soft tissue deformation and realtime interactivity. Steering and Manipulation: The concepts of needle steering and manipulation have been proposed, and a Needle Manipulation Jacobian has been introduced for this purpose. In clinical applications, needle steering may be able to avoid obstacles, and reach regions that are currently inaccessible using straight-line trajectories. Needle motion manipulation and steering are illustrated using the simulation system that has been developed in this work. This is the first time that such concepts have been demonstrated using physically-based needle motion models. Motion Planning: Geometric constraints and sensitivity analysis for needle insertions have been considered in prior work; however, this thesis presents the first needle trajectory motion planning technique that incorporates soft tissue motion, needle flexibility and a physically-based contact model. The high-degree-of-freedom motion planning problem is solved using a task-space potential field approach. The algorithm has been demonstrated in simulation, and validated for simple trajectories in tissue phantoms. 8.2 Future Work 8.2 115 Future Work The limitations of this work can be addressed as follows: Linearity: As a first approximation, linear elastostatic tissue models were used to model tissue deformation, but these may not be representative of tissues that exhibit more complex behaviour. Many of the methods presented in this work are independent of a specific tissue model Q, and may be generalisable to other—more realistic—tissue models. The interactive simulation techniques rely on the assumption of linearity in order to achieve satisfactory real-time performance; therefore, significant contributions will be required for non-linear systems. Scalability: The two-dimensional material models that are used throughout this thesis are quite easily extended to three-dimensions by using mesh geometries containing solid elements (e.g., tetrahedra and hexahedra). 3-D models were used for model identification and force computations in Sections 3.6 and 4.2; however, model simulations were developed and illustrated in 2-D. For 3-D models, the simulation algorithms will be similar; however, for practical problems, the size of the model data-sets may exceed available memory resources, requiring hierarchical representations, data caching and predictive caching. Interactive simulations of complex tissue models in three dimensions will benefit from the multi-resolution and mesh adaptation techniques described in [150-155]. Tissue Phantoms: Soft plastic tissue phantoms were used for modelling and validation because of the high level of control and repeatability that they afforded during experiments. Further experiments are required to characterise insertions into second generation tissue phantoms that are representative of typical soft tissue structures that might contain several distinct layers, inclusions, rigid structures, fluid pockets, etc. This will set the stage for ex vivo tissue samples, and finally in vivo measurements. Accurate 3-D deformation measurement of tissue, or even tissue phantoms is challenging [144]; however, several medical imaging modalities may be used, e.g., three dimensional ultrasound, Magnetic Resonance Imaging (MRI), TaggedM R I and Computed Tomography (CT). Measurement and elastography methods based on such imaging modalities may be suitable for extending the modelling methodology presented in this work [156,157]. Tractable physically-based, non-linear, anisotropic, inhomogeneous, multi-phase tissue models are also not yet available. 8.2 Future Work 116 Contact Model: The stick-slip contact model used in simulation has not been validated for needle insertions that do not occur at a constant velocity. The contact model, working in conjunction with the force-velocity relationship measured in experiments, should be evaluated. The effects of needfe lubrication can also be investigated using the experimental system presented in this work. Dynamic tissue models may also be useful. Clinical Application: This work provides new approaches for solving problems that have clinical relevance; however, several advancements are required in order to apply these fundamental concepts to clinical applications. Imaging systems that can accurately measure the complex deformation of tissue volumes are required. Methods for determining tissue parameters in vivo, such as ultrasound and M R I elastography [f56], are still in development. In addition, the needle force estimation technique will need to be applied differently in a clinical setting. For example, an iterative approach might proceed as follows: (1) advance the needle, (2) observe tissue deformation, (3) estimate shaft forces, (4) predict motion based on estimated shaft force, (5) select appropriate motion, (6) advance needle, repeat. In cases where tissue models and force estimates are not feasible, it may still be possible to estimate the Needle Manipulation Jacobian (e.g., by observing needle tip motion following small perturbations applied at the base) and to guide the needle based on the potential field approach. fn addition, the concepts presented in this work can be extended in a number of interesting ways: Human Factors: Needle insertion simulations for medical training will be based on anatomical models and clinical procedures. The effectiveness of these models, and the level of realism required, are yet to be determined. Prototype systems should be developed with help from the medical community in order to investigate these human factors. Haptic Simulation: The prototype haptic system can be extended to include mesh adaptation and flexible needle models. Some further work is required to implement the mesh adaptation algorithm at haptic rates, and this must be done prior to incorporating a flexible needle model, due to issues of stability (the coupled needle-tissue system can become unstable when nodes are "snapped" to the needle shaft, particularly in stiff tissues). Trajectory Analysis: Applications involving trajectory analysis (trajectory feasibiiity, planning and optimisation), navigation and automatic control are extremely exciting, and may lead to sig- 8.3 The Big Picture 117 nificant improvements in tool placement for medical applications in the future. In addition, these methods may spawn new treatment and interventional techniques. For example, insertion models may be used to test several surgical scenarios and needle trajectories (straight-line or steered) in order to identify insertion points and puncture trajectories that maximise the probability of targeting success in minimally invasive surgical procedures. The potential function approach may also be applicable to straight-line anatomical needle insertion planning. Robotic Needle Manipulation: Robotic needle manipulation systems that are based on the steering and manipulation concepts introduced in this work promise to provide more precise needle placement. New steering mechanisms, model-based controllers, and robust imaging and image tracking systems will be required in order to close the loop on needle insertion in clinical applications. Such systems will complement the surgical mechanisms and telemanipulators that have recently started entering the operating room. 8.3 The Big Picture In addition to contributing models and simulation methods that will lead to advanced training systems for medical students, this thesis addresses a fundamentally important problem in percutaneous therapy, namely that of manipulating a flexible needle to hit a target in a deforming organ. Although some of the assumptions in this work, concerning tissue elasticity and linearity, may not be realistic for real organs, the approach makes an important contribution to knowledge, and lays the basis for future work. The concept of steering a flexible needle by manipulating its base breaks away from the way that we think about linear needle insertions in percutaneous therapy, and may lead to systems that are capable of high targeting accuracy in regions of the anatomy that were formerly unreachable by conventional means. The elastic tissue models and plastic phantoms that were used to approximate soft tissues provided a tractable means for developing these ideas. The Needle Manipulation Jacobian and the resulting steering and motion planning methodologies stand independent of the specific models used here, and will continue to be developed using more realistic environments in the future. In conjunction with real-time closed-loop compensation methods, based on medical imaging systems, the shortcomings of these models will become less relevant, and the methods more robust in the future. Bibliography [1] M . B . Harm, P. M . McQuillan, and G . J . Sheplock, Regional anesthesia : an atlas of anatomy and techniques. Mosby, 1996. [2] S. Nag, ed., Principles and Practice of Brachytherapy. Futura Publishing Company, 1997. [3] F . S. Azar, D . N . Metaxas, and M . D . Schnall, " A Finite Element Model of the Breast for Predicting Mechanical Deformations during Biopsy Procedures," in Proceedings of the IEEE Workshop on Mathematical Methods in Biomedical Image Analysis, pp. 38-45, 2000. [4] J . J . Bauer, J . Zeng, I. A . Sesterhenn, J . W . Moul, and S. K . M u n , "Comparison of Prostate Biopsy Protocols using 3-D Computer Simulation," in Proceedings of the Pacific Medical Technology Symposium, pp. 109-114, Aug. 1998. [5] S. Nath, Z. Chen, N . Yue, S. Trumpore, and R. Peschel, "Dosimetric effects of needle divergence in prostate seed implant using 1 2 5 I and 1 0 3 P d radioactive seeds," in Medical Physics, pp. 1058-1066, 2000. [6] J . De Andres, M . A . Reina, and A . Lopez-Garcia, "Risks of regional anaesthesia: Role of equipment - needle design, catheters," in The 17 th Annual European Society of Regional Anaesthesia Congress, Available Online: http://www.esraeurope.org/andresl4.htm, Sept. 1998. [7] P. Gorman, T . Krummel, R. Webster, M . Smith, and D . Hutchens, " A Prototype Haptic Lumbar Puncture Simulator," in Proceedings of Medicine Meets Virtual Reality, pp. 106-109, 2000. [8] S. Datta, Annual "Complications of European Society of regional Regional analgesia Anaesthesia http://www.esraeurope.org/dattal.htm, Sept. 1998. 118 and anaesthesia," Congress, in Available The 17 th Online: BIBLIOGRAPHY [9] B . 119 Finucane, nual "Complications European Society of of brachial Regional http://www.esraeurope.org/finucanl.htm, [10] D . L. Annual Brown, "What European do Society we of anaesthesia," Anaesthesia Congress, in The 17 th Available An- Online: Sept. 1998. need to Regional http://www.esraeurope.org/brownl2.htm, plexus advance regional Anaesthesia anaesthesia?," Congress, in Available The Online: 1999. [11] L . T. Kohn, J . M . Corrigan, and M . S. Donaldson, eds., To err is human : building a safer health system. National Academy Press, 2000. [12] A . Zorcolo, E . Gobbetti, P. P i l i , and M . Tuveri, "Catheter insertion simulation with combined visual and haptic feedback," in Proceedings of the First PHANToM (PURS'99), Users Research Symposium May 1999. [13] E . Gobbetti, M . Tuveri, G . Zanetti, and A . Zorcolo, "Catheter insertion simulation with coregistered direct volume rendering and haptic feedback," in Medicine Meets Virtual Reality 2000, pp. 96-98, Jan. 2000. [14] Immersion Medical, "CathSim Vascular Access Simulator." Available Online: http://www.immersion.com/products/medical/vascular.shtml. [15] M . Bro-Nielsen, J . L . Tasto, R. Cunningham, and G . L . Merril, " P r e O p r M Endoscopic Simu- lator - A PC-Based Immersive Training System for Bronchoscopy," in Medicine Meets Virtual Reality 7 (MMVR-7), 1999. [16] K . Ikuta, K . Iritani, and J . Fukuyama, "Mobile Virtual Endoscope System with Haptic and Visual Information for Non-invasive Inspection Training," in IEEE International on Robotics and Automation, Conference May 2001. [17] P. N . Brett, T. J . Parker, A . J . Harrison, T. A . Thomas, and A . Carr, "Simulation of resistance forces acting on surgical needles," in Journal of Engineering in Medicine, Institute of Mechanical Engineers, Part H, vol. 211, pp. 335-347, 1997. [18] T. Dang, T. M . Annaswamy, and M . A . Srinivasan, "Development and Evaluation of an Epidural Injection Simulator with Force Feedback for Medical Training," in Proceedings of Medicine Meets Virtual Reality, pp. 97-102, 2001. BIBLIOGRAPHY 120 [19] L . Hiemenz, D . J . McDonald, D . Stredney, and D . Sessanna, " A Physiologically Valid Simulator for Training Residents to Perform an Epidural Block," in Proceedings of the 15 th Southern Biomedical Engineering Conference, pp. 170-173, Mar. 1996. [20] L . Hiemenz, A . Litsky, and P. Schmalbrock, "Puncture Mechanics for the fnsertion of an Epidural Needle," in 21 st Annual Meeting of the American Society of Biomechanics, Available Online: http://asb-biomech.org/onlineabs/abstracts97/34/index.html, 1997. [21] L . Hiemenz, D . Stredney, and P. Schmalbrock, "Development of the Force-Feedback Model for an Epidural Needle fnsertion Simulator," in Medicine Meets Virtual Reality, pp. 272-277, 1998. [22] R. S. Haluck, W . B . Murraya, R. Webster, N . Mohler, and M . Melkonian, " A Haptic Lumbar Puncture Simulator," in Proceedings of Medicine Meets Virtual Reality, Jan. 2000. [23] S. K . Singh, M . Bostrom, D . O. Popa, and C. W . Wiley, "Design of an interactive lumbar puncture simulator with tactile feedback," in Proceedings of the IEEE International Confer- ence on Robotics and Automation, vol. 2, pp. f734-f739, f994. [24] D . O. Popa and S. K . Singh, "Creating Realistic Force Sensations in a Virtual Environment: Experimental System, Fundamental Issues and Results," in Proceedings of the IEEE International Conference on Robotics and Automation, pp. 59-63, May f 998. [25] D . Kwon, J . Kyung, S. M . Kwon, J . B . Ra, H . W . Park, H . S. Kang, J . Zeng, and K . R. Cleary, "Realistic Force Reflection in a Spine Biopsy Simulator," in Proceedings of the IEEE International Conference on Robotics and Automation, pp. 1358-1363, 2001. [26] K . B . Shimoga and P. K . Khosla, "Visual and Force Feedback to A i d Neurosurgical Probe Insertion," in 16 th International Conference of the IEEE Engineering in Medicine and Biology Society, vol. 2, pp. 1051-1052, 1994. [27] S. Miller, C. Jeffrey, J . Bews, and W . Kinsner, "Advances in the Virtual Reality Interstitial Brachytherapy System," in Proceedings of the Canadian Conference on Electrical and Computer Engineering, pp. 349-354, May 1999. BIBLIOGRAPHY 121 [28] J. Zeng, C. Kaplan, J. Bauer, J. Xuan, I. A. Sesterhenn, J. H. Lynch, M. T. Freedman, and S. K. Mun, "Optimizing prostate needle biopsy through 3-D simulation," in Proceedings of SPIE Medical Imaging, vol. 3335, pp. 488-497, Feb. 1998. [29] J. Moline, "Virtual Environments for Health Care," tech. rep., A White Paper for the Advanced Technology Program (ATP) National Institute of Standards and Technology, Oct. 1995. [30] R. M. Satava and S. B. Jones, "Current and Future Applications of Virtual Reality for Medicine," in Proceedings of the IEEE, vol. 86(3), pp. 484-489, 1998. [31] D. Sorid and S. K. Moore, "The Virtual Surgeon," in IEEE Spectrum Magazine, pp. 26-31, July 2000. [32] M. Teschner, S. Girod, and B. Girod, "Optimization Approaches for Soft-Tissue Prediction in Craniofacial Surgery Simulation," in Medical Intervention, Image Computing and Computer-Assisted pp. 1183-1190, 1999. [33] L. M. Auer, A. Radetzky, C. Wimmer, G. Kleinszig, F. Schroecker, D. P. Auer, H. Delingette, B. Davies, and D. P. Pretschner, "Visualization for Planning and Simulation of Minimally Invasive Neurosurgical Procedures," in Medical tervention, Image Computing and Computer-Assisted In- pp. 1199-1209, 1999. [34] A. Zivanovic and B. L. Davies, "A Robotic System for Blood Sampling," in IEEE on Information Technology in Biomedicine, Transactions vol. 4(1), pp. 8-14, Mar. 2000. [35] K. Khodabandehloo, P. N. Brett, and R. O. Buckingham, "Special-Purpose Actuators and Architectures for Surgery Robots," in Computer-Integrated Applications, Surgery: Technology and Clinical The MIT Press, 1996. [36] H. Saito and T. Togawa, "Detection of Puncturing Vessel Wall for Automatic Blood Sampling," in Proceedings of the First Joint BMES/EMBS Conference. Serving Humanity, Ad- vancing Technology, p. 866, Oct. 1999. [37] S. B. Solomon, A. Patriciu, M. E. Bohlman, L. R. Kavoussi, and D. Stoianovici, "Robotically Driven Interventions: A Method of Using CT Fluoroscopy without Radiation Exposure to the Physician," in Radiology, pp. 277-282, 2002. BIBLIOGRAPHY 122 [38] M . H . Loser, A New Robotic System for Visually Controlled Percutaneous Interventions under X-ray- or CT-Fluoroscopy. P h D thesis, University of Freiburg, Germany, 2002. [39] C . Simone and A . Okamura, "Haptic Modeling of Needle Insertion for Robot-Assisted Percutaneous Therapy," in Proceedings of the IEEE International Automation, Conference on Robotics and May 2002. [40] R. Alterovitz, J . Pouliot, R. Taschereau, I. J . Hsu, and K . Goldberg, "Needle Insertion and Radioactive Seed Implantation in Human Tissues: Simulation and Sensitivity Analysis," in IEEE International Conference on Robotics and Automation, 2003. [41] K . Miller, K . Chinzei, G . Orssengo, and P. Bednarz, "Mechanical properties of brain tissue in-vivo: experiment and computer simulation," in Journal of Biomechanics, vol. 33, pp. 13691376, 2000. [42] M . P. Ottensmeyer and J . K . Salisbury Jr., "7n Vivo Data Acquisition Instrument for Solid Organ Mechanical Property Measurement," in Medical Image Computing and Computer Aided Intervention (MICCAI), pp. 975-982, 2001. [43] I. Brouwer, J . Ustin, L . Bentley, A . Sherman, N . Dhruv, and F . Tendick, "Measuring In Vivo Animal Soft Tissue Properties for Haptic Mmodeling in Surgical Simulation," in Proceedings of Medicine Meets Virtual Reality (J. D . Westwood et al., ed.), pp. 69-74, IOS Press, 2001. [44] H . Nienhuys, Cutting in deformable objects. P h D thesis, Utrecht University, 2003. [45] T . Mineta, T . Mitsui, Y . Watanabe, S. Kobayashi, Y . Haga, and M . Esashi, "Batch fabricated flat meandering shape memory alloy actuator for active catheter," in Sensors and Actuators A: Physical, vol. 88, pp. 112-120, Feb. 2001. [46] S. Chin, D . Berube, D . Mody, and N . Norris, "Medical Instrument Positioning Tool and Method." US Patent Application 2002/0128636, Sept. 2002. [47] R. Slater, "Steerable Drug Delivery Device." US Patent App. 2002/0143291, Oct. 2002. [48] R. Ebrahimi, S. Okazawa, R. Rohling, and S. Salcudean, "Hand-held steerable needle device," in Medical Image Computing and Computer Aided Intervention (MICCAI), Springer-Verlag, 2003. BIBLIOGRAPHY 123 [49] S . D . Pathak, P. D. Grimm, and Y . K i m , "Pubic Arch Detection and Interference Assessment in Transrectal Ultrasound Guided Prostate Cancer Therapy." US Patent 6,027,446, Feb. 2000. [50] S. P. DiMaio and S. E . Salcudean, "Simulated Interactive Needle Insertion," in Proceedings of the 10 th IEEE Symposium on Haptic Interfaces for Virtual Environments and Teleoperator Systems, Virtual Reality, pp. 344-351, Mar. 2002. [51] S. P. DiMaio and S. E . Salcudean, "Needle Insertion Modelling and Simulation," in IEEE International Conference on Robotics and Automation, May 2002. [52] S. P. DiMaio and S. E . Salcudean, "Needle Insertion Modelling for the Interactive Simulation of Percutaneous Procedures," in Medical Image Computing and Computer Aided Intervention (MICCAI), vol. 2489, pp. 253-260, 2002. [53] S. P. DiMaio and S. E . Salcudean, "Needle Insertion Modelling and Simulation," in IEEE Transactions on Robotics and Automation: Special Issue on Medical Robotics., Oct. 2003. [54] S. P. DiMaio and S. E . Salcudean, "Needle Steering and Model-Based Trajectory Planning," in Medical Image Computing and Computer Aided Intervention (MICCAI), 2003. [55] D . Terzopoulos and K . Fleischer, "Modeling Inelastic Deformation: Viscoelasticity, Plasticity, Fracture," in Computer Graphics, vol. 22(4), pp. 269-278, A u g . 1988. [56] D . Terzopoulos and K . Fleischer, "Deformable Models," in The Visual Computer, vol. 4, pp. 306-331, Springer-Verlag, 1988. [57] D . T. Chen and D. Zeltzer, "Pump It Up: Computer Animation of a Biomechanically Based Model of Muscle Using the Finite Element Method," in Computer Graphics, vol. 26, pp. 89-98, July 1992. [58] M . Paule, C. Gascuel, and M . Desbrun, "Animation of Deformable Models Using Implicit Surfaces," in IEEE Transactions on Visualization and Computer Graphics, vol. 3(1), pp. 3950, 1997. [59] S. Coquillart, "Extended Free-From Deformation: A Sculpturing Tool for 3D Geometric Modeling," in Computer Graphics, vol. 24(4), pp. 187-196, Aug. 1990. BIBLIOGRAPHY 124 [60] S. Sclaroff and A . Pentland, "Generalized Implicit Functions for Computer Graphics," in Computer Graphics, vol. 25(4), pp. 247-250, July 1991. [61] X . Provot, "Deformation Constraints in a Mass-Spring Model to Describe Rigid Cloth Behavior," in Proceedings of the Graphics Interface Conference, pp. 147-154, 1995. [62] D . L . James and D . K . Pai, "ArtDefo: Accurate Real Time Deformable Objects," in Computer Graphics - Proceedings of SIGGRAPH '99, pp. 65-72, 1999. [63] D . L . James and D . K . Pai, "Multiresolution Green's Function Methods for Interactive Simulation of Large-scale Elastostatic Objects," in A CM Transactions on Graphics, vol. 22, pp. 4782, 2003. [64] D . L . James and D. K . Pai, "Real Time Simulation of Multizone Elastokinematic Models," in IEEE International Conference on Robotics and Automation, May 2003. [65] D . L . James and D. K . Pai, "DyRT: Dynamic Response Textures for Real Time Deformation Simulation with Graphics Hardware," in ACM Transactions on Graphics, vol. 21, pp. 582-585, 2002. [66] H . Delingette, "Initialization of Deformable Models from 3D Data," in Proceedings of the IEEE International Conference on Computer Vision, pp. 311-316, 1998. [67] D . E . Altobelli, R. Kikinis, J . B . Mulliken, H . C. W . Lorensen, and F . Jolesz, "ComputerAssisted Three-Dimensional Planning in Craniofacial Surgery," in Plastic and Reconstructive Surgery, vol. 92(4), pp. 576-585, 1993. [68] H . Delingette, G . Subsol, S. Cotin, and J . Pignon, " A Craniofacial Surgery Simulation Testbed," Tech. Rep. 2199, Institut National de Recherche en informatique et en automatique, Feb. 1994. [69] J . X i a , F . Q i , W . Yuan, D . Wang, W . Qiu, Y . Sun, Y . Huang, G . Shen, and H . Wu, "Computer Aided Simulation System for Orthognatic Surgery," in Eighth IEEE Symposium on ComputerBased Medical Systems, pp. 237-244, 1995. [70] E . Keeve, S. Girod, P. Pfeifle, and B . Girod, "Anatomy-Based Facial Tissue Modeling Using the Finite Element Method," in Proceedings of the IEEE Visualization Conference, pp. 21-28, 1996. BIBLIOGRAPHY 125 [71] A . Sarti, R. Gori, and C. Lamberti, " A Physically Based Model to Simulate Maxillo-Facial Surgery from 3D C T Images," in Future Generations Computer Systems, pp. 217-221, 1999. [72] F . Schutuser, J . V . Cleyenbreugel, J . Schoenaers, G . Marchal, and P. Suetens, " A Simulation Environment for Maxillofacial Surgery Including Soft Tissue Implications," in Medical Image Computing and Computer-Assisted Intervention, pp. 1210-1217, 1999. [73] A . Jodicke, W . Deinsberger, H.Erbe, A . Kriete, and D . K . Boker, "Intraoperative ThreeDimensional Ultrasonography: A n Approach to Register Brain Shift using Multidimensional Image Processing," in Minimally Invasive Neurosurgery, vol. 41, pp. 13-19, 1998. [74] H . Dickhaus, K . Ganser, A . Staubert, M . M . Bonsanto, C. R. Wirtz, V . M . Tronier, and S. Kunze, "Quantification of brain shift effects by MR-imaging," in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology, pp. 491-494, 1997. [75] D . G . Gobbi, R. M . Comeau, and T. M . Peters, "Ultrasound Probe Tracking for Real-Time Ultrasound/MRI Overlay and Visualization of Brain Shift," in Medical Image Computing and Computer-Assisted Intervention, pp. 920-927, Springer-Verlag, 1999. [76] R. M . Comeau, A . Fenster, and T. M . Peters, "Integrated M R and Ultrasound Imaging for Improved Image Guidance in Neurosurgery," in SPIE Conference on Image Processing, pp. 747-754, Feb. 1998. [77] N . Hata, William M . Wells III, M . Halle, S. Nakajima, P. Viola, R. Kikinis, and F. A . Jolesz, "Image Guided Microscopic Surgery System Using Mutual-Information Based Registration," in Visualization in Biomedical Computing, vol. 1131, pp. 317-326, Springer-Verlag, 1996. [78] A . Hagemann, K . Rohr, H . S. Stiehl, U . Spetzger, and J . M . Gilsbach, "Intraoperative Image Correction using a Biomechanical Model of the Human Head with Different Material Properties," in Mustererkennung (W. Forstner, J . M . Buhmann, A . Faber, and P. Faber, eds.), pp. 223-231, Springer-Verlag, 1999. [79] M . I. Miga, K . D . Paulsen, F . E . Kennedy, A . Hartov, and D . W . Roberts, "Model-Updated Image Guided Neurosurgery Using the Finite Element Method: Incorporation of the Falx BIBLIOGRAPHY 126 Cerebri," in Medical Image Computing and Computer-Assisted Intervention, pp. 900-909, 1999. [80] K . Miller, "Constitutive model of brain tissue suitable for finite element analysis of surgical procedures," in Journal of Biomechanics, vol. 32, pp. 531-537, 99. [81] A . A . Young, D . L . Kraitchman, and L . Axel, "Deformable Models for Tagged M R Images: Reconstruction of Two- and Three-Dimensional Heart Wall Motion," in Proceedings of the IEEE Workshop on Biomedical Image Analysis, pp. 317-323, 1994. [82] A . A . Amini, Y . Chen, R. W . Curwen, V . Mani, and J . Sun, "Coupled B-snake grids and constrained thin-plate splines for analysis of 2-D tissue deformations from tagged M R I , " in IEEE Transactions on Medical Imaging, vol. 17, pp. 344-356, 1998. [83] A . A . Amini, Y . Chen, and D . Abendschein, "Comparison of Land-Mark-Based and CurveBased Thin-Plate Warps for Analysis of Left-Ventricular Motion from Tagged M R I , " in Medical Image Computing and Computer-Assisted Intervention, pp. 498-507, 1999. [84] J . Park, D . Metaxas, and L . Axel, "Volumetric Deformable Model with Parameter Functions: A New Approach to the 3D Motion Analysis of the L V from M R I - S P A M M , " in IEEE International Conference on Computer Vision, pp. 700-705, 1995. [85] E . Bardinet, L . D . Cohen, and N . Ayache, "Analyzing the deformation of the left ventricle of the heart with a parametric deformable model," Tech. Rep. 2797, Institut National de Recherche en informatique et en automatique, Feb. 1996. [86] X . Papademetris, A . J . Sinusas, D . P. Dione, and J . S. Duncan, "3D Cardiac Deformation from Ultrasound Images," in Medical Image Computing and Computer-Assisted Intervention, pp. 420-429, 1999. [87] J . M . Guccione and A . D . McCulloch, "Finite element models of ventricular mechanics," in Theory of Heart: Biomechanics, Biophysics and Nonlinear Dynamics of Cardiac Function (L. Glass and A . D . M . P. J . Hunter, eds.), pp. 121-145, Springer-Verlag, 1991. [88] J . Calisse, A . Rohlmann, and G . Bergmann, "Estimation of trunk muscle forces using the finite element method and in vivo loads measured by telemetrized internal spinal fixation devices," in Journal of Biomechanics, vol. 32, pp. 727-731, Elsevier, 1999. BIBLIOGRAPHY 127 [89] M . A . Schill, C. Wagner, M . Hennen, Hans-Joachim Bender, and R. Manner, "EyeSi - A Simulator for Intra-ocular Surgery," in Medical Image Computing and Computer-Assisted Intervention, pp. 1166-1174, 1999. [90] E . Chen and B . Marcus, "Force feedback for surgical simulation," in Proceedings of the IEEE, vol. 86, pp. 524-530, Mar. 1998. [91] G . Burdea, G . Patounakis, V . Popescu, and R. E . Weiss, "Virtual Reality Training for the Diagnosis of Prostate Cancer," in IEEE International cations (VRAIS Symposium on Virtual Reality Appli- '98), pp. 190-197, Mar. 1998. [92] http://www.sensable.com. SensAble Technologies. Cambridge, M A . [93] V . Popescu, G . Burdea, and M . Bouzit, "Virtual Reality Simulation Modeling for a Haptic Glove," in Computer Animation, May 1999. [94] M . Dinsmore, N . Langrana, G . Burdea, and J . Ladeji, "Virtual Reality Simulation for Palpation of Subsurface Tumors," in Proceedings of the IEEE Virtual Reality Annual International Symposium, 1997. [95] N . A . Langrana, G . Burdea, K . Lange, D . Gomez, and S. Deshpande, "Dynamic Force Feedback in a Virtual Knee Palpation," in Artificial Intelligence in Medicine, pp. 321-333, 1994. [96] P. S. Wellman and R. D . Howe, "Modeling Probe and Tissue Interaction for Tumor Feature Extraction," in ASME Summer Bioengineering Conference, June 1997. [97] P. S. Wellman and R. D . Howe, "Extracting Features from Tactile Maps," in Medical Image Computing and Computer-Assisted Intervention, p. 1133:1142, 1999. [98] J . Yan, P. K . Scott, and R. S. Fearing, "Inclusion probing: Signal detection and haptic playback of 2D F E M and experimental data," in 8 th Virtual Environment Symposium on Haptic Interfaces for and Teleoperator Systems, Proceedings of the ASME IMECE DSC-5, Nov. 1999. [99] J . C. Berrios and P. C. Pedersen, "Deep Vein Thrombosis Detection," in Ultrasonics Symposium, pp. 1567-1570, 1990. BIBLIOGRAPHY 128 [100] T. J . Barloon, G . R. Bergus, and J . E . Seabold, "Diagnostic Imaging of Lower Limb Deep Venous Thrombosis," in American Family Physician, vol. 56(3), Sept. 1997. [101] J . R. Harding, "Investigating Deep Venous Thrombosis with Infrared Imaging," in IEEE Engineering in Medicine and Biology, pp. 43-46, July 1998. [102] J . Guerrero, S. Salcudean, J . A . McEwen, B . A . Masri, and S. Nicolaou, "Measurement-Based Deep Venous Thrombosis Screening System," in Medical Image Computing and Computer Aided Intervention, 2003. [103] D . d'Aulignac, M . C . Cavu§oglu, and C . Laugier, "Modeling the Dynamics of the Human Thigh for Realistic Echographic Simulator with Force Feedback," in Medical Image Computing and Computer-Assisted Intervention, pp. 1191-1198, 1999. [104] A . J . Hodgson, J . G . Person, S. E . Salcudean, and A . G . Nagy, "Assessing Potential Benefits of Enhanced Dexterity in Laparoscopic Surgery," in Medical Image Analysis (Oxford University Press), 1999. [105] S. Cotin, H . Delingette, M . Bro-Nielsen, N . Ayache, J . M . Clement, V . Tassetti, and J. Marescaux, "Geometric and Physical Representations for a Simulator of Hepatic Surgery," in Medicine Meets Virtual Reality IV, Interactive Technology and the New Paradigm for Healthcare, pp. 139-151, 10S Press, Jan. 1986. [106] S. Cotin and H . Delingette, "Real-time Surgery Simulation with Haptic Feedback using F i nite Elements," in Proceedings of the 1998 IEEE International Automation, Conference on Robotics and pp. 3739-3744, May 1998. [107] S. Cotin, H . Delingette, and N . Ayache, "Real-time elastic deformations of soft tissues for surgery simulation," Tech. Rep. 351 f, fnstitut National de Recherche en informatique et en automatique, Oct. 1998. [108] S. Cotin, H . Delingette, and N . Ayache, "Efficient Linear Elastic Models of Soft Tissues for Real-Time Surgery Simulation," Tech. Rep. 3510, fnstitut National de Recherche en informatique et en automatique, Oct. 1998. BIBLIOGRAPHY 129 [109] S. Cotin, H . Delingette, and N . Ayache, "Real-Time Elastic Deformations of Soft Tissues for Surgery Simulation," in IEEE Transactions on Visualization and Computer Graphics, vol. 5(f), pp. 62-73, 1999. [110] F . Boux de Casson and C. Laugier, "Modeling the Dynamics of a Human Liver for a Minimally Invasive Surgery Simulator," in Medical Image Computing and Computer-Assisted Interven- tion, pp. 1156-1165, 1999. [Ill] J . B . A . Maintz and M . A . Viergever, " A Survey of Medical Image Registration," in Medical Image Analysis, vol. 2(1), pp. 1-36, 1998. [112] P. A . van den Elsen, E . D . Pol, and M . A . Viergever, "Medical Image Matching - A Review with Classification," in IEEE Engineering in Medicine and Biology, pp. 26-39, Mar. 1993. [113] Lisa Gottesfeld Brown, " A Survey of fmage Registration Techniques," in ACM Computing Surveys, vol. 24(4), pp. 325-376, Dec. 1992. [114] J . West, J . M . Fitzpatrick, M . Y . Wang, B . M . Dawant, C. R. Maurer Jr., R. M . Kessler, R. J . Maciunas, C. Barillot, D . Lemoine, A . Collignon, F . Maes, P. Suetens, D . Vandermeulen, P. A . van den Elsen, S. Napel, T. S. Sumanaweera, B . Harkness, P. F . Hemler, D . L . G . H i l l , D. J . Hawkes, C. Studholme, J . B . A . Maintz, M . A . Viergever, G . Malandain, X . Pennec, M . E . Noz, G . Q. Maguire Jr., M . Pllack, C. A . Pelizzari, R. A . Robb, D . Hanson, and R. P. Woods, "Comparison and Evaluation of Retrospective Intermodality Brain Image Registration Techniques," in Journal of Computer Assisted Tomography, vol. 21, pp. 554-566, 1997. [115] D . L . G . H i l l , C. Studholme, and D . J . Hawkes, "Voxel Similarity Measures for Automated fmage Registration," in Proceedings of SPIE, vol. 2359, pp. 205-216, 1994. [116] J . C. Gee, C. Barillot, L . L . Briquer, D . R. Haynor, and R. Bajcsy, "Matching Structural Images ofthe Human Brain Using Statistical and Geometrical Image Features," in Proceedings of SPIE, vol. 2359, pp. 191-204, 1994. [117] A . Gueziec, K . Wu, B . Williamson, P. Kazanzides, R. V . Vorhis, and A . Kalvin, "Exploiting 2-D to 3-D Intra-operative Image Registration for Qualitative Evaluations and Post-operative Simulations," in Medical Image Computing and Computer-Assisted Springer-Verlag, 1999. Intervention, pp. 820-83f, BIBLIOGRAPHY 130 [118] T. Gaens, F . Maes, D . Vandermeulen, and P. Seutens, "Non-rigid Multimodal Image Registration Using Mutual Information," in Medical Image Computing and Computer-Assisted Intervention, vol. 1496, pp. 1099-1106, Springer-Verlag, 1998. [119] J . Feldmar, J . Declerck, G . Malandain, and N . Ayache, "Extension of the I C P Algorithm to Nonrigid Intensity-Based Registration of 3D Volumes," in Computer Vision and Image Understanding, vol. 66, pp. 193-206, May 1997. [120] R. Bajcsy and S. Kovacic, "Multiresolution Elastic Matching," in Computer Vision, Graphics and Image Processing, vol. 46, pp. 1-21, 1989. [121] B . M . Dawant, S. L . Hartmann, and S. Gadamsetty, "Brain Atlas Deformation in the Presence of Large Space-occupying Tumors," in Medical Image Computing and Computer-Assisted Intervention, pp. 589-596, Springer-Verlag, 1999. [122] T. Mclnerney and D . Terzopoulos, "Deformable Models in Medical Image Analysis," in Proceedings of the Workshop on Mathematical Methods in Biomedical Image Analysis, pp. 171— 180, 1996. [123] X . Pennec, P. Cachier, and N . Ayache, "Understanding the 'Demon's Algorithm': 3D Nonrigid Registration by Gradient Descent," in Medical Image Computing and Computer-Assisted Intervention, pp. 597-605, Springer-Verlag, 1999. [124] M . Ferrant, S. Warfield, C. R. G . Guttmann, R. V . Mulkern, F . A . Jolesz, and R. Kikinis, "3D Image Matching Using a Finite Element Based Elastic Deformation Model," in Medical Image Computing and Computer-Assisted Intervention, pp. 202-209, 1999. [125] A . Hagemann, K . Rohr, H . S. Stiehl, U . Spetzger, and J . M . Gilsbach, " A Biomechanical Model of the Human Head for Elastic Registration of MR-Images," in Bildverarbeitung fuer fie Medizin (H. Evers, G . G . dn T. Lehmann, and H . P. Meinzer, eds.), pp. 44-48, SpringerVerlag, 1999. [126] A . Hagemann, K . Rohr, H . S. Stiehl, U . Spetzger, and J . M . Gilsbach, "Nonrigid Matching of Tomographic Images Based on a Biomechanical Model of the Human Head," in Medical Imaging - Image Processing, pp. 583-592, 1999. BIBLIOGRAPHY 131 [127] C . Davatzikos, D . Shen, A . Mohamed, and S. K . Kyriacou, " A Framework for Predictive Modeling of Anatomical Deformations," in IEEE Transactions on Medical Imaging, vol. 20, pp. 836-843, 2001. [128] S. F . Gibson, "3D ChainMail: a Fast Algorithm for Deforming Volumetric Objects," in Proceedings of the Symposium on Interactive 3D Graphics, Apr. 1997. [129] S. Gibson, J . Samosky, A . Mor, C . Fyock, E . Grimson, T . Kanade, R. Kikinis, H . Lauer, N . McKenzie, S. Nakajima, H . Ohkami, R. Osborne, and A . Sawada, "Simulating Arthroscopic Knee Surgery using Volumetric Object Representations, Real-Time Volume Rendering and Haptic Feedback," in Proceedings of CVRMed II and MRCAS III, Mar. 1997. [130] P. Meseure and C . Chaillou, "Deformable Body Simulation with adaptive subdivision and cuttings," in Proceedings of the WSCG, pp. 361-370, Feb. 1997. [131] S. E . Salcudean and T . D . Vlaar, "On the Emulation of Stiff Walls and Static Friction W i t h a Magnetically Levitated Input/Output Device," in Journal of Dynamic Systems, Measurement and Control, vol. 119, pp. 127-132, Mar. 1997. [132] A . Radetzky, A . Nurnberger, M . Teistler, and D . P. Pretscher, "Elastodynamic shape modeling in virtual medicine," in Shape Modelling and Applications, Proceedings of the International Conference on Shape Modeling, pp. 172-178, Mar. 1999. [133] R. J . Lapeer and R. W . Prager, "Finite Element Model of a Fetal Skull Subjected to Labour Forces," in Medical Image Computing and Computer-Assisted Intervention, pp. 1143-1155, 1999. [134] M . Bro-Nielsen, "Finite Element Modeling in Surgery Simulation," in Proceedings of the IEEE, vol. 86, pp. 490-503, 1998. [135] S. De and M . A . Srinivasan, "Thin Walled Models for Haptic and Graphical Rendering of Soft Tissues in Surgical Simulations," in Medicine Meets Virtual Reality (J. D . W . et al., ed.), pp. 94-99, fOS Press, 1999. [136] Y . C . Fung, Biomechanics - Mechanical Properties of Living Tissues. Springer-Verlag, 1993. 132 BIBLIOGRAPHY [137] H . Kataoka, T . Washio, M . Audette, and K . Mizuhara, " A Model for Relations between Needle Deflection, Force, and Thickness on Needle Penetration," in Medical Image Computing and Computer Aided Intervention, pp. 966-974, 2001. [138] H . Kataoka, T . Washio, K . Chinzei, K . Mizuhara, C . Simone, and A . M . Okamura, "Measurement of the T i p and Friction force Acting on a Needle during Penetration," in Medical Image Computing and Computer Aided Intervention, pp. 216-223, 2002. [139] J . Latombe, Robot Motion Planning. Kluwer Academic, 1991. [140] J . Barraquand and J . C . Latombe, "On Non-Holonomic Mobile Robots and Optimal Maneuvering," in Revue dTntelligence Artificielle, vol. 3(2), pp. 77-103, 1989. [141] O. Khatib, "Real-Time Obstacle Avoidance for Manipulators and Mobile Robots," in International Journal of Robotics Research, pp. 90-98, 1996. [142] M . R. Sirouspour, S. P. DiMaio, S. E . Salcudean, P. Abolmaesumi, and C . Jones, "Haptic Interface Control - Design Issues and Experiments with a Planar Device," in Proceedings of the IEEE International Conference on Robotics and Automation, Apr. 2000. [143] J . Heikkila, and O. Silven, " A Four-step Camera Calibration Procedure with Implicit Image Correction," in Proceedings of the IEEE Computer Society Conference on Computer and Pattern Recognition (CVPR'97), Vision pp. 1106-1112, 1997. [144] A . E . Kerdok, S. M . Cotin, M . P. Ottensmeyer, A . Galea, R. D . Howe, and S. L . Dawson, "Truth Cube: Establishing Physical Standards for Real Time Soft Tissue Simulation," in International Workshop on Deformable Modeling and Soft Tissue Simulation, pp. 65-72, Nov. 2001. [145] K . R. Castleman, Digital Image Processing. Prentice Hall, 1996. [146] D . H . Ballard and C . M . Brown, Computer Vision. Prentice Hall, 1982. [147] Y . Zhuang, Real-time Simulation of Physically Realistic Global Deformations. P h D thesis, University of California at Berkeley, Fall 2000. [148] Y . Zhuang and J . Canny, "Haptic Interaction with Global Deformations," in Proceedings of the IEEE International Conference on Robotics and Automation, pp. 2428-2433, A p r . 2000. BIBLIOGRAPHY 133 [149] G . H . Golub and C . F . V . Loan, Matrix computations. Johns Hopkins University Press, 3 rrf ed., 1996. [150] H . Nienhuys and A . F . van der Stappen, "Interactive needle insertions in 3D nonlinear material," Tech. Rep. UU-CS-2003-019, Institute of Information and Computing Sciences, Utrecht University, 2003. www.cs.uu.nl. [151] D . L . James, Multiresolution Green's Function Methods for Interactive Simulation of Large- scale Elastostatic Objects and other Physical Systems in Equilibrium. P h D thesis, Institute of Applied Mathematics, University of British Columbia, 200L [f 52] G . Debunne, M . Desbrun, M.-P. Cani, and A . H . Barr, "Dynamic Real-Time Deformations Using Space k, Time Adaptive Sampling," in SIGGRAPH 2001, Computer Graphics Proceedings (E. Fiume, ed.), pp. 31-36, A C M Press / A C M S I G G R A P H , 2001. [153] P. L . Tallec and M . Vidrascu, "Efficient solution of mechanical and biomechanical problems by domain decomposition," Numerical linear algebra with applications, vol. 6, no. 7, pp. 599-616, 1999. [154] X . Wu, M . S. Downes, T . Goktekin, and F . Tendick, "Adaptive Nonlinear Finite Elements for Deformable Body Simulation Using Dynamic Progressive Meshes," in Proceedings of EUROGRAPHICS (A. Chalmers and T . - M . Rhyne, eds.), vol. 20, pp. 349-358, Blackwell Publishing, 2001. [155] A . B . Mor and T . Kanade, "Modifying Soft Tissue Models: Progressive Cutting with Minimal New Element Creation," in Medical Image Computing and Computer-Assisted (MICCAI), Intervention pp. 598-607, Springer-Verlag, 2000. [156] J . Ophir, S. K . Alam, B . Garra, F . Kallel, E . Konofagou, T. Krouskop, and T . Varghese, "Elastography: ultrasonic estimation and imaging of the elastic properties of tissues," in Proceedings of the Institute of Mechanical Engineers, Journal of Engineering in Medicine, vol. 213, pp. 203-233, 1999. [157] K . R. Nightingale, M . L . Palmeri, R. W . Nightingale, and G . E . Trahey, "On the feasibility of remote palpation using acoustic radiation force," in The Journal of the Acoustical Society of America, vol. ff0(f), pp. 625-634, July 2001. BIBLIOGRAPHY 134 [158] S. P. DiMaio, S. E . Salcudean, and M . R. Sirouspour, "Haptic interaction with a planar environment," in 9 th Symposium on Haptic Interfaces for Virtual Environments Systems, International and Teleoperator Mechanical Engineering Congress and Exposition, vol. DSC-Vol. 69-2, pp. 1223-1230, 2000. [159] S. Salcudean and L . Stocco, "Isotropy and Actuator Optimization in Haptic Interface Design," in Proceedings of the IEEE International Conference on Robotics and Automation, Apr. 2000. [160] D . Constantinescu, I. Chau, S. P. DiMaio, L . Filipozzi, S. E . Salcudean, and F . Ghassemi, "Haptic rendering of planar rigid-body motion using a redundant parallel mechanism ," in Proceedings of the IEEE International Conference on Robotics and Automation, A p r . 2000. [161] M . W . Spong and M . Vidyasagar, Robot Dynamics and Control. Wiley, 1989. [162] P. Hacksel and S. Salcudean, "Estimation of Environment Forces and Rigid-Body Velocities using Observers," in Proceedings of the IEEE International Automation, Conference on Robotics and (San Diego, C A ) , pp. 931-936, May 1994. [163] S. Nicosia and P. Tomei, "Robot control by using only joint position measurements," in IEEE Transactions on Automatic Control, pp. 1058-1061, Sept. 1990. [164] S. Salcudean and T. Vlaar, "On the emulation of stiff walls and static friction with a magnetically levitated input/output device," in Proceedings of the 1994 International Mechanical Engineering Congress and Exposition, pp. 303-309, 1994. [165] S. E . Salcudean, "Control for Teleoperation and Haptic Interfaces," in Lecture Notes in Control and Information Sciences 230 - Control Problems in Robotics and Automation (B. Siciliano and K . P. Valavanis, eds.), pp. 51-65, Springer-Verlag, 1997. [166] R. J . Adams, M . R. Moreyra, and B . Hannaford, "Stability and Performance of Haptic Displays: Theory and Experiments," in Proceedings of the ASME, Dynamic Systems and Control Division, (Anaheim, C A ) , pp. 227-234, 1998. [167] R. Adams and B . Hannaford, " A Two-Port Framework for the Design of Unconditionally Stable Haptic Interfaces," November 1998. in Proceedings of IROS 98, (Victoria, Canada), pp. 1254-1259, BIBLIOGRAPHY 135 [168] T. Yoshikawa and H . Uead, "Construction of Virtual World Using Dynamics Modules and Interaction Modules," in Proceedings of the 1996 IEEE International and Automation, Conference on Robotics (Minneapolis, Minnesota), pp. 2358-2364, 1996. [169] A . Nahvi, D . Nelson, J . Hollerbach, and D . Johnson, "Haptic Manipulation of Virtual Mechanisms from Mechanical C A D Designs," in Proceedings of the 1998 IEEE International ference on Robotics and Automation, (Leuven, Belgium), pp. 375-380, May 1998. [170] D . Lawrence, "Stability and Transparency in Bilateral Teleoperation," in IEEE on Robotics and Automation, Con- Transactions vol. 9(5), pp. 624-637, Oct. 1993. [171] S. Salcudean, M . Wong, and R. Hollis, "Design and Control of a Force-Reflecting Teleoperation System with Manetically Leviated Master and Wrist," in IEEE Transactions on Robotics and Automation, vol. 11, Dec. 1995. [172] D . C . Brown, "Close Range Camera Calibration," in Photogrammetric Engineering and Re- mote Sensing, vol. 37(8), pp. 855-866, 1971. [173] J . Bouguet, "Camera Calibration Toolbox for Matlab." Available Online: http://www.vision.caltech.edu/bouguetj/calib_doc/index.html. [174] P. L . Gould, Introduction to Linear Elasticity. Springer-Verlag, 1994. [175] Y . C . Fung, A First Course in Continuum Mechanics. Prentice Hall, 1977. [176] M . Bro-Nielsen, Medical Image Registration and Surgery Simulation. P h D thesis, Technical University of Denmark, July 1997. [177] C . Cuvelier, A . Segal, and A . A . van Steenhoven, Finite Element Methods and Navier-Stokes Equations. D . Reidel Publishing Company, 1986. [178] E . P. Popov, Engineering Mechanics of Solids. Prentice Hall, 1990. [179] R. E . White, An Introduction to the Finite Element Method with Applications to Nonlinear Problems. John Wiley and Sons, 1985. [180] M . D . Gunzburger, "Mathematical Aspects of Finite Element Methods for Incompressible Viscous Flows," in Finite Elements - Theory and Application, Proceedings of the ICASE BIBLIOGRAPHY Finite Element Theory and Application 136 Workshop (D. L . Dwoyer, M . Y . Hussaini, and R . G . Voigt, eds.), ch. 6, pp. 124-150, Springer-Verlag, 1986. Appendix A The Planar Pantograph for M a n i p u l a t i o n a n d The Mechanism Haptics planar pantograph mechanism is used as both a needle manipulator and haptic interface for needle insertion experiments and interactive simulation, as shown in Figure A.l. Figure A.l: This appendix The planar pantograph mechanism for 3-DOF manipulation and haptics. describes the mechanism, as well as sensing, actuation and control issues that are relevant to this work. For further resources related to this mechanism, see [142,158-160]. The manipulator mechanism has three degrees of freedom allowing for planar translation and unlimited rotation about a single axis. Two identical five-bar linkages (pantographs) are mounted 137 A.l Mechanism Actuation and Dynamics 138 parallel to one another, and their endpoints are coupled by a linkage that allows the end-effector to rotate. Each pantograph comprises four rigid carbon fibre links and five joints, as illustrated in Figure A.2. Joint angles 6\ and 62 are measured by digital optical encoders with a resolution y d Figure A.2: Pantograph dimensions and parameters. of 0.0225 degrees (16,000 counts per revolution), and are actuated by D C motors. The remaining joints are passive. A.l Mechanism Actuation and Dynamics Four 90W D C motors drive the active pantograph joints, and are considered to be torque sources. Each motor is driven by a constant current source, with the assumption that motor torque r is proportional to armature current i : a T = ki t a , (A.l) where kt is called the motor torque constant, which is specified by the manufacturer, and verified to be 0 . 0 5 2 5 N m A . The D C motors have a continuous current rating of 2A, and a peak rating _1 exceeding 10A. Mechanism dynamics are required for control purposes. The equations of motion for a single pantograph mechanism are derived using the Euler-Lagrange approach [161], and are expressed in A.2 Force Observer 139 joint variables 0 = [9\ # 2 ] as follows: r D {0)0 p C (0,0)0 = + p T P - J (A.2) ; where r is the vector of applied actuator torques, J is the pantograph manipulator Jacobian and p p i is the external force applied to the pantograph end-point P . Mass and Christoffel matrices, D p e p and C , are also present. p The equations of motion describing the workspace dynamics of two coupled pantographs are expressed in the coordinate space: M x c + C x = { c c c h + T = C f/j + u , (A.3) where M and C are mass and Christoffel matrices of the coupled mechanism, x c c c is a vector of interface handle coordinates [x y cx ], f^ is the external force applied at P , r is a vector of torques c c c e c applied at the active joints and J is the manipulator Jacobian. c A.2 Force Observer External forces applied at P e are estimated from joint angle measurements and applied motor torques by means of a force observer (there is no physical force sensor present). Given an accurate dynamic model, as well as measured joint angles and applied actuator forces, the system states (angular velocities at the joints) and unknown external disturbances (external force applied at P ) e can be computed, as indicated in Figure A.3. This approach has been demonstrated for a single free body [162], but is shown here to be applicable to a parallel mechanism (a single pantograph mechanism), by using a simplified Nicosia Observer [163]: 0 p — Q + k0 0 V — Dp{6) (~Cp(6,0 )0p v v p 1 p + kpOp + T ) , P (A.4) where 0 and 6 are angular position and velocity states, 0 is the position state estimation error P V P (d - 0 ), k and k are state feedback gains. B y substituting (A.4) into the pantograph equations P v p of motion (A.2), the error dynamics become: 140 A.3 Closed-Loop Position Control Controller -6 Pantograph Force/Velocity Observer (O. Figure A.3: Force observation using only applied actuator torques and measured joint angles. D (0)0 p p + [C {0,0) p + + C {0,0 ) p P k D {0)\ v p 0 + k 0p p = p J f = r T e At steady-state, the effective joint torques due to applied hand forces angle estimation errors by a simple stiffness relationship, T env = k0. p p e T env env . (A.5) are related to joint The observer also provides estimates of joint velocities: ^estimated = @p @v + k 0 = v p . (A.6) The selection of observer state feedback gains k and k is made such that the error 0 , and p v p consequently the force estimate, converge as quickly as possible. Their magnitudes are limited by observer stability considerations [164], convergence rate and noise sensitivity. Figure A.4 shows a comparison of external forces predicted off-line by an inverse dynamic model, with those estimated by the on-line force observer. The observer clearly tracks the applied force closely. Force estimates for each pantograph mechanism are combined in order to determine forces and torques applied to the interface end-effector. A.3 Closed-Loop Position Control During probing and needle insertion experiments, the planar pantograph mechanism is controlled by a simple coordinate-space PID controller, as illustrated in Figure A.5. Each axis is controlled A.3 Closed-Loop Position Control I 141 1 1 i 0.1 0.05 \ / \ / \ /~ 0 I 0.05 I I ! i i 0.5 i 4.5 Figure A.4: Force observer performance. Planar Pantograph Interface 9- s ! PID Controller Figure A.5: PID controller for needle trajectory control. independently. The controller gains for translational degrees of freedom are: Kp = 800Nm , KD = lOiVsm -1 K and Ki = 7.8 x lO~ Nm~ s~ , while for end-effector orientation, they are: 3 = INmrad- , K 1 P 1 = 0.005JVmsracT and K, = 3.9 x 10~ Nmrad' s~ . 1 D l tween the desired end-effector position 6 1 l The error be- and the resolved end-effector position Xh drives the PID A.4 Haptic Control Architecture 142 controller, which produces a command force/torque u in the coordinate space. Coordinate-space wrench commands are transformed to joint-space torques, as described in Section A.5. A.4 Haptic Control Architecture Haptically rendered interactive virtual environments resemble teleoperation systems in which the haptic interface and the virtual slave are almost always both kinematically and dynamically dissimilar, resulting in some difficulty realizing transparent interaction between the user and the virtual environment [f65]. For haptic displays, impedance and admittance simulations have been proposed [166—169]. Impedance displays pass sensed hand positions to the dynamic simulator, and return forces from the simulated environment—a two-channel approach. The transmission of positions and forces in both directions between master and slave has been found to be important for achieving high performance in teleoperation systems [165,170]. The general teleoperation architecture is shown in Figure A.6, where Fh and F are hand and e environment forces; V/j and V are master and slave velocities; and Z , e m Z , Z^ and Z are master, s e slave, hand and environment impedances, respectively [f70]. C\ and C 4 are position channel controllers; whereas, C2 and C 3 are force channel controllers. Finally, C m c Operator Master and C are local master and s 1 Communication Channel Slave i Environment Figure A.6: Four-channel architecture with local force feedback. slave position controllers. In a haptic simulation, the slave and the environment are virtual, and A.4 Haptic Control Architecture 143 the dynamic simulator replaces Z and Z . s C is a local controller for the virtual slave, designed e s independently of the dynamic simulator. A four-channel coupling between the haptic interface and dynamic simulation allows the interface to behave as a force sensor, or as a position sensor, depending upon the impedance of the virtual environment, and is therefore a hybrid of the two traditionally adopted approaches [165]. A n impedance controller shapes the pantograph interface dynamics Z m to match that of the virtual needle. The desired equations of motion have the following form: M d x + c Bk where M is the desired mass matrix, B d d m + i_dx c (A.7) = t + f, h m is the desired damping matrix, K d matrix and f c d is the desired stiffness is an external control command i.e., the sum of control commands by C , m C2 and C4, as shown in Figure A.6. B y combining (A.7) and (A.3), the following impedance control law is established: u = (M MJ l c - I){ h + M M^{ c + (C - M MJ B )x l m c c d c - M M^K x c d c , (A.8) where the hand force, or equivalently the acceleration x , is required in order to alter the apparent c mass of the device (M d ^ M ). c The observed hand force (Section A.2) is used to synthesise the control law. The following controller parameters were used: = 35 + C = 35 + Ci = c = 1 = 1 C m s 2 C7 4 = 100 s C s ~C m (A.9) A.4 Haptic Control 144 Architecture The environment and hand forces are fed forward with unity gains C2 and C 3 , and second-order position controllers act through the position channels (C\ and C4) and local controllers ( C and m C ). s An adaptive damping term B d a (K~h = 60s/m, Bi m n = K~h\f e p \ + Bmin is added to C in order to reduce chattering s = Okg/s and f is the environment force) [171]. e The advantage of the four-channel approach is that it provides a clear and general methodology for interactive virtual environment design. The explicit modelling of a virtual slave environment means that complex needle insertion mechanics can be incorporated independently of the haptic interface or its associated control system. This is in contrast with traditional approaches in which slave and environment dynamics are often implied within the haptic interface control system itself. The resulting control architecture for the interactive needle simulator is illustrated in Figure A.7. Low-level Control velocity & force Force and Velocity Observer Teleoperation Control Block Low-level Sensing and Actuation velocity & force Master Interface Virtual World State Collision Detection Contact and Slave Dynamics Slave Environment Figure A.7: The virtual environment system architecture. A.5 Redundant Actuation A.5 145 Redundant Actuation The motor torques r are computed from the desired control command u in the coordinate space; however, since there is redundancy in the actuation, r is not unique. This redundancy may be used to minimize the internal force applied along the linkage that connects the two pantograph end-points. The following kinematic relationships are known: Jel 0 X2 e 0 621 Je2 #12 Ve2 and (A.10) #22 (A.ll) where X g is the vector of pantograph endpoint velocities, q are the angular velocities of the active joint angles, and J is a block matrix containing the Manipulator Jacobians of both pantographs. e This implies the following static force transformations: T ": fxl T21 fyl J f T = e t i J fx2 T12 and (A.12) J 2_ 2 fxh u = fyh = J f T (A.13) Th where r are the torques applied at the active joints, f is a vector of forces applied by the two pantograph end-points, and u are the wrenches applied at the handle by the combined mechanism. To obtain the minimum norm solution for r in the coordinate space, the right pseudo-inverse is A.5 Redundant Actuation 146 computed in (A. 13): • f = Uo ) u f = MJ^Jo)-^. T f (A.14) Substitute (A.14) into (A. 12): r = jfJoiflJo)- !! • 1 (A.15) The minimum norm solution for T in joint space starts from (A. 12): l=J~ T. (A.16) T Substitute (A.16) into (A. 13): u = r = T = jTj~ r T (JZJ- )*U T J- J {J%J- J- J )- u l T 0 l l Q . (A.17) If the torque capabilities of the actuators are not identical, or if the actuators have different gear ratios, then a weighted inverse may be necessary. The torque weighting is expressed as follows: r = WT' 0 0 Wi 2 0 0 0 0 0 0 0 Wi2 o. 0 0 (A.18) W2 2 where r ' are the actuator torques before any gearing. The minimum norm solution in "actuator" A.5 Redundant Actuation 147 joint space is computed from (A. 17) as follows: jTj- T = U TT T-T jTjWr' =u T T r' = WJ~ J {j7j- W Jl T Q 2 1 J ) 0 _ 1 U . (A.19) Appendix B Motion Imaging and Camera Calibration Tissue phantom deformation is measured from video images taken during needle insertion experiments, as described in Chapter 3. This appendix details hardware specifications, camera calibration and marker tracking. B.1 Camera Specifications and Image Acquisition A Silicon Imaging SI-3170RGB, 3.2 Megapixel digital camera is used to image the tissue phantom for motion analysis. This is a high-definition C M O S camera that is capable of acquiring 30 frames per second at its full 2056 x f 560 resolution. For higher frame rate applications, a windowed area can 6mm, be measured at over 500 frames per second. Detailed specifications are given in Table B.1. A f/f.2 Fujinon D F 6 H A - 1 C C T V lens is attached to the SI-3170 camera. The camera system is shown in Figure 3.2. Table B.1: SI-3170 Camera Specifications (Silicon Imaging) Sensor: Active Pixels Optical Imaging Format 2056 (H) x 1544 (V) I" Pixel Size (pitch) 3.3/im x 3.3^m Pixel Type Aspect Ratio 2 C M O S Active Pixel 1:1 continued on next page 148 149 B.1 Camera Specifications and Image Acquisition continued from previous page Spectral Response 350 ~ HOOnm Dynamic Range 66dB (Vsat/Read Noise) Fill Factor 38% Fill Factor with MicroLens 60% Sensitivity 0.5 Lux at F1.0 Linearity (5-70%) ±2.5% SAT QE at 540nm 0.58e~/photon Read Noise 20eDark Current Noise at 295K < 3e~ Saturation Capacity 35,000e" Conversion Gain 36.0^ Vsat 2.7Ve Shutter Rolling Shutter Shutter Speed Variable, 4 to 4091 line times Long Integration n-Frame Time Readout Progressive Scan, Windowed Colour Filter Bayer RGB A/D Conversion and Digital Clock Synthesizer: A/D Conversion 2 channels at 50MHz (Nominal) Vertical Resolution 12-bit Pixel Clock Frequency 40 ~ 100MHz Programmable Adjustment Method Serial command protocol A/D SNR 67.5dB Output Noise 0.2 LSB rms Digital Video Output, 12 Bit Multiplex LVDS (CameraLink): Readout Rate 100MHz @ 12-bit Readout Format 12-bit Dual Channel (Ports A, B, C) Frame Rate 2056 x 1544 at 30fps 1600 x 1200 at 38fps UXGA 1920 x 1080 at 48fps (16:9) 1280 x 1024 at 48fps SXGA 1280 x 720 at 60fps (16:9) 1024 x 768 at 60fps XGA 640 x 480 at lOOfps VGA 640 x 240 at 200fps 256 x 128 at 360fps 24 x 6 at 2500fps Signal-to-Noise > 60dB (fc=20MHz, Gains=1.0) continued on next page B.1 Camera Specifications and Image Acquisition 150 continued from previous page Connector M D R 26-pin connector (3M 10226-6212VC) CameraLink Control and Communication: Serial Communication RS-232 Protocol 9.6Kbps Signaling T X and R X (LVDS) G P Asynch Triggers 3 (LVDS) General Purpose I / O 4 T T L (on camera P C B header) High Speed Shutter 40nsec ~ 2msec Long Integration 1 sec ~ 5 min Region-of Interest 24 x 6 to 2056 x 1544 in 24 x 6 steps R , G , B Independent Gains 4 Settings each ( l x , 1.2x, 1.5x, 2x) Overall Gain 3 Settings ( l x , 2 x , 4 x ) Setting Timing Next top of frame Data, Power, Trigger, RS-232 MDR-26 (optional + 5 V D C Power jack) Serial Control Mode Full/window, Live/Single, Shutter Speed, Frame rate, Long Integration, Clock Speed and Gains. Power: Input Voltage +5VDC ±10% Power 2.5 Watts Power Connection CameraLink Connector Using C C 2,3,4 (3 pairs) Mechanical: Lens Mount C-Mount Enclosure Size 1.8" W x 2.0" H x 1.5" L Weight 12 oz. Camera Mount | " x 20 standard tripod mount The SI-3170 signal interface conforms to the digital Camera Link standard, and is connected to an E P I X PIXCI® C L 1 Camera Link Interface. The PIXCI® C L 1 operates in a 32-bit, 33MHz P C I bus slot, and transfers video data to the P C I bus at a sustained rate of up to 100Mbps. For full frames of 2056 x 1560 x 8-bit pixels, this permits a frame rate of 30fps. A pixel depth of 12bits would require almost 200Mbps data transfer, which is infeasible and unnecessary in this work. During experiments, images are transferred to memory and subsequently stored to disk for offline processing that includes camera calibration and motion tracking. B.2 Camera B.2 Calibration 151 Camera Calibration A video camera, with its associated optics, maps three dimensional scene coordinates to two dimensional image coordinates. This mapping differs from camera to camera and the purpose of camera calibration is to describe the mapping for a particular camera system. Parametric camera models have been studied extensively (see e.g., [172]). Physical camera parameters include both intrinsic and extrinsic parameters, where intrinsic parameters typically include the effective focal length, aspect ratio, scale factor, skew, principal point (centre of the image) and lens distortion parameters. The extrinsic parameters describe the transformation from object coordinates to a camera-centred coordinate frame. This section describes a camera model, as well as an approach for determining its physical intrinsic parameters experimentally. B.2.1 Camera Model Consider a pinhole camera model, as shown in Figure B.1. Each point in a scene projects along a straight line through the projection centre onto the image plane. The origin of the camera coordinate Figure B.1: The pinhole model projects scene points along a straight line through a projection centre (the pin hole), onto the image plane. frame lies at the projection centre, and the z-axis is perpendicular to the image plane. Expressed in the camera frame, a point P has coordinates (x , y , z ), and projects onto the image plane at s s s s B.2 Camera 152 Calibration point P with image coordinates (x , y ): n n n = !_ (B.1) Vn where / is the focal length. For a focal length / = 1, P = [x y ] T n s s jz . s This simple approximation does not account for optical aberrations that are observed for physical lens systems, as shown in Figure B.2. Most C C T V lenses exhibit some distortion in the image Figure B.2: Image projection through a lens system. coordinates, which is typically expressed as a combination of radial and decentering distortions as follows: d +d r = c (1 + hr + 2 kr+ A 2 kr 6 3 + ...)P + 2p x y 1 n + p (r + 2x1) 2 n 2 n (r 2 Pl (B.2) + 2y ) + 2p x 2 n 2 nUn where P& is the image point coordinate after being displaced by lens distortions and r is the distance between P and the principal point, corresponding to a radius ( r = x + y ) [143]. Radial distortion 2 2 2 n causes points to be displaced radially in the image plane and is expressed as a power series in the radius, with coefficients Ari, k and k% etc. Decentering distortions arise when the lens centres are 2 B.2 Camera Calibration 153 not strictly collinear, giving rise to further radial and tangential distortion. The distortion vector d is parameterised by p\ and c Finally, the image is mapped to pixel coordinates by the image sensor: P = fP p d +P , (B.3) c where / is the focal length in pixels (/ includes the scale factor), P is the location of the principal c point on the image plane (in pixel coordinates), and P is expressed as pixel coordinates (u , p p v ). p Therefore, a point P in the scene is measured at P in the camera image. The SI-3170 has square s p pixel cells; therefore, the aspect ratio is 1:1 and the skew factor is zero. Sample images are shown in Figure B.3. The effect of radial distortion is clearly evident in the images—straight lines within the scene become curved in the images. B.2.2 E s t i m a t i n g the Intrinsic C a m e r a Parameters The intrinsic camera parameters are estimated by observing images of known 3-D geometry in the scene. Figure B.3 shows a planar checkerboard pattern that is imaged at several different positions and orientations in the scene. A n object coordinate frame is attached to the checkerboard so that there exists a transformation between object coordinates and camera coordinates in each image. Each transformation is described by a rotation vector ctj and a displacement vector dj; extrinsic parameters that are not known. Both the intrinsic and extrinsic parameters are estimated simultaneously during camera calibration in order to minimise the reprojection error between point features in the object space and point features in the image space. Corner points in the checkerboard pattern are chosen as candidate feature points, since they are easily and accurately detected in the image space, as shown in Figure B.4. Because the dimensions of the checkerboard pattern are known precisely, the feature point locations are also known in the object space. Projection from object space to image space using the camera model and estimated intrinsic and extrinsic parameters is compared to the feature points detected in the images. The reprojection error is minimised, in the B.2 Camera Calibration 154 Figure B.3: Sample images ofthe planar calibration pattern, acquired using the SI-3170RGB camera. least squares sense, over the intrinsic and extrinsic parameters $ j and $ : e N mm J2 [(«« - <W) 2 + & - ^)) } v ' 2 where, $ = = $ , U $ C {/, P , fci, k , fa, p i , p } U { a i , a , c 2 2 2 a, n d u d, 2 d} , n (B.5) B.2 Camera Calibration 155 Extracted corners Xc (in camera frame) Figure B.4: Corner feature points on the checkerboard pattern are located in the calibration images, and used to evaluate the reprojection error. [Hi Vi] is the i th the i th image feature point extracted from the image, w^($)] is the location of image feature point derived by reprojecting using intrinsic and extrinsic parameters N is the total number of image feature points (corners) and n is the number of calibration images. Algorithms for implementing this procedure are available in the Camera Calibration Toolbox for M A T L A B ® [173]. Given a sufficiently large set of images, each showing the checkerboard pattern in a distinct pose, the minimisation in (B.4) is shown to converge using a gradient-descent algorithm [173]. The parameters $ that minimise (B.4) for a set of twenty-one images, including those in Figure B . 3 , are given in Table B . 2 . In this example, the standard deviation ofthe reprojection error along each image coordinate axis was approximately 0.4 pixels. B.2 Camera Calibration 156 Table B.2: Intrinsic parameters estimated from camera calibration images. Principal Point [1044, 737] Focal Length 1884 pixels Radial Distortion Coefficients [k\, k , £3] [-0.268, 0.183, 0.000] Tangential Distortion Coefficients [pi, p ] [-0.00096, -0.00650] 2 2 B.2.3 Image Rectification The radial and tangential distortion fields described by these parameters are illustrated in Figure B.5. Radial distortions are clearly significant for this lens, resulting in displacements of up to Radial Component of the Distortion Model 0 200 400 600 BOO 1000 1200 1400 Tangential Component of the Distortion Model 1G0O 1BO0 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Figure B.5: Radial and tangential distortion fields. Image coordinate displacements due to radial distortions are significantly larger than those due to tangential distortions. Contour annotations are measured in pixels. 80 pixels. For image-based motion tracking it is desirable to compensate for this displacement by rectifying image features using an inverse of the distortion model (B.2). Therefore, given an image pixel coordinate P , it is necessary to estimate fP p n + P , the image coordinate that would have c B.2 Camera Calibration 157 been observed in the absence of lens distortion. This is achieved by estimating P . From (B.3), n Pd = P (Pn)o = Pd, p j ° P , and (B.6) where (P )o is an initial estimate for P , initially set to P^. n n r = \J{Pn)J{Pn)i (B.7) Equation (B.8) is iterated so that (PnWi tends toward P . n rectified image pixel coordinate is P = f{P )M r n After M iterations the estimated + P - In practice, approximately 20 iterations are c required for convergence to within one pixel. Figure B.6 shows an image both before and after rectification. Figure B.6: A sample image shown before and after rectification. Straight lines and edges are restored in the corrected image. Appendix C Linear Elasticity and the Finite Element Method This appendix details the modelling of linear elastic bodies from constitutive equations through discretisation by means of the Finite Element Method. The soft plastic tissue phantoms and steel needles that are used in the experimental system (described in Chapter 3) are modelled and simulated using the methods described here. C.l State of Stress An elastic body under external load experiences internal stresses, as shown in Figure C.l. A small area element AA, shown oriented with unit normal n, experiences a resultant force Af and moment Am. At the limit, where AA is an infinitesimal area, Af h m h m A A A df - ° AA = d A ^ ° AA = d A = ' t = C - a n d ( C - 1 } Thefirst-ordertensor t is called a traction or stress vector, while c is called the couple-stress vector. From the elementary theory of elasticity c = 0, and the stress vector is of primary interest [174]. t is the stress vector at a point, with respect to a particular area element specified by the unit normal n; whereas, a complete description of stress at a point requires the state of stress in all directions. 158 C.1 State of Stress 159 Figure C.1: Force and moment vectors at a point, with respect to an area element, within a deformable body under load. Consider an infinitesimal rectangular parallelepiped located at a point inside the elastic body. A coordinate frame is chosen with orthonormal basis vectors {x, y, z} parallel to the faces of the parallelepiped, as shown in Figure C.2. Each face corresponds to an area element similar to AA Figure C.2: Components of stress. above, and the stress vectors at each face are shown as t , tj,, t . Each stress vector can be written x 2 in terms of its Cartesian components, e.g., t = o x + o y + <r z, where the coefficients x xx xy xz are called the the components of stress. The second-order tensor _> containing all of these components C.2 Strain and 160 Deformation defines the state of stress at a point: ®xy °~xz &yx °yy 0~yz °~zx 0~zy 0~zz (C.2) where the stress vector t can be computed with respect to an infinitesimal area element with unit normal vector n , i.e., t = \I>n. Under the assumption that body-moment and couple-stress do not exist, it can be shown that _> is symmetric, thus CTJJ = CTJJ [175]. The state of stress at a point inside an elastic body is therefore described by six components of stress: (7 T where o, u xx yy and o zz are normal &yz yy stress ®z components, and (C.3) >xy o, yz a z x and a x y are the shear stress components. C.2 Strain and Deformation The deformation of a material body under the action of stress is characterised by strain. The second-order strain tensor that describes the state of strain at a point is defined as: €xx ^xy €xz T = (CA) tyx (-yy £yz €zx £zy €zz Similarly to Vl/, it can be shown that T is symmetric for elastic bodies, so that e-jj = tji [175]. The resulting strain tensor is written as: ^ tyy &zz £yz ^zx ^xy (C.5) A point P in a non-rigid body is illustrated in Figure C.3, both before and after the body is deformed by external loads. It experiences a displacement u due to the deformation, and is centred at the C.2 Strain and Deformation 161 X Figure C.3: A non-rigid body before and after deformation. reference coordinate frame Co m the undeformed state, written as ° P . fn the deformed body state, the point is expressed as P, and is centred at the new coordinate frame C_ . X x Point-wise displacements in the deforming body are written as u = [it(x) v(x) w ( x ) ] , where T x = (x, y, z) is a point in the reference body frame C . The engineering strain tensor is expressed 0 C.2 Strain and Deformation 162 in terms of u [174,175], and becomes: du I dx ^ dv dy i + 1 2 I 2 (t) +(g) +(f) ] 2 2 2 yy c dw I 1 dz ~t~ 2 ; e(u) lyz Izx 9«) dv I du I f2l£ 9u i 92 ( i ) + 2 + ( i ) + (^) ] 2 2 du du i dv dv dy dz ' dy dz + du du dv dv dz dx + dv dv dx dy + dw dx 9a; + [ dz dx cH> + i dz dx. Jxy dy dx du du dx dy d_ dx 0 o 4z + d_ dy 1 2 0 0 0 T 0 e Ou T z o z 0 dx d_ dx 0 °z w d_ dy dw dy 0 0 dy o i ^y el 0 (C.6) where 0 r = du dx dv dx dw dx and 6 and 6 are defined similarly [147]. This is called a Lagrangian y Z formulation, as opposed to a Eulerian formulation, because e(u) is computed with respect to the reference coordinate frame C , rather than the frame C i in the deformed body. The vector e(u) is 0 called the Green-Lagrange strain tensor. C.3 Constitutive Equations 163 If u is small, then the quadratic terms in (C.6) can be neglected: du dx d_ dx 0 0 dv dy 0 JL dy 0 dw dz 0 0 d dz lyz dv _|_ djv dz dy 0 d dz d_ dy Izx du dz dw_ dx d dz 0 d_ dx dv_ dx d_ dy d_ dx 0 _ - €xx yy € e(u) = lxy_ du dy i i = Bu. (C.7) This is called Cauchy's infinitesimal strain tensor, and is the basis of the small-strain analysis of elastic materials. The Cauchy strain tensor is not rotation invariant [147], and must be used with care. In this thesis, the strains experienced by the tissue phantoms and tissue models are assumed to be small; therefore, the linear strain tensor is used in order to obtain a linear system of equations. This is shown later. The needle deflection model is based on the Green-Lagrange strain tensor that contains quadratic terms that become significant due to large strains. C.3 Constitutive Equations A material that is loaded to a level below the yield strain, and that unloads along the same path is U called elastic. If the path is also linear, then the material is said to be l i n e a r l y elastic" [174]. The generalised Hooke's law states that: a = Ce , (C.8) where C is a material matrix that characterises the Hookean elastic material. For isotropic materials, C requires two parameters, namely Young's modulus E and Poisson's ratio u, or the Lame constants CA The Discrete Solution 164 u. and A: A = A* = Ev (l + u)(l-2u) E 2(1 + 22/) Young's modulus relates normal stress and strain a xx = Ee , xx (C.9) while Poisson's ratio v is the ratio of lateral contraction to longitudinal strain (it is also related to volume change or dilation). The material matrix C is defined as: A + 2/x A A 0 0 0 A A + 2,u A 0 0 0 A A 0 0 0 • fi 0 0 0 0 0 0 n 0 0 0 0 0 0 fi C = A+ 2 / x O O O (CIO) and is positive definite; therefore, (C.8) is an elliptic partial differential equation that can be used to relate forces and deformations applied to a body. C.4 The Discrete Solution Solutions to boundary-value problems involving (C.8) may be computed analytically only for simple cases. Complex body geometries can be divided into simple shapes that can each be handled analytically. This is the basis of Finite Element Analysis and the Finite Element Method. The geometry of a continuous body is approximated by a discrete representation. In order to solve (C.8), it is convenient to use an energy formulation [134], where the strain energy E strain of an elastic body is written as: E train = \ s 1 f e ( x ) cr(x) dx . r Jn (Cll) CA The Discrete Solution 165 After substituting the material matrix: Estrazn = \ I e (x)Ce(x) dx , (C.12) T where cr and e are stress and strain vectors respectively, and the strain energy is integrated over the body domain f_. The total energy of the body under load can be rewritten in terms of body displacements u(x) at each point x £ f_, as follows: E(u) = I [{BufCBudx: - J /f (x)u«ix, (C.13) T where the first term is strain energy and the second term is work done by external forces applied to the body (surface or internal) [134,176]. For Finite Element Analysis, the displacement field over the body domain Q, is discretised into a finite number of elements Q , as illustrated in Figure C.4, in two dimensions. The elements are e Figure C.4: The domain __ of the elastic body is discretised into elements fie- assembled in a mesh-like lattice of nodes/vertices. The displacement at each node, due to body deformation, is written Uj = [u{ vi Wi] , and the compound vector of all nodal displacements over T the body is written u = [uj the e th ... u^] , where there are n nodes in the mesh. The vertices for T element lie at the points Pf, where i = 1,2,..., n and n is the number of vertices that e e make up the element. Examples of triangular and tetrahedral elements are shown in Figure C.5. The compound vector of element node displacements is written as u = [(uf ) e T (uf ) ... (u^ ) ] , r r e T CA The Discrete Solution 166 (a) (b) Figure C.5: Triangular and tetrahedral elements. where uf are the node displacement vectors at the vertices of element e. The displacement at a point x lying in fl , can be expressed in terms of element node displacements: e «(x) £ iV7(x)u: = i=l = N u e (C.14) e where A^(x) are called shape, basis, or weighting functions. For linear shape functions, iV?(x) are the natural or barycentric coordinates ofthe element [177], which have the following form: Nf (x) = of + b\x + dy + d e iZ i = 1,2,...,n e (C.15) The solution to the deformation problem lies at the equilibrium at which the total energy of the elastic system is minimised, where the first variation of the energy functional vanishes: SE(u) = J2 6E<i( u) = ° (C.16) C.4 The Discrete Solution 167 The energy differential for a single element is: SE {u) = VE -8u e e e 6E{u) = J2( VEe- S^ e) = 0 e • ( C- 17) The nodal displacement vector u contains the independent variables, therefore \7E e e = 0, which implies that: / {BN f e C{B~N )u e [ d x - f = 0 e B CB u eT e u If e B CB u eT e dx = f e e dx == i e 0.' /B e T C5 u e Ku c e where f K e e e = f = i e (C.18) is a compound vector containing the applied force components at each element vertex, is the element stiffness matrix. Equation (C.18) may be computed for each element in the discretised body, and combined to form the entire set of system equations: K u = f, (C.19) where K is the global stiffness matrix, u is the compound node displacement vector and f is the compound node force vector, that are assembled as follows: K = ^global(i_" ) e e u = global (u ) e e f = ^global(f) , (C.20) e where the function "globalQ" maps element node numbers to global node numbers [134]. The global stiffness matrix K is singular as it stands in (C.19), since there are no geometrical C.5 Large Strain Analysis 168 constraints on the body mesh. Consider a single displacement constraint on the i th degree of freedom in Equation (C.19): 11 •• • fcii • 21 ' ' • k • • 2i kin h k h 2n U 'il •• = kii kin fi nl • • kni knn fn kin fl k u .. o • kl 2 .. o • • 0 • • 1 •• 0 kni .. o • diku h - k 2n dik 2i u = (C.21) di fn knn where di is the displacement constraint applied to the i degree of freedom. In Euclidean-space, at th least six independent displacement constraints are required in order to uniquely define the location of the body in space, and hence ensure that K is invertible. C.5 Large Strain Analysis Analysis using the quadratic strain terms proceeds in a similar fashion with the Green-Lagrange strain tensor defined in Equation (C.6) as follows: e(u) =5u+ |A(u)6(u) , (C.22) C.5 Large Strain Analysis 169 where 0 ( u ) is defined as: ° T dx 0(u) = A (3x3) u = Tu C? u e (C.23) e 9 T djj (3x3) A The strain energy is expressed as d follows: T Wz (3x3) A E.strain — . - /f eTc dx 1 \ I (Bu + ]-AT\i) C{B\i 2 JQ 2 + ]:ATu) dy 2 T 2 ii J / T + -u T A CBu 2 T T T + -u B CATu 2 T T u £ ? C B u + \U B CATVL T 2 Jnl T T T J 2 + T -u T A CATu d x T 4 T T (C.24) dx , = 0. B y substituting (C.23) into (C.24), the element strain energy is expressed as: where T A T u B CBu T zpe (N u ) B CB N u e strain \2 f e T T e e u B CB u Jn L eT eT e + -(NV) B CiTN u e T + e e T \u B CAG u e eT 2 eT e e e dx (C.25) dx . From (C.18): 5E e = f B CB u dx eT e + \ f e v B CB u e eT e e B CAG u dx eT e + \v eT e + e Ku e e f e B CA{u )G u e Ku where K - e e e e e K (u )u e e e q 0 f f (C.26) + f^(u ) e is the element stiffness matrix resulting from the linear strain terms, and Kg(u ) e is the element stiffness matrix resulting from the quadratic strain terms, f^ is is a compound force vector containing "correction" forces due to quadratic strain terms. The global linear stiffness matrix, compound displacement and force vectors are assembled according to Equation (C.20). C. 6 Linear Elastic Models in Two C.6 Dimensions 170 Linear Elastic Models in Two Dimensions This section contains details for computing Equations (C.19) and (C.26) in two dimensions, using the barycentric coordinates of triangular elements, as shown in Figure C.5(a). The elasticity problem can be reduced from three to two dimensions if there is no traction in one dimension through the body; therefore, if the is stress-free then o z-axis xz = o vz = o zz = 0. In addition, it is assumed that all remaining stresses are functions of x and y only. Elementary plane stress theory also neglects e , so that both the stress and strain tensors are reduced to 3-tuples 22 [174]. In two dimensions x = [x y] , u = [u v] ', and the displacement vector at the i node on T T th element e is uf = [uf vf] . T C.6.1 Linear Strain In order to form the elementary stiffness matrices K in the linear strain case, it is necessary to e compute C, v and B in Equation (C.18). The (6 x 6) material matrix in (C.10) is reduced to the e e following (3 x 3) matrix in two dimensions: C = X + 2p A 0 A A + 2/i 0 0 0 ix (C.27) The scalar value v comes from the integral [J dx] in Equation (C.f8), which is the area of element e Qe e. If the plain stress model approximates a body having one dimension much smaller than the other two dimensions, such as a thin sheet or the tissue phantom shown in Chapter 3, then this area must be scaled by the thickness of the body, resulting in the element volume v . e The matrix B is the product of B and N (see Equation C.f8), where B is defined in (C.6) e e and (C.7). N is defined as follows: e N e = N l 1(2x2) 2 N 1(2x2) $ 1(2x2) N (C.28) where Nf are the barycentric coordinates of element e, which are functions of x and y. Therefore C.6 Linear Elastic Models in Two Dimensions 171 u(x) = N u . The barycentric coordinates are derived as follows: e e 1 i i i iVf(x) X xi x xz ATf(x) 2/1 2/2 2/3_ 2 (C.29) _/V|(x)_ where (xi, j/i), (x , y ) and (0:3, j/3) are the element node coordinates in the reference frame C . 2 2 0 This implies the following: -1 A/f(x) 7V|(x) _/V|(x)_ = 1 1 1 1 »i ^2 £3 x 2/1 2/2 2/3 af — e n cf 6f ue 2 "2 «3 *>3 u e X 2 c c af + 1 = y. 3_ b\ x + c\y (C.30) a | + b\x + c\y a% + b\x + cfy Finally, 0 b\ e(x) = B u ( x ) = 5 N ( x ) u e 0 bi e S U c C.6.2 6§ «e l Le °l c e 2 °2 e r c 3 e u (C.31) e ue °3 Quadratic Strain For large strain analysis using the Green-Lagrange strain tensor, two additional matrices A(u ) and e G are required in order to compute K (u ), e e q and hence f^(u ) in (C.26). A(u ) e du dx A(u) = dv dx 0 dv dy du 0 ^ dv ^ dy dy du dv dx dx 0 du dy e follows from (C.6): 0 (C.32) 172 C. 7 Beyond Elastostatic Models where: du dx + 63 u b\ui + b\ui b\vi + b%V2 + b\v3 c\u\ + C%U2 + C3U3 3 dv dx du dy dv (C.33) C-\vi + C%V2 + C3V3 . dy G follows from (C.23), and is defined as follows: e d_ dx G = TN e e 0 a_ 0 = a_ 0 0 d_ dy. dy w 0 0 0 bl 0 °2 0 °3 0 0 K 0 °2 0 °3 l 0 2 0 3 0 c 0 C.7 Nf dx c l c 0 e r 2 c c 0 C NI (C.34) 3. Beyond Elastostatic Models This section provides a brief outline of dynamic, plastic and fluid models that extend the linear elastostatic formulations detailed in this appendix. C.7.1 Dynamic Elastic Materials Dynamic behaviour may be incorporated into the linear elastic model by means of the Formulation, Newtonian which has the following form: Mu + Bu + Ku = f , where M and B are mass and damping (C.35) matrices respectively. Although Equation (C.35) is shown to have mass, damping and stiffness matrices, it is important to note that this formulation is not C. 7 Beyond Elastostatic Models 173 equivalent to a mass-spring-damper model. A continuum model is discretised using the Finite Element Method, and as such, no spring lengths are computed. Therefore, numerical simulations of these models will remain stable if semi-implicit integration techniques are used. For a discussion of the form of M and B, as well as issues related to mass-lumping, see [147]. Equation (C.35) may be integrated using both explicit and implicit numerical methods, with a trade-off between speed and numerical stability [134]. C.7.2 Plasticity Many materials exhibit plastic behaviour once their yield strength is exceeded [178], resulting in permanent inelastic deformation. This phenomenon can be modelled by displacing the node positions in the reference mesh when local strains exceed a material dependent plastic or yield limit, as shown by Terzopoulos et al. [55]. A similar approach has been used by Terzopoulos et al. to simulate visco-elastic materials. C.7.3 Stokes F l u i d s The behaviour of viscous fluids is distinctly different from that of elastic materials. This section briefly outlines Stoke's Fluid models. Stoke's equation for incompressible viscous fluids [125,177] is derived from the Navier-Stokes equation by neglecting thefluidinertia term pv • Vv. This is a valid simplification forfluidswith a low Reynolds Number and yields Stoke's equation: - Vp + pV v + f = 0 2 V-v = 0 (C.36) where p is unknownfluidpressure field, \x a viscosity parameter, v the unknown velocity field and f the force field applied to the fluid body. Incompressibility is achieved by imposing a zero divergence constraint on the velocity field. A weak formulation is employed: - / Vp-wd__ + / pV v-wd__ + / fwd__ 2 Jn Jn Jn I qV-vdtt Jn = 0 = 0, (C.37) Models 174 test functions with piecewise-continuous derivatives over the fluid body domain C. 7 Beyond Elastostatic where w and q are Q. In R,, Equation (C.37) is discretised by means ofthe 2 Galerkin method [179]: TV P TV = i=i M ^PiQiix) i=l (C.38) ,and -/ / Jn TV .Jn • Vw = / /» Wj dfi , / -E« Jn 0, ^ ) fix^ '("«'i') ydy 1 V J i = 1,2,N, 1 1,2, ...,M (C.39) This can be expressed in matrix form as follows: C.7.4 V f p 0 (C.40) Composite Models This section demonstrates the implementation of a composite model containing both linear elastic and fluid elements. A rectangular fluid region trapped within a linear elastic body was simulated using Equations (C.19) and (C.39). Different fluid velocity and pressure grids were selected with a linear-linear element pair in order to satisfy the div-stability criterion [180]. Both the velocity and pressure domains are triangulated; however, the velocity field is given afinertriangulation resulting in four triangular velocity elements per pressure element. Therefore w, and qi are simply the barycentric coordinates of the triangular velocity and pressure elements. The elastic and fluid bodies are discretised as described C. 7 Beyond Elastostatic Models 175 in Section C.7.3, yielding systems of the form: Elastic Body u ==> 1 (C.41) = £+ Kn Fluid Body => Arr Ar2 Aiv Ar A22 Lv2 Aty LIT Lr 2 2 2 2 0 = £+ g (C.42) 2 P 2 where T denotes the shared boundary between elastic and fluid domains. Note that the Stoke's formulation has been converted from an Eulerian formulation to a Lagrangian formulation by scaling the velocity field by the discrete time step interval A i . This transformation assumes that deformations between time steps are infinitesimal, but allows for the coupling of boundary conditions. Therefore, for continuity and equilibrium at the fluid-elastic body interface, u = A i v . = u . The 2 r r composite system becomes: K Kn ir 0 0 u K-pr + Apr A2 Ln Mr F 0 A A22 Lr2 0 Lir Lr 2r 2 0 1 Aiv P 2 f + g 1 f 2 f + g (C.43) 2 0 Given a suitable set of constraining boundary conditions, as well as the global force field f and external forces g and g acting upon the elastic and fluid bodies, the deformation vector [u 1 2 1 u Aiv ] 2 r r can be computed, as shown in Figure (C.6). The incompressible nature of the fluid body is evident in the simulation result. C. 7 Beyond Elastostatic Models Figure C.6: A composite fluid-elastic body under compressive and tensile loads. The grey areas are regions of fluid, while the remainder of the body is simulated using a linear elastic model. 176
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Modelling, simulation and planning of needle motion...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Modelling, simulation and planning of needle motion in soft tissues Dimaio, Simon P. 2003
pdf
Page Metadata
Item Metadata
Title | Modelling, simulation and planning of needle motion in soft tissues |
Creator |
Dimaio, Simon P. |
Date Issued | 2003 |
Description | Precise needle placement is required for the success of a wide variety of percutaneous interventions in medicine. Insertions into soft tissues can be difficult to learn and to perform, due to tissue deformation, needle deflection and limited visual feedback. Little quantitative information is known about the interaction between needles and soft tissues during puncture, and no effective physically-based training, planning and guidance systems exist for such procedures. This work aims to characterise needle-tissue interaction by measuring contact forces and deformations that are applied during insertions into soft tissue phantoms. A new methodology for estimating the forces that occur along the needle shaft during insertion is described. The approach is based on physical experiments, as well as on linear elastic phantom models that are discretised by traditional Finite Element Methods. Shaft force distributions are derived from insertions into homogeneous and simple layered inhomogeneous tissue phantoms at several driving velocities, and are applied as boundary conditions to tissue models for physically-based simulations of needle insertion trajectories. A large-strain elastic needle model is coupled to the tissue models to account for needle deflection and bending during simulated insertion. Since the force-displacement relationship is only of interest along the needle shaft, a condensation technique is shown to reduce the computational complexity of linear simulation models significantly. The boundary conditions that determine the tissue and needle motion change as the needle penetrates, or is withdrawn from the tissue model. Boundary condition and local material coordinate changes are facilitated by fast low-rank matrix updates. Such numerical schemes have been seen in prior work involving point and surface interaction; however, in this work the condensation state, boundary conditions and material coordinates evolve as the needle penetrates the tissue volume, and as internal contact states change. These novel interactive simulation techniques allow users to manipulate a three-degree-of-freedom virtual needle as it penetrates virtual tissue models, while experiencing steering torques and forces through a planar haptic interface. Models and simulations are also used to formulate needle insertion as a trajectory planning and control problem. The concept of needle steering is developed, and a Needle Manipulation Jacobian is defined to express the relationship between the needle base and tip velocities. This concept is used in conjunction with a potential-field-based path planning technique to demonstrate needle tip placement and obstacle avoidance. Results from open loop insertion experiments are also provided. |
Extent | 31554356 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065587 |
URI | http://hdl.handle.net/2429/15012 |
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 |
Graduation Date | 2003-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2003-860108.pdf [ 30.09MB ]
- Metadata
- JSON: 831-1.0065587.json
- JSON-LD: 831-1.0065587-ld.json
- RDF/XML (Pretty): 831-1.0065587-rdf.xml
- RDF/JSON: 831-1.0065587-rdf.json
- Turtle: 831-1.0065587-turtle.txt
- N-Triples: 831-1.0065587-rdf-ntriples.txt
- Original Record: 831-1.0065587-source.json
- Full Text
- 831-1.0065587-fulltext.txt
- Citation
- 831-1.0065587.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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065587/manifest