"Applied Science, Faculty of"@en . "Electrical and Computer Engineering, Department of"@en . "DSpace"@en . "UBCV"@en . "Goksel, Orcun"@en . "2009-12-02T21:34:33Z"@en . "2004"@en . "Master of Applied Science - MASc"@en . "University of British Columbia"@en . "Brachytherapy is a prostate cancer treatment where radioactive pellets (seeds) are implanted\r\ninto prostate region through long flexible needles. Steering them to preplanned locations require\r\nthe use of tip bevel deflection effect and a thorough comprehension of tissue deformation and\r\nneedle flexibility. A simulator is proposed to help medical residents in acquiring these procedural\r\nskills. In prostate brachytherapy, two important feedback modalities during needle insertions are\r\nthe transrectal ultrasound (TRUS) imaging and the kinesthetic feedback received through the\r\nneedle base. Thus, these have been incorporated in the simulator, the subsystems of which are\r\npresented here.\r\nIn this thesis, a needle-tissue interaction model from prior work which was based on finite\r\nelement method is extended into 3D and adapted for custom meshes. The issues related to\r\nneedle-tissue coupling are addressed. A meshing scheme which generates soft tissue models with\r\nsmall number of nodes for haptic implementation from contours on parallel slices is proposed.\r\nIt is utilized to generate a prostate volume mesh which consists of well-conditioned elements for\r\nfinite element analysis. A simulation interface is developed where a 3D rendering of the model\r\ndeformation is displayed while the needle can be manipulated by entering target positions for its\r\nbase.\r\nA method for B-mode TRUS image synthesis out of a deformed tissue volume for which\r\nthe patient data was collected at the BC Cancer Agency is also presented. Some interpolation\r\nschemes for remotely spaced ultrasound data are compared on a sample case and a method for\r\noverlaying needle shaft reflection onto TRUS images is described."@en . "https://circle.library.ubc.ca/rest/handle/2429/16194?expand=metadata"@en . "14690589 bytes"@en . "application/pdf"@en . "Ultrasound Image and 3D Finite Element based Tissue Deformation Simulator for Prostate Brachytherapy by Orcun Goksel B.S., Middle East Technical University, Ankara, T U R K E Y , 2002 A THESIS S U B M I T T E D IN 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 O F 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 (Electrical and Computer Engineering) 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 2004 \u00C2\u00A9 Orcun Goksel, 2004 A b s t r a c t Brachytherapy is a prostate cancer treatment where radioactive pellets (seeds) are implanted into prostate region through long flexible needles. Steering them to preplanned locations require the use of tip bevel deflection effect and a thorough comprehension of tissue deformation and needle flexibility. A simulator is proposed to help medical residents in acquiring these procedural skills. In prostate brachytherapy, two important feedback modalities during needle insertions are the transrectal ultrasound (TRUS) imaging and the kinesthetic feedback received through the needle base. Thus, these have been incorporated in the simulator, the subsystems of which are presented here. In this thesis, a needle-tissue interaction model from prior work which was based on finite element method is extended into 3D and adapted for custom meshes. The issues related to needle-tissue coupling are addressed. A meshing scheme which generates soft tissue models with small number of nodes for haptic implementation from contours on parallel slices is proposed. It is utilized to generate a prostate volume mesh which consists of well-conditioned elements for finite element analysis. A simulation interface is developed where a 3D rendering of the model deformation is displayed while the needle can be manipulated by entering target positions for its base. A method for B-mode TRUS image synthesis out of a deformed tissue volume for which the patient data was collected at the B C Cancer Agency is also presented. Some interpolation schemes for remotely spaced ultrasound data are compared on a sample case and a method for overlaying needle shaft reflection onto T R U S images is described. u Table of Contents Abstract ii Table of Contents iii List of Tables vi List of Figures vii Glossary ix Acknowledgements xi 1 Introduction 1 1.1 Prostate Cancer 1 1.2 Brachytherapy 2 1.3 Need for Training 5 1.4 Related Work 7 1.4.1 Tissue Modeling 8 1.4.2 Needle Insertion 10 1.4.3 Meshing Techniques 11 1.4.4 Ultrasound Imaging and 3D Reconstruction 13 1.5 Proposed Simulator 16 1.6 Thesis Organization 16 2 Needle-Tissue Interaction Model 17 2.1 Stiffness Matrix Inversion and Condensation 17 2.1.1 Inversion 18 2.1.2 Condensed Subsystem 19 2.2 Contact Behaviour 19 2.3 Boundary Condition Change 21 2.4 Local Frame Rotation 22 2.4.1 Stick State 23 iii TABLE OF CONTENTS iv 2.4.2 Slip State 24 2.5 Decision of Remeshing 27 2.6 Mesh Modification 29 2.6.1 0(mn2) Approach 30 2.6.2 0(m2n) Approach 31 2.7 Needle Model 32 2.8 Projecting Forces onto the Needle 32 2.9 Integrating Needle Forces 33 2.10 Discussion 34 3 Mesh Generation 35 3.1 Prostate Surface Construction 35 3.1.1 Treatment Planning . 36 3.1.2 Contour Extraction via D I C O M - R T 36 3.1.3 Culling the Contour Nodes 36 3.1.4 Obtaining the Surface Definition 37 3.2 Meshing the Region 39 3.3 Results 40 3.4 Local Element Numbering 42 3.5 Mesh for the Needle Model 44 3.6 Discussion 46 4 Ultrasound Synthesis 47 4.1 Patient Data Collection 47 4.2 The Shift and the Torque 48 4.3 Problem Definition 51 4.4 A Direct Approach 52 4.5 Inverting the Problem . 54 4.6 Image Interpolation Results 56 4.7 Needle Overlay 57 4.8 Discussion 60 5 Implementation 62 5.1 Program Overview 62 5.2 Simulation Structure 63 5.3 Tissue Model Initialization 63 5.4 Needle Model Initialization 65 5.5 A n Iteration 66 5.6 Notes 69 5.7 Discussion 70 TABLE OF CONTENTS v 6 Conclusion 71 6.1 Summary and Contributions 71 6.2 Future Work 72 Bibliography 73 A Finite Element Method 79 A . l Definitions of Stress and Strain 79 A.2 Deformation 81 A.3 Linear Elastic Behaviour 82 A.4 Discretization 83 A.5 Element energy 86 B Boundary Condition Change 89 L i s t o f T a b l e s 3.1 M e s h quali ty 43 5.1 Tissue F E A variables 65 5.2 Needle F E A variables 66 v i L i s t o f F i g u r e s 1.1 The anatomy of the prostate region 1 1.2 Prostate brachytherapy operation 3 1.3 Radioactive seeds implanted during brachytherapy 3 1.4 A representation of the 3D volume and two T R U S B-scans 4 1.5 Fluoroscopy images from a brachytherapy procedure 5 1.6 The brachytherapy needle 6 2.1 Automaton describing the stick-slip friction model of a contact node 20 2.2 Tissue and needle-base axes in our implementation 20 2.3 A n instance of contact node local frames lying on a bent needle 2.1 2.4 Subregion remeshing for a new element penetration 28 2.5 Four elements neighbouring to node P which is to be repositioned 29 2.6 Projecting the working node forces onto the needle knots 33 2.7 Integration of knot forces and torques on the base 34 3.1 Node culling example 37 3.2 Node culling with a threshold of 4 mm on a sample data 38 3.3 Nuages surface output 39 3.4 Mesh visualization 41 3.5 Mesh quality analysis 42 3.6 A right-handed tetrahedron 43 3.7 Brachytherapy needle deflection under a lateral tip force 44 3.8 Needle mesh topology 44 3.9 Mesh topology of a section for two different designs 45 3.10 Lateral tip deflection with respect to number of sections 45 4.1 C T images after a brachytherapy procedure 48 4.2 MRI after a brachytherapy procedure 48 4.3 Prostate shifting into the bladder during a needle insertion 49 4.4 The prostate region and the T R U S images in 2D 50 4.5 Demonstration of the structural deformation during needle insertion 50 vii GLOSSARY vi i i 4.6 A sketch representing the 3D deformation of the initial T R U S images 51 4.7 Demonstration of the structural deformation during needle insertion 52 4.8 Determining whether a point is inside a tetrahedron 53 4.9 Ultrasound imaging plane transformed into the nominal mesh 55 4.10 Ultrasound plane with deformed tetrahedra 56 4.11 Interpolation results and needle overlay 58 4.12 Comet tail and multiple reflections artifacts of the needle. . . 59 4.13 Needle shaft 60 5.1 Matlab interface running the simulation 64 5.2 A main iteration flowchart of needle interacting with the tissue 66 5.3 The simulation forces and timings 69 A . l Resultant force and moment vectors 79 A.2 Components of stress on an infinitesimal cube 80 A.3 A 2-D discretization example of a body Q, into finite elements Q,1 84 A.4 A tetrahedral element and its corners in local numbering 84 G l o s s a r y Roll-pitch-yaw angles 2D Two-dimension(al) 3D Three-dimension(al) B C C Body Centered Cubic B C C A British Columbia Cancer Agency B E M Boundary Element Method bNED biochemical no evidence of disease C D T Constrained Delaunay Tetrahedralization C G Conjugate Gradients C T Computed Tomography D I C O M - R T Digital Imaging and Communications in Medicine - Radiation Therapy D O F Degree of Freedom DW Distance Weighted F E A Finite Element Analysis F E M Finite Element Method H D R High-dose-rate (brachytherapy) MRI Magnetic Resonance Imaging NURBS Nonuniform Rational B-Splines P N N Pixel Nearest Neighbour PSA Prostate Specific Antigen GLOSSARY x P T V Planning target volume R B F Radial Basis Functions T R U S Transrectal Ultrasound US Ultrasound V N N Voxel Nearest Neighbour V R Virtual-Reality x Mixed force-displacement outputs of working subsystem y Mixed force-displacement inputs of working subsystem / Force vector K~l Inverse global stiffness matrix KH Condensed system matrix Kw Working matrix Kw' Updated working matrix u Displacement vector m The number of neighbouring nodes of a node to be relocated n The total number of nodes in a tissue mesh w Number of working nodes Acknowledgements \u00E2\u0080\u00A2 First of all, I would like to thank my supervisors T im Salcudean, Robert Rohling, and Simon DiMaio for their guidance. In this journey of two years, their support has always been a positive motivation for me. \u00E2\u0080\u00A2 I would also like to thank all my lab mates: Daniela, for her always being kind and helpful in whatever the topic of my problem is; Danny, for making our hospital visits enjoyable; Julian, for his endless \"literal\" help; Giselle, for making the lab environment fun; to Emre, for our valuable discussions; Ehsan, for always sharing his vast algebraic knowledge with me; and to all others including but not only Nicu, Hani, Xu, Sara, Hassan, Reza, ... \u00E2\u0080\u00A2 I am also very grateful to Dr. Jim Morris and Dr. Mira Keyes. Their expertise in the clinical setting was invaluable and necessary in order to realize this project. And, without the help of the friendly staff of the B C Cancer Agency, the data collection would not be easy at all. \u00E2\u0080\u00A2 I appreciate the logistic and psychological support of my fellow Johannians through the process of (especially writing of) this thesis. \u00E2\u0080\u00A2 I definitely owe my being here to my parents whom I should acknowledge the most. Their endless support has always been with me. x i C h a p t e r 1 I n t r o d u c t i o n 1.1 Prostate Cancer The prostate (see Figure 1.1) is a walnut-sized gland. It lies below the urinary bladder and surrounds the urethra. It adds important nutrients to the semen, which is the fluid containing sperms. The rectum wall lies just a few millimeters behind the gland and lying on each side of their interface are nerves and blood vessels. The urinary sphincter is a muscle just beyond the prostate and it functions to prevent involuntary leakage of urine. Finally, the vas deferens, a tube originating from each testicle, conducts the sperm to the urethra inside the prostate. Figure 1.1: The anatomy of the prostate region. Prostate cancer is the most common cancer among men in Canada, with 20,100 new cases and more than 4,200 deaths estimated in 2004 [1]. Despite the high rate of the disease (it will afflict 1 in 8 Canadian men in their lifetime), it is a relatively less aggressive and slow processing cancer type especially in the elderly. That is why the decision of treatment method has to be 1 1.2 Brachytherapy 2 made carefully considering involved risks and life expectancy, and in many cases an expectant management (\"watchful waiting\") may be opted for. Some treatment options are: prostatectomy, the surgical removal of the prostate; chemother-apy, the use of drugs; cryotherapy, a new technique of freezing the gland in situ; hormone therapy, which downsizes the prostate and suppresses the cancer; and radiation therapy. First type of radiation therapy is the external beam radiation which focuses beams on the region. The second type, called brachytherapy, applies radiation internally. Both radiation and cryo approaches kill cancerous cells by basically destroying the prostate glands which are later replaced by scar tis-sue. Usually combined treatment modalities are employed rather than a monotherapy in order to eliminate all cancerous tissue within and around the prostate also eradicating extracapsular disease and seminal vesicle involvement. A n important tool for assessing the level of the disease is the serum prostate specific antigen (PSA) which can be tested in a blood sample. 1.2 Brachytherapy Brachy is a Greek term for short distance. High-dose-rate (HDR) brachytherapy refers to the temporary application of radiation intra-operatively using a powerful source in an inserted nee-dle. For early stage locally confined prostate cancer, low-dose-rate brachytherapy which is the permanent implantation of low energy radioactive pellets (seeds) is often the treatment of choice. Its main advantage over other methods is that the delivered dose is very conformal to the prostate and it spares normal tissue. Furthermore, this cost-effective solution maintains urinary conti-nence and it has a short convalescence time. Long term results from the pioneering Seattle group shows the effectiveness of this treatment. They reported a promising 10 year biochemical no evidence of disease (bNED) at 84-93% for patients with no more than one favourable risk factor ( P S A M O ng/mL or Gleason 1 ^) [2]. In spite of the various treatment options, permanent seed implantation (from now on referred as brachytherapy) combined with other modalities has become common because of the flexibility in individual customization of source placement and target volume. A right choice of plan with proper application is very likely to eradicate the cancer while preserving the healthy tissue. How-ever, despite the low risk of brachytherapy, seed placement errors, which may lead to subsequent complications such as urinary incontinence and impotence, are still common. In [3], an experi-enced physician implanting seeds in 20 patients achieved an average placement error of 0.63 cm, a substantial error of 21% of average prostate diameter (3 cm). In addition, a suboptimal ap-plication resulting in an undesired dose distribution at the target volume may in turn decrease the effectiveness of the treatment. Therefore, while not as invasive as surgery, brachytherapy requires significant skill on the part of the performing physician. The implant procedure has numerous preoperative and postoperative components. Most patients are managed preoperatively with external beam radiation to treat locally disseminated 1Gleason score is a cancer aggressiveness grade determined via biopsy. 1.2 Brachytherapy 3 disease and/or w i t h neoadjuvant hormonal therapy to downsize the gland and inhibi t the spread. For planning the ind iv idua l seed arrangement a preimplant ultrasound (US) volume study is done several weeks before the procedure. B o t h the volume study and the procedure are performed wi th the patient i n the dorsal l i thotomy posi t ion (lying on his back). A sketch of the procedure setup is shown i n Figure 1.2(c). A transrectal ul trasound ( T R U S ) probe which is posit ioned using a stepper unit mounted on a stabilizer captures images at 5 m m increments (see Figure 1.2(b)). The template gr id (which is described below) attached to the stepper also discretizes this volume transversally at 5 m m intervals. Stabilizer .Stylet TRUS * probe (a) The operation room. (b) The TRUS setup. (c) A sketch of the procedure. Figure 1.2: Prostate brachytherapy operation. For planning, the ultrasound prostate volume is expanded wi th margins to a planning tar-get volume ( P T V ) , which encompasses likely extracapsular disease. The decision of P T V seed arrangement is made by a radiat ion oncologist using special software to evaluate parameters, e.g., dose dis tr ibut ion. P r io r to the procedure, the sterilized seeds (see \"standard seeds\" i n F i g -ure 1.3) are loaded into needles. The B C Cancer Agency ( B C C A ) uses OncoSeed\u00E2\u0084\u00A2 that is a welded t i t an ium capsule containing radioactive iodine isotope 1 2 5 I (with a half-life of 60 days) absorbed onto a silver rod. ' H III flSH Corrugated seeds I\u00E2\u0080\u0094\u00E2\u0080\u0094 Seeds In vicryl sutures j Figure 1.3: Radioact ive seeds implanted during brachytherapy (reprinted from [4]). Fluoroscopy is a medical imaging modality, which is also uti l ized in brachytherapy, to view 1.2 Brachytherapy 4 the anatomy by means of X-ray . The operation room w i t h the a C-a rm (for fluoroscopy) and U S setup can be seen i n Figure 1.2(a). Al though the C - a r m can rotate for different imaging directions, fluoro-images are only taken along the posterior-anterior axis i n brachytherapy. The intra-operative T R U S probe is set up i n the same way as i n the preimplant volume study so that the prostate can be registered to the implant plan in the beginning of the procedure. Loaded needles are inserted v i a a template gr id (see Figure 1.2(c)) the coordinates of which are given by the preplan. Needle posi t ioning is confirmed using the T R U S (see Figure 1.4) and fluoroscopy (see Figure 1.5). For guidance, white dots denoting template gr id holes are overlayed onto the ultrasound images by the T R U S hardware as i n Figure 1.4(b). Since there is a risk of X -ray overdose, the number of fluoro-images taken during a procedure must be kept low (20-50). Consequently, T R U S becomes the only modal i ty supplying continuous feedback to the physician. Before the procedure the patient is anaesthetized . If anaesthesia is not safe, the spinal blockade is applied so that the patient is temporari ly paralyzed below waist. F i rs t the physician fixes the stabilizer and the precision stepper which are used to posi t ion the ultrasound probe. The stepper is calibrated by matching the T R U S images w i t h the ones that were taken i n advance. These manual adjustments are done by visually comparing the printed images wi th the intra-operative T R U S imaging. It takes relatively long time and it is subject to the experience of the physician. T h e n the fluoroscope is quickly adjusted and the radioactive seeds are taken into the operation room. The brachytherapy needles are loaded wi th up to six seeds each prior to the procedure. If a l l seeds i n a needle are spaced by 1 cm, then a special type of seed packaging (the coated ones at the bo t tom in Figure 1.3) can be used. The B C C A uses the brand R A P I D S t rand\u00E2\u0084\u00A2 - R i g i d Absorbable Permanent Implant Device). Th is manufacturer wrapped strands are preferable since they are less likely to j a m i n the needle. They also stay at the implanted location and orientation whereas the indiv idual seeds may rotate during retraction or can travel through blood vessels. Cus tom spacing is achieved by adding spacer between ind iv idua l seeds. The needles are inserted to their targets through a gr id template. The grid hole and insertion depth are double-checked by a nurse who is reading the therapy plan to the doctor. The stepper may need to be readjusted (a) Needle and a T R U S image. (b) a B-scan in the beginning. (c) a B-scan in the end. Figure 1.4: A representation of the 3D volume wi th T R U S layers and a needle (a), and two B-scans from a procedure. 1.3 Need for Training 5 (a) Prior to the procedure (b) After the 5th needle (c) Towards the end of the procedure Figure 1.5: Fluoroscopy images from a brachytherapy procedure. dur ing the procedure because of the prostate swelling. The doctor may also modify the procedure plan slightly and/or use extra seeds at the end of the operation to correct and improve the dose dis tr ibut ion. Dur ing a procedure wlOO seeds are implanted using 20 to 25 needles. The needles are 20 cm i n length and quite flexible w i th a significant t ip bevel. A needle consists of an outer cannula (that holds the seeds) w i th black stripes marking every centimeter and a stylet that t ightly fits i n this cannula (see Figure 1.6). So, the seeds can be released by pul l ing the needle cannula while holding the stylet steady. Just before that, measuring the distance from the cannula base to the template and to the stylet holder ensures the depth of the current implant and the queueing of the first seed to the needle t ip , respectively. A n arrow on the cannula base marks the direction of t ip bevel. Us ing this bevel and the needle flexibility, the physician can steer the needle. Th i s can be done to refine the grid discretization by reaching in-between template holes, i f the physician wishes to make slight modifications to the preplan. Moreover, a direct insertion cannot reach behind the pubic arch where steering becomes necessary. Proper steering places seeds away from certain anatomical landmarks, e.g. the urethra and the rectum, while not overdosing or underdosing a region. However, this requires the experience of many procedures. Successful implants rely upon not only a thorough comprehension of the relative 3D seed locations and surrounding anatomical structures, but also taking their deforma-t ion during needle insertion into account. These challenges reveal the need and importance of brachytherapy training. 1.3 Need for Training The lack of direct visual feedback makes the comprehension of the needle path wi th respect to the anatomical layout rely only upon the medical imaging modalities, ul trasound and fluoroscopy, 1.3 Need for Training 6 (k(i,i) + 1); k (r,*)J (2.9) [fc(l,i)> [%,1)> (k(i,i) ~~ 1); \u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2) %,r)] So, via an inexpensive low-rank update, the boundary condition of a node for a given axis can be switched between force and displacement. This is similar to the capacitance matrix strategy of James and Pai in [28] and is also implemented in [40]. The CiTi multiplication results in 0(w2) computations for this operation. Throughout this thesis, computations refer only to multiplications and divisions. 0(w2) makes the speed of this algorithm greatly dependant on the number of contact nodes. Although this is a relatively fast operation, its time requirement is not constant per iteration. Its biggest computational burden comes as one boundary change per contact node, when the needle penetration/retraction stops so all contact nodes switch from slip to stick state in a single iteration. Consequently, this sets a limit on how coarse the mesh can be. In the current implementation and meshing, there are usually no more than 10 nodes on the needle. 2.4 Local Frame Rotation Although this section is based on the work by DiMaio [40], the increase in matrix dimensions for 3D needed some modifications on the method, especially for the slip state, which required a new algebraic solution. Furthermore, the need for more than one angle in 3D resulted in extra steps. For example, whereas in 2D an angle-wise addition of two consecutive rotations can solve 2.4 Local Frame Rotation 23 the total rotation, in 3D consecutive rotations needed transformation matrices to be employed. The local frames of the contact nodes seen in Figure 2.3 rotated as the needle twists or flexes. Depending on its state, this rotation is handled correspondingly. 2.4.1 Stick State For a node i in this state, all three axes are under displacement constraints. So, the vectors and \u00C2\u00B0 / j in the initial coordinate frame can be expressed as follows: Rw R12 Rl3 = R21 R22 R23 R31 R32 R33 . \ z . \"\u00C2\u00B0/f\" R11 R\2 Rl3 \" Vf \" OfV J i = R21 R22 R23 J i Ofz . J i . R31 R32 R33 .'ft . (2.10) where \ii and 1fi are the vectors after rotating the coordinate frame and the matrix i?3X3 is the transformation from the rotated frame to the initial frame. This expression is chosen for simplicity in the equations below. (2.4) can then be rewritten in the rotated frame as: 1x= [ATKwA)ly (2.11) \D= Kw' \ where K w is the updated working matrix and: 1 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0 0 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0 0 1 0 0 0 0 A = 0 0 Ru R\2 Riz 0 0 0 R21 R22 R23 0 0 0 R31 R32 -R33 0 0 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0 0 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 1 (2.12) 2.4 Local Frame Rotation 24 If the rotation angles (, 9, ip) are given in roll-pitch-yaw form, then R becomes: R Rll R\2 7?i3 R21 R22 -R23 R3I Rz2 -R33 c cip ccf) s9 cip + scj) sip S(f)c9 S(f> SO S1p + C(f)C1p S(f) sO C1p \u00E2\u0080\u0094 C(f) sip \u00E2\u0080\u0094s9 c9 sip c9 cip (2.13) where s and c denotes sin and cos functions, respectively. This is an 0(1) update per contact node in every iteration. 2.4.2 Slip State While a contact node is slipping, the constraint along the shaft, which is the local y axis, becomes a friction force that has a magnitude equal to the needle-tissue slippage. Thus, from (2.4) the mixed force-displacement constraints are written in the form: J i OfV J i \u00C2\u00B0/f \u00C2\u00B0X3w \ c?3w which can then be separated into: 0 I 0 Ofx J i 0 0 0 + = Kw \u00C2\u00B0 7 ? + 0 \u00C2\u00B0fi 0 0 ^ Zw 0 \ 0 (2.15) Using (2.10) each term above can be denoted in terms of the rotated frame vectors lx and (2.14) Mlx + N \ =KW [JV^ + A f V (2.16) 2.4 Local Frame Rotation 25 where M and N are sparse transformation matrices of the form: 1 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0 0 0 . . . 0 0 1 0 0 0 0 M = 0 0 Ru 0 ^13 0 0 0 0 R22 0 0 0 0 Rz\ 0 -R33 0 0 0 \u00E2\u0080\u00A2 0 0 0 . . . 1 0 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0 0 0 . . . 0 0 0 0 0 0 0 N = 0 0 0 R12 0 0 0 0 R21 0 R23 0 0 0 0 R32 0 0 0 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0 0 0 . . . 0 (2.17) (2.18) (2.16) leads to the expression of xx in terms of \j as: *x= [M - K w N}'1 [Kw M - N] ly (2.19) ke = Kw' \ which gives the updated working matrix K w . The implementation of this formula requires the inversion of [M \u00E2\u0080\u0094 KWN}. Since it is sparse as shown below, its inverse can be computed symbolically and then be hard-coded in the C file as an algebraic set of equations. N has only 4 non-zero values, so KfwxSwN3^x3\u00E2\u0084\u00A2 takes the form [03u)x3(j-i) -^3u)x3 \u00C2\u00AE3wx3(w-j)] ^na^ n a s 3 non-zero columns, where j is the working node to be 2.4 Local Frame Rotation 26 updated. Thus, the matrix to be inverted is expressed as: 1 0 \u00E2\u0080\u00A2 a(l,3j-2) a(l,3j-l) \u00C2\u00B0(l,3i) \u00E2\u0080\u00A2 0 0 1 \u00E2\u0080\u00A2 a(2,3;-2) a(2,3j-l) a(2,3j) \u00E2\u0080\u00A2 0 M - KWN = 0 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0(3j-2,3j-2) a(3j-2,3j-l) a(3j-2,3j) \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0 0 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 a(3j-l,3j-2) a(3j-l,3j-l) a(3j-l,3j) ' \u00E2\u0080\u00A2 0 0 0 \u00E2\u0080\u00A2 ' a(3j,3j-2) a(3j,3j-l) Q(3j,3j) \u00E2\u0080\u00A2 0 0 0 \u00E2\u0080\u00A2 a(3u>,3j-2) a(3w,3j-l) a(3w,3j) \u00E2\u0080\u00A2 1 For convenience, elements of the 3 x 3 central matrix will be renamed as: On 0-12 \u00C2\u00AB13 \u00C2\u00AB21 \u00C2\u00AB22 \u00C2\u00AB23 3^1 \u00C2\u00AB32 \u00C2\u00AB33 a(3j-2,3j-2) \u00C2\u00AB(3 i -2 ,3j - l ) G(3j-2,3j) a(3j-l,3j-2) \u00C2\u00AB ( 3 i - l , 3 j - l ) a(3j-l,3j) \u00C2\u00AB(3j,3j-2) C(3i,3j-1) a(3j,3i) Symbolic toolbox of Matlab 1 gives the inverted matrix in the following form: 1 0 \u00E2\u0080\u00A2 a(l,3j-2) a(l,3j-l) a(1.3j) ' \u00E2\u0080\u00A2 0 0 1 \u00E2\u0080\u00A2 \u00C2\u00B0(2,3j-2) a'(2,3j-l) 0(2,3j) \u00E2\u0080\u00A2 \u00E2\u0080\u00A2\u00E2\u0080\u00A2 0 [M-KWN]~1 = 0 0 0 \u00E2\u0080\u00A2 0 \u00E2\u0080\u00A2 a'l2 a22 a'l3 a23 \u00E2\u0080\u00A2 0 \u00E2\u0080\u00A2 0 0 0 \u00E2\u0080\u00A2 a31 a32 a33 \u00E2\u0080\u00A2 0 0 0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 n' a(3w,3j-2) a'(3w,3j-l) a'(3w,3j) ' \u00E2\u0080\u00A2 1 (2.20) (2.21) (2.22) where the elements are: for the rows 3 J - 2 , 3 J - 1 , 3 J : for the rest of the rows a e{l,...,3j'-3,3j+l,...,3u;} 11 a21 = a31 = [ a - c n -Cl2 -C13 *12 -C21 a13 = -C 3 1 a23 = -C32 a33 = -C33 a22 = ~C22 a32 = _ c 23 *(3ij-) - CnO( 3 i j i) + Ci 2 a ( 3 l , j ) + C l30(3i , j ) a'(3z+lJ) = c21a(3i+lj) + C 2 2 0 ( 3 j + i i j ) + C23a(3i+l,j) [ a ( 3 i + 2 , j ) = C31a(3i+2,j) + C32a(3i+2,j) + C33a(3z+2,j) Matlab\u00C2\u00AE is a registered trademark of The MathWorks, Inc. (2.23) 2.5 Decision of Remeshing 27 where the c values are found in the following form: Cll = 0 2 3 0 3 2 -022033 r C12 = 0 3 3 0 2 1 -023031 r Cl3 = 0 2 2 0 3 1 -032 021 r C21 = 0 3 3 0 1 2 -013032 r C22 = 013031 -011 ,033 r C23 = O l i a 3 2 -031012 r (2.24) C31 = 0 1 3 0 2 2 -012023 r C32 = 023 0 U - O 1 3 O 2 1 r C33 = O12O21 - 0 2 2 0 1 1 r while r = 0 1 3 0 2 1 0 3 2 - 0 1 3 0 2 2 0 3 1 - 0 2 3 0 1 1 0 3 2 - 0 3 3 0 1 2 0 2 1 + 0 2 2 0 3 3 0 1 1 + 0 3 1 0 1 2 0 2 3 -In every iteration, this update requires 0(w) computations per slipping node. 2.5 Decision of Remeshing Due to the simulation design and implementation, the manipulated tissue nodes can only be on the needle. To relax this restriction a virtual coupling could have been achieved by interpolating the forces of the surrounding nodes and linearly distributing the needle displacements onto them similarly to the approach of Berkley et al. [39]. However, this would make the application of boundary constraints difficult. Also a larger working matrix would then be necessary which would make the boundary condition changes slower. Thus, the approach of constraining the tissue nodes onto the needle is preferred. Alterovitz et al. [9] and DiMaio [40] use similar techniques as parts of their applications. In brief, as the needle is advancing when its tip enters into a new tissue element, the closest corner of this newly intercepted element is rediscretized in order to have a mesh node coinciding with the needle. Figure 2.4 depicts a 2D case. Note that all the element geometry shown in here is from a deformed mesh. Although that geometry's nominal structure is also referred in this section, it is not drawn for representation simplicity. While the needle is advancing in the tissue, if its tip penetrates into a \"new\" (having none of its corners already on the needle) element as in Figure 2.4(b), then the node addition steps are taken as follows: 1. The point that P should move to is determined (it is P' in our implementation). 2. This point (in ours, P') is mapped to the nominal mesh. 3. The region is remeshed as described in mesh modification (Section 2.6): The nominal position of P is moved to the nominal position of the desired point that is found in step 2. 4. The condensed matrix is updated to accommodate this new working node. 5. Along with the other working node constraints, this node is also constrained on the needle (in our implementation, confined to the tip at stuck state initially). Although from a general perspective, not all the intercepted elements need to have a corner as a working node (let us say, some of the nodes could be skipped), this was not preferred since 2.5 Decision of Remeshing 28 (a) Advancing towards a new ele- (b) Tip penetrated into new ele- (c) Remeshing the region, ment. ment. Figure 2.4: When a new element is punctured, the subregion is remeshed in order to have a node around the needle tip. In (c), solid lines represent the new elements to be formed and the dashed lines denote the elements to be removed during remeshing. more nodes would achieve a tighter coupling between tissue and the needle. Moreover, the above mentioned approach of adding all the nodes leads to some useful assumptions (such as \"the tip can only be in an element connected to the last working node on the shaft\") which can speed up other algorithms throughout the implementation. The tip penetration check is performed using the Algorithm 4.2 in Section 4.4. Each element connected to the last contact node is checked for tip inclusion. For step 1 of node addition, [40] preferred the point of edge intersection, P\" in Figure 2.4. First this was implemented in 3D as well. It required plane-line intersection calculations, as opposed to the line-line case of 2D. Later this was changed to P' considering this node is to be stuck at the tip in step 5. This also removed the intersection calculations. In step 2, nominal position calculation is easily performed by [40] in 2D for P\". The ratio of |PP\" | to I-PP2I is applied on the nominal mesh. However, the 3D case required a slightly longer analytical solution. The solution for an inner node, P' was even harder analytically. Thus, the author made use of the basis functions to correctly map this point to the nominal mesh in fewer well-structured calculations. First the basis function from the deformed element (PP2P3) to its nominal geometry is found, and then using this basis function P' is identified in the nominal mesh. This is the same method (Algorithm 4.3) applied for the pixel positions of the ultrasound imaging plane in Section 4.5. Unfortunately, using flexible needle the above scheme caused instabilities when constraining the tip to either P' or P\". After confining the node on the needle in step 5, there has been a residual force on that node of the mesh, which only settled in many iterations, if it converged at all. Further investigation revealed the fact that finding the exact point which would deform onto P' is theoretically impossible. When a mesh node is relocated, since the global stiffness matrix changes, even the point P will not be at its previous location with the same constraints applied on the contact nodes. This is due to the different response of the new tissue discretization. Nevertheless, preventing this instability caused by the slight position miscalculation is only a 2.6 Mesh Modification 29 matter of using a higher damping for the tissue-needle coupling described in Section 2.7. Briefly, since both models are elastic, their coupled system has a tendency to resonate and it needs to damped. The above mentioned tip residual force is a case where this coupled system is pushed to its stability limits. Therefore, the damping was determined empirically where the needle does not resonate too much during insertion but still responds quickly. 2.6 Mesh Modification Since the global stiffness matrix K is obtained by summing up all individual element stiffness matrices (Ke), subtracting them will effectively remove the elements' stiffnesses from this global matrix. For example, consider a part of the tissue that is formed by the elements e\ = P2P4P3P, e2 = P1P3P4P, e3 = P1P4P2P, and e4 = P1P2P3P as seen in Figure 2.5. Removing those four elements from the global stiffness matrix K decouples the node P from the rest of the tissue. In other words, it inhibits all the strain effects caused by P on other nodes and caused on P by other nodes. Consequently, we can safely relocate it to anywhere (let us say P') inside the tetrahedron P1P2P3P4 and add the new tissue elements ex' = P2P4P3P', e2' = P1P3P4P', e3' = P i P 4 P 2 P ' , and c\ = P1P2P3P' onto the global stiffness matrix K to have an intact mesh. Here, note that in most mesh geometries the nodes usually have more than four neighbouring elements. Yet, the descriptions below are made without loss of generality, so they are independent from the number of surrounding elements. P 2 Figure 2.5: Four elements neighbouring to node P which is to be repositioned. Thus, the general form of the matrix update for a displacement it of a node P which relocates it to the position P' \u00E2\u0080\u0094 P + u is written as follows: K' = K - global (Kei (P)) + ^ global (Kei (P')) (2.25) where ej are the elements connected to P. The global function maps the element matrices to 2.6 Mesh Modification 30 the global mat r ix numbering (Appendix A . 5 ) . Note that this modifies only a submatr ix of K corresponding to the corner vertices of elements ej. For simplicity, throughout the rest of this section, let us assume that the matrices are reordered i n order to have the nodes that are directly connected to P at the end and call those nodes P j . Let m be the number of such nodes (namely P ' s valence or order). 2.6.1 0(mn2) Approach The real-time code uses the inverse stiffness mat r ix K^1 as in (2.1). However, Ke defines a relationship from the displacements to the forces. Therefore, a direct operation on this K - 1 is not possible. One approach can apply the operations of (2.25) on the or iginal global stiffness matr ix K and then invert the entire full-rank matr ix . A s discussed earlier, a real-time inversion of rank 3n is infeasible. A fast approach is the extended version of Sherman-Morr ison formula, namely the Woodbury formula [90]. If A K = - X)e. K & I (P) + 2 J e i K & I ( P ' )> then the inversion of the perturbed mat r ix is given by: (K + VAKVT) 1 = K L - \K-1VAK(1 + VTK'1VAK) 1 VTK~1 (2.26) Here V is a 3n x 3m matr ix having the last row as an identity matr ix of the form Ismx3m- Th i s requires 27nm(n + 2m) computations, i n addi t ion to a 3m x 3m matr ix inversion which can be done i n 0(m3) i n the worst case. A l though n > m is a val id assumption, m w i l l also be included in the b ig O notat ion to express the effect of node valence. Instead of this 0(mn2) approach, D i M a i o [40] employed another 0(mn2) method w i t h sl ightly more computations (162mn 2 ) . He converts a part of K~L to displacement constraints, and after applying the Ke' operations on this region, converts it back. The conversion is done by swapping the 3m rows (corresponding to Pi) at the end using boundary switches i n Section 2.3. Tha t yields to: Up1 UPn-n F P r n-L\2 l21 l22 Pn-UPn-(2.27) Here u and / values are triplets containing a l l three axes. Naming the vectors conveniently gives: L N L\2 l21 l22 u7 (2.28) Note that none of the submatrices L holds values of either K or K 1 directly. Instead they 2.6 Mesh Modification 31 define a mixed boundary condition similar to (2.4) such that: U\ = LuFi + liyiUl F2 = L21F1 + L22U2 (2.29) The advantage of this conversion is that now Kei subtractions and additions can be applied on L 2 2 , since all the relationships of the displacements to the forces of the nodes Pj are now imbedded in it. Thus, finally subtracting the previous elements' stiffnesses and adding the new elements' stiffnesses, and then applying the boundary condition switches once more onto those nodes (to convert them back to original force constraints) results in the (K~1)', the inverse global stiffness matrix of the new mesh. Note that finding and adding/subtracting the individual element stiffness matrices, Kei, take negligible time compared to the other operations. 2.6.2 0(m2n) Approach We propose an approximation to this matrix update in 0(m2n). Before describing this, let us examine how the remeshing of a node P can affect the inverse global stiffness matrix K~l. Firstly, changing the position of P with respect to the other nodes will alter two factors, the displacement caused on P by other nodes' forces and the displacements of other nodes caused by the forces on P . In other words, a node's position in the mesh definitely determines the values of the columns and rows of K'1 corresponding to P . Note that if the tissue deformation was not modeled using discrete elements, this effect of P on itself would be the only required change for a position update. However, since the effects diffuse in finite element model through elements, surrounding nodes of P affect each other also through these elemenets. For instance considering two neighbours of P that are P i and P2, the displacement of P 2 is also affected from a force on P i through the element PPiP2 so that the position of P can affect their mutual behaviour. We will call these the secondary effects. In fact, there are tertiary effects and so on, meaning that even the most distal nodes reciprocal behaviour is slightly affected by where exactly this node is located. As mentioned earlier, all effects except for the primary are side-effects of F E M . Obviously in reality, the mutual behaviour of the distal nodes should have nothing to do with the virtual partitioning of the in-between material. To see this, we analyzed an instance of K.~x matrix for a regular mesh with 240 nodes and 810 tetrahedra. This example was remeshed at a node and the difference matrix (K\"1)' \u00E2\u0080\u0094 K~l was found. The sum of squares of the elements of this matrix is called the total power oj difference. Then, it was seen that 21% of this power was concentrated on the rows/columns of P for this instance. When we also included the secondary nodes, this ratio reached to 96%. Experimenting with similar matrices revealed that including only primary and secondary effects can mostly approximate the final matrix with an error less than \u00C2\u00AB 5 % in terms of power of difference. We implemented this with a modified version of the boundary condition change function which only calculates Li2,L2i, and L22 of (2.29). In there, Ln is the effect of other nodes' forces 2.7 Needle Model 32 (Fi) on other node's displacements (U\). Neglecting Ln in the boundary condition change reduces the complexity of each switch to 0(mn). Consequently, this eases the mesh adaptation to 0(m2n). However, this method needs a better validation and verification. Extreme case with strange meshes (e.g., having slivers) should also be studied in depth. 2.7 Needle Model In the simulation, the needle is represented by line segments which are connected to each other at the needle knots. In Figure 2.6, each is a needle knot and each Pf P(\_i is a needle segment. As described in Section 3.5, each knot is defined as the controid of a needle section interface so that the sections are dual for segments in the F E A . The needle flexibility is modeled using 3D F E M which was derived from [40]. A quadratic strain solution is applied on this model so that the relationship between the needle node dis-placements and forces are expressed as follows: => f = K ( u ) u K, + K,(u) u K<0) then return False 6: end for 7: return True Figure 4.8: Determining which side of the plane PiPmPn the point X is. This is one of the four checks before X \u00E2\u0082\u00AC e conclusion. After calculating the positions of the known intensity values in the deformed mesh using Algorithm 4.1, the intensities on L can then be interpolated. This interpolation could be defined in numerous ways. A n inherent nature of ultrasound is the speckle. Its intensity and distribution inside a tissue region is caused by the superpositioning of the sound reflecting from unidentifiable scatterers in a complex manner. Therefore, these B-mode speckles (the dark/bright spots) do not really convey much information about the underlying tissue. Thus, their spatial interpolation from coarsely sampled data using most common techniques may just give too blurred images where black and white speckles gray out each other. To properly interpolate the data at point A, the points on which the interpolation will work on > For each face of the element > P's ought to be in right-hand order > Find the surface normal > Dot product finds the projection 4.5 Inverting the Problem 54 and their weightings on the output have to be determined. For such problems there are scattered data interpolation techniques available [77]. Pixel nearest neighbourhood (PNN), voxel nearest neighbourhood (VNN), radial basis functions (RBF) [80] are some of the methods that were attempted in literature. A valid approach is finding the closest points, Bi and Bj+i, on each curve (see Figure 4.7). However, this would bring an extreme computational burden. A n approximation could be the use of points Ci and Cj+i. They are the transversally (height-wise in 2D) closest points of L\ and L*+1 to A. Once the points are determined, the intensity of A can be found by linear or quadratic interpolation. For instance, linear interpolation can be applied as follows: = IBi\Bi+1A\+IBi+1\BjA\ A \BiA\ + \Bi+1A\ [ ' ' where Ip denotes the intensity of P. One drawback of this direct approach is the algebraic complexity in finding the proper points (ideally the closest ones, e.g., B{ and -Bj+i for A) on the neighbouring curves. Further disadvan-tages arise when L cuts one or more of L*, which happens if the tissue shifts more than 5 mm. In this case, more than two L* need to be found by Algorithm 4.1. Dealing with all separate cases complicates the problem algorithmically. Therefore, the following simpler yet accurate approach was followed. 4.5 Inverting the Problem This ultrasound synthesis problem can be tackled by inverting the steps above, so that the points to be interpolated from would be in a grid structure letting us use the regular-grid interpolation techniques which are fast and well-known. Considering the input T R U S slices already have regularly-sampled voxels on them, leaving those as they are might be beneficial. Thus, let us imagine the deformation of the tissue as a backward deformation of the imaging plane L as seen in Figure 4.9. This is equal to finding the set of initial (nominal) mesh points, which the given deformation would relocate onto the imaging plane. Note that \"a surface that is fit on those points should look similar to the mirror image of the deformed ultrasound input data planes\". The positions of the imaging plane points can be transformed from the deformed mesh to the nominal mesh using a similar algorithm to Algorithm 4.1. However, this time the basis functions have to be calculated at every run, meaning every iteration. That is simply because the transformation is now from the deformed elements to the nominal ones and the basis functions of these deformed elements change with their changing geometry as the tissue is manipulated by the needle. Figure 4.10 demonstrates an instance, where P? are the nominal corner positions of the element e. The dotted lines denote the nominal mesh whereas the solid lines depict the deformed geometry with vertices Pf. 4.5 Inverting the Problem 55 Lj.i L Lj L i + ] Figure 4.9: Ultrasound imaging plane transformed into the nominal mesh. Algorithm 4.3 presents the pseudo-code for this deformed-to-nominal transformation. For each point X* on the imaging plane L, it first finds the deformed element e* that confines X* (lines 2-3). Then using the basis functions of this element e* and the negated displacements, \u00E2\u0080\u0094 u | , the corresponding nominal position X is calculated. In our implementation, each calculated basis function was stored in an array to speed up further references to the same element. A l g o r i t h m 4.3 Finding the nominal positions of the imaging plane pixels, l: for all X* E L do 2: for e* using Pf as x in (A. 18) 5: X <- X* + YA=I Nf (x)(~ui) > Find nominal position for US data point using (A.17) 6: end for A drawback of this scheme is the low speed of the ISPOINTINSIDEELEMENT function, which bottlenecks the current synthesis process. Although it was implemented in Matlab, even a C implementation with optimizations would be slow since checking all mesh elements is not an efficient approach. First a simple yet effective addition was made before Line 2 to limit this search only to the elements cutting the imaging plane. Elements having nodes on both sides of the imaging plane were found in 0(n) computations. For example, in our current mesh (given in Chapter 3) the imaging plane passing through uBase+2.0\" cuts 337 of the total 2801 elements and this can be found in relatively negligible time. This technique is open to many other further optimizations, some of which are: \u00E2\u0080\u00A2 For each voxel, first the element of a neighbouring voxel (the parent element of which has already been found) can be checked considering that X* is most probably inside the same element with its neighbours \u00E2\u0080\u00A2 Before traversing the elements, sorting the elements with respect to their proximity of X* 4.6 Image Interpolation Results 56 ,P!I L Figure 4.10: A n instance showing the imaging plane, the deformed mesh (solid lines), and its nominal geometry (dotted lines). to their centroids increase the chance and speed of finding the right element. (This sorting could just be done once in a while) \u00E2\u0080\u00A2 Element connectivity graphs can be used to make an educated guess of the next element trial. To avoid element traversal, analytical methods could be employed as well. For example, finding the intersections of the element edges with the imaging plane gives the triangular patches of L spanned by each element. We can then easily tag the points inside these patches. Real-time consecutive image generation brings even more options of optimization, with the main assumption of relatively slow deforming tissue with respect to the iteration rate. T R U S is working at 12 Hz imaging rate during the operation. Although the current implementation cannot achieve this speed, it is anticipated that a C code using a combination of above mentioned methods can easily reach this rate. 4 .6 I m a g e I n t e r p o l a t i o n R e s u l t s For the inverted problem, the ultrasound was synthesized using the fast interpS function of Mat-lab. The amount of input data is not noticeably changing the speed of this function. Therefore, instead of only the adjacent planes, the voxel locations and intensity values of all the initial data 4.7 Needle Overlay 57 planes were sent to the function input. As the output, the intensity values at the transformed positions X are asked for. For instance, in a deformed mesh an ultrasound synthesis at uBase+2.0\" is an interpolation probably between the initial images of uBase+2.0\" and uBase+2.5\", considering the shift of the latter towards the base during needle insertion. Those initial images of a patient can be seen in Figure 4.11(a) and 4.11(b), respectively. To see the interpolation effects caused by the deformation, the simulation implementation in Chapter 5 was used. For the given pubic arch-prostate bonding, an insertion at the lower left corner of the prostate twists the region coronally as well as sagittally. This caused half of the ultrasound image to be close to one input data plane and the other half to be closer to the other, which generated a good case to see interpolation artifacts. First the nearest-neighbourhood interpolation was tried. In Figure 4.11(e) the fracture line separating the regions can easily be seen around the column E. This input plane switch arti-fact and other discontinuities caused by skipped image lines completely disappeared when linear interpolation was used as seen in Figure 4.11(f). Furthermore, the blurring caused by the in-terpolation was minimized using cubic interpolation as demonstrated in Figure 4.11(d). On the average interpS finishes the linear interpolation in 0.6 sec and the cubic interpolation in 3 sec. For all three interpolation schemes, it can easily be observed that the left side of the images bear features from \"Base+2.0\" and on the right of the column E they represent \"Base+2.5\" more. To give examples on Figure 4.11, the hospital naming convention will be used. It names the row 1, 1.5, 6, 6.5, and the columns A, a, f, G. So, let us say that the white spot around e6 and the seed at Fl on \"Base+2.5'\" in contrast do not exist on \"Base+2.0\". As expected, they are clearly visible in the generated images where this region of the imaging plane is still close enough to \"Base+2.5\". Similarly, the seed at a2 on the plane \"Base+2.0\", which is one of its distinguishable features, is observed on the interpolated image, whereas the seed of \"Base+2.5\" at C5 does not show up at all as expected. 4.7 Needle Overlay Intra-operatively, the main use of the ultrasound is seeing the placement of the needle and the seeds, rather than observing the deformation of the underlying tissue. That is why the needle appearance is also needed on the ultrasound images. As the seeds are implanted, they are also have to be added on the images. However, to show the concept of overlaying these info, only the needle shaft case is presented here and a discussion of seeds and the needle tip is included at the end. Since during the insertions, the needle shaft is mostly parallel to the probe, the rest of this section assumes that they are always aligned. In the ultrasound image, the needle shaft looks basically as a bean-like bright spot followed by its comet tail. This bean-shape is caused by the orbital elongation of features under T R U S and it can be approximated with an ellipsoid curved around the probe. From here on, its slight 4.7 Needle Overlay 58 (c) needle overlayed on cubic (d) cubic (e) nearest-neighbour (f) linear Figure 4.11: T R U S slices (a,b), interpolation results (e,f,d), and needle overlay (c). 4.7 Needle Overlay 59 curvature will be neglected and it will be referred as the ellipsoid, which is a valid assumption especially distant from the probe. The comet tail artifact is a narrow pie-like region, which appears behind the shaft (see Figure 4.12). It is mainly caused by the repetitive sound reflections inside the needle. Any sound wave that enters in the shaft, starts reflecting back and forth inside the needle while loosing some fraction of itself at every bounce. Since these consecutively escaped sound in time is interpreted as coming from different depths by the ultrasound hardware, the comet tail forms and it consists of repetitive bright and dark stripes. If the speed of sound were the same in steel and tissue, the bright stripe separation would be twice the needle thickness. However, because of the higher speed of sound in steel, their separation becomes less than that. Figure 4.12: Comet tail and multiple reflections artifacts of the needle. Another artifact on the comet tail arises from the high reflectivity of the shaft, which is caused by the relatively high acoustic impedance of the needle with respect to the tissue (ztiSSue ~ 1.54 x 106 Rayles, zsteel ~ 47 x 106 Rayles, where Rayle is kg/m 2s). That creates multiple reflections artifacts, where obscure copies of the shaft are seen at the integer multiples of the actual shaft depth. This is caused by the sound that reflects back and forth between the needle and the probe. When the initial needle reflection hits the probe surface, most of it reflects back to the tissue medium perpendicularly, and it is then reflected from the needle surface once more. This second reflection and possible further ones generate the multiple less-bright shaft looks behind the actual one (see C4.5 in Figure 4.12). Considering the above phenomena, the factors affecting the look of the needle shaft on the image can be categorized in two. The first is the distance from the probe to the shaft. This mainly determines the brightness and the radius of the ellipsoid denoting the shaft, and the separation between multiple reflections. It also indirectly affects the brightness of the comet tail and multiple reflections. The second factor is the angle between the horizon and a virtual line that is perpendicular to the probe surface and passing through the needle shaft (for instance, any point on column D has this angle of instance equal to 90\u00C2\u00B0) . This angle is defined on the US image plane and it determines the orientation of the comet tail and the ellipsoid axis. Thus, the comet tail extends on this virtual line and the longer axis of the ellipsoid lies perpendicular to 4.8 Discussion 60 this line. For needle shaft synthesis, real needle ultrasound images were used. They were acquired at the hospital i n a phantom using the same ultrasound system and probe settings that are used during the procedure. A n acquired image can be seen i n Figure 4.13(a). F i r s t its background was removed by subtracting the image of the region before needle insertion, then the remaining noise was el iminated using an adaptive thresholding filter. T h e processed \"clean\" version is shown i n Figure 4.13(b). (a) raw (b) filtered (c) overlayed (d) real Figure 4.13: Needle shaft. In the code, we are only simulat ing the ellipsoid overlay at the proper posit ion and orientation. A l though mult iple reflections are relatively easy to implement, comet ta i l synthesis needed further in-depth study. Since this thesis is not focused on U S artifact modeling, better artifact synthesis schemes are yet to be researched in future work. For the shaft image, the posit ion where the needle intersects the imaging plane is easily found using the needle knots just before and just after the imaging plane. Th i s intersection gives where the needle should be positioned on the T R U S . W i t h this information and using the pr ior i -known coordinates of the probe, the angle of incidence can be calculated. The cleaned shaft image is then rotated accordingly and its brightness is adjusted to match w i t h the T R U S image. Since adaptive thresholding filter leaves sharp edges, it is also smoothed wi th a 7 x 7 Gaussian filter. Final ly , it is added on to the synthesized ultrasound image where its brightest p ixe l would correspond to the desired shaft center coordinate. The final generated image can be seen i n Figure 4.11(c) w i th the overlayed needle between 2B and 2.5b. 4.8 Discussion In this chapter, two methods for generating T R U S images of soft tissue, the deformation of which was given by F E M , have been presented. The first one (in Section 4.4) has an advantage of precomputing basis functions and tagging the in i t ia l ly known voxels w i th the F E M elements they reside in . However, its interpolat ion requires many geometrical computations which would make its implementation slow. Therefore, the second one (in Section 4.4) was opted for. It makes use of the well-known regular-grid interpolations by transforming the plane of interest into the in i t i a l undeformed mesh. 4.8 Discussion 61 Spatial optimizations (e.g., the displacements of some are computed, and the rest is ap-proximated) or temporal estimations (e.g., a voxel stays mostly in the same element in-between iterations) can be applied. For the first presented method (namely direct method), there exists approximations for the interpolation. These should be addressed in the future work to decrease the running time of this synthesis process down to \u00C2\u00AB 0 . 1 s to be able to generate real-time T R U S images (original rate is 12 Hz). Using 3D ultrasound would also give better resolution between slices. The computational times in the end showed that the FEM-based transformations take in the order of minutes whereas interpS function could interpolate in seconds (for linear ?s0.6 s, for cubic \u00C2\u00AB 3 s). The cause of the slow FEM-based transformation was its running on all the points (~ 1,400 pixels). Approximations can be used for increasing this performance. The negative points on the speed side were compensated on the quality side. The main problem could be the discontinuities in the image. The qualitatively comparisons assured the smoothness and realism of these generated images. However, a user study with sonographers would give a better idea. The needle shaft appearance has also been simulated on a generated T R U S image. This was a preliminary attempt to show the possibility of synthesizing the needle tip, brachytherapy seeds, and their artifacts. It was seen that there were no previous studies on this. In this work, the feasibility of overlaying artificial images, such as the shaft reflection, without disturbing the realism of the underlying image has been studied. It has given promising results, however, once more modalities (e.g., artifacts) are added, the issue of validation has to be addressed. The current needle overlay can be validated through the following scheme: in the same phantom which the initial needle image was taken, other images at different angles can be taken and these can be compared with the needle overlayed ones generated by the algorithm above. Chapter 5 Implementation In this chapter, the implementation of the concepts given throughout the thesis is summarized. Although the F E A and ultrasound synthesis implementations below are currently not working in real-time, during their design and programming, their portability to a real-time system was given the utmost priority. 5.1 Program Overview The tissue interaction is handled by a C code. This C - M E X structured file can be compiled in Matlab so that when it is called as a function with the state variables such as the needle base position, it will return the feedback values such as the needle base forces. Thus, a simulation of the system is possible such that the above function can be called repetitively while modifying the variables representing the scene change or disturbance. Each call represents an iteration cycle of a real-time haptic control loop. At each iteration of this code, first the passed values and pointers to arrays and matrices are checked for any discrepancy among them. For an average iteration where the needle is in the tissue, the following operations are performed in the given order (they are detailed later in the chapter): 1. Transform the needle knot positions from the local needle frame into the world frame; find the orientations of the needle segments; and calculate the rotations from the previous iteration. 2. If a new tissue element is intercepted, modify the nominal mesh so that the closest node would be on the needle and add this node to the working system matrix. 3. Considering the needle direction, insertion/retraction, and the forces on the working nodes, decide the finite state on each working node and if the node slips off the tip, remove it from the working system matrix. 4. Set the constraints of each working node depending on its stick-slip state. 62 5.2 Simulation Structure 63 5. Solve the working subsystem to find the displacements/forces on the working nodes. 6. Calculate the relative position of each working node on the needle and determine the tissue contact forces (shaft friction and tip puncture) on them to be used in the next iteration. 7. Project the tissue forces of the working nodes onto the knots of the needle model which is in a separate frame by itself. 8. Sum up the working node forces to find out the force feedback to be applied at the needle base on the user's hand. 9. Solve the needle model for its knot displacements in its own local frame. 10. Multiply the working node forces (since forces at other nodes are zero) by the nominal inverse stiffness matrix to find the displacements of all the nodes. 5.2 Simulation Structure For simulation of the insertion the C - M E X file is compiled into a D L L executable. Through an interface designed in Matlab, the user can input the next destination of the needle base. Then the Matlab program calculates a linear trajectory and runs a loop, in which first the needle base position is updated and then the executable is called. The tissue and the needle vertices are used by the interface to draw them on screen. Since Matlab drawing functions are relatively slow, the drawing is done at every Nth iteration, where TV can be set through the interface. Figure 5.1 shows screenshots of the interface. At the beginning of the simulation, the matrices modeling the tissue and the needle are initialized. Since the tissue model takes up to 30 sec to compute, when tissue matrices are generated once, they are saved in a file to be loaded in further runs. Then the initial values of the variables are defined according to the initial setup of the volume. Finally the iterative scheme in Section 5.5 is started. Currently the TRUS image synthesis is implemented separately and it is initiated by a button in the interface. This generates a slice at a constant depth in order of minutes. That is why the current implementation is not suitable for a visual TRUS rate of 12 Hz. Further approximations and optimizations have been suggested in the related chapter (Chapter 4). 5.3 Tissue Model Initialization The variables are named here to make them easier to follow on Table 5.1 where n is the number of nodes, e is the number of tetrahedra, t is the number of materials, b is the number of constrained nodes, and v is the maximum node valence. First, node numbers (NodeList) and material types (ElemType) of tetrahedra elements, and vertex positions (Xn) are read from the mesh file generated by GiD. Using the element types, the nodes inside the regions and on the interfaces are also set a type (NodeType). Then each node 5.3 Tissue Model Initialization 64 (c) A perspective view. Figure 5.1: Matlab interface running the simulation. is assigned a friction constant and a puncture threshold (Fd). Tissue parameters (TissuePar) define Young's modulus E and Poisson's ratio u of prostate and surrounding soft tissue. During the simulation, at many points the elements connected to a node are needed. There-fore, a node-element connectivity list (NEconn) is generated offline to avoid traversing the element list each time. Before starting the stiffness matrix creation, the nodes (Nfixed) which are constrained in space are defined. In the current version of the simulation, the entire volume is fixed in space only from a couple of nodes on the surface of the prostate above its apex (denoted by the colored patch on the prostate surface in Figure 5 .1). This is to represent the stiff tissue connecting it to the pubic arch and to display its effect on the torquing of the prostate as the needle is inserted. If the pubic arch were segmented, all of its nodes would have been used as fixed displacement boundary conditions since a major skeletal bone can be assumed a mechanical ground. Segmenting that other anatomy would make the simulation more realistic. 5.4 Needle Model Initialization 65 N a m e Dimension Content X n n x 3 Nominal positions of tissue nodes, (its order sets the global num-bering) NodeList e x 4 Indices of each element's corner nodes, (sets element numbering) ElemType e The tissue content of each element. Nodetype n The tissue type at each node. F d n x 2 Shaft friction and tip puncture forces for each node. TissuePar t x 2 Parameters (E, v) of each tissue type. NEconn n x v The indices of elements surrounding each node. Nfixed b The indices of the nodes with displacement boundary condition. Kinv 3n x 3n The inverse global stiffness matrix for tissue F E M . Table 5.1: Tissue F E A variables. One of the last offline initialization steps is the computation of the inverse global stiffness matrix K~l as in Algorithm 5.1. On Line 11, this huge matrix is inverted using the Matlab's inv function. A l g o r i t h m 5.1 Creating the inverse stiffness matrix offline. 1: function CREATElNvSTlFFMAT(TissuePar, ElemType, Xn , NodeList, Nfixed) 2: TissuePar(S, v) \u00E2\u0080\u0094> Lame constants (A, u) \u00E2\u0080\u0094> C > (A.12) and (A. l l ) 3: Kglobal (3-1) 6: Find Be > (A.18) and (A.24) 7: #12x12 <\" vBeTCBe > (A.30) 8: Add Ki2xi2 o n * \u00C2\u00B0 t n e corresponding elements of Kglobal > (A.32) 9: for all i \u00E2\u0082\u00AC Nfixed do 10: the rows/columns {3i \u00E2\u0080\u0094 2,3i \u00E2\u0080\u0094 1,3i} of Kglobal ^ hX3 > (A.33) 11: K'1 <- (Kglobal)' 5.4 Needle Model Initialization The needle model is also created before the simulation starts. For this, the variables in Table 5.2 are set. First the mesh is generated by setting node positions (NXn) and elements (NNodeList) in a loop, according to the mesh structure desired. A simple or a symmetric mesh can be chosen as presented in Section 3.5. For the final simulation, the symmetric structure is used. Then its stiffness matrix is calculated using the same function in Algorithm 5.1, the needle's base end is positionally constrained, and the inverse stiffness matrix (NKinv) is calculated. 5.5 An Iteration 66 N a m e Dimension Content N X n n x 3 Nominal positions of needle nodes. NNodeList e x 4 Indices of each needle element's corner nodes. NKinv 3n x 3n The inverse global stiffness matrix for needle F E M . Table 5.2: Needle F E A variables. 5.5 An Iteration The flowchart of a code iteration is given in Figure 5.2. Its functions are detailed as follows: L O C A T E N E E D L E S E T C O N S T R A I N T A N G L E S M O D I F Y M E S H x , 9, *p)segment > found geometrically 2: for all j G contact nodes do Find R3TN { j j ( j \~X 3 : ^TiTi-i * ~ ^TN,iy^TN,i-lJ \u00E2\u0080\u00A2 C H E C K N E W E L E M E N T I N T E R C E P T : If the element that the tip is in has no (working) nodes on the needle, then its closest corner is needed to be added on the needle after a remeshing stage. Otherwise the flow jumps to contact node servicing. \u00E2\u0080\u00A2 M O D I F Y M E S H : This implements the steps 1 to 3 in Section 2.5. First, the (approximate) nominal coordinate for the tip is found, then the boundary conditions are switched on K\"1 for the surrounding nodes. After that, the current surrounding elements are removed by subtracting their Ke and the new elements are added similarly. Finally, the boundary conditions are switched back. In the end, X n is updated with P' to reflect the coordinate change. \u00E2\u0080\u00A2 A D D N O D E : The new node j is added to the list of contact nodes and its variables (e.g., the rotation matrices j and - f t ^ T - i ) a r e initialized. Rows/columns corresponding to the contact nodes are copied from K^1 to the condensed system matrix KH. Finally, adjusting the L C F s using cumulative rotations RTN.I a n d switching the boundary conditions {x,y,z for stuck; only x and z for slip), the working matrix Kw is obtained. \u00E2\u0080\u00A2 S E R V I C E F I N I T E S T A T E M A C H I N E : If a working node slips off the tip, its corresponding rows/columns are removed from Kw- The stick-slip states are serviced as given in Al -gorithm 5.3. \u00E2\u0080\u00A2 S E T C O N S T R A I N T S : Each contact node's world-frame position is found using the orienta-tions and the knots of the segment that it resides on. For stuck nodes, corresponding displacements applied on three y rows/columns; for slipping nodes, two displacements and one force (friction threshold) are applied. Then, the working subsystem is solved and the forces/displacements are arranged into separate vectors. For each contact node, an axial friction force is defined according to the tissue material and the length of needle that this contact node is representing. In other words, the friction on 5.5 An Iteration 68 A l g o r i t h m 5.3 Finite state machine servicing for a contact node. 1: if stuck then 2: Update L C F for stuck using R3T,T, > (2.11) 3: if \"force along y\" > \"friction threshold\" then 4: Boundary condition change along y 5: Set its state to \"slip\" 6: if slip then 7: if TipMotionj ^ TipMotiom_i then 8: Boundary condition change along y 9: Update L C F for stuck using RJTT. t> (2.11) 10: Set its state to \"stuck\" 11: else 12: Update L C F for slip using R^,.T. > (2.19) the needle surface between this node and the next contact node (or the tip) is accumulated and applied as a friction force on this contact node. It is used both as a force boundary condition of a slipping node and a static friction threshold for a stuck node (such that the node will start slipping past that limit). \u00E2\u0080\u00A2 P R O J E C T F O R C E S : This first transforms the contact node forces into needle frame (using RTN,I) and then distributes them onto the needle knots linearly as in Section 2.8. \u00E2\u0080\u00A2 I N T E G R A T E F O R C E S : This sums up the forces and the torques of contact nodes on the base (see Section 2.9). This is to be applied by the haptic device on the user hand. If a valid model for tip bevel deflection is found, it will probably be in the form of a pure tip force. If this is the case, this force can be added onto the needle model at this point of the iteration cycle. \u00E2\u0080\u00A2 S O L V E N E E D L E : After distributing the knot forces (dividing them by 3), the needle model is iteratively solved (FEA) for the needle displacements using the applied forces as shown in Section 2.1. At each iteration a velocity gain is applied on the displacements such that Ui = Ui-i + kAu. For stability, this gain needs to be adjusted according to the Young's modulus and Poisson's ratio of the needle. This solution converges within 1% margin in less than 10 iterations for proper values and acceptable forces/deflections. Finally, the displacements and positions of all the nodes are calculated as mentioned in Section 2.7. Although it is not crucial to compute them at haptic refresh rate, the current implementation is also using some of the deformed node positions in various operations (e.g., penetration check in C H E C K N E W E L E M E N T I N T E R C E P T ) . Moving this step to a second (visualization) computer and its drawbacks are discussed in Section 5.7. 5.6 Notes 69 5.6 Notes A n advantage of the C - M E X file is that it can easily be compiled and loaded onto a target machine in the future without the need for many changes. Here target machine refers to a computer that can handle fast haptic loops while communicating with optical encoders and motor driving amplifiers at a fast and reliable rate by usually running a real-time operating system (e.g., V x W o r k s \u00C2\u00AE ) . That structure allows the use of WindRiver's Tornado\u00C2\u00AEdeve lopment environment and the control of runtime variables in Matlab S imul ink\u00C2\u00AEthrough the Real-Time Workshop\u00C2\u00AEtoo lbox , in which the author has experience with. For matrix operations, the functions such as matrix multiplication and rank 3-4 inversion are included in the code. However, F E A matrix operations were not optimized, i.e., although there are many zero elements in the matrices Be and C, they are actually multiplied element by element in the code. Skipping zero terms would could speed up the algorithms, but the current bottleneck is the 0(n2) operations (mesh adaptation), therefore all the other optimizations are given a lower priority. The presented (D(m2n) approach decreased the 0(mn2) mesh adaptation time of two sample meshes with 836 nodes and 240 nodes from \u00C2\u00AB 1 3 s and wl s down to \u00C2\u00AB 1 . 2 5 s and \u00C2\u00AB 0 . 0 2 s, respectively. The memory requirements for K~l storage are 50 M B and 4 M B , respectively. Figure 5.3 shows the needle base reaction force and the time requirement of each iteration during the insertion shown in Figure 5.1. The needle is modeled using a symmetric mesh (shown in Section 3.5) having 32 elements and it is pushed at a constant velocity of 0.1 mm/iteration. This would correspond to a 5 cm/s speed in a haptic implementation running at 500 Hz. In Fig-ure 5.3(b), note that the computations take 2-3 s during mesh modification using 0(mn2). The effect of the different stiffness and the friction of prostate tissue can be observed in Figure 5.3(a) around the 380 t / l iteration. 100 150 200 250 300 350 400 450 Iteration number (a) Needle base forces (from top - x, z, y). ?00 150 200 250 300 350 400 450 Iteration number (b) Time taken by each iteration. Figure 5.3: The simulation forces and timings. 5.7 Discussion 70 5.7 D i s c u s s i o n In this chapter, the implementation of extensions to DiMaio's work [40] have been described. One major change was extending the algorithms to 3D, which generated challenges in local frame handling and meshing (needle asymmetry discussed in Chapter 3). The other main improvement was the custom meshing. The data structures and algorithms were set to accommodate any number of nodes and elements having arbitrary connectivity (without dependency on maximum node valence) so that high-level meshing can be used in order to have a more accurate model with high-quality elements and a small number of nodes. Although the current code is not fast enough, employing a node snap (where the closest node is constrained on the needle without remeshing, which is also implemented in DiMaio's needle insertion simulator) instead of remeshing would let it run at a haptic rate (>200Hz). Although the needle flexibility causes problems due to force spikes, these issues can be attempted using force smoothing techniques. In order to make this implementation faster, two burdens, namely the mesh adaptation and visualization calculations, can be removed from the haptic loop together. They can be both moved to a separate computer (let us call it TRUS computer), which holds the K~1 matrix. When a node addition is required, the T R U S computer can send the condensed matrix (KH) to the haptic controller which processes it to obtain Kw- The remeshing can also be implemented such that the new node P' found by the haptic controller can be sent to the T R U S computer, which can then update the K~x matrix (and sends back the new condensed matrix). This structure has two flaws: the time lag of the remeshing and the node positions for new element interception check. The former one refers to the time passes on the haptic side between the element interception and the return of KH- Nonetheless, the haptic loop can continue handling the interaction and apply the changes when the matrix is available (maybe using some additional force smoothing). The latter refers to the need of the penetration check algorithm for the deformed positions of the surrounding nodes. Since there only a few of those nodes, their positions can be sent back by the TRUS computer. Chapter 6 Conclusion 6.1 Summary and Contributions A prostate brachytherapy simulator design has been presented. The design consists of two major subsystems, a needle-tissue interaction model and a T R U S image synthesizer. The specific contributions of this thesis can be summarized as follows: 3 D Extension: In this thesis, the concepts which were proposed in prior work [40] are extended to 3D. The presented system is the first physically-based 3D interaction model of a flexible needle insertion into a soft deformable body. Existing literature has presented heuristic approaches [20, 21], 2D F E M implementations [9, 40, 36], and a sample insertion of a rigid needle into nonlinear 3D F E M [97]. Nevertheless, the only physically-based 3D work by [97] utilized a non-conforming F E M , whereas the proposed method successfully couples a conforming quasi-static linear F E M of the tissue with the flexible needle. \"Effective\" Meshes: Using the presented scheme, a 3D volume can be meshed with conform-ing region boundaries on surfaces generated from contours on parallel slices. This method achieves relatively high-quality meshes which are also suitable for fast haptic implementa-tions of the F E M . Brachytherapy Simulation: A prostate mesh is generated using the above mentioned method. Incorporation of this mesh into the proposed interaction model has resulted in the first physically-based 3D needle insertion simulation of a brachytherapy procedure. The issues related to the haptic implementation have been addressed through this simulation which is currently running in a Matlab interface. T R U S Synthesis: The feasibility of TRUS generation as a visual feedback modality in the above mentioned simulation has been studied. To the best of our knowledge, the presented technique for ultrasound synthesis from a deformed volume with deformations defined by a finite element model is a novel contribution to this field. The speed considerations have also been addressed and a sample case of needle overlay has been presented. 71 6.2 Future Work 72 6.2 Future Work Segmenting more anatomical structures of the region and incorporating them into the simulator will increase the realism. Most importantly the pubic arch is needed as a fundamental boundary condition. The bladder is also important as a fluid pocket. The fatty and muscle tissue and the urethra can also be added to increase the realism. Moreover, a finer meshing of the tissue regions, where the needle passes frequently or seeds are implemented (namely P T V ) , might increase the simulation accuracy. On the other hand, the values which were set empirically for the presented mesh (e.g., preferred element sizes) can be studied for their effects on the final performance so that these parameters can be theoretically bounded by means of expected accuracy. The linear quasi-static model should be verified to see how well it represents the reaction forces and deformations of the prostate region. Before an in vivo validation, 3D deformations of tissue-like phantoms and in vitro tissue samples can be used. Depending on the performance of these deformation models, nonlinear and dynamic models may be attempted next. In particular dynamic characteristics may bring more realism to the simulator (e.g., the inertia of the tissue can then be felt). Other elasticity models may be attempted as well. The demonstrated iterative solution of the quadratic-strain needle model does not converge in some ranges of its elasticity parameters. Other solution methods which were shown to con-verge require large matrix inversions ( K \" 1 which is rank \u00C2\u00AB 5 0 - 1 0 0 ) at each iteration. Thus, before a haptic implementation, this may need to be addressed in order to prevent instability. The feasibility of other approaches (such as mass-spring models) for needle bending can also be investigated. The needle tip bevel effect should also be studied and incorporated into the simulator. Finding the bevel deflection through model fitting and parameter optimization using experimental insertions may be a promising approach. The T R U S synthesis rate has to be brought up to 12 Hz for real-time implementation. A l -though the current computation of T R U S generation is rather slow, there are numerous ap-proximation schemes that can be applied to improve performance. For instance, calculating deformation for one every 10 pixels in each direction and then interpolating the deformations of the rest would increase speed by 2 orders of magnitude and drop the computation time down to the order of seconds. A combination of different approaches can be utilized for interpolation. For example, different regions of the images can be interpolated using different methods. However, these patches might be hard to fuse realistically. The proposed simulator or a simplified version of it can be implemented on a haptic device in order to address possible implementation problems. Two computer or dual processor implemen-tations can be considered. Priority scheduling of VxWorks can also be utilized for background processing of low-priority tasks. The literature addresses similar problems and implementations extensively and a suitable solution can be extracted from these related studies. Bibliography [1] National Cancer Institute of Canada, \"Canadian cancer statistics 2004,\" 2004. http://www.cancer.ca. [2] W. M. Butler, \"Review of modern prostate brachytherapy,\" in Proceedings of the 22nd Annual EMBS International, (Chicago, IL), pp. 748-752, 2000. [3] R. Taschereau, J. Pouliot, J. Roy, and D. Tremblay, \"Seed misplacement and stabilizing needles in transperineal permanent prostate implants,\" Radiotherapy and Oncology, vol. 55, no. 1, pp. 59-63, 2000. [4] B. J. Davis, R. R. Kinnick, M. Fatemi, E. P. Lief, R. A. Robb, and J. F. Greenleaf, \"Measure-ment of the ultrasound backscatter signal from three seed types as a function of incidence an-gle: Application to permanent prostate brachytherapy,\" International Journal of Radiation On-cology* Biology* Physics, vol. 57, no. 4, pp. 1174-1182, 2003. [5] P. J. Gorman, T. M. Krummel, R. Webster, M. Smith, and D. Hutchens, \"A prototype haptic lumbar puncture simulator,\" in Proceedings of 8th Medicine Meets Virtual Reality (MMVR), pp. 263-274, 2000. [6] F. S. Azar, D. N. Metaxas, and M. D. Schnall, \"Methods for modeling and predicting mechan-ical deformations of the breast under external perturbations,\" in Medical Image Computing and Computer-Assisted Intervention (MICCAI), pp. 1267-1270, 2001. [7] S. Miller, C. Jeffrey, J. Bews, and W. Kinsner, \"Advances in the virtual reality interstitial brachyther-apy system,\" in IEEE Canadian Conference on Electrical and Computer Engineering, (Edmonton, Alberta), 1999. [8] J. K. Hahn, M. J. Manyak, G. Jin, D. Kim, J. Rewcastle, S. Kim, and R. J. Walsh, \"Cryother-apy simulator for localized prostate cancer,\" in Proceedings of ltfh Medicine Meets Virtual Reality (MMVR), 2002. [9] R. Alterovitz, K. Goldberg, J. Pouliot, R. Taschereau, and I.-C. Hsu, \"Sensorless planning for med-ical needle insertion procedures,\" in IEEE/RSJ International Conference on Intelligent Robots and Systems, (Las Vegas, NV), 2003. [10] G. Picinbono, H. Delingette, and N. Ayache, \"Non-linear and anisotropic elastic soft tissue mod-els for medical simulation,\" in Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), pp. 1370-1375, 2001. [11] Immersion, \"Endoscopy, laparoscopy, intra-vascular simulators,\" 2004. http://www.immersion.com. [12] Simbionix, \"Endoscopy, laparoscopy, intra-vascular simulators,\" 2004. http://www.simbionix.com. [13] Mentice, \"Laparoscopy simulator,\" 2004. http://www.mentice.com. [14] D. Terzopoulos and K. Waters, \"Physically-based facial modeling, analysis and animation,\" Visual-ization and Computer Animation, vol. 1, no. 2, pp. 73-80, 1990. [15] M. Downes, M. C. Cavusoglu, W. Gantert, L. W. Way, and F. Tendick, \"Virtual environments for training critical skills in laparoscopic surgery,\" in Proceedings of 6th Medicine Meets Virtual Reality (MMVR), (San Diego, CA), pp. 316-322, 1998. 73 BIBLIOGRAPHY 74 [16] S. P. DiMaio and S. E. Salcudean, \"Needle insertion modeling and simulation,\" IEEE Transactions on Robotics and Automation, vol. 19, no. 5, pp. 864-875, 2003. [17] L. Hiemenz, J. McDonald, D. Stredney, and D. Sessana, \"A physiologically valid simulator for training residents to perform an epidural block,\" in Proceedings of the 15th Southern Biomedical Engineering Conference, pp. 170-173, 1996. [18] G. Hirota, S. Fisher, and A. State, \"An improved finite element contact model for anatomical simu-lations,\" The Visual Computer, vol. 19, no. 5, pp. 291-309, 2003. [19] H. Fuchs, A. State, E. D. Pisano, W. F. Garrett, G. Hirota, M. Livingston, M. C. Whitton, and S. M. Pizer, \"Towards performing ultrasound-guided needle biopsies from within a head-mounted display,\" in Proceedings of Visualization in Biomedical Computing, (Hamburg, Germany), pp. 591-600, 1996. [20] A. N. Kimura, J. Camp, R. Robb, and B. Davis, \"A prostate brachytherapy training rehearsal system - simulation of deformable needle insertion,\" in Proceedings of Medical Image Computing and Computer-Assisted Intervention (MICCAI), pp. 264-271, 2002. [21] X. Wang and A. Fenster, \"A haptic-enhanced 3D real-time interactive needle insertion simulation for prostate brachytherapy,\" in Prooceedings of SPIE Medical Imaging 2004: Visualization, Image-Guided Procedures, and Display, vol. 5367, pp. 781-789, 2004. [22] C. M. Schneider, A. M. Okamura, and G. Fichtinger, \"A robotic system for transrectal needle insertion into the prostate with integrated ultrasound,\" in Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), (New Orleans, LA), pp. 365-370, 2004. [23] Z. Wei, G. Wan, L. Gardi, D. B. Downey, and A. Fenster, \"Robotic aided 3D TRUS guided intra-operative prostate brachytherapy,\" in the Proceedings of SPIE Medical Imaging 2004: Visualization, Image-Guided Procedures, and Display, vol. 5367, pp. 361-370, 2004. [24] S. F. F. Gibson and B. Mirtich, \"A survey of deformable modeling in computer graphics,\" Tech. Rep. TR-97-19, MERL - Mitsubishi Electric Research Laboratory, 1997. [25] S. F. F. Gibson, \"3D ChainMail: A fast algorithm for deforming volumetric objects,\" in Proceedings of 1997 Symposium on Interactive 3D Graphics, pp. 149-154, 1997. [26] K. Waters, \"A muscle model for animating three-dimensional facial expression,\" Computer Graphics, vol. 21, no. 4, pp. 17-24, 1987. [27] E. Promayon, P. Baconnier, and C. Puech, \"Physically-based deformations constrained in displace-ments and volume,\" Computer Graphics Forum (EUROGRAPHICS'96), vol. 15, no. 3, pp. 155-164, 1996. [28] D. L. James and D. K. Pai, \"Artdefo: Accurate real time deformable objects,\" in Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques, pp. 65-72, 1999. [29] D. L. James and D. K. Pai, \"Dyrt: Dynamic response textures for real time deformation simulation with graphics hardware,\" ACM Transactions on Graphics (SIGGRAPH Proceedings), vol. 21, no. 3, pp. 582-585, 2002. [30] D. Terzopoulos, J. Piatt, A. Barr, and K. Fleischer, \"Elastically deformable models,\" Computer Graphics, vol. 21, no. 4, pp. 205-214, 1987. [31] S. Cotin, H. Delingette, and N. Ayache, \"Real-time elastic deformations of soft tissues for surgery simulation,\" IEEE Transactions On Visualization and Computer Graphics, vol. 5, no. 1, pp. 62-73, 1999. [32] Y. Zhuang and J. Canny, \"Real-time simulation of physically realistic global deformation,\" in SIG-GRAPH99 Sketches and Applications, (Los Angeles, CA), 1999. [33] Y. Zhuang, Real-time simulation of physically realistic global deformations. PhD thesis, CS, University of California Berkeley, 2000. BIBLIOGRAPHY 75 G. Debunne, M. Desbrun, M.-P. Cani, and A. H. Barr, \"Dynamical real-time deformations using space & time adaptive sampling,\" ACM Transactions on Graphics (SIGGRAPH Proceedings), pp. 31-36, 2001. X. Wu, M. S. Downes, T. Goktekin, and F. Tendick, \"Adaptive nonlinear finite elements for de-formable body simulation using dynamic progressive meshes,\" Computer Graphics Forum (EURO-GRAPHICS'01), vol. 20, no. 3, 2001. H. -W. Nienhuys, Cutting in deformable objects. PhD thesis, Utrecht University, 2003. A. E. Kerdok, S. M. Cotin, M. P. Ottensmeyer, A. M. Galea, R. D. Howe, and S. L. Dawson, \"Truth cube: Establishing physical standards for soft tissue simulation,\" Medical Image Analysis, vol. 7, pp. 283-291, 2003. S. P. DiMaio and S. E. Salcudean, \"Simulated interactive needle insertion,\" in Proceedings of the 10th Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems (HAPTICS), pp. 344-351, 2002. J. Berkley, G. Turkiyyah, D. Berg, M. A. Ganter, and S. Weghorst, \"Real-time finite element modeling for surgery simulation: An application to virtual suturing,\" IEEE Transactions on Visualization and Computer Graphics, vol. 10, no. 3, pp. 314-325, 2004. S. DiMaio, Modelling, Simulation and Planning of Needle Motion in Soft Tissues. PhD thesis, E C E , University of British Columbia, 2003. M. Bro-Nielsen and S. Cotin, \"Real-time volumetric deformable models for surgery simulation using finite elements and condensation,\" Computer Graphics Forum (EUROGRAPHICS'96), vol. 15, no. 3, pp. 57-66, 1996. A. H. Gosline, S. E. Salcudean, and J. Yan, \"Haptic simulation of linear elastic media with fluid pockets,\" in Proceedings of the 12th Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems (HAPTICS), pp. 266-271, 2004. S. Cotin, H. Delingette, and N. Ayache, \"A hybrid elastic model for real-time cutting, deformations, and force feedback for surgery training and simulation,\" The Visual Computer, vol. 16, pp. 437-452, 2000. C. Simone and A. M. Okamura, \"Modeling of needle insertion forces for robot-assisted percutaneous therapy,\" in Proceedings of IEEE International Conference on Robotics and Automation, (Washing-ton, DC), 2002. H. Kataoka, T. Washio, K. Chinzei, K. Mizuhara, C. Simone, and A. M. Okamura, \"Measurement of the tip and friction force acting on a needle during penetration,\" in Medical Image Computing and Computer-Assisted Intervention (MICCAI), pp. 216-223, 2000. M. D. O'Leary, C. Simone, T. Washio, K. Yoshinaka, and A. M. Okamura, \"Robotic needle insertion: Effects of friction and needle geometry,\" in Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), (Taipei, Taiwan), pp. 1774-1780, 2003. A. M. Okamura, C. Simone, and M. D. OLeary, \"Force modeling for needle insertion into soft tissue,\" IEEE Transactions on Biomedical Engineering, vol. 51, no. 10, pp. 1707-1716, 2004. S. P. DiMaio and S. E. Salcudean, \"Needle steering and model-based trajectory planning,\" in Pro-ceedings of Medical Image Computing and Computer-Assisted Intervention (MICCAI), (Montreal, QC), 2003. R. J. Websterlll, N. J. Cowan, G. Chirikjian, and A. M. Okamura, \"Nonholonomic modeling of needle steering.\" to appear at 9th International Symposium on Experimental Robotics, June 2004. R. Alterovitz, K. Goldberg, and A. Okamura, \"Planning for steerable bevel-tip needle insertion through 2D soft tissue with obstacles.\" to appear at \"IEEE International Conference on Robotics and Automation (ICRA) 2005\", April 2005. BIBLIOGRAPHY 76 [51] B. Cutler, J. Dorsey, L. McMillan, M. Miiller, and R. Jagnow, \"A procedural approach to authoring solid models,\" ACM Transactions on Graphics (SIGGRAPH Proceedings), vol. 21, no. 3, pp. 302-311, 2002. [52] J. F. O'Brien, A. W. Bargteil, and J. K. Hodgins, \"Graphical modeling and animation of ductile fracture,\" ACM Transactions on Graphics (SIGGRAPH Proceedings), vol. 21, no. 3, pp. 291-294, 2002. [53] K. Dhinsa, C. Bailey, and K. Pericleous, \"Investigation into the performance of turbulence mod-els for fluid flow and heat transfer phenomena in electronic applications,\" in 2,0th Annual IEEE Semiconductor Thermal Measurement and Management (SEMI-THERM), pp. 278-285, 2004. [54] C. Singh Verma and V. C. V. Rao, \"Parallelization of finite volume computations for heat transfer application using unstructured mesh partitioning algorithms,\" in 4th International Conference on High Performance Computing, 1997. [55] F. Ganovelli, P. Cignoni, C. Montani, and R. Scopigno, \"Enabling cuts on multiresolution represen-tation,\" The Visual Computer, vol. 17, no. 5, pp. 274-286, 2001. [56] E. Grinspun, P. Krysl, and P. Schroder, \"Charms: a simple framework for adaptive simulation,\" ACM Transactions on Graphics (SIGGRAPH Proceedings), vol. 21, no. 3, pp. 281-290, 2002. [57] M. Miiller and M. Teschner, \"Volumetric meshes for real-time medical simulations,\" in Proceedings of 7th Workshop on Bildverarbeitung fur die Medizin, (Erlangen, Germany), pp. 279-283, Springer-Verlag, 2003. [58] P. Fleischmann, R. Kosik, and S. Selberherr, \"Simple mesh examples to illustrate specific finite element mesh requirements,\" in 8th International Meshing Roundtable, pp. 241-246, 1999. [59] R. V. Garimella and M. S. Shephard, \"Boundary layer meshing for viscous flows in complex domains,\" in 7th International Meshing Roundtable, pp. 107-118, 1998. [60] M. A. Yerry and M. S. Shephard, \"Automatic three-dimensional mesh generation by the modified octree technique,\" International Journal For Numerical Methods in Engineering, no. 20, pp. 1965-1990, 1984. [61] R. Lohner and P. Parikh, \"Three-dimensional grid generation by the advancing front method,\" In-ternational Journal for Numerical Methods in Fluids, vol. 8, pp. 1135-1149, 1988. [62] S. H. Lo, \"Volume discretization into tetrahedra - II: 3D triangulation by advancing front approach,\" Computers and Structures, vol. 39, no. 5, pp. 501-511, 1991. [63] B. N. Delaunay, \"Sur la Sphere - vide, izvestia akademia nauk sssr,\" VII Seria, Otdelenie Matem-aticheskii i Estestvennyka Nauk, vol. 7, pp. 793-800, 1934. [64] D. F. Watson, \"Computing the delaunay tesselation with application to voronoi polytopes,\" The Computer Journal, vol. 24, no. 2, pp. 167-172, 1981. [65] J. R. Shewchuk, \"Tetrahedral mesh generation by delaunay refinement,\" in Symposium on Compu-tational Geometry (SCG), (Minneapolis, MN), pp. 86-95, 1998. [66] J. R. Shewchuk, \"Triangle: Engineering a 2D quality mesh generator and delaunay triangulator,\" 1996. http://www.cs.cmu.edu/~quake/triangle.html. [67] H. Borouchaki, F. Hecht, E . Saltel, and P. L. George, \"Reasonably efficient delaunay based mesh generator in three dimensions,\" in 4th International Meshing Roundtable, pp. 3-14, 1995. [68] P. L. George, F. Hecht, and E. Saltel, \"Automatic mesh generator with specified boundary,\" Computer Methods in Applied Mechanics and Engineering, vol. 92, pp. 269-288, 1991. [69] M. Bern and P. Plassmann, \"Mesh generation,\" in Handbook of Computational Geometry (J.-R. Sack and J. Urrutia, eds.), Elsevier Science, 2000. BIBLIOGRAPHY 77 [70] S. J. Owen, \"A survey of unstructured mesh generation technology,\" in 7th Inter-national Meshing Roundtable, (Dearborn, MI), pp. 239-267, 1998. also available at http://www. andrew. emu. edu/user/sowen /survey / . [71] N. Molino, R. Bridson, J. Teran, and R. Fedkiw, \"A crystalline, red green strategy for meshing highly deformable objects with tetrahedra,\" in 12th International Meshing Roundtable, pp. 103-114, 2003. [72] T. R. Nelson, D. B. Downey, D. H. Pretorious, and A. Fenster, Three Dimensional Ultrasound. Lippincott Williams & Wilkins, 1999. US. [73] R. N. Rohling, A. H. Gee, and L. Berman, \"A comparison of freehand three-dimensional ultrasound reconstruction techniques,\" Medical Image Analysis, vol. 3, no. 4, pp. 339-359, 1999. [74] P. R. Detmer, G. Bashein, and T. Hodges, \"3-D ultrasonic image feature location based on magnetic scan head tracking: In vitro calibration and validation/' Ultrasound in Medicine and Biology, vol. 20, pp. 923-936, 1994. [75] R. N. Rohling, A. H. Gee, and L. Berman, \"Three-dimensional spatial compounding of ultrasound images,\" Medical Image Analysis, vol. 1, no. 3, pp. 177-193, 1997. [76] D. F. Leotta and R. W. Martin, \"Three-dimensional ultrasound imaging of the rotator cuff: Spatial compounding and tendon thickness measurement,\" Ultrasound in Medicine and Biology, vol. 25, no. 4, pp. 509-525, 2000. [77] R. Franke, \"Scattered data interpolation: Tests of some methods,\" Math Computation, vol. 38, no. 157, pp. 181-200, 2000. [78] J. M. Thijssen and B. J. Oosterveld, \"Speckle and texture in echography: Artifact or information?,\" in Proceedings of Ultrasonics Symposium, pp. 803-809, 1986. [79] C. D. Barry, C. P. Allott, N. W. John, P. M. Mellor, P. A. Arundel, D. S. Thomson, and J. C. Wa-terton, \"Scattered data interpolation: Tests of some methods,\" Ultrasound in Medicine and Biology, vol. 23, no. 8, pp. 1209-1224, 1997. [80] R. N. Rohling, A. H. Gee, L. Berman, and G. Treece, \"Radial basis function interpolation for freehand 3D ultrasound,\" in Proceedings of the International Conference on Image Processing in Medical Imaging, (Visegrad, Hungary), 1999. [81] T. Roxborough and G. Nielson, \"Tetrahedron based, least square, progressive volume models with applications to freehand ultrasound data,\" in Proceedings of IEEE Visualization, pp. 93-100, 2000. [82] S. Meairs, J. Beyer, and M. Hennerici, \"Reconstruction and visualization of irregularly sampled three-and four-dimensional ultrasound data for cerebrovascular applications,\" Ultrasound in Medicine and Biology, vol. 26, no. 2, pp. 263-272, 2000. [83] H. Watanabe, H. Kaiho, M. Tanaka, and Y. Tersawa, \"Diagnostic application of ultrasonography to the prostate,\" Invest. Ural, vol. 8, pp. 548-549, 1971. [84] S. Tong, D. B. Downey, H. N. Cardinal, and A. Fenster, \"A three-dimensional ultrasound prostate imaging system,\" Ultrasound in Medicine and Biology, vol. 22, pp. 735-746, 1996. [85] R. W. Prager, A. Gee, and L. Berman, \"Stradx: Real-time acquisition and visualisation of freehand 3D ultrasound,\" Medical Image Analysis, vol. 3, no. 2, pp. 129-140, 1999. [86] D. Aiger and D. Cohen-Or, \"Real-time ultrasound imaging simulation,\" Real-Time Imaging, vol. 4, no. 4, pp. 263-274, 1998. [87] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method. Butterworth Heinemann, 2000. [88] T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Mineola, New York: Dover Publications, 2000. [89] M. Bro-Nielsen, \"Finite element modeling in medical VR,\" Journal of the IEEE, vol. 86, no. 3, pp. 490-503, 1998. BIBLIOGRAPHY 78 [90] I. Wolfram Research, \"Mathworld.\" http://mathworld.wolfram.com/WoodburyFormula.html. [91] E. Keppel, \"Approximating complex surfaces by triangulation of contour lines,\" IBM Journal of Research Developments, vol. 19, pp. 2-11, 1975. [92] J.-D. Boissonnat, \"Shape reconstruction from planar cross-sections,\" Computer Vision, Graphics and Image Processing, vol. 44, pp. 1-29, 1988. [93] Institut National de Recherche en Informatique et en Automatique (INRIA), \"Nuages,\" 2000-2004. http: / / www-sop.inria.fr / prisme/logiciel / nuages.html. [94] E. J. Nicolson, Tactile sensing and control of a planar manipulator. PhD thesis, EECS, University of California Berkeley, 1987. [95] S. Hang, TetGen, a quality tetrahedral mesh generator and three-dimensional delaunay triangulator, vl.3 user's manual. Weierstrass Institute for Applied Analysis and Stochastics (WIAS), Berlin, Germany, 9 ed., 2004. http://tetgen.berlios.de/. [96] International Center for Numerical Methods in Engineering (CIMNE), \"Gid: The personal pre and postprocessor,\" 2000-2004. http://gid.cimne.upc.es/. [97] H.-W. Nienhuys and A. F. van der Stappen, \"Interactive needle insertions in 3D nonlinear material,\" Tech. Rep. 2003-019, UU-CS: Utrecht University - Information and Computing Sciences, Utrecht, 2003. Appendix A Finite Element Method A . l Definitions of Stress and Strain Consider a small area element AA, oriented with unit normal vector n, on the surface or within a deformable body under external forces as shown in Figure A . l . Let A F and A M be the resultant force and moment vectors, respectively, on A A . Their intensities on the area can be found at the limit as follows: i - A , A M _ dM _ ft ^ ' l i m A A ^ o \" A X ~~ HA ~ where T is called stress vector or traction, and C is called couple-stress vector. The elementary theory of elasticity proceeds on the assumption that C = 0 [87]. T is the stress intensity at a point with respect to a particular area element specified by the unit normal n. A complete description of stress at a point requires the traction for all directions. This can be found by traction in principle directions. Figure A . l : Resultant force and moment vectors on a small area element within a deformable body. Consider an infinitesimal cube centered at a point inside a deformable body. Let the faces of the cube be oriented with the orthonormal basis vectors x, y, z of the coordinate frame as 79 A.l Definitions of Stress and Strain 80 shown in Figure A.2. Then the tractions on the faces with surface normals along the principle directions are shown as Tx,Ty,Tz. If they are written in terms of their Cartesian coordinates, e.g., Tx = oxxx + oxyy + oxzz, the scalar stress components form a second-order tensor ^ defining the state of stress at a point as: * = 0~xx 0~xy O^x (A.2) where the traction of an area element with unit normal n is given by \17 h X Figure A.2: Components of stress on an infinitesimal cube. Under the assumption that couple-stress does not exist, the equilibrium of zero moment implies Oij = oji. Therefore the symmetric ^ has only 6 independent variables, which can be expressed in the vector representation of stress as: 0~xx 0~yy 0~zz \u00C2\u00AE~yz 0~zx 0~xy (A.3) where the first three components are called normal stresses, and the other three are called shear stresses. Strain characterizes the deformation of a body under stress. Similar to the stress tensor a second-order strain tensor can completely describe the state of strain at a point: T = (A.4) A.2 Deformation 81 A more complicated mathematical treatment [87] shows the symmetry of T for elastic bodies. Thus, \u00E2\u0082\u00ACij = eji, and the vector representation of strain becomes: e = (A.5) A . 2 Deformation Displacements of a point x = (x, y, z) can be written as u \u00E2\u0080\u0094 [u(x) v(x) w(x)]T in the reference coordinate frame of the undeformed body. Then the strain components are given by: _ du 1 S x x ~ dx' + 2 du du dx + dv dx _ dv \u00E2\u0082\u00ACxy ~ ~7\ h ~x h dy dx du du dv dv dw dw dx dy dx dy dx dy (A.6) (A.7) The other four strain components are defined similarly. Note that the square-bracketed second-order terms are relatively small for small displacements. The matrix expression of the engineering strain tensor in terms of u becomes: \u00E2\u0082\u00AC(\u00C2\u00AB) = eyy eyz \u00E2\u0082\u00ACzx exy dx 0 0 0 d_ dz d_ dy d_ dy d_ dz du , 1 dx ^ 2 dv_ _i_ 1 dy ^ 2 dw i 1 dz ^ 2 dv I dw_ dz dy du I djv dz dx du _,_ dji dy dx du\^ I I dv_\^ I i dw \* dx' ~*~ V dx ) \dx> / 9 u \ 2 , / dv> \dy> \"rAdy> ^ du d-z) + dz 0 0 d_ dz d_ dy d_ dx 0 \u00C2\u00A3 d_ dx Bu + -A\u00C2\u00AE - ) 2 dz > dw dw du du I 9\u00C2\u00BB3\u00C2\u00BB j dy dz dy dz y z du du I dv dv_ , dw_ dw dz dx dz dx dz dx du du i dv_ dv_ , dw dw dx dy dx dy dx dy U V W + du dx dv dx dw dx 0 0 0 0 0 0 0 0 0 du dy dv dy dw dy 0 0 0 0 0 0 0 0 0 du dz dv dz dw dz 0 0 0 du dz dv dz dw dz du dy dv dy dw dy du dz dv dz dw dz 0 0 0 du dx dv dx dw dx du . dy dv dy dw dy du dx dv dx dw dx 0 0 0 du dx dv dx dw dx du dy dv dy dw dy du dz dv dz dw dz (A.8) This vector cr(u) is called the Green-Lagrange strain tensor. Neglecting the second-order A.3 Linear Elastic Behaviour 82 terms under small deformations leads to the linear strain approximation as follows: Ui= - ^ - [ f c i i u i + . . . + 0\" + . . . + fcirur-/i] (B.3) ^ j e H i : fi = kj\u\ + . . . + 0Jl + . . . + kjrur k-- J 1 [kilUi + ...+0U + ...+ kirUr - fi] To be able to write it as a difference from the initial K~x matrix like in (2.8) we add and subtract the term kjifi and reorder the terms: V j e N ; : fj = kjiUi + . . . + kjifi + \u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2 + kjrUr k-~r~~ [fca\u00C2\u00ABi + \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 + (ku -l)fi + ... + kirur] (B.4) To make the U{ in (B.3) also look like the form in (2.7) we should add and subtract all the 1The terms 0\" and 0J' are used to indicate the absence of the terms with the factors ku and kji, respectively. 89 90 missing terms and regroup them properly: Uj = J = ki\U\ + . . . + kafi + ... + kiruT \u00E2\u0080\u0094 \ki\U\ + . . . + kafi + ... + kirur] 1 1 \u00E2\u0080\u0094 7\u00E2\u0080\u0094 [knui + ... + 011 + . . . + kirur] + \u00E2\u0080\u0094fi hum + ... + kafi + ... + k irur kii + 1 1 - k-2 [knUi + ...+0U + ... + kirUr] + J-^-fi k%\U\ + . . . + kafi + ... + k{rur ka H~ 1 kj' [knu-i + ... + (ku - + \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 + kirur] (B.5) In (B.4) and (B.5) the upper line terms are the initial coefficients of the K 1 matrix. We can rewrite the common square bracketed terms in matrix form as: hiUi + \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 + (ka ~ l ) / i + \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 + kirur \u00E2\u0080\u0094 \kn \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 (ka ~ 1) \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 h (B-6) And the coefficients of these bracketed terms become: hi ka I T hi kii 1 ka + 1 (B.7) kf kj-"@en . "Thesis/Dissertation"@en . "2004-11"@en . "10.14288/1.0065481"@en . "eng"@en . "Electrical and Computer Engineering"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en . "Graduate"@en . "Ultrasound image and 3D finite element based tissue deformation simulator for prostate brachytherapy"@en . "Text"@en . "http://hdl.handle.net/2429/16194"@en .