UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Ultrasound image and 3D finite element based tissue deformation simulator for prostate brachytherapy Goksel, Orcun 2004

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

Item Metadata


831-ubc_2005-0045.pdf [ 14.01MB ]
JSON: 831-1.0065481.json
JSON-LD: 831-1.0065481-ld.json
RDF/XML (Pretty): 831-1.0065481-rdf.xml
RDF/JSON: 831-1.0065481-rdf.json
Turtle: 831-1.0065481-turtle.txt
N-Triples: 831-1.0065481-rdf-ntriples.txt
Original Record: 831-1.0065481-source.json
Full Text

Full Text

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 REQUIREMENTS FOR T H E DEGREE OF MASTER OF APPLIED  SCIENCE  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 UNIVERSITY  O F BRITISH C O L U M B I A  December 2004 ©  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 T R U S 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  T a b l e of C o n t e n t s  iii  List of Tables  vi  List of F i g u r e s  vii  Glossary  ix  Acknowledgements  xi  1  2  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  Needle-Tissue Interaction M o d e l  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  23  Stick State iii  TABLE  3  4  5  OF CONTENTS  iv  2.4.2  24  Slip State  2.5  Decision of Remeshing  27  2.6  Mesh Modification  29  2.6.1  0(mn )  Approach  30  2.6.2  0(m n)  Approach  31  2  2  2.7  Needle Model  32  2.8  Projecting Forces onto the Needle  32  2.9  Integrating Needle Forces  33  2.10 Discussion  34  M e s h 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  U l t r a s o u n d 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  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 6  A  B  v  Conclusion  71  6.1  Summary and Contributions  71  6.2  Future Work  72  Bibliography  73  Finite Element M e t h o d  79  A.l A.2 A.3 A.4 A.5  79 81 82 83 86  Definitions of Stress and Strain Deformation Linear Elastic Behaviour Discretization Element energy  Boundary Condition Change  89  L i s t  o f  T a b l e s  3.1  M e s h quality  43  5.1  Tissue F E A variables  65  5.2  Needle F E A variables  66  vi  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  M R I 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  viii  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  A.3 A 2-D discretization example of a body Q, into finite elements Q,  80 1  A.4 A tetrahedral element and its corners in local numbering  84 84  G l o s s a r y  Roll-pitch-yaw angles 2D  Two-dimension(al)  3D  Three-dimension(al)  BCC  Body Centered Cubic  BCCA  British Columbia Cancer Agency  BEM  Boundary Element Method  bNED  biochemical no evidence of disease  CDT  Constrained Delaunay Tetrahedralization  CG  Conjugate Gradients  CT  Computed Tomography  D I C O M - R T Digital Imaging and Communications in Medicine - Radiation Therapy DOF  Degree of Freedom  DW  Distance Weighted  FEA  Finite Element Analysis  FEM  Finite Element Method  HDR  High-dose-rate (brachytherapy)  MRI  Magnetic Resonance Imaging  NURBS  Nonuniform Rational B-Splines  PNN  Pixel Nearest Neighbour  PSA  Prostate Specific Antigen  GLOSSARY  x  PTV  Planning target volume  RBF  Radial Basis Functions  TRUS  Transrectal Ultrasound  US  Ultrasound  VNN  Voxel Nearest Neighbour  VR  Virtual-Reality  x  Mixed force-displacement outputs of working subsystem  y  Mixed force-displacement inputs of working subsystem  /  Force vector  K~  Inverse global stiffness matrix  KH  Condensed system matrix  Kw  Working matrix  K'  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  l  w  Acknowledgements • First of all, I would like to thank my supervisors T i m 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. • 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, X u , Sara, Hassan, Reza, ... • 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. A n d , without the help of the friendly staff of the B C Cancer Agency, the data collection would not be easy at all. • I appreciate the logistic and psychological support of my fellow Johannians through the process of (especially writing of) this thesis. • I definitely owe my being here to my parents whom I should acknowledge the most. Their endless support has always been with me.  xi  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; chemotherapy, 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 tissue. 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 needle. 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 continence 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 G l e a s o n ^ ) [2]. 1  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. However, 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 experienced 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 application 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 1  Gleason score is a cancer aggressiveness grade determined via biopsy.  1.2 Brachytherapy  3  disease a n d / o r w i t h neoadjuvant h o r m o n a l therapy to downsize the gland a n d i n h i b i t the spread. For p l a n n i n g the i n d i v i d u a 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  w i t h the patient i n the dorsal l i t h o t o m y p o s i t i o n (lying on his back). A sketch of the procedure setup is shown i n F i g u r e 1.2(c). A transrectal u l t r a s o u n d ( T R U S ) probe w h i c h is positioned using a stepper unit mounted on a stabilizer captures images at 5 m m increments (see F i g u r e 1.2(b)). T h e template g r i d (which is described below) attached to the stepper also discretizes this volume transversally at 5 m m intervals.  .Stylet  TRUS * probe  Stabilizer  (a) The operation room.  (b) The TRUS setup.  (c) A sketch of the procedure.  Figure 1.2: Prostate brachytherapy operation. For p l a n n i n g , the ultrasound prostate volume is expanded w i t h margins to a p l a n n i n g target volume ( P T V ) , w h i c h encompasses likely extracapsular disease. T h e decision of P T V seed arrangement is made by a r a d i a t i o n oncologist using special software to evaluate parameters, e.g., dose d i s t r i b u t i o n . P r i o r to the procedure, the sterilized seeds (see "standard seeds" i n F i g ure 1.3) are loaded into needles.  T h e B C Cancer A g e n c y ( B C C A ) uses O n c o S e e d ™ that is a  welded t i t a n i u m capsule containing radioactive iodine isotope  1 2 5  I (with a half-life of 60 days)  absorbed onto a silver r o d .  ' H III flSH Corrugated seeds I—— Seeds In vicryl sutures  j  Figure 1.3: R a d i o a c t i v e seeds i m p l a n t e d d u r i n g brachytherapy (reprinted from [4]). Fluoroscopy is a m e d i c a l imaging modality, w h i c h is also utilized i n brachytherapy, to view  1.2 Brachytherapy  4  the anatomy by means of X - r a y . T h e operation r o o m w i t h the a C - a r m (for fluoroscopy) a n d U S setup can be seen i n F i g u r e 1.2(a).  A l t h o u g h the C - a r m can rotate for different i m a g i n g  directions, fluoro-images are o n l y taken along the posterior-anterior axis i n brachytherapy. T h e intra-operative T R U S probe is set up i n the same way as i n the p r e i m p l a n t volume study so that the prostate can be registered to the i m p l a n t p l a n i n the b e g i n n i n g of the procedure.  Loaded  needles are inserted v i a a template g r i d (see F i g u r e 1.2(c)) the coordinates of which are given by the p r e p l a n . Needle p o s i t i o n i n g is confirmed using the T R U S (see F i g u r e 1.4) a n d fluoroscopy (see F i g u r e 1.5).  For guidance, w h i t e dots denoting t e m p l a t e g r i d holes are overlayed onto  the u l t r a s o u n d images by the T R U S hardware as i n F i g u r e 1.4(b).  Since there is a risk of X -  ray overdose, the number of fluoro-images taken d u r i n g a procedure must be kept low (20-50). Consequently, T R U S becomes the only m o d a l i t y s u p p l y i n g continuous feedback to the physician. Before the procedure the patient is anaesthetized .  If anaesthesia is not safe, the s p i n a l  blockade is applied so that the patient is t e m p o r a r i l y p a r a l y z e d below waist. F i r s t the p h y s i c i a n fixes the stabilizer a n d the precision stepper w h i c h are used to p o s i t i o n the u l t r a s o u n d probe. T h e stepper is calibrated by m a t c h i n g the T R U S images w i t h the ones that were taken i n advance. These m a n u a l adjustments are done by visually c o m p a r i n g the p r i n t e d images w i t h the i n t r a operative T R U S imaging. It takes relatively long time a n d it is subject to the experience of the physician. T h e n the fluoroscope is quickly adjusted a n d the radioactive seeds are taken into the operation r o o m . T h e brachytherapy needles are loaded w i t h up to six seeds each prior to the procedure. If a l l seeds i n a needle are spaced by 1 c m , then a special type of seed packaging (the coated ones at the b o t t o m i n F i g u r e 1.3) can be used. T h e B C C A uses the b r a n d R A P I D S t r a n d ™ - R i g i d A b s o r b a b l e Permanent Implant Device). T h i s manufacturer w r a p p e d strands are preferable since they are less likely to j a m i n the needle. T h e y also stay at the i m p l a n t e d location a n d o r i e n t a t i o n whereas the i n d i v i d u a l seeds may rotate d u r i n g r e t r a c t i o n or can travel t h r o u g h b l o o d vessels. C u s t o m spacing is achieved by a d d i n g spacer between i n d i v i d u a l seeds. T h e needles are inserted to their targets t h r o u g h a g r i d template. T h e g r i d hole a n d insertion depth are double-checked b y a nurse who is reading the therapy p l a n to the doctor. T h e stepper may need to be readjusted  (a) Needle and a T R U S image. Figure 1.4:  (b) a B-scan in the beginning.  (c) a B-scan in the end.  A representation of the 3 D volume w i t h T R U S layers a n d a needle (a), a n d two B-scans from a procedure.  1.3 Need for  5  Training  (a) Prior to the procedure  (b) After the 5  th  needle  (c) Towards the end of the procedure  Figure 1.5: F l u o r o s c o p y images from a brachytherapy procedure.  d u r i n g the procedure because of the prostate swelling. T h e doctor may also modify the procedure p l a n slightly a n d / o r use e x t r a seeds at the end of the operation to correct a n d improve the dose distribution. D u r i n g a procedure wlOO seeds are i m p l a n t e d using 20 to 25 needles. T h e needles are 20 c m i n length and quite flexible w i t h a significant t i p bevel. A needle consists of an outer  cannula  (that holds the seeds) w i t h black stripes m a r k i n g every centimeter and a stylet that t i g h t l y fits i n this c a n n u l a (see F i g u r e 1.6). So, the seeds can be released b y p u l l i n g the needle cannula w h i l e 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 i m p l a n t and the queueing of the first seed to the needle t i p , respectively. A n arrow on the cannula base marks the direction of t i p bevel. U s i n g this bevel and the needle flexibility,  the physician can steer the needle. T h i s can be done to refine the grid discretization  b y reaching in-between template holes, i f the p h y s i c i a n wishes to make slight modifications to the preplan.  Moreover, a direct insertion cannot reach b e h i n d the p u b i c arch where steering  becomes necessary. P r o p e r steering places seeds away from certain a n a t o m i c a l landmarks, e.g. the urethra a n d the rectum, while not overdosing or underdosing a region. However, this requires the experience of m a n y procedures.  Successful i m p l a n t s rely u p o n not o n l y a thorough comprehension of the  relative 3 D seed locations a n d surrounding a n a t o m i c a l structures, but also t a k i n g their deformat i o n d u r i n g needle insertion into account. These challenges reveal the need and importance of brachytherapy training.  1.3  Need for Training  T h e lack of direct visual feedback makes the comprehension of the needle p a t h w i t h respect to the a n a t o m i c a l layout rely only u p o n the medical i m a g i n g modalities, ultrasound and  fluoroscopy,  _  1.3 Need for Training  6  <n—i  11::::::::  Tip bevel  Cannula (shaft)  V»  Needle (cannula) base i  stylet  "Stylet /  holder / I  Figure 1.6: T h e brachytherapy needle.  a n d the kinesthetic feedback o n the needle base.  T h e fluoroscopy gives a clear view of seeds  a n d the T R U S i m a g i n g supplies continuous information about the region. In a d d i t i o n to those imaging modalities, h a p t i c  2  feedback also gives important information to the physician about  w h i c h m e d i u m the t i p is i n , w h e n the prostate membrane is punctured, or whether the pubic arch is hit. However, correctly processing the images, the kinesthetic feedback, and their combined information i n terms of the anatomy of interest requires t r a i n i n g . A s mentioned before, there is a l i m i t on the number of fluoroscopy images that can be taken d u r i n g brachytherapy because of the X - r a y r a d i a t i o n . A s i d e from that l i m i t a t i o n , fluoroscopy does not convey a l l 3 D geometry (e.g., the height info) either. O n the other hand, unfortunately, the prostate, the pathology, the seeds, a n d the needle are not always clearly visible i n T R U S images. Furthermore, the procedure must be completed i n a relatively short period of time to m i n i m i z e edema and anesthesia side-effects, while the best possible seed d i s t r i b u t i o n may require intra-operative p l a n m o d i f i c a t i o n w h i c h also adds up to the complexity. In order to develop a l l those skills, m e d i c a l residents use mannequins, animals, or cadavers, w h i c h are either not very realistic or may involve ethical issues. T h e t r a i n i n g of a resident also involves observation and practice on real patients, w h i c h may lead to irreversible consequences even t h o u g h procedures are supervised by an experienced physician. T h e identification of the need for better t r a i n i n g schemes led us to consider alternative noninvasive solutions. Based on this m o t i v a t i o n , this thesis proposes a simulator of the brachytherapy procedure, on w h i c h the m e d i c a l residents can t r a i n and the expert physicians can excel their skills while experimenting new techniques. A s outlined above, a brachytherapy simulator should convey the haptic and imaging feedback to the user. T h u s , a tissue m o d e l interacting w i t h a needle, w h i c h is m a n i p u l a t e d by the user holding on its base, is needed so that the reaction forces of this m o d e l can be applied back Haptic comes from a Greek term for touch. Haptics is the science concerned with the tactile or kinesthetic sense of touch. 2  1.4 Related Work  7  onto the user's hand via the base. The imaging feedback may be ultrasound and/or fluoroscopy. Due to its more frequent referral by the physician during the procedure, ultrasound has been incorporated first in our simulator. Its continuous representation of tissue deformation was also coherent with our goal of tissue deformation modeling. For realistic T R U S feedback, images need to be generated according to the current placement of the tissue and the needle. This required a model for the deformed tissue image generation and needle representation on it.  1.4  Related Work  Just as a flight simulator is used for training pilots and it helps in identifying the complications of the current and proposed procedures, a medical simulator is built to recreate the realm of a surgery or examination to be used in training and testing. In this area there exists significant amount of previous work including epidural lumbar puncture [5], breast biopsy [6], interstitial brachytherapy [7], prostate cryotherapy [8], prostate brachytherapy [9], hepatic surgery [10]. There are also current commercialized products on laparoscopy [11, 12, 13], endoscopy [11, 12], and intra-vascular access [11, 12]. Our proposed needle insertion simulator for prostate brachytherapy includes a tissue deformation model. Generating a tissue mesh for this model and the interaction of the model with a brachytherapy needle are two issues that are addressed in this thesis.  Finally, ultrasound  synthesis as a visual feedback is discussed. Some related studies in those areas are mentioned throughout the rest of this chapter. The realism of a medical simulator is drastically affected by the tissue model of interest. Although the interconnection of muscle fibers, veins, nerves, etc. makes the prediction and generalization of tissue behaviour almost impossible, still some major characteristics (e.g., deformation) can be approximated for certain cases with available techniques [9, 10, 14, 15, 16]. For the simulation of these characteristics first a tissue model of the region of interest is needed. Section 1.4.1 presents some of them. Earlier attempts for haptics in medicine were simple approaches such as look-up tables [17]. In the course of time technological advances have enabled the use of haptics for physically-based interaction in medical simulations with complex control loops requiring fast update rates. Haptic implementations including palpation, incision, and needle insertion require a stable and realistic contact handling which has been a challenge in surgical tool simulations.  Hirota et al. [18]  recently proposed a contact method for finite element analysis ( F E A ) of non-linearly elastic objects and applied it to a knee bending simulation. Friction models for tissue interaction is another concern of medical simulation. Section 1.4.2 mentions methods that were applied in needle-tissue interaction. Owing to vision being our primary sense, most medical simulators have focused on graphics. Many created a virtual-reality (VR) environment that renders the region of interest from the surgeon's eyes [10]. Fuchs et al. [19] overlayed oriented ultrasound images onto the patient using  1.4 Related Work  8  V R goggles during US examinations.  Some researchers preferred rendering the 3D anatomy  including the obscured internal structures in their simulations [20].  The implementation and  realism of rendering methods were largely dealt in computer graphics. However, during a needle insertion nothing more than the tissue surface can be seen. This makes it completely dependent on medical imaging techniques such as ultrasound. Hahn et al. [8] developed a graphical interface for cryotherapy without haptics. They cut the desired US slices from a 3D volume data acquired by a 3D US device and overlayed the freezing probes as white dots. Kimura et al. [20] designed a brachytherapy simulation with US, fluoroscopy, and a 3D computer graphics rendering of the region. They added the needle as a black line on fluoro image and as a white disk on US. They also supplied a simple haptic feedback, a peak for puncture. Unlike the above mentioned ones, Wang and Fenster [21] proposed a modelbased approach to brachytherapy simulation. For a 3D graphical rendering of the region they used a modified version of the 3D ChainMail method described in the next section. They loosely coupled this deformation model with haptic feedback of only the prostate. Some recent research on improving brachytherapy procedures are by Schneider et al. [22] and Wei et al. [23]. The former proposes a new technique for inserting the needles through the rectum via a partial sheath surrounding the T R U S probe in order to maximize the seed positioning and the comfort to the patient.  The latter introduces a robot arm guiding the orientation of the  template in order allow oblique insertions. A major part of their work focused on finding the needle on US volume data gathered by a rotational fan scan using a side-firing transrectal probe.  1.4.1  Tissue M o d e l i n g  The trend toward computer-based medical simulators brought a high demand for modeling the complex characteristics of the tissue. Neglecting the inhomogeneities, anisotropy, relaxation, and hysteresis, soft tissue is generally considered to be a non-linearly deformable elastic material. In computer graphics, deformation modeling has been studied for noisy data fitting, clothing, facial expressions, and human/animal characters. [24] gives a survey of deformable modeling in computer graphics. In tissue modeling, one of the proposed heuristic approaches is the ChainMail algorithm introduced by Gibson [25].  It estimates the object as small patches which only transfer their  motion to their neighbour when their displacement is higher than a preset value. ChainMail does not supply a model for haptic force feedback calculations and the deformations are not realistic. A physically-based approach, the mass-spring model, approximates the body as a set of masses connected via springs (and possibly dampers). Waters [26] was the first to model the muscles of the face and generate gestures using this model.  Later, Terzopoulos and Waters  constructed a 3-layer dynamic mass-spring mesh for face in [14].  Although some degree of  realistic visual and haptic feedback can be achieved as shown in cloths made of double cross springs by Promayon [27] and laparoscopic simulator by Downes [15], the performance of massspring models highly depend on their spring network structure. Even after proper adjustment  1.4 Related Work  9  of the parameters, the models can only accommodate certain interactions and the instability is a major problem. Unlike the above mentioned approaches, the finite element method ( F E M ) and boundary element method (BEM) are based on continuum mechanics.  James and Pai introduced B E M  in [28] by solving the boundary integration of Navier's equation formulating the linear elastostatic objects. However, this method can only accommodate one material in a model and only solves for the boundary deformations/forces. In their later work [29], the same authors achieved realtime deformable textures using graphics hardware. They stored precomputed modal vibration models in graphics card memory and in real-time excited them according to rigid body motions. Continuum mechanics defines the deformation of an elastic body such that the equilibrium of the internal strain energy and the work done by external forces is preserved.  Although a  spatially continuous model is possible, because of the computational limitations in solving them on computers the implemented forms are ultimately discrete. Thus, the object is subdivided into smaller pieces and numerical interpolation/integration techniques are employed to approximate the system at discrete nodes.  This is called the finite element method.  performance dependent on this discretization (mesh).  It makes the solution  Meshing techniques are summarized in  Section 1.4.3 and our approach is given in Chapter 3. The earlier work of Terzopoulos et al. [30] showed the potential of F E M and finite difference method, which a special form of F E M . Cotin et al. [31] demonstrated a real-time visual performance in an FEM-based surgery simulator using linearly elastic materials with quasi-static deformations. Zhuang and Canny [32] reached real-time rates for large meshes using quadratic strain and dynamic F E M , which incorporates mass and damping effects onto the basic stiffness equation. They used a lumped mass matrix and exploited the small stiffness of the soft material to use explicit Newmark integration scheme for large time steps ([33]). Picinbono et al. [10] introduced an incompressibility constraint on a nonlinear dynamic mass-tensor model with lumped mass and damping and they simulated anisotropy. Both above mentioned models by Zhuang and Picinbono et al.  use St.  Venant-Kirchoff  elasticity which is the non-linear (non-simplified) version of the linear elasticity described in Appendix A . Debunne et al. [34] implemented a space-time adaptive mesh refinement on top of this same non-linear dynamic simulation scheme. Wu et al. [35] and Nienhuys [36] attempted other elasticity formulations (e.g., neo-Hookean). [35] implemented an adaptive mesh refinement using offline progressively coarsened meshes. [36] also compared some iterative solution techniques such as conjugate gradients (CG) and applied some 2D mesh modification and 3D mesh refinement techniques on incision and needle insertion examples. Conceptual hypotheses can be easier experimented in predefined media (phantoms)  than  complex tissue structures. The concepts laid out and validated through these phantom studies are then considered to be applicable to various tissues using their biomechanical properties reported in literature. In a recent study, Kerdok et al. [37] validated their finite element deformation  1.4 Related  Work  10  using the volumetric displacement data of markers embedded in a rubber cube and tracked by 3  C T (Computed Tomography).  1.4.2  Needle Insertion  The simulation problems related to needle insertion slightly differs from the other tool interactions in a couple of ways: it does not manipulate only the organ surface; the friction is an effective force during its interaction; and it is flexible unlike mostly-rigid surgical tools. The pioneering work of DiMaio and Salcudean [38], solved the above mentioned issues in 2D using a condensation technique, stick-slip friction with low-rank boundary condition switch, and local frame low-rank matrix update, respectively. The related work on those issues are presented below. In conforming F E M the needle surface forms a boundary condition for the tissue. Although Berkley et al. [39] showed a virtual coupling is possible by interpolating the forces/displacements, this interpolation brings an overhead and also the boundary criteria to be applied become vague. Moreover, the increase in the number of manipulated nodes negatively affect the speed of condensation techniques mentioned in Section 2.1.2. Therefore, a conforming tissue model with nodes falling on the needle shaft was chosen. To retain the tissue mesh nodes on the needle three approaches were implemented: Alterovitz et al. [9] and Nienhuys [36] insert new nodes while subdividing elements; Alterovitz et al. [9] and DiMaio [40] relocate nodes onto the needle path; and the haptic implementation of [38] forces mesh nodes to be on the needle which locally deforms the tissue. The tissue nodes on the needle can efficiently be solved using the condensation technique introduced by Bro-Nielsen and Cotin [41].  In a recent work, Gosline et al. [42] introduce fluid pockets onto DiMaio's  simulation [40] by adding their boundary force constraints onto this condensed working matrix. Generally, the online mesh modification is a major problem of F E M also required during incision [36]. Nonetheless, more direct approaches such as simply removing elements still exist. For instance, Cotin et al. [43] implemented this on their hybrid elastic model with partly quasistatic and partly dynamic elements. Needle insertion forces were studied by many researchers. Some earlier methods employed look-up tables such as the 1-DOF (degree of freedom) epidural block simulator of Hiemenz et al. [17]. Simone and Okamura [44] used a modified Karnopp model for needle friction estimation during in vitro insertions into the liver. [45] attempted a cantilever-beam-like needle flexibility to model the bevel deflection from in vitro muscle insertions, however in their work this model did not perform well. It makes an assumption of a distributed lateral force along the shaft dependent only on the insertion depth and the needle diameter as opposed to DiMaio's proposal of one deflection force at the tip as a function of velocity in [40]. O'Leary et al. [46] also attempted an insertion force modeling of bending and friction, however their approach of using the combined force on the base did not 3  This data set is available on the internet.  1.4 Related Work  11  lead to significant results. Although Okamura et al. [47] modeled the in-vivo tip cutting force and the insertion-retraction hysteresis effect, this and many other studies only considered the needle base forces.  DiMaio  and Salcudean [16] decoupled the shaft friction from tip cutting force and validated velocity dependent models physically. They were the first to model a velocity dependent shaft friction and tip cutting force through F E A of phantoms during needle insertions. In their later work [48], they present a needle steering scheme using potential fields. If needle steering is divided in two, steering using tip bevel and redirecting the tip by base manipulation, they explored the latter one. For needle flexibility, they employed a large-strain linear elastic model, that uses the GreenLagrange strain tensor. Webster et al. [49] recently showed that the highly flexible needles follow a path of constantcurvature for a fixed bevel rotation. They used a 3-DOF nonholonomic bicycle model for the tip bevel and found its parameters experimentally via insertions into semi-transparent and stiff phantoms. Alterovitz et al. [50] incorporated a simplified version of this concept into their needle insertion simulator.  1.4.3  M e s h i n g Techniques  Tetrahedral meshes are commonly used in many areas including solid modeling [51], haptics [34], fracture [52], fluid flow [53], heat transfer [54], surgical simulations [55], and biological tissue modeling such as the brain [56], the knee [18], the liver [35], and the uterus [57]. Therefore, generation of meshes with desirable characteristics has always been an interesting research topic. However, in this broad field mesh requirements highly depend on the application. The optimum mesh for soft tissue modeling with large deformation, solid modeling with small deformation, and fluid flow and heat transfer modeling with no deformation differ significantly.  Likewise,  the quality measures for mesh/element evaluation also depend on the problem type [58]. Many solutions were proposed for specific areas, which would not be considered even a suboptimal solution in another field. For example, it is desirable to have anisotropically compressed elements in boundary layers of a mesh for an Eulerian fluid flow simulation [59], since the solution gradient is obviously smaller in the direction of the flow. However, those elongated elements would result in anisotropy in a deformation simulation. Therefore, we first defined the desired characteristics of a deformable soft tissue mesh for prostate brachytherapy in Chapter 3. Although there are various algorithms for tetrahedral mesh generation, most of them are based on and classified in one of the following three methods. The algorithmically simplest one is the octree technique, which was originally developed at the Scientific Computation Research Center of Rensselaer Polytechnic Institute [60]. It is the recursive division of the model space into smaller elements until the desired resolution is reached for the elements that are cut by the object boundary. Another meshing technique is the advancing fronts [61] which progressively builds the tetrahedra starting from the triangulated surface and maintaining an active front that moves inwards by the formation of new elements.  In 3D, for each triangular facet on  1.4 Related  Work  12  the moving front a fourth node position for a well-shaped tetrahedron is calculated and the corresponding element is added. Before this node addition, the current nodes are also considered for a well-shaped element formation. This technique initiated with a poor surface triangulation will engender ill-conditioned tetrahedra. Moreover, its performance degrades when the fronts are merging which may as well occur near the surface with high curvatures. The popular F E A software A N S Y S ® i n c l u d e s a version of Lo's advancing fronts implementation [62]. The third and the most exploited tetrahedral meshing technique is based on the criterion set by Delaunay [63] in 1934. Until the work [64] in 1981, this well-known criterion was not utilized for mesh generation. The Delaunay tessellation of a set of vertices V in 3D has the property that the circumsphere of any tetrahedron encloses no other vertex of V. Since this criterion is just a connectivity definition, not a meshing algorithm itself, extra tetrahedralization methods are required for mesh generation. These methods basically differ by their decision mechanisms on the position of the new point that is to be introduced. For example, a simple approach can be filling the space with a regular grid. More complicated approaches usually first tessellate the triangulated surface nodes and then considering the quality and relative positions of the current elements, they insert nodes to divide the ill-conditioned tetrahedra. Shewchuk [65] introduced a guaranteed quality triangulation based on the node insertion at the element circumcircle centers and he implemented this in a 2D mesher . 4  The method [67] developed at INRIA employs a  incremental node insertion from a list of pregenerated candidate vertices until a background sizing function is satisfied. This was also implemented in their meshing software. Mesh postprocessing schemes such as smoothing, cleanup, and refinement can be applied to improve element qualities. Smoothing refers to the slight relocation of the mesh nodes.  Some  of the smoothing techniques use averaging, optimization-based, or physically-based (e.g. spring interaction between neighbours) methods. optimized using edge/element swaps (flips).  Cleanup is the step where element connectivity is It is used either for shape improvement, which  leads to better-quality elements, or for topological improvement, which optimizes the number of edges sharing a single node (node valence or degree). Refinement is the reduction of the local element sizes by node insertion for the purpose of capturing a local physical phenomenon or simply increasing the mesh quality. One drawback of the postprocessing may be that most of the methods above violate the Delaunay criterion of a mesh. Successful adaptive approaches were proposed based on the combination of mesh postprocessing steps and an iterative application of the Delaunay tessellation iteratively. Nevertheless, Delaunay techniques do not perform well in 3D, since many straightforward 2D methods such as edge flips become excessively complicated when adapted to 3D and may be intractable due to the non-polynomial complexities of the algorithms involved. In 3D, iterative edge flips may even not terminate. Moreover, slivers (very flat elements with negligible volume) become another source of problem in 3D despite the sliver removal steps.  [65] pro-  vides an overview of these methods including a discussion on why theoretical results often fail 4  It is available on the internet at [66]  1.4 Related  Work  13  in practice. It is very common that objects to be meshed are not convex and they contain holes or subregions (e.g., different materials) separated by internal boundaries. as boundary mesh generation or boundary conformity. tetrahedralization  This problem is referred  One solution is a constrained  Delaunay  ( C D T ) , which can generate a mesh that has no facets cutting the internal  boundary. A general approach is the facet recovery by edge swaps. However, it is extremely complicated and there are even cases where no valid tessellation is guaranteed without new point insertion [68]. During the Delaunay criterion check of a C D T , circumspheres passing to the other side of a constraining boundary are ignored. Therefore, elements of C D T s are not necessarily Delaunay. Nonuniform Rational B-Splines (NURBS) is a powerful way of geometrical representation. In meshing, N U R B S can represent the boundaries of a region as a set of parametrical curves instead of line segments between nodes. Since the prostate boundary has low curvature and our segmented nodes are already finely spaced, we did not employ any N U R B S techniques. A recent study [57] achieves an almost uniform tetrahedralization as follows: first the boundary cutting elements on an initial uniform mesh are subdivided; then the tetrahedra located outside are removed; and finally the positions of the remaining vertices are iteratively optimized considering the element quality measures. However, Chapter 3 mentions the F E A drawbacks of uniform meshing. Detailed surveys of mesh generation can be found in [69] and online at [70]. In general, the initial surface mesh is preserved and the number of elements is comparatively large. The high demand by the recent soft tissue studies lead to a lot of research in the area. There are methods promising well-conditioned meshes for subsequent large deformations, e.g., in [71] the refinement of an initial mesh with body-centered-cubic (BCC) crystalline structure towards the surface. However, in those mainly the surface representation and interaction are focused.  Therefore,  these coarse inner structure is likely to fail during a haptic interaction with inner nodes. O n the other hand, there is also related prior work (e.g., Debunne et al. [34]) generating good-quality elements locally from coarse meshes using adaptive meshing techniques. The types of elements are named according to the interpolating basis functions used. For instance, in our simulation T E T 4 elements are used. This name refers to a tetrahedron where inner displacements are approximated linearly using its four corners. Higher order shape approximations (e.g., T E T 1 0 with 6 more nodes on the edges) obviously perform better with a computational complexity drawback.  1.4.4  U l t r a s o u n d Imaging and 3D Reconstruction  Ultrasound offers tomographic imaging capability in a rapid (10-100 images/sec), noninvasive manner [72]. Tremendous progress over the past few years has improved its image quality and visualization clarity significantly. During brachytherapy the only continuous imaging feedback is the T R U S . It gives information on the positions of the needle tip and the already implanted seeds  1.4 Related  Work  14  with respect to the volume. That is why this imaging modality is crucial in decision making and so it should be incorporated in our simulator. The problem to be tackled is the generation of T R U S images in a deformed volume, the initial undeformed state of which is known in terms of regularly sampled B-mode T R U S data. If we were to divide our problem into two parts, volume deformation and the synthesis of T R U S images out of this deformed volume, the second part becomes a slicing problem of 3D (reconstructed) ultrasound data. The first part is attempted using FEM-based techniques (basis functions) in Chapter 4. In order to gain more insight into the second part, let us present a summary of prior work on 3D ultrasound volume reconstruction. Traditional single 2D B-mode ultrasound mostly relies on the operator's interactive ability of acquiring a properly oriented image of the target structure for diagnosis. This highly depends on 3D comprehension of the anatomy and its 2D ultrasound representation. O n the other hand, 3D ultrasound can be sliced in any direction, can be rendered for easier visualization, or can be stored and reviewed at a later time. In this concept, reconstruction  refers to the generation  of 3D volume data on computer using the 2D slices acquired from the region of interest.  So,  basically this 3D data can then be displayed using computer graphics techniques or resliced in order to generate 2D B-mode images. Nevertheless, note that in our system this 3D data changes dynamically as the tissue deforms and this necessitates additional steps before reslicing if the same approach is followed. For volume reconstruction, scanning can be done in various ways. [72] classifies the major ones as: M e c h a n i c a l : The transducer array or the entire probe is moved mechanically. Regular sampling achieves an accurate geometry, however, the bulky mechanisms are generally not practical. P o s i t i o n sensing: The probe position is tracked using sensors (e.g., magnetic, infrared). The probe is easier to manipulate, but the data is harder to reconstruct due to sensing errors. T r a j e c t o r y a s s u m p t i o n : The operator moves the probe along a predetermined trajectory (e.g., spatially regular data sampling generated from slices sampled at regular temporal intervals during a constant velocity sweep). Distortions and inaccuracies are common in this method. 2 D t r a n s d u c e r arrays: A pyramid-like volume is scanned using dynamic beam steering.  A  relatively small region is scanned, but at a fast rate (suitable for echocardiography). The reconstruction from 2D images can be done using feature-based or voxel-based approaches.  The former segments the surfaces or the volumes in the acquired data, and only  displays the simple anatomy. Although it is fast, subtle tissue information is lost and the segmented boundary can easily be erroneous.  The latter puts each pixel at the corresponding  3D position. A l l information is retained, however, large data files must be manipulated and fast algorithms are needed. Some of the factors that play a role in 3D reconstruction quality are the accuracy of the initial data position information and the selection of a coordinate system and a bounding box.  1.4 Related Work  15  These problems, such as preprocessing of the data to minimize patient motion and breathing artifacts, have been tackled by several authors. However, since this thesis does not focus on the above mentioned issues, only the reconstruction phase (i.e., the resampling into a regular grid) is addressed here. Reconstruction can briefly be divided into two stages, a bin-filling and a hole-filling, as suggested by Rohling et al. [73]. The former is also referred as compounding. It is the process of combining multiple pixel information for a voxel.  This necessity arises either because of  the voxel/pixel ratios greater than 1, or more commonly because of the overlapping B-scans. Detmer et al. [74], Rohling et al. [75], and Leotta and Matin [76] achieved some of the significant contributions in 3D freehand ultrasound compounding. Specifically, [75] introduces an imagebased registration approach of the compounded data, to avoid common problems (e.g., blurring and misregistered structures) in spatial compounding. Bin-filling is not applicable to our problem, since in our proposed initial data does not have more than one value at a coordinate. However, since the T R U S probe is moving through the volume many times during the whole procedure, effective methods that utilize compounding techniques are also applicable. The hole-filling stage is the interpolation, where the values of voxels which have not yet been filled are to be determined.  The interpolation can be done in several ways.  a review and comparison of scattered data interpolation methods in [77].  Franke presents  However, they are  not very helpful in tracking features for our case of having a much greater distance between adjacent B-scans (5 mm) than their internal resolution (~0.17 mm). Another challenge is the complex characteristics of speckle on individual B-scans.  The basics of speckle formation is  given in the early work [78] of Thijssen and Oosterveld.  Regarding the high dimensionality  of the problem, simpler approaches are usually preferred.  In addition to the common ones  like voxel nearest neighbour (VNN) and pixel nearest neighbour (PNN), various others such as distance weighted (DW) [79], radial basis functions (RBF) [80], tetrahedral grids [81], and ellipsoid Gaussian convolution kernel [82] were also proposed.  [73] provides a survey of well-  known reconstruction methods. Visualization, one of the major challenges of 3D ultrasound, is approached by slice projection, surface fitting, segmentation, volume rendering and depth cue enhancement techniques. However, our problem is generating only a (2D) B-mode slice of the volume, therefore, visualization of 3D is not one of our concerns. On the prostate side, 3D ultrasound shows some evidence of better prostate disease evaluation and more accurate volume estimation. It has the potential to improve the planning and evaluation of brachytherapy procedures. After its introduction by Watanabe [83] over 30 years ago, transrectal ultrasound had a very important role in diagnosis and intervention of various prostate diseases. Tong et al. [84] produced a 3D prostate volume rotating a side-firing transducer, although the modern prostate brachytherapy is performed using end-firing transducers. One of the major difference of our problem from the above mentioned ones is the spacing between individual data slices. T R U S images have a 5 mm inter-slice spacing as opposed to an  1.5 Proposed Simulator  16  only « 0 . 1 7 mm pixel spacing on a slice. This prevents an efficient use of the typical interpolation approaches, since the features such as speckle change a lot in 5 mm and the speckle shape and distribution are intractable to model [78]. Although studies showed the potential of real-time visualization based on raw B-scans [85], the construction of a regular grid is still the preferred method. It enables further processing of the data with off-the-shelf algorithms. However our proposed method synthesizes the images directly from raw B-scans, since it takes much less time compared to the rest of the algorithm and the volume reconstruction method would bring huge memory requirements.  Nevertheless,  for a shorter computation time this method can also easily be applied on a volume which is reconstructed by interpolation. Simple interpolation algorithms are successfully used on reconstructed volume sets to achieve B-mode slices. Aiger and Cohen-Or [86] implemented this idea in a commercialized application called UltraSim™ ignoring the directional effects of ultrasound.  [8] included a similar image  generator as the ultrasound modality in their cryosurgery simulator. There are available commercialized 3D ultrasound systems doing such volume generation and slicing.  1.5  Proposed Simulator  A 3D needle-tissue interaction model is proposed for the specific application to prostate brachytherapy. In this interaction simulation, the coupling of a flexible needle with a linearly elastic tissue is achieved where the tissue is modeled as a soft deformable body by extending the quasistatic and linearly elastic F E A techniques presented in [40]. In this simulation implementation, the reaction forces on the base are computed while the needle base is being manipulated. A T R U S imaging method of the deformed tissue model is also proposed as a visual feedback modality.  1.6  Thesis Organization  Chapter 2 describes the F E A techniques used in simulation. The reader is referred to Appendix A and the cited books for a further insight on F E M . Then, Chapter 3 details the mesh generation methods attempted and the final mesh applied. In Chapter 4 the deformation in 3D ultrasound and its reslicing for intra-operative T R U S synthesis are discussed. The overall structure of the simulation with the data structures and implementation details are presented in the following chapter, Chapter 5. Finally, our achievements, conclusions, and possible future directions are laid out in Chapter 6.  Chapter  2  Needle-Tissue Interaction M o d e l Both the tissue and the needle are modeled using the finite element method.  The reader is  referred to Appendix A and related books [87, 88] for an introduction to finite element analysis ( F E A ) . Many of the following approaches were extended to 3D from the work of DiMaio [40]. In brief, the tissue and the needle are modeled using F E M , the displacement of the needle is applied to a set of tissue model nodes, and the force feedback of this model is applied back to the needle F E M model. As the needle advances, the set of tissue model nodes that the forces are applied are updated when needle tip penetrates a new tissue element.  Considering that  these nodes slide along the needle their positions are updated using a stick-slip  friction model.  This model basically either fixes the node onto the needle shaft (by setting three displacement constraints) or lets it slide along (by setting two displacement constraints in lateral axes and a force constraint along the needle shaft axis). A l l of the above are discussed in detail throughout this chapter. In here, our ultimate goal is to lay the grounds for a stable and realistic model of flexible needle and 3D soft tissue. To achieve this the following techniques are applied.  2.1  Stiffness Matrix Inversion and Condensation  For visualization of an F E M model, the displacements of some (e.g., the surface) or all of the model nodes are needed. In the proposed simulation, US is the visualization modality. So, the generation of the T R U S images requires the displacements of a subset of model nodes. Thus, u = K~ f  (2.1)  1  needs to be solved, where K~  l  is the inverse global stiffness  matrix.  The displacements are  needed at 12 Hz considering the T R U S screen update rate during brachytherapy procedure. Although visualization is not an issue in haptic interaction and it can mostly be handled by a separate computer, in the proposed implementation, (2.1) is handled inside the haptic loop. This is the scheme followed by DiMaio in [40] and it was proven to work efficiently in his 17  2.1 Stiffness Matrix Inversion and Condensation  18  implementation although for very large meshes speed issues may arise. Various numerical methods can invert this K matrix relatively fast. For example, L U decomposition is a powerful tool considering the very sparse characteristic of the stiffness matrix. Reordering the rows/columns using Cholesky Decomposition or reverse Cuthill-McKee method can also improve the speed. However, for large 3D meshes, offline precomputation of this inversion is necessary for real-time performance. Even then, an online inversion might be needed if the mesh structure is modified such as an element is removed (e.g., incision) or a node is relocated (e.g., mesh adaptation in Section 2.6).  2.1.1  Inversion  Considering only a subset of displacements (e.g., the anatomical boundaries, seed locations, or elements on the T R U S plane) are needed, multiplying only the rows of K~  l  corresponding to  the these displacements would be enough. This would decrease the computations from to 0(n).  0(n ) 2  However, in this thesis this is implemented for all the node displacements since the  required nodes for ultrasound synthesis are hard to predetermine and the positions of some element corners are also needed at haptic rates (e.g., the corners close to the needle tip are used for element penetration check described in Section 2.5). The above mentioned method is half of the condensation approach by Bro-Nielsen [89]. The other half suggests using only the columns of K~  1  which correspond to the applied external  forces on the system. The fact that external forces on the tissue mesh exist only at the tissue mesh nodes, which are constrained on the needle, led to the exploitation of that other half. These nodes are called the contact nodes or the working nodes. In other words, while the needle base is being manipulated, the interaction with the tissue model only occurs at the contact nodes, which is a subset of tissue mesh vertices lying on the needle. Thus, if (2.1) is reordered so that contact nodes are at the end, the formulation becomes:  fo  oo  fi  /  0  = 0  u  fi )n (2.2)  where for the vectors u and /  the subscripts 1 and 0 denote the rows corresponding to the  contact nodes and the rest of the nodes, respectively.  2.2 Contact  2.1.2  Behaviour  19  Condensed Subsystem  The solution to (2.2) requires the non-zero forces, /  l  5  on the contact nodes. During the simulation  these contact boundary conditions are set by the needle. However, they are mostly applied as pure displacement boundary conditions or as mixed displacement-force boundary conditions (these cases correspond to the stick and the slip states, respectively, as described in the following section). Fortunately, the entire relationship between the displacements u\ and the forces f  1  of  the contact nodes can be described by the set of equations on the last line of (2.2). This smaller system of equations is called the condensed subsystem. Its matrix representation (K )^, -1  is a submatrix of the inverted global stiffness matrix K matrix and denoted by Kjj  - 1  , is named the condensed  which (system)  basically to refer further in the implementation chapter. Note that  this definition attributes only to a special (the initial - before any update) form of the working matrix that will be defined later in this chapter. As an example, when all contact nodes are stuck, there will be displacement constraints on all of them and the condensed subsystem is then expressed as: f ^ K ^ m .  (2.3)  Here if the number of working nodes is w, the dimension of the system will become 3w. This set of equations, consequently its matrix definition, need to be updated at each iteration according to the contact node states (Section 2.3) and needle orientation (Section 2.4).  2.2  Contact Behaviour  After the successful implementation of the stick-slip approach by DiMaio [40], we integrated the same method into our simulator. It had been proven to be computationally fast and yet visually realistic. It is based on physics and makes use of the fact that two surfaces start slipping only if the interaction force exceeds their friction and they continue slipping as long as they are in motion with respect to each other. This shapes the quasi-static behaviour of our simulation. At an iteration, a node either slips along or sticks onto the needle shaft. This decision for each contact node is made at every iteration depending on its current state. Our stick-slip automaton is shown in Figure 2.1. Here F is the static friction threshold at the needle-tissue interface and s  Ft ex  is the magnitude of the axial force exerted by the contact node on the needle surface, and  Vi and Vi-i are the velocities in the current and previous iterations of a simulation, respectively. So, the sign of Vi •  shows whether the direction of the motion has changed since the previous  iteration. A direction switch implies a zero-crossing of the velocity which consequently denotes the end of slipping due to the above mentioned kinematic fact. After being constrained onto the needle, a contact node stays on it until it slips off the tip. Therefore, its lateral displacement with respect to the local orientation of the needle is set by the needle shaft. On the other hand, along the shaft axis the contact is switched between force and  2.2 Contact Behaviour  F  20  <F  ext —  x  s  Figure 2.1: Automaton describing the stick-slip friction model of a contact node.  displacement constraints such that in the stick state the displacement is set by the needle, and in the slip state an axial friction force is exerted by the shaft. These axial constraint switches are applied using the boundary condition change technique described in Section 2.3, where the local axis along the shaft is chosen to be y in our particular implementation shown in Figure 2.2.  Figure 2.2: Tissue and needle-base axes in our implementation. Since the unknowns in (2.3) switch between displacement and force, depending on the boundary condition, this system of equations can be rearranged in order to have all the unknown axis components on the left-hand side as: x = where Kyv  Ky  (2.4)  w  is the working matrix, which is the boundary condition changed and local frame  adjusted version of the condensed matrix KH at any simulation iteration. In order to correctly define a relationship between displacements and forces of the contact nodes, the matrix needs to be updated if there is a stick-slip state change at an iteration. For example, at any given time if the contact node at the tip is stuck while two other contact nodes are slipping along the shaft, x and y should look like:  *  =  [/f  u\  fl  /|  u  V  =  u\  f\  u\  uf  f  y  v 2  2  fi u\  /£  fl  u% u\  fl u\  T  (2.5)  2.3 Boundary  Condition  Change  21  The above scheme necessitates the constraint application on a local frame that is axially oriented with the needle shaft at the contact points. For the three contact nodes given above, an instance of local frames on a bent needle can be oriented as in Figure 2.3. When the needle base is manipulated and/or the needle shaft flexes, these local frames also need to be updated (rotated) accordingly as described in Section 2.4. z  Figure  2.3  2  2.3: A n instance of contact node local frames lying on a bent needle.  Boundary Condition Change  When a contact state changes, in (2.4) one row of the knowns y has to be swapped with the corresponding row of the unknowns x in order to accommodate the new boundary condition. The procedure for this operation is easier to show on a sample case. Nonetheless, note that the following operations are applied without loss of generality. So, this sample case can later be generalized. First let us assume that all the contact nodes are stuck. Thus, the condensed subsystem is denoted as:  = fl  fc(l,l)  y  w  (l,r)  k  Ul  Ui-l  fi-1 fi  K  —  fi+1  fr  (2.6)  Ui Ui+1  u  r  Here a generic variable r is used just because this same technique is utilized also for mesh adaptation with a full-rank matrix. So, the method and the related appendix can be applied to this other case using a proper r value. In here r is obviously 3io, where w is the number of working nodes. For any given node axis, whenever the force is known and the displacement is sought, the set  2.4 Local Frame  22  Rotation  of linear equations to be solved becomes: Ul  h  Ui-l  fi-1 Ui  = K  t  (2.7)  fi  w  Ui+1  fi+1  u  r  This updated working subsystem which is denoted by this new matrix K  w  consequently has  its unknowns in the vector on the left, which is the new x, and the known variables on the right, which is the new y vector. In Appendix B , K K w where Cj and  is proven to be of the following form:  w  K w  are the modified forms of the i  th  1  (2.1  ~CjT' i k(i,i) column and i  th  row of K , w  respectively, such  that: Ci  [ (l,i)> fc  i  fyi-l,i)>  [%,1)>  (k(i,i) + 1);  k  (k(i,i) ~~ 1);  •••) %,r)]  (r,*)J  (2.9)  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(w ) 2  computations for this operation. Throughout this  thesis, computations refer only to multiplications and divisions. 0(w ) 2  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  23  2.4 Local Frame Rotation  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 ° / j in the initial coordinate frame can be expressed as follows:  =  "°/f" OfV Ji Ofz . Ji . where \ii and f 1  i  =  Rw  R12  Rl3  R21  R22  R23  R31  R32  R33  . \  R11  R\2  Rl3  " Vf "  R21  R22  R23  R31  R32  R33  z  .  (2.10)  Ji .'ft  .  are the vectors after rotating the coordinate frame and the matrix i?3 3 is X  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:  <k= => x= 1  \D=  where K  w  K  w  \  [A K A) y T  (2.11)  l  w  K'  \  w  is the updated working matrix and:  A =  1  0 • •  0  0  0  0  1  0  0  0  0  0  0  Ru  R\2  Riz  0  0  0  R21  R22  R23  0  0  0  R31  R32  -R33  0  0  0 • •  0  0  0  • •  • •  0  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  R3I  Rz2  c<j)c9 c<p sd sip — s(j> cip  -R23 -R33  ccf) s9 cip + scj) sip (2.13)  S(f)c9 S(f> SO S1p + C(f)C1p S(f) sO C1p — C(f) sip —s9  c9 sip  c9 cip  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:  Ji OfV Ji  (2.14)  °/f  \  ° 3w X  c?3w  which can then be separated into: I  0  Ofx Ji 0  0  +  0  =  °fi  0  ^ Zw  0  0  Kw  °7?  +  (2.15)  0  0  \  0  Using (2.10) each term above can be denoted in terms of the rotated frame vectors x and l  Mx l  + N \  =K  W  [JV^ + A f V  (2.16)  2.4 Local Frame Rotation  25  where M and N are sparse transformation matrices of the form:  M =  N  =  1  0  0  •  ...  0  0  0  0  1  0  0  0  0  0  0  Ru  0  ^13  0  0  0  0  R22  0  0  0  0  Rz\  0  -R33  0  0  0  •  0  0  0  ...  1  0  0  •  0  0  0  ...  0  0  0  0  0  0  0  0  0  0  R12  0  0  0  0  R21  0  R23  0  0  0  0  R32  0  0  0  0  0  0  0  •  •  •  •  ...  (2.17)  (2.18)  0  (2.16) leads to the expression of x in terms of \j as: x  *x=  [M - K  ke =  K'  N}'  1  w  [K  M - N] y  (2.19)  l  w  \  w  which gives the updated working matrix K  w  .  The implementation of this formula requires the inversion of [M — K N}. W  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 Kf N3^x3™ wxSw  [03u)x3(j-i) -^3u)x3 ®3wx3(w-j)]  ^ ^ na  n  a  s  takes the form  3 non-zero columns, where j is the working node to be  2.4 Local  Frame  26  Rotation  updated. Thus, the matrix to be inverted is expressed as:  M  - K N  =  W  1  0 •  a  0  1 •  a  0  0 • •  0  0 • •  0  0 • '  0  0 •  (l,3j-2) (2,3;-2)  a  a  (l,3j-l)  °(l,3i)  • 0  a  (2,3j-l)  a  (2,3j)  • 0  0(3j-2,3j-2)  a  (3j-2,3j-l)  a  (3j-2,3j)  • • 0  (3j-l,3j-2)  a  (3j-l,3j-l)  a  (3j-l,3j)  ' • 0  a  (3j,3j-2)  a  (3u>,3j-2)  a  (3j,3j-l)  (3j,3j)  • 0  (3w,3j)  • 1  Q  (3w,3j-l)  a  a  (2.20)  For convenience, elements of the 3 x 3 central matrix will be renamed as:  On  0-12  «13  «21  «22  «23  ^31  «32  «33  a  (3j-2,3j-2)  «(3i-2,3j-l)  a  (3j-l,3j-2)  «(3i-l,3j-l)  «(3j,3j-2)  C(3i,3j-1)  G(3j-2,3j) (2.21)  (3j-l,3j)  a  (3j,3i)  a  Symbolic toolbox of Matlab gives the inverted matrix in the following form: 1  [M-K N]~ W  1  =  1  0 •  0  1 •  0  0 •  0  0 •  0  0 •  0  0 • • n'  a  (l,3j-2)  °(2,3j-2)  a  (l,3j-l)  a  '(2,3j-l)  '  (2,3j)  • ••  0  31  33  • 0  22 32  a  '(3w,3j-l)  '(3w,3j)  a  (3w,3j-2)  23  • 0  a  a  a  • 0  a  a  0  'l3  'l2  a  a  • 0  (1.3j)  a  a  '  (2.22)  • 1  where the elements are:  for the rows  *12  -C21  21  =  -Cl2  a  22  = ~ 22  a  31  =  -C13  a  32  =  [ a *( -)  -  3ij  '(3z+lJ)  =  (3i+2,j)  =  a  C  CnO( ) 3iji  c  _ c  13  =  -C 1  a  23  =  -C32  a  33  =  -C33  a  a  3J-2,3J-1,3J :  for the rest of the rows  -cn  a 11  23  + Ci a ,j) 2  (3l  21 (3i+lj) + C 2 2 0 ( j i j ) a  3  +  i  3  + Cl30(3i,j) + C 3a(3i+l,j) 2  e{l,...,3j'-3,3j+l,...,3u;} [  a  31 (3i+2,j)  C  a  Matlab® is a registered trademark of The MathWorks, Inc.  + C32a(3i+2,j) + C33a(3z+2,j)  (2.23)  2.5 Decision of Remeshing  27  where the c values are found in the following form:  while r =  -022033  Cll =  023032  C21 =  033012  -013032 r  C31 =  013022  -012023  013021032 -  r  r  013022031 -  033021  Cl3 =  =  C22  013031 - 0 1 1 , 0 3 3 = r  C23 =  C32  =  023 0 U - O 1 3 O 2 1 r  C33 =  023011032 -  r  022031  - 0 3 2 021  Olia  -031012  3 2  r  (2.24)  r  O12O21 - 0 2 2 0 1 1  r  033012021 + 022033011 + 031012023-  In every iteration, this update requires 0(w)  2.5  -023031  C12  computations per slipping node.  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]. boundary constraints difficult.  However, this would make the application of  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  (a) Advancing towards a new ele- (b) Tip penetrated into new element. ment.  28  (c) Remeshing the region,  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 | P P " | 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 (K ), e  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, e  2  = P1P3P4P, e3 = P1P4P2P, and e = P1P2P3P as seen in Figure 2.5. Removing those four 4  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 e ' = P P P P', e ' = P1P3P4P', e ' = P i P P P ' , and c\ = P1P2P3P' onto the global stiffness matrix K to have an intact mesh. Here, note that x  2  4  3  2  3  4  2  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 ' — P + u is written as follows:  K' = K -  global (K  ei  (P)) + ^  global (K  ei  (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 m a t r i x numbering ( A p p e n d i x A . 5 ) . N o t e that this modifies o n l y a s u b m a t r i x of K corresponding to the corner vertices of elements ej. F o r 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 e n d a n d call those nodes P j . L e t m be the number of such nodes (namely P ' s valence or order).  2.6.1  0(mn ) 2  Approach  T h e real-time code uses the inverse stiffness m a t r i x K^  as i n (2.1).  1  However, K  e  defines a  relationship from the displacements to the forces. Therefore, a direct o p e r a t i o n o n this K  -  1  is  not possible. O n e approach c a n apply the operations of (2.25) o n the o r i g i n a l g l o b a l stiffness m a t r i x K a n d t h e n invert the entire full-rank m a t r i x . A s discussed earlier, a real-time inversion of rank 3 n is infeasible. A fast approach is the extended version of S h e r m a n - M o r r i s o n formula, namely the W o o d b u r y formula [90]. If A K = - X) .  K  &  e  I  ( ) + 2J P  K e i  &  I  ( ')> t h e n the inversion of P  the p e r t u r b e d m a t r i x is given by: (K + VAKV ) T  1  = K  L  -  + V K' VAK)  \K- VAK(1 1  T  1  (2.26)  V K~  1  T  1  Here V is a 3 n x 3 m m a t r i x having the last row as a n identity m a t r i x of the form Is x3m-  This  m  requires 27nm(n + 2m) computations, i n a d d i t i o n to a 3 m x 3 m m a t r i x inversion w h i c h c a n be done i n 0(m ) 3  i n the worst case. A l t h o u g h n > m is a v a l i d assumption, m w i l l also be i n c l u d e d  i n the b i g O n o t a t i o n to express the effect o f node valence. Instead of this 0(mn ) 2  D i M a i o [40] employed another 0(mn ) 2  converts a part of K~  L  approach,  m e t h o d w i t h slightly more c o m p u t a t i o n s ( 1 6 2 m n ) . H e 2  to displacement constraints, a n d after a p p l y i n g the K ' e  operations o n  this region, converts it back. T h e conversion is done b y swapping the 3 m rows (corresponding to Pi) at the e n d using b o u n d a r y switches i n Section 2.3. T h a t yields to:  Up  1  L\2  Pn-n  U  F  P n-  r  l21  l22  Pn-  (2.27)  Pn-  U  Here u a n d / values are triplets containing a l l three axes.  N a m i n g the vectors conveniently  gives: L  N  l21  L\2 l22  (2.28)  u  7  Note that none of the submatrices L holds values of either K or K  1  directly. Instead they  31  2.6 Mesh Modification define a mixed boundary condition similar to (2.4) such that: U\  =  LuFi  F  =  L21F1 + L22U2  2  +  The advantage of this conversion is that now K  ei  liyiUl (2.29)  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  the inverse  (K~ )', 1  global stiffness matrix of the new mesh. Note that finding and adding/subtracting the individual element stiffness matrices, K , take negligible time compared to the other operations. ei  2.6.2  0(m n) Approach 2  We propose an approximation to this matrix update in 0(m n). 2  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 P , the displacement of P 2 is also affected from a force 2  on P i through the element PPiP  2  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.~ matrix for a regular mesh with 240 nodes and x  810 tetrahedra. This example was remeshed at a node and the difference matrix  (K" )' 1  —  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 « 5 % in terms of power of difference. We implemented this with a modified version of the boundary condition change function which only calculates Li ,L i, 2  2  and L 2 of (2.29). In there, Ln is the effect of other nodes' forces 2  2.7 Needle Model (Fi)  32  on other node's displacements (U\).  Neglecting Ln  reduces the complexity of each switch to 0(mn). to 0(m n). 2  in the boundary condition change  Consequently, this eases the mesh adaptation  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 displacements and forces are expressed as follows: K(u)u =>  f =  K, + K,(u) u K<u + f J u )  (2.30)  where the subscript £ and q denotes the linear and quadratic terms, respectively. The hat (") is used to denote the needle vectors and matrices. For force boundary conditions, the equation above can be solved iteratively as shown by DiMaio [40] as follows: (2-31) Thus, if the forces on the needle model is known, its quadratic model displacements can be approximated iteratively. Next section shows how those needle mesh forces are found. However, since both needle and tissue are modeled using quasi-static elastic material, their coupling obviously needs to be damped for a stable behaviour. The maximum needle base manipulation speed and any internal step input such as the force during the node repositioning are the two dominant factors determining the minimum damping value to be used. The needle damping values are found empirically and applied on the needle forces of the next section before putting in (2.31).  2.8  Projecting Forces onto the Needle  The working nodes are the only points at which needle constraints are applied onto the tissue. Moreover, the F E A of the tissue solves for only those nodes' forces which are then applied back onto the needle to achieve needle-tissue coupling. However, similarly the needle model can accept constraints only at its mesh nodes, which are represented by the needle knots in the segmented  2.9 Integrating Needle Forces  33  approximation. Thus, while the displacements of the contact nodes are derived from the needle knots, the reaction forces of the tissue mesh are also to be projected onto these knots. Figure 2.6 shows an example with two working nodes  and P^f on the needle shaft. They  can be either stuck or sliding on the needle segments P^P^  and P^ Pf.  Assume that the forces  1  caused by tissue deformation on those working nodes are solved by F E M and found to be F\ and F . Then, the force on each working node is linearly distributed onto two closest knots as: 2  Figure 2.6: Projecting the working node forces onto the needle knots.  2.9  Integrating Needle Forces  For a haptic implementation of the simulation, the force feedback at the needle base has to be calculated in order to be applied on the user's hand. This feedback is equal to the sum of the forces along the needle, in addition to their wrench effects as seen in Figure 2.7. Since it is very unlikely to have an axial rotation torque around the shaft, r  y  is neglected.  The forces can easily be added up for all s needle segments to find the total base force as follows: s  = Y. ? F  F*  i=l  s  FY  = J2 ? F  i=l  FZ  s  = T. i F  (- ) 2  i=l  Calculation of the torque needs the distance from the base, which can easily be found in terms of the knot position in the local needle frame. Thus, the total torque becomes:  (2.34) i=l  i=l  33  34  2.10 Discussion  Figure 2.7: Integration of knot forces and torques on the base.  2.10  Discussion  The bottleneck of the presented approach is the time complexity of the mesh adaptation algorithm. Its dependency on node valence is also a concern, which constrains the mesh that can be used. Section 5.7 discusses implementation solutions regarding the speed of remeshing. Those proposed approaches also address the issue of current visualization calculations inside the haptic loop. Instead of mesh adaptation, node addition techniques can also be considered. For instance, in Figure 2.4, P' could be added as a new node to the mesh. This would require the removal of two elements (PP1P2, PP2P3)  and addition of four new ones (PPiP', P i P P ' , P 2 P 3 P ' , 2  P3PP').  Although in 3D six elements are needed instead of four, this method at least ensures a constant running time regardless of node valence as opposed to the current adaptation technique. Using this proposed addition with a Woodbury formula update could also make use of a hard-coded algebraic matrix inversion of rank 18 (6 involved nodes in 3D). Finally, the issue related to the undeterministic nominal position of the desired new node (P ) can be dealt with force interpolation techniques over iteration cycles so that the user will 1  not feel the discontinuity. DiMaio [40] reported these as force spikes without discussing their causes. Nevertheless, they have not resulted in poor haptic performance even without employing any remeshing (just snapping the closest node onto the needle). Although these spikes rendered some problems with our 3D flexible needle, they were overcome by damping (Chapter 5). The iterative solution of the needle in (2.31) has some stability issues for certain ranges of Young's modulus and Poisson's ratio. This may need to be addressed using different solution techniques or some approximations that can be validated using experimental needle deflections. Some preliminary work which is not presented in here has been done on this issue. This work has shown considering other needle bending models could be a promising (fast and accurate enough) approximation. This has to be extensively addressed in a future work.  Chapter 3  M e s h Generation The performance of F E M depends on the element discretization geometry which is called the mesh. In [58] Fleischmann et al. state that "for the finite element method, geometrical quality criteria suffice for most applications". There are many techniques and applications for creating those meshes.  However, for a (fast) haptic implementation also a small number of nodes is  desirable. This chapter shows how traditional Delaunay tessellation may perform poorly in those aspects due to our irregularly-spaced relatively-dense prostate surface definition. As a result, an advancing fronts implementation which also enabled us to set element sizes for different regions is described later in the chapter. The above mentioned geometrical quality criteria can be defined in various ways. Usually aspect ratio or dihedral angle are considered.  They are discussed in Section 3.3 on our final  generated mesh. Throughout this chapter the tackled problem is converting the contour polygons to a surface definition and generating a quality mesh with small number of nodes which are concentrated around our region of interest (being the common needle insertion paths) with a conforming boundary on this surface. After presenting our results local element renumbering which is needed as a final tissue meshing step for our implementation is described briefly. Finally our needle mesh is introduced.  3.1  Prostate Surface Construction  The contour data of the prostate slices are stored in the database of the B C C A . For an insight to the characteristics of the available contour data the treatment planning is first described briefly. This data was to be obtained for further processing. Then using these points which are lying on parallel slices a prostate surface definition was to be generated. Knowing this prostate boundary the entire region could be meshed accordingly.  35  3.1 Prostate Surface Construction 3.1.1  36  Treatment P l a n n i n g  The brachytherapy operation planning starts with a T R U S volume study of the prostate weeks prior to the procedure.  This ultrasound examination is recorded on V C R . 10 to 12 parallel  prostate slices are captured from this V C R and saved as 640x480 images in the B C C A patient database. During the following couple weeks, prostate segmentation, procedure planing, and its approval phases are carried out on a radiation therapy planning software called VariSeed™. First a doctor segments the prostate on each captured ultrasound image by identifying its contour. During this segmentation the pixels of the image are clicked using the mouse so that these clicked pixels will become the corners of a closed contour polygon. Since the prostate boundary is not clear especially around its apex, the physician also relies on the anatomical knowledge. In the second phase of planning, different radiation therapists prepare implantation plan proposals for the same patient.  Then the physician that will perform the procedure chooses one of them  considering dosimetry and the cancer level and location.  3.1.2  C o n t o u r E x t r a c t i o n via D I C O M - R T  The easiest way of extracting data from the database was exporting in D I C O M - R T (Digital Imaging and Communications in Medicine - Radiation Therapy) format in VariSeed.  This  format stores the images with their relative spatial information. Two additional files contain contour and dosimetry data. The contour file was read by the dicomread function of Matlab and a script was written to extract the contour corner points of the requested anatomy into an array. The anatomical structures segmented in the T R U S volume studies (of 20 study cases mentioned in next chapter) are the urethra, the prostate, and the rectum. Dicom contour file stores contour polygons (sets), volume image slices, and the anatomical structures in separate object structures since not all the images contain a contour polygon and not all the anatomical structures of interest span exactly the same slices. This information is parsed and rearranged using the Dicom manual from N E M A (National Electrical Manufacturers Association).  3.1.3  C u l l i n g the C o n t o u r Nodes  The resolution of ultrasound images is determined as 55.3 pixels/cm using the grid in the image, where the spacing is known to be 0.5 cm. When the contours are drawn over these images, it is seen that there are very close nodes (<0.5 mm) probably caused by the physician clicking the mouse button twice. There is also a node abundance such that many nodes are located almost linearly. This superfluity of nodes remains in the final mesh, since the meshing algorithms often add nodes on top of the current ones. Therefore, it was decided that some of the nodes could be eliminated initially without the loss of much anatomical information. Considering the prostate is a smooth structure having ellipsoid-like cross-sections, very closely spaced nodes are not carrying much spatial information. A Matlab code was written to iteratively remove the contour nodes, which are closer to the rest of the nodes than a set threshold. When  3.1 Prostate Surface Construction  37  two nodes are close, the one which would leave less gap after its removal is chosen to be removed. For example, if the minimum gap is between the nodes Pi and Pi+\ at an iteration of this algorithm, then the distances | P j _ i P j | and | P j i P j 2 | are compared and depending on which is +  +  smaller, either Pj or P j i is removed from the node list of that slice. So, for the given example in +  Figure 3.1 since dj_i < dj+i, Pi will be removed (otherwise Pj+i would be removed) leaving the list of remaining nodes as {...,  Pj_2, P t - i , P i + i , Pi+2, • • •}• This is applied iteratively until there  is no more nodes spaced less than the set threshold on a slice. The same procedure is repeated for each slice separately.  Figure 3.1: Node culling example. A sample case having 262 nodes is downsized to 194 nodes using a threshold of 3 mm, to 160 with 4 mm, and to 142 nodes with 5 mm. Figure 3.2(a) shows a slice where the circles indicate the removed nodes and the stars are the culled ones after the application of our algorithm with 4 mm threshold. For the same case the nodes of all slices are shown in Figure 3.2(b). that most of the closely spaced nodes are around the base and the apex.  Note  This is due to the  physician inputting almost as many contour nodes for those regions as the other slices, although the contours have shorter perimeter here. Also these figures show the nodes which are probably erroneously added due to the physician double-clicking or clicking twice (instead of single-click). As seen those are also removed by this algorithm (an example is at the coordinate [50,30] in Figure 3.2(a)). A 5 mm culling threshold was chosen. One reason was the inter-planar distance being 5 mm which sets our minimum sampling in one dimension. The other was the qualitative observation of the final mesh. Although the number of nodes increased significantly with decreasing threshold, the surface and the mesh did not change much visually.  3.1.4  O b t a i n i n g the Surface Definition  The culled nodes are points lying on parallel slices. A first attempt was direct meshing from those surface points. However, simple Delaunay tessellation of all these nodes in Matlab gave us tetrahedra outside this region since the prostate surface was apparently not convex. Although there are algorithms to remove the elements outside an object, for any algorithm to work on this data first the surface definition was needed in order to describe the connectivity of those nodes  3.1 Prostate Surface Construction  38  °"*»-»  »  50  55  60  65  70  o .  ;s  75  X[mm]  (a) node culling on slice "Base+1.0"  (b) A l l surface nodes on 11 contours  Figure 3.2: Node culling with a threshold of 4 mm on a sample data. Circles are the removed nodes and the stars are the remaining ones.  on the surface. Thus, simple tessellation did not work and the problem is divided in two parts, surface generation and volume meshing. This section presents the solution to the first part. We attempted to surface construction by tessellating the nodes of each two consecutive planes separate from the others and then we were to accumulate the triangles that would be on the prostate surface. Although identifying those triangles was theoretically possible, it did not have a straight-forward solution either. Thus, the solution to this was searched in the literature. Keppel [91] has one of the earliest publications on surface reconstruction from planar slices. For each two neighbour nodes on a slice, he looked for a third one in the adjacent slice and formed the surface triangles this way. Various researchers have proposed methods dealing with contour tiling and branching problems, where the former refers to finding better-conditioned surface triangles and the latter addresses the differences in number of contour polygons on adjacent slices. Boissonat [92] proposed a method based on the Delaunay tessellation of each pair of slices. The method then extracts the boundaries of each slab, once their volumes are tetrahedralized. This (5)  is very similar to what we attempted above and its code (Nuages  [93] by INRIA) is available  on the internet. Nuages gets a set of simple closed polygons on parallel cross-sectional planes and outputs the set of surface triangles or volume tetrahedra representing the 3D polyhedra. The tetrahedral elements generated by simple Delaunay methods are generally poor-conditioned for F E M . Moreover, this software cannot generate C D T s and earlier trials of meshing the two regions separately faced the problem of non-conforming surface triangulation of the tissue tessellation inside and outside the prostate. Therefore, only the surface output is obtained using this software and this triangulated surface is meshed in a further step. The culled nodes are written in a file format that Nuages can read and a dxf file is requested as the output format for the triangular surface representation of the object, since this native  3.2 Meshing the Region  39  A u t o C a d ® f i l e format can be imported by many popular meshing software. For the sample case mentioned above with 142 culled nodes, a surface definition is generated consisting of 280 triangles as shown in Figure 3.3.  Figure 3.3: Nuages surface output (280 triangles) for 140 input nodes on 11 parallel slices.  3.2  Meshing the Region  First, element sizing was to be determined. Many meshing studies employ larger elements inside the objects based on Nicolson's empirical observation from his simulations in [94] that the spatial cut-off frequency of an elastically deformable body in response to external loads decreases faster than 1/d in terms of the depth d from the surface of the body. This result suggests that element sizes that are proportional to the depth can be accurate enough for many surface interaction simulations such as in [33]. However, since our model is manipulated via inner nodes, the same statement does not necessarily apply to our case. That is why a low-resolution mesh with similar element sizes regardless of the depth was aimed during the mesh generation for our simulation. Various software packages were considered for meshing the surface definition. Octree-based ones were neglected at the first place since they would have small elements on the prostate boundary. Many programs which cannot generate a boundary constrained mesh are also neglected since at least the prostate needed a conforming boundary separating it from the surrounding tissue. A good and easy to use option seemed the TetGen [95] package since it offers C D T s with high-quality tetrahedra for F E M . The C + + source code compiled on Linux accepts a polyhedral definition of the object and can output the mesh in various popular file formats.  The mesh  can be constrained on the minimum radius-edge ratio and the maximum element volume, and a post-processing option can attempt to remove slivers with high dihedral angles. The surface definition (dxf file) is translated into TetGen input format by a Matlab script. When its output is visualized, it is seen that with any quality setting many tiny tetrahedra on the prostate surface are generated. Moreover, the non-optimal surface triangulation generated by  3.3 Results  40  Nuages is almost not touched. It is concluded that this package is intended for small deformation solid modeling, since its underlying technique aims for feature preservation and gradual element sizing outwards to the surface. Thus, it ends up with ?s3000 nodes just for the interior of the prostate. In addition, the surrounding tissue requires mesh nodes as well. This made us aware that not many packages are intended to mesh for low number of elements. Owing to its extensive F E M analysis options, A N S Y S ® w a s as a strong candidate which was tried next. However, its complicated quality settings, relatively large number of generated nodes, and individual installation and stability problems (experienced with version 7.0) led us to consider other solutions. At this point, a package that was tried was found to satisfy the shortcomings of the above mentioned ones. G i D from C I M N E (International Center for Numerical Methods in Engineer1  ing) in Spain is a geometrical modeling, meshing, and visualization software. It uses advancing fronts so that a preferred element size setting can define the rate that the front advances (which in turn determines an average element height). It can also generate conforming boundaries allowing separate preferred element size settings for different material subregions. The prostate was imported from the dxf file. It was enclosed in a rectangular prism and different materials were attributed to the interior and exterior of the prostate. The width and the height of this rectangular prism were set to the dimensions of the ultrasound.  A larger  depth than the span of the volume study is meshed in order to account for the offset of the prostate underneath the skin. Considering the 5 mm culling threshold and slice separation, after experimenting with various values while observing the size of the resultant mesh, the tetrahedron edge lengths are chosen 7 mm for the prostate and 14 mm for the surrounding soft tissue. The two ends of the prostate, the base and the apex, have relatively fewer surface nodes since the large polygon segmentations on these surfaces directly passes from Nuages and only one or two nodes are added by G i D . During the further steps, when this mesh is used for the needle insertion simulation, it is concluded that more nodes on these surfaces would create a better path for the needle.  Since it especially passes through the apex quite frequently (for almost  every insertion), having a finer mesh on it would improve the simulation. Therefore, a further step is added before volume construction, where the first and the last contours are displayed and the user may add extra surface nodes on them. For the specific example presented, 5 and 4 nodes are added on the base and the apex, respectively. This increased the total number of nodes from 152 to 161.  3.3  Results  The final generated mesh has 570 nodes and 2801 tetrahedra. Figure 3.4(a) shows the bounding box surface and the prostate surface triangulations. A rendered prostate image with the 4 added vertices at the apex can be seen in Figure 3.4(b). In this mesh, the smallest element volume 1  Its free academic version [96] allows up to 700 nodes with 3000 tetrahedra.  3.3 Results  41  was 6.67 m m , where the largest one is 719.41 m m . T h e shortest a n d the longest edges were 3  3  3.06 m m a n d 24.04 m m , respectively. C u m u l a t i v e histograms of these quality measures for a l l the mesh elements are presented i n F i g u r e 3.5. I n those graphs the height difference of a bar from its left adjoining bar indicates how many elements that value applies. T h u s , the higher the gradient is, the more elements are concentrated i n that range of values.  For instance, F i g u r e 3.5(a)  has a steep ascent from « 3 0 0 to ssl600 between 4.9 a n d 6.5 a n d almost flat above 14, w h i c h consequently mean there are w l 3 0 0 elemens h a v i n g a m i n i m u m edge about 5-6 m m a n d only a couple elements w i t h a n edge longer t h a n 14 m m .  (a) Region boundaries  (b) Rendered prostate  Figure 3.4: M e s h visualization. Dihedral angle is the angle between two planes. For a n element, m i n i m u m and m a x i m u m d i h e d r a l angles give a general idea of its flatness. A too s m a l l m i n i m u m or a too large m a x i m u m w o u l d indicate the presence of slivers. I n our mesh, the smallest m i n i m u m is 9.94° where the largest m a x i m u m is 159.49°.  I n T e t G e n meshing package, a sliver is defined for a m a x i m u m  dihedral angle over 175°. O u r m a x i m u m value is well below this quality l i m i t . Moreover, the low rate of change for values higher t h a n 135° i n F i g u r e 3.5(b) a n d less t h a n 19° i n F i g u r e 3.5(d) show that there are very few low-quality elements. Well-conditioned elements are very i m p o r t a n t i n F E A , since the formulas a n d calculations supply better a p p r o x i m a t i o n s for elements w i t h geometries closer to an o p t i m a l shape. tetrahedra, this shape is a  For  regular tetrahedron w h i c h has four equilateral triangular faces. T h e  performance of F E M degrades as tetrahedra deviate from this o p t i m a l structure. various c r i t e r i a to measure this d e v i a t i o n .  T h e r e are  Aspect ratio is the c o m m o n t e r m used for this as-  sessment factor, however it is definition changes at different resources. Generally it is a ratio of the circumsphere radius or the longest edge length to the inscribed sphere radius or the shortest edge length. W e use the longest edge length d i v i d e d by the diameter of the inscribed sphere. Table 3.1(a) presents the number of elements w i t h certain aspect ratios. O u r mesh h a v i n g no  3.4 Local Element Numbering  42  # Elements  # Elements  (c) Tetrahedron volume  (d) Minimum dihedral angle  Figure 3.5: Mesh quality analysis.  elements with an aspect value larger than 15 is considerably good. Another element criterion is the radius-edge ratio which is the circumsphere radius divided by the shortest edge length. It is a weaker quality measure since it cannot identify many slivers. The values below 2 are preferable and our mesh has only 20 elements out of this limit as seen in Table 3.1(b).  3.4  Local Element Numbering  Some parts of the code requires a special local numbering of the element corners in order to make faster run-time operations. This corner ordering would make the decision of whether a point is inside an element or not. This basically needs to be checked for tip penetration in Section 2.5 and during ultrasound synthesis in Section 4.4. It also saves us from taking absolute values of the determinant, where the volume of an element is calculated using: Xi  2/1  Zl  1  X2  V2  Z2  1  X3  V3  Z3  1  X4  V4  ZA  1  (3.1)  3.4 Local Element Numbering  43  (a) Aspect ratio A s p e c t ratio 1.1547 - 2 2 - 2.5  (b) Radius-edge ratio  elements  Radius-edge  elements  0  0 - 0.707  253  31  0.707 - 1  1566  2.5 - 3  331  1 - 1.1  298  3-4  1284  1.1 - 1.2  222  4-6  989  1.2 - 1.4  240  6-10  156  1.4 - 1.6  113  10 - 15  10  1.6 - 1.8  67  15 - oo  0  1.8 - 2  22  2 - 2.5  Table 3.1:  14  2.5 - 3  1  3-10  5  10 - oo  0  quality.  If the corners are in right-hand ordering, (3.1) gives a positive volume. Right-handed local node numbering is geometrically equal to having the node P | on the same side with the surface normal, P{P\  x P{P\,  of the plane passing through P£, P | , and P | as seen in Figure 3.6. So,  the local numbering of element corners generated by G i D is processed by a Matlab script to ensure right-handedness before simulation starts. For each element, first its volume is calculated using (3.1). If it is positive, the order is not touched.  Otherwise, the nodes P | and P | are  switched, which puts the element in right-handed form, since it changes the side of P respect to P P P e  1  e  2  e 3  plane.  P4  6  r  Figure 3.6: A right-handed tetrahedron.  e 4  with  3.5 Mesh for the Needle Model  3.5  44  Mesh for the Needle Model  T h e needles used i n prostate brachytherapy procedures are long a n d very flexible. For instance, F i g u r e 3.7 shows the last 16 c m of a needle straight and bent.  For a 0.6 N orthogonal force  the lateral deflection is more t h a n 6 c m while the deformation is s t i l l elastic. T h e i r flexibility complicates the insertion m a k i n g the prediction of the t i p m o t i o n harder.  Nonetheless,  the  flexibility is a desirable property that also makes the needles steerable to reach occluded targets. T h i s steering challenge is incorporated into our simulator by modeling the needle  flexibility.  Figure 3.7: B r a c h y t h e r a p y needle deflection under a lateral t i p force. Different models i n c l u d i n g b e a m elements, an a n a l y t i c a l cantilever b e a m model, and spring models have been considered. Because of its p h y s i c a l validity a n d previous successful implement a t i o n as a needle m o d e l [40], finally F E M was chosen. T h e needle mesh is prepared m a n u a l l y using the same type of ( T E T 4 ) elements as the tissue m o d e l . T h i s enables us to use the same basis functions, solution techniques, and p r o g r a m m i n g procedures for the needle and the tissue.  x 10 J  Figure 3.8: Needle mesh topology. T h e criteria we chose as the guidelines for the needle mesh construction were filling the needle volume using the least possible number of flexibility-representative elements, i n order to decrease c o m p u t a t i o n a l time, a n d keeping at least one surface between adjacent elements, i n order to r o t a t i o n a l l y fix them. A simple uniform structure w h i c h comes to m i n d at first is shown i n F i g u r e 3.8. However, experimenting w i t h this i m p l e m e n t a t i o n revealed how susceptible F E A is to the a x i a l asymmetry of these elements and their connectivity. E s p e c i a l l y w i t h few elements,  3.5 Mesh for the Needle Model  45  the ratio of t i p deflection caused by equal-magnitude forces i n perpendicular axes raises up to 2:3. T h a t means a force a p p l i e d i n pure + Z direction causes a n asymmetric t i p deflection of 2 units i n - X , i n a d d i t i o n to the expected 3 units i n + Z . T o avoid this a new structure was designed i n w h i c h the elements are connected i n a way to cancel out each other's asymmetric effect. F i g u r e 3.10 shows the t i p deflection due to a constant lateral force rotating around the t i p for b o t h models. F i g u r e 3.9 presents sections of b o t h mesh designs. T h r o u g h o u t this thesis the t e r m section denotes a b u i l d i n g block of needle mesh so that similar blocks can be stitched at the end of each other to form the entire needle shaft. Needle knots are v i r t u a l points (for interaction model) l y i n g at the interfaces of these mesh sections. T h u s , each needle segment, w h i c h connect consecutive needle knots, represent a line passing t h r o u g h the center a needle section.  (a) Simple asymmetric  (b) Symmetric  Figure 3.9: M e s h topology of a section for two different designs.  Lateral axis Z Lum]  Lateral axis Z [fim]  // / 0  /  \\ V -2.5  -2  -1.5  -1  -0.5  0  0.5  1  1.5  2  2.5  Lateral axis X [^m]  (a) Simple asymmetric (for 20,60,...,180 sections)  -2  -1 Lateral axis X [^m]  (b) Symmetric (for 45,65,...,225 sections)  Figure 3.10: L a t e r a l t i p deflection w i t h respect to number of sections.  3.6 Discussion  3.6  46  Discussion  A method for generating a volume mesh from contours on parallel slices has been presented. Considering many seeds are implanted in the prostate, a finer mesh is preferred in it. Although only the prostate is defined for fine meshing in the current process, the entire P T V and the soft tissue in front of it can also be meshed fine considering the frequent needle interaction. The values for node culling threshold and preferred element size were found empirically. To achieve a small number of nodes for high speed in haptics, various applications and techniques were attempted. The final prostate mesh generated which has a reasonable number of nodes was shown to have good-quality elements for F E A . Both the F E A in Chapter 2 and the mesh generation scheme presented in this chapter are applicable with more than two tissue regions. Thus, the segmentations from different medical imaging modalities can also be used to generate the meshes for other anatomical structures inside the volume of interest. These meshes can be incorporated into the simulator. This could possibly bring issues especially for the bladder and the pubic arch. The bladder can be modeled as a fluid pocket [42]. Generally, the bones can generate problems because they have relatively high stiffness and represent fundamental boundary conditions. The high stiffness results in large force discontinuities when a node is snapped (see Section 2.6). Therefore, the pubic arch may require special treatment (such as solid contact models) out of the F E A . The presented F E M asymmetry problem on the needle mesh was overcome by a new design. It was shown to be symmetrically modeling the bending of a cantilever structure using a small number of optimally structured elements throughout its volume. This finding also made us aware that a simple uniform tetrahedralization of the tissue could result in asymmetric behaviour. That was why uniform meshes were not considered as the first tissue meshing solution.  Chapter 4  Ultrasound Synthesis According to the instantaneous deformed layout of the tissue, T R U S images are to be generated and displayed to the user. This requires proper deformation and interpolation models for fast realistic image synthesis, accounting for the characteristics of ultrasound (e.g., speckle).  Seeing  the needle reflection on these synthesized images is also important for the physician in judging where the needle path lies. Therefore, in the following sections we propose a T R U S slicing method with an integrated needle reflection model. A typical brachytherapy procedure was summarized in Section 1.2. Below the B C C A study during which our dataset was acquired is described.  4.1  Patient Data Collection  For a research study at the B C Cancer Agency, 20 prostate brachytherapy cases were to be scanned in C T and M R I right after each procedure. The data collected from these cases would be beneficial for us in segmenting all the tissue types around the prostate.  These regions can  be incorporated into our simulator. Moreover, the segmented seed locations from these cases could be used to validate the designed simulator. Grabbing the ultrasound data for the same cases would also enable us to validate other aspects of our work. For example, the seed and the needle shaft overlays onto T R U S could be evaluated. The T R U S image of the deformed region during an insertion could also be compared with the synthesized B-modes. Thus, we captured intra-operative ultrasound data using high resolution video hardware.  In order to know the  instantaneous stepper position (consequently, the slice depth) of a T R U S image and the instantaneous motion of the needle base which are needed for the above mentioned validations, we recorded a frame including the precision stepper, the needle, and the physician's hand holding onto it using a digital camcorder. This data could also help us verify needle flexibility and deflection models. Furthermore, tissue elasticity parameterization (elastography) and/or verification of its known biomechanical values can be performed on this captured data. The study conducted in the B C C A aimed to compare the dosimetry calculations using USC T (Computed Tomography) fusion with the ones obtained from M R I - C T (Magnetic Resonance  47  4.2 The Shift and the  Torque  48  Imaging) fusion. T h u s , unlike t y p i c a l procedures for these 20 cases prostate ultrasound volume studies were conducted at the end of each procedure. T h e y segmented the urethra i n a l l m o d a l ities and registered it using VariSeed for b o t h fusion types above. A catheter w i t h contrast fluid is used to make u r e t h r a visible i n T R U S images. B y means of this registration, the seed locations easily found i n C T (see F i g u r e 4.1) can be projected on M R I where soft tissue can more easily be segmented (see F i g u r e 4.2). We also recorded this volume study, w h i c h gave us T R U S slices scanned i n 5 m m increments w i t h a l l i m p l a n t e d seeds from the base to the apex.  (a)  (b)  (c)  Figure 4.1: C T images after a brachytherapy procedure.  4.2  The Shift and the Torque  T h e forces applied by the needle cause the tissue to displace a n d deform. T h i s is m a i n l y observed d u r i n g the insertion w h e n the forces on the tissue are higher because of the t i p c u t t i n g force added on the shaft friction. A l t h o u g h the insertion and the retraction m a i n l y deforms the region parallel  4.2 The Shift and the Torque  49  to the needle axis, the mesh geometry a n d the b o n d i n g between a n a t o m i c a l structures m a y also cause more complex behaviour (e.g., sagittal torquing), details of which are discussed later i n this section. A l t h o u g h a sagittal i m a g i n g o p t i o n is available i n the T R U S hardware, the images are mostly acquired transversally to observe the prostate cross-sections while the stepper that positions the probe is moved w i t h 5 m m increments.  T h e instantaneous ultrasound image represents  tissue content residing m o m e n t a r i l y i n a slab. M e a n w h i l e , when the tissue deforms d u r i n g needle insertion, a n anatomy i n a separate region can possibly shift into this slab. T h i s effect is easily observed especially at the base of the prostate, where the T R U S image before the needle insertion is m o s t l y covered w i t h the bladder a n d it has almost no prostate tissue i n it (see F i g u r e 4 . 3 ( a ) ) . A s the needle penetrates into the prostate or while it is advancing through, the prostate base w h i c h basically just shifts i n 3 D suddenly appears i n the image (see F i g u r e 4 . 3 ( b ) ) . F i g u r e 4 . 3 shows two images of the same d e p t h w i t h «1 s t e m p o r a l separation.  T h i s means the stepper  (consequently the probe) is not moved between these two images, so their difference comes from the shift of the tissue i n this one second. I n the first image, the physician is not a p p l y i n g any force o n the needle. T h e n , when he pushes the needle, the prostate base can be easily seen i n Figure 4.3(b) (which is a c t u a l l y the same slice but just the prostate shifts into i t ) . Note that the needle t i p deflects upwards w i t h the torquing prostate. T h i s torquing effect is described below.  Bladder  Soft Tissue  (a) Relaxed tissue.  (b) Force applied on the needle.  Figure 4.3: Prostate shifting into the bladder d u r i n g a needle insertion. A n o t h e r observable deformation phenomenon occurs when T R U S is viewing a slice near the m i d d l e or the base of the prostate.  If the needle is pushed meanwhile, the observed prostate  cross-section moves upwards up to 1 cm. T h i s is again because of the torque, which the advancing needle generates around the stiffer tissue connecting the prostate to the pubic arch. For simplicity, the p r o b l e m can be studied i n 2 D sagittal direction without loss of generality. In 2 D , the T R U S imaging planes look like lines as shown i n the simplified sketch i n F i g u r e 4 . 4 . T h e physician can adjust the stepper to move the T R U S back a n d forth to see different slices. T h e medical n a m i n g convention for these T R U S slice depths are "5ase", Base+0.5", u  etc. L e t these lines (which are planes i n 3 D ) be Z / i , Z / 2 , • • •, L  n  "Base+l.O",  before any deformation (needle  insertion). Note that after deformation the T R U S probe w i l l obviously still be acquiring planar  4.2 The Shift and the Torque  50  (linear in 2D) data whereas the tissue contents of the voxels on those planes will not be the same as Li anymore. 5  mm  I  I  k  Bladder  #  1  I  ,  ,  |  |  ;  /  Prostate N,^  ,  y <o  1  '  f «J  ;  1  :  s\!  IX)  ' L, Figure  1  L  i  2  L  | Rectum  3  •••  L  n  4.4: A 2D sketch of the region and the T R U S images before needle intervention.  When the needle is being inserted into the prostate, the structural deformation can be demonstrated as in Figure 4.5.  With the friction and the cutting forces of the needle, the prostate,  which is relatively stiffer than the surrounding tissue, shifts axially and rotates sagittally (due to the stiff tissue above its apex), such that Q will shift and occupy f T at a time instant.  Needle  L Figure  4.5: Demonstration of the structural deformation during needle insertion.  The first observed phenomenon, the base shift into the bladder, can be easily visualized considering the L\ in Figure 4.4: while we do not normally see the prostate in L i , after the shift the base will be visible in it. Let us also analyze the second observed phenomenon, the upward image shift due to the torquing effect. Let L be the slice imaged by the T R U S probe. Assume that before the needle insertion the pixel P of L is inside the bladder. Consequently, let us say that the physician sees P almost black (due to the fluid filled bladder). While the needle is being inserted, the sagittal torquing may have the prostate occupy the coordinate of P. If this happens, P will then show prostate tissue (since P e l l * ) changing the T R U S image look. In  4.3 Problem  Definition  51  other words, the cross-section of Q is then be shifted upwards showing the cross-section of Q*.  4.3  Problem Definition  Due to the above mentioned deformations, the i n i t i a l l y observed lines L , that denote the tissue content at those locations change their shape d u r i n g the insertion. a 3 D sketch  1  T o visualize this better,  of the tissue w i t h a n advancing flexible needle is shown i n F i g u r e 4.6.  Here a  representation of how the T R U S planes may deform is given. W e w i l l refer to the deformed geometries of these T R U S slices as L J , L\,...,  L * . I n this chapter, the superscript "*" is used  for the geometries i n deformed mesh layouts.  Figure 4.6: A sketch representing the 3 D deformation of the i n i t i a l T R U S images. N o t e that whereas L * can be curved because of the deformation, any instantaneous  TRUS  imaging plane L w i l l always be linear (and similarly planar i n 3D) since the probe acquires planar images regardless of the tissue content. T h e ultrasound d a t a o n slices Li w i t h a 5 m m separation is acquired d u r i n g the preoperative T R U S volume study. T h u s , the voxel intensity values of the i n i t i a l mesh is k n o w n on those planes (Li).  Let us assume that "the voxel intensities of L j are representative for the corresponding  voxels of L * i n the deformed state". T h i s is not always true, since U S image pixels are i n fact a cumulative effect of scattering and reflections.  T h i s refer to two aspects of the u l t r a s o u n d  imaging, the speckle f o r m a t i o n a n d the i m a g i n g directionality. T h e former affects us i n a way that the same voxel (tissue content) shifted to another part of the volume is not necessarily be denoted by the same image intensity on an U S image.  For example, comet t a i l w h i c h is  observed b e h i n d most of the metallic surgical tools is formed along a line o r i g i n a t i n g from the probe surface, therefore it is orientation dependent. A l s o acoustic shadowing of different density anatomical structure depends on the angle of U S incidence. A s a result, the intensity of a point may not convey m u c h information about the tissue. Nevertheless, p r i o r work [86, 84] achieved success using interpolation for reslicing.  However,  This figure does not depict the mesh. The lines on it are only for representing the TRUS planes and the general structure of the tissue.  4.4 A Direct Approach  52  similar studies in literature have dense input datasets as opposed to our 5 mm inter-slice spacing. In brief, the problem can be formalized as "determining the ultrasound image on L at any instance of time with the ultrasound data given on the initial slices L j which are deformed into L* at that instance".  Two approaches and their comparison are presented below.  The  reader is referred to Appendix A.4 for information on basis functions and discretization which are mentioned in this chapter.  4.4  A Direct Approach  A simple instance of the problem is sketched in Figure 4.7, where the imaged plane L falls completely in-between two adjacent deformed data slices. intensity on any point A on L, the geometry of L* and L*  Thus, in order to interpolate the +l  are needed. The points of these  slices can be found using Algorithm 4.1.  Figure 4.7: Demonstration of the structural deformation during needle insertion.  A l g o r i t h m 4.1 Finding the locations of the input data after deformation. 1: for x 4- each voxel on Li do 2:  Determine the mesh element Q where x € Cl  3:  Find the basis functions Nf for the element ti  4:  Using the corner displacements u\, find deformed position x*  e  e  e  This algorithm is suitable for various optimizations. For instance, on line 2 of Algorithm 4.1 the elements containing the data points are needed, and this can be calculated offline and put in a look-up table.  In other words, in an offline step all input data points, which are on the  planes Li with depths known prior to the simulation, can be tagged by the element numbers  4.4 A Direct  53  Approach  they reside in. In addition to this, the basis functions already found in offline K  e  computations  can be stored in arrays for fast access on line 3 of Algorithm 4.1. Determination of whether a point is inside an element or not is performed using the function in Algorithm 4.2. A n instance for one loop of this algorithm, where X is in the element can be seen in Figure 4.8. Note that it requires the local numbering of the elements to be in  right-hand  order as described in Section 3.4. A l g o r i t h m 4.2 The function for deciding whether a point is inside a tetrahedron, l: function I S P O I N T I N S I D E E L E M E N T ( X , e) 2:  > For each face of the element  for all Face G e do {Pi, P , P }  4:  N* <—  5:  if ( P )t •N <0)  m  n  X pTti >  t  6:  e n d for  7:  r e t u r n True  > P's ought to be in right-hand  <r- corners of Face  3:  order  > Find the surface normal > Dot product finds the projection  t h e n r e t u r n False  Figure 4.8: Determining which side of the plane PiP P m  n  the point X is. This is one of the four  checks before X € 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  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:  I \B A\+I \BjA\  =  Bi  i+1  Bi+1  \BiA\ + \B A\  A  i+1  [  ''  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 disadvantages 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, — 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: 3:  for e* <r- 1 to #e do if I S P O I N T I N S I D E E L E M E N T ( X * , e*) = True t h e n break  4:  Find basis functions Nf  5:  X <- X* +  YA=I  f  N  for the deformed element e*  ( )(~ i) x  u  >  t> using Pf  as x in (A. 18)  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 Base+2.0" u  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: • 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 • 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) • 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  Image Interpolation  Results  For the inverted problem, the ultrasound was synthesized using the fast interpS function of Matlab. 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 Base+2.0" u  probably between the initial images of Base+2.0" u  is an interpolation  and Base+2.5", considering the shift of the u  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 artifact 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 interpolation was minimized using cubic interpolation as demonstrated in Figure 4.11(d). O n 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), a n d 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 (zti ue SS  ~  1.54 x 10 Rayles, 6  steel  z  ~ 47 x 10  6  Rayles, where Rayle is kg/m s). 2  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°). 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 u l t r a s o u n d images were used. T h e y were acquired at the h o s p i t a l i n a p h a n t o m using the same u l t r a s o u n d system a n d probe settings that are used d u r i n g the procedure. A n acquired image can be seen i n F i g u r e 4.13(a). F i r s t its b a c k g r o u n d was removed by s u b t r a c t i n g the image of the region before needle insertion, then the r e m a i n i n g noise was e l i m i n a t e d using a n adaptive thresholding filter. T h e processed "clean" version is shown i n F i g u r e 4.13(b).  (a) raw  (b) filtered  (c) overlayed  (d) real  Figure 4.13: Needle shaft. In the code, we are only s i m u l a t i n g the ellipsoid overlay at the proper position a n d orientation. A l t h o u g h m u l t i p l e reflections are relatively easy to implement, comet t a i l synthesis needed further in-depth study. Since this thesis is not focused on U S artifact m o d e l i n g , better artifact synthesis schemes are yet to be researched i n future work. For the shaft image, the p o s i t i o n where the needle intersects the i m a g i n g plane is easily found using the needle knots just before a n d j u s t after the i m a g i n g plane. T h i s intersection gives where the needle should be positioned o n the T R U S . W i t h this i n f o r m a t i o n a n d using the p r i o r i - k n o w n coordinates of the probe, the angle of incidence can be calculated.  T h e cleaned shaft image is then rotated accordingly a n d its  brightness is adjusted to m a t c h w i t h the T R U S image. Since adaptive thresholding filter leaves sharp edges, it is also smoothed w i t h a 7 x 7 G a u s s i a n  filter.  F i n a l l y , it is added on to the  synthesized u l t r a s o u n d image where its brightest p i x e l w o u l d correspond to the desired shaft center coordinate. T h e final generated image can be seen i n F i g u r e 4.11(c) w i t h the overlayed needle between 2B a n d 2.5b.  4.8  Discussion  In this chapter, two methods for generating T R U S images of soft tissue, the deformation of w h i c h was given by F E M , have been presented. T h e first one (in Section 4.4) has a n advantage of p r e c o m p u t i n g basis functions a n d tagging the i n i t i a l l y k n o w n voxels w i t h the F E M elements they reside i n . However, its i n t e r p o l a t i o n requires m a n y geometrical computations w h i c h w o u l d make its i m p l e m e n t a t i o n 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 i n 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 approximated) 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 « 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 « 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 N  th  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 T R U S 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 T R U S 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 G i D . 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  64  Initialization  (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. Therefore, 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 Name  Initialization  65  Dimension  Content  Xn  n x 3  Nominal positions of tissue nodes, (its order sets the global numbering)  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.  x 2  Shaft friction and tip puncture forces for each node.  Fd  n  TissuePar  t x 2  NEconn  n x  Parameters (E, v) of each tissue type.  v  The indices of elements surrounding each node.  b  Nfixed Kinv  The indices of the nodes with displacement boundary condition.  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~  as in Algorithm 5.1. On Line 11, this huge matrix is inverted using the Matlab's  l  inv function. A l g o r i t h m 5.1 Creating the inverse stiffness matrix offline.  1: function CREATElNvSTlFFMAT(TissuePar, ElemType, X n , NodeList, 2: TissuePar(S, v) —> Lame constants (A, u) —> C 3:  Kglobal <r- 0 „  4:  for  3  i i—  > (A.12) and ( A . l l )  x 3 n  1 to e do  5:  Find element volume v  6:  Find B  > (3-1) > (A.18) and (A.24)  e  7:  #12x12 <"  8:  A d d Ki2xi2  9:  Nfixed)  vB CB eT  o n  *° t  e  n e  > (A.30) corresponding elements of Kglobal  > (A.32)  for all i € Nfixed do  10: 11:  K'  5.4  Needle Model Initialization  the rows/columns {3i — 2,3i — 1,3i} of Kglobal ^ 1  h3 X  > (A.33)  <- (Kglobal)'  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 Name  66  Dimension  NXn  nx  NNodeList  e x 4  3  Nominal positions of needle nodes. Indices of each needle element's corner nodes.  3n x 3n  NKinv  Content  The inverse global stiffness matrix for needle F E M . Table 5.2: Needle F E A variables.  5.5  A n Iteration  The flowchart of a code iteration is given in Figure 5.2. Its functions are detailed as follows:  LOCATE NEEDLE  x  <r-  \if ,u.w) w  SETCONSTRAINT ANGLES  MODIFYMESH  Ky  SERVICEFINITE STATEMACHINE  w  <- (x,y)  Set contact node friction thresholds  SETCONSTRAINTS  PROJECTFORCES INTEGRATEFORCES  SOLVENEEDLE  ADDNODE  I  f-  Xn +u  A d d tip bevel force Figure 5.2: A main iteration flowchart of needle interacting with the tissue.  • L O C A T E N E E D L E : The needle base position and orientation are taken from the graphical interface (for haptics they will be read from the haptic simulator device). Adding Uk, the needle knot world-frame positions are computed. Comparing the tip coordinate with the previous iteration's, the needle motion direction (advance/retract) is set as well.  (This  direction is used for the stick-slip automatons.) • SETCONSTRAINTANGLES:  A S also given in Algorithm 5.2, first using the needle knot  positions, the needle segment orientations are calculated geometrically (atan2 function 1  was employed). These are necessary for local frame updates in Section 2.4. Then, for each The needle geometry is rotationally symmetric. The tip bevel disobeys this axial symmetry, however, the bevel is not a concern of F E A either. Thus, 8 is taken constant which saves many operations during transformations. lr  5.5 An Iteration  67  contact node j the cumulative rotation matrix Rj- j from the needle frame to the tissue N  frame is found, where i denotes the current iteration. matrix Rx x _ i  i  ii^^TN i-i)  is set by calculating  1  Finally, an incremental rotation  • The incremental rotation matrices  are used for Kw update at each iteration and the cumulative one is used when a new node is added. To add a node, KH matrix is reformed from K~  x  which resets the orientations  to the initial geometry. Therefore, the initial rotated working matrix Kw is regenerated by updating KH using cumulative rotations. A l g o r i t h m 5.2 Constraint angle calculation for a contact node. > found geometrically  1: F i n d {(f>, 9, *p)segment  2: for all j G contact nodes d o Find R  3  3 :  j ^TiTi-i  •  *~  j  TN  {  ( j \~ ^TN,iy^TN,i-lJ  X  CHECKNEWELEMENTINTERCEPT:  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. •  MODIFYMESH:  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 K  e  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. • ADDNODE:  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^ to the condensed system matrix K . 1  H  the L C F s using cumulative rotations RTN.I  a n  Finally, adjusting  d switching the boundary conditions {x,y,z  for stuck; only x and z for slip), the working matrix Kw is obtained.  •  SERVICEFINITESTATEMACHINE:  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 A l -  gorithm 5.3.  •  SETCONSTRAINTS:  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 t h e n 2:  Update L C F for stuck using  3:  if "force along y" > "friction threshold" t h e n  R ,,  > (2.11)  3  T T  4:  Boundary condition change along y  5:  Set its state to "slip"  6: if slip t h e n 7:  if TipMotionj ^ TipMotiom_i t h e n  8:  Boundary condition change along y  9:  Update L C F for stuck using  10:  R .  t> (2.11)  J  TT  Set its state to "stuck" else  11: 12:  Update L C F for slip using  R^,. .  > (2.19)  T  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). •  PROJECTFORCES:  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. •  INTEGRATEFORCES:  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. •  SOLVENEEDLE:  After distributing the knot forces (dividing them by 3), the needle model  is iteratively solved ( F E A ) 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  CHECKNEWELEMENTINTERCEPT).  Moving this step to a second  (visualization) computer and its drawbacks are discussed in Section 5.7.  5.6 Notes  5.6  69  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 ® ) . That structure allows the use of WindRiver's T o r n a d o ® d e v e l o p m e n t environment and the control of runtime variables in Matlab S i m u l i n k ® t h r o u g h the Real-Time W o r k s h o p ® t o o l b o x , 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 B  e  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(n ) 2  operations (mesh adaptation), therefore all the other optimizations are  given a lower priority. The presented (D(m n) approach decreased the 0(mn ) 2  2  mesh adaptation  time of two sample meshes with 836 nodes and 240 nodes from « 1 3 s and w l s down to « 1 . 2 5 s and « 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 Figure 5.3(b), note that the computations take 2-3 s during mesh modification using 0(mn ). 2  The  effect of the different stiffness and the friction of prostate tissue can be observed in Figure 5.3(a) around the 380  100  150  t/l  iteration.  200  250  300  350  400  450  ?00  150  200  250  300  350  400  Iteration number  (a) Needle base forces (from top - x, z, y).  Iteration number  (b) Time taken by each iteration.  Figure 5.3: The simulation forces and timings.  450  70  5.7 Discussion  5.7  Discussion  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. both moved to a separate computer (let us call it  They can be  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. positions can be sent back by the T R U S computer.  Since there only a few of those nodes, their  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 E x t e n s i o n : 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" M e s h e s : Using the presented scheme, a 3D volume can be meshed with conforming region boundaries on surfaces generated from contours on parallel slices. This method achieves relatively high-quality meshes which are also suitable for fast haptic implementations of the F E M . B r a c h y t h e r a p y S i m u l a t i o n : 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 T R U S 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  6.2  72  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. O n 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 converge require large matrix inversions ( K " which is rank « 5 0 - 1 0 0 ) at each iteration. 1  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 approximation 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 implementations 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 International, (Chicago, IL), pp. 748-752, 2000.  22 Annual EMBS nd  [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, "Measurement of the ultrasound backscatter signal from three seed types as a function of incidence angle: Application to permanent prostate brachytherapy," International Journal of Radiation Oncology*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 8 Medicine Meets Virtual Reality (MMVR), pp. 263-274, 2000. th  [6] F. S. Azar, D. N. Metaxas, and M . D. Schnall, "Methods for modeling and predicting mechanical 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 brachytherapy 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, "Cryotherapy simulator for localized prostate cancer," in Proceedings of ltf Medicine Meets Virtual Reality (MMVR), 2002. h  [9] R. Alterovitz, K. Goldberg, J. Pouliot, R. Taschereau, and I.-C. Hsu, "Sensorless planning for medical 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 models 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," Visualization 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 6 Medicine Meets Virtual Reality (MMVR), (San Diego, CA), pp. 316-322, 1998. th  73  BIBLIOGRAPHY  74  [16] S. P. DiMaio and S. E . Salcudean, "Needle insertion modeling and simulation," IEEE on Robotics and Automation, vol. 19, no. 5, pp. 864-875, 2003.  Transactions  [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 15 Southern Biomedical Engineering Conference, pp. 170-173, 1996. th  [18] G. Hirota, S. Fisher, and A. State, "An improved finite element contact model for anatomical simulations," 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, ImageGuided 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 intraoperative 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, M E R L - 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 displacements 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 26 Annual Conference on Computer Graphics and Interactive Techniques, pp. 65-72, 1999. th  [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 SIGGRAPH99 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 deformable body simulation using dynamic progressive meshes," Computer Graphics Forum (EUROGRAPHICS'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 10 Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems (HAPTICS), pp. 344-351, 2002. th  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 12 Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems (HAPTICS), pp. 266-271, 2004. th  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, (Washington, 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 Proceedings 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 9 International Symposium on Experimental Robotics, June 2004. th  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 models for fluid flow and heat transfer phenomena in electronic applications," in 2,0 Annual IEEE Semiconductor Thermal Measurement and Management (SEMI-THERM), pp. 278-285, 2004. th  [54] C. Singh Verma and V. C. V. Rao, "Parallelization of finite volume computations for heat transfer application using unstructured mesh partitioning algorithms," in 4 International Conference on High Performance Computing, 1997. th  [55] F. Ganovelli, P. Cignoni, C. Montani, and R. Scopigno, "Enabling cuts on multiresolution representation," 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 7 Workshop on Bildverarbeitung fur die Medizin, (Erlangen, Germany), pp. 279-283, SpringerVerlag, 2003. th  [58] P. Fleischmann, R. Kosik, and S. Selberherr, "Simple mesh examples to illustrate specific finite element mesh requirements," in 8 International Meshing Roundtable, pp. 241-246, 1999. th  [59] R. V. Garimella and M . S. Shephard, "Boundary layer meshing for viscousflowsin complex domains," in 7 International Meshing Roundtable, pp. 107-118, 1998. th  [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. 19651990, 1984. [61] R. Lohner and P. Parikh, "Three-dimensional grid generation by the advancing front method," International 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 Matematicheskii 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 Computational 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 4 International Meshing Roundtable, pp. 3-14, 1995. th  [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 7 International Meshing Roundtable, (Dearborn, MI), pp. 239-267, 1998. also available at http://www. andrew. emu. edu/user/sowen /survey / . th  [71] N. Molino, R. Bridson, J. Teran, and R. Fedkiw, "A crystalline, red green strategy for meshing highly deformable objects with tetrahedra," in 12 International Meshing Roundtable, pp. 103-114, 2003. th  [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. Waterton, "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 threeand 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 . resultant force and moment vectors, respectively, on A A .  Let A F and A M be the  Their intensities on the area can be  found at the limit as follows:  i, lim A^o A  A  where T is called  stress vector  or  AM  _  " A X ~~  traction,  and  dM  HA  C is  _  ft  ^  '  ~ called  couple-stress vector.  theory of elasticity proceeds on the assumption that C = 0 [87].  The elementary  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 T ,T ,T . x  e.g., T  x  y  = o x + o y + o z, xx  xy  xz  z  If they are written in terms of their Cartesian coordinates, 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 ®~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, €ij = eji, and the vector representation of strain becomes:  (A.5)  e =  A.2  Deformation  Displacements of a point x = (x, y, z) can be written as u — [u(x) v(x) w(x)] in the reference T  coordinate frame of the undeformed body. Then the strain components are given by:  _ du 1 ~ dx' 2  S x x  du dx  +  _ du dv ~7\ h dy dx  €xy ~  +  dv dx  du du h dx dy  ~x  (A.6)  dv dv dx dy  dw dw dx dy  (A.7)  The other four strain components are defined similarly. Note that the square-bracketed secondorder terms are relatively small for small displacements. The matrix expression of the  engineering strain tensor in terms of u becomes:  du , 1 dx ^ 2 dv_ _i_ dy ^  2 i1  yy  e  dw dz ^ 2  €(«) = yz  dv dz  e  I  1/ 9 u \  2  dw_ dy  i dw \* \dx>  , / dv>  \dy>  ^ du d-z)  I  " Ady> r  +  - )  dz  du du I 9»3» j dy dz  dy dz  2  dz > dy dw dz dw  €zx  du dz  djv dx  du du dz dx  I  dv dv_ dz dx  xy  du _,_ dji dy dx  du du dx dy  i  dv_ dv_ , dw dw dx dy dx dy  e  du dx  dv dx  dw dx  0  0  0  0  0  0  d_ dy  0  0  0  0  du dy  dv dy  dw dy  0  0  0  0  0  0  0  0  0  d_ dz  d_ dz d_ dy d_ £ dx  0  0  0  du dz  dv dz  dw dz  du dz du dy  dw dz dw dy  du dz du dy  dv dz dv dy  dw dz dw dy  0  0  0  du dx  dv dz dv dy dv dx  dw dx  du dx  dv dx  dw dx  0  0  0  0 0 d_ dz d_ dy  , dw_ dw dz dx  0  dx  0  I  du\^ I I dv_\^ dx' ~*~ V dx )  0  U  V  +  W  d_ dx  Bu + -A® This vector cr(u) is called the  .  du dx dv dx dw dx du dy dv dy dw dy du dz dv dz dw dz  (A.8)  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:  <r(«) =  £dx  o  o  0  d_ dy  o U  d_ dz  0 0  d_ dz  -2-  0  dz d_ dy  V  a_  dy d_ dx  m  W  0  = Bu This is called Cauchy's infinitesimal  (A.9)  strain tensor. It is rotation variant, therefore must be  used with care. This linear approximation will be used for the tissue model in this thesis since the tissue deformation and displacements during prostate brachytherapy is assumed to be small.  A.3  Linear Elastic Behaviour  Elastic materials load and unload along the same stress-strain curve. If they follow a linear path, they are called linear elastic materials. A Hookean elastic material obeys the generalized Hooke's law: a = Ce  (A.10)  where C is the material matrix relating the strain to the stress. For homogeneous and isotropic materials C is defined as:  C =  X + 2p  A  A  0  0  0  A  X + 2p  A  0  0  0  A  A  X + 2p  0  0  0  0  0  0  i  0  0  0  0  0  0  0  0  0  0  1  (A.ll)  0 0  where the Lame constants A and p are expressed as follows: A =  Ev  E  (l + i / ) ( l - 2 i / )  2(1 + 1/)  (A.12)  Here E and v are called Young's modulus and Poisson's ratio, respectively. The first one is the relation between the normal stress, o , and the normal strain, e . The latter is the incompressx  x  ibility of the material. In other words Poisson's ratio, u, is the lateral contraction per longitudinal expansion.  A.4  Discretization  83  A boundary-value problem of a linear elastic object Q consisting of particles x = (x, y, z) can be defined as setting force or displacement constraints at some points and trying to solve for the displacements and forces of all other element particles. For the solution of such a problem, the strain energy of the elastic body is a convenient point to start from. It can be expressed as:  E.strain  1 I e o dx Jn  (A.13)  T  2  where e is the engineering stress vector and a is the engineering strain vector. Substituting the Hooke's law in (A. 10):  Egtrain  ^J f e Ce dx  (A.14)  T  Jn  2  Then the strain energy stored in the body for the displacements u(x) = [u v w]  T  can be  written using e(tt) = Bu from (A.9) as: E train(u) s  = \ f (Bu(x)f Jn  C\Bu(x))  dx  (A.15)  1  Assuming the potential energy of an object at rest is E(0), its strain energy and the work done on the body particles by the external/internal forces, f{x), can be related using the conservation of energy as follows: E(u) = \ f {Bu) C{Bu) Jn  \ f f udx Jn  dx-  T  Note that {Bu)  T  A.4  ^ u B T  T  (A.16)  T  2  2  since B is an operator matrix.  Discretization  Since an analytical continuous solution of (A. 10) over a complex shape is computationally impossible, one possible approximation can be the discretization of the object into simpler geometrical structures. The behaviour of the body, Vt, can be modeled by dividing it into the elements Q  l  where i = 1 , . . . , s as in Figure A.3 having the nodes P where q = 1 , . . . , n. Although those q  elements ought not necessarily have the same structure, assume that they are all tetrahedrons in our implementation. Considering an element, Q , shown in Figure A.4 its corners are expressed e  as Pf where i is the local node number within the element and it is completely unrelated to the global node numbering. While the i  th  node displacement is itj = [ui V{ Wi] , the compound displacement vector can  be written as u — [uj u^...  T  u^] '. T  For the tetrahedron element, fl , the compound element e  displacement vector becomes u = [{ul) e  T  ( « 2 ) { t) T  u  T  {ul) ] '. T  T  The displacement of a point inside an element can be estimated as a weighted average of its corner displacements: 4  u(x) = Y,Nt(x)nt  (A.17)  A.4  Discretization  84  Figure A.4: A tetrahedral element and its corners in local numbering.  where Nf are called the basis functions. A linear displacement interpolation satisfies zero-order continuity between elements. Although it was reported that a quadratic interpolation gives much better F E A results, for our purposes the accuracy is given up for the speed and the simplicity of the linear estimation. The linearly interpolating basis functions are called natural coordinates. Natural coordinates should be a function of the position, x, since the weighting of displacements should depend on where the internal node is. For example, if the node is closer to a corner, then this corner's displacement will determine the node's displacement the most. A n d also the total displacement shouldn't be magnified which brings the constraint  Nf = 1. Consider the relocation of the  corners onto the origin which requires negative displacements of their 3D coordinates. This would confine any internal point to the origin as well, resulting in a displacement equal to the negative of its coordinate. Using this fact, given the element corner positions [xi yi Zi], the position of an internal point [x y z] can be estimated by the natural coordinates as follows: 1  1  i  X  Xl  y  yi  z  mix)  1  1  x  2  •'•'3  £4  1)2  :'/s  P4  N§(x)  ~-2  ~3  Z4,  _ Nl(x) _  Nf(x)  1  Nl(x)  X  N§(x)  ^3  4  y z  (A.18)  A.4  Discretization  85  where:  b\  1  d\  d%  b\  1  1  1  Xi  X2  Xz  yi  yi  ys yi  Zl  Z  Zz  2  X  4  (A.19)  Z4  So the linearly interpolating basis functions for a tetrahedron is expressed as:  Nt{x) = a +b x e  e  i  To calculate Bu(x)  i  + c y+ d z e  e  i  (A.20)  i  in (A. 16) inside an element e, the partial derivatives of the discretized  displacements in (A. 17) are needed:  du dx  i=l  du  E  dy  i=l  du ~¥z dw  A n d similarly | ^ and ^  i=l N  i( ) x  V^' e e  e  (A.21)  i=l  y  s—^  i=l  1=1  are obtained:  4  dv dx  i=i  dv  4  dy dv dz~  <9a;  £  - E* 4 -E i=i  (A.22)  i=l  i=l  4  1=1  i=l  In (A.21) and (A.22) all the partial derivatives are constants. Therefore this constant  Bu(x)  can be rewritten inside an element as:  n  0  0  Bu(x)  = J ^ - Ri i=i  i  a  where B\  0  0  0  0  <  dt  0  c!  0  dt  e  0  (A.23)  A.5 Element energy  86  If this is expressed in one matrix notation using the 1 2 x 1 compound element displacement vector, u , Bu(x)  becomes B u  e  e  for any x inside a given element e where:  e  bl  0  0  b\  l  0  0  0  0  dl  0  dl  dl cl  0  0  0  b% 0  0  2  0  0  3  0  0  0  dl  0  0  dl  cl  0  dl  4  0  bl  dl  0  b  b\  0  4  b\  0  c  C  e 2  0  0  0  c  0  0  0  0  d% c% 0  d\  dl  0  0  C  bl c  0  4  d% e  61  4  (A.24)  r4 c  0  Now the strain energy in (A.15) can be rewritten as a sum of element energies using the property of definite integrals, J  ( )dV — Y2 fv ( e  v  E train(u) s  =JT  [  a s :  e  (B u ) e  e  T  dx = J V  C (B u ) e  e  Note that whereas B is an operator matrix, B  e  (B u ) e  e  = u  T  e T  B  e T  e  =  e  s  \ j  A.5  eT  eT  e  e  dx  (A.25)  is a constant matrix. That's how the equality  U B CB U eT  Ju  2  e  u B CB u  holds. So the strain energy of an element can be expressed as:  E train(u )  where v  f  e  eT  e  =  -u B CB u  =  ^u B CB u v  eT  eT  eT  e  eT  f  e  e  dx  e  e  dx (A.26)  e  is the element volume.  Element energy  Considering the tendency towards the minimum potential energy, (A.16), in the tissue, the solution of the strain problem occurs when the variation of the energy is zero, 5E(u) = 0. Then the discretized form becomes:  SE(u) =  SE (u ) e  e  =0  (A.27)  e  Considering J2 ^E (u ) e  e  e  = Y^, f^i"^ e  6  a  n  d each node displacement, uf, is independent from  the others, the following is concluded:  dE V j e n e :  As the element energy definition, E , e  e  du^  =  0  ( A  -  2 8 )  the equation (A.16) can be written using the element  A.5 Element energy  87  strain energy expression in (A.26) as follows: E {u ) e  where f  = \u  e  e T  B  e T  CB u v e  e  e  - f  (A.29)  e  is the 12 x 1 compound element force vector. Applying vanishing criteria in (A.28) on  e  this expression using the algebraic fact that -J^x Ax  = 2Ax gives:  T  3  ^u B CB u v \2 d (u B CB u ^j 2 du v B CB u eT  eT  du v  e  e  - f]  e  6  e  =  0  =  0  =  0  =  f  e  eT  eT  e  - f  e  e  (  e  eT  e  - f  e  e  Ku e  where K  e  e  (A.30)  e  is the 12 x 12 element stiffness matrix. Calculating these for each element and then  combining them to form the entire set of system stiffness equations yields: Ku where K [fl  = f  (A.31)  is the global stiffness matrix, u is the compound node displacements, and /  f l • • • fn]  T  i * s  n e  3n x 1 compound force vector such that fi — [/f ff f?]  T  =  is the force  components at the node i. Combining the set of system equations of each element necessitates the mapping of the local element node numbering to the global node numbers: K  =  global (K ) e  (A.32)  e e  u  =  global (u ) e  Note that the displacements are not summed not to accumulate the displacement of the same node shared by the neighbouring elements, whereas their force contributions and thus the individual stiffness matrices are summed up. Since K is singular, the linear system in (A.31) cannot be solved. One way of solving is  A.5 Element energy  88  modifying K so that the displacements of some nodes are constrained as follows: K ku  ku  • ••  k  •  ki  k22  • ••  ki  • •  kn  ki2  • •  k--  knl  k2  • ••  2  u  2  u  kin  fl  k  h  2n  U ---  n  h • •  •  "•ni  k  ku  • ••  0  k2i  k  22  • ••  0 •  0  •••  1  • •  n2  • ••  0  •  n  •  kin  h  k r^nn  fn  '  kin •  k  knl  k  0  knn  diku  h~ • dik  2i  2 n  u  0  fl~-  (A.33)  =  di  fn ~~  Here setting dj = 0 fixes the node to its initial position. The necessity of this node constraining process can also be visualized intuitively: A body floating in space without any constraints can move to anywhere even with no force exerted. This would result in infinitely many solutions of the system. It can easily be seen that fixing 3 non-collinear nodes in space, constrains the body (both translations and rotations).  Although 3 nodes bring 9 displacement constraints,  it can analytically be proven that six independent constraints are enough to ensure that K is invertible.  Appendix B  Boundary Condition Change (2.6) gives us a set of equations which we can write the row i and the rest seperately as: hum + ... + kuUi + . . . + k u  fi =  ^jew '• fi  ir  kjiui + . . . +  =  kjiUi  + ...  +  (B.l)  r  kj u r  (B.2)  r  where N = {1, . . . , i — 1, i + 1, . . . , r}. J  The form in (2.7) requires fj and ui in terms of Uj and fi. If we find Ui from (B.l) and substitute it in (B.2) it becomes : l  => Ui= ^jeH  i  :  fi  =  - ^ - [ f c i u i + . . . + 0" + . . . +  fc u -/ ]  i  kj\u\ + . . . + 0  Jl  + ... +  ir  r  i  (B.3)  kj u r  r  k- J  1  [kilUi  +  ...+0  U  + ...+  kU ir  To be able to write it as a difference from the initial K~  jeN  fi]  and reorder the terms:  kjifi  V  -  matrix like in (2.8) we add and  x  subtract the term  r  ;  :  fj  =  kjiUi  + . . . + kjifi  + ••• +  kj U r  r  k-  ~r~~ [fca«i + • • • + (ku -l)fi  1  + ... + k u ] ir  r  (B.4)  To make theand U{ 0in' are (B.3) look like form should add ku and subtract all the The terms 0" usedalso to indicate the the absence of in the(2.7) termswe with the factors and kji, respectively. J  89  90 missing terms and regroup them properly:  Uj =  ki\U\ + . . . + kafi + ... + ki u r  T  — \ki\U\ + . . . + kafi + ... + ki u ] r  1 — 7— [knui + ... + 0  11  r  + . . . + k u] ir  +  r  1 —fi  hum + ... + kafi + ... + k ir r u  kii  J  =  + 1  [knUi + ...+0  U  1 -  + ... + k U ] ir  r  +  kJ-^-fi 2  k%\U\ + . . . + kafi + ... + k{ u r  ka H~ 1 kj'  r  [knu-i + ... + (ku -  (B.5)  + ••• + k u] ir  r  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  + • • • + (ka ~ l ) / i + • • • + ki u r  r  — \kn  • • • (ka ~ 1) • • •  h  (B-6)  And the coefficients of these bracketed terms become: hi  I  ka  kii  kf kj-  1  T  hi  ka + 1  (B.7)  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items