Prostate Registration UsingMagnetic Resonance Elastography forCancer LocalizationbyGuy NirB.Sc., Technion – Israel Institute of Technology, 2006M.Sc., Technion – Israel Institute of Technology, 2010A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2014c© Guy Nir 2014AbstractNoninvasive detection and localization of prostate cancer in medical imaging is an impor-tant, yet difficult task. Benefits range from diagnosis of cancer, to planning and guidance ofits treatment. In order to characterize cancer and evaluate its localization in volumetric im-ages, such as ultrasound or magnetic resonance imaging (MRI), their spatial correspondencewith the “gold standard” provided by histopathology must be established.In this thesis, we propose a general framework for a multi-slice to volume registrationthat is applied to register a stack of sparse, unaligned two-dimensional histological slices toa three-dimensional volumetric imaging of the prostate. The approach uses particle filteringthat allows deriving optimal pose parameters of the slices in a Bayesian approach.We then propose a novel registration method between in vivo and ex vivo MRI of theprostate to facilitate its registration to histopathology. The method incorporates elasticityinformation, acquired by magnetic resonance elastography (MRE), to generate a patient-specific biomechanical model of the prostate and periprostatic tissue.Next, we propose a registration method between preoperative MRI and intraoperativetransrectal-ultrasound. The method can be incorporated with a robotic surgical system toaugment the surgeon’s visualization during robot-assisted prostatectomy. We also studythe use of elasticity-based registration of ultrasound elastography and MRE.We then present an image processing approach for enhancing MRE data. The approachemploys registration to compensate for motion of patients during the scan to improve theaccuracy of the reconstructed elastogram. A super-resolution technique is employed toincrease the resolution of the acquired images by utilizing unique properties of MRE.Finally, we develop a theory for optimization-based design of motion encoding in MREthat allows reducing scanning time and increasing signal-to-noise ratio of elasticity recon-struction. We formulate the displacement estimation of the mechanical wave as an ex-perimental design problem, by which we quantify performance of sequences, and optimizemultidirectional designs.The proposed methods have been evaluated in simulations and on a diverse set of clinicaldata. Results may pave the way for a broader clinical deployment of elastography andelastography-based image processing.iiPrefaceThis thesis is based on several manuscripts, resulting from collaboration between multipleresearchers. All publications have been modified to make the thesis a coherent document.Ethical approval for conducting this research has been provided by the Clinical ResearchEthics Board, certificate numbers: H09-03163, H08-02696.• A version of Chapter 2 has been published as:G. Nir, R. S. Sahebjavaher, P. Kozlowski, S. D. Chang, E. C. Jones, S. L. Goldenberg,and S. E. Salcudean, “Registration of whole-mount histology and volumetric imagingof the prostate using particle filtering,” IEEE Transactions on Medical Imaging, 2014,vol. 33, no. 8, pp. 1601–1613.The author’s contribution in that article was proposing and developing the idea, for-mulating and implementing the algorithm, conducting numerical simulations, col-lecting the data, conducing the experiments, evaluating the results, and writing themanuscript. Dr. Sahebjavaher contributed to the MRI data collection, sectioning thespecimens, and scanning the histological slides. Dr. Kozlowski contributed to the7-Tesla MRI data collection and provided the histology sectioning device. Dr. Changassisted in segmentation of the data and identification of corresponding landmarks.Dr. Jones processed the whole-mount histological slices and provided the histopatho-logical analysis. Dr. Goldenberg contributed to the US data collection, performedthe radical prostatectomy procedures, and provided the postoperative specimens. Dr.Salcudean provided guidance, support, and assisted in editing the manuscript.• A version of Chapter 3 has been published as:G. Nir, R. S. Sahebjavaher, P. Kozlowski, S. D. Chang, R. Sinkus, S. L. Goldenberg,and S. E. Salcudean, “Model-based registration of ex vivo and in vivo MRI of theprostate using elastography,” IEEE Transactions on Medical Imaging, 2013, vol. 32,no. 7, pp. 1349–1361.The author’s contribution in that article was proposing and developing the idea, for-mulating and implementing the algorithm, conducting numerical simulations, col-iiilecting the data, conducing the experiments, evaluating the results, and writing themanuscript. Dr. Sahebjavaher contributed to the article through discussions and col-lected the MRI data. Dr. Kozlowski contributed to the article through discussionsand assisted with the 7-Tesla MRI data collection. Dr. Chang assisted in segmenta-tion of the data and identification of corresponding landmarks. Dr. Sinkus helped toset up the MRE system at UBC. Dr. Goldenberg performed the radical prostatec-tomy procedures and provided the postoperative specimens. Dr. Salcudean providedguidance, support, and assisted in editing the manuscript.• A version of Chapter 6 is in press for publication as:G. Nir, R. S. Sahebjavaher, R. Sinkus, and S. E. Salcudean, “A framework foroptimization-based design of motion encoding in magnetic resonance elastography,”Magnetic Resonance in Medicine, 2014, DOI: 10.1002/mrm.25280.The author’s contribution in that article was proposing and developing the idea, pro-viding mathematical formulation and analysis, conducting numerical simulations, col-lecting the phantom data, processing the data, evaluating the results, and writing themanuscript. Dr. Sahebjavaher contributed to the article through discussions, imple-mented the proposed MRI pulse sequences, collected the phantom data, provided thecode for phase unwrapping, and assisted in editing the manuscript. Dr. Sinkus helpedto set up the MRE system at UBC, contributed to the article through discussions,and assisted in editing the manuscript. Dr. Salcudean provided guidance, support,and assisted in editing the manuscript.• A version of Chapter 4 is to be submitted as:G. Nir, R. S. Sahebjavaher, J. Lobo, O. Mohareri, S. S. Mahdavi, and S. E. Salcud-ean, “Registration of Ultrasound and Magnetic Resonance Imaging for Image GuidedRobot-Assisted Prostatectomy,”The author’s contribution in this work was developing and implementing the regis-tration algorithms, conducting numerical simulations, collecting the phantom data,assisting in collection of patient data, processing the data, evaluating the results, andwriting the manuscript. Dr. Sahebjavaher collected the phantom and patient MRIand MRE data, contributed to the work through discussions, and assisted in editingthe manuscript. Mr. Lobo and Mr. Mohareri collected the phantom and patient USand USE data, assisted in integration of the algorithm with the da Vinci system, con-tributed to the work through discussions, and assisted in editing the manuscript. Dr.S. S. Mahdavi provided the code and support for the semi-automatic segmentationivalgorithm. Dr. Salcudean provided guidance, support, and assisted in editing themanuscript. In addition, Ms. S. Rashid and Ms. A. Ruszkowski provided evaluationof the surface-to-surface registration algorithm.• A version of Chapter 5 is to be submitted as:G. Nir, R. S. Sahebjavaher, M. Honarvar, and S. E. Salcudean, “Motion Compensa-tion and Super-Resolution in Magnetic Resonance Elastography,”The author’s contribution in this work was proposing, developing and implementingthe motion compensation and super-resolution idea, conducting numerical simula-tions, collecting the phantom data, assisting in collection of patient data, processingthe data, evaluating the results, and writing the manuscript. Dr. Sahebjavaher col-lected the phantom and patient data, contributed to the work through discussions,and assisted in editing the manuscript. Mr. Honarvar provided reconstruction of theelastograms for the motion compensated data, contributed to the work through dis-cussions, and assisted in editing the manuscript. Dr. Salcudean provided guidance,support, and assisted in editing the manuscript.Portions of the above papers, mainly introductory text, figures and literature review,also appear in Chapter 1.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 Prostate Anatomy and Prostate Cancer Interventions . . . . . . . . 11.1.2 Imaging of the Prostate . . . . . . . . . . . . . . . . . . . . . . . . . 31.1.3 Elastography of the Prostate . . . . . . . . . . . . . . . . . . . . . . 41.1.4 Whole-Mount Histopathology of the Prostate . . . . . . . . . . . . . 61.1.5 Prostate Registration . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Registration of Whole-Mount Histology and Volumetric Imaging of theProstate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.2 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2.1 Slice Transformation and Volume Re-Slicing . . . . . . . . . . . . . 15vi2.2.2 Multi-Slice to Volume Registration . . . . . . . . . . . . . . . . . . 162.3 Particle Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3.1 Prediction Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3.2 Observation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.3.3 Initial State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.3.4 Particle Filtering Algorithm . . . . . . . . . . . . . . . . . . . . . . 202.3.5 Similarity Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.4 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.4.1 Synthetic Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.4.2 Clinical Data Acquisition . . . . . . . . . . . . . . . . . . . . . . . . 272.4.3 Parameter Selection and Implementation Details . . . . . . . . . . . 302.4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 343 Model-Based Registration of the Prostate Using Elastography . . . . . 373.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.2.1 Elastic Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.2.2 Region-Based Active Contours . . . . . . . . . . . . . . . . . . . . . 403.2.3 Model-Based Alignment . . . . . . . . . . . . . . . . . . . . . . . . . 423.3 Model-Based Registration Using Inhomogeneous Elasticity . . . . . . . . . 433.3.1 Model-Based Force . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.3.2 Inhomogeneous Regularization . . . . . . . . . . . . . . . . . . . . . 443.3.3 Energy Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . 443.3.4 Synthetic Example . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.3.5 Clinical Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.4 Elastic Homogeneity Versus Inhomogeneity . . . . . . . . . . . . . . . . . . 483.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.5.1 Data Acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.5.2 Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . . 533.5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.6 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 574 Registration of Ultrasound and Magnetic Resonance Imaging for ImageGuided Interventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60vii4.2 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.3 Segmentation-Based Registration . . . . . . . . . . . . . . . . . . . . . . . 634.3.1 Semi-Automatic Prostate Segmentation in TRUS . . . . . . . . . . 634.3.2 Registration Process . . . . . . . . . . . . . . . . . . . . . . . . . . . 644.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.3.4 Application in MR-Guided RALP . . . . . . . . . . . . . . . . . . . 684.4 Elasticity-Based Registration . . . . . . . . . . . . . . . . . . . . . . . . . . 684.4.1 Registration Process . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 725 Motion Compensation and Super-Resolution in Magnetic Resonance Elas-tography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 755.2.1 Motion Compensation Using Registration . . . . . . . . . . . . . . . 755.2.2 Super-Resolution in MRE . . . . . . . . . . . . . . . . . . . . . . . . 765.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 785.3.1 Phantom Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . 785.3.2 Patients’ Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 805.4 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 816 Optimization-Based Design of Motion Encoding in Magnetic ResonanceElastography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 846.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 846.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 876.2.1 Parameterized Motion Encoding Gradients . . . . . . . . . . . . . . 876.2.2 Parameterized Phase Accumulation . . . . . . . . . . . . . . . . . . 876.2.3 Optimal Estimation of Wave Properties . . . . . . . . . . . . . . . . 896.2.4 Optimal Design of Motion Encoding . . . . . . . . . . . . . . . . . . 906.2.5 Unidirectional Designs . . . . . . . . . . . . . . . . . . . . . . . . . 916.2.6 Multidirectional RME Design . . . . . . . . . . . . . . . . . . . . . 926.2.7 Optimal Design with Varying Start Phases . . . . . . . . . . . . . . 946.2.8 Optimal Design with Varying Phase Shifts . . . . . . . . . . . . . . 966.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 976.3.1 Simulation of a Plane Wave . . . . . . . . . . . . . . . . . . . . . . 976.3.2 Phantom Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . 99viii6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1016.4.1 Simulation of a Plane Wave . . . . . . . . . . . . . . . . . . . . . . 1016.4.2 Phantom Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . 1016.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 1027 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1077.1 Contributions and Limitations . . . . . . . . . . . . . . . . . . . . . . . . . 1087.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114AppendixA Derivation of the Parameterized Phase Accumulation . . . . . . . . . . . 131ixList of Tables2.1 Summary of the parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.2 Registration results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.1 MRE results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.2 Registration results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.1 Segmentation-based registration results. . . . . . . . . . . . . . . . . . . . . 665.1 Super-resolution results on clinical data . . . . . . . . . . . . . . . . . . . . 826.1 Summary and description of the notations . . . . . . . . . . . . . . . . . . . 866.2 Phantom results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102xList of Figures1.1 Prostate anatomy with major anatomical landmarks . . . . . . . . . . . . . 21.2 Imaging of the prostate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 MRE of the prostate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.4 TRUS-VE of the prostate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.5 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Problem illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Random particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3 Similarity metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.4 Registration of synthetic data: illustration . . . . . . . . . . . . . . . . . . . 252.5 Registration of synthetic data: sensitivity analysis . . . . . . . . . . . . . . 272.6 Registration results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.1 Problem illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.2 Registration of synthetic images . . . . . . . . . . . . . . . . . . . . . . . . . 453.3 Registration of noisy synthetic images . . . . . . . . . . . . . . . . . . . . . 463.4 Mean displacement error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.5 Registration process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.6 Inhomogeneity parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.7 Difference maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.8 Inhomogeneity parameter performance . . . . . . . . . . . . . . . . . . . . . 513.9 Inhomogeneous registration results . . . . . . . . . . . . . . . . . . . . . . . 553.10 Inhomogeneous registration results . . . . . . . . . . . . . . . . . . . . . . . 574.1 Problem illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.2 Semi-automatic segmentation. . . . . . . . . . . . . . . . . . . . . . . . . . . 644.3 Segmentation-based registration . . . . . . . . . . . . . . . . . . . . . . . . . 654.4 Proposed augmented visualization . . . . . . . . . . . . . . . . . . . . . . . 674.5 Elasticity-based registration of a phantom . . . . . . . . . . . . . . . . . . . 70xi4.6 Elasticity-based registration of patient data . . . . . . . . . . . . . . . . . . 715.1 Motion compensation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 775.2 Super-resolution of a phantom MRE . . . . . . . . . . . . . . . . . . . . . . 795.3 Phantom errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 805.4 Super-resolution of a patient’s prostate MRE . . . . . . . . . . . . . . . . . 816.1 Example of MEGs sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . 886.2 RME design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 936.3 MD-RME design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 956.4 Optimal design with varying start phases . . . . . . . . . . . . . . . . . . . 966.5 Optimal design with varying phase shifts . . . . . . . . . . . . . . . . . . . . 986.6 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 996.7 Phantom data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1006.8 Simulation errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1016.9 Phantom results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1036.10 Phantom errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1047.1 Clinical flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108xiiList of Abbreviations2-D two-dimensional3-D three-dimensionalCaP carcinoma of the prostateDFT discrete Fourier transformDRE digital rectal examinationEPI echo-planar imagingFEM finite element methodFSE fast spin-echoFOV field of viewi.i.d. identically and independently distributedLFE local frequency estimationMAP maximum a posterioriMD-RME multidirectional reduced motion encodingMEG motion encoding gradientMI mutual informationMR magnetic resonanceMRE magnetic resonance elastographyMRI magnetic resonance imagingMSE mean squared errorxiiiNCC normalized cross-correlationPDE partial differential equationpdf probability density functionPF particle filterPSA prostate-specific antigenRALP robot-assisted laparoscopic radical prostatectomyRF radio-frequencyRME reduced motion encodingSE spin-echoSNR signal-to-noise ratioSTDV standard deviationSSD sum of squared differencesT2w T2-weightedTE echo timeTR pulse repetition timeTRE target registration errorTRUS transrectal ultrasoundUS ultrasoundUSE ultrasound elastographyVE vibro-elastographyxivAcknowledgementsI sincerely thank my supervisor, Professor Septimiu (Tim) E. Salcudean, for his supportthroughout the course of this work. His knowledge, experience, and constructive criticism,both on the professional and personal levels, are invaluable.I thank my fellow members of the Robotics and Control Laboratory for their collaborationand friendship during these years at UBC.I further thank my collaborators at Vancouver General Hospital, BC Cancer Agency, andthe UBC MRI Research Centre for their time, help and fruitful discussions.Last but not least, I thank my family for their love and support, with a special gratitudeto Moran.xvTo Carmel –the evergreen mountain in my lifexviChapter 1Introduction1.1 Background and Motivation1.1.1 Prostate Anatomy and Prostate Cancer InterventionsThe prostate is a fibromuscular gland about the size and shape of a walnut that lies at alow level in the lesser pelvis, and is traversed by the urethra and ejaculatory ducts [48].Being somewhat pyramidal, it presents a base superiorly, an apex inferiorly, and surfacesanteriorly, posteriorly and laterally. The prostate surrounds the prostatic urethra and is en-compassed superiorly by the bladder neck, inferiorly by the urogenital diaphragm, anteriorlyby the symphysis pubis and pubic arch, and posteriorly by the rectum wall (Figure 1.1). Theglandular tissue may be subdivided into three distinct zones: peripheral (70% by volume),central (25% by volume), and transition (5% by volume).With an estimated number of 233,000 new cases and 29,480 deaths in 2014, prostatecancer is the most commonly diagnosed cancer and second most common cause of cancerdeath among American men [6]. Prostate cancer, also known as carcinoma of the prostate(CaP), occurs mostly in the peripheral zone of the prostate, followed by the transition zone,then by the central zone. In its early stages, CaP has no symptoms and is usually suspectedwhen a prostate-specific antigen (PSA) test has shown a significant increase and/or whena digital rectal examination (DRE) shows stiff nodules or an asymmetric gland. Diagnosiscannot be confirmed until a biopsy test is returned positive.Clinical management of CaP is based on combination of the patient’s physical condition(e.g., age, weight, other health issues), tumor staging, tumor grading, and PSA [114]. In theGleason grading system, a grade is assigned by a pathologist to the most- and second-mostdominant microscopic patterns of a tumor based on histopathology analysis. Each graderanges from one to five, representing well- to poorly-differentiated carcinoma, respectively.The final Gleason score is the sum of the primary and secondary grades.There are various treatment methods available to CaP patients. These include watchfulwaiting, hormone therapy, chemotherapy, external beam radiation therapy, low/high doserate brachytherapy, and high-intensity focused ultrasound. Radical prostatectomy, i.e., the11.1. Background and MotivationVerumontanumBladder neckEjaculatory ductsSphincter musclesPeripheral zoneCentral zoneProstatic urethraTransition zoneSeminal vesicles Bladder(a) Coronal cross-sectionPeripheral zoneCentral zoneFibromuscular stromaTransition zoneSeminal vesiclesEjaculatory ductsRectal wallBladderUrethra(b) Sagittal cross-sectionFigure 1.1: Prostate anatomy with major anatomical landmarks.surgical extirpation of the prostate, is viewed by many as the “gold standard” treatment forclinically localized CaP. The long term prognosis of patients who undergo prostatectomyhas improved significantly over the past few decades. However, there remain significant ratesof disease recurrence and complications. In addition to regular risks of a surgical procedure,common complications include incontinence and impotence. Therefore, the primary goalsof prostatectomy are (i) cancer control, (ii) preservation of urinary continence, and (iii)preservation of erectile function.The state-of-the-art radical prostatectomy technique is the minimally invasive robot-assisted laparoscopic radical prostatectomy (RALP), which involves a robotic surgery sys-tem (da Vinci Surgical System, Intuitive Surgical, Sunnyvale, CA). Compared with theopen prostatectomy and standard laparoscopic prostatectomy techniques, RALP providesthe surgeon with excellent three-dimensional (3-D) visualization of the surgical site and im-proved dexterity. Studies show comparable [3, 89, 115], or superior [99] outcomes of RALPcompared with other prostatectomy techniques. The consensus [39] is that while RALPresults in lower blood loss, hospitalization duration, and complication rates, it suffers fromloss of haptic feedback, longer operating time, and cost [81].The usage of the da Vinci System is rapidly increasing worldwide, with RALP performedin as many as four out of five radical prostatectomies procedures in the United States [119].In addition to the advantages above, the RALP platform has been the subject of researchand development projects to further improve performance and surgery outcomes. Such toolsinclude real-time ultrasound (US) data [80, 94] that is displayed on the surgeon’s consolefor image guidance and improved visualization of anatomy and cancer (Figure 4.4(a)).21.1. Background and Motivation(a) In vivo T2w (b) In vivo B-mode (c) Ex vivo T2wFigure 1.2: Imaging of the Prostate. Examples of modalities used in our data collection.1.1.2 Imaging of the ProstateNoninvasive detection and localization of CaP in medical imaging is an important, yetdifficult task. Benefits range from diagnosis of cancer, to preoperative planning and in-traoperative carrying out of its treatment. An accurate localization of CaP may lead toless conservative treatments that reduce possible side effects such as impotence and urinaryincontinence.The most common modalities used for imaging of the prostate and CaP are transrectalultrasound (TRUS) B-mode images and magnetic resonance imaging (MRI) T2w images(Figure 1.2). The advantages of TRUS imaging are mainly that it is widely accessible,simple to operate, fast (real-time), non-ionizing, inexpensive and can be easily used intra-operatively. However, B-mode images present poor quality and are not accurate enough toreliably detect tumors and delineate the prostate [152].On the other hand, high resolution T2w images (typically acquired with an endorec-tal coil) provide a very detailed view of the prostate and periprostatic anatomy. Multi-parametric MRI presents the most accurate image based characterization of CaP (see [75]and references therein), and can be used for detection, localization and grading of CaP.Yet, magnetic resonance (MR) acquisition is time consuming, resource costly and, in manytimes, inaccessible, which makes it inapplicable for real-time guidance.1It is therefore a common practice for health care practitioners to employ MR preop-eratively for CaP detection and treatment planning, and to employ US intraoperativelyfor guidance. There are ongoing efforts to integrate preoperative MR data into US-guidedprocedures for better intraoperative localization of anatomy and cancer (see Chapter 4).Imaging of the postoperative specimens using MRI is also of interest, as these images canbe used to further improve CaP characterization by correlating data with confirmed tumors.1MR-guidance for prostate biopsy and treatment has been reported, e.g., in [153, 172], but is still cum-bersome and not widely used.31.1. Background and Motivation1.1.3 Elastography of the ProstateThe characterization of mechanical properties of biological tissues is important in the de-tection and diagnosis of diseases such as CaP, in which affected tissues typically altertheir viscoelastic properties. Manual palpation, e.g., DRE in case of CaP, provides anon-quantitative estimation of such properties that is limited to the surface of the body.Elastography, first introduced in [108], refers to noninvasive techniques for characterizingmechanical properties of tissues.An elastography system typically comprises an actuator for delivering mechanical ex-citations to the tissue, an imaging hardware for capturing and computing the resultingdisplacements, and an algorithm for solving the inverse problem of extracting an elasticityimage from the displacements. The elasticity image, also known as elastogram, usuallycorresponds to the elastic (Young’s) modulus of the imaged tissues. Elastography has beenshown to aid in CaP detection [118] and interventions [74].Elastography volumes are typically acquired at the same session along with “regular”imaging with respect to the modality, i.e., US B-mode and MR T2w images. These imagescan be aligned with their corresponding elastograms using simple or no processing. In vibro-elastography (VE), the radio-frequency (RF) data collection to compute the elastogram, isusually followed by a B-mode scan of the same volume under the same conditions. In mag-netic resonance elastography (MRE), the elastogram is computed from the phase componentof the acquired signal, while its magnitude is available as low resolution T2w images.Magnetic Resonance ElastographyIntroduced in [100], MRE uses a sinusoidal mechanical excitation and a synchronized MRpulse sequence that encodes the tissue motion in the phase of T2w (akin to) images. Thetechnique have been applied to study tissue changes associated with diseases of, e.g., theliver, breast, prostate, brain, and heart [46]. A detailed description of the acquisition processcan be found, e.g., in [150].In [135], an MRE system with an actively shielded electromagnetic transducer is pro-posed for in vivo prostate imaging (Figure 1.3(a)). Excitations are applied to the perinealregion of the subject for efficient transmission of mechanical waves to the prostate. MREelastograms have been shown to be promising in improving the visibility of the prostategland. However, to speed up the acquisition process, so as to minimize inconvenience andmovements of the patient, the images typically have low resolution. Also, the boundary ofthe prostate may be blurred and assimilated with the background.41.1. Background and Motivation(a) Transperineal MRE setup([135])(b) Sagittal T2w (c) MRE magnitude(d) MRE phase 3/8 (e) MRE phase 5/8 (f) MRE phase 7/8 (g) MRE elastogramFigure 1.3: MRE of the prostate. (a) Illustration of the transperineal driver (green) posi-tioning. (b) The prostate (cyan) overlaid on a sagittal scout image of the patient. (c) Aslice of the magnitude image. (d)-(f) Three phases of displacements along the Y-axis, pi/2apart, depict wave propagation in the tissue. (g) The corresponding slice of the elastogram.Ultrasound ElastographyIn ultrasound elastography (USE), the resulting displacements of the vibrated tissue arerecorded using time domain normalized cross-correlation (NCC) between sequences of highframe rate RF data obtained from US imaging. These techniques can be classified intostatic and dynamic elastography methods. In static elastography, an axial force is appliedto the tissue, and the pre- and post-compressed sets of images are used to determine themotion along the direction of the applied force. Dynamic elastography consists of applyingvibration to the tissue, during which RF frames are acquired continuously and elastogramsare generated in real-time by estimating strain between sequential frames.The authors in [139] present a TRUS-VE system, in which transfer function analysisof tissue motion over a frequency range of simultaneous vibrations is used to acquire 3-Delastograms of the prostate and periprostatic tissue (Figure 1.4(a)). Compared to B-modeimages, with MR T2w as “gold standard,” the produced elastograms have shown, bothqualitatively (visually) and quantitatively, to have superior object-background contrast,especially around the base and apex [143]. The authors concluded that the VE imaging51.1. Background and MotivationTRUS cradleRotation motorVibration motorTRUS probeSagittal plane Transverse plane(a) TRUS-VE system ([139]) (b) TRUS B-mode (c) TRUS-VEFigure 1.4: TRUS-VE of the prostate.is a promising modality for prostate interventions. However, edge continuity of prostateboundaries in the elastogram was generally inferior to B-mode.In [9], a USE method using free-hand US is proposed. The technique, referred to asabsolute elastography, involves a specific excitation frequency of the tissue, and trackingits motion using the algorithms proposed in [9, 37]. Assuming linearity, the motion hasthe same temporal frequency, and the phase and amplitude of the displacement at eachpixel can be estimated similar to MRE. The inverse problem is then solved, e.g., by localfrequency estimation (LFE), to reconstruct the elastogram.1.1.4 Whole-Mount Histopathology of the ProstateWhole-mount histopathology refers to placing an entire organism or structure directly ontoa microscope slide in order to study the manifestations of disease by a pathologist. Whenspecimen size allows and longer processing time is available, whole-mount histopathologypreserves the structural relationships between cells, and facilitates the comparison of thespecimen to imaging data. In whole-mount histopathology of the prostate, analysis ofcomplete transversal cross-sections, placed on histological slices, help in predicting the longterm disease outcome of the CaP patients, and provide a ground truth to determine theability of imaging to detect cancer.The process required for pathologists to prepare scans of whole-mount histological slicesin our study is summarized as follows (see also Figure 2.1). The specimen should firstbe fixed in 10% formalin solution for about 72 hours, followed by its transverse sectioninginto 4 mm thick parallel sections using a single blade or a multi-bladed cutting device [32].Each section is embedded in a paraffin block, followed by slicing a thin, 4 µm, layer fromeach block using a microtome. Each slice is then glided onto a glass slide, followed by itsstaining in hematoxylin and eosin (H&E). Finally, the stained slides are digitally scannedto generate a series of digital images that represent the different levels of the specimen.61.1. Background and Motivation1.1.5 Prostate RegistrationRegistration is the process of establishing a spatial correspondence between two or moreimage data sets. The images may be acquired at different times, from different perspectivesor by different imaging techniques. In medical image processing, registration is a criticaltask for integration of information acquired from pre-, intra-, and postoperative imaging. InCaP diagnosis and treatment, registration can, for example, enable US-guided proceduresto benefit from the accuracy of cancer detection and anatomical detail in MR.In general, given a reference image and a (set of) template image, the registration prob-lem is to find a geometric transformation or mapping, such that the transformed (deformed)template is close to the reference, in some sense. A detailed classification of the numerousregistration methods in medical imaging is provided in [84]. The main classification criteriarelevant to the algorithms discussed in this thesis are the dimensionality of the registra-tion, the nature of transformation the registration is confined to, and the nature of theregistration basis.The spatial dimensionality we consider is registration of a set of two-dimensional (2-D)images to a 3-D volume (Chapter 2), and registration of a 3-D volume to a 3-D volume(Chapters 3-5). The nature of the spatial transformations we propose is rigid, i.e., trans-lation and rotation, with an additional isotropic or anisotropic scaling (Chapters 2-4), anddeformable (also known as curved or nonrigid) transformations (Chapters 3-5).In terms of the registration basis, the focus is on intrinsic registration methods ofintensity-based (Chapters 2, 4-5), model-based (Chapters 2-3), and segmentation-based(Chapters 2, 4) nature. Intensity-based methods utilize solely image intensities to computethe mapping. In model-based methods, a template model of the target structure is ex-tracted (preoperatively, usually) and registered to the reference image. Segmentation-basedmethods register a segmented template into a segmented reference.The registration problem is often formulated as an optimization of a similarity metricbetween the deformed template and the reference images, with respect to a set of param-eters (parametric techniques), or to the mapping displacements directly (nonparametrictechniques) [154]. In the former case, in order to correct for the problem’s ill-posednessand produce more probable maps, the transformation is confined to a function space andthe parameters to be optimized are the weights of the basis functions. In the latter case, aregularization metric is defined and combined with the similarity metric in order to smooththe displacements. Registration algorithms differ in the choice of the similarity metric, pa-rameters or regularization metric, and in the numerical method for solving the optimizationproblem.71.1. Background and MotivationCommon choices for the similarity measure are the sum of squared differences (SSD),and NCC. While simple and, in many cases, provide plausible results, these methods aresensitive to noise and require the reference and template to be of the same intensity scale.Therefore, these measures are not suitable for multi-modality registration. An extensivelyused similarity measure, that is applicable to multi-modality registration, is the (normalized)mutual information (MI) [26, 165].Many parametric registration techniques employ geometric transformations derived frominterpolation theory, such as polynomial or B-spline function spaces. In nonparametrictechniques, a variety of metrics have been proposed to regularize the deformations of thetemplate according to geometrical transformations derived from physical models. A pop-ular physical regularizer, which we also employ in this thesis (Chapters 3-5), is the elasticpotential of the displacements [11, 20]. Other physical metrics are based on diffusion [58],fluid [25] and curvature [41] models. Another metric, the earth mover’s distance, preservesthe “mass” associated with the images, and yields the optimal mass transport [52].An extensive survey of deformable medical image registration methods is presentedin [154]. A general framework for many of these techniques is presented in [92]. Theframework is based on a variational formulation and the corresponding Euler-Lagrangeequations that characterize a minimizer. Numerical schemes to solve the resulting partialdifferential equations (PDEs) can be formulated and discretized using the finite differencemethod. Other schemes employ the finite element method (FEM) [176].Biomechanical models for registration are segmentation-based methods, in which thereference and/or template segmented models are associated with mechanical properties,e.g., elasticity, viscosity. The template model can then be deformed in a physical manner,in order for it to match the reference model, by applying forces. Typically, FEM is employedto numerically model the deformation.In [38], a linear elastic biomechanical model is proposed for registration of MR brainimages. The method was adapted in [16] for registration of the prostate. In [157], theauthors propose a global surface difference to register linear elastic models extracted fromprostate TRUS and MR images. A nonlinear elastic model that is driven by a distance fieldrepresentation of the models, was used for fast registration of prostate data in [106].Biomechanical models produce plausible and physical deformations. However, accuracyof the registration highly depends on accuracy of the segmentation to extract the models.An intensity-based registration of the prostate in TRUS, that does not require segmentation,was presented in [13] for biopsy tracking. The technique involves a generalized correlationratio as a similarity metric, with elastic regularization and an inverse consistency constraint.81.2. Thesis Objectives1.2 Thesis ObjectivesThe objectives of this thesis are as follows:• Developing and validating methods for registration between histopathology slices andvolumetric imaging of the prostate that allow both training and evaluation of CaPclassifiers.• Developing and validating methods for registration between preoperative MRI andintraoperative TRUS during RALP, in order to facilitate surgery and minimize sideeffects while avoiding positive margins.• Incorporating elastography information into the registration framework, and investi-gating its effect on registration performance.• Improving the quality and/or shortening acquisition time of MRE.1.3 Thesis OutlineThe overall goal of the research undertaken in this thesis is to develop algorithms for reg-istration of the prostate among different modalities. The motivation for this research is toprovide tools that allow for an automatic preoperative and intraoperative localization ofCaP. In addition, the thesis studies the reciprocity between registration and elastography.The diagram in Figure 1.5 depicts the general structure and flow of the thesis. We beginwith a “pure” registration algorithm, proceed to incorporation of elastography in registra-tion that can be applied to detection of CaP and its treatment, then we incorporate imageregistration to enhance MRE reconstruction, and end with a general scheme to improveMRE acquisition.This chapter introduced the research topic, motivation for performing this research andthe thesis objectives. The remainder of this thesis is organized as follows. In Chapter 2, weconsider a multi-slice to volume registration method in which a stack of sparse, unaligned2-D whole-mount histological slices is registered to a 3-D volumetric imaging of the prostate.A particle filtering framework is proposed in order to derive optimal registration parametersin a Bayesian approach. We demonstrate and evaluate our method on a diverse set of datathat includes a synthetic volume, ex vivo and in vivo MRI, and in vivo US.In Chapter 3, we propose a novel registration method that uses a patient-specific biome-chanical model acquired using MRE to deform the in vivo volume and match it to the surface91.3. Thesis OutlineChap. 5Chap. 6Chap. 2Chap. 3Chap. 4IntroductionConclusionMRERegistrationElastography in registrationRegistration in MREFigure 1.5: Thesis outline.of the ex vivo specimen. The incorporation of elastography data into the registration frame-work allows the measured inhomogeneous elasticity to be assigned to the in vivo volume.The method is demonstrated and evaluated on simulation and clinical cases.In Chapter 4, we propose a segmentation-based registration method between preopera-tive MRI and intraoperative TRUS. The resulting corresponding images may be displayedon the surgeon’s console during RALP. Incorporation of elastography is proposed as aninter-modality to facilitate MRI-TRUS alignment, by elasticity-based registration of thecorresponding MRE and USE.In Chapter 5, we employ registration and image processing techniques to improve MREreconstruction. Registration is used to compensate for patient’s motion during acquisition,while super-resolution technique is used to enhance the resolution of the acquired images.In Chapter 6, we present an optimization-based approach for encoding motion in MRE.Multidirectional sequences are derived by setting the problem in an experimental designframework. Such designs are implemented and evaluated on simulation and phantom data.Finally, in Chapter 7, we summarize the major contribution of this thesis and describepotential research directions in the future.10Chapter 2Registration of Whole-MountHistology and Volumetric Imagingof the Prostate2.1 IntroductionIn order to characterize cancer and to evaluate its localization in volumetric images, suchas those obtained by an US or MR imaging, the images must be compared to the “goldstandard” provided by a histopathological analysis of the excised tissue following biopsy orprostatectomy. Thus, registration of histology to the imaging is required so as to find theirspatial correspondence.The overall objective of the work described in this chapter is to develop and validatea method for automatic registration between 2-D digital scans of whole-mount histologicalslices and 3-D volumes of the prostate. Such multi-slice to volume registration is challengingdue to the dissimilarity between the modalities’ intensity values, the similarity in geome-try of different cross-sections of the volume, the sparsity of the histological slices and theinaccuracies that occur during their preparation.Each step in the preparation process of whole-mount histopathology, as described inSection 1.1.4 above, introduces uncertainty to the geometric origin of the histological slices inthe volumetric images. We distinguish between the different mechanisms for misalignments:• In-plane misalignment – gliding the slice and scanning the slide are prone to 2-Drotation and translation of the histology image. These misalignments are statisticallyindependent among slices, but distributed around the correct pose, since efforts aremade to align the slices visually during preparation.• Out-of-plane misalignment – sectioning and slicing are susceptible to 3-D rotationand translation of the specimen with respect to the imaging volume. These misalign-ments may be caused directly by errors in position and angulation of the blade duringsectioning, and/or indirectly by elastic deformation of the specimen such as bending.112.1. Introduction(a) MR volume (b) Prostate specimen (c) Prostate sections (d) Histological slicesFigure 2.1: Problem illustration. (a) Volumetric imaging of the prostate, such as an in vivoMR, is being acquired. (b) The excised fixed prostate is sectioned, during which out-of-plane misalignments occur. (c) Each 4 mm section is embedded in a paraffin block, fromwhich a 4 µm thin layer is sliced. (d) Each slice is glided onto a slide, stained and scanned,during which in-plane misalignments occur. These scans should be then registered to thevolumetric imaging.While each slice may have a different pose, they are typically correlated, especiallywhen using a multi-bladed device, as described later, that maintains orientation andspacing between the sections.• Scaling – fixation of the specimen causes change in its volume. The scaling is anisotropicand typically manifested as shrinkage due to tissue dehydration [145]. In addition,implicit scaling can be caused by errors in measuring histology resolution, which isoften done manually. Scaling is highly correlated among the slices.• Residual deformation – the shape of the prostate may be deformed after imaging.An algorithm for registration of histological slices to a volumetric image should takethese mechanisms into account. Overall, excluding deformation, the number of degrees offreedom in registration of N histological slices to a volume is 9N , where typically N∼10.In addition, the transformation space is constrained, e.g., such that slices do not intersect.Figure 2.1 illustrates the registration problem.Previously, different apparatuses were proposed to minimize misalignment between his-tological slices. In [12], the authors have inserted needles into the specimen before sectioningto provide landmarks for correcting in-plane misalignment. In [32], a multi-bladed cuttingdevice that constrains the specimen mechanically during sectioning allows its slicing intoparallel cuts of fixed intervals, thus minimizing out-of-plane misalignment between slices.However, in both approaches, the sectioning axis might not coincide with the scanning axisof the imaging and further registration is needed.122.1. IntroductionThe authors in [129] proposed a rotating platform for imaging a desired plane in ex vivoMR, and then cut the specimen in parallel to that plane. The works [167, 168] involveplacement of fiducial markers into the ex vivo specimen, in order to guide its sectioningalong the imaging plane with a magnetically tracked probe, while [45] employs such markersto fit affine transformations between histological slices and the specimen. These methodsachieve good registration results, but require a special setup, and additional clinical timeand approvals that may not be available in practice.For registration of in vivo imaging, the authors in [148, 162] employ a 3-D printer togenerate a patient-specific prostate mold, in which the fixed specimen is placed and cut inorder to preserve its orientation. However, the postoperative specimen after fixation mightnot fit into the mold that is based on preoperative in vivo imaging. In [60], the authorspropose an apparatus for insertion of three needles from apex to base at distinct anglesinto fresh prostate specimens. The resulting punctures serve as fiducial markers, based onwhich the slices are aligned to each other. This alignment method is intended to allow a3-D registration to in vivo imaging. Although requiring a special setup, as well as invasivemanipulation of the specimen, the authors report that their protocol is employed in a rapidand standardized manner, which does not compromise the histopathology information.A 3-D deformable registration between histology and ex vivo MR of the prostate wasproposed in [175]. The method involves extraction of landmarks on, and inside the sur-face of the gland. However, in order to acquire landmarks on the histological surface, thesectioning should be dense, which may not be feasible in practice. The authors in [24]propose to register segmented histology with segmented multi-parametric in vivo MR. Acorresponding MR slice is chosen for each histological slice by an expert, followed by a 2-Dto 2-D registration. Such manual search for corresponding slices may be a difficult and/ortime consuming task for an expert.In [111], a 3-D histological volume is reconstructed from aligned slices, and a 3-D to 3-Dregistration method is employed to match it to the imaging volume. However, as discussedin [24], a volume reconstruction from the sparsely sectioned histology is inaccurate and maylead to registration errors. Thus, the authors in [171] consider a more practical registrationbetween multiple histology slices and an in vivo MR volume. They propose an iterativealgorithm that consists of three modules: finding a subset of MR slices that correspondsto the histological slices, 2-D registration of each histological slice to its corresponding MRslice, and 3-D registration of the MR volume to the registered histological slices.In summary, the approaches proposed in the literature typically require special hardware[32, 129], involve insertion of markers to the specimen [12, 45, 60, 167, 168], are 2-D in nature[24] or assume a dense sectioning for a 3-D registration [111], do not take into account all132.1. Introductiondegrees of freedom and/or do not utilize their interdependence, are sensitive to segmentationerrors, and/or are prone to local extrema [171]. While fusing several methods and runningthem multiple times with different initializations may allow for an accurate registration, itwould result in a cumbersome and long process that might not be suitable in some scenarios.Some works in the literature, e.g., [7, 43, 73, 82, 88, 93, 105, 117, 141], have appliedfiltering techniques to image registration. This approach allows derivation of an optimalestimation of the registration parameters in a Bayesian fashion to overcome issues suchas susceptibility to initialization and local extrema. Specifically, particle filtering was em-ployed for rigid/affine registration of points-to-surface [82] and two point clouds [43, 141]using a point-based metric, a model-to-slice registration using a region-based metric [105],an elastic/affine registration of two images using an intensity-based metric [7, 88], andoptimization of a parameterized deformable registration field [73].We propose a unified approach for registration of multiple slices to a volume in 3-Dwithout intervening with the standard histopathological processing. The method employsparticle filtering that allows derivation of an optimal estimation of the registration param-eters in a Bayesian fashion to overcome issues such as susceptibility to initialization andlocal extrema. The framework allows incorporation of prior knowledge on admissible mis-alignments according to the misalignment mechanisms discussed above, and knowledge ofthe imaging noise or segmentation error.A particle in our method represents a combination of the slices in various 3-D poses.Each such particle is assigned a likelihood based on the similarity of its slices to the volume,and a prior probability based on the admissibility of their pose misalignment. This allowsto implicitly neglect non-admissible pose of slices, associated with lower probabilities, anddiffuse the rest towards an optimal estimation of the true pose. This approach contends wellwith the multimodal nature of the optimization, and, as the number of particles increases,converges to the global optimum.Our approach generalizes the method in [171], in the sense that the histological slicesare not assumed to be parallel to each other, and the volumetric image does not need tobe segmented if it contains sufficient intensity information. As in [171], we do not handledeformations of the histological slices to avoid over fitting the sparse data. However, theaffine model we consider provides a first order approximation of the deformations, and adeformable registration algorithm may be employed as postprocessing to compensate forthe residual deformation. We evaluate the method on a diverse set of data that includes asynthetic volume, ex vivo MR, in vivo MR, and in vivo US.The remainder of this chapter is organized as follows. In Section 2.2, we formulatethe multi-slice to volume registration problem. In Section 2.3, we model the problem in142.2. Problem Formulationstate-space, and present the particle filtering algorithm for registration. In Section 2.4, weprovide details on implementation and the clinical experiments conducted to evaluate theproposed method, and summarize the results. Finally, we conclude in Section 2.5 with asummary and future research directions.2.2 Problem Formulation2.2.1 Slice Transformation and Volume Re-SlicingLet r¯ = (x, y, z)T ∈ R3 be the 3-D spatial coordinates vector. A rigid transformation ofr¯ can be parameterized by six pose parameters, namely, three translation values ∆i alongeach axis i = x, y, z, and three rotation angles θi about each axis. In addition, anisotropicscaling is parameterized by three scaling factors ξi of each axis. We denote the pose vectoras s¯ = ({∆i}, {θi}, {ξi})T , i = x, y, z.We represent the spatial transformation by a concatenation of matrices using homo-geneous coordinates, such that the transformed coordinates vector ˜¯r ∈ R3 is written as(˜¯r1)= T (s¯)(r¯1),whereT (s¯) = Tout-of-plane(∆z, ξz, θx, θy)Tin-plane(∆x,∆y, ξx, ξy, θz)= TθyTθxTξzT∆zTθzTξyTξxT∆yT∆x . (2.1)The matrix T ∈ R4×4 is the affine transformation matrix with respect to the posevector. Explicit expressions of the matrices can be found, e.g., in [53]. The specific orderof the matrix multiplication in (2.1) was set to distinguish between “in-plane” and “out-of-plane” transformations, such that the former will act prior to the latter, so as to model themisalignment mechanisms described in Section 2.1. This order also guarantees that slicessharing a same out-of-plane pose are parallel, even though their in-plane pose may vary.For notation purposes we will omit the “1” and simply write ˜¯r = T (s¯)r¯.Let {r¯0 = (x, y, z0)T : (x, y) ∈ R2} define a planar region or slice on z = z0. We mayapply an affine transformation with respect to a pose s¯0 to each point on the slice, suchthat {˜¯r0 = T (s¯0)r¯0} is the shifted, rotated and scaled slice.Suppose we are given a volumetric image, J : R3 → R. In order to find the values of theimage on the transformed slice, we can use Euler coordinates as the interpolation reference152.2. Problem Formulationframe and writeJ˜0(x, y; s¯0) = J(T (s¯0)r¯0), ∀(x, y) ∈ R2. (2.2)We refer to this operation as re-slicing of the volume J , and to J˜0 as the re-sliced image.2.2.2 Multi-Slice to Volume RegistrationAssume we are given N scans of 2-D histological slices with the approximated spacingbetween consecutive pairs. In addition, assume that the region inside the specimen in eachslice is segmented. We can represent this information as the set of triplets {In, Hn, zn}N−1n=0 ,where In : R2 → R is a grayscale image of the slice, Hn : R2 → {0, 1} is a binary imagewith values of “1” on the region inside the specimen, and zn ∈ R is the position or “height”of the slice in 3-D space. Without loss of generality, we set z0 = 0 and assign values for{zn}N−1n=1 with respect to the given spacing.Let J be a volumetric image of the same specimen the histology was processed from. Wedenote the 3-D pose vector of the nth slice as s¯n. The registration of histology to volumetricimaging problem is finding N poses s¯0:N−1 = {s¯n}N−1n=0 , such that the transformed segmentedhistological slices optimally match, in some sense, the corresponding re-sliced images of thevolume. The admissible space of 3-D poses is constrained such that the transformed slicesdo not intersect.In general, the multi-slice to volume registration can be formulated as the optimizationproblemmins¯0:N−1(admissible)N−1∑n=0αnD[In(·), J(·);Hn(·), s¯n], (2.3)where D[In(·), J(·);Hn(·), s¯n] is a distance functional that measures (dis-)similarity betweenthe slice In and the volumetric image J , with respect to the slice segmentation Hn and itspose s¯n. The weights {αn}N−1n=0 can be used to prioritize or neglect certain slices, as discussedin Section 2.4.3 below.The formulation (2.3) supports any choice of metric for measuring similarity. We con-sider three types of metrics, as discussed in Section 2.3.5. Specifically, we employ anintensity-based metric for the ex vivo MR data sets, a region-based metric for the in vivoMR, and a point-based metric for US data sets (in which we assume that the volumetricimage J is segmented as well).162.3. Particle Filtering(a) (b) (c)(d) (e)Figure 2.2: Random particles. Each particle represents a combination of 3-D poses assignedto the histological slices (in blue). The surface of the prostate in the volumetric image isshown (in red) as reference. Some particles are more likely to be correct than others.Particle (e) corresponds to the MAP estimation.2.3 Particle FilteringIntroduced in [47], particle filters (PFs) are sequential Monte Carlo methods based on pointmass (or “particle”) representations of probability density functions (pdfs). As the num-ber of particles becomes larger, this representation becomes closer to the usual functionaldescription of pdf, and the PF approaches the optimal Bayesian estimate. Unlike the tradi-tional Kalman filter and its extensions, PFs are capable in handling nonlinear state-spacemodels with non-Gaussian noise. The tradeoff is that PFs are computationally expensive,although, with the increase in computational power in recent years, PFs were successfullyemployed even in real-time applications.In order to set the registration problem in a filtering framework, we represent it as astate-space model. We define the state vector as S¯ = s¯0:N−1, i.e., a column vector comprisingthe 9N pose and scaling parameters of all slices. In this formulation, a point in state-space,or a particle, represents a combination of slices in various 3-D poses. Figure 2.2 showsexamples of different particles that were generated for one case.172.3. Particle Filtering2.3.1 Prediction ModelSince our problem is static, we assume a process in which the state S¯(k) at a time step (afiltering iteration) k ≥ 1, can be predicted from the previous state S¯(k−1) simply byS¯(k) = S¯(k−1) +√γk−1U (k), (2.4)where U is an identically and independently distributed (i.i.d.) system noise that reflectsthe uncertainty in the pose parameters, and 0 < γ ≤ 1 is a discount factor that reflects thedecrease of this uncertainty with every time step.In order to model the statistical dependency of the pose parameters among the slices,we define U to have a multivariate normal distribution with a specific form of covariance.We write s¯i,0:N−1 to denote a vector that comprises N values of the pose parameter ·i ofall slices, e.g., (θx,0, . . . , θx,N−1)T . We assume that s¯i,0:N−1 is normally distributed arounda zero mean with the N ×N covariance matrixΣs¯i(σs¯i , ηs¯i) =σ2s¯i σ2s¯i −η2s¯i2 · · · σ2s¯i −η2s¯i2σ2s¯i −η2s¯i2 σ2s¯i · · ·...... ... . . . ...σ2s¯i −η2s¯i2 · · · · · · σ2s¯i,where σ2s¯i is the variance of the pose parameter, and 0 < η2s¯i ≤ 3σ2s¯i is the variance of thedifference of this parameter between each pair of slices, i.e., η2s¯i = Var(s¯i,n1−s¯i,n2), ∀n1 6= n2.The smaller the value of η2s¯i is, the more (positively) correlated the slices become, withη2s¯i → 0 yielding the same parameter value for all the slices (as if the value was drawnfrom a univariate normal distribution with a variance σ2s¯i). Specifically, setting η2θx , η2θy → 0yields transformed slices that are parallel, and η2∆z → 0 yields equally spaced slices. Settingη2s¯i = 2σ2s¯i yields a diagonal covariance matrix, thus an independent parameter for each slice(with a variance σ2s¯i).According to the discussion on the misalignment mechanisms for each type of poseparameter in Section 2.1, low values of η2s¯i (high correlation among slices) are typicallyassigned to the out-of-plane pose and scaling parameters, and a value of η2s¯i = 2σ2s¯i (nocorrelation among slices) is assigned to the in-plane pose parameters.The system noise is therefore defined as U∼N (0,ΣU ), where ΣU is a 9N×9N covariancematrix that comprises the nine covariance matrices Σs¯i , one for each pose parameter, in itsdiagonal. The system (prediction) model (2.4), with the modeled noise U , can be employed182.3. Particle Filteringto formulate the prior pdf, i.e., the probability of the current state given the previous statevector, asp(S¯(k)|S¯(k−1)) ∝ exp{−12(S¯(k) − S¯(k−1))T1γk−1 Σ−1U(S¯(k) − S¯(k−1))}. (2.5)2.3.2 Observation ModelIn a state-space model, at each time k ≥ 1, a measurement (observation) becomes availableand is related to the state vector via the observation model. In case the intensity- or region-based metrics are being used, we consider that observation to be the static volumetric imageJ with an additive i.i.d. measurement noise V that represents the intensity noise of theimage. In the point-based metric, the points of the segmented surface are the observation,and the measurement noise represents the errors in segmentation and reconstruction of thesurface. In both cases, the noise is modeled as V∼N(0, σ2v), i.e., a zero mean Gaussiannoise with a variance σ2v .The relation between the observation and the state vector is a nonlinear, complex andunknown function. We can define the observation model implicitly by defining the likelihoodpdf, i.e., the probability of the observation given the state vector, asp(J |S¯(k)) ∝ exp{−N−1∑n=0αnD[In(·), J(·);Hn(·), s¯(k)n ]/σ2v}. (2.6)Note that if the state S¯(k) provides a good alignment between the histological slices Inand the volume J , the functional D takes small values which results in high likelihood. Inturn, high likelihood states minimize the objective functional (2.3). Such interpretation of ametric or an energy functional as log-likelihood has been used, e.g., as in [21], for switchingbetween variational and probabilistic formulation of optimization problems.2.3.3 Initial StateIn order to get an initial estimate of the state vector, we try to explicitly minimize theobjective functional (2.3). Of course, due to the multimodal nature of the problem, theprocess is most likely to converge to a local minimum. Nevertheless, this local minimumtypically provides an initialization to the particle filtering algorithm that allows a fasterconvergence to the global minimum.192.3. Particle FilteringMoreover, in order to simplify the minimization, we reduce the search space in this stepby “locking” the pose parameters of all the slices, i.e. s¯n1 = s¯n2 , ∀n1, n2. Thus, the slices inthis step move together in 3-D as a “sparse” volume. This minimization process brings theslices into the coordinate system of the volume in order to correct for global misalignment.It may also be seen as analogous to the volume alignment step in [171], where the entirevolume is transformed to maximize similarity with the histological slices on their pre-definedplanes.We minimize (2.3) with respect to the “locked” pose vector s¯ until it converges, suchthat the difference between the pose parameters of two consecutive iterations is less thansome threshold. Section 2.4.3 below provides additional details on the implementation andoptimization methods. We then apply the pose vector s¯ to all slices, and denote the resultingstate vector as S¯(0). With no better estimate about the pose, we define the initial state pdfas a deterministic probability distribution about the resulting state. We denote this pdf asp(S¯(0)) = δ(S¯(0) − S¯(0)), where δ is the Dirac delta function.2.3.4 Particle Filtering AlgorithmThe filtering problem is to recursively construct the posterior pdf of the state given theobservation, i.e., p(S¯(k)|J). A PF algorithm produces at each time step k, a cloud of Mstates {S¯m,(k)}Mm=1, known as particles, that provides a discrete weighted approximation ofthe posterior distribution.Several particle filtering schemes have been proposed in the past, but can be shown tobe derived from a generic algorithm [8]. For a comprehensive discussion on PFs, we refer thereader to [30, 31] and references therein. In this chapter we employ a sampling-importance-resampling filter. The algorithm, described below, is summarized as a pseudocode in Algo-rithm 2.1.The algorithm starts with sampling M times from the initial state distribution p(S¯(0)),and then employing the Bayes’ recursion at each time step. Assuming that one can samplefrom the posterior, an empirical estimate of this distribution is given byp(S¯(k)|J) ≈M∑m=1ωm,(k)δ(S¯(k) − S¯m,(k)), (2.7)where ωm,(k) is a weight associated with the mth particle, such that∑Mm=1 ωm,(k) = 1. Ingeneral, it is difficult to sample from the posterior. One solution is to employ importancesampling, which consists of using a distribution known as the proposal that is related to theposterior but easier to generate samples from, and then normalizing the weights.202.3. Particle FilteringWe follow an approach proposed in [125] for tracking. In this approach, predictedparticles { ˆ¯Sm,(k)}Mm=1 are first sampled from the prior distribution, and then propagatedtowards higher likelihood states based on observations via an update stage, which canbe interpreted as an importance sampling of the particles {S¯m,(k)}Mm=1. The weights arerecursively updated byωm,(k) ∝ ωm,(k−1) p(J |S¯m,(k)) p(S¯m,(k)|S¯m,(k−1)). (2.8)After executing the alignment process in Section 2.3.3 above, we generate M particlesaccording to p(S¯(0)), i.e., M state vectors that equal to the initial state. We set the particles’weights {ωm,(0)}Mm=1 uniformly to 1M and proceed to time step k = 1. Then, and at everynew time step, we predict a state { ˆ¯Sm,(k)}Mm=1 for each particle m, according to the systemmodel (2.4).The predicted state of each particle is propagated towards higher likelihood states, asdefined by (2.6), based on the “observed” volume J . We do this by minimizing the objectivefunctional for each particle, with its predicted state ˆ¯Sm,(k) taken as the initial state. Thisupdate stage is denoted asS¯m,(k) = FLmin( ˆ¯Sm,(k)), (2.9)where FLmin stands for performing L iterations of the optimization process in order to min-imize (2.3) with respect to the entire (“unlocked”) state of the particle s¯m,(k)0:N−1.Only a fixed number of L iterations is performed in order to avoid overfitting the stateof each particle to the observation. Large values of L lead to convergence of the slices tothe volume, and thus express trust in the observations, while small values of L keep theslices close to their predicted state, and thus express high trust in the initial alignment andsystem model. Too large values of L would result in slices reaching unfeasible states thatmaximize their similarity to the volume, perhaps even violating the non-intersecting slicesconstraint, and causing degeneration of particle diversity. On the other hand, too smallvalues would result in the particles staggering around the initial alignment state, keepingthem in low likelihood states.After all particles are updated, we assign an importance weight to each one, based on(2.8). To this end, we calculate the prior and likelihood for each particle S¯m,(k) using (2.5)and (2.6), respectively. The posterior pdf of the state given past and current observationsis then approximated by p(S¯(k)|J) ≈ ∑Mm=1 ωm,(k)δ(S¯(k) − S¯m,(k)). The posterior allowsan optimal estimate of the state to be obtained, with respect to any criterion. In ourexperiments, we used the maximum a posteriori (MAP) state, which corresponds to the212.3. Particle Filteringxy(a) intensity-, region-based: Initialxy(b) intensity-, region-based: Final(c) point-based: Initial (d) point-based: FinalFigure 2.3: Minimization of the similarity metrics. (a)-(b) The segmented region of ahistological slice is overlaid on the cross-section of the volumetric MR image on the re-slicing plane. In the intensity-based case, the re-slicing planes move in order to maximizethe MI metric between the cross-sections of the image and the corresponding histologicalslices. In the region-based case, the re-slicing planes move in order to minimize the intensityvariances of the image’s cross-sections inside and outside the segmentation contour of thecorresponding histological slices. (c)-(d) The segmented histological contours (filled in blue)are interwoven with the surface (in red) that was segmented from an US volume. The pointsrepresent the closest counterparts of the contour points on the surface. The minimizationprocess moves each slice to minimize the distance between the pairs of closest points.state of the highest weighted particle, i.e.,S¯MAP,(k) = S¯m?,(k), with m? = arg max1≤m≤Mωm,(k). (2.10)Finally, resampling is performed. We resample M times according to the approximatedposterior distribution above, to have new particles {S¯m,(k)}Mm=1 that replace the currentones. This step eliminates particles that have very low weights and concentrates on higherweighted particles. The importance weights for the resampled particles are reset to ωm,(k) =1M (thus, ωm,(k−1) can be omitted from (2.8)). The particle filtering algorithm iterates untilconvergence of the MAP state.222.3. Particle FilteringAlgorithm 2.1: Histological slices to volumetric image registration using particlefiltering (see Table 2.1 for a list of parameters).Input: {In}N−1n=0 , {Hn}N−1n=0 , {zn}N−1n=0 , J ,{σ2i }, {η2i }, σ2vM , L, γ, α, p(S¯(0)).Output: p(S¯(k)|J), S¯MAP,(k) (at each time step k).Initialization: generate {S¯m,(0)}Mm=1 according to p(S¯(0)),set weights {ωm,(0)}Mm=1 to 1M ,k = 1.repeatfor m = 1 to M doobtain prediction:ˆ¯Sm,(k) = S¯m,(k−1) +√γk−1Um,(k).perform updates: S¯m,(k) = FLmin( ˆ¯Sm,(k), J).calculate prior: p(S¯m,(k)|S¯m,(k−1)) using (2.5).calculate likelihood: p(J |S¯m,(k)) using (2.6).set weight:ωm,(k) ∝ p(J |S¯m,(k))p(S¯m,(k)|S¯m,(k−1)).end fornormalize weights such that:∑Mm=1 ωm,(k) = 1.approximate posterior:p(S¯(k)|J) = ∑Mm=1 ωm,(k)δ(S¯(k) − S¯m,(k)).resample {S¯m,(k)}Mm=1 according to the posterior.reset weights {ωm,(k)}Mm=1 to 1M .k := k + 1.until convergence.232.3. Particle Filtering2.3.5 Similarity MetricsWe consider three types of metrics according to the imaging modality of the volume. Onthe ex vivo MR data sets, for which intensity information has been shown to be useful inregistration to histology [24, 171], we employ the intensity-based MI metric.For the in vivo MR data sets, we found that the MI metric fails to match the histologicalslices to the volume on some cases, and we apply a region-based metric instead. Such metricsemploy intensity statistics to describe the prostate region, and have been shown to be usefulfor automatic segmentation of the prostate in T2w MR [163].For the US data sets, we found that the metrics above fail to measure similarity tohistology. Therefore, we assume that the volumetric image J is segmented as well and usea point-based metric, namely, the closest point metric. An illustration of the minimizationprocess of the two similarity metrics is shown in Figure 2.3.Intensity-BasedWe employ the inverse of the (normalized) MI [26, 165] between the histological slice Inand the re-sliced image J˜n with respect to the pose s¯n. In this case, the distance functionalin (2.3) isD[In(·), J(·);Hn(·), s¯n] =E[In(·), J˜n(·; s¯n)]E[In(·)] + E[J˜n(·; s¯n)], (2.11)where E[In] and E[J˜n] are the intensity entropies of the 2-D slice In and re-sliced image J˜n,respectively, while E[In, J˜n] is their joint entropy [122]. The entropies are calculated usingthe marginal and joint histograms of the region inside the segmented histological slice, andits corresponding region in the re-sliced volumetric image.Region-BasedWe employ the Chan-Vese metric [23] that is typically used as an active contours model forsegmentation. The metric minimizes the intensity variances of the re-sliced image J˜n onthe regions inside and outside the segmentation contour of the corresponding histologicalslice In. The metric is minimized with respect to the pose s¯n, and thus may be seen asfitting the re-sliced image to the contour of the histology, rather than evolving the contourto segment the image as was done in [163].242.4. Experiments(a) Initial state (b) After initial alignment (c) Final MAP state (d) True stateFigure 2.4: Registration of synthetic data: illustration of the algorithm stages. The ground-truth, (d), is shown in comparison with the MAP state resulting from the proposed algo-rithm, (c).In this case, the distance functional in (2.3) isD[In(·),J(·);Hn(·), s¯n] =∫R2(J˜n(x, y; s¯n)− cin)2Hn(x, y) dxdy+∫R2(J˜n(x, y; s¯n)− cout)2(1−Hn(x, y)) dxdy, (2.12)where cin and cout are the mean intensities of the re-sliced image J˜n inside and outside thehistological contour, respectively.Point-BasedHere, we assume that the contour of each histological slice In is represented by Qn points,{r¯Inq }Qnq=1, and the surface in the segmented volumetric image J is represented by a cloud ofP points, {r¯Jp }Pp=1. We employ the SSD between pairs of closest points in the point cloudsof the volume and in transformed histology slice ˜¯rInq = T (s¯n)r¯Inq with respect to the poses¯n. Thus, the distance functional in (2.3) isD[In(·), J(·);Hn(·), s¯n] =1QQ∑q=1‖˜¯rInq − r¯Jpq‖2, (2.13)where we employ the `2 norm ‖ · ‖ (Euclidean distance), and r¯Jpq is the point in J that isclosest to ˜¯rInq in In with respect to the norm.252.4. Experiments2.4 Experiments2.4.1 Synthetic CaseWe have tested and evaluated the proposed multi-slice to surface registration algorithm ona synthetic data set for which we know the ground truth. Specifically, we used the 3-D pointdata set of the Stanford Bunny [155] as the segmented volumetric image, and generated nineslices of 256 × 256 pixels from cross-sections of the volume. Rather than transforming thesurface, we applied random transformations to each of the slices by sampling an initial statefrom the prior distribution (2.5) around some mean state (different from the true state).We allowed both different in-plane poses among slices, in order to model the misalignmentduring slide scanning, and different out-of-plane poses that yield non-parallel slices (withoutviolating the non-intersecting slices constraint) that implicitly model bending of the surface.Figures 2.4(a)-2.4(c) show the registration process between the contours of the slicesand the surface after employing the proposed algorithm for 10 time steps with 100 par-ticles. Since the data is given as point cloud, we used the point-based metric describedin Section 2.3.5 below. Figure 2.4(d) shows the true state of the slices for a qualitativeevaluation of the results.A quantitative evaluation of the results is provided in Figure 2.5. In Figure 2.5(a),the mean squared distance between contour points of the MAP state’s slices and theircounterparts on the surface (the objective metric), as well as their mean squared distanceto the corresponding contour points of the ground truth state’s slices, obtained by using 100particles are plotted versus the number of filtering iterations. Indeed, the distances fromthe MAP to surface points, and from the MAP to ground truth points are minimized andconverge to 1.67± 0.65 pixels (0.7± 0.3%) and 2.98± 1.16 (1.1± 0.4%) pixels, respectivelyIn order to evaluate the sensitivity of the results to the number of particles, we reranthe experiment with five filtering iterations for different number of particles (the rest of theparameters were fixed). Figure 2.5(b) shows that indeed the trend is that the objectivemetric decreases as the number of particle increases. We note that comparison of theresults achieved by using fewer particles (a single particle in particular), to those achievedby using many particles, can be regarded as a comparison of our stochastic approach toa deterministic algorithm such as [171]. Also, we note that with 150 particles and fivefiltering iterations, the obtained distances from MAP to surface points, and from MAP toground truth points are 1.38± 0.86 pixels (0.5± 0.3%) and 2.59± 0.98 pixels (1.0± 0.4%),respectively, which are lower than those obtained for 100 particles and 10 filtering iterations.This is due to the quick convergence of the algorithm with respect to the filtering iterations,and the fact that more particles allow further “exploration” of the search space.262.4. Experiments0 1 2 3 4 5 6 7 8 9 100102030405060Filtering iterationsDistance(pixels) MAP to surfaceMAP to ground truth(a)0 50 100 1500246810121416ParticlesDistance(pixels) MAP to surfaceMAP to ground truth(b)Figure 2.5: Registration of synthetic data: sensitivity analysis. (a) sensitivity of the al-gorithm to the number of filtering iterations (with 100 particles). (b) sensitivity of thealgorithm to the number of particles (with 5 filtering iterations). The mean squared dis-tance between contour points of the MAP state and the surface/ground truth are used asperformance indicators (see text).2.4.2 Clinical Data AcquisitionAfter approval of the institutional ethics board and signed consents, ten patients (meanage 65, range 57–76 years old) who were scheduled for radical prostatectomy, underwentpreoperative MR and intraoperative US. A postoperative MR was acquired for the prostatespecimen as well. We have tested the proposed method in order to register the histologyscans to the different imaging modalities of each case.In Vivo MRA T2w scan in a 3-Tesla system (Achieva 3.0T, Philips, The Netherlands) was performed. Astandard 6-channel cardiac coil with acceleration factor (SENSE) 2 was used. A standardaxial T2w fast spin-echo (FSE) sequence was acquired (echo time (TE)/pulse repetitiontime (TR) = 86/2500 ms, field of view (FOV) = 320×320×70 mm3, with 0.5 mm in-planeresolution and 2 mm slice thickness).272.4.ExperimentsNotation Description Typical values RemarksK Number of time steps (filtering iterations) 10 More steps allow convergence of the particles, but increaseruntime; Increased with the degree of uncertainty.M Number of particles 250 More particles estimate the posterior better, but are compu-tationally expensive; Increased with the number of slices N .L Number of update iterations 5 (intensity-based) High number of iterations lead to (noisy) data overfitting;10 (region-based) Low number of iterations lead to slow convergence.25 (point-based)αn Slice weighting 0.5 (n = 0, N − 1), High values prioritize slices;1 (otherwise) Is used to neglect the base and apex.γ Pose discount factor 0.85 Low values reflect decrease of uncertainty (pose STDV)with every filtering iteration.σ∆x,y In-plane translation STDV 3 mm In-plane misalignment occurs during sliding the slicesσθz In-plane rotation STDV 3◦ and scanning the slides;η∆x,y In-plane translation slice difference STDV√2 · σ∆x,y Statistically independent among slices.ηθz In-plane rotation slice difference STDV√2 · σθzσ∆z Out-of-plane translation STDV 1 mm (ex vivo MR) Out-of-plane misalignment occurs during sectioning3 mm (in vivo MR, US) of the specimen;σθx,y Out-of-plane rotation STDV 1◦ (ex vivo MR)3◦ (in vivo MR, US)η∆z Out-of-plane translation slice difference STDV 0.1 · σ∆z Highly dependent among slices when using a multi-ηθx,y Out-of-plane rotation slice difference STDV 0.1 · σθx,y bladed cutting device (equally spaced and parallel slices).σξx,y,z Scaling STDV 0.1 (ex vivo MR) Occurs during excision and fixation of the specimen;0.3 (in vivo MR, US)ηξx,y,z Scaling slice difference STDV 0.01 · σξx,y,z Highly dependent among slices. Values are logarithmic.σv imaging noise STDV (intensity-, region-based) 3% dynamic range Affects the likelihood of states; High values favour the/ surface error STDV (point-based metric) 1 mm updated state, lower values favour the predicted state.Table 2.1: Summary of the parameters.282.4. ExperimentsIn Vivo USIntraoperatively, prior to the procedure, B-mode images were collected using an US machine(Sonix RP, Ultrasonix Medical, Richmond, BC, Canada) with a biplane transrectal probe(BPL9-5/55, 5–9 MHz) that was mounted on a motorized system [143], to collect a contin-uous sagittal sweep from −45◦ to 45◦ that spans a 3-D volume (FOV = 84× 50× 56 mm3).The volume was then interpolated in order to produce conventional transverse images ofthe prostate (with 0.3 and 0.42 mm in-plane and axial resolution, respectively).The prostate in the transverse volume was segmented by a radiologist on key slices (every∼5 − 7 slices in a volume of 128 slices), in which the prostate was visualized better. Theradiologist also used the corresponding sagittal view for assistance, and relied on symmetryof the gland to guide its contouring. We have used Stradwin (Cambridge University, UK),[160], which employs [161] in order to interpolate the sparse contours on intermediate slicesand generate a triangulated surface of the prostate. We set the surface resolution andsmoothing strength of the algorithm to medium.Ex Vivo MRAfter fixation and right before sectioning, the postoperative prostate specimen was wrappedin a plastic bag, taped to a dedicated holder and scanned in a 7-Tesla system (Biospec,Bruker, Germany). A rapid acquisition with relaxation enhancement (RARE) pulse se-quence was used to acquire axial T2w images (TE/TR = 70/5000 ms, FOV = 70× 70× 42mm3, with 0.2734 mm in-plane resolution and 2 mm slice thickness).Whole-Mount HistologyThe multi-bladed cutting device in [32] was used for whole-mount sectioning of the specimenin order to minimize orientation and spacing variance among the sections. Following thehistology preparation process described in Section 2.1, cancerous regions were marked bya pathologist on the histological slides. A typical number of 8-10 slides were then digitallyscanned using a commercial flatbed Scanner with 600 DPI that yields an in-plane resolu-tion of about 0.17 mm. The prostate boundary was segmented manually on each of thehistological slices. Although the segmentation was relatively fast (< 1 min per slice) andstraightforward due to the small number of slices and homogeneous background, the clearboundaries may also allow employment of a (semi-) automatic segmentation algorithm.292.4. Experiments2.4.3 Parameter Selection and Implementation DetailsWe have tuned the registration parameters by first running the algorithm multiple times ona study data set from each modality, and then applying these values to the rest of the datasets of the same modality. A summary of the different parameters and their assigned valuesis presented in Table 2.1. Below we explain the reasoning behind some of the choices.As a rule of thumb, the number of particles should be increased exponentially with thedimensionality of the state-space [30] (9N in our case). However, as in [125] and [141], due tothe importance sampling we employ, the actual growth in particles is less than exponential,and for a typical number of 10 slices, we found that 100− 250 particles are required. Withno time limitation, we have chosen the upper value of that range to increase the chances ofconverging to the global minimum.From our observations, we have also estimated typical values for the uncertainty inthe pose parameters. We set the variance of each pose parameter using the fact thatabout 99.7% of the samples drawn from a normal distribution are within three standarddeviations (STDVs) from the mean, which is initially taken as zero (before initial alignment).Since we have used a multi-bladed cutting device, we set low values to the variances of theout-of-plane pose difference among slices to reflect the high probability of equally spacedand parallel slices. Without prior knowledge about the in-plane pose, we set their differencevariances to allow a statistical independence. We have also set the imaging noise to berelative to the dynamic range of the images, when the intensity- or region-based metricswere employed, and estimated the surface segmentation error for the point-based metric.We note that, while keeping the same number of particles and assuming equally spacedand parallel slices, the algorithm is robust with respect to the choice of the pose variances,with similar results achieved for a range of STDVs in the order of a few mm or degrees.This is due to a fast convergence (2− 3 filtering iterations) of the particles to a state-spaceregion around the true state.To minimize the metrics with respect to the state vector we employed the Nelder-Meadsimplex method for the intensity-based metric, a gradient descent for the region-basedmetric, and the iterative closest point (ICP) algorithm [15] for the point-based metric. Asdescribed in Section 2.3.4, the number of iterations in the update stage (2.9) was chosento be relatively low to avoid overfitting of the histological slices to the volumetric image.We note that, since the Nelder-Mead optimization process is not smooth, the algorithm issensitive to the number of update iterations when employing the intensity-based metric, andthis number should be kept low to avoid “drifting” of the state. Consequently, more particlesand filtering iterations are usually required in comparison to the point-based metric.302.4. ExperimentsWe used a discount factor of γ = 0.85 to reflect the decrease in uncertainty with everytime step. We also set the weights α0 = αN−1 = 0.5 and αn = 1, n = 1, . . . , N − 2, in orderto prioritize the mid-gland slices, for which, typically, prostate imaging is more reliable thanfor the top and bottom slices (base and apex). This is due to ambiguity in boundaries ofthe gland on extreme slices, where its apex and base blend into the pelvic floor and thebladder neck tissue, respectively, and thus we would like the algorithm to be less affectedby them. When using a sufficient number of particles, the algorithm is robust with respectto the choice of values for both the discount factor and weights. However, employing theseparameters may speed up the convergence of the algorithm, and in our experiments theywere set by a trial-and-error process on example cases.We have implemented our registration method and tested it under MATLAB on a 2.93GHz Quad Core CPU machine with 8 GB of RAM. The average runtime per registrationwas 25 minutes for the point-based metric cases, 8 hours for the region-based metric cases,and 3 hours for the intensity-based metric cases.2.4.4 ResultsTo evaluate the registration performance on the clinical data, we compare the registeredhistological slices with the corresponding re-sliced images of the volume. For each slice, wecompare its histological segmentation with the corresponding cross-section of the prostatesurface, as obtained from a manual segmentation of the volumetric image by a radiologist.We measure the area overlap of the regions inside the two segmentation contours, in thesense of Dice’s coefficient, i.e., twice the size of the intersected area divided by the sum ofthe sizes of the areas.In addition, for the MR cases, the radiologist has identified matching landmarks in-side the prostate on both the histology and volumetric images after manually selectingcorresponding slices. These landmarks include the urethra, nodules, scars (from previousbiopsies), calcifications and other general distinguished anatomical features.We measure the Euclidean distance between the registered landmarks on histology andtheir counterparts in the volume. These distances provide us a good approximation ofthe target registration error. Since, typically, the 3-D coordinates of the landmarks in thevolume do not lie exactly on the registered slice, we show the projections of the landmarkson that plane (however, the error is calculated in 3-D).Selected histological slices and their corresponding re-sliced volume images according tothe registration results are shown in Figure 2.6 for each type of imaging. The quantitativeresults for each modality are also summarized in Table 2.2.312.4.Experimentsxy10 20 30 40 50 600510152025303540 xy0 5 10 15 20 25 30 35 40 45 5005101520253035404550xy10 20 30 40 50 600510152025303540xy20 30 40 50 60 70 8005101520253035404550 xy0 20 40 60 80 100 120020406080100120xy20 30 40 50 60 70 8005101520253035404550xy0 10 20 30 40 50 60051015202530354045(a) Histological slicexy0 10 20 30 40 50 60 70 80051015202530354045(b) Imaged slice (selected manually) (c) Registration re-sultsxy0 10 20 30 40 50 60051015202530354045(d) Registered slice (auto-matic)Figure 2.6: Registration results. Example slices from ex vivo MR (top row), in vivo MR (middle row) and in vivo US (bottomrow) imaging. A mid-gland histological scan (left column), a manually selected corresponding volumetric slice (second left column),visualization of the registration results (second right column), and the corresponding registered volumetric slice. The segmentation(contour) and landmarks (markers) of the histology (blue) and volumetric imaging (red) are overlaid on the images.322.4.ExperimentsHistology to ex vivo MR Histology to in vivo MR Histology to in vivo USArea overlap Mean error Min-max error Area overlap Mean error Min-max error Area overlap(%) (mm) (mm) [# of landmarks] (%) (mm) (mm) [# of landmarks] (%)Case 1 93.0± 6.5 1.8± 1.0 0.8− 5.0 [11] 89.2± 5.7 3.4± 1.6 0.9− 5.6 [11] 91.1± 4.1Case 2 94.7± 1.4 1.5± 1.0 0.5− 3.5 [9] 91.0± 4.1 2.9± 1.6 0.3− 5.3 [9] 93.0± 2.7Case 3 93.1± 2.2 1.4± 0.5 0.4− 2.5 [14] 89.2± 6.4 3.2± 1.3 0.7− 5.2 [14] 92.7± 1.2Case 4 94.4± 3.2 1.6± 0.5 0.2− 2.6 [13] 86.1± 7.1 4.3± 1.6 1.0− 7.1 [13] 94.2± 2.1Case 5 90.2± 4.1 1.4± 0.5 0.6− 2.2 [19] 81.5± 6.3 4.1± 1.5 1.5− 6.7 [19] 89.1± 4.9Case 6 91.4± 5.0 1.6± 0.7 0.5− 3.0 [14] 87.1± 8.0 3.3± 1.1 1.1− 4.9 [14] 92.4± 1.8Case 7 94.1± 1.3 1.4± 0.3 0.8− 1.9 [10] 92.2± 5.3 2.8± 1.3 0.8− 4.6 [10] 90.5± 2.9Case 8 91.9± 5.5 1.6± 0.6 0.6− 2.6 [12] 84.4± 10.3 3.9± 1.2 1.8− 5.4 [12] 87.1± 8.4Case 9 89.9± 23.8 1.7± 0.6 1.0− 2.7 [12] 79.1± 8.9 6.1± 2.5 2.0− 9.4 [12] 85.1± 14.1Case 10 94.5± 1.8 4.4± 1.2 2.5− 6.5 [20] 84.7± 12.6 3.9± 1.0 1.8− 5.3 [20] 93.8± 0.9Mean 92.0± 5.6 1.7± 1.2 0.8− 3.3 [14] 86.9± 8.2 3.8± 1.5 1.3− 5.9 [14] 90.1± 5.8Table 2.2: Registration results. Quantitative summary.332.5. Discussion and Conclusion2.5 Discussion and ConclusionWe presented a method for a multi-slice to volume registration in order to solve the chal-lenging clinical problem of alignment between histological slices and volumetric imaging ofthe prostate. The proposed algorithm does not alter the clinical flow of histopathologicalprocessing, and can be employed to characterize CaP and evaluate its localization in vol-umetric images by allowing the mapping of lesions, which are marked by a pathologist onthe histological slices, onto the volume.The proposed algorithm has been evaluated on both synthetic and clinical data. Theexperiments on synthetic data, for which the ground truth of slices’ pose with respect tothe surface is known, showed a fast convergence of the algorithm and its ability to recoverthe pose while avoiding local optima. The registration of histological slices to ex vivoMR, in vivo MR, and in vivo US achieved relatively high area overlaps and registrationerrors below the diameter of typical lesions that are marked during histopathology analysis(∼10± 5 mm [34]).Lower error values were obtained in the ex vivo MR experiments. This is due to thefact that the ex vivo MR scan is performed after excision and fixation, and just beforesectioning, and therefore captures the shape of the deformed specimen, between which thehistological slices are more likely to have a rigid transformation. Surprisingly, the resultingarea overlaps in the registered in vivo US cases are comparable to those in the ex vivo MRcases. We believe this is due to the high resolution of the US scans that allowed a goodreconstruction of the imaged prostate’s surface, to which the slices are directly registeredusing the point-based metric.In most cases, the error ranges in Table 2.2 are below the clinical significance canceroustumor threshold of 0.2 cm3 in volume, as reported in [36], which translates to a sphericaltumor of 7.2 mm in diameter. Thus, for the intended application of cancer characterization,such errors should allow feature extraction and training of classification algorithms, e.g.,support vector machines [146] or random forests [19], with techniques to discard possibleoutliers caused by larger registration errors.The employment of a particle filtering framework allows us to contend with the highdimensionality of the constrained search space, and with the multimodal nature of theoptimization problem. The uncertainties in the alignment, imaging noise and surface errorare modeled stochastically, which allows derivation of optimal registration parameters in aBayesian approach.Other methods that employ PFs for registration typically assume a system noise associ-ated with a diagonal covariance matrix, i.e., such that the state parameters are statistically342.5. Discussion and Conclusionindependent. Here, in addition to the novel application of the algorithm to a multi-sliceto volume registration, we have incorporated a misalignment model into the PF frameworkthrough a parameterized structure of the covariance matrix. This covariance matrix allowsthe generation of slice poses that are correlated according to the misalignment mechanismsin histopathological processing.The incorporation of intensity-, region- and point-based similarity metrics into the reg-istration algorithm, in order to account for the different imaging modalities, demonstratesthe generality and flexibility of the proposed method, with which any similarity metric canbe employed. Specifically, we have employed the MI and Chan-Vese metrics on the ex vivoand in vivo MR data sets, respectively, based on their previously reported successful appli-cations to these modalities (see Section 2.3.5). In case these metrics do not perform welldue to corrupted or noisy data, the simple SSD between closest points that was employedon the US data sets can still be used as the metric for these modalities, at the expense ofan added segmentation or other processing of the volumetric image.The main limitation of the proposed method, which stems from its generality, is therelative high number of parameters. Tuning these parameters requires a preliminary studyof the variances of the misalignment, and/or a trial-and-error process. However, as discussedin Section 2.4.3, we found the algorithm to be robust within a range of parameters’ values,and in most scenarios they need to be set only once for processing data sets that wereacquired under the same clinical settings.Note that, due to its stochastic nature, running the algorithm on the same data set un-der the same parameters may still yield different results. However, in case the parametersare tuned properly, the variance among such results is not significant as they should all bearound the global minimum. To ensure a convergence to a global minima, a high number ofparticles should be employed, which makes the algorithm computationally expensive. How-ever, in a typical setting, the algorithm is running offline, after histopathology processing,and therefore registration time may be expendable.Another drawback of the current implementation is that the rigid transformation withanisotropic scaling we assume is unable to fully capture distortion caused by tear, shear,fold, or stretch, as assumed in [60] (a small percentage of slices in our data, on which glandsegmentation was approximated by the expert). A straightforward extension of our methodis to consider a full affine transformation by adding three additional shear parameters, oneper axis, with corresponding standard variations. However, in that case, additional particlesand/or iteration are expected to be required for convergence.The proposed framework can be also extended to accommodate a deformable registrationstep. Such a step can further improve the results by capturing any residual deformations,352.5. Discussion and Conclusionand is required, e.g., when using an endorectal coil for acquiring in vivo MR. Any deformableregistration algorithm can be employed based on the data, e.g., an intensity-based algorithmsuch as free-form deformations [130], or a point-based algorithm such as landmark-basedthin-plate splines [128]. We note that, due to the sparsity of slices, a 2-D registrationshould be applied between each histological slice and its corresponding registered slice onthe volumetric image.In case a 3-D deformation map between the specimen and in vivo imaging is required, apossible approach is to first register the histology to an ex vivo MR with minimal error usingthe proposed algorithm, and then compose the resulting transformation with a deformablemapping that was obtained from a 3-D elastic registration between the ex vivo MR and thein vivo MR, e.g., using the method described in Chapter 3. Future work may include theemployment of such nonrigid registration algorithms, either as postprocessing, or as part ofthe particle filtering framework.With the proposed similarity metrics, the approach is either susceptible to segmentationerrors, or to lack of intensity features. Future work may improve the usage of the MI metricby employing a Parzen window approximation [165] to allow a smooth minimization of theintensity-based metric. This will increase the robustness of the algorithm to the choice ofthe number of update iterations when using the intensity-based metric.Other future work includes investigation of more metrics that can describe the simi-larity between histological slices and volumetric imaging. Specifically, we are interested inemploying a point-based metric between descriptors of intensity-based extracted features.Incorporated into our PF framework, such metric will allow avoiding segmentation, as whenintensity-based metrics are used, while benefiting from the shorter runtime of a point-basedmetric. An automatic extraction of salient features can be based, e.g., on a scale invariantapproach, as proposed in [175], or on the modality independent neighbourhood descriptor[55] that has recently gained popularity for multimodality registration of medical imaging,including the prostate [156]. Another research direction is to employ a linear combinationof metrics as the similarity, and either train or manually tune their weights in order toprovide a better shaping of the target functional according to the data.36Chapter 3Model-Based Registration of theProstate Using Elastography3.1 IntroductionAs discussed in Chapter 2, a typical sparse sectioning of histological slices does not allow adirect construction of a 3-D deformation map between the specimen and in vivo imaging.Recently, researchers have proposed the use of an ex vivo MRI of the prostate specimen[167, 168] as an intermediate modality between histopathology and in vivo MRI. Performedafter excision and fixation, and just before sectioning, such ex vivo scan captures the shape ofthe deformed specimen, and can therefore be better aligned with the histological slices usingonly a rigid registration. Compared to the histological slices, the superior axial resolutionof the ex vivo scan provides true 3-D modeling of the specimen that facilitates a nonrigidregistration to the in vivo MRI.Thus, the slice-to-volume registration of histological sections to in vivo MRI can bedivided into two more tractable registrations of histopathology to ex vivo MRI, and ex vivoto in vivo MRI. The former registration can be solved with minimal error by approachesproposed, e.g., in [45] and [70], or by employing a multi-bladed cutting device [32] thatmechanically constrains the prostate specimen during sectioning. For the latter registration,however, the authors in [167] propose only an interactive rigid registration between thein vivo and ex vivo volumes, which has been recognized to be a shortcoming [168].It is the objective of the work described in this chapter to develop and validate a methodfor an automatic and accurate 3-D registration between ex vivo and in vivo MRI. Figure 3.1illustrates the problem with typical images to be registered. We propose a novel deformableregistration approach, in which the ex vivo MRI volume of the prostate is segmented andtreated as a reference model, and the in vivo MRI volume is treated as a deformable templatewith an underlying patient-specific biomechanical model. First, a similarity transformation(translation, rotation and isotropic scaling) is employed to align the ex vivo and in vivoscans and correct for the change in volume. Following this initial alignment, an elasticregistration is performed to match the template with the reference. Similar to [76], the373.1. Introduction(a) In vivo T2w (b) MRE elastogram⇔(c) Ex vivo T2w (d) Prostate specimenFigure 3.1: Problem illustration. (a)-(b) The in vivo MRE scan provides measured elasticityfor each imaged voxel. (c)-(d) A 3-D model of the postoperative specimen, acquired fromex vivo MRI, should be registered to the in vivo MRI. Cross-sections of the unregisteredex vivo model (red) and the unknown in vivo prostate surface (cyan) are overlaid on slices.deformation is driven by a region-based energy as the distance measure, and the elasticpotential as a regularizer. Unlike [76], we deform the image in order to fit the model, anduse an inhomogeneous biomechanical model, acquired from MRE [100], in our regularizationterm.While potentially faster because of the much coarser discretization of images, biome-chanical models derived from the FEM, such as those employed in [16] and [5], requiresegmentation and meshing of the internal structure of the prostate and its surroundinganatomy. These steps are time consuming, often require user intervention and may gener-ate errors that affect registration performance. Unlike FEM-based approaches, our methodallows to solve the registration problem on a regular grid by using a variational approachand eliminates volume meshing. Such an approach was successfully employed in [13] for USimage registration in prostate biopsy tracking.The MRE technique that we use integrates relatively easily in a standard multi-parametricin vivo MRI session. It involves the transperineal application of single- or multiple-toneharmonic mechanical excitation that generates shear waves in the prostate. The waves aremeasured by encoding the resulting displacements in the phase component of a synchro-nized MRI pulse sequence. Images, known as elastograms, that depict tissue elasticity orstiffness are generated from the imaged tissue displacements. With MRE, each voxel in thein vivo MRI is associated with the measured elasticity of the corresponding voxel in theelastogram. Additional details about the technique are given in [134] and Section 3.5.1.With few exceptions, e.g., [29, 64, 78], biomechanical models have been assigned arbi-trary and typically constant elastic properties that are optimized to produce low registrationerrors. Even methods that do consider inhomogeneous elasticity for registration, express383.2. Preliminariestissue inhomogeneity by setting values that are empirically found to work well. This is notalways realistic and may result in inaccurate deformation maps.In general, the current use of elastography in image processing is limited (see [140]for a review on applications in prostate imaging), and, in particular, the use of patient-specific elasticity maps for registration is new. We justify our approach by providing anovel quantitative analysis that shows the advantage in utilizing inhomogeneous elasticityfor registration. We validate the method by presenting the results from synthetic and sixpatient data sets.The remainder of this chapter is organized as follows. Section 3.2 describes the theoryand basic ideas behind the proposed methods. Section 3.3 presents a method for registeringthe ex vivo model to an in vivo volume using elastography. In Section 3.4 we investigatethe registration performance with respect to elastic inhomogeneity. In Section 3.5, weprovide experimental results of the proposed registration method on clinical data. Finally,we conclude in Section 3.6 with a discussion and future research directions.3.2 Preliminaries3.2.1 Elastic RegistrationWe denote the 3-D spatial coordinate vector as r¯ = (x, y, z)T ∈ R3. Let IR, IT : D → Rdenote the reference and template volumes, respectively, on a 3-D domain D ⊆ R3. Also,let u¯ = (ux, uy, uz)T : D → R3 be a displacement field, such that the deformed templatewith respect to u¯ isI˜T (r¯; u¯) = IT (r¯ − u¯(r¯)). (3.1)The elastic registration problem is finding an admissible displacement field u¯, such thatthe deformed template I˜T is similar, in some sense, to the reference IR. This problem isformulated as the minimization of the joint functional, or energy,J [u¯; IR, IT ] = F [IR, I˜T (·; u¯)] + S[u¯] , (3.2)where F represents a distance measure between the reference and the deformed template,and S is a metric that determines the smoothness or regularization of the displacements.Common choices for a distance measure are SSD, NCC and (normalized) MI [165].In [76], the authors proposed a region-based energy (see Section 3.2.2) as the distancemeasure for a combined segmentation and registration. Their method maps a contour,which is extracted from a segmented template image, onto the reference image to provide393.2. Preliminariesthe segmentation of the reference and its spatial correspondence with the template.A variety of metrics have been proposed to regularize the deformations based on physicalmodels, e.g., diffusion [58], fluid [25] and curvature [41] models. The authors in [20] and[11] proposed the linear elastic potential of the displacements as a regularizer, namelyS[u¯] =∫Dµ43∑j=13∑k=1(∂u¯j∂r¯k+ ∂u¯k∂r¯j)2 + λ2 (div u¯)2 dr¯ , (3.3)where, from continuum mechanics, λ and µ are the Lame´ parameters that reflect the elasticproperties (µ is also known as the shear modulus).We follow the general framework presented in [92], and provide a variational formulationof the nonrigid registration problem and the corresponding Euler-Lagrange equations thatcharacterize a minimizer. In the variational approach, in order to minimize the functional(3.2), the Gaˆteaux derivative is computed with respect to u¯. This yields the Euler-LagrangeequationA[u¯](r¯)− f(r¯, u¯(r¯)) = 0 , (3.4)where −f , also known as the (external) force, corresponds to the derivative of F , and A[u¯],also known as the internal force, corresponds to the derivative of S and typically representsa partial differential operator. Thus, a minimizing displacement map should satisfy thisequilibrium equation.In the case of elastic potential as a regularizer, the Gaˆteaux derivative of (3.3) yieldsthe Navier-Lame´ operatorA[u¯](r¯) = µ∆u¯(r¯) + (λ+ µ)∇div u¯(r¯) , (3.5)for which the Euler-Lagrange equation, (3.4), is a second-order nonlinear PDE, also knownas the elastostatic or Navier-Cauchy equation.3.2.2 Region-Based Active ContoursActive contour models [66] have been extensively employed for image segmentation. In thesemethods, an initial segmenting surface is evolved in order to minimize an energy functional.Edge-based (classical) active contours, [22, 67], drive the surface to lock onto local maximaof image gradient magnitude values that typically characterize edges. The evolution of thesurface in this case depends strictly on nearby pixels, and is therefore local. Region-basedactive contours, [23, 174], rely on regional statistics, such as sample mean and variance403.2. Preliminarieswithin the image, for the evolution of the segmentation surface, and direct its movementtoward boundaries that are not necessarily defined by clear edges. Region-based models areglobal and tend to be robust to noise, since they take into account intensities within entireregions, rather than intensities of individual pixels.The level set method, [109, 147], is an effective technique for implicit representationof evolving surfaces, frequently employed to implement active contours algorithms, sinceit inherently allows for cusps, corners and automatic topological changes. In the level setmethod, a surface (or contour in 2-D) is embedded as the zero level set of a function Φ.We will use the convention that Φ takes negative values in the region inside the surfaceand positive values in the region outside. A common choice for a level set function is thesigned distance function, of which the value Φ(r¯) at each point r¯ equals the Euclidean (`2)distance to its nearest point on the surface, multiplied by 1 or −1 according to the side ofthe surface r¯ lies in.Formulating energy functionals and their variations in the level set method often requiresusage of the Heaviside function and the Dirac delta function, defined respectively byH(z) =1, z ≥ 0;0, z < 0,and δ(z) = ddzH(z).In practice, regularized (smoothed) versions of H and δ are used, as suggested in [23], inorder for the minimization process to be less local, thus increasing the chances of reachinga global minimizer, independent of the initial surface.A popular region-based active contour model for image segmentation was proposed byChan and Vese in [23]. In its level set formulation, the Chan-Vese algorithm minimizesthe following energy functional with respect to the level set function Φ that embeds thesegmenting surface,FCV [Φ; IT ] =∫D(IT − c−)2H(−Φ) dr¯+∫D(IT − c+)2H(Φ) dr¯ , (3.6)where c− and c+ are the mean intensities of the image to be segmented, IT , in the regionsinside and outside the segmenting surface (i.e., the zero level set of Φ), respectively, givenbyc∓[Φ; IT ] =∫D ITH(∓Φ) dr¯∫DH(∓Φ) dr¯. (3.7)413.2. PreliminariesThe minimization of (3.6) can be understood as partitioning the image into two homo-geneous regions, by minimizing the variances of the intensities on the regions inside andoutside the segmenting surface.3.2.3 Model-Based AlignmentA nonrigid registration typically follows a rigid registration that corrects for the misalign-ment and change of coordinate systems between the two images. In this work we employa model-based alignment based on an algorithm that was originally used in [163] for ashape-based segmentation.Let ΦR be a level set function that embeds the surface of a reference model, which wasextracted by a segmentation of the reference image IR. A similarity transformation in 3-Dis defined by seven pose parameters: three translation values along each axis, three rotationangles about the axes, and one isotropic scale factor. Using a linear transformation of thespatial coordinates, the transformed model may be written asΦ˜R(r¯; ∆, θ, σ) = ΦR(σ−1Θ−1(r¯ −∆)), (3.8)where ∆ and θ are vectors that contain the translation and rotation parameters, respectively,σ is the scale, and Θ is the rotation matrix, which is determined by θ.The authors in [163] propose to segment medical images by optimally fitting a modelonto the image. They extract the principal modes of a training set that comprises differentshapes of the target object, and obtain a parametric shape model by tuning the shapeparameters, namely the coefficients in a linear combination of the principal modes. Theirmethod minimizes region-based energy functionals in a level set formulation, among themthe Chan-Vese energy in (3.6), with respect to the pose and shape parameters.We consider a single element training set, in which the shape is the surface of thereference model. Since the minimization in this case is performed only with respect to thesimilarity pose parameters, the method is degenerated into alignment of the model. Thefunction to be minimized is thereforeFCV (∆, θ, σ) = FCV [Φ˜R(·; ∆, θ, σ); IT ]. (3.9)The alignment algorithm employs gradient descent to minimize (3.9) by computing itsderivatives with respect to the pose parameters. We refer the interested reader to [163] fordetails.423.3. Model-Based Registration Using Inhomogeneous Elasticity3.3 Model-Based Registration Using InhomogeneousElasticityIn our framework, high quality ex vivo images with a homogeneous background (as seen inFigure 3.1(c)) allow a 3-D model of the prostate to be constructed easily by using eithermanual, semi- or fully automatic segmentation. In contrast, the in vivo images containsurrounding anatomy and tissue, with which the prostate blends. Thus, in our model-basedregistration, we take the ex vivo model as the reference, and the in vivo image as thetemplate.3.3.1 Model-Based ForceMotivated by the adequacy of region-based energies to delineate the prostate in MRI, asdemonstrated in [163], we follow a model-based registration approach that incorporates aregion-based energy as the distance measure between the template image and the referencemodel. Specifically, we propose a distance measure which is based on the Chan-Vese energyin (3.6),F [u¯] = 12FCV [ΦR; I˜T (·; u¯)], (3.10)where, as in Section 3.2.3, ΦR is a level set function that embeds the surface of the referencemodel. Note that in this formulation, the “segmenting” surface ΦR is constant and theminimization is with respect to the displacement field.Thus, in contrast to active contours segmentation methods and the registration approachin [76], here the image is deformed in order to fit the model, rather than the model surfaceto fit the image. This formulation allows for the inhomogeneous elasticity of the voxels inthe template to regularize its deformations, as explained below.By computing the Gaˆteaux derivative of (3.10), we have that the force associated withthis energy isfCV (r¯, u¯(r¯)) = − [(IT (r¯ − u¯(r¯))− c˜−)H(−ΦR)+ (IT (r¯ − u¯(r¯))− c˜+)H(ΦR)]∇IT (r¯ − u¯(r¯)) , (3.11)where, analogous to the discussion in Section 3.2.2, c˜− and c˜+ are the mean intensities ofthe deformed template image in the regions inside and outside the reference model.Being derived mathematically, the magnitude of the forces acting on the template couldbe arbitrarily high and may force unrealistic deformation to bring it into the reference.However, using a scaling factor, we can normalize the forces such that their maximum433.3. Model-Based Registration Using Inhomogeneous Elasticitymagnitude inside the model is in the order of 0.1 N, which is physically plausible andwithin the same order of magnitude as the gravitational force on the prostate [157].3.3.2 Inhomogeneous RegularizationAs discussed in Section 3.1, in most elastic registration algorithms, the Lame´ parametersλ and µ are constant and chosen arbitrary. In elastography, we obtain measurements ofthe Young’s modulus that is associated with each voxel of the template, ET : D → R. Thecorresponding Lame´ parameters in this case are the functions given byλT =ET · ν(1 + ν)(1− 2ν) and µT =ET2(1 + ν) , (3.12)with an (almost) incompressible Poisson ratio ν = 0.499.Thus, we modify (3.5) to contend with the inhomogeneous Lame´ parameters,AT [u¯](r¯) = µT (r¯ − u¯(r¯))∆u¯(r¯)+ (λT (r¯ − u¯(r¯)) + µT (r¯ − u¯(r¯)))∇div u¯(r¯) . (3.13)Note that λT and µT are deformed according to the template so as to “follow” its vox-els. Also note that (3.13) is only a piecewise constant approximation of the problem. Thenon-approximated operator, derived by taking the Gaˆteaux derivative of (3.3) with inho-mogeneous Lame´ parameters (see, e.g., [29] and [64]), contains terms with gradients of theLame´ functions. Due to the noisy nature of elastography, numerical calculation of thesegradients is unreliable and the minimization is unstable. However, the approximation isjustified by a typical “slow” spatial variation of the true elasticity (another option is toemploy a low-pass filtered elastogram and use the non-approximated operator).3.3.3 Energy MinimizationIn order to solve the PDE (3.4) for the displacements u¯, with A = AT and f = fCV , wefollow a fixed-point iteration scheme [92], in which an initial (e.g., zero) displacement u¯(0)is assumed. We then solveAT [u¯(k)](r¯) = fCV (r¯, u¯(k−1)(r¯)), k = 1, 2, . . . , (3.14)for u¯(k), until convergence. Since the elastic model relates to small (infinitesimal) deforma-tions, we chain the resulting displacements at each iteration, as if the deformed templateat each iteration is to be registered to the reference.443.3. Model-Based Registration Using Inhomogeneous Elasticityx2(a) Initialx2(b) Final (homogeneous)x2(c) Final (inhomogeneous)Figure 3.2: Registration of synthetic images. The reference contour (red) is overlaid on(a) the initial binary template, and (b)-(c) the warped template after registration usinghomogeneous/inhomogeneous elasticity. The grid (blue) illustrates the displacements.On a discrete grid, the technique is to first discretize AT and fCV using finite differences.Next, we generate a sparse convolution matrix AT that corresponds to the operator AT afterdiscretization, such that (3.14) readsATu(k) = f (k−1)CV , (3.15)where fCVand u are column vectors that contain samples of the forces and (unknown)displacements on the discrete grid, respectively. The samples are ordered lexicographicallyby stacking the voxels of each component on top of one other to form one large column.Zero boundary conditions are assumed to compute the convolution matrix.3.3.4 Synthetic ExampleThe proposed method is illustrated on a synthetic example in 2-D. Here, two binary imagesare being registered. The template is an ellipse with intensity values of “1” assigned to itsinside and “0” outside. The reference is a similar but dented ellipse, of which the contouris used as the model.In Figure 3.2(a), the reference contour is overlaid on the initial template image. First,homogeneous elasticity was assigned to the template. Figure 3.2(b) shows the final warpedtemplate (after 50 iterations). Notice the symmetry of the warped grid around the horizontalaxis. In the second example, the bottom half of the template was assigned three times stifferelasticity values. Figure 3.2(c) shows the final warped template in this case. As expected,the displacements are larger in the top half of the image.453.3. Model-Based Registration Using Inhomogeneous Elasticityx2(a) Initialx2(b) Final (homogeneous)x2(c) Final (inhomogeneous)Figure 3.3: Registration of noisy synthetic images. The reference contour (red) is overlaidon (a) the initial binary template with a Gaussian noise of 0.25 STDV, and (b)-(c) thewarped template after registration using homogeneous/inhomogeneous elasticity. The grid(blue) illustrates the displacements.Due to the employed region-based force, which takes into account region statistics ratherthan specific pixel values, the method is robust to noise. In Figure 3.3 we see the sameregistration examples as before, but with a zero mean Gaussian noise added to the binarytemplate image. The corresponding results (after 100 iterations) are shown in Figures 3.3(a)-(c). Indeed, the results are very similar to the noiseless experiments (a slight grid asymmetrycan be noticed in Figure 3.3(b)).In order to quantify the robustness to noise, we repeat the experiments for a Gaussiannoise with different values of STDV. We evaluate both robustness to intensity variations, byadding the noise to the template, and robustness to elasticity variations, by adding the noiseto the elasticity map. The resulting pixels’ displacement map for each noisy registrationis compared to the noiseless registration map. Figure 3.4 plots the mean Euclidean norm(over all pixels) of the difference between the displacement maps against the correspondingSTDV values of the noise.The intensity noise and elasticity noise affect the results to a small extent (less thanone pixel of error on average), up to STDV values of 0.1 (10% of noiseless intensity andelasticity). We see that intensity noise causes a constant error in the displacement. Thisis mostly due to small displacement errors in the background, and not due to errors in thelarge displacements around the object. The elasticity noise however, causes a logarithmicincrease of the error. This is due to violation of the “slow” varying assumption in theapproximation of the Navier-Cauchy equation, (3.13).For STDVs above 0.1, the error caused by intensity noise increases, while elasticity463.3. Model-Based Registration Using Inhomogeneous Elasticity10−3 10−2 10−1 10010−210−1100Noise Standard DeviationMeanDisplacementError(pixels) Intensity NoiseElasticity NoiseFigure 3.4: Robustness to noise. The error is the mean Euclidean distance between thedisplacement maps resulting from noiseless and noisy registration of synthetic images, asseen in Figures 3.2 and 3.3. The STDV of the noise represents the (inverse) SNR withrespect to the ellipse (intensity and elasticity values of “1”).noise yields convolution matrices AT that are singular (to working precision), due to pixelswith (almost) zero elasticity, causing the algorithm to fail. We also found a logarithmicrelationship between the noise and the required number of iterations for convergence.3.3.5 Clinical ExampleHere, we describe the utilization of the proposed method for registration between an ex vivomodel and the in vivo MRI volume and demonstrate the process on clinical data. As apreprocessing step, we translate the ex vivo model to match its center with the center ofthe in vivo volume. Next, we employ the model-based approach described in Section 3.2.3for fine alignment with respect to translations, scaling and rotations. Finally, we computethe residual nonrigid map between the aligned model and the image using the proposedmodel-based registration with inhomogeneous elasticity.The registration process is illustrated in Figure 3.5. Note that in Figure 3.5(c), thetemplate image is warped to fit the prostate inside the reference model. The inverse map isthen applied to the rigidly aligned model in order to produce the registered ex vivo modelon the coordinate system of the in vivo image, as seen in Figure 3.5(d).We note that in [163], the authors transformed the T2w images by applying the gradient473.4. Elastic Homogeneity Versus Inhomogeneity(a) Initial (b) Aligned (c) Warped template (d) Registered referenceFigure 3.5: Registration process. (a) The initial ex vivo reference model (red) is overlaid ona slice of the in vivo T2w template image. (b) Translated, rotated and scaled model afterusing the model-based alignment. (c) Warped template after using the proposed model-based elastic registration. (d) Registered model after applying the inverse map. A 3-D view(top row) of the model, and a 2-D cross-section (bottom row) are shown.operator in order to increase intensity homogeneity inside the prostate and accentuate itsboundaries for their region-based segmentation to work. In our case, we found the registra-tion to perform better on the MRE magnitude images (see Section 3.5.1), and superimposedthe results on the T2w images.3.4 Elastic Homogeneity Versus InhomogeneityIn this section we provide analysis to show that inhomogeneous elasticity improves theregistration performance. In order to do so, we degenerate the MRE measured elastogram,ET , by employing the k-means clustering algorithm [33] with different values of k’s.The feature space to be clustered comprises the elasticity value, with the addition ofthe spatial coordinates r¯, in order to create geometrically connected clusters. The “points”to be clustered are therefore the feature vectors generated for each voxel. A clusteredelastogram ET,k is created by applying the elasticity value of the centroid to each voxel inthe corresponding cluster.The extreme case of k = 1 yields a homogeneous elastogram, ET,1 = const, in whichthe intensity is the mean elasticity value. The other extreme, which we denote by k →∞,yields the original inhomogeneous elastogram, ET,∞ = ET . We therefore refer to k as theinhomogeneity parameter.483.4. Elastic Homogeneity Versus Inhomogeneity(a) ET,1 (homogeneous) (b) ET,3 (c) ET,10 (d) ET,50(e) ET,100 (f) ET,250 (g) ET,1000 5101520253035404550(h) ET,∞ (totalinhomogeneous)Figure 3.6: Elastograms for different values of the inhomogeneity parameter. Intensity scaleis in kPa.In fact, any value k ≥ N , where N is the total number of voxels, and, in practice, dueto repetitions in the quantized data, even lower k’s, yield the same clustering as k → ∞.Figure 3.6 depicts a mid-gland slice of the elastograms, which correspond to the examplecase in Figure 3.5, for selected values of the inhomogeneity parameter in the range betweenhomogeneity and total inhomogeneity.We repeatedly perform the registration described in Section 3.3.5 for different values ofthe inhomogeneity parameter, with respect to the resulting clustered elastogram. We runeach registration until convergence to the same energy value that we have achieved in theoriginal experiment. Therefore, the resulting warped templates look similar to each otherand to the result in Figure 3.5(c). However, as expected, the displacements that map eachof the templates are different.Figure 3.7 shows the Euclidean norms of the displacement difference between the result-ing maps that correspond to the inhomogeneity parameter values shown before. Specifically,Figure 3.7(h) shows the difference between the total inhomogeneous and homogeneous reg-istration maps on the mid-gland slice, and we may notice a difference of ∼ 3.5 mm insidethe prostate region.The performance is evaluated for each registration through the volume overlap, in the493.4. Elastic Homogeneity Versus Inhomogeneity(a) ‖u¯k=3 − u¯k=1‖ (b) ‖u¯k=10 − u¯k=3‖ (c) ‖u¯k=50 − u¯k=10‖ (d) ‖u¯k=100 − u¯k=50‖(e) ‖u¯k=250 − u¯k=100‖ (f) ‖u¯k=1000 − u¯k=250‖ (g) ‖u¯k→∞ − u¯k=1000‖ 0.511.522.533.544.55(h) ‖u¯k→∞ − u¯k=1‖Figure 3.7: Euclidean distance of the displacement maps resulting from registrations usingdifferent values of the inhomogeneity parameter. Intensity scale is in mm.sense of Dice’s coefficient (3.16), between the registered ex vivo model and a manuallysegmented in vivo model, and through the distance between manually selected landmarkson the two modalities (see Section 3.5.3). Figure 3.8 shows these measures for differentvalues of the inhomogeneity parameter.Indeed, we identify the general trends that the volume overlap and the landmark errorare, respectively, in direct and inverse proportion with the inhomogeneity parameter. Wetherefore conclude that performance improves when an inhomogeneous model is used. Wenote, however, that for two cases, the maximum volume overlap and minimum error wereattained at k = 250 and k = 500 rather than k → ∞. The reason for that is the noisynature of the inhomogeneous elastogram that violates the “slow” varying assumption, asdiscussed in Section 3.3.4.Thus, setting a value for k that is too small leads to elastic homogeneity, while setting itsvalue too high may lead to inaccurate regularization. In our experiments we take k = 250,since for greater values, when applicable, the improvement is not significant.503.5. Experiments100 101 102 103 1047075808590VolumeOverlap(%) Volume OverlapMean Landmark Error 2.42.62.833.2MeanLandmarkError(mm)InhomogeneityFigure 3.8: Registration performance for different values of the inhomogeneity parameter.3.5 ExperimentsIn this section we describe the acquisition of the data sets that were used in our study. Weprovide details on the imaging protocol, as well as on the implementation of the registrationalgorithm, its execution time and evaluation. Finally, we summarize the results of theregistration experiments.3.5.1 Data AcquisitionIn Vivo ScanAfter approval of the institutional ethics board and signed consents, six patients (mean age65, range 59–76 years old) who were scheduled for radical prostatectomy, underwent a preop-erative in vivo T2w scan followed by MRE scan in a 3-Tesla system (Achieva 3.0T, Philips,The Netherlands). A standard 6-channel cardiac coil with acceleration factor (SENSE) of2 was used.First, a sagittal scout image was acquired to ensure that the transperineal MRE driveris properly positioned at the beginning of the examination, as illustrated in Figure 1.3(a).Then, a standard axial T2w FSE sequence was acquired (TE/TR = 86/2500 ms, FOV= 320× 320× 70 mm3, with 0.5 mm in-plane resolution and 2 mm slice thickness).Next, a transperineal prostate MRE scan was acquired using a gradient echo sequencenamed eXpresso [136]. Eight phases were encoded with a mechanical excitation of 70 Hz513.5. ExperimentsMin–max inside Mean inside Min–max outside Mean outsideCase 1 4.8− 50.3 19.5± 9.8 1.6− 75.3 13.1± 6.0Case 2 3.0− 36.6 17.7± 5.3 1.4− 79.6 11.1± 5.2Case 3 2.0− 38.0 13.7± 6.8 1.0− 36.4 8.5± 4.5Case 4 1.6− 26.0 11.6± 3.6 1.3− 53.1 9.7± 5.3Case 5 6.6− 37.0 21.0± 6.7 1.9− 79.2 15.8± 9.0Case 6 10.3− 50.8 24.0± 6.5 2.8− 54.5 15.0± 8.2Table 3.1: MRE results. Quantitative summary of elasticity inside and outside the prostate.Values are in kPa.applied to the perineum of the patient [134]. Wave images were acquired on a 128×128×24matrix with 2 mm isotropic voxel size. Figure 1.3(d)-(f) shows an example slice of these waveimages. The entire MRE scan lasted about 8 min for a 3-D displacement field acquisition.Each slice was cropped to accommodate 64 × 64 pixels around the prostate region forcomputational efficiency, and then processed offline to generate an elastogram similarly tothe approach described in [150].The quantitative elasticity values that were obtained are summarized in Table 3.1. Therange and mean elasticity were computed over the regions inside and outside (excluding thepubic bone) of the prostate boundary (manually chosen by an expert for the registrationevaluation). The results demonstrate the inhomogeneity in elasticity and the fact that, ingeneral, the prostate is stiffer (higher values) than its surrounding tissue. For other reportedelasticity values of the prostate, see [134] and references therein.In addition to the elastogram, an image was generated from the magnitude componentof the MRE signal. This image, referred to as the magnitude image, is intrinsically alignedwith the elastogram and provides complementary intensity information. An example slice ofthe magnitude image in Figure 1.3(c) shows the intensity homogeneity inside the prostate,which is ideal for a region-based algorithm to be employed on. The scanning geometry,as given in the scanner output, determines the alignment of both the elastogram and themagnitude image with the T2w volume.Ex Vivo ScanThe postoperative ex vivo prostate specimen was fixed in 10% buffered formalin, typicallyfor 72 hours. It was then wrapped in a plastic bag, taped to a dedicated holder and scannedin a 7-Tesla system (Biospec, Bruker, Germany).A rapid acquisition with relaxation enhancement (RARE) pulse sequence was used to523.5. Experimentsacquire axial T2w images (TE/TR = 70/5000 ms, FOV = 70× 70× 42 mm3, with 0.2734mm in-plane resolution and 2 mm slice thickness).To extract the ex vivo model, the images were semi-automatically segmented by em-ploying the region-based active contours algorithm described in Section 3.2.2, with manualcorrections by a radiologist that were required mainly due to artifacts caused by imagingof the wrapping bag around the specimen. We have then used Stradwin (Cambridge Uni-versity, U.K.) [160], which employs [161], in order to generate a surface from the segmentedex vivo slices. We set the surface resolution and smoothing strength to “medium.”3.5.2 Implementation DetailsWe have implemented our registration method and tested it under MATLAB on a 2.93GHz Quad Core CPU machine with 8 GB of RAM. It typically takes 100 iterations of rigidregistration, and additional 100 iterations of nonrigid registration to reach convergence ofthe energy. The bottleneck of the algorithm is solving (3.15) for the displacement values. Itinvolves finding 3×N unknowns in each iteration, where N is the total number of voxels.This is usually a large number and makes the algorithm computationally expensive.In general, the inhomogeneity of the Lame´ parameters yields sparse, but non-symmetricand non-periodic convolution matrices AT , which do not allow numerical schemes to takeadvantage of their structure, e.g., as in [40] where the authors employ the fast Fouriertransform (FFT) for fast inversion. Thus, an approximate solution of the problem is foundby employing the biconjugate gradient method [123] for a fixed number of iterations, whilekeeping execution times within reason. We found 50 iterations to provide plausible results,which resulted in a total execution time of about 12 min for registering 64×64×24 volumeson an un-optimized code.3.5.3 ResultsTo evaluate the registration performance, we compare the final registered ex vivo model toan in vivo prostate model that was segmented manually by an expert on the T2w slices. Asa general performance indicator, we measure V O, the volume overlap between the registeredand manual 3-D models in the sense of Dice’s coefficient, i.e., twice the size of the intersectedvolume divided by the sum of the sizes of the volumes. In the level set method, this can beformulated asV O = 2∫DH(−Φ˜R)H(−ΦT ) dr¯∫D(H(−Φ˜R) +H(−ΦT ))dr¯, (3.16)533.5. Experimentswhere Φ˜R and ΦT are the level set functions embedding the registered (warped) model andthe manual model, respectively. Similarly, in order to evaluate the performance at differentregions of the prostate, we measure AO, which is the area overlap between cross-sections ofthe registered and manual models on slices from the base, mid-gland, and apex regions ofthe gland.In addition, the expert has identified matching landmarks inside the prostate regionon both the ex vivo and in vivo images. We measure the Euclidean distance between theregistered ex vivo landmarks and their in vivo counterparts. Since our goal is to register(unknown) cancerous regions inside the prostate, these distances between landmarks provideus a good approximation of the target registration error (TRE).Finally, in order to evaluate the improvement that was achieved by employing MRE,we repeated the experiments for all cases after assigning the voxels a homogeneous elas-ticity value that is equal to the mean elasticity value inside the prostate (the elastogramregion inside the expert’s manual segmentation). Note that this experiment is different thanthe simulation with k = 1 in Section 3.4, where the mean elasticity of the entire field ofview was used, and is expected to perform better since the elasticity inside the prostate ismore relevant. We compare these homogeneous registration results to the inhomogeneousregistration results.Figure 3.9 visualizes the inhomogeneous registration results. In Figure 3.9(a), the sur-faces of the registered and manual models of each case are shown in 3-D, together with theregistered and manual landmarks. Figures 3.9(b)-(d) show cross-sections of the registeredand manual models overlaid on selected slices from the base, mid-gland and apex regions ofthe prostate, respectively. Quantitative and statistical summaries of the registration resultsare presented in Table 3.2 and Figure 3.10, respectively.543.5. Experiments(a) 3-D (b) Base (c) Mid-gland (d) ApexFigure 3.9: Inhomogeneous registration results. The registered ex vivo reference model (red) anda manually segmented model (cyan) of the in vivo T2w template image. Spheres in column (a)represent registered and manually marked landmarks. Columns (b)-(d) show cross-sections of themodels overlaid on selected T2w slices from the base, mid-gland and apex regions of the prostate.553.5.ExperimentsInhomogeneous Elasticity Homogeneous ElasticityBase AO (%) Mid-gland AO (%) Apex AO (%) VO (%) Error (mm) VO (%) Error (mm)Case 1 77.7± 4.5 91.0± 2.3 88.3± 0.8 87.5 2.5± 0.8 85.3 3.1± 1.1Case 2 76.6± 10.9 92.4± 2.6 78.1± 7.8 85.9 1.5± 0.6 84.0 2.8± 0.7Case 3 90.3± 0.6 92.1± 3.0 78.9± 2.4 87.3 3.2± 1.0 82.9 3.1± 0.8Case 4 83.5± 6.0 83.5± 4.7 76.0± 2.9 81.6 4.3± 2.5 76.9 5.2± 3.6Case 5 90.0± 4.9 93.5± 2.5 89.9± 1.0 88.1 2.3± 1.3 84.1 2.5± 1.2Case 6 85.7± 2.3 91.2± 2.4 79.4± 9.3 86.9 4.1± 1.4 82.9 5.5± 2.3Mean 83.4± 7.7 90.8± 4.0 81.8± 7.1 86.2± 2.4 3.1± 1.4 82.7± 3.0 3.72± 1.9Table 3.2: Registration results. Quantitative summary.563.6. Discussion and Conclusion5075100AreaOverlap(%) AO (All slices)Base AOMid-gland AOApex AOLandmark Error 012345678910LandmarkError(mm)Case 1 Case 2 Case 3 Case 4 Case 5 Case 6Figure 3.10: Inhomogeneous registration results. Statistical summary of the resulting areaoverlap (AO) among slices, and the landmark error.3.6 Discussion and ConclusionWe have presented a method for an automatic model-based registration of ex vivo and in vivoprostate MRI. The method can be fused with established methods that consider registrationof ex vivo MRI and histological sections, in order to recover the spatial correspondencebetween histopathology and the in vivo scan. Since histopathology of the prostate providesthe ground truth cancerous and noncancerous regions, a correct mapping of these regionsonto an MRI volume allows characterization of CaP on co-registered multi-parametric MRI,and may provide training and validation sets for a classifier that produces cancer probabilitymaps on preoperative images. Typical training sets may include 10-20 cases, and validationsets may include 50-100 cases. Focusing on regions with high probability of cancer couldimpact treatments such as brachytherapy and radical prostatectomy, and reduce possibleside effects such as impotence and urinary incontinence.The proposed method aligns a model of the prostate, constructed from the ex vivo MRI,onto the in vivo volume, and then warps the in vivo volume so as to fit the aligned model.A variational approach has been used to derive the registration algorithm using a region-based active contour model in a level set formulation. The method was first demonstratedon synthetic images to evaluate its robustness to noise in the data, and was further testedand evaluated on six clinical cases by comparing the registered ex vivo model to manualsegmentations of the in vivo volume, and by computing the distances between registeredlandmarks on the two modalities.573.6. Discussion and ConclusionOur approach employs an MRE scan that allows measurement and assignment of elas-tic properties of the tissue to corresponding voxels in the image. In turn, this providesan implicit modeling of the prostate and its surrounding anatomical structures that allowswarping the in vivo volumetric image onto the ex vivo model in a physical manner. More-over, we found the MRE magnitude image to be an excellent candidate for the region-basedapproach due to its intensity homogeneity within the prostate.While the most valuable information is pertaining to the prostate itself, utilizing theelasticity of both the prostate and its surrounding tissue is more likely to result in a mappingthat portrays the actual displacement of the voxels within the prostate. This is attributableto the fact that most deformation occurs in vivo (ex vivo deformation is typically manifestedas shrinkage and corrected by scaling), and is therefore governed by the in vivo elasticityand interaction of the prostate with its surrounding anatomy. Indeed, we have showed thatthis approach yields improved registration as opposed to applying a constant elasticity valueto the volume, and that, in general, the performance improves with the inhomogeneity ofelasticity inside and outside of the prostate. To the best of our knowledge, no such study hasbeen conducted thus far. We note that an ex vivo measurement of the elasticity, althoughmight be technically simpler, reflects the much stiffer elasticity of the post-fixation specimenwhich is irrelevant to the deformation.In addition, our approach does not require segmentation and meshing of the in vivoanatomy, as opposed to explicit FEM-based modeling approaches. Although, as discussedabove, the designated application of classification requires a limited number of cases, seg-mentation and meshing can still consume a significant amount of time. On the other hand,MRE is fast and easily integrated into the in vivo MRI protocol. Even in case that FEMis to be employed instead of our approach, MRE could benefit the chosen algorithm byproviding the model with patient-specific elasticity.The performance of the registration, and the significance of its improvement using MRE,should be evaluated with respect to the designated application. The mean registrationerror that was achieved during the experiments is within the tumor diameter that MRI isreported to detect and characterize CaP [71], and within the size of most clinically significantlesions that are marked during histopathology analysis (∼10 ± 5 mm [34]). Thus, in ourcase, the method is suitable for mapping of these tumors on the in vivo MRI. However,the registration might not be accurate enough to characterize early stage tumors that areconsidered in, e.g., training a classifier for an image guided biopsy application. Since MREdepends on fairly complex modulus reconstruction, we expect the registration results toimprove with the developments in our reconstruction methods.As other registration techniques, the proposed algorithm performs better on the mid-583.6. Discussion and Conclusiongland region of the prostate. This finding is due to the ambiguity in contouring the prostatearound the base and apex on the ex vivo and in vivo images, that impacts the fidelity ofboth the model and evaluation in these regions. In future work, we may utilize the intensitysimilarity between in vivo and ex vivo T2w MRI to derive forces from intensity-baseddistance measures, e.g., MI. This may eliminate the required segmentation of the ex vivovolume, and improve speed and accuracy.A disadvantage of the proposed approach is that the ex vivo MRI scan is time consuming,expensive and not accessible. However, as discussed above, the proposed method shouldbe employed on a limited number of cases for the designated application. Moreover, asdiscussed in Section 3.1, methods that propose direct registration between in vivo MRI tohistopathology require either dense sectioning of histological slices [111], or manual selectionof corresponding slices [24]. Another direct approach, proposed in [148] and [162], generatesa patient-specific mold, based on the (non-deformed) in vivo prostate, to preserve sectioningorientation. Thus, current direct registration methods can be just as time consuming as theex vivo scan, while being more user-dependent and error-prone.We note, however, that since the ex vivo scan in our approach is used only to acquirea model of the specimen, a direct registration to histopathology is still possible using ourmethod by using a 3-D model reconstruction based on the histological slices, provided thattheir sectioning is dense enough. In addition, other ex vivo imaging modalities, which aremore accessible than MRI, can be employed. In fact, an US generated ex vivo model maybe combined with methods such as [158] to register the ex vivo scan to histology followingour method. Thus, our framework is general and can be fitted into different workflows.The current implementation does not allow real-time performance, which is not a re-quirement in the case of cancer characterization. The execution times are in the order ofminutes, so multiple registrations can be processed on the background in a typical clinicalscenario to provide a reasonable throughput. Another current limitation is the assump-tion of small deformations. Although the residual deformations are small after a similaritytransformation, we may still consider nonlinear elasticity for future work in order to contendwith large deformations of the prostate, e.g., in case of employing a transrectal coil duringimage acquisition.59Chapter 4Registration of Ultrasound andMagnetic Resonance Imaging forImage Guided Interventions4.1 IntroductionIntraoperative TRUS imaging has been shown to be valuable during CaP interventions suchas RALP [80], but is of poor quality and not reliable for detecting cancer. The ability ofmulti-parametric MRI to generate the most accurate characterization of CaP [75], has led tothe development of methods for MRI-guided treatment, mainly biopsy and brachytherapy[159]. Indeed, a systematic literature review in [97] suggests that MRI-guided biopsy detectsclinically significant CaP using fewer samples compared to a standard biopsy that employsTRUS alone.While direct MRI-guided methods have been reported, e.g., in [153], for prostate biopsyand brachytherapy, intraoperative MRI is still cumbersome, time consuming, resource costlyand not widely used. Cognitive fusion, in which the clinician estimates the lesion’s locationon the intraoperative TRUS based on a preoperative MRI, varies greatly with expertise.A more feasible approach to allow integration of MRI data in the operating room involvesregistration of the preoperative MRI to the intraoperative TRUS, and visualization of thecorresponding images to assist the clinician during treatment. Previously, such an approachwas successfully demonstrated, e.g., in [51, 86, 91, 120, 149, 164], for prostate biopsy, andin [28, 126] for prostate brachytherapy.Due to the difference in intensity information, contrast, and anatomical visualizationpresented by the two modalities, techniques for registration of MRI and TRUS typicallyinvolve preprocessing of the images. In [101, 157] both modalities are segmented and thesurfaces are registered using deformable registration for targeted prostate biopsy. In [90],corresponding 2-D MRI and TRUS images are manually selected, segmented, and thenregistered by minimizing the Bhattacharyya distance between their shape representations,604.1. Introductionwith thin-plate splines as regularization. In [59, 87], a patient-specific biomechanical modelis generated from segmented T2w images, and a model-to-image deformable registrationis performed based on features extracted from 3-D TRUS images. In [156], features areextracted by employing [55], and the distance between modality independent neighbourhooddescriptors (MIND) is minimized with a smoothness regularization.In this chapter, we propose to introduce MRI information during RALP. To this end,in the first part of this chapter, we propose a segmentation-based registration method thatincorporates a semi-automatic prostate segmentation algorithm, developed in our lab, tosegment a TRUS volume, acquired from an intraoperative sweep of the probe. A preopera-tive segmented MRI is then registered to the segmented intraoperative TRUS. The methodcan be incorporated with a tracking algorithm, developed in our lab as well, that allowsthe TRUS probe, that is mounted on a motorized system [143], to track the robot’s tooltip.Corresponding TRUS and MRI planes that follow the tooltip may then be visualized inreal-time on the surgeon console.In the second part of this chapter, we propose to employ elastography for intensity-basedregistration of MRI and US. Elastograms have been shown to be promising in improving thevisibility of the prostate gland [143]. The produced images typically have superior prostate-to-background contrast compared to T2-weighted images, especially at the base and apex.This is due to the fact that prostate tissue is generally stiffer than the surrounding tissue.The challenges in employing elastography for image processing of the prostate are theirlow resolution (typically of MRE volumes) and poor signal-to-noise ratio. Also, the bound-ary of the prostate may be blurred and assimilated with the background at some regions.Nevertheless, previously in [83], the authors showed that fusion of region information fromTRUS-VE images and edge information from B-mode images improved the performance ofan active shape model for prostate segmentation. In [102], MRE is employed for prostatesegmentation by using the elastogram and corresponding magnitude image for driving aregion- and edge-based active contours, respectively. In Chapter 3, we described the uti-lization of elastography as a regularizer for prostate registration.In our data collection, both MRE and USE images are acquired and should ideally rep-resent the same tissue elasticity in the scanned anatomy. We investigate the use of the elas-ticity values in a similarity metric to drive the registration. Since the elastography volumesare inherently registered with their corresponding T2w and B-mode volumes, MRE-USEregistration may serve as an inter-modality to facilitate such multimodality registration.614.2. Problem Formulation(a) IUS,SAGx [mm]y[mm]−50 −40 −30 −20 −10 0 10 20 30 40 5010203040506070(b) IUS,TRV⇔x [mm]y[mm]0 20 40 60 80 100 120020406080100120(c) IMR,TRVFigure 4.1: Problem illustration. (a) Intraoperative sagittal B-mode images acquired by asweep of the TRUS probe span a sector in cylindrical coordinates. (b) A transverse B-modeimage after conversion to cartesian spatial coordinates. (c) A preoperative transverse T2wMR slice.4.2 Problem FormulationSuppose we are given IUS,SAG(ρ, θ, z), the intraoperative B-mode volume spanned by sagittalimages that are acquired by a sweep of a TRUS probe. The cylindrical coordinate vectoris denoted by (ρ, θ, z)T , with ρ ∈ [ρ0, ρ1] as the radial distance from the center of theprobe (ρ0 is the probe’s radius), θ ∈ [θ0, θ1] as the probe’s angle, and z ∈ [0, z1] as theaxial coordinate on the sensor. We can then use the following conversion to obtain thecorresponding transverse TRUS volumeIUS,TRV(r¯) = IUS,SAG(x2 + y2, arctan 2(y, x)−pi2 , z), (4.1)where r¯ = (x, y, z)T ∈ R3 denotes the 3-D cartesian spatial coordinate vector, with its originat the probe’s center. Figure 4.1 illustrates the TRUS B-mode volumes.Suppose we are also given IMR,TRV(r¯), the preoperative T2w MR volume in cartesiancoordinates. We are interested in a displacement field u¯ = (ux, uy, uz)T , such that thedeformed transverse MR volume,I˜MR,TRV(r¯; u¯) = IMR,TRV(r¯ − u¯(r¯)), (4.2)is similar, in some sense, to the transverse TRUS, volume IUS,TRV. Similar to Chapter 3,the registration problem is formulated as the minimization of the energy,J [u¯; IUS,TRV, IMR,TRV] = F [IUS,TRV, I˜MR,TRV(·; u¯)] + S[u¯] , (4.3)where F and S are the similarity metric and regularization functionals, respectively.624.3. Segmentation-Based RegistrationWe can then use the following conversion to obtain the corresponding registered sagittalMR volumeI˜MR,SAG(ρ, θ, z) = I˜MR,TRV(−ρ sin(θ), ρ cos(θ), z), (4.4)in order to visualize the TRUS and MR volumes together.4.3 Segmentation-Based RegistrationTypically, as part of a volume study, the prostate gland on each transverse slice in thepreoperative T2w MR volume is segmented manually by an expert before surgery. Inaddition, we employ a real-time semi-automatic algorithm for 3-D segmentation of theprostate on the intraoperative TRUS B-mode. The algorithm, described in details in [142],has been routinely employed during brachytherapy, and found to be a fast, consistent andaccurate tool for the delineation of the prostate gland in TRUS images. Below, we providea short summary of the method.4.3.1 Semi-Automatic Prostate Segmentation in TRUSThe user interaction part of the algorithm is required for initialization and comprising twosteps. In the first step, the user selects two points that correspond to the base and apex ofthe prostate on a sagittal image. This selection of base and apex set the center as mid-glandplane, on which the prostate is the largest and most visible.The second step is to select seven key geometrical points on a transverse image fromthe mid-gland. These points consist of the center of the TRUS probe, four points onthe prostate’s boundary – the lowest posterior lateral, extreme right, mid-posterior, mid-anterior, and two points that intersect the boundary with perpendicular lines that connectsthe extreme right point with the lowest posterior lateral and mid-anterior points. Thesesteps of the algorithm and the resulting segmentation contour on the mid-gland image aredepicted in Figure 4.2.After user initialization, the 3-D segmentation is performed automatically. The algo-rithm un-warps the mid-gland image to reduce the deformation caused by the probe, thenuses a tapered ellipse, fitted to the selected points, to guide an edge detector algorithm, [1]that evolves the contour edge on each 2-D transverse image using an interacting multiplemodel estimator. The images are then un-tapered, and the contour estimations are cor-rected to allow for a tapered ellipsoid fitting for smoothness and continuity of the contours.Finally, the contours are tapered and warped to match the original images.634.3. Segmentation-Based Registration(a) Base-apex selection (b) Mid-gland points selectionx(c) Segmented sliceFigure 4.2: Semi-automatic segmentation. (a) The user selects the base and the apex ona sagittal image. (b) The user selects the center of the probe, and six boundary points onthe transverse mid-gland image. (c) The resulting segmentation contour.4.3.2 Registration ProcessBased on the segmented surfaces of the prostate in the TRUS and MR volumes, we constructthe binary volumes HUS,TRV and HMR,TRV, respectively, each with values of “1” on theregion inside the prostate, and “0” outside. We may then register the two binary volumesand apply the resulting displacement map back to the MR volume IMR,TRV.First, the two binary volumes are rigidly aligned (and scaled) using the principal axestransformation [4]. This is a fast, one-step registration that was found to provide a good ini-tial alignment for the followed deformable registration. To this end, the binary volumes aretreated as pdfs, and the corresponding “centers of mass”, m¯US, m¯MR ∈ R3, and covariancesMUS,MMR ∈ R3×3, are calculated asm¯i =∫r¯Hi,TRV dr¯∫Hi,TRV dr¯Mi =∫(r¯ − m¯i)(r¯ − m¯i)THi,TRV dr¯∫Hi,TRV dr¯= UiΣ2iUTi , i = US,MR (4.5)where UUS,ΣUS and UMR,ΣMR ∈ R3×3 are the eigendecomposition of the correspondingcovariance matrices such that Ui is comprising the eigenvectors in its columns, and Σicontaining the eigenvalues in its diagonal.The aligned MR binary volume can be then computed using a linear coordinate trans-formation,HrigidMR,TRV(r¯) = HMR,TRV(Ar¯ + b), (4.6)with b = m¯MR −Am¯US, and A = UMRΣMRΣ−1USUTUS.644.3. Segmentation-Based Registration(a) Initial (b) Rigidly aligned (c) RegisteredFigure 4.3: Segmentation-based registration. The surface of the segmented prostate on theTRUS (red) and MR (blue) volumes at different stages of the registration process.Next, we apply (4.3) to the binary volumes, and minimize J [u¯;HUS,TRV, HrigidMR,TRV]. WetakeF [u¯] =∫ (HrigidMR,TRV(r¯ − u¯(r¯))−HUS,TRV(r¯))2dr¯ ,S[u¯] =∫µ43∑j=13∑k=1(∂u¯j∂r¯k+ ∂u¯k∂r¯j)2 + λ2 (div u¯)2 dr¯ , (4.7)i.e., the SSD metric, and an homogenous elastic regularizer (the Lame´ parameters λ andµ are constant and set arbitrarily with an almost incompressible Poisson ratio ν = 0.499).Similar to Chapter 3, we minimize (4.7) in a variational approach by computation of thecorresponding Euler-Lagrange equation and its discretization (see also [92]). The proposedregistration process is illustrated in Figure 4.34.3.3 ResultsThe proposed method was tested on six data sets of patients who were scheduled for radicalprostatectomy and underwent preoperative MR and intraoperative US. The MR volumeswere acquired by a 3-Tesla system (Achieva 3.0T, Philips, The Netherlands) using a standard6-channel cardiac coil with acceleration factor (SENSE) 2. A standard axial T2w FSEsequence was acquired (TE/TR = 86/2500 ms, FOV = 320× 320× 70 mm3, with 0.5 mmin-plane resolution and 2 mm slice thickness).The B-mode images were collected using an US machine (Sonix RP, Ultrasonix Medical,Richmond, BC, Canada) with a biplane transrectal probe (BPL9-5/55, 5–9 MHz, withradius ρ0 = 10 mm) mounted on a motorized cradle [143]. A continuous sweep of sagittal654.3. Segmentation-Based RegistrationRigid registration FEM method ([157]) Proposed methodVO Mean error VO Mean error VO Mean error(%) (mm) (%) (mm) (%) (mm)Case 1 85.8 1.39 95.7 1.19 97.3 1.57Case 2 70.1 3.03 95.0 1.52 98.0 2.06Case 3 75.1 1.67 94.6 0.87 97.6 1.11Case 4 76.6 2.64 95.6 1.06 97.9 1.77Case 5 72.9 2.73 94.3 1.41 98.0 1.10Case 6 80.6 1.24 95.6 0.93 97.1 1.07Mean 76.8± 5.1 2.11± 0.77 95.1± 0.5 1.16± 0.26 97.7± 0.3 1.44± 0.42Table 4.1: Segmentation-based registration results.images was acquired (radial distance ρ1 = 79 mm, angular range θ0,1 = ±pi4 , with 0.5 mmradial and axial resolution, and 0.2◦ angular resolution).Both of the MR and TRUS transverse volumes were segmented using Stradwin (Cam-bridge University, UK), [160], and the surfaces of the prostate were generated with mediumresolution and smoothing strength settings. The urethra on both modalities was segmentedmanually for evaluation as well. In order to maintain consistent processing times, we ranthe deformable registration algorithm a fixed number of 30 iterations that takes about 184sec, and was found to converge sufficiently. Including the times, as reported in [142], formanual selection of initialization points (32 ± 14 sec) and segmentation (14 ± 1 sec), thetotal runtime of the entire segmentation-based registration process is about 4 min.We evaluated the registration performance using the volume overlap, in the sense ofDice’s coefficient (3.16), between the MR and the TRUS volumes after rigid and deformableregistration. The mean distance between splines on both modalities that were fitted throughthe center points of the urethra was also measured for the rigid and deformable registeredvolumes. In addition, our proposed variational approach was compared using the sameevaluation metrics to a surface-based FEM registration approach proposed in [157]. Thequantitative results are summarized in Table 4.1.664.3.Segmentation-BasedRegistrationz [mm]ρ[mm]−30 −20 −10 0 10 20 3010203040506070z [mm]ρ[mm]−30 −20 −10 0 10 20 3010203040506070(a) Current console displayx [mm]y[mm]−50 −40 −30 −20 −10 0 10 20 30 40 5010203040506070(b) B-modex [mm]y[mm]−50 −40 −30 −20 −10 0 10 20 30 40 5010203040506070(c) T2wFigure 4.4: Proposed augmented visualization. (a) The current surgeon’s console display augmented with intraoperativeTRUS using TilePro, as proposed in [94]. (b) Intraoperative sagittal (top) and transverse (bottom) TRUS. (c) The proposedadditional data of the registered sagittal (top) and transverse (bottom) MR. The line (red) represents the real-time probe’ssagittal angle and axial position to facilitate navigation and orientation.674.4. Elasticity-Based Registration4.3.4 Application in MR-Guided RALPAs discussed in Section 4.1 above, the purpose of the proposed registration is to allowintraoperative MR-guidance during CaP interventions. Specifically, we intend to employ theregistration during RALP and augment the display of the surgeon’s console with relevantpreoperative MR data.In [94], the authors introduced a system of an intraoperative motorized TRUS that iscalibrated to the da Vinci system and tracks the robot’s tooltip. This allows the surgeon tocontrol the TRUS imaging plane automatically by moving the surgical instrument, and towatch the corresponding image of the prostate and periprostatic anatomy from the surgeon’sconsole in real-time during RALP. The system has been tested clinically and was founduseful in certain stages of the surgery, e.g., in defining the prostate-bladder interface duringposterior and anterior bladder neck dissection, and identifying the prostate boundary at theapex.The tracking method involves a rigid registration between the da Vinci surgical ma-nipulators coordinates and the motor of the TRUS probe manipulator. This registration,described in [2], is performed as initialization and requires localization of the surgical tooltipsin the 3-D TRUS. The da Vinci application programming interface (API) is then used toprovide the locations of the tool tips in the da Vinci kinematic frame, and to update theselocations in real-time.The registered TRUS images are streamed to the surgeon’s console and displayed belowthe standard camera view of the surgical site using the TilePro feature of the robot. Wepropose to employ the segmentation-based MR-TRUS registration above, and to display, inaddition to the sagittal B-mode, the corresponding registered T2w image, and the transversecounterparts of both the B-mode and T2w. Figure 4.4 illustrates the current, and ourproposed augmented visualization.We note, however, that while the sagittal TRUS image can be streamed and displayedin real-time, the rest of the images are frozen and only valid with respect to the time of theinitial sweep of the probe, when the data was acquired. In case time and clinical conditionsallow, additional sweeps followed by the proposed registration process may be acquired inorder to update the images.4.4 Elasticity-Based RegistrationAs discussed in Chapter 3, utilizing elasticity information can facilitate registration and im-prove its accuracy. Indeed, one approach to utilize the acquired elastography information684.4. Elasticity-Based Registrationfor TRUS-MR registration is to employ the proposed model-based algorithm from Chap-ter 3, with the segmented intraoperative TRUS volume as the reference model, and MREmagnitude and elastogram volumes as the template and its regularizer. Alternatively, onemay employ this algorithm with the segmented MRI as the reference model, and the TRUSB-mode and USE volumes as the template. However, this intensity-driven and regularizedapproach may require more iterations with smaller step sizes, which will yield longer runningtimes that may not be suitable for intraoperative deployment of the algorithm.Another approach to employ elastography is to use [143] and [102] to segment the USEand MRE volumes, and employ the proposed segmentation-based method above for TRUS-MR registration. However, this approach too requires preprocessing steps that may generateerrors. In this section, we are taking a different approach and investigate the utilization ofthe “pure” elastography data, of both modalities, as an inter-modality to allow intensity-based registration of the corresponding MR and US volumes.4.4.1 Registration ProcessSuppose that we have EUS,SAG(ρ, θ, z), the sagittal intraoperative USE volume. Using theconversion in (4.1), we obtain the transverse USE volume, EUS,TRV(r¯). In addition, we aregiven EMR,TRV, the transverse MRE volume, and assume that the center of the rectum(that corresponds to the TRUS probe in USE) is found, either manually or automatically.The volume can be then converted to the sagittal EMR,SAG volume using (4.4) with respectto that center.The USE method is based on estimation of displacements along the radial direction onthe sagittal plane, on which the data is in its native resolution. As shown in Figure 4.6, theelastograms are blurred and noisy. Thus, we limit the transformation model to translationand rotation about the z-axis, and formulate the problem in sagittal coordinates asminδθ,δzF [EUS,SAG, EMR,SAG(ρ, θ + δθ, z + δz)], (4.8)where F is the distance metric between the two volumes.Ideally, the elastograms of both MRE and USE represent the same tissue elasticity, andtherefore an SSD or NCC metric would be a proper choice. However, as can be seen by atypical joint histogram of two registered elastograms in Figure 4.6(c), there is a nonlinearand non-injective relationship between them. Therefore, we chose to minimize the inverseof the (normalized) MI as the metric F .694.4. Elasticity-Based Registrationz−translation[mm]z−rotation[◦ ]−10 −8 −6 −4 −2 0 2 4 6 8 10−10−8−6−4−20246810(a) MI metric(b) USE (c) Registered MREFigure 4.5: Elasticity-based registration of a phantom. (a) The MI similarity map withrespect to axial translation and rotation. The maximum value is at (0,0), where the vol-umes are registered. (b)-(c) Cross-sections of the registered volumes. Intensity scale of theelastograms is 0-50 kPa.Phantom ExperimentWe first study registration of a prostate elastography phantom (Model 066, CIRS, VA,USA). The USE system, described in [79], includes an US machine (BK Medical, Herlev,Denmark) with a biplane transrectal transducer (Endocavity 8848, 4–12 MHz) mountedon a motorized cradle, similar to Section 4.3.3. The algorithm in [10] was implementedto provide real-time absolute elastography using a transperineal voice coil as the shaker(with 75 Hz excitation). We also used the transperineal MRE system with similar scanningparameters and processing to the description in Section 3.5.1.Here, we study the use of the MI metric as a similarity measure. The MRE and USEvolumes were registered manually, based on the prostate and inner inclusions. Then, axialtranslations and rotations were applied, with the MRE volume as the template, and theMI value was recorded for each pose. In Figure 4.5(a), the resulting MI similarity mapis depicted with respect to the axial translation and rotation parameters. As desired, the704.4. Elasticity-Based Registrationz[mm]ρ[mm]−30 −20 −10 0 10 20 3010203040506070(a) USEz[mm]ρ[mm]−30 −20 −10 0 10 20 3010203040506070(b) Registered MREEUSE MR0 5 10 15 20 25 30 35 40 45 5001020304050(c) Joint histogramFigure 4.6: Elasticity-based registration of patient data. (a) A sagittal slice of the USEreference. (b) A sagittal slice of the registered MRE template. (c) The joint histogram ofthe registered volumes.global maximum value is at (0,0), where the volumes are registered. However, we note thatthere are many local maxima as well, and that the function changes slowly with respect torotations. Thus, a gradient descent based algorithm fails to optimize the metric, and weemploy search-based methods, e.g., the Nelder-Mead simplex method.Patient DataThe registration was performed offline on four patient’s data sets that comprise both pre-operative MRE and intraoperative USE. First, the volumes were smoothed using 3× 3× 3median filtering. Due to off-axis rotation of the two volumes, the center of the rectum in theMRE volume was selected on both the base and apex of the transverse magnitude image,and the sagittal volume was rotated accordingly. Then, the Nelder-Mead simplex methodwas used to minimize the inverse of the normalized MI between the volumes, with respectto axial rotation and translation.Illustration of the resulting registration of one data set is presented in Figure 4.6. On theselected sagittal slice, we see that there are matching landmarks on both volumes. However,some landmarks significantly differ on the two modalities. To evaluate the registration, theresulting transformation was applied to the binary volume of the segmented MR, and wasthen compared to the rigid registered volume that was obtained using the principal axestransformation, as described above in Section 4.3.2. In two of the cases, depending on theinitial pose, we noticed a convergence to local minima. Still, the resulting mean volumeoverlap of 81.2%± 6.4% is relatively high.714.5. Discussion and Conclusion4.5 Discussion and ConclusionIn this chapter we presented methods for registration of intraoperative TRUS and preopera-tive MRI. With MRI correctly registered, multi-parametric imaging and cancer probabilitymaps can be displayed on the surgeon’s console to augment visualization of the prostate,periprostatic anatomy, and CaP for facilitating orientation, planning and navigation duringRALP.We proposed a segmentation-based approach that incorporates an existing semi-automaticsegmentation of the prostate in the B-mode volume, and then deforms the segmented T2wvolume to match the two surfaces. The registration performance was found to be plausible,with comparable values to other TRUS-MR registration methods. Specifically, in terms ofvolume overlap, the proposed method was found to be slightly superior to a similar methodthat uses FEM. The registration errors for the proposed method, based on the urethra,were slightly inferior to the FEM method. However, since the urethra can be significantlydeformed within the prostate, especially in the presence of a Foley catheter, as is the caseduring intraoperative TRUS imaging, it is not necessarily a good anatomical landmark.In order to evaluate the registration alone, i.e., differentiate its errors from the segmen-tation errors, we used the manual segmentations of the TRUS in our experiments ratherthan the resulting semi-automatic segmentations. Thus, to evaluate the expected intra-operative errors, one should aggregate the reported segmentation error reported in [142].Similarly, for using the tooltip tracking system, the errors in calibration and during track-ing should be considered. Further evaluation of the algorithm may employ histopathologyimages, registered to both MR and TRUS by using, e.g., the algorithm in Chapter 2, as aninter-modality for comparing the MR volumes before and after registration to TRUS.In addition, compared with the FEM method, the proposed method does not requiremeshing, and its runtime is predictable (constant, for a fixed number of iterations) andsignificantly shorter (3 min as opposed to 7–9 min). Including the setting up the equipment,the initial sweep of the probe, segmentation, registration and calibration to the robot, weexpect the entire process to consume 12-15 minutes of operating room time prior to surgery.From our past experience during RALP in [94, 98], this is a feasible time frame at thebeginning of the procedure, with which clinicians are comfortable.We note, however, that even with the most accurate imaging and registration, relyingsolely on preoperative information to define a surgical plane is not recommended, due tohigh risk of positive margins CaP. Therefore, the proposed MR-TRUS registration is moresuitable for surgical planning, at the beginning of the procedure. A clinical study thatemploys the proposed system is to be conducted in order to determine its role in achieving724.5. Discussion and Conclusionsurgery goals and improving outcomes. The proposed registration can also be employedoffline for learning 2-D TRUS features that correspond to findings in MR, which may thenassist clinicians in real-time guidance using TRUS.In the second part of this chapter, we investigated intensity-based registration of MREand USE. Such registration can be used for aligning the corresponding T2w and B-modevolumes without preprocessing such as segmentation. We found that the discrepancy be-tween the elasticity values in the two elastograms does not allow employment of unimodalsimilarity measures. Nevertheless, by employing a proper metric, e.g., MI, the intensitiesof the elastograms may allow a rigid alignment.The difference in elasticity values of the two volumes may stem from the noisy character-istics of the elastograms, different acquisition methodologies, different spatial sampling thatmay introduce aliasing of different higher frequency components, difference in displacementestimation (full 3-D in MRE, and radial in USE, errors caused by the pressure of the TRUSprobe, and a slightly different excitation frequency (70 Hz in MRE and 75 Hz in USE).Other future work in elasticity-based registration may include extraction of intensityinvariant features that are stable and unique to the elastograms, thus reducing the problemto point-based registration. These features can be stiff inclusions that are clearly visibleon both modalities, or edges that can be used as surface landmarks. Another approach isto utilize the MRI and MRE volumes in order to simulate images that are closer to USor USE, similar to the approach in [169], where US images were simulated from computedtomography volume to facilitate multi-modality registration. Also, improving data acqui-sition, e.g., by increasing MRE resolution, may allow for a better registration. However,such improvements may not be feasible in practice, either preoperatively or intraoperatively,due to a longer scanning time. In the following chapters we take different approaches forimproving MRE, without prolonging acquisition time.73Chapter 5Motion Compensation andSuper-Resolution in MagneticResonance Elastography5.1 IntroductionIn conventional MRE, each spatial component of the 3-D displacements is encoded sepa-rately by applying motion encoding gradients (MEGs) along the corresponding direction.The wave is sampled at different stages during its cycle by varying the phase shifts be-tween the mechanical excitation and the MEGs in each acquisition. Although the Nyquistrate requires only two samples per direction to estimate the amplitude and phase of thedisplacements at each voxel [104, 166], typically more samples are acquired to improvesignal-to-noise ratio (SNR). A temporal discrete Fourier transform (DFT) is then per-formed to extract the displacements’ amplitude and phase from the fundamental frequencycomponent of the samples.Due to the oversampling of the wave, an MRE scan with a typical field of view is alengthy process that may cause inconvenience to the patient, and increases the possibilityof both voluntary and involuntary movements. As in MRI, patient’s movements in MREintroduce motion artifacts such as ghosting in the reconstructed images. Furthermore, inMRE, such motion corrupts displacement estimation, which in turn degenerates elastogramreconstruction. Thus, there is a tradeoff in acquiring more samples to improve SNR, orfewer samples to reduce patient’s movements. For a clinical usage of MRE, a compromiseis to allow reasonable scanning times by acquiring four or eight samples of lower spatialresolution images.Previously, several approaches were considered to suppress patient’s movement duringMRI scans. These range from physical constraints such as a belt wrapped around thepatient in free-breathing abdominal MRI [62], through specialized filling of k-space [121],to image processing techniques [170].745.2. MethodsIn this work, we utilize the redundancy of images acquired in MRE, and apply imageregistration to compensate for patient’s motion during the scan. We also employ a super-resolution technique to enhance the resolution of the acquired images to allow for a betterestimation of tissue displacements.In the MRE acquisition process, the phase component of the complex imaging signal, alsoknown as the phase image or volume, encodes the displacement, while its magnitude com-ponent is still available as a low-resolution T2w image, also known as the magnitude imageor volume, that provides similar intensity information and features as high-resolution one[102]. Since the corresponding phase images are derived from the same signal as the magni-tude images, each such pair of phase-magnitude images is inherently registered. Ideally, allthe magnitude images, acquired in different directions and wave samples, are identical. Inpractice, however, the magnitude images vary due to noise and patient’s movements. Thus,motion among different acquisitions of phase images can be compensated by registering thecorresponding magnitude images.In case of multiple images of the same scene or object, super-resolution methods canbe used to increase the native size of images [113]. Rather than basic interpolation meth-ods, which introduce artifacts such as blur and contrast loss, super-resolution techniquesefficiently improve resolution by incorporating information from the different images. Super-resolution methods have been applied to medical imaging to enhance the acquired data andincrease diagnosis [49].In MRE, the DFT processing of the oversampled phase images coincides with a linearleast-squares estimation of the displacements [104]. Here, we consider pairs of phase imagesthat represent the same displacement’s phase and amplitude (up to a sign flip and noise),and should complement each other in other in terms of intensity information after motioncompensation. In addition, due to the mechanical nature of MRE sequence, althoughactual data acquisitions are at rest, small displacements of the scanned object occur on asub-voxel scale, and the images may contain different information. Thus, we propose toemploy a super-resolution method to utilize this complementary property and increase theresolution of the acquired data.5.2 Methods5.2.1 Motion Compensation Using RegistrationLet I(i,k), Φ(i,k) : R3 → R denote the kth acquisition of the magnitude and phase volumes,respectively, in direction i ∈ {x, y, z} on the 3-D domain D ⊆ R3. We consider a conven-755.2. Methodstional MRE sequence that comprises K acquisitions in each direction, with phase shifts ofthe MEGs that are evenly distributed over [0, 2pi).As discussed above, we compensate for patient’s motion by registering the magnitudevolumes. Without loss of generality, we consider I(x,1) as the reference volume, and registereach template volume I(i,k) by minimizing the following energy functionalJ [u¯(i,k);I(i,k), I(x,1)] =∫D(I(i,k)(r¯ − u¯(i,k))− I(x,1)(r¯))2dr¯ +∫Dµ43∑j=13∑l=1(∂u¯(i,k)j∂r¯l+ ∂u¯(i,k)l∂r¯j)2 + λ2 (div u¯)2 dr¯ ,∀ i ∈ {x, y, z}, k ∈ {1, . . . ,K}, (5.1)where r¯ = (x, y, z)T is the spatial coordinate vector, and λ and µ are the Lame´ parametersthat reflect the elastic properties (µ is also known as the shear modulus). The minimizationof the functional is with respect to the deformation field, u¯(i,k) = (u(i,k)x , u(i,k)y , u(i,k)z )T : D →R3, which maps I˜(i,k) = I(i,k)(r¯ − u¯(i,k)) to I(x,1).The first component of the functional (5.1) is a similarity measure between the templateand reference volumes, while the second component is a regularizer of the deformation map.Since the template and reference are both magnitude volumes of the same intensity range,we employ SSD as an intramodality similarity measure. Due to the elastic nature of thedeformation, the regularizer was chosen as the linear elastic potential.The minimization of (5.1) is performed in a variational approach. The algorithm isimplemented on a fixed grid using finite differences, similar to the description provided in[92]. The registration of one template volume to the reference is demonstrated in Figure 5.1.5.2.2 Super-Resolution in MREAs discussed in Section 5.1, the phase and magnitude images are inherently registered. Thus,we apply the resulting deformation map of each magnitude image to the correspondingphase image, Φ˜(i,k) = Φ(i,k)(r¯ − u¯(i,k)), to compensate for motion in the phase images.Considering a conventional MRE sequence of K phase shifts evenly distributed over [0, 2pi)in each direction, the phase images Φ˜(i,k) and −Φ˜(i,k+K2 ) should be similar and differ on asub-voxel scale, since they sample a certain phase and anti-phase of the displacement wave.As in most super-resolution approaches, we attempt to recover a high-resolution imageof double the size in each axis of the native resolution volume. We employ an iterativeback-projection super-resolution method, first introduced in [61], with a steepest-descent765.2. Methods[mm][mm]−60 −40 −20 0 20 40 60−60−40−200204060(a) I(x,1) (the reference)[mm][mm]−60 −40 −20 0 20 40 60−60−40−200204060(b) I(z,8) (a template)[mm][mm]−60 −40 −20 0 20 40 60−60−40−200204060(c) I˜(z,8) (registered temp.)Figure 5.1: Motion compensation using registration. Cross-section of the MRE magnitudevolumes from the mid-gland region of a patient’s prostate scan. The blue grid representsthe sampled deformation field.choice of the back-projection kernel [35].The iterative algorithm, adapted to the context of phase images in MRE, starts with aninitial guess of the super-resolution image, Φˆ(i,k),0, which in our case is taken as the meanphase image 12(Φ˜(i,k) − Φ˜(i,k+K2 )), interpolated on a high-resolution volume. Then, in eachiteration j = 1, 2, . . . , the high-resolution volume is updated byΦˆ(i,k),j+1 = Φˆ(i,k),j +ABP(Φ˜(i,k) − Φ˜(i,k+K2 ) − 2ATBP Φˆ(i,k),j),∀ i ∈ {x, y, z}, k ∈ {1, . . . , K2 }, (5.2)where the back-projection “tall” matrix ABP comprises translation, blurring and downsam-pling of the high-resolution image. Note that the 3-D volumes are rearranged as columnvectors in (5.2), and similarly, the resulting long column vector is rearranged back into thesuper-resolution volume.In each iteration, the algorithm minimizes the SSD between the two native resolutionphase images and the translated, blurred and downsampled estimation of the current super-resolution volume. Note that since the images are motion compensated, zero translationsare chosen. Also, without prior information on the scanner, a 3×3×3 kernel of a Gaussianblurring is used. A temporal DFT is then performed on the super-resolution phase imagesΦˆ(i,k), k = 1, . . . , K2 to extract the super-resolution displacements’ amplitude and phasefrom the fundamental frequency component of the samples, followed by reconstruction ofthe super-resolution elastogram.775.3. Results5.3 Results5.3.1 Phantom ExperimentWe have acquired an MRE scan of an elasticity phantom (Model 049, CIRS, VA, USA)using a custom-made shielded electromagnetic transducer synchronized to the 3-Tesla scan-ner (Achieva 3.0T, Philips, The Netherlands). We employed a spin-echo (SE)-echo-planarimaging (EPI) pulse sequence with a TE/TR of 30/1600 ms. A 200 Hz vibration was ap-plied to the rectangular shaped phantom. As seen in Figure 5.2(e), the phantom containsspherical inclusions of different radii and stiffness values relative to the background mate-rial that mimic different tissue types. The entire phantom was scanned in a field of viewof 128 × 128 × 20 grid, with an isotropic 2 × 2 × 2mm3 resolution. We have oversampledthe scan of each acquisition with eight phase shifts of the MEGs ranging evenly over [0, 2pi)relative to the mechanical excitation. The scanning time was 364s per acquisition.A 80 × 60 × 20 region of interest on the acquired axial phase accumulation imageswas selected, compensated and phase unwrapped. Since the phantom is static, the motioncompensation algorithm was not used. We applied the super-resolution algorithm describedin Section 5.2.2 to double the resolution of the phase images. In addition, for comparison,a high-resolution phase images was computed using bilinear interpolation of the nativeresolution phase image. The native and high-resolution elastograms were reconstructedfrom the corresponding displacement maps using a LFE technique [85].The elastograms are illustrated in Figure 5.2 for qualitative comparison. A correspond-ing mid-gland T2w image with an in-plane resolution of 0.732×0.732mm2 is also presentedto reveal the geometry of the phantom’s inclusions. Further quantitative evaluation, pre-sented in Figure 5.3, was performed by comparing the elastograms to the factory reportedelasticity values of the inclusions. The root-mean-square errors between the factory values,and the native resolution and the super-resolution elasticity values were 6.57 ± 5.55 kPaand 4.73± 4.34 kPa, respectively.785.3.Resultsx [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(a) Native resolutionx [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(b) Bilinear interpolationx [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(c) Super-resolutionx [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(d) T2-weightedx [mm]y[mm] Type 4 Type 3 Type 2 Type 1Type 0−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050 102030405060(e) True elastogram (in kPa)Figure 5.2: Super-resolution of a phantom MRE. (a)-(c) Corresponding cross-sections of the reconstructed elastograms.Intensity scale is in kPa and same as in (e). (d) The T2w image. (e) The factory reported elasticity values of tissue types.795.3. Results5152535455565Tissue typeElasticity[kPa] 0 1 2 3 4Native resolutionSuper−resolutionTrueFigure 5.3: Phantom elasticity errors. The errors between the native and super-resolutionelastograms in comparison with the true elasticity on each tissue type.5.3.2 Patients’ DataThe patients’ data was acquired as part of a project for evaluating CaP detection usingMRE. After approval of the institutional ethics board and signed consents, twelve patients(mean age 65, range 59–76 years old) who were scheduled for radical prostatectomy, under-went a preoperative scan in a 3-Tesla MRI scanner (Achieva 3.0T, Philips, The Netherlands)with a standard 6-channel cardiac coil with acceleration factor (SENSE) of 2.A standard axial T2w FSE sequence was acquired (TE/TR = 86/2500 ms, FOV =320 × 320 × 70 mm3, with 0.5 mm in-plane resolution and 2 mm slice thickness). Next, aprostate MRE scan was acquired using a transperineal shielded electromagnetic transducer[135], and a gradient echo sequence named eXpresso [136]. Eight phases were encoded witha mechanical excitation of 70 Hz applied to the perineum of the patient. Wave images wereacquired on a 128× 128× 24 matrix with 2 mm isotropic voxel size. The MRE scan lasted8 min for a 3-D displacement field acquisition. Each slice was cropped to accommodate64× 64 pixels around the prostate region for computational efficiency.We applied the motion compensation algorithm proposed in Section 5.2.1 to the MREmagnitude images. We note that in the MRE pulse sequence we used, the different ac-quisitions in each direction are multiplexed, and therefore, in the motion compensationalgorithm, we have registered the mean volumes over all acquisitions for each direction,rather than each individual acquisition.805.4. Discussion and Conclusionx [mm]y[mm]−60 −40 −20 0 20 40 60−60−40−200204060(a) Native resolutionx [mm]y[mm]−60 −40 −20 0 20 40 60−60−40−200204060(b) Super-resolutionx [mm]y[mm]−60 −40 −20 0 20 40 60−60−40−200204060(c) T2-weightedFigure 5.4: Super-resolution of a patient’s prostate MRE. Cross-section of the elastogramsand T2w image. Manual segmentation (red) of the prostate, and the inside and outsidebands that were used for contrast evaluation. Intensity scale of the elastograms is 0-30 kPa.We then employed the proposed super-resolution technique on the patients’ data. Qual-itative results on a mid-gland slice of one case are presented in Figure 5.4. We may noticethat on the super-resolution image, the prostate is more separable from its surroundingtissue than on the native resolution image. In order to quantify the results, we utilizethe corresponding segmented T2w volume. For evaluation of the gland contrast and edgestrength in the elastograms, we measure the contrast-to-noise ratio [17, 143] and the meangradient on the surface, i.e.,CNR =2(cin − cout)2/(σ2in + σ2out)MeanEdge =∮Γ‖∇E‖ dr¯/∮Γdr¯ (5.3)where cin, σin and cout, σout are the mean and STDV values of the elastogram E on twobands, inside and outside of the segmentation surface Γ, as illustrated in Figure 5.4(c). Inthe computation of the edge strength, a 3× 3× 3 Gaussian filtering was performed on theelastograms for noise reduction in the gradients.5.4 Discussion and ConclusionWe have presented an image processing approach for enhancing MRE data. The approachemploys image registration to compensate for motion of patients during the, typicallylengthy, acquisition process, in order to improve the accuracy of the reconstructed elas-togram. Then, a super-resolution technique is employed to increase the resolution of the815.4. Discussion and ConclusionNCC Difference (%) CNR MeanEdge (kPa/mm)Patient y z Native Super-res. Native Super-res.P01 1.12± 0.08 2.27± 0.11 1.1 1.2 3.1 3.8P04 5.34± 0.58 1.21± 0.61 5.5 7.8 2.4 3.5P05 0.93± 0.15 4.76± 0.16 8.9 13.2 1.3 2.5P06 1.21± 0.38 4.88± 0.51 1.8 2.5 1.5 2.7P08 2.22± 0.12 4.28± 0.15 9.3 12.5 2.3 2.8P10 0.23± 0.05 0.14± 0.06 23.1 20.9 3.9 3.3P12 1.62± 0.06 3.55± 0.10 12.9 7.3 1.2 1.6P15 0.69± 0.01 0.11± 0.01 0.6 4.3 2.3 3.2P17 0.20± 0.06 1.51± 0.07 9.9 10.4 2.7 3.6P18 0.13± 0.03 0.71± 0.09 7.0 7.7 2.6 2.5P19 1.03± 0.12 3.36± 0.09 10.1 13.6 2.3 3.0P20 0.52± 0.40 2.44± 0.43 1.9 3.7 1.2 1.8Mean 1.33± 0.55 2.84± 0.74 7.7± 6.4 8.8± 5.7 2.2± 0.8 2.9± 0.7Table 5.1: Super-resolution results on clinical data. The increase in NCC of the magnitudeimages in directions y, z (direction x is the reference) resulting from motion compensation.The contrast-to-noise ratio and edge strength of the prostate in the native and super-resolution elastograms (higher is better).registered phase images. The proposed approach utilizes unique properties of the MREacquisition process, namely, the similarity of the “by-product” magnitude images used forcompensating motion in the phase images, and the redundancy of the phase images usedfor super-resolution.We note that, in MRI, super-resolution has been proven to be of limited benefit in-planedue to the inherent bandlimited nature of the Fourier encoded images [144]. However, in ourapproach, we employ super-resolution on the phase component of the signal, and assumethat the nonlinear residual physical movements of the object (rather than shifts in FOV),and synchronization errors in MRE allow for higher frequency components to be recovered.Indeed, in [132], the authors proposed similar nonlinear deformations for in-plane super-resolution in diffusion MRI.We have tested the proposed approach on phantom and clinical data of CaP patients.Indeed, the NCC was increased after registration in our motion compensation approach.However, in some cases, the improvement is not significant, which implies that patient’s825.4. Discussion and Conclusionmotion was low. In a recent study by our group, [137], the effect of motion compensationon elasticity-based CaP detection was examined. The results are mixed, with classificationof lesions being improved on average in the peripheral zone, where the majority of cancerlesions are, but worsened in the central zone. This may be due to the fact that the gland isstiffer in the central zone, while we assumed a homogeneous elasticity regularization. Futurework may include an iterative registration process that uses the inhomogeneous elasticityof the reconstructed elastogram in each iteration, similar to [103], for improving estimationaccuracy of the deformations.The evaluation of the super-resolution method on an elasticity phantom showed that,on average, the super-resolution elasticity values were closer to the reported factory values.However, the values of the stiffest inclusion are lower in the super-resolution elastogram.This result might be due to a higher wavelength-to-voxel size ratio that may not be suitablefor optimal reconstruction. In general, on both phantom and clinical data experiments, thegeometry of the inclusions/prostate in the super-resolution volumes is shown better than innative resolution. In some cases, however, the improvement is not significant, and on somesuper-resolution images there are image artifacts and patterns, typically around the baseand apex, that might be caused by registration error. Nevertheless, on average, the improvedcontrast and edge strength around the prostate allow better separability of the gland fromits surrounding tissue, and may facilitate automatic segmentation algorithms. We note thatin the displacement estimation for reconstruction of the super-resolution elastogram, onlyfour samples are used (pairs of registered phase and anti-phase images), rather than eightsamples used in the native resolution elastogram. The SNR in the super-resolution imageis therefore improved due to the specific utilization of the information.For simplicity and due to a lack of prior information on the scanner imaging process,we have employed a Gaussian blur in the proposed super-resolution algorithm. However, asfuture work, the scanner’s actual point spread function of phase might be calculated similarto [127] in order to accurately model the blurring. Future work may also include improvedreconstruction methods using FEM combined with sparsity regularization [57]. In addition,increasing the vibration frequency to reduce the wavelength-to-voxel size ratio, or acqui-sition of multi-frequency MRE, may allow for better elasticity reconstruction. However,transperineal excitation beyond 100 Hz may not be feasible due to the high attenuation ofthe waves near the perineum.83Chapter 6Optimization-Based Design ofMotion Encoding in MagneticResonance Elastography6.1 IntroductionThe most prevalent dynamic MRE technique involves the application of a harmonic me-chanical excitation that generates shear waves in the tissue. The resulting displacementsare estimated by employing an MRI pulse sequence with MEGs. These MEGs are bipolargradients (typically in shape of a sinusoid or trapezoid) that are synchronized to the me-chanical excitation. The phase accumulated by each isochromat is encoded in the phasecomponent, also known as the phase image, of the MRI signal.Conventional motion encoding in MRE comprises unidirectional MEGs, in which eachcomponent of the displacement vector is encoded separately by applying the gradients inthe corresponding direction. The measured phase accumulation, i.e., each voxel of the phaseimages, is a sinusoid with respect to a phase shift between the MEG and the mechanicalwave. Thus, a sufficient number of consecutive acquisitions with varying phase shifts arerequired in order to sample the temporal characteristics of the wave. The amplitude andphase of the displacement at each voxel are then estimated through a DFT of the samples.Elastograms, i.e., images that depict tissue elasticity or stiffness can be reconstructed fromthe displacements by solving an inverse problem using, e.g., [57, 85, 107, 112, 151].In [166], a generalized sampling method is proposed in order to achieve a sub-Nyquistsampling rate of two acquisitions per direction using a least-squares approach, thus reducingacquisition time of unidirectional sequences. In practice, however, due to noise in themeasured phase images, MRE sequences typically require four or eight acquisitions in eachdirection in order to provide a displacement map and, in turn, an elastogram, with areasonable SNR. This yields a longer acquisition time that hinders usage of MRE in somecases due to time constraints and limits on image quality. Longer acquisition times also846.1. Introductionintroduce artifacts in the images due to patient movement.There are approaches that have been proposed to improve SNR or, alternatively, re-duce the acquisition time of MRE [46]. Most approaches propose improved imaging pulsesequences, such as EPI [85] and gradient-echo [44], as opposed to the basic SE sequence.Other approaches deal directly with the design of the MEGs, e.g., fractional encoding[56, 131], in which the mechanical frequency is taken as a fraction of the MEGs oscillationfrequency. In [54, 116], sequence designs with combination of MEGs were applied to phase-contrast MRI, a field related to MRE, in order to reduce scan time and increase sensitivityfor estimating the magnitude of flow in MR angiography.Recently, [72, 173] have proposed approaches for encoding the 3-D displacements in MREby applying the MEGs in all directions simultaneously. These multidirectional sequenceshave the potential of reducing the number of scans by three, while achieving the same SNRas in unidirectional sequences. In [173], MEGs are applied simultaneously with differentfrequencies, and the displacement is encoded for each direction by utilizing the filteringcondition in [133]. In [72], the authors propose a sample interval modulation MRE, inwhich MEGs are applied simultaneously with different phase shift intervals, such that thephase accumulation in each direction can be recovered from the total accumulation by usingDFT.In this chapter we focus on optimization of the MEGs design. We parameterize generalmultidirectional sequences, and formulate the estimation of the 3-D displacement as a linearsystem in the time-domain. This allows derivation of a least-squares estimator that isoptimal under a given design and noise model. We then characterize and set the problemin an experimental design framework, by which we quantify the performance of MEGssequences, and derive parameters that lead to optimal multidirectional designs under a givenoptimality criterion. We fit the conventional MEGs sequence in the proposed framework,and compare the sequences on simulation and phantom data.The remainder of this chapter is organized as follows. First, we present the generalizedMEGs and corresponding phase accumulation, and propose a least-squares estimation ofwave properties. We then present the framework for optimal design of motion encoding,and derive optimal sequences. We continue with a description of the experiments that wereperformed to evaluate our approach, and summarize the results. Finally, we conclude withdiscussion and future research directions. A summary of the notations used throughout thiswork is provided in Table 6.1.856.1.IntroductionNotation Description Typical valuesr¯ = (x, y, z)T spatial coordinate vector in mt time in sγ/(2pi) the gyromagnetic ratio 42.58MHz/Ti spatial index x, y, zωm/(2pi) mechanical frequency 200Hzξi displacement amplitude in direction i To be estimated (in the order of µm)ϕi displacement phase in direction i To be estimated ([−pi, pi))ωi/(2pi) MEG oscillation frequency 200HzTi MEG period time (Ti = 2pi/ωi) 5msqi fractional encoding ratio qi = ωm/ωi qi ≤ 1|Gi| MEG amplitude 40mT/mαi phase shift between MEG and vibration [0, 2pi)βi start phase of the MEGs [0, 2pi)δi net phase (δi = αi − βi) (−2pi, 2pi)Ni number of periods MEGs are applied 1WNiTi(t) a window function with a [0, NiTi] support 1 for t ∈ [0, NiTi], 0 otherwiseConventional a sequence comprising 24 unidirectional acquisitions3-D design (8 in each direction)RME a sequence comprising a minimal number of 6unidirectional acquisitions (2 in each direction)MD-RME a sequence comprising a minimal number of 6multidirectional acquisitions with sign-flipsOptimal design a sequence satisfying an optimality criterion(comprising 8 multidirectional acquisitions in our examples)Table 6.1: Summary and description of the notations.866.2. Theory6.2 Theory6.2.1 Parameterized Motion Encoding GradientsWe denote the spatial coordinate vector as r¯ = (x, y, z)T and the time as t. We representa steady-state harmonic mechanical wave through tissue by its 3-D displacement vectorΨ¯ = (Ψx,Ψy,Ψz)T . Without loss of generality, we writeΨi(r¯, t) = Im{ |ξi(r¯)|ejϕi(r¯) ejωmt } = |ξi(r¯)| sin (ωmt+ ϕi(r¯)) , i = x, y, z (6.1)The wave frequency ωm is the same as the known mechanical excitation frequency. Thespatially varying amplitudes and phases in each direction, |ξi| and ϕi, are unknown andhave to be determined at each imaged voxel based on a number of phase images acquiredin different acquisitions. Note that this representation can accommodate locally any typeof wave (e.g., traveling or standing), wavefront (e.g., plane or spherical), polarization (e.g.,linear or circular) and attenuation.We parameterize general multidirectional sequences of sinusoidal MEGs, G¯ = (Gx, Gy, Gz)T .Specifically, for each direction, we allow control of the (angular) frequency ωi, the amplitude|Gi|, the start phase βi, the phase shift relative to the mechanical wave αi, and the numberof periods the gradients are applied Ni ∈ N. We therefore haveGi(t;ωi, Ni, |Gi|, αi, βi) = |Gi| sin (ωi(t− αi/ωi) + βi)WNiTi(t− αi/ωi) , i = x, y, z (6.2)where WNiTi(t) is a window function with a [0, NiTi] support, and Ti = 2pi/ωi is the timeperiod of each MEG.The effect of the parameters αi and βi on the MEGs is demonstrated on a single acqui-sition in Figure 6.1. Geometrically, this may be seen as if the gradients are rotating in a(piecewise) helical fashion as a function of time in gradient space.6.2.2 Parameterized Phase AccumulationIn this chapter, we limit our discussion to MEGs with the same frequency as the mechanicalfrequency, i.e., ωi = ωm, ∀i = x, y, z. However, the theory can be extended to fractionalencoding in each direction.The phases accumulated by the precession of each moving isochromat during the ac-quisition sequence, as result of the MEGs, are encoded in the measured phase images ofthe MRI signal. As derived in the Appendix, the phase accumulation at each voxel as a876.2. Theory−10 −5 0 5 10 15 20t [ms]MechGxGyGz(a) Infinite slew rate−40−20 0 20 40 −40 −200 2040−40−30−20−10010203040Gy [mT/m]t=0Gx [mT/m]Gz[mT/m](b) Infinite slew rate−10 −5 0 5 10 15 20t [ms]MechGxGyGz(c) Finite slew rate−40−20 0 20 40 −40 −200 2040−40−30−20−10010203040Gy [mT/m]t=0Gx [mT/m]Gz[mT/m](d) Finite slew rateFigure 6.1: Example of MEGs sequence with α¯ = (0, pi2 , 0), β¯ = (0, pi, 3pi2 ). (a),(c) Plots ofthe MEGs in time with infinite and finite slew rate. (b),(d) Plots of the sequence in gradientspace with infinite and finite slew rate.function of the MEGs parameters is given byΦ(r¯; N¯ , |G¯|, α¯, β¯) =∑i=x,y,zγpiNi|Gi|ωm|ξi(r¯)| cos (ϕi(r¯) + αi − βi) (6.3)where γ is the gyromagnetic ratio. This equation is a generalization of the phase accumu-lation that was originally presented in [100].In order to formulate the problem as a linear one, for every voxel, we define Ai, Bi tobe such that|ξi(r¯)| =√A2i (r¯) +B2i (r¯) , ϕi(r¯) = arctan(−Bi(r¯)Ai(r¯)), i = x, y, z (6.4)886.2. Theoryand, using the harmonic addition theorem, substitute in (6.3) to obtainΦ(r¯; N¯ , |G¯|, α¯, β¯) =∑i=x,y,zγpiNi|Gi|ωm[Ai(r¯) cos (αi − βi) +Bi(r¯) sin (αi − βi)] (6.5)Thus, for every voxel, estimation of the amplitudes |ξi| and phases ϕi in each direction, isgiven by estimation of the six unknowns Ai and Bi, i = x, y, z. This requires a minimumof six acquisitions, i.e., six acquisitions of the phase image with different MEG parameters.6.2.3 Optimal Estimation of Wave PropertiesIn general, we acquire K ≥ 6 acquisitions with, possibly, different α(k)i , β(k)i , N(k)i and|G(k)i | in each direction i = x, y, z, and for each acquisition k = 1, 2, . . . ,K. We denoteδ(k)i = α(k)i − β(k)i as the net phase of each acquisition, and define the K × 6 (“tall”)acquisition design matrixM = γpiωmN (1)x |G(1)x | cos (δ(1)x ) N (2)x |G(2)x | cos (δ(2)x ) · · · N (K)x |G(K)x | cos (δ(K)x )N (1)x |G(1)x | sin (δ(1)x ) N (2)x |G(2)x | sin (δ(2)x ) · · · N (K)x |G(K)x | sin (δ(K)x )N (1)y |G(1)y | cos (δ(1)y ) N (2)y |G(2)y | cos (δ(2)y ) · · · N (K)y |G(K)y | cos (δ(K)y )N (1)y |G(1)y | sin (δ(1)y ) N (2)y |G(2)y | sin (δ(2)y ) · · · N (K)y |G(K)y | sin (δ(K)y )N (1)z |G(1)z | cos (δ(1)z ) N (2)z |G(2)z | cos (δ(2)z ) · · · N (K)z |G(K)z | cos (δ(K)z )N (1)z |G(1)z | sin (δ(1)z ) N (2)z |G(2)z | sin (δ(2)z ) · · · N (K)z |G(K)z | sin (δ(K)z )T(6.6)By (6.5), for each voxel r¯ we have a (typically overdetermined) linear system of equationsΦ¯ =MA¯+ w¯ (6.7)where Φ¯ = (Φ(1), . . . ,Φ(K))T is a vector that comprises the measured phase accumulationsin each acquisition, A¯ = (Ax, Bx, Ay, By, Az, Bz)T is the unknown coefficient vector that weneed to solve the system for, and w¯ is a random equation error that depends on the typeof noise in the system.We can re-write (6.7) as the estimation problemˆ¯A = arg minA¯‖Φ¯−MA¯‖ , given M, Φ¯ (6.8)The distribution of w¯, which can be interpreted as the noise introduced into the measuredphase images during each acquisition, determines the type of norm in (6.8), and the choice896.2. Theoryof estimator ˆ¯A. We model w¯ as a white Gaussian noise (see Discussion section) witha variance of σ2, and assume that M is full column rank. Under these assumptions, theGauss-Markov theorem [65] states that the ordinary least-squares estimator is optimal (bestlinear unbiased estimator) in the sense of minimizing the mean squared error (MSE) in theestimation E‖ ˆ¯A− A¯‖2. The least-squares estimator is defined byˆ¯A = (MTM)−1MT Φ¯ =M†Φ¯ (6.9)where M† is the Moore-Penrose pseudoinverse of the matrix M.6.2.4 Optimal Design of Motion EncodingThe proposed MEGs parametrization provides several degrees of freedom that allow us tomodify the acquisition design matrix in (6.6). The parameters should be diverse amongacquisitions to yield linearly independent rows, such that M is full column rank. Clearly,the more acquisitions we acquire (larger K), the lower is the estimation error that theleast-squares estimator obtains. However, this would increase acquisition time, which isundesirable as explained above.Another option to reduce error is to increase the MEGs amplitude |Gi|, the number ofperiods that they are applied Ni, or to decrease their frequency ωi. However, the allowablerange for the parameters above is rather small due to hardware limitations, acquisition time,and the size of the scanned target (such that enough wavelengths are imaged). Therefore,we focus on finding the parameters αi and βi that optimize the acquisition design matrixin the sense of minimizing the estimation error by the least-squares estimator in (6.9).The problem above can be considered as a special case of the experimental design problemin statistics [69, 124]. A popular criterion for an optimal design is the so called D-optimality[18, 63], in which the determinant of MTM, also known as the information matrix of thedesign, is maximized. Such D-optimum design minimizes the (generalized) variance of thebest linear estimator [68]. Thus, maximizing the determinant of the information matrix isequivalent to minimizing the estimation error.We can therefore formulate the D-optimal design problem in MRE as(α(k)i , β(k)i ) = arg maxδ(k)idet (MTM) , i = x, y, z , k = 1, 2, . . . ,Kgiven ωm, N (k)i , |G(k)i |, and possible conditions on α(k)i , β(k)i (6.10)As we typically do not favor any particular direction or acquisition, unless stated otherwise,906.2. Theorywe choose the same N (k)i = N and |G(k)i | = |G|, ∀i = x, y, z, k = 1, . . . ,K.While designs that yield the same determinant are theoretically equivalent in theirestimation error, in practice, issues such as TE and TR increase, finite slew-rate, anddifferent flow compensation characteristics affect SNR. Thus, in some cases, to minimizethe TE of the sequence in each acquisition, we set the same α(k)i = α(k), ∀i = x, y, z. Inother cases, to prevent discontinuities (“jumps”) in the MEGs, we set β(k)i = {0, pi}.Looking at (6.6), we have a constant multiplied by a matrix with elements cos (δ(k)i ) andsin (δ(k)i ). In general, finding an L × L matrix that maximizes the (absolute) determinantamong all L × L matrices with elements with an absolute value less or equal to one, isthe well-known Hadamard’s maximum determinant problem [50]. It turns out that such(non-unique) matrix HL, called a Hadamard matrix, is known to exist on the real field forL = 1, 2, L ≡ 0(mod 4), and its elements take the extreme values ±1.For cases of L in which a Hadamard matrix does not exist, there may still be a D-optimum design matrix [42]. However, it is typically of more value to employ the closestH-matrix [27], which is a non-square matrix K × L with K > L that comprises elementsof ±1, and satisfies HTK×LHK×L = KIL, with IL being the L × L identity matrix. Suchmatrices are D-optimum with respect to a K × L design with det (HTK×LHK×L) = KL.6.2.5 Unidirectional DesignsHere we fit a unidirectional MEGs sequence designs into the proposed framework. Thisis done in order to compare and evaluate the multidirectional designs that are presentednext. This special case is also important as it shows that our framework allows estimationof the displacement map using a unidirectional design with a minimum of two acquisitionsfor each direction, as suggested in [166].In order to apply only one MEG direction in each acquisition, we toggle the MEGsamplitudes as |G(k)i | = {0, |Gi|}. We follow the result of [166], which states that a minimalerror with respect to sampling of the phase images is attained with evenly distributedacquisitions (in each direction) over the [0, pi) range for two acquisitions per direction.Thus, we limit our choice to α(k)i = {0, pi2 }. Finally, we choose zero start phases β(k)i = 0.916.2. TheoryWe therefore define a unidirectional design with a minimum of six acquisitions as|G(1,...,6)x | = (|G|, |G|, 0, 0, 0, 0)|G(1,...,6)y | = (0, 0, |G|, |G|, 0, 0)|G(1,...,6)z | = (0, 0, 0, 0, |G|, |G|)α(1,...,6)i = (0, pi/2, 0, pi/2, 0, pi/2) , i = x, y, zβ(1,...,6)i = 0 , i = x, y, z (6.11)This sequence design, referred to as reduced motion encoding (RME) MRE in [166], isillustrated in Figure 6.2, and by (6.6) it yields a simple diagonal design matrixMrme =γpiN |G|ωmI6 , with[det (MTrmeMrme)] 12 =(γpiN |G|ωm)6(6.12)For unidirectional designs that comprise more than two acquisitions per direction, [166]shows that αi should be evenly distributed over the [0, 2pi) range in each direction. Similarto (6.11), we can construct such sequences and compute their design matrices. We considera conventional design of 24 acquisitions (8 in each direction) for our evaluation and denoteits design as Mcon. We note that by explicitly writing the proposed least-squares solutionin (6.9) forMcon, it can be shown that the estimated wave properties, |ξi| and ϕi, using ourtime-domain approach coincide with the results achieved by a standard DFT processing ofthe phase images, in which |ξi| and ϕi are estimated by the amplitude and phase of thefundamental frequency.6.2.6 Multidirectional RME DesignHere we are interested in generating a multidirectional MEGs sequence that takes advantageof the parameterized approach in order to reduce estimation error. We note that employingdifferent phase shifts αi 6= αj , in two directions or more, increases the TE of the acquisitionso as to cover the entire MEGs sequence, which decreases SNR due to a weaker signal. TheTE is measured based on the union of the time the MEGs are applied in a single acquisition,rather than from t = 0, therefore TE is still minimal even if αi 6= 0, as long as it is the samein all directions.In addition, employing a start phase βi 6= {0, pi} produces discontinuities in the MEGswaveform. In practice, gradient amplifiers impose a finite slew rate, and limitations onspecific absorption rate restrict excessive dB/dt. This introduces errors to the phase ac-cumulation calculation that lead to estimation errors in the displacement. The effect of a926.2. Theory−10 −5 0 5 10 15 20t [ms]MechGxGyGz(a) Acquisition 1−10 −5 0 5 10 15 20t [ms]MechGxGyGz(b) Acquisition 2−10 −5 0 5 10 15 20t [ms]MechGxGyGz(c) Acquisition 3−10 −5 0 5 10 15 20t [ms]MechGxGyGz(d) Acquisition 4−10 −5 0 5 10 15 20t [ms]MechGxGyGz(e) Acquisition 5−10 −5 0 5 10 15 20t [ms]MechGxGyGz(f) Acquisition 6Figure 6.2: RME design. Comprises a minimum number of unidirectional acquisitions.finite slew rate on the sequence can be observed in Figure 6.1.Thus, in this example we add to the optimal design problem in (6.10) the constraintsβi = {0, pi} and αi = α, ∀i = x, y, z. We propose the following design with sign-flipsα(1,...,6)i = (0, pi/2, 0, pi/2, 0, pi/2) , i = x, y, zβ(1,...,6)x = (pi, pi, 0, 0, 0, 0)β(1,...,6)y = (0, 0, pi, pi, 0, 0)β(1,...,6)z = (0, 0, 0, 0, pi, pi) (6.13)Since it comprises the minimum number of acquisitions, we refer to it as a multidirectionalreduced motion encoding (MD-RME) design. The sequence is illustrated in Figure 6.3, and936.2. Theoryyields a Toeplitz design matrixMmd-rme =γpiN |G|ωm−1 0 1 0 1 00 −1 0 1 0 11 0 −1 0 1 00 1 0 −1 0 11 0 1 0 −1 00 1 0 1 0 −1,with[det (MTmd-rmeMmd-rme)] 12 = 16(γpiN |G|ωm)6(6.14)Thus, for the same minimum number of acquisitions, the determinant value of the estimationcovariance is reduced by a factor of 16 in comparison with the unidirectional RME sequence.We note that, since we are interested in the estimation error of each couple of coefficients,Ai, Bi, we should look at the elements in the diagonal of the inverse information matrix(MTmd-rmeMmd-rme)−1 that correspond to the variances of each coefficient. We find thatthese elements are half of their value in the unidirectional RME case, thus we expect theSTDV of the estimation error to be reduced by√2 using the MD-RME sequence.6.2.7 Optimal Design with Varying Start PhasesHere, rather than proposing a sequence and evaluating its performance, we take the inverseapproach where we find a design that produces an optimal acquisition matrix. In our caseL = 6 ≡ 2(mod 4), a Hadamard matrix does not exist, and the closest H-matrix is forK = 8. Any such (non-unique) matrix can be chosen, and, using the technique in [42], weconstructedH8,6 =1 1 −1 1 −1 −1 1 −11 1 1 −1 −1 −1 −1 11 −1 1 1 −1 1 −1 −11 1 −1 1 1 1 −1 11 1 1 −1 1 1 1 −11 −1 1 1 1 −1 1 1T(6.15)Given constant gradient amplitude and number of periods, the matrixH8,6 suggests thatδ(k)i = {−3pi/4,−pi/4, pi/4, 3pi/4}. The technique to determine δ(k)i , for each acquisition kand direction i, is to observe the pair of elements (k, 2i − 1), (k, 2i). Based on (6.6),these elements are proportional to cos (δ(k)i ) and sin (δ(k)i ), respectively. Thus, the value of946.2. Theory−10 −5 0 5 10 15 20t [ms]MechGxGyGz(a) Acquisition 1−10 −5 0 5 10 15 20t [ms]MechGxGyGz(b) Acquisition 2−10 −5 0 5 10 15 20t [ms]MechGxGyGz(c) Acquisition 3−10 −5 0 5 10 15 20t [ms]MechGxGyGz(d) Acquisition 4−10 −5 0 5 10 15 20t [ms]MechGxGyGz(e) Acquisition 5−10 −5 0 5 10 15 20t [ms]MechGxGyGz(f) Acquisition 6Figure 6.3: MD-RME design. Comprises a minimum number of multidirectional acquisi-tions, without discontinuities in the MEGs or increase in TE.δ(k)i is determined by the quadrant defined by these elements (e.g., the pair (1,−1) yieldsδ = −pi/4).We can derive several sequences that produce H8,6. Here, in order to minimize TE, weare interested in having no phase shift between directions, i.e., α(k)i = α(k), ∀i = x, y, z. Itcan be proved that no combination with only β(k)i = {0, pi} can produce H8,6. Thus, we areforced to have discontinuities in the MEGs in this case, and in order to minimize the slewrate effect we use β(k)i = {0, pi} when possible.We propose the following sequenceα(1,...,8)i = (pi4 ,3pi4 ,3pi4 ,3pi4 ,5pi4 ,5pi4 ,5pi4 ,3pi4 ) , i = x, y, zβ(1,...,8)x = (0,pi2 , 0, pi, 0, 0,3pi2 , 0)β(1,...,8)y = (0, 0, pi,pi2 ,pi2 , pi, 0, 0)β(1,...,8)z = (0, pi,pi2 , 0, pi,3pi2 , pi, 0) (6.16)956.2. Theory−10 −5 0 5 10 15 20t [ms]MechGxGyGz(a) Acquisition 1−10 −5 0 5 10 15 20t [ms]MechGxGyGz(b) Acquisition 2−10 −5 0 5 10 15 20t [ms]MechGxGyGz(c) Acquisition 3−10 −5 0 5 10 15 20t [ms]MechGxGyGz(d) Acquisition 4−10 −5 0 5 10 15 20t [ms]MechGxGyGz(e) Acquisition 5−10 −5 0 5 10 15 20t [ms]MechGxGyGz(f) Acquisition 6−10 −5 0 5 10 15 20t [ms]MechGxGyGz(g) Acquisition 7−10 −5 0 5 10 15 20t [ms]MechGxGyGz(h) Acquisition 8Figure 6.4: Optimal design with varying start phases. Notice the discontinuities in theMEGs. Using phase shifts we achieve discontinuity of only one direction per acquisition,without increase in TE.The sequence is illustrated in Figure 6.4, and yields the design matrixMopt =γpiN |G|√2ωmH8,6 ,with[det (MToptMopt)] 12 = 64(γpiN |G|ωm)6(6.17)Note we have a 1/√2 loss (cosine and sine of multiples of pi/4), but the overall gain of 64of the optimal design compared to the RME design is still significant. By looking at thediagonal of the inverse information matrix, we find that its elements are 0.25 of their valuein the RME case, thus we expect the STDV of the estimation error to be reduced by 2 bythe optimal design.Remember that we have used 8 acquisitions for this design, rather than 6 as in theprevious designs. In case we employ 8 acquisitions of an RME sequence (by replicating twoacquisitions), we have that the optimal design has a gain of 32 relative to such a design.In fact, the gain of this optimal design with only 8 acquisitions equals the gain of theconventional design Mcon, which comprises 24 acquisitions.6.2.8 Optimal Design with Varying Phase ShiftsIn the previous optimal design, the single discontinuity in some acquisitions may hinder itsusage in practice due to a limited slew rate. In addition, as pointed out in [72], different966.3. Methodsstart phases may affect the SNR in in vivo environments, since MEGs with different startphase have different flow compensation characteristics. Thus, here we propose a sequencewith a start phase that is limited to β(k)i = {0, pi}, but with different phase shifts α(k)ifor each direction. Such multidirectional approach was proposed in [72] as sample intervalmodulation MRE to modulate the sampling interval such that the phase accumulation ineach direction has a different frequency.As described above, to yield the same optimal design as in (6.17), we derive the valuesof α(k)i such that δ(k)i are in the corresponding quadrants determined by H8,6. We proposethe sequenceα(1,...,8)x = (pi4 ,pi4 ,3pi4 ,3pi4 ,pi4 ,pi4 ,3pi4 ,3pi4 )α(1,...,8)y = (pi4 ,3pi4 ,3pi4 ,pi4 ,3pi4 ,pi4 ,pi4 ,3pi4 )α(1,...,8)z = (pi4 ,3pi4 ,pi4 ,3pi4 ,pi4 ,3pi4 ,pi4 ,3pi4 )β(1,...,8)x = (0, 0, 0, pi, pi, pi, pi, 0)β(1,...,8)y = (0, 0, pi, 0, 0, 0, pi, 0)β(1,...,8)z = (0, pi, 0, 0, 0, pi, 0, 0)(6.18)The sequence is illustrated in Figure 6.5. As discussed in [72], a disadvantage of employ-ing different phase shifts for each direction is an increase of TE. Using the combinations ofα(k)i and β(k)i above, the phase shifts between the directions in each acquisition are mini-mized such that the maximum phase shift is pi/2, which translates into a TE increase by aquarter of the time period.6.3 Methods6.3.1 Simulation of a Plane WaveWe first test our approach by simulating a linearly polarized wave propagating through a3-D volume of size 64×64×64 that comprises a background with a spherical inclusion. Thestiffer inclusion is associated with a higher propagation speed (lower wave number), and alower displacement amplitude than the background. We have used a shear wave frequencyas in the phantom experiment (200Hz), and chosen the propagation speeds (∼ 10m/s),and amplitudes (∼ 10µm) to be of the same order of magnitude as in the phantom exper-iment below. The goal of this simulation is to validate the least-squares wave estimation976.3. Methods−10 −5 0 5 10 15 20t [ms]MechGxGyGz(a) Acquisition 1−10 −5 0 5 10 15 20t [ms]MechGxGyGz(b) Acquisition 2−10 −5 0 5 10 15 20t [ms]MechGxGyGz(c) Acquisition 3−10 −5 0 5 10 15 20t [ms]MechGxGyGz(d) Acquisition 4−10 −5 0 5 10 15 20t [ms]MechGxGyGz(e) Acquisition 5−10 −5 0 5 10 15 20t [ms]MechGxGyGz(f) Acquisition 6−10 −5 0 5 10 15 20t [ms]MechGxGyGz(g) Acquisition 7−10 −5 0 5 10 15 20t [ms]MechGxGyGz(h) Acquisition 8Figure 6.5: Optimal design with varying phase shifts. The maximum phase shift at eachacquisition dictates the increase in TE. Using sign-flips, this phase shift is at most pi2 .framework, and theoretical performance analysis of the proposed designs.We define a gradient sequence with the proposed parametrization of the MEGs, andcompute the associated phase images for each acquisition by integrating the product of thewave and MEGs over time. Note that, in order to validate our analytical expression for thephase accumulation, we compute the integral numerically. To model the imaging noise, weadd a zero mean Gaussian noise to each voxel of the phase images.We repeat the simulations for the unidirectional RME, MD-RME and optimal designsproposed above, each with a range of noise variances. We calculate the coefficient vector ˆ¯Aby using the least-squares estimator, and derive the estimated amplitude and phase of thewave, |ξˆi|, ϕˆi in each direction and for each voxel. Slices from the estimated volumes areshown in Figure 6.6.Since we know the ground truth of the wave phases and amplitudes, we may computethe estimation error as the MSE between the estimations and the ground truth, i.e.,ε2ξ =1|Ω|∫Ω√∑i=x,y,z(|ξˆi(r¯)| − |ξi(r¯)|)2dr¯ , ε2ϕ =1|Ω|∫Ω√ ∑i=x,y,z(ϕˆi(r¯)− ϕi(r¯))2 dr¯(6.19)where Ω and |Ω| are the volume and its size, respectively.986.3. MethodsAmplitudex [mm]y[mm]−30 −20 −10 0 10 20 30−30−20−100102030 x [mm]y[mm]−30 −20 −10 0 10 20 30−30−20−100102030 x [mm]y[mm]−30 −20 −10 0 10 20 30−30−20−100102030 x [mm]y[mm] −30 −20 −10 0 10 20 30−30−20−100102030 05101520253035404550Phasex [mm]y[mm]−30 −20 −10 0 10 20 30−30−20−100102030(e) RME (6 acquisi-tions)x [mm]y[mm]−30 −20 −10 0 10 20 30−30−20−100102030(f) MD-RME (6 acq.)x [mm]y[mm]−30 −20 −10 0 10 20 30−30−20−100102030(g) Optimal (8 acq.) orConventional (24 acq.)x [mm]y[mm] −30 −20 −10 0 10 20 30−30−20−100102030−3−2−10123(h) TrueFigure 6.6: Simulation results. The estimated amplitude and phase of the wave using theunidirectional RME, MD-RME and optimal design with minimum acquisitions for a STDVof σ = pi/5 of the noise added to the phase accumulation images. The conventional designyields similar results as the optimal design, but with three times the number of acquisitions.Top: The amplitude in the z-axis, intensity scale is in µm. Bottom: The phase in the z-axis(wrapped in [−pi, pi)), intensity scale is in rad.6.3.2 Phantom ExperimentWe have implemented the RME and the MD-RME sequences on a 3-Tesla system (Achieva3.0T, Philips, The Netherlands) using a SE-EPI pulse sequence with a TE/TR of 30/1600 ms.A 200 Hz vibration was applied to a rectangular shaped elasticity phantom (Model 049,CIRS, VA, USA) using a custom-made shielded electromagnetic transducer synchronizedto the scanner [135]. As seen in Figure 5.2(e), the phantom contains spherical inclusions ofdifferent radii and stiffness values relative to the background material that mimic differenttissue types. The entire phantom was scanned in a field of view of 128× 128× 20 grid, withan isotropic 2× 2× 2mm3 resolution. The scanning time was 364s per acquisition.In order to evaluate our approach, we have oversampled the scan of each acquisitionin both designs with the phase shifts of the MEGs ranging evenly over [0, 2pi) relative tothe mechanical excitation. A 80 × 60 × 20 region of interest on the acquired axial phaseaccumulation images was selected, compensated and phase unwrapped. We estimated thedisplacements using the proposed approach with a corresponding design. We have also gen-996.3. Methodsx [mm]y[mm] −80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050 05101520253035404550(a) Estimated wave amplitude (in µm)x [mm]y[mm] −80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050 −3−2−10123(b) Estimated wave phase (in rad)Figure 6.7: Phantom data. A cross-section of the best estimated amplitude and phase ofthe wave in the z-axis, as obtained by the least-squares approach over the entire set ofacquisitions, and used as the “ground truth” for evaluation.erated an elastogram using a LFE technique [85]. We note that these offline processing stepsare performed in the same manner as in standard processing of MRE, with the exception ofthe displacement estimation that involves multiplications of matrices and inversion of theinformation matrix instead of the standard DFT computation (total processing times aresimilar).We first estimated the displacements using the unidirectional and multidirectional de-signs with a minimum of 6 acquisitions each (RME and MD-RME). We then repeated forthe 12 and 24 acquisition versions of the designs (where we refer to the unidirectional de-sign with 24 acquisitions as the conventional MRE sequence), and evaluated the estimationerrors in the amplitude, phase and elasticity. Since we do not have the ground truth of thedisplacements, we have employed the best estimation obtained by using our least-squaresestimation approach over the entire set of acquisitions (i.e., a total of 8×6 acquisitions thatinclude both unidirectional and multidirectional designs). A cross-section of this best esti-mation’s amplitude and phase is shown in Figure 6.7. The MSE of each design is computedusing (6.19).Further evaluation of the approach is performed by comparing the elastograms to thereported elasticity values of the phantom. We segmented the inclusions, and calculated theMSE between the reconstructed elastograms and the elasticity of each type of tissue overthe corresponding region.1006.4. Results0 0.5 1 1.5102030405060708090σ [rad]Error[µm] RME (6 acquisitions)MD−RME (6 acq.)Optimal (8 acq.)Conventional (24 acq.)(a) Amplitude Error0 0.5 1 1.50.511.522.533.5σ [rad]Error[rad] RME (6 acquisitions)MD−RME (6 acq.)Optimal (8 acq.)Conventional (24 acq.)(b) Phase ErrorFigure 6.8: Simulation Errors. The MSE of the estimated amplitude and phase as a functionof the STDV σ of the noise added to the phase accumulation images. The lines and barsrepresent the mean and half of the STDV of the error over all voxels, respectively.6.4 Results6.4.1 Simulation of a Plane WaveThe error plots as a function of the STDV σ of the noise are shown in Figure 6.8. Wenotice linear and nonlinear relationships between the STDV and the MSE of the amplitudeand phase, respectively. As expected, the results using the conventional sequence with 24unidirectional acquisitions coincide with those of the optimal sequence with 8 acquisitions.By dividing the amplitude error plots, we find that the estimation errors of the RMEdesign are reduced by 1.51± 0.04 and 2.03± 0.06 using the MD-RME and optimal designs,respectively. Similar slope ratios were computed for the phase error plots in the linear regionthat corresponds to small values of STDVs. These results fit the theoretical error reductionsof√2 and 2 in the estimation of the coefficients, to which, by (6.4), the amplitude and phaseare proportional linearly and nonlinearly.6.4.2 Phantom ExperimentThe quantitative results are summarized in Table 6.2. The ratio between the amplitudeestimation errors of the unidirectional and the multidirectional with sign-flips designs forthe 6, 12 and 24 acquisitions experiments were 4.05, 1.19 and 1.47, respectively. Thus, theerror reductions by the multidirectional designs are close to the theoretical√2 value thatwas calculated from the information matrices, and in the 6 acquisitions case significantly1016.5. Discussion and ConclusionUnidirectional designs Multidirectional designs6 Acq. 12 Acq. 24 Acq. 6 Acq. 12 Acq. 24 Acq.(RME) (Conventional) (MD-RME)Amplitude Error (µm) 46.6± 33.9 2.2± 1.8 2.1± 1.7 11.5± 8.3 1.8± 1.5 1.4± 1.1Phase Error (rad) 1.4± 0.9 0.4± 0.4 0.3± 0.4 0.9± 0.8 0.2± 0.3 0.2± 0.3Elasticity Error (kPa) 29.6± 22.6 3.81± 3.3 3.7± 3.1 14.9± 13.1 2.6± 2.2 1.7± 1.5Table 6.2: Phantom results. Estimation errors in comparison with the best estimate ob-tained by the entire set of acquisitions.surpasses it.The elastograms that were reconstructed based on the unidirectional and multidirec-tional designs for 6, 12 and 24 acquisitions are illustrated in Figure 6.9. The correspondingerrors in comparison with the factory reported elasticity are presented in Figure 6.10. Wenotice that, as expected, the more acquisitions we have, the closer are the elasticity valuesto their factory values. Also, on average, the values generated from the multidirectionaldesign are closer to their reported values.6.5 Discussion and ConclusionIn this chapter we have presented a framework for motion encoding design in MRE. Theapproach was tested and evaluated using simulation and experiments on a phantom. Thesimulation confirms that optimal sequences have the potential of reducing the number ofacquisitions required by unidirectional MRE up to a factor of three while acquiring imageswith the same SNR. Alternatively, for the same number of acquisitions, the estimation errorof displacement amplitude can be reduced up to a factor of two. The produced elastogramsare also more accurate, since reconstruction is positively correlated with the accuracy ofthe displacement input.A simplified analogy of the problem is the Yates-Hotelling weighing problem, in whichit was shown that by using the same number of weighings, individual weights of objects canbe estimated better by weighing them in combinations, rather than separately [96]. Here,each acquisition is considered a “weighing,” the measured phase accumulation are the “scalereadings,” and using the MEGs parameters we assign each direction a “pan on the scale”or exclude it.1026.5.DiscussionandConclusionU ni di r ec ti on alx [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(a) RME (6 acquisitions)x [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(b) 12 Acquisitionsx [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(c) Conventional (24 Acq.)Mu lt i di r ec ti on alx [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(d) MD-RME (6 acquisitions)x [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(e) 12 Acquisitionsx [mm]y[mm]−80 −60 −40 −20 0 20 40 60 80−50−40−30−20−1001020304050(f) 24 AcquisitionsFigure 6.9: Phantom results. The elastograms (mean over three central slices) that were reconstructed from the estimatedamplitude and phase of the wave with a different number of acquisitions. Top: Using the unidirectional design. Bottom:Using the multidirectional design. Intensity scale is in kPa and same as in Figure 5.2(e).1036.5.DiscussionandConclusion5305580105130155180205Tissue typeElasticity[kPa] 0 1 2 3 4RME (6 acq.)MD−RME (6 acq.)True(a) 6 Acquisitions5152535455565Tissue typeElasticity[kPa] 0 1 2 3 4Unidirectional (12 acq.)Multidirectional (12 acq.)True(b) 12 Acquisitions5152535455565Tissue typeElasticity[kPa] 0 1 2 3 4Conventional (24 acq.)Multidirectional (24 acq.)True(c) 24 AcquisitionsFigure 6.10: Phantom elasticity errors. Errors between the elastograms of the unidirectional and multidirectional designs incomparison with the true elasticity on each tissue type.1046.5. Discussion and ConclusionWith scanner limitations on the slew rate of the gradients, an optimal design withvarying start phases may cause errors in phase accumulation by encoding higher frequencycomponents. Alternative implementations of optimal designs that are not affected by theslew rate, but require longer TE, can be achieved with varying phase shifts as proposed in[72]. Indeed, our framework shows that the sequence in [72] yields an optimal design underthe noise model and optimality criterion we assumed, and their DFT approach coincideswith the proposed least-squares solution.Thus, in practice, we have a trade-off of SNR gain by using an optimal design versusSNR loss by lower gradient amplitudes or longer TE. Currently, the MD-RME design andits extension to 12 and 24 acquisitions with sign-flips provide the best of both worlds andare straightforward to implement. Specifically, the 12 acquisition sign-flips sequence attainsthe same estimation errors compared to a conventional sequence of 24 acquisitions, withhalf of the acquisition time.We are also looking into additional parametrization of the acquisition design withinour framework. A potential approach is to employ multidirectional fractional encodingratios [173]. In the Appendix, we extend the parameterized phase accumulation for thiscase, in which the MEG frequency is higher than the mechanical frequency, and may varyin each direction and acquisition. This may allow derivation of optimal designs withoutdiscontinuities in the MEGs and an increase of TE. However, this also affects the encodingefficiency and therefore requires changing the gradient amplitudes in (possibly) large factorsto compensate each direction, which may hinder the use of this approach in practice.We limited our experiments to employ SE-EPI sequences. However, other pulse se-quences in the literature, e.g., [44], can be employed to further improve SNR and/or shortenacquisition time. Trapezoidal waveform can also be applied rather than sinusoidal, to max-imize encoding efficiency under limitation of the slew rate. Note that applying multidirec-tional MEGs introduces concomitant-field terms that should be corrected, e.g., by usinggradient symmetrization [14]. In the sequences we use, the MEGs are mirrored after theecho pulse and cancel such effects.We note that modeling the noise in the phase images as white Gaussian noise is justifiedby the fact that the least-squares solution for the conventional sequence using our approachcoincides with the standard DFT processing. Thus, such noise has been implicitly assumedso far in MRE imaging. However, in practice, it might not be the case due to, e.g., change intemperature of the coils between measurements. An advantage of our framework is that wemay employ a different noise model with its proper optimal estimator (e.g., weighted least-squares [18]). In addition, optimality criteria other than D-optimality can be considered forevaluating design performance. Our approach may allow to adapt the acquisition matrix1056.5. Discussion and Conclusionaccordingly by, e.g., using non-uniform sampling of the acquisitions. The flexibility inshaping the matrix may also be important in incorporating compressed sensing in MRE.In conclusion, the formulation of motion encoding as both linear estimation and op-timal design problems allows employing tools from fields such as optimization, statisticsand information theory. A careful design of parameterized multidirectional MEGs yieldsdisplacement maps with a better SNR and/or shorter acquisition times. Along with theongoing improvement in transducers and reconstruction algorithms, accurate elastogramsthat characterize the mechanical properties of in vivo tissue can be derived. In turn, thesemay pave the way for a broader clinical deployment of MRE.106Chapter 7ConclusionIn this thesis we have developed algorithms for registration of the prostate in differentmodalities. Specifically, we considered registration of histopathology of excised prostatespecimens with their volumetric MR and US imaging, and registration of preoperative MRIto intraoperative TRUS imaging of the prostate during RALP. The former registration canbe applied for training classifiers in order to allow an automatic localization of CaP. Thelatter registration can be applied to augment the surgeon’s console with the preoperativeMR image that corresponds to the da Vinci instrument’s tooltip, and, combined with aclassifier, a corresponding CaP probability map.We have also studied the reciprocity between registration and MRE, and showed that byutilizing the properties of the registration and elastography acquisition process, MRE canfacilitate registration, and vice versa. The results of the research presented in the variouschapters complement each other, and are intended to be jointly applied as depicted by thediagram in Figure 7.1 and described below.First, by employing the optimization-based designed motion encoding sequences duringa preoperative MRI scanning session, as proposed in Chapter 6, the quality of the acquiredMRE data is increased, and/or scanning time (and therefore costs and discomfort to thepatient) is reduced. Estimation of the MRE elastograms can be further enhanced or cor-rected by employing the motion compensation and super-resolution technique discussed inChapter 5.In turn, the enhanced elastograms provide accurate biomechanical models of the prostateand periprostatic tissue in preoperative in vivo MRI, which may improve the accuracyof its registration to the corresponding postoperative ex vivo MRI using the deformableregistration method developed in Chapter 3. By composing the resulting deformation mapwith the transformation between the histological slices and ex vivo MRI, as computed usingthe particle filtering algorithm developed in Chapter 2, we obtain an accurate mappingbetween histopathology and in vivo MRI that can train a classifier for localization of CaP.Finally, during future RALP procedures, a reliable classifier will allow a CaP probabilitymap to be overlaid on MR images that correspond to the robot’s tooltip by employing theMRI to TRUS registration algorithm proposed in Chapter 4.1077.1. Contributions and LimitationsChap. 6 Chap. 5 Chap. 3 Chap. 2 Chap. 4In vivo MRI + MRE DataEnhanced MREIn vivo to ex vivoMRI mapTrained CaPclassifierEx vivo MRI Histopathology + in vivo TRUSFigure 7.1: Proposed clinical flow.We note that the data collection that was employed in this thesis was acquired eitherprior or in parallel to development of the methods, and follows protocols that were set andapproved in advance. In retrospect, the multi-parametric MRI settings could have been setto produce images that are, e.g., better correlated with histology, or of higher resolution.Nevertheless, our methods were designed and adapted to contend with the typical dataacquired by practical and clinical protocols.The clinical significance of the research conducted throughout this thesis is that, withMRI correctly displayed on the surgeon’s console and automatic localization of CaP, re-section of the prostate during RALP can be carried out with minimal damage to adjacentstructures. Due to the increasing usage of RALP worldwide, even a slight improvement byusing the proposed methods can impact treatment outcomes of thousands of CaP patients.The proposed methods can also benefit other CaP treatments such as image guided biopsy,brachytherapy, and watchful waiting. In addition, the time reduction in acquisition of MREdata increases its potential in clinical deployment.The main scientific significance lies in the novel registration methods, the usage andstudy of elastography in registration, and the established interface between MRE acquisitionand information theory. Below, we summarize the main contributions that were made inthe course of achieving the thesis objectives, describe the limitations, and discuss possiblefuture research directions.7.1 Contributions and LimitationsRegistration of Whole-Mount Histology and Volumetric Imaging:In Chapter 2, we developed a general framework for multi-slice to volume registration that isapplied to align histopathology with volumetric MR and US imaging of the prostate. Theregistration is formulated as a filtering problem, which allows prior knowledge regardingthe physical misalignment mechanisms to be incorporated into the framework in a Bayesianfashion. We also implemented an accompanying visualization tool allows comparison of1087.1. Contributions and Limitationseach histological slice side-by-side with corresponding slices from two other modalities, andyield the intensity histograms of these slices on user selected points, lines and regions.The registration of histological slices to ex vivo MR, in vivo MR, and in vivo US achievedrelatively high area overlaps and registration errors below the diameter of typical lesions thatare marked during histopathology analysis (∼10± 5 mm [34]). The registration is thereforebeing utilized in ongoing studies to optimize and evaluate CaP detection on a range ofimaging modalities. These studies have yielded several publications ([95, 98, 137, 138]) thatcharacterize cancer on US, TRUS-VE, and on both in vivo and ex vivo MRI and MRE.The main limitation of the proposed method, which stems from its generality, is therelative high number of parameters that need to be set. These parameters can be chosen bya trial-and-error process on a few training cases. In addition, the registration time may belong when using intensity- or region-based metrics, whereas the point-based metric requiresan additional segmentation of the volumetric images. Nevertheless, longer running timesare tolerable in histology processing as the algorithm is running offline, and segmentationerrors can be taken into account by modeling of the noise in the filtering formulation.Model-Based Registration Using Elastography:In Chapter 3, we developed and implemented an algorithm for deformable registration ofex vivo and in vivo MRI. By allowing a two-step registration with the ex vivo MR asan intermediate modality, this algorithm can be integrated with the multi-slice to volumealignment above in order to capture the residual deformations between histological slicesand in vivo MR.The method presents the first utilization of measured in vivo elasticity by MRE for im-proving deformable registration. The MRE elastogram is incorporated into the registrationframework as a biomechanical model that regularizes the deformation according to physicalinhomogeneous elasticity. By studying the effect of inhomogeneity in elasticity on registra-tion errors, we have showed that this approach yields improved registration as opposed toapplying a constant elasticity value to the volume, and that, in general, the performanceimproves with the inhomogeneity of elasticity inside and outside of the prostate. We havealso proposed a novel similarity metric that does not require segmentation and meshing ofthe in vivo anatomy, as opposed to explicit FEM-based modeling approaches.As in other studies that concern deformable registration methods, the ground-truthdeformation map is unknown and the results are being evaluated on a sparse set of landmarksand volume overlaps. Another limitation of our study is that the choice of landmarks andsegmentations were performed once by one expert, thus we do not have intra- or inter-1097.1. Contributions and Limitationsobserver variability measurements. Such information would have allowed evaluation ofwhether the lower accuracy for base and apex stems from the registration itself or fromsegmentation inaccuracy. A disadvantage of the proposed approach is that the ex vivoMRI scan is time consuming and expensive. In addition, current in vivo MRE protocolsalso involve longer scanning time, inconvenience to the patient, and are still not widelyaccessible. This limitation may be mitigated by the results of Chapter 6.Registration of Ultrasound and Magnetic Resonance Imaging for ImageGuided Interventions:In Chapter 4, we developed a feasible method for introducing MRI data in the operatingroom, by registering preoperative MRI to intraoperative TRUS. The proposed approachcan employ existing building blocks to enable tracking and display of the MR image thatcorresponds to the da Vinci instrument’s tooltip during RALP. The segmentation-basedapproach deforms the segmented T2w volume to match its surface with that of an intraop-eratively segmented B-mode. The registration is fast and robust, while its performance wasfound to be plausible, with results comparable to other TRUS-MR registration methods.A limitation of the proposed approach is that, as in other segmentation-based methods,the quality of the registration depends on the quality of the segmentation. Another limi-tation is that the registration is only valid with respect to the time the TRUS volume wasacquired, thus rendering its usage during later stages of the surgery challenging.We have also studied the usage of elastography as an inter-modality for registeringTRUS and MRI. While we found the elasticity values to be nondescriptive for a deformableintensity-based registration, the phantom and four clinical data sets we used are statisticallyinsufficient to conclude that. Typically, in studies that evaluate new registration methods,data sets in the range of 10-20 patients are being used. However, the required sample sizefor drawing a statistically significant conclusion in our case should be determined by settingthe desired statistical power and significance level, and using power-analysis tools [77].Motion Compensation and Super-Resolution in MRE:In Chapter 5, we turned to employ registration for enhancing MRE data. Magnitudeimages, which are often regarded as a by-product of the MRE acquisition process andignored, are utilized to compensate for motion of patients during the, typically lengthy,acquisition process. The inherent alignment of the magnitude and phase images allows forthe enhancement of elastogram reconstruction by employing a super-resolution techniquethat increases the resolution of the corresponding registered phase images.1107.1. Contributions and LimitationsThe main contribution in this work is the improvement of existing MRE data by usingimage processing techniques without altering the current data acquisition protocol. This ispossible due to observation and utilization of the redundancy in the data acquired in MRE.We noticed a better separability of objects from their background on the reconstructedhigher resolution elastograms in comparison with the native resolution ones.We note, however, that assuming a Gaussian blurring kernel in the super-resolutiontechnique is inaccurate, and a model of the actual point spread function of the scanneris required in order to yield accurate reconstructions of the high resolution images. Inaddition, [35] suggested that the highest resolution increase in each axis using two images(each pair of phase images in our case) is by√2. Thus, doubling the resolution based on [61]may introduce further inaccuracies to the super-resolution phase images and elastograms.We also note that other MRI sequences that encode motion were not considered. Specif-ically, the harmonic phase (HARP) algorithm [110] is popular for detection and trackingcardiac motion using tagged MRI. However, tagged MRI sequences, which employ RFto encode motion, are slower than MEGs and typically do not allow capturing motion infrequencies high enough for dynamic elastography.Optimization-Based Design of Motion Encoding in MRE:In Chapter 6, we proposed a general framework for improving MRE acquisition by opti-mizing designs of the motion encoding process. The framework lays the groundwork forformulating the MRE acquisition problem within fields such as optimization, statistics andinformation theory, and thus allows employing established tools from these fields for quan-tification of design performance and derivation of optimal sequences by different criteria.Generalization and parametrization of the MEGs were introduced, and by using theproposed framework we derived and implemented a multidirectional sequence with a min-imum number of acquisitions that does not require an increase in TE, nor discontinuitiesin the MEGs. The reduced acquisition time and improved data quality that is achieved bythe proposed approach may assist in bringing MRE to broader use.While multidirectional MEGs allow minimization of the variance in motion estimationerror, estimations among the different directions become correlated. Thus, if the scannerhas inhomogeneity along one direction, the motion estimations in all directions are affected,rather than just the corresponding direction. Physical limitations such as concomitant-fieldterms should be taken into account when employing multidirectional gradients. In addition,hardware issues or overheating of the coils may limit the maximum gradient amplitude orrequire cooling down periods between acquisitions of multidirectional sequences.1117.2. Future Work7.2 Future WorkBased on the results of this thesis, the main future work is to conduct a clinical study forevaluating the performance of the proposed MR-TRUS registration intraoperatively, andperform a cost-benefit analysis of using this approach during RALP, based on evidence andfeedback from clinicians. In addition, the two-step registration between histological slicesand in vivo MR that employs the ex vivo MR as an intermediate modality is to be usedfor training classifiers of CaP. Other possible future research directions that may furtherimprove or extend the proposed methods are summarized below.Registration of Whole-Mount Histology and Volumetric Imaging:• Incorporation of a nonrigid registration step in order to improve registration of his-tology and in vivo imaging. Due to the sparsity of slices, a 2-D registration shouldbe applied between each histological slice and its corresponding registered slice on thevolumetric image.• Investigation of other metrics that can describe the similarity between histologicalslices and volumetric imaging. Specifically, employing a point-based metric betweendescriptors of intensity-based extracted features.Model-Based Registration Using Elastography:• Derivation of forces through intensity-based similarity metrics between in vivo andex vivo T2w MRI to eliminate the required segmentation of the ex vivo volume.• Consider nonlinear elasticity to contend with large deformations of the prostate, e.g.,in case of employing a transrectal coil during image acquisition.Registration of Ultrasound and Magnetic Resonance Imaging for ImageGuided Interventions:• Allow a manual selection of corresponding points on MR and TRUS images that areused as anchor points to constrain the registration. This will provide the user withsome control over the registration, and the ability to correct it in case the automaticresults are not plausible.• Allow switching the display between different MR images, e.g., diffusion weightedimaging or MRE, which were acquired during the multi-parametric preoperative MRIsession, and can be co-registered to the T2w volume (and thus to the TRUS volume).1127.2. Future Work• Optimization of the code and its implementation on C/C++ to reduce running time,thus making the proposed approach more feasible in a clinical setting.• Develop a tracking algorithm that allows real-time deformation of the registered MRaccording to the TRUS volume.• Study other metrics and the extraction of features that may be better suitable forintensity-based registration of MRE and USE.Motion Compensation and Super-Resolution in MRE:• Iterative optimization of motion compensation and elasticity reconstruction, by ap-plying the reconstructed elastogram in the current iteration as regularization for themotion compensation displacements in the next iteration.• Incorporation of multi-frequency information for adding additional data to the super-resolution technique and improving displacement estimation.• Considering blurring kernels in the super-resolution technique that model the pointspread function of the scanner in order to improve the accuracy of the high resolutionimage.Optimization-Based Design of Motion Encoding in MRE:• Additional parametrization of the MEGs, e.g., including varying fractional encodingratios in each direction, in order to increase the degrees of freedom in designs withinour framework.• Employment of the proposed framework to extend existing pulse sequences, e.g., [44],into multidirectional in order to maximize encoding efficiency under limitation of theslew rate.• Employment of different noise models with corresponding optimal estimators.• Employment of different optimality criteria for evaluating design performance.113Bibliography[1] P. Abolmaesumi and M. R. Sirouspour. An interacting multiple model probabilisticdata association filter for cavity boundary extraction from ultrasound images. IEEETransactions on Medical Imaging, 23(6):772–784, 2004.[2] T. K. Adebar, M. C. Yip, S. E. Salcudean, R. N. Rohling, C. Y. Nguan, and S. L.Goldenberg. Registration of 3D ultrasound through an air–tissue boundary. IEEETransactions on Medical Imaging, 31(11):2133–2142, 2012.[3] T. E. Ahlering, D. Woo, L. Eichel, D. I. Lee, R. Edwards, and D. W. Skarecky. Robot-assisted versus open radical prostatectomy: a comparison of one surgeon’s outcomes.Urology, 63(5):819–822, 2004.[4] N. Alpert, J. Bradshaw, D. Kennedy, and J. Correia. The principal axes transforma-tion: a method for image registration. Journal of Nuclear Medicine, 31(10):1717–1722,1990.[5] R. Alterovitz, K. Goldberg, J. Pouliot, I. Hsu, Y. Kim, S. Noworolski, andJ. Kurhanewicz. Registration of MR prostate images with biomechanical modelingand nonlinear parameter estimation. Medical Physics, 33:446, 2006.[6] American Cancer Society. Cancer facts & figures 2014. http://www.cancer.org/(retrieved on April 18, 2014).[7] E. R. Arce-Santana, D. U. Campos-Delgado, and A. Alba. Affine image registrationguided by particle filter. IET Image Processing, 6(5):455–462, 2012.[8] M. S. Arulampalam and et al. A tutorial on particle filters for online nonlinear/non-gaussian Bayesian tracking. IEEE Transactions on Signal Processing, 50(2):174–188,2002.[9] A. Baghani, A. Brant, S. Salcudean, and R. Rohling. A high-frame-rate ultrasoundsystem for the study of tissue motions. IEEE Transactions on Ultrasonics, Ferro-electrics and Frequency Control, 57(7):1535–1547, 2010.114[10] A. Baghani, H. Eskandari, W. Wang, D. Da Costa, M. N. Lathiff, R. Sahebjavaher,S. Salcudean, and R. Rohling. Real-time quantitative elasticity imaging of deep tissueusing free-hand conventional ultrasound. Medical Image Computing and ComputerAssisted Intervention (MICCAI), pages 617–624, 2012.[11] R. Bajcsy and S. Kovacic. Multiresolution elastic matching. Computer Vision Graph-ics and Image Processing, 46:1–21, 1989.[12] S. Bart, P. Mozer, P. Hemar, G. Lenaour, E. Comperat, R. Renard-Penna, E. Chartier-Kastler, and J. Troccaz. MRI-histology registration in prostate cancer. In Surgetica,pages 361–367, 2005.[13] M. Baumann, P. Mozer, V. Daanen, and J. Troccaz. Prostate biopsy tracking withdeformation estimation. Medical Image Analysis, 2011.[14] M. A. Bernstein, X. J. Zhou, J. A. Polzin, K. F. King, A. Ganin, N. J. Pelc, and G. H.Glover. Concomitant gradient terms in phase contrast MR: analysis and correction.Magnetic Resonance in Medicine, 39(2):300–308, 1998.[15] P. J. Besl and N. D. McKay. A method for registration of 3-D shapes. IEEE Trans-actions on Pattern Analysis and Machine Intelligence, 14(2):239–256, 1992.[16] A. Bharatha et al. Evaluation of three-dimensional finite element-based deformableregistration of pre-and intraoperative prostate imaging. Medical Physics, 28:2551,2001.[17] M. Bilgen and M. Insana. Predicting target detectability in acoustic elastography.IEEE Ultrasonics Symposium, 2:1427–1430, 1997.[18] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press,Cambridge UK, 2004.[19] L. Breiman. Random forests. Machine learning, 45(1):5–32, 2001.[20] C. Broit. Optimal registration of deformed images. PhD thesis, Computer and Infor-mation Science Dept., University of Pennsylvania, Philadelphia, PA, 1981.[21] T. Brox and D. Cremers. On local region models and a statistical interpretation ofthe piecewise smooth mumford-shah functional. International Journal of ComputerVision, 84(2):184–193, 2009.115[22] V. Caselles, R. Kimmel, and G. Sapiro. Geodesic active contours. InternationalJournal of Computer Vision, 22(1):61–79, 1997.[23] T. F. Chan and L. A. Vese. Active contours without edges. IEEE Transactions onImage Processing, 10(2):266–277, 2001.[24] J. Chappelow, B. Bloch, N. Rofsky, E. Genega, R. Lenkinski, W. DeWolf, and A. Mad-abhushi. Elastic registration of multimodal prostate MRI and histology via multiat-tribute combined mutual information. Medical Physics, 38(4):2005, 2011.[25] G. E. Christensen, R. D. Rabbit, and M. I. Miller. Deformable templates using largedeformation kinematics. IEEE Transactions on Medical Imaging, 5:1435–1447, 1996.[26] A. Collignon, F. Maes, D. Delaere, D. Vandermeulen, P. Suetens, and G. Marchal. Au-tomated multi-modality image registration based on information theory. In Y. Bizais,C. Barillot, and R. Di Paola, editors, Information Processing in Medical Imaging,pages 263–274. Kluwer Academic Publishers, Dordrecht Netherlands, 1995.[27] W. A. Coppel. Hadamard’s determinant problem. In Number Theory, pages 223–259.Springer, 2009.[28] V. Daanen, J. Gastaldo, J.-Y. Giraud, P. Fourneret, J.-L. Descotes, M. Bolla, D. Col-lomb, and J. Troccaz. MRI/TRUS data fusion for brachytherapy. The InternationalJournal of Medical Robotics and Computer Assisted Surgery, 2(3):256–261, 2006.[29] C. Davatzikos. Spatial transformation and registration of brain images using elas-tically deformable models. Computer Vision and Image Understanding, 66:207–222,1997.[30] A. Doucet, N. de Freitas, and N. Gordon. Sequential Monte Carlo Methods in Practice.Springer, New York, 2001.[31] A. Doucet and A. Johansen. A tutorial on particle filtering and smoothing: fifteenyears later. Handbook of Nonlinear Filtering, pages 656–704, 2009.[32] B. Drew, E. Jones, S. Reinsberg, A. Yung, S. Goldenberg, and P. Kozlowski. Devicefor sectioning prostatectomy specimens to facilitate comparison between histology andin vivo MRI. Journal of Magnetic Resonance Imaging, 32(4):992–996, 2010.[33] R. O. Duda, P. E. Hart, and D. G. Stork. Pattern Classification. Wiley-Intersciance,New York, second edition, 2001.116[34] L. Eichelberger, M. Koch, J. Daggy, T. Ulbright, J. Eble, and L. Cheng. Predictingtumor volume in radical prostatectomy specimens from patients with prostate cancer.American Journal of Clinical Pathology, 120(3):386–391, 2003.[35] M. Elad and A. Feuer. Restoration of a single superresolution image from severalblurred, noisy, and undersampled measured images. IEEE Transactions on ImageProcessing, 6(12):1646–1658, 1997.[36] J. I. Epstein, P. C. Walsh, M. Carmichael, and C. B. Brendler. Pathologic andclinical findings to predict tumor extent of nonpalpable (stage t1 c) prostate cancer.The Journal of the American Medical Association (JAMA), 271(5):368–374, 1994.[37] H. Eskandari, O. Goksel, S. E. Salcudean, and R. Rohling. Bandpass sampling ofhigh-frequency tissue motion. IEEE Transactions on Ultrasonics, Ferroelectrics andFrequency Control, 58(7):1332–1343, 2011.[38] M. Ferrant, A. Nabavi, B. Macq, F. Jolesz, R. Kikinis, and S. Warfield. Registrationof 3-D intraoperative MR images of the brain using a finite-element biomechanicalmodel. IEEE Transactions on Medical Imaging, 20(12):1384–1397, 2001.[39] V. Ficarra, G. Novara, W. Artibani, A. Cestari, A. Galfano, M. Graefen, G. Guaz-zoni, B. Guillonneau, M. Menon, F. Montorsi, et al. Retropubic, laparoscopic, androbot-assisted radical prostatectomy: a systematic review and cumulative analysis ofcomparative studies. European Urology, 55(5):1037–1063, 2009.[40] B. Fischer and J. Modersitzki. Fast inversion of matrices arising in image processing.Numerical Algorithms, 22(1):1–11, 1999.[41] B. Fischer and J. Modersitzki. A unified approach to fast image registration anda new curvature based registration technique. Linear Algebra and its Applications,380:107–124, 2004.[42] Z. Galil and J. Kiefer. D-optimum weighing designs. The Annals of Statistics, pages1293–1306, 1980.[43] Y. Gao, R. Sandhu, G. Fichtinger, and A. Tannenbaum. A coupled global registrationand segmentation framework with application to magnetic resonance prostate imagery.IEEE Transactions on Medical Imaging, 29(10):1781–1794, 2010.117[44] P. Garteiser, R. S. Sahebjavaher, L. C. Ter Beek, S. Salcudean, V. Vilgrain, B. E.Van Beers, and R. Sinkus. Rapid acquisition of multifrequency, multislice and multi-directional MR elastography data with a fractionally encoded gradient echo sequence.NMR in Biomedicine, 26(10):1326–1335, 2013.[45] E. Gibson, C. Crukley, M. Gaed, J. A. Go´mez, M. Moussa, J. L. Chin, G. S. Bauman,A. Fenster, and A. D. Ward. Registration of prostate histology images to ex vivoMR images via strand-shaped fiducials. Journal of Magnetic Resonance Imaging,36(6):1402–1412, 2012.[46] K. J. Glaser, A. Manduca, and R. L. Ehman. Review of MR elastography applicationsand recent developments. Journal of Magnetic Resonance Imaging, 36(4):757–774,2012.[47] N. J. Gordon, D. J. Salmond, and A. F. Smith. Novel approach to nonlinear/non-gaussian bayesian state estimation. In IEE Proceedings F (Radar and Signal Process-ing), volume 140, pages 107–113. IET, 1993.[48] H. Gray, S. Standring, H. Ellis, and B. Berkovitz. Gray’s anatomy: the anatomicalbasis of clinical practice. Gray’s Anatomy: The Anatomical Basis of Clinical Practice.Elsevier Churchill Livingstone, 2005.[49] H. Greenspan. Super-resolution in medical imaging. The Computer Journal, 52(1):43–63, 2009.[50] J. Hadamard. Re´solution d’une question relative aux de´terminants. Bulletin desSciences Mathe´matiques, 17(1):240–246, 1893.[51] B. A. Hadaschik, T. H. Kuru, C. Tulea, P. Rieker, I. V. Popeneciu, T. Simpfendo¨rfer,J. Huber, P. Zogal, D. Teber, S. Pahernik, et al. A novel stereotactic prostate biopsysystem integrating pre-interventional magnetic resonance imaging and live ultrasoundfusion. The Journal of Urology, 186(6):2214–2220, 2011.[52] S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent. Optimal mass transport forregistration and warping. International Journal of Computer Vision, 60(3):225–240,2004.[53] R. Hartley and A. Zisserman. Multiple View Geometry in Computer Vision. Cam-bridge University Press, 2nd edition, 2004.118[54] R. Hausmann, J. S. Lewin, and G. Laub. Phase-contrast MR angiography withreduced acquisition time: New concepts in sequence design. Journal of MagneticResonance Imaging, 1(4):415–422, 1991.[55] M. P. Heinrich, M. Jenkinson, M. Bhushan, T. Matin, F. V. Gleeson, S. M. Brady, andJ. A. Schnabel. MIND: Modality independent neighbourhood descriptor for multi-modal deformable registration. Medical Image Analysis, 16(7):1423–1435, 2012.[56] D. Herzka, M. Kotys, R. Sinkus, R. Pettigrew, and A. Gharib. Magnetic resonanceelastography in the liver at 3 Tesla using a second harmonic approach. MagneticResonance in Medicine, 62(2):284–291, 2009.[57] M. Honarvar, R. S. Sahebjavaher, S. E. Salcudean, and R. Rohling. Sparsity reg-ularization in dynamic elastography. Physics in Medicine and Biology, 57(19):5909,2012.[58] B. K. P. Horn and B. G. Schunck. Determining optical flow. Artificial Intelligence,17:185–203, 1981.[59] Y. Hu, H. U. Ahmed, Z. Taylor, C. Allen, M. Emberton, D. Hawkes, and D. Barratt.MR to ultrasound registration for image-guided prostate interventions. Medical ImageAnalysis, 16(3):687–703, 2012.[60] C. Hughes, O. Rouviere, F. Mege-Lechevallier, R. Souchon, and R. Prost. Robustalignment of prostate histology slices with quantified accuracy. IEEE Transactionson Biomedical Engineering, 60(2):281–291, 2013.[61] M. Irani and S. Peleg. Improving resolution by image registration. CVGIP: GraphicalModels and Image Processing, 53(3):231–239, 1991.[62] M. Ishida, A. Schuster, S. Takase, G. Morton, A. Chiribiri, B. Bigalke, T. Schaeffter,H. Sakuma, and E. Nagel. Impact of an abdominal belt on breathing patterns andscan efficiency in whole-heart coronary magnetic resonance angiography: comparisonbetween the uk and japan. Journal of Cardiovascular Magnetic Resonance, 13:71,2011.[63] R. S. John and N. R. Draper. D-optimality for regression designs: a review. Techno-metrics, 17(1):15–23, 1975.[64] S. Kabus, A. Franz, and B. Fischer. Spatially varying elasticity in image registration.Methods of Information in Medicine, 46(3):287–291, 2007.119[65] T. Kailath, A. H. Sayed, and B. Hassibi. Linear estimation, volume 1. Prentice HallNew Jersey, 2000.[66] M. Kass, A. Witkin, and D. Terzopoulos. Snakes: Active contour models. Interna-tional Journal of Computer Vision, 1(4):321–331, 1988.[67] S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi. Gradient flowsand geometric active contour models. International Conference on Computer Vision,pages 810–815, 1995.[68] J. Kiefer. Optimum designs in regression problems, ii. The Annals of MathematicalStatistics, pages 298–325, 1961.[69] J. Kiefer and J. Wolfowitz. Optimum designs in regression problems. The Annals ofMathematical Statistics, pages 271–294, 1959.[70] S. Kimm, T. Tarin, J. Lee, B. Hu, K. Jensen, D. Nishimura, and J. Brooks. Methodsfor registration of magnetic resonance images of ex vivo prostate specimens withhistology. Journal of Magnetic Resonance Imaging, 2012.[71] A. Kirkham, M. Emberton, and C. Allen. How good is MRI at detecting and charac-terising cancer within the prostate? European Urology, 50(6):1163–1175, 2006.[72] D. Klatt, T. K. Yasar, T. J. Royston, and R. L. Magin. Sample interval modulationfor the simultaneous acquisition of displacement vector data in magnetic resonanceelastography: theory and application. Physics in Medicine and Biology, 58(24):8663–8675, 2013.[73] I. Kolesov, J. Lee, P. Vela, and A. Tannenbaum. A stochastic approach for non-rigidimage registration. SPIE Medical Imaging, 86550:U1–U15, 2013.[74] K. Konig, U. Scheipers, A. Pesavento, A. Lorenz, H. Ermert, and T. Senge. Initialexperiences with real-time elastography guided biopsies of the prostate. The Journalof Urology, 174(1):115–117, 2005.[75] P. Kozlowski et al. Combined prostate diffusion tensor imaging and dynamic con-trast enhanced MRI at 3T–quantitative correlation with biopsy. Magnetic ResonanceImaging, 28(5):621–628, 2010.[76] C. Le Guyader and L. Vese. A combined segmentation and registration frameworkwith a nonlinear elasticity smoother. Computer Vision and Image Understanding,2011.120[77] R. V. Lenth. Some practical guidelines for effective sample size determination. TheAmerican Statistician, 55(3):187–193, 2001.[78] H. Lester, S. Arridge, K. Jansons, L. Lemieux, J. Hajnal, and A. Oatridge. Non-linearregistration with the variable viscosity fluid algorithm. In Information Processing inMedical Imaging, pages 238–251. Springer, 1999.[79] J. Lobo. Towards intra-operative dosimetry for prostate brachytherapy: improved seeddetection and registration to ultrasound using needle detection. PhD thesis, Electricaland Computer Engineering, University of British Columbia, Vancouver, BC, 2012.[80] J.-A. Long, B. H. Lee, J. Guillotreau, R. Autorino, H. Laydner, R. Yakoubi,E. Rizkala, R. J. Stein, J. H. Kaouk, and G.-P. Haber. Real-time robotic tran-srectal ultrasound navigation during robotic radical prostatectomy: initial clinicalexperience. Urology, 80(3):608–613, 2012.[81] Y. Lotan, J. A. Cadeddu, and M. T. Gettman. The new economics of radical prosta-tectomy: cost comparison of open, laparoscopic and robot assisted techniques. TheJournal of Urology, 172(4):1431–1435, 2004.[82] B. Ma and R. Ellis. Surface-based registration with a particle filter. Medical ImageComputing and Computer Assisted Intervention (MICCAI), pages 566–573, 2004.[83] S. S. Mahdavi, M. Moradi, W. J. Morris, and S. E. Salcudean. Automatic prostatesegmentation using fused ultrasound B-mode and elastography images. Medical ImageComputing and Computer Assisted Intervention (MICCAI), pages 76–83, 2010.[84] J. Maintz and M. Viergever. A survey of medical image registration. Medical ImageAnalysis, 2(1):1–36, 1998.[85] A. Manduca, T. E. Oliphant, M. Dresner, J. Mahowald, S. Kruse, E. Amromin, J. P.Felmlee, J. F. Greenleaf, and R. L. Ehman. Magnetic resonance elastography: non-invasive mapping of tissue elasticity. Medical Image Analysis, 5(4):237–254, 2001.[86] L. Marks, S. Young, and S. Natarajan. MRI–ultrasound fusion for guidance of targetedprostate biopsy. Current opinion in urology, 23(1):43, 2013.[87] S. Martin, M. Baumann, V. Daanen, and J. Troccaz. MR prior based automaticsegmentation of the prostate in TRUS images for MR/TRUS data fusion. In IEEEInternational Symposium on Biomedical Imaging, pages 640–643. IEEE, 2010.121[88] A. Mejia-Rodriguez, E. Arce-Santana, E. Scalco, D. Tresoldi, M. Mendez, A. Bianchi,G. Cattaneo, and G. Rizzo. Elastic registration based on particle filter in radiotherapyimages with brain deformations. In Engineering in Medicine and Biology Society(EMBC), pages 8049–8052. IEEE, 2011.[89] M. Menon, A. Tewari, B. Baize, B. Guillonneau, and G. Vallancien. Prospectivecomparison of radical retropubic prostatectomy and robot-assisted anatomic prosta-tectomy: the vattikuti urology institute experience. Urology, 60(5):864–868, 2002.[90] J. Mitra, Z. Kato, R. Mart´ı, A. Oliver, X. Llado´, D. Sidibe´, S. Ghose, J. C. Vi-lanova, J. Comet, and F. Meriaudeau. A spline-based non-linear diffeomorphism formultimodal prostate registration. Medical Image Analysis, 16(6):1259–1279, 2012.[91] T. Miyagawa, S. Ishikawa, T. Kimura, T. Suetomi, M. Tsutsumi, T. Irie, M. Kondoh,and T. Mitake. Real-time virtual sonography for navigation during targeted prostatebiopsy using magnetic resonance imaging data. International Journal of Urology,17(10):855–860, 2010.[92] J. Modersitzki. Numerical methods for image registration. Oxford University Press,USA, 2004.[93] H. Moghari, M. and P. Abolmaesumi. Point-based rigid-body registration using anunscented Kalman filter. IEEE Transactions on Image Processing, 26(12):1708–1728,Dec. 2007.[94] O. Mohareri, M. Ramezani, T. K. Adebar, P. Abolmaesumi, and S. E. Salcudean.Automatic localization of the da vinci surgical instrument tips in 3-d transrectalultrasound. IEEE Transactions on Biomedical Engineering, 60(9):2663–2672, 2013.[95] O. Mohareri, A. Ruszkowski, J. Lobo, J. Ischia, A. Baghani, G. Nir, H. Eskandari,E. Jones, L. Fazli, L. Goldenberg, et al. Multi-parametric 3D quantitative ultrasoundvibro-elastography imaging for detecting palpable prostate tumors. Medical ImageComputing and Computer Assisted Intervention (MICCAI), pages 561–568, 2014.[96] A. M. Mood. On hotelling’s weighing problem. The Annals of Mathematical Statistics,17(4):432–446, 1946.[97] C. M. Moore, N. L. Robertson, N. Arsanious, T. Middleton, A. Villers, L. Klotz, S. S.Taneja, and M. Emberton. Image-guided prostate biopsy using magnetic resonanceimaging–derived targets: a systematic review. European Urology, 63(1):125–140, 2013.122[98] M. Moradi, S. S. Mahdavi, G. Nir, O. Mohareri, A. Koupparis, L.-O. Gagnon, L. Fazli,R. G. Casey, J. Ischia, E. C. Jones, et al. Multiparametric 3d in vivo ultrasoundvibroelastography imaging of prostate cancer: Preliminary results. Medical Physics,41(7):073505, 2014.[99] P. S. Moran, M. O’Neill, C. Teljeur, M. Flattery, L. A. Murphy, G. Smyth, andM. Ryan. Robot-assisted radical prostatectomy compared with open and laparoscopicapproaches: A systematic review and meta-analysis. International Journal of Urology,20(3):312–321, 2013.[100] R. Muthupillai, D. Lomas, P. Rossman, J. Greenleaf, A. Manduca, and R. Ehman.Magnetic resonance elastography by direct visualization of propagating acoustic strainwaves. Science, 269(5232):1854–1857, 1995.[101] R. Narayanan, J. Kurhanewicz, K. Shinohara, E. Crawford, A. Simoneau, and J. Suri.MRI-ultrasound registration for targeted prostate biopsy. IEEE International Sym-posium on Biomedical Imaging, pages 991–994, 2009.[102] G. Nir, R. S. Sahebjavaher, A. Baghani, R. Sinkus, and S. E. Salcudean. Prostatesegmentation in MRI using fused T2-weighted and elastography images. SPIE MedicalImaging, 90340:C1–C6, 2014.[103] G. Nir, R. S. Sahebjavaher, P. Kozlowski, S. D. Chang, R. Sinkus, S. L. Goldenberg,and S. E. Salcudean. Model-based registration of ex vivo and in vivo MRI of theprostate using elastography. IEEE Transactions on Medical Imaging, 32(7):1349–1361, 2013.[104] G. Nir, R. S. Sahebjavaher, R. Sinkus, and S. E. Salcudean. A framework foroptimization-based design of motion encoding in magnetic resonance elastography.Magnetic Resonance in Medicine, 2014.[105] G. Nir and A. Tannenbaum. Temporal registration of partial data using particlefiltering. In IEEE International Conference on Image Processing, pages 2177–2180.IEEE, 2011.[106] K. Noe and T. Sørensen. Solid mesh registration for radiotherapy treatment planning.Biomedical Simulation, pages 59–70, 2010.[107] T. E. Oliphant, A. Manduca, R. L. Ehman, and J. F. Greenleaf. Complex-valuedstiffness reconstruction for magnetic resonance elastography by algebraic inversion ofthe differential equation. Magnetic Resonance in Medicine, 45(2):299–310, 2001.123[108] J. Ophir, I. Ce´spedes, H. Ponnekanti, Y. Yazdi, and X. Li. Elastography: a quan-titative method for imaging the elasticity of biological tissues. Ultrasonic Imaging,13(2):111–134, 1991.[109] S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed:algorithms based on hamilton-jacobi formulations. Journal of Computational Physics,79(1):12–49, 1988.[110] N. F. Osman, E. R. McVeigh, and J. L. Prince. Imaging heart motion using harmonicphase MRI. IEEE Transactions on Medical Imaging, 19(3):186–202, 2000.[111] Y. Ou, D. Shen, M. Feldman, J. Tomaszewski, and C. Davatzikos. Non-rigid reg-istration between histological and MR images of the prostate: a joint segmentationand registration framework. In IEEE Conference on Computer Vision and PatternRecognition Workshops, pages 125–132. IEEE, 2009.[112] S. Papazoglou, J. Rump, J. Braun, and I. Sack. Shear wave group velocity inversionin MR elastography of human skeletal muscle. Magnetic Resonance in Medicine,56(3):489–497, 2006.[113] S. C. Park, M. K. Park, and M. G. Kang. Super-resolution image reconstruction: atechnical overview. Signal Processing Magazine, IEEE, 20(3):21–36, 2003.[114] A. W. Partin, M. W. Kattan, E. N. Subong, P. C. Walsh, K. J. Wojno, J. E. Oester-ling, P. T. Scardino, and J. Pearson. Combination of prostate-specific antigen, clin-ical stage, and gleason score to predict pathological stage of localized prostate can-cer: a multi-institutional update. The Journal of the American Medical Association(JAMA), 277(18):1445–1451, 1997.[115] V. R. Patel, R. Thaly, and K. Shah. Robotic radical prostatectomy: outcomes of 500cases. BJU international, 99(5):1109–1112, 2007.[116] N. J. Pelc, M. A. Bernstein, A. Shimakawa, and G. H. Glover. Encoding strategies forthree-direction phase-contrast MR imaging of flow. Journal of Magnetic ResonanceImaging, 1(4):405–413, 1991.[117] X. Pennec and J. Thirion. A framework for uncertainty and validation of 3-d registra-tion methods based on points and frames. International Journal of Computer Vision,25(3):203–229, 1997.124[118] A. Pesavento and A. Lorenz. Real time strain imaging and in-vivo applications inprostate cancer. IEEE Ultrasonics Symposium, 2:1647–1652, 2001.[119] C. Phillips. Tracking the rise of robotic surgery for prostate cancer. NCI CancerBulletin, 8(16), 2011.[120] P. A. Pinto, P. H. Chung, A. R. Rastinehad, A. A. Baccala Jr, J. Kruecker, C. J.Benjamin, S. Xu, P. Yan, S. Kadoury, C. Chua, et al. Magnetic resonance imag-ing/ultrasound fusion guided prostate biopsy improves cancer detection followingtransrectal ultrasound biopsy and correlates with multiparametric magnetic resonanceimaging. The Journal of Urology, 186(4):1281–1285, 2011.[121] J. G. Pipe et al. Motion correction with PROPELLER MRI: application to head mo-tion and free-breathing cardiac imaging. Magnetic Resonance in Medicine, 42(5):963–969, 1999.[122] J. P. W. Pluim, J. B. A. Maintz, and M. A. Viergever. Mutual information basedregistration of medical images: a survey. IEEE Transactions on Medical Imaging,22:986–1004, 2003.[123] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. NumericalRecipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press,New York, NY, USA, 3 edition, 2007.[124] F. Pukelsheim. Optimal design of experiments, volume 50. SIAM, 1993.[125] Y. Rathi, N. Vaswani, and A. Tannenbaum. A generic framework for tracking usingparticle filter with dynamic shape prior. IEEE Transactions on Image Processing,16(5):1370–1382, 2007.[126] C. Reynier, J. Troccaz, P. Fourneret, A. Dusserre, C. Gay-Jeune, J.-L. Descotes,M. Bolla, and J.-Y. Giraud. MRI/TRUS data fusion for prostate brachytherapy.preliminary results. Medical Physics, 31(6):1568–1575, 2004.[127] M. D. Robson, J. C. Gore, and R. T. Constable. Measurement of the point spreadfunction in MRI using constant time imaging. Magnetic Resonance in Medicine,38(5):733–740, 1997.[128] K. Rohr, H. Stiehl, R. Sprengel, T. Buzug, J. Weese, and M. Kuhn. Landmark-basedelastic registration using approximating thin-plate splines. IEEE Transactions onMedical Imaging, 20(6):526–534, 2001.125[129] O. Rouvie`re, C. Reynolds, T. Hulshizer, P. Rossman, Y. Le, J. P. Felmlee, and R. L.Ehman. MR histological correlation: a method for cutting specimens along the imag-ing plane in animal or ex vivo experiments. Journal of Magnetic Resonance Imaging,23(1):60–69, 2006.[130] D. Rueckert, L. Sonoda, C. Hayes, D. Hill, M. Leach, and D. Hawkes. Nonrigidregistration using free-form deformations: application to breast MR images. IEEETransactions on Medical Imaging, 18(8):712–721, 1999.[131] J. Rump, D. Klatt, J. Braun, C. Warmuth, and I. Sack. Fractional encoding ofharmonic motions in MR elastography. Magnetic Resonance in Medicine, 57(2):388–395, 2007.[132] L. Ruthotto, S. Mohammadi, and N. Weiskopf. A new method for joint susceptibilityartefact correction and super-resolution for dMRI. SPIE Medical Imaging, 90340:P1–P4, 2014.[133] I. Sack, C. K. Mcgowan, A. Samani, C. Luginbuhl, W. Oakden, and D. B. Plewes. Ob-servation of nonlinear shear wave propagation using magnetic resonance elastography.Magnetic Resonance in Medicine, 52(4):842–850, 2004.[134] R. S. Sahebjavaher, A. Baghani, M. Honarvar, R. Sinkus, and S. E. Salcudean.Transperineal prostate MR elastography: initial in vivo results. Magnetic Resonancein Medicine, 69(2):411–420, 2013.[135] R. S. Sahebjavaher, S. Frew, A. Bylinskii, L. Beek, P. Garteiser, M. Honarvar,R. Sinkus, and S. E. Salcudean. Prostate MR elastography with transperineal elec-tromagnetic actuation and a fast fractionally encoded steady-state gradient echo se-quence. NMR in Biomedicine, 27(7):784–794, 2014.[136] R. S. Sahebjavaher, P. Garteiser, B. Van Beers, S. E. Salcudean, and R. Sinkus.Rapid 3D motion-encoding using spoiled FFE: application towards multi-frequencyMR rheology in liver. In European Society for Magnetic Resonance in Medicine andBiology, page 124, 2011.[137] R. S. Sahebjavaher, G. Nir, M. Honarvar, L. O. Gagnon, J. Ischa, E. C. Jones, S. D.Chang, L. Fazli, S. L. Goldenberg, R. Rohling, P. Kozlowski, R. Sinkus, and S. E. Sal-cudean. MR elastography of prostate cancer: quantitative comparison to histopathol-ogy, repeatability of methods. NMR in Biomedicine, 28(1):124–139, 2015.126[138] R. S. Sahebjavaher, G. Nir, M. Honarvar, E. C. Jones, S. D. Chang, L. O. Gagnon,J. Ischa, A. Yung, L. Fazli, S. L. Goldenberg, R. Rohling, R. Sinkus, P. Kozlowski,and S. E. Salcudean. MR elastography and diffusion-weighted imaging of ex vivoprostate cancer: quantitative comparison to histopathology. NMR in Biomedicine,28(1):89–100, 2015.[139] S. Salcudean, D. French, S. Bachmann, R. Zahiri-Azar, X. Wen, and W. Morris. Vis-coelasticity modeling of the prostate region using vibro-elastography. Medical ImageComputing and Computer Assisted Intervention (MICCAI), pages 389–396, 2006.[140] S. E. Salcudean, R. S. Sahebjavaher, O. Goksel, A. Baghani, S. S. Mahdavi, G. Nir,R. Sinkus, and M. Moradi. Biomechanical modeling of the prostate for procedureguidance and simulation. In Y. Payan, editor, Soft Tissue Biomechanical Modelingfor Computer Assisted Surgery, volume 11 of Studies in Mechanobiology, Tissue En-gineering and Biomaterials, pages 169–198. Springer Berlin Heidelberg, 2012.[141] R. Sandhu, S. Dambreville, and A. Tannenbaum. Point set registration via parti-cle filtering and stochastic dynamics. IEEE Transactions on Pattern Analysis andMachine Intelligence, 32(8):1459–1473, 2010.[142] S. Sara Mahdavi, N. Chng, I. Spadinger, W. J. Morris, and S. E. Salcudean.Semi-automatic segmentation for prostate interventions. Medical Image Analysis,15(2):226–237, 2011.[143] S. Sara Mahdavi, M. Moradi, X. Wen, W. Morris, and S. Salcudean. Evaluationof visualization of the prostate gland in vibro-elastography images. Medical ImageAnalysis, 15(4):589–600, 2011.[144] K. Scheffler. Superresolution in MRI? Magnetic Resonance in Medicine, 48(2):408–408, 2002.[145] A. Schned, K. Wheeler, C. Hodorowski, J. Heaney, M. Ernstoff, R. Amdur, andR. Harris. Tissue-shrinkage correction factor in the calculation of prostate cancervolume. The American Journal of Surgical Pathology, 20(12):1501–1506, 1996.[146] B. Scho¨lkopf and A. J. Smola. Learning with kernels: Support vector machines, regu-larization, optimization, and beyond. MIT press, 2002.[147] J. A. Sethian. Level Set Methods and Fast Marching Methods Evolving Interfaces inComputational Geometry, Fluid Mechanics, Computer Vision, and Materials Science.Cambridge University Press, New York, second edition, 1999.127[148] V. Shah, T. Pohida, B. Turkbey, H. Mani, M. Merino, P. Pinto, P. Choyke, andM. Bernardo. A method for correlating in vivo prostate magnetic resonance imagingand histopathology using individualized magnetic resonance-based molds. Review ofScientific Instruments, 80:104301, 2009.[149] M. M. Siddiqui, S. Rais-Bahrami, H. Truong, L. Stamatakis, S. Vourganti, J. Nix,A. N. Hoang, A. Walton-Diaz, B. Shuch, M. Weintraub, et al. Magnetic resonanceimaging/ultrasound–fusion biopsy significantly upgrades prostate cancer versus sys-tematic 12-core transrectal ultrasound biopsy. European Urology, 64(5):713–719, 2013.[150] R. Sinkus, J. Lorenzen, D. Schrader, M. Lorenzen, M. Dargatz, and D. Holz. High-resolution tensor MR elastography for breast tumour detection. Physics in Medicineand Biology, 45:1649–1664, 2000.[151] R. Sinkus, M. Tanter, S. Catheline, J. Lorenzen, C. Kuhl, E. Sondermann, andM. Fink. Imaging anisotropic and viscous properties of breast tissue by magneticresonance-elastography. Magnetic Resonance in Medicine, 53(2):372–387, 2005.[152] S. Smith, K. Wallner, G. Merrick, W. Butler, S. Sutlief, and P. Grimm. Interpretationof pre- versus postimplant TRUS images. Medical Physics, 30(5):920–924, 2003.[153] S.-E. Song, N. B. Cho, G. Fischer, N. Hata, C. Tempany, G. Fichtinger, and I. Ior-dachita. Development of a pneumatic robot for MRI-guided transperineal prostatebiopsy and brachytherapy: New approaches. In IEEE International Conference onRobotics and Automation (ICRA), pages 2580–2585. IEEE, 2010.[154] A. Sotiras, C. Davatzikos, and N. Paragios. Deformable medical image registration:A survey. IEEE Transactions on Medical Imaging, 32(7):1153–1190, 2013.[155] Stanford University Computer Graphics Laboratory. The Stanford 3D scanning repos-itory. http://graphics.stanford.edu/data/3Dscanrep/ (retrieved on November16, 2012).[156] Y. Sun, J. Yuan, M. Rajchl, W. Qiu, C. Romagnoli, and A. Fenster. Efficient con-vex optimization approach to 3D non-rigid MR-TRUS registration. Medical ImageComputing and Computer Assisted Intervention (MICCAI), pages 195–202, 2013.[157] F. Taquee, O. Goksel, S. Mahdavi, M. Keyes, W. Morris, I. Spadinger, and S. Sal-cudean. Deformable prostate registration from MR and TRUS images using surfaceerror driven FEM models. SPIE Medical Imaging, 831612:1–10, 2012.128[158] L. Taylor, B. Porter, G. Nadasdy, P. di Sant’Agnese, D. Pasternack, Z. Wu, R. Baggs,D. Rubens, and K. Parker. Three-dimensional registration of prostate images fromhistology and ultrasound. Ultrasound in Medicine & Biology, 30(2):161–168, 2004.[159] C. Tempany, S. Straus, N. Hata, and S. Haker. MR-guided prostate interventions.Journal of Magnetic Resonance Imaging, 27(2):356–367, 2008.[160] G. M. Treece, R. W. Prager, and A. H. Gee. Stradwin. http://mi.eng.cam.ac.uk/~rwp/stradwin/ (retrieved on July 12, 2011). Medical Imaging Group, MachineIntelligence Laboratory, Cambridge University, UK.[161] G. M. Treece, R. W. Prager, and A. H. Gee. Regularised marching tetrahedra: im-proved iso-surface extraction. Computers & Graphics, 23(4):583–598, 1999.[162] H. Trivedi, B. Turkbey, A. Rastinehad, C. Benjamin, M. Bernardo, T. Pohida,V. Shah, M. Merino, B. Wood, W. Linehan, et al. Use of patient-specific MRI-based prostate mold for validation of multiparametric MRI in localization of prostatecancer. Urology, 79(1):233–239, 2012.[163] A. Tsai, A. Yezzi Jr, W. Wells, C. Tempany, D. Tucker, A. Fan, W. Grimson, andA. Willsky. A shape-based approach to the segmentation of medical imagery usinglevel sets. IEEE Transactions on Medical Imaging, 22(2):137–154, 2003.[164] O. Ukimura, M. M. Desai, S. Palmer, S. Valencerina, M. Gross, A. L. Abreu, M. Aron,and I. S. Gill. 3-dimensional elastic registration system of prostate biopsy loca-tion by real-time 3-dimensional transrectal ultrasound guidance with magnetic reso-nance/transrectal ultrasound image fusion. The Journal of Urology, 187(3):1080–1086,2012.[165] P. Viola and W. M. Wells III. Alignment by maximization of mutual information.International Journal of Computer Vision, 24(2):137–154, 1997.[166] H. Wang, J. B. Weaver, M. M. Doyley, F. E. Kennedy, and K. D. Paulsen. Optimizedmotion estimation for MRE data with reduced motion encodes. Physics in Medicineand Biology, 53(8):2181–2196, 2008.[167] A. Ward, C. Crukley, C. McKenzie, J. Montreuil, E. Gibson, J. Gomez, M. Moussa,G. Bauman, and A. Fenster. Registration of in vivo prostate magnetic resonanceimages to digital histopathology images. Prostate Cancer Imaging. Computer-AidedDiagnosis, Prognosis, and Intervention, pages 66–76, 2010.129[168] A. Ward, C. Crukley, C. McKenzie, J. Montreuil, E. Gibson, C. Romagnoli, J. Gomez,M. Moussa, J. Chin, G. Bauman, et al. Prostate: Registration of digital histopatho-logic images to in vivo MR images acquired by using endorectal receive coil. Radiology,263(3):856–864, 2012.[169] W. Wein, S. Brunke, A. Khamene, M. R. Callstrom, and N. Navab. Automatic CT-ultrasound registration for diagnostic imaging and image-guided intervention. MedicalImage Analysis, 12(5):577–585, 2008.[170] M. White, D. Hawkes, A. Melbourne, D. Collins, C. Coolens, M. Hawkins, M. Leach,and D. Atkinson. Motion artifact correction in free-breathing abdominal MRI usingoverlapping partial samples to recover image deformations. Magnetic Resonance inMedicine, 62(2):440–449, 2009.[171] G. Xiao, B. Bloch, J. Chappelow, E. Genega, N. Rofsky, R. Lenkinski, J. Tomaszewski,M. Feldman, M. Rosen, and A. Madabhushi. Determining histology-MRI slice corre-spondences for defining MRI-based disease signatures of prostate cancer. Computer-ized Medical Imaging and Graphics, 35(7):568–578, 2011.[172] H. Xu, A. Lasso, S. Vikal, P. Guion, A. Krieger, A. Kaushal, L. Whitcomb, andG. Fichtinger. MRI-guided robotic prostate biopsy: A clinical accuracy validation.Medical Image Computing and Computer Assisted Intervention (MICCAI), pages 383–391, 2010.[173] T. K. Yasar, D. Klatt, R. L. Magin, and T. J. Royston. Selective spectral displacementprojection for multifrequency MRE. Physics in Medicine and Biology, 58(16):5771,2013.[174] A. Yezzi, A. Tsai, and A. Willsky. A statistical approach to snakes for bimodal andtrimodal imagery. International Conference on Computer Vision, 2:898–903, 1999.[175] Y. Zhan, Y. Ou, M. Feldman, J. Tomaszeweski, C. Davatzikos, and D. Shen. Reg-istering histological and MR images of prostate for image-based cancer detection.Academic Radiology, 14(11):1367, 2007.[176] O. Zienkiewicz and R. Taylor. The finite element method for solid and structuralmechanics, volume 2. Butterworth-Heinemann, 2005.130Appendix ADerivation of the ParameterizedPhase AccumulationHere we derive expressions of the parameterized phase accumulation. Based on the Lar-mor precession, we calculate the phase accumulated by each isochromat at r¯ that movesaccording to the wave displacement Ψ¯ under the MEGs G¯ (excluding the constant phaseaccumulated by the static magnetic field B¯0). We have thatΦ(r¯;ωi, N¯ , |G¯|, α¯, β¯) = γ∫ ∞−∞G¯(t; ω¯, N¯ , |G¯|, α¯, β¯) · r¯(t) dt= γ∫ ∞−∞G¯(t; ω¯, N¯ , |G¯|, α¯, β¯) ·(r¯ + Ψ¯(r¯, t))dt=∑i=x,y,zΦi(r¯;ωi, Ni, |Gi|, αi, βi)By using (6.1) and (6.2), the elements in G¯ · r¯ are cancelled (integration of a constantmultiplied by sinusoid over its time period). For i = x, y, z, we have the individual phaseaccumulationsΦi(r¯;ωi, Ni, |Gi|, αi, βi) = γ∫ αi/ωi+NiTiαi/ωi|ξi(r¯)| sin (ωmt+ ϕi(r¯))|Gi| sin (ωit− αi + βi) dt= γ2 |Gi||ξi(r¯)|∫ αi/ωi+NiTiαi/ωi[cos ((ωm − ωi)t+ ϕi(r¯) + αi − βi)− cos ((ωm + ωi)t+ ϕi(r¯)− αi + βi)] dtIn case ωi = ωm, the second cosine is cancelled in the integration and we haveΦi(r¯;Ni, |Gi|, αi, βi) =γpiNi|Gi|ωm|ξi(r¯)| cos (ϕi(r¯) + αi − βi)131Appendix A. Derivation of the Parameterized Phase AccumulationOtherwise, for ωi > ωm, i.e., qi = ωm/ωi < 1, i = x, y, z, we haveΦi(r¯; qi < 1, Ni, |Gi|, αi, βi) =γ2 |Gi||ξi(r¯)|[ 1(ωm − ωi)sin ((ωm − ωi)t+ ϕi(r¯) + αi − βi)− 1(ωm + ωi)sin ((ωm + ωi)t+ ϕi(r¯)− αi + βi)]αi/ωi+NiTit=αi/ωi= γ2ωi|Gi||ξi(r¯)|[ 1(qi − 1)[sin (qi(αi + 2piNi) + ϕi(r¯)− βi)− sin (qiαi + ϕi(r¯)− βi)]− 1(qi + 1)[sin (qi(αi + 2piNi) + ϕi(r¯) + βi)− sin (qiαi + ϕi(r¯) + βi)]]= γωi|Gi||ξi(r¯)| sin (piqiNi)[ 1(qi − 1)cos (qiαi + ϕi(r¯)− βi + piqiNi)− 1(qi + 1)cos (qiαi + ϕi(r¯) + βi + piqiNi)]= γωm|Gi||ξi(r¯)| sin (piqiNi)−qi(1− q2i )[(1 + qi) cos (qiαi + ϕi(r¯)− βi + piqiNi)+ (1− qi) cos (qiαi + ϕi(r¯) + βi + piqiNi)]= γωm|Gi||ξi(r¯)| sin (piqiNi)−2qi(1− q2i )[cos (βi) cos (qiαi + ϕi(r¯) + piqiNi)+ qi sin (βi) sin (qiαi + ϕi(r¯) + piqiNi)]= γωm|Gi| sin (piqiNi)−2qisign(cos (βi))(1− q2i )√cos2 (βi) + q2i sin2 (βi)· |ξi(r¯)| cos (qiαi + ϕi(r¯) + piqiNi − arctan (qi tan (βi)))= −2γqi cos (βi) sin (piqiNi)|Gi|ωm(1− q2i )√1 + q2i tan2 (βi)· |ξi(r¯)| cos (qiαi + ϕi(r¯) + piqiNi − arctan (qi tan (βi))) , i = x, y, zwhere we have used trigonometric identities and the harmonic addition theorem.Specifically, the phase accumulation for qi < 1 when βi = {0, pi} (no discontinuities) isΦ(r¯; q¯ < 1, N¯ , |G¯|, α¯, β¯ = {0, pi}) =∑i=x,y,z∓2γqi sin (piqiNi)ωm(1− q2i )|Gi||ξi(r¯)| cos (qiαi + ϕi(r¯) + piqiNi)with a ’−’ sign when βi = 0, and a ’+’ sign when βi = pi. The term∓2γqi sin (piqiNi)ωm(1− q2i )|Gi|is also known as the encoding efficiency [131], and we note that it depends on qi.132
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Prostate registration using magnetic resonance elastography...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Prostate registration using magnetic resonance elastography for cancer localization Nir, Guy 2014
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Prostate registration using magnetic resonance elastography for cancer localization |
Creator |
Nir, Guy |
Publisher | University of British Columbia |
Date Issued | 2014 |
Description | Noninvasive detection and localization of prostate cancer in medical imaging is an important, yet difficult task. Benefits range from diagnosis of cancer, to planning and guidance of its treatment. In order to characterize cancer and evaluate its localization in volumetric images, such as ultrasound or magnetic resonance imaging (MRI), their spatial correspondence with the "gold standard" provided by histopathology must be established. In this thesis, we propose a general framework for a multi-slice to volume registration that is applied to register a stack of sparse, unaligned two-dimensional histological slices to a three-dimensional volumetric imaging of the prostate. The approach uses particle filtering that allows deriving optimal pose parameters of the slices in a Bayesian approach. We then propose a novel registration method between in vivo and ex vivo MRI of the prostate to facilitate its registration to histopathology. The method incorporates elasticity information, acquired by magnetic resonance elastography (MRE), to generate a patient-specific biomechanical model of the prostate and periprostatic tissue. Next, we propose a registration method between preoperative MRI and intraoperative transrectal-ultrasound. The method can be incorporated with a robotic surgical system to augment the surgeon's visualization during robot-assisted prostatectomy. We also study the use of elasticity-based registration of ultrasound elastography and MRE. We then present an image processing approach for enhancing MRE data. The approach employs registration to compensate for motion of patients during the scan to improve the accuracy of the reconstructed elastogram. A super-resolution technique is employed to increase the resolution of the acquired images by utilizing unique properties of MRE. Finally, we develop a theory for optimization-based design of motion encoding in MRE that allows reducing scanning time and increasing signal-to-noise ratio of elasticity reconstruction. We formulate the displacement estimation of the mechanical wave as an experimental design problem, by which we quantify performance of sequences, and optimize multidirectional designs. The proposed methods have been evaluated in simulations and on a diverse set of clinical data. Results may pave the way for a broader clinical deployment of elastography and elastography-based image processing. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2015-01-06 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0167097 |
URI | http://hdl.handle.net/2429/51769 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2015-02 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2015_february_nir_guy.pdf [ 5.83MB ]
- Metadata
- JSON: 24-1.0167097.json
- JSON-LD: 24-1.0167097-ld.json
- RDF/XML (Pretty): 24-1.0167097-rdf.xml
- RDF/JSON: 24-1.0167097-rdf.json
- Turtle: 24-1.0167097-turtle.txt
- N-Triples: 24-1.0167097-rdf-ntriples.txt
- Original Record: 24-1.0167097-source.json
- Full Text
- 24-1.0167097-fulltext.txt
- Citation
- 24-1.0167097.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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0167097/manifest