Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Simulation of linear elastic media with fluid inclusions Gosline, Andrew Havens 2003

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_2004-0031.pdf [ 8.22MB ]
JSON: 831-1.0065176.json
JSON-LD: 831-1.0065176-ld.json
RDF/XML (Pretty): 831-1.0065176-rdf.xml
RDF/JSON: 831-1.0065176-rdf.json
Turtle: 831-1.0065176-turtle.txt
N-Triples: 831-1.0065176-rdf-ntriples.txt
Original Record: 831-1.0065176-source.json
Full Text

Full Text

Simulation of Linear Elastic Media with Fluid Inclusions by Andrew Havens Gosline B . S c , Queen's University, Kingston, Ontario, Canada, 2001 A THESIS S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R OF A P P L I E D S C I E N C E in T H E F A C U L T Y O F G R A D U A T E STUDIES (Department of Electrical and Computer Engineering) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A December 2003 © Andrew Havens Gosline, 2003 Library Authorization In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Name of Author (please print) Date (dd/mm/yyyy) Title of Thesis: ^ liMpi Ctjiort 0^ L/hea r (~'/ps'fyt Degree: /[AA^C Y e a r : lOO^ Department of Q1-CctrtQoA -i CouU^Ljicr (n^ioJ^C,^ The University of British Columbia & Vancouver, B C Canada Abstract This work presents a methodology for haptic rendering of fluid filled structures that are enclosed in linear elastic media using the Finite Element Method. Haptic medical simulation is a growing field of research motivated by creating risk-free virtual envi-ronments for medical students to learn and practise surgical procedures. One of the challenges in creat-ing medical simulators is modeling the deformation of living tissues. Due to minimum haptic update rate requirements, the deformable methods are simplified and precomputed. Most medical simulators model anatomy as elastic material with constant or varying stiffness. However, human anatomy includes a variety of fluid filled structures. To improve the realism of these simulators, fluid filled bodies should be modeled in addition to elastic media. This thesis presents a method for simulating fluid effects by adding hydrostatic fluid pressure to a body of elastic material modeled with the Finite Element Method. By distributing fluid pressure across an interior cavity surface, the fluid can be modeled using a force boundary condition. Proportional feedback is used to solve for an incompressible fluid relationship between the cavity pressure and volume. Linear finite elements are used so the stiffness matrix can be condensed to achieve real-time haptic rates. To validate that this method predicts deformation of a fluid filled cavity in a realistic manner, the defor-mation a fluid filled phantom is tracked and compared to a Finite Element simulation of the same phantom. The data is found to agree well with the simulation. A real-time haptic simulation of elastic media enclosing incompressible fluid, based on an existing 2D needle insertion simulation, is presented. Numerical tests show that simulation of a relatively large 3D fluid body will be possible at haptic rates. i i Contents A b s t r a c t i i Table of Contents i i i L i s t of Tables v i i L i s t of Figures v i i i G l o s s a r y x A c k n o w l e d g e m e n t s x i i 1 I n t r o d u c t i o n 1 1.1 Background 1 1.1.1 Medical Training Simulation 2 1.1.2 Deformable Models for Tissue Simulation 4 1.1.3 Fluid Simulation 5 1.2 Motivation 5 1.3 Research Objectives 6 1.4 Thesis Outline 6 2 P r o b l e m F o r m u l a t i o n 8 2.1 Deformable Methods in Haptics 8 2.1.1 Finite Element Method in Haptics 10 2.1.2 Fluid Modeling in Haptics Applications 12 iii CONTENTS iv 2.1.3 Justification for choice of Finite Element Method 13 2.2 Fluid-Structure Interaction with the Finite Element Method 13 2.2.1 Eulerian-Lagrangian Frame 14 2.2.2 Lagrangian Frame . . . 15 2.2.3 On Incompressible Finite Element Analysis . ." 16 2.3 Problem Statement 17 2.4 Proposed Solution 18 2.4.1 Drawbacks and Limitations 18 2.5 Summary 19 3 Methodology 20 3.1 Finite Element Method for Linear Elasticity 20 3.1.1 Volume Discretization 22 3.1.2 Constraints and Boundary Conditions 23 3.2 Fluid Model 23 3:2.1 Hydrostatic Fluid Pressure as a FEM Boundary Condition 25 3.2.2 Relationships between Fluid Pressure, Volume, and Flow 27 3.2.3 Confined Incompressible Fluid 28 3.3 Numerical Method . 29 3.3.1 Proportional Controller 29 3.3.2 System Identification and Controller Tuning 30 3.3.3 Performance and Robustness of Controller 31 3.3.4 Volume Calculation 32 3.4 Summary • 34 4 Experimental Validation 35 4.1 Phantom 35 4.1.1 Ultrasound Imaging Requirements 36 4.1.2 Phantom Material Properties 37 4.2 Surface Tracking with Digitizing Pen 38 4.2.1 Calibration of Pen 39 CONTENTS v 4.3 Experimental Procedure 40 4.4 Finite Element Analysis 41 4.4.1 Mesh Generation 41 4.4.2 Boundary Condition Changes and Constraints 43 4.5 Results 44 4.5.1 Ultrasound Contour Results 44 4.5.2 Surface Tracked Results 46 4.6 Summary 48 5 Haptic Simulation 50 5.1 Haptic Algorithm 50 5.1.1 Real-time Haptic Algorithm for Incompressible Fluid 51 5.2 Interactive Haptic Simulation in 2D 52 5.2.1 3DOF Twin Pantograph 52 5.2.2 Modifications Required for Fluid Cavity Simulation 53 5.2.3 Results from 2D Interactive Haptic Simulation 53 5.3 3D Performance 55 5.3.1 Timing of an Off-line 3D Simulation 56 5.4 Summary 56 6 Conclusions and Recommendations 58 6.1 Contributions of this Thesis 58 6.1.1 Fluid Pocket Model 58 6.1.2 Validation of Fluid Structure Model using 3D Fluid Filled Phantom 59 6.1.3 Real-time Haptic Simulation of Fluid Pockets 60 6.2 Recommendations for Future Work 60 Bibliography 62 Appendices 68 A Finite Element Modeling ' 68 CONTENTS vi A. 1 Linear Elastostatic Methods 68 A. 1.1 Linear Strain Triangular Element 69 A. 1.2 Linear Strain Tetrahedral Element 71 A. 1.3 Stiffness Matrix Generation by Node Numbering 72 A. 1.4 Constraints 73 A. 1.5 Boundary Condition Changes by Low Rank Updates 74 A. 1.6 Matrix Condensation 75 A.2 3D Mesh Generation 77 A.2.1 Removal of A Node from the Connectivity List 78 A. 2.2 Procedure for Creating a 3D Tetrahedral FEM stiffness matrix 80 A. 3 Centroid of a Triangle 80 Appendices 82 B On Graphics Accelerators for Real-Time Processing 82 B. l Introduction and Previous Work 82 B.2 Graphics Accelerator Architecture 83 B.3 Larsen's Matrix Multiplaction Algorithm 84 B. 3.1 OpenGL Code Details for Matrix Multiplication Using Larsen's Method 85 B.4 Discussion of Results 87 B.4.1 Trade-offs 87 B.4.2 Conclusions 88 List of Tables 3.1 Performance of Controller 31 5.1 CPU Times for an 80 Node Fluid Structure 56 A- l Matrix Multiply Comparison Between GPU and Optimized CPU. Size = 1024 x 1024 . . . 87 vii List of Figures 1.1 Examples of Medical Training Simulators with Haptic and Graphic Feedback1 2 1.2 Data Flow in a Typical Haptic Virtual Environment 4 2.1 Truss Elements to Decouple Tangential Displacements 15 2.2 Fluid Filled Structures Enclosed in Elastic Media 17 3.1 Elastic Body Discretized with 2D Triangular Finite Elements 22 3.2 Fluid Model using Hydrostatic Pressure Normals 24 3.3 Finite Element Solution of Elastic Object with Static Fluid Pressure Load 24 3.4 Pressure Distribution over one Edge of a 2D Element 25 3.5 Center of Pressure of a Triangle 26 3.6 Fluid Relationships shown with FEM 28 3.7 Proportional Controller 30 3.8 Pressure vs. Volume of an Elastic Body 31 3.9 2D Polyhedral Area Calculation by Summing the Areas of Trapezoids 32 4.1 Ultrasound Phantom Construction 36 4.2 Ultrasound Phantom Material Properties 37 4.3 OptoTrak® Digitizing Pen 39 4.4 Complete Phantom Compression Apparatus 40 4.5 Schematic Diagram 40 4.6 Sequence of Ultrasound Images used to Re-create Phantom Geometry 42 4.7 Hand Segmented Ultrasound Image of Fluid Pocket 42 viii LIST OF FIGURES ix 4.8 Finite Element Mesh Showing Fluid Pocket Shape 43 4.9 Ultrasound vs. FEM Contours 45 4.10 Error Plot from Top Surface Tracking 47 4.11 Top Surface Motion, Tracked and FEM 48 4.12 Vertical Expansion Plot, Phantom Markers vs. FEM Simulation 49 5.1 Haptic Algorithm Data Flow 51 5.2 3DOF Twin Pantograph Haptic Interface 52 5.3 Haptic Simulation of an Incompressible Fluid inclusion in a 2D Environment 54 A. l Triangular Element 69 A.2 Tetrahedral Element 71 A.3 Tetrahedral Elements 73 A.4 Simple ID Truss System 74 A.5 3D Mesh of Parallel Slices with Tetrahedral Elements 77 A.6 Node Removal 79 A.7 Right Angle Triangle Geometry 80 A. 8 Triangle Composed of Two Right Angle Triangles 81 B. 1 A Typical Graphics Pipeline based on the OpenGL API 84 B.2 Graphical Description of Larsen's Algorithm 85 1 Glossary BEM Boundary Element Method CT Computed Tomography DOF Degree of Freedom FEM Finite Element Method FLOPS Floating Point Operations GPU Graphics Processing Unit IRED Infra Red Emitting Diode LEM Long Element Method MRI Magnetic Resnonance Imaging NURBS Nonuniform Rational B-splines SMS Spring Mass System SMVM Selective Matrix Vector Multiplication X GLOSSARY xi E Young's Modulus [Pa] V Poisson's Rado k,fi Lame Parameters a Stress Vector [Pa] e Strain Vector D Elasticity Matrix B Strain Tensor K Stiffness Matrix KP Controller Gain h Normal Vector u FEM Displacement Vector f FEM Force Vector K Updated Stiffness Matrix X Composite FEM Vector X y Composite FEM Vector Y p Fluid Density [^] p Fluid Pressure [Pa] 8 Gravitational Constant, 9.81 [ ^ Acknowledgements • First of all I thank my co-supervisors Tim Salcudean and Joseph Yan for their guidance, support, and continous feedback. I am also very grateful for the excellent motivation that they gave to encourage me to work hard towards good project goals. • I acknowledge Rob Rohling for letting me use the OptoTrak and Ultrasonix machine in my experi-ments. • I thank all of my labmates, each of whom helped me complete this work. Purang and Julian for ultra-sound and signal processing, Parvaneh for image processing, Danny and Orcun for help understanding prostate brachytherapy, Daniela for haptics, Emre for ultrasound support, Qi for graphics cards, Scott for phantoms, You Wei for OptoTrak, and Simon DiMaio for FEM. • I could not have completed this work without the support of my Mother and Father, who feed me whenever I'm in need. Also much gramatical and emotional support:) • And finally, I thank Sara for the loving support and friendship that she has given to me for the duration of my work at UBC. xii Chapter 1 Introduction 1.1 Background Haptic Simulation of medical procedures is an active area of research in biomedical engineering. Anal-ogous to flight simulators for pilots, medical simulators allow students to learn and practice delicate proce-dures that are performed on soft tissues in a risk-free virtual environment with visual and force feedback. Such procedures as suturing [15], bone dissection [3], laparoscopy [10,24], and needle insertion [34] have been developed and implemented using a variety of architectures. One of the challenges in developing medical training simulations is the real-time simulation of soft tissue deformation. Most present surgical simulations [10,15,16,25,28,34,57] assume that soft tissue is comprised of elastic material with either uniform or varying stiffness. With this assumption, anatomical structures such as fluid filled organs, blood vessels, or cysts are ignored even though their presence has a considerable impact on how the surrounding tissue deforms. By simulating fluid filled structures enclosed in elastic media, this anatomy can be included in surgical simulations, increasing the realism of the simulation. Also, procedures that specifically involve either interaction with, or avoidance of, fluid filled structures could be simulated if fluid bodies enclosed in the tissue are included. Before describing the deformable methods, a brief introduction to haptic simulation of medical proce-dures is given to acquaint the reader with the general area of research. 1 1.1 Background 2 1.1.1 Medical Training Simulation The goal of medical training simulation is to provide a realistic virtual environment within which stu-dents can practice medical procedures. As almost all medical procedures involve hand held instruments like needles, scalpels, or laparoscopic tools, these training simulator systems must allow students to inter-act with similar tools and feel their effect on the tissue, as well as see it. Figure 1.1 shows two examples of medical training simulations in current medical research today. Figure 1.1(a) shows a needle insertion simulation [34], and Figure 1.1(b) shows a laparoscopic training simulation [17]. (a) Needle Insertion (b) Laparoscopy Figure 1.1: Examples of Medical Training Simulators with Haptic and Graphic Feedback1 In general, medical simulation systems have three basic components: 1. Haptic Interface: The haptic device is the hardware with which the user physically interacts. These devices are instrumented so that the position and orientation of the tool can be calculated, and some or all of the joints are actuated so that force feedback can be displayed as the user interacts with the virtual environment. Because human touch is very sensitive, the force feedback must be updated at high frequency, with a minimum of 300 Hz [24] for free arm movement. Several different haptic interfaces have been used in a variety of medical training simulators. For example, a 3 degree of freedom (DOF) planar pantograph is used in a simulation of needle insertion as shown in Figure 'Needle insertion photo is reproduced with permission from Simon DiMaio, Laparoscopic photo is reproduced with permission from Iman Brouwer, 1.1 Background 3 1.1(a). For minimally invasive abdominal procedures, a Laparoscopic Impulse Engine® interface is used for the laparoscopic simulator as shown in Figure 1.1(b). For bone dissection [3] or suturing [15] a PHANToM® [51] haptic pen is used (Not shown in Figure). 2. Real-time Simulation of the Virtual Environment: The real-time simulation is responsible for com-puting the interaction physics between virtual tools and virtual objects as the user moves the input device. Computational algorithms for collision detection, deformation, and dynamics are carried out by the simulation for display to the user. 3. Graphic Display: Where applicable, visual feedback is shown to the user. The visual feedback can come in a variety of forms. For the example of minimally invasive abdominal surgery, a camera mounted to a flexible optical line is used to show the surgeon what is happening in the surgery site on a monitor. To mimic this procedure, a monitor is used to display visual feedback to the student [10]. In removal of air cavities under the skin behind the ear, a procedure called mastoidectomy, the surgeon looks at the surgery site through a stereoscopic microscope. To mimic this method of graphical feedback, a set of video goggles is used [3]. Both of these examples use a mono 2D display. An alternative method for display is to use stereo glasses to provide 3D visual feedback, as used in the Zeus [49] surgical robotics suite. According to [24], 24Hz is sufficient for the visual update rate. These three components work together to allow the user to interact with a virtual object, feeling the simulated physical forces through the haptic interface, and seeing the motion and deformation of the object on the graphic display. Although graphics and haptics are updated at different rates, they are synchronized so that the touch and visual cues combine to yield a convincing feeling that the user is actually deforming the object on screen. The flow of information in a virtual environment forms a closed loop [25], as illustrated in Figure 1.2. The model deforms according to the user's input; this deformation results from a deformation force, and the loop is closed when this force is exerted on the user by the haptic interface. The loop is also closed with a visual representation of the deformation. 1.1 Background 4 f \ -Sample Virtual Tool Position -Collision Detection -Compute Penetration -Compute Deformation -Integrate Force -Output Deformation to Graphics Display -Output Force to Haptic Interface Figure 1.2: Data Flow in a Typical Haptic Virtual Environment 1.1.2 Deformable Models for Tissue Simulation The work presented in this thesis focuses on modeling deformable objects in the virtual environment for tissue simulation. Closed-form solutions exist for the deformation of simple geometric shapes such as beams or cylinders. However, anatomy can rarely be modeled using simple geometry because shapes are typically quite complex. As a result, the anatomy must be discretized into a set of volumetric elements, and the deformation solved numerically. Several methods for computing deformations are available from the literature. Elasticity based meth-ods such as the Finite Element Method (FEM) [16, 34] and Boundary Element Method (BEM) [53] use continuum mechanics to describe the deformation of an elastic body that has been discretized using either volumetric (FEM) or surface (BEM) elements. Due to it's popularity in structural engineering, the FEM has a wealth of research to build upon. These methods are well suited for accurate simulation of tissue because elasticity parameters from experiments can be incorporated into the model easily. However, these methods are very computationally expensive and require precomputation and numerical optimization for real-time applications [41]. Alternatively, solid objects can be modeled as a mesh of interconnected springs/damper assemblies and point masses. These Spring Mass System (SMS) based methods are quite efficient because of the simplicity of the equations of motion and are very prevalent in computer animation [41]. Unfortunately, the deformations are very sensitive to the spring parameters, and it is difficult to achieve tissue like-behavior because the deformation is not based on accurate physical modeling [31]. 1.2 Motivation 5 1.1.3 Fluid Simulation It is common in medical procedures that one has to interact with fluid as well as solid materials. To model fluid-filled structures, such as blood vessels, cysts, or glands, fluid physics must be computed using fluid mechanics analysis. Fluid mechanics can be solved using either static or dynamic analysis. Fluid statics is often used in civil engineering for the static analysis of structures like dams or pressure vessels. Static fluid analysis does not include any time varying quantities, so the only quantity used is hydrostatic pressure. In fluid dynamics, however, fluid motion is considered, and quantities such as velocity and hydrodynamic pressure are used to describe the flow of fluid [62]. Both static and dynamic fluid analysis have been used in medical simulations. Any fluid simulation that requires flow, such as bleeding [9] or removal of water and bone paste [2], can be modeled using Navier-Stokes based fluid dynamics. The drawback of using fluid dynamics is that it is very computation-ally expensive. Due to the computational expense, fluid dynamics are solved only for graphical presentation, which requires a much slower update rate than haptics does. However, it has been suggested that the de-formation of fluid-filled structures, such as blood vessels or organs, could be modeled using fluid statics to obtain the high update rates necessary for haptics [30]. This thesis will demonstrate that fluid statics can indeed be used to model the deformation of fluid-filled structures for haptic rendering. 1.2 Motivation The motivation of this work is to increase the realism and scope of medical haptic simulation by mod-eling fluid structures enclosed in bodies of elastic material. Since many medical procedures occur either in close proximity to, or directly with, fluid filled structures, the contribution of fluid forces to the surrounding elastic media should be modeled to make medical simulations as realistic as possible. Prostate brachytherapy is an excellent example of a medical procedure that is performed in close proxim-ity to a fluid-filled structure. This procedure is a treatment for prostate cancer where many small radioactive seeds are implanted with a long needle into the prostate gland to irradiate and kill the prostate tissue while avoiding the surrounding organs. The prostate gland is located inferior to the bladder, and at the time of the procedure the bladder is filled with a dilute contrast agent so that it will be visible in CT images [18,45]. 1.3 Research Objectives 6 Because the bladder is filled with fluid, the deformation of the bladder should be modeled as a fluid filled cavity surrounded by elastic tissue. If an extension of [34] to train for prostate brachytherapy is to be created, a fast model for simulating elastic bodies filled with incompressible fluid is required. To further motivate this work, procedures that take place near large blood vessels will have effects from the typically low pressure and easily deformable veins or, from high pressure arteries that exhibit very stiff behavior [1]. By modeling fluid structures enclosed in elastic media, blood vessels could be included in future medical simulations. 1.3 Research Objectives The first objective of this work is to develop a method for modeling fluid enclosures in elastic media for haptics that: 1. is fast enough for real-time simulation. 2. can accommodate any relationship between the fluid pressure and volume. 3. can accommodate a volume of elastic material in addition to the fluid body (ie. not a thin-walled model). The second objective is to show that the model predicts deformations in a realistic manner. To ac-complish this, a set of experiments is required to compare real deformations of a fluid-filled phantom to simulated deformations using the model. The third objective is to implement the model in a real-time haptic simulation to show that it is indeed fast enough to yield convincing results when deformed interactively. 1.4 Thesis Outline 1. Chapter 1 provides a brief description of haptic virtual environments for medical training, deformable methods, and fluid mechanics. The motivation for and research objectives of this thesis are discussed. 1.4 Thesis Outline 7 2. Chapter 2 describes the problem that this thesis attempts to solve. A literature review of the real-time modeling of deformable objects, haptic rendering of fluids, and fluid-structure interaction with the Finite Element Method are presented along with findings from early experiments. 3. Chapter 3 provides a brief introduction to the linear elastostatic Finite Element Method, which is used as the deformable model in this work. Static fluid analysis and its application to pressure boundary conditions are discussed, and the numerical method used to compute the desired pressure/volume relationship for an incompressible fluid inclusion is described in detail. 4. Chapter 4 describes the experimental validation of the incompressible fluid structure model. The phantom, data acquisition, and modeling of the experiments are discussed in detail. 5. Chapter 5 describes matrix condensation, the haptic algorithm, and a real-time haptic simulation of an incompressible fluid inclusion that was created based on the 2D needle insertion simulation described by [34]. Issues in 3D implementation of the haptic algorithm are also discussed. 6. Chapter 6 discusses the conclusions from the work presented in this thesis as well as recommendations for future work. Chapter 2 Problem Formulation In this chapter, a literature survey in the areas of deformable methods for haptics, fluid modeling for haptics, and fluid-structure interaction using the FEM is presented and discussed. From the conclusions of this survey and results from preliminary work, a problem statement is formulated, and a solution is proposed. 2.1 Deformable Methods in Haptics The literature shows a multitude of methods for describing the deformation of objects in haptic simula-tions. In this section, they are broken up into continuum mechanics based and spring mass based. Spring Mass Systems (SMS) are very prominent in computer animation, particularly for dynamic cloth simulation [8]. Because equations of motion are quite simple, yielding a relatively low computational cost [31], they are also an attractive method for haptic rendering of deformable bodies. d'Aulignac [28] uses a SMS system to model the deformation of the human thigh in a virtual echographic simulator. Cotin et al. [26] use a combination of FEM and SMS to simulate laparoscopic surgery of the human liver. The authors note that SMS simulation is faster to compute and simpler to modify mesh topology than with the FEM. The major disadvantage of SMS methods is that deformations are very sensitive to mesh topology and spring parameters. Much of the SMS related research is directed at computing the parameters so that the 8 2.1 Deformable Methods in Haptics 9 deformation of the SMS model follows a set of training data or elasticity parameters. For example, Deussen et al. [32] used simulated annealing to obtain spring parameters such that the SMS would mimic elastic deformations based on a FEM solution. d'Aulignac et al. [29] used a two step nonlinear least squares approach to match SMS parameters with a set of human thigh deformations taken with a PUMA robot arm. Crochemore et al. [27] used an evolutionary genetic algorithm to train SMS parameters for a simulation of nonlinear elasticity. Continuum mechanics based methods, like the Finite Element Method (FEM) and Boundary Element Method (BEM), are prominent because of their accuracy and long history in engineering analysis of mate-rials. However, these methods are computationally very intensive and require considerable precomputation and optimization for real-time display. James and Pai [52, 53] used the BEM to precompute the Green's Functions of linear elastic material, which relate a field of displacement vectors to traction vectors on the boundary [56]. In this method, the Green's Functions were used to render real-time interaction of a trac-tion vector with an elastic body using linear superposition and low-rank matrix update techniques. This method works well for non-invasive haptics because after precomputation over the domain, only a boundary representation is used. As a result, it would be inefficient to simulate procedures that require sub-surface interaction because local domain decompositions in the area of interest would be required. In a similar pre-computation approach, Cotin and Delingette [24] precomputed elementary deformations of a linear FEM analysis of a human liver. Section 2.1.1 is dedicated to the use of FEM in haptics. Two survey papers by Delingette [31], and Gibson and Mirtich [41] give an excellent overview of de-formable methods for computer graphics and medical simulation respectively. Delingette [31] observes that a different rate of computing is required for different medical applications. Planning applications that are not time critical can use larger, more accurate models that include the nonlinear nature of living tissues [39]. Surgical training with haptics, however, requires very fast update rates, and linear models with considerable precomputation and optimization are used. Delingette provides a comprehensive trade-off list between two very prominent deformable methods in haptic interaction, FEM and SMS. Delingette concludes that FEM provides an accurate description of deformable objects, but for real-time deformation, only linear and pre-computed models are fast enough. Topology changes from suturing or cutting are difficult with FEM in real-time because these changes require alteration to the often condensed or precomputed stiffness matrix. With regards to SMS models, Delingette concludes that SMS methods can solve quickly and accommodate 2.1 Deformable Methods in Haptics 10 topology changes easily at the cost of physical accuracy. Also, the 3D "ChainMail" algorithm for fast com-putation of volumetric deformation is discussed [40]. Although "ChainMail" is well-suited to accomodate topology changes, it is unclear how realistic the deformation is. Gibson and Mirtich [41] survey a wide range of deformable methods for use in computer graphics that extends from purely geometric representations like NURBS and B-splines to physically based methods like FEM, SMS, and approximate continuum models. Because of space constraints, only the physically based models are summarized in this survey. As in [31], the authors concentrate on the trade-offs between FEM. and SMS and come to similar conclusions. In addition to FEM and SMS, however, approximate continuum models like Snakes [55] and Discretized Deformation Energy [77], and low degree of freedom models like Modal Analysis [64] and Dynamic Global Deformations [78] are discussed. Although there exist many methods that are faster than the FEM, the authors conclude that for applications like surgery simulation, FEM is the best choice because it can accurately describe the deformation of nonlinear tissue. To improve the real-time performance of FEM, multi-resolution, hierarchical, and combined kinematic-dynamic [40], techniques are suggested. 2.1.1 Finite Element Method in Haptics One of the acceleration techniques that allows FEM analysis to run in real-time haptics is matrix conden-sation, described by Bro-Nielsen and Cotin [16]. It is important to emphasize that condensation techniques are only applicable to linear systems, hence linear elastic, linear strain analysis is the only method that can be accelerated using condensation. Matrix condensation allows for the typically large and sparse system matrix to be reduced to a much smaller and denser sub-matrix that contains only the entries of the nodes of interest. Once condensed, this smaller system can be used to describe the deformation and force of the condensed nodes very quickly. Selective Matrix Vector Multiplication (SMVM) can be used to compute the deformation of the entire body very quickly. Typically only a select few nodes actually have force ap-plied, so SMVM allows for only those nodes that are active to be multiplied through, eliminating all zero multiplies. Many researchers have used matrix condensation with FEM in haptic simulations of medical procedures. DiMaio and Salcudean [34] used the FEM and condensation techniques to simulate needle insertion interac-tively. In their work, the global stiffness matrix is inverted (precomputation stage), then the required matrix 2.1 Deformable Methods in Haptics 11 entries are put into the condensed matrix, which has only the nodes on the needle. Boundary condition changes are used to allow for force and displacement boundary conditions to co-exist without a full matrix inversion. Also, in [35], a method for changing mesh topology in the condensed matrix without a full matrix inversion is presented. This simulation was implemented on a 3DOF twin pantograph [36]. Berkley et al. [15,42] used banded matrix algorithms with Gaussian Elimination to precompute and condense linear FEM analysis of soft tissue for a suturing simulation. In their work, the global stiffness matrix is condensed to include only the nodes that are on the boundary of the suturing site. The author notes that this method is not capable of mesh topology changes without doing the precomputation steps again. The simulation was implemented using a PHANToM® interface, and initial subject tests by surgeons revealed that the linear FEM "felt" real. Basdogan and Srinivasan [10] used a combination of FEM and SMS to simulate laparoscopic common bile duct exploration. Linear, triangular FEM elements were used to model the cystic duct, while SMS was used to model the catheter. Matrix condensation was used to reduce the FEM system matrix to include only the translational degrees of freedom from a mesh of 3D triangular elements. Triangular elements require rotational constraints at each node when used in 3D, and the authors note that the rotational variables are unnecessary for the deformation in the simulation. This simulation was implemented on a laparoscopic interface. ' As an alternate method to condensation, it is possible to precompute FEM deformations in a method similar to that used by James and Pai [53] so that the solution can be rendered for real-time haptic display. Cotin et al. [24,25] used linear FEM to precompute and save the elementary deformations at each node. At run-time, the force and deformation are computed using superposition of the elementary deformations. This work was then extended to include quasi-nonlinear deformation by adding radial and tangential deformation corrections based on the experimental data found by Chinsei and Miller [20] It is possible to render deformable objects for haptic display without solving for the deformation at a high haptic update rate. Instead, the results from the model can be interpolated to achieve an update rate that is faster than the solution time of the model. Zhuang and Canny [80,81] describe how nonlinear strain analysis of elastic bodies can be used to model deformation in a haptic simulation. In their work, the haptic and graphic display are delayed by one timestep so that a simple linear interpolation between states tn and r n + i can be used. The authors note that the several millisecond delay is within the tolerance of human 2.1 Deformable Methods in Haptics 12 perception. While this approach may work for deformable objects, it would be difficult to simulate stiff objects like walls or bones because the maximum stiffness for stable haptic interaction is limited by the update frequency [70]. 2.1.2 Fluid Modeling in Haptics Applications Most of the present medical haptic simulators, including the commercial products from Reachin1 and Surgical Science2 include fluid simulation only for graphical feedback, to visually indicate bleeding or motion of fluid generated by surgical tools. In this case, the fluid simulation does not contribute to the haptic force, and is therefore only required to be updated at typical graphical rates of order 20Hz. For a laparoscopy simulator, Basdogan and Srinivasan [9] describe a method of incorporating the simpli-fied 2D Navier-Stokes based wave equations from [54] and [63] with deformable objects using the method of auxiliary surfaces. In this method, a 2D fluid grid is projected onto a deformable object, and the fluid simulation is carried out using the resulting depths from the projection. Because fluid motion is only fed back graphically, it is updated every 0.1s. In a simulation of mastoidectomy, Agus et al. [3] used a combination of the solid mechanics based animation methods described by [58] and [75] to show the motion of water and bone paste. This simulator used a decoupled simulation architecture where haptic force computation was performed at 1000 Hz in the "fast" subsystem, and fluid simulation and bone removal were performed at 20 Hz in the "slow" subsystem. In contrast to the previously mentioned simulation techniques that used fluid only for graphic feedback, De and Srinivasan [30] propose a nonlinear FEM formulation to model human organs as thin elastic walls enclosing incompressible fluid. This method was not implemented in real-time, and it's applicability to real-time simulation is based solely on the reduction in dimensionality from a volume to a surface. That is, to represent a volumetric deformable object using only a thin surface reduces the complexity from 0(n3) to 0(n2). If a body of elastic solid is required in addition to a thin walled organ, this reduction in dimensionality is no longer achieved, and the solution time will likely be too large for real-time applications. The authors note that a solution using 64 triangular elements and 41 nodes take the order of 2M floating point operations 1 2.2 Fluid-Structure Interaction with the Finite Element Method 13 (FLOPs) to solve. This is, in fact, a considerable computational cost for a system with only 41 nodes, and this result will be discussed in Chapter 5. Laugier et al. [57] used the Long Element Method (LEM), developed by Balaniuk et al. [7] in a surgical simulation. Despite using physically based principles like incompressible fluid pressure and density, the LEM does not describe any physical phenomena and is unsuitable for haptic fluid simulation. 2.1.3 Justification for choice of Finite Element Method FEM is chosen as the deformable model for this work for the following reasons: • It is physically based on the continuum mechanics of elasticity. This allows for elastic parameters like Young's Modulus (E) and Poisson's Ratio (v) to describe the object accurately. • Linear FEM can be condensed and run at typical haptic rates without the need for interpolation, even for fairly large and complicated models. • FEM uses volumetric elements, like tetrahedra in 3D or triangles in 2D, to discretize elastic bodies rather than only a surface representation as with the BEM. This is important because sub surface in-teraction is often required in surgical simulation. Even if the interior nodes are condensed out to begin with, they can be always be added to the condensed matrix without any considerable computational cost. • There is a large body of literature in fluid-structure interaction using FEM from which to draw ideas and methods. • As computational power increases, more accurate and/or complex nonlinear models will be able to run at haptic rates. 2.2 Fluid-Structure Interaction with the Finite Element Method As previously mentioned, there is a large body of engineering literature that is focused on the interaction of elastic structures with fluid bodies, and the general research area is referred to as fluid-structure interac-2.2 Fluid-Structure Interaction with the Finite Element Method 14 tion. There are two methods with which the fluid can be modeled, using the mixed Euler-Lagrangian Frame or the Lagrangian Frame. 2.2.1 E u l e r i a n - L a g r a n g i a n F r a m e The Eulerian approach to finite element analysis assumes that mesh nodes are static. This differs from the Lagrangian approach which is based on mesh node displacement. When fluid flow is modeled with the Lagrangian approach, the mesh moves according to the flow and needs to be re-zoned to alleviate mesh distortion, adding considerable computation and momentum conservation errors that prevent the accurate prediction of flow around comers. Consequently, fluid flow is usually modeled with an Eulerian approach, using velocity as the nodal quantity and pressure as driving potential energy [14]. To accurately model dynamic fluid flow coupled with structural mechanics, a mixed Eulerian-Lagrangian approach is a prominent method in current FEM research. In this approach, the fluid-structure coupling is achieved by matching normal velocities over a time-step (in effect a displacement) with the structural displacements at the contact nodes. The resulting matrices are asymmetric and are poorly banded. Conse-quently, special solvers are required to solve the system, introducing a large computational effort [43]. Bathe et al. [13,76] modeled the fluid-structure interaction of several common engineering systems with the ADINA 3 finite element package. The authors notes that to fully couple the fluid analysis with a structural analysis, the Eulerian-Lagrangian approach is needed for the fluid, while the Lagrangian is used for the structure. According to [11] and [48], iterative solution methods such as the Gauss-Siedel or Preconditioned Conjugate Gradient are used to solve the poorly conditioned system of equations. DiMaio [35] implemented an Eulerian-Lagrangian fluid-structure simulation. In his work, the stiffness matrix, which contains viscosity and pressure coefficients for the fluid and stiffness coefficients for the solid, is near singular, so an implicit iterative algorithm is used to solve the system. The MATLAB® Optimization Toolbox was used to solve the set of equations with the Levenberg-Marquardt solution method, which entailed considerable computational cost. This type of simulation could not mn in real-time on present hardware. 2.2 Fluid-Structure Interaction with the Finite Element Method 15 2.2.2 Lagrangian Frame The Lagrangian approach to fluid analysis has been used in civil/mechanical engineering for many years to solve hydrostatic and hydrodynamic problems. A fluid body can be modeled using modified elastic elements and easily integrated into a standard structural analysis formulation. FEM in structural analysis is done assuming the material is isotropic, meaning it has the same elasticity parameters in all directions. This allows the deformation of a body to be characterized by two parameters. The most common sets of parameters are: 1. Young's Modulus, E, and Poisson's Ratio, v. Young's Modulus is physically interpreted as the stiff-ness of the material, while Poisson's ratio is a measure of the compressibility. 2. Lame Parameters X and fi, where /i is the Shear Modulus that characterizes the ability to resist shear strain. According to [4], an irrotational fluid element can be created by setting the Shear Modulus to zero. This results in an element that cannot resist shear strain, and therefore the fluid's elasticity is defined solely by the Bulk Modulus, which characterizes the ability to resist volumetric strain. Several authors have reported good results from the use of irrotational fluid elements in static and/or dynamic analysis of common engineering structures such as water tanks [37], dams [19], and pipes [12]. In each work, a different method for releasing the tangential displacements at the fluid-solid interface was used, enforcing the 'slip' condition. [19] used constraint equations, [12] used a thin layer of the irrotational elements, and [37] used short, stiff truss elements. Figure 2.1 illustrates the truss method. Figure 2.1: Truss Elements to Decouple Tangential Displacements 2.2 Fluid-Structure Interaction with the Finite Element Method 16 Preliminary tests with irrotational fluid elements and truss elements to enforce the slip condition showed that the trass elements only decouple the tangential component for displacements as large as the length of the truss. This is not a problem for stiff structures like dams or pipes that experience small deformations, but it is a problem for soft tissue where deformations are much larger. The same would be true if thin fluid elements were used to decouple tangential displacements. The multipoint constraints were not attempted because the resulting system would require penalty methods. Another problem with using irrotational fluid elements that was uncovered in preliminary tests is that the stiffness matrix entries are poorly scaled and the matrix is ill-conditioned when it contains entries from both soft tissue and fluids like water because water has a Bulk Modulus roughly 5 orders of magnitude larger than the Young's Modulus of soft tissue. This effect worsens as the fluid tends to incompressibility. Due to problems enforcing the slip condition for large deformations and the inability to model near incompress-ible fluid in contact with soft tissue, irrotational elastic fluid elements are not suitable for modeling fluid structures in haptic medical simulation. In an alternate approach to model fluid-structure interaction with FEM, Shiraz-Adl and Shrivastava [72] used a nonlinear formulation similar to that used by De and Srinivasan [30] in a static analysis of various structures enclosing fluid. The effects of fluid pressure are added to the FEM formulation in a tangent stiffness matrix, and the virtual work done on the fluid creates an extra term in the FEM formulation. A constraint is imposed such that the fluid cavity exhibits a constant volume. The authors note that it would be possible to use the bulk modulus of a compressible fluid in the formulation as an alternative to the constant volume constraint. The authors also note that changes in fluid pressure could be accommodated if they are known a priori. Their method is well suited to model tissue that encloses fluid structures because it does not require slip condition enforcement nor does it exhibit numerical problems due to incompressibility. However, because it requires a nonlinear formulation, computation time will be too high to solve at haptic update rates in excess of 300 Hz, and an interpolation scheme similar to that described in [80,81] would be required. 2.2.3 On Incompressible Finite Element Analysis It is often necessary in engineering applications to model incompressible materials. In fluid flow, low mach number flows can usually be assumed incompressible, hence removing the compressibility terms from the 2.3 Problem Statement 17 Navier-Stokes equations [62]. This poses no problem for finite element modeling of fluid flow in the Eulerian Frame. However, it is well discussed in literature that the deformation of incompressible materials like rubber or plastically deforming metals cannot be solved using Lagrangian linear finite element analysis. As the Poisson's ratio tends towards 0.5, the stiffness matrix becomes over constrained, and the solution exhibits "mesh locking" [65]. Higher order elements can be used to compute incompressible deformation [59,65]. 2.3 Problem Statement The literature has not shown a method for solving the deformation of an elastic body that encloses fluid structures that is fast enough for real-time haptic simulation. Motivated by improving medical simulation technology, the focus of this thesis is to develop, test, and implement a method that is capable of simulating soft tissue with fluid filled anatomical structures at real-time haptic rates. Furthermore, the literature review concerning incompressible fluid has shown that the simulation of incompressible media with linear FEM yields inaccurate results. Consider an arbitrarily shaped elastic body Qe, as illustrated in Figure 2.2. This body may contain pockets of fluid, Q,f\ and Q/2, each of which has a quasi-static relationship between the pressure and volume that is known. The exterior boundary, may be constrained (dQ,d), and traction or displacement boundary conditions may be imposed (dQ.n). The problem is to compute the resulting body deformations and forces at a rate that is fast enough for haptic simulation which requires force updates to be of order 300 9Q, Figure 2.2: Fluid Filled Structures Enclosed in Elastic Media Hz. 2.4 Proposed Solution 18 2.4 Proposed Solution The proposed solution is to use linear FEM to simulate an elastic body with a cavity, and to use hy-drostatic fluid pressure to simulate the presence of fluid inside the cavity. Hydrostatic fluid pressure is the simplest of the methods for modeling fluid described in the literature review, but despite being simple it yields accurate results for fluid-structure interaction problems. The work described in [30] uses the hy-drostatic pressure method for simulating fluid with a nonlinear FEM formulation. The major drawback of including the fluid pressure in a non-linear FEM formulation is that the solution cannot be condensed. In-stead of including work done on the fluid in a nonlinear FEM formulation, pressure force will be applied as a boundary condition to a linear elastic FEM formulation. This completely abstracts the fluid simulation from the FEM, and because the stiffness matrix of a linear FEM system can be condensed, fluid computa-tions can be done only on the fluid exposed nodes. This reduction in the system size will allow the fluid simulation to run at real-time haptic update rates. A further advantage to using condensed linear FEM is that a large volume of elastic solid can be included in the system without affecting the real-time performance. The nonlinear FEM based thin-walled model described in [30] requires iteration on each elastic node, so if a volume of elastic solid is required in addition to the fluid structure, the solution time will likely be too long for haptic display. 2.4.1 Drawbacks and Limitations There are several drawbacks to modeling fluid filled structures with linear FEM and hydrostatic pressure boundary conditions, the first of which stems from the linearity of the FEM methods. Linear elastic FEM is a coarse first approximation to the highly nonlinear deformation of living tissues [39]. This effect is worsened by the fact that to maintain a linear FEM formulation, the linear strain tensor is required where it is assumed that strains are small. The primary justification for using linear FEM is that it can be condensed and solved very quickly. A second limitation is that the fluid-elastic interaction is assumed to be quasi-static. Although the simulation runs in real-time, at any update step only a static force balance is maintained. For example, the pulsing of an artery could be modeled using a time varying pressure, but the pressure would be applied to a static solid. As a result, sloshing of fluid or vibration of the elastic material caused either by user interaction 2.5 Summary 19 or biological forces is not considered. Again, this assumption is made to allow for fast computation time. 2.5 Summary A review of the literature concerning deformable methods for haptics, fluid modeling for haptics, and fluid-structure interaction with FEM was presented. No method was found that will allow for a body of elastic material that encloses fluid structures to be computed at haptic update rates. A method that uses linear FEM in conjunction with an iterative solver to enforce a fluid volume/pressure relationship has been proposed. The limitations of the proposed method were discussed. Chapter 3 Methodology The fluid inclusion model presented in this thesis is based on the FEM for linear elastic materials. This chapter will describe the use of the FEM by briefly describing the continuum mechanics based derivation, volume discretization using finite elements, boundary conditions and constraints. Once the FEM has been described, the fluid model is presented by describing the distribution of static pressure over the FEM mesh, and the numerical method that solves the fluid model. The method's performance and robustness are then discussed. Finally, volume calculation from a closed polyhedral surface is described. A detailed description of all FEM methods used, including mesh and stiffness matrix generation, is provided in Appendix A. 3.1 Finite Element Method for Linear Elasticity The Finite Element Method for a homogeneous, isotropic elastic material is based on generalized Hooke's law: a = De (3.1) where a and e are the stress and strain vectors on a body, respectively, and D is the elasticity matrix which is parametrized by two independent elasticity variables. In this thesis, the elasticity variables that are used are Young's Modulus, E, and Poisson's Ratio, v. Young's Modulus is a measure of the materials stiffness, and can be determined experimentally from the slope of a stress vs. strain curve. Poisson's Ratio describes the compressibility of the material and ranges from 0, for very compressible materials, to 0.5, which represents 20 3.1 Finite Element Method for Linear Elasticity 21 incompressible materials. From 3.1, the strain energy of a linear elastic body, Q, can be computed by: Estrain = \ f £TOdx (3.2) 2 JQ Static equilibrium between deformation energy and external force is achieved when the first variation of E(u)strain vanishes. This occurs when SE(u)strain = 0 and corresponds to the potential energy of the system reaching a minimum value [16]. &E(u)slrain = 0 = / BTeDBeudx-f (3.3) JQ. where / is the discretized force vector and Be is a constant matrix defined in Appendix A. (3.3) is integrated numerically and interpolated using shape functions over each discrete volumetric element to generate a stiffness matrix for each element. Each stiffness entry is added to the global stiffness matrix, K, which creates a system of 3n (for a 3D object) linear equations that describes the deformation of each of the n nodes by the linear system: f = Ku (3.4) where u and / are the displacement and force vectors, ordered based on node number such that: f\,x f\,z fn,x fn,y fn,z M l , z (3.5) 3.1 Finite Element Method for Linear Elasticity 22 3.1.1 Volume Discretization In order to compute the global stiffness matrix for an elastic body, it first must be discretized into finite volume elements. Figure 3.1 shows a 2D example of an elastic body with a hole in the middle, that has been discretized using triangular finite elements into a FEM "mesh". The mesh presented in Figure 3.1(b) is created by first generating a set of regularly spaced nodes, then generating a connectivity list that describes which nodes appear in each triangle. After the regular mesh is created, nodes are removed from the middle of the mesh to create the hole. When nodes are removed, triangles that include these nodes must be removed as well (see Appendix A.2 for a detailed description of this algorithm). Finally, the nodes on the boundary of the cavity are moved to lie on a circular perimeter. The FEM meshing done in this thesis is done using the same general algorithm of generating a regularly spaced mesh and then removing and moving nodes to match the desired geometry. All of the FEM work done in this thesis was done using MATLAB® scripts. : Empty ' Cavity Elastic Uodv I (a) Elastic Body with Fluid Inclu- (b) Discretization of Elastic Body sion with Triangular Elements Figure 3.1: Elastic Body Discretized with 2D Triangular Finite Elements A regularly spaced mesh works well for some forms of medical imaging like 3D ultrasound, MRI or CT. Because these modalities generate regularly spaced planes of image data, it is possible to attach mesh to each plane, move nodes to match up with anatomy, and discretize the volume using predefined connectivity. This approach is used to generate 3D mesh for the FEM analysis in this thesis. However, it is important to note that the use of a regular mesh for complex 3D objects could be difficult. As an alternative to regular / / / / / / / / / / / 3.2 Fluid Model 23 spacing, a variety of automatic mesh generation programs like NETGEN 1 or GRUUMP 2 and many volume reconstruction algorithms, [6], that are capable of generating 3D tetrahedral mesh could be used. 3.1.2 Constraints and Boundary Conditions The FEM described in this chapter uses equilibrium equations that are based upon elastic deformation. If the object is not constrained in each DOF, an external force will cause the object to translate, and elastic FEM is not capable of solving for motion. As a result, the elastic object must be constrained in all DOFs to ensure that the object cannot translate or rotate. Appendix A. 1.4 describes the procedure and requirements for implementing appropriate constraints. Once the stiffness matrix has been sufficiently constrained, boundary conditions can be applied to com-pute deformation. In most engineering applications, the load vector, / , is known from the application, and the displacement vector, u, is unknown. The implicit system Ku = f can be solved for the unknown displace-ments either by a direct method like Gaussian Elimination or by an iterative matrix solver, like Conjugate Gradient. An alternate approach is to invert the stiffness matrix and and solve the system explicitly. This is usually avoided due to numerical errors from inversion, long inversion time, and much larger storage requirements for the inverted matrix [16,22]. 3.2 Fluid Model The fluid model presented in this thesis is based on static fluid analysis, where the only contribution from the fluid to the elastic body is hydrostatic force [62]. Consider a fluid cavity that is enclosed by elastic media and at rest, as shown in Figure 3.2. The pressure inside the cavity manifests itself as a distributed force normal to the fluid elastic boundary (shown in Figure 3.2 by arrows). Fluid obeys the laws of gravity, and hydrostatic pressure increases with depth. The pressure difference, AP, varies linearly with the height of the fluid, h, the fluid density, p, and the gravitational constant, g, by: 1 2http.7/ 3.2 Fluid Model 24 A P = pgAh ( 3 . 6 ) If the fluid object is very small, or the fluid has a low density, the effects of gravitational pull can be neglected. For example, cysts, glands and fat pockets are of the order I0~2m in dimension, and have a negligible pressure difference across their volume. If the fluid is at high pressure, such as the blood pressure of arteries, the pressure difference across the vessel is also negligible. However, for larger organs that are filled with dense fluid, like the bladder, the gravitational effects should not be neglected. / • / / / / / / / < • -• / > (a) F E M Mesh of Elastic Object with Computed Pressure Normals Figure 3 . 3 : Finite Element Solution of (b) F E M Solution with Static Fluid Pressure Load Object with Static Fluid Pressure Load To illustrate the effect of fluid pressure applied to an elastic body, Figure 3.3(a) shows an elastic body that has been constrained on the left vertical edge with pressure normal vectors for each of the elements 3.2 Fluid Model 25 that is on the fluid-solid boundary. Figure 3.3(b) shows the deformation that results from a static pressure force being applied to the elastic body. Note how the cavity has enlarged, and the elements on the perimeter appear stretched. 3.2.1 Hydrostatic Fluid Pressure as a F E M Boundary Condition The fluid model used in this thesis is based entirely on the principle that fluid pressure can be treated as a force boundary condition on an FEM discretization of an elastic body. Figure 3.4 illustrates how pressure is distributed along one edge of a 2D finite element. Even though 2D FEM is planar, each element is computed with a thickness, r, so each line segment is actually a rectangular surface. The effective force that an element on the fluid boundary is exposed to is solved by Fe = JA^PdA, where Ae = Let is the area of the element. This pressure force is applied along its normal by: F. = PA„he (3.7) Where he is the element normal. The uniform pressure shown in Figure 3.4(a) is equivalent to F_g applied at the center of pressure of the element, as shown in Figure 3.4(b) [62]. Using static force analysis, the two nodes at opposite ends of the line segment each recieve half the force over the element, =f, which agrees with intuition. F. I l l l l i l l i l i l l (a) Pressure Distribution Le. 2 (b) Effective Pressure Force 2 2 (c) Pressure Force at Nodes Figure 3.4: Pressure Distribution over one Edge of a 2D Element 3.2 Fluid Model 26 The 3D FEM elements used in this thesis are tetrahedral elements. When an elastic body with a fluid enclosure is discretized with tetrahedra, the inner cavity surface is defined by a closed contour of triangular surfaces. It is therefore necessary to determine how a uniform pressure is distributed amongst the three nodes of a triangular surface. By using a sum of forces in the surface normal direction, and a sum of moments about either a point or two edges of the triangle, the force on each node can be determined [47]. A constant pressure on a triangular surface can be replaced with the an effective pressure force that acts at the centoid of the triangle [62]. (a) Pressure Distribution on a Triangle (b) Center of Pressure of a Triangle Figure 3.5: Center of Pressure of a Triangle The sum of forces in the normal direction yields: S£„ = £ < - £ i - £ 2 - £ 3 = 0 (3.8) The sum of moments along two edges: = E 3 h - ^ = 0 (3.9) S M 1 3 = £ 2 ^ - ^ = 0 (3.10) Solving these three equations yields: (3.11) A triangle has the special case that it's center of pressure is one third of its height away from any edge regardless of the geometry of the triangle, as shown in Figure 3.5(b). An analytical proof of this is given 3.2 Fluid Model 27 in Appendix A.3. Therefore the support forces at each node are one third of the pressure force provided the pressure is spatially constant. If the pressure is not constant, as it would not be if gravitational effects were considered, the center of pressure is located below the centroid in the plane of the triangle, and the distribution of force over the nodes will not be equal. Distributing hydrostatic pressure over the fluid boundary requires one pressure computation per element. One third of F^ is added to each of the three nodes on each tetrahedral surface triangle in 3D, or half of F\_ is added to the nodes of each triangle edge in 2D. 3.2.2 Relationships between Fluid Pressure, Volume, and Flow The fluid model described in this work reduces the effects of fluid to a pressure force that is applied over the fluid-elastic boundary. Although the fluid is assumed static in time, this does not limit this approach to modeling purely static phenomena. Because a static equilibrium is used for each update step, and each update step is computed over a finite time, this type of simulation is quasi-static, and it could be used to approximate dynamic fluid-elastic interaction for haptic rendering. The majority of the fluid modeling in this thesis deals with constrained incompressible fluid, in which flow is not a important feature. However, the fluid model can be applied to a variety of fluid systems, such as the flow of blood through the vascular system. As an approximation, vascular flow could be simulated as steady, inviscid fluid confined within a flexible pipe using (3.12) for changes in height, h, pressure p, and velocity V. The pressure effects from this flow could be added to a linear FEM simulation of an artery to create a fast simulation of the fluid-elastic interaction that occurs near arteries using human blood pressure and heart rate. Bernoulli's Equation, 1 , P + 2 (>V + Pgh = constant along streamline (3.12) holds for steady, incompressible flow where viscous effects are negligible [62]. Blood flow is neither steady nor are viscous effects negligible [39], but within the limits of a real-time haptic simulator, this simplified model could provide a reasonable approximation. 3.2 Fluid Model 28 3.2.3 Confined Incompressible Fluid Because incompressible fluid enclosures are both challenging to compute in real-time and have been iden-tified as an appropriate model for anatomical structures such as the bladder, it is discussed in detail in the following sections. An incompressible fluid will resist volumetric strain with hydrostatic pressure applied on the boundary. To simulate the presence of an incompressible fluid inside a cavity, hydrostatic force can be applied to the boundary to conserve the cavity volume. The problem is solving for the hydrostatic fluid pressure such that the cavity does not change volume. 0.01 [ml (a) At Rest 0.01 0.04 -0.02 0 0.02 0.04 [m] (b) Deformed, Cavity at External Pressure 0.01 -0.02 0 0.02 0.04 [m] (c) Deformed, Confined Incompress-ible Fluid Inclusion Figure 3.6: Fluid Relationships shown with FEM Figure 3.6 illustrates the relationship between volume and pressure for a fluid filled cavity in three configurations. Figure 3.6(a) shows an elastic body with a square cavity in its undeformed state. Figure 3.6(b) shows the same elastic body after an external load is applied by the arrow. This corresponds to the case where the cavity is filled with a fluid that is connected to the surroundings or is filled with a very compressible fluid, and hence there is no pressure difference across the elastic material. Because there is no pressure difference, no fluid pressure is applied to the boundary and the volume of the cavity is significantly reduced. Figure 3.6(c) illustrates the case where the same elastic object is filled with incompressible fluid that is confined withing the elastic object. Hydrostatic fluid pressure has been added to compensate for the volumetric strain on the fluid cavity that the external load imposes. Note that the volume of the cavity in 3.6(c) is the same as Figure 3.6(a) and the cavity has a slight bulge in the vertical direction from the fluid pressure. 3.3 Numerical Method 29 The incompressible fluid pressure in this example was set by trial and error such that the volume of the deformed cavity is within 5% of the cavity at rest. For simulation of a fluid structure, a numerical method that solves for the hydrostatic pressure quickly is required. 3.3 Numerical Method Given an external force, shown by an arrow in Figure 3.6(b), the fluid cavity will deform from an original volume, V0, to a new volume V, at sample i, as illustrated between Figures 3.6(a) and 3.6(b). The relative error at sample i is defined by: The error should be minimized to within a pre-defined tolerance from zero to solve for the incompress-ible fluid pressure. It is important to keep in mind that this method for modeling fluid is designed specifically for use in haptic applications where update rates are on the order of 300 Hz at typical hand speeds on the order of 1 cm/s. With a high sampling rate, the change in deformation between each sample should be small. For example, a 5 x 5cm block of tissue, and a hand speed of 1 cm/s would yield a strain rate of 20 %/s. At an update rate of 300 samples per second, the strain rate per sample would be approximately 0.06 %/Sample. With strain rates at this level, the change in incompressible fluid pressure required between samples is small, and an appropriate numerical method should be able to converge to the correct pressure quickly. The numerical method that is used to essentially "track" the correct incompressible fluid pressure is a proportional controller. 3.3.1 Proportional Controller Figure 3.7 illustrates a digital feedback loop controlled by a proportional controller that results in a pressure update equation of: Errorj =1 (3.13) o PI+X=PI + KpErrori (3.14) The FEM block starts with a solution to the elastic system given the pressure and external forces. This 3.3 Numerical Method 30 » Errori Controller + ) Pi. i+l Figure 3.7: Proportional Controller is followed by a computation of the volume of the cavity, V;, which is fed back to form a closed loop control structure. Because the performance of a proportional controller is dependent on the value of the controller gain, Kp, a system identification analysis is performed to estimate an optimal value of Kp. 3.3.2 System Identification and Controller Tuning To choose an appropriate value for the controller gain, a basic understanding of the system response is required. In this case, the system response is the relationship between the volume and pressure of a cavity in an elastic body. Figure 3.8(a) shows a square elastic body with a square internal cavity. The body is 5 x 5cm square with a thickness of 1.2cm, a Young's Modulus of 34kPa, and a Poisson's Ratio of 0.34. Figure 3.8(b) shows the pressure to volume response of this body. The data in this figure were generated by applying a static pressure to the cavity and computing the subsequent volume change. The pressure vs. volume curve is shown from approximately 30% reduction in volume to a 30% increase in volume, which is an estimated region of interest for haptic interaction with the fluid inclusion. These data were generated assuming that there is no pressure difference across the body when it is at rest (Figure 3.8(a)). The data in Figure 3.8(b) exhibits a nonlinear progression. However, over the region of interest, a good estimation to the pressure vs. volume relationship can be made with a linear fit. Figure 3.8(b) shows a linear fit to the data, which yields a slope of approximately 238 Pa/%. This slope should be used as the gain of the controller because over the region of interest it describes the pressure to volume relationship of this hollow elastic object quite well. 3.3 Numerical Method 31 [m] Volume Ratio (a) Elastic Body at Rest (b) Pressure vs Volume Curve Figure 3.8: Pressure vs. Volume of an Elastic Body 3.3.3 Performance and Robustness of Controller When the controller gain is set to the approximate slope of the region of interest from the pressure vs. volume curve, Figure 3.8(b), the performance of the controller is quite good. Table 3.1 lists different convergence requirements for different levels of strain. These values are generated from a simulation of the elastic body described above that allows the user to apply forces to the object interactively (though not in real-time) using the arrow keys. Approximate Strain (%) Iterations for 1% tolerance 0.3 Oor 1 3 1 15 1 30 1 or 2 Table 3.1: Performance of Controller The results demonstrate that when the controller is tuned to the approximate slope of the pressure vs. volume curve, it converges quickly to the correct incompressible fluid pressure, and is robust to very large deformations. 3.3 Numerical Method 32 3.3.4 Volume Calculation It is necessary to compute the volume of the cavity for each iteration numerically. The computation of the area of a 2D polyhedral shape can be computed by summing the areas of all of the trapezoids in order of connectivity. Figure 3.9 shows the calculation of an arbitrarily shaped polyhedron using a counter-clockwise order. y Figure 3.9: 2D Polyhedral Area Calculation by Summing the Areas of Trapezoids (3.15) expresses the sum of each trapezoid around the contour. It is important that the order be consistent so that trapezoids on the upper surface contribute positively to the area, and trapezoids on the lower surface contribute negatively. " n n y~^~y+i l " i=l i=l i=l 1 1 i=\ Computing the volume of an arbitrary 3D object that is defined by a polyhedral surface, and is not necessarily convex is not as trivial as in 2D. Brian Mirtich [60] developed a fast method for computing polyhedral mass properties based on a surface representation for graphics applications. This method uses three steps to reduce a volume integration to a set of line integrals over each face. The method first uses the Divergence Theorem to reduce a volume integral of a vector field F over a 3.3 Numerical Method 33 volume V, to an integral over a closed surface, dl/ by: / V-FdV = f F-hdA (3.16) J<V Jd<V If F = [x, 0,0], then V • F = 1, and the left hand side of 3.16 reduces to f^, 1 dV = V. The volume of the closed surface composed of discrete faces, f can be computed by: V= X hx f xdA (3.17) Next, the integral over a face, !F, in (x,y,z) is reduced to a projection integral over the bounded region, II, in the (a, P) plane where (a, P,y) is a right handed permutation of (x,y,z) chosen so that the projection onto the plane is maximized to minimize numerical errors. A plane in (a, P,y) space is defined by: n a a + «pP + nyY + w = 0 (3.18) So the face integral is reduced to a projection integral by: / / ( a , p , Y ) d A = ^ - / / ( a , p , / i ( a , P ) ) £ / a d P (3.19) where h(a, P) = (n a« + wpP + w) Finally, the projection integral is reduced to a set of line integrals by Green's Theorem in a Plane, which is a 2D version of Divergence theorem. f VHdA=<( Hrhds (3.20) Jn Jdn where m is the outward normal of each line segment, length s, that make up the bounded region, n, in a counter clockwise order. It is interesting to note that Equation 3.20 reduces exactly to Equation 3.15, so the 2D area calculation presented previously is an intuitive implementation of Green's Theorum in a Plane. 3.4 Summary 34 3.4 Summary This chapter briefly described the FEM for linear elastostatic objects and the mechanics of using FEM to simulate deformable objects. A method for modeling fluid enclosed in an elastic object using force contributions hydrostatic fluid pressure was described. A numerical method for enforcing a given pressure-volume relationship inside the fluid pocket was presented, and it's robustness and tuning were discussed. Chapter 4 Experimental Validation To validate the method for modeling a fluid pocket in an elastic body presented in this thesis, a tissue phantom that encloses a pocket filled with fluid is deformed with a known load while the contour of the pocket is imaged using ultrasound, and the top surface of the phantom is tracked by hand with a digitizing pen. The data from these experiments are then compared with results from a FEM simulation that mimics the experiment using the incompressible fluid model presented in Chapter 3. This chapter discusses each of the components in detail, and presents the results from the experiments. Because these experiments are meant to validate the hydrostatic pressure model of fluid structures for haptics, the FEM analysis that is used to model the experiment is done with the same methods that are used to create a real-time haptic simulation (as will be described in Chapter 5). The use of these methods has two impacts on the accuracy of the analysis: 1. Linear FEM is used despite its limited accuracy beyond approximately 1% strain. 2. The elastic system is solved explicitly, so numerical errors and long computational times are present due to the inversion of the stiffness matrix (the reason for this is described in Section A. 1.6). 4.1 Phantom A linear elastic phantom that encloses a pocket filled with incompressible fluid is required for validation. 35 4.1 Phantom 36 The elastic solid for this phantom is made by dissolving 13% B-Type Gelatin from bovine skin in de-ionized water. Gelatin is used as the phantom material because it is non-toxic, easy to work with, and linear elastic to approximately 12 % strain [46]. The fluid pocket is made from the fingertip of a latex glove filled with a 10% glycerol-water mixture. When the glove tip is being filled, care is taken to make sure that the latex is not very taut so the pocket is as compliant as possible, and its effect can be neglected from the FEM solution. The fluid pocket is suspended by threads in a 6.5 x 5.5 x 3.5cm3 mold when the molten gelatin is poured in, as shown in Figure 4.1(a). These threads are carefully removed prior to the experiments. Figure 4.1(b) shows the phantom after removal from the mold with marker dots drawn on the top surface. There is one major drawback to using gelatin as the phantom material. Gelatin has a tendency to fracture and crumble when it is deformed with high local stress. If a small probe or needle is used to deform the phantom, the gelatin in local area around the probe will crumble very easily. Because gelatin is used, the deformation is limited to force applied over a large surface to prevent the destruction of the phantom. (a) Phantom Mold and Suspended Fluid (b) Phantom After Removal from Mold Pocket Figure 4.1: Ultrasound Phantom Construction 4.1.1 Ultrasound Imaging Requirements To capture the shape of the fluid pocket under the various loads that the phantom is exposed to, it is imaged using an Ultrasonix 500 ultrasound imaging system [23]. A 9 MHz 40mm linear array probe is used to image the phantom both for mesh generation and for experimental data collection, and can be seen rigidly 4.1 Phantom 37 constrained to the plate on the left hand side of Figure 4.4 In order to image the contour of the fluid pocket with ultrasound, there needs to be contrast between the fluid and the elastic solid. The contrast texture that is available in ultrasound imaging is referred to as speckle, which is a result of the sound waves bouncing off of small scattering objects in random patterns. When there is sufficient speckle, the difference between a fluid body and a solid body is clear, because the fluid does not generate any speckle. 3% SigmaCell cellulose [68] is added to the gelatin so that sufficient tissue-like speckle is generated. Images from the ultrasound phantom are shown in Section 4.4.1. The pocket is clearly evident as a black void in the center. 4.1.2 Phantom Material Properties The elastic material properties are required as input parameters to the FEM simulation of experiments. As such, the Young's Modulus and Poisson's Ratio must be determined experimentally. The Young's Mod-ulus of a material can be determined by the slope of a stress vs. strain curve. An unconfined parallel plate compression is performed to generate the data set for the stress vs. strain curve. Figure 4.2(a) shows a sam-ple of the gelatin/cellulose in the compression apparatus. An ATI Nano-17 Force/Torque sensor is mounted vertically on a Newport 3DOF motion stage, and the sample is compressed between flat aluminum plates. 5000 (a) Parallel Plate Compression Appa-ratus C o m p r e s s i o n Exper iment Y o u n g ' s Modu lus = 1 5 . 2 k P a 10 15 20 25 30 Strain [%] (b) Compression Results Figure 4.2: Ultrasound Phantom Material Properties 4.2 Surface Tracking with Digitizing Pen 38 Stress and strain are related by Young's Modulus by: oax = Eeax (4.1) Gax — where Fc is the compressive force, A is the exposed area, is the axial stress, and is the axial strain. E is the Young's Modulus of the Material and the slope of the stress vs. strain graph, e = f; where d is the compressive displacement and h is the original height of the sample. According to results from [71] gelatin is nearly incompressible. It is assumed that volume is conserved in the compression test analysis such that: V n — d where A; is the area at the i'h displacement step, V is the volume of the sample, h is the original height of the sample, dj is the ith displacement. Force readings are taken at 0.1mm increments until approximately 25% strain, generating a set of force vs. displacement data. Using Equation 4.2, this data can be converted into a stress vs. strain data set which is show in Figure 4.2(b). A least squares fit to the linear region, to approximately 15 % strain, produced a slope of approximately 15.2 kPa, also shown in Figure 4.2(b). 4.2 Surface Tracking with Digitizing Pen The OptoTrak® 3020 optical tracking system is used to track the surface deformation of the phantom in 3D by using a digitizing pen. The OptoTrak® system has an accuracy of 0.1mm in the viewing plane, 0.15mm in depth, and a resolution of 0.01mm [50]. The digitizing pen consists of 5 IRED markers that are rigidly attached in a cross pattern to a sharp pointer. Figure 4.2 shows the digitizing pen in the foreground and the OptoTrak® 3020 camera system in the background. 4.2 Surface Tracking with Digitizing Pen 39 Figure 4.3: OptoTrak® Digitizing Pen 4.2.1 Calibration of Pen The pen is calibrated using a combination of the methods described by Rohling et al. [69] and Youwei Zhang [79]. Calibration is necessary so that at any instant, the position of the tip of the digitizing pen can be determined using the positions of the 5 IREDs. The process by which the pen is calibrated: 1. Rotate the probe around a stationary contact point S to collect several hundred samples of the IRED positions. 2. Use a least squares tit of the centroid from each sample to determine the point °S, where the superscript 0 means that the point is given in the world or OptoTrak® co-ordinate frame. 3. Use the non-linear least squares method described by [5] to determine a robust co-ordinate system with which PS is computed. The point of contact in the probe co-ordinate system. 4. Average the several hundred values of PS. Once the value of PS is determined, the position of the needle point can be determined using the position of the 5 IRED markings by: °S = °pTpS (4.3) where °T is the homogeneous transformation matrix from probe co-ordinates to world co-ordinates. 4.3 Experimental Procedure This matrix is generated using the centroid and orientation from [5] for every set of IRED positions. 40 4.3 Experimental Procedure Figure 4.4 shows the experimental setup used to measure deformation, and Figure 4.5 illustrates the setup in a schematic diagram. Figure 4.4: Complete Phantom Compression Apparatus Surface Markers-Ultrasound Probe Fluid Pop Plane Phantom O ( ) O Force Sensor Motion Stage Lubricated Surfaces ' 0 1 z Figure 4.5: Schematic Diagram The experimental procedure is as follows: 1. Make a fresh phantom. Cast, cool, and allow the phantom to come to room temperature. Mark the top surface of the phantom where FEM nodes are known to lie while the phantom is in the mold. At the same time, make a small sample of the same mixture that will be used in the parallel plate test. 4.4 Finite Element Analysis 41 2. Perform an unconfined parallel plate compression test on the small gelatin sample to acquire the Young's Modulus of the gelatin mixture. 3. Take ultrasound images of the phantom every 2mm using the motion stage. 4. Calibrate the digitizing pen. 5. Place the phantom in the compression apparatus with all surfaces lubricated with ultrasound gel. 6. Move the motion stage until it just touches the phantom, this is achieved by watching the force sensor until it reads a non-zero force. 7. Capture an ultrasound image of the fluid pocket. 8. Capture the position of the surface markers with the digitizing pen. 9. Move the motion stage to deform the phantom, first 3mm, then 6mm, finally 9mm. For each displace-ment, repeat steps 7 and 8. 4.4 Finite Element Analysis A FEM simulation using the incompressible fluid pocket model is created to compare experimental results with simulated results. The first step in re-creating the experiment in simulation is generating a FEM mesh of the phantom. 4.4.1 Mesh Generation v The sequence of ultrasound images of the phantom taken at 2mm increments is used to determine the size, shape, and location of the fluid pocket inside the phantom. The images that show the fluid pocket, such as Figure 4.6(d), are segmented by hand to determine the radius and location of the fluid pocket. The fluid pocket is assumed to have circular contours. Figure 4.7 shows typical-results from a hand segmented ultrasound image. The top, bottom, and vertical lines correspond to the edges of the phantom, while the radius of the pocket is denoted by the arrow extending from the center of the pocket to the edge. Using the 4.4 Finite Element Analysis 42 (a) (b) (c) (d) (e) (f) ( g) Figure 4.6: Sequence of Ultrasound Images used to Re-create Phantom Geometry known dimensions of the phantom, and the relative distances from the ultrasound images, the center point and radius are computed. Figure 4.7: Hand Segmented Ultrasound Image of Fluid Pocket The 3D FEM mesh is created in 4 steps: 1. Create a regular grid with nodes at 5mm spacing using the known dimensions of the phantom mold. For a 65 x 55 x 35mm phantom, the grid size is 14 x 12 x 8, or 1344 nodes. 2. Create a connectivity list to fill the grid with tetrahedral elements. 3. Remove nodes in the center of the block that correspond to the location and approximate radius of the pocket. All tetrahedra that lie inside the fluid pocket are removed from the connectivity list and the list is re-ordered to compensate for the new node numbers. After removal of nodes, the grid size is 1328 nodes. 4.4 Finite Element Analysis 43 4. Use the hand segmented radii and centerpoint locations to move the fluid pocket surface nodes to their appropriate location in the phantom. As the mesh has a 5mm spacing and the ultrasound images have a 2mm spacing, linear interpolation is used to interpolate between the ultrasound images. Figure 4.8 shows the shape and location of the fluid pocket inside the phantom after mesh generation. (a) FEM Mesh Cut away to show Fluid Pocket (b) FEM Surface of Fluid Pocket Figure 4.8: Finite Element Mesh Showing Fluid Pocket Shape Using the node locations and the tetrahedral connectivity list, the stiffness matrix for the phantom is created using the isotropic, linear strain tetrahedral elements. Appendix A gives detailed information on how the stiffness matrix is generated for tetrahedral elements, and how the mesh is created and manipulated. The Young's Modulus for the stiffness matrix is 15.2 kPa, as found from the previous parallel com-pression experiment. Poisson's Ratio for gelatin is taken from the literature. Sciarretta et al performed a Magnetic Resonance Imaging (MRI) based validation of gelatin deformation using the FEM. In their work, the best results were obtained at a Poisson's Ratio of 0.485. 4.4.2 Boundary Condition Changes and Constraints The motion stage imparts a known displacement to the phantom. As such, the boundary condition that 4.5 Results 44 must be applied to the FEM model should also be a displacement. To achieve this, a boundary condition change as described in Section A. 1.5 is necessary on all of the nodes that are in contact with the motion stage. The phantom in the FEM simulation is constrained to slide horizontally along the bottom surface, and to slide vertically between the ultrasound and motion platform. This sufficiently constrains the phantom in the X and Y directions, according to Figure 4.5, and mimics the lubricated surfaces that the phantom is in contact with. A constraint in Z is required, so two nodes at the on the left vertical surface are constrained in Z. These nodes are chosen to be in the center of the phantom, on the bottom surface, so that this constraint has a negligible effect on the deformation of the top surface and fluid pocket. 4.5 Resu l t s 4.5.1 Ultrasound Contour Results The ultrasound images of the phantom are used to compare the shape of the fluid pocket from the phan-tom with the predicted shape from FEM analysis under the various loads imparted during the experiments. In order to overlay the FEM mesh results and ultrasound images, the images need to be scaled, and the fluid pocket has to be registered so that the FEM mesh agrees with the image of the phantom. Figure 4.9 shows the FEM analysis results overlayed on the ultrasound image for the rest case and three deformations. In the rest case, the nodes from the FEM that correspond to those contacting the ultrasound probe were aligned with the top of the ultrasound image, and the nodes that were in contact with the motion stage were aligned with the bright horizontal line across the bottom of the ultrasound image. This bright line is a reflection from the end of the phantom. The FEM nodes were then moved left and right visibly until the nodes matched the contour as well as possible. The interpolation scheme used to display the FEM contour results is nearest-neighbor, which selects a plane of FEM mesh nodes that was purposely placed at approximately the same height as the centerline of the ultrasound probe. The FEM solution predicts the contour of the fluid pocket well. The largest error occurs in the 9mm de-4.5 Results 45 Figure 4.9: Ultrasound vs. FEM Contours 4.5 Results 46 formation case where the top left of the contour was predicted by FEM to move approximately 3.3mm while i the ultrasound contour moved 3.71mm. It is important to note that the accuracy of linear FEM is limited to approximately 1% strain [41]. The three deformation cases shown in this analysis are at approximately 5%, 10% and 15% strain. The assumption of small strain is definitely not true in this analysis. The results from the FEM are not expected to predict the deformations very accurately at these large strains. It is more important that the trends of deformation are consistent. The justification for this accuracy trade-off is the high update rates that linear FEM can achieve with condensation. 4.5.2 Surface Tracked Results Figure 4.10 shows the locations of the phantom markers compared with the locations of the correspond-ing nodes from the FEM simulation. The surface tracked data is fit using a 2 step least squares approach. First, the non-linear least squares method described by An et al. is used to find the best plane that fits the marker positions when the experiment is at rest. The locations of the markers in OptoTrak co-ordinates are then projected onto this plane for comparison. Once projected, a linear least squares fit, as described by [74] is performed to align the marker positions for the various deformation cases with the corresponding FEM node positions. Results in Figure 4.10 show that despite the errors introduced by hand tracking the markers, the FEM simulation agrees well qualitatively with the experimental data in the top plane. The "o" symbols that represent marker positions are approximately 1mm in diameter, which is representative of the estimated 1mm error in the plane from hand tracking. The majority of positions do fall within the "o", indicating that within the error of this experiment, the data agree well with the simulation. Deformation trends are consistent between the FEM simulation and the experimental data. Markers adjacent to the displaced surface exhibit large vertical displacements, while markers adjacent to the fixed surface displace only slightly in the vertical direction. Also, markers in the center of the phantom show only slight expansion horizontally, where markers closer to the edge expand considerably. Figure 4.10 shows a comparison between the surface displacements from the FEM simulation with the corresponding marker displacements. Arrows are used to illustrate the magnitude and direction of each subsequent displacement. As in the previous figure, the markers at the top of the figure are adjacent to the 4.5 Results 47 0.06 r 0.05 0.04 £ , 0 . 0 3 0.02 0.01 Displaced Surface i Fixed Surface • F E M Node Positions O Tracked Markers o o o o 0 O .o o © o o © .o o o o o- O p Q •O D p o © o Q 0 © 0 o 6 G © o-G G O o. o o o o © © © O o-o o D o © c? Q G © G % Q o 0 0 o © o © G o-Q ®° O i o 0 o © \ o °8 0.01 0.02 0.03 x[m] 0.04 0.05 Figure 4.10: Error Plot from Top Surface Tracking 0.06 displaced surface, while the markers near the bottom are adjacent to the fixed surface. Although significant position errors are incurred from hand tracking markers on a soft, deformable object, the trends are consistent between the experimental and simulated data. The markers near the displaced surface move considerably, while the markers near the fixed surface move only slightly. Markers in the center of the phantom compress almost vertically, with very little horizontal expansion, while markers near the edges show considerable horizontal expansion. Figure 4.12 shows the height of the phantom along two sets of markers located above the fluid pocket. The FEM simulation predicts a noticeable bulge in the phantom in this location due to the fluid pockets deformation. Unfortunately, the experimental height data is has considerable error from hand placing the sharp tip of the digitizing pen as close as possible to the phantom markers. The tip penetrated and hovered above the marker with an estimated 1mm error. This error causes the height readings to be being very 4.6 Summary 48 0.05 0.04 £.0.03 N 0.02 0.01 r Displaced Surface M i l ) I \ \ I i i l I \ I \ \ s 0.05 0.04 • 0.03 0.02 0.01 Displaced Surface / / I! / / ! / / / / i / / / t \ \ \ \ \ \ 0.02 0.04 x[m] 0.06 0.02 0.04 x[m] 0.06 Displacements from Experiments Displacements from FEM Simulation Figure 4.11: Top Surface Motion, Tracked and FEM scattered, and as a result this bulge is not visible in the experimental data. It is, however, promising that the average height of the phantom markers is similar to the average height of the FEM simulation. 4.6 S u m m a r y In this chapter, the experimental validation of the fluid pocket model was presented. Technical aspects of the experiment were described, including the use and calibration of an OptoTrak® based digitizing pen, the construction of a fluid filled phantom for ultrasound viewing, and the experimental measurement of Young's Modulus with an unconfined parallel plate compression test. The method for modeling the experiment with a FEM analysis was described, and experimental results were compared with simulated results. The simulated deformation shows similar trends to the experimental data. 4.6 Summary 49 | Displaced Surface Fluid Pocket • • • • • • < • Fixed Surface Marker Line Shown 0.04 0.039 0.038 ^0.037 0.036 0.035 0.034 0 0.04 0.039 0.038 * No Load • 3mm Disp * 6mm Disp 9mm Disp £-0.037 >> « o * 0.036 • • • • • • 0.035 * * * * * * 0.034 0.06 0.02 0.04 0.06 0 0.02 0.04 x [m] x [m] Top Surface Height from Experiments Top Surface Height from F E M Simulation (a) Location (b) Displaced Surface | Fluid Pocket Fixed Surface Marker Line Shown 0.04 0.04 0.039 0.039 0.038 0.038 i.0.037 o * •> * • • •i-0.037 * •> 0.036 0.035 0.034 • • D • * * * * 0.036 0.035 0.034 o o * * * No Load a 3mm Disp <• 6mm Disp 9mm Disp 0.06 0.02 0.04 0.06 0 0.02 0.04 x [m] x [m] Top Surface Height from Experiments Top Surface Height from FEM Simulation (c) Location (d) Figure 4.12: Vertical Expansion Plot, Phantom Markers vs. FEM Simulation Chapter 5 Haptic Simulation This chapter first describes how FEM simulations can be accelerated to haptic update rates using matrix condensation. The step by step haptic algorithm is given to illustrate how the fluid pocket model is imple-mented in real time. A 2D implementation of an incompressible fluid structure based on the haptic needle insertion simulation from [34] is presented. Issues concerning 3D simulation are discussed. 5.1 Haptic Algorithm The haptic algorithm for creating a fluid pocket simulation in real time has several components, the most important of which is creating a FEM solution of a given geometry that can be solved at real-time haptic rates. The pre-processing step that allows for real-time simulation is condensation, which is described in detail in Appendix A. 1.6. Without going into too much detail, condensation allows for the complete, linear FEM system to be reduced to a smaller subset of nodes with the knowledge that many of the nodes do not have any load applied on them. As a result, these "free" nodes can be eliminated from the system by either Gaussian Elimination or Selective Matrix Vector Multiplication (SMVM), depending whether the implicit or explicit system is being solved. When a FEM system is reduced with condensation, the notation is changed from Ku = / to Kwy = x where x and y are the composite displacement and load vectors. This notation change is done because the nodes in contact with the surgical tool often require displacement and force boundary conditions to co-exist 5 , 50 5.1 Haptic Algorithm 51 in the load vector, y. Kw is referred to as the "Working" matrix because it is the matrix that is used for real-time solution of the FEM system. 5.1.1 Real-time Haptic Algorithm for Incompressible Fluid In general, haptic simulators have two separate update loops for graphics and haptics because the update rate requirements for the each differ by more than an order of magnitude. Because required haptic update rates are of the order 300Hz, an absolute minimum of computations should be performed in the haptic loop. The haptic algorithm for simulating an incompressible fluid filled structure is illustrated in Figure 5.1. Initially, Kw is condensed to include only the nodes on the fluid boundary. As the simulation runs, nodes that are in contact with the surgical tool should be added to Kw. Tool Position Compute x = K~ly with P and tool position Compute V Error < Tolerance Haptic Force YES NO Update = P + KpError Figure 5.1: Haptic Algorithm Data Flow The algorithm in the flowchart is presented with the assumptions that the surgical tool is in contact with the elastic body, and the appropriate force and/or displacements have been added to the load vector. It is also assumed that the proportional gain, Kp, has been chosen for the elastic object's pressure/volume response as described in Section 3.3.2. Because Kw is condensed to use only the "active" nodes, the haptic loop described above can be computed at real-time haptic rates. 5.2 Interactive Haptic Simulation in 2D 52 5.2 Interactive Haptic Simulation in 2D A 2D haptic simulation of an incompressible fluid filled pocket enclosed in elastic media was imple-mented using the 3DOF Twin Pantograph, using the software framework provided by the needle insertion simulation found in [33,34]. This addition allows the user to feel the difference between the needle response to a compressible cavity and the needle response to a cavity filled with incompressible fluid. 5.2.1 3DOF Twin Pantograph The 3DOF twin pantograph uses a redundant parallel linkage design to allow for planar interaction (x and y-axis translation, and rotation around z-axis) as shown in Figure 5.2. Each link is actuated by a DC motor in a direct drive configuration, and each motor is monitored by an optical 12 bit shaft encoder. The motors are powered by a 4 channel current driver. Haptic interaction with the twin pantograph uses two independent PC's that communicate over a UDP socket connection. One PC, that runs VxWorks®, contains a DAC board that reads the pantograph's en-coders and outputs control voltages to the current driver. This PC runs all of the pantograph control, dy-namics, and simulation software, which is written using the Matlab's Realtime Workshop and the Tornado® integrated development environment. The second PC runs Windows®, and is responsible for displaying the graphical output of the virtual environment and compiling the real-time code for the VxWorks® PC. The Figure 5.2: 3DOF Twin Pantograph Haptic Interface 5.2 Interactive Haptic Simulation in 2D 53 real-time software is run with a hard real-time update rate of 512Hz [21,36,73]. 5.2.2 Modifications Required for Fluid Cavity Simulation The original needle insertion simulation from [34] uses a condensed FEM solution in combination with ex-perimentally based needle insertion force profiling for a solid elastic body. To implement the incompressible fluid simulation the following software components were added: • Precomputed pressure normal directions and the undeformed cavity volume were added to the initial-ization stage. • The FEM system was initially condensed to include the 24 fluid boundary nodes. As the needle simulation progresses, needle nodes are added to the system. • The haptic algorithm as outlined in Section 5.1.1 was implemented in the inner haptic loop. The optimized gain from Section 3.3.2 was used because the geometry and stiffness is identical. The maximum number of iterations was fixed at 2 to ensure a finite solution time. This is done as a precaution to avoid cut-off by the hard real-time operating system. • The needle puncture and friction force thresholds are raised to increase the compression force that is transmitted to the cavity so that deformations are more visible. 5.2.3 Results from 2D Interactive Haptic Simulation Figure 5.3 shows the visual output from the interactive haptic simulation of an incompressible fluid inclu-sion. Figure 5.3(a) shows the incompressible fluid filled structure before the needle is penetrated into the elastic material. Figures 5.3(b) and 5.3(c) show approximately the same needle insertion depth to illustrate the difference in cavity shape. Figure 5.3(b) illustrates the response to the needle when the cavity is empty and there is no pressure difference across the elastic body. The cavity exhibits an indentation from needle compressive force, but there is no outward bulging of cavity. Figure 5.3(c) shows the response to the needle when the cavity is filled with incompressible fluid. Note how the cavity bulges to maintain it's volume. It is important to note that the elastic object is constrained on three sides, which causes the elastic material to 5.2 Interactive Haptic Simulation in 2D 54 (c) Incompressible Fluid Pocket Figure 5.3: Haptic Simulation of an Incompressible Fluid inclusion in a 2D Environment 5.3 3D Performance 55 bulge slightly to the left to accommodate the bulging of the pressurized fluid cavity. This bulging will have an impact on the perceived stiffness of the object. A quantitive comparison between the incompressible fluid filled and compressible fluid filled objects has not been done. However, there is a noticeable difference in the perceived stiffness between the two configurations. Because the incompressible fluid pocket resists volumetric strain by bulging, the needle displaces a smaller distance before penetrating the elastic boundary than it would in the compressible fluid filled pocket case. This causes the incompressible fluid filled object to feel stiffer, and the compressible filled object to feel more compliant. As previously mentioned, the fluid pressure algorithm was set to iterate 3 times to ensure a fixed solution time. This is, in fact, a conservative number because the analysis in Section 3.3.2 showed that convergence to a 1% tolerance is usually achieved in only 1 iteration. To investigate how many iterations would be possible with this size of a system, the number of iterations was increased until it was too large for the hard real-time cut-off. 7 iterations of a 24 node fluid pocket with 3 or 4 nodes on the needle was found to exceed the cut-off. Because it was possible to do 6 iterations, and generally only 1 is necessary to achieve a 1% tolerance, a larger fluid system could be simulated. 5.3 3D Performance The complexity of the method for simulating fluid structures in an elastic body presented in this work scales with 0(n2) on the number of nodes on the fluid-elastic boundary. That is: 1. The system product x = K~1 y requires 3M 2 + (3n — 1) (3n) FLOPs. 2. The volume calculation in 3D requires approximately 128n FLOPs. For comparison, the method described by De and Srinivasan [30] required approximately 2M FLOPs for a thin-walled system with 41 nodes in 3D. With the linear FEM approach presented in this work, a similar sized pocket would require approximately 90k FLOPs, assuming 2 iterations were necessary for convergence, thus reducing the processing requirement by almost 2 orders of magnitude. Also, because the system can be condensed to include only active nodes, the linear FEM approach can accommodate for a 5.4 Summary 56 large volume of elastic solid to be simulated in addition to the fluid structure with a minimal increase in computational cost. 5.3.1 Timing of an Off-line 3D Simulation A real-time haptic simulation that runs at 512 Hz has been created using the fluid model presented in this work. However, to illustrate that the method will be able to simulate a relatively large fluid inclusion in 3D using present hardware, an off-line simulation is timed to assess the speed. Using the same 3D phantom mesh created in Section 4.4, an off-line fluid simulation was implemented using C++. The full phantom mesh consists of 1328 nodes, which is clearly too large for simulation in real-time on present hardware. Using condensation, the mesh is reduced to include the 80 nodes that are on the closed fluid boundary and 2 nodes that are on the exterior surface. The reduced system, consisting of 82 nodes, is used as an example of a FEM system that would include fluid nodes and nodes from a surgical tool, like a needle, in the haptic update loop of a surgical simulator. Table 5.1 lists approximate update rates for the 82 node, 3D FEM simulation using a Pentium 4 2.4 GHz PC with 512Mb ram. The times correspond to a simulation that required 2 iterations to achieve desired tol-erance with single precision floating point variables. 2 iterations of a 82 node system requires approximately 370k FLOPs. The algorithm presented in Section 5.1.1 is used to generate these times. Compiler Approx. Update Rate (Hz) Approx. CPU Performance (FLOPs per sec) Microsoft C++ 6.0, optimized 1060 0.36 Giga Intel C++ 7.0, optimized 3120 1.06 Giga Table 5.1: CPU Times for an 80 Node Fluid Structure In both cases, the maximum optimization compiler flag is used to achieve the highest possible perfor-mance. The C++ implementation is also written with many common optimization practices such as passing arrays by reference and minimizing array evaluations with pointer arithmetic. 5.4 Summary This chapter first described matrix condensation, then discussed the use of the hydrostatic pressure based 5.4 Summary 57 fluid pocket model with condensed FEM matrices for real-time haptic simulations. A 2D interactive haptic simulation that runs at 512 Hz using the 3DOF twin pantograph was described, and the performance of an off-line 3D a simulation was discussed. Chapter 6 Conclusions and Recommendations 6.1 Contributions of this Thesis This thesis has presented contributions to the development of new and more life-like haptic medical simulations by introducing a method for modeling fluid structures enclosed in an elastic body with the FEM that is fast enough for haptic display. To ensure that this method predicts deformation in a realistic manner, an experimental validation was performed. A 2D fluid pocket simulation that runs at 512 Hz was implemented and an offline 3D simulation was timed to illustrate that the method is viable for real-time simulation. 6.1.1 Fluid Pocket Model The fluid pocket model developed in this thesis is based on physical principles of fluid and solid me-chanics. Hydrostatic fluid pressure is applied as a force boundary condition to a linear strain FEM solution of an elastic body to simulate the presence of fluid. A fast numerical method, which employs a proportional feedback control algorithm, is tuned with system identification analysis and is used to enforce a volume to pressure relationship inside the fluid inclusion. In this thesis, the fluid was simulated as incompressible for the case of fluid filled organs like the bladder, or various small anatomical structures like cysts or glands. Although only incompressible fluid was simulated in the thesis, the methods can accommodate any pressure 58 6.1 Contributions of this Thesis 59 to volume relationship. The fluid pocket model described is based on a static analysis of fluid, and uses a quasi-static implemen-tation in real-time. Only a static force balance is done at every iteration, so the acceleration and velocity of fluid and elastic material are not considered. This method could be readily used, however, to simulate dynamic phenomena such as the pulsing of arteries. As an approximation, the pulsing of an artery could be modeled using a fluid pressure that changes in time and magnitude according to a humans blood pressure and heart rate. Although the model was developed and implemented using the FEM, the method employs a general principle, and could work with other deformable methods such as SMS, BEM, or any other method that can accept a force boundary condition from hydrostatic fluid pressure. 6.1.2 Validation of Fluid Structure Model using 3D Fluid Filled Phantom To validate that the hydrostatic pressure based fluid pocket model predicts the deformation of an incom-pressible fluid in a realistic manner, a fluid filled phantom was constructed and deformed while recording displacements so the data could be compared with a FEM solution of the same phantom. The surface of the phantom was tracked by hand with a digitizing pen, and the contour of the fluid pocket was imaged using ultrasound. A compression test was performed to experimentally compute the elastic properties of the phantom, and ultrasound images were used to create a 3D FEM mesh of the phantom. The ultrasound imaged results agreed well with the FEM simulation. Results from the surface tracking and contour imaging are compared with the corresponding deforma-tions from the FEM simulation. The displacements from the top plane agree well between the simulation and experiment despite the error introduced by hand tracking. However, it is difficult to correlate the data displaying the height of the phantom due to varying heights that the pen tip penetrates into the phantom. Although experimental data agrees well with the FEM simulation, the validation work done presented in this thesis should be regarded as preliminary. Validation with more complicated structures and more complicated loadings are suggested to further validate the method. 6.2 Recommendations for Future Work 60 6.1.3 Real-time Haptic Simulation of Fluid Pockets The hydrostatic fluid pressure model has been added to the needle insertion simulation described in [34]. This implementation is computed in 2D, and simulates the deformation of an incompressible fluid pocket at an update rate of 512 Hz using a 3DOF haptic interface. To achieve real-time rates, the FEM stiffness matrix is pre-inverted and condensed to include only working nodes that are either on the fluid-solid boundary or on the needle. With regards to 3D real-time simulation, offline CPU times for an 80 node, 3D fluid pocket were recorded to illustrate that this method will be able to compute for a relatively large 3D fluid body on cur-rent hardware. Using a 2.4GHz Pentium 4 with the Intel C++ compiler with maximum optimization, an update rate of approximately 3000Hz was recorded. It is important to note that the 3000 Hz update rate included only the FEM simulation. Haptic interface control, graphics and other related computations are not considered. 6.2 Recommendations for Future Work There are several opportunities for the work presented in this thesis to be continued and improved. Some of the possibilities include: 1. The surface tracking done in this thesis was done by hand with a digitizing pen. Because the object being tracked is a soft, deformable object, and the pen has a sharp tip for calibration, it is very difficult to accurately touch the surface marker with the pen. The use of rigidly attached IREDs could increase the accuracy of the surface tracking if a way of fastening IREDs to the deformable surface is found. Although rigidly fixing the IREDs to the surface would remove hand tracking error, it would be possible to collect a much more accurate and comprehensive set of deformation data if experiments were imaged using MRI or CT. The surface of the phantom could be outfitted with small markers, and the fluid pocket could be filled with the same contrast agent used in [18]. A set of data collected in this manner would have almost no error introduced by the user, and the 3D deformation of the entire surface and shape of fluid pocket could be collected quickly and accurately. 6.2 Recommendations for Future Work 61 2. As a further extension to the simulation presented by [34], a needle insertion simulation in 3D could be created with a variety of fluid inclusions. Easily deformable, low pressure veins, stiff high pressure arteries, glands, and cysts could all be added to increase the realism and difficulty of needle placement. 3. Application to prostate brachytherapy. The bladder is filled with a fluid contrast agent during prostate brachytherapy so that it is visible in CT images [18]. The incompressible fluid model presented in this thesis could be used to model the bladder and surrounding tissues if a haptic simulation of prostate brachytherapy is to be implemented. 4. Applications to vascular related research. Using a phantom similar to that described in [44], this method for simulating fluid structures could be validated for long cylindrical structures with appro-priate pressure/volume relationships taken from literature regarding the vascular system. 5. Application to surgical planning. Many surgical procedures are planned using medical images and deformable modeling [31]. Any planning procedure that involves fluid filled organs could be improved by the addition of fluid-elastic interaction as described in this thesis. 6. Fluid pressure prediction. It would be very interesting to outfit the fluid filled section of a phantom with a pressure transducer so that the pressure of the fluid can be monitored and compared to the value computed in simulation. 7. Psychophysics experiments. The presence of a fluid pocket affects the deformation of the elastic body, especially in close proximity to the fluid pocket. It would be interesting to do a variety of user tests with simple geometric models to see if users can tell the difference between bodies that have, for example, an incompressible fluid pocket, an empty pocket, or a softer elastic material inside. Bibliography [1] P. Abolmaesumi. Image-Guided Robot-Assisted Diagnostic Ultrasound. PhD thesis, University of British Columbia, 2002. [2] M. Agus, A. Giachetti, E. Gobbetti, G. Zanetti, N. W. John, and R. J. Stone. Mastoidectomy Simulation with Combined Visual and Haptic Feedback. In Medicine Meets Virtual Reality, pages 17-23, 2002. [3] M. Agus, A. Giachetti, E. Gobbetti, G. Zanetti, and A. Zorcol. Real-time Haptic and visual Simulation of Bone Dissection. In IEEE Virtual Reality, pages 209-216, 2002. [4] N. Akkas, H. U. Akay, and C. Yilmaz. Applicability of General-purpose Finite Element Programs in Solid-fluid Interaction Problems. In Computers and Structures, volume 10, pages 773-783, 1979. [5] C. H. An, C. G. Atkeson, and J. M. Hollerbach. "Model-Based Control of a Robot Manipulator ". MIT Press, 1988. [6] N. Archip and R. Rohling. Fast Reconstruction of Volumetric Models of Anatomical Structures. In IEEE: Engineering in Medicine and Biology, 2003. [7] R. Balaniuk and K. Salisbury. Dynamic Simulation of Deformable Objects using the Long Elements Method. In 10th Symposium On Haptic Interfaces for Virtual Environments and Teleoperator Systems, 2002. [8] D. Baraff and A. Witkin. Large Steps in Cloth Simulation. In Proceedings of Computer Graphics, pages 44-54, 1998. [9] C. Basdogan, C. H. Ho, and M. A. Srinivasan. Simulation of Tissue Cutting and Bleeding for Laparo-scopic Surgery Using Auxiliary Surfaces. In Medicine Meets Virtual Reality, pages 38^14, 1999. [10] C. Basdogan, C. H. Ho, and M. A. Srinivasan. Virtual Environments for Medical Training: Graphical and Haptic Simulation of Laparoscopic Common Bile Duct Exploration. IEEE/ASME Transactions on Mechatronics, 6(3):269-285, 2001. [11] J. K. Bathe. "Finite Element Procedures ". Prentice Hall Inc., 1996. [12] J.K. Bathe and W. F. Hahn. On Transient Analysis of Fluid-Structure Analysis. Computers And Structures, 10:383-391, 1979. [13] J.K. Bathe, H. Shang, and S. Ji. Finite Element Analysis of Fluid Flows Fully Coupled with Structural Interactions. Computers And Structures, 72:1-16, 1999. 62 BIBLIOGRAPHY 63 [14] T. Belytschko. Fluid-structure interaction. Computers and Structures, 12:459^-69, 1980. [15] J. Berkley, S. Weghosrt, H. Gladstone, G. Raugi, D. Berg, and M. Ganter. Banded Matrix Approach to Finite Element Modeling for Soft Tissue Simulation. In Virtual Reality: Research, Development and Application, volume 4, pages 203-212, 1999. [16] M. Bro-Nielsen and S. Cotin. Real-time Volumetric Deformable Models for Surgery Simulation using Finite Elements and Condensation. Proceedings of Eurographics, 15(3):57-66, 2002. [17] I. Brouwer, K. E. MacLean, A. J. Hodgson, A.G. Nagy, K.A. Qayumi, and L.M. Rucker. Price-Quality Trade-offs In Haptic Interfaces for Simulation of Laparoscopic Surgery. In Poster Presented at Medicine Meets Virtual Reality, 2003. [18] W.M.Butler. Review of Modern Prostate Brachytherapy. In Proceedings of the 22nd Annual EMBS International Conference, pages 748-752, 2000. [19] Y. Calayir and A. Donangun. Static and Dynamic analysis of Fluid and Fluid-structure Systems by the Lagrangian Method. Computers and Structures, 49(4):625-632, 1993. [20] K. Chinsei and K. Miller. Compression of Swine Brain Tissue - Experiment in Vitro. Journal of Mechanical Engineering Laboratory, 5(4): 106—115, 1996. [21] D. Constantinescu, I. Chau, S.P. DiMaio, L. Filipozzi, S.E. Salcudean, and F. Ghassemi. Haptic Ren-dering of Planar Rigid-Body Motion using a Redundant Parallel Mechanism. In IEEE International Conference on Robotics and Automation, 2000. [22] R. D. Cook, Malkus D. S., and Pelsha M. E. Concepts and Applications of Finite Element Analysis. John Wiley and Sons, 1989. [23] Ultrasonix Medical Corp. World Wide Web Page. [24] S. Cotin and H. Delingette. Real-time Surgery Simulation with Haptic Feedback using Finite Ele-ments. In IEEE International Conference on Robotics and Automation, 1998. [25] S. Cotin, H. Delingette, and N. Ayache. Real-Time Elastic Deformations of Soft Tissues for Surgery Simulation. IEEE Transactions on Visualization and Computer Graphics, 5(l):52-73, 1999. [26] Stephane Cotin, Herve Delingette, and Nicholas Ayache. A Hybrid Elastic Model for Real-time Cut-ting, Deformations, and Force Feedback for Surgery Training and Simulation. The Visual Computer, 16(7):437^152, 2000. [27] D. Crochemore, J Louchet, and X. Provot. Evolutionary Identification of Cloth Animation Model. In ln Workshop on Computer Animation and Simulation (Eurographics 95), pages 44—54, 1995. [28] D. d'Aulignac, R. Balaniuk, and C. Laugier. A Haptic Interface for a Virtual Examination of the Human Thigh. In IEEE International Conference on Robotics and Automation, 2000. [29] D. d'Aulignac, M. C. Cavusoglu, and C. Laugier. Modeling the Dynamics of the Human Thigh for a Realistic Echographic Simulator with Force Feedback. In Medical Image Computing and Computer-Assisted Intervention, 1999. BIBLIOGRAPHY 64 [30] 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, pages 94-99, 1999. [31] H. Delingette. Towards Realistic Soft Tissue Modeling in Medical Simulations. In Proceedings of the IEEE: Special Issue on Surgery Simulation, pages 512-523, 1998. [32] O. Deussen, L. Kobbelt, and P. Tucke. Using Simulated Annealing to Obtain good Approximations of Deformable Bodies. In In Proc. Eurographics Workshop on Animation and Simulation, 1995. [33] S. P. DiMaio and S. E. Salcudean. "Needle Insertion Modelling and Simulation". In Proceedings of the IEEE International Conference on Robotics and Automation, pages 2098—2105, May 2002. [34] S. P. DiMaio and S. E. Salcudean. Simulated Interactive Needle Insertion. In 10th Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems, 2002. [35] S.P. DiMaio. Physical Modeling of Tissue for Needle Steering. PhD thesis, University of British Columbia, 2003. [36] S.P. DiMaio, S.E. Salcudean, and M.R. Sirosupour. Haptic Interaction with a Planar Environment. In 9th Symposium on Haptic Interfaces for Virtual Environments and Teleoperator Systems, pages 1223-1230. International Mechanical Engineering Congress and Exposition (ASME Winter Annual Meeting), 2000. [37] A. Donangun, A. Durmus, and Y. Ayvaz. Static and Dynamic Analysis of Rectangular Tanks by using the Lagrangian Fluid Finite Element. Computers and Structures, 59(3):547-552, 1996. [38] Larsen E.S. and McAllister D. Fast Matrix Multiplies using Graphics Hardware. In Proceedings ofSC 2001. ACM. [39] Y. C. Fung. "Biomechanics: Mechanical Properties of Living Tissues ". Springer-Verlag, second edition, 1993. [40] S. F. F. Gibson. 3D ChainMail: A Fast Algorithm for Deforming Volumetric Objects. Tech-nical Report TR-96-22, MERL - A MITSUBISHI ELECTRIC RESEARCH LABORATORY,, 1996. [41] S. F. F. Gibson and B. Mirtich. A Survey of Deformable Modeling in Computer Graphics. Technical Report TR-97-19, MERL - A MITSUBISHI ELECTRIC RESEARCH LABORATORY,, 1997. [42] H. B. Gladstone, G. J. Raugi, D. Berg, J. Berkley, S. Weghorst, and M. Ganter. Virtual Reality for Dermatologic Surgery: Virtually a Reality in the 21st Century. Journal of the American Academy of Dermatology, 42(1), 2000. [43] S. Gopalakrishnan and G. Devi. A Lagrangian Quadratic Triagnular Fluid Finite Element for Fluid-Structure Interaction Problems. International Journal of Computational Engineering Science, 1(2):257-272, 2000. [44] J. Guerrero, S. E. Salcudean, J. A. McEwen, B. A. Basri, and S. Nicolaou. Mesaurement-Based Deep Venous Thrombosis Screening System. In Proceedings of Medical Image Computing and Computer Aided Interventions, pages 214—221, 2003. BIBLIOGRAPHY 65 [45] K. Haberman, S. D. Pathak, and Y. Kim. Effects of Video Digitization in Pubic Arch Interfer-ence Assessment for Prostate Brachytherapy. In IEEE Transactions on Information Technology in BioMedicine, volume 7, pages 8-15, 2003. [46] T.J. Hall, M. F. I. Bilgen, and T. A. Krouskop. Phantom Materials for Elastography. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency control, 44(6): 1355-1365, 1997. [47] R. C. Hibbeler. Statics and Mechanics of Materials. Macmillan Publishing Company, 1993. [48] T. J. R. Hughes, W. K. Liu, and T. K. Zimmermann. Lagrangian-Eulerian Finite Element Formula-tion for Incompressible Viscous Flows. Computer Methods in Applied Mechanics and Engineering, 29:329-349, 1981. [49] Computer Motion Inc. World Wide Web Page. [50] Northern Digital Inc. World Wide Web Page. [51] Sensable Technologies Inc. World Wide Web Page. [52] D. L. James and D. K. Pai. A Unified Treatment of Elastostatic Contact Simulation for Real Time Haptics. Haptics-E, 2(1): 1-13, 2001. [53] D.L. James and D. Pai. ARTDEFO: Accurate Real Time Deformable Objects. In Proceedings of Siggraph, 1999. [54] M. Kass and G. Miller. Rapid, Stable Fluid Dynamics for Computer Graphics. Computer Graphics, 24(4):49-57, 1990. [55] M. Kass, A. Witkin, and D. Terzopoulos. Snakes: Active Contour Models. International Journal of Computer Vision, l(4):321-3324, 1988. [56] J. Lang. Deformable Model Acquisition and Validation. PhD thesis, University of British Columbia, 2001. [57] C. Laugier, C. Mendoza, and K. Sundaraj. Faithfull Haptic Feedback in Medical Simulators. In 8th International Symposium on Experimental Robotics, 2002. [58] X. L i and J. Moshell. Modeling Soil: Realtime Dynamic Models for Soil Slippage and Manipulation. In Computer Graphics, 1993. [59] D. Malkis. A finite element displacement model valid for any value of the compressibility. In Interna-tional Journal of Solids and Structures, volume 12, pages 731-738, 1976. [60] B. Mirtich. Fast and Accurate Computation of Polyhedral Mass Properties. Journal of Graphics Tools, 1(2), 1996. [61] A. Moravanzky. Dense Matrix Algebra on the GPU. In ShaderX 2, 2003. [62] B. R. Munson, D. F. Young, and T. H. Okiishi. "Fundamentals of Fluid Mechanics ". John Wiley and Sons, third edition, 1998. BIBLIOGRAPHY 66 [63] J. F. O'Brien and J. K. Hodgins. Dynamic Simulation of Splashing Fluids. In Computer Animation, 1995. [64] A. Pentland and J. Williams. Good Vibrations: Modal Dynamics for Graphics and Animation. In Proceedings of SIGGRAPH, pages 215-222, 1989. [65] H. H. T. Pian and S. W. Lee. Notes on Finite Elements for Nearly Incompressible Materials. AIAA Journal, 14(6):824-826, 1976. [66] K. Proudfoot, W.R. Mark, S. Tzvetkov, and Hanrahan P. A Real-Time Procedural Shading System for Programmable Graphics Hardware. In Proceedings of SIGGRAPH. ACM, 2001. [67] T.J. Purcell, I. Buck, R.M. William, and P. Hanrahan. Ray Tracing with Programable Graphics Hard-ware. In Proceedings of SIGGRAPH. ACM, 2002. [68] D.W. Rickey, P. A. Picot, D. A. Christopher, and A. Fenster. A Wall-Less Vessel Phantom for Doppler Ultrasound Studies. Ultrasound in Medicine and Biology, 21(9): 1163-1176, 1995. [69] R. Rohling, P. Munger, J. M. Hollerbach, and T. Peters. Comparison of Relative Accuracy Between a Mechanical and an Optical Position Tracker for Image-Guided Neurosurgery. Journal of Image Guided Surgery, 1:30-34, 1995. [70] S.E. Salcudean and T.D. Vlaar. On the Emulation of Stiff Walls with Static Friction with a Magnetically Levitated Input/Output Device. ASME Journal of Dynamics, Measurement, and Control, 119:127-132, 1997. [71] J. Sciarretta, A. Samani, J. Bishop, and D. B. Plewes. MR Validation of Soft Tissue Mimicing Phantom Deformation as Modeled by Nonlinear Finite Element Analysis. Medical Physics, 29(l):65-72, 2002. [72] A. Shiraz-Adl and S. C. Shrivastava. Large Deformation Finite Element Treatment of Changes in the Volume of Fluid-Filled Cavities Enclosed in a Structure. Computers and Structures, 34(2):225-230, 1990. [73] M.R. Sirouspour, S.P. DiMaio, S.E. Salcudean, P. Abolmaesumi, and C. Jones. Haptic Interface Con-trol - Design Issues and Experiments with a Planar Device. IEEE Intl Conf on Robotics and Automa-tion, 2000. [74] A. Soroushi. Real-Time Motion Tracking of Mobile Robots via Image Registration. Master's thesis, University of British Columbia, 2003. [75] R. W. Summer, J. F. O'Brien, and J. K. Hodgins. Animating Sand, Mud and Snow. Eurographics Computer Graphics Forum, 18(1), 2002. [76] D. Tang, C. Yan, and J.K Bathe. A 3-D Thin-walled Model with Fluid-structure Interactions for Blood Flow in Carodit Arteries with Symmetric and Asymmetric Stenoses. Computers And Structures, 72:357-377, 1999. [77] D. Terzopoulos, J. Piatt, A. Barr, and K. Fleischer. Elastically Deformable Models. In Proceedings of SIGGRAPH, pages 205-214, 1987. BIBLIOGRAPHY 67 [78] A. Witkin and W. Welch. Fast Animation and Control of Nonrigid Structures. In Proceedings of SIGGRAPH, pages 243-252, 1990. [79] Y.W. Zhang. Direct Surface Extraction from 3D Freehand Ultrasound Images. Master's thesis, Uni-versity of British Columbia, 2003. [80] Y. Zhuang and J. Canny. Real-time Simulation of Physically Realistic Global Deformation. In IEEE Visualization, 1999. [81] Y. Zhuang and J. Canny. Haptic Interaction with Global Deformations. In IEEE International Confer-ence on Robotics and Automation, 2000. Appendix A Finite Element Modeling A . l Linear Elastostatic Methods This appendix describes the procedure used for FEM analysis in this thesis. Linear strain analysis using triangular and tetrahedral elements are discussed. Creation of the global stiffness matrix and issues regarding 3D meshing are described. This appendix covers only the methods used in this thesis. If any further information regarding FEM analysis is desired, readers are refered to [22] and [11]. The strain energy of an elastic body is defined by: where a is the body stress vector and e is the body strain vector. (A-l) can be expressed using only displacements by substituting rj = De (where D is the material elasticity matrix) and e = Bu. B is the Cauchy Strain Tensor: (A-l) d 0 0 0 9 0 B = 0 0 9. 3z (A-2) 3. dy d Tx 0 3 3i 0 dx 0 9z -2-dy J 68 A.l Linear Elastostatic Methods 69 The Cauchy Strain Tensor is often referred to as linear strain because it is assumed that strains are small and therefore higher order terms are neglected. Also, because only tensile and shear strain terms appear in the Cauchy Strain Tensor, it is not invariant to rotation, and must be used accordingly. The after substitution, (A-l) becomes: After discretizing into volumetric elements and interpolating with the natural co-ordinates, the strain tensor B becomes a constant matrix Be. Static equilibrium between deformation energy and external force, / , is achieved when the first variation of E(u)strain vanishes. This occurs when 5£ ,(«)w r a,„ — 0 and corresponds to the potential energy of the system reaching a minimum value [16]. Everything inside the integral is constant, so the integral is reduced to define the stiffness matrix Ke: where Ve is the volume of the element. Be is defined in the following sections for triangular and tetrahedral elements. A . l . l Linear Strain Triangular Element The FEM stiffness matrix for a linear strain triangular finite element is a function of the three node positions [JCI ,yi], [x2,y2], [xi,yi], and the material elasticity matrix, D (A-3) (A-4) Ke = BTeDBeVe (A-5) 2 Figure A . l : Triangular Element A.l Linear Elastostatic Methods 70 The material elasticity matrix for an isotropic material, in 2D, has two possible forms: Plane Strain or Plane Stress. 1. Plane Stress is the case where a flat plate is loaded only within it's plane, and it is unconstrained in the out of plane direction. For this case: D = ( 1 - v 2 ) 1 v 0 v 1 0 0 0 l - v (A-6) 2. Plane Strain is the case where the planar object is constrained in the out of plane direction. For this case: D = ( l + v ) ( l - 2 v ) l - v V 0 V l - v 0 0 0 l -2v 2 1 X = _ y _ Be is defined by the Natural Co-ordinates, L\, L2 and L3 from: 1 1 1 x\ x2 x3 Vl V2 V3 which, as long as nodes are not co-linear, can be inverted to: a\ b\ c\ a2 b2 C2 a 3 b-i c-i L\ L2 = . L 3 Lx L2 L3 1 X y (A-7) (A-8) (A-9) A. I Linear Elastostatic Methods 71 then Be is: b{ 0 b2 0 b3 0 0 c i 0 c2 0 C3 cj bi c2 b2 c 3 c73 (A-10) For the 2D triangular element case, Be is 3 x 6, D is 3 x 3, and is 6 x 3. So the product BTeDB^/e .yeilds a 6 x 6 FEM stiffness matrix Ke for that element. A.1.2 L inear St ra in Tetrahedral Element The FEM stiffness matrix for a linear strain tetrahedral finite element is a function of the four node positions [xi,vi,zi], [x2,y2,z2], [x3,y3,z3], [x4,y4,Z4] and the material elasticity matrix, D Figure A.2: Tetrahedral Element The material elasticity matrix for an isotropic material, in 3D is: + 2fi X X 0 0 0 X X + 2fi X 0 0 0 X X X + 2[i 0 0 0 0 0 0 M 0 0 0 0 0 0 M 0 0 0 0 0 0 M (A-l l) where X = ( 1 + v ^ _ 2 v ) and fi= ^ A A Linear Elastostatic Methods 72 Be is defined by the Natural Co-ordinates, Lx,L2, L3 and L4 from: 1 1 1 1 1 L\ X X4 L2 y yi yi y^ yn L3 z Z3 Z4 LA (A-12) which is inverted to: Li a\ 6l c\ 1 L2 a2 b2 c2 X L3 «3 C3 <*3 v L 4 «4 b4 c 4 rf4 z (A-13) And Be is: 61 0 0 b2 0 0 h 0 0 4^ 0 0 0 c\ 0 0 Q 0 0 c3 0 0 CA 0 0 0 dx 0 0 d2 0 0 0 0 ti4 c\ 61 0 c2 b2 0 C3 *>3 0 c 4 0 0 rfi c\ 0 d2 Q 0 d3 C3 0 J4 c 4 dl 0 bx <*2 0 ^ 2 d3 0 &3 0 b4 (A-14) For the 2D trianglular element case, Be is 6 x 12, D is 6 x 6, and is 12 x 6. So the product BTeDBeVe yeilds a 12 x 12 FEM stiffness matrix Ke for that element. A.1.3 Stiffness Matrix Generation by Node Numbering The formula Ke = BTe,DB^/e gives the local stiffness matrix for a single finite element. By generating a local stiffness for each element in a system of nodes, the global stiffness matrix can be created. Figure A.3 illustrates a simple 2D elastic system consisting of 4 nodes. For illustration, each node has only one DOE The linear system Ku = f is created by: A. 1 Linear Elastostatic Methods 73 a\ 0 a2 a 3 0 0 0 0 0:4 0 a5 as a7 0 a% a<) bi bi bi 0 bA b5 be 0 bi bs b9 0 0 0 0 0 (A-15) (A-16) Kghbal =Ka + Kb (A-17) The resulting linear system, with node numbers according to figure A.3, that describes the deformation of the elastic body is: a\ + b\ 2>2 #2 + hi a-i b\ b5 b6 0 a^-j-bj fc8 a5+bc> a6 a-i 0 a% ag A.1.4 Constraints U\ fx U2 fl W3 h W4 (A-18) The FEM formulation used in this work yields a set of linear equations that describes the deformation of A. 1 Linear Elastostatic Methods 74 the interconnected system of elements and nodes. (A-19) is the stiffness matrix for a simple elastic system made of two ID truss elements, shown by Figure A.4. U\ U2 «3 L\ L2 Figure A.4: Simple ID Truss System K k\ -k\ 0 -Jti ki+k2 -k2 (A"19) 0 -k2 k2 were is the stiffness entry for each truss, respectively, and kj = It is clear that K is singular because column two is a linear combination of columns one and three. As result, there is no unique solution to the system unless one of the nodes is constrained. Any node can be constrained by replacing the corresponding column and row in the stiffness matrix with the the same column and row from the identity matrix, and replacing the corresponding force and displacement entries with zero. It is necessary to constrain a FEM system in each direction to ensure that a unique solution exists [22]. (A-20) shows a suitably constrained system were the first node is fixed. 1 0 0 0 k\+k2 -k2 0 -k2 h 0 0 u2 = h "3 h (A-20) A.1.5 Boundary Condition Changes by Low Rank Updates It is often necessary to have both displacement and force entries in the load vector of a FEM system. A.l Linear Elastostatic Methods 75 Any linear system can be processed using a low-rank matrix update to rewrite: as, A • ^1,3 B = C ^3,1 • • • ^3,3 A • • #1,3 P = C • • • ^3,3 a P . Y . (A-21) a B . Y . (A-22) This low-rank update can be used on FEM stiffness matrices to exchange displacement boundary condi-tions for force boundary conditions, or vice versa, by: (A-23) where K~1 is the original inverted stiffness matrix, K~x is an updated composite matrix, c, and r, are the ith row and column of K~x with their ith entry set to p, — 1 and pi is the pivot value, K(ity After boundary condition changes have been performed on K~x the FEM notation changes to x = K~ly where x replaces u and y replaces / because what used to be homogeneous displacement and force vectors now contain a mix of displacement and force entries. A.1.6 Matrix Condensation Matrix Condensation is an optimization technique that reduces the number of DOF of the FEM solution by substitution using Gaussian Elimination [22]. The linear system K u = f can be re-ordered as: Kss K-si Us L Hi (A-24) A.l Linear Elastostatic Methods 76 where A^w^anci/^ are stiffness matrix entries, nodal displacements and nodal forces for all of the nodes on the surface respectively. From (A-24), a condensed system containing only the surface nodes can be produced: KsMs = C (A-25) Where: = Kss-KsiKrlKis (A-26) £ = f^-K^f. (A-27) The results from the condensed system are exactly the same as those from the full system as long as force is only applied on the surface nodes, hence = / \ This technique can be used optimize FEM solutions for real-time, as typically very few nodes are ac-tive, and it is possible to condense the system for an arbitrary set of nodes. The only drawback using the condensed implicit system K*^ = is that any change in the required nodes or mesh topology requires that the system be re-condensed. In some applications, it is possible predict which nodes will be necessary for the simulation. Palpation, for example, would require only the surface nodes. For needle insertion or laparoscopy, however, the user should be able to interact with nodes that are under the surface, making it difficult to know how wide around the site, or how deep into the tissue the condensed system needs to extend. The explicit system, u = K~lf, can be reduced using Selective Matrix Vector Multiplication (SMVM) [16], yielding the same system size as if it had been condensed using (A-25), (A-26) and (A-27). The drawbacks with this approach are inversion time, and numerical errors due to inversion. These drawbacks are far outweighed by the advantage of SMVM, that it requires no computation to reduce the system. Reducing the explicit system with SMVM consists only of re-arranging and removing entries from K~x, which allows for the condensed matrix to change on the fly to accomodate for the addition or removal of nodes to the active system [34]. Equation A-29 shows a simple linear FEM system that is to be reduced by removing nodes, for example, 2 and 4 because they are not active, and hence f2,fn = 0. Equation A-29 shows the resulting linear FEM A.2 3D Mesh Generation 77 system after it is reduced, which is analogous to a condensed system. u2 U4 U\ M3 K. 21 31 K, 41 K-31 K 1 2 Kn K 1 4 K 2 2 K 2 3 K 2 4 K 3 2 K 3 4 KA2 K 4 3 K 4 4 K 1 3 K 4 4 fl h h fl h (A-28) (A-29) A.2 3D Mesh Generation The 3D mesh presented in the Chapter 4 was created assuming a regular mesh spacing. This results in a series of parallel planes in space, each containing a a uniform grid of nodes as illustrated in Figure A.5. Each node has a number that is used to generate the connectivity list, and eventually the stiffness matrix. Figure A.5: 3D Mesh of Parallel Slices with Tetrahedral Elements Using the node numbering in Figure A.5 for demonstration, the six tetrahedra required to fill the cube A.2 3D Mesh Generation 78 are: NodeList (A-30) 1 3 7 4 1 2 5 4 1 4 7 5 2 5 8 6 4 5 7 8 4 5 2 8 The connectivity list is refered to as the NodeList, and for each cube in the elastic body, equivalent entries are entered into the NodeList until the entire object has been filled with tetrahedra using a triple nested loop. A.2 .1 Removal of A Node from the Connectivity List Once the NodeList that describes the entire solid object has been created, which in the case of our experiment was a rectangular volume, nodes can be removed from the interior of the elastic object to create a cavity. Using the node locations to determine which nodes should be removed to correspond with the desired cavity geometry, nodes can be removed from the mesh and NodeList. For demonstration, (A-31) describes a 2D system consisting of 9 nodes at regular spacing shown in Figure A.6. Nodes = x\ V i * 2 yi M y4 X5 V5 X6 V6 xi yi x% y% 9 y9 NodeList = 1 2 5 1 5 4 2 3 6 2 6 5 4 5 8 4 8 7 5 6 9 5 9 8 (A-31) A.2 3D Mesh Generation 79 Figure A.6: Node Removal To remove, for example, node 7. Al l triangles that include node seven are removed from the NodeList. The algorithm for removing a node, delete-node, from both the Nodes and NodeList is: 1. Search for elements in NodeList that contain the delete-node 2. Remove that element from the NodeList 3. Subtract 1 from every node number larger than the delete-node in NodeList 4. Remove the delete-node from the mesh It is important to note that when more than one node is removed, the nodes should be removed in descending order to maintain proper connectivity. (A-32) shows the appropriate mesh and NodeList after node 7 is removed. x\ yi X2 yi * 3 V3 X4 V4 *5 ^5 X6 ye Xi V8 X<) V9 1 2 5 1 5 4 2 3 6 Nodes = I " I NodeList = I 2 6 5 I (A"32) 4 5 7 5 6 8 5 8 7 Although a 2D example was shown, the procedure for removing nodes and tetrahedra from a 3D mesh is identical. A.3 Centroid of a Triangle 80 A.2.2 Procedure for Creating a 3D Tetrahedral F E M stiffness matrix 1. Generate 3D node locations. These may be a regular grid, as described in this work, or it may be from a set of contours or a surface representation generated by any number of segmentation algorithms. 2. Create a connectivity list that fills the volume with tetrahedra. In this thesis, the connectivity list was created by packing each of the regularly spaced cubes with six tetrahedra. The use of various automatic mesh generators can be used to create a volumetric mesh from a surface representation, for example. Also, CAD programs like ProEngineer can create mesh from a solid model of the geometry which can be exported in various industry standard data files that define node location and connectivity (such as IGES). 3. Use the mesh node locations and connectivity list as input arguements to a script that creates a stiffness matrix for each tetrahedra, then populates the global stiffness matrix as outlined in the above sections. A.3 Centroid of a Triangle h rJ Figure A.7: Right Angle Triangle Geometry The centriod of a 2D object is defined by: x y IAdA IAydA IAdA • (A-33) (A-34) (A-35) A.3 Centroid of a Triangle 81 Where (x,y) = (x, \y) for an integral in x [47]. Using the triangle in Figure A.7, x is: x = Ipxydx &x{\x)dx \bh Jo 'l^dx \bh A v 3 \ b 3bx |0 2 2b bh (A-36) (A-37) (A-38) (A-39) (A-40) And y is: Soy{\y)dx \bh Ioi(y2)dx \bh fb ]f_Jl J O 2b2 x dx \bh h2 ?, \b lo \bh h 3 (A-41) (A-42) (A-43) (A-44) (A-45) (A-46) This proves that for an arbitrarily shaped right angle triangle, the centriod is one third of the height away from the base. Because any triangle can be represented by two right angle triangles, each with the same height, this relationship holds true for all triangles. Figure A. 8: Triangle Composed of Two Right Angle Triangles Appendix B On Graphics Accelerators for Real-Time Processing B.l Introduction and Previous Work As a result of ongoing competition amongst PC graphics accelerator manufacturers, commodity graphics hardware continues to increase in complexity and computational power. Each design cycle, more of the 3D graphics pipeline is been optimized and implemented in hardware, allowing for higher frame-rates and polygon counts in real-time 3D applications. Researchers have begun to use the computational power of PC graphics accelerators to speed up and create new algorithms and techniques. Purcell et al. [67] have proposed a method of utilizing "Near-Term" programmable graphics hardware to aid in the computation of ray tracing and ray casting. This work makes assumptions about what features future hardware will contain. These assumptions include increased accu-racy in the frame buffer and improved fragment shader flexibility. Proudfoot et al. [66] have designed and implemented a shading system that will take advantage of the programmable features of modern day graphics hardware, while abstracting the hardware from the graphics code. Presently, the only way to access what programmable FPU processing is available is to use proprietary assembly. It is the goal of this work to increase programmer productivity and portability by adding a layer of abstraction. In 2001, Larsen and McAllister [38] showed that it is possible to compute large matrix multiplications 82 B.2 Graphics Accelerator Architecture 83 very quickly on a present day graphics cards due to its high memory bandwidth and optimization for floating point texture operations. Larsen's matrix multiplication algorithm uses two OpenGL techniques. 1. Multi-texturing. Multi-texturing allows developers to render multiple textures onto a polygon. This operation is achieved by a pixel to pixel multiplication in RGB space performed with at least 16 bit floating point accuracy (on a Nvidia GeForce 3 GPU according to [38]). 2. Frame Buffer Blending. When enabled, frame buffer blending will perform many mathematical func-tions on consecutive frames and accumulate the answer in the frame buffer. The primary use of this feature in graphics programming is alpha blending, where the alpha values of multiple frames are averaged to give an enhanced 3D effect. Unfortunately, the PC graphics subsusystem networks with the CPU and main memory over the PCI or AGP bus. This causes read-back of rendered frames into main memory to be very slow. Additionally, at the time, PC graphics accelerators had a fixed output accuracy of 8 bits per channel. According to Larsen and McAllister, these hardware restrictions compromise the effectiveness of the method since present CPU's can achieve similar overall speeds. In 2002, GPUs improved considerabley when 24 or 32-bit floating point pixel formats were introduced. With the improvement of floating point texture processing, GPU performance can be directly compared to CPU performance on a level playing field. Adam Moravanzky [61] presents a comprehensive comparison between the performance of an opti-mized C++ library, ATLAS (Automatically Tuned Linear Algebra Software), with optimized GPU imple-mentations of dense linear algebra operations using the most current programmable features of the graphics hardware. The computation times show that for matrix multiplies, the CPU achieves higher processing speeds, while for linear system solvers like Conjugate Gradient and the Linear Complementarity Problem the GPU achieves higher processing speed. B.2 Graphics Accelerator Architecture Graphics accelerators are designed to take the graphics processing load off of the CPU with special purpose chipsets that are accessible to programmers through graphics API's such as OpenGL. Figure B.l illustrates a typical graphics pipeline. B.3 Larsen's Matrix Multiplaction Algorithm 84 1. Model-View Transform 2. Lighting —• 3. Projection Transform —*• 4. Clipping 5. Primitive Assembly 6. Rasterization Figure B. l : A Typical Graphics Pipeline based on the OpenGL API. Current graphics hardware allows for user defined shader programs to interact in: 1. Steps 1 - 5 by using Vertex Shaders 2. Step 6 with Fragment or Pixel Shaders With the advent of floating point pixel formats, Fragment Shaders can be used to perform pixel by pixel operations such as add, multiply, sin, cos and square root on on values from multiple floating point textures to perform an almost limitless variety of functions on the GPU. All PC graphics accelerator boards have a frame buffer, which is where the fully rendered frame is stored before it is sent out to a display. The frame buffer allows a great deal of flexibility for graphics programmers to blend frames, render frames to textures for later use, or read-back the current frame into main memory for further processing. Unfortunately, most vendors do not implement optimized read-back hardware and software routines, so the read-back is very slow. B.3 Larsen's Matrix Multiplaction Algorithm In order to map a multiplication onto a graphics card using standard OpenGL calls, it has to be broken down into basic texture operations. Presently, this algorithm works only for square matrices. The operations are as follows for two nxn matrices A and B: 1. Use texture addressing to address the ith column of A and the ith row of B. Each respective column and row are copied n times onto a square polygon. 2. Use multitexture to perform element to element multiply of the ith row and columns together (o^ Application^ Framebuffer 7. Texture Fetching & Other Per-fragment Operations B.3 Larsen's Matrix Multiplaction Algorithm 85 3. Perform 1 and 2 for each of the n columns and rows, each time using the blend function to sum each subsequent pair into the framebuffer. Figure B.2 illustrates the algorithm graphically. a l l al2 al3 b l l bl2 bl3 a21 a22 a23 X b21 b22 b23 a31 a32 a33 b31 b32 b33 a l l b l l a l lb l2 a l lb l3 a21bll a21bl2 a21bl3 a31bll a31bl2 a31bl3 + al2b21 al2b22 al2b23 a22b21 a22b22 a22b23 a32b21 a32b22 a32b23 + al3b31 al3b32 al3b33 a23b31 a23b32 a23b33 a33b31 a33b32 a33b33 First column of A multitextured with first row of B Second column of A multitextured with second row of B Third column of A multitextured with third row of B Figure B.2: Graphical Description of Larsen's Algorithm B.3.1 OpenGL Code Details for Matrix Multiplication-Using Larsen's Method Larsen's method for matrix multiplication uses 6 steps 1. bind the matrices that are to be multiplied as textures glBindTexture(GL_TEXTURE_2D, texNames[0]); glTexImage2D(GL_TEXTURE_2D, 0, 4, t e s t M a t r i x D i m , t e s t M a t r i x D i m , 0, GL_RGBA, GL_UNSIGNED_BYTE, t e s t M a t r i x O ) ; 2. enable the modulation (or multiplicative) environment for multitexturing B.3 Larsen's Matrix Multiplaction Algorithm 86 glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE); 3. enable the sum mode of frame buffer blending glEnable(GL_BLEND); glBlendEquation(GL_FUNC_ADD); 4. render a rectangular polygon in orthogonal mode g l O r t h o ( 0 . 0 , t e s t M a t r i x D i m , 0.0, t e s t M a t r i x D i m , -1.0, 1.0); 5. multitexture the polygon repeatedly with columns of the first texture, and the rows of the second (if the matrix is nxn, then the polygon will be rendered n times and each time the result is added to what is in the frame buffer) f o r ( i = 0; i < t e s t M a t r i x D i m ; i++){ dim = ( ( f l o a t ) i ) / t e s t M a t r i x D i m + p a r t ; glBegin(GL_QUADS); glMultiTexCoord2fARB(GL_TEXTURE0_ARB, dim , 1.0); glMultiTexCoord2fARB(GL_TEXTURE1_ARB,0.0, d i m ) ; g l V e r t e x 3 f ( 0 . 0 , 0.0, 0.0); glMultiTexCoord2fARB(GL_ T E X T U R E 0_ARB, dim , 1 . 0 ) ; glMultiTexCoord2fARB(GL_TEXTURE1_ARB,1.0, dim ); g l V e r t e x 3 f ( t e s t M a t r i x D i m , 0.0, 0.0); glMultiTexCoord2fARB(GL_TEXTURE0_ARB, dim , 0.0); glMultiTexCoord2fARB(GL_TEXTURE1_ARB,1.0, dim ); g l V e r t e x 3 f ( t e s t M a t r i x D i m , t e s t M a t r i x D i m , 0.0); glMultiTexCoord2fARB(GL_TEXTURE0_ARB, dim , 0.0); glMultiTexCoord2fARB(GL_TEXTURE1_ARB,0.0, d i m ) ; g l V e r t e x 3 f ( 0 . 0 , t e s t M a t r i x D i m , 0.0); g l E n d ( ) ; ( 6. read the answer back through the frame buffer g l R e a d P i x e l s ( 0 , 0 , w i n _ w i d t h , win_height,GL_RGBA,GL_UNSIGNED_BYTE,buffer); B. 4 Discussion of Results 87 B.4 Discussion of Results The literature has shown that GPUs can be used to perform linear algebra and other useful computations. Table A- l lists approximate cpu times for matrix multiplies. Processor Approximate CPU Time (ms) Radeon 9700 Pro [61] 550 Pentium 4 1.6 GHz (ATLAS) [61] 530 Pentium 4 2.5 GHz (Intel MKL 6.0) 325 Table A - l : Matrix Multiply Comparison Between GPU and Optimized CPU. Size = 1024 x 1024 The times listed are taken either from the results of Moravanzky [61] or from in-house testing using the Intel Math Kernel Library (MKL) Version 6.0. B.4.1 Trade-offs The choice to use a GPU as a cost effective parallel processing unit has several trade-offs. Advantages: 1. GPU processors are heavily pipelined and parallelized to achieve very high processing speeds. 2. GPU processors are a relatively cost effective upgrade to a PC. 3. With current shaders, the user can implement almost any mathematical functions desired to be com-puted entirely on the GPU. 4. GPUs have been shown to process conjugate gradient and linear complimentarity computations faster than an optimized C++ library. Disadvantages: 1. Read-back of data from the framebuffer is relatively slow. 2. The use of shaders entails considerable software development time. 3. CPUs have been shown to process matrix multiplies faster than GPUs using optimized CPU libraries. B.4 Discussion of Results 88 B.4.2 Conclusions Despite the advantages, the use of GPUs for parallel processing of real-time applications is not recom-mended for several reasons. The first and foremost reason is that it is time consuming to transfer large portions of data from GPU memory to main memory. If the data is meant to stay on the GPU for display to the screen, this is not a problem. But, if the data is needed for any other reason, the read-back time is a limiting factor. The second reason is that GPUs do not achieve any processing improvement over optimized CPU libraries. The third reason is that to take advantage of floating point pixel formats that current GPUs offer, programmers must become familiar with shader programming. This requires the programmer to be confined to a limited number of operations per program, and all operations must fit into the context of texture operations. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items