University of Southern California Dissertations and Theses

Variational techniques for cardiac image analysis: algorithms and applications

(USC Thesis Other)

Variational techniques for cardiac image analysis: algorithms and applications

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content V ARIATIONAL TECHNIQUES FOR CARDIAC IMAGE ANALYSIS: ALGORITHMS AND APPLICATIONS by Jonghye Woo A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ELECTRICAL ENGINEERING) December 2009 Copyright 2009 Jonghye Woo Dedication To my family ii Acknowledgements This dissertation would not have been possible without the help and encouragement of a number of people. First of all, my deepest gratitude to my advisor, Prof. C.-C. Jay Kuo, for his continued support and guidance throughout my Ph.D. research. Not only has he shown the passion for research, but also he has been a good mentor to me. I would like to thank Prof. Piotr Slomka who generously supported me, for his kindness and guidance. Also I am indebted to Prof. Byung-Woo Hong who has been incredibly helpful with enthusiastic discussions about research. He always put time and effort into sharing his knowledge with me. I would also thank Prof. Krishna Nayak and Prof. K. Kirk Shung for being the members in my dissertation committee, and Prof. Richard Leahy for being the member in my qualifying exam committee. It is indeed a great privilege to have their valuable comments and feedback on my work. Over the last three years, I’ve had the opportunity to work at two excellent labs filled with great people. First, the USC Media Communications Lab is a perfect place to learn a variety of fields. Discussions with people about work and about anything but work have greatly con- tributed to the overall enjoyment of my experience. I am grateful to my friends and colleagues: Do-kyoung Kwon, Sungwon Lee, Taehoon Shin, In suk Chong, Jeong-Yoon Lee, Woo-shik Kim, Byung-Ho Cha, Byung Tae Oh, Namgook Cho, Joohyun Cho, Jong Dae Oh, SeongHo Cho, Sungje Cho, Yongjin Cho, Youngmin Kwak, Dong-Woo Kang, Jewon Kang, and Shahryar iii Karimi-Ashtiani. Also the Artifical Intelligence in Medicine (AIM) Lab at Cedars-Sinai Medi- cal Center is a tremendous resource that brings together researchers and clinicians. Through the AIM, I was able to work with people from radiology and software engineering. I would like to thank my collaborators: Amit Ramesh, Manisha Pal, Mithun Prasad, Yuan Xu, Sherry Casanova, Paul Kavanagh, Geoff Pollard, Parker Waechter, Chih-Chun Wei, Dr. Serge Van Kriekinge, Dr. Damini Dey, Dr. Victor Cheng, Dr. Guido Germano, and Dr. Daniel Berman. Last but not least, I am deeply grateful to my beloved wife, Seon and my family for their encouragement, motivation and love along this journey. iv Table of Contents Dedication ii Acknowledgements iii List of Tables viii List of Figures ix Abstract xiv Chapter 1 : Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Clinical Background and Significance . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Technical Background and Challenges . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Contributions of the Research . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.5 Outline of the Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Chapter 2 : Curve Evolution with A Dual Shape Similarity and Its Application to Segmenta- tion of Left Ventricle 14 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Variational Formulation and Bayesian Inference for Coupled Segmentation . . . 17 2.3 Proposed Shape Constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.1 Synthetic Example for the Shape Constraint . . . . . . . . . . . . . . . . 20 2.4 Energy Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.1 Energy Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.5 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6 Conclusion and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Chapter 3 : Multi-modal Data Integration for Computer-Aided Ablation of Atrial Fibrillation 27 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Review of Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3 Multi-modal Data Integration Model for Simultaneous Surface Reconstruction and Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3.1 Surface Reconstruction and Registration . . . . . . . . . . . . . . . . . . 33 3.3.2 Derivation of Energy Terms . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3.3 Numerical Implementation . . . . . . . . . . . . . . . . . . . . . . . . . 36 v 3.3.4 Relationship between Variational Formulation and Bayesian Inference . . 38 3.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.4.1 Synthetic data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.4.2 Patient data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.5 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 Chapter 4 : Geometric Feature-based Multi-modal Image Registration of Contrast-enhanced Cardiac CT with Gated Myocardial Perfusion SPECT 46 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2.1 Description of the Method . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2.1.1 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2.1.2 Registration Model . . . . . . . . . . . . . . . . . . . . . . . 53 4.2.1.3 Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.2.1.4 Cardiac Phase Matching with Thin-Plate-Spline Warping . . . 56 4.2.2 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2.2.1 Patient selection . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2.2.2 Coronary CTA: acquisition and reconstruction . . . . . . . . . 59 4.2.2.3 Stress CT perfusion: acquisition and reconstruction . . . . . . 60 4.2.2.4 Myocardial perfusion SPECT protocol . . . . . . . . . . . . . 61 4.2.2.5 Evaluation of automatic registration . . . . . . . . . . . . . . . 61 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Chapter 5 : Nonlinear Ultrasound Image Registration Based on Intensity and Local Phase Information 75 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.2 Review of Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.3 Local Phase Computation and Properties . . . . . . . . . . . . . . . . . . . . . . 78 5.3.1 Local Phase and Local Energy . . . . . . . . . . . . . . . . . . . . . . . 78 5.3.2 Monogenic Signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.4 Registration Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.4.1 Variational Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.4.2 Euler-Langrange Optimization . . . . . . . . . . . . . . . . . . . . . . . 84 5.4.3 Numerical Approximation . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.5 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.5.1 Synthetic Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.5.2 In vivo Image . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Chapter 6 : Nonlinear Registration of Serial Coronary CT Angiography (CCTA) for Plaque Development Monitoring 94 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.2.1 Description of Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 vi 6.2.2 Feature Mask . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.2.3 Global Displacement Model . . . . . . . . . . . . . . . . . . . . . . . . 100 6.2.4 Local Deformation Model . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.2.5 Constrained V olume Preserving Flow . . . . . . . . . . . . . . . . . . . 103 6.2.6 Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.2.7 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.3 Experimental Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.3.1 Clinical Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.3.2 Evaluation of Registration Quality . . . . . . . . . . . . . . . . . . . . . 109 6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 Chapter 7 : Conclusion and Future Work 121 7.1 Summary of the Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 7.2 Future Research Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 References 125 vii List of Tables 2.1 Comparison of the manual segmentation and proposed algorithm . . . . . . . . . . . . 24 3.1 Performance comparison of ICP and the proposed method . . . . . . . . . . . . 45 4.1 Registration error and inter-observer variability . . . . . . . . . . . . . . . . . . 67 4.2 Evaluation of best matching phase . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.1 Comparison of three different data terms (unit:pixel). . . . . . . . . . . . . . . . 89 5.2 Comparison of three different similarity measures for the mice image registration 92 5.3 Comparison of three different similarity measures for the human subject image registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.1 Patient characteristics and image parameters of baseline and follow-up scans (for 5 cases only) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.2 TRE comparison with different registration methods (Mean SD (mm)) . . . . . 112 6.3 Inverse consistency error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.4 Average landmark distances between two observers for each landmark (n=56) . . 116 viii List of Figures 1.1 (a) Example of coronary plaque and (b) fusion of MPS-CTA . . . . . . . . . . . 3 1.2 (a) A CARTO map of the left ventricle of an animal and (b) an image of the left atrium and pulmonary veins. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 Illustration of notion of distance between two shapes. (a) Examples of two shapes, a circle S 1 in red and a rectangle S 2 in blue, represented by f 1 and f 2 respectively. (b) Profiles of distances that are measured by the shortest Euclidean distance from S 1 to S 2 at each point on S 1 (red line) and the shortest Euclidean distance from S 2 to S 1 at each point on S 2 (blue line) along the curves parameter- ized by arc length. (c) the difference between two signed distance functions that represents the two shapes. The outlines of the two shapes are superimposed and the values on the shape outlines represent the shortest Euclidean distance from one shape to the other. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2 Examples of the shape constraint with different s. Outer structure of interest is circle in both examples but the shape suffers from fuzzy boundary (top row) and noise (bottom row). The first column shows segmentation results of multiphase Chan-Vese model. Our shape constraint can capture the desired outer shape by adjusting s. As s becomes bigger, the allowable distance between two shapes becomes larger. In the second column, the outer shape is restricted to the circle due to the small s whereas the variation of distance becomes bigger in last two columns. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Bland-Altman plots for EDV and LVEF . . . . . . . . . . . . . . . . . . . . . . 24 2.4 Experimental results with ED and ES phases of three subjects. Manual segmen- tation results of endocardium (green) and epicardium (yellow) are shown in first and third column and results of proposed algorithm are presented with endo- cardium (red) and epicardium (white) using same images in second and fourth column. 3D volume is presented in the rightmost column. Parameters were se- lected empirically same for all (a = 400,b = 30,l = 400, ands = 2). . . . . . . 25 ix 3.1 The shape reconstruction results for a synthetic 2D star shape, (a) the original shape, (b)-(d) ICP results using different Gaussian noise levels, (e) the deformed shape, and (f)-(h) results of the proposed method using different Gaussian noise levels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2 Comparison of mean distances between the reconstructed surface and measured data points for the 2D star example at different noise levels using ICP and the proposed algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3 The shape reconstruction results for a synthetic 3D image: (a) the original shape, (b)-(e) ICP results using different Gaussian noise levels, (f) the deformed shape, and (g)-(j) results of the proposed method using different Gaussian noise levels. . 42 3.4 Comparison of mean distances between the reconstructed surface and measured data points for the 3D jar example at different noise levels using ICP and the proposed algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.5 The 3D patient data: (a) MRA of LV and (b) 3D reconstruction result of the given MRA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.6 The surface registration results for the patient data. . . . . . . . . . . . . . . . . 44 4.1 Examples of CT (top) and MPS (bottom) volume cross-sections illustrating image variations between patients: (a) and (b) show coronal views for two different patients. (c) and (d) show the axial/transverse views for the respective patients. . 47 4.2 An example of myocardial segmentation of MPS: (a) shows transverse views for an abnormal case (top) and a normal case (bottom); (b) shows respective con- tours from segmentation overlaid on the images from (a); (c) illustrates the blood pool (D 1 ), myocardial (D 2 ) and extracardiac (D 3 ) regions; (d) shows 3D contours obtained from the myocardial segmentation, endocardium (gray surface) and epi- cardium (orange mesh). Note the robust segmentation result despite abnormality in MPS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.3 Overview of the registration procedure . . . . . . . . . . . . . . . . . . . . . . 53 4.4 An example of finding the best matching phase of gated MPS by minimizing the registration cost functional: (a) plot of cost functional value across different phases of MPS; (b) overlaid registration result of phase 9 MPS on CT perfusion; (c) overlaid registration result of best matching phase (16) with CT perfusion. Note that the best matching MPS phase is visually well correlated with CT. . . . 56 x 4.5 “Motion-frozen” technique: All MPS frames are registered to a common refer- ence frame, which is the best matching MPS phase for the CT volume. TPS-based warping ( f i ) is used to deform and map each frame to I N . The control points from the segmentation contours serve as landmark points for TPS warping, and cor- responding points from the source and target contours determine the mapping. Final “motion-frozen” image is obtained by averaging over all the warped vol- umes. Note that perfusion defects are clearly seen in the final image (final column). 58 4.6 An example of registration results: original CT image (top), images before (mid- dle) and after (bottom) registration are shown. Errors were (2, 1, 0) mm for translation and (1, 1, 0) degrees for rotation, as compared to visual alignment. . . 64 4.7 An example of registration results (CT Perfusion): original CT image (top), im- ages before (middle) and after (bottom) registration are shown. Errors were (1, 0, 4) mm for translation and (5, 2, 3) degrees for the rotation, as compared to visual alignment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.8 Comparison of NMI based method and proposed method: (a) registration result using NMI based method; (b) the NMI values for translational offsets from the ground truth for each axis; (c) registration result using the proposed method; (d) the cost values of proposed method for translational offsets from the ground truth for each axis. Note that the anatomical features are very different and the image segmentation will aid the registration of MPS to CT. Additionally cost value of the proposed method is optimal at the origin whereas NMI does not provide any such robust optimum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.9 Registration errors: (a) translational errors of stress/rest MPS, and CT obtained from site A; (b) rotational errors of stress/rest MPS, and CTA obtained from site A; (c) translational errors of stress/rest MPS, and CT obtained from site B; (d) rotational errors of stress/rest MPS, and CT obtained from site B; (e) translational errors of stress/rest MPS and perfusion CT; (f) rotational errors of stress/rest MPS and perfusion CT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.10 Capture ranges of the proposed method: (a) translation errors for different mis- alignment initializations along X axis; (b) translation errors for different mis- alignment initializations along Y axis; (c) translation errors for different mis- alignment initializations along Z axis. . . . . . . . . . . . . . . . . . . . . . . . 69 4.11 Capture ranges of the NMI based method: (a) translation errors for different mis- alignment initializations along X axis; (b) translation errors for different mis- alignment initializations along Y axis; (c) translation errors for different mis- alignment initializations along Z axis. . . . . . . . . . . . . . . . . . . . . . . . 69 xi 5.1 (a) An examplary cardiac ultrasound image, and its local phases using two DoG bandpass filters with (b) a small standard deviation (s 1 = 2 p 2;s 2 = 2) and (c) a large standard deviation (s 1 = 4 p 2;s 2 = 4). . . . . . . . . . . . . . . . . . . . . 78 5.2 (a) The first synthetic image and (b) its deformed image. . . . . . . . . . . . . . 88 5.3 (a) The second synthetic image and (b) its deformed image. . . . . . . . . . . . . 89 5.4 In vivo mice ultrasound images: (a) the source image, (b) the target image, (c) the difference image between source and target images, (d) results obtained by the proposed method, (e) the displacement vector from the boundary in (a) (arrows being scaled 1.5 times for the visualization purpose), and (f) the deformation field. 91 5.5 In vivo human subject ultrasound images: (a) the source image, (b) the target image, (c) the difference image between source and target images, (d) results ob- tained by the proposed method, (e) the displacement vector from the myocardium boundary in (a) (arrows being scaled 1.5 times for the visualization purpose), and (f) the deformation field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.1 The flowchart of the proposed method. . . . . . . . . . . . . . . . . . . . . . . . 98 6.2 Illustration of the coronary tree mask used in image registration: (a) a representa- tive 2D slice (yellow plane) of the pre-segmented 3D binary mask that covers the coronary tree region; (b) the representative slice from a 3D image overlaid with the pre-segmented mask boundary (black); (c) a 2D pre-segmented binary mask for this slice, and (d) the soft mask generated by convolving with the Gaussian kernel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.3 Anatomical landmarks used for validation shown in one example case with im- ages from a baseline scan in the top row, and images from a follow-up scan in the bottom row. We have obtained independent sets of anatomical landmarks from two clinical observers. Red crosses indicate the positions of landmarks, includ- ing the left main coronary artery origin (a), bifurcation of left anterior descending and left circumflex arteries (b), right coronary artery origin (c), and branch point of the first diagonal artery from the left anterior descending artery (d). . . . . . . 109 6.4 Comparison between total energy minimization without and with the volume- preserving constraint. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.5 Comparison of global displacement and local deformation registration techniques: (a) the global displacement model, (b) the result of subsequent local deformation, and (c) the target image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 xii 6.6 Example of registration results. The baseline image is registered to the subse- quent image. The transverse orientation is shown in the top row while the coronal orientation is shown in the bottom row. The original baseline images before and after registration are shown in (a) and (b), respectively. Panel (c) shows the sub- sequent image. As shown in (a), the original image is significantly misaligned. The “roving window” (yellow) technique in (d) using a portion of the registered image from (b) to the target image in (c) shows that the result from our method is visually accurate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 6.7 Example of plaque lesion progression (worsening) in an 82-year old man. The top row shows a follow-up study while the bottom row shows a registered baseline study. The time difference between the two scans was 11 months. The pixel sizes of baseline and follow-up studies were 0.40.40.29 mm and 0.450.450.3 mm, respectively. Green arrows show increased, non-calcified plaque and steno- sis in the follow-up study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.8 Example of a plaque without significant changes in a 59-year old man. The top row shows the follow-up study while the bottom row shows the registered base- line study. The lesion is marked with green arrows. The time difference between the two scans was 12 months. The pixel sizes of baseline and follow-up study were 0.310.310.3 mm and 0.340.340.3 mm, respectively. . . . . . . . . . 115 xiii Abstract In this dissertation we investigate several image segmentation and registration techniques based on the variational formulation for medical imaging applications. The five main research results are summarized below. First, a novel segmentation approach is proposed to jointly delineate the boundaries of epi- and endocardium of the left ventricle on Magnetic Resonance Imaging (MRI) under a variational framework using level sets. While most left ventricle segmentation approaches incorporate a shape prior obtained by a training process from an ensemble of examples, we exploit a novel shape constraint using an implicit shape prior knowledge, which assumes shape similarity be- tween epi- and endocardium allowing a variation under the Gaussian distribution. Second, we examine multi-modal data integration with an electroanatomic mapping (EAM) data and MRI images for computer-aided catheter ablation of atrial fibrillation accurately. Specif- ically, we propose a variational formulation for surface reconstruction and incorporate the prior shape knowledge, which results in a level set method. The proposed method enables simultaneous reconstruction and registration under nonlinear deformation. Third, a fully automated registration method is presented utilizing geometric features from a reliable segmentation of gated myocardial perfusion SPECT (MPS) volumes, where regions of myocardium and blood pools are extracted and used as an anatomical mask to de-emphasize xiv the inhomogeneities of intensity distribution caused by perfusion defects and physiological vari- ations. A multi-resolution approach is employed to represent coarse-to-fine details of both vol- umes. The extracted voxels from each level are aligned using a similarity measure with a piece- wise constant image model and minimized using a gradient descent method. We then perform limited nonlinear registration of gated MPS to adjust for phase differences by automatic cardiac phase matching between CT and MPS. For phase matching, we incorporate nonlinear registration using thin-plate-spline-based warping Fourth, a nonlinear ultrasound image registration method is proposed using the intensity and the local phase information under a variational framework. One application of this technique is to register two consecutive images in an ultrasound image sequence. Although intensity is the most widely used feature in traditional ultrasound image registration algorithms, speckle noise and lower image resolution make the registration process difficult. By integrating the intensity and the local phase information, we can find and track the nonlinear transformation of each pixel under diffeomorphism between the source and target images. Finally, we develop a fully automatic and accurate nonlinear volume registration for longitu- dinal Coronary CT angiography (CCTA) scan pairs. The proposed algorithm combines global dis- placement and local deformation using nonlinear volume co-registration with a volume-preserving constraint. Histogram matching of intensities between two serial scans is performed prior to non- linear co-registration with dense nonparametric diffeomorphism in which sum of squared differ- ence is used as a similarity measure. The segmented coronary artery trees provide initial anatom- ical landmarks for the co-registration algorithm that help localize and emphasize the structure of interest. To avoid possible bias caused by incorrect segmentation, we convolve the Gaussian kernel with the segmented binary coronary tree mask and define an extended weighted region xv of interest. A multi-resolution approach is employed to represent coarse-to-fine details of both volumes. The energy functional is optimized using a gradient descent method. Extensive computer simulations have been conducted and clinical validations have been per- formed to demonstrated the improved accuracy of the proposed techniques. xvi Chapter 1 : Introduction 1.1 Motivation Medical image analysis aims at extracting useful information from medical images with im- proved accuracy and consistency so as to facilitate doctors to make sound clinical decision. Academia, hospitals and the industry have made a substantial amount of efforts in the devel- opment of computer-aided diagnosis tools. While such a field has been around for decades, it has become more popular and research results have been applied to real world applications and products since early eighties. This is partly due to recent advances in imaging hardware and mathematical tools for image analysis. However, there still exist a variety of unsolved problems. In order to process a large number of high-resolution acquired images, automated algo- rithms, which demands as little human intervention as possible, are increasingly in great demand. Medical image analysis covers a wide range of research topics, including acquisition, forma- tion/reconstruction, enhancement, compression/storage, segmentation, registration and visualiza- tion. Among them, image/video segmentation and registration are two of the most challenging problems. 1 Image/video segmentation in medical applications can be cast as a constrained optimization problem. The objective of segmentation is to partition an image into a finite number of seman- tically meaningful regions according to anatomical or functional structures. There is no general definition of a region of interest since it may depend on the image type and the application context. However, it can be defined as relatively homogeneous regions with respect to a given measure or regions that fit some prior knowledge the best. Image/video registration (also known as image/video matching, morphing and alignment) refers to the process of estimating the linear/nonlinear mapping function that transforms each point of one image/video to a point in the target image/video. It can be viewed as the process of finding the correspondence or aligning two or more images of the same scene taken at different times, from different perspective, or using different sensors. For instance, multi-modal image registration or integration can provide complementary information in the sense that images with functional, physiological and anatomical information could be fused into one display in an anal- ysis platform [35]. On the other hand, image registration is mono-modal registration used for monitoring disease progression or evaluation of therapies by longitudinal image alignment [122] or estimating motion across different images. In this dissertation, we focus on cardiac applications associated with life-threatening coro- nary artery disease and cardiac arrhythmias, and study the related segmentation and registration problems. Medical image analysis for cardiac applications is a rapidly growing interdisciplinary research due to its high death rate, where disciplines of engineering, applied mathematics and statistics meet. Seminal works were heuristic but recent studies are benefited by sophisticated techniques from partial differential equations (PDE), differential geometry, statistical inference and machine learning theory. 2 1.2 Clinical Background and Significance (a) (b) Figure 1.1: (a) Example of coronary plaque and (b) fusion of MPS-CTA Sudden cardiac death (SCD) is the single largest cause of death in the developed world. In the United States (US) alone more than 400,000 persons die annually. In addition, more than 1 million people in the US and 19 million people worldwide experience SCD or myocardial infarction each year [141]. Nearly 80% of the SCD deaths are estimated to be a consequence of coronary artery disease (CAD) [145]. 75% of these acute coronary episodes is due to rupture of a vulnerable atherosclerotic plaque [40] as shown in Fig. 1.1 (a). Although the current gold standard for evaluating coronary arteries is invasive coronary angiography (ICA) [38], Coronary CT angiography (CCTA) with multi-slice helical scanners has become the integral part of major diagnostic pathways for CAD [9]. This is because ICA requires cardiac catheterization and it is occasionally associated with complications. In contrast, noninvasive cardiac imaging techniques of coronary arteries provide minimal risks in detecting CAD and serially monitor response to 3 therapy [38]. Thus, detection of changes in plaque parameters is one of the urgent issues, which can be estimated by image registration of serial data and subsequent direct estimation of plaque development from registered serial images on a voxel basis. In addition, cardiac computed tomography (CT) and single photon emission computed to- mography (SPECT) provide clinically complementary information in the diagnosis of CAD. Fused anatomical and physiological data acquired sequentially on separate scanners can be co- registered to accurately diagnose CAD in specific coronary vessels as shown in Fig. 1.1 (b). (a) (b) Figure 1.2: (a) A CARTO map of the left ventricle of an animal and (b) an image of the left atrium and pulmonary veins. The most common mechanism of SCD in the presence of CAD is the development of lethal arrhythmias such as atrial fibrillation (AF) and ventricular tachycardia (VT). Arrhythmias are characterized by abnormalities in electrical impulse formation or conduction within the heart. Treatment of arrhythmias has dramatically evolved over past decades, and minimally-invasive catheter-based therapy is the preferred method of eliminating arrhythmias today [88]. With an 4 electroanatomical mapping (EAM) system, which precisely tracks the position of catheters in- side patient’s body, it is possible to construct three-dimensional maps of the ventricular and atrial chambers of the heart as shown in Fig. 1.2 (a), which gives a typical CARTO map of the left ven- tricle in an animal with a scar created by inflating a balloon in left anterior descending coronary artery. The map was created using a catheter which physically touches the ventricular inner wall and records bipolar voltages at a minimum of 300 points. The procedure is laborious yet only a representation of the anatomy (not the real anatomy) is given. Colors represent bipolar volt- ages (V), such that red indicates V 0.5mV (scarred tissue) and magenta indicates V1.5mV (normal tissue). A voltage level in between indicates a myocardium that is neither scarred nor fully functioning, i.e. the arrhythmic substrate indicated by blue-green. Red dots indicate sites of ablation based on abnormal electrograms. The left atrium and pulmonary veins obtained using the Cardiac MR Angiography (CMRA) are shown in Fig. 1.2(b). CMRA is a useful method of producing detailed pictures of blood vessels. Each point in these maps is annotated with bio- electrical signals recorded from electrodes located at the tip of the catheter. These maps are then used to guide catheter ablation within the heart. However, a major limitation of the EAM system is its inability to create the true cardiac anatomy. Instead, a virtual shell is created to represent the inner wall of the left ventricle or atrial wall based on intracardiac electrical signals. Also, the electroanatomical mapping procedure results in relatively sparse sampling of the heart and a significant amount of time and skill are require to generate these maps. Thus, real anatomy using a patient-specific 3D model obtained from MRI/CT, e.g. Fig. 1.2 (b), is required for accurate localization and further ablation of the target site. Furthermore, to locate lesions and ablate the target site accurately, it would be desirable to have an optimal surface model using a registration 5 or a reconstruction method that integrates EAM points and a patient-specific anatomical heart model in a complementary manner. 1.3 Technical Background and Challenges Image segmentation and registration models can be formulated in a variational framework, where an objective function (also called the energy functional) is defined by incorporating appropriate cues available in applications of interest. The objective function consists of an appropriate simi- larity measure (i.e. the data fidelity term) that indicates the degree the model is fitted to the given object of interest or the closeness the transformed source image is aligned with the target image. Additionally, it consists of a regularizing term, which is included to guarantee the well-posedness of the optimization problem or provide a deformation that is diffeomorphic [25]. Mathematically, it can be written as inf u fE(u)= F(u)+lR(u)g; (1.1) where E is the combined energy functional, F is the data fidelity functional, R is the regularization functional,l2R + is a weighting parameter, and u is the desired solution. To solve the minimization problem in (1.1), several well-known global minimization tech- niques such as simulated annealing [127] and dynamic programming [6] may not be practical due to their high computational complexity. Graph cuts, which offers the global minimum or a strong local minimum, is an emerging and promising technique. Since the seminal work [15], this technique has been widely applied to a variety of computer vision problems, including stereo matching [71] and image segmentation [14], due to its efficiency. However, it may not be appli- cable to some energy functionals [72]. 6 A suboptimal strategy based on the calculus of variation has also been developed to minimize the energy functional in (1.1). A necessary condition for the minimum, u, of the energy functional to exist is that its Gateaux derivative dE(u;d) vanishes for all purtabations ofd. This derivative, known as the first variation of E in the direction ofd, can be written as dE(u;d)= lim e!0 E(u;ed) E(u) e = ¶E(u+ed) ¶e e=0 : (1.2) If u is a local extremum of E, its first variation at u must vanish, which yields the Euler-Lagrange equation: ¶ u(x) E(u(x)) = 0: (1.3) Thus, the energy functional can be minimized using the gradient descent of the Euler-Lagrange equation by introducing an artificial time t as ¶u(x,t) ¶t =¶ u E(u): (1.4) The variational formulation based on partial differential equations and level set represen- tation is adopted in this research. This formulation is a state-of-the-art tool in medical image analysis [99] owing to several important advantages [102]. First, it is less empirical since it in- corporates a better knowledge of the model (i.e. the energy functional). This is particularly true when the prior knowledge is built upon a statistical model (e.g. the maximum a posteriori (MAP) estimate in the Bayesian formulation [129]). Second, the existence of a global minimum and the well-posedness of the minimization procedure can be mathematically proven in some cases. Finally, once the variational problem is defined, it can be tackled with a variety of minimization 7 procedures depending on the available prior information. Usually, it is important to design a con- vex energy functional or deal with a non-convex energy functional by incorporating appropriate prior knowledge to avoid the local minima. Several technical challenges and constraints imposed by applications of consideration are discussed below. Segmentation of left ventricle remains to be a challenging problem since image intensities of cardiac chambers vary due to differences in the blood velocity. In particular, delineat- ing epicardium is more challenging since there is no clear boundary due to occlusions, illumination variations and other unmodeled phenomena, which make robust and accurate segmentation difficult. In the integration process of MR/CT with EAM data, one obvious difficulty stems from noise or outliers that are inevitably associated with the MR/CT imaging process and the EAM data collection procedure. Unlike other organs in the body, human heart under- goes contractile motion, apart from the respiratory motion, thus making it unique and very challenging to register and integrate with other modalities. These motions, in addition to physiological variations including changes in the heart beat rate, the heart rhythm and the respiratory effect, are sources of outliers. Accurate registration of contrast-enhanced cardiac CT with MPS is challenging as different anatomical features are visible in different images. For example, although the myocardium is visible on both cardiac CT and MPS, there could be severe perfusion defects present on MPS images which do not have equivalent appearance on CT. Additionally, variation in 8 the intravenous contrast distribution for different patients causes intensity variations on CT images within the blood pool region. Registration of ultrasound images is challenging task to find correspondence since ultra- sound images have a low signal-to-noise (SNR) ratio due to artifacts and attenuation. Fur- thermore, speckle noise caused by sub-resolution tissue structures is a big foe in analyzing lesions and performing various processes such as segmentation and registration by limiting the detectability of low contrast lesions. Serial Coronary CT angiography registration is challenging since the structure of interest (i.e. the plaque area) is very small and there is no clear definition of anatomical detail such as plaque shape change. Specifically, in images with a plaque, the fundamental assumption of topological equivalence between the baseline and follow-up scans may be violated due to anatomical changes caused by plaque development. Moreover, different acquisition time with two serial scans may cause changes in image intensities, which makes the task of finding the correspondence difficult. 1.4 Contributions of the Research The main contributions of the research in this dissertation are summarized below. Curve evolution with a dual shape similarity for segmentation of left ventricle In Chapter 2, we present a novel segmentation approach to jointly delineate the boundaries of epi- and endocardium of the left ventricle on MRI images under a variational framework using level sets. This technique is in great demand in clinical cardiological applications. 9 One strategy to tackle segmentation under undesirable conditions such as subtle bound- aries and occlusions is to exploit prior knowledge specific to the object to segment, e.g., the knowledge about heart anatomy in this case. While most left ventricle segmentation approaches incorporate a shape prior obtained by a training process from an ensemble of examples, we exploit a novel shape constraint using an implicit shape prior knowledge, which assumes shape similarity between epi- and endocardium allowing a variation under the Gaussian distribution. Our approach does not demand a training procedure which is usually subject to the training examples and is also laborious and time-consuming in gen- erating the shape prior. Instead, we model a shape constraint by a statistical distance be- tween the shape of epi- and endocardium employing signed distance functions. We applied this technique to cardiac MRI data with quantitative evaluations performed on 10 subjects. Experimental results show robustness and effectiveness of our shape constraint within a Mumford-Shah segmentation model in the segmentation of left ventricle from cardiac MRI images in comparison with manual segmentation results. Multi-modal data integration for computer-aided ablation of atrial fibrillation Image-guided percutaneous interventions have successfully replaced invasive surgical meth- ods in some cardiologic practice, where the use of 3D reconstructed cardiac images gener- ated by MRI and CT plays an important role. To conduct computer-aided catheter ablation of atrial fibrillation accurately, multi-modal information integration with EAM data and MRI/CT images is considered in Chapter 3. Specifically, we propose a variational formu- lation for surface reconstruction and incorporate the prior shape knowledge, which results 10 in a level set method. The proposed method enables simultaneous reconstruction and reg- istration under nonlinear deformation. Preliminary experimental results show the potential of the proposed approach. Geometric Feature-based Multi-modal Image Registration of Contrast-enhanced Cardiac CT with Gated Myocardial Perfusion SPECT CT and SPECT provide clinically complementary information in the diagnosis of CAD. Fused anatomical and physiological data acquired sequentially on separate scanners can be co-registered to accurately diagnose CAD in specific coronary vessels. A fully automated registration method is presented utilizing geometric features from a reliable segmentation of gated MPS volumes, where regions of myocardium and blood pools are extracted and used as an anatomical mask to de-emphasize the inhomogeneities of intensity distribution caused by perfusion defects and physiological variations. A multi-resolution approach is employed to represent coarse-to-fine details of both volumes. The extracted voxels from each level are aligned using a similarity measure with a piecewise constant image model and optimized using a gradient descent method. We then perform limited nonlinear regis- tration of gated MPS to adjust for phase differences by automatic cardiac phase matching between CT and MPS. For phase matching, we incorporate nonlinear registration using thin-plate-spline-based warping. Nonlinear Ultrasound Image Registration Based on Intensity and Local Phase Information A nonlinear ultrasound image registration method is proposed using the intensity as well as the local phase information under a variational framework. One application of this tech- nique is to register two consecutive images in an ultrasound image sequence. Although 11 intensity is the most widely used feature in traditional ultrasound image registration algo- rithms, speckle noise and lower image resolution make the registration process difficult. By integrating the intensity and the local phase information, we can find and track the nonlinear transformation of each pixel under diffeomorphism between the source and tar- get images. Experiments using synthetic and cardiac images of in vivo mice and human subjects are conducted to demonstrate the advantages of the proposed method. Nonlinear registration of serial coronary CT angiography for plaque development moni- toring CCTA is a high-resolution 3D imaging technique for the evaluation of coronary arteries. Registration of serial CCTA scans can aid in recognition of subtle changes and monitoring of coronary plaque development. We aim to develop a fully automatic and accurate nonlin- ear volume registration for longitudinal CCTA scan pairs in Chapter 5. Our algorithm com- bines global displacement and local deformation using nonlinear volume co-registration with a volume-preserving constraint. Histogram matching of intensities between two serial scans is performed prior to nonlinear co-registration with dense nonparametric diffeomor- phism in which sum of squared difference is used as a similarity measure. The segmented coronary artery trees provide initial anatomical landmarks for the co-registration algorithm that help localize and emphasize the structure of interest. To avoid possible bias caused by incorrect segmentation, we convolve the Gaussian kernel with the segmented binary coronary tree mask and define an extended weighted region of interest. A multi-resolution approach is employed to represent coarse-to-fine details of both volumes. The energy func- tional is optimized using a gradient descent method. The proposed algorithm is applied to 12 data from 10 patients with serial CCTA scans obtained within 10.75.7 months from each other on a 64-slice CT scanner. It can register serial CCTA scans accurately and allow di- rect comparison of calcified and non-calcified atherosclerotic plaque change between two scans. 1.5 Outline of the Dissertation The rest of the dissertation is organized as follows. The left ventricle segmentation problem using a mutual shape constraint is examined in Chapter 2. Multi-modal data integration for computer- aided ablation of atrial fibrillation is studied in Chapter 3. Multi-modal image registration of CT with MPS is studied in Chapter 4. A nonlinear ultrasound image registration based on in- tensity and local phase information is presented in Chapter 5. In addition, nonlinear registration technique of serial coronary CT angiography in plaque development monitoring is presented in Chapter 6. Finally, research ideas along with two future research topics are given in Chapter 7. 13 Chapter 2 : Curve Evolution with A Dual Shape Similarity and Its Application to Segmentation of Left Ventricle 2.1 Introduction Cardiac image segmentation has been one of crucial prerequisites in a wide variety of cardiac applications, which however still remains one of the most active and challenging areas of research in cardiac image analysis. The objective of segmentation is to partition an image into a finite number of important regions according to anatomical or functional structures. By now, their scope encompasses not only delineating chambers but also other parts of heart. Cardiac MRI is widely recognized as the gold standard for volumetric analysis of the left ventricle (LV). Accurate segmentation of LV , particularly delineating the epi- and endocardium, is a fundamental and crucial prerequisite for the derivation of heart parameters such as wall mo- tion, end-diastolic volume (EDV), end-systolic volume (ESV) and ejection fraction (EF) [48]. In general, segmentation in medical images remains to be a challenging problem since not only should the unique imaging properties be captured but there also exist intrinsic problems such as 14 low contrast, occlusions, illumination variations and other unmodeled phenomena which would make robust and accurate segmentation difficult. Several general segmentation approaches have been proposed so far. One of the widely used segmentation methods is the edge-based segmentation method, which includes the use of edge operators and active contours called snakes [70, 92]. The snake technique is to evolve a curve defined in the image domain under the influence of internal and external forces. In classical active contour models, an edge detector is used to derive the external force. They detect objects with edges defined by gradients. This approach has several limitations. First, the underlying image may suffer from low contrast and noise and the segmentation accuracy can be affected accordingly. Second, it is difficult to deal with topological change where curve parameterization is required. On the other hand, the region-based approach adopts a more global partitioning method that integrates intensity over the whole image domain [21, 95]. The region-based level set method represents curves implicitly as the zero level set of a scalar function. It can deal with topological change easily. Also, the region-based approach is more robust than the edge-based approach when the contrast between the object and the background is low. Recently, the model-based approach [29, 111, 132] has gained a lot of attention as a solu- tion to image segmentation problems with incomplete image information. It incorporates prior knowledge to improve the segmentation robustness and accuracy caused by poor image contrast, noise and missing boundaries. The prior shape model can be obtained using either explicit or implicit contour representations. The parametric point distribution model [28] with landmarks is an example of the explicit representation. However, the explicit boundary representation has sev- eral drawbacks such as fixed topology or pairwise correspondence. To overcome them, the level set method is employed for implicit representation. Segmentation with prior information was 15 adopted in the level set method before [29,111,132]. The implicit contour representation has sev- eral advantages over the explicit one in that it does not depend on a specific parameterization or a shape dissimilarity measure. Thus, it can handle topological change more easily. In either explicit or implicit contour representation, principal component analysis (PCA) has been widely used to get the shape prior. For example, Leventon et al. [81] proposed a shape prior model to restrict the flow of the geodesic active contour, where the prior was derived by performing PCA on a col- lection of signed distance function of the training shape. Most previous work [100, 132] tackled ill-posedness due to noise or subtle boundary by incorporating the explicit statistical shape prior. Given a set of training shapes, the prior information can be incorporated in the segmentation pro- cess. Although the shape-based method has led to successful segmentation results in [100, 132], they demand a manual segmented training data set, which is labor intensive. It is worthwhile to point out that the detection of epicardium is more challenging than that of endocardium since the background contrast is lower than the blood pool contrast, resulting in a partially missing epicardium boundary. Besides, it has a heterogeneous characteristics due to various structures immediately adjacent to the heart. Here, we propose a segmentation algorithm by introducing a novel shape constraint and apply it to cardiac MRI segmentation of LV; namely, the implicit shape prior, without a training pro- cedure based on the assumption that the endocardium and epicardium shapes are anatomically similar. Specifically, we formulate this problem as a coupling level set segmentation. Cardiac MRI data are used in our experiment and quantitative tests of 10 subjects in end-diastolic (ED) and end-systolic (ES) phase demonstrate the accuracy and robustness of the proposed approach. The rest of this chapter is organized as follows. Variational formulation and bayesian infer- ence for segmentation is presented in Sec. 2.2 followed by proposed shape constraint in Sec. 2.3. 16 Overall energy formulation and its minimization are given in Sec. 2.4 and real datasets are tested to demonstrate the efficiency and accuracy of the proposed method in Sec. 2.5. Concluding re- marks and future work are given in Sec. 2.6. 2.2 Variational Formulation and Bayesian Inference for Coupled Segmentation In our problem we need to extract two structures simultaneously in a single image. We propose a multi-phase segmentation scheme where the boundary of most probable structures of interest are jointly estimated. We model the boundaries of the epi- and endocardium as the zero level sets of two signed distance functions,f 1 (x) andf 2 (x), wherefxjf 1 (x)= 0g represents the endocardium boundary andfxjf 2 (x)= 0g represents the epicardium boundary, respectively. The observed image I can then be divided into three disjoint regions as follows: 8 > > > > > > < > > > > > > : endocardium(blood pool)=W 1 =fx2Wjf 1 (x)> 0g; epicardium(myocardium)=W 2 =fx2Wjf 1 (x)< 0 andf 2 (x)> 0g; extracardiac background=W b =fx2Wjf 2 (x)< 0g The region-based segmentation method can be obtained from a Bayesian formulation of the image partition problem. Given image I, the optimal segmentation boundary can be obtained by maximizing the following posteriori probability: p(f 1 ;f 2 jI)µ p(Ijf 1 ;f 2 )p(f 1 ;f 2 ); (2.1) 17 where p(Ijf 1 ;f 2 ) is the likelihood function and p(f 1 ;f 2 ) is the prior probability of the segment- ing curves. Maximizing this conditional probability with image I for the segmenting functionsf 1 andf 2 is equivalent to minimizing its negative logarithm and this gives the following minimizing process: E(f 1 ;f 2 )=log p(Ijf 1 ;f 2 ) log p(f 1 ;f 2 ); (2.2) where the first term represents data fidelity and the second term represents a shape constraint prior, respectively. 2.3 Proposed Shape Constraint The level set representation provides topology free interface and allows deformation of surface without any specific parameterization. LetWR 2 denote the image domain and a shape SW can be represented by zero level set of a higher dimensional embedding functionf(x):W!R as given by: 8 > > > > > > < > > > > > > : S=fx2Wjf(x)= 0g; interior(S)=fx2Wjf(x)> 0g; exterior(S)=fx2Wjf(x)< 0g; wheref(x) is a signed distance function that impliesjÑfj = 1 almost everywhere. Given two shapes S 1 and S 2 represented by signed distance functions f 1 (x) and f 2 (x) re- spectively, let S 1 (t) :[0;1]!R 2 be a parameterized curve that is an explicit representation for the interface of f 1 (x) (e.g. circle in Figure 2.1 (a)) and S 2 (s) :[0;1]!R 2 be a parameterized curve that is an explicit representation for the interface of f 2 (x) (e.g. rectangle in Figure 2.1 (a)). Then, we can measure a distance between two shapes in different aspects: the distance 18 P1 P2 P3 10 20 30 40 50 60 70 80 90 100 30 35 40 45 50 55 60 65 70 75 arc−length distance 35 40 45 50 55 60 65 70 (a) Two shapes (b) Distance along the arc length (c)f 2 (x) -f 1 (x) Figure 2.1: Illustration of notion of distance between two shapes. (a) Examples of two shapes, a circle S 1 in red and a rectangle S 2 in blue, represented by f 1 and f 2 respectively. (b) Profiles of distances that are measured by the shortest Euclidean distance from S 1 to S 2 at each point on S 1 (red line) and the shortest Euclidean distance from S 2 to S 1 at each point on S 2 (blue line) along the curves parameterized by arc length. (c) the difference between two signed distance functions that represents the two shapes. The outlines of the two shapes are superimposed and the values on the shape outlines represent the shortest Euclidean distance from one shape to the other. d 1 (t) can be measured by the Euclidean distance from the point S 1 (t) on S 1 to the closest point S 2 (s) on S 2 , which implies that S 1 (t) corresponds to S 2 (s) in a Euclidean sense. In the same way, the distance d 2 (s) can be measured by the Euclidean distance from the point S 2 (s) on S 2 to the closest point S 1 (t 0 ) on S 1 , which implies that S 2 (s) corresponds to S 1 (t 0 ). Utilizing the com- putational advantage with level set functions allows to havefd 1 (t)g=d(f 1 (x))(f 2 (x)f 1 (x)) andfd 2 (s)g = d(f 2 (x))(f 2 (x)f 1 (x)). Note that S 1 (t) is not necessarily identical to S 1 (t 0 ) as shown in Figure 2.1 where P3 is a corresponding point on S 1 from P2 on S 2 , however P1 is a corresponding point on S 2 from P3 on S 1 in measuring d 2 and d 1 . This implies that a distance between S 1 and S 2 cannot be uniquely determined based on either d 1 (t) or d 2 (s) un- less two shapes are isotropic such as concentric circles. Therefore, we introduce a distance D(x)=(f 2 (x)f 1 (x))(1H(f 1 (x))H(f 2 (x)) that can be interpreted as a quasi-average distance which is obtained by interpolated value between d 1 and d 2 and use this distance D(x) in defining a shape constraint. 19 We propose a novel shape prior knowledge which allows variation of Euclidean distances as well as geometrical changes between two shapes in a statistical distribution sense. We propose that D(x) can be approximated by Distance Gaussian Distribution in the following form: p(Djm;s) = 1 p 2ps 2 e (f 2 (x)f 1 (x)m) 2 2s 2 (1 H(f 1 (x)))H(f 2 (x)) (2.3) At fixed f 1 (x) and f 2 (x), the closed form of optimal statistical parameters can be obtained from the derivation: 8 > > < > > : m = R W (f 2 (x)f 1 (x))(1H(f 1 (x))H(f 2 (x))dx R W (1H(f 1 (x))H(f 2 (x))dx s = R W (f 2 (x)f 1 (x)m) 2 (1H(f 1 (x))H(f 2 (x))dx R W (1H(f 1 (x))H(f 2 (x))dx (2.4) By fixing s, distance variation, eventually the variation of shape can be controlled in a sta- tistical manner. This means that a larges allows a large deformation between two shapes and a smalls restricts the geometrical deviation from one shape to the other. 2.3.1 Synthetic Example for the Shape Constraint To illustrate the role of our shape constraint, simple yet demonstrative synthetic examples have been employed as shown in Figure 2.2 where an example of one circular shape (inside) and another deformed circular shape with a rectangular corner (outside) is presented in top row and an example of one circular shape (inside) and its noisy version (outside) is presented in the bottom row. Assuming that the inner circle and outer circle have similar shape, frequent phenomena, unusual local deformation due to subtle boundary and noise, are simulated in this experiment. By unusual local deformation, we mean that the boundary is missing because it does not exhibit features that can help characterize the boundary compared to the adjacent region. 20 (a) Chan-Vese model (b)s = 2 (c)s = 4 (d)s = 6 (e) Chan-Vese model (f)s = 2 (g)s = 4 (h)s = 6 Figure 2.2: Examples of the shape constraint with different s. Outer structure of interest is circle in both examples but the shape suffers from fuzzy boundary (top row) and noise (bottom row). The first column shows segmentation results of multiphase Chan-Vese model. Our shape constraint can capture the desired outer shape by adjustings. Ass becomes bigger, the allowable distance between two shapes becomes larger. In the second column, the outer shape is restricted to the circle due to the small s whereas the variation of distance becomes bigger in last two columns. By setting different values for standard deviation s, we can control the degree of deviation from one shape to the other, which results in constraining the allowable space of the mutual shape deformation. A small s restricts the both shapes to be similar and a large s allows a high variation of the distance between two shapes. The segmentation results using Chan-Vese segmentation model [21] (e.g. (a) and (e) in Figure 2.2) are compared to the results with our shape constraint prior for differents. 2.4 Energy Formulation Following the Mumford-Shah energy formulation [95] and Chan-Vese region-based segmentation model [21], we have an image model I= c 1 c 1 \c 2 +c 2 (1c 1 )c 2 +c b (1c 1 )(1c 2 ) where 21 c 1 and c 2 are characteristic functions for the regions W 1 ;W 2 respectively. Assuming intensity values to be homogeneous in each region, the data fidelity term is defined from in the following form and the intensitiesfc 1 ;c 2 , and c b g are approximated by the average within each region fW 1 ;W 2 , andW b g, respectively. E data (f 1 ;f 2 )= Z W H(f 1 (x))jI c 1 j 2 dx+ Z W (1 H(f 1 (x)))H(f 2 (x))jI c 2 j 2 dx + Z W (1H(f 2 (x))jI c b j 2 dx (2.5) In the optimization, keepingf 1 andf 2 fixed, the average constant intensity, c 1 ;c 2 , and c b , can be updated as 8 > > > > > > < > > > > > > : c 1 = 1 jW 1 j R W IH(f 1 )dx; c 2 = 1 jW 2 j R W I(1 H(f 1 ))H(f 2 )dx c b = 1 jW b j R W I(1 H(f 2 ))dx In addition, we impose length regularization for each segmented boundary using the following energy functional: E reg = Lengthf(f 1 = 0)+(f 2 = 0)g= Z W jÑH(f 1 (x))jdx+ Z W jÑH(f 2 (x))jdx (2.6) The data fidelity term may be insufficient to extract the structures of interest in the presence of missing or subtle boundaries. Thus shape prior knowledge needs to be incorporated to deal with those obstacles and avoid local minima. Instead, we propose to use a shape constraint prior, implicit shape prior knowledge, between structures of interest without a need of training data. 22 Then the shape constraint energy can be defined using Eq. 2.3 as E shape (f 1 ;f 2 )= log p(Djm;s); (2.7) where m is updated at every iteration using the closed form formula described in Eq. 2.4 and s can be chosen as constant empirically based on the deformation properties of the regions of interest to segment. Finally, we define the external energy using the edge indicator function, which speeds up the curve evolution. g= 1 1+jÑG s Ij 2 ; (2.8) where G s is the Gaussian kernel with standard deviations. E ext (f 2 )= Z W gH(f 2 )dx (2.9) When the function g is constant 1, the external energy functional is the area of the region W f 2 =fxjf 2 (x) < 0). Thus the external energy can be considered as the weighted area ofW f 2 . 2.4.1 Energy Minimization The optimal solution to our problem should minimize the following energy: E total (f 1 ;f 2 )= E data (f 1 ;f 2 )+aE reg (f 1; f 2 )+bE shape (f 1; f 2 )+lE ext (f 2 ) (2.10) 23 Table 2.1: Comparison of the manual segmentation and proposed algorithm Manual Segmentation Proposed Segmentation Max EDV (ml) 140 44 139 41 222/216 ESV (ml) 69 45 68 49 159/170 LVEF (%) 55 16 55 19 68/74 wherea controls the smoothness of segmenting shapes,b weights the significance of the shape prior knowledge that specifies an allowable variation in distance between two shapes, andl is a coefficient of the external energy which helps speed up the curve evolution. A gradient descent approach is used as a minimization scheme and alternative minimization with respect tof 1 andf 2 is performed. We omit the detailed derivation of the associated Euler- Lagrange equations. 2.5 Experimental Results 80 100 120 140 160 180 200 220 -15 -10 -5 0 5 10 15 Average of manual and proposed algorithm (ml) Difference between manual and proposed (ml) 20 30 40 50 60 70 80 -10 -8 -6 -4 -2 0 2 4 6 8 Average of manual and proposed algorithm (%) Difference between manual and proposed (%) (a) EDV (b) LVEF Figure 2.3: Bland-Altman plots for EDV and LVEF The proposed algorithm is implemented and tested on cardiac MRI at ED and ES phases of 10 subjects, followed by quantitative comparison to expert manual segmentation previously 24 obtained. The images were acquired on a clinical 1.5T MRI scanner (Siemens Sonata, Erlangen, Germany). The imaging parameters were as follows: field of view = 350mm, matrix = 255 192 to 255 255, slice thickness = 8mm, slice gap = 2mm, TE/TR = 1.6/3.1 msec, pixel bandwidth = 930 kHz, flip angle = 60 and segment size = 5-9 lines depending on heart rate. We used two circles from the endocardium center as manual initialization. Figure 2.4: Experimental results with ED and ES phases of three subjects. Manual segmentation results of endocardium (green) and epicardium (yellow) are shown in first and third column and results of proposed algorithm are presented with endocardium (red) and epicardium (white) using same images in second and fourth column. 3D volume is presented in the rightmost column. Parameters were selected empirically same for all (a = 400,b = 30,l = 400, ands = 2). Our approach provides statistical modeling for the margin between two shapes based on the fact that the distance between epi- and endocardium is not fixed but variable, and the epicardium delineation often suffers from vague boundary. Assuming fixed distance between two shapes, curve evolving can be robust to vague boundary but it is not flexible enough to capture the actual boundary. On the other hand, given no constraint, curve evolution depends only on the data 25 fidelity term, which results in passing through the boundary. Thus we achieved balance between them by adopting statistical modeling for the margin between two shapes. The mean values of EDV , ESV and LVEF calculated from 10 studies by using manual seg- mentation and proposed algorithm are presented in Table 2.1 and Bland-Altman plot is provided in Figure 2.3. Segmentation results of 2D and 3D volume for three cases are also presented in Figure 2.4. 2.6 Conclusion and Discussion We have proposed a novel shape constraint, implicit shape prior knowledge, to jointly segment the epi- and endocardium of LV simultaneously. With this, curve evolution is proposed using Bayesian inference in a variational framework. This shape constraint has an alternative role of substituting the shape prior from training. It is defined on the domain between endo- and epicardium and it allows to evolve the contour using statistical modeling for the margin between two shapes. We achieved a robust and efficient segmentation algorithm that performed well on cardiac MRI data. Currently, the initialization involves manual interaction to specify the size of the circle from the central point. However, this procedure can be replaced by an automatic initialization scheme that approximates a fitting circle by analyzing line histograms in a future work. 26 Chapter 3 : Multi-modal Data Integration for Computer-Aided Ablation of Atrial Fibrillation 3.1 Introduction Current treatment of cardiac arrhythmias ranges from non-invasive strategies such as pharmaco- logical therapy, to minimally invasive techniques such as catheter-based ablation, and to open surgical techniques. While medical therapy can mitigate the occurrence of arrhythmias, these treatments may have significant side effects since most drugs used have some toxicity that is not suitable for long-term therapy. The catheter-based procedure is proven to be an effective method in treating patients with certain cardiac arrhythmias [146]. It is much less invasive and more established. It also demands shorter recovery time than the surgical approach. Thus, catheter- based radio frequency (RF) ablation has become a widely accepted method in the treatment of cardiac arrhythmias, including Atrial Fibrillation (AF) and Ventricular Tachycardia (VT). These arrhythmias affect a large number of people and result in significant morbidity and mortality. AF is the most common sustained cardiac arrhythmia encountered in clinical practice. In the United States alone, there are over 3.5 million patients with this disorder [104]. AF can 27 result in serious complications, including congestive heart failure and thromboembolism. Despite recent advances, drug therapy to control this disease is still unsatisfactory. As an alternative, a non-pharmacological, interventional approach based on creating percutaneous catheter-based lesion inside the heart has been developed. Lesions are delivered in the left atrial-pulmonary vein junction with an aim to electrically isolate these veins from the rest of the atrium. This protects the atrium from fast heart beating that is originated in the veins, which initiate and perpetuate AF. The procedure of interventional AF treatment entails mapping the left atrium and the attached pulmonary veins using an electroanatomic mapping (EAM) system. This mapping information can be used to deliver lesions as well. This electrical approach is suitable for the heart since it is an electromechanical organ, where mechanical contractions are driven by electrical stimulus. However, there is a serious limitation of the EAM system in that it is not able to provide an accurate anatomical information of heart. Typically, a virtual shell is used to represent the atrial wall and the vein. The points on the atrial wall, where the catheter is manually touched, are used to create this shell [98]. The catheter-based ablation process can be greatly improved if a real anatomy is used instead of the virtual shell. To ensure safe catheter maneuverability and enable delivery of effective le- sions with minimal collateral damage and complications, it is critical to have both the anatomical information and the electrical information available to the operator. This is particularly impor- tant for performing ablation in a complex structure such as the left atrium that is surrounded by important organs, which are vulnerable to damage if lesions are not appropriately directed with close anatomical guidance. Furthermore, even the pulmonary veins themselves are liable to be damaged with grave long-term consequences if the lesions extend deeply into the veins instead of being restricted to the ostia. 28 The use of a multi-modal data integration process can provide an anatomical, physiological and functional representation simultaneously. In practice, this can be achieved by combining an anatomical surface model acquired by MRI/CT images and the localized electrical information measured by an EAM system. In the registration process, one obvious difficulty stems from the noise and/or outliers that are inevitably associated with the MRI/CT imaging process and the EAM data collection procedure. Unlike other organs in the body, heart undergoes contractile motion, apart from respiratory motion, thus making it unique and very challenging to register and integrate data of different modalities. In addition to physiological variations such as changes in the heart rate, the heart rhythm and the respiratory effect, various types of heart motion are the source of outliers. The main contribution of this chapter is to provide enhanced imaging of the anatomical heart surface from sparse and noisy EAM data in combination with a heart shape model obtained from MRI/CT reconstruction as a prior knowledge. For 3D surface reconstruction, we propose a variational formulation in the level set framework that is an efficient numerical scheme. The level set method is particularly of great use in representing a shape due to its topology free and implicit characteristics. By leveraging the 3D heart shape model, we can compensate incomplete EAM data, thereby representing the anatomical heart more accurately. The proposed method has two important advantages as follows. First, the proposed method is robust against non-rigid deformation caused by cardiac motion and noise and thus it can provide more accurate result. Second, it can construct the optimal surface without an explicit correspondence between the MR/CT surface and EAM data due to the implicit surface representation. The rest of this chapter is organized as follows. Previous related work is reviewed in Sec. 3.2 followed by the proposed multi-modal data integration method for computer-aided ablation of 29 atrial fibrillation in Sec. 3.3. Both synthetic and real data are tested to demonstrate the efficiency of the proposed method in Sec. 3.4. Concluding remarks and future work are given in Sec. 3.5. 3.2 Review of Related Work Developing a computer-guided system for ablative heart surgery involves image registration or integration techniques. In practice, this system requires fast converging algorithms for planning, monitoring due to time constraint imposed by real-time nature of cardiac intervention application. Consequently, they are usually performed under a rigid transformation between pre-operative MRI/CT reconstruction and intra-operative EAM data points [34, 105]. Among various registra- tion algorithms, the Iterative Closest Point (ICP) method and its variants have been widely used for this application due to their computational efficiency and rubustness [88]. The ICP algorithm begins with two meshes and an initial guess for their relative rigid-body transform. It refines the transform iteratively by generating pairs of corresponding points on the meshes and minimizing an error metric repeatedly [10]. However, the standard ICP algorithm does not take noise and outliers into account. Since noise and outliers may affect the ICP perfor- mance substantially, several ICP variants have been proposed in [112] to mitigate this problem. One popular approach to identify outliers is to use a threshold, including a certain constant, a frac- tion of a sorted distance and some multiple of the standard deviation of a distance [23, 89, 103]. Even with these variants, however, it is still challenging to deal with non-rigid deformation and differentiate inliers from outliers. Most of previous schemes used an ICP-based method without addressing the above-mentioned problem. Instead, they focused on clinical registration due to real-time nature of this application. 30 For example, Reddy et al. [88,105] showed the feasibility of combining MRI with CARTO-XP in a porcine model of myocardial infarction (MI). They used the mICP (modified Iterative Closest Point) scheme for registration, but did not address the outlier problem. The modification is to adopt hierarchical registration by adding the class information in the algorithm. A clinical regis- tration strategy that combines landmarks and surface registration was proposed in [34]. This study assessed the accuracy for each cardiac chamber using a different clinical registration method. It was observed in [62] that the size of the left atrium affects the accuracy. The patient who has a bigger chamber volume tends to have more ablation errors. The rigid transformation assumption made by existing schemes is simple yet insufficient in most cases. It often yields unsatisfactory results since a non-rigid deformation is involved be- tween the anatomical heart model reconstructed by MRI/CT images and temporal instances of the heart at the collection of EAM data points. This physiological and anatomical variation that occurs in the formation of the heart surface model and the collection of EAM data points de- mands a non-rigid transformation (or equivalently diffeomorphism) between the model and the data. Woo et al. proposed a novel image integration technique by incorporating non-rigid defor- mation using the level set method in [136]. To overcome the limitation of the traditional registration approach based on the rigid-transformation assumption, we formulate this problem as a 3D surface reconstruction problem from EAM data points with a given surface prior. A similar context arises in surface reconstruction from point clouds in a scanned noisy image. Surface reconstruction using an explicit representation has been considered by researchers, e.g., [26, 107]. Typically, this approach needs to parameterize a large point set that could be difficult to manipulate. Another approach was proposed in [5, 36] to construct triangulated surfaces using Delaunay triangulations and V oronoi diagrams. It has 31 to determine the right connection among points in the point set, which could be challenging in handling noisy and unorganized point data. Surface reconstruction based on an implicit shape representation using the level set tech- nique has been studied for almost two decades by applied mathematicians, e.g., [117]. The non-parametric (or implicit) surface representation has an advantage in dealing with arbitrary topology change and deformation. Hoppe et al. [67] proposed an algorithm in reconstructing a surface using the signed distance function from unorganized points. Zhao et al. [143] proposed another algorithm using the unsigned distance function and the weighted minimal energy to re- construct the surface. These algorithms are however restricted to situations where the population of points is dense enough to characterize the target surface. They are not applicable to our appli- cation where EAM data points are sparse and insufficient. The problem of insufficient EAM data encountered in the 3D heart surface construction (including the left Atrium and its pulmonary veins) can be mitigated by incorporating a heart shape prior provided by MRI/CT imaging. Then, the optimal surface can be obtained by minimizing the energy functional that consists of a data fitting term and a prior knowledge term as detailed in the next section. 3.3 Multi-modal Data Integration Model for Simultaneous Surface Reconstruction and Registration In this section, we present a multi-modal data integration algorithm for simultaneous surface reconstruction and registration. This algorithm reconstructs the heart surface from measured EAM data points and a heart shape prior obtained by MRI/CT imaging using the level set method. 32 3.3.1 Surface Reconstruction and Registration Under the level set framework [117], a surface S(x) inR 3 (a curve inR 2 ) to reconstruct can be represented by the zero level set of a higher dimensional embedding function f(x) :W!R as given by: 8 > > > > > > > > > < > > > > > > > > > : S=fx2Wjf(x)= 0g; interior(S)=fx2Wjf(x)> 0g; exterior(S)=fx2Wjf(x)< 0g; (3.1) whereW is the domain off(x). This implicit representation can provide the geometric shape of surface with the region defined by the union of S and its interior, denoted by ¯ S. Then, the shape of S is given by ¯ S(x)= H(f(x)) where H(x) is the heaviside function as defined by: H(x)= 8 > > < > > : 1; x 0 0; x< 0 (3.2) For the representation of a surface model M(x) obtained from MRI/CT images, we define its shape ¯ M(x) in a similar way using another level set function psi(x) as given by H(y(x)). Now We denote the set of measured EAM data points by D=fp 1 ; p 2 ;:::; p n gR 3 ; (3.3) where n is the number of data points. The essential assumption in this surface reconstruction application is that surface S to recon- struct is an equivalent class to a given prior surface model M under small non-rigid deformation 33 and the surface is close to the data points in D. Thus, surface S can be obtained by minimizing an energy functional that consists of a data fitting term and a prior knowledge term. Our goal is to find the embedding functionf associated with surface S that minimizes the following energy functional: E(f)= E point (f;D)+aE prior (f): (3.4) where detail of each term will be described in the following section. 3.3.2 Derivation of Energy Terms The energy functional in (3.4) consists of two terms. The first term measures how well the surface is fit to measured points based on the distance between the surface and these points. The second term measures how plausible the surface is in terms of the prior knowledge of the target surface. Parametera 0 is a weight that adjusts the importance of these two factors. The data fitting term can be written as E point (fjD)= n å i=1 Z W jf(x)d(x p i )j 2 dx; (3.5) where d(x)= d dx H(x) is Dirac measure and this term sums up the Euclidean distance between point p i and f. Recall that f is a signed distance function where the value of each point gives the euclidean distance between the point and the interface. To impose the prior knowledge on the target surface, S should be close to prior surface model M with a smooth surface assumption. In 34 other words, we penalize any abrupt change of the surface gradient. Here, we use two terms to represent the prior knowledge; i.e., E prior (f)= E reg (f)+ E shape (fjy); (3.6) where E reg (f) is the smoothness regularization term in the form of E reg (f)= Z W jÑH(f(x))jdx: (3.7) and E shape (fjy) is the shape dissimilarity term of the following form E shape (fjy)= Z W jH(f(x))H(y(T(x)))j 2 dx; (3.8) and where T(x) is a rigid transformation resulting from scaling, rotation and translation. Note that the smoothness regularization term measures the length for the 2D curve and the area of for a 3D surface while the shape dissimilarity term measures the symmetric difference between H(f) and H(y) under a rigid transformation. Combining (3.4)-(3.8), the total energy functional E(f) is given by: E(f)= E point (fjD)+a (E reg (f)+ E shape (fjy)): (3.9) Then, surface reconstruction can be achieved using the energy minimization principle as f = argmin f E(f) (3.10) 35 under constraintjÑfj= 1, which is a property of a signed distance function. Instead of applying the same weight a for both E reg and E shape as shown in (3.9), different weights can be used for each term. 3.3.3 Numerical Implementation For numerical implementation, we use the following approximations for the heaviside function and the Dirac delta measure as defined in [21, 142]: d(z)= 8 > > < > > : 0; ifjzj>e 1 2e h 1+ cos( pz e ) i ; ifjzje (3.11) and H(z)= 8 > > > > > > > < > > > > > > > : 1; ifz>e 0; if z<e 1 2 1+ z e + 1 p sin( pz e ) ; ifjzje (3.12) The energy functional in (3.10) can be minimized with respect to f(x) using the Euler- Lagrange equation with f(t = 0;x) = f 0 (x) defining the initial surface. Finally, the gradient descent method is applied to the resultant Euler-Lagrange equation, which leads to ¶f ¶t =2 n å i=1 f(x)d(x p i )+a[d(f)div( Ñf jÑfj ) 2 Z W (H(f) H(y(T(x))))d(f)dx]; f(t= 0;x)=f 0 (x) in W; d(f) jÑfj ¶f ¶~ n = 0 on¶W (3.13) where ! n denotes the exterior normal to the boundary¶W, and¶W=¶ ! n denotes the normal derivative off at the boundary. 36 To discretize the equation in f, we adopt a finite differences explicit scheme. The usual notations are as follows: let h be the space step andDt be the time step. The finite differences are expressed as D x f= f i+1; j;k f i1; j;k 2Dx ; D y f = f i; j+1;k f i; j1;k 2Dy ; D z f = f i; j;k+1 f i; j;k1 2Dz ; D xx f = f i+1; j;k +f i1; j;k 2f i; j;k Dx 2 ; D yy f = f i; j+1;k +f i; j1;k 2f i; j;k Dy 2 ; D zz f = f i; j;k+1 +f i; j;k1 2f i; j;k Dz 2 ; D xy f = f i+1; j+1;k f i+1; j1;k f i1; j+1;k +f i1; j1;k 4DxDy ; D xz f = f i+1; j;k+1 f i+1; j;k1 f i1; j;k+1 +f i1; j;k1 4DxDz ; D yz f = f i; j+1;k+1 f i; j+1;k1 f i; j1;k+1 +f i; j1;k1 4DyDz We first setf n as the initial surface and then updatef n+1 using following discretization. f n+1 i; j;k f n i; j;k Dt =2 n å m=1 f n i; j;k d(x p m )+ a h 2 d(f n i; j;k )k 2 Z W (H(f n i; j;k ) H(y(T(x))))d(f n i; j;k )dx; (3.14) where k =(D x f 2 D yy f 2D x fD y fD xy f+D y f 2 D xx f +D x f 2 D zz f 2D x fD z fD xz f+D z f 2 D xx f +D y f 2 D zz f 2D y fD z fD yz f): To solve the above partial differential equations numerically is challenging since the time step should be constrained to a small value in maintaining numerical stability. In addition, it 37 is computationally expensive to find a high dimensional surface. Thus, it is desired to employ an efficient numerical scheme and we naturally use multi-grid method that adopts a hierarchical representation of the data in multiple scales and propagates the solution from the coarse scale to the fine scale to achieve computational efficiency. 3.3.4 Relationship between Variational Formulation and Bayesian Inference This variational approach presented in Secs. 3.3.1 and 3.3.2 can be interpreted from the interpre- tation of Bayesian inference under a probabilistic framework. This relationship is presented in this subsection. The target surface S can obtained by maximizing the following posterior proba- bility: P(SjD)= P(DjS)P(S) P(D) (3.15) where P(DjS) is the likelihood function, P(S) is the prior probability of the surface. Maximiz- ing this conditional probability with data points D for surface S is equivalent to minimizing its negative logarithm: log(P(SjD))=log(P(DjS)) log(P(S))+ c; (3.16) where c is a constant. Thus, by setting E point (f;D)=log(P(DjS)); and aE prior (f)=log(P(S)); 38 we can convert Eq. (3.16) to Eq. (3.4). 3.4 Results and Discussion We begin with simple yet illustrative examples to demonstrate the efficiency and robustness of the proposed algorithm. We use synthetic geometric objects in 2D and 3D which have geomet- ric features that aim to be preserved under reconstruction process. Then, the evaluation of the algorithm is performed based on real patient data. 3.4.1 Synthetic data (a) Original shape (b) ICP (2% noise) (c) ICP (6% noise) (d) ICP (10% noise) (e) Deformed shape (f) Proposed (2% noise) (g) Proposed (6% noise) (h) Proposed (10% noise) Figure 3.1: The shape reconstruction results for a synthetic 2D star shape, (a) the original shape, (b)-(d) ICP results using different Gaussian noise levels, (e) the deformed shape, and (f)-(h) results of the proposed method using different Gaussian noise levels. We first compare the proposed scheme with the ICP scheme in registration accuracy using a synthetic 2D star shape as shown in Figure 3.1. The original 2D star shape (image size: 200 200 pixels) is shown in Figure 3.1(a). It is deformed as shown in Figure 3.1(e). Furthermore, 39 33 noisy contour points are generated by adding Gaussian noise (2%, 6% and 10% standard deviation of contour points, respectively) to original points obtained by sampling the original shape in Figure 3.1(a). 2% noise 6% noise 10% noise 0 1 2 3 4 5 6 7 Noise Level Mean Distance ICP Our Method Figure 3.2: Comparison of mean distances between the reconstructed surface and measured data points for the 2D star example at different noise levels using ICP and the proposed algorithm. The deformed shape stands for the reconstructed heart shape and data points with different noise levels represent the EAM data points. EAM data points can have errors from sensor, heart movement, and patient breathing. In synthetic experiments we used Gaussian noise to represent noises introduced in EAM data points. Graphical illustration of the registration results between deformed shape of different noise level and noisy contour points are presented in Figure 3.1(b)-(d) for the ICP scheme. In Fig- ure 3.1(f)-(h), shape reconstruction using both shape prior Figure 3.1(e) and noisy contour points 40 are presented using the proposed scheme. It is clear from Figure 3.2 that the proposed scheme outperforms ICP. For quantitative error analysis, we measure the mean Euclidean distance be- tween the reconstructed surface and measured data points with varying degree of Gaussian noise. The result is shown in Figure 3.2. Again, the proposed algorithm outperforms ICP significantly. This is especially true when the noise level is higher. The proposed algorithm produces more stable mean Euclidean distance than ICP as well. Next, we compare the proposed scheme with ICP using a 3D synthetic jar example. Experi- mental results are shown in Figure 3.3, where the original and the deformed shapes are shown in (a) and (f), respectively. Points extracted from the corrupted surfaces with various Gaussian noise levels (3%, 6%, 9% and 12%) are used for visual evaluation. Registration between deformed surface and data points with different noise levels using ICP is shown in Figure 3.3 (b)-(e). Re- constructed surfaces based on data points at different noise levels with the proposed algorithm are presented in Figure 3.3 (g)-(j). The proposed method generated the reconstructed surface which incorporates the data point as well as deformed prior shape in Figure 3.3 (e). The mean distances are also measured for accuracy comparison as shown in Figure 3.4. Again, the proposed algo- rithm is significantly better than the ICP scheme in terms of stability and the distance, which is especially obvious at higher noise levels, as shown in Figure 3.4. We can adjust the importance of both a data fitting term and a prior knowledge term that includes a regularization and shape similarity term by tuning the parametera. For the choice of the parametera, we usea= 1 3 that gives reasonable results. In general, ifa is too small, then the importance of the data fitting term increases, which results in overfitting to the data point. On the other hand, too biga produces the reconstructed shape which is almost same as the prior shape. Thus parameter is desired to be carefully chosen according to the application. 41 (a) Original shape (b) ICP (3% noise) (c) ICP (6% noise) (d) ICP (9% noise) (e) ICP (12% noise) (f) Deformed shape (g) Proposed (3% noise) (h) Proposed (6% noise) (i) Proposed (9% noise) (j) Proposed (12% noise) Figure 3.3: The shape reconstruction results for a synthetic 3D image: (a) the original shape, (b)- (e) ICP results using different Gaussian noise levels, (f) the deformed shape, and (g)-(j) results of the proposed method using different Gaussian noise levels. 3.4.2 Patient data The final example is a set of real patient data. 3D pre-operative contrast-enhanced MR angiogra- phy (MRA) was performed to delineate endocardial boundaries of the left atrium and pulmonary veins. The voxel size was 0.78125 0.78125 1.5mm and 45 slices were used in the exper- iment. We obtained MRA and 250 EAM data points from the same patient. The EAM data consists of the CARTO points imported from the CARTO-XP, including measurement points as well as ablation points. After delineating and removing unwanted regions such as the left ventricle (LV) and other small veins, we reconstruct the 3D model as shown in Figure 3.5 using ITK-SNAP [140] and Matlab software. Afterwards, a two-step registration process are applied for the ICP scheme. First, we perform the landmark registration using three junctions between LA and pulmonary 42 3% noise 6% noise 9% noise 12% noise 0 0.2 0.4 0.6 0.8 1 1.2 Noise Level Mean Distance ICP Our Method Figure 3.4: Comparison of mean distances between the reconstructed surface and measured data points for the 3D jar example at different noise levels using ICP and the proposed algorithm. veins: LA-LIPV , LA-LSPV , and LA-RSPV . These points are used for the initial pose of sub- sequent registration. Second, surface registration using the ICP scheme is performed to refine accuracy furthermore. The resulting image is shown in Figure 3.6(a). To validate the proposed algorithm, the optimal surface is reconstructed using 250 EAM data points by incorporating a heart shape prior from pre-operative MRA. By minimizing the energy functional, the final result is shown in Figure 3.6(b), where diamonds (in blue) represent EAM data points and circles (in red) represent ablation points. Blurred points are located inside. A quantitative evaluation result can be obtained in terms of the mean Euclidean distances of EAM and ablation points from the surface of the left atrium. They are reported in Table 3.1, which shows that the proposed approach outperforms the ICP method significantly. 43 (a) MRA of LA (b) 3D Reconstruction Figure 3.5: The 3D patient data: (a) MRA of LV and (b) 3D reconstruction result of the given MRA. (a) ICP (b) Proposed method Figure 3.6: The surface registration results for the patient data. 3.5 Conclusion and Future Work A novel multi-modal data integration technique using the level set method for catheter ablation of AF was presented in this chapter. This technique enables reconstruction and registration si- multaneously using data fitting, regularization, and shape prior energy terms. It provides better performance than the existing ICP method in accuracy. In the proposed framework, the heart shape model from MRA reconstruction is used as a prior shape knowledge. Thus, we can use the shape information to compensate for insufficient EAM data. Clinically, this technique can improve efficacy and safety of AF ablation by integrating EAM data and 3D imaging data. 44 Table 3.1: Performance comparison of ICP and the proposed method ICP Proposed EAM point mean distance 4.5087mm 2.4113mm Ablation point mean distance 3.2046mm 2.0921mm Dynamic cardiac shape analysis will make the current integration method more precise and meaningful. We plan to incorporate a richer set of spatio-temporal shape models using dynamic shape information in the future. Besides, we may consider a localized regularization method around the point data to obtain more precise reconstruction. 45 Chapter 4 : Geometric Feature-based Multi-modal Image Registration of Contrast-enhanced Cardiac CT with Gated Myocardial Perfusion SPECT 4.1 Introduction The common use of multiple imaging techniques in the same patient poses a great demand for multi-modal image registration to show functional and anatomical information on a single in- teractive fused display. Integration of complementary information provides the foundation for a variety of clinical applications [122, 123], including localization and boundary definition of organs and lesions, comparison of anatomical information with function, planning of radiation therapy and image-guided surgery. Computed tomography (CT) and myocardial perfusion single photon emission computed to- mography (MPS) convey different characteristics. This is due to the fact that the CT image is created by the attenuation of external photons throughout the thoracic volume, while MPS is created by emissions of photons primarily in the myocardium, with poorly perfused regions ex- hibiting a reduced photon signal. The advent of multi-slice CT scanners with rapid rotation has 46 increased the availability of noninvasive coronary CT angiography (CTA), a promising method for cardiac imaging. Contrast-enhanced coronary CTA allows depiction of coronary anatomy and non-invasive assessment of coronary artery stenosis and coronary artery plaque [2]. Recently, it has been suggested that contrast-enhanced cardiac CT can also be used to detect ischemic heart disease by assessment of rest and stress myocardial perfusion during CTA, but this method of perfusion assessment is not yet in common use [11]. Additionally, MPS is a mainstream imag- ing technique used to estimate myocardial hypoperfusion due to coronary stenosis [9], which however cannot be observed by MPS directly. Consequently, it is essential to combine different imaging modalities that provide clinically complementary information in the diagnosis, which can be achieved by registration. (a) (b) (c) (d) Figure 4.1: Examples of CT (top) and MPS (bottom) volume cross-sections illustrating image variations between patients: (a) and (b) show coronal views for two different patients. (c) and (d) show the axial/transverse views for the respective patients. It has been proposed that the combination of CTA and MPS modalities can provide clini- cally complementary information and offer better understanding of perfusion defects [49, 50, 96, 47 137], and most applications to date utilize manual or semi-automated image fusion for this pur- pose [49]. For example, visual analysis of fused MPS and coronary CTA images can improve the diagnostic value of sequential combined imaging [49], and manual tools for coronary CTA-MPS fusion have been developed [50]. However, manual alignment is tedious and observer-dependent, which complicates clinical protocols and reduces the practicality of such tools. Nakajo et al. [96] proposed a semi-automated co-registration of coronary CTA with gated MPS using a Left Ven- tricle (LV) binary model of coronary CTA. In their work, gated SPECT data were registered and mapped to a left ventricle binary model extracted from CTA data using manual, rigid and affine registration methods with a cost functional of the SPECT count value. However, manual pro- cedures were involved and the potential phase mismatch problem was not addressed although gated SPECT data were used for registration. On the other hand, our group was the first to pro- pose a preliminary automated multi-modal combination of MPS with coronary CTA scans using pre-segmented MPS volumes. A preliminary approach from our group [137] has dealt with the phase mismatch problem by assuming that reconstructed CT cardiac phases were obtained from the ED phase and MPS data could be integrated to the ED phase using a “motion-frozen” tech- nique [125]. However, one limitation of this approach was the rigid body assumption, which may be significantly inaccurate when reconstructed cardiac phases differ in the two modalities. In general, there have been other related studies on the registration of MPS with other cardiac modalities. Faber et al. [39] proposed a point-based registration of the MPS surface with 3D coronary anatomy reconstructed from invasive 2D coronary angiography, and it did not use vol- ume registration. Guetter et al. [56] developed a registration method for MPS and non-contrast CT data to improve attenuation correction. However, in their studies, images were obtained by the hybrid SPECT-CT scanner (in contrast with different stand-alone scanners) where images are 48 already in approximate alignment, requiring only small correction. Guetter et al. [56] proposed a learning-based method that exploits the mutual information and an intensity co-occurrence prior. However, this approach used an off-line training process. Besides, all data sets have to be pro- vided a priori; otherwise, the prior model may not be incrementally extended as new image data (e.g. different vendors or abnormal data) are available. In related development, Aladl et al. [4] proposed an MRI-MPS volume registration technique which utilized the MRI motion to pre- segment the heart for registration with MPS. However, this approach is not practical for CTA since typically only a few cardiac phases with acceptable image quality may be available for analysis due to radiation dose considerations. Registration methods can be categorized into two types.“Feature-based” methods rely on lo- cal image features such as edges or relevant anatomical structures, whereas “intensity-based” methods directly use image intensity and do not require any additional function for characteriz- ing structures in image. A number of different measures have been applied for the computation of dissimilarity. MI and its variant are the most popular approaches, especially for multi-modal image registration [86]. However, in multi-modal image registration problems, intensity alone is often insufficient, and additional information needs to be incorporated. This is especially true for images from abnormal subjects who have inherent physiological variations between differ- ent modalities. For this reason, we employed the regional feature- and intensity-based energy functional. This is mainly due to the fact that both MPS and coronary CTA have large variations caused by the different injection contrast and perfusion defects, resulting in highly variant in- tensity distribution, even in the same region. Also, the standard dissimilarity measure for multi- modal image registration such as mutual information [133] may not be very effective without specialized pre-processing. 49 In this chapter, we aim to develop a fast and fully automated multi-modal image registration method for aligning functional MPS with anatomical CT scans by extending our preliminary ap- proach [137]. Accurate registration of contrast-enhanced cardiac CT with MPS is challenging as different anatomical features are visible in different images. For example, although the my- ocardium is visible on both cardiac CT and MPS, there could be severe perfusion defects present on MPS images which do not have equivalent appearance on CT as shown in Fig. 4.1. Addi- tionally, variation in the intravenous contrast distribution for different patients causes intensity variations on CT images within the blood pool region. To alleviate the problems and provide robust registration results regardless of image varia- tions, we proposed to utilize the explicit and geometrically characteristic information (i.e. seg- mented boundary of MPS) to find the correspondence and estimate intensity values in the my- ocardium and blood pool regions simultaneously using a variational framework, where a piece- wise constant image model of each segmented region is incorporated as a similarity measure. Both segmentation and registration are interrelated [139] and heavily subject to image charac- teristics that provide discriminative power between the regions of interest and the background. However, a joint segmentation and registration scheme may not be applicable to our problem since the characteristic features that define the boundary of regions of interest in MPS are usu- ally subtle (especially in abnormal datasets) and this will severely degrade the performance of both segmentation and registration. Consequently, we have proposed an algorithm that applies segmentation as an initial step and subsequently registration procedure. We first perform rigid registration of the segmented boundary of gated MPS with cardiac CT. Subsequently, the best matching cardiac phase is automatically found based on the cost functional, followed by nonlinear landmark-based warping of gated MPS to the specific phase 50 matching to the CT phase. The nonlinear registration is applied within gated MPS for better visualization of defects to improve spatial resolution and to remove blurring due to motion ar- tifact. To our knowledge, our technique is the first to address the issues of multi-modal image registration and cardiac phase matching problem at the same time. In our application, we only use segmentation of MPS images, since cardiac CT has relatively homogeneous intensity distri- bution compared with MPS volume. Although the degree of contrast between the blood pool and myocardium of cardiac CT varies according to the dataset obtained, the proposed scheme can estimate corresponding intensity values in an alternative manner. 4.2 Materials and Methods 4.2.1 Description of the Method 4.2.1.1 Pre-processing A fully automated segmentation algorithm [53] for the left ventricle (LV) of gated SPECT pre- viously developed in our laboratory is adopted to define specific regions of interest, including myocardium, blood pool, and extracardiac region, in the MPS volume. In this step, all gated MPS volumes are segmented. This technique has been previously validated and provides high accuracy of MPS segmentation [51, 52]. In brief, for each interval of a gated MPS volume, an asymmetric Gaussian is fitted to each profile from which a maximal count midmyocardial surface is determined. The inflection points of the Gaussian are taken to be the surface points. Geometric features of boundary of left ventricle obtained from the segmentation of MPS vol- umes are of great importance as they can provide the anatomical information to guide registration. Those regions to be extracted from MPS volumes usually do not present distinguished boundary 51 (a) (b) (c) (d) Figure 4.2: An example of myocardial segmentation of MPS: (a) shows transverse views for an abnormal case (top) and a normal case (bottom); (b) shows respective contours from segmen- tation overlaid on the images from (a); (c) illustrates the blood pool (D 1 ), myocardial (D 2 ) and extracardiac (D 3 ) regions; (d) shows 3D contours obtained from the myocardial segmentation, endocardium (gray surface) and epicardium (orange mesh). Note the robust segmentation result despite abnormality in MPS. features, resulting in an inhomogeneous intensity distribution within the same anatomical region as shown in Fig. 4.2. Let I and J be segmented MPS and cardiac CT volumes, respectively, inW, which correspond to an open and bounded domain inR 3 . Both volumes have a different image size but the trans- formation is calculated on the same space. The segmented MPS volume, I, is the union of three disjoint regions; namely, the blood pool, myocardium, and extracardiac structures. We denote regions of blood pool, myocardium, and extracardiac structures by D 1 , D 2 , and D 3 respectively, whereW= D 1 [ D 2 [ D 3 as shown in Fig. 4.2 (c). We assume that CTA has relatively homoge- neous intensity values in blood pool and myocardium and apply a Gaussian-shaped filter of 3mm standard deviation to the CTA volume to reduce the intensity variation due to noise and lower the CTA image resolution closer to the MPS image. 52 Figure 4.3: Overview of the registration procedure 4.2.1.2 Registration Model Given pre-processed MPS and CT volumes, we automatically align the centers of the end-diastolic (ED) phase of segmented MPS and the reconstructed CT volume as an initialization step, and then apply the proposed multi-modal registration method. An overview of the registration method is presented in Fig. 4.3. Our goal in multi-modal registration is to estimate the transformation, T Q (x), that maps the MPS gated phase to the cardiac CT volume. The transformation T Q (x), is a rigid transformation of parametersQ=(t x ;t y ;t z ;q x ;q y ;q z ) that denote translations with respect to the x, y, and z directions and rotations in the x, y, and z directions, respectively. The optimal transformation ˆ T Q (x) was found by iteratively shifting and rotating the MPS until the energy functional defined in Eq. (4.1) computed with the two volumes is minimized. The energy functional based on the segmented MPS volume and CT volume can be defined as E(c 1 ;c 2 ;c 3 ;Q)= 3 å i=1 Z D i jc i G s (J T Q )(x)j 2 dx = 3 å i=1 Z W c D i (x)jc i G s (J T Q )(x)j 2 dx; (4.1) 53 where * denotes convolution operator, c 1 , c 2 and c 3 are constants that approximate the average intensities of regions D 1 , D 2 , and D 3 , respectively, G s is a Gaussian kernel with standard devia- tions, and T Q is a rigid transformation of 6 parametersQ and J T q is transformed CT volume. The characteristic functionc D is defined by c D (x)= 8 > > < > > : 1; x2 D 0; x = 2 D: (4.2) For both CT and MPS volumes, transverse image orientation is used for the registration pro- cess. 4.2.1.3 Minimization The optimal values for the parameters (c 1 ;c 2 ;c 3 ;Q) that minimize the given energy functional are obtained by solving the associated Euler-Lagrange equations, and a gradient descent method is performed. Parameterizing the descent direction by an artificial time t> 0, the derivation of the Euler-Lagrange equations is given below: ¶Q i ¶t = ¶E ¶Q i (i= 1;2;;6): (4.3) In the derivation, we consider only one term for a single region in the energy functional for simplicity, which is given by 1 2 ¶E ¶Q i = Z W fc D (x)(c G s (J T Q )(x))g ¶ ¶Q i (c G s (J T Q )(x) dx (4.4) 54 By setting T Q (x)= T(Q;x), we get ¶ ¶Q i (G s (J T Q )(x)) = Z R 2 G s (y) ¶ ¶Q i fJ(T(Q;x y))gdy = Z R 2 G s (y)ÑJ(T(Q;x y)) ¶T(Q;x y) ¶Q i dy: (4.5) Furthermore, by setting ¶JT Q (x) ¶Q i = L i (Q;x)=ÑJ(T(Q;x)) ¶T(Q;x) ¶Q i , we obtain ¶ ¶Q i G s (J T Q )(x) = Z R 2 G s (y)L i (Q;x y)dy= G s L i (Q;x): (4.6) Hence, Eq. (4.3) becomes 1 2 ¶E ¶Q i = Z W fc D (x)(c G s (J T Q )(x))gG s L i (Q;x)dx: (4.7) By keepingQ fixed, the minimization of the energy E with respect to constants c 1 , c 2 and c 3 leads to 8 > > > > > > > > > < > > > > > > > > > : c 1 = R W c D 1 (x)(G s (J T Q )(x))dx R W c D 1 dx ; c 2 = R W c D 2 (x)(G s (J T Q )(x))dx R W c D 2 dx ; c 3 = R W c D 3 (x)(G s (J T Q )(x))dx R W c D 3 dx : (4.8) The constant intensity values and transformation parameters can be minimized using Eqs. (4.7) and (4.8) in an alternative manner. We employ a coarse-to-fine multi-resolution scheme for computational efficiency and robustness. 55 4.2.1.4 Cardiac Phase Matching with Thin-Plate-Spline Warping (a) (b) (c) Figure 4.4: An example of finding the best matching phase of gated MPS by minimizing the registration cost functional: (a) plot of cost functional value across different phases of MPS; (b) overlaid registration result of phase 9 MPS on CT perfusion; (c) overlaid registration result of best matching phase (16) with CT perfusion. Note that the best matching MPS phase is visually well correlated with CT. It is necessary to match cardiac phases between MPS and CT, since the reconstructed CT cardiac phase may vary. Consequently there could be a mismatch between cardiac phases of two modalities, which might affect further visual analysis and image registration. To match the cardiac phase between these two modalities, we perform rigid registration of gated MPS to CT, using the ED phase of MPS (since cardiac CT is usually reconstructed in diastole), and then evaluate the cost functional value for all the other phases of gated MPS. Subsequently, the best matching phase of MPS can be selected according to the cost value where the minimum value is found: k = argmin k2f1;:::;ng E k (c 1 ;c 2 ;c 3 ;Q); (4.9) where k denotes phase number, n is total number of phase, and E k (c 1 ;c 2 ;c 3 ;Q) denotes the cost functional defined in Eq. (4.1) with phase k of MPS. In Fig. 4.4, we illustrate one example of car- diac phase matching, where phase 16 is the best matching phase of MPS (defined as the smallest 56 cost value). By using cardiac phase matching, it can be clearly seen that both the modalities are aligned more accurately as shown in Fig. 4.4 (c). In what follows, derived LV contours are used in combination with nonlinear registration to compensate for phase mismatch, thereby creating the average MPS volume, warped to the best matched cardiac phase that best corresponds to CT. For the nonlinear registration, 3D thin- plate-spline (TPS) [12] is used to model the nonrigid deformation between all image phases of gated MPS and the best matching cardiac phase position. This is similar to the “motion-frozen” technique [125] which is also previously developed by our group. This nonlinear registration step is applied only within gated MPS data for better visualization of defects. Once detection of epi- and endocardial surfaces on the gated MPS images is accomplished using the technique previously developed [53], the contour points are positioned on vectors nor- mal to the myocardial surface. These corresponding points are used as template and target control points to be registered where the control points are kept unchanged. Let a template shape surface denote P and a target shape surface be Q. We assume that P L , the landmarks in the template shape P, are known and fixed. Given n control points in a template volume ˆ p=( ˆ p x ; ˆ p y ; ˆ p z )2R 3 , the TPS finds a mapping f =( f x ; f y ; f z ) from P L to corresponding landmarks Q L , ˆ q=( ˆ q x ; ˆ q y ; ˆ q z )2R 3 , on the target surface Q. The thin-plate bending energy which measures the required energy to de- form a volume to match these two sets of landmarks is then characterized by d(P L ;Q L )= ZZZ ¥ ¥ (L( f x )+ L( f y )+ L( f z ))dxdydz; (4.10) where L()= ¶ 2 ¶x 2 2 + ¶ 2 ¶y 2 2 + ¶ 2 ¶z 2 2 + 2 ¶ 2 ¶x¶y 2 + 2 ¶ 2 ¶y¶z 2 + 2 ¶ 2 ¶z¶x 2 . The thin-plate bending energy is invariant to any affine transformation. The TPS model is a transformation 57 Figure 4.5: “Motion-frozen” technique: All MPS frames are registered to a common reference frame, which is the best matching MPS phase for the CT volume. TPS-based warping ( f i ) is used to deform and map each frame to I N . The control points from the segmentation contours serve as landmark points for TPS warping, and corresponding points from the source and target contours determine the mapping. Final “motion-frozen” image is obtained by averaging over all the warped volumes. Note that perfusion defects are clearly seen in the final image (final column). that is able to represent elastic deformation. Successive TPS-based transformations f 1 through f N1 are computed that map I 1 through I N1 , respectively, to I N which is the best matching phase corresponding to CT phase as shown in Fig. 4.5. Final “motion-frozen” volume (MFV) can be given as MFV = 1 N N1 å i=1 I i ( f i )+ I N ! (4.11) Note that perfusion defects are observed more clearly in Fig. 4.5 (final column) using the cardiac phase matching and TPS warping that are applied individually to each phase, for all the 58 cardiac phases. All the registered image frames are then integrated to represent a final “motion- frozen” volume in the cardiac phase, which corresponds to the CT phase. For example, for a 16-phase gated data set, 15 image phases are registered to the best matching phase. 4.2.2 Validation 4.2.2.1 Patient selection Between October 1, 2005 and May 31, 2007, we retrospectively identified 35 consecutive patients (26 men and 9 women; mean age 6712 years) who underwent myocardial MPS and coronary CTA within a 90-day period at Cedars-Sinai Medical Center (site A) or a neighboring outpatient imaging center (Cardiovascular Medical Group of Southern California, Beverly Hills, California) (site B). All had full datasets available for processing. Additionally, we also retrospectively iden- tified 5 patients who underwent MPS and stress CT angiography for assessment of myocardial perfusion within 30 days at site A. This study was approved by the Institutional Review Board at Cedars-Sinai Medical Cente 4.2.2.2 Coronary CTA: acquisition and reconstruction Coronary CTA was performed on the SOMATOM Definition dual-source CT (DSCT) scanner (Siemens Medical Systems, Forchheim, Germany) at site A and by a Lightspeed VCT 64-slice CT scanner (GE Healthcare, Waukesha, WI) at site B. Pre-scan beta-blockade was performed at both sites to attain a heart rate <70 beats/min (<60 beats/min at site B). At site A, ECG-gated helical CTA was performed during a 10-12 second breath-hold by injecting 92 ml of intravenous contrast and using scan parameters of: heart-rate dependent pitch (range 0.2 to 0.45), 330 ms gantry rotation time, 83 ms temporal resolution, 0.6 mm slice thickness, 120 kVp tube voltage 59 and 600 mAs tube current. At site B, ECG-gated helical CTA was performed during a 9-12 second breath-hold by injecting 100 ml of intravenous contrast and using scan parameters of: heart-rate dependent pitch, 350 ms gantry rotation time, 0.625 mm slice thickness, 120 kVp tube voltage, 300-700 mA tube current depending on patient size. At both sites, ECG-based dose modulation was used for all patients to limit radiation dose. The specifications of reconstruction at both sites are as follows: Site A. Retrospectively-gated reconstruction of raw CTA data was routinely performed at 40%, 65%, 70%, 75% and 80% of the R-R interval using the following parameters: 0.6 mm slice thickness (0.75 mm if BMI> 35 kg/m2), 0.3 mm slice increment, 250 mm field-of- view, 512x512 matrix, and B26f medium smooth kernel. Reconstructed data for the cardiac phase with the best quality were transferred to a separate Windows workstation via DICOM protocol for coronary CTA-SPECT fusion. Site B. Raw data from 5% phases of the cardiac cycle representing 40-80% of the R-R interval were reconstructed with 0.625 mm slice thickness, and 250 mm field- of-view, using single-segment reconstruction and a standard convolution kernel. The matrix size was 512x512 with a transverse pixel spacing of 0.49 mm/pixel. Reconstructed data for the 75% phase were transferred to a separate Windows workstation via DICOM protocol for coronary CTA-SPECT fusion. 4.2.2.3 Stress CT perfusion: acquisition and reconstruction Patients enrolled in the CT perfusion study were part of an IRB-approved pilot study at the Cedars-Sinai Medical Center (site A). Under this protocol, patients with perfusion defects in their MPS study were recruited for a stress CT perfusion study within 1 month of their MPS scan. Adenosine (Adenosine, Astellas Pharma, Illinois, US, Inc.) was infused over 5 min at a constant rate of 0.14 mg:min 1 kg 1 . At the end of 3 minutes of adenosine infusion, contrast-enhanced 60 CT was performed as described in II.B.2. From this contrast-enhanced stress CT scan, the best- quality cardiac phase was identified (either 70% or 40% for the patients in this study, depending on the heart-rate at the time of the scan). Stress CT perfusion series was retrospectively recon- structed from this best-quality cardiac phase with the following parameters: 512x512 matrix, 250 mm field-of-view 1 mm slice thickness, 0.5 mm slice increment and a smooth reconstruction ker- nel (B10f). The reconstructed stress CT perfusion data was transferred to a separate Windows workstation via DICOM protocol for CTA-SPECT fusion. 4.2.2.4 Myocardial perfusion SPECT protocol Rest stress MPS studies were acquired with dual isotope (Tl-Tc) protocol [8] or a low dose-high dose same day Tc protocol [60]. MPS acquisitions were performed with non-circular orbits, obtaining 64 projections over 180 degrees (45 RAO to 45 LPO). Acquisitions were performed on Philips CardioMD or Forte cameras or on Siemens e.cam camera. The stress was performed with exercise, adenosine injection or adeno-walk protocol [7]. No attenuation or scatter correction was used. Gated and summed data were reconstructed with filtered back projection (FBP) and with Butterworth filter (cut-off 0.83 cycles/cm, order 5) to original transverse orientation without any short axis reorientation. Site A had 16-bin gated acquisition and site B had 8-bin gated acquisition. The matrix and reconstructed voxel size of MPS are 64x64x24 and 6.4x6.4x6.4 mm, respectively. 4.2.2.5 Evaluation of automatic registration To evaluate the performance of the proposed method in terms of the registration accuracy and robustness, we have performed a series of experiments on datasets of 20 stress/rest MPS images 61 and coronary CTA datasets acquired from two different sites (as mentioned above), and addition- ally 5 stress MPS images and CT perfusion datasets. The proposed method was also compared to the conventional intensity-based normalized mutual information (NMI) method without any specialized pre-processing. Registration results using the proposed registration method and NMI based method were compared with manual alignment by two expert observers. Each image S : D7! [0;1] is associated with a probabilistic model by introducing a local random variable X, which is uniformly distributed in D, and its related intensity random variable S(X). For registration, MPS images were chosen as the reference images M : D M 7![0;1], and CT images were chosen as the floating images N : D N 7![0;1]. The transformation is represented by a mapping T Q : D N 7! D M . The joint intensity distribution of M(X) and N(T Q (X)) in the overlap region V = D N \ T 1 Q (D M ) can be estimated from pixel samples using Parzen window given by p(m;n;T)= 1 jVj Z V j m M(x) r j n N(x) r dx; (4.12) wherej is a Gaussian kernel andr determines the width of the window. The joint histogram and NMI can be defined to represent the global statistics between two images and have been widely employed for multi-modal image registration, and proven to be less sensitive than mutual information to changes in the region of overlap [128]. The joint histogram, H, and NMI, N, are defined as follows: H(M;N;T)= ZZ p(m;n;T)log p(m;n;T)dmdn (4.13) N(M;N;T)= H(M;T)+ H(N;T) H(M;N;T) (4.14) 62 For our experiment, we set the number of histogram as 10, kernel bandwidthr as 0.4 and num- ber of sample size as 10% of the total pixel number. Additionally 2 levels of multi-resolution optimization were performed for both NMI based method and the proposed method. As for con- vergence criteria, the algorithm stops when translation is less than 0.01 mm or iteration reaches the predefined iteration number 100 for both methods. The registration errors were measured by absolute difference of 6 parameters(t x ;t y ;t z ;q x ;q y ;q z ) between manual alignment and results obtained from the two methods. We tested the data sets under the same initialization conditions. For all the data sets, multiple registrations were done with different initializations (translations from -100 mm to 100 mm in each direction with step size 10 mm) away from the ground truth to measure the capture range. 4.3 Results The registrations were performed on an Intel Core2 Duo CPU with a clock speed of 2.5 GHz and 4 GB memory. The mean computation time for the whole process of the proposed method was 7.30.4 sec. The NMI based method had a larger error than the proposed technique and the average registration time was 14.88.6 sec. Our registration of MPS with cardiac CT achieved a success rate of 87% (6 failed out of 45 data sets). Examples of the coronary CTA and MPS pair and registration results of CT perfusion-MPS registration are shown in Fig. 4.6 and Fig. 4.7, respectively. Fig. 4.8 plots the translational probes for registering one representative example case of CT perfusion using the proposed method and NMI based method. Fig. 4.8 (a) plots the result using NMI based method and its NMI values with different offsets from the ground truth as visually 63 Figure 4.6: An example of registration results: original CT image (top), images before (middle) and after (bottom) registration are shown. Errors were (2, 1, 0) mm for translation and (1, 1, 0) degrees for rotation, as compared to visual alignment. determined are presented in (b) and Fig. 4.8 (c) plots the result using the proposed method with different offsets and its cost values are presented in (d). The translational and rotational registration results were not much different in inter-observer variability and in the stress and rest MPS (p=NS) as shown in Table 4.1. No significant differences were observed between the 2 different coronary CTA/MPS protocols, as shown in Fig. 4.9. That is, (4.33.7, 3.63.6, 3.62.1) mm for translation and (2.13.2, 0.40.9, 0.71.2) degrees for rotation for site A and (3.82.6, 3.53.4, 2.62.3) mm for translation and (1.22.0, 1.73.1, 64 Figure 4.7: An example of registration results (CT Perfusion): original CT image (top), images before (middle) and after (bottom) registration are shown. Errors were (1, 0, 4) mm for translation and (5, 2, 3) degrees for the rotation, as compared to visual alignment. 2.03.8) degrees for rotation for site B (p=NS). The results of CT perfusion were (3.02.9, 3.52.4, 2.81.0) mm for translation and (3.02.4, 0.60.9, 1.21.3) degrees for rotation. In Fig. 4.9 and Table 4.1, we have compared the different registration results with different data sets. The result shows that the proposed method performed better than the NMI based method in translation, and similarly in terms of rotation. In these results, we excluded the 6 failed cases. In Table 4.2, we present validation results of finding the best cardiac matching phase by the algorithm as compared with the visual ground truth obtained from an expert observer. It can 65 −80 −60 −40 −20 0 20 40 60 80 1.02 1.025 1.03 1.035 1.04 1.045 1.05 1.055 Misliagnment along each axis Cost value of NMI based method X axis Y axis Z axis (a) (b) −80 −60 −40 −20 0 20 40 60 80 2.095 2.1 2.105 2.11 2.115 2.12 2.125 2.13 x 10 5 Misliagnment along each axis Cost value of proposed cost functional X axis Y axis Z axis (c) (d) Figure 4.8: Comparison of NMI based method and proposed method: (a) registration result using NMI based method; (b) the NMI values for translational offsets from the ground truth for each axis; (c) registration result using the proposed method; (d) the cost values of proposed method for translational offsets from the ground truth for each axis. Note that the anatomical features are very different and the image segmentation will aid the registration of MPS to CT. Additionally cost value of the proposed method is optimal at the origin whereas NMI does not provide any such robust optimum. 66 (a) (b) (c) (d) (e) (f) Figure 4.9: Registration errors: (a) translational errors of stress/rest MPS, and CT obtained from site A; (b) rotational errors of stress/rest MPS, and CTA obtained from site A; (c) translational errors of stress/rest MPS, and CT obtained from site B; (d) rotational errors of stress/rest MPS, and CT obtained from site B; (e) translational errors of stress/rest MPS and perfusion CT; (f) rotational errors of stress/rest MPS and perfusion CT. Table 4.1: Registration error and inter-observer variability Translation (mm) Rotation (degree) X Y Z X Y Z NMI based Stress 9.58.2 8.86.7 8.56.2 1.82.7 1.12.3 1.63.0 Method Rest 7.35.6 10.68.6 12.76.5 2.03.8 1.12.5 3.11.5 Proposed Stress 4.63.4 4.03.0 2.92.1 1.72.6 0.92.3 1.43.1 Method Rest 4.33.2 3.82.4 3.82.4 1.83.3 0.10.3 0.30.7 Inter-observer Stress 3.72.4 3.52.8 2.21.8 3.53.1 1.93.0 3.12.8 Variability Rest 4.32.6 3.83.1 3.22.9 3.03.5 1.83.1 2.94.7 67 Table 4.2: Evaluation of best matching phase Pation No. Stress/Rest Visual ground truth Our method (phase/total phase) (phase/total phase) 1 Stress 5/8 7/8 Rest 6/8 7/8 2 Stress 5/8 7/8 Rest 6/8 7/8 3 Stress 5/8 7/8 Rest 7/8 8/8 4 Stress 6/8 7/8 Rest 6/8 7/8 5 Stress 8/8 8/8 Rest 1/8 1/8 6 Stress 7/8 7/8 Rest 1/8 7/8 7 Stress 5/8 6/8 Rest 5/8 7/8 8 Stress 8/8 7/8 Rest 8/8 7/8 9 Stress 14/16 15/16 Rest 16/16 16/16 10 Stress 1/16 15/16 Rest 15/16 16/16 11 Stress 14/16 14/16 Rest 12/16 11/16 12 Stress 2/16 1/16 Rest 2/16 1/16 13 Stress 2/16 1/16 Rest 2/16 1/16 14 Stress 14/16 16/16 Rest 15/16 16/16 15 Stress 1/16 1/16 Rest 2/16 1/16 16 Stress 16/16 16/16 Rest 16/16 1/16 17 Stress 16/16 1/16 Rest 14/16 16/16 18 Rest 6/16 6/16 19 Stress 1/16 15/16 20 Stress 1/16 15/16 21 Stress 1/16 15/16 22 Stress 8/16 8/16 68 −100 −80 −60 −40 −20 0 20 40 60 80 100 0 10 20 30 40 50 60 70 80 90 Misalignment along X axis Translation error (mm) X axis Y axis Z axis −100 −80 −60 −40 −20 0 20 40 60 80 100 0 10 20 30 40 50 60 70 Misalignment along Y axis Translation error (mm) X axis Y axis Z axis −100 −80 −60 −40 −20 0 20 40 60 80 100 0 10 20 30 40 50 60 Misalignment along Z axis Translation error (mm) X axis Y axis Z axis (a) (b) (c) Figure 4.10: Capture ranges of the proposed method: (a) translation errors for different misalign- ment initializations along X axis; (b) translation errors for different misalignment initializations along Y axis; (c) translation errors for different misalignment initializations along Z axis. −40 −30 −20 −10 0 10 20 30 40 0 10 20 30 40 50 60 70 Misalignment along X axis Translation error (mm) X axis Y axis Z axis −40 −30 −20 −10 0 10 20 30 40 0 10 20 30 40 50 60 70 80 Misalignment along Y axis Translation error (mm) X axis Y axis Z axis −20 −10 0 10 20 30 40 50 0 10 20 30 40 50 60 70 80 Misalignment along Z axis Translation error (mm) X axis Y axis Z axis (a) (b) (c) Figure 4.11: Capture ranges of the NMI based method: (a) translation errors for different mis- alignment initializations along X axis; (b) translation errors for different misalignment initializa- tions along Y axis; (c) translation errors for different misalignment initializations along Z axis. be observed that best matching phases were not different by more than 2 phases, as visually determined. In Figs. 4.10 and 4.11, we have presented the plots of capture range of each axis with different misalignment, which shows the proposed method is robust with different initializations compared to NMI based method. Average of 5 cases was computed according to the offset misalignment introduced along each axis. The proposed method provided robust result using different initial- izations of translation in the range of -80 to 70 mm, -80 to 70 mm and -50 to 50 mm from manual 69 ground truth in x, y, and z axes, respectively. As shown in Fig. 4.11, the capture range of the NMI based method was much smaller than the proposed method in the range of -20 to 20 mm, -20 to 20 mm and -10 to 40 mm in x, y, and z axes, respectively. 4.4 Discussion In this study, we proposed and evaluated a fully automated method to register cardiac CT (anatom- ical information) to MPS (physiological information). Anatomical features, segmentations of LV in the MPS volume, were automatically found in advance to guide the registration process. With a piecewise constant image model, both volumes were then registered using the region-based sum of squared difference as the similarity measure. Subsequently, the best matching phase of MPS was selected based on the minimum cost value, followed by the TPS warping using derived LV contours to match all image phases of gated MPS to the best matching cardiac phase position (i.e., the best corresponding to the CT phase). As compared with existing methods, we provided a cor- rect solution to the compensation of the moving heart using the TPS-based warping within gated MPS. This is the first report describing both the registration and the phase mismatch problems. The proposed methods were validated on MPS and cardiac CT (coronary CTA and CT perfusion) obtained from multiple MPS cameras and two CT scanners. The evaluation of our algorithm was performed with the rigid registration scheme. The subsequent nonlinear registration was applied only within gated MPS data for the visualization purpose. Thus, the registration performance should be measured by the registration scheme with rigid registration while the nonlinear regis- tration scheme is solely aimed for better visualization in the observation of defects. Robust and 70 accurate registration accuracy was obtained with respect to the ground truth, as visually deter- mined by two expert observers. The translational error was less than the size of a MPS voxel, and thus is feasible to be used for the clinical practice. For rotations in rigid registration, a small time step in the minimization step was adopted to constrain the angular movement to -2 to 2 de- grees, since both cardiac CT and MPS were acquired in a similar posture. Furthermore, the best matching phase was also validated using visual assessment from an expert observer. Registration of the patient data achieved 87% success rate with accuracy< 10 mm. As for the robustness test, the capture range was used. The capture range represents the range of positions from which a registration method can converge to the correct minima/maxima. The proposed method provided a more robust capture range as compared with that of the NMI based method. We additionally have presented validation of best matching phase method compared to visual assessment obtained from an expert observer. Best matching phases found by our proposed method were not different by more than 2 phases, as visually determined. These results are acceptable since visual assess- ment is not absolute and in some cases it was equivocal to find only 1 phase that corresponds to CT phase, as visually determined. In general, as compared with brain image registration, both cardiac image registration and its validation are inherently more difficult because of cardiac and respiratory motion [87]. Con- sequently, rigid registration alone may not be sufficient, and subsequent nonlinear component should be used to compensate for changes in organs or tissue due to the breathing pattern or movement of internal organs. Since the exact ground truth is unknown and only visual judgment is available, quantitative evaluation of the registration is challenging and subject to inter-observer variability. No significant differences were observed between inter-observer variability and the results of the proposed method. 71 Our proposed method might have potential limitations. First, we applied surface-based regis- tration using segmentation of left ventricle in MPS. However, this contour-driven approach will admittedly affect the registration accuracy when it comes to the poor segmentation result of MPS. In our system, if the MPS contours are not correct, they are manually adjusted for the perfusion and fusion analysis of MPS; therefore from the practical standpoint correct segmentation will be available in all cases, even if adjusted manually for purposes of MPS analysis. However, in our experiment, no manual segmentation adjustment aided by user was observed. Secondly, there are several possible image artifacts such as windmill artifact that might affect the performance of the proposed method. However, the patients in this study were without known coronary artery disease, and patients with existing coronary artery bypass grafts were excluded from the study. There were no datasets that had any metallic implants. For CTA data, we performed Gaussian smoothing to suppress noise and artifacts. We have previously published systematically catego- rized image artifacts in our patient population, and percentage of patients with artifacts due to metallic implants was 3% [33], suggesting that our proposed method will not be subject to this difficulty in the vast majority of applicable patients. In our experiment, there were 6 failed cases. The reason for failure was that all 6 datasets have large atria with different blood pool size be- tween CT and MPS rather than the artifacts present. These datasets can be considered as outliers of our assumption that each region should be corresponding in terms of size and constant intensity level. Although SPECT/CT scanners become more widely available, there are two reasons why it will be preferable to perform register MPS and CTA images acquired on standalone machines rather than perform a hybrid MPS/CTA scan. Firstly, current SPECT/CT scanners are configured with lower performance CT scanners. Typically it is a CT machine capable of attenuation cor- rection or calcium scoring but not of high-speed resolution CT angiography. SPECT/CT is not 72 offered in with 64-slice scanners which are preferred for CT angiography. Secondly, from the clinical point of view it is most valuable to know the results of one test before deciding if addi- tional test (and additional radiation dose) is required. For example, when the results are equivocal for CTA in case of heavily calcified arteries, a second test such as MPS can be ordered to resolve the diagnostic question - but this can not be predicted before the test. On a hybrid SPECT/CT ma- chine it would be difficult to modify the protocol on-the-fly based on the intermediate results of one test and therefore the use of standalone modalities is more practical clinically. The selective application of one of these tests or a combination as presented here will be most cost-effective as it will allow minimizing the ionizing radiation dose. In addition, even if MPS and CTA were per- formed in the same scanning session on the hybrid scan there will likely be residual registration errors during the respiratory differences, since CT angiography is performed during a breath-hold while MPS is performed during normal breathing and software image registration as presented here may still be required. The proposed technique can also be extended to PET-CT fusion since the segmentation technique has been applied and use in for PET imaging [114]. Consequently, after PET segmentation, the registration problem would be essentially the same as for MPS, since no raw MPS data are used in the registration. 4.5 Conclusions A multi-modal image registration method of cardiac CT with SPECT in a variational framework was presented. A segmentation method for extracting the blood pool, myocardium and extrac- ardiac structures from MPS was adopted prior to registration. In the registration process, these 73 segmented regions in MPS were aligned to the corresponding regions in contrast CT by utiliz- ing geometrically characteristic features obtained by MPS segmentation and the homogeneity property of intensities within regions of interest in contrast CT. A specialized cost functional for matching segmented heart image to cardiac CT was designed. We demonstrated by experimental results that the proposed method outperformed a standard registration method based on normal- ized mutual information. The results show that the proposed approach provides fast, accurate and robust registration with error less than the size of a MPS voxel as tested in datasets from two dif- ferent sites obtained with three different acquisition protocols. This fully automated method may facilitate integrated gated SPECT and cardiac CT analysis to resolve borderline cases in clinical practice, while reducing operator time and variability needed for manual adjustments. 74 Chapter 5 : Nonlinear Ultrasound Image Registration Based on Intensity and Local Phase Information 5.1 Introduction Ultrasound imaging has been widely used as an medical imaging modality during last several decades due to its non-invasiveness, safety and economical efficiency [119, 120]. For instance, the use of ultrasound is suitable for fetus imaging and the diagnosis and assessment of arterial and cardiac diseases [97]. Left ventricular contractile function can also be evaluated effectively by echocardiography imaging [101]. However, ultrasound images have a low signal-to-noise (SNR) ratio due to artifacts and attenuation. Furthermore, speckle noise caused by sub-resolution tissue structures is a big foe in analyzing lesions and performing various processes such as segmentation and registration by limiting the detectability of low contrast lesions. Due to these limitations, validation of estimated of motion is currently performed in a subjective visual manner. Registration of ultrasound images can be used in estimating cardiac motion and/or velocities to find abnormalities in segmental wall motion and analyzing tissue mechanic properties such as 75 quantification of elasticity and contractility of the myocardium [78]. It could be also used for spa- tial compounding of images to improve image quality and contrast [73]. In general, registration of ultrasound images is a challenging task as compared with that of CT (Computed Tomography) or MRI due to special properties of ultrasound images, including speckle noise artifacts and low images contrast. Broadly speaking, registration algorithms can be categorized to two classes: rigid registration using linear transformation (translation, rotation and scaling) and nonlinear registration using nonlinear transformation to compensate for changes in organs or tissue due to the breathing pat- tern or movement of internal organs [121]. Nonlinear registration using variational minimization and/or various similarity measures is well studied in the context of computer vision. For example, to determine optical flow in temporal images can be viewed as a registration problem [16, 17]. In this chapter, we propose a novel registration method that exploits both intensity and local phase information by incorporating them in a variational framework and present some preliminary results when this technique is applied to test ultrasound images. The local phase feature is well- suited for ultrasound image registration since it is invariant to image brightness, contrast and noise [90]. We combined this feature with image intensity to register different frames in more accurate manner. To the best of our knowledge, this is the first effort in deriving an energy functional using the intensity and the local phase information under the variational framework to solve the ultrasound image registration problem. The rest of this paper is organized as follows. Some previous work is reviewed in Sec. 5.2. The local phase and its properties are discussed in Sec. 5.3. The registration method derived based on the variational framework is presented in Sec. 5.4. Experimental results are shown in Sec. 5.5. Finally, concluding remarks and future research directions are given in Sec. 5.6. 76 5.2 Review of Previous Work Several registration methods have been proposed before. They are usually classified into two types: the feature-based approach using extracted image features and the volume-based method using statistical voxel dependencies [87]. Although intensity values contain all available informa- tion in images, other derived features may reveal the underlying image structure more clearly [90]. Thus, we may need to search for additional relevant features and then study how to apply these features jointly for efficient nonlinear ultrasound image registration [121]. Different features have been incorporated to tackle the ultrasound image registration problem. Quite a few registration methods are based on the intensity information in the image registration [20, 54, 64] under the assumption that mono-modal images are identical when correctly aligned. An intensity-based deformable registration technique for 3D ultrasound images was presented in [144] under a vari- ational framework. However, this method may not work well since the assumption of gray value constancy may not hold in ultrasound images. For example, it is difficult to detect the heart muscle consistently, thus resulting in missing and/or false positives. Gradient-based methods are conducted based on the assumption that the contrast of a pattern is invariant over time. However, computed gradient values are highly sensitive to noise since the differentiation of gray levels tends to magnify noise, leading to inaccurate results. Intensity and its derivatives were used to generate an attribute vector in [47]. However, in general, the gradient-based approach may not be suitable for the processing of noisy ultrasound images. Another class of registration methods is developed based on finding the corresponding landmarks and, furthermore, two images are aligned with approximating thin-plate splines [110]. However, it demands manual selection of landmark points, and finding landmark points is a difficult task especially for ultrasound image 77 because of its limited field-of-view. The local phase feature has been developed by Morrone et al. [93] and has become increasingly popular. Its robustness and accuracy for ultrasound images were studied in [55, 91, 94]. The local phase information can provide more local structural in- formation than the intensity feature. We demonstrate in this chapter that, when the intensity and the local phase information can be combined, they can supplement each other, thereby achieving better registration results. 5.3 Local Phase Computation and Properties 5.3.1 Local Phase and Local Energy The local phase function provides a qualitative and contrast invariant local description of an image. It can be interpreted as an amplitude-weighted phase of a windowed (bandpass filtered) Fourier component of a signal [91]. We will show in this section that the local phase serves as an effective feature for ultrasound image registration. The concept of local phase and local energy were originally developed for analyzing one dimensional signals. The analytic signal is most easily interpreted in terms of Fourier analysis [90]. (a) (b) (c) Figure 5.1: (a) An examplary cardiac ultrasound image, and its local phases using two DoG bandpass filters with (b) a small standard deviation (s 1 = 2 p 2;s 2 = 2) and (c) a large standard deviation (s 1 = 4 p 2;s 2 = 4). 78 The analytic signal, f A (x), of a real signal, f(x), in the space domain can be written as f A (x)= f(x) j f H (x); (5.1) where f H (x) is the Hilbert transform of f(x), defined by f H (x)= 1 p Z ¥ ¥ f(t)dt t x , F H (w)= F(w) j sign(w); (5.2) where F(w) is the Fourier transform of f(x) and sign(w)= 8 > > < > > : 1; w < 0; +1; w> 0: (5.3) As given in Eq. (5.1), F A (w) is equal to zero at the negative frequency and is twice of the original signal at the positive frequency, which results in an analytic signal given by [13, 41] F A (w)= F(w) j F H (w) F A (w)= F(w) j F(w) j sign(w) F A (w)= F(w)[1+ sign(w)] (5.4) Then, the local phasef(x) and local energy A(x) of signal f(x) is defined as f(x)= tan 1 ( f(x)= f H (x)); A(x)=k f A (x)k= q f 2 (x)+ f 2 H (x): (5.5) 79 However, it is difficult to compute the local phase based on the above definition directly. To extract the proper feature, both time and frequency localization is desired. In other words, we may apply a window function of bandpass filtering characteristics to f(x) with an objective to reduce the fine scale noise effect and localize the desired signal. The window function has two properties. First, it should be even and symmetric to maintain the phase information of the original signal. Second, it should have the zero DC value so that steady-state values will not be affected. There are several choices in the design of bandpass filters. We adopt the Difference of Gaussian (DoG) in our chapter. 5.3.2 Monogenic Signal It is difficult to extend the Hilbert transform from 1D to 2D directly. For several years it is found that there is no straightforward extension of the Hilbert transform to higher dimensions. Instead, the monogenic signal was proposed in [41] to estimate and extend the phase notion to 2D effi- ciently. This representation preserves core properties such as isotropic properties of 1D analytic signals and decomposes a signal into the energy and the structure (local phase) features. The even filter is rotationally symmetric (DoG), and the single odd filter used previously is replaced by two odd isotropic filters (vector-valued) in the frequency domain as H 1 = u p u 2 + v 2 ; H 2 = v p u 2 + v 2 ; (5.6) 80 where u and v are frequency-domain variables. Then, the local phase,f, and local energy can be calculated using a bandpass filtered image I b and its filter responses as f(x;y)= tan 1 ( I b p (h 1 I b )+(h 2 I b ) ); E(x;y)= q I 2 +(h 1 I b ) 2 +(h 2 I b ) 2 ; (5.7) where denotes the convolution operator. The local phase can be interpreted as a qualitative description of salient regions in images such as edges or ridges. An examplary cardiac ultrasound image, and its local phases using two DoG bandpass filters with a small and a large standard deviations are shown in Fig. 5.1. We see from this figure that the local phase provides valuable structure information of the underlying image. 5.4 Registration Method The problem of ultrasound image registration is equivalent to finding the displacement vector un- der diffeomorphism between the source and the target images. One application of this technique is to register two consecutive images in an ultrasound image sequence. Then, the source and the target images are the previous and the current image frames, respectively, and the displacement vector denotes the motion field of an image pixel. 5.4.1 Variational Framework Before deriving a variational formulation for the registration problem, we provide some intu- itive explanations on constraints and assumptions used in this model. The proposed registration method is derived from a variational framework based on the following observations. 81 The intensity and its derivatives, which are commonly used as features in mono-modal image registration and optical flow estimation, have similar values over two consecutive image frames. The property, called gray value constancy, can be written mathematically as I(x;y;t) I(x+ u;y+ v;t+ 1) ÑI(x;y;t)ÑI(x+ u;y+ v;t+ 1); (5.8) where I :WR 3 !R denotes an image sequence,(u;v) is the local displacement vector of a pixel over two adjacent frames, t is the frame index number, andÑ=(¶ x ;¶ y ) is the 2D spatial gradient vector. The linearized formula of Eq. (5.8) yields the following optical flow constraint [68]: I x u+ I y v+ I t = 0; (5.9) where I x =¶ x I and I y =¶ y I. However, Eq. (5.9) alone does not carry sufficient and accurate information for ultrasound image registration due to low contrast and low resolution nature of the underlying images. In contrast, as explained in Sec. 5.3, the local phase value is a more accurate feature than the intensity value in ultrasound imaging. It is assumed that local phase values are also similar along their temporal trajectory curves: LP(x;y;t)= LP(x+ u;y+ v;t+ 1); (5.10) where LP is the local phase of a given image. The registration problem is however an ill-posed one. That is, the model given in Eqs. (5.8) and (5.10) estimates the displacement of a pixel only 82 without taking any neighborhood relation into account. Thus, a piecewise smooth flow is needed for regularization since a problem occurs in the model when the gradient vanishes. Let x=(x;y) T and w(x)= w(x;y)=(u(x;y);v(x;y)) T . Then, we can incorporate the intensity and the local phase in the data term as E data (u;v)= Z W Y(jS I (x+ w(x)) T I (x)j 2 +gjS LP (x+ w(x)) T LP (x)j 2 )dx; (5.11) where g is a weight, and Y is a quadrature penalizer (e.g., Y(s 2 )= p s 2 +e 2 , which results in modified L 1 norm), and S I , S LP , T I and T LP are the intensity and the local phase of source image S and target image T , respectively. The smoothness term is derived under the assumption of a piecewise smooth flow field. This can be achieved by penalizing the total variation of the flow field and expressed as E smooth (u;v)= Z W Y(jÑuj 2 +jÑvj 2 )dx (5.12) with the sameY function as defined above. The total energy functional can be expressed as the weighted sum of the data term and the smoothness term: E(u;v)= E data +aE smooth : (5.13) By adjustinga, we can control the balance between the data term and the smoothness term. The smoothness term can prevent overfitting to the data term and thus produce a piecewise continuous solution. 83 5.4.2 Euler-Langrange Optimization The energy functional in Eq. (5.13) can be minimized for the desired gradient descent using the Euler-Lagrange method. With the calculus of variation, the proposed energy functional satisfies the following Euler-Lagrange equation with respect to u and v as ˜ Y(I 2 d +gLP 2 d )(I d I x + LP d LP x )a div( ˜ Y(jÑuj 2 +jÑvj 2 )Ñu)= 0; ˜ Y(I 2 d +gLP 2 d )(I d I y + LP d LP y )a div( ˜ Y(jÑuj 2 +jÑvj 2 )Ñv)= 0; (5.14) where I d := S I (x+ w) T I (x);LP d := S LP (x+ w) T LP (x);I x :=¶ x S I (x+ w); I y :=¶ y S I (x+ w);LP x :=¶ x S LP (x+ w);LP y :=¶ y S LP (x+ w): (5.15) and ˜ Y(x)=Y 0 (x)= 1 2 p x+e 2 and wheree is a small number. In our experiment, we choosee = 0:001. 5.4.3 Numerical Approximation To solve the Euler-lagrange equations in (5.14), we adopt an alternating minimization method. That is, we first minimize the energy functional with respect to variable u with fixed variable v and, then, minimize the energy functional with respect to variable v with fixed variable u. This alternating optimization process can be performed iteratively. 84 Let us consider the minimization process with respect to variable u. We define F N (e)= E(u+eN;v)= Z W Y(jS I (x+ u+eN;y+ v) T I (x;y)j 2 +gjS LP (x+ u+eN;y+ v) T LP (x;y)j 2 )dx +a Z W Y(jÑ(u+eN)j 2 +jÑvj 2 )dx (5.16) wheree is a small parameter. If u is the minimizer, we have F N (0)6 F N (e);8e. Thus, we have derived the following condition: dF N (e) de u=u ;e=0 = 0: (5.17) Let F 0 N (e)= dF N (e) de : It is straightforward to find that F 0 N (e)= Z W 2Y 0 (jS I (x+ u;y+ v) T I (x;y)j 2 +gjS LP (x+ u;y+ v)T LP (x;y)j 2 ) ¶S I (x+ u+eN;y+ v) ¶(x+ u+eN) +g ¶S LP (x+ u+eN;y+ v) ¶(x+ u+eN) Ndx +a Z W 2Y 0 (jÑuj 2 +jÑvj 2 ) ¶u ¶x +e ¶N ¶x ¶N ¶x + ¶u ¶y +e ¶N ¶y ¶N ¶y dx (5.18) 85 By settinge = 0, we obtain F 0 N (0)= Z W 2Y 0 (jS I (x+ u;y+ v) T I (x;y)j 2 +gjS LP (x+ u;y+ v) T LP (x;y)j 2 ) ¶S I (x+ u;y+ v) ¶(x+ u) +g ¶S LP (x+ u;y+ v) ¶(x+ u) Ndx+a Z W 2Y 0 (jÑuj 2 +jÑvj 2 ) ¶u ¶x ¶N ¶x + ¶u ¶y ¶N ¶y dx (5.19) Then, we can apply integration by parts to the last integral as Z W (¶ i u) vdx= Z W u(¶ i v)dx+ Z ¶W (u n i ) vds; (5.20) where n i is a component along the normal direction. By setting F 0 N (0)= 0;8N, we obtain the Euler-Lagrange equation in u as Y 0 (jS I (x+ u;y+ v) T I (x;y)j 2 +gjS LP (x+ u;y+ v) T LP (x;y)j 2 ) ¶S I (x+ u;y+ v) ¶(x+ u) +g ¶S LP (x+ u;y+ v) ¶(x+ u) =aY 0 (jÑuj 2 +jÑvj 2 )(Du) (5.21) The above equation can be put in time-dependent form using the gradient descent as ¶u ¶t = dE du = dF N (e) de e=0 ; (5.22) 86 where t is an artificial time. Then, we have ¶u ¶t =Y 0 (jS I (x+ u;y+ v) T I (x;y)j 2 +gjS LP (x+ u;y+ v) T LP (x;y)j 2 ) ¶S I (x+ u;y+ v) ¶(x+ u) +g ¶S LP (x+ u;y+ v) ¶(x+ u) +aY 0 (jÑuj 2 +jÑvj 2 )(Du) (5.23) We are able to solve Eq. (5.23) for u numerically. In the numerical implementation, we adopt the following approximation for time differencing: ¶u ¶t u (n+1) (i; j) u (n) (i; j) Dt ; (5.24) and the following 5-point approximation for the Laplacian operator: Du= u(i+ 1; j) 2u(i; j)+ u(i 1; j) h 2 + u(i; j+ 1) 2u(i; j)+ u(i; j 1) h 2 : (5.25) As a result, we adopt the explicit scheme to solve the above equation with smallDt of the follow- ing form: u (n+1) (i; j) u (n) (i; j) Dt =Y 0 (jS I (x+ u;y+ v) T I (x;y)j 2 +gjS LP (x+ u;y+ v) T LP (x;y)j 2 ) ¶S I (x+ u;y+ v) ¶(x+ u) +g ¶S LP (x+ u;y+ v) ¶(x+ u) +aY 0 (jÑuj 2 +jÑvj 2 )(Du) (5.26) Following a similar procedure, the energy functional can be minimized with respect to v. 87 5.5 Experimental Results 5.5.1 Synthetic Images We begin with some simple yet illustrative examples to demonstrate the efficiency and robustness of the proposed algorithm. Two synthetic images are used in our quantitative validation as shown in Fig. 5.2 (a) and 5.3 (a). Fig. 5.2 (a) consists of simple geometrical shapes while Fig. 5.3 (a) is the cyst phantom that consists of a collection of point targets, five cyst regions and five highly scattering regions with SNR 40 dB, followed by adding artificial speckle noise. Then, these two synthetic images are deformed by randomly generated displacement fields as shown in Fig. 5.2 (b) and Fig. 5.3 (b), respectively . Our goal is to evaluate registration accuracy under nonlinear deformation. The displacement vectors that incur deformations are generated by random fields with pixel difference of 15 and 30 followed by Gaussian smoothing with s = 6 and s = 15 in Fig. 5.2 (b) and Fig. 5.3 (b), respectively. (a) Example 1 (b) Deformed image Figure 5.2: (a) The first synthetic image and (b) its deformed image. Experimental results with these two synthetic examples are shown in Table 5.1, where w and w 0 denote the ground truth displacement vector and the computed one using different methods, 88 (a) Example 2 (b) Deformed image Figure 5.3: (a) The second synthetic image and (b) its deformed image. Table 5.1: Comparison of three different data terms (unit:pixel). Method I I+ G Proposed Example1 m ww 0 0.648 0.657 0.602 s ww 0 0.995 0.918 0.850 SSD 201.028 203.943 190.403 Example2 m ww 0 4.150 4.015 3.776 s ww 0 6.747 6.606 6.773 SSD 23.183 22.731 21.719 respectively. The mean (m ww 0) and the standard deviation (s ww 0) of the error and the mean of the Sum of Squared Differences (SSD) between target and registered images are used as quantitative performance metrics. Three methods using different features are compared in Table 5.1. First, we only use the intensity term [144]; namely, Z W Y(jS I (x+ w(x)) T I (x)j 2 dx; 89 in the data term, which denoted by I in Table 5.1. Next, both the intensity and its gradient features [16] are used in the data term. Mathematically, we have Z W Y(jS I (x+ w(x)) T I (x)j 2 +g(jÑS I (x+ w(x))ÑT I (x)j 2 )dx; which is denoted by I+ G in the table. The third one is the proposed method that uses both the intensity and the local phase information. As to the smoothness term, the first two methods adopt the same one as the proposed method. We see from Table 5.1 the proposed method has a smaller registration error in terms of the mean and the standard deviation as well as smaller mean SSD than the other two methods. We can adjust the importance of the data term and the smoothness term by tuning parameter a. Likewise, by choosingg, we can adjust the importance between different data term including an intensity, gradient and local phase term. We find the parameter values heuristically which minimize the energy functional most. 5.5.2 Invivo Image Furthermore, we use the in vivo cardiac ultrasound images of mice and human subjects for per- formance comparison. Gaussian smoothing is first performed on these ultrasound images as a pre-processing step since the intensity and its gradient can be sensitive to noise, resulting in un- desired matching. The smoothed source and the target ultrasound images of mice and human subject are shown in Figs. 5.4 and Figs. 5.5 (a) and (b) along with their intensity difference in Fig. 5.4 and Fig. 5.5 (c). For visual evaluation of the proposed algorithm, the difference image between the target and the registered images is presented in Figs. 5.4 and Fig. 5.5 (d) and the 90 displacement vector from the source image and the deformation field are presented in (e) and (f), respectively. (a) (b) (c) (d) (e) (f) Figure 5.4: In vivo mice ultrasound images: (a) the source image, (b) the target image, (c) the difference image between source and target images, (d) results obtained by the proposed method, (e) the displacement vector from the boundary in (a) (arrows being scaled 1.5 times for the visu- alization purpose), and (f) the deformation field. As shown in the obtained deformation fields, the proposed method is able to find the displace- ment of pixels in the dark region that has similar intensity levels. Since there is no ground truth in the in vivo data, we compare registration accuracy using various similarity measures such as the mean SSD, the Mutual Information (MI) and the Normalized Correlation (NC) in Tables 5.2 and 5.3. Since we used the L1 norm which is similar to the mean SSD, comparison using more generic metrics such as MI and NC demonstrates the validity of the mean SSD as a similarity measure in ultrasound image registration. 91 (a) (b) (c) (d) (e) (f) Figure 5.5: In vivo human subject ultrasound images: (a) the source image, (b) the target image, (c) the difference image between source and target images, (d) results obtained by the proposed method, (e) the displacement vector from the myocardium boundary in (a) (arrows being scaled 1.5 times for the visualization purpose), and (f) the deformation field. Table 5.2: Comparison of three different similarity measures for the mice image registration Method I I+ G Proposed Mice Data SSD 26.145 25.99 25.885 MI 1.182 1.182 1.183 NC 0.995 0.995 0.996 We see that there is a general trend between mean SSD, MI and NC in Tables 5.2 and 5.3. That is, when SSD is lower, both MI and NC become higher. This demonstrates the usefulness of the similarity measure (i.e., the mean SSD) as well as the energy functional. 5.6 Conclusion We presented a novel nonlinear registration method for ultrasound images in this chapter. Specif- ically, we proposed to use the local phase function as a geometric feature to find correspondences 92 Table 5.3: Comparison of three different similarity measures for the human subject image regis- tration Method I I+ G Proposed Human Data SSD 37.377 37.126 34.768 MI 1.203 1.21 1.221 NC 0.957 0.959 0.963 of pixels in two ultrasound images that suffer from speckle, artifact and occlusion. By combining both the local phase and the intensity information under a variational framework, the proposed method outperforms others using the intensity and its gradient. Different similarity measures, including mean SSD, NC and MI, were used to compare the performance of different algorithms. Clinically, the proposed method can be applied to the analysis of tissue mechanic properties and object motion (e.g. beating heart) for the treatment of the myocardial perfusion and other diseases. In the future, we would like to extend this method so as to find motion of tissues in multiple frames, which is similar to optical flow and other vision related problems such as mo- tion segmentation and tracking. Also the local phase can be used for velocity estimation or other tracking purposes as well. 93 Chapter 6 : Nonlinear Registration of Serial Coronary CT Angiography (CCTA) for Plaque Development Monitoring 6.1 Introduction Coronary artery disease (CAD) is a major cause of mortality in the industrial world. Coronary CT angiography (CCTA) has become an increasingly effective clinical tool for non-invasive assess- ment of the coronary arteries following the introduction of 64-slice CT scanners with fast gantry rotation times [2, 9]. CCTA has demonstrated a great potential for plaque detection, quantifica- tion, and characterization [76, 77], and assessment of positive coronary artery remodeling [3], which is associated with an unstable plaque. Recent studies have suggested that serial CCTA may be used to monitor untreated plaque progression over time [115], or assess plaque changes in response associated with lipid-lowering therapy [18]. Currently, serial CCTA scans of the same patient, acquired at two separate times, are as- sessed separately and manually with multi-planar orientation and/or curvilinear displays [42]. This process is time-consuming and it is difficult to ensure that the identical views of the artery are presented to the clinician in both scans. CT plaque imaging is challenging in general because 94 of cardiac motion and the small size of coronary arteries (2-3 mm in diameter) [1]. Detection of plaque progression is particularly challenging because areas exhibiting changes within the coronary arteries but displayed over several slices. Therefore for visual comparison it is not guar- anteed that 2D images (either multi planar reconstructions or curvilinear displays) refer to the same image position. Due to the tortuous anatomy of the coronary vessels, even small differ- ences in image orientation and curvilinear image surface position may result in apparent changes in perceived plaque characteristics, introducing significant subjectivity. Recently studies [22,126] on manual plaque quantification of serial scans shows that variability between different observers in measuring plaque volume and composition quantification can be quite significant, especially in small plaques. Accurate and automated co-registration of serial CCTA scans would ensure that multipla- nar/curvilinear views created for lesion comparison are visualized in a common frame of refer- ence. It would allow voxel-wise comparison of coronary plaque composition, arterial stenosis and remodeling at two time-points, and potentially allow direct quantitative monitoring of plaque change. However, co-registration for assessment of serial CCTA poses several challenges. First, there could be significant intensity differences between the two datasets due to acquisition at two separate times, with two different contrast injections and parameters occurring at potentially dif- ferent physiological states. Second, the algorithm must deform the source image in the coronary artery region but preserve true plaque characteristics of source volume. Finally, there could be missing portions or occlusion due to plaque development between the two datasets, which might cause intensity differences and makes establishing correspondence challenging. With regard to serial CCTA co-registration, rigid registration alone is not sufficient to com- pensate for inherent cardiac and respiratory motion [87]. Consequently, nonlinear component 95 should be incorporated to describe the detailed deformation. Nonlinear registration generally refers to the process of estimating the nonlinear mapping function, which transforms each point in one image to a point in the target image. It can be categorized according to the following three aspects, (i) similarity measure, (ii) nature of the transformation, and (iii) optimization pro- cedure. In principle, a similarity measure is required to be convex, a transformation to be smooth and invertible, and an optimization method to guarantee the globally optimal solution. Common similarity measures are landmark matching [109], sum of squared differences (SSD) [65], cross correlation (CC) [27], and mutual information (MI) [124]. SSD and CC are well-suited when the source and target images have been acquired similarly, and are expected to present similarly in both intensity range and distribution [66]. The MI measure has been adopted successfully for multi-modal image registration without specific assumptions about the two images. With respect to transformation, most registration techniques can be divided into two types. In the first type, a nonparametric model is adopted, and the mapping of each voxel is estimated using an inter- polation method [19, 83]. In the second type, the deformation is determined by using a function of select parameters, which are estimated by using polynomial transformation or spline based methods [74, 116]. Other approaches to regularizing deformation have also been proposed to compensate for the ill-posedness of correspondence problems. They include biharmonic [44,45], diffusion [43], and physical continuum models [31, 79], such as linear elastic and viscous fluid regularizations. In addition to regularization techniques, there are also methods that impose con- straints on transformation such as inverse consistency [24], topology preservation [66], and dif- feomorphism [25]. Although this problem has not been addressed before, some related techniques have been proposed previously. Registration was proposed for application of measuring volume or shape 96 changes in MRI scans [30, 61, 80, 106]; in addition, related serial nonlinear CT registration tech- niques have been proposed for estimating interval changes for the application of dose tracking, CT radiotherapy planning, tracking of CT, and temporal subtraction of thoracic CT [85,118,122, 130, 134]. Nonetheless, CCTA application is uniquely demanding due to the high resolution im- ages, complex three-dimensional anatomy, and the requirement for high accuracy to preserve the plaque characteristics after registration. In this study, we propose a novel and fully automatic serial co-registration method for CCTA scans using a dense nonparametric co-registration model via diffeomorphism by extending our preliminary approach which combines segmentation of the vessels with image co-registration [135]. Aiming to solve the difficulties mentioned above, our method utilizes a binary feature mask, ob- tained by segmentation of the arterial lumen, with global displacement. The global displacement refers to rigid registration as an initialization. To avoid possible bias caused by incorrect segmentation, we introduce a soft mask to con- sider the coronary artery and surrounding region; a volume-preserving constraint is also used to guarantee that the volume remains constant before and after the co-registration. Co-registration is a two-step procedure which consists of global displacement using binary feature masks and subsequent local deformation model under a variational framework. Prior to local deformation, histogram matching is performed to compensate for potential intensity differences, resulting in similar image characteristics (i.e. intensity) between the two datasets. The proposed approach in- corporates combined information of intensity and geometrical features thereby establishing more accurate and robust correspondences. The rest of this chapter is organized as follows. A hierarchical transformation model that captures both global displacement and local deformation is presented in Sec. 6.2. Experimental 97 validation is described in Sec. 6.3. Experimental results are provided to support theoretical as- sumptions, and the performance of the proposed registration algorithm is demonstrated using real patient datasets in Sec. 6.4. Discussion and conclusions are given in Secs. 6.5 and 6.6, respec- tively. 6.2 Methods 6.2.1 Description of Method Global Displacement Baseline CCTA Volume Baseline Feature Mask Follow-up CCTA Volume Follow-up Feature Mask Histogram Matching Local Deformation Figure 6.1: The flowchart of the proposed method. We classify the problem as the same-modality and intra-subject co-registration, since serial CCTA scans are obtained from the same patient at different time instances. The goal of our serial 98 co-registration is to estimate transformation that aligns coronary arteries in 3D, thereby allow- ing assessment of the regression and progression of coronary plaques in the arteries. To sum- marize, our co-registration scheme consists of two sequential steps: (1) a rigid co-registration scheme for global displacement based on a pre-segmented binary mask, and (2) a final nonlin- ear co-registration that accounts for local vessel deformation after histogram matching of inten- sity distribution in the two datasets. Global displacement is initially obtained by the optimal rigid co-registration between the two pre-segmented binary masks, and the local deformation is subsequently obtained by optimal nonlinear co-registration via diffeomorphism, with a volume- preserving constraint. In the co-registration process, we utilize coronary trees that are initially obtained by the existing coronary artery segmentation software and use them as initial feature masks. The flowchart of the proposed co-registration method is given in Fig. 6.1. 6.2.2 Feature Mask (a) (b) (c) (d) Figure 6.2: Illustration of the coronary tree mask used in image registration: (a) a representative 2D slice (yellow plane) of the pre-segmented 3D binary mask that covers the coronary tree region; (b) the representative slice from a 3D image overlaid with the pre-segmented mask boundary (black); (c) a 2D pre-segmented binary mask for this slice, and (d) the soft mask generated by convolving with the Gaussian kernel. We utilize binary coronary trees as initial feature masks. The binary 3D map containing vox- els of the coronary arterial lumen may not be accurate and complete due to initial segmentation 99 errors (missing portions, inclusion of noise and extra coronary structures), which could lead to subsequent bias in the registration process. However, it does provide a useful reference for the initial region of interest as long as plaques of interest are contained within. We utilize a 3D soft mask, generated by convolving the Gaussian kernel with the binary coronary tree mask, to account for neighboring vessel region; this can be viewed as a Gaussian weighting function to emphasize the central area and weaken the bias of pre-segmented binary maps. Examples of the binary and soft mask are shown in Fig. 6.2. Binary masks of coronary trees in target volume I and source volume J are M I ;M J :W!R + , respectively. They were obtained automatically with the CCTA segmentation software (Circula- tion, Syngo MMWP Version VE31A, Siemens Medical Solutions) [63] available on the CCTA workstation. The smoothed binary masks, referred here as soft mask, favor certain locations of domain x2W for source and target volume and are given by m I = M I G s ;m J = M J G s ; (6.1) where * denotes a convolutional operator, and G s is a Gaussian kernel with standard deviations. 6.2.3 Global Displacement Model The global displacement model is described by a rigid registration which includes rotation and translation (6 degrees of freedom in 3D) to account for large displacements. Since the 3D voxel masks with pre-segmented coronary tree regions are routinely and automatically obtained from CCTA scans, we can use them for source and target shapes in the initial registration without reducing practicality of the approach. We employ explicit parameters to account for translation 100 (by displacement parameter m) and rotation (by angular parameter q) from source image M J to target image M I , which are defined in the image domain, W. The energy functional for rigid transformation is given by E global (x;m;q)= Z W (M I (x) M J (R q x+m)) 2 dx (6.2) Furthermore, we employ a coarse-to-fine multi-resolution scheme for computational efficiency and robustness. The optimal values for the 6 parameters that minimize the given energy functional are obtained by solving the associated Euler-Lagrange equations, and a gradient descent method is performed. A global displacement model is used as the initialization for feature-based nonlinear volume registration in the later stage. 6.2.4 Local Deformation Model Rigid transformation only captures the global displacement of the coronary anatomy between two scans. Further local deformation should be considered to achieve more accurate co-registration results due to myocardial tissue deformations and residual cardiac phase differences. The basic assumption of the proposed energy functional is similar to that of optical flow [68] that maintains the brightness constancy assumption in the two datasets. To compensate for poten- tial intensity differences between them, contrast correction is performed using histogram match- ing. Additionally, to preserve the true plaque morphology of plaque lesions in the process of serial co-registration, several constraints are imposed in the nonlinear co-registration process. First, the primary structures of interest in both co-registration and comparison are the plaque 101 regions within the coronary artery trees in the two volumes. Second, the total volume of the coro- nary vessels should be preserved in the transformation. We therefore exploit the coronary tree binary masks and obtain local displacement by an optimal nonlinear transformation, represented by diffeomorphism with a volume-preserving constraint as detailed below. Finally, we incorpo- rate regularization term to establish correspondences with spatial coherence even in the presence of potential plaque development. Let I and J be the source volume and the target volume contained inW, which corresponds to an open and bounded domain in R n . Without loss of generality, we consider the 3D case. Registration between two volumes is equivalent to finding a smooth diffeomorphism defined on the image domainW: h(x) :W!W; (6.3) where h(x) denotes a displacement vector at each voxel that allows infinite dimensional diffeo- morphism. The optimal solution is found by minimizing an energy functional which consists of a data fidelity term that measures the dissimilarity of source and target volume within the soft mask region, and a smoothness term for volume matching as follows: E(x)= E data (x)+aE smooth (x); (6.4) wherea is a weighting parameter. For the data fidelity term, we exploit the normalized intensity difference within the soft mask region as a similarity measure given by E data = Z W m I (x)jI(x) J h(x)j 2 dx: (6.5) 102 We can emphasize the coronary tree region using the soft mask as a weighting function that restricts the domain of the similarity measure. Recovering the optimal potential of this objective function is not straightforward. Due to the lack of information, the complete recovery of the local deformation field is an ill-posed problem. The use of regularization is a common practice to overcome this limitation [43] with smoothness of the pixel-wise motion field as a natural registration assumption. The regularization term penalizes large variations in the transformation to ensure that dis- placement h(x) is continuous. We adopt the so-called diffusion regularization [43] that penalizes the total variation of the flow field, which can be expressed as E smooth = Z W jÑh(x)j 2 dx: (6.6) Now, the objective is to find h(x) that minimizes this energy functional while preserving the volume. We use the soft mask as a weighting function, which provides a mechanism to include the local neighborhood of regions of interest in the registration process. In our registration scheme, we applied Gaussian kernel with standard deviation s (e.g. 10 pixels) to include neighboring region and to weaken the bias caused by the segmentation algorithm. 6.2.5 Constrained Volume Preserving Flow Nonlinear co-registration can potentially change the original volume of structures, which in our case is highly undesirable. As stated before, the volume-preservation constraint is necessary to monitor changes of the target plaque region. The idea of using the Jacobian determinant 103 in preserving volume when registering two volumes has been widely exploited in the past [37, 57, 108] and this volume preserving regularization term could reduce the volume change of the total mass between source and target images. This regularization restricts the deformation to be incompressible. In this study, however, we will exploit different approach to preserve the volume of the coro- nary artery tree explicitly so that the volume of soft mask J remains constant while the flow is minimized using the projection technique [138]. First, we compute the variation of E(h(x))= Z W m I (x)jI(x) J(h(x))j 2 dx+a Z W jÑh(x)j 2 dx: (6.7) The variation of E (if E is perturbed infinitesimally in the direction ˜ h from h) is E(h) ˜ h= Z W [2m I (x)(I(x) J(h(x)))ÑJ(h(x))aDh(x)] ˜ hdx Z ¶W (Ñh 1 (x)+Ñh 2 (x)+Ñh 3 (x)) ˜ h(x)ds(x): (6.8) This leads to ÑE(h(x))=2m I (x)(I(x) J h(x))ÑJ h(x)aDh(x): (6.9) We then construct a flow that minimizes E and also preserves the volume of m J h(x): V(h)= Z W m J (h(x))dx: (6.10) The gradient of V is given by ÑV(h)=Ñm J h(x): (6.11) 104 The flow that preserves volume and minimizes E becomes ÑE+lÑV; (6.12) where we choosel so that dV(h)(ÑE+lÑV)= 0; (6.13) because the volume does not change. Therefore, we get Z W (ÑE+lÑV)ÑV dx= 0: (6.14) By solving (6.14) forl, we have l = R W ÑEÑV dx R W ÑVÑV dx = R W [2m I (I J h)ÑJ h+aDh]ÑJ hdx R W jÑJ(h(x))j 2 dx : (6.15) This derivation is based on the projection technique that imposes a condition that the total flow in Eq. (6.12) should be orthogonal to the volume flow in Eq. (6.11) for preserving the volume. V olume constraint parameter,l, is updated according to Eq. (6.15) while the energy is minimized. 6.2.6 Minimization Given the initialization of global displacement, we can find the optimal transformation by mini- mizing the flow associated with data fidelity, smoothness and the volume-preserving constraint. 105 The Euler-Lagrange equation, associated with a flow that minimizes E and preserves volume V , is given by h t+1 = h t +dt(ÑE(h t )+lÑV(h t )) = h t +dt(2m I (I J h t )ÑJ h t +aDh t +lÑm J h t ); (6.16) where l is given in Eq. (6.15), dt is the time step, t is the time index, and a is a weighting parameter that defines the trade off between the alignment of two volumes and the smoothness of the transformation. We found experimentally thata = 0.1 provides a good compromise between the two terms. For computational efficiency, global displacement is computed using an iterative multi-resolution technique. In the subsequent stage, the energy functional of local deformation is minimized nu- merically using an iterative gradient descent technique. The algorithm stops if a local minimum of the energy functional is found. 6.2.7 Implementation The proposed co-registration algorithm is fully automatic and the global displacement model and subsequent local deformation model as described above are implemented using rigid co- registration in ITK (Insight Segmentation and registration Toolkit) [69]. In addition, we used pub- licly available implementation [32] of dense nonparametric MI based nonlinear co-registration and ITK [69] implementation of demons algorithm [131] for comparison with the proposed method. We employed a multi-resolution scheme for our global displacement model. That is, a multi- scale pyramid is built for both the source and the target images, and the energy functional is minimized using a coarse-to-fine approach. The registration parameters are first estimated at the 106 coarsest level, and these are then used as the starting parameters for the next level. Transforma- tions at each level of the pyramid are accumulated and propagated, thus generating a single final transformation. The multi-resolution scheme helps avoid local minima while achieving compu- tational efficiency and robustness. 6.3 Experimental Validation 6.3.1 Clinical Data 10 paired CCTA datasets (20 scans in total) obtained on a dual source CT (DSCT) scanner from 10 patients were evaluated in this retrospective analysis. The imaging protocol has previously been described in detail by Dey et al. [33]. All CCTA scans were obtained with the same 64- slice DSCT scanner (Definition; Siemens Medical Solutions, Erlangen, Germany) with gantry- rotation time of 330 milliseconds and a standard detector collimation of 0.6 mm. Eighty ml of intravenous contrast was administered before each scan. Raw data were reconstructed using 0.6 mm slice thickness, 0.3 mm slice increment, a 250250 mm field of view, a transverse field of view encompassing the whole heart (with z-axis heart dimensions ranging from 159-241 mm in the tested datasets), single-segment reconstruction and a medium-smooth reconstruction kernel (B26f). The time difference between baseline and follow-up scans ranged from 1 to 19 months. Patient characteristics and image parameters are described in Table 6.1. The best phase of the cardiac cycle for visualization of coronary arteries was determined by an expert reader at the time of clinical assessment, which ranged from 65% to 80%. The study was conducted according to the guidelines of the Cedars-Sinai Medical Center Institutional Review Board, and all patients gave written informed consent to the retrospective use of their data. 107 Table 6.1: Patient characteristics and image parameters of baseline and follow-up scans (for 5 cases only) Subject 1 Subject 2 Subject 3 Subject 4 Subject 5 Gender Male Male Male Female Male Age 79 64 82 79 62 Time difference (months) 4 15 11 5 1 Matrix size of baseline (pixels) 512512417 512512306 512512684 512512591 512512472 Matrix size of follow-up (pixels) 512512406 512512404 512512404 512512584 512512447 Pixel size of baseline (mm) 0.430.430.3 0.390.390.4 0.40.40.29 0.360.360.3 0.350.350.3 Pixel size of follow-up (mm) 0.470.470.3 0.450.450.3 0.450.450.3 0.360.360.3 0.360.360.3 Cardiac phase of baseline 70% 70% 70% 65% 65% Cardiac phase of follow-up 80% 70% 70% 65% 70% Image quality of baseline Good Good Excellent Excellent Excellent Image quality of follow-up Good Good Excellent Excellent Good Estimated maximal stenosis severity of baseline 25-49% 25-49% 25-49% 0 1-24% Plaque type at maximal stenosis of baseline Mixed Mixed Mixed 0 Calcified Plaque type at maximal stenosis of follow-up Mixed Mixed Non-calcified 0 Non-calcified 108 6.3.2 Evaluation of Registration Quality We compared the results obtained with our proposed method to the dense nonparametric MI based co-registration with same diffusion regularization and demons algorithm without any pre- processing. In our experiment, for MI based co-registration we set the number of histogram 10 and number of sample size to 10% of the total pixel number. We found the parameters that produce best results heuristically. The algorithm stops when iteration reaches the predefined iteration number 200 for all methods. (a) (b) (c) (d) Figure 6.3: Anatomical landmarks used for validation shown in one example case with images from a baseline scan in the top row, and images from a follow-up scan in the bottom row. We have obtained independent sets of anatomical landmarks from two clinical observers. Red crosses indicate the positions of landmarks, including the left main coronary artery origin (a), bifurcation of left anterior descending and left circumflex arteries (b), right coronary artery origin (c), and branch point of the first diagonal artery from the left anterior descending artery (d). Evaluation of the proposed co-registration algorithm was carried out with quantitative and qualitative analysis. To assess the accuracy of the co-registration algorithm quantitatively, we 109 calculated target registration error (TRE), as defined by Fitzpatrick et al. [46], which is the dis- tance between a certain point in the target volume and the corresponding point in the source volume after registration, and can be described as T RE(r)= 1 m m å i=1 D(p i t (r); p i s (T r)); (6.17) where p t (r) and p s (Tr) indicate coordinates of the landmark point from target volume and coor- dinate of the corresponding registered point from the source volume, respectively. D indicates the Euclidean distance between two points, and m denotes the number of points used for evaluations. Two expert observers independently selected several corresponding anatomical landmarks from each scan, and we calculated the mean and standard deviation of the pair-wise distances of these landmarks after automatic alignment after co-registration. We obtained several 3D landmarks at specific anatomical locations from each scan. These locations included the following anatomical locations: Left Main (LM) artery origin, bifurcation of the Left Anterior Descending (LAD) and Left Circumflex (LCX) arteries, Right Coronary Artery (RCA) origin and branch point of the First Diagonal (D1) from the LAD. Examples of those selected landmarks for one case are presented in Fig. 6.3 In addition, the alignment of serial scans was visually assessed qualitatively with the use of the ”roving window” technique, where with a moving window one could interactively “rove” and reveal the registered image at the same position. 6.4 Results Serial CCTA registrations were performed on an Intel Core2 Duo CPU with a clock speed of 2.5 GHz and a 4 GB memory. The mean computation time for the whole process (both global 110 0 5 10 15 20 25 30 35 40 45 50 2 4 6 8 10 12 14 x 10 -3 Iteration Energy Total Energy Volume Constraint Energy Figure 6.4: Comparison between total energy minimization without and with the volume- preserving constraint. displacement and local deformation) was 5.20.5 minutes. The mean computation times for MI based method and demons algorithm were 5.30.5 and 4.40.1 minutes, respectively. Our proposed method, MI based method and demons algorithm achieved a success rate of 90% as assessed visually (1 failed out of 10 datasets). The error of failed case was 5.4 mm. In this patient, there was no change in the coronary arteries, which were assessed to be normal. However, there was severe aneurysmal dilation of the aortic root and the ascending aorta which caused a significant change. We have therefore excluded the failed case in the following results. 111 (a) (b) (c) Figure 6.5: Comparison of global displacement and local deformation registration techniques: (a) the global displacement model, (b) the result of subsequent local deformation, and (c) the target image. Table 6.2: TRE comparison with different registration methods (Mean SD (mm)) TRE (mm) LM artery origin Bif. of LM and LCX RCA Origin Bif. of LAD and D1 Average Before registration 49.5141.71 49.9 41.01 50.5440.53 49.6441.42 49.5141.71 Global displacement 2.401.10 1.88 1.02 2.731.53 1.890.65 2.221.15 MI based method 1.930.89 1.46 0.84 2.361.53 1.570.58 1.831.06 Demons method 1.970.82 1.41 1.13 2.441.16 1.550.45 1.841.00 Proposed method without volume constraint 1.670.65 1.24 0.68 1.961.13 1.540.43 1.600.79 Proposed method without mask 1.950.89 1.44 1.12 2.231.23 1.520.47 1.781.00 Proposed method 1.650.59 1.21 0.68 1.751.02 1.640.51 1.560.74 In Fig. 6.4, we present graphs of the average total energy and the average energy derived by the volume preserving for 5 consecutive cases as a function of the number of iterations. In Table 6.2, we presented the TRE comparison of original distance before co-registration, global displacement and local deformation. In addition, we performed experiments without soft mask and volume constraint and with the MI-based method and with demons algorithm. To ob- tain these measurements, we considered each anatomical landmark as a 3D point and then trans- formed it using each of the co-registration methods. In all cases the original data misalignments 112 (a) (b) (c) (d) Figure 6.6: Example of registration results. The baseline image is registered to the subsequent image. The transverse orientation is shown in the top row while the coronal orientation is shown in the bottom row. The original baseline images before and after registration are shown in (a) and (b), respectively. Panel (c) shows the subsequent image. As shown in (a), the original image is significantly misaligned. The “roving window” (yellow) technique in (d) using a portion of the registered image from (b) to the target image in (c) shows that the result from our method is visually accurate. were larger than 40 mm. Our proposed method resulted in mean TRE of global displacement was 2.221.15 mm and further improved to 1.560.74 mm with subsequent local deforma- tion (p< 0:05) compared to 1.831.06 mm and 1.841.00 mm for MI based co-registration and demons algorithm, respectively (p< 0:05). The proposed method also had slightly better accuracy compared to the same methods without using soft mask. The importance of local de- formation is shown in Fig. 6.5, where images of global displacement correction and additional local deformation correction are presented. Although global displacement registration using a pre-segmented mask performed well in terms of and TRE, we still observed misalignment when visually comparing source and target images in Fig. 6.5. This visual misalignment was signifi- cantly reduced with local deformation correction. 113 Figure 6.7: Example of plaque lesion progression (worsening) in an 82-year old man. The top row shows a follow-up study while the bottom row shows a registered baseline study. The time difference between the two scans was 11 months. The pixel sizes of baseline and follow-up stud- ies were 0.40.40.29 mm and 0.450.450.3 mm, respectively. Green arrows show increased, non-calcified plaque and stenosis in the follow-up study. In addition, although we did not impose the symmetric consistency constraint, which guar- antees that the estimated forward transformation (source to target) equals that of inverse trans- formation (target to source) due to computational efficiency, we performed inverse consistency experiments using the landmarks obtained. To show the inverse consistency error, registered landmarks obtained using forward transformation (source to target) were transformed using in- verse mapping by replacing the source and target volumes. We then calculated inverse target registration errors between source landmarks and the inversely registered points. As shown in Table 6.3, the error was less than the voxel size of 0.15 mm, confirming that the transformation was consistent. 114 Figure 6.8: Example of a plaque without significant changes in a 59-year old man. The top row shows the follow-up study while the bottom row shows the registered baseline study. The lesion is marked with green arrows. The time difference between the two scans was 12 months. The pixel sizes of baseline and follow-up study were 0.310.310.3 mm and 0.340.340.3 mm, respectively. We also calculated the landmark distances of 100 point pairs (5 for each scan) between two observers to measure the inter-observer variability of the reference standard. Table 6.4 shows the statistics of observer variability. For qualitative assessments of image registration, the “roving window” technique was utilized as shown in Fig. 6.6. In Figs. 6.7 and 6.8, we show two examples of the final registration results for one case with observed change and for the other case without observed changes, respectively. In Fig. 6.7, we show images that have significant plaque progression in proximal LAD, where luminal diameter stenosis from mixed plaque increases from approximately 45% to90%, as visually assessed by an expert reader. In Fig. 6.8, images from the baseline and the follow-up CCTAs show no 115 Table 6.3: Inverse consistency error Inverse Consistency Error Mean SD (mm) LM artery origin 0.13 0.02 Bifurcation of LM and LCX 0.14 0.04 RCA Origin 0.12 0.10 Bifurcation of LAD and D1 0.17 0.06 Average 0.14 0.06 Table 6.4: Average landmark distances between two observers for each landmark (n=56) Observer Variability Mean SD (mm) Max (mm) Min (mm) LM artery origin 1.20 0.41 1.80 0.41 Bifurcation of LM and LCX 1.30 0.67 2.31 0.63 RCA Origin 1.33 0.54 2.12 0.42 Bifurcation of LAD and D1 1.42 0.75 2.45 0.44 Average 1.31 0.91 2.42 0.41 significant change in the proximal calcified plaque, which appeared to cause an estimated 25- 49% stenosis in both studies. Visual comparison can be made easily on co-registered images. 6.5 Discussion In this study, we proposed and evaluated a fully automated method to register serial CCTA scans for plaque evaluation. Although there have been some other studies which measure shape/volume changes or longitudinal changes [85, 118, 122, 130, 134] with different modalities, to the best of our knowledge, this is the first report applied to measuring longitudinal changes with CCTA. We have demonstrated that our two-step co-registration using feature mask is an effective approach for monitoring plaque development obtained from serial high-resolution CCTA scans. Combi- nation of an initial global displacement model and subsequent local deformation model is an effective approach that can register small structure of interest. A global displacement model was performed to find global transformation based on the segmented coronary tree structures from 116 CCTA scans. Subsequently, nonlinear transformation in a local deformation model using the dif- feomorphism approach was performed to find optimal local deformations, where a 3D feature mask (i.e. the soft mask) was introduced as a weighting function approximating the coronary vessel and surrounding region. The mask was also utilized for a volume-preserving constraint. The proposed method was validated on 10 paired CCTA datasets (20 scans in total) obtained on a dual source CT (DSCT) scanner. Robust and accurate target registration error was obtained with respect to the ground truth, as visually determined by two expert CCTA readers. The over- all co-registration errors were 2.221.15 mm for global displacement and 1.540.76 mm after further local deformation. In Comparison, the inter-observer variability was 1.310.91 mm. In addition, we have measured inverse consistency error and obtained 0.140.06 mm. Registration of the patient data achieved 90% success rate by visual confirmation. Sum of squared difference (similar to optical flow) was utilized as a similarity measure, since serial scans from the same modality were considered. We should note that in general, the SSD measure has been widely used for serial MR registration [58, 59] and CT registration [82, 84, 113]. To correct for the pos- sible image intensity discrepancies between the two serial scans, we performed additional step of histogram matching prior to local deformation model to ensure that both volumes have similar intensity distribution and satisfy the assumptions of SSD method. Additionally energy functional of the proposed co-registration algorithm is dependent on both intensity-based data driven term, SSD, and a shape-based term using soft mask. Thus shape-based term can compensate for the limitation using SSD, producing robust and accurate results. The MI measure might be suitable for this application as well since the intensities of corresponding voxel intensity may differ in 117 time. However we performed additional experiments using dense nonparametric MI-based co- registration and compared with the proposed method. MI approach (1.831.06 mm) appears to be slightly less accurate compared to the proposed method (1.560.74 mm). Serial co-registration is challenging in this case since the region of interest (coronary arteries) is small, coronary plaque size is small compared to CT spatial resolution, volume preservation is crucial for lesion comparison, and CCTA scans are acquired at two separate times with varia- tions in coverage, contrast intensity level, cardiac phase selected for optimal reconstruction, and different imaging artifacts. In addition, because of potential progression in coronary plaque, topo- logical equivalence between baseline and follow-up scans may be violated but volume must be preserved. To our knowledge, these challenges have not been previously studied for this appli- cation. To address these challenges, we exploited the binary mask obtained from pre-segmented coronary tree to further guide co-registration. However, this segmented feature-based approach will admittedly affect the co-registration accuracy, if the coronary tree segmentation is incorrect. In our algorithm, Gaussian kernel was applied to weaken the bias of segmentation in local defor- mation model. We observed that visual assessment of our method showed accurate co-registration of serial volumetric CCTA scans with one case failing due to severe aneurysmal dilation of the ascending aorta. Mean increase of 22% in non-calcified plaque volume measured from serial CCTA has been previously reported for the left main coronary artery [115]. A recent study [75] from our group re- ported mean measured plaque lengths of 10.56.2 mm in the proximal segments of the coronary arteries. The TRE after our proposed registration was 1.560.74 mm. If we consider proximal plaques with lengths similar to our reported study, this error is around 15% of the plaque length. However, we should note that even though visual assessment showed accurate registration, the 118 inter-observer variability in determining landmarks to determine TRE was 1.31 0.91 mm. Thus, it is most likely the TRE is exaggerated by the observer variability. Our approach has several potential limitations. First, quantitative evaluation of the registration is challenging, since there is no true rigid gold standard other than visual judgment, which is associated with inter-observer variability. To account for this, we measured target registration error using representative anatomical landmarks. Still, measured co-registration error can provide a relatively objective measure in mm but remains subject to variability between human observers. Second, the correspondences for the area of plaque at occlusion in the source volume are established mainly based on the regularization term since there is no corresponding region in the target volume due to the occlusion. In our variational framework, the regularization term is driven by diffusion regularization. However, such diffusion regularization may not be well suited for the phenomena of plaque progression, resulting in distortion at the vicinity of the plaque boundary. In addition to diffusion regularization, we introduced volume preserving term to guide desirable deformation field. However their roles are not completely independent and balance is desirable in our co-registration framework. Finally, there could be several physiologic factors that cause the change of volume. However, the objective of serial registration is to align same structure of interest (i.e. plaque area) that is located in lumen of artery without any artificial bending or deformation and thus it is reasonable to preserve original region of interest (i.e. entire coronary artery) in order to measure and monitor the progression or regression of plaque area in that region. The basic assumption made, therefore, is to preserve the volume of the entire coronary artery, in the source volume. To de-emphasize the strong bias of the binary mask obtained from segmenting the coronary lumen, we applied Gaussian kernel to the binary luminal mask. 119 Serial CCTA is clinically attractive because it holds promise for assessing changes in coronary artery plaque burden and response to therapy. With improvement in multi-detector CT technology and reduction in associated radiation dose, we expect efforts for improved assessment of change in coronary plaque to significantly increase in the near future. Rapid and effective CCTA co- registration can help facilitate improved evaluation of plaque change. 6.6 Conclusion We have developed a novel, fully automatic nonlinear feature-based volume co-registration algo- rithm for serial CCTA scans. Serial co-registration is challenging because of small target site size and differences in the field-of-view and coverage between the serial scans. Thus we incorporated segmented coronary tree structures as anatomical and geometric features for the co-registration that helps localize the structure of interest. Our proposed feature-based method provides accurate and robust co-registration with accuracy of 1.56 mm and 90% success rate as tested in datasets from 10 paired scans. Our results demonstrated the feasibility of the proposed co-registration algorithm for the analysis of serial plaque changes. 120 Chapter 7 : Conclusion and Future Work 7.1 Summary of the Research The segmentation and registration problems in clinical applications associated with cardiac dis- eases were studied in this research. We have presented novel segmentation and registration algo- rithms which were derived under a variational framework where these algorithms can incorporate cues available such as shape constraints, prior shape knowledge and intrinsic features. Especially, the level set method with implicit shape representation played a crucial role in solving segmenta- tion and reconstruction problems in Chapters 2 and 3, where the implicit shape representation is used to handle the complex anatomy of human hearts and their topological changes. In Chapter 2, we presented a novel segmentation method that jointly delineates the boundaries of epi- and endocardium of the left ventricle on the MR images under a variational framework using level sets. While most of traditional left ventricle segmentation methods incorporate the shape prior obtained by a training process to delineate the target object with a subtle boundary, the proposed method exploits a novel shape constraint as the implicit shape prior, which allows distance variation between epi- and endocardium in the Gaussian distribution. This prior offers a 121 statistical shape constraint by adjusting the distance distribution. The accuracy of this algorithm was demonstrated by experimental results. In Chapter 3, we proposed novel image integration technique for the use of catheter ablation. It enables reconstruction and registration simultaneously and achieves better performance than existing methods in terms of accuracy. Clinically, this technology has the potential to improve efficacy and safety of Atrial Fibrillation ablation by creating a more precise integration of EAM data and 3D imaging. Further studies using multiple data sets are needed to gauge its actual utility. In Chapter 4, we tackled another multi-modal registraiton problem where CT needs to be co- registered with MPS. Both modalities provide clinically important information for diagnosis of coronary artery disease. We presented a fully automated registration method utilizing geometric features from a reliable segmentation of gated MPS volumes, where regions of myocardium and blood pools were extracted and used as an anatomical mask to de-emphasize the inhomogeneities of intensity distribution caused by perfusion defects and physiological variations. In what follows, the best matching phase and nonlinear warping were performed for better defect visualization. In Chapter 5, a nonlinear ultrasound image registration method was proposed using the inten- sity and the local phase information under a variational framework. By integrating the intensity and the local phase information, we can find and track the nonlinear transformation of each pixel under diffeomorphism between the source and target images. In Chapter 6, we examined the problem of CCTA volume registration based on paired CCTA scans with nonlinear registration. We adopted a combined registration strategy to achieve this task. That is, we used the global rigid registration as an initialization, followed by the local regis- tration using nonlinear volume registration to meet a volume preserving constraint. We exploited 122 extracted coronary trees, as approximate binary masks and blur masks by convolving the Gaus- sian, to help localize and emphasize the region of interest. A volume preserving constraint was imposed so that the total volume of the binary mask before and after co-registration is preserved. Experimental results were given to demonstrate the performance of this technique. Throughout these studies, we were concentrated on the variational principle and its applica- tions to real-world cardiac applications. Sucessful results were provided by incorporating inher- ent cues such as shape constraint, geometric features, intrinsic features, etc. and we were able to obtain better results than numerous recent works with a lower complexity. 7.2 Future Research Directions To make our work more complete, we would like to extend current results along the following two directions. Detection and quantitation of soft plaques We would like to examine a larger patient group and perform direct voxel-wise quantifica- tion of plaque changes on CCTA volumes registered using our proposed technique. Since the soft plaques lie within the arterial wall, the lumen and the aterial wall should be delin- eated first and further detection is then required. This is however a challenging issue since the soft plaque has no salient feature but a similar intensity distribution to its surrounding areas. Multi-modal registration of CT/CTA with nuclear imaging obtained from hybrid scanner Multi-modal image integration of anatomical and functional data can be performed us- ing dedicated hybrid scanners or software registration. The respiratory effect and patient 123 movement may degrade the performance of hybrid scanners, which demand additional la- borious manual adjustment. In contrast, the software approach based on multi-modal image registration allows flexible, patient-specific and radiation dose-effective application. This problem is challenging since nuclear imaging such as SPECT and PET has low resolution without salient features. 124 References [1] S. Achenbach, “Computed tomography coronary angiography,” J Am Coll Cardiol, vol. 48, pp. 1919–1928, 2006. [2] ——, “Cardiac ct: State of the art for the detection of coronary arterial stenosis,” J Car- diovasc Comput Tomogr, vol. 1, no. 1, pp. 3–20, 2007. [3] S. Achenbach, D. Ropers, U. Hoffmann, B. MacNeill, U. Baum, K. Pohle, T. Brady, E. Pomerantsev, J. Ludwig, F. Flachskampf, S. Wicky, I. Jang, and W. Daniel, “Assess- ment of coronary remodeling in stenotic and nonstenotic coronary atherosclerotic lesions by multidetector spiral computed tomography,” J Am Coll Cardiol., vol. 43, pp. 842–847, 2004. [4] U. E. Aladl, G. A. Hurwitz, D. Dey, D. Levin, M. Drangova, and P. J. Slomka, “Auto- mated image registration of gated cardiac single-photon emission computed tomography and magnetic resonance imaging,” J Magn Reson Imaging, vol. 19, no. 3, pp. 283–90, 2004. [5] N. Amenta and M. Bern, “Surface reconstruction by voronoi filtering,” Discrete and Com- putational Geometry, vol. 22, no. 4, pp. 481–504, 1999. [6] A. A. Amini, T. E. Weymouth, and R. Jain, “Using dynamic programming for solving variational problems in vision,” IEEE Trans Pattern Anal Mach Intell, vol. 12, no. 9, pp. 855–867, 1990. [7] D. S. Berman, X. Kang, S. W. Hayes, J. D. Friedman, I. Cohen, A. Abidov, L. J. Shaw, A. M. Amanullah, G. Germano, and R. Hachamovitch, “Adenosine myocardial perfusion single-photon emission computed tomography in women compared with men. impact of diabetes mellitus on incremental prognostic value and effect on patient management,” J Am Coll Cardiol, vol. 41, pp. 1125–1133, 2003. [8] D. S. Berman, H. Kiat, J. D. Friedman, F. P. Wang, K. van Train, L. Matzer, J. Maddahi, and G. Germano, “Separate acquisition rest thallium-201/stress technetium-99m sestamibi dual-isotope myocardial perfusion single-photon emission computed tomography: a clini- cal validation study,” J Am Coll Cardiol, vol. 22, no. 5, pp. 1455–1464, 1993. [9] D. S. Berman, L. J. Shaw, R. Hachamovitch, J. D. Friedman, D. M. Polk, S. W. Hayes, L. E. J. Thomson, G. Germano, N. D. Wong, K. Xingping, and A. Rozanski, “Compara- tive use of radionuclide stress testing, coronary artery calcium scanning, and noninvasive coronary angiography for diagnostic and prognostic cardiac assessment,” Seminars in Nu- clear Medicine, vol. 37, no. 1, pp. 2–16, 2007. 125 [10] P. J. Besl and N. D. McKay, “A method for registration of 3-d shapes,” IEEE Trans Pattern Anal Mach Intell, vol. 14, no. 2, pp. 239–256, 1992. [11] R. Blankstein, D. R. Okada, J. A. Rocha-Filho, F. J. Rybicki, T. J. Brady, and R. C. Cury, “Cardiac myocardial perfusion imaging using dual source computed tomography,” Int J Cardiovasc Imaging, vol. 25, no. Suppl 2, pp. 209–216, 2009. [12] F. Bookstein, “Principal warps: Thin-plate splines and the decomposition of deforma- tions,” IEEE Trans Pattern Anal Mach Intell, vol. 11, no. 6, 1989. [13] D. Boukerroui, A. Noble, and M. Brady, “On the choice of band-pass quadrature filters,” Journal of Mathematical Imaging and Vision, vol. 21, no. 1, pp. 53–80, 2004. [14] Y . Boykov and V . Kolmogorov, “Computing geodesic and minimal surfaces via graph cuts,” in International Conference on Computer Vision, 2003, pp. 26–33. [15] Y . Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” IEEE Trans Pattern Anal Mach Intell, vol. 23, no. 11, pp. 1222–1239, 2001. [16] T. Brox, A. Bruhn, N. Papenberg, and J. Weickert, “High accuracy optical flow estimation based on a theory for warping,” in European Conference on Computer Vision, ser. LNCS, vol. 3024. Prague, Czech Republic: Springer, May 2004, pp. 25–36. [17] A. Bruhn, J. Weickert, C. Feddern, T. Kohlberger, and C. Schnorr, “Variational optical flow computation in real time,” IEEE Trans Image Process, vol. 14, pp. 608– 615, 2005. [18] C. Burgstahler, A. Reimann, T. Beck, A. Kuettner, D. Baumann, M. Heuschmid, H. Bro- doefel, C. D. Claussen, A. F. Kopp, and S. Schroeder, “Influence of a lipid-lowering ther- apy on calcified and noncalcified coronary plaques monitored by multislice detector com- puted tomography: results of the new age ii pilot study,” Invest Radiol., vol. 42, pp. 189–95, 2007. [19] Y . Cao, M. Miller, R. Winslow, and L. Younes, “Large deformation diffeomorphic metric mapping of vector fields,” IEEE Trans Med Imaging, vol. 24, no. 9, pp. 1216–1230, 2005. [20] V . Chalana, D. Linker, D. Haynor, and Y . Kim, “A multiple active contour model for cardiac boundary detection on echocardiographic sequences,” IEEE Trans Med Imaging, vol. 15, no. 3, pp. 290–298, 1996. [21] T. Chan and L. Vese, “Active contours without edges,” IEEE Trans Image Process, vol. 10, no. 2, pp. 266–277, 2001. [22] V . Cheng, R. Nakazato, D. Dey, S. Gurudevan, J. Tabak, M. Budoff, R. Karlsberg, J. Min, and D. Berman, “Reproducibility of coronary artery plaque volume and composition quan- tification by 64-detector row coronary computed tomographic angiography: An intraob- server, interobserver, and interscan variability study,” J Cardiovasc Comput Tomogr, 2009 (in press). [23] D. Chetverikov, D. Svirko, D. Stepanov, and P. Krsek, “The trimmed iterative closest point algorithm,” International Conference on Pattern Recognition, vol. 03, p. 30545, 2002. 126 [24] G. Christensen and H. Johnson, “Consistent image registration,” IEEE Trans Med Imaging, vol. 20, pp. 1384–1397, 2001. [25] G. Christensen, R. Rabbitt, and M. Miller, “Deformable templates using large deformation kinematics,” IEEE Trans on Image Process, vol. 5, no. 10, pp. 1435–1447, 1996. [26] C.-H. Chuang and C.-C. J. Kuo, “A wavelet descriptor of planar curves: theory and appli- cations,” IEEE Trans on Image Process, vol. 5, no. 1, pp. 56–70, 1996. [27] D. Collins, T. Peters, and A. Evans, “Automated 3d nonlinear deformation procedure for determination of gross morphometric variability in human brain,” in SPIE, 1994, pp. 180– 190. [28] T. F. Cootes, C. J. Taylor, D. H. Cooper, and J. Graham, “Active shape models-their training and application,” Computer Vision and Image Understanding, vol. 61, no. 1, pp. 38–59, 1995. [29] D. Cremers, S. J. Osher, and S. Soatto, “Kernel density estimation and intrinsic align- ment for shape priors in level set segmentation,” International Journal of Computer Vision, vol. 69, no. 3, 2006. [30] W. R. Crum, R. I. Scahill, and N. C. Fox, “Automated hippocampal segmentation by re- gional fluid registration of serial mri: validation and application in alzheimer’s disease,” Neuroimage, vol. 13, no. 5, pp. 847–855, 2001. [31] E. D’Agostino, F. Maes, D. Vandermeulen, and P. Suetens, “A viscous fluid model for mul- timodal non-rigid image registration using mutual information,” Medical Image Analysis, vol. 7, no. 4, pp. 565–575, 2003. [32] M. De Craene, A. du Bois dAische, B. Macq, and S. Warfield, “Incorporating metric flows and sparse jacobian transformations in itk,” The Insight Journal. URL http://hdl. handle. net/1926/183. [33] D. Dey, C. J. Lee, M. Ohba, A. Gutstein, P. J. Slomka, V . Cheng, Y . Suzuki, S. Suzuki, A. Wolak, L. L. Meunier, L. E. J. Thomson, I. Cohen, J. D. Friedman, G. Germano, and D. S. Berman, “Image quality and artifacts in coronary ct angiography with dual-source ct: Initial clinical experience,” J Cardiovasc Comput Tomogr, vol. 2, pp. 105–114, 2008. [34] J. Dong, H. Calkins, S. B. Solomon, S. Lai, D. Dalal, A. Lardo, E. Brem, A. Preiss, R. D. Berger, H. Halperin, and T. Dickfeld, “Integrated Electroanatomic Mapping With Three- Dimensional Computed Tomographic Images for Real-Time Guided Ablations,” Circula- tion, vol. 113, no. 2, pp. 186–194, 2006. [35] J. S. Duncan and N. Ayache, “Medical image analysis: progress over two decades and thechallenges ahead,” IEEE Trans Pattern Anal Mach Intell, vol. 22, no. 1, pp. 85–106, 2000. [36] H. Edelsbrunner, “Shape reconstruction with delaunay complex,” in LATIN ’98: Proceed- ings of the Third Latin American Symposium on Theoretical Informatics. London, UK: Springer-Verlag, 1998, pp. 119–132. 127 [37] P. J. Edwards, D. L. G. Hill, J. A. Little, and D. J. Hawkes, “A three-component deforma- tion model for image-guided surgery,” Medical Image Analysis, vol. 2, no. 4, pp. 355–367, 1998. [38] H. K. Eltzschig and R. Ehlers, “Noninvasive tests in patients with stable coronary artery disease,” N Engl J Med., vol. 345, no. 18, p. 1351, 2001. [39] T. L. Faber, C. A. Santana, E. V . Garcia, J. Candell-Riera, R. D. Folks, J. W. Peifer, A. Hop- per, S. Aguade, J. Angel, and J. L. Klein, “Three-dimensional fusion of coronary arteries with myocardial perfusion distributions: Clinical validation,” J Nucl Med, vol. 45, no. 5, 2004. [40] E. Falk, “Pathogenesis of atherosclerosis,” Journal of the American College of Cardiology, vol. 47, no. 8, Supplement 1, pp. C7–C12, 2006. [41] M. Felsberg and G. Sommer, “A new extension of linear signal processing for estimating local properties and detecting features,” in DAGM, 2000, pp. 195–202. [42] M. Ferencik, D. Ropers, S. Abbara, R. C. Cury, U. Hoffmann, K. Nieman, T. J. Brady, F. Moselewski, W. G. Daniel, and S. Achenbach, “Diagnostic accuracy of image postpro- cessing methods for the detection of coronary artery stenoses by using multidetector ct,” Radiology, vol. 243, pp. 696–702, 2007. [43] B. Fischer and J. Modersitzki, “Fast diffusion registration,” AMS Contemporary Mathe- matics, Inverse Problems, Image Analysis, and Medical Imaging, vol. 313, pp. 117–129, 2002. [44] ——, “Curvature based image registration,” Journal of Mathematical Imaging and Vision, vol. 18, no. 1, pp. 81–85, 2003. [45] ——, “A unified approach to fast image registration and a new curvature based registration technique,” Linear Algebra and its Applications, vol. 23, no. 7, pp. 107–124, 2004. [46] J. M. Fitzpatrick, J. B. West, and C. R. Maurer, “Predicting error in rigid-body point-based registration,” IEEE Trans Med Imaging, vol. 17, pp. 694–702, 1998. [47] P. Foroughi and P. Abolmaesumi, “Elastic registration of 3d ultrasound images,” in Medical Image Computing and Computer Assisted Intervention (MICCAI), vol. 3749. Springer, 2005, pp. 83–90. [48] A. F. Frangi, W. J. Niessen, and M. A. Viergever, “Three-dimensional modeling for func- tional analysis of cardiacimages, a review,” IEEE Trans on Med Imaging, vol. 20, no. 1, pp. 2–26, 2001. [49] O. Gaemperli, T. Schepis, V . Kalff, M. Namdar, I. Valenta, L. Stefani, L. Desbiolles, S. Leschka, L. Husmann, H. Alkadhi, and P. A. Kaufmann, “Validation of a new cardiac image fusion software for three-dimensional integration of myocardial perfusion spect and stand-alone 64-slice ct angiography,” Eur J Nucl Med Mol Imaging, vol. 34, no. 7, pp. 1097–1106, 2007. 128 [50] O. Gaemperli, T. Schepis, I. Valenta, L. Husmann, H. Scheffel, V . Duerst, F. R. Eberli, T. F. Luscher, H. Alkadhi, and P. A. Kaufmann, “Cardiac image fusion from stand-alone spect and ct: Clinical experience,” J Nucl Med, vol. 48, no. 5, pp. 696–703, 2007. [51] G. Germano, J. Erel, H. Lewin, P. B. Kavanagh, and D. S. Berman, “Automatic quantitation of regional myocardial wall motion and thickening from gated technetium-99m sestamibi myocardial perfusion single-photon emission computed tomography,” J Am Coll Cardiol, vol. 30, pp. 1360–1367, 1997. [52] G. Germano, P. B. Kavanagh, and D. S. Berman, “An automatic approach to the analy- sis, quantitation and review of perfusion and function from myocardial perfusion spect images,” Int J Card Imaging, vol. 13, pp. 337–346, 1997. [53] G. Germano, H. Kiat, P. B. Kavanagh, M. Moriel, M. Mazzanti, H. T. Su, K. F. V . Train, and D. S. Berman, “Automatic quantification of ejection fraction from gated myocardial perfusion spect,” J Nucl Med, vol. 36, no. 11, pp. 2138–47, 1995. [54] A. Giachetti, “On-line analysis of echocardiographic image sequences,” Medical Image Analysis, vol. 2, no. 3, pp. 261–284, 1998. [55] V . Grau, H. Becher, and J. A. Noble, “Phase-based registration of multi-view real-time three-dimensional echocardiographic sequences,” in Medical Image Computing and Com- puter Assisted Intervention (MICCAI), 2006, pp. 612–619. [56] C. Guetter, M. Wacker, C. Xu, and J. Hornegger, “Registration of cardiac spect/ct data through weighted intensity co-occurrence priors,” Medical Image Computing and Com- puter Assisted Intervention (MICCAI), vol. 4791, 2007. [57] E. Haber and J. Modersitzki, “Numerical methods for volume preserving image registra- tion,” Inverse Problems, vol. 20, pp. 1621 – 1638, 2004. [58] J. Hajnal, N. Saeed, A. Oatridge, E. Williams, I. Young, and G. Bydder, “Detection of subtle brain changes using subvoxel registration and subtraction of serial MR images,” Journal of computer assisted tomography, vol. 19, no. 5, pp. 677–691, 1995. [59] J. Hajnal, N. Saeed, E. Soar, A. Oatridge, I. Young, and G. Bydder, “A registration and interpolation procedure for subvoxel matching of serially acquired MR images,” Journal of computer assisted tomography, vol. 19, no. 2, p. 289, 1995. [60] C. L. Hansen, R. A. Goldstein, O. O. Akinboboye, D. S. Berman, E. H. Botvinick, K. B. Churchwell, C. D. Cooke, J. R. Corbett, S. J. Cullom, S. T. Dahlberg, R. S. Druz, E. P. Ficaro, J. R. Galt, R. K. Garg, G. Germano, G. V . Heller, M. J. Henzlova, M. C. Hyun, L. L. Johnson, A. Mann, B. D. McCallister, R. A. Q. Jr., T. D. Ruddy, S. N. Sundaram, R. Taillefer, R. P. Ward, and J. J. Mahmarian, “Myocardial perfusion and function: single photon emission computed tomography,” J Nucl Cardiol, vol. 14, pp. e39–60, 2007. [61] T. Hartkens, D. L. G. Hill, A. D. Castellano-Smith, C. R. Maurer, A. J. Martin, W. A. Hall, H. Liu, and C. L. Truwit, “Measurement and analysis of brain deformation during neurosurgery,” IEEE Trans Med Imaging, vol. 22, no. 1, pp. 82–92, 2003. 129 [62] E. Heist, J. Chevalier, G. Holmvang, and et al., “Factors affecting error in integration of electroanatomic mapping with ct and mr imaging during catheter ablation of atrial fibril- lation,” Journal of Interventional Cardiac Electrophysiology, vol. 17, no. 1, pp. 21–27, October 2006. [63] A. Hennemutha, T. Boskampa, D. Fritzc, C. Kuhnela, S. Bocka, D. Rinckb, M. Scheuer- ingb, and H.-O. Peitgen, “One-click coronary tree segmentation in ct angiographic im- ages,” Computer Assisted Radiology and Surgery, vol. 1281, pp. 317–321, 2005. [64] I. Herlin and N. Ayache, “Features extraction and analysis methods for sequences of ultra- sound images,” Image and Vision Computing, vol. 10, no. 10, pp. 673–682, 1992. [65] C. K. Hoh, M. Dahlbom, G. Harris, Y . Choi, R. A. Hawkins, M. E. Phelps, and J. Mad- dahi, “Automated iterative three-dimensional registration of positron emission tomography images,” J Nucl Med, vol. 34, pp. 2009–2018, 1993. [66] M. Holden, “A review of geometric transformations for nonrigid body registration,” IEEE Trans Med Imaging, vol. 27, no. 1, pp. 111–128, 2008. [67] H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, and W. Stuetzle, “Surface reconstruction from unorganized points,” Computer Graphics, vol. 26, no. 2, pp. 71–78, 1992. [Online]. Available: citeseer.ist.psu.edu/article/hoppe92surface.html [68] B. Horn and B. Schunck, “Determining optical flow,” Artificial Intelligence, vol. 17, pp. 185–203, 1981. [69] L. Ibanez, W. Schroeder, L. Ng, and J. Cates, The ITK Software Guide: The Insight Seg- mentation and Registration Toolkit. Kitware, 2003. [70] M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: Active contour models,” International Journal of Computer Vision, vol. 1, no. 4, pp. 321–331, 1998. [71] V . Kolmogorov, “Graph-based algorithms for multi-camera reconstruction problem,” Ph.D. dissertation, Cornell University, 2003. [72] V . Kolmogorov and R. Zabih, “What energy functions can be minimized via graph cuts?” IEEE Trans Pattern Anal Mach Intell, vol. 26, no. 2, pp. 147–159, 2004. [73] J. Krucker, G. LeCarpentier, J. Fowlkes, and P. Carson, “3d spatial compounding of ul- trasound images using image-based nonrigid registration,” Ultrasound Med Biol., vol. 26, no. 9, pp. 1475–1488, 2000. [74] J. Kybic and M. Unser, “Fast parametric elastic image registration,” IEEE Trans Med Imag- ing, vol. 12, no. 11, pp. 1427–1442, 2003. [75] L. Le Meunier, P. Slomka, D. Dey, V . Cheng, L. Thomson, S. Hayes, J. Friedman, G. Ger- mano, and D. Berman, “Feasibility of coronary plaque imaging with FDG high-definition PET: Phantom study,” in Society of Nuclear Medicine Annual Meeting Abstracts, vol. 49, no. Supplement 1. Soc Nuclear Med, 2008, p. 188P. 130 [76] A. W. Leber, A. Becker, A. Knez, F. V . Ziegler, M. Sirol, K. Nikolaou, B. Ohnesorge, Z. A. Fayad, C. R. Becker, M. Reiser, G. Steinbeck, and P. Boekstegers, “Accuracy of 64-slice computed tomography to classify and quantify plaque volumes in the proximal coronary system: A comparative study using intravascular ultrasound,” J Am Coll Cardiol., vol. 47, no. 3, pp. 672–677, 2006. [77] A. W. Leber, A. Knez, A. Becker, C. Becker, F. V . Ziegler, K. Nikolaou, C. Rist, M. Reiser, C. White, G. Steinbeck, and P. Boekstegers, “Accuracy of multidetector spiral computed tomography in identifying and differentiating the composition of coronary atherosclerotic plaques: a comparative study with intracoronary ultrasound,” J Am Coll Cardiol., vol. 43, no. 7, pp. 1241–1247, 2004. [78] M. J. Ledesma-Carbayo, J. Kybic, M. Desco, and et al., “Spatio-temporal nonrigid regis- tration for ultrasound cardiac motion estimation,” IEEE Trans Med Imaging, vol. 24, no. 9, pp. 1113–1126, September 2005. [79] H. Lester and S. R. Arridge, “A survey of hierarchical non-linear medical image registra- tion,” Pattern Recognition, vol. 32, no. 1, pp. 129–149, 1999. [80] K. K. Leung, M. Holden, N. Saeed, K. J. Brooks, J. B. Buckton, A. A. Williams, S. P. Campbell, K. Changani, D. G. Reid, Y . Zhao, M. Wilde, D. Rueckert, J. V . Hajnal, and D. L. Hill, “Automatic quantification of changes in bone in serial mr images of joints,” IEEE Trans Med Imaging, vol. 25, no. 12, pp. 1617–1626, 2006. [81] M. Leventon, W. Grimson, and O. Faugeras, “Statistical shape influence in geodesic active contours,” in International Conference on Computer Vision, 2000, pp. 316–323. [82] B. Li, G. Christensen, E. Hoffman, G. McLennan, and J. Reinhardt, “Establishing a nor- mative atlas of the human lung intersubject warping and registration of volumetric CT images,” Academic radiology, vol. 10, no. 3, pp. 255–265, 2003. [83] T. Liu, D. Shen, and C. Davatzikos, “Deformable registration of cortical structures via hybrid volumetric and surface warping,” NeuroImage, vol. 22, pp. 1790–1801, 2004. [84] W. Lu, M. Chen, G. Olivera, K. Ruchala, and T. Mackie, “Fast free-form deformable registration via calculus of variations,” Physics in Medicine and Biology, vol. 49, no. 14, pp. 3067–3087, 2004. [85] T. R. Mackie, J. Kapatoes, K. Ruchala, W. Lu, C. Wu, G. Olivera, L. Forrest, W. Tome, J. Welsh, R. Jeraj, P. Harari, P. Reckwerdt, B. Paliwal, M. Ritter, H. Keller, J. Fowler, and M. Mehta, “Image guidance for precise conformal radiotherapy,” Int J Radiat Oncol Biol Phys, vol. 56, pp. 89–105, 2003. [86] F. Maes, A. Collignon, D. Vandermeulen, and et al., “Multimodality image registration by maximization of mutual information,” IEEE Trans Med Imaging, vol. 16, pp. 187–198, 1997. [87] T. Makela, P. Clarysse, O. Sipila, N. Pauna, Q. C. Pham, T. Katila, and I. E. Magnin, “A review of cardiac image registration methods,” IEEE Trans Med Imaging, vol. 21, no. 9, pp. 1011–1021, 2002. 131 [88] Z. Malchano, “Image guidance in cardiac electrophysiology,” Master’s thesis, Mas- sachusette Institute of Technology, Boston, MA, June 2006. [89] T. Masuda, K. Sakaue, and N. Yokoya, “Registration and integration of multiple range images for 3-d model construction,” in Pattern Recognition, 1996., Proceedings of the 13th International Conference on, 1996, pp. 879–883. [90] M. Mellor, “Phase methods for non-rigid medical image registration,” Ph.D. dissertation, Oxford University, 2004. [91] M. Mellor and M. Brady, “Non-rigid multimodal image registration using local phase,” in Medical Image Computing and Computer Assisted Intervention (MICCAI), vol. 3216. Springer, 2005, pp. 789–796. [92] A. Montillo, D. Metaxas, and L. Axel, “Automated segmentation of the left and right ven- tricles in 4d cardiac spamm images,” Medical Image Computing and Computer Assisted Intervention (MICCAI), vol. 2488, pp. 620–633, 2002. [93] M. C. Morrone and R. A. Owens, “Feature detection from local energy,” Pattern Recogn. Lett., vol. 6, no. 5, pp. 303–313, 1987. [94] M. Mulet-Parada and A. Noble, “2d+t acoustic boundary detection in echocardiography,” in Medical Image Computing and Computer Assisted Intervention (MICCAI), vol. 1496. Springer, 1998, pp. 806–813. [95] D. Mumford and J. Shah, “Optimal approximation by piecewise smooth functions and associated variational problems,” Commun. Pure Appl. Math, vol. 42, pp. 577–685, 1989. [96] H. Nakajo, S. Kumita, K. Cho, and T. Kumazaki, “Three-dimensional registration of my- ocardial perfusion spect and ct coronary angiography,” Ann Nucl Med, vol. 19, pp. 207– 215, 2005. [97] T. R. Nelson, D. B. Downay, D. H. Pretorius, and A. Fenster, Three-Dimensional Ultra- sound. Philadelphia, PA: Lippincott Williams and Wilkins, 1990. [98] C. Pappone, S. Rosanio, G. Orte, and et al., “Circumferential radiofrequency ablation of pulmonary vein ostia,” circulation, vol. 102, pp. 2619–2628, November 2000. [99] N. Paragios, “Variational methods and partial differential equations in cardiac image anal- ysis,” in IEEE International Symposium on Biomedical Imaging, 2004, pp. 17–20. [100] ——, “A level set approach for shape-driven segmentation and tracking of the left ventri- cle,” IEEE Trans Med Imaging, vol. 22, no. 6, pp. 773–776, 2003. [101] M. Picard, R. Popp, and A. Weyman, “Assessment of left ventricular function by echocar- diography: a technique in evolution,” Journal of the American Society of Echocardiogra- phy, vol. 21, no. 3, pp. 14–21, 2002. [102] J.-P. Pons, “Methodological and applied contributions to the deformable models frame- work,” Ph.D. dissertation, ENPC, 2005. 132 [103] K. Pulli, “Multiview registration for large data sets,” in Inlt. Conf. on 3D Digital Imaging and Modeling, 1999, pp. 160–168. [104] I. B. Ray and E. Heist, “Treating atrial fibrillation. what is the consensus now?” Postgrad- uate medicine, vol. 118, no. 4, October 2005. [105] V . Reddy, Z. Malchano, G. Holmvang, and et al., “Integration of cardiac magnetic reso- nance imaging with three-dimensional electroanatomic mapping to guide left ventricular catheter manipulation: feasibility in a porcine model of healed myocardial infarction,” J Am Coll Cardiol, vol. 44, no. 11, pp. 2202–2213, 2004. [106] D. Rey, G. Subsol, H. Delingette, and N. Ayache, “Automatic detection and segmentation of evolving processes in 3d medical images: Application to multiple sclerosis,” Medical Image Analysis, vol. 6, no. 2, pp. 163–179, 2002. [107] D. Rogers, An Introduction to NURBS: With Historic Perspective. San Francisco, CA: Morgan Kaufmann Publisher Inc., 2001. [108] T. Rohlfing, C. R. Maurer, D. A. Bluemke, and M. A. Jacobs, “V olume-preserving nonrigid registration of mr breast images using free-form deformation with an incompressibility constraint,” IEEE Trans Med Imaging, vol. 22, pp. 730–741, 2003. [109] K. Rohr, Landmark-Based Image Analysis: Using Geometric and Intensity Models. Boston, MA: Kluwer, 2001. [110] K. Rohr and H. S. et al., “Landmark-based elastic registration using approximating thin- plate splines,” IEEE Trans Med Imaging, vol. 20, no. 6, pp. 526–534, 2001. [111] M. Rousson, N. Paragios, and R. Deriche, “Implicit active shape models for 3d segmen- tation in mri imaging,” Medical Image Computing and Computer Assisted Intervention (MICCAI), vol. 2217, pp. 209–216, 2007. [112] S. Rusinkiewicz and M. Levoy, “Efficient variants of the ICP algorithm,” in Proceedings of the Third Intl. Conf. on 3D Digital Imaging and Modeling, 2001, pp. 145–152. [113] D. Sarrut, V . Boldea, S. Miguet, and C. Ginestet, “Simulation of four-dimensional CT images from deformable registration between inhale and exhale breath-hold CT scans,” Medical physics, vol. 33, p. 605, 2006. [114] W. M. Schaefer, C. S. Lipke, B. Nowak, H. J. Kaiser, P. Reinartz, A. Buecker, G. A. Krombach, U. Buell, and H. P. Kuhl, “Validation of qgs and 4d-mspect for quantification of left ventricular volumes and ejection fraction from gated 18f-fdg pet: comparison with cardiac mri,” J Nucl Med, vol. 45, pp. 74–79, 2004. [115] M. Schmid, S. Achenbach, D. Ropers, S. Komatsu, U. Ropers, W. G. Daniel, and T. Pfled- erer, “Assessment of changes in non-calcified atherosclerotic plaque volume in the left main and left anterior descending coronary arteries over time by 64-slice computed to- mography,” J Am Coll Cardiol., vol. 101, pp. 579–584, 2008. 133 [116] M. Sdika, “A fast nonrigid image registration with constraints on the jacobian using large scale constrained optimization,” IEEE Trans Med Imaging, vol. 27, no. 2, pp. 271–281, 2008. [117] J. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Compu- tational Geometry. Cambridge University Press, 1998. [118] J. Shi, B. Sahiner, H. P. Chan, L. Hadjiiski, C. Zhou, P. N. Cascade, N. Bogot, E. A. Kazerooni, Y . T. Wu, and J. Wei, “Pulmonary nodule registration in serial ct scans based on rib anatomy and nodule template matching,” Med Phys, vol. 34, pp. 1336–1347, 2007. [119] K. K. Shung, Diagnostic Ultrasound: Imaging and Blood Flow Measurements. Boca Raton, FL: Francis Taylor: CRC Press, 2005. [120] K. K. Shung, M. B. Smith, and B. M. W. Tsui, Principles of Medical Imaging. San Diego, CA: Academic, 1992. [121] P. J. Slomka, “Software approach to merging molecular with anatomic information,” Jour- nal of Nuclear Medicine, vol. 45, no. 1, pp. 36–45, 2004. [122] P. J. Slomka and R. P. Baum, “Multimodality image registration with software: state-of- the-art,” Eur J Nucl Med Mol Imaging, vol. 36, no. Suppl 1, pp. 44–55, 2009. [123] P. J. Slomka, D. S. Berman, and G. Germano, “Applications and software techniques for integrated cardiac multimodality imaging,” Expert Review of Cardiovascular Therapy, vol. 6, no. 1, pp. 27–41, 2008. [124] P. J. Slomka, J. Mandel, D. Downey, and A. Fenster, “Evaluation of voxel-based registra- tion of 3-d power doppler ultrasound and 3-d magnetic resonance angiographic images of carotid arteries,” Ultrasound in medicine & biology, vol. 27, no. 7, pp. 875–892, 2001. [125] P. J. Slomka, H. Nishina, D. S. Berman, X. Kang, C. Akincioglu, J. D. Friedman, S. W. Hayes, U. E. Aladl, and G. Germano, ““motion-frozen” display and quantification of my- ocardial perfusion,” J Nucl Med, vol. 45, pp. 1128–34, 2004. [126] S. Sola and M. Y . Desai, “Measuring coronary artery plaque volume and composition by mdct: Too much, too soon? (in press),” J Cardiovasc Comput Tomogr, 2009. [127] G. Storvik, “A bayesian approach to dynamic contours through stochastic sampling and simulated annealing,” IEEE Trans Pattern Anal Mach Intell, vol. 16, no. 10, pp. 976–986, 1994. [128] C. Studholme, D. L. G. Hill, and D. J. Hawkes, “An overlap invariant entropy measure of 3d medical image alignment,” Pattern Recognition, vol. 32, pp. 71–86, 1999. [129] R. Szeliski, “Bayesian modeling of uncertainty in low-level vision,” International Journal of Computer Vision, vol. 5, no. 3, pp. 271–301, 1990. [130] H. Takao, I. Doi, and M. Tateno, “Evaluation of an automated system for temporal sub- traction of thin-section thoracic ct,” Br J Radiol, vol. 80, pp. 85–89, 2007. 134 [131] J. P. Thirion, “Image matching as a diffusion process: an analogy with maxwell’s demons,” Medical Image Analysis, vol. 2, no. 3, pp. 243–260, 1998. [132] A. Tsai, A. Yezzi, W. Wells, C. Tempany, D. Tucker, A. Fan, W. E. Grimson, and A. Will- sky, “A shape-based approach to the segmentation of medical imagery using level sets,” IEEE Trans Med Imaging, vol. 22, no. 2, pp. 137–154, 2003. [133] P. Viola and W. M. Wells, “Alignment by maximization of mutual information,” Interna- tional Journal of Computer Vision, vol. 24, no. 2, pp. 137–154, 1997. [134] H. Wang, L. Dong, J. O’Daniel, R. Mohan, A. S. Garden, K. K. Ang, D. A. Kuban, M. Bon- nen, J. Y . Chang, and R. Cheung, “Validation of an accelerated ’demons’ algorithm for de- formable image registration in radiation therapy,” Phys Med Biol, vol. 50, pp. 2887–2905, 2005. [135] J. Woo, B.-W. Hong, D. Dey, V . Cheng, A. Ramesh, G. Sundaramoorthi, C.-C. J. Kuo, D. S. Berman, G. Germano, and P. J. Slomka, “Feature-based non-rigid volume registration of serial coronary ct angiography,” in SPIE, Orlando, FL, Feb. 2009. [136] J. Woo, B. Hong, S. Kumar, I. B. Ray, and C.-C. J. Kuo, “Joint reconstruction and regis- tration using level sets: Application to the computer-aided ablation of atrial fibrillation,” in International Conference Frontiers in the Convergence of Bioscience and Information Technologies, Jeju, Korea, Oct 2007. [137] J. Woo, P. J. Slomka, D. Dey, V . Cheng, A. Ramesh, B.-W. Hong, D. S. Berman, C.- C. J. Kuo, and G. Germano, “Automated multi-modality registration of 64-slice coronary ct angiography with myocardial perfusion spect,” in IEEE International Symposium on Biomedical Imaging, Boston, MA, June 2009. [138] A. Yezzi, A. Tsai, and A. Willsky, “A fully global approach to image segmentation via coupled curve evolution equations,” Journal of Visual Communication and Image Repre- sentation, vol. 13, pp. 195–216, 2002. [139] A. Yezzi, L. Zollei, and T. Kapur, “A variational framework for integrating segmentation and registration through active contours,” Medical Image Analysis, vol. 7, no. 2, pp. 171– 185, 2003. [140] P. A. Yushkevich, J. Piven, H. Cody, S. Ho, J. C. Gee, and G. Gerig, “User-guided level set segmentation of anatomical structures with ITK-SNAP,” Insight Jounral, vol. 1, 2005, special Issue on ISC/NA-MIC/MICCAI Workshop on Open-Source Software. [141] S. Yusuf, S. Reddy, S. Ounpuu, and S. Anand, “Global burden of cardiovascular diseases: part i: general considerations, the epidemiologic transition, risk factors, and impact of urbanization,” Circulation, vol. 104, no. 22, pp. 2746–2753, 2001. [142] H. Zhao, T. Chan, B. Merriman, and S. Osher, “A variational level set approach to multi- phase motion,” J. Comput. Phys., vol. 127, pp. 179–195, 1996. 135 [143] H. Zhao, S. Osher, and R. Fedkiw, “Fast surface reconstruction using the level set method,” 1st IEEE Workshop on Variational and Level Set Methods, 8th ICCV , Vancouver, pp. 194– 202, 2001. [144] D. Zikic, W. Wein, A. Khamene, D. Clevert, and N. Navab, “Fast deformable registration of 3d-ultrasound data using a variational approach,” in Medical Image Computing and Computer Assisted Intervention (MICCAI), vol. 3216. Springer, 2006, pp. 915–923. [145] D. Zipes and J. Jalife, Cardiac electrophysiology: from cell to bedside, 4th ed. Philadel- phia, Pa.: W B SAUNDERS COMPANY , 2004. [146] D. P. Zipes and H. J. J. Wellens, “What have we learned about cardiac arrhythmias?” Circulation, vol. 102, no. 4, November 2000. 136

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Geometric methods of image registration and analysis

PDF

Techniques for efficient cloud modeling, simulation and rendering

PDF

Digital signal processing techniques for music structure analysis

PDF

Advanced coronary CT angiography image processing techniques

PDF

Advanced liquid simulation techniques for computer graphics applications

PDF

Image and video enhancement through motion based interpolation and nonlocal-means denoising techniques

PDF

Radio-frequency non-uniformity in cardiac magnetic resonance imaging

PDF

Advanced intra prediction techniques for image and video coding

PDF

Image registration with applications to multimodal small animal imaging

PDF

Signal processing methods for interaction analysis of functional brain imaging data

PDF

Wideband steady-state free precession for cardiac MRI

PDF

Advanced techniques for high fidelity video coding

PDF

Mocap data compression: algorithms and performance evaluation

PDF

Robust speaker clustering under variation in data characteristics

PDF

Timing analysis of coupled interconnect and CMOS logic cells in the presence of crosstalk noise

PDF

Visualizing and modeling vocal production dynamics

PDF

Texture processing for image/video coding and super-resolution applications

PDF

Efficient management techniques for large video collections

PDF

Distributed edge and contour line detection for environmental monitoring with wireless sensor networks

PDF

Resource allocation in OFDM/OFDMA cellular networks: protocol design and performance analysis

Asset Metadata

Creator Woo, Jonghye (author)

Core Title Variational techniques for cardiac image analysis: algorithms and applications

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 10/09/2009

Defense Date 08/27/2009

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag cardiac image analysis,image registration,image segmentation,OAI-PMH Harvest,variational method

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Kuo, C.-C. Jay (committee chair), Nayak, Krishna S. (committee member), Shung, Kirk K. (committee member), Slomka, Piotr (committee member)

Creator Email jonghyew@usc.edu,jschant@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m2654

Unique identifier UC1470181

Identifier etd-Woo-3316 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-268111 (legacy record id),usctheses-m2654 (legacy record id)

Legacy Identifier etd-Woo-3316.pdf

Dmrecord 268111

Document Type Dissertation

Rights Woo, Jonghye

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu