In-vivo 3T and Ex-vivo 7T Diffusion Tensor Imaging of Prostate Cancer: Correlation With Histology by Carlos Felipe Uribe Mu˜noz B.Sc. Physics, Universidad de los Andes, Bogota, Colombia,2008 BA.Sc. Mechanical Engineering, Universidad de los Andes, Bogota, Colombia,2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Physics) The University Of British Columbia (Vancouver) April 2012 c Carlos Felipe Uribe Mu˜noz, 2012 To my Parents and my Sister. To my dear Lila. To my dog Homero. For the battle against Cancer. ii Abstract Diffusion Tensor Imaging has been successfully applied in prostate cancer diagnosis [27]. It has been well established that the water Apparent Diffusion Coefficient has a lower value in the prostate carcinomas when compared to normal prostatic tissue [2, 16, 27, 31, 41, 46, 56]. However, fractional anisotropy values in prostatic carcinoma have been reported to be higher [16], lower [31], and unchanged [56] when compared to the prostate’s normal peripheral zone. Preliminary data from a study involving diffusion tensor imaging measurements in prostate glands, in-vivo and ex-vivo following radical prostatectomy, is presented. Histology whole mount slides were registered to T2 weighted images and diffusion parametric maps using a mutual information voxel intensity registration algorithm using software developed in-house. Regions of interest which included the normal peripheral zone, the normal peripheral zone with enlarged glands, and tumours were taken into account for this study. The tumours were highlighted and graded with the Gleason score grading system by a specialized pathologist. Values of the apparent diffusion coefficient and the fractional anisotropy parameters were calculated. Monte-Carlo simulations of the behaviour of the fractional anisotropy with respect to the value of the apparent diffusion coefficient, the signal to noise ratio, and the b-value of the Stejskal-Tanner equation were performed. The results show lower values of the apparent diffusion coefficient for regions of tumours for the ex-vivo and in-vivo cases. Values of the fractional anisotropy in the prostate carcinomas are slightly higher ex-vivo than in-vivo, which may be explained by the dependence of the fractional anisotropy on partial volume effects and noise. These preliminary results show that the fractional anisotropy does not show significant differences between normal and cancerous tissue, strongly suggesting that it is not likely to contribute iii significantly to the diagnostic capabilities of diffusion tensor imaging in prostate cancer. iv Preface This research was approved by the University of British Columbia Research Ethics Board (H02-70400) and all participants gave signed consent prior to entering the study. v Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi List of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I Literature Review 2 The Prostate Gland . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1 Anatomy and Physiology . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Prostate Cancer Pathology . . . . . . . . . . . . . . . . . . . . . 11 2.3 Histology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3.1 15 3 Gleason Grading . . . . . . . . . . . . . . . . . . . . . . 1 Magnetic Resonance Imaging (MRI) and Diffusion Tensor Imaging (DTI) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.1 19 Nuclear Magnetic Resonance (NMR) . . . . . . . . . . . . . . . . vi 3.1.1 Generation of Macroscopic Magnetization . . . . . . . . . 23 3.1.2 The Rotating Reference Frame . . . . . . . . . . . . . . . 25 3.1.3 The RF Pulse . . . . . . . . . . . . . . . . . . . . . . . . 27 Relaxation Times . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2.1 Detecting the Induction Signal . . . . . . . . . . . . . . . 31 3.3 Spin-echo Experiments . . . . . . . . . . . . . . . . . . . . . . . 31 3.4 Spatial Localization . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.4.1 Slice Selection . . . . . . . . . . . . . . . . . . . . . . . 35 3.4.2 Frequency Encoding . . . . . . . . . . . . . . . . . . . . 36 3.4.3 Phase Encoding . . . . . . . . . . . . . . . . . . . . . . . 37 3.5 Signal to Noise Ratio and Partial Volume Effects . . . . . . . . . 39 3.6 Some Pulse Sequences . . . . . . . . . . . . . . . . . . . . . . . 41 3.6.1 Spin-Echo . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.6.2 Gradient Echo . . . . . . . . . . . . . . . . . . . . . . . 44 3.6.3 Rapid Acquisition with Relaxation Enhancement (RARE), 3.2 Fast Spin Echo (FSE), or Turbo Spin-echo (TSE) . . . . . 45 Echo Planar Imaging . . . . . . . . . . . . . . . . . . . . 46 Measuring Diffusion . . . . . . . . . . . . . . . . . . . . . . . . 47 3.7.1 Characterizing Diffusion . . . . . . . . . . . . . . . . . . 47 3.7.2 Stejskal-Tanner Experiment . . . . . . . . . . . . . . . . 49 3.7.3 Diffusion Anisotropy Indices . . . . . . . . . . . . . . . . 54 3.7.4 Diffusion Echo Planar Imaging (EPI) . . . . . . . . . . . 54 Image Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1 Types of Transformations . . . . . . . . . . . . . . . . . . . . . . 58 4.1.1 Dimensionality Transformations and Degrees of Freedom 58 4.1.2 Linear Transformations . . . . . . . . . . . . . . . . . . . 60 4.1.3 One to One Transformations . . . . . . . . . . . . . . . . 61 4.1.4 Nonlinear or Non-rigid Transformations . . . . . . . . . . 61 4.1.5 Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Image Registration Algorithms . . . . . . . . . . . . . . . . . . . 65 4.2.1 Landmark Based Registration . . . . . . . . . . . . . . . 66 4.2.2 Surface Matching . . . . . . . . . . . . . . . . . . . . . . 67 3.6.4 3.7 4 4.2 vii 5 4.2.3 Voxel Intensities . . . . . . . . . . . . . . . . . . . . . . 69 4.2.4 Optimization . . . . . . . . . . . . . . . . . . . . . . . . 79 Previous Work in Prostate DTI . . . . . . . . . . . . . . . . . . . . . 81 5.1 Bashar’s Study . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.2 Sato’s Group Study . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.3 Pickles’ Group Study . . . . . . . . . . . . . . . . . . . . . . . . 85 5.4 Manenti’s Group Study . . . . . . . . . . . . . . . . . . . . . . . 85 5.5 G¨urses’ Group Study . . . . . . . . . . . . . . . . . . . . . . . . 87 5.6 Xu’s Group Study . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.6.1 In-vivo MRI . . . . . . . . . . . . . . . . . . . . . . . . 89 5.6.2 Ex-vivo . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.6.3 MRI and Histology Coregistration . . . . . . . . . . . . . 90 Kozlowski’s Group Study . . . . . . . . . . . . . . . . . . . . . . 90 5.7 II 6 7 Methodology, Results, and Discussion Methods and Procedure . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.1 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . 95 6.2.1 Patients . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.2.2 In-vivo MRI . . . . . . . . . . . . . . . . . . . . . . . . 95 6.2.3 Ex-vivo MRI . . . . . . . . . . . . . . . . . . . . . . . . 97 6.2.4 Histology . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.2.5 Histology and MRI Coregistration . . . . . . . . . . . . . 101 6.2.6 Launching the Program . . . . . . . . . . . . . . . . . . . 102 Results and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 114 7.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 7.1.1 Validation of the Registration Software . . . . . . . . . . 114 7.1.2 Values of the Diffusion Parameters for the Different Types of Tissue in the Prostate . . . . . . . . . . . . . . . . . . 117 7.1.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 118 viii 7.2 7.3 Future Studies and Suggestions . . . . . . . . . . . . . . . . . . . 122 7.2.1 Digitizing the Whole Mount Sections (WMS) . . . . . . . 122 7.2.2 Registration Software . . . . . . . . . . . . . . . . . . . . 123 7.2.3 Correspondence of WMS and MRI . . . . . . . . . . . . . 124 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 III Appendix A Parameter Files for the Registration Algorithm . . . . . . . . . . . . 132 A.1 Affine Transformation Parameter Files (Initial Parameter Files) . . 132 A.1.1 Ex-vivo T2 Weighted (T 2 W) Images and WMS Registration 132 A.1.2 In-vivo T 2 W Images and WMS Registration . . . . . . . . 134 A.1.3 In-vivo DTI to in-vivo T 2 W Images Registration . . . . . . 136 A.2 Bsplines Transformation Parameter Files (Final Parameter Files) . 138 A.2.1 Ex-vivo T 2 W Images and WMS Registration . . . . . . . . 138 A.2.2 In-vivo T 2 W Images and WMS Registration . . . . . . . . 140 A.2.3 In-vivo DTI to in-vivo T 2 W Images Registration . . . . . . 142 A.3 Example of an Inversion Parameter File . . . . . . . . . . . . . . 144 ix List of Tables Table 1.1 Prostate cancer statistics. . . . . . . . . . . . . . . . . . . . . Table 3.1 Optimum parameters for the different types of contrast in spin- 1 echo images. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Table 5.1 Apparent diffusion coefficient values in Bashar’s experiments. . 83 Table 5.2 G¨urses’ results. . . . . . . . . . . . . . . . . . . . . . . . . . 88 Table 5.3 Kozlowski’s results. . . . . . . . . . . . . . . . . . . . . . . . 92 Table 6.1 TSE sequence details for the acquisition of the in-vivo T 2 W im- ages. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Table 6.2 Diffusion EPI details for the acquisition of the in-vivo DTI data. 96 Table 6.3 RARE sequence details for the acquisition of the ex-vivo T 2 W images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table 6.4 Spin-echo sequence details for the acquisition of the ex-vivo DTI Table 6.5 97 images. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Scaling factors required to load and save the different maps of the DTI reconstruction function dtiUI . . . . . . . . . . . . . . 110 Table 7.1 Percentage differences between the mean Apparent Diffusion Coefficient (ADC) and Fractional Anisotropy (FA) values. . . . 116 x List of Figures Figure 2.1 Location of the prostate gland within the body. . . . . . . . . 7 Figure 2.2 Sagittal view of the prostate gland. . . . . . . . . . . . . . . . 10 Figure 2.3 Zonal anatomy of the prostate gland. . . . . . . . . . . . . . . 10 Figure 2.4 Digital rectal examination. . . . . . . . . . . . . . . . . . . . 12 Figure 2.5 Functions of the prostate specific antigen. . . . . . . . . . . . 13 Figure 2.6 Patient’s screening procedure. . . . . . . . . . . . . . . . . . 14 Figure 2.7 Prostate needle biopsy. . . . . . . . . . . . . . . . . . . . . . 15 Figure 2.8 Gleason grades in prostate adenocarcinomas. . . . . . . . . . 17 Figure 3.1 NMR experiments, basis of MRI signal collection. . . . . . . 20 Figure 3.2 Precession in an external magnetic field. . . . . . . . . . . . . 22 Figure 3.3 Energy states in Zeeman effect. . . . . . . . . . . . . . . . . 24 Figure 3.4 Magnetization vector in a laboratory reference frame. . . . . . 26 Figure 3.5 Path of a walk on two reference frames. . . . . . . . . . . . . 26 Figure 3.6 Magnetization being tipped from the z axis. . . . . . . . . . . 27 Figure 3.7 View in the rotating reference frame. . . . . . . . . . . . . . . 28 Figure 3.8 Plot of the longitudinal magnetization as a function of time. . 30 Figure 3.9 FID signal detected in the coil. . . . . . . . . . . . . . . . . . 32 Figure 3.10 Spin behavior in the Carr and Purcell sequence. . . . . . . . . 33 Figure 3.11 Modification of the CP sequence by Meiboom and Gill. . . . . 34 Figure 3.12 Spin-echo diagram. . . . . . . . . . . . . . . . . . . . . . . . 34 Figure 3.13 Phase encoding example. . . . . . . . . . . . . . . . . . . . . 38 Figure 3.14 Spin-echo pulse diagram. . . . . . . . . . . . . . . . . . . . . 42 Figure 3.15 Spin-echo alternate pulse diagram. . . . . . . . . . . . . . . . 43 Figure 3.16 Gradient-echo pulse sequence. . . . . . . . . . . . . . . . . . 45 xi Figure 3.17 RARE pulse sequence. . . . . . . . . . . . . . . . . . . . . . 46 Figure 3.18 EPI pulse sequence. . . . . . . . . . . . . . . . . . . . . . . . 47 Figure 3.19 Diagram of isotropic and anisotropic diffusion. . . . . . . . . 48 Figure 3.20 Diffusion spin-echo sequence. . . . . . . . . . . . . . . . . . 50 Figure 3.21 Stejskal-Tanner sequence. . . . . . . . . . . . . . . . . . . . 51 Figure 3.22 Diffusion EPI pulse sequence. . . . . . . . . . . . . . . . . . 55 Figure 4.1 Sketch of the image registration process. . . . . . . . . . . . . 57 Figure 4.2 Different types of transformations. . . . . . . . . . . . . . . . 62 Figure 4.3 Misalignment in registered images. . . . . . . . . . . . . . . 75 Figure 4.4 Two-dimensional histograms. . . . . . . . . . . . . . . . . . 77 Figure 5.1 Bashar’s results. . . . . . . . . . . . . . . . . . . . . . . . . 83 Figure 5.2 Sato’s results. . . . . . . . . . . . . . . . . . . . . . . . . . . 84 Figure 5.3 Pickles’ results. . . . . . . . . . . . . . . . . . . . . . . . . . 86 Figure 5.4 Manenti’s results. . . . . . . . . . . . . . . . . . . . . . . . . 87 Figure 5.5 Xu’s results. . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Figure 6.1 Experimental set up of the specimen for the collection of the ex-vivo data. . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Figure 6.2 Ex-vivo pulse sequences used to acquire the data. . . . . . . . 99 Figure 6.3 Changes in the slicer device. . . . . . . . . . . . . . . . . . . 101 Figure 6.4 Sectioned specimen. . . . . . . . . . . . . . . . . . . . . . . 102 Figure 6.5 in ex reg function which launches the registration program. . 103 Figure 6.6 Main window for the ex-vivo to WMS registration procedure. . 103 Figure 6.7 slreg loadData window. . . . . . . . . . . . . . . . . . . . . 104 Figure 6.8 slreg pairData windows with two paired slices. . . . . . . . . 105 Figure 6.9 slreg segment specimen segmentation. . . . . . . . . . . . . . 106 Figure 6.10 slreg register window. . . . . . . . . . . . . . . . . . . . . . 106 Figure 6.11 Information displayed at the end of the ex-vivo registration process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 Figure 6.12 In-vivo registration procedure main window. . . . . . . . . . 108 Figure 6.13 dtiUI main window. . . . . . . . . . . . . . . . . . . . . . . . 109 Figure 6.14 Final window of the in-vivo registration process. . . . . . . . 111 xii Figure 6.15 Registration of in-vivo DTI to histology. . . . . . . . . . . . . 112 Figure 7.1 Ex-vivo registered slice and overlay. . . . . . . . . . . . . . . 115 Figure 7.2 A histology slice with outlined regions of interest. Figure 7.3 Median values of distance between corresponding landmarks . . . . . . 115 in the registered and fixed images for both ex-vivo and in-vivo. 116 Figure 7.4 ADC and FA values for different types of tissue in the prostate gland. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 Figure 7.5 Magnified images of healthy peripheral zone and prostate cancer.119 Figure 7.6 Monte-Carlo simulations of FA and its dependence on noise. . 120 Figure 7.7 Monte-Carlo simulations of the behavior of FA as a function of ADC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Figure 7.8 A new WMS digitization scanner. . . . . . . . . . . . . . . . . 122 Figure 7.9 Leica scn400 slide scanner specifications. . . . . . . . . . . . 123 xiii List of Acronyms AFMS Anterior fibromuscular stroma ADC Apparent Diffusion Coefficient BPH Benign Prostatic Hyperplasia BUN Blood Urea Nitrogen CCS Canadian Cancer Society CG Central Gland CZ Central Zone CT Computed Tomography CC Correlation Coefficient DSC Dice Similarity Coefficient DAI Diffusion Anisotropy Indices DWI Diffusion Weighted Images DTI Diffusion Tensor Imaging DRE Digital Rectal Examination DCE Dynamic Contrast Enhancement EPI Echo Planar Imaging xiv TE Echo Time FSE Fast Spin Echo FOV Field Of View FA Fractional Anisotropy FID Free Induction Decay FFDS Freeform Deformations GS 33 Gleason Score of 3 + 3 GS 34 Gleason Score of 3 + 4 ICP Iterative Closest Point MRI Magnetic Resonance Imaging NPZ Normal Peripheral Zone NPZEG Normal Peripheral Zone With Enlarged Glands NMR Nuclear Magnetic Resonance PIU Partition Intensity Uniformity PET Positron Emission Tomography PDF Probability Density Function PCA Prostate Cancer PSA Prostate Specific Antigen PDW Proton Density Weighted RARE Rapid Acquisition with Relaxation Enhancement RTI Rate of Transmission of Information RIU Ratio Image Uniformity xv ROI Regions of Interest TR Repetition Time SNR Signal to Noise Ratio SAD Sum of Absolute Differences SSD Sum of Squared Intensity Differences T1W T1 Weighted T2W T2 Weighted TZ Transitional Zone TSE Turbo Spin-echo WMS Whole Mount Sections xvi Acknowledgements I want to thank my supervisor Dr. Piotr Kozlowski for giving me the opportunity to do research on this field which I am truly passionate about. Without his support and patience this project would not have been possible. Thanks for pushing me to take control of my project and helping to improve my skills every day. I am grateful for always encouraging me to develop my career, and motivating me to attend conferences and symposiums. Special thanks to Dr. Stefan Reinsberg for spending time reviewing and commenting on my work in order to make it better as time went by. The discussions on how to explain my results and the help with the simulations were essential for the project’s success. Thanks to the staff from the Vancouver General Hospital, our pathologist Dr. Edward Jones (Ted), the histology technician Margaret Luk, and the resident Dr. Brian Skinnider. All the care required to successfully store, transport, slice, and annotate the slices could not have been better. Dr. Larry Goldenberg, head of the department of Urology, and Dr. Silvia D. Chang, part of the 3T MRI Facility at UBC hospital, thank you very much for your dedication in recruiting patients and the high quality data obtained from their scans. I would also like to thank Dr. Stefan Klein and Kyzyl Herzog for their assistance developing the image registration software. Although I never got to meet them, I express my gratitude to the patients who xvii decided to enrol in this research project. We are working hard to keep future generations away from the difficult times caused by cancer. All this would be impossible without that altruistic contribution. Thanks to the Canadian Institutes of Health Research for believing and providing the funds for this project. My special gratitude to the Government of Canada and the University of British Columbia for giving me the opportunity to be here and work with this wonderful group of people. To my parents and my sister, if I have managed to be here is due to all the effort that you have made throughout my life and all the lessons you have taught me. I miss you all around here, and I am pretty sure that you miss me too, but Homero is definitely capable of protecting you all and sharing its affection keeping us together no matter how far we are; he is a very good dog, not to say the best one that has ever existed. Lila, this paragraph is probably the hardest one that I have written in this thesis. There are just not enough words to express my feelings and appreciation for all those wonderful moments, company, and support that you have given me. Thanks for showing your interest in reading this project and helping me put it together. You have always been inspiring me on the difficult moments and sharing happiness on the good ones. Libia and Lucy, is sad that you left this world and are not physically sharing this milestone with me. I am however grateful for all the care you gave me through all the moments of my life. I am confident that you will think that everything has been worth it. This will not just end here, it still requires many more years of work and I won’t stop doing it while I’m still alive. I just wanted to let you know that you will always have a very special place in my heart on all the things I get to do. Thanks to my friends and colleagues from the 7T facility in the MRI Research Centre. Thanks to my friends back home, friends from Canada, and friends from all around the world for being by my side in this new achievement. xviii Chapter 1 Introduction According to the global burden of disease published by the World Health Organization [40], prostate cancer is a type of non-communicable disease which is expected to cause one of the highest number of deaths in men. In particular, a projected pessimistic scenario for the year 2015 shows that after trachea, bronchus, and lung, prostate cancer will be the cause of death for most men in the Americas (north, central, and south). The Canadian Cancer Society (CCS), whose statistics are shown in table 1.1, states that 27% of all new cancer cases in men and 10% of all cancer deaths were caused by prostate cancer in just 2011 [48]. Movements like the Canadian Movemeber, in which thousands of men in Canada grow moustaches to collect funds during the month of November, has been created in order to raise the community awareness for men’s health; especially regarding prostate cancer. Cases 25,500 Incidence rate per 100,000 Incidence Rank Deaths Death Rate per 100,000 5-year survival (2004-2006) 122 1 4,100 21 96% Table 1.1: Prostate Cancer Statistics. Data from the Canadian Cancer Society 2011. [48] 1 This project begins with a literature review starting with chapter 2, showing a description of the prostate’s gland anatomy and physiology, and it focuses on the way in which prostate cancer is diagnosed today. Additionally, the zone characterization of the gland proposed by McNeal is shown in detail, and its typical diseases are explained. Following this, methods of identifying prostate cancer as a pathology are presented, and an emphasis in the prostate specific antigen in blood serum as well as the digital rectal examination is made. The chapter closes with an explanation of the histological analysis of the cancer, where the Gleason score grading system for the different types of prostate carcinomas is detailed, and is the grading scheme used throughout this study. Chapter 3 is about magnetic resonance imaging that provides a way of acquiring images in a non-invasive or at least a less invasive way that how is usually done in current days. A description of the theory behind nuclear magnetic resonance, and the way in which a macroscopic magnetization is generated is presented. Concepts such as the rotating reference frame, the RF pulse, as well as how natural relaxation times change the behaviour of magnetization are described. Moreover, signal detection, the way in which noise affects the signal, and spatial localization, including an explanation of several typical pulse sequences, are then shown. Finally, the chapter closes with the theory behind the measurement of diffusion properties with magnetic resonance, and two pulse sequences for this are also explained. In order to be able to analyze data, the project involved the development of a registration software. Chapter 4 presents the theory behind it starting with an explanation of different types of transformations. Several image registration algorithms are mentioned, emphasizing the one that involves voxel intensity values. Different metrics can be used to compare the images being registered, but special attention is given to the mutual information technique. In practice, the non-invasive localization and grading of prostate cancer remains a challenge. Different researchers have been trying to develop new techniques that lead to the improvement of prostate cancer diagnosis. Chapter 5 dis2 cusses several studies from the first decade of this century, in which diffusion techniques and parameters are used to grade and determine the existence of prostate carcinomas. The studies agree that the apparent diffusion coefficient is reduced in carcinomas as compared to healthy tissue, but fractional anisotropy still needs to be validated; all type of values have been presented by different research groups. Finally, the project concludes with a description of the methods, results, and discussion of the current study performed at the UBC MRI Research Centre. The purpose of the study was to determine if the fractional anisotropy parameter provides useful information for the localization and grading of prostate cancer. The description of the software developed for the image registration is detailed, and the way images were collected for in-vivo and ex-vivo cases is presented. A discussion of the results has been made suggesting that fractional anisotropy does not add any useful information to the diagnosis of prostate cancer. 3 Part I Literature Review 4 Chapter 2 The Prostate Gland 2.1 Anatomy and Physiology The prostate gland is only present in males and it is part of the reproductive and urinary system. It is developed during the third month of gestation influenced by the mesenchyme, which is a type of connective tissue. The gland remains almost identical until puberty, with its base measuring about 4 cm transversely, 2 cm anterior to posterior, and 3 cm in its vertical diameters and it weights around 8 grams [49]. During puberty, it follows some morphological changes and it will enlarge until it reaches a weight of approximately 20 grams at 25 to 30 years of age [19]. It usually and inevitably enlarges as a man ages with the development of a Benign Prostatic Hyperplasia (BPH) getting to weight 40 g and even 150 g in some men older than 50. People usually describe the gland without BPH as having the shape of a croissant while an enlarged gland would have the shape of a doughnut.[49] The human prostate gland has a pyramidal shape formed of fibromuscular tissue, and made up of a mixture of glandular and non-glandular components kept together inside a capsule. The capsule is not really a fibrous one but rather enclosed by visceral fascia 1 that contains neurovascular tissue. By slicing the gland in well selected planes, the anatomical features can be clearly observed. The gland 1A fascia is a fibrous connective tissue that wraps around muscle bundles to hold them in place[45] 5 can be seen as being convex towards the anterior due basically to non-glandular tissue. The prostatic urethra is a primary reference point as four glandular regions arise from different segments of it. [11, 45] The muscular tissue of the gland is mostly smooth muscle. Anterior to it, skeletal muscle continuous with the external urethral sphincter, and supplied by the pudendal nerve, compresses the urethra pulling the urethral crest to the back and the prostatic sinuses to the front. This causes the urethra to be dilated and glandular contents, simultaneously released into it, contains between 3 to 5 ml of seminal fluid before the ejaculation [11] making around 80% of the fluid released during ejaculation. These fluids are drained to the urethra by ducts originating from the tiny glands and are the ones that give the signal to the sperm cells to begin to “swim”. The seminal vesicles located above the gland store sperm and ejaculatory fluids and are considered an extension of the prostate gland. The gland’s muscle fibre that squeeze the urethra helps control the urinary flow.[14] Due to its pyramidal shape, the prostate has a base located superiorly and an apex in the inferior. Figure 2.1 shows the location of the prostate gland in a human male. It can be found at a low level in the lesser pelvis, behind the symphysis pubis and pubic arch, and anterior to the rectourethralis and rectal ampulla. When men visit their urologist for the yearly prostate exam, the gland is palpated through the rectal ampulla in order to identify any malignant growth of it. [11] The gland’s base is located at the bladder neck while the apex reaches the urogenital diaphragm. It is separated from the rectum by a small layer of connective tissue called Denonvilliers’ fascia. The transverse measure of the base of the gland is about 2 cm, the anterior to posterior length is of approximately 2 cm and close to 3 cm in its vertical length [49]. The sphincter at the base of the bladder closes and prevents urine from entering the urethra at the moment of ejaculation. [14] During the last century, studies of animal prostate glands showed lobes and suggested that human glands should be described by different lobes [19]. The human prostate gland has five lobes only in the first 20 weeks of gestation and then only three lobes can be recognized; two lateral ones and one median [49]. It was 6 (a) Coronal Slice (b) Sagittal Slice Figure 2.1: Coronal and sagittal views of the location of the prostate gland within a man. ( c Elsevier taken from Coakley and Hricak [8] with permission.) 7 McNeal who established the accepted concept of zones that is used today [1, 33– 35]. The gland is divided into three glandular zones and one non-glandular tissue as it can be seen in figure 2.3. • Peripheral Zone: It covers around 70% of the gland’s volume. It involves all the glandular tissue at the apex and the one posteriorly close to the capsule. Its shape is similar to a cup and it encloses the other two glandular zones and the preprostatic urethra, except anteriorly. Most carcinomas occur here as well as chronic prostatitis 2 and post inflammatory atrophy. [19, 49] • Central Zone: It covers about 25% of the prostate’s volume. This zone has a cone-shape that surrounds the ejaculatory ducts posterior to the preprostatic urethra and its apex is located at the verumontanum 3 . This zone is rarely involved in any prostatic disease. [19, 49] • Transition zone: It covers about 5% of the gland’s volume. It involves two equal parts of glandular tissue lateral to the distal urethra and is located proximal to the apex of the central zone and ejaculatory ducts. The ducts enter the urethra below the preprostatic sphincter and above the ducts from the peripheral zone. BPH usually develops here and though less common, cases of adenocarcinoma have also been detected. BPH can grow to form a bulk in the gland. It begins as micronodules that grow and combine together forming macronodules around the inferior preprostatic urethra just above the verunontanum.[19] BPH has a high incidence of 50% in 50 year old men and no evidence has been found that it leads to cancer. The enlarged gland squeezes the urethra tighter than in normal situations obstructing the urine flow from the bladder. The muscles from the bladder then become thicker and stronger in order to compensate for this effect by contracting harder to move the urine through the obstacle. If the prostate enlarges so much that 2 Prostatitis is an infection of the gland that is usually caused by no apparent reasons. It may be present as a chronic and acute infections. It’s symptoms are urgency to urinate, burning sensation, and difficulty to pass urine. Antibiotics are used to treat prostatitis. [14] 3 Location where the seminal ducts enter the urethra. 8 the urethra becomes very narrow, the bladder might not be able to empty completely and residual urine will remain in it increasing the frequency of urination needs.[14, 49] • Anterior fibromuscular stroma (AFMS): The anterior external surface of the prostate gland has a convexity formed by this stromal tissue. The region close to the apex of the gland contains striated muscle. On the other hand, the half portion close to the base, is rich in smooth muscle. Figure 2.2 shows a sagittal slice of the prostate in which the urethra can be seen having an abrupt angulation, between prostate apex and bladder neck, dividing it in proximal and distal segments. Though a typical angle of deviation of about 35◦ is usually seen, it can be bigger in men with nodular hyperplasia. The ejaculatory ducts join the distal urethral segment while around 95% of the glandular prostate ducts end at this point as well making this segment the only one involved in the ejaculation process. Voluntary sphincter functions and involuntary sphincter functions are controlled by the distal and proximal portions of the AFMS respectively.[11, 19] The zonal anatomy of the prostate gland is important due to the fact that most carcinomas originate in the peripheral zone; zonal anatomy can be easily identified with Magnetic Resonance Imaging (MRI). High signal intensity on T2 Weighted (T 2 W) images is obtained for the Normal Peripheral Zone (NPZ) as well as for fluids within the seminal vesicles. The Central Zone (CZ) and Transitional Zone (TZ) are commonly grouped into what is known as the Central Gland (CG) and they show a relatively low signal on the T 2 W images. The verumontanum shows a high signal within the CG. It has to be kept in mind that the zones change as the person ages. The CG atrophies and the TZ enlarges. A low signal band, the surgical pseudo capsule, can be seen in the limit of the atrophied TZ and compressed peripheral zone on T2-weighted images. [49] Details regarding MRI are given in chapter 3. 9 Figure 2.2: Sagittal view of the prostate gland. The angulation of the urethra divides the prostate in proximal and distal segments. ( c Elsevier taken from Foster and Bostwick [11] with permission). Figure 2.3: Zonal anatomy of the prostate gland. ( c Nature Publishing Group taken from Marzo et al. [32] with permission.) 10 2.2 Prostate Cancer Pathology Out of all prostatic neoplasms, adenocarcinomas 4 account for 95% of them. Most of the times they are silent diseases, but in some cases they present obstructive symptoms similar to the ones caused by BPH. Cases of patients being diagnosed with an advanced prostate cancer with metastases without any symptoms related to the prostate area are common. Early detection programs have been developed during the last two decades and techniques involving screening with Prostate Specific Antigen (PSA) and the Digital Rectal Examination (DRE) are the most typically used in present days.[19] DRE consists of a palpation of the prostate gland by using a lubricated gloved finger in the patient’s rectum and feeling towards the front. Figure 2.4 illustrates the process in which the physician might be able to feel a hard lump caused by the cancer. This exam usually lasts just seconds and, though painless, it does cause an unpleasant need to urinate or defecate. Several positions are used to carry out the exam including the patient lying on his side, on his back, on his stomach, or standing. BPH feels like a soft rubber and it is usually symmetric, while the cancer lump is felt as hard plastic or wood.[14] Blood tests are also used to detect the cancer in the prostate. Hemoglobin, white blood cell counts, and platelet count are routine exams. Other exams to measure concentration of salts in blood, and kidney functioning like the Blood Urea Nitrogen (BUN) and serum creatinine, help to determine any blockage of the ureters by a prostate growth. However, the most used blood test is the one that measures the amount of PSA, which is produced by the prostate gland and is present in the blood. The PSA test indicates whether the disease might exist, but it does not diagnose it.[14] PSA was detected in seminal fluid during the early 1970’s and isolated from the prostatic tissue by the end of that decade. There has been well documented correlation between the higher risk of undetected prostate cancer and its stage with 4 Cancer that originates within the tiny glands of the prostate. [14] 11 Figure 2.4: Illustration of the digital rectal examination. (This image is in the public domain and can be freely reused. Taken from NCI [38].) increased serum levels of PSA. Its primary biological function is to cause semen liquefaction that leads to fertilization. Normal prostatic epithelial cells highly express this antigen. Some evidence has lead to think that PSA might accelerate carcinogenesis in the prostate. It can also stimulate the cell detachment, and then cause tumour cell invasion and metastases. Treated patients have shown a decrease in PSA levels in correlation with time to pain progression and survival. Figure 2.5 shows the functions of PSA in both the healthy and cancerous tissue. [7] The normal value of PSA circulating in the blood is considered to be below 4.0ng/mL. In order to better interpret the PSA results, researchers have been proposing an age related scale in which older men are allowed higher levels of the antigen in the blood compared to younger men. [14] In general to determine the type of follow up procedure that the patient requires, physicians follow the diagram shown in figure 2.6. The patient follows a DRE and PSA blood test. If the DRE is negative and the levels of PSA are normal then no action is needed and the screening is repeated every one or two years. If on the 12 Figure 2.5: Functions of PSA. It helps in liquefaction and fertilization in the healthy prostate while causing cell detachment, proliferation and tumour cell invasion and metastasis when carcinogenesis is present. ( c Nature Publishing Group taken from Bok and Small [7] with permission.) other hand, the PSA levels are high even if the DRE is negative, or if the DRE is mildly to highly suspicious regardless of the PSA levels then the patient will be followed up by the physician. 13 Figure 2.6: Patient’s screening procedure. 2.3 Histology Needle biopsy samples and radical prostatectomy are used to determine the histological grading of the cancer. In the former case, a sample of the lump is obtained for microscopic analysis. Using ultrasound to guide the needle, the physician inserts a needle into the area where the cancer is suspected to be. As illustrated in figure 2.7, a very thin needle passes through the wall of the rectum without damaging it into the prostate, and a small portion of the tissue is removed. The patient might feel a sharp jab but the pain is gone by the end of the procedure.[14] Some patients require radical prostatectomy, in which the gland and some surrounding tissues are completely removed, in order to provide a good chance of cure. The physician makes the decision based on [14]: • The tumour must be confined to the prostate. If that is not the case, chances of the patient being cured are small and the risks of the prostatectomy are not worth it. • The levels of PSA must not be very high. 14 Figure 2.7: Illustration of the prostate needle biopsy guided with ultrasound. ( c 2005 Terese Winslow, U.S. Govt. has certain rights. Taken from Winslow [54] with permission.) • The patient should be in good physical condition in order to resist the anaesthetic and the surgery. Age is not necessarily a factor involved unless the patient’s life expectancy is at least 10 years. 2.3.1 Gleason Grading Grading methods in prostate cancer are used just for adenocarcinomas, which represent the majority of prostate malignancies. It helps the pathologist to predict the tumours biological behavior. Though several methods have been proposed to determine the morphological features of the adenocarcinomas, only the Gleason Score grading method, proposed by Gleason in 1966 [13], will be described as it is the one relevant for the current project. [11] Gleason participated in a study performed by the Veterans Administration Cooperative Urological Research Group (VACURG) that analyzed histomorphologic appearances of more than 4000 prostate specimens. A low magnification was used 15 to categorize the histological patterns of the glands and the pattern of growth of the tumour in the stroma. Figure 2.8 shows the pattern for the different Gleason grades.[11] The different Gleason grades can be described as follows: [11] • Grade 1: Architecture of the gland and tumours is almost not distorted. The neoplastic glands are uniform in shape, round, and single and clearly delineated from the fibrovascular stroma. To differentiate from some type of BPH , grade 1 tumour cells contain enlarged nucleoli. • Grade 2: This type of tumours are mild but a separation of the tumour glands by the stroma can be observed. In comparison to grade 1, the variation of the size and the shape of the glands is higher. The tumour is circumscribed, but neoplastic gland separation can be observed at the periphery. • Grade 3: Carcinomas described by this grade are divided in three histological forms depending on the architectural dispersion of the neoplastic glands. Compared to grade 2 tumours, grades 3A and 3B show more variation in size, shape, and separation between the individual glands which is usually greater than one average gland diameter. An extension of glands into the surrounding stroma is also present. The nuclei of these cells are bigger than the ones from grades 1 and 2. Grades 3A and 3B tumours have different average sizes of the glands. The former case shows medium to large sizes while for the latter one the glands are small. Grade 3C show a mix of masses and cords with cribriform tumours. They don’t show invasive edges and are rather rounded and smooth. • Grade 4: This grade of cancer can take a cribriform pattern in which solid masses of tumour cells are formed. It is associated with big volume tumours that contain ragged invading edges (grade 3C show smooth edges). • Grade 5: No glandular differentiation can be seen, and the tumours are formed by solid masses of cells. It also includes anaplastic tumours that grow in nonglandular tumour masses. 16 (a) Grade 1 tumour at 50X magnification. ( c Elsevier taken from Foster and Bostwick [11] with permission.) (b) Grade 2 tumour at 50X magnification. ( c Elsevier taken from Foster and Bostwick [11] with permission.) (c) Grade 3 tumour at 100X magnification. ( c Elsevier taken from Foster and Bostwick [11] with permission.) (d) Grade 4 tumour at 50X magnification. ( c Elsevier taken from Foster and Bostwick [11] with permission.) (e) Grade 5 tumour at 100X magnification. ( c Elsevier taken from Foster and Bostwick [11] with permission.) (f) Diagram of the Gleason score grading system. (Egevad et al. [10] by permission.) Figure 2.8: Magnified images and diagram of the different Gleason grades in prostate adenocarcinomas. 17 Usually, prostate cancers present more than one Gleason grade in their histological patterns. A primary grade is then assigned to the predominant pattern and a second grade is assigned to the second most dominant one, with the restriction that it must include at least 5% of the cells of this type. The grades are scored using integer numbers ranging from 1 to 5, according to the patterns shown in figure 2.8. Grades are then added providing the histological score. Scores of 6 and 7 are the most common diagnosed tumours, and it is known that patients with score 7 have a higher risk of suffering from lymph node metastases.[11] 18 Chapter 3 and Diffusion Tensor Imaging (DTI) MRI 3.1 Nuclear Magnetic Resonance (NMR) In 1938, Rabi and his group were able to measure the nuclear magnetic moment by using atomic beams [43] and won the Nobel prize in 1944. The work gave origin to nuclear magnetic resonance experiments. In 1952, the Nobel prize in Physics was shared by F. Bloch, from Stanford University, and E.M. Purcell from Harvard university. The former (figure 3.1a), used an induction method [5, 6] in which two orthogonal coils were placed close to a sample of water. A static magnetic field was present and its direction was perpendicular to the axes of the coils. An alternating current was passed through one of the coils and the induction signal was collected in the other one. Purcell’s method (figure 3.1b), called the absorption method, consisted in placing a sample in a strong magnetic field perpendicular to the axes of a coil. Power loss was noticed when a certain frequency of the alternating current was set. [42] This NMR experiments provide the basis for the detection of the signal in MRI. Today, induction methods are preferred than absorption ones.[55] To be able to understand how the nuclear magnetic moments interact with the external magnetic fields, a closer look to the atoms has to be seen first. Protons have an intrinsic property called spin and it has a quantum mechanical origin. Charged 19 (a) NMR induction experiment performed by Bloch and his group. ( c Xiang taken from Xiang [55] with permission.) (b) NMR absorption experiment performed by Purcell and his group. ( c Xiang taken from Xiang [55] with permission.) Figure 3.1: NMR experiments, basis of MRI signal collection. 20 particles generate a magnetic moment µ which, for the case of a proton, is given by equation 3.1. The gyro-magnetic ratio γ is a constant of the nucleus defined as the ratio of the magnetic moment and the total angular momentum J. In order to have a non-zero magnetic moment, the total angular momentum of the nuclei must not be null. This occurs in atoms with an odd number of protons or neutrons. The Hydrogen atom 1 H is of particular interest as it is present in high quantities in the human body, and are the used nuclei to obtain the images, as it will be described later in this chapter. It contains one proton and from now on there will be no distinction made between a proton and the hydrogen atom.[26, 55] µ = γJ (3.1) When a proton is located in an external magnetic field B, the magnetic moment µ tends to align with the field. However, as angular momentum is present, a torque N shown in equation 3.2 causes it to precess around the magnetic field. Precession can be observed classically in a spinning top and the axis of the Earth. Figure 3.2 shows the precession on the three mentioned cases. [26, 55] N = µ ×B (3.2) Newton’s second law, equation 3.3, gives the relation between the torque and the change in angular momentum. Equation 3.4 is obtained by reorganizing the terms in 3.3, providing the mathematical description of the magnetic moment precession around B.[26] N= dJ 1 d µ = dt γ dt dµ = γµ ×B dt 21 (3.3) (3.4) (a) Precession in a spinning top. (b) Precession of the Earth’s axis. (c) Precession of a proton around an external magnetic field. ( c Kozlowski taken from Kozlowski [26] with permission.) Figure 3.2: Diagram of a precessing top, the earth, and a proton in an external magnetic field. 22 It can be seen from equation 3.4 that the frequency ω of precession is proportional to B. For this case, ω receives the name of Larmor frequency and equation 3.5 is known as the Larmor equation. The importance in signal localization will become clear in a moment.[26, 55] ω = −γ · B 3.1.1 (3.5) Generation of Macroscopic Magnetization A sample of a solid or liquid contains a huge amount of nuclei whose spins are oriented randomly in all directions. The total angular momentum of the sample is statistically (vectors in opposite directions cancel out) equal to zero, and thus the macroscopic magnetic moment is null as well. On the other hand, according to the quantum Zeeman effect in which an external magnetic field is present at the location of a spin, its magnetic moment interacts with the field in such a way that the energy states are given by equation 3.6. [26, 55] E = −µ · B (3.6) As a reference frame can always be chosen as desired, it is common to pick up one such that the Z direction points along B. With this in mind, the component in the z direction, µz = γ jz , and the angular momentum in that direction, jz = h¯ s, where h¯ is Planck’s constant divided by 2π and s is the spin value of the particle. Equation 3.6 can then be rewritten as shown in 3.7. For the case of the hydrogen atom, the spin value is 1/2 leaving just two possible energy states when an external magnetic field is present; one parallel to the field and one anti-parallel to the field. Figure 3.3 represents the energy split when the external magnetic field is applied to spin 1/2 particles.[26, 55] E = −γ h¯ sB (3.7) When the sample, which includes many hydrogen nuclei, is located in an external magnetic field, they will now be in one of the two allowed states. One might 23 Figure 3.3: Zeeman effect in spin 1/2 particles. The energy states split in two, one corresponding to the magnetic moment being aligned parallel to the field and the other one anti-parallel to it. ( c Xiang taken from Xiang [55] with permission.) be tempted to think that in this situation the spins pointing in the direction of the magnetic field N − and the ones opposite to it N + will cancel each other. Luckily, the probabilities of finding a nucleus in one of the states is given by the Boltzmann distribution, which is dependent on the energy of the state. The ratio between the population in each energy state is shown in equation 3.8 where ∆E is the difference in energy for the two states, k is the Boltzmann constant and T the temperature of the system in Kelvin. In a common water sample, at room temperature and an external magnetic field of 1.5 Tesla, the ratio is of around 1.0000102 (10 for every 106 protons) [55] meaning that a slightly higher population of spins pointing in the direction of the field is present. This causes a macroscopic magnetization M0 in the sample, with a magnetic moment that will precess around the external magnetic field in the same way as individual spins do when tipped from its equilibrium position; that is with the Larmor frequency.[26, 55] γ h¯ B ∆E N− = e kT = e− kT + N (3.8) Figure 3.4 shows the magnetization being a vector with a longitudinal component M|| = Mz zˆ and a transverse component M⊥ = Mx xˆ + My y. ˆ Quantum mechanically, the expectation value of the magnetic moment in which the Hamiltonian for 24 the system is Hˆ = −µ · B can be calculated as shown:[55] d i ˆ µ] >= i < Hˆ µ − µ Hˆ > < µ > = < [H, dt h h i = < (−µ · B)µ − µ(−µ · B) > h (3.9) By definition, µ = γ j and making use of the commutation relation of angular momentum j × j = ih j, equation 3.9 can be modified to obtain. [55] d < µ >= −γ 2 < B × j >= γ < µ > ×B dt (3.10) When N nuclei are present, the total magnetization M = N < µ > and thus: [55] d M = γM × B dt (3.11) From figure 3.4 it should be clear that the longitudinal magnetization is stationary, while the transverse component is the one that precesses around B. This can be expressed mathematically as dMz dt = 0 and dM⊥ dt = γ · M⊥ × B. The signal detected in NMR involves the precession of the macroscopic magnetization around B. [55] 3.1.2 The Rotating Reference Frame To simplify the description of the motion of the magnetization vector, a reference frame that rotates with the Larmor frequency is chosen. The concept is similar to the one shown in figure 3.5, in which the path of a man that travels from the north pole to the equator is shown in a stationary and a rotating frame. Figure 3.6 shows the behavior of the magnetization when it is being flipped in both the stationary and rotating frames. In the stationary reference frames the trajectory is a spiral, while in the rotating frame it is just a straight line. The magnetization does not precess when seen in the frame rotating with the Larmor frequency, and thus, will remain constant as if there was no magnetic field at all. This means that the effective magnetic field B0e f f = 0 when this frame rotates at the Larmor frequency.[26, 55] 25 Figure 3.4: Magnetization vector in a laboratory reference ( c Kozlowski taken from Kozlowski [26] with permission). frame Figure 3.5: Representation of a walk from the north pole as seen in a stationary (left) and rotating (right) reference frame. ( c Xiang taken from Xiang [55] with permission.) 26 Figure 3.6: Magnetization being tipped from the z axis as seen in a stationary reference frame (left) and a rotating one whose angular frequency is equal to the Larmor frequency. ( c Xiang taken from Xiang [55] with permission.) 3.1.3 The RF Pulse The thermal equilibrium magnetization M0 , which points in the z axis, can be tipped. From now on, the fixed external magnetic field will be noted as B0 and a second magnetic field B1 , which is perpendicular to B0 is used to change the direction of the magnetization. By leaving the amplitude of B1 constant but making it rotate around the z axis it is possible to tip M0 . It turns out that hydrogen atoms, at the external magnetic fields typical in NMR and MRI precess at a Larmor frequency in the order of radio frequency. The magnetic field B1 is applied in this frequency range and the process is called RF excitation.[26, 55] Seen in the rotating frame (figure 3.7), the effective field Be f f is then B1 and the magnetization will precess around it with a frequency ω1 = γB1 . Depending on the time τ for which the RF pulse is applied, a flip angle α by which the magnetization will be tipped from the z axis is given by equation 3.12.[26, 55] α = γB1 τ 27 (3.12) Figure 3.7: View in the rotating reference frame in which Be f f = B1 and the flip angle α depends on the time for which the RF pulse is applied. ( c Xiang taken from Xiang [55] with permission.) If the field B1 does not rotate close to the Larmor frequency, the effect of tipping the magnetization is very low and will only occur for an instant every time the magnetization and the RF pulse are in phase. On the other hand, even if B1 has a small magnitude but it rotates at a frequency similar to the Larmor frequency, the magnetization will be tipped. This is what is known as the resonance condition, in which at the Larmor frequency energy will be absorbed by the spins, and transition between energy states will occur. This does not occur at other frequencies. [26, 55] Although basically any flip angle can be achieved, α = 90◦ and α = 180◦ are commonly used; they receive the names of 90 and 180 degree pulses. The way this is done is by applying a known B1 pulse for the corresponding time to achieve the desired flip angle.[26, 55] 3.2 Relaxation Times Equation 3.11 suggests that when the magnetization is tipped from the z axis it will precess forever around it. However, in real life, the magnetization returns to its thermal equilibrium value M0 after some time in a process called relaxation. When being tipped, protons exchange energy with the environment in order to minimize 28 its potential energy. This process causes a relaxation of the longitudinal magnetization and is called spin-lattice relaxation. The transverse magnetization is also affected by interactions between spins known as spin-spin relaxation. The relaxation times are characterized by T1 and T2 for spin-lattice and spin-spin relaxation respectively.[26, 55] Equation 3.11 was modified by Bloch in order to account for the relaxation processes that occur, and a set of equations known as Bloch equations are shown in equation 3.13.[26, 55] Mx xˆ + My yˆ (Mz − M0 )ˆz dM = γM × B − − dt T2 T1 (3.13) This can be separated into components which include just the spin-lattice relaxation or the spin-spin relaxation. Equation 3.14 shows the differential equation for the longitudinal component of the magnetization Mz and whose solution is given by equation 3.15, taking into account an equilibrium magnetization M0 . A plot of the solution is shown in figure 3.8, in which it can be seen that T1 is the time constant of the exponential growth. In other words, T1 is the time required to reach 63% of the equilibrium magnetization in the longitudinal direction.[26, 55] dMz M0 − Mz = dt T1 (3.14) Mz (t) = Mz (t = 0)e−t/T1 + M0 (1 − e−t/T1 ) (3.15) For the transverse components, the differential equations are given in 3.16. Local differences in the magnetic field cause individual spins to precess at slightly different frequencies. The total transverse magnetization is reduced due to this. The solutions of these equations are presented in equation 3.17.[26, 55] 29 Figure 3.8: Plot of the longitudinal magnetization as a function of time after being tipped. Due to spin-lattice relaxation process, the magnetization reaches equilibrium after some time. T1 is the time required to recover 63% of the equilibrium magnetization. ( c Kozlowski taken from Kozlowski [26] with permission.) dMx Mx = γMy B0 − dt T2 dMy My = −γMx B0 − dt T2 Mx (t) = [Mx (t = 0)Cos(ω0t) + My (t = 0)Sin(ω0t)]e−t/T2 My (t) = [−Mx (t = 0)Sin(ω0t) + My (t = 0)Cos(ω0t)]e−t/T2 (3.16) (3.17) It has to be clear that ω0 = γB0 with the external magnetic field being constant in time, and that B1 has been turned off before the relaxation process described by these equations takes place. The solutions to the transverse magnetization differential equations show how the components in the xy plane precess with the Larmor frequency. It is simpler then to describe the transverse magnetization in a similar way as a phasor in which the magnitude and the phase are used instead as shown in equation 3.18.[26, 55] M⊥ (t) = |M⊥ (t)|e−t/T2 eiφt 30 (3.18) As the magnitude of the magnetization vector cannot become bigger than the equilibrium value, a condition in which T1 is longer than T2 always applies.[26] 3.2.1 Detecting the Induction Signal According to Faraday’s law shown in equation 3.19, when a magnetic flux Φ across a closed circuit changes in time, it generates a potential difference ε, according to Lens’ law, the direction of the current produced by this voltage is such, that the resulting magnetic field compensates for the flux change.[26, 55] ε =− dΦ dt (3.19) As discussed before, in the case of NMR the transverse magnetization precesses around B0 and this precession is the one that causes a change in the magnetic flux across a closed circuit. For NMR purposes, a coil is used as the closed circuit and the current generated around it is detected. Figure 3.9 shows a diagram of the coil and the signal detected for the current around it. This signal is known as the Free Induction Decay (FID) generated by the transverse magnetization and whose amplitude changes in time due to the relaxation effects.[26, 55] 3.3 Spin-echo Experiments Due to the spin-spin relaxation effects, the transverse magnetization is expected to decay proportional to e−t/T2 . However, when measured it was observed that the decay rate was higher than expected. The difference in Larmor frequencies due to spin generated magnetic fields inside the sample generate a field inhomogeneity in it. The spins then dephase (some precess slightly faster and some others slightly slower) and this creates a reduced magnetization thus decreasing the FID. This decay is characterized by T2∗ given by equation 3.20, where ∆B provides a measure of the inhomogeneities in the magnetic field within the sample.[26, 55] 1 1 = + γ∆B ∗ T2 T2 (3.20) Hahn [18] showed that the field inhomogeneity is recoverable by applying two 31 Figure 3.9: FID signal detected in the coil. ( c Kozlowski taken from Kozlowski [26] with permission.) 90◦ RF pulses separated by a certain time τ. Hahn’s spin echo method was modified by Carr and Purcell and whose diagram is shown in figure 3.10. A 90◦ pulse is applied at t = 0. After a certain time τ, a 180◦ pulse is applied in the same direction in the rotating reference frame. By specifying the coordinates in the rotating frame as (x , y , z ) this sequence is denoted as 90◦x − τ − 180◦x − τ − echo, which means that the pulses are applied in the x direction on that frame. When the 90◦ pulse is applied, the magnetization is tipped to the x y plane and points in the y direction. Field inhomogeneities cause the spins to dephase that reduce the magnetization as it can be seen for t = τ − in the figure. When the 180◦ pulse is applied at time t = τ, all the spins are rotated around the x axis causing them to refocus as it can be appreciated in the diagram for time t = τ + . Finally, when the spins are refocused at time t = 2τ a large induction signal is detected, called the echo. It has to be clear that the 180◦ pulse does not change the frequency at which 32 Figure 3.10: Spin behavior in the Carr and Purcell sequence. ( c Xiang taken fromXiang [55] with permission.) the individual spins are rotating and that’s why they are refocused after this.[26, 55] Meiboom and Gill made a second modification to the spin-echo sequence shown in figure 3.11 and named the CPMG sequence. This sequence is denoted as 90◦x − τ − 180◦y − τ − echo in which the difference with respect to Carr and Purcell’s sequence is the direction at which the 180◦ pulse is applied, i.e. the 180◦ pulse flips the individual spins around the y axis. In this case, the refocused spins will generate a magnetization that points in the positive y direction in contrast to the CP sequence which points in the negative direction. By doing this, the CPMG sequence is better if RF flip angles errors exist. If the refocusing pulses are not exactly 180◦ but are rather slightly greater, then at t = τ + the flipped magnetization will have a longitudinal component. In the case of the CP sequence, when a second refocusing pulse is applied along the x axis, the overshooting will increase the longitudinal component even more. The useful transverse magnetization is reduced due to the RF errors. In the case of the CPMG sequence, by applying the 180◦ pulse along the y axis, the overshooting cancels every second echoe.[26, 55] The spin-echo sequences just described can be extended to produce trains of echoes by applying 180◦ pulses at times t = τ, 3τ, 5τ, .. with the echoes occurring at t = 2τ, 4τ, 6τ, .... The signal follows equation 3.21, where S(0) is the signal at t = 0 and n is the echo number which is of course an integer. Figure 3.12 shows an echo from a spin-echo sequence in which now the signal intensity decays with the T2 value rather than T2∗ . The time at which the echo occurs is referred as the Echo Time (TE).[26, 55] 33 Figure 3.11: Modification of the CP sequence by Meiboom and Gill. ( c Xiang taken from Xiang [55] with permission.) Figure 3.12: Spin-echo diagram in which the signal decays with the relaxation time constant T2 as described by equation 3.21. ( c Kozlowski taken from Kozlowski [26] with permission.) S(2nτ) = S(0)exp 2nτ T2 (3.21) By acquiring the signal of the different echoes it is then possible to calculate the value of T2 . 3.4 Spatial Localization The localization of the NMR signal is performed in MRI by using three sets of coils that produce field gradients in the three orthogonal directions x, y, z. In general, the magnetic field perturbation ∆B, follows equation 3.22, and is zero at the isocenter of the coils and varies linearly with space coordinates. The superpo- 34 sition of the static magnetic field B0 and the perturbation, causes a linear dependence of the magnetic field with the position that is desired as presented in equation 3.23.[26, 55] 3.4.1 ∆B = r · G = xGx + yGy + zGz (3.22) B(r) = B0 + ∆B (3.23) Slice Selection When a magnetic field gradient is turned on spins located in the plane perpendicular to it, will have the same Larmor frequency. If the RF pulse with frequency ω given by equation 3.24 is applied, it will excite the spins in that plane. MRI acquisitions are usually made by acquiring data for many different slices. The way this is achieved, is by modifying the RF frequency while keeping the gradient. The gradient can be applied in any direction, providing the option to acquire oblique images which is a great advantage of MRI in comparison to other imaging techniques. [26, 55] ω = γ(B0 + r · G) (3.24) The gradient is applied in the z direction, which is used as the notation for slice selection coordinate. The slice thickness depends on the bandwidth ∆ω as shown in equation 3.25. The RF pulse is applied at the same time as the the slice selection gradient and it lasts on the order of milliseconds. This way, only magnetization within the desired slice and thickness will be excited.[26, 55] d= ∆ω γG 35 (3.25) 3.4.2 Frequency Encoding Frequency encoding is achieved when the Larmor frequency is made position dependent by using another magnetic field gradient. The distribution of magnetization in space is determined from the echo signal using frequency analysis with the help of the Fourier transformation. After the RF pulse is applied, the frequency encoding or reading gradient is turned on causing the difference in Larmor frequencies for different positions within the sample. Lines perpendicular to the gradient direction can be thought as having magnetization vectors precessing at the same frequency, but different from spins belonging to a different line. This means that the FID signal of each line individually has an intensity proportional to the number of spins in that line oscillating at the same Larmor frequency. However, the FID signal detected includes information from all the lines (i.e. all the sample).[26, 55] The acquired signal S is a function of time given by equation 3.26 in which the transverse magnetization takes the form shown in equation 3.18. As the gradient is present during signal acquisition, a change of variable can be made in such a way that the signal depends on k = γGt; i.e. the angular frequency ω from equation 3.21 is just k. This signal can then be written as shown in equation 3.27 in which P(x), known as a projection, is a function that depends only on x. A projection is a spatial distribution of the magnetization along the axis perpendicular to the direction of the frequency encoding gradient. This equation states that the signal detected is a Fourier transform of the projection, meaning that the P(x) can be determined by calculating the inverse Fourier transform of the acquired signal S(k) as presented in equation 3.28. [26, 55] s(t) ∝ M⊥ (r) · e−iω(r)t (3.26) +∞ S(k) = P(x)exp(ikx)dx (3.27) −∞ P(x) = 1 2π +∞ S(k)exp(−ikx)dk (3.28) −∞ In other words, applying the inverse Fourier transform generates the spectrum 36 of frequencies involved in the different perpendicular lines to the gradient. The intensity of a certain frequency depends on the amount of protons present and this then generates the profile of the sample across the gradient. The spectrum represents the projection of the sample in the direction perpendicular to the one of the frequency encoding gradient.[26, 55] 3.4.3 Phase Encoding Methods applied to image reconstruction in nuclear medicine imaging techniques that involve the acquisition of a series of projections to generate a two dimensional image could be applied to MRI with data collected as mentioned in the previous section. Projections in MRI are projections of a vector field which differentiates it from the ones obtained in nuclear medicine. This means that each projection profile has a complex value obtained by adding magnetization vectors in a line perpendicular to the gradient. The magnitude of the vector projection shows the scalar quantity of the spin density, but only if those spins are in phase. Magnetic field inhomogeneities, cause phase dispersion in the object making the reconstruction by projections very sensitive to the field inhomogeneities. [26, 55] Taking advantage of the vector property of the magnetization, allows a better data encoding. Using another gradient, the phase encoding gradient, a linear phase distribution can be created in the magnetization that contains spatial information. All this is then used to generate the vector projection reconstruction in which frequency encoding projections performed at different angles are avoided and replaced by just one frequency encoding projection with many different phase dispersion. This then allows to solve for the 2D distribution of the magnetization in the sample. [26, 55] Back again to equation 3.18, the phase of the magnetization depends on time. To determine the phase at the beginning of the acquisition, the components of the magnetization in the x and y axis are measured and the phase is calculated. A second encoding gradient in the direction orthogonal to both the slice selection and frequency encoding ones, and usually referred as the y direction, is used to localize 37 Figure 3.13: Example of a phase encoding procedure in which a gradient is applied with different amplitudes for three RF excitations. The magnetizations a, b, and c can be determined from the projections. ( c Xiang taken from Xiang [55] with permission.) the signal. The RF pulse is applied as usual and after the pulse the phase encoding gradient is turned on for a certain amount of time. A point located in the y axis will have a relative angular frequency ω = γyGy , which will cause a phase φ = γyGy T when the gradient is applied at time T . The vector projections P can be known from the frequency encoding, and the phase values of these projections, are also known from the time and magnitude of the phase encoding gradient.[55] Considering a case in which three different amplitudes for the phase encoding gradient are applied after the slice selection, but before the frequency encoding one, is a good example. As seen in figure 3.13, the magnetization of the three shown pixels, whose locations are known, wants to be determined and is given by the values a,b, and c. After applying the gradient for a time interval T , the vector projections p− ,p0 , and p+ can be written as follows: p− = aexp(−iθ ) + b + cexp(iθ ) 38 (3.29) p0 = a + b + c (3.30) p− = aexp(iθ ) + b + cexp(−iθ ) (3.31) The phase rotation of the magnetization is included in the exponential. The frequency encoding projections provide all the information needed to determine a, b, and c. If the set of equations 3.29 to 3.31 are expressed in matrix form, where F is a matrix that includes the phase information (i.e. the exponentials), m the magnetization coefficients, and P the projections, then equation 3.32 summarizes the process.[55] Fm = P (3.32) Usually a number n of points are collected for the y direction. The phase of the spins is increased linearly from negative to positive values by equally spaced steps, thus providing n different phases symmetrical around zero. In this sense, matrix F takes the form given in equation 3.33 where −n/2 ≤ j ≤ n/2 and −n/2 ≤ k ≤ n/2. [55] 2π jk 1 Fjk = √ exp −i n n (3.33) Every point in P is a weighted sum due to vector rotation of magnetizations in a line parallel to the y axis, in which the weighting factor is then an exponential that represents the vector rotation. F and F−1 are actually a Fourier transform and its inverse. It means that the 2D image of a sample can be obtained by two inverse Fourier transforms; one in the frequency encoding direction and the other one in the phase encoding one. [55] 3.5 Signal to Noise Ratio and Partial Volume Effects Signal to Noise Ratio (SNR) is defined as: 39 SNR = I σ (3.34) with I representing the signal intensity and σ the standard deviation of the noise in the image. The signal in an MRI image changes between different points and the noise is spatially invariant and receives the name of (i.e. white noise). The SNR is usually determined by using a phantom in which a Regions of Interest (ROI) is selected and the mean and standard deviation of the values in that region are used in equation 3.34, but it can also be determined in an image of a real sample by selecting an ROI in the background of the image and another one in the sample. Better methods though, measure noise from a scan without RF excitations or subtract images after a series of repeated scans. [55] SNR can be written as shown in equation 3.35 in which p is a constant depend- ing on the imaging equipment known as the noise power, Nx , Ny , Nz are the number of samples taken in the three encoding directions, ∆t is the time between samples of the signal, and M is the local transverse magnetization. [55] 1 SNR = √ M∆x∆y∆z Na Nx Ny Nz ∆t p (3.35) The MRI machine equipment and settings can then be modified in order to obtain a better SNR (the higher SNR the better quality images and data) taking into account the following:[55] • The higher the value of the static field, the higher the transverse magnetization which increases the SNR. • Reducing the size of the coil increases the SNR as the noise comes from the whole volume and if the volume is smaller then the noise power p is reduced. • If the slice thickness is big, then more signal is detected and the SNR is better. • Leaving the number of pixels constant, if the Field Of View (FOV) is increased then the voxel size is bigger and more signal is detected which increases the SNR. 40 • Increasing the matrix size changes the number of samples in each direction but at the same time reduces the voxel size so in general the collected signal for each voxel drops down and SNR is reduced. In general, to obtain a better SNR the size of the voxel should be increased. Another effect, known as partial volume effect, happens when two different tissue types are present in a sample. If tissue A is smaller than a certain voxel, that voxel is shared with another background tissue. The signal from the voxel is then a combined signal of both types of tissue. Its contrast in an MRI image is reduced according to the fraction of volume that occupies in the voxel. The size of the voxel should then be optimized to maximize SNR while minimizing partial volume effects. [17] 3.6 3.6.1 Some Pulse Sequences Spin-Echo Figure 3.14 shows the pulse sequence diagram of the spin-echo technique. The horizontal axis represents time where, as mentioned before, T E is the echo time and T R indicates repetition time, in which the magnetization is recovered, before the sequence is repeated.[26] A 90◦ RF pulse, applied at the same time as the slice selection gradient, usually has a shape of a Gaussian or a Sinc function in order to excite a small range of frequencies to provide better slice selection. Due to the spread of Larmor frequencies across the slice, caused by the presence of the slice selective gradient, there is a dephasing of the spins. The spins are flipped to the transverse plane at the time where the center of the pulse is reached, thus between the center of the pulse and the end of the pulse the spins at different locations across the slice precess with different frequencies. The negative lobe of the slice selecting gradient is applied in order to refocus the spins, whose area is then equal to half the total area of the positive lobe in order to achieve the perfect refocusing. The 180◦ RF pulse is applied after a time T E/2 with the slice selection gradient on in order to refocus the transverse magnetization in the selected slice and generate the echo.[26] 41 Figure 3.14: Spin-echo pulse diagram. ( c Kozlowski taken from Kozlowski [26] with permission.) The phase encoding gradient is applied before the echo acquisition to encode the spins in the y direction. This phase encoding gradient can be applied at any time before the echo is detected and while the RF pulse is turned off as the phase accumulation, due to phase encoding gradient, will remain the same. Finally, the frequency encoding gradient is turned on to encode the spins in the x direction. A negative frequency encoding lobe is used in the frequency encoding direction and then switched polarity in order to rephase the magnetization and detect the echo. Due to this, the frequency encoding direction is also called the reading direction as this gradient is on during signal acquisition. The magnetization will be in phase again when the area of the negative lobe equals half the area of the positive lobe. The sequence is repeated with a different value of the phase encoding gradient, the data is stored in the computer and later used to reconstruct the image. The frequency and phase encoding gradients should be turned on when the RF pulse is off as the gradients should affect the the whole transverse magnetization.[26] The spin-echo sequence can also be presented as shown in figure 3.15. It can be seen that the phase encoding gradient and the frequency encoding refocusing 42 Figure 3.15: Spin-echo alternate pulse diagram. ( c Kozlowski taken from Kozlowski [26] with permission.) lobe are now placed in another period of time, before the 180◦ pulse. However, the polarity has been changed as the 180◦ RF pulse flips the magnetization. Weighted Images The spin-echo signal is a function of the relaxation times as well as the amount of spins (more spins, more magnetization) present in the sample. If T E and T R are modified, different contrasts can be obtained in the images. Its signal can be described as shown in equation 3.36, where M0 is the equilibrium magnetization and F1 and F2 are functions depending on T1 and T2 respectively and given by equations 3.37 and 3.38. [55] S ∝ M0 F1 F2 −T R T1 F1 = 1 − exp F2 = exp −T E T2 The three type of images are as follows: 43 (3.36) (3.37) (3.38) • Proton Density Weighted (PDW): In this types of images, water content is a highly dominant source of contrast. The idea is to reduce the effect of spin-lattice and spin-spin relaxation processes in the signal detected. A repetition time long enough to allow the magnetization to return back to its equilibrium value is used; usually T R ∼ 5T1 . On the other hand, in order to reduce the decay in signal strength caused by the spin-spin relaxation process, a very short echo time, compared to T2 is required; T E << T2 . By doing this, equation 3.36 is basically M0 as the other terms will be close to unity.[26, 55] • T 2 W: These images are the ones in which the contrast highly depends on the spin-spin relaxation process but in which spin-lattice relaxation can be neglected. To take out the effect of T1 , according to equation 3.37, a high value of T R should be selected. To make the effect of T2 have a high contribution, F2 in equation 3.38 should have a high value, and this is achieved by selecting a long T E. In general T R ∼ 5T1 and T E > T2 . [26, 55] • T1 Weighted (T 1 W): In a similar way as T 2 W, this images provide more information from the spin-lattice relaxation while T2 effects are neglected. In this situation, F2 in equation 3.38 should be close to unity and that is achieved when T E << T2 . F1 in equation 3.37 on the other hand should have a high contribution and that is achieved when T R < T1 . [26, 55] T1 and T2 contrast is defined as the derivative of 3.37 and 3.38 with respect to their relaxation times respectively. With this in mind, it can be shown that the optimal parameters for the different types of images are the ones shown in table 3.1. [55] 3.6.2 Gradient Echo Another simple and common pulse sequence, shown in figure 3.16, is the gradient echo sequence. A 90◦ RF pulse is applied with the slice selection gradient as in spin-echo. A negative gradient in the reading direction is applied and then its polarity is changed and the signal is detected. In this case there is no 180◦ pulse used and the static field inhomogeneities will affect the signal. [26] 44 Contrast Parameters PDW T R >> T1 and T E << T2 T R >> T1 and T E ∼ T2 T E << T2 and T R ∼ T1 T2W T1W Table 3.1: Optimum parameters for the different types of contrast in spinecho images. Figure 3.16: Gradient-echo pulse sequence. The shaded regions in the slice selecting gradient indicate that the highlighted areas are equal. A similar situation occurs in the frequency encoding gradient lobe. When these areas are the equal, the spins are perfectly rephased. ( c Kozlowski taken froom Kozlowski [26] with permission.) 3.6.3 Rapid Acquisition with Relaxation Enhancement (RARE), Fast Spin Echo (FSE), or Turbo Spin-echo (TSE) The pulse sequence for this acquisition technique is shown in figure 3.17, in which the sequence begins with a 90◦ RF pulse, followed by a train of refocusing pulses. It is common to apply the refocusing pulses at time intervals of 2τ in which τ is the time between the initial RF excitation and the first refocusing pulse. When the refocusing pulses have a flip angle of 180◦ the spin echo signals are maximized. B1 field nonuniformity and imperfections in the slice profile, create limitations even 45 Figure 3.17: RARE pulse sequence. ( c Elsevier taken from Bernstein et al. [4] with permission). when the 180◦ flip angles are used. Some practical considerations need refocusing pulses with smaller flip angles. The refocusing pulses do not need to have the same flip angle throughout the complete train. The flip angles can be adjusted individually to regulate signal variations. The slice selection gradient is turned on for all the RF pulses in the train. The frequency encoding gradient is always the same for each echo detection. A phase encoding gradient is turned on between refocusing pulses but before the reading gradient is applied. The phase encoding lobes have the same pulse width but different amplitudes throughout the echo train, providing information for different phase values in a single Repetition Time (TR), and thus reducing the scanning time.[4] 3.6.4 Echo Planar Imaging Echo Planar Imaging (EPI) is one of the fastest pulse sequences available in MRI. It differentiates from common pulse sequences in the way in which the frequency and phase encoding gradients are applied. Bipolar readout gradients are used to 46 Figure 3.18: EPI pulse sequence. ( c Xiang taken from Xiang [55] with permission.) generate a train of echoes. A spin-echo EPI pulse sequence is shown in figure 3.18 in which a 90◦ RF pulse is applied followed by a 180◦ refocusing pulse. The slice selection gradient is identical to the one in the conventional spin-echo sequence. [4] 3.7 3.7.1 Measuring Diffusion Characterizing Diffusion Diffusion is a process of random Brownian motion of molecules throughout a medium. Einstein’s diffusion equation shown in 3.39 determines a way of characterizing this random movement for free water particles. The mean square distance travelled by a water molecule R2 in a time interval τ and in an N dimensional space has a diffusion coefficient D. Molecular weight, viscosity, and temperature influence diffusion in the solution. In a living organism, cells provide additional barriers by creating compartments in between. [3, 21] R2 = 2NDτ 47 (3.39) Figure 3.19: Diagram of isotropic and anisotropic diffusion. ( c John Wiley and Sons taken from Beaulieu [3] with permission.) Molecules might not necessarily diffuse the same distance in all directions. Isotropic diffusion is referred to pure liquids and other mediums in which no barriers are present to limit the motion of the particles or in which those barriers are randomly oriented. The particles will diffuse and will be able to travel the same distance in all directions. There are some other cases in which the barriers limit diffusion in certain directions and this is referred to, as anisotropic diffusion. A representation of these two is shown in figure 3.19 Diffusion can also be characterized using the Fick’s laws shown in 3.40 with boundary condition P(y, 0), where J represents the flux, D is again the diffusion coefficient, y is taken to be the direction of diffusion, and P is known as the propagator.[29] J = −D P δP =D δt 2 (3.40) P By solving Fick’s laws, the propagator is given as follows: 48 P(y,t) = exp(y2 /2σ 2 ) √ σ 2π (3.41) P(y,t) is basically a Gaussian probability density distribution in which σ represents the root mean square distance (in one dimension) travelled by the particle. Substituting Einstein’s relation 3.39 into 3.41, the propagator can then be expressed in terms of the diffusion coefficient D as given by 3.42. [29] P(y,t) = 3.7.2 exp(y2 /4Dt) √ 4πDt (3.42) Stejskal-Tanner Experiment Stejskal and Tanner introduced a time dependent gradient in the spin echo sequence, and derived an equation that allows measuring diffusion using the NMR signal obtained [50]. Figure 3.20 shows the suggested pulse sequence in which two gradients G, each of them applied for a time δ and with a time separation of ∆ are introduced in the phase encoding direction. Assuming that no motion is present in the particles of the sample then the two new gradients will have no effect on the signal detected as, due to the 180◦ pulse, all phases accumulated due to the first gradient pulse will be rephased by the second gradient pulse. If however, a spin moves a distance y during time ∆, a phase accumulation of γGδ y will exist by the end of the second pulse. [29] The detected signal will be given by equation 3.43 where S0 is the signal intensity detected when no diffusion is present (i.e. no motion). S = S0 exp(iγGδ y) (3.43) To account for diffusion in the whole region, the propagator is used and an integration is made as given by equation 3.44. S= P(y,t)exp(iγGδ y)dy (3.44) This equation can be modified by changing variable as q = γGδ and can be 49 Figure 3.20: Spin-echo sequence with two time dependent gradients used to acquire diffusion information. ( c MacKay taken from MacKay [29] with permission.) written as follows for time ∆. S∆ (q) = P(y, ∆)exp(iqy)dy (3.45) Equation 3.45 basically states that the signal detected is a Fourier transform of the propagator and thus the propagator can be obtained by detecting the signal and applying an inverse Fourier transform as follows. [29] P(y, ∆) = S∆ (q)exp(−iqy)dq (3.46) Substituting the propagator from equation 3.42 into equation 3.45, results in equation 3.47, which is basically the signal detected as a function of y, D, ∆, and δ. 1 S= √ 4πD∆ exp(−y2 /4D∆)exp(iγGδ y)dy The result of the integration is then 50 (3.47) Figure 3.21: Stejskal-Tanner sequence in which the lapse of time in which the time dependent gradients are applied is longer and is of the same order of magnitude than the time interval between them. ( c MacKay taken from MacKay [29] with permission.) S = S0 exp(−γ 2 G2 δ 2 D∆) (3.48) It has to be noted that equation 3.48 and figure 3.20 suggest an infinitely small δ compared to ∆. In real life it is not possible to have an infinitely small δ value on existing MRI scanners and the sequence is modified as shown in figure 3.21. [29] To take into account for the longer duration of the gradient pulses, equation 3.47 is modified as shown: 1 S= √ 4πD∆ exp(−y2 /4D∆)exp iγ G(t)tdt y dy (3.49) The solution is now S = S0 exp[−γ 2 G2 δ 2 (∆ − δ /3)D] (3.50) Equation 3.50 is the Stejskal-Tanner relation in which a variable b, called the gradient factor, is introduced as: 51 TE b = γ2 2 t G(t )dt 0 dt (3.51) 0 In the case of the sequence shown in figure 3.21, b = γ 2 G2 δ 2 (∆ − δ /3), and equation 3.50 can be rewritten as shown in equation 3.52. In isotropic diffusion the value of D is independent of the direction in space, and thus diffusion measurements can be done using only one direction of the diffusion sensitizing gradient. In an anisotropic medium, water diffusion is characterized by a second rank tensor shown in equation 3.53. [29] S = S0 e−bD (3.52) Dxx Dxy Dxz D = Dyx Dyy Dyz Dzx Dzy Dzz (3.53) Turning the attention back again at figure 3.19, the value of D in the isotropic case provides the same mean square distance travelled by a particle in any direction. This means that the possible space covered by this particle can be enclosed by a sphere whose radius is that mean square distance. When diffusion is anisotropic, different values are obtained for different directions and the sphere is deformed into an ellipsoid. By diagonalizing the tensor from equation 3.53, the three eigenvectors provide the directions of the axis of this ellipsoid, and the eigenvalues are related to the squares of the axis lengths. [21] The diagonalized tensor is represented by a diagonal matrix given by equation 3.54 in which λi represents the ith eigenvalue. [29] λ1 0 0 D = 0 λ2 0 0 0 λ3 (3.54) The tensor follows the following properties when rotated by any angle around 52 any axis: [21] • The tensor is symmetric Di j = D ji • A rotation by θ is equivalent to a rotation by θ + 180◦ around an ellipsoid axis. • A rotation by 90◦ around the x, y, or z axis exchanges the other two diagonal elements (Dzz exchanges for Dyy and viceversa when rotated about the x axis) and the corresponding off-diagonal element is multiplied by −1 (Dyz for the example given). • The diagonal elements Dii fall in the range λmin ≤ Dii ≤ λmax . When the eigenvalues are positive and real the following properties also apply: • The diagonal elements of the tensor are all positive. • The determinant of the tensor is positive. • ∀ i and j Dii + D j j > 2Di j • Any off-diagonal element is always smaller than the largest diagonal element. In order to measure the tensor properly, the diffusion sensitizing gradients must be applied in at least six non-collinear directions and the b value of equation 3.52 will take the form bi j = γ 2 δ 2 gi g j (∆ − δ /3) with i and j being x, y, or z. Using the symmetry property of the tensor the Stejskal-Tanner relation for anisotropic diffusion is then: [29] S = exp[−bxx Dxx − byy Dyy − bzz Dzz − 2bxy Dxy − 2bxz Dxz − 2byz Dyz ] S0 (3.55) Equation 3.55 has six unknowns, the different components of the diffusion tensor, and thus in order to generate the tensor at least six gradient orientations should be set and the ratio S/S0 is measured and S0 is measured when b = 0. [29] 53 3.7.3 Diffusion Anisotropy Indices In the case of anisotropic diffusion, a value of the diffusion coefficient is usually calculated, given by equation 3.56 in terms of both the eigenvalues and the tensor elements. It is known as the Apparent Diffusion Coefficient (ADC) or mean diffusivity. This value represents an average diffusion coefficient throughout the anisotropy ellipsoid, which is an indirect measure of the distance that a water molecule will travel when diffusing. [21] ADC = λ1 + λ2 + λ3 Dxx + Dyy + Dzz = 3 3 (3.56) In order to determine the degree of anisotropy of the diffusion ellipsoid, a quantification procedure has been generated. Several Diffusion Anisotropy Indices (DAI) have been proposed with different properties. Most commonly used (and useful) indices are single numbers that vary from 0 (isotropic) to 1 (anisotropic). The most widely used is the Fractional Anisotropy (FA) given in equation 3.57 in terms of the eigenvalues of the tensor, and in equation 3.58 in terms of the tensor elements.[21] 3 ∑(λi − λ )2 FA = i 2 ∑ λi2 (3.57) i FA = 3.7.4 3[(Dxx − ADC)2 + (Dyy − ADC)2 + (Dzz − ADC)2 + 2(D2xy + D2xz + D2yz )] 2[D2xx + D2yy + D2zz + 2(D2xy + D2xz + D2yz )] (3.58) Diffusion EPI Figure 3.22 shows a single-shot spin-echo diffusion EPI pulse sequence. This sequence is the most commonly used to measure diffusion, due to its high acquisition speed and motion intensitivity. Identical diffusion weighed gradient lobes are located on both sides of the 180◦ RF pulse. The direction is modified by changing the amplitude of each of the gradients in the three orthogonal directions and is being 54 Figure 3.22: Diffusion EPI pulse sequence. ( c Elsevier taken from Bernstein et al. [4] with permission.) shown as the shaded regions in the figure.[4] 55 Chapter 4 Image Registration In order to gain as much information as possible, different medical imaging modalities (e.g. X-ray, MRI, Computed Tomography (CT), etc) can be used to obtain images from a patient. These images provide information regarding shape, structure, size, spatial relationships of anatomical structures within the patient and the function and pathology or any abnormality [39]. Image registration, is the process in which two or more images of the same object that are taken at different times, different viewpoints, or by different sensors are overlayed. Usually, there will be differences in the images due to different imaging conditions. It is important to register images in which the information provided comes from the combination of all data sources. [57] According to Zitov´a [57] the applications of image registration differ but they can be divided in four main groups: • Different viewpoints for multi view analysis. • Different times for multi temporal analysis. • Different sensors for multi modal analysis. • Scene to model registration. This thesis focuses on multi modal analysis, in which images are acquired with different sensors in order to integrate the information from all of them to gain more 56 Figure 4.1: Sketch of the image registration process in which the voxel p from the image in the left is mapped into voxel q of the image on the right. The image has been taken from the manual given in reference ( c Klein and Staring [23] with permission.) detailed representations [57]. Image registration establishes a correspondence between different features from different modalities. It is becoming common to image patients several times either using the same modality or different ones. [39] Two images, the moving image IM (x) which is deformed to match the fixed image IF (x) are needed for the registration process. Both images are of dimension d defined on their own spatial domains. All the images that will be registered on a computer must be available in digital form. These images are made up of pixels that have a certain intensity associated to it. When two dimensional images are stacked together, a three dimensional volume is formed and each pixel is now a small volume called a voxel. Physical attributes are measured and the intensity that corresponds to each voxel is an average of that value in that small volume [39]. The registration is the process of finding a transformation T (x) = x + u(x) that aligns IM (T (x)) with IF (x). The transformation T (x) is a mapping from the fixed image to the moving image whose quality of alignment is defined by a similarity measure S. Figure 4.1 shows how voxels in one image are mapped to voxels in another image [23]. 57 4.1 4.1.1 Types of Transformations Dimensionality Transformations and Degrees of Freedom Due to the fact that we are three dimensional beings that move around in space, registration should be carried out in a four-dimensional space. However, approximations can be made in order to represent a body with fewer dimensions. The number of parameters that are required to describe a registration transformation is called degrees of freedom. [39] • 2D-2D: When the acquisition of the images is very well controlled, 2D images can be registered using a rotation and two translations which are in fact orthogonal. • 3D-3D: Assumptions that the internal anatomy of the imaged part of the body has not changed, leads to take the body as a rigid body in which three translations and three rotations will register the images. • 2D-3D: This is required when volumes and projection images must be registered. • Time: When sequences of images change with time, as for example images of a beating heart or other dynamic processes, time has to be considered. For the cases in which the object can be treated as a rigid body, the registration can be achieved by using translations and rotations. In 3D, a total of 6 parameters are needed; i.e. the transformation has 6 degrees of freedom, defined as translations in the x, y, and z directions and rotations α, β , and γ around those axes respectively. A transformation matrix Trigid can be found which maps the points from one image into the other. The transformation can be represented as a rotation R composed with a translation t = (tx , ty , tz )T . This can be applied to any point x = (x, y, z) in the image. [39] Trigid = Rx + t (4.1) The rotation matrix is formed by the rotation angles and is shown in equation 4.2. 58 cos(β )cos(γ) R = −cos(β )sin(γ) sin(β ) cos(α)sin(γ) + sin(α)sin(β )cos(γ) cos(α)cos(γ) − sin(α)sin(β )sin(γ) −sin(α)cos(β ) sin(α)sin(γ) − cos(α)sin(β )cos(γ) sin(α)cos(γ) + cos(α)sin(β )sin(γ) (4.2) cos(α)cos(β ) When a 2D to 3D rigid body registration is required, a projection of the 3D image into a plane is also needed. Combining the projection with the rigid body translations and rotations generates a useful 4 × 4 matrix, given in equation 4.3, that uses homogeneous coordinates.[39] cos(β )cos(γ) −cos(β )sin(γ) Trigid (x) = sin(β ) 0 cos(α)sin(γ) + sin(α)sin(β )cos(γ) cos(α)cos(γ) − sin(α)sin(β )sin(γ) −sin(α)cos(β ) 0 sin(α)sin(γ) − cos(α)sin(β )cos(γ) sin(α)cos(γ) + cos(α)sin(β )sin(γ) cos(α)cos(β ) 0 tx ty tz 1 x y z 1 (4.3) The projection of an object defined in the (x, y, z) space is taken along the z axis onto the u, v plane. The imaging system intrinsic parameters (u0 , v0 , ku , kv ) are used to characterize the projection. As an example, in the case of an x-ray, u0 and v0 define a point from which a vector normal to the u − v imaging plane points in the direction that passes through the x-ray source. The ku and kv parameters are the pixel sizes divided by the distance between imaging plane and the focal spot in the horizontal an vertical directions. The parameters can be determined by the calibration of the equipment or they can be taken as unknowns, in which case, 4 degrees of freedom are added to the registration algorithm. In the latter case, the transformation Tprojection shown in equation 4.4 takes the form of a 4 × 3 matrix that projects the 3D object along the z axis. When a Tprojection is applied to a point (x, y, z, 1)T , a vector (λ u, λ v, λ )T is obtained. To generate a 2D point (u, v)T both the first and second component of the vector are divided by the scaling factor λ . [39] ku Tprojection = 0 0 59 0 u0 0 kv v0 0 0 1 0 (4.4) To perform the 2D-3D rigid body registration, the composition of a Tprojection and Trigid has to be carried out as shown. T2D−3D = Tprojection Trigid (4.5) Bone structures are very well aligned even if the object is in different positions. However, when patients are repositioned, this changes the position of objects as well and the alignment cannot be described by translations and rotations only. Additional changes that include stretches and skews to extend the transformation are used to solve the problem; this generates what is known as an affine transformation. The rigid transformation conserves the distance between points on the object, while an affine transformation preserves parallel lines. [39] An affine transformation has twelve degrees of freedom and is described as follows: x a00 a01 a02 a03 y a10 a11 a12 a13 = T(x, y, z) = z a20 a21 a22 a23 1 0 0 0 1 x y z 1 (4.6) The affine transformations are used when some of the image acquisition parameters are unknown or when a small amount of shape variability occurs. [39] 4.1.2 Linear Transformations A linear transformation must satisfy: T(αxA + β x A ) = αT(xA ) + β T(x A ) ∀xA , x A ∈ RD (4.7) Affine transformations are usually confused and with the linear transformations; however, the translational part of it does not follow equation 4.7. An affine transformation is a composition of linear transformations with translations. Reflections do follow equation 4.7, but they are not desirable in medical imaging. Reflections are detected when the determinant of the transformation matrix has a 60 negative value. [39] 4.1.3 One to One Transformations Though different modalities can be used to generate images, a registration algorithm in medical imaging is usually applied to several images from the same patient. This leads to the assumption that the transformation T should take all the points in one image and transform them to points on the second image in a one to one correspondence. However, there are situations in which this is not the case. When the images have different dimensions or when different field of views and sampling are used between the images, parts of the patient that can be seen on one image can be missing on the other one. Registration of images from different patients or from the same patient, but before and after surgery, results in some structures in the images being absent in one of them. In these cases, a one to one transformation is either not possible to find or not desirable. [39] 4.1.4 Nonlinear or Non-rigid Transformations Non rigid registration is described by a transformation relating the two images, a similarity measurement between the images, and an optimization that determines optimal parameters. The optimization is a function of the similarity measurement. By adding additional degrees of freedom, linear or affine transformation models can be extended to nonlinear ones. The quadratic model of equation 4.8 for example, is defined by second order polynomials and has 30 degrees of freedom. Extending the model to third, fourth, and fifth order, 60, 105, and 168 degrees of freedom are obtained. Higher order polynomials though, introduce artifacts and it is uncommon to use them in nonrigid registration problems. They model global shape changes and have problems with local shape changes in anatomy. Figure 4.2 shows the behavior of the deformation field when different types of transformations are applied.[39] 61 Figure 4.2: Different types of transformations on a grid. ( c Klein and Staring [23] with permission.) a00 ... a08 a09 x2 y a10 ... a18 a19 Tquad (x, y, z) = z = a 20 ... a28 a29 1 0 0 0 1 y2 .. . x (4.8) 1 The Use of Basis Functions In order to describe the deformation, a linear combination of basis functions denoted by θi can be used as shown in equation 4.9. Orthonormal basis functions are usually used to describe the deformation field including trigonometric or wavelet basis functions. When trigonometric functions are used, a spectral representation of the deformation field is made, in which each frequency corresponds to a different basis function.[39] 62 x y a10 = Tquad (x, y, z) = a z 20 0 1 4.1.5 θ1 (x, y, z) .. ... a1n . ... a2n θn (x, y, z) 0 1 1 a00 ... a0n (4.9) Splines Many registration techniques use splines and they take into account that some landmarks, usually called control points, can be identified between images. To determine the control points, anatomical and geometrical landmarks that can be seen in both images are usually used. Control points arranged with a constant separation across the image form a regular mesh forming a parametrization of the transformation; these are not anatomical or geometrical landmarks and are called pseudo-landmarks. The spline-based transformations interpolate or approximate displacements at the control points in order to map the location of the point in the target image into its corresponding point in the source image. A smoothly varying displacement field is generated between the control points in which the interpolation condition is given by equation 4.10. φi and φi represent the location of the control point and its equivalent in the target and source images respectively. [39] T(φi ) = φi i = 1, ..., n (4.10) Thin Plate Splines Thin plate splines, equation 4.12, are defined as a linear combination of n radial basis functions θ (s) given by equation 4.11. Multi quadratics and Gaussian are other choices for radial basis functions.[39] |s|2 log(|s|) in 2D θ (s) = |s| in 3D 63 (4.11) n t(x, y, z) = a1 + a2 x + a3 y + a4 z + ∑ b j θ (|φ j − (x, y, z)|) (4.12) j=1 When the transformation is defined by three separate thin-plate splines, T(t1 , t2 , t3 )T the coefficients ai describe the affine part of the transformation while the b j coefficients describe the non-affine part. A total of 3n linear equations are obtained from equation 4.10. A total of 3(n + 4) coefficients must be determined, thus requiring 12 additional equations that guarantee that the b j coefficients add to zero and that the cross products with the control points coordinates is also zero. By taking a as the 4 × 3 vector of a coefficients, b as the n × 3 vector of the non-affine coefficients, and defining Θ as a kernel matrix whose elements Θi j = θ (|φi − φ j |), equation 4.13 can be solved algebraically to find a and b. The thin-plate spline transformation that interpolates the displacements at the control points is then known.[39] Θ Φ b ΦT 0 a = Φ 0 (4.13) Thin-plate splines are widely used to incorporate constraints into the transformation model. B-splines Freeform Deformations (FFDS) have been widely used in computer animations and graphics. This tool works really well when modelling 3D deformable objects and are useful in image registration. The idea is to manipulate a mesh of control points in order to deform the object. This deformation controls the shape of the 3D object and produces a smooth and continuous transformation. In the case of thin-plate splines, the control points could be arbitrarily configured, however FFDs require a regular mesh of control points that are uniformly spaced.[39] FFDs based on splines are defined on Ω = {(x, y, z)|0 ≤ x < X, 0 ≤ y < Y, 0 ≤ z < Z} and are called the domain of the image. The mesh of control points Φ, is formed by nx × ny × nz points with spacing δ between them. The displacement 64 field u in this case can be expressed as: [39] 3 3 u(x, y, z) = ∑ 3 ∑ ∑ θl (u)θm (v)θn (w)φi+l,j+m,k+n l=0 m=0 n=0 x δ i= u= x δ − 1, j = − x δ , v= y δ − 1, y y δ − δ k= w= z δ −1 z z δ − δ (4.14) θl represents the l-th basis function of the B-splines: θ0 (s) = (1 − s)3 /6 θ1 (s) = (3s3 − 6s2 + 4)/6 θ2 (s) = (−3s3 + 3s2 + 3s + 1)/6 θ (s) = s3 /6 (4.15) 3 In the case of thin-plate splines, each radial basis function contributes to the transformation and each control point has a global influence on it. When trying to model local deformations, thin-plate splines do not provide the optimal way of doing it. When the number of control points is large, the complexity of radial functions limits their use. On the other hand, FFDs are locally controlled and computationally efficient, even for a high number of control points. As an example, basis functions for the cubic B-splines have the property that when a control point is changed, the transformation is only affected in the local neighborhood of that point. [39] 4.2 Image Registration Algorithms Image registration algorithms’ purpose is to find the transformation T when the fixed and moving images are given. These algorithms are divided into three categories: [39] • Identified points in the image. • Delineated surfaces from the images. • Voxel intensity values. 65 4.2.1 Landmark Based Registration This type of procedure is based on the identification of corresponding points in the images to be aligned. The registration is carried out for those points and the transformation for the whole image is inferred from it. If the subject can be thought of as a rigid body, three corresponding non-collinear points, called fiducial markers, are sufficient. In practice it is common to use more than three points as errors become smaller with higher number of selected points. The algorithm consists of finding the centroid of each set of points (in fixed and moving images). Calculating the distance difference between the centroids specifies the required translation needed for the set of points. The set of points is then rotated around the new centroid and the sum of square distances between corresponding points is minimized. The root mean square error is usually stored by the algorithm. [39] The mathematical procedure to find this transformation is called the Procrustes Problem, which is a fitting problem of the least square type. Two configurations of N points, P = {pi } and Q = {qi }, for D dimensions are known. The goal is to find the transformation T that minimizes G(T) = |T(P) − Q|2 , in which P and Q are N × D matrices whose rows specify the coordinates of points pi and qi respectively, and T(p) is the matrix of transformed points.[39] The algorithm steps are as follows: 1. The mean values of pi and qi are calculated and subtracted from each point in order to reduce the problem to an orthogonal Procrustes problem. pi → pi − p qi → qi − q (4.16) 2. The orthogonal rotation R needs to be found. The correlation matrix K := PT Q, whose dimensions are D × D, states how points in Q are predicted by the ones in P. P = [pT1 , ..., pTN ]T and Q = [qT1 , ..., qTN ]T are matrices of row vectors, K = Σi Ki , taking into account that Ki := pi qTi . With this in mind, 66 K = UDV T ⇒ R = V ∆U T ∆ := diag(1, 1, det(VU T )) (4.17) K as shown in equation 4.17 is the Singular Value Decomposition of K. 3. When calculating the determinant of VU T two values are possible: +1 will show that the rotation has no reflections and −1 will show that the rotation includes reflections. As stated before, reflections are not desired in medical imaging registration. 4. The translation t = q − Rp is then applied. This Procrustes algorithm can be used to find rigid-body or affine transformations. The use of point landmarks is also present in non-affine transformations as the thin-plate spline.[39] 4.2.2 Surface Matching The boundaries of surfaces can be more easily identified in medical images than landmarks, and some segmentation algorithms locate the high contrast surfaces. When equivalent surfaces can be automatically segmented from the two images, rigid-body registration can be made fitting the surfaces together. The first algorithm used for this is called the head and hat algorithm. In this algorithm, the contours of the surfaces are drawn on different slices of one modality (the head). Points that correspond to the same surface but are present in the other modality are also identified (the hat). The computer attempts to fit the points from the hat to the head contours and is refined by iterations. The sum of squares of the distance between each point on the hat with respect to the head is minimized. The algorithm seems to fail when symmetry to rotation is present on the surfaces. To measure how close the points are to the head, the direction towards the centroid of the head is chosen. The Powell optimization algorithm is used in which a one dimensional optimizations are performed followed by one another. This finds the best solution around the six degrees of freedom after which it returns to the first 67 degree of freedom, and it stops when it fails to find a new solution that has a lower tolerance factor than the one it has found at that point. [39] Precomputing the distance from every point in space to one of the surfaces to be registered using a distance transform improved the performance of the head and hat algorithm. The distance transform speeds up the computation at each iteration. The transform is applied to a binary image and pixels are labelled with the distance from the surface of the object.[39] Another algorithm used in surface matching registration is the Iterative Closest Point (ICP). One surface is represented by a set of points, the collected data P, while the other one, the model data χ, could be represented by point sets, line segment sets, implicit surface, parametric curves, triangle sets, or parametric surfaces. In medical imaging, surfaces represented by triangular patches called facets are the most common ones due to the availability of algorithms to delineate them from the medical images. The ICP finds the closest model point for each data point and it then calculates the least-square rigid-body transformation using those points. A tolerance threshold is set and the algorithm iterates to find the local minimum. [39] The algorithm works as follows: [39] 1. The data surface P is converted to a set of points {pi }, while the model data χ remains in its original representation. In medical imaging both sets of data are obtained from radiological images and the model might be complemented with stereo video from surgery. For each point pi on the surface P, a point x is found on the model surface χ such that the distance between them is the minimum, and a set of closest points {qi } is generated. [39] d(pi , χ) = min ||x − pi || x∈χ (4.18) For a triangulated model representation χ, most common in medical imaging, a set of triangles {ti } is also obtained. 2. Using a linear interpolation across the facets ti whose vertices are r1 , r2 , and 68 r3 , the smallest distance between the point pi and the triangle ti is given by equation 4.19 in which u ∈ [0, 1], v ∈ [0, 1], and w ∈ [0, 1]. d(pi , ti ) = min u+v+w=1 ||ur1 + vr2 + wr3 − pi || (4.19) 3. qi = (ur1 , vr2 , wr3 ) is then the closest model point to the data point pi . 4. The points {pi } and {qi } are registered using a least square registration method, like the Procrustes method explained before. This gives a known rigid-body registration transformation. 5. The transformation is applied to the points {pi } to find points pi . 6. The process is repeated until the mean square error is smaller than a specified tolerance value. The algorithm should be started several times with different estimates of rotation alignments and the minima of the minimum obtained should be chosen at the end. [39] 4.2.3 Voxel Intensities This type of algorithms are becoming popular as they use the intensities in the two images alone, avoiding the segmentation or delineation of corresponding structures. As they use almost all of the data in each image, errors caused by noise or intensity fluctuations are averaged. The registration transformation τ 1 is calculated using a measure from the voxel values in the images instead of geometrical structures and optimizing it. In contrast with the other registration algorithms, voxel similarity measures need special care when images are obtained from different modalities. In general, there is no relationship in the intensities of images A and B, and a simple arithmetic operation to quantify for misregistration is inexistent in registration of intermodality images. However, though images taken using different modalities give complementary information, they also have a big amount of 1 Notation changed as for this type of registration the voxel values of the transformed image depend on the type of interpolation required for the registration. 69 shared information. People have been trying to solve this issue by preprocessing one of the images and making it look similar to the other. As an example, a remapping of the high CT intensities to lower ones can help the image to look more like an MRI one. A second technique, in which scale-space derivatives are applied to the images to make it easier to identify intensity changes that should be similar in both images at certain scales is also used. The newer algorithms include techniques that can be applied to intra and inter modality registration without the need of preprocessing.[39] An image can be thought of as a mapping of points with a certain domain Ω to intensity values. A : xA ∈ ΩA → A(xA ) B : xB ∈ ΩB → B(xB ) (4.20) Images can have a different FOV, in which case the domains ΩA and ΩB will be different. The registration, as stated before, finds the transformation that maps a point in image B onto a point in image A, from ΩA to ΩB within the overlap of both of them ΩTA,B that can be defined as: ΩTA,B = {xA ∈ ΩA |T −1 (xA ) ∈ ΩB } (4.21) The voxel similarity measures require the analysis of level sets in the images. For a given image A, an isointensity set with value A is a set of voxels in A in which Ωa = {xa ∈ ΩA |A((xA ) = a} (4.22) Not all algorithms use isointensity sets but rather use small bins of intensities; here a will be used to describe both cases accordingly. The level set in the overlap domain when voxel similarity measures are used are functions of the transformation T. The intensity set is defined as ΩTa = {xA ∈ ΩTA,B |A(xA ) = a} (4.23) The transformed image B on the other hand has a different definition for the 70 intensity set of value b in the transformed image Bτ . ΩTb = {xA ∈ ΩTA,B |Bτ (xA ) = b} (4.24) It has to be clear that the domain over which the transformation τ is valid is ΩTA,B , which is smaller than ΩA or ΩB and is also a function of the transformation T . Isointensity sets used by the voxel intensity algorithms are the sets of fixed intensities in ΩTA,B , and if the algorithms are too sensitive to changes in the overlap domain they are not reliable. [39] Medical images are also discrete, sampling the object with a finite number of voxels. The sampling is different for images A and B and it may or may not be isotropic. This discretization is important when the registration is made. The domain is then defined as: Ω =: Ω ∩ Γζ (4.25) In this case Ω is a bounded continuous set that defines the volume of the image and Γ is an infinite discrete sampling grid, characterized by a spacing ζ = (ζ x , ζ y , ζ z ) in the three dimensions. This grid is usually different for images A and B and this is denoted by ΓζA and ΓζB for ΩA and ΩB respectively. For a given T, the intersection of discrete ΩA and ΩB is very likely to be empty as sampling points might not overlap exactly. So in order to compare A and B with any estimate of T, an interpolation between sample positions taking into account sampling spacings ζA and ζB is necessary. Fast interpolation algorithms have errors that introduce blurring and ringing in the image, changing the histograms that will change the isointensity sets. When the image that is going to be transformed, B, has a higher sampling resolution than the reference, A, aliasing can occur when B is resampled from ΩB to Bτ in ΩA . B should be blurred first with a filter whose resolution is ζA or lower before resampling. τ has to take into account discrete sampling as it maps position and intensities at those positions.[39] 71 Voxel Similarity Measures Intensity Difference The Sum of Squared Intensity Differences (SSD) is the simplest voxel similarity measure applied to the images. SSD = 1 Nx ∑ |A(xA − Bτ (xA )|2 (4.26) T A ∈ΩA,B This equation takes image A, as composed of N voxels whose locations are given by xA on the overlapped domain ΩTA,B . Viola [53] showed that when images differ only by Gaussian noise (which is not the case for intermodality registration) the SSD, shown in equation 4.26, is the optimum measure. SSD is sensitive to large intensity differences in voxels between the images A and B. To reduce the effect of those pixels the Sum of Absolute Differences (SAD) is used instead of SSD as shown.[39] SAD = 1 Nx ∑ |A(xA − Bτ (xA )| (4.27) T A ∈ΩA,B Correlation Instead of thinking that the registration will leave images that differ only by Gaussian noise, the intensity values between the images could be thought as having a linear relationship. If this is the case, the similarity measure that should be used is the Correlation Coefficient (CC). If A and B are the averaged voxel values for images A and Bτ respectively in the domain ΩTA,B then the CC is given by ∑ CC = (A(xA ) − A) · (Bτ (xA ) − B) xA ∈ΩTA,B 12 ∑ (A(xA ) − A)2 · xA ∈ΩTA,B ∑ (Bτ (sA ) − B)2 xA ∈ΩTA,B which is a normalized version of the cross correlation measure 72 (4.28) 1 Nx C= ∑ A(xA ) · Bτ (xA ) (4.29) T A ∈ΩA,B The correlation technique can be applied both in the spatial and frequency domains.[39] Ratio Image Uniformity (RIU) RIU method uses a derived ratio image that is calculated from the images A and B. The transformation τ is calculated by iterating and maximizing the uniformity of the ratio. Quantification of it uses the normalized standard deviation of the voxels in the ratio image. The ratio is calculated as R(xA ) = A(xA ) Bτ (xA ) ∀xA ∈ ΩTA,B (4.30) It is useful to take the intermediate ratio R taking N voxels within ΩTA,B R= 1 Nx ∑ R(xA ) (4.31) T A ∈ΩA,B The RIU is then calculated as: [39] 1 N ∑ (R(xA ) − R)2 xA ∈ΩTA,B RIU = R (4.32) Partition Intensity Uniformity (PIU) PIU method was first used for intermodality registration between Positron Emission Tomography (PET) and MRI. This algorithm assumes that there is a small corresponding range of intensities in the PET image for each intensity in the MRI one. When registering images A and B, the PIU can be calculated by adding the normalized standard deviation of the voxel values in B for the value a from image A or by calculating the sum of the normalized standard deviations of the voxel values in image A for intensity value b from image B. [39] 73 PIUB = ∑ a na σB (a) N µB (a) PIUA = ∑ b nb σA (b) N µA (b) (4.33) Where na and nb represent the number of pixels with intensity a and b in images A and B respectively. The terms in equation 4.33 follow the following properties: µB (a) = σb (a) = 1 na 1 na ∑ Bτ (xA ) ΩTa τ µA (b) = 2 ∑(B (xA ) − µB (a)) σA (b) = ΩTa 1 nb 1 nb ∑ A(xA ) ΩTb ∑(A(xA ) − µA (b))2 ΩTb (4.34) Information Techniques When registration by voxel similarity is made, it can be said that the goal of it is to maximize the shared information from both images A and B. If the images are well registered, duplicates from known structures should not occur (e.g. two noses or four ears). With this in mind, registration can be thought of as a reduction of information in the combined image, in which a registration metric is used to measure it. Figure 4.3 shows a diagram of this, where the nonregistered images show the presence of more information. The Shannon-Wiener entropy measure H is the most commonly used in signal and image processing, where H represents the average information given by a set of symbols whose probabilities are Pi . [39] H = − ∑ pi logpi (4.35) i Equation 4.35 assumes that: 1. The functional is continuous in pi . 2. Taking n as the number of symbols i, if the probability pi = 1/n, H is monotonically increasing in n. 3. If a choice is broken down into a series of choices, the original value for H should be H(p1 , p2 , p3 ) = H(p1 , p2 + p3 ) + (p2 + p3 )H which is the weighted sum of the constituent H 74 p2 p3 p2 +p3 p2 +p3 , Figure 4.3: Diagram of two images being registered in which it can be seen how the misalignment of the images contains more information that the correctly registered one. The only functional that fulfils all the three conditions is H as stated in equation 4.35. The entropy value will be maximum when all symbols have exactly the same probability, and it will be null when one of the probabilities of occurrence is 1 while the others are zero. When the probabilities of the different symbols is made more uniform, the entropy is increased. [39] When talking about image registration, two symbols (one for each image) at voxel locations are available to find the estimate of τ. The joint entropy is what measures the amount of the information when the images are combined. The joint entropy is given by: [39] H(A, B) ≤ H(A) + H(B) (4.36) When the images are not related, the joint entropy is the sum of the entropies of both images independently. As the images are more and more related, the joint entropy becomes smaller in comparison to the sum of the independent entropies of the images. Two-dimensional histograms, made by plotting the intensity of corresponding voxels in the overlapping domain one against the other, are commonly used to evaluate the registration. The joint histogram is normalized dividing by the total number of voxels N in the overlap domain, in which case can be used as a joint Probability Density Function (PDF) pτAB (which depends on the estimate τ of the transformation) of the two images. The joint entropy for the images is given 75 by: H(A, B) = − ∑ ∑ pτAB (a, b)logpTAB (a, b) a (4.37) b The number of elements in the PDF is determined by the range of intensity values or a reduced number of bins from that range. Figure 4.4, taken from the book of Neuman [39], shows 2D histograms. Bright regions in the histograms get darker and some dark regions become lighter when misregistration occurs. The misregistration reduces both the highest values and the number of zeroes in the PDF, thus increasing the entropy. So in general, the registration should produce a small number of PDF elements with values as high as possible, while giving as many zero probability elements in order to minimize the joint entropy. [39] The joint entropy, that depends on τ and pTAB , has a high dependence on the overlap domain ΩTA,B and on the interpolation algorithm used to generate Bτ on each iteration. To solve the overlap problem, the information that is being contributed to the overlapping volume and the joint information from the registered image is also taken into account. The contribution is the entropy of the part of the image that overlaps with the other image volume. [39] H(A) = − ∑ pTA (a)logpTA (a) (4.38) H(B) = − ∑ pτB (b)logpτB (b) (4.39) a b pTA and pτB are the projection of the joint PDF onto the intensity axes of images A and B, and are known as marginal probabilities. The marginal entropies are not constant during the registration process, as the portion of each image that overlaps the other changes with each estimate of T. B is also transformed to Bτ in each iteration and the required interpolation also changes the probabilities. The Rate of Transmission of Information (RTI) was stated by Shannon as part of communication theory and gives tools to calculate the joint entropy from the marginal ones. This is better known as mutual information I(A, B) used for intermodality medical imaging registration.[39] 76 Figure 4.4: Two-dimensional histograms. (a) MRI images of the head. (b) MRI and CT images of the head. (c) MRI and PET images of the head. The left column is generated when the images are aligned. The center column shows the histogram when the images are translated by 2mm. The right column shows the images when translated by 5 mm. ( c SPIE taken from Hill et al. [20] with permission from the publisher.) I(A, B) = H(A) + H(B) − H(A, B) = ∑ ∑ pτAB (a, b)log a b pτAB (a, b) pTA (a) · pτB (b) (4.40) Mutual information is maximized in order to obtain the best possible registration. Thinking in terms of conditional probability p(b|a), in which B takes the value b when it is known that A has value a, the conditional entropy is given by: 77 H(B|A) = − ∑ pτAB (a, b)logpτ (b|a) = H(A, B) − H(A) (4.41) a,b Knowing this, equation 4.40 can be written as: I(A, B) = H(A) − H(B|A) = H(B) − H(A|B) (4.42) When the intensity of A(XA ) is known, the term H(B|A) is zero and the value corresponding to the intensity of Bτ can be predicted. Mutual information registration consists of finding the transformation that takes A to predict Bτ as best as possible in the overlap. If the images of the same object are correctly aligned, and the value of the intensity of a voxel in A is known, then the uncertainty for the corresponding value in image B is reduced. [39] Another still remaining problem is that changes in low intensity regions in the overlap (e.g. noise) contribute to mutual information in an uncontrollable way. Normalizing the joint entropy helps to solve this problem. The following equations are ways of normalizing this: I˜1 (A, B) = 2I(A, B) H(A) + H(B) (4.43) I˜2 (A, B) = H(A, B) − I(A, B) (4.44) H(A) + H(B) I˜3 (A, B) = H(A, B) (4.45) Equations 4.43 and 4.44 were suggested by Maes et. al. [30] and equation 4.45 was proposed by Studholme [51], showing that it works much better than mutual information without normalization. Equation 4.43 can be related to equation 4.45 as follows: [39] H(A) + H(B) I(A, B) 1 I˜3 (A, B) = = +1 = ˜ H(A, B) H(A, B) I1 (A, B) − 2 78 (4.46) 4.2.4 Optimization When the Procrustes algorithm is used or when two images that have very similar intensities are registered, the transformation that establishes the correspondence can be calculated in a straightforward way; these are the only two cases in which this happens. The others require iterations to improve the initial estimate of the transformation. Each iteration uses the actual estimate of the transformation and calculates a similarity measure. The algorithm then generates a new estimate and calculates the similarity measure again and again, until it converges and the similarity measure cannot be improved (above a certain tolerance level). Sometimes the algorithm might converge to what is known as a local optimum. [39] Thinking in terms of the parameters of the similarity measure, each point in that parameter space is a different estimate of the transformation. Non-rigid registration algorithms have more degrees of freedom and thus their parameter spaces have more dimensions. This parameter space can be seen as a high dimensionality image, in which the value of the similarity measure that gives the estimate for the transformation is the intensity at each location. If for example, low intensities mean that a good similarity value is obtained and high intensities mean that a poor similarity value is obtained, the ideal parameter space is one such that it has low intensity values that monotonically increases with distance from the optimum position. The purpose of the optimization is to determine the optimum location when an estimate is given. [39] Multiple optimum points can be present in the parameter space and the algorithm can converge to the wrong one. Some small optima can be caused by interpolation artifacts or a good localized alignment of the intensities. Blurring the image before the registration is performed helps reduce the local optima from the parameter space. In general, registering the images at low resolutions to find the transformation, and then using that result as the initial estimate of registration at higher resolutions is very common to avoid this problem. Although this helps, it does not necessarily solve the problem of multiple optima and starting the optimization with several estimates that lead to multiple solutions is required. Usually 79 the solution with the most appropriate similarity measure is chosen, but in voxel similarity measures that usually leads to a local optima instead of the global one. To show this point, Neuman [39] gives an example in which, when registering images using the joint entropy, a very good value of the similarity measure is found when the overlap of the images is just air. This has a very low entropy that is lower than in the case of the correct alignment.[39] 80 Chapter 5 Previous Work in Prostate DTI Several research groups have been trying to identify Prostate Cancer (PCA) and its tumour grade by using diffusion techniques and identifying if the ADC and FA values in different regions of the tissue provide any information that can lead to a successful diagnosis of the disease. Ultrasound guided biopsy (section 2.3), which is the standard method for PCA diagnosis, misses between 20 to 30% of the tumours [56]. MRI has gained acceptance in the evaluation of prostatic carcinoma before surgery in terms of extracapsular spread, seminal vesicle and neuro-vascular bundle invasion, and lymph node metastases [16]. MRI performed with an endorectal coil together with a phased-array pelvic coil has been shown to be one of the best techniques to localize and stage PCA. Local metastases can be studied by acquiring and combining high resolution images using T 2 W images, spectroscopy, and dynamic techniques [31]. T 2 W images however, have low sensitivity and specificity (e.g. prostatitis and tumoural infiltration both present signal decrease on this type of images) [16]. Since 2007, reports show that the use of DTI in combination with T2W images accurately detect PCA [52]. Miao and his group demonstrated that Diffusion Weighted Images (DWI) were much better in detecting PCA than normal T2W images [37]. The information regarding the microscopic organization of tis- sue is what makes DWI so successful. DTI was initially applied to brain studies but it has been gaining importance in other applications such as PCA when Sinha et. al. performed a preliminary study on 6 healthy volunteers. They concluded that the prostate has an architectural anisotropy, in which early detection of changes in its 81 structure can be identified with DTI. [31, 47] 5.1 Bashar’s Study The objective of Bashar’s work goal was to measure, for the first time, the ADC values in certain anatomical regions of the prostate gland. He recruited 7 healthy volunteers and 19 patients with high PSA and positive DRE (section 2.2) and scanned them on a 1.5 T MRI scanner. T 2 W and DTI data were acquired using the following settings in which no endorectal coil was used: [2] • Slice thickness of 4 mm with a 1 mm intersection gap was used for the T 2 W images. • Slice thickness of 7 mm with no information regarding any gaps provided was used for the DTI images. • FOV of 16 cm. • Matrix size for the T 2 W is not specified. • Matrix size for the DTI images is of 96 × 96. • Six b-values, b = 64, 144, 257, 401, 578, 786s/mm2 were selected. ROI ’s were selected that included cancer in the peripheral zone and healthy tis- sue in the NPZ and CG and the mean ADC value was calculated for each region. Figure 5.1 shows the behavior of the ADC in the peripheral zone of the prostate for the patients, while table 5.1 summarizes Bashar’s results. This work was published in 2002 and the details are available in reference [2]. 82 Figure 5.1: Behavior of the ADC value in the patients recruited by Bashar. ( c John Wiley and Sons taken from Bashar [2] with permission.) Volunteer Patient Tumour (in peripheral zone) NPZ CG 1.91 ± 0.46 1.82 ± 0.53 1.63 ± 0.30 —– —– 1.38 ± 0.52 BPH —– 1.62 ± 0.41 Table 5.1: ADC values (×10− 3mm2 /s) for the different zones in volunteers and patients from Bashar’s experiments. ( c John Wiley and Sons taken from Bashar [2] with permission.) 5.2 Sato’s Group Study In this study, Sato et. al recruited 29 patients with suspected PCA due to its elevated PSA serum levels during 2003 and early 2004. Ultrasound guided biopsy of the prostate was performed after the imaging of the patients. A 1.5 T scanner was used in which T 1 W, T 2 W, and DTI images were acquired with the following settings:[46] • Slice thickness of 7 mm with a 1.4 mm interslice gap. • FOV of 250 mm. • Matrix size of 253 × 384 was used for the T 1 W and T 2 W images. • Matrix size of 119 × 192 was used for the DTI images. 83 (a) ADC values of all ROI’s with cancer. (b) ADC values for peripheral zone. (c) ADC values for transition zone. (d) ADC values for ROI’s without cancer in both the trasition an peripheral zones ROI ’s in the ROI ’s in the Figure 5.2: ADC values obtained from Sato et. al. experiments. ( c John Wiley and Sons taken from Sato et al. [46] with permission.) • Three b-values of 0, 300, and 600 s/mm2 were selected. After obtaining the images, a urologist drew the biopsy sites on the T 2 W images, and ROI’s were then specified by taking into account landmarks such as the urethra, limit between the peripheral zone and the transition zone, and capsule for that purpose. The placement of the ROI’s on the DTI images was manually done by a radiologist and mean values of ADC and FA were calculated. Figure 5.2 summarizes the results. This work was published in 2005 and all the details can be found in reference [46]. 84 5.3 Pickles’ Group Study Pickles mentions how, by the date of his publication, no 3 T studies had been performed in the prostate cancer diagnosis using diffusion data. They recruited 53 patients during 2004 and selected according to elevated serum PSA levels or if histology proved PCA. Nine asymptomatic volunteers were recruited. T 2 W and data were collected using a 3 T scanner without any endorectal coil. Though DTI the details for the T 2 W images acquisition are not specified, the DTI was collected as follows: [41] • Slice thickness of 5 mm. • FOV of 26 cm. • Matrix size of 224 × 224. • Two b-values of 0 and 500 s/mm2 were selected. The manufacturer’s software was used and ROI’s were defined around tumours and NPZ for the patient images, while NPZ and CG were selected for the volunteers. The ROI’s were drawn by consensus of a doctor, two researchers, and a urologist on the T 2 W images, taking into account the radiology report for the examination. The mean ADC values for the different regions were calculated and the results are summarized in figure 5.3.[41] Study by Pickles et. al. study was published in 2006 and the complete details can be found in reference [41]. 5.4 Manenti’s Group Study Manenti et. al. present a in which DTI is performed at 3T to discriminate between normal and cancerous tissue in the prostate gland. The increase in SNR at 3T improved the data quality in this study. His study involved 30 patients with elevated levels of PSA (Section 2.2) and biopsy proven cancer. The tumours identified in this study had an average Gleason score of 7, ranging from 5-9. [31] The setup of his experiments can be summarized as follows: 85 Figure 5.3: ADC values for tumour and normal tissue in Pickles et. al. study. ( c John Wiley and Sons taken from Pickles et al. [41] with permission.) • Slice thickness of 3 mm with 0.3 mm gaps were selected when imaging. • FOV of 160 mm. • Matrix size of 720 × 720 for T 2 W images. • Matrix size of 80 × 120 for DTI images. • Two b-values (equation 3.52) of 0 and 1000 s/mm2 were selected. After an average of 10.4 days after the MRI procedure was done, the patients underwent radical prostatectomy (section 2.3), their prostate specimens were fixed in 10% formalin, and sliced into 3 mm thick slices trying to fit the plane corresponding to the T 2 W images. Histology slices, 4µm thick were analyzed by a pathologist with no knowledge about the MRI results, and the tumour edges were drawn on the slices. The tumours were finally graded using the Gleason score system. After digitizing the slices using a flatbed scanner and using a DTI analysis software a pathologist traced an ROI on the b = 0 images correlating the histology slides with the MRI slices. Another two identical ROI’s were also selected; one in the peripheral zone and another one in the central gland. Mean values for FA and ADC were calculated for the different ROI’s (section 3.7.3). Figures 5.4a and 5.4b summarize the results of this study. The results of this study show unusually high values of FA in NPZ and CG compared to other papers published on the subject, such as the ones detailed in the following sections. This behavior is hard to explain 86 (a) ADC values in Cancer, Peripheral Zone, and Central Gland selected ROI’s (b) FA values in Cancer, Peripheral Zone, and Central Gland selected ROI’s Figure 5.4: DTI results from Manenti’s experiments. ( c Wolters Kluwer Health taken from Manenti et al. [31] with permission.) based on the prostate gland morphology. [31] Manenti however, did not correlate the ADC or FA values with the Gleason score of the different tumours of his study. The patients were recruited between March and November 2005 and his work was published in the year 2007. The full paper is found in reference [31] 5.5 ¨ Gurses’ Group Study G¨urses’ et. al. study was performed on 30 volunteers without any history of prostate disease. The scans were performed on a Philips Achieva 3T scanner in which T 2 W and DTI images were acquired.[16] The details were given as follows: • Slice thickness of 2.7 mm without gaps. • FOV of 430 mm. • Matrix size of 243 × 512 for T 2 W images. • Matrix size of 128 × 126 for DTI images. 87 • Two b-values of 0 and 700s/mm2 were selected. The b-value was chosen based on the ADC value of the prostate. G¨urses mentions that the optimal b value in a diffusion experiment is one such that b × ADC = 1. [16] DTI data was processed using the manufacturer’s software. The ROI’s were drawn on the b = 0 images, following the T 2 W images as an anatomical guide. Four ROI’s of identical size were selected; CG left, CG right, NPZ left, and NPZ right. FA and ADC mean values for the selected regions were calculated. G¨urses concludes that statistically significant differences are observed between the central and peripheral zones. The ADC being lower in the central gland than in the peripheral zone, and the FA value of the peripheral zone was significantly lower than the one in the central zone. Table 5.2 shows the mean values obtained by this group.[16] ADC (×10− 3mm2 /s) FA 1.220 ± 0.271 1.610 ± 0.347 0.260 ± 0.060 0.160 ± 0.030 Central Zone Peripheral Zone Table 5.2: Results from G¨urses’ et. al. experiments. ( c Springer taken from G¨urses et al. [16] with permission.) G¨urses’ study was published in 2008 and the complete details can be found in reference [16]. 5.6 Xu’s Group Study Xu states that for the first time ever the diffusion characteristics of histologically defined PCA were determined both in-vivo and ex-vivo for the same patient. This group analyzed 24 patients, all of which were scheduled to undergo a radical prostatectomy procedure (section 2.3). Any preoperative treatment such as radiotherapy was excluded from the study. Fourteen prostate glands were examined, both in-vivo and ex-vivo, while ten of them were studied only ex-vivo after the 88 extraction of the specimen. The study was then divided into three steps, one for the in-vivo scan, another for the ex-vivo case, and a histology image coregistration. [56] 5.6.1 In-vivo MRI A 1.5 T Siemens scanner was used for the scans. Both DWI and T 2 W images were generated using the details that follow: [56] • 2.5 mm thick transverse slices. • FOV of 256 × 256 mm2 . • Matrix size of 128 × 128 for DTI images. • Matrix size of 256 × 256 for T 2 W images. • Two b-values of 0 and 500s/mm2 . The DWI were post processed to correct for any eddy currents and motion.[56] 5.6.2 Ex-vivo After the radical prostatectomy procedure, the prostate specimens were fixed in 10% formalin for a minimum time of 48 hours. The specimens were then sliced from base to apex at 4mm intervals (slice thickness) with a constructed prostate slicer built at their laboratory. In order to prevent tissue dehydration, the slices were then placed in formalin again. A 4.7 T scanner was then used to image the ex-vivo slices. The settings used to acquire this images were:[56] • Slice thickness of 0.5 mm. • FOV of 6.5cm × 6.5cm • Matrix size of 128 × 128 for DTI images. • Matrix size of 128 × 128 for T 2 W images. • Two b-values of 45 and 1130s/mm2 . 89 5.6.3 MRI and Histology Coregistration By using the software ImageJ the coordinates for the ex-vivo MRI were corrected to match the histology image coordinates using a set of control points. A 3D rigid body registration was then performed with 9 degrees of freedom. The ex-vivo ADC and FA maps were also modified to the coordinates of the corresponding in-vivo T2W images using manual alignment of intraglandular structures provided by the ITK module in the Analyze software. A final 12 degrees of freedom affine transformation was applied to coregister the in-vivo FA and ADC maps with the T 2 W images. This concluded in a match between the histology slides and the diffusion maps with the standard in-vivo T 2 W images.[56] Having registered both ex-vivo and in-vivo ADC and FA maps and the histology slides to the in-vivo T 2 W images, regions of PCA, NPZ, stromal and epithelial BPH (see section 2.1) were selected, and the mean values of the diffusion parameters were calculated. Figure 5.5a shows the results obtained by Xu for the ADC values in both ex-vivo and in-vivo, while figure 5.5b shows the FA values. [56] The work from this group was published in 2009 and the details can be found in reference [56]. 5.7 Kozlowski’s Group Study Kozlowski et. al. present a study in which DTI and Dynamic Contrast Enhancement (DCE) techniques are combined quantitatively to determine if the specificity and sensitivity in detecting PCA is improved. Twenty-five patients were recruited due to its elevated levels of PSA serum and/or palpable nodules, which had not received previous treatment. Biopsies were performed under local anesthetic and its number depended on the volume of the gland. MRI data was acquired on a 3 T Philips scanner combining an endorectal coil and a cardiac phased-array coil. A FSE sequence in order to generate T 2 W images and an EPI diffusion one were used to collect both anatomical and diffusion data. The details of the setup are as follows: [27] • Slice thickness of 4 mm with no gaps. 90 (a) Ex-vivo (top) and invivo (bottom) ADC values obtained by Xu (b) Ex-vivo (top) and invivo (bottom) FA values obtained by Xu Figure 5.5: Ex-vivo and in-vivo results from Xu et. al. experiments. ( c John Wiley and Sons taken from Xu et al. [56] with permission.) • FOV of 14 cm for the T 2 W images. • FOV of 24 cm for the DTI images. • Matrix size of 284 × 225 for the T 2 W. • Matrix size for the DTI images is 128 × 115. • Two b-values, b = 0, 600s/mm2 . • Diffusion was measured in six different directions. Diffusion weighted images were registered to the non-weighted b = 0 image using a mutual information registration algorithm. The diffusion tensor was then calculated and maps of the ADC and FA parameters were generated. ROI’s were manually selected according to biopsy results and according to low values of ADC . The mean values of those parameters in the tumours were compared with 91 the healthy peripheral zone and central gland. Table 5.3 summarizes the obtained results.[27] ADC (×10− 3mm2 /s) FA 1.61 ± 0.12 1.76 ± 0.23 1.07 ± 0.12 0.20 ± 0.03 0.20 ± 0.06 0.20 ± 0.05 Central Gland Peripheral Zone Prostate Carcinomas Table 5.3: Results from Kozlowski’s et. al. experiment. ( c Elsevier taken from Kozlowski et al. [27] with permission.) The work from this group was published in 2010 and the details can be found in reference [27]. 92 Part II Methodology, Results, and Discussion 93 Chapter 6 Methods and Procedure 6.1 Purpose As it was shown in chapter 5, some studies have been performed in order to determine if diffusion data acquired with MRI provides information that leads to the correct diagnosis of PCA. All of the works presented from section 5.1 to section 5.6 demonstrate that the ADC values (section 3.7.3) in carcinomas are reduced in comparison with healthy tissue. Of those studies, three of them included the FA value as an interesting parameter that may provide even more information for the correct grading and identification of the cancer, but lead to contradictory results. • In 2007, Manenti suggested that FA values in the prostate carcinomas is lower than in the NPZ. At the same time, FA values of the CG tissue is lower than the one in the NPZ but higher than the one in the carcinomas.[31] • One year later, in 2008, G¨urses, though his studies were performed in healthy volunteers, his results suggested that the FA values for tissue in the CG and NPZ are the opposite of Manenti’s study; CZ having a higher value. [16] • In 2009 and 2010 Xu and Kozlowski respectively, suggested that the FA value of tissue in NPZ and PCA is basically unchanged.[27, 56] With these inconsistencies in mind, the purpose of the project is to identify and quantify the role of FA values in the different types of prostate tissue, and determine 94 if it provides reliable and useful information in the PCA grading and diagnosis. As far as it is known, a correlation between ADC and FA values with the Gleason score has not been studied before. The study was approved by the institutional human ethics board and all participants gave signed consent prior to entering the study. 6.2 6.2.1 Materials and Methods Patients Thirteen patients with biopsy proven carcinoma (section 2.3) were recruited during the period of August 2010 to October 2011. The mean age of these patients was 64 years and range 56-73 years. The prostate glands of all the patients were examined both in-vivo (before surgery) and ex-vivo (after resection). 6.2.2 In-vivo MRI The in-vivo examination was performed using a 3T Philips Achieva MR scanner with the use of a combined cardiac phased array/endorectal coil, and T 2 W and DTI images were obtained. Axial T 2 W images were obtained using a TSE pulse sequence (section 3.6), whose details are shown in table 6.1. In-vivo T 2 W Matrix Size TR FOV Slice Thickness 284 × 225 1851.190ms 140 cm 4 mm without gaps. Table 6.1: TSE sequence details for the acquisition of the in-vivo T 2 W images. 95 In-vivo DTI Matrix Size TR FOV Slice Thickness Number of directions b-values Diffusion Time 128 × 115 2100.000ms 24 cm 4 mm without gaps 6 0 and 600 s/mm2 36.4 ms Table 6.2: Diffusion EPI details for the acquisition of the in-vivo DTI data. 96 The diffusion images were obtained using an EPI diffusion sequence (section 3.7.4), and its details are shown in table 6.2. The six directions of the diffusion √ √ √ √ √ √ sensitizing gradient were: ( 2/2, 0, 2/2), (− 2/2, 0, 2/2), (0, 2/2, 2/2), √ √ √ √ √ √ (0, 2/2, − 2/2), ( 2/2, 2/2, 0), (− 2/2, 2/2, 0). 6.2.3 Ex-vivo MRI After the robotic surgery was performed, the excised glands were fixed with 10% buffered formalin for at least 24 hours. A 1.8mm diameter plastic rod was inserted through the urethra, to facilitate the alignment of the MRI slices with histology sections. The specimens were placed on a small plastic platform and covered with a cotton towel soaked in formalin as shown in figure 6.1. The specimens were placed inside a volume coil and data were collected using a 7 T Bruker Biospin 30 cm bore MRI scanner. T 2 W images were acquired using a RARE sequence, and DTI data were obtained using a spin-echo sequence with diffu- sion sensitizing gradients (section 3.7.2). The details are given in tables 6.3 and 6.4 and the diagrams of these pulse sequences, generated with the Bruker ParaVision software, are shown in figures 6.2a and 6.2b. The six directions of the diffusion √ √ √ √ √ √ sensitizing gradient were: ( 2/2, 0, 2/2), (− 2/2, 0, 2/2), (0, 2/2, 2/2), √ √ √ √ √ √ (0, 2/2, − 2/2), ( 2/2, 2/2, 0), (− 2/2, 2/2, 0). Ex-vivo T 2 W Matrix Size TR Effective TE FOV Slice Thickness 256 × 256 5000.000ms 46.98ms 6 cm 4 mm without gaps Table 6.3: RARE sequence details for the acquisition of the ex-vivo T 2 W images. Using the lack of signal due to the presence of the plastic rod in the urethra, the slice orientation was selected as close as possible to the orientation in which the specimens were then sliced, as it is described in the next section. The number of 97 (a) Specimen placed on a plastic platform. (b) The gland is covered with a cotton towel soaked in formalin (c) It was then well taped to the platform and wrapped with plastic to avoid any leakage of the formalin Figure 6.1: Experimental set up of the specimen for the collection of the ex-vivo data. 98 (a) RARE pulse sequence used to acquire the ex-vivo T 2 W images. (b) Spin-echo pulse sequence used to acquire the ex-vivo DTI data. Figure 6.2: Ex-vivo pulse sequences used to acquire the data. 99 Ex-vivo DTI Matrix Size TR TE FOV Slice Thickness Number of directions b-values Diffusion Time 128 × 128 1500.000ms 22.277ms 6 cm 4 mm without gaps 6 0 and 750 s/mm2 9.46 ms Table 6.4: Spin-echo sequence details for the acquisition of the ex-vivo DTI images. slices was chosen to cover the entire specimen. A total time of approximately 1.5 hours was required for the DTI scan. 6.2.4 Histology After the ex-vivo scan, a device [9] designed at the 7T facility of the UBC MRI Research Centre was used to section the specimens into 4 mm thick slices. Some small modifications, shown in figure 6.3, were made to the original slicer device in order to facilitate easier setup of the blades. The openings for the blades were extended to the top, and were made slightly wider, allowing easier insertion of the blades. The piece to hold the blades on the opposite side of the handle was similarly modified, by opening the bottom part of the slots and making them slightly wider to simplify the process of assembling the device. The seminal vesicles, bladder neck, and the small part of the apex were removed by the pathologist and the specimen was then sliced using the slicer device shown in figure 6.4a. Four millimeter sections, shown in figure 6.4b, were obtained and were then immersed in paraffin. Whole Mount Sections (WMS) were cut using the Leica RM2245 microtome, mounted on oversized glass slides and stained with hematoxylin and eosin (H&E) stain. 100 (a) The blade openings were extended to the top and widened at the ends. (b) The blade openings were also extended on this piece whose purpose is to hold the blades at the end opposite to the handle. Figure 6.3: Two modified parts of the prostate slicer device [9] in order to make assembling the device easier at the moment of slicing the specimen. Tumour areas were manually outlined and graded by a pathologist with the use of a microscope, and digital images were generated from the annotated glass slides using a flatbed scanner. The tumours were graded following the Gleason score grading system (section 2.3.1). 6.2.5 Histology and MRI Coregistration In order to register the histology slides to the MRI images, a software program was written in Matlab (Mathworks, Natick, MA, USA) that incorporated the ITK li101 (a) Slicer device.[9] (b) Slices obtained after the sectioning process. Figure 6.4: Four millimeter thick sections were obtained after slicing the specimen. The seminal vesicles, bladder neck, and the small part of the apex were removed before the sectioning. braries available with elastix [24]. Due to different deformations of the prostate, as well as the way in which the diffusion parameters are calculated for the in-vivo and ex-vivo data acquisition, the software was divided into two branches corresponding to either of the scanning modalities. 6.2.6 Launching the Program The main function of the software is called in ex reg, in which the user selects one of the two type of registrations available; ex-vivo to WMS or in-vivo to WMS. The displayed windows is shown in figure 6.5. 102 Figure 6.5: in ex reg function which launches the registration program. Figure 6.6: Main window for the ex-vivo to WMS registration procedure. Ex-vivo to WMS The ADC and FA maps obtained from the ex-vivo MRI scan procedure were calculated using Bruker’s provided software, Paravision 3.0.2 (Biospec R system). After selecting the ex-vivo to WMS registration, a new window (figure 6.6) opens up providing the option to “Run Registration”, which should be clicked to begin the process. The user is asked to select the folder where the ex-vivo DTI data is stored. The user is then notified of how many slices are available for the corresponding specimen (this depends on the size of each particular specimen) and the function slreg loadData opens a new window. 103 Figure 6.7: slreg loadData window with two sets of images selected. The fixed image set is made up of the T 2 W images while the WMS are the moving images. slreg loadData is a function in which the user loads the two sets of images which will be registered. As anatomical details are better seen in T 2 W images and the diffusion maps overlap precisely with these images, the histology slides are selected as the moving images and the T 2 W ones as the fixed ones (chapter 4). Figure 6.7 shows how the window generated by slreg loadData looks like. The user has three options to load the images: • Load Volume: Loads a binary file that contains the information of different images in the scanned volume. When the volume is loaded, the user must specify the number of rows and columns, the number of images per slice, and the number of slices in the binary file, as well as the byte order (little endian or big endian). • Add Folder: Loads all the images located in a folder. The images’ extension can be jpeg, jpg, tiff, tif, or png. • Add Image: Loads only one image with extension jpeg, jpg, tiff, tif, or png. The user selects the file from a new menu. The next function is called slreg pairData (figure 6.8) where the user must pair the fixed image that corresponds to a moving image. Clicking on the “Show ADC 104 Figure 6.8: slreg pairData windows with two paired slices. Image” button displays the ADC map of the current T 2 W image in order to help with the pairing process. After pairing the slices and clicking continue, a new window slreg segment (figure 6.9), opens up and the user is required to segment all the paired slices. This is done to make the registration algorithm work better, as only the desired information of each image is given as an input. This is done by manually defining an ROI around the prostate specimen. After having segmented the prostate specimens, the registration algorithm can be run. The function that runs it is called slreg register (figure 6.10) in which an initial and a final set of parameters are required. In order to perform an initial alignment of the images, an affine transformation is carried out, whose parameters are given in the “initial parameter” file and can be found in appendix A. The parameter files have been specified according to the elastix manual [24] and follow a 2D, multiresolution algorithm with a B-spline interpolator. Though it follows the mutual information theory described in section 4.2.3, the metric, called in elastix as the advanced mattes mutual information, is given by equation 4.40. The “final parameters” specify a second non-rigid b-spline transformation (section 4.1.5), composed of the affine transformation described before. Though the parameter 105 Figure 6.9: slreg segment allows segmentation of the specimen in which non useful data is removed in order to obtain a better registration. Figure 6.10: slreg register runs the registration algorithm by calling the elastix [24] software with the parameter files specified in appendix A. files contain more information than the one explained in chapter 4, they are explained in the elastix software manual [24]. After clicking on the “Register Slices” button the whole algorithm to register the sets of paired images is run. Following the registration, different ROI’s can be selected, and their diffusion 106 Figure 6.11: Information displayed at the end of the ex-vivo registration process. parameters calculated. The number of ROI’s in each of the registered slices needs to be set before clicking on the “Draw ROIs” button. A window in which the different ROI’s are specified opens and the user outlines them manually. Finally, the software displays the registered images in a new window. The user must manually select the contour of the registered prostate for each slice. These contours are used to determine the quality of the registration by calculating the Dice Similarity Coefficient (DSC) value, explained in detail in the next section. When the process is concluded, the ex vivo registration window displays the slices, through which the user can browse using the buttons at the bottom, and for which the ADC and FA values, as well as the DSC value are displayed. This can be seen in figure 6.11. It has to be clear that though the registration involved the T 2 W images, they are already registered with the FA and ADC maps as they are collected without moving the specimen. In-vivo to WMS Unlike the ex-vivo case, the in-vivo ADC and FA maps are calculated in Matlab using the developed software. Due to low resolution provided by the ADC and FA 107 Figure 6.12: In-vivo registration procedure main window. maps, the WMS are registered to the T 2 W images. The procedure is repeated to register the diffusion maps to the T 2 W images selected before. In order to leave the in-vivo ADC and FA maps unchanged, the inverse transformation of the registration, which aligned the diffusion images with the T 2 W images, is applied to the registered histology images, and WMS registered to the original diffusion parameter maps are obtained. When the user selects the in-vivo registration from the main menu, the window shown in figure 6.12, will open up. From this window, the user can select the following functions: DTI reconstruction, T 2 W-Histology registration, and T 2 W-DTI registration. The first step is to generate the diffusion maps from the data obtained from the 3 T MRI scanner. After selecting the DTI reconstruction, the window shown in figure 6.13 will open up. The function in charge of this is called dtiUI, and the binary data has to be loaded from the file menu. A scaling factor has to be applied due to the way Philips stores the data. Values of the scaling factors are provided in table 6.5. 108 Figure 6.13: dtiUI main window in which the DTI reconstruction is made and diffusion maps are generated. A floating point FP, which is the correct value to use for the calculations, depends on the displayed value DV in the console, a rescale slope RS, and a scale slope SS as given by equation 6.1. FP = DV RS · SS (6.1) The displayed value is given by equation 6.2 where PV represents the pixel value and is the one stored in the Philips data file, the so called REC file, RI is the rescale intercept. DV = PV × RS × RI (6.2) After inserting 6.2 in equation 6.1 in which RI = 0 for the cases treated here FP = PV SS (6.3) When the data is loaded, it is important to load the parameter file provided by Philips, so that the software knows the magnitudes and directions of the diffusion sensitizing gradients, in order to diagonalize the diffusion tensor (section 3.7.2). This can be done in the menu “File/set DTI parameters” . After loading the 109 data and parameter files, the software calculates the tensor components as given by equation 3.55. The user can then generate the ADC and FA parametric maps (section 3.7.3). The generated parametric maps can be saved by “Save Current Data” from the “File” menu. An appropriate value of the scaling factor, which needs to be applied to the parametric maps at the time of saving, has to be input by the user (see table 6.5). Note that every time the file is loaded in any part of the software, in order to obtain the correct data, the numbers have to be divided by the corresponding value of the table 6.5 as well. When the ADC map is loaded and scaled by the factor shown in the table, the numbers obtained will be in units of 10−3 mm2 /s. DTI Original File Eigenvalue FA map ADC map Load (÷) Save (×) 1.976 × 10−2 4.095 × 102 4.095 × 103 4.095 × 102 1.976 × 10−2 4.095 × 105 4.095 × 103 4.095 × 105 Table 6.5: Scaling factors required to load and save the different maps of the DTI reconstruction function dtiUI Having obtained the diffusion parametric maps, the next step is to run the histology WMS to T 2 W registration procedure. This process runs exactly the same as for the ex-vivo registration; however the user is not required to select the contour of the prostate specimens at the end of the process. The parameter files for the transformations are slightly different as the endorectal coil deforms the prostate gland in a much bigger way than it is deformed during the ex-vivo scan (which is basically negligible). These parameters can be found in appendix A. The registration of DTI to T 2 W images is almost the same as in the ex-vivo case, however after the registration is completed, the generated ADC and FA parameter maps have to be paired to the corresponding registered T 2 W images. Following that, the already registered WMS are inverse transformed and aligned to the original (unregistered) parametric maps. The way this is done was proposed and evaluated by Metz et. al. [36] and explained in the elastix manual [24]. It consists 110 Figure 6.14: dtiGUI Final window of the in-vivo registration process in which the mean values of the diffusion parameters are shown at the left. of changing the mutual information metric for the displacement magnitude penalty shown in equation 6.4 where Tµ is the transformation wanted to be inverted and that has parameters specified in the vector µ for each point x. ||Tµ (x) − x||2 (6.4) The inversion is performed as another registration problem, in which the transformation to be inverted is used as the initial transform but elastix is run with the magnitude penalty metric. The images, both fixed and moving, are set to be the initial fixed images which for this case are the T 2 W images. [24] This generates the parameters that are applied to the registered histology images. The WMS are finally registered to the original parametric diffusion maps, after which the DTIGUI window opens up (see figure 6.14) and the average ADC and FA values for different regions are shown on the left side of the window. A diagram of the in-vivo registration procedure is shown in figure 6.15. 111 Figure 6.15: Registration of in-vivo DTI to histology. Both histology (left) and in-vivo DTI (right) are registered to an in-vivo T2W image (top). The inverse of the transformation found for DTI registration is applied to histology and resized in order to match the original ADC and FA maps obtained from the DTI data. Evaluation of the Registration In order to determine if the registration software was working as required, 34 registered slices for both in-vivo and ex-vivo cases (68 in total), were studied to quantitatively determine if the software was reliable. This validation was based on the alignment of the contour and the correspondence of the anatomical details in the interior of the gland explained as follows: • Contour Alignment: Overlap is usually measured by using the DSC shown in equation 6.5 in which X and Y represent the fixed and the registered image respectively. The numerator gives the cardinality (number of pixels) in the intersection of both images and the denominator is the sum of cardinalities of the images independently. The closer the resulting value is to unity the better the overlap (registration) of the images.[24] 112 DSC = 2|X ∩Y | |X| + |Y | (6.5) • Anatomical Details: The internal anatomical details of the registered images were qualitatively and a quantitatively analyzed. An overlay of the registered and fixed images, where one image is displayed on top of the other, with the possibility of changing the transparency of the top image, was used. This is available by pressing the “Show Overlay” button in the registration software. The quantitative evaluation involved selecting corresponding landmarks in both images and, as the coordinates are known, determining the distance between them. In addition, similar ROI’s to the ones annotated by the pathologist were manually selected on the ADC maps, based on low ADC values (5), and the average ADC and FA values were calculated and compared to those from the registered ROI ’s. 113 Chapter 7 Results and Conclusions 7.1 Results The results presented here include data obtained for four of the mentioned patients as histology had not been received for the remaining ones at the moment of writing this document. 7.1.1 Validation of the Registration Software Thirty-four slices from the four data sets were used for the evaluation of the registration software mentioned in the last section. Figure 7.1 shows an example of a registered slice for the ex-vivo case. The matching of the contour of the images is delineated in green and red. The bottom of the image shows the overlay explained in the last section, with a transparency level change in order to evaluate the registration method qualitatively. The ROI’s of the desired areas selected in this study included tumours with Gleason Score of 3 + 3 (GS 33) and Gleason Score of 3 + 4 (GS 34), NPZ, and Normal Peripheral Zone With Enlarged Glands (NPZEG). Figure 7.2 shows a histology slide in which the regions have been marked. A mean DSC of 0.97 ± 0.01, 0.950.04, and 0.940.05 for ex-vivo DTI to WMS, in-vivo WMS to T 2 W, and in-vivo DTI to T 2 W was obtained respectively for the 34 slices. Five slices showed a poor match of the fine anatomical details in the overlay and were not included in this analysis. Percentage differences between the average 114 Figure 7.1: Registration of ex-vivo DTI (top left) with histology (top middle) shows a very good overlap measured with the DSC value (top right) in which the red and green contours represent the fixed and registered images respectively. On the bottom, an overlay of the original fixed image with the histology shown on top of it; the level of transparency is changed in steps of 0.2 from 0 (DTI image bottom left) to 1 (registered histology bottom right). Figure 7.2: A histology slice with outlined regions of interest. values from ROIs selected manually and with the registration algorithm are shown in table 7.1. Sixty-five and fifty landmarks were selected for the in-vivo and ex-vivo registration respectively. The distribution of the distances between landmarks showed a different behavior than the one of a normal distribution, and thus the median value was used instead of the mean. Figure 7.3 shows a plot of the results in which a 115 % diff ADC ex-vivo % diff FA ex-vivo % diff ADC in-vivo % diff FA in-vivo GS 33 GS 34 NPZ NPZEG Mean NPZ and NPZEG 2.4 4.12 13.9 11.6 8.2 5.2 10.4 10.1 1.1 1.1 2.7 3.3 0.1 5.1 4.4 8.1 0.7 1.6 1.1 6.3 Table 7.1: Percentage differences obtained between the averages of ADC and FA on manually selected regions of interest in comparison with the ones generated by the registration algorithm. Five slices showed a poor match of the fine anatomical details in the overlay and were excluded from this analysis. Figure 7.3: Median values of distance between corresponding landmarks in the registered and fixed images for both ex-vivo and in-vivo. The data is displayed as median value and the bars correspond to the range of values. median distance of 0.66 mm (range 0.23, 2.90 mm) and 1.55 mm (range 0.3 to 3.1 mm) was obtained for the ex-vivo and in-vivo registrations respectively. 116 Figure 7.4: Ex-vivo (left) and in-vivo (right) ADC and FA values for the different types of tissue in the patient’s prostate gland. 7.1.2 Values of the Diffusion Parameters for the Different Types of Tissue in the Prostate The 34 registered slices were evaluated to determine the ADC and FA values of the GS 33, GS 34, NPZ , NPZEG . A mean value corresponding to the peripheral zone in which no distinction was made between normal and enlarged glands was also calculated as the average of the NPZ and NPZEG regions of the gland. Central gland tumours were not taken into account in this procedure. Figure 7.4 shows the distribution of the results for both ex-vivo and in-vivo. The results are presented as mean ± standard deviation. No statistical significance test has been done yet due to the current low number of patients. 117 7.1.3 Discussion Evaluation of the Registration Software The evaluation of the registration software suggests that, an accurate registration between DTI and WMS has been achieved in this study. Differences in the values obtained from the manual selection of similar ROI’s annotated by the pathologist and the ones marked with the registration software are believed to be caused by tumour heterogeneity and cancer infiltrations that cannot be detected in MRI but that are observed through the microscope and included in the annotated tumour region. Higher differences between the manually selected ROI’s and the annotated ones invivo are explained by the lower resolution provided in this images compared to the ex-vivo case. The ADC and FA Parameters Based on the presented results it has been seen that the ADC values in tumours is lower, when compared to the average peripheral zone in both the in-vivo and the ex-vivo case. This lower values show that the mean square distance travelled by water molecules is shorter than in normal tissue. Tumour tissue cells reduce the water volume fraction and generate randomly oriented restrictions for the molecules to move [41]. Figure 7.5 shows the micro structure of healthy peripheral zone and PCA . The high cellularity of the prostate carcinomas that cause the restrictions to the water molecules act on a distribution of length scales on a range that goes from sub-microns to tens of microns. DTI is sensitive to those length scales. Cancer infiltrations does not change the cellularity of the prostate and ADC changes d not necessarily occur.[56] Ex-vivo ADC values are seen to be lower than in the in-vivo cases, but follow a similar trend. Diffusion is highly dependent on different factors, for example temperature. The ex-vivo scan is done at room temperature which is in fact around 17◦C lower than the in-vivo scan and this behavior was expected. The fixation process is also involved in the reduction of the ADC values as compared to the in-vivo 118 (a) Healthy Peripheral Zone (b) Prostate Cancer Figure 7.5: Magnified images of healthy peripheral zone and PCA.( c John Wiley and Sons taken from Xu et al. [56] with permission.) case. The ADC value of the NPZEG is higher than the one in the NPZ. The enlarged glands contain less barriers to diffusion compared to NPZ making the value larger. In the ex-vivo case, the NPZ values are similar to the tumour ones, in which dehydration effects and fixation might be playing an important role; this is not the case in-vivo. Partial volume effects present in-vivo might also be present in the reduction of this heterogeneity of the ADC values in the peripheral zone. As mentioned in chapter 2, the human prostate gland is formed by a mixture of fibromuscular tissue with glandular and non-glandular components. The fibromuscular tissue is what causes the observed diffusion anisotropy within the gland. When determining the FA values for the different prostate zones and PCA, ex-vivo FA values were found to be slightly higher in the carcinomas compared to the pe- ripheral zone, while the opposite effect was observed in-vivo. To explain this discrepancy between the in-vivo and ex-vivo trend of FA values in the NPZ and PCA, the sources of errors in estimating FA in DTI can be evaluated. Figure 7.6 shows a Monte-Carlo simulation performed by Reinsberg et. al. [44] where it can be seen that the measured value of FA rises as the noise increases. The 119 Figure 7.6: Monte-Carlo simulations have been made in order to study the behavior of the measured value of FA and its dependence on noise. ( c ISMRM taken from Reinsberg et al. [44] with permission.) healthy peripheral zone seems to be showing this effect as FA values are larger for the noisier in-vivo scans. Monte-Carlo simulations similar to the ones presented by Reinsberg [44] were carried out to explain the FA values in the tumour. Figure 7.7 shows the behavior of FA as a function of the ADC value for different realistic values of FA and SNR. The simulations suggest that as the value of ADC is lowered, FA increases, while as the SNR increases the FA value gets smaller. As the real value of FA grows, the different curves seem to get closer and the effect of the SNR is reduced. The plots do not show a big effect of the b-value (section 3.7.2) on the FA behavior. [44] Due to the tissue structure of the prostate it is hard imagine that the real FA value can be equal to 0, i.e. that diffusion can be fully isotropic. It was concluded that though there might be an effect of FA depending on the value of ADC, for the acquired range of ADC values in this study this effect was negligible. Lower tumour FA values may also result from partial volume effects. The spatial resolution of the in-vivo DTI data was 5 × 6.2 × 4mm3 compared to the ex-vivo case, where the spatial resolution was 0.23 × 0.23 × 4mm3 . The study shows very 120 (a) Real FA = 0 (b) Real FA = 0.1 (c) Real FA = 0.2 (d) Real FA = 0.3 Figure 7.7: Monte-Carlo simulations of the behavior of FA as a function of ADC values for different realistic values of FA and SNR . The scale has been kept the same for all the cases. similar results to Xu et.al. work [56] (figure 5.5), described in section 5.6. They also concluded that the discrepancy between in-vivo and ex-vivo is very likely due to the partial volume effects. Our preliminary results show that the proposed registration algorithm provides accurate registration between DTI and histology data and may be useful in validating the diagnostic accuracy of the DTI technique. Though probably too soon to state a definite conclusion, the results shown suggest that there are no significant differences in FA between normal peripheral zone and prostatic carcinoma. These results together with the results of the studies that show the FA values being much less repeatable than the ADC values [12], strongly suggest that FA is not likely to 121 (a) Leica scn4000 scanner. (b) Digital image generated by the Leica scn400 scanner. Figure 7.8: A new WMS digitization scanner will be available soon. ( c Leica Microsystems taken from Leica Microsystems [28] with permission). contribute significantly to diagnostic capabilities of DTI in prostate cancer. 7.2 7.2.1 Future Studies and Suggestions Digitizing the WMS Histology of the remaining specimens still need to be analyzed. A new technology for the digitization of the WMS will be available soon in the Vancouver Prostate Centre: a new scanner, the Leica scn400 (Leica Microsystems CMS GmbH, Wetzlar, Germany) shown in figure 7.8 will be used to digitize the histology slides. Among the outstanding characteristics that the new technology will bring, whose details are shown in figure 7.9, is the short time required to scan the slide at up to a 40× magnification.[28] The manufacturer provides a software that can be used to build a database of the scanned slices. Through an online server, the pathologist will be able to annotate any slices, marking all the desired ROIs and commenting on them. This will be done from any internet browser the pathologist is registered to as an authorized user. Digital magnification up to twice the optical one can be achieved by 122 Figure 7.9: Leica scn400 slide scanner specifications. ( c Leica Microsystems taken from Leica Microsystems [28] with permission). the software; in addition the areas of the selected ROIs are provided as well. The technology will bring better quality images, with more accurately specified ROIs and easy access to relevant information.[28] 7.2.2 Registration Software The most current elastix version, released after completing the registration software, now includes an option to input corresponding points in the images in order to help the registration algorithm. The software then will rely not only on the intensity values but will also take into account positions of known landmarks. This is a way of selecting the anatomical landmarks described in section 6.2.6 before running the algorithm. This type of option can be added to the software to solve 123 cases where the registration is not satisfying. [23] Some of the windows that the user encounters when running the registration software might be given the option to go back to the previous function without losing the data that has been stored until that moment. As an example, currently a user cannot go back to select another slice or another pair of slices without having to start the whole process again. 7.2.3 Correspondence of WMS and MRI Though hard work has been done in order to reliably slice the prostate specimen in such a way that the slides correspond to the locations of the corresponding MRI slices, this can still be improved. It has been suggested that immersing the prostate specimen in a gel, that solidifies after a while, before the ex-vivo scanning, will make it easier to position the gland in the scanner in the same position as it will be sliced in the hospital. The procedure only works for the ex-vivo case as for in-vivo the problem is almost impossible to solve. 7.3 Conclusions Several groups of researchers have been studying the diffusion properties of the cancerous prostate gland, and DTI has been successfully applied in prostate cancer diagnosis [27]. It has been well established that the water ADC has a lower value in tumours compared to normal prostatic tissue. However, FA values in prostatic carcinoma have been reported to be higher [16], lower [31], and unchanged [56] when compared to the prostate’s normal peripheral zone. The purpose of this project was to to identify and quantify the role of FA values in the different types of prostate tissue, and determine if it provides reliable and useful information in the PCA grading and diagnosis. Histology WMS were registered to both in-vivo 3T and ex-vivo 7T T 2 W images, and the DTI parameters were calculated for the tumours, normal peripheral zone, and normal peripheral zone with enlarged glands. The results shown suggest that there are no significant differences in FA between normal peripheral zone and prostatic carcinoma, strongly 124 suggesting that FA is not likely to contribute significantly to diagnostic capabilities of DTI in prostate cancer. 125 Bibliography [1] M. Amin, A. Khalid, N. Tazeen, and M. Yasoob. Zonal anatomy of prostate. Annals of KEMU, 16(3):138–142, 2010. → pages 8 [2] I. Bashar. In vivo measurement of the apparent diffusion coefficient in normal and malignant prostatic tissues using echo-planar imaging. Magn. Reson. Imaging, 16:196–200, 2002. → pages iii, 82, 83 [3] C. Beaulieu. The basis of anisotropic water diffusion in the nervous system a technical review. NMR in Biomedicine, 15:435–455, 2002. → pages 47, 48 [4] M. A. Bernstein, K. F. King, and X. J. Zhou. Handbook of MRI Pulse Sequences. Elsevier Academic Press, 1 edition, 2004. → pages 46, 47, 55 [5] F. Bloch. Nuclear induction. Phys. Rev., 70:460–474, 1946. → pages 19 [6] F. Bloch, W. Hansen, and M. Packard. Nuclear induction. Phys. Rev., 69: 127, 1946. → pages 19 [7] R. A. Bok and E. J. Small. Bloodborne biomolecular markers in prostate cancer development and progression. Nature Reviews Cancer, 2:918–926, 2002. → pages 12, 13 [8] F. V. Coakley and H. Hricak. Radiologic anatomy of the prostate gland: A clinical approach. Radiologic Clinics of North America, 38(1):15–30, 2000. → pages 7 [9] B. Drew, E. C. Jones, S. Reinsberg, A. C. Yung, S. L. Goldenberg, and P. Kozlowski. Device for sectioning prostatectomy specimens to facilitate comparison between histology and in vivo mri. Magn. Reson. Imaging, 32: 992–996, 2010. → pages 100, 101, 102 [10] L. Egevad, F. Montorsi, and R. Montironi. Obituary donald f. gleason, 1920-2008. European Urology, 55(5):1247–1249, 2009. → pages 17 126 [11] C. S. Foster and D. G. Bostwick. Pathology of the Prostate. W.B. Saunders Company, 1 edition, 1998. → pages 6, 9, 10, 15, 16, 17, 18 [12] P. Gibbs, M. D. Pickles, and L. W. Turnbull. Repeatability of echo-planar-based diffusion measurements of the human prostate at 3 t. Magn. Reson. Imaging, 25:1423–1429, 2007. → pages 121 [13] D. Gleason. Classification of prostatic carcinomas. Cancer Chemother Rep, 50:125–128, 1966. → pages 15 [14] S. Goldenberg and C. Glegg. The intelligent patient guide to prostate cancer: all you need to know to take an active part in your treatment. Intelligent patient guide. Intelligent Patient Guide, 1997. ISBN 9780969612537. → pages 6, 8, 9, 11, 12, 14 [15] A. Grove and E. Lopez-Gunn. Uncertainty in climate change (wp), 2010. URL http://www.realinstitutoelcano.org/wps/portal/rielcano eng/Content? WCM% GLOBAL CONTEXT=/elcano/elcano in/zonas in/dt25-2010. → pages [16] B. G¨urses, N. Kabakci, A. Kovanlikaya, Z. Firat, A. Bayram, A. M. Uluo, and I. Kovanlikaya. Diffusion tensor imaging of the normal prostate at 3 tesla. Eur Radiol, 18:716–721, 2008. → pages iii, 81, 87, 88, 94, 124 [17] E. M. Haacke, R. W. Brown, M. R. Thompson, and R. Venkatesan. Magnetic Resonance Imaging: Physical Principles and Sequence Design. Wiley-Liss A John Wiley & sons, Inc, Publication, 1 edition, 1999. → pages 41 [18] E. Hahn. Spin echoes. Phys. Rev., 80(4):580–594, 1950. → pages 31 [19] K. H. Hammerich, G. E. Ayala, and T. M. Wheeler. Anatomy of the prostate gland and surgical pathology of prostate cancer. Prostate Cancer. Cambridge University Press, 2008. ISBN 9780521887045. → pages 5, 6, 8, 9, 11 [20] D. L. Hill, C. Studhome, and D. J. Hawkes. Voxel similarity measures for automated image registration. SPIE, 2359(205):205–216, 1994. → pages 77 [21] P. B. Kingsley. Introduction to diffusion tensor imaging mathematics: Part i. tensors, rotations, and eigenvectors. Concepts in Magnetic Resonance Part A, 28A:155179, 2006. → pages 47, 52, 53, 54 [22] P. B. Kingsley. Introduction to diffusion tensor imaging mathematics: Part i. anisotropy, diffusion-weighting factors, and gradient encoding schemes. Concepts in Magnetic Resonance Part A, 28A:155179, 2006. → pages 127 [23] S. Klein and M. Staring. elastix, the Manual, September 2011. → pages 57, 62, 124 [24] S. Klein, M. Staring, K. Murphy, M. A. Viergever, and J. P. W. Pluim. elastix: a toolbox for intensity based medical image registration. IEEE Transactions on Medical Imaging, 29(1):196–205, 2010. → pages 102, 105, 106, 110, 111, 112 [25] P. Kozlowski. Magnetic Resonance Imaging Principles in Medical Imaging: Principles, Detectors, and Electronics. John Wiley and Sons Inc., 2009. ISBN 0470391642. → pages [26] P. Kozlowski. Ubc physics 540: Image reconstruction, 2010. → pages 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 41, 42, 43, 44, 45 [27] P. Kozlowski, S. D. Chang, R. Meng, B. M¨adler, R. Bell, E. C. Jones, and S. L. Goldenberg. Combined prostate diffusion tensor imaging and dynamic contrast enhanced mri at 3t - quantitative correlation with biopsy. Magn Reson Imaging, 28:621–628, 2010. → pages iii, 90, 92, 94, 124 [28] Leica Microsystems. Leica SCN400: The fast, reliable and flexible way to scan and digitize your slides, September 2011. → pages 122, 123 [29] A. L. MacKay. Diffusion and magnetic resonance: For ubc physics 542, 2010. → pages 48, 49, 50, 51, 52, 53 [30] F. Maes, A. Collingnon, D. Vandermeulen, G. Marchal, and P. Suetens. Multimodality image registration by maximization of mutual information. IEEE Trans. Med. Imaging, 16:187–198, 1997. → pages 78 [31] G. Manenti, M. Carlani, S. Mancino, vittorio Colangelo, M. D. Roma, E. Squillaci, and G. Simonetti. Diffusion tensor magnetic resonance imaging of prostate cancer. Investigative Radiology, 42(6):412–419, 2007. → pages iii, 81, 82, 85, 87, 94, 124 [32] A. M. D. Marzo, E. A. Platz, S. Sutcliffe, J. Xu, and H. Gronberg. Radiologic anatomy of the prostate gland: A clinical approach. Nature Reviews Cancer, 7(4):256–269, 2007. → pages 10 [33] J. E. McNeal. Anatomy of the prostate:an historicalsurvey of diverget views. The Prostate, 1(1):3–13, 1980. → pages 8 [34] J. E. McNeal. Normal and pathologic anatomy of prostate. Urology, 17(3): 11–16, 1981. → pages 128 [35] J. E. McNeal. Normal histology of the prosate. The American Journal of Surgical Pathology, 12(8):619–633, 1988. → pages 8 [36] C. Metz, S. Klein, M. Schaap, T. van walsum, and W. Niessen. nonrigid registration of dynamic medical imaging data using nd+t b-splines and a groupwise optimization approach. Medical Image Analysis. → pages 110 [37] H. Miao, H. Fukatsu, and T. Ishigaki. Prostate cancer detection with 3t-mri comparison of diffusion-weighted imaging. Eur J radiol, 61:297–302, 2007. → pages 81 [38] NCI. Digital rectal exam, Aug 2008. URL http://visualsonline.cancer.gov/details.cfm?imageid=7136. → pages 12 [39] M. Neuman. Medical Image Registration. CRC Press, 2001. ISBN 0849300649. → pages 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 78, 79, 80 [40] W. H. Organization. The Global Burden of Disease. WHO Library Cataloguing, 2004. ISBN 9789241563710. → pages 1 [41] M. D. Pickles, P. Gibbs, M. Sreenivas, and L. W. Turnbull. Diffusion-weighted imaging of normal and malignant prostate tissue at 3.0t. Magn. Reson. Imaging, 23:130–134, 2005. → pages iii, 85, 86, 118 [42] E. Purcell, H. Torrey, and R. Pound. Resonance absorption by nuclear magnetic moments in a solid. Phys. Rev., 69:37–38, 1946. → pages 19 [43] I. Rabi, J. R. Zacharias, S. Millman, and P. Kusch. A new method of measuring nuclear magnetic moment (letter). Phys. Rev., 53:318, 1938. → pages 19 [44] S. A. Reinsberg, J. M. Brewester, G. S. Payne, M. O. Leach, and N. M. deSouza. Anisotropic diffusion in prostate cancer: Fact of artefact? volume 13, page 269. Intl. Soc. Mag. Reson. Med., 2005. → pages 119, 120 [45] D. C. Rizzo. Fundamentals of Anatomy and Physiology. Delmar Publishers Inc, 3 edition, 2009. → pages 5, 6 [46] C. Sato, S. Naganawa, T. Nakamura, H. Kumada, S. Miura, O. Takizawa, and T. Ishigaki. Differentiation of noncancerous tissue and cancer lesions by apparent diffusion coefficient values in transition and peripheral zones of the prostate. Magn. Reson. Imaging, 21:258–262, 2005. → pages iii, 83, 84 129 [47] S. Sinha and U. sinha. In vivo diffusion tensor imaging of the human prostate. Magn Reson Med, 52:530–537, 2004. → pages 82 [48] C. C. Society. Prostate Cancer Statistics at a Glance. 2011. URL http://www.cancer.ca/Canada-wide/About%20cancer/Cancer%20statistics/ Stats%20at%20a%20glance/Prostate%20cancer.aspx?sc lang=en. → pages 1 [49] S. Standring. Gray’s Anatomy. The Anatomical Basis of Clinical Practice. Elsevier, 40 edition, 2008. → pages 5, 6, 8, 9 [50] E. O. Stejskal and J. E. Tanner. Spin diffusion measurements: Spin echoes in the presence of a time-dependent field gradient. The Journal of Chemical Physics, 42(1):288, 1965. → pages 49 [51] C. Studholme, D. Hill, and D. Hawkes. An overlap invariant entropy measure of 3d medical image alignment. Pattern Recognition, 32:71–86, 1999. → pages 78 [52] A. Tanimoto, J. Nakashima, and S. Kuribayashi. Prostate cancer screening: the clinical value of diffusion-weighted imaging and dynamic mr imaging in combination with t2-weighted imaging. Magn Reson Imaging, 25: 146–4252, 2007. → pages 81 [53] P. Viola. Alignment by Maximization of Mutual Information. PhD thesis, Massachusetts Institute of Technology, 1995. → pages 72 [54] T. Winslow. Transrectal ultrasound, Apr 2008. URL http://visualsonline.cancer.gov/details.cfm?imageid=7224. → pages 15 [55] Q.-S. Xiang. Ubc physics 542: Magnetic resonance imaging, 2010. → pages 19, 20, 21, 23, 24, 25, 26, 27, 28, 29, 30, 31, 33, 34, 35, 36, 37, 38, 39, 40, 43, 44, 47 [56] J. Xu, P. A. Humphrey, A. S. Kibel, A. Z. snyder, V. R. Narra, J. J. Ackerman, and S.-K. Song. Magnetic resonance diffusion characteristics of histologically defined prostate cancer in humans. Magnetic Resonance in Medicine, 61:842–850, 2009. → pages iii, 81, 89, 90, 91, 94, 118, 119, 121, 124 [57] B. Zitova and J. Flussera. Image registration methods: a survey. Image and Vision Computing, 21(11):977–1000, 2003. → pages 56, 57 130 Part III Appendix 131 Appendix A Parameter Files for the Registration Algorithm A.1 Affine Transformation Parameter Files (Initial Parameter Files) A.1.1 Ex-vivo T 2 W Images and WMS Registration // Description: affine, MI, ASGD //ImageTypes (FixedInternalImagePixelType ”float”) (FixedImageDimension 2) (MovingInternalImagePixelType ”float”) (MovingImageDimension 2) //Components (Registration ”MultiResolutionRegistration”) (FixedImagePyramid ”FixedSmoothingImagePyramid”) (MovingImagePyramid ”MovingSmoothingImagePyramid”) (Interpolator ”BSplineInterpolator”) (Metric ”AdvancedMattesMutualInformation”) 132 (Optimizer ”AdaptiveStochasticGradientDescent”) //(Optimizer ”FiniteDifferenceGradientDescent”) (ResampleInterpolator ”FinalBSplineInterpolator”) (Resampler ”DefaultResampler”) (Transform ”AffineTransform”) (ErodeMask ”false” ) (NumberOfResolutions 4) (HowToCombineTransforms ”Compose”) (AutomaticTransformInitialization ”true”) (AutomaticScalesEstimation ”true”) (WriteTransformParametersEachIteration ”false”) (WriteResultImage ”true”) (ResultImageFormat ”png”) (ResultImagePixelType ”unsigned short”) (UseDirectionCosines ”true”) (CompressResultImage ”false”) (WriteResultImageAfterEachResolution ”true”) (ShowExactMetricValue ”false”) //Maximum number of iterations in each resolution level: //(MaximumNumberOfIterations 975) (MaximumNumberOfIterations 300) //Number of grey level bins in each resolution level: (NumberOfHistogramBins 32 ) (FixedLimitRangeRatio 0.0) (MovingLimitRangeRatio 0.0) (FixedKernelBSplineOrder 3) (MovingKernelBSplineOrder 3) //Number of spatial samples used to compute the mutual information in each resolution level: (ImageSampler ”RandomCoordinate”) 133 (FixedImageBSplineInterpolationOrder 1 ) (UseRandomSampleRegion ”false”) (NumberOfSpatialSamples 2000) (NewSamplesEveryIteration ”true”) (CheckNumberOfSamples ”true”) (MaximumNumberOfSamplingAttempts 10) //Order of B-Spline interpolation used in each resolution level: (BSplineInterpolationOrder 1) //Order of B-Spline interpolation used for applying the final deformation: (FinalBSplineInterpolationOrder 1) //Default pixel value for pixels that come from outside the picture: (DefaultPixelValue 63500) A.1.2 In-vivo T 2 W Images and WMS Registration // Description: affine, MI, ASGD //ImageTypes (FixedInternalImagePixelType ”float”) (FixedImageDimension 2) (MovingInternalImagePixelType ”float”) (MovingImageDimension 2) //Components (Registration ”MultiResolutionRegistration”) (FixedImagePyramid ”FixedSmoothingImagePyramid”) (MovingImagePyramid ”MovingSmoothingImagePyramid”) (Interpolator ”BSplineInterpolator”) (Metric ”AdvancedMattesMutualInformation”) (Optimizer ”AdaptiveStochasticGradientDescent”) //(Optimizer ”FiniteDifferenceGradientDescent”) (ResampleInterpolator ”FinalBSplineInterpolator”) (Resampler ”DefaultResampler”) (Transform ”AffineTransform”) 134 (ErodeMask ”false” ) (NumberOfResolutions 4) (HowToCombineTransforms ”Compose”) (AutomaticTransformInitialization ”true”) (AutomaticScalesEstimation ”true”) (WriteTransformParametersEachIteration ”false”) (WriteResultImage ”true”) (ResultImageFormat ”png”) (ResultImagePixelType ”unsigned short”) (UseDirectionCosines ”true”) (CompressResultImage ”false”) (WriteResultImageAfterEachResolution ”true”) (ShowExactMetricValue ”false”) //Maximum number of iterations in each resolution level: (MaximumNumberOfIterations 975) //Number of grey level bins in each resolution level: (NumberOfHistogramBins 32 ) (FixedLimitRangeRatio 0.0) (MovingLimitRangeRatio 0.0) (FixedKernelBSplineOrder 3) (MovingKernelBSplineOrder 3) //Number of spatial samples used to compute the mutual information in each resolution level: (ImageSampler ”RandomCoordinate”) (FixedImageBSplineInterpolationOrder 1 ) (UseRandomSampleRegion ”false”) (NumberOfSpatialSamples 2000) (NewSamplesEveryIteration ”true”) (CheckNumberOfSamples ”true”) (MaximumNumberOfSamplingAttempts 10) //Order of B-Spline interpolation used in each resolution level: (BSplineInterpolationOrder 1) //Order of B-Spline interpolation used for applying the final deformation: 135 (FinalBSplineInterpolationOrder 1) //Default pixel value for pixels that come from outside the picture: (DefaultPixelValue 63500) A.1.3 In-vivo DTI to in-vivo T 2 W Images Registration // Description: affine, MI, ASGD //ImageTypes (FixedInternalImagePixelType ”float”) (FixedImageDimension 2) (MovingInternalImagePixelType ”float”) (MovingImageDimension 2) //Components (Registration ”MultiResolutionRegistration”) (FixedImagePyramid ”FixedSmoothingImagePyramid”) (MovingImagePyramid ”MovingSmoothingImagePyramid”) (Interpolator ”BSplineInterpolator”) (Metric ”AdvancedMattesMutualInformation”) (Optimizer ”AdaptiveStochasticGradientDescent”) //(Optimizer ”FiniteDifferenceGradientDescent”) (ResampleInterpolator ”FinalBSplineInterpolator”) (Resampler ”DefaultResampler”) (Transform ”AffineTransform”) (ErodeMask ”false” ) (NumberOfResolutions 4) (HowToCombineTransforms ”Compose”) (AutomaticTransformInitialization ”true”) (AutomaticScalesEstimation ”true”) (WriteTransformParametersEachIteration ”false”) (WriteResultImage ”true”) (ResultImageFormat ”png”) (ResultImagePixelType ”unsigned short”) 136 (UseDirectionCosines ”true”) (CompressResultImage ”false”) (WriteResultImageAfterEachResolution ”true”) (ShowExactMetricValue ”false”) //Maximum number of iterations in each resolution level: (MaximumNumberOfIterations 400) //Number of grey level bins in each resolution level: (NumberOfHistogramBins 32 ) (FixedLimitRangeRatio 0.0) (MovingLimitRangeRatio 0.0) (FixedKernelBSplineOrder 3) (MovingKernelBSplineOrder 3) //Number of spatial samples used to compute the mutual information in each resolution level: (ImageSampler ”RandomCoordinate”) (FixedImageBSplineInterpolationOrder 1 ) (UseRandomSampleRegion ”false”) (NumberOfSpatialSamples 2000) (NewSamplesEveryIteration ”true”) (CheckNumberOfSamples ”true”) (MaximumNumberOfSamplingAttempts 10) //Order of B-Spline interpolation used in each resolution level: (BSplineInterpolationOrder 1) //Order of B-Spline interpolation used for applying the final deformation: (FinalBSplineInterpolationOrder 1) //Default pixel value for pixels that come from outside the picture: (DefaultPixelValue 63500) 137 A.2 Bsplines Transformation Parameter Files (Final Parameter Files) A.2.1 Ex-vivo T 2 W Images and WMS Registration // ImageTypes (FixedInternalImagePixelType ”float”) (FixedImageDimension 2) (MovingInternalImagePixelType ”float”) (MovingImageDimension 2) //Components (Registration ”MultiResolutionRegistration”) (FixedImagePyramid ”FixedSmoothingImagePyramid”) (MovingImagePyramid ”MovingSmoothingImagePyramid”) (Interpolator ”BSplineInterpolator”) (Metric ”AdvancedMattesMutualInformation”)// ”TransformBendingEnergy”) (Optimizer ”AdaptiveStochasticGradientDescent”) (ResampleInterpolator ”FinalBSplineInterpolator”) (Resampler ”DefaultResampler”) (Transform ”BSplineTransform”) // ********** Pyramid // Total number of resolutions (NumberOfResolutions 4) // ********** Transform (AutomaticTransformInitialization ”true”) (AutomaticScalesEstimation ”true”) (HowToCombineTransforms ”Compose”) // ********** Optimizer //(MaximumNumberOfIterations 200 50 30 10) (MaximumNumberOfIterations 200 50 40 25) (FinalGridSpacingInVoxels 32) 138 (GridSpacingSchedule 8 8 2 2 2 2 2 2) // ********** Metric //Number of grey level bins in each resolution level: (NumberOfHistogramBins 32) (FixedKernelBSplineOrder 3) (MovingKernelBSplineOrder 3) // ********** Several (WriteResultImageAfterEachResolution ”true”) (WriteTransformParametersEachIteration ”false”) (WriteTransformParametersEachResolution ”false”) (ShowExactMetricValue ”false”) // ********** ImageSampler // Number of spatial samples used to compute the // mutual information in each resolution level: (ImageSampler ”RandomCoordinate”) (NumberOfSpatialSamples 10000) (NewSamplesEveryIteration ”true”) // ********** Interpolator and Resampler //Order of B-Spline interpolation used in each resolution level: (BSplineInterpolationOrder 1) //Order of B-Spline interpolation used for applying the final deformation: (FinalBSplineInterpolationOrder 1) //Default pixel value for pixels (DefaultPixelValue 63500) (WriteResultImage ”true”) (ResultImageFormat ”png”) (ResultImagePixelType ”unsigned short”) (UseDirectionCosines ”true”) (CompressResultImage ”false”) (Metric0Weight 1.0) (Metric1Weight 0.2) 139 A.2.2 In-vivo T 2 W Images and WMS Registration // ImageTypes (FixedInternalImagePixelType ”float”) (FixedImageDimension 2) (MovingInternalImagePixelType ”float”) (MovingImageDimension 2) //Components (Registration ”MultiResolutionRegistration”) (FixedImagePyramid ”FixedSmoothingImagePyramid”) (MovingImagePyramid ”MovingSmoothingImagePyramid”) (Interpolator ”BSplineInterpolator”) (Metric ”AdvancedMattesMutualInformation”)// ”TransformBendingEnergy”) (Optimizer ”AdaptiveStochasticGradientDescent”) (ResampleInterpolator ”FinalBSplineInterpolator”) (Resampler ”DefaultResampler”) (Transform ”BSplineTransform”) // ********** Pyramid // Total number of resolutions (NumberOfResolutions 4) // ********** Transform (AutomaticTransformInitialization ”true”) (AutomaticScalesEstimation ”true”) (HowToCombineTransforms ”Compose”) // ********** Optimizer (MaximumNumberOfIterations 200 50 30 100) (FinalGridSpacingInVoxels 32) (GridSpacingSchedule 8 8 2 2 2 2 2 2) // ********** Metric 140 //Number of grey level bins in each resolution level: (NumberOfHistogramBins 32) (FixedKernelBSplineOrder 3) (MovingKernelBSplineOrder 3) // ********** Several (WriteResultImageAfterEachResolution ”true”) (WriteTransformParametersEachIteration ”false”) (WriteTransformParametersEachResolution ”false”) (ShowExactMetricValue ”false”) // ********** ImageSampler // Number of spatial samples used to compute the // mutual information in each resolution level: (ImageSampler ”RandomCoordinate”) (NumberOfSpatialSamples 10000) (NewSamplesEveryIteration ”true”) // ********** Interpolator and Resampler //Order of B-Spline interpolation used in each resolution level: (BSplineInterpolationOrder 1) //Order of B-Spline interpolation used for applying the final deformation: (FinalBSplineInterpolationOrder 1) //Default pixel value for pixels (DefaultPixelValue 63500) (WriteResultImage ”true”) (ResultImageFormat ”png”) (ResultImagePixelType ”unsigned short”) (UseDirectionCosines ”true”) (CompressResultImage ”false”) (Metric0Weight 1.0) (Metric1Weight 0.2) 141 A.2.3 In-vivo DTI to in-vivo T 2 W Images Registration // ImageTypes (FixedInternalImagePixelType ”float”) (FixedImageDimension 2) (MovingInternalImagePixelType ”float”) (MovingImageDimension 2) //Components (Registration ”MultiResolutionRegistration”) (FixedImagePyramid ”FixedSmoothingImagePyramid”) (MovingImagePyramid ”MovingSmoothingImagePyramid”) (Interpolator ”BSplineInterpolator”) (Metric ”AdvancedMattesMutualInformation”) (Optimizer ”AdaptiveStochasticGradientDescent”) (ResampleInterpolator ”FinalBSplineInterpolator”) (Resampler ”DefaultResampler”) (Transform ”BSplineTransform”) // ********** Pyramid // Total number of resolutions (NumberOfResolutions 4) // ********** Transform (AutomaticTransformInitialization ”true”) (AutomaticScalesEstimation ”true”) (HowToCombineTransforms ”Compose”) // ********** Optimizer (MaximumNumberOfIterations 200 50 30 100) (FinalGridSpacingInVoxels 32) (GridSpacingSchedule 8 8 2 2 2 2 2 2) // ********** Metric //Number of grey level bins in each resolution level: (NumberOfHistogramBins 32) 142 (FixedKernelBSplineOrder 3) (MovingKernelBSplineOrder 3) // ********** Several (WriteResultImageAfterEachResolution ”true”) (WriteTransformParametersEachIteration ”false”) (WriteTransformParametersEachResolution ”false”) (ShowExactMetricValue ”false”) // ********** ImageSampler // Number of spatial samples used to compute the // mutual information in each resolution level: (ImageSampler ”RandomCoordinate”) (NumberOfSpatialSamples 10000) (NewSamplesEveryIteration ”true”) // ********** Interpolator and Resampler //Order of B-Spline interpolation used in each resolution level: (BSplineInterpolationOrder 1) //Order of B-Spline interpolation used for applying the final deformation: (FinalBSplineInterpolationOrder 1) //Default pixel value for pixels (DefaultPixelValue 63500) (WriteResultImage ”true”) (ResultImageFormat ”png”) (ResultImagePixelType ”unsigned short”) (UseDirectionCosines ”true”) (CompressResultImage ”false”) (Metric0Weight 1.0) (Metric1Weight 0.2) 143 A.3 Example of an Inversion Parameter File // ImageTypes (FixedInternalImagePixelType ”float”) (FixedImageDimension 2) (MovingInternalImagePixelType ”float”) (MovingImageDimension 2) //Components (Registration ”MultiResolutionRegistration”) (FixedImagePyramid ”FixedSmoothingImagePyramid”) (MovingImagePyramid ”MovingSmoothingImagePyramid”) (Interpolator ”BSplineInterpolator”) (Metric ”DisplacementMagnitudePenalty”)// ”TransformBendingEnergy”) (Optimizer ”AdaptiveStochasticGradientDescent”) (ResampleInterpolator ”FinalBSplineInterpolator”) (Resampler ”DefaultResampler”) (Transform ”BSplineTransform”) // ********** Pyramid // Total number of resolutions (NumberOfResolutions 4) // ********** Transform (AutomaticTransformInitialization ”true”) (AutomaticScalesEstimation ”true”) (HowToCombineTransforms ”Compose”) // ********** Optimizer (MaximumNumberOfIterations 200 50 30 100) (FinalGridSpacingInVoxels 32) (GridSpacingSchedule 8 8 2 2 2 2 2 2) // ********** Metric //Number of grey level bins in each resolution level: (NumberOfHistogramBins 32) 144 (FixedKernelBSplineOrder 3) (MovingKernelBSplineOrder 3) // ********** Several (WriteResultImageAfterEachResolution ”true”) (WriteTransformParametersEachIteration ”false”) (WriteTransformParametersEachResolution ”false”) (ShowExactMetricValue ”false”) // ********** ImageSampler // Number of spatial samples used to compute the // mutual information in each resolution level: (ImageSampler ”RandomCoordinate”) (NumberOfSpatialSamples 10000) (NewSamplesEveryIteration ”true”) // ********** Interpolator and Resampler //Order of B-Spline interpolation used in each resolution level: (BSplineInterpolationOrder 1) //Order of B-Spline interpolation used for applying the final deformation: (FinalBSplineInterpolationOrder 1) //Default pixel value for pixels (DefaultPixelValue 63500) (WriteResultImage ”true”) (ResultImageFormat ”png”) (ResultImagePixelType ”unsigned short”) (UseDirectionCosines ”true”) (CompressResultImage ”false”) (Metric0Weight 1.0) (Metric1Weight 0.2) 145
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- In-vivo 3T and ex-vivo 7T diffusion tensor imaging...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
In-vivo 3T and ex-vivo 7T diffusion tensor imaging of prostate cancer : correlation with histology Uribe Muñoz, Carlos Felipe 2012
pdf
Page Metadata
Item Metadata
Title | In-vivo 3T and ex-vivo 7T diffusion tensor imaging of prostate cancer : correlation with histology |
Creator |
Uribe Muñoz, Carlos Felipe |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | Diffusion Tensor Imaging has been successfully applied in prostate cancer diagnosis (Kozlowski et al., 2010). It has been well established that the water Apparent Diffusion Coefficient has a lower value in the prostate carcinomas when compared to normal prostatic tissue (Bashar, 2002; Gürses et al., 2008; Kozlowski et al., 2010; Manenti et al., 2007; Pickles et al., 2005; Sato et al., 2005; Xu et al., 2009). However, fractional anisotropy values in prostatic carcinoma have been reported to be higher (Gürses et al., 2008), lower (Manenti et al., 2007), and unchanged(Xu et al., 2009) when compared to the prostate’s normal peripheral zone. Preliminary data from a study involving diffusion tensor imaging measurements in prostate glands, in-vivo and ex-vivo following radical prostatectomy, is presented. Histology whole mount slides were registered to T2 weighted images and diffusion parametric maps using a mutual information voxel intensity registration algorithm using software developed in-house. Regions of interest which included the normal peripheral zone, the normal peripheral zone with enlarged glands, and tumours were taken into account for this study. The tumours were highlighted and graded with the Gleason score grading system by a specialized pathologist. Values of the apparent diffusion coefficient and the fractional anisotropy parameters were calculated. Monte-Carlo simulations of the behaviour of the fractional anisotropy with respect to the value of the apparent diffusion coefficient, the signal to noise ratio, and the b-value of the Stejskal-Tanner equation were performed. The results show lower values of the apparent diffusion coefficient for regions of tumours for the ex-vivo and in-vivo cases. Values of the fractional anisotropy in the prostate carcinomas are slightly higher ex-vivo than in-vivo, which may be explained by the dependence of the fractional anisotropy on partial volume effects and noise. These preliminary results show that the fractional anisotropy does not show significant differences between normal and cancerous tissue, strongly suggesting that it is not likely to contribute significantly to the diagnostic capabilities of diffusion tensor imaging in prostate cancer. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-04-02 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0072658 |
URI | http://hdl.handle.net/2429/41915 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2012-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2012_spring_uribe_munoz_carlos.pdf [ 18.26MB ]
- Metadata
- JSON: 24-1.0072658.json
- JSON-LD: 24-1.0072658-ld.json
- RDF/XML (Pretty): 24-1.0072658-rdf.xml
- RDF/JSON: 24-1.0072658-rdf.json
- Turtle: 24-1.0072658-turtle.txt
- N-Triples: 24-1.0072658-rdf-ntriples.txt
- Original Record: 24-1.0072658-source.json
- Full Text
- 24-1.0072658-fulltext.txt
- Citation
- 24-1.0072658.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0072658/manifest