Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Deformable prostate registration from MR and TRUS images using surface error driven FEM models Taquee, Farheen 2011-12-31

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
24-ubc_2011_fall_taquee_farheen.pdf [ 17.5MB ]
Metadata
JSON: 24-1.0072359.json
JSON-LD: 24-1.0072359-ld.json
RDF/XML (Pretty): 24-1.0072359-rdf.xml
RDF/JSON: 24-1.0072359-rdf.json
Turtle: 24-1.0072359-turtle.txt
N-Triples: 24-1.0072359-rdf-ntriples.txt
Original Record: 24-1.0072359-source.json
Full Text
24-1.0072359-fulltext.txt
Citation
24-1.0072359.ris

Full Text

Deformable Prostate Registration from MR and TRUS images Using Surface Error Driven FEM Models  by Farheen Taquee B.Sc, Sir Syed University of Engineering and Technology, 2002  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  Master of Applied Science in THE FACULTY OF GRADUATE STUDIES (Electrical and Computer Engineering)  The University Of British Columbia (Vancouver) October 2011 c Farheen Taquee, 2011  Abstract TransRectal Ultrasound (TRUS) is used for image guidance during prostate biopsy and for treatment planning of brachytherapy due to low cost and accessibility in operating room. However, tumors have better visibility in Magnetic Resonance (MR) images. The fusion of TRUS and MR images of the prostate can aid with the diagnosis and treatment planning for prostate cancer and with post-brachytheraphy quality assurance. We developed a 3D deformable registration method using the segmentations obtained from TRUS and MR images and a biomechanical model that employs stiffness values derived from elastography. The segmented source volume is meshed and a linear finite element model is created for it. This volume is deformed to the target image volume by applying surface forces computed by assuming a negative relative pressure between the non-overlapping and the overlapping regions of the volumes. This pressure drives the model to increase the volume overlap until the surfaces are aligned. We tested our algorithm on prostate surfaces extracted from postoperative MR and TRUS images for 14 patients and pre-operative MR and TRUS images for 4 patients, using a model with elasticity within the range reported in the literature for the prostate. We used three evaluation metrics for validation: the Dice Similarity Coefficient (DSC) (ideally equal to 1.0), the volume change in source surface during registration, and the Target Registration Error (TRE) defined as the mean distance between landmarks such as urethrae and calcifications. For post-operative images, we obtained a DSC of 0.96±0.02 and a TRE of 1.5±1.4mm. The change in the volume of the source surface was 1.5±1.4%. For pre-operative images, we obtained the DSC of 0.96±0.01 and a TRE of 1.3±0.8mm. The change in the volume of the source surface was -0.9±0.2%. Our results show that this ii  method is a promising tool for physically-based deformable surface registration. We also used our technique to register ultrasound strain images to free mount histo-pathology images with the goal of correlating cancer with areas of low strain. This was done using relative stiffness values derived from vibroelastography data. We also performed Computed Tomography (CT) and Ultrasound (US) kidney surface registration using this technique.  iii  Preface The idea behind the work presented in this thesis was developed by Septimiu Salcudean and Orcun Goksel. It was brought to its present stage by the author under the supervision of Septimiu Salcudean. The registration technique presented in this thesis uses Finite Element Modeling (FEM), finding intersection and union of surfaces and non-overlapping volume. The codes for these steps were written by Orcun Goksel. The principal axis alignment was done by using part of a code written by S. Sara Mahdavi. The rest of the programming was done by the author. Septimiu Salcudean wrote parts of abstract and parts of chapters 2 and 3. He also reviewed and suggested corrections on the thesis manuscript. The work presented in this thesis has been accepted for presentation in SPIE (International Society of Optics and Photonics) Medical Imaging Conference 2012 - Image-Guided Procedures, Robotic Interventions, and Modeling and will appear in Proceedings of SPIE Volume 8316. Orcun Goksel and Sara Mahdavi also reviewed and suggested corrections to the said paper. The registration technique presented in this thesis has been evaluated using four sets of data. All the data sets were obtained as part of studies approved by Research Ethics Board (REB). The pre-brachytherapy prostate data was collected as part of a study for visualization and real-time dosimetry for prostate brachytherapy. The UBC CREB number of this study is H06-70146. The principal investigator for this study is Septimiu Salcudean. The images were obtained by Xu Wen and S. Sara Mahdavi and Mehdi Moradi at Vancouver Cancer Center. The post-operative prostate data was collected as part of a study done at British Columbia Cancer Agency for comparison of post-implant US/CT and MRI/CT  iv  image fusion in determining I125 prostate brachytherapy post implant dosimetry 1 . The REB number for this study was H03-60010. The primary investigator for this study was Mira Keyes. The prostatectomy data was collected as a part of study for Optimization of Elastography Imaging of the Prostate. The UBC CREB number of this study is H08-02696. The principal investigator for this study is Septimiu Salcudean. The images were collected by S. Sara Mahdavi, Mehdi Moradi and Guy Nir at Vancouver General Hospital. The kidney data was collected as part of a study for Real-time Image Guidance for Robot-Assisted Laparoscopic Partial Nephrectomy. The UBC CREB number of this study is H08-02798. The principal investigator for this study is Christopher Nguan. The images were collected by Caitlin Schneider, John Bartlett, Raoul Kingma and Michael Yip.  1 Patanjali, N.;  Keyes, M.; Morris, W.; Liu, M.; Harrison, R.; Spadinger, I. & Moravan, V. A comparison of post-implant US/CT image fusion and MRI/CT image fusion for 125I prostate brachytherapy post implant dosimetry Clinics in Occupational and Environmental Medicine, Elsevier, 2009, 8, 124-124  v  Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  x  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xi  Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xiii  Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xv  1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  1.1  Image Registration . . . . . . . . . . . . . . . . . . . . . . . . .  2  1.2  Prostate Cancer . . . . . . . . . . . . . . . . . . . . . . . . . . .  3  1.2.1  Prostate Anatomy . . . . . . . . . . . . . . . . . . . . . .  4  1.2.2  Diagnosis . . . . . . . . . . . . . . . . . . . . . . . . . .  5  1.2.3  Treatment Options . . . . . . . . . . . . . . . . . . . . .  6  1.2.4  Need of Registration . . . . . . . . . . . . . . . . . . . .  9  Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . .  10  1.3.1  Registration Model . . . . . . . . . . . . . . . . . . . . .  11  1.3.2  Boundary Conditions . . . . . . . . . . . . . . . . . . . .  12  1.4  Proposed Method . . . . . . . . . . . . . . . . . . . . . . . . . .  13  1.5  Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . .  14  1.3  vi  2 Deformable Registration Model . . . . . . . . . . . . . . . . . . . .  15  2.1  Deformation and Elasticity . . . . . . . . . . . . . . . . . . . . .  15  2.2  The Linear Elastic Model . . . . . . . . . . . . . . . . . . . . . .  16  2.2.1  The Finite Element Method . . . . . . . . . . . . . . . .  18  Factors Affecting the Stiffness Matrix . . . . . . . . . . . . . . .  21  2.3.1  Inhomogeneity of Material . . . . . . . . . . . . . . . . .  21  2.3.2  Stiffness from Elastography . . . . . . . . . . . . . . . .  22  Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . .  22  2.4.1  Overlapping Volume . . . . . . . . . . . . . . . . . . . .  23  2.4.2  Intersection and Union of the Surfaces . . . . . . . . . . .  24  2.4.3  Calculation of Surface Force . . . . . . . . . . . . . . . .  24  Computational Model . . . . . . . . . . . . . . . . . . . . . . . .  28  2.5.1  Chattering of Forces . . . . . . . . . . . . . . . . . . . .  29  2.5.2  Convergence Criteria . . . . . . . . . . . . . . . . . . . .  30  2.5.3  Time Step Reduction . . . . . . . . . . . . . . . . . . . .  30  Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  30  3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  32  2.3  2.4  2.5  2.6  3.1  Surface Generation . . . . . . . . . . . . . . . . . . . . . . . . .  32  3.2  Mesh Generation . . . . . . . . . . . . . . . . . . . . . . . . . .  33  3.3  Stiffness Parameters and Amount of Pressure . . . . . . . . . . .  34  3.3.1  Stiffness from Transfer Function Images . . . . . . . . . .  35  3.3.2  Amount of Pressure . . . . . . . . . . . . . . . . . . . . .  37  3.4  Registration Algorithm . . . . . . . . . . . . . . . . . . . . . . .  38  3.5  Finding the Position of Landmarks in the Deformed Mesh . . . .  40  3.6  Displaying the Result in 2-D Images . . . . . . . . . . . . . . . .  40  3.7  Validation Measures . . . . . . . . . . . . . . . . . . . . . . . . .  41  3.8  Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  42  3.8.1  Error in Marking Boundary and Landmarks . . . . . . . .  42  3.8.2  Error Due to Difference in Slice Thickness . . . . . . . .  42  3.8.3  Difference in Volume . . . . . . . . . . . . . . . . . . . .  43  3.8.4  Discretization Error . . . . . . . . . . . . . . . . . . . . .  43  3.8.5  Computational Performance . . . . . . . . . . . . . . . .  43  vii  4 Evaluation and Application . . . . . . . . . . . . . . . . . . . . . . . 4.1  Advantage of Using Non - Homogeneous Stiffness . . . . . . . .  4.2  Prostate Magnetic Resonance (MR) and TransRectal UltraSound  44 44  (TRUS) Surface Registration . . . . . . . . . . . . . . . . . . . .  45  4.2.1  Patient Data . . . . . . . . . . . . . . . . . . . . . . . . .  47  4.2.2  Landmarks . . . . . . . . . . . . . . . . . . . . . . . . .  47  4.2.3  Results . . . . . . . . . . . . . . . . . . . . . . . . . . .  48  4.2.4  Discussion . . . . . . . . . . . . . . . . . . . . . . . . .  50  Pathology and Vibroelastography Registration . . . . . . . . . . .  51  4.3.1  Patient Data . . . . . . . . . . . . . . . . . . . . . . . . .  53  4.3.2  Landmarks . . . . . . . . . . . . . . . . . . . . . . . . .  53  4.3.3  Results . . . . . . . . . . . . . . . . . . . . . . . . . . .  53  4.3.4  Discussion . . . . . . . . . . . . . . . . . . . . . . . . .  54  Kidney CT and US Registration . . . . . . . . . . . . . . . . . .  54  4.4.1  Patient Data . . . . . . . . . . . . . . . . . . . . . . . . .  54  4.4.2  Landmarks . . . . . . . . . . . . . . . . . . . . . . . . .  55  4.4.3  Results . . . . . . . . . . . . . . . . . . . . . . . . . . .  55  5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  57  4.3  4.4  5.1  Summary of Contributions . . . . . . . . . . . . . . . . . . . . .  57  5.2  Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . . .  59  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  61  Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  70  A Stress and Strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  70  A.1 Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  70  A.2 Strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  71  B Prostate Registration Results . . . . . . . . . . . . . . . . . . . . . .  74  B.1 Results of Pre-Operative Prostate Surfaces . . . . . . . . . . . . .  74  B.2 Results of Post-Operative Prostate Surfaces . . . . . . . . . . . .  79  C Pathology and VE Registration Results . . . . . . . . . . . . . . . .  94  viii  D Kidney Registration Results . . . . . . . . . . . . . . . . . . . . . . 100 E Amount of Pressure and TRE . . . . . . . . . . . . . . . . . . . . . . 104  ix  List of Tables Table 4.1  Results for deformation recovery test using homogeneous and non-homogeneous stiffness . . . . . . . . . . . . . . . . . . .  46  Table 4.2  Results for the four pre-brachytherapy images. . . . . . . . . .  48  Table 4.3  Results for the fourteen post brachytherapy images. . . . . . .  49  Table 4.4  Results for the fourteen post brachytherapy images using distance gradients based boundary conditions. . . . . . . . . . . .  52  Table 4.5  Results of kidney surface registration . . . . . . . . . . . . . .  55  Table E.1  Results for first five cases with two levels of pressure. . . . . . 105  x  List of Figures Figure 1.1  Prostate zones . . . . . . . . . . . . . . . . . . . . . . . . . .  4  Figure 2.1  Illustration of discretization . . . . . . . . . . . . . . . . . .  19  Figure 2.2  Intersection and union surfaces . . . . . . . . . . . . . . . . .  24  Figure 2.3  Overlapping and non-overlapping volume . . . . . . . . . . .  25  Figure 2.4  Illustration of connected surface triangles . . . . . . . . . . .  26  Figure 2.5  Force distribution in a mesh . . . . . . . . . . . . . . . . . .  27  Figure 3.1  Block diagram showing deformable registration algorithm . .  32  Figure 3.2  Prostate imaging using TRUS probe . . . . . . . . . . . . . .  34  Figure 3.3  A screen shot of Stradwin . . . . . . . . . . . . . . . . . . .  35  Figure 3.4  Triangulated surface in Matlab . . . . . . . . . . . . . . . . .  36  Figure 3.5  Exterior and interior of a mesh . . . . . . . . . . . . . . . . .  37  Figure 3.6  Young’s modulus from transfer function images . . . . . . . .  38  Figure 4.1  Deformably registered post-op MR and TRUS surfaces . . . .  50  Figure 4.2  Comparison of a result with surface distance based registration  51  Figure 4.3  2-D images showing VE and pathology registration . . . . . .  54  Figure 4.4  Deformably registered CT and US kidney surfaces . . . . . .  56  Figure A.1  Illustration of deformation and displacement between an at rest and deformed body . . . . . . . . . . . . . . . . . . . . . . .  71  Figure B.1  Registration result for pre-op case 1. . . . . . . . . . . . . . .  75  Figure B.2  Registration result for pre-op case 2. . . . . . . . . . . . . . .  76  Figure B.3  Registration result for pre-op case 3. . . . . . . . . . . . . . .  77  xi  Figure B.4  Registration result for pre-op case 4. . . . . . . . . . . . . . .  78  Figure B.5  Registration result for post-op case 1. . . . . . . . . . . . . .  80  Figure B.6  Registration result for post-op case 2. . . . . . . . . . . . . .  81  Figure B.7  Registration result for post-op case 3. . . . . . . . . . . . . .  82  Figure B.8  Registration result for post-op case 4. . . . . . . . . . . . . .  83  Figure B.9  Registration result for post-op case 5. . . . . . . . . . . . . .  84  Figure B.10 Registration result for post-op case 6. . . . . . . . . . . . . .  85  Figure B.11 Registration result for post-op case 7. . . . . . . . . . . . . .  86  Figure B.12 Registration result for post-op case 8. . . . . . . . . . . . . .  87  Figure B.13 Registration result for post-op case 9. . . . . . . . . . . . . .  88  Figure B.14 Registration result for post-op case 10. . . . . . . . . . . . . .  89  Figure B.15 Registration result for post-op case 11. . . . . . . . . . . . . .  90  Figure B.16 Registration result for post-op case 12. . . . . . . . . . . . . .  91  Figure B.17 Registration result for post-op case 13. . . . . . . . . . . . . .  92  Figure B.18 Registration result for post-op case 14. . . . . . . . . . . . . .  93  Figure C.1  Registered transfer function images for case 1 . . . . . . . . .  95  Figure C.2  Registered transfer function images for case 2 . . . . . . . . .  96  Figure C.3  Registered transfer function images for case 3 . . . . . . . . .  97  Figure C.4  Registered transfer function images for case 4 . . . . . . . . .  98  Figure C.5  Registered transfer function images for case 5 . . . . . . . . .  99  Figure D.1  Registration result for case 1. . . . . . . . . . . . . . . . . . . 101  Figure D.2  Registration result for case 2. . . . . . . . . . . . . . . . . . . 101  Figure D.3  Registration result for case 3. . . . . . . . . . . . . . . . . . . 102  Figure D.4  Registration result for case 4. . . . . . . . . . . . . . . . . . . 102  Figure D.5  Registration result for case 5. . . . . . . . . . . . . . . . . . . 103  xii  Glossary 2- D  Two-Dimensional  3- D  Three-Dimensional  3 D - CRT 3-Dimensional Conformal Radiotherapy CGAL  Computational Geometry Algorithms Library  CT  Computed Tomography  CTV  Clinical Target Volume  DICOM  Digital Imaging and Communications in Medicine, a standard for handling, storing, printing, and transmitting information in medical imaging.  DSC  Dice Similarity Coefficient  EBRT  External Beam Radiation Therapy  FEM  Finite Element Method  GTV  Gross Tumor Volume  HDR  High Dose Radiation  IGRT  Image-Guided Radiotherapy  IMRT  Intensity-Modulated Radiotherapy  LDR  Low Dose Radiation xiii  LRP  Laparoscopic Radical Prostatectomy  MR  Magnetic Resonance  MRS  Magnetic Resonance Spectroscopy A set of polyhedra obtained from Boolean operations of set intersection and set complement of half spaces  NOV  Non-Overlapping Volume  OV  Overlapping Volume  PAA  Principal Axis Alignment  PTV  Planning Target Volume  SNC  Selective Nef Complex  SI  Superior-Inferior  TF  Transfer Function  TRE  Target Registration Error  TRUS  TransRectal UltraSound  US  Ultrasound  VC  Volume Change  VD  Volume Difference  VE  Vibro-elastography  VRML  Virtual Reality Modeling Language  xiv  Acknowledgments This research project would not have been possible without the support of many people. First and foremost I want to express my deep and sincere gratitude towards my supervisor Dr. Tim Salcudean for his support, patience and immense knowledge. His guidance and enthusiasm kept me motivated through out the research. I would not have imagined having a better supervisor. I want to take this opportunity to thank Orcun Goksel and Sara Mahdavi for all the discussions, ideas and help. Thanks for being so patient with me. My work is based on finite element model and volume overlap calculation codes written and provided by Orcun. Thanks Sara for sharing little pieces of codes and literature that helped me in my. I am lucky to have a lab mate and friend like you. My research would not have been completed with out the support of my husband Husain Taquee. He kept me going. I am grateful to staff at Liliput for taking care of my daughter so that I could study and do research. I would like to thank my parents and siblings for their encouragement, advice and positive words. The financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) is gratefully acknowledged.  xv  Chapter 1  Introduction Prostate cancer is among the most prevalent non-cutaneous cancer in men: in 2011 it is estimated that 25,500 men will be diagnosed with and 4,100 will die from it in Canada. Due to better treatment options and early diagnosis, survival for prostate cancer is consistently high (>95%) among men aged 40 to 79 for the past five years [1]. Prostate cancer can be treated using minimally non-invasive surgical procedures which reduce hospital stay and patient trauma. TransRectal UltraSound (TRUS) has become the standard image guidance tool for most of the prostate related minimally invasive surgical procedures and treatment planning due to its low cost, safety and availability in operating room. It is a method of imaging pelvic organs using an ultrasound probe inserted in the rectum. However, due to low contrast in Ultrasound (US) B-mode images, tumors are generally not visible [57]. Magnetic Resonance (MR) images, although expensive, provide better visibility of the tumors [4].  MR  images also show patient specific  pelvic anatomy and give information on the spatial spread of cancer which can be useful in keeping trauma to minimum and avoiding injuries.  MR  can be brought in  the operating room virtually by fusing MR acquired before the procedure to the US images acquired during it [49, 53]. This fusion, a product of image registration, can be useful for biopsy guidance, treatment planning for brachytherapy and focal therapy. Fusion of Computed Tomography (CT) or  MR  with  TRUS  which shows the seed distribution at the end of  the brachytherapy procedure can also be used for quality assurance after brachy1  therapy.  1.1 Image Registration “Image registration is the process for determining the correspondence of features between images collected at different times or using different imaging modalities. Once this correspondence is known, it can be used to change the spatial mapping of one image so it more closely resembles another so the pair can be directly compared, combined or analyzed” [16]. Medical image registration is commonly used to correct for different patient positions between the scans [42]. It is also used to watch progress of a disease over time in a patient [24, 37] and to create anatomical atlases [35]. It basically allows viewing and analyzing of two different types of images in one coordinate frame. The correspondence between two images or image volumes is established by finding the spatial transformation between the pixel positions that correspond to the same structural or functional areas in two images. Sometimes, this transformation can be achieved by translating and rotating one image so that it matches the other. This is called rigid registration. This type of registration is limited to six degrees of freedom for Three-Dimensional (3- D) to 3- D images and may not be enough to accurately match the two images. In cases where it is not enough, a more complex transformation involving stretching and shrinking of areas of images is required. Again, this transformation tells how one image can be changed so that it matches the other. This is called non-rigid or deformable registration. Images of soft tissues of human body usually require deformable registration. The image that is deformed or registered is called is the Source and the image which it tries to match is called the Target. In addition to a transformation model, a registration algorithm uses a similarity measure to determine how well two images match. An optimization process is employed to find the maximum similarity by varying different parameters of the transformation model. Medical image registration can be divided into two broad categories: model based and intensity based. The former builds explicit models of surfaces, curves and point landmarks corresponding to anatomical structures in one image that can  2  be matched with their counterparts in the other image. These correspondences are used to find the transformation from one image to the other. As model based registration matches the structural information, the registration has better anatomical correctness and can be interpreted in terms of underlying anatomy [16]. The models most commonly used for this type of registration are splines (including thin-plate splines and B-splines), elastic models [6, 9], viscous fluid model [14], demons [10, 62] and active shape models [15]. Intensity based registration algorithms assume a statistical relationship between intensities of two images and define an intensity based measure of similarity. The optimization process adjusts the transformation until the similarity measure is maximized [16]. The similarity measures commonly used are the squared differences in intensities, the correlation coefficient, measures based on optical flow, and mutual information. Intensity based registration matches intensity patterns over the whole image but do not use anatomical knowledge. More details of these measures can be found [29]. Hybrid models combine intensity based and model based techniques to create better correspondence [18, 50]. We have chosen a model based approach for registration of prostate TRUS  MR  and  images because due to the poor contrast and presence of noise in B-mode  images it is very hard to establish a statistical relationship between the intensities in two modalities. Another reason for this choice is that for every prostate cancer intervention prostate surfaces are segmented from TRUS images. This can be used in model based registration.  1.2 Prostate Cancer Prostate cancer is the malignant growth of cells in the prostate gland. The cancerous cells divide at a harmful rate resulting in a large number of cells with possibly abnormal numbers of chromosomes. Cancerous cells do not die and get replaced like healthy cells do and therefore grow into malignant tumors, which can spread through the tissue of the prostate gland and metastasize into surrounding tissues [3].  3  Figure 1.1: A diagram showing the different zones of prostate  1.2.1  Prostate Anatomy  The prostate is a gland in the male reproductive system. Its function is to store and secrete a slightly alkaline fluid which usually constitutes 20-30% of the volume of the semen. In healthy adult males its size is slightly larger than a walnut. The weight of a healthy prostate in adult males ranges from 7 to 16 gms with an average weight about 11 grams [34]. The prostate sits above the base of the penis below the urinary bladder and backs onto the front wall of the rectum. The apex of the prostate is pointed down to the perineum as opposed to the base which is wider and located next to the bladder. The prostatic urethra is the portion of urethra that runs from the urinary bladder through the prostate and exits from the apex via the urinary sphincter which is a group of muscles that prevents involuntary leakage of urine. The prostate is surrounded by a membrane called the prostatic capsule. In pathology, the regions of the prostate are classified as zones Figure 1.1. The prostate gland has four distinct glandular regions: 1. Peripheral zone (PZ): This zone occupies approximately 70% of the volume of gland. 70-80% of prostatic cancers originate in peripheral zone [3]. 2. Central zone (CZ): This zone constitutes approximately 25% of the prostate 4  gland. About 5% of prostate cancer cases originate in the central zone [3]. 3. Transition zone (TZ): The transition zone is the innermost part of the prostate gland and surrounds the urethra. It makes up about 5% of the prostate volume. About 10% of prostate cancers occur in this zone. This zone also enlarges with age and can result in benign prostatic enlargement [3]. 4. Anterior fibro-muscular zone (or stroma): The anterior zone is located close to the abdomen (away from the rectum). This zone is constitutes 5% of the gland volume and is composed mostly of muscular tissue. The elastography studies of the prostate have shown that varying stiffnesses occur in these zones. For example, [32] reports that the central zone is less stiff than the peripheral zone. Moreover, the cancerous tissues have higher stiffness than healthy tissue [33, 66]. Zhang et al. [66] has also reported lower than healthy stiffness for tissue with benign prostatic hyperplasia (BPH). This difference in stiffness affects the way the prostate deforms. This should be taken into account when modeling the prostate as an elastic object for simulation or registration.  1.2.2  Diagnosis  Prostate cancer is usually suspected when a high level of prostate specific antigen (PSA) is detected in blood tests. A digital rectal examination (DRE), in which the physician palpates the prostate through the rectum is then performed to detect any abnormalities within the prostate. A biopsy is performed to confirm the presence of cancer. Prostate Biopsy A prostate gland biopsy is a diagnostic procedure which involves removal and examination of small samples of tissue. To remove tissue samples a needle is inserted either through the rectum (transrectal biopsy), through the urethra, or through the area between the anus and scrotum (perineum). A transrectal biopsy is the most common method used.  TRUS  is commonly used to guide the placement of the nee-  dle during a prostate biopsy. However, due to low contrast in B-mode images, it is impossible to localize cancer directly from the 5  TRUS  images. To detect cancer  the prostate is usually sampled uniformly using a sextant biopsy template, a fiveregion biopsy template, an eight-biopsy systematic template, or an the eleven-core biopsy multi-site directed template. The number of biopsy cores, even increased, is still limited and these cores may miss the cancer region, leading to inaccurate or false negative results [22, 45]. To avoid false negatives, localized biopsy, where samples are taken from suspected cancer regions, seems an obvious solution. This requires knowledge of the location and spatial spread of cancer from MR and Magnetic Resonance Spectroscopy (MRS)images [4, 8, 54]. This is explored more in Section 1.2.4.  1.2.3  Treatment Options  Modern health care gives many treatment options for prostate cancer. In cases of cancers confined to the prostate gland, radiation is often used to destroy cancer cells either through non invasive External Beam Radiation Therapy (EBRT) or minimally invasive brachytherapy.  EBRT  uses electron, proton, or neutron beams directed  at the prostate (see Section 1.2.3) while brachytherapy involves small radioactive seeds implantation inside the tissue using needles. In advanced stages of prostate cancer, chemotherapy is used to destroy cancer cells via drugs. Another minimally invasive therapy known as cryotherapy freezes up cancerous tissue using freezing gases and needles. Hormone therapy is used to cut down the growth of prostatic tissue and thus cancer. Another option is a surgical procedure that removes the entire prostate gland called the radical prostatectomy. In early stages of cancer, watchful waiting is often recommended. This involves carefully monitoring the development of tumor before starting advanced therapy. We examine some of the treatment options in more detail in the following section. External Beam Radiation Therapy The  EBRT  uses high energy radiation beams to irradiate target cancerous tissues  while sparing the neighboring healthy tissue. This includes conventional EBRT, 3Dimensional Conformal Radiotherapy (3 D - CRT), Intensity-Modulated Radiotherapy (IMRT), and Image-Guided Radiotherapy (IGRT). Localization of the target and  6  adjacent normal tissue is critical in the planning of all types of EBRT treatments. Conventional  EBRT  entails irradiation of the entire pelvis. Flouroscopy and  plain radiographs are used to plan treatment. The four field (anteroposterior, posteroanterior, left lateral, right lateral) technique which targets the prostate, seminal vesicles, and regional lymphatics is used. The dose is delivered over a period of 5 -6 weeks. An additional dose, called boost, is delivered to the target volume (prostate and seminal vesicles plus 1-2 cm margin) separately. Complications due to exposure to radiation can arise in skin, bladder, intestines, pelvic bone etc. 3 D - CRT delivers a precisely shaped field to the 3- D anatomical map of target volume (prostate, seminal vesicals etc.) obtained from  CT  or  MR .  The visibly  affected (cancerous) areas are delineated from images to form Gross Tumor Volume (GTV). This volume is expanded to include areas possibly affected but not visible in on the images to form a Clinical Target Volume (CTV). This is further expanded to include a Planning Target Volume (PTV) to include any position errors due to patient positioning and setup.  IMRT  can be considered an advanced form of  3 D - CRT. It delivers highly conformal doses by using non-uniform radiation beams directed at target volume from different directions. It should be noted again that the dose is delivered over a period of 5 to 6 weeks. Locking devices are used in 3 D - CRT and  IMRT  to reproduce the patient pose to  match the one at which images used for planning were obtained. These locking devices however do not guarantee that the target volume obtained at the beginning of treatment conforms to the orientation and shape of prostate at the start of or during each treatment session and can introduce systematic errors [61]. There is possibility of deformation in the prostate due to differences in bladder and rectal filling. Bylund et al. report significant average prostate motion of 6.7 mm for 24 patients [12]. Oncologists address this issue by using PTV which includes margins that allow for these differences. Adjustment to the initially obtained target volume is an intuitive solution.  IGRT  provides this solution by imaging the target area  during the treatment. Abdominal Ultrasound and x-ray cameras are being used to image the pelvic organs during treatment.  7  Brachytherapy Brachytherapy, also known as interstitial radiation therapy, seed therapy, or seed treatment is a minimally invasive radiation therapy. There are two types of brachytherapy: permanent Low Dose Radiation (LDR) and temporary High Dose Radiation (HDR). In  LDR  brachytherapy, brachytherapy seeds consisting of iodine-  125 and palladium-103 inside titanium pallets are implanted permanently in the prostate gland. These emit a low level of radiation over the course of their radioactive lives. In HDR brachytherapy, high dose from a single radioactive seed made of iridium-194 (also called the wire) is delivered through plastic catheters inserted in the gland. A patient undergoes two or three dose treatments during an overnight stay for HDR brachytherapy [3]. Before brachytherapy, TRUS is used to obtain transverse B-mode images from base to apex. The images are segmented and a prostate volume or surface is obtained. Margins are added to the prostate volume to form a  PTV .  This is done to  encompass any spread of prostate cancer out of the prostatic capsule. A radiation oncologist devises a plan for seed arrangement after calculating dose distribution. Radioactive seeds are implanted in the prostate according to the plan using  TRUS  for image guidance. It is important to cover the entire prostate according to the plan and spare the healthy tissue surrounding it. If this plan is not executed well, some parts of prostate may miss seeds and/or rectum and urethra can an additional extra dose. This can cause complications. Also due to lack of visibility, in some cases, injury to blood vessels and nerves and seed drift can give rise to complications. Radical Prostatectomy Radical prostatectomy is the surgical removal of the entire prostate along with prostatic urethra. The two detached ends of urethra are joined using in a connection called the anastomosis. The seminal vesicles and surrounding blood vessels and nerves may also be removed. This is oldest treatment for prostate cancer. The prostate surgery can be performed open i.e., through an incision in the lower abdomen, called Retropubic Radical Prostatectomy (RRP), or through the perineum, called Perineal Radical Prostatectomy (PRP), or through a laparoscope, called Laparoscopic Radical Prostatectomy (LRP). A laparoscope is a slender tube like in-  8  strument that is inserted inside the body through a small incision to show surgeons the view of the anatomy. Other small incisions for surgical tools are also made through which the surgeon performs the surgery. In Robot Assisted Laparoscopic Prostatectomy (RALP) surgeons use robotic arms to perform the surgery remotely. Laparoscopic surgeries result in less blood loss and shorter convalescence time. Sexual dysfunction is one of the most common side effects of prostatectomy. Even though surgeons try to spare the nerves to preserve sexual potency, the location of cancer and skill and experience of the surgeon determine the outcome.  1.2.4  Need of Registration  Even though the the location of tumors can be known from pre-operative  MR  im-  ages, due to differences in patient position, bladder fullness and the presence of the TRUS  probe, the prostate is deformed in intra-operative  TRUS  images and cancer  can not be correctly localized in B-mode images unless a registration is performed. As rigid registration cannot provide accurate registration [56], deformable registration is needed. For Biopsy To guide biopsy, the  MR  images acquired prior to the procedure can be registered  to the TRUS images obtained at the start of or during the procedure. This will give a 3- D map of cancers within the prostate volume. This can lead to the possibility of localized biopsy which eliminates false negatives and can reduce the number of samples taken from prostatic tissue [51]. For External Beam Radiation The  PTV  used to initially plan treatment can be registered to the  obtained during each treatment session of  IGRT.  US  or  CT  images  This will result in dose delivery  that is truly in agreement with the original plan and further decreases morbidity in the surrounding tissues and structures. Ling et al. recognizes that result in re-examination of current practice of expanding that to  PTV  GTV  to  IGRT  CTV  could  and from  [36]. Registration of during-the-treatment images to those used for  planning can help create target volumes with tighter margins which will in effect 9  reduce dose to normal tissues . For Brachytherapy In brachytherapy, after seeds are implanted, it is important to know how close the seed distribution is to the original plan. Deviation from the plan can affect the planned dose significantly [19].  CT  and  MR  images show seeds very well.  Registering prostate surfaces from post-brachytherapy images of TRUS  CT  or  MR  to the  surface on which the original plan was made, can show how well the actual  positions of seeds match the ones in the plan. This is essential to study dose-effect relations of prostate implants and long term quality assurance [55]. For Prostatectomy In order to avoid recurrence of cancer, the removal of organ and associated structures and tissue should been done in a manner as to include all the cancer affected areas. This is achieved by using wider surgical margins. Most of the cancers occur in the peripheral zone and wider margins in that zone may damage the neurovascular bundles and the urethral sphincter. Ukimuraa et al. showed that using TRUS  guidance during  LRP  can significantly reduce the incidence of positive sur-  gical margins (i.e., when part of cancerous tissue is left inside the body) [60]. For nerve-sparing prostatectomy, a pre-operative MR or TRUS which identifies the nervovascular bundle can be registered with intra-operative TRUS to effectively avoid damage to nerves. Watchful Waiting Registration of images is very useful to monitor disease progression. The images obtained periodically when registered show the changes in the different areas of the organ explicitly. This can help physicians assess the progress and plan the treatment accordingly.  1.3 Literature Review We have divided the literature review for this project in two categories: registration models and boundary conditions. The former discusses registration models used 10  in the literature for TRUS and MR prostate registration. The later gives an overview of techniques available for calculating boundary conditions for our chosen model: linear elastic finite element mesh.  1.3.1  Registration Model  There is extensive literature available on the general subject of medical image registration; however, prior work in the specific category of TRUS to MR prostate image registration is more limited. In the past rigid registration methods [63, 64] were used for quick registration allowing for real-time use. These methods, however, ignored the possibility of significant deformations due to the presence of the TRUS  probe, change in patient position, bladder fullness etc. Zaider et al. [65] use  a very simplified registration technique using the assumption that the position of points within the gland with respect to the axial contours remains unchanged, and no displacement occurs in the Superior-Inferior (SI) direction. For deformable registration a few spline based models were proposed. Shao et al. [51] used thin plate splines to integrate  MR  or  MRS  with  TRUS  for targeted  robotic prostate biopsy. The deformable registration carried out Two-Dimensional (2- D) to 2- D after the initial alignment and therefore significant deformation in the  SI  direction was ignored. More recently, Mitra et al. [40] also used thin plate  splines for deformable registration of 2- D  TRUS  and  MR  prostate images with au-  tomatic generation of control points from principal axes. This method relies on finding the corresponding slice in two modalities which can only be approximate due to difference in slice spacing. Again, deformation in the  SI  direction is ig-  nored. The paper reports a Dice Similarity Coefficient (DSC) to measure overlap of contours but does not report a Target Registration Error (TRE) validating internal registration. Reynier et al. [49] and Daanen et al. [17] use an octree-spline for deformable registration after performing rigid registration. Reynier et al. [49] reported TRE for four patients. Another popular technique for deformable registration is to use a biomechanical model. A biomechanical model represents organs of interest as elastic meshed objects. Bharatha et al. [7] used an finite element linear elastic model with different material properties for the central and peripheral zones of the prostate driven  11  by surface distance based forces. Alterovitz et al. [5] also use a 2- D linear elastic Finite Element Method (FEM) model to register MRS/MR images with and without endo-rectal coil for treatment planning of  HDR  Therapy. They argue that a 2- D  model is sufficient as out of plane deformations are less then the slice thickness. Their FEM model uses a non-linear optimizer to estimate tissue stiffness and external forces causing deformation. Hensel et al. [30] used a multi-organ linear elastic  FEM  model of the pelvic  anatomy with different stiffness values assigned to prostate, bladder, rectum etc to register MR with endo-rectal coil to one without for radiation therapy planning. Noe et al. [44] used non-linear elastic models to perform  MR  to  MR  prostate reg-  istration. They evaluated their model with prostate data by validating inverse consistency of the transformation. Hu et al. [31] perform deformable registration between MR and TRUS using a patient specific finite element statistical model. They test their algorithm on eight patient data sets and report a mean TRE of 2.4mm.  1.3.2  Boundary Conditions  As mentioned in Section 1.4 and explained in Chapter 2 we use a linear elastic FEM  mesh as our registration model in this thesis. Therefore it is important that we  examine the boundary conditions that are used to derive the FEM model. Boundary conditions refer to the external force required to deform one surface to match the other. In most part of the literature, these forces are derived from some metric representing the distance between the two surface [7, 13, 44]. The goal of the registration is then to minimize this metric and thus have perfect alignment between the two surfaces. Bharatha et al. [7] used active surfaces and distance fields. This method was presented by Ferrant et al. [23] for surface matching of brain meshes. Noe et al. [44] use a surface distance term as well as a surface normal term to calculate the boundary conditions allowing a surface to fill in narrow extrusions while registering. Brock et al. [11] find correspondences based of curvature and use that for driving their boundary conditions. Choi et al. [13] also use a surface distance term with an iterative mesh fitting. Zhang et al. [67] use contact impact analysis for boundary conditions. The surfaces that are or could be in contact are with the organ of interest are modeled using  12  the FEM. The contact is determined by determining the penetration of elements and applying appropriate pressure. Alterovitz et al. [5] find boundary conditions and tissue stiffness parameters for their 2- D FEM model by an optimization routine that maximizes  DSC  and minimizes strain energy. Hensel et al. [30] use surface node  projection from source to target surface guided by the curvature of the curve to which the node belongs.  1.4 Proposed Method We propose a registration technique for registering the TRUS prostate images with the  MR  prostate images. In both the  TRUS  and the  MR  images, the prostate can  be segmented easily and therefore surfaces are available. Taking advantage of this fact, we propose a surface based registration method to deformably register surfaces obtained from the segmented contours using a biomechanical model. The model is solid linear elastic body which deforms when subjected to external forces on the surface. We use global surface difference to find boundary conditions instead of more popular local surface distances (see Section 1.3.2) and therefore produce a more realistic distribution of forces driving the model to the target. Our registration algorithm assumes the presence of a negative pressure in non-overlapping regions of the surface. This pressure increases the overlap so that the surfaces match. We extend the goal of our deformable registration method to the following, so that the transformation is physically realistic: 1. We wish to deform the source surface to match the target surface 2. We wish to use parameters (Young’s Modulus, Poisson’s ratio and forces) that are close to the values reported in literature. 3. We wish to preserve the source surface volume i.e., keep it constant during registration We validate our approach with elasticity parameters typical of those reported in the literature. The results are evaluated in terms of the improvement achieved  13  in the matching the prostate volumes. In addition, the distances between the urethra segmentations in the two modalities are used as an independent measure of registration accuracy.  1.5 Thesis Organization We start with Chapter 2 explaining our deformable registration model. In addition to providing details on our linear elastic model, this chapter discusses parameters such as tissue stiffness and boundary conditions. The Chapter 3 describes the implementation details of our algorithm. We move on to evaluate our algorithm in Chapter 4 and include some applications outside  TRUS  and  MR  prostate surface  registration. Finally conclusions, achievements and future directions of our work are laid out in Chapter 5.  14  Chapter 2  Deformable Registration Model When performing the deformable registration of two surfaces, we fix one surface during the registration and let the other move in thousands of degrees of freedom. The former is called the target and the later is called the source. It is the source that deforms to match the target. The source is modeled as a linear elastic solid and subjected to boundary conditions that aim for a perfect fit between the two surfaces. This is a dynamic process and therefore uses time integration. We use an iterative algorithm that converges when the cost of registration is minimum.  2.1 Deformation and Elasticity In continuum mechanics, deformation is defined as the transformation of the positions of constituent particles of a body from one configuration to another. Deformation is usually caused by external loads, body forces (such as gravity or electromagnetic forces), or temperature changes within the body. These factors induce internal body forces called the stress. The deformation defined in terms of relative displacement of the constituent particles of a body is called the strain. These terms are defined in more detail in Appendix A. The deformations in a continuous body can be elastic or plastic. Elastic deformations are characterized by removal of strain (i.e., reverting back of body to its pre-stress configuration) when stress is removed. In cases where strain is permanent and stays even after the removal of stress, the deformation is called plastic.  15  Plasticity is exhibited when the stress exceeds the elastic limit. Most materials, including biological tissues, are normally exposed to stresses much lower than the elastic limit and thus exhibit elasticity. Hyper-elasticity, a model for ideally elastic materials, can be used to describe the complex deformations of biological tissues. However, it is often computationally efficient to assume linear elasticity for modeling and simulating biological tissues. As all hyper-elastic materials can be approximated by linear elasticity when small deformations are applied, this assumption is quite accurate [25, 39]. For these reasons, we have used linear elasticity to model the prostate for now. In future more complex models like Neo-Hookean incompressible model can be used. Linear elasticity is a simplified version of non-linear elasticity that assumes infinitesimal strains or “small” deformations and linearly relates stress and strain.  2.2 The Linear Elastic Model If a 3- D linear elastic solid body Ω composed of particles p(x, y, z), initially under no stress, is subjected to a force F(x, y, z) and, as a result, deforms so that each of its particles undergo displacement u(x, y, z), the deformation energy is given by [9]: Estrain =  1 2  Ω  ε T σ dΩ −  Ω  F u dΩ  (2.1)  where ε and σ are the strain and stress vectors respectively. The first term presents strain energy and second term presents work done by external forces. As described in Appendix A, for small deformations, the strain vector can be written as:  δu δu δu δu δu δu δu δu δu , , , + , + , + ε= δx δy δz δy δx δz δx δz δy where B is strain-displacement matrix:  16  T  = Bu  (2.2)    δ δx  0 δ δy   0  0  B=δ δy  δ  δz 0  0 δ δx  0 δ δz  0     0  δ  δz   0  δ  δx  (2.3)  δ δy  The constitutive equation that describes the relationship between stress and strain for linear elastics is based on Hook’s law:  σ = Cε  (2.4)  where C is material stiffness matrix that characterizes Hookean elastic material. It is a symmetric 6x6 matrix for 3- D objects. For a homogenous and isotropic material this matrix is defined by the two Lam´e material constants λ and µ [9]:  λ + 2µ λ λ   λ λ + 2µ λ   λ λ λ + 2µ  C=  0 0 0   0 0 0  0  0  0    0  0  0  0  0  0  0  µ 0 0  0   0  0   0  0  µ  µ 0  (2.5)  The Lam´e’s constants can be derived from Young’s Modulus (E) and Poisson’s ratio (ν ) using these relationships:  λ=  Eν (1 + ν )(1 − 2ν )  (2.6)  E 2(1 + ν )  (2.7)  µ=  Now Equation 2.1 can be written as: EStrain =  1 2  Ω  uT BT CBu dΩ −  17  Ω  F u dΩ  (2.8)  This continuous function of energy can be discretized over displacement using the finite element model.  2.2.1  The Finite Element Method  As C is positive definite, Equation 2.4 is an elliptic partial differential equation which can be solved for simple shapes. The complex shapes can be dealt with by considering them as combinations of some simple shapes and then solving the boundary-value problem for each simple shape separately. This is the idea behind the  FEM .  These simple shapes are called the elements. Common choices for ele-  ments are squares, triangles and hexagons for 2- D surfaces and cubes, tetrahedra, and hexahedra for 3- D objects. In our finite element model we use 4-node tetrahedra connected with each other in a mesh like configuration to model a 3- D object. A tetrahedron is a four faced polyhedron. The vertices of the tetrahedron are called nodes. The surface of the object is formed by the outward facing triangular facets of the tetrahedra near the surface. The nodes that belong to the triangular faces on the surface are called the surface nodes. The nodes that exist inside the surface are called the interior nodes. Please see Figure 2.1 for an illustration of tetrahedra used as elements to model an ellipsoid. The  FEM  discretizes the displacement field resulting from deformation in the  3- D body. If a 3- D body Ω is divided in elements Ωe with a total number of n nodes and the displacement at a node i, due to deformation, is ui = [ui vi wi ]T , the compound vector of all the nodal displacements is given by u = [uT1 uT2 . . . uTn ]T . This vector represents the deformation of the entire 3- D object. Our chosen element for FEM, the tetrahedron, had four nodes. If the four nodes of an element have displacements uem where m = 1 . . . 4, the compound vector of element node displacements is given by ue = [(ue1 )T (ue2 )T (ue3 )T (ue4 )T ]T . The displacement of a point x lying somewhere in an element Ωe can be determined from the interpolation of element node displacements.  18  Figure 2.1: A 3-D object made up of small tetrahedral elements. One of the tetrahedral elements is magnified and shown with its nodes.  4  ∑ Nie (x)uei  u(x) =  (2.9)  i=1 e e  = Nu  (2.10)  The Nie (x) is element basis function. For linear interpolation, it consists of barycentric coordinates of the element nodes: Nie (x) = aei + bei x + cei y + die z  i = 1...4  (2.11)  Revisiting the problem of seeking equilibrium at which total energy (2.8) is minimized, consider when first variation of total energy is zero: e δ EStrain (u) = ∑ δ EStrain (u) = 0  (2.12)  e  The differential energy for an element is: e e δ EStrain (u) = ∇EStrain · δ ue  The Equation 2.12 can be written as: 19  (2.13)  δ E(u) = ∑(∇E e · δ ue ) = 0  (2.14)  e  The term on the side of Equation 2.14 can be zero for arbitrary displacement from equilibrium only when ∇E e is zero. This is because ue contains independent variables. This implies that  Ωe  (BN e )T C(BN e )ue dΩ − Fe = 0  ⇒  Ωe  BeT CBe ue dΩ = Fe  BeT CBe ue  ⇒  Ωe  dΩ = Fe  ⇒  V e BeT CBe ue = Fe  ⇒  Ke ue = Fe  (2.15)  where V e is the volume of element. Ke , the product of V e , BeT , C and Be , is the element stiffness matrix. Fe is the compound vector of the applied element nodal forces. The volume of a tetrahedron can be calculated using: x1 y1 z1 1 Ve =  1 x2 y2 z2 1 6 x3 y3 z3 1 x4 y4 z4 1  The element stiffness matrices for all the elements of the body can be combined to assemble a global stiffness matrix K for the body. This is done so that every constituting point of the elastic body receives its contribution from every Ke that its part of. Now we have a linear equation that describes the whole system: Ku = F  (2.16)  where size of K is 3n × 3n. The external force F and nodal displacements u are represented by a vectors of size 3n × 1. It should be noted here that for a certain 20  amount of stress F, the resulting deformation (displacement of nodes- u) depends on the stiffness matrix K which characterizes the stiffness of the whole object. Dynamics can be included in the model by using classical Newtonian formulation: Mu¨ + Bu˙ + Ku = F  (2.17)  where M and B are mass and damping matrices respectively.  2.3 Factors Affecting the Stiffness Matrix The stiffness of an object, in part, determines the degree of deformation it will undergo when a certain amount of force is applied to it. This depends on the elasticity of its constituent material as well as its geometric shape. Therefore, elasticity moduli can be used to measure relative stiffness for objects of approximately the same geometric shape. Another factor that affects the amount of deformation taking place irrespective of the shape of the object is the compressibility of the constituent material. Compressibility determines the volume change in an object as a result of application of stress. A perfectly incompressible object retains its volume in the state of deformation. We saw in Section 2.2, Young’s Modulus and Poisson’s ratio are used to calculate material stiffness matrix C and thus elemental stiffness matrix Ke . A perfectly incompressible object has a Poisson’s ratio of 0.5.  2.3.1  Inhomogeneity of Material  The elasticity properties of the constituent material of tissue is rarely homogeneous. Elasticity mapping from elastography show a wide range [32, 33, 46] within a organ, depending on the tissue type (muscular vs. glandular [46]) and the presence of a tumor or calcification [33, 48]. If the elasticity map of an organ is present, the stiffness matrix for each element of the mesh can be formed by using the average material properties of the volume of tissue occupied by that element. The global stiffness matrix can then be constructed in the manner described in Section 2.2. In case the elasticity map is not available, then homogeneous stiffness can be assumed for simplicity.  21  2.3.2  Stiffness from Elastography  As mentioned above, if an elasticity map of an organ is available, the stiffness matrix can be constructed using this map. Elastography is the technique used for determining the stiffness of a tissue from strain images obtained by application of some form of stress. The stress is usually mechanical compression or vibration with the application of which the softer regions of the tissue deform more and stiffer regions deform less. The strain images show the extent of deformation at each point (limited by the resolution of device used) and can be used to determine relative elasticity at each pixel in a strain image. Using realistic relative elasticity values is important for registration accuracy especially for areas inside the mesh because for the same amount of force, two regions of different stiffness will deform differently. This can be explained by the following example. Consider a piece of soft tissue with a stiff tumor in it. If pressure is exerted on this tissue, the soft tissue will deform and probably displace the tumor from its position but might not change the shape of the tumor. If this tissue is modeled as an elastic mesh using a homogeneous stiffness and the same amount of pressure is applied to the mesh, the result will be a deformed tissue with a deformed tumor. A more accurate displacement map can be obtained by using relative elasticity values derived from elastography (please see Section 4.1). Once the elastic model for the source surface has been created, the next task is to determine boundary conditions that explain the deformation from the source to the target surface.  2.4 Boundary Conditions By performing registration of two surfaces using the linear elastic FEM model, we try to find boundary conditions for which one object deforms to match the other. These boundary conditions are surface forces F in Equation 2.16 that produce deformation indicated by the nodal displacement vector u. Most of the work done in the literature uses some sort of distance measure between the surfaces for driving surface forces (see Section 1.3.2). The goal of these methods is to minimize distance between the surfaces. We have proposed a method of registration that takes into account the overall surface misalignment and maximizes the overlapping vol22  ume, or in other words, minimizes the non overlapping volume. While at a glance our method of force calculation might not look significantly different from the distance based methods, there is a big difference. Consider the state of two surfaces that exist so that some parts of them are close to each other and some parts are quite far. In the case of distance based methods, the parts of surface with small misalignment will experience small surface force while the parts with large misalignment will experience large forces. In our method, all areas of the surfaces get similar forces. This can result in the development of a more significant torsional effect. As the method uses overlapping volume to drive forces, the first step in our calculation is finding the overlapping volume which is also used for calculating our registration error.  2.4.1  Overlapping Volume  The term Overlapping Volume (OV) as used from now onwards refers to the volume of the region overlapped by two surfaces. It is calculated by forming an intersection surface from the two surfaces (see Section 2.4.2). Once the intersection surface is formed, the Non-Overlapping Volume (NOV) is calculated by subtracting OV  from the union of the source surface volume (VolSource ) and the target surface  volume (VolTarget ). This OV is our similarity measure (see Section 1.1) that we try to maximize.  NOV  = (VolSource − OV ) + (VolTarget − OV ) = VolSource +VolTarget − 2 · OV  (2.18)  The volume of source, target and intersection surfaces are calculated using the divergence theorem. The NOV can also be represented as a percentage: NOV (%) =  NOV × 100 VolSource +VolTarget  (2.19)  We use Equation 2.19 as cost of registration and minimize it through an iterative process. To find the OV it is essential to find the intersection of the two surfaces.  23  (a) Two surfaces  (b) Intersection surface  (c) Union surface  Figure 2.2: Intersection and union surfaces shown for two prostate surfaces  2.4.2  Intersection and Union of the Surfaces  The programming routines based on Computational Geometry Algorithms Library (CGAL)[2] provided by Orcun Goksel are used to calculate the union and intersection surfaces. The  CGAL  routines create a Selective Nef Complex (SNC) from the  input triangulated surfaces [28]. The SNC are generalized Nef-polyhedra, obtained from partitioning space by various planes, with labels attached to them. These labels are boolean (in, out) and are called set selection marks. The SNC data structure can be used to perform binary boolean set operations such as complement and intersection. These operations can be used to deduce results for other set operations such as union, difference, and symmetric difference. Details of this procedure can be found in [28]. If A and B are two surfaces, and their intersection surface SA∩B and relative complement surfaces SB\A and SA\B are known then union is given by: SA∪B = SA∩B + SB\A + SB\A  (2.20)  The intersection and union surfaces for two prostate surfaces are shown in Figure 2.2.  2.4.3  Calculation of Surface Force  For the calculation of forces, we assume the presence of a negative pressure (called the pressure from now onwards) that exists in the non-overlapping regions of the surfaces relative to the overlapping volume and the complement of the union the  24  Figure 2.3: Overlapping volume and non-overlapping region are shown in white and gray colors respectively. Consider A as the target surface and B as the source. The forces produced as the result of the negative pressure P -ve in the gray regions are shown by the arrows. two surfaces. This pressure results in surface forces on the source surface that increases the overlapping region while reducing non-overlapping region. Please see Figure 2.3 for an illustration of direction of surface forces due to this pressure. As described in Section 2.2.1, the surface of the mesh is formed by the outward facing triangles of the tetrahedra. We call them surface triangles. Each surface triangle experiences a force due to the presence of the pressure. The magnitude of this force depends on the area of the triangle. The direction of the force is along the normal to the face and can be outwards or inwards with respect to the mesh depending on whether that face lies in the region enclosed by the target surface or not. These directions will be such that the force increases OV and decreases NOV. This force will be distributed equally to the surface nodes that form that surface triangle. Every node receives contribution from its neighboring surface triangles. If n(f) is a vector representing the face normal of a triangular face ’f’ in the outward direction to the surface, and A(f) is the area of the face and P is the pressure, then the force on the face Ft ( f ) is given by: Ft ( f ) = D × A( f ) × P × n(f)  (2.21)  where direction D is determined by whether the face is part of intersection surface Sint or union surface Suni : 25  Figure 2.4: A number of connected surface triangles. The tetrahedral structure is ignored for ease of visualization. The pressure results in force on a triangular face with direction normal to it (direction shown with dashed lines). The force at the node (direction shown with solid line) connected to all these faces is the sum of these forces.  D=  +1 if t ∈ Sint −1 if t ∈ Suni  This force is then divided equally to the nodes forming the triangular face. Every node is connected to several faces and therefore the total force on the node is the sum of forces on all the faces. See Figure 2.4 for an illustration. If a node is surrounded by n faces then: 1 Fn(p) = ∑ × Ft ( fn ) n 3  (2.22)  The vector of nodal forces can be formed from Equation 2.22 such that: F = [ f nx1 , f ny1 , f nz1 , ...., f nxm , f nym , f nzm ]T  (2.23)  where m is the total number of nodes and f nx, f ny and f nz are x, y and z components of Fn. Figure 2.5 shows distribution of nodal forces in a mesh.  26  Figure 2.5: Force distribution in a mesh from real patient data. The mesh is color-coded with the magnitude of nodal forces. The direction of forces can be seen from the arrows originating from the surface node of the mesh.  27  2.5 Computational Model The solution of a differential equation is found by integration. The term time integration is used to discretize time in order to solve an equation such as 2.17. There are two popular time integration methods used for simulation of elastic models: explicit and implicit. The explicit and the implicit methods differ in whether the deformation is calculated by applying internal forces estimated at the previous or at the current iteration. Thus, at iteration k + 1 the displacement uk+1 is estimated by applying internal forces calculated at iteration k for explicit iterative methods and at k + 1 for implicit or semi implicit methods. The implicit and explicit methods use forward and backward Euler finite differ˙ Implicit and semi-implicit methods ence methods respectively to estimate u¨ and u. offer better stability at the expense of high computation time. The high computation time is due to the calculation of the inverse of the stiffness matrix. However, implicit methods can handle larger time steps as they are unconditionally stable. At this point in our research we have chosen a semi-implicit method for time discretization of our dynamic model which is described below. The advantage of semi-implicit over standard Euler is that it conserves energy. Consider the Equation 2.17 again: Mu¨ + Bu˙ + Ku = F For simplicity’s sake we ignore the mass and set damping to identity we get: u˙ + Ku = F  (2.24)  Using backwards Euler finite difference method, if τ is the time step:  ⇒ ⇒  uk+1 − uk + K uk+1 = Fk τ uk+1 − uk + τ K uk+1 = τ Fk (I + τ K)uk+1 = uk + τ Fk  28  (2.25)  Use a backward time step: (I + τ K)uk = (uk−1 + τ Fk−1 )  (2.26)  where uk is the array of nodal displacements at iteration k, Fk−1 is the driving force calculated at iteration k − 1, I is identity matrix and τ is the time step. The time step should be small enough to realize the system through finite differences. We assume no initial strain and therefore the initial condition u( 0) is a zero vector. The iterations end when convergence criteria is met.  2.5.1  Chattering of Forces  Our registration scheme involves repeated iterations as explained above. During the registration process, as the source and the target surfaces become close to each other, chattering of forces can be observed as the surfaces cross over and the forces change directions abruptly. This was noticed in our preliminary experiments with registration of two cylindrical. The surfaces had equal volumes and were aligned in the manner similar to a piston inside a cylinder. The goal of the registration was to bring the piston, or the moving cylinder, inside the fixed cylinder so that the two surfaces overlap completely. High stiffness parameters (Young’s Modulus of 100KPa) were used to simulate rigid registration. If the sides of the cylinders are completely aligned, the net amount of force that moves the piston depends on the area of the circular face of the piston inside the cylinder and therefore remains the same no matter how near the surfaces are to equilibrium. As our model has no collision detection and the surface is penetrable, if the time step is not small enough, this results in surface crossing over so that part of the piston is outside the fixed cylinder. Now the amount of force on the piston will be the same as before, however the direction would have changed to opposite. The next iteration will bring the piston back inside the fixed cylinder and the one after it back outside. These sudden changes in direction of force and resulting oscillations were addressed by filtering. Low pass filtering of the forces produces sufficient averaging to insure oscillationfree convergence of our dynamic model. The following finite impulse response filter is used: 29  Fk = Fk−1 + 0.9 × Fk  (2.27)  where k is the iteration at which force is being calculated.  2.5.2  Convergence Criteria  Ideally the algorithm converges when equilibrium is achieved i.e., when the mesh stops deforming. The change in  NOV(%)  is an indication of any changes in the  mesh and therefore can be used as a convergence criteria. Our algorithm stops when the change in NOV(%) is less than a pre-specified value we call the limit.  2.5.3  Time Step Reduction  As the system is iterative, the time step should be small enough for stability and large enough for speed. The NOV(%) decreases almost linearly at the start and therefore a larger time step can be used. However near convergence, as the surfaces come close to each other, this time step can be too large to find the equilibrium and therefore a smaller time step should be used. We use a scheme to reduce the initially large time step that follows the following simple rule: The NOV(%) should always decrease. With a large time step, as surfaces come close and cross over, NOV(%) can increase and therefore indicate that the time step was too large. In this case, the algorithm steps back to the last iteration and decreases the time step by half. The time step can be further decreased if the NOV(%) still increases.  2.6 Discussion This chapter provided the details of the registration model and algorithm used for registration. We have chosen a linear elastic mesh instead of a nonlinear mesh due to consideration of speed. In future non-linear elastic models can be evaluated for any improvement in the registration. Time integration is an important issue. Implicit methods have an inherent global truncation error which makes them inaccurate especially for large time steps. Better numerical methods can be used in future. It should be noted that the pressure and stiffness are relative to each other. This  30  means that no matter how soft a body is, if the force is too low it will not be able to deform it. Therefore the right amount of pressure is important and should be chosen carefully. This topic is covered in Section 3.3. The basis of linear elastic model is Cauchy’s infinitesimal strain tensor (see Section A.2) which is not rotation invariant. The assumption of small deformations is therefore key to the accuracy of the model. In cases where large deformations are involved, rotation compensation should be performed [13, 43]. Not doing so results in ghost forces which can inflate the elements and have an overall effect of increased volume of the surface [41]. To understand ghost forces, one must realize that pure rotation of an element with no change in shape should not produce strain. However, as our model right now has no way of knowing whether displacement of nodes is due to rotation or change in shape of an element therefore, internal forces based on position change of the nodes are produced not taking into account any rotation involved. These forces are unrealistic and are called ghost forces. This is certainly a shortcoming of our model in cases where the mesh undergoes large rotations during registration.  31  Chapter 3  Implementation There are various implementation considerations for the model described in the Chapter 2. The image data is available usually in the form of a stack of uniformly spaced 2- D images. These images are segmented to mark the boundary of the organ of interest and landmarks. The contours resulting from segmentation are used to create surfaces. These surfaces are then used by the registration algorithm which creates a mesh from the source surface. An overview of our algorithm is given in Figure 3.1. This chapter describes these steps in detail.  3.1 Surface Generation Physicians form a 3- D image volume of the prostate anatomy by scanning 2- D images equally spaced along a third dimension. This forms a 3- D array of 2- D images. Please see Figure 3.2 for an illustration of how the prostate image vol-  Moving Surface  Find Intersection and Union of Surfaces  Model as Elastic Solid  Calculate External Force from Pressure  No Deformed Surface  Yes  Stop?  Calculate Non Overlapped Volume  Deform Elastic Solid  Figure 3.1: Block diagram showing deformable registration algorithm 32  ume is obtained using a  TRUS  probe. Physicians then segment the organs and  areas of interest in medical images for treatment planning and surgery purposes. Once the segmentation is done, the contours can be interpolated to form 3- D surfaces. We have used the software “Stradwin” [58, 59] for this purpose. The initial image data is available as a set of 2- D images in Digital Imaging and Communications in Medicine (DICOM) format. The DICOM images are converted to Stradwin compatible format using a tool developed in house by Orcun Goksel. Stradwin opens the images as a 3- D array and can create surfaces with various resolutions and smoothing strength from the contours. We generally use the lowest resolution option available with low smoothing strength. The created surfaces are made up of triangular elements and can be visualized in Stradwin through surface rendering. The surfaces are saved on the hard drive of the computer in Virtual Reality Modeling Language (VRML) file format directly from Stradwin. A screen shot of Stradwin with images containing the segmented contours of the prostate, the urethra and the rectum is shown in Figure 3.3. The surfaces created by Stradwin are also shown in the figure.  3.2 Mesh Generation The saved surfaces can be opened in “Matlab” using another function developed inhouse by Orcun Goksel. The surface is represented in Matlab by two parameters: a list of positions in Cartesian coordinates of nodes or vertices of the triangles that make up the surface and a list of arrangement showing how they are connected to form triangular faces. The surface can be visualized using the “patch” command in Matlab (see Figure 3.4). The volume enclosed by the surface designated as source is then meshed using a software called “TetGen” [52]. Tetgen uses Delaunay tetrahedralization to divide the volume of the surface in tetrahedral elements. The Delaunay tetrahedralization ensures that the mesh is conformal i.e., that the elements do not intersect each other. This process introduces internal nodes which are distributed in the volume of the mesh. The outward appearance of the surface is the same as before the meshing. The triangular elements of the surface become the outward faces of the tetrahedra that lie close to the surface. Please see Figure 3.5 for an illustration. The created  33  Figure 3.2: Prostate imaging using TRUS probe. This image shows acquisition of transverse images spaced along axial direction. mesh is represented in Matlab by a list of nodes presented by their positions in Cartesian coordinates. The list consists of surface nodes and internal nodes, and a connection matrix that describes how nodes are connected to form tetrahedra.  3.3 Stiffness Parameters and Amount of Pressure Once the mesh is created, suitable stiffness parameters are use to create a corresponding stiffness matrix K. As described in Section 2.3, the stiffness parameters and the boundary forces are used determine the deformation. We have used realistic values for both. For prostate surface registration, the mean shear modulus for peripheral zone reported in [32] was used for homogeneous stiffness. In cases where Vibro-elastography (VE) data was available, we used Young’s modulus estimated from Transfer Function (TF) estimated from motion of tissue.  34  Figure 3.3: A screen shot of Stradwin. The screen shows four different view. The top left view shows a transverse image pelvic anatomy with segmented contours. This image is part of a series of images that are spaced in the axial direction. The top right view shows surfaces created from the segmented contours. The bottom left shows the interpolated coronal plane view and bottom right shows the interpolated sagittal plane view  3.3.1  Stiffness from Transfer Function Images  The transfer functions are available as series of 2-D images in which intensity is assumed to be directly proportional to the stiffness. As absolute values of Young’s modulus from elastography can not be extracted at present, we used a linear relationship that gives the Young’s modulus in the range of that reported [32] for the peripheral and central zones. In  VE ,  the axial strain is computed from one radio frequency (RF) data frame  to the next using a correlation based method. In cases where the correlation coefficient is low, the data is deemed unreliable. This unreliable or “bad” data appears as bright patches in the transfer function image which, when mapped linearly to Young’s Modulus can result in false low stiffness areas in the mesh. We remove the TF image data corresponding to low correlation coefficients and replace it with values based on inverse-distance weighted interpolation. In addition to that we also 35  Figure 3.4: Triangulated surface in Matlab smooth the 3- D  TF  image array with 5x5x5 box filter.  The prostate has been segmented in  TF  images [38]. From the contours we  first create a surface and then a mesh using the techniques described in Section 3.1 and Section 3.2. We consider the mesh to be embedded in the 3- D TF image array. We assign to each pixel within each tetrahedron of the mesh the average value of intensities inside the tetrahedron. We find the average value of intensities of pixels present inside each tetrahedron of the mesh and assign them that value. The Young’s modulus for the tetrahedra is estimated using these average values. We visualize the range of stiffness in the prostate by color-coding the mesh with the magnitude of estimated Young’s moduli. Such a color-coded mesh is shown in Figure 3.6b. We find the mapping which shows the color in the peripheral and central zones corresponding to the values from [32] in the color bar. The following relationship was found to give the Young’s modulus values in the range found in literature: 36  (a) Exterior of a mesh  (b) Interior of a mesh  Figure 3.5: Exterior and interior of a mesh (a) shows an opaque mesh with triangular facets of tetrahedra forming the surface (b) shows a mid plane cutting the mesh to reveal the tetrahedra inside  Young′ s Modulus = 106 × Trans f er Function Image Intensity  (3.1)  This mapping is applied to average intensity calculated for each tetrahedra. We assumed that the tissue is nearly incompressible and therefore used a Poisson’s ratio of .499 for all cases. This might not be the case if there is a fluid outflow from the organ in the time between the two sets of images are taken. However, as volume preservation is one of the goals of registration, the high incompressibility of the mesh is desirable.  3.3.2  Amount of Pressure  The amount of pressure is set so that the net force resulting from it on the mesh is of the same order of magnitude as the weight of the organ being modeled. An average prostate has a mass of 11gms [34] and a weight of 0.11N. For all of the cases of prostate surface registration we have performed, a pressure of 1KPa was enough to produce a force in the order of 10−1 N. An increase in TRE was observed 37  (a) VE transfer function image  (b) Corresponding Young’s Modulus  Figure 3.6: Young’s modulus from transfer function images (a)VE transfer function prostate image re-sliced at coronal plane. The plane passes through the center of mass of the surface. (b) The mesh is color-coded to Young’s modulus to show the variation of stiffness inside the prostate. Please note the range of Young’s modulus on the colorbar. when for higher pressure in our early results as shown in Appendix E.  3.4 Registration Algorithm We start the registration by aligning the centers of mass of the two surfaces to bring them in the same Cartesian frame. Next we perform Principal Axis Alignment (PAA) to recover the deformation due to rotation. The overview of algorithm is given in Algorithm 1. We used homogeneous stiffness where elastography data was not available. As the computation of the inverse term in Equation 2.26 is a time consuming task, we compute it only once for every time step τ . The limit (see Section 2.5.2) for change in NOV (%) for convergence was kept at 1 × 10−4 . Once the surfaces are registered, the positions of the landmarks in the deformed mesh are obtained by using the natural coordinates of points on segmented landmarks.  38  Algorithm 1 The registration algorithm Let the target organ shape be represented as a triangular surface St with its center of mass at the origin Let the source organ shape be represented as a triangular surface Ss with its center of mass at the origin Align principal axis of Ss to that of St Let the Ss be represented as a tetrahedral mesh with nodes at positions p0 Calculate stiffness matrix K for the source mesh. Let time step τ , pressure P and standard for convergence limit be known. Calculate the intersection and union of surfaces St and Ss Calculate NOV0 Calculate external force F0 A = (I + τ K)−1 k←0 repeat k ← k+1 uk−1 = pk−1 − p0 uk = A · (I uk−1 + τ FK−1 ) pk = pk−1 + uk Get surface Ss k from pk Calculate intersection and union of surfaces St and Ss k Calculate NOVk (%) Calculate external force Fk Fk = 0.9 × Fk−1 + Fk if NOVk > NOVk−1 then k ← k−1 τ ← τ /2 A = (I + τ K)−1 end if until (NOVk−1 (%) − NOVk (%)) ≤ limit  39  3.5 Finding the Position of Landmarks in the Deformed Mesh As the source surface deforms, the volume interior to the surface deforms too. Thus, in the registered surface, the landmark has moved and has a new position. For determining this position, we use the natural coordinates of the points contouring the landmark. For every point po = (xo , yo , zo ) on the landmark contour(s), we find the tetrahedral element e from the original source surface mesh that contains that point and find its natural coordinates by using:    −1   N1e 1 1 1 1 1  e  o o o o  o N2  x1 x2 x3 x4  x   =    N e  yo yo yo yo  yo   3   1 2 3 4   N4e zo1 zo2 zo3 zo4 zo  (3.2)  where (x1o , yo1 , zo1 ), (x2o , yo2 , zo2 ), (x3o , yo3 , zo3 ), (x4o , yo4 , zo4 ) are the positions of the four nodes forming the tetrahedra inside which the point po lies. We use the same relationship to find the Cartesian coordinates of the landmark point pd = (xd , yd , zd ) in the deformed mesh:      1 1 1 1 Ne 1  d   d d d d   1e  x  x1 x2 x3 x4  N2     = yd  yd yd yd yd  N e     1 2 3 4 3 zd zd1 zd2 zd3 zd4 N4e  (3.3)  where (x1d , yd1 , zd1 ), (x2d , yd2 , zd2 ), (x3d , yd3 , zd3 ), (x4d , yd4 , zd4 ) are the new positions of the corners of the same tetrahedra. Once the positions of all the landmark contour points in deformed mesh is known, the registration error (see Section 3.7) can be calculated.  3.6 Displaying the Result in 2-D Images While the registered surface can be enough for use in the volume studies like quality assurance for brachytherapy, 2- D images are often better suited for observing the results of registration. In order to obtain 2- D images, the images from which 40  the source surface was derived can be deformed according to the displacement field obtained from registration. This is done in a manner very similar to the one described in Section 3.5. The undeformed source mesh is brought into the image coordinate frame and the natural coordinates of all the pixels of the images with respect to the tetrahedron they are present inside are calculated. As the mesh exists only inside the surface, all the pixels outside the surface are ignored and assigned a constant value. Using 3.2 and 3.3, the new position of each pixel in the deformed image is calculated. These new scattered pixel positions are linearly interpolated to find a 3- D grid which is a stack of 2- D deformed images.  3.7 Validation Measures The validation measures for determining the accuracy of our registration are chosen in accordance with our goals as described in Section 1.4: 1. Dice Similarity Coefficient (DSC) 2. Target Registration Error (TRE) 3. Volume Change (%) (VC) The DSC can be used to measure overlap of the two surfaces. In context of our problem, it can be defined by: DSC =  2 × NOV VolSource +Voltarget  (3.4)  A DSC of 1 indicates complete overlap. This measure tells how well the boundaries of the two surfaces match. A surface with thousands of degrees of freedom can match another surface in various ways. To determine accuracy of registration, the concept of  TRE  is very  commonly used. This is done by finding a landmark i.e., a feature that is visible in both sets of images and calculating the distance between the landmarks as visible in target images and registered source images. The lower the  TRE ,  the better a  registration is. Volume preservation is often desired in registration specially for cases where the volume of an organ is not known to change in between the time it is imaged by 41  the modalities involved. In our registration we use a nearly incompressible elastic model and therefore the volume of the source surface should not change after the registration. To test if it has changed we use the Volume Change (VC)(%) defined simply as: Volume Change(%) =  Volsource(a f ter) −Volsource(be f ore) × 100 Volsource(be f ore)  (3.5)  where Volsource(be f ore) and Volsource(a f ter) are volumes of the surfaces before and after registration. For a perfectly preserved volume, this value should be zero.  3.8 Discussion Image and mesh quality and resolution can influence registration in terms of accuracy as well as computational performance. These issues are discussed in detail below.  3.8.1  Error in Marking Boundary and Landmarks  The segmentation errors in outlining the organ can be caused by low visibility in the images or inexperience of the person performing the segmentation. The boundary of the organ thus marked inaccurately results in a surface that is not true to the surface of organ. This directly affects the accuracy of any surface based registration method. As described in Section 3.7, the accuracy of registration is measured with TRE which is a distance measure between the marked landmarks after registration. The segmentation error in marking the landmarks can results in a value not representative of the accuracy of registration.  3.8.2  Error Due to Difference in Slice Thickness  Slice thickness or image spacing in the z direction (which is axial direction for transverse images) is different in every modality. Therefore there cannot be exact correspondence in the position of a specific landmark seen in 2- D images of two different imaging modalities. The error resulting from these is related to the similarity of the slice position along z axis. To avoid this 3- D structural landmarks 42  (surfaces, curves etc.) can be used.  3.8.3  Difference in Volume  One of the goals of our deformable registration is volume preservation. If accomplished along with the goal of accurate registration, this requires surface volumes to be similar. Errors in segmentation and/or loss of fluids from the tissue can result in difference in volumes of surface obtained from different modalities at different times. If such difference exists, the two goals cannot be achieved at the same time and some kind of compromise is necessary to obtain a result.  3.8.4  Discretization Error  When a continuous object is approximated by discrete elements, an error called discretization error arises which accounts for the differences between the actual shape and the  FEM  model. The discretization error can be reduced by using a  higher number of elements. However using highly refined meshes increases computational burden and therefore a compromise must be found.  3.8.5  Computational Performance  We implementated of our algorithm in Matlab where it took 20-45 minutes for a registration to complete. The computation speed depends on the resolution of mesh. It can be made quicker using optimized coding in C++. This might also allow the use of high resolution meshes for increased registration accuracy.  43  Chapter 4  Evaluation and Application In this chapter, we first study impact of using non-homogeneous stiffness on our registration method. Then we evaluate our registration technique using prostate surfaces derived from  TRUS  and  MR .  After we establish it as an effective de-  formable registration method, we use it to register prostate VE images to pathology images using stiffness from the former. We find another application of our algorithm in registration of Kidney surfaces which is covered in Section 4.4.  4.1 Advantage of Using Non - Homogeneous Stiffness We performed a series of simulations to evaluate the importance of using the stiffness values close to that of tissue for registration accuracy. The tissue deforms under an external force according to stiffness distribution in the tissue. We create this scenario by deforming an ellipsoid by registering it with an ellipsoid of different shape. The registration is carried out by using an elasticity map created for the ellipsoid which contains the average elasticity value for the peripheral zone [32] everywhere inside the surface except for a blob which was four times stiffer than its surroundings. The registered or deformed ellipsoid surface is then registered to its undeformed version in an attempt to remove the deformation. This process is carried out first with homogeneous stiffness and then with stiffness from the elasticity map (one with stiff blob). The results are evaluated on three basis: mean distance be-  44  tween surface nodes, mean distance between internal nodes corresponding to the stiff blob and the NOV(%). For a perfect deformation recovery, all three should be zero. We performed this test for three cases: 1. Ellipsoid surface containing stiff blob in the center of the surface deformed to a different shape. 2. Ellipsoid surface containing stiff blob off the center of the surface deformed to a different shape. 3. Ellipsoid surface containing stiff blob off the center of the surface rotated and deformed to a different shape. The results are shown in Table 4.1. The results clearly show that non-homogeneous stiffness performs better in deformation recovery specially for internal stiffer nodes.  4.2 Prostate MR and TRUS Surface Registration We tested our registration algorithm on four pre-operative and fourteen post-operative MR  and TRUS transverse prostate image data sets. The pre-operative images have segmented prostate. We were able to identify  and segment common structures such as calcifications and the urethra on the images. These were used as landmarks for validation of the registration. All the post-operative images have the prostate and urethra segmented. As  VE  data was  not available for these images, we used homogeneous stiffness (see Section 2.3.1) with Young’s modulus of 10KPa. This is the mean Young’s modulus for the peripheral zone of prostate reported in [32]. The Poisson’s ratio of ideally incompressible material is 0.5. Using this value however causes instability in our model and therefore we use a Poisson’s ratio of 0.499. The MR surface was used as target and the US  surface was used as the source. Once the surfaces were registered, the land-  mark or urethra position in the deformed US was obtained in the manner described in Section 3.5. We used a pressure of 1KPa which resulted in a force in the order of 10−1 N.  45  Pressure  5000 46 1000  Case 1 2 3 1 2 3  Non-homogeneous Stiffness Surface Distance(mm) Blob Distance(mm) 0.015 ± 0.008 0±0 0.015 ± 0.008 0.000 ± 0.000 0.043 ± 0.018 0.000 ± 0.000 0.002 ± 0.001 0±0 0.002 ± 0.001 0.000 ± 0.000 0.146 ± 0.027 0.002 ± 0.001  NOV(%) 0.25% 0.25% 0.15% 0.11% 0.12% 0.43%  Homogeneous Stiffness Surface Distance(mm) Blob Distance(mm) 0.502 ± 0.738 0.009 ± 0.013 0.502 ± 0.738 0.015 ± 0.023 0.717 ± 0.619 0.003 ± 0.003 0.484 ± 0.79 0.009 ± 0.014 0.491 ± 0.794 0.014 ± 0.023 0.135 ± 0.019 0.002 ± 0.001  Table 4.1: Results for deformation recovery test using homogeneous and non-homogeneous stiffness  NOV(%) 6.39% 6.39% 0.48% 6.82% 6.81% 0.40%  4.2.1  Patient Data  We have two separate sets of prostate MR and TRUS images obtained for two different studies. The images in the two sets belong to different patients going through different procedures. The pre-operative images are taken from patients going for LDR  brachytherapy. The images were taken a few weeks prior to the procedure.  In these images, the slice thickness in  MR  and  US  is 4mm and and 5mm respec-  tively and the in-plane resolution of MR and US images is 0.27 mm/pixel and 0.16 mm/pixel respectively. We have fourteen sets of post-brachytherapy prostate  MR  and  TRUS  images.  These were collected from patients that went through prostate brachytherapy at the BC Cancer Agency. The original purpose of data collection was to study the feasibility and comparability of post-implant  US / CT  image fusion with  MR / CT  image  fusion for the purpose of post implant dosimetry [47]. The TRUS data was collected at the end of the treatment and the MR images were collected the same day. A Foley catheter was placed in the urethra to ease image fusion. The prostate and the urethra were segmented in all of the images. In these images, the slice thickness in MR and US is 3mm and and 5mm respectively and the in-plane resolution of MR and US images is 0.55-0.70 mm/pixel and 0.17-0.18 mm/pixel respectively. In US, the prostate is observed to be deformed due to the presence of the TRUS probe and change in the patient position. Our registration attempts to recover the prostate surface from this deformation.  4.2.2  Landmarks  In pre-operative images, we found common intra-prostatic structures like the urethra and calcifications visible in both the modalities as landmarks. The landmarks were segmented and the position of landmarks after registration was determined using the method described in Section 3.5. It must be noted that due to the slice difference, the position of landmarks seen in two image planes from cannot be the same along the z-axis. Considering this, the  TRE  MR  and  US  reported is the  distance in the x-y plane between the centers of the landmarks. In post-operative images, due to the presence of the Foley catheter, the urethra is visible in both  MR  and  US .  The segmented urethra in both sets of images was  47  used as the landmark. The segmented contours of the urethra from both sets of images were read into a Matlab function. The points on these contours are often sparse, therefore, we fit ellipses on each of the contours and found the centers of these ellipses. These centers represent the center of the urethra in the corresponding image planes. We fit a 3- D spline passing through the centers of the ellipses for both the  MR  and the  surfaces to form the urethral path. After the  TRUS  surface is registered to the  MR  TRUS  surface, the position of each point on the ellipse  corresponding to the TRUS urethra is found in the deformed TRUS mesh using the method described in Section 3.5. This gives “deformed” ellipses corresponding to the urethra in the deformed mesh. The center of each deformed ellipse is found and a 3- D spline is fit to find the urethral path. The  TRE  reported is the mean and  maximum distance between urethral paths calculated for MR images and registered TRUS  4.2.3  mesh.  Results  We performed registration for all the sets first without and then with  PAA  as the  initial step before performing deformable registration to test the sensitivity of our registration technique to initial alignment. Table 4.2 shows the results of registration for the four pre-operative image sets. The TRE improved from 2.1±0.9mm to 1.7±0.8mm when PAA was used and from 2.9±2.3mm to 1.3±0.8 when PAA was not used. In both cases the DSC improved from 0.90±0.02 to 0.96 ±0.01. The maximum VC in the source surface was 1%. Case  VD(%)  1  0.3  2  -8.3  3  -4.1  4  -3.6  PAA Yes No Yes No Yes No Yes No  DSC Before After 0.88 0.95 0.89 0.96 0.90 0.95 0.90 0.95 0.92 0.97 0.93 0.97 0.90 0.97 0.89 0.97  TRE(mm) Before After 3.4 2.5 6.2 2.5 2.0 1.2 2.5 1.2 1.8 2.3 1.6 0.7 1.2 0.9 1.2 1.0  VC(%) -0.9 -1.1 -1.0 -1.0 -0.2 -0.7 -1.0 -0.9  Table 4.2: Results for the four pre-brachytherapy images. 48  Case  VD(%)  1  2.7  2  2.5  3  0.8  4  2.7  5  3.2  6  3.8  7  1.9  8  2.6  9  4.8  10  -2.1  11  -2.6  12  -4.1  13  2.6  14  -4.4  PAA No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes  DSC Before After 0.81 0.95 0.81 0.95 0.89 0.95 0.89 0.95 0.89 0.94 0.89 0.94 0.94 0.98 0.93 0.98 0.85 0.95 0.88 0.94 0.92 0.96 0.92 0.96 0.92 0.98 0.93 0.97 0.81 0.91 0.80 0.91 0.86 0.95 0.86 0.94 0.90 0.96 0.90 0.96 0.88 0.96 0.85 0.95 0.94 0.97 0.93 0.96 0.93 0.97 0.94 0.98 0.91 0.96 0.91 0.95  TRE Before (mm) Mean Max 0.8±0.5 2.7 1.6±1.4 4.9 1.0±0.7 2.6 1.0±0.7 2.6 0.7±0.5 1.7 0.7±0.5 1.7 1.1±0.9 3.0 1.4±1.2 3.4 2.3±1.6 6.4 0.7±0.5 2.2 1.3±0.9 3.0 1.2±0.8 3.4 4.4±2.4 7.7 4.5±1.5 6.3 2.9±2.3 8.8 2.0±2.4 7.6 1.9±1.5 4.7 3.5±1.9 7.5 0.9±0.7 2.5 0.9±0.7 2.5 1.6±0.9 3.4 2.6±1.7 6.0 1.6±1.0 4.4 2.6±1.2 4.9 2.3±1.2 3.7 1.9±0.9 3.0 2.9±2.2 6.6 3.3±2.4 7.4  TRE After (mm) Mean Max 0.7±0.5 1.9 0.8±0.6 1.9 1.7±1.4 4.5 1.8±1.5 4.8 0.9±0.8 2.6 0.9±0.8 2.6 0.9±0.8 2.9 0.8±0.9 2.9 0.7±0.6 2.0 0.5±0.4 1.6 0.5±0.4 1.8 0.6±0.5 2.3 3.8±1.0 5.5 3.7±1.0 5.2 1.0±1.1 4.4 1.2±1.2 5.0 3.0±0.4 3.9 3.5±0.7 4.7 1.1±0.7 2.0 1.0±0.6 1.9 1.4±0.8 2.9 1.8±1.0 3.7 1.7±0.7 3.4 1.9±1.0 4.2 1.3±0.7 2.1 1.6±0.8 2.6 3.0±2.1 6.3 2.9±2.1 6.2  Table 4.3: Results for the fourteen post brachytherapy images. The results for post-operative data sets are presented in Table 4.3. When was used, after registration, the  DSC  PAA  improved from 0.89±0.05 to 0.95±0.02 and  TRE improved from 2.0±1.8mm to 1.6±1.4mm and the change in volume of source  surface was 1.9±1.4%. When improvement in  DSC  PAA  was not used, after registration, we saw an  from 0.89±0.04 to 0.96±0.02 and an improvement in 49  TRE  VC(%) 1.4 2.0 0.0 0.5 2.4 2.5 1.2 2.3 3.7 -0.9 -0.3 -0.3 2.3 1.5 1.2 4.6 4.8 5.7 -1.2 -1.2 -0.2 3.4 -0.3 0.3 -0.3 1.1 0.4 1.0  (a) Before Deformable registration  (b) After deformable registration  Figure 4.1: Result of our deformable registration on patient data case number 5. The pink and green meshes represents the target and the source surfaces respectively. The urethra from both the volumes are shown by the fitted splines of the corresponding colors. PAA was not performed in this case. from 1.8±1.7mm to 1.5±1.4mm. The  VC  in source surface during registration  was 1.5±1.4%. The Volume Difference (VD) between source and target surfaces is shown as a reference. Please see Figure 4.1 for registration result of one case.  4.2.4  Discussion  In three cases out of fourteen post-operative images and in one case out of four preoperative images, the (maximum) PAA  TRE  changes by more then a millimeter when  is not used. This shows that our registration method can be sensitive to initial  orientation. The main contribution of our model is the novel determination of boundary conditions. We, therefore, compare our registration method with a method presented in [23] that also uses linear elastic model and semi implicit time integration. This method has been evaluated by registering ten pairs of MR prostate surfaces [7]. It uses gradients on distance map created for two surfaces to drive its  FEM  elastic  model. We implemented this method by using the same computational model as 50  (a)  (b)  (c)  Figure 4.2: Case 9. The pink and green meshes represent the target and source surfaces respectively. The urethras from both volumes is shown by the splines of the corresponding colors. (a) shows unregistered surfaces, (b) shows surfaces registered using our approach, and (c) shows surfaces registered using [23] our method. The boundary conditions were calculated using a Matlab function created by Orcun Goksel that used Matlab’s built in binary distance map and gradient functions. We did not use low pass filtering as there was no need for it. For the fourteen patients, we obtained a 1.8±1.7mm when  PAA  TRE  of 1.7±1.6 mm when  was used. The  volume was 2.0±1.6% when  PAA  DSC  PAA  was not used and  was 0.97±0.01 and change in source  was used. These values were 0.97±0.01 and  1.7±1.4% respectively when it was not used. Please see Figure 4.2 for differences in registration results for one case. The detailed results are presented in Table 4.4. In terms of TRE we see that our method performs equally or slightly better than the distance based method.  4.3 Pathology and Vibroelastography Registration We used our registration algorithm to register five pairs of VE transfer function images obtained from vibroelastography of the prostate to pathology images. The prostate was segmented in both. In pathology, the cancer contours were also segmented. The prostate surface from transfer function images was used as the source surface and the surface from pathology was used as the target surface. The  51  Case  VD(%)  1  2.7  2  2.5  3  0.8  4  2.7  5  3.2  6  3.8  7  1.9  8  2.6  9  4.8  10  -2.1  11  -2.6  12  -4.1  13  2.6  14  -4.4  PAA No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes No Yes  DSC Before After 0.81 0.96 0.81 0.95 0.89 0.98 0.89 0.98 0.89 0.98 0.89 0.98 0.94 0.98 0.93 0.98 0.85 0.97 0.88 0.97 0.92 0.98 0.92 0.97 0.92 0.98 0.93 0.98 0.81 0.95 0.80 0.95 0.86 0.96 0.86 0.96 0.90 0.98 0.89 0.98 0.88 0.96 0.85 0.96 0.94 0.98 0.93 0.98 0.93 0.99 0.94 0.98 0.91 0.98 0.91 0.97  TRE Before (mm) Mean Max 0.8±0.5 2.7 1.6±1.4 4.9 1.0±0.7 2.6 1.0±0.7 2.6 0.7±0.5 1.7 0.7±0.5 1.7 1.1±0.9 3.0 1.4±1.2 3.4 2.3±1.6 6.4 0.7±0.5 2.2 1.3±0.9 3.0 1.2±0.8 3.4 4.4±2.4 7.7 4.5±1.5 6.3 2.9±2.3 8.8 2.0±2.4 7.6 1.9±1.5 4.7 3.5±1.9 7.5 0.9±0.7 2.5 0.9±0.8 2.5 1.6±0.9 3.4 2.6±1.7 6.0 1.6±1.0 4.4 2.6±1.2 4.9 2.3±1.2 3.7 1.9±0.9 3.0 2.9±2.2 6.6 3.3±2.4 7.4  TRE After (mm) Mean Max 1.3±1.1 3.3 1.2±1.0 2.9 0.7±0.4 1.4 1.0±0.8 2.4 0.9±0.4 1.6 0.9±0.4 1.6 1.1±1.0 3.2 1.6±1.4 4.1 1.1±0.5 1.8 0.5±0.3 1.0 0.5±0.4 2.3 0.8±0.5 2.4 4.2±1.7 6.8 3.8±1.2 5.4 1.9±1.7 6.1 1.0±1.2 4.8 3.2±2.0 6.0 3.8±2.2 7.0 1.2±0.8 2.5 1.2±0.8 2.2 1.7±0.7 3.1 2.7±1.7 5.6 1.6±0.6 2.9 1.9±0.9 3.8 1.0±0.5 1.9 1.1±0.5 2.1 3.3±2.3 7.1 3.2±2.3 6.9  Table 4.4: Results for the fourteen post brachytherapy images using distance gradients based boundary conditions. badly correlated data was removed and replaced by interpolated data and stiffness matrix was created from the transfer function data using methods described in Section 3.3.1. We used a higher amount of pressure (5KPa) than that used in prostate  MR  and  TRUS  surface registration. After the surface registration, the 52  VC(%) -4.3 -5.2 -0.1 -0.9 -1.2 -1.2 -0.2 -0.3 -0.3 -1.9 0.1 -0.1 -0.5 -0.7 -4.3 -4.4 -1.6 -2.4 -3.0 -2.5 -3.0 -3.8 -1.9 -1.9 0.7 -0.2 -2.4 -2.5  source images were deformed according to the displacement field obtained using method described in Section 3.6. Slices are cut from 3- D deformed image matrix at the same depths as pathology images and cancer contours from pathology images were overlaid on them.  4.3.1  Patient Data  The VE image data used here were acquired from patients going through the radical prostatectomy at the Vancouver General Hospital. The B-mode  US  and  VE  were  collected intra-operatively, just prior to the actual prostatectomy intervention. Detail on the data acquisition and processing of this patient data set can be found in [38]. The  VE  transfer function images have in-plane resolution of 5mm/pixel and  slice spacing of 0.4mm. The pathology slices were obtained from the prostate extracted from radical prostatectomy. Using a multi-bladed device, as descibed in [21], 4mm thick sections are cut from prostate. From these sections, a fine slice is cut and mounted on a glass slide for hematoxylin and eosin staining. From each prostate, 9-13 slides, depending on the size of the gland, are produced on which cancer regions are detected and marked by a certified pathologist. The photographic images of these slices are used to create 3- D image volume which is segmented using Stradwin. The cutting of fine slice from the section can be done from anywhere in the 4mm section therefore exact slice position or slice spacing is unknown. However, we used slice spacing of 4mm to produce the 3- D images.  4.3.2  Landmarks  No visible common landmarks could be seen in pathology and strain images.  4.3.3  Results  We obtained 2- D deformed transfer function images for five pairs of patient data. The cancer contours were overlaid on the deformed images to study if the intensity differences can be used to predict the presence of cancers. Please see Figure 4.3 for example of a deformed 2- D transfer function image. Results are laid out in Appendix C. 53  (a)  (b)  (c)  Figure 4.3: A typical result of VE and pathology registration (a) shows a pathology image with segmented prostate boundary and cancer contours (b) shows a VE transfer function image at approximately the same depth with segmented prostate boundary and (c) shows a registered VE transfer image at the same depth with cancer contours from pathology  4.3.4  Discussion  These are preliminary results. As we can not find any common landmarks, the results cannot be validated for internal accuracy. We saw an improvement in  DSC  from 0.89±0.0 to 0.98±0.01.  4.4 Kidney CT and US Registration We used our registration method to register Kidney surfaces. We have five patient Computed Tomography (CT) and US image data sets with the kidney segmented in all of them. The registration was carried out in the same way as for the prostate problems. Due to the absence of elastography data, we use homogeneous stiffness using the same parameters as described in Section 4.2. The amount of pressure used results in a force in the order of 1 N.  4.4.1  Patient Data  The  and  CT  tomy.  US  US  data was collected from the patients about to go through nephrec-  was taken from the patient’s side with the patient in the flank position on  the CT table (affected kidney up). Patients were told to hold their breath during the acquisition of the US volume. Patients were then told to hold their position during CT  imaging. The  CT  image was taken approximately 5 minutes after the 54  US  vol-  ume. A contrast agent was injected prior to CT imaging, but this was not observed to affect kidney shape/volume.  4.4.2  Landmarks  Due to the low contrast and resolution in  US,  no common landmarks were identi-  fied. We were unable to calculate TRE in this case.  4.4.3  Results  We used our algorithm on five sets of patient images. Due to the absence of landmarks,  DSC ,  Hausdorff distance and  are reported. Please see Table 4.5 for  VC  results. The registered kidney surfaces are shown in Appendix D. Case  VD(%)  1 2 3 4 5  -1.2 1.7 -1.0 -4.8 3.0  Initial 0.63 0.66 0.60 0.68 0.65  DSC After PAA 0.81 0.83 0.77 0.81 0.78  Final 0.94 0.96 0.94 0.94 0.91  Hausdorff Distance (mm) Initial After PAA Final 24.9 17 16.6 18.7 9.6 6.4 27.3 14.5 10.8 21.5 14.6 13.9 27.4 22.7 14.5  VC(%) 3.6 3.1 1.4 -1.3 11.41  Table 4.5: Results of kidney surface registration For the five patient data sets, the DSC improved from 0.80±0.02 to 0.95±0.01 after deformable registration. The change in source surface volume during the registration was 2.1±2.3%. Please see Figure 4.4 for registered surfaces of case 1.  55  (a) Before Deformable registration  (b) After Deformable registration  Figure 4.4: Result of our deformable registration on kidney patient data (case 1). The pink and green meshes represents the target (CT) and the source (US) surfaces respectively.  56  Chapter 5  Conclusions 5.1 Summary of Contributions We proposed a new 3- D surface based deformable registration technique that was evaluated for eighteen patient data sets and showed reasonable results. The specific contributions and achievements of this project are: • Use of boundary conditions derived from a global non-overlap in volumes: We presented an innovative method of calculating boundary conditions which maximizes overlap between the surfaces by assuming a negative pressure in non-overlapping regions of the surfaces. We observed a torsional effect in some cases which could eliminate the need to perform a principal axis alignment prior to performing deformable registration. • Comparison with distance derived finite element based registration method: We implemented a distance based boundary conditions function and carried out registration for fourteen post-operative prostate data sets first with initial PAA  and then with out initial  PAA.  We compared the results to those from  our registration. • Achieved a compromise between high surface fit and low  TRE  to com-  pensate for segmentation errors: Our registration method found a compromise between the surface match and internal accuracy in the presence of  57  image quality related errors. This was demonstrated by high TRE .  DSC  In comparison, the other registration method found lower  and low DSC  and  higher TRE especially for the cases where rotation was involved. • First reported use of stiffness parameters from the elastography for registration: In literature, the bio-mechanical models used for prostate registration use either fixed stiffness values [11, 44] that remain same for all types of registrations or stiffness values obtained through an optimization process [5, 7, 31]. To our knowledge, we are first to report the use of stiffness value reported in literature for the prostate in a biomechanical model used for registration. In addition, we used relative stiffness parameters derived directly from elastography transfer function images for registering prostate surface from transfer function images to prostate surface from pathology images. To the best of our knowledge, this is the first time this has been reported. This shows the potential of our technique for using real stiffness values derived from elastography. • Use of 3- D curve distance between urethra for the purpose of evaluation: We used 3- D continuous curves representing urethras as landmarks. Use of such 3- D landmarks evaluate internal registration accuracy better despite of the differences in slice thickness. This unique method for calculating  TRE  has been reported for the first time. • Registration of elastography transfer function images to histo-pathology images: We performed prostate surface registration with surfaces from elastography transfer function images to histo-pathology images for five patient image sets. This work can be used for studying the relationship of stiffness in a region (obtained from elastography) to the presence of tumors in that region. • Finding deformed image from a deformed mesh: We obtained 2- D deformed transfer function images after the registration of 3- D prostate surfaces from the elastography transfer function and histo-pathology. This was done by using our finite element mesh which was used to determine position of every pixel of the underformed images in the deformed mesh to form 58  warped images. Using this method, warped images can be produced after any 3- D  FEM  registration.  Perhaps its most promising aspect of our registration method is the ability to combine the approach with elastography, so the Young’s modulus used in our registration is the actual one obtained in-vivo in one or both of the imaging modalities. Due to anatomical correctness of this registration, this work has immense potential in quality assurance of brachytherapy, treatment planning of  HDR  therapy  and localized biopsy.  5.2 Future Works • Use with other surfaces and modalities: This registration technique is not limited to a specific modality and anatomy. As long as surfaces of organ of interest are available, this method can be used. Therefore there is potential of using this technique for registration of surfaces of other organs. • Implementation of rotation compensation: As discussed in Section 2.6 our linear elastic model might be improved by addressing the issue of rotation compensation. The impact of rotation compensation and difference in efficiency of algorithm resulting from it should be explored in future. • Better time integration: The implicit time integration methods accumulate error over time. Therefore, computational performance might be improved in future by implementing better time integration and optimization methods. • Use of optimized C++ code: One draw back of our registration method is the computation speed which can be improved in future through optimizing the algorithm for speed and potentially using C++ or other programming languages. Achieving better computational performance can enable us to use higher resolution meshes and therefore increase accuracy of registration. • Presence of FEM mesh outside the surface:  Another drawback when  creating a deformed 2- D image from deformed 3- D volume is that due to the presence of finite element mesh only on the inside of the prostate surface we cannot determine the positions of pixels outside it. While this might be 59  adequate for most of the applications, the use of a multi-organ mesh might be better to create a deformed 2- D image. • Mathematical basis for comparison: We know that our method performs better when rotation is involved as compared to the other registration method we implemented. This comparison needs further analyzing. • Use of stiffness from elastography: geneous stiffness for registering  MR  In this thesis we have used homo-  and  TRUS  prostate surfaces. In future,  elasticity values from elastography can be used to perform registration with better results. • Registration of histo-pathology images to MR images: We leave this registration to future work. This can be done with stiffness derived from elastography transfer function images. The result will show the correspondence between cancer tumors as they appear in MR images and pathology slices. • Volume scaling for pathology surfaces: Prostate volume shrinkage after it is mounted for pathology is a well known issue. We did not notice much difference in volumes of the surfaces from elastography transfer function images and from histo-pathology images. This is probably due to tighter margins of segmentation in elastography. However, this should be investigated in future.  60  Bibliography [1] Canadian cancer societys steering committee on cancer statistics. canadian cancer statistics 2011. toronto, on: Canadian cancer society. → pages 1 [2] C GAL, Computational Geometry Algorithms Library. http://www.cgal.org. → pages 24 [3] Prostate cancer information from the foundation of the prostate gland. prostate cancer treatment guide. Web, August 2011. URL http://www.prostate-cancer.com. → pages 3, 4, 5, 8 [4] H. Ahmed, A. Kirkham, M. Arya, R. Illing, A. Freeman, C. Allen, and M. Emberton. Is it time to consider a role for mri before prostate biopsy? Nature Reviews Clinical Oncology, 6(4):197–206, 2009. → pages 1, 6 [5] R. Alterovitz, K. Goldberg, J. Pouliot, I. Hsu, Y. Kim, S. Noworolski, and J. Kurhanewicz. Registration of mr prostate images with biomechanical modeling and nonlinear parameter estimation. Medical physics, 33:446, 2006. → pages 12, 13, 58 [6] R. Bajcsy and S. Kovacic. Multiresolution elastic matching*. Computer Vision, Graphics, and Image Processing, 46(1):1–21, 1989. → pages 3 [7] A. Bharatha, M. Hirose, N. Hata, S. Warfield, M. Ferrant, K. Zou, E. Suarez-Santana, J. Ruiz-Alzola, A. DAmico, R. Cormack, et al. Evaluation of three-dimensional finite element-based deformable registration of pre-and intraoperative prostate imaging. Medical Physics, 28: 2551, 2001. → pages 11, 12, 50, 58 [8] R. Bourne, P. Katelaris, S. Danieletto, T. Dzendrowskyj, P. Stanwell, and C. Mountford. Detection of prostate cancer by magnetic resonance imaging and spectroscopy in vivo. ANZ Journal of Surgery, 73(8):666–668, 2003. → pages 6 61  [9] M. Bro-Nielsen and S. Cotin. Real-time Volumetric Deformable Models for Surgery Simulation using Finite Elements and Condensation. In Computer graphics forum, volume 15, pages 57–66. Wiley Online Library, 1996. → pages 3, 16, 17 [10] M. Bro-Nielsen and C. Gramkow. Fast fluid registration of medical images. In Visualization in Biomedical Computing, pages 265–276. Springer, 1996. → pages 3 [11] K. Brock, M. Sharpe, L. Dawson, S. Kim, and D. Jaffray. Accuracy of finite element model-based multi-organ deformable image registration. Medical physics, 32:1647, 2005. → pages 12, 58 [12] K. Bylund, J. Bayouth, M. Smith, A. Hass, S. Bhatia, and J. Buatti. Analysis of interfraction prostate motion using megavoltage cone beam computed tomography. International Journal of Radiation Oncology* Biology* Physics, 72(3):949–956, 2008. → pages 7 [13] J. Choi and A. Szymczak. Fitting solid meshes to animated surfaces using linear elasticity. ACM Transactions on Graphics (TOG), 28(1):1–10, 2009. → pages 12, 31 [14] G. Christensen, R. Rabbitt, and M. Miller. Deformable templates using large deformation kinematics. Image Processing, IEEE Transactions on, 5(10): 1435–1447, 1996. → pages 3 [15] T. Cootes, C. Taylor, D. Cooper, J. Graham, et al. Active shape models-their training and application. Computer vision and image understanding, 61(1): 38–59, 1995. → pages 3 [16] W. Crum, T. Hartkens, and D. Hill. Non-rigid image registration: theory and practice. British journal of radiology, 77(Special Issue 2):S140, 2004. → pages 2, 3 [17] V. Daanen, J. Gastaldo, J. Giraud, P. Fourneret, J. Descotes, M. Bolla, D. Collomb, and J. Troccaz. Mri/trus data fusion for brachytherapy. The International Journal of Medical Robotics and Computer Assisted Surgery, 2(3):256–261, 2006. → pages 11 [18] E. D’Agostino, F. Maes, D. Vandermeulen, and P. Suetens. A viscous fluid model for multimodal non-rigid image registration using mutual information. Medical image analysis, 7(4):565–575, 2003. → pages 3  62  [19] J. Dawson, T. Wu, T. Roy, J. Gu, and H. Kim. Dose effects of seeds placement deviations from pre-planned positions in ultrasound guided prostate implants. Radiotherapy and Oncology, 32(3):268 – 270, 1994. ISSN 0167-8140. doi:DOI:10.1016/0167-8140(94)90027-2. URL http://www.sciencedirect.com/science/article/pii/0167814094900272. → pages 10 [20] S. DiMaio. Modeling, Simulation and Planning of Needle Motion in Soft Tissues. PhD thesis, University of British Columbia, 2003. → pages 70, 72 [21] B. Drew, E. Jones, S. Reinsberg, A. Yung, S. Goldenberg, and P. Kozlowski. Device for sectioning prostatectomy specimens to facilitate comparison between histology and in vivo mri. Journal of Magnetic Resonance Imaging, 32(4):992–996, 2010. → pages 53 [22] K. Eichler, S. Hempel, J. Wilby, L. Myers, L. Bachmann, and J. Kleijnen. Diagnostic value of systematic biopsy methods in the investigation of prostate cancer: a systematic review. The Journal of urology, 175(5): 1605–1612, 2006. → pages 6 [23] M. Ferrant, A. Nabavi, B. Macq, F. Jolesz, R. Kikinis, and S. Warfield. Registration of 3-D intraoperative MR images of the brain using a finite-element biomechanical model. Medical Imaging, IEEE Transactions on, 20(12):1384–1397, 2002. ISSN 0278-0062. → pages 12, 50, 51 [24] P. Freeborough and N. Fox. Modeling brain deformations in alzheimer disease by fluid registration of serial 3d mr images. Journal of Computer Assisted Tomography, 22(5):838, 1998. → pages 2 [25] Y. Fung. Biomechanics: Mechanical Properties of Living Tissues. Springer-Verlag, Berlin, Germany, 1993. → pages 16 [26] O. Gonzalez and A. Stuart. A first course in continuum mechanics, volume 42. Cambridge Univ Pr, 2008. → pages 72 [27] P. Gould. Introduction to linear elasticity. Springer, 1994. → pages 72 [28] P. Hachenberger, L. Kettner, and K. Mehlhorn. Boolean operations on 3d selective nef complexes: Data structure, algorithms, optimized implementation and experiments. Computational Geometry, 38(1-2):64–99, 2007. → pages 24 [29] J. Hajnal. Medical image registration. 2001. → pages 3 63  [30] J. Hensel, C. M´enard, P. Chung, M. Milosevic, A. Kirilova, J. Moseley, M. Haider, and K. Brock. Development of multiorgan finite element-based prostate deformation model enabling registration of endorectal coil magnetic resonance imaging for radiotherapy planning. International Journal of Radiation Oncology* Biology* Physics, 68(5):1522–1528, 2007. → pages 12, 13 [31] Y. Hu, H. Ahmed, et al. MR to ultrasound registration for image-guided prostate interventions. Medical Image Analysis, 2010. ISSN 1361-8415. → pages 12, 58 [32] J. Kemper, R. Sinkus, J. Lorenzen, C. Nolte-Ernsting, A. Stork, and G. Adam. MR elastography of the prostate: Initial in-vivo application. In Rofo. Fortschritte auf dem Gebiet der Rontgenstrahlen und der bildgebenden Verfahren, volume 176(8), pages 1094–1099. Thieme, 2004. → pages 5, 21, 34, 35, 36, 44, 45 [33] T. Krouskop, T. Wheeler, F. Kallel, B. Garra, and T. Hall. Elastic moduli of breast and prostate tissues under compression. Ultrasonic imaging, 20(4): 260–274, 1998. → pages 5, 21 [34] K. Leissner and L. Tisell. The weight of the human prostate. Scandinavian Journal of Urology and Nephrology, 13(2):137–142, 1979. → pages 4, 37 [35] B. Li, G. Christensen, E. Hoffman, G. McLennan, and J. Reinhardt. Establishing a normative atlas of the human lung* 1:: Intersubject warping and registration of volumetric ct images. Academic radiology, 10(3): 255–265, 2003. → pages 2 [36] C. C. Ling, E. Yorke, and Z. Fuks. From imrt to igrt: Frontierland or neverland? Radiotherapy and Oncology, 78(2):119 – 122, 2006. ISSN 0167-8140. doi:DOI:10.1016/j.radonc.2005.12.005. URL http://www.sciencedirect.com/science/article/pii/S0167814005005694. → pages 9 [37] H. Livyatan, Z. Yaniv, and L. Joskowicz. Gradient-based 2-d/3-d rigid registration of fluoroscopic x-ray to ct. Medical Imaging, IEEE Transactions on, 22(11):1395 –1406, nov. 2003. ISSN 0278-0062. doi:10.1109/TMI.2003.819288. → pages 2 [38] S. S. Mahdavi, M. Moradi, X. Wen, W. J. Morris, and S. E. Salcudean. Evaluation of visualization of the prostate gland in vibro-elastography images. Medical Image Analysis, 15(4):589 – 600, 2011. ISSN 1361-8415. 64  doi:DOI:10.1016/j.media.2011.03.004. URL http://www.sciencedirect.com/science/article/pii/S1361841511000466.  Special section on IPMI 2009. → pages 36, 53 [39] W. Maurel, D. Thalmann, Y. Wu, and N. Thalmann. Biomechanical models for soft tissue simulation. Springer-Verlag, 1998. → pages 16 [40] J. Mitra, A. Oliver, R. Mart´ı, X. Llad´o, J. Vilanova, and F. Meriaudeau. Multimodal Prostate Registration using Thin-Plate Splines from Automatic Correspondences. In 2010 International Conference on Digital Image Computing: Techniques and Applications, pages 587–592. IEEE, 2010. → pages 11 [41] M. M¨uller and M. Gross. Interactive virtual materials. In Proceedings of Graphics Interface 2004, pages 239–246. Canadian Human-Computer Communications Society, 2004. → pages 31 [42] M. Murphy. An automatic six-degree-of-freedom image registration algorithm for image-guided frameless stereotaxic radiosurgery. Medical Physics, 24:857, 1997. → pages 2 [43] M. Nesme, M. Marchal, E. Promayon, M. Chabanas, Y. Payan, and F. Faure. Physically realistic interactive simulation for biological soft tissues. Recent Research Developments in Biomechanics, 2, 2005. → pages 31 [44] K. Noe and T. Sørensen. Solid Mesh Registration for Radiotherapy Treatment Planning. Biomedical Simulation, pages 59–70, 2010. → pages 12, 58 [45] M. Norberg, L. Egevad, L. Holmberg, P. Sparen, B. Norlen, and C. Busch. The sextant protocol for ultrasound-guided core biopsies of the prostate underestimates the presence of cancer*. Urology, 50(4):562–566, 1997. → pages 6 [46] K. Parker, S. Huang, R. Lerner, F. Lee Jr, D. Rubens, and D. Roach. Elastic and ultrasonic properties of the prostate. In Ultrasonics Symposium, 1993. Proceedings., IEEE 1993, pages 1035–1038. IEEE, 1993. → pages 21 [47] N. Patanjali, M. Keyes, W. Morris, M. Liu, R. Harrison, I. Spadinger, and V. Moravan. A comparison of post-implant us/ct image fusion and mri/ct image fusion for 125i prostate brachytherapy post implant dosimetry. Clinics in Occupational and Environmental Medicine, 8(2):124–124, 2009. → pages 47 65  [48] S. Phipps, T. Yang, F. Habib, R. Reuben, and S. McNeill. Measurement of tissue mechanical characteristics to distinguish between benign and malignant prostatic disease. Urology, 66(2):447–450, 2005. → pages 21 [49] C. Reynier, J. Troccaz, et al. MRI/TRUS data fusion for prostate brachytherapy. Preliminary results. Medical physics, 31:1568, 2004. → pages 1, 11 [50] A. Roche, G. Malandain, X. Pennec, and N. Ayache. The correlation ratio as a new similarity measure for multimodal image registration. Medical Image Computing and Computer-Assisted InterventationMICCAI98, pages 1115–1124, 1998. → pages 3 [51] W. Shao, R. Wu, C. Thng, K. Ling, H. Ho, C. Cheng, and W. Ng. Deformable registration for integration of MRI/MRSI information in TRUS-guided prostate biopsy. In Proceedings of SPIE, volume 5747, page 1263, 2005. → pages 9, 11 [52] H. Si. Tetgen, a quality tetrahedral mesh generator and three-dimensional delaunay triangulator. Technical report, Technical Report 9, WIAS, Berlin, Germany, 2004. → pages 33 [53] A. Singh, J. Kruecker, et al. Initial clinical experience with real-time transrectal ultrasonography-magnetic resonance imaging fusion-guided prostate biopsy. BJU international, 101(7):841–845, 2008. ISSN 1464-410X. → pages 1 [54] A. K. Singh, J. Kruecker, S. Xu, N. Glossop, P. Guion, K. Ullman, P. L. Choyke, and B. J. Wood. Initial clinical experience with real-time transrectal ultrasonography-magnetic resonance imaging fusion-guided prostate biopsy. BJU International, 101(7):841–845, 2008. ISSN 1464-410X. doi:10.1111/j.1464-410X.2007.07348.x. URL http://dx.doi.org/10.1111/j.1464-410X.2007.07348.x. → pages 6 [55] M. Steggerda, C. Schneider, M. van Herk, L. Zijp, L. Moonen, and H. van der Poel. The applicability of simultaneous trus-ct imaging for the evaluation of prostate seed implants. Medical physics, 32:2262, 2005. → pages 10 [56] R. Ten Haken, J. Forman, et al. Treatment planning issues related to prostate movement in response to differential filling of the rectum and bladder. International Journal of Radiation Oncology* Biology* Physics, 20(6): 1317–1324, 1991. ISSN 0360-3016. → pages 9 66  [57] M. Terris, E. Wallen, and T. Stamey. Comparison of mid-lobe versus lateral systematic sextant biopsies in the detection of prostate cancer. Urologia internationalis, 59(4):239–242, 1997. → pages 1 [58] G. Treece, R. Prager, A. Gee, and L. Berman. Surface interpolation from sparse cross sections using region correspondence. IEEE Transactions on Medical Imaging, 19(11):11061114, 2000. → pages 33 [59] G. Treece, R. Prager, A. Gee, and L. Berman. Surface interpolation from sparse cross sections using region correspondence. Medical Imaging, IEEE Transactions on, 19(11):1106–1114, 2000. → pages 33 [60] O. Ukimura, C. Magi-Galluzzi, and I. S. Gill. Real-time transrectal ultrasound guidance during laparoscopic radical prostatectomy: Impact on surgical margins. The Journal of Urology, 175(4):1304 – 1310, 2006. ISSN 0022-5347. doi:DOI:10.1016/S0022-5347(05)00688-9. URL http://www.sciencedirect.com/science/article/pii/S0022534705006889. → pages 10 [61] M. van Herk. Errors and margins in radiotherapy. Seminars in Radiation Oncology, 14(1):52 – 64, 2004. ISSN 1053-4296. doi:DOI:10.1053/j.semradonc.2003.10.003. URL http://www.sciencedirect.com/science/article/pii/S1053429603000845. High-Precision Radiation Therapy of Moving Targets. → pages 7 [62] H. Wang, L. Dong, J. O’Daniel, R. Mohan, A. Garden, K. Ang, D. Kuban, M. Bonnen, J. Chang, and R. Cheung. Validation of an accelerated’demons’ algorithm for deformable image registration in radiation therapy. Physics in Medicine and Biology, 50:2887, 2005. → pages 3 [63] R. Wu, K. Ling, W. Shao, and W. Ng. Registration of organ surface with intra-operative 3d ultrasound image using genetic algorithm. Medical Image Computing and Computer-Assisted Intervention-MICCAI 2003, pages 383–390, 2003. → pages 11 [64] S. Xu, J. Kruecker, B. Turkbey, N. Glossop, A. Singh, P. Choyke, P. Pinto, and B. Wood. Real-time mri-trus fusion for guidance of targeted prostate biopsies. Computer Aided Surgery, 13(5):255–264, 2008. → pages 11 [65] M. Zaider, M. Zelefsky, E. Lee, K. Zakian, H. Amols, J. Dyke, G. Cohen, Y. Hu, A. Endi, C. Chui, et al. Treatment planning for prostate implants using magnetic-resonance spectroscopy imaging1* 1. International Journal 67  of Radiation Oncology* Biology* Physics, 47(4):1085–1096, 2000. → pages 11 [66] M. Zhang, P. Nigwekar, B. Castaneda, K. Hoyt, J. Joseph, A. di Sant’Agnese, E. Messing, J. Strang, D. Rubens, and K. Parker. Quantitative characterization of viscoelastic properties of human prostate correlated with histology. Ultrasound in medicine & biology, 34(7): 1033–1042, 2008. → pages 5 [67] T. Zhang, N. Orton, T. Mackie, and B. Paliwal. Technical note: A novel boundary condition using contact elements for finite element based deformable image registration. Medical physics, 31:2412, 2004. → pages 12 [68] Y. Zhuang and J. Canny. Real-time simulation of physically realistic global deformation, 1999. → pages 72, 73  68  Appendices  69  Appendix A  Stress and Strain The linear elastic model used in this thesis is based on infinitesimal stress and strain. This appendix describes these terms in details. Consider a non-rigid 3- D elastic body at rest. The body is comprised of particles with position coordinates X = (x, y, z)T with respect to Coordinate frame C0 . This body deforms with the application of boundary conditions such that the positions of constituent particles change to Φ(X) where Φ is the deformation function. This deformation can be explained by a displacement vector U = (u(X), v(X), w(X))T which states the variation between the position of particles at rest and deformed state of the body. Please see Figure A.1 for an illustration.  A.1 Stress When a body is deformed by applying external forces, as a reaction, internal forces develop inside the elastic material. A body under these internal forces is said to be in the state of stress. A first order stress tensor t is described as the force f acting on an infinitesimal area ∆A oriented in a certain direction with a surface normal n. f ∆A→0 ∆A  t = lim  (A.1)  The state of stress of a particle inside the 3- D elastic body can be defined using three first order stress tensors [tx , ty , tz ] oriented in three orthogonal directions aligned with the coordinate frame [20]. The second-order stress tensor describing 70  Figure A.1: Illustration of deformation and displacement between an at rest and deformed body the state of stress of that particle is given by:   αxx αxy αxz   Σ(X) = αyx αyy αyz  αzx αzy αzz  (A.2)  where αxx , αxy and αxz are the components of tx . The components of ty and tz are defined similarly. Assuming the couple stress and body moment do not exist, it can be shown that Σ(X) is symmetric and therefore, the state of stress on a particle inside an elastic body:  α = [αxx  αyy  αzz  αyz  αzx  αxy ]T  (A.3)  where αxx , αyy and αzz are normal and αyz , αzx and αxy are shear components of stress.  A.2 Strain A normalized measure of the deformation produced as a result of application forces is called strain. Just as stress described in Section A.1, a strain is in general a tensor quantity and can be decomposed into normal and shear components. The 71  normal components correspond to stretch and compression and shear components correspond deformation due to sliding of planes against each other. The state of strain of a point on a continuum body can be described by these components on a set of three mutually independent directions[20]. A strain tensor can be thus written as:  ε = [εxx  εyy  εzz  εyz  εzx  εxy ]T  (A.4)  where εxx , εyy and εzz are normal and εyz , εzx and εxy are shear components. Expressing engineering strain tensor in terms of u [26, 27] while applying quadratic strain [68], we have:    δu  δx  1 2  δu δx  2  δv δx  2  δw δx  2    + + +      εxx 2 2   δv 1 δu 2 δ v δ w     +  εyy   δ y 2 δ y + δ y + δ y     2 2 2  ε   δ w 1 δ u δv δw   zz   + ε =   =  δz 2 δz + δz + δz   εyz      δv δw δu δu δv δv δw δw  ε   δ z + δ y + δ y δ z + δ y δ z + δ y δ z    zx    δu δw δu δu δv δv δw δw  + + + +  δz δx εxy δz δx δz δx δz δx    δu δv δu δu δv δv δw δw + + + + δy δx δy δx δy δx δy δx   δ  0 0 ΘTx 0 0   δx  0 ΘTy   0     δ   Θx u  0 δy 0    T   0 Θz    1 0 δ    =  Θy   0 0 δz  v + 2   0 ΘTz ΘTy   δ δ    w  0 δz δy ΘT 0 ΘT  Θz   z x δ 0 δδx T T δz Θy Θx 0 1 = Bu + AΘ 2 where Θx =  δu δv δw δx , δx, δx  T  (A.5)  and Θy and Θz are defined similarly [20]. The strain  tensor in the above form is called Green-Lagrange strain tensor. It is non-linear and rotation invariant. The linear form can be obtained if u is small enough for the  72  quadratic terms in Equation A.5 to be neglected:     εxx   εyy    ε   zz  ε = = εyz    ε   zx  εxy  δ  δx   0  0   0  δ δz  = Bu  0 δ δy  0 δ δz  0  0     0  δ  δz  δ  δy  δ δx    u   v w  (A.6)  This is called Cauchy’s infinitesimal strain tensor. This linear strain tensor enables us to use linear system of equations assuming that the deformations are small. However, unlike Green-Lagrange strain tensor, Cauchy’s strain tensor is not rotation invariant [68].  73  Appendix B  Prostate Registration Results We present the results of our deformable registration on prostate patient data in this appendix. The pink and green meshes represents the target (MR) and the source (TRUS) surfaces respectively. The meshes color-coded to the surface displacement (in meters) resulting from registration are also shown for every case.  B.1 Results of Pre-Operative Prostate Surfaces In this section we show images of unregistered and registered surfaces from the post-operative images. The landmarks from both the volumes are shown by the filled points. Red points represent the centers of segmented landmark from  MR  images. The black and blue filled markers represent center of landmark contour in unregistered and registered TRUS images respectively.  74  −3  x 10 11 10 0.015 9 0.01 8 0.005  7  0 Z  6  −0.005  75  5  −0.01  4  −0.015  3  −0.02  2 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.1: Registration result for pre-op case 1.  −0.01  0 Y  0.01  0.02  (c) Deformation field  1  −3  x 10 0.02  3  0.015 2.5 0.01  Z  0.005  2  0 1.5  76  −0.005 −0.01  1 −0.015 0.5  −0.02 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.2: Registration result for pre-op case 2.  −0.01  0 Y  0.01  (c) Deformation field  0.02  −3  x 10 0.02  3  0.015 2.5 0.01  Z  0.005  2  0 1.5  77  −0.005 −0.01  1  −0.015 0.5 −0.02  −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.3: Registration result for pre-op case 3.  −0.01  0 Y  0.01  (c) Deformation field  0.02  −3  x 10 0.02 4.5 0.015 4 0.01 3.5  Z  0.005  3  0  2.5  78  −0.005  2  −0.01  1.5  −0.015  1 0.5  −0.02 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.4: Registration result for pre-op case 4.  −0.01  0 Y  0.01  (c) Deformation field  0.02  B.2 Results of Post-Operative Prostate Surfaces In this section we show images of unregistered and registered surfaces from the pre-operative images. The urethra from both the volumes are shown by the fitted splines of the corresponding colors.  79  −3  x 10 12 0.015 10  0.01 0.005  8 Z  0 6  −0.005  80  −0.01  4  −0.015 2  −0.02 −0.025 −0.02  −0.01  0  0.01  Y  (a) Before Deformable registration  (b) After Deformable registration  Figure B.5: Registration result for post-op case 1.  (c) Deformation field  −3  x 10 6 0.02 0.015  5  0.01 4  Z  0.005 0  3  81  −0.005 −0.01  2  −0.015 1  −0.02 −0.025 −0.02  −0.01  0  0.01  Y  (a) Before Deformable registration  (b) After Deformable registration  Figure B.6: Registration result for post-op case 2.  (c) Deformation field  −3  x 10 8 0.015 7 0.01 6 0.005 5 Z  0 4  82  −0.005 3  −0.01  2  −0.015  1  −0.02 −0.015 −0.01 −0.005  0  0.005 0.01 0.015  Y  (a) Before Deformable registration  (b) After Deformable registration  Figure B.7: Registration result for post-op case 3.  (c) Deformation field  −3  x 10 5 0.02 4.5  0.015  4  0.01  3.5  0.005  3  Z  0  2.5  −0.005  83  −0.01  2  −0.015  1.5  −0.02  1  −0.025  0.5 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.8: Registration result for post-op case 4.  −0.01  0 Y  0.01  (c) Deformation field  0.02  −3  x 10 0.02  9  0.015 8 0.01 7 0.005 6 Z  0 5 −0.005  84  4 −0.01 3 −0.015 2 −0.02 −0.025 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.9: Registration result for post-op case 5.  1 −0.01  0 Y  0.01  (c) Deformation field  −3  x 10 5.5  0.02  5 4.5  0.01  4 3.5  Z  0  3  85  2.5  −0.01  2 1.5  −0.02  1 0.5  −0.03 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.10: Registration result for post-op case 6.  −0.01  0 Y  0.01  0.02  (c) Deformation field  −3  x 10 10 0.02 9 8  0.01  7 0 Z  6 5  86  −0.01 4 3  −0.02  2 −0.03  1 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.11: Registration result for post-op case 7.  −0.01  0 Y  0.01  0.02  (c) Deformation field  −3  x 10 12  0.015 0.01  10  0.005 8  Z  0 −0.005  6  87  −0.01 4 −0.015 −0.02  2  −0.025 −0.02  −0.01  0  0.01  Y  (a) Before Deformable registration  (b) After Deformable registration  Figure B.12: Registration result for post-op case 8.  (c) Deformation field  0.02  −3  x 10 10  0.02  9  0.015  8  0.01  7  0  6  −0.005  5  −0.01  4  −0.015  3  −0.02  2  Z  0.005  88  −0.025  1 −0.02  −0.01  0  0.01  Y  (a) Before Deformable registration  (b) After Deformable registration  Figure B.13: Registration result for post-op case 9.  (c) Deformation field  −3  x 10 7  0.01  6 0.005 5 0 Z  4 −0.005  89  3 −0.01  2  −0.015  1  −0.015  −0.01  −0.005  0  0.005  Y  (a) Before Deformable registration  (b) After Deformable registration  Figure B.14: Registration result for post-op case 10.  (c) Deformation field  0.01  −3  x 10 7  0.02 0.015  6  0.01 5  0.005 0 Z  4  −0.005  90  3  −0.01 −0.015  2  −0.02 1  −0.025 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.15: Registration result for post-op case 11.  −0.01  0 Y  0.01  (c) Deformation field  0.02  −3  x 10 0.015  3  0.01 2.5 0.005 2 Z  0 1.5  91  −0.005 −0.01  1  −0.015 0.5 −0.02 −0.015 −0.01 −0.005  (a) Before Deformable registration  (b) After Deformable registration  Figure B.16: Registration result for post-op case 12.  0 Y  0.005 0.01 0.015  (c) Deformation field  −3  x 10 4 0.02 3.5 0.015 3  0.01 0.005  2.5  Z  0 2  −0.005  92  −0.01  1.5  −0.015 1 −0.02 0.5  −0.025 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.17: Registration result for post-op case 13.  −0.01  0 Y  0.01  (c) Deformation field  0.02  −3  x 10 5.5 0.02  5 4.5  0.01  4 3.5  0 Z  3  93  2.5  −0.01  2 1.5  −0.02  1 −0.03 −0.02  (a) Before Deformable registration  (b) After Deformable registration  Figure B.18: Registration result for post-op case 14.  0.5 −0.01  0 Y  0.01  (c) Deformation field  Appendix C  Pathology and VE Registration Results We present the results of our deformable registration on prostate surfaces from pathology and elastography fransfer function images in this appendix. We obtained deformed 2- D images to study the results which are shown below.  94  Figure C.1: Case 1. Registered transfer function images. Colored contours are cancers contours marked on histo-pathology images overlaid on registered transfer function images.  95  Figure C.2: Case 2. Registered transfer function images. Colored contours are cancers contours marked on histo-pathology images overlaid on registered transfer function images.  96  Figure C.3: Case 3. Registered transfer function images. Colored contours are cancers contours marked on histo-pathology images overlaid on registered transfer function images.  97  Figure C.4: Case 4. Registered transfer function images. Colored contours are cancers contours marked on histo-pathology images overlaid on registered transfer function images.  98  Figure C.5: Case 5. Registered transfer function images. Colored contours are cancers contours marked on histo-pathology images overlaid on registered transfer function images.  99  Appendix D  Kidney Registration Results We present the results of our deformable registration on kidney patient data in this appendix. The pink and green meshes represents the target (CT) and the source (US) surfaces respectively.  100  (a) Before Deformable registration  (b) After Deformable registration  Figure D.1: Registration result for case 1.  (a) Before Deformable registration  (b) After Deformable registration  Figure D.2: Registration result for case 2.  101  (a) Before Deformable registration  (b) After Deformable registration  Figure D.3: Registration result for case 3.  (a) Before Deformable registration  (b) After Deformable registration  Figure D.4: Registration result for case 4.  102  (a) Before Deformable registration  (b) After Deformable registration  Figure D.5: Registration result for case 5.  103  Appendix E  Amount of Pressure and TRE We perform registration for all the cases except for  VE  and pathology registration  with the amount of pressure 1KPa. It can be argued that higher pressure can result in better surface fit. We claim that increased pressure can reduce the tortional effect observed in certain cases. To test this claim, we performed registration with higher pressure without PAA for five cases where rotation was involved in deformation and observed a greater  TRE  when higher pressure was used. Please see Table E.1 for  results. The results show that a good surface match does not translate into equally good internal registration.  104  Case  105  1 4 5 8 9  DSC 0.95 0.98 0.95 0.91 0.95  Pressure 1KPa Mean TRE (mm) 0.7±0.5 0.9±0.8 0.7±0.6 1.0±1.1 3.0±0.4  VD(%) 1.4 1.2 3.7 1.2 4.8  DSC 0.97 0.98 0.96 0.93 0.96  Pressure 10KPa Mean TRE (mm) 1.4±1.2 1.2±1.2 0.9±0.6 2.1±1.7 3.2±1.9  Table E.1: Results for first five cases with two levels of pressure.  VD(%) -4.3 -0.4 -2.2 -5.2 -1.5  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0072359/manifest

Comment

Related Items