Global and Local Cardiac Functional Analysis with Cine MR Imaging: A Non-Rigid Image Registration Approach Except where reference is made to the work of others, the work described in this dissertation is my own or was done in collaboration with my advisory committee. This dissertation does not include proprietary or classified information. Wei Feng Certificate of Approval: Stanley J. Reeves Professor Electrical and Computer Engineering Thomas S. Denney Jr. Professor Electrical and Computer Engineering Jitendra K. Tugnait James B. Davis Professor Electrical and Computer Engineering George T. Flowers Dean Graduate School Global and Local Cardiac Functional Analysis with Cine MR Imaging: A Non-Rigid Image Registration Approach Wei Feng A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Auburn, Alabama August 10, 2009 Global and Local Cardiac Functional Analysis with Cine MR Imaging: A Non-Rigid Image Registration Approach Wei Feng Permission is granted to Auburn University to make copies of this dissertation at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iii Vita Wei Feng, only son and the third child of Zhiguo Feng and Shulan He, was born in Tangshan city, Hebei Province, China on Nov. 28, 1976. He graduated from the First High School of Luanxian in 1996 and was enrolled in Shanghai Jiao Tong University the same year, where he earned both his Bachelor of Engineering and Master of Science degrees in electrical engineering in 2000 and 2003, respectively. He then worked for Magima Inc. in 2003 for several months before he went to the United States with his beloved wife, Shuang Feng. He then entered the PhD program in electrical engineering of Auburn University in Aug. 2004, where he is expected to acquire the degree in Aug. 2009. He was accepted in the Master?s program in mathematics of Auburn University in Feb. 2008. He and his wife shared the same high school, the same university for their undergraduate and graduate years in China, and the same university in the United States in Auburn. They are blessed with a lovely daughter, Grace T. Feng, who was born on May 19, 2005. iv Dissertation Abstract Global and Local Cardiac Functional Analysis with Cine MR Imaging: A Non-Rigid Image Registration Approach Wei Feng Doctor of Philosophy, August 10, 2009 (M.S. (Mathematics), Auburn University, 2008) (M.E. (Electrical Engineering), Shanghai Jiao Tong University, 2003) (B.S. (Electrical Engineering), Shanghai Jiao Tong University, 2000) 143 Typed Pages Directed by Thomas S. Denney Jr. and Stanley J. Reeves According to the World Health Report 2003, cardiovascular disease (CVD) made up 29.2% of total global deaths, which highlights the importance of clinical cardiovascular research. Magnetic Resonance Imaging (MRI) has become an important technology to assist clinical diagnosis and treatment of cardiovascular disease that is non-invasive and radiation-free. Cine MRI can provide high-quality images of the beating heart with a good time resolution. Tagged MRI can be used to image the myocardium deformation by altering the magnetization spatially, which deforms with the myocardium during the cardiac cycle. The quantitative evaluation of cardiac functions can be divided into three categories: volumetric analysis, geometry analysis and deformation analysis. Both volumetric and geometry analysis requires the segmentation of the myocardium in the MR images. The myocardium segmentation identifies the shape and size of the ventricles that are used to compute ventricular volumes and derive geometric parameters, such as curvature. The deformation analysis measures the myocardium strains, which reflect the contractibility and stretchability of the myocardium on a local scale. v Manual myocardium segmentation can be extremely labor-intensive due to the large number of images that need to be processed in a limited amount of time. A large num- ber of fully automatic contouring algorithms have been proposed. But they are generally unreliable, and manual corrections are usually needed. In this dissertation, we propose a dual-contour propagation algorithm that propagates a small number of manual contours at two critical frames of the cardiac cycle to all the other time frames. Since manual contours are usually drawn at the two critical frames in most institutions, no extra labor is needed to perform the propagation. We validate our contour propagation algorithm on patient data, and it is shown to be statistically as accurate as manual contours. Although myocardium deformation analysis is usually performed with tagged MRI, there are several disadvantages. First, the tags in tagged MR images fade quickly and can only be reliable over about half of the cardiac cycle. Second, the spatial resolution of the tag lines are limited and sparse inside the myocardium. Furthermore, tagged MR imaging is a more complicated protocol and is not as commonly available as cine MR imaging. So it will be very beneficial for clinical purposes to be able to measure myocardium strains through cine MR imaging. In this dissertation, we propose a comprehensive algorithm with several regularization schemes to measure myocardium strain for both left and right ventricles. Both the contour propagation and myocardium strain analysis are based on non-rigid image registration. In this dissertation, we propose a comprehensive non-rigid image reg- istration algorithm that is adaptive, topology preserving and consistent. This algorithm is computationally more efficient due to its adaptive optimization scheme. In addition, the inverse consistency and topology preservation are achieved by regularization. vi Acknowledgments The author would like to thank Dr. Thomas S. Denney, Jr. and Dr. Stanley J. Reeves not only for their extensive teachings in the academic field, but also their guidance and kindness through the years. Without them, this work would never have existed. I also want to thank Dr. Jitendra K. Tugnait for his support and assistance in the course of my doctoral years. In addition, I want to thank Dr. Soo-Young Lee for introducing me to the community of Auburn University. The author want to thank Dr. Thomas H. Pate and Dr. Gary M. Sampson for their clear, concise and devoted teachings and guidance in the field of mathematics, which is a major source of enjoyment in my academic life. In addition, I want to thank Dr. Amnon J. Meir, Dr. Jo W. Heath, Dr. Narendra K. Govil, Dr. Stephen E. Stuckwisch and Dr. Jerzy Szulga for their teachings. The author also want to thank Dr. Louis J. Dell?Italia, Dr. Himanshu Gupta and Dr. Steven G. Lloyd for their kindness and accommodation for the last two years. I also want to thank Dr. Hosakote Nagaraj and Dr. Ravi Desai for their cooperation. Finally, I want to thank all my family for the endless love they gave me. Special thanks to my mom Shulan He and dad Zhiguo Feng for giving me a life and raising me through the hard times. I want to thank my two older sisters, Hui Feng and Bo Feng, who like my parents, loved me and supported me for my whole life. Lastly, I want to thank my wife, Shuang Feng, for sharing my adult life and sticking with me despite my numerous flaws, and most importantly, for giving birth to my beautiful and precious baby, Grace T. Feng. They would always be in my heart. vii Style manual or journal used Journal of Approximation Theory (together with the style known as ?aums?). Bibliograpy follows van Leunen?s A Handbook for Scholars. Computer software used The document preparation package TEX (specifically LATEX) together with the departmental style-file aums.sty. viii Table of Contents List of Figures xi 1 Introduction 1 1.1 Physiology of the Human Heart . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Cardiac Magnetic Resonance Imaging . . . . . . . . . . . . . . . . . . . . . 3 1.3 Overview of Cardiac Functional Analysis . . . . . . . . . . . . . . . . . . . . 5 1.4 Cardiac Functional Analysis: A Non-Rigid Registration Approach . . . . . 9 1.5 Adaptive and Topology Preserving Consistent Image Registration . . . . . . 11 1.6 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.7 Overview of the Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2 Adaptive and Topology-Preserving Consistent Image Registration 15 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 The B-spline Deformation Model . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Traditional Non-Rigid Image Registration . . . . . . . . . . . . . . . . . . . 22 2.3.1 The Cost Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3.2 Gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3.3 Hessian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3.4 Multi-resolution Optimization . . . . . . . . . . . . . . . . . . . . . . 25 2.3.5 Optimization Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.4 Topology-Preserving Consistent Image Registration . . . . . . . . . . . . . . 29 2.4.1 Consistent Image Registration with Topology Preservation . . . . . 30 2.4.2 Gradient and Hessian of Inverse Consistency Regularization . . . . . 33 2.4.3 Gradient and Hessian of Topology Preservation Regularization . . . 42 2.5 Adaptive Optimization Strategy . . . . . . . . . . . . . . . . . . . . . . . . 45 2.5.1 Selecting Adaptive Set of Control Points . . . . . . . . . . . . . . . . 45 2.5.2 Divide and Conquer for Large DOF . . . . . . . . . . . . . . . . . . 48 2.5.3 Dynamic Removal of Stagnant Control Points . . . . . . . . . . . . . 50 2.6 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.6.1 Inverse Consistency . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.6.2 Topology Preservation . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.6.3 Adaptive Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3 Dual Contour Propagation for Global Cardiac Volumetric Analysis 60 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2 Review of Existing Myocardium Segmentation Methods . . . . . . . . . . . 62 3.2.1 Automatic Segmentation Methods . . . . . . . . . . . . . . . . . . . 62 ix 3.2.2 Contour Propagation Methods . . . . . . . . . . . . . . . . . . . . . 67 3.3 Dual Contour Propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.3.1 Subjects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.3.2 Image Acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.3.3 Image Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.3.4 Contour Propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.3.5 Volumetric Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.3.6 Parameter Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . 72 3.3.7 Validation of Functional Parameters . . . . . . . . . . . . . . . . . . 72 3.3.8 Inter-User Variability . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.3.9 Comparison Between Single and Dual-Contour Propagation . . . . . 73 3.3.10 Comparison of PER and PFR Values in Normals and Hypertensives 73 3.3.11 Statistical Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.4.1 Parameter Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . 74 3.4.2 Validation of Functional Parameters . . . . . . . . . . . . . . . . . . 75 3.4.3 Inter-User Variability . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.4.4 Comparison Between Single- and Dual-Contour Propagation . . . . 78 3.4.5 Peak Ejection and Filling Rates in HTN . . . . . . . . . . . . . . . . 80 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4 Shape Regularized Myocardial Strain Analysis with Cine MRI 86 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.2 The Registration Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.2.1 Contour Regularization for Left Ventricular Registration . . . . . . . 90 4.2.2 Polar Regularization for Left Ventricular Registration . . . . . . . . 91 4.2.3 Custom Shape Regularization for Right Ventricular Registration . . 99 4.3 Strain Computation from Myocardial Deformation . . . . . . . . . . . . . . 101 4.4 Left and Right Ventricular Strain Results . . . . . . . . . . . . . . . . . . . 105 4.4.1 Contour Regularization for LV Strains . . . . . . . . . . . . . . . . . 105 4.4.2 Preliminary Study of LV Strains with Hypertensive Patients . . . . . 107 4.4.3 Polar Regularization for LV Strains . . . . . . . . . . . . . . . . . . . 108 4.4.4 RV Strain Measurement . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.4.5 RV Strain Results for Patient Population . . . . . . . . . . . . . . . 115 5 Conclusion 120 Bibliography 123 x List of Figures 1.1 Illustration of a cut section of a normal human heart . . . . . . . . . . . . . 2 1.2 Slice prescriptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 cine and tag MR images at ED and ES . . . . . . . . . . . . . . . . . . . . . 6 2.1 Traditional inverse consistency constraint illustration . . . . . . . . . . . . . 18 2.2 B-spline bases functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Embeded subspaces formed with multi-resolution B-spline bases . . . . . . . 26 2.4 CIR inverse consistency illustration . . . . . . . . . . . . . . . . . . . . . . . 31 2.5 Topology preservation penalizing function . . . . . . . . . . . . . . . . . . . 33 2.6 Flow chart for adaptive optimization algorithm. . . . . . . . . . . . . . . . . 46 2.7 Registration results with different inverse consistency parameters . . . . . . 52 2.8 Cost vs. iteration plots for different inverse consistency parameters . . . . . 52 2.9 LV registrations with different topology penalizing functions . . . . . . . . . 53 2.10 Cost vs. iteration plots for LV registration with different topology parameters 54 2.11 Brain registration with different topology penalizing functions . . . . . . . . 55 2.12 Cost vs. iteration plots for different topology penalizing functions: brain . . 56 2.13 Evolution of adaptive control points: forward direction . . . . . . . . . . . . 57 2.14 Evolution of adaptive control points: backward direction . . . . . . . . . . . 57 2.15 Cost vs. iteration plots for non-adaptive and adaptive brain registration . . 58 2.16 Time vs. iteration plots for non-adaptive and adaptive brain registration . . 59 3.1 Short axis cine MR images at ED and ES . . . . . . . . . . . . . . . . . . . 61 xi 3.2 Dual contour propagation scheme . . . . . . . . . . . . . . . . . . . . . . . . 70 3.3 Dual-contour propagated contours compared to manual contours . . . . . . 74 3.4 VTCs computed from manual and propagated contours . . . . . . . . . . . 76 3.5 Scatter and Bland-Altman plots of LV PER, ePFR and aPFR . . . . . . . . 77 3.6 Scatter and Bland-Altman plots of LV PER, ePFR and aPFR . . . . . . . . 79 3.7 VTC curves of LV for a normal volunteer . . . . . . . . . . . . . . . . . . . 80 3.8 VTC plots for a normal volunteer and a hypertensive patient . . . . . . . . 81 3.9 LV PER, ePFR and aPFR in normals and hypertensives . . . . . . . . . . . 81 4.1 ED and ES images of RV free wall . . . . . . . . . . . . . . . . . . . . . . . 89 4.2 Original image and contour image . . . . . . . . . . . . . . . . . . . . . . . 91 4.3 Polar decomposition of deformation . . . . . . . . . . . . . . . . . . . . . . . 93 4.4 Custom shape regularization: circle fitting . . . . . . . . . . . . . . . . . . . 100 4.5 Curvature illustration of custom shape regularization. . . . . . . . . . . . . 101 4.6 Normal directions of custom shape regularization . . . . . . . . . . . . . . . 102 4.7 Correlation between circumferential strains from cine and tagged MRI . . . 106 4.8 Mid-ventricular circumferential and radial strains . . . . . . . . . . . . . . . 108 4.9 Circumferential expansion rates in early and atrial diastole . . . . . . . . . 109 4.10 Circumferential strains with polar regularized registration . . . . . . . . . . 110 4.11 Estimated deformations without and with polar regularization . . . . . . . 111 4.12 Circumferential strains overlaid on images . . . . . . . . . . . . . . . . . . . 112 4.13 Registration results with different topology penalizing function . . . . . . . 114 4.14 Registration results with different topology weights . . . . . . . . . . . . . . 115 4.15 Registration results without and with custom shape regularization . . . . . 116 4.16 RV circumferential strain vs. time plot for a normal volunteer . . . . . . . . 118 4.17 RV circumferential strain vs. time plot for 34 studies . . . . . . . . . . . . . 119 xii Chapter 1 Introduction The American Heart Association?s 2009 update of heart disease and stroke statistics shows that coronary heart disease is the single leading cause of death in America for the year 2005, followed by stroke, which is still higher than lung cancer and breast cancer. This justifies the importance of clinical research in cardiovascular disease. Magnetic resonance imaging (MRI) is an important technological tool for assisting the clinical assessment of cardiovascular disease. MRI provides a non-invasive and radiation-free venue to visualize and measure the function of the beating heart. Through MRI, important cardiac func- tional parameters can be measured quantitatively to help both diagnosis and treatment of cardiovascular diseases. The human heart contains four chambers, left and right atrium and left and right ventricle (LV/RV). Both ventricles contract and relax about 72 times per minute throughout a person?s whole life. The left ventricle is responsible for pumping oxygenated blood from the lung to the whole body, while the right ventricle takes the deoxygenated blood from the whole body and delivers it through the pulmonary arteries to the lung to re-oxygentate. This leads to stronger myocardium for the left ventricle than for the right ventricle. Thus for a normal heart, the functional parameters of the left ventricle are the most important in evaluating the health of the heart. However, the functional parameters of the right ventricle can be of significant importance for assessing certain diseases, such as pulmonary hypertension, where the pulmonary artery connecting the right ventricle and the lung suffers high blood pressure due to dysfunction of the lung. In this case, the right ventricle will get bigger and stronger to compensate for the high pulmonary pressure. In this dissertation, our research focuses on functional analysis for both LV and RV, while putting a slightly heavier emphasis on LV. 1 1.1 Physiology of the Human Heart A whole period of contraction and relaxation of the LV is called a cardiac cycle, which can be divided into three major stages: systole, early diastole and atrial diastole. During diastole, the LV relaxes and blood flows from the left atrium into the LV. Before the end of diastole, the left atrium contracts and more blood is pushed into the LV (called atrial kick). Then, during systole, the LV contracts and pumps blood into the aorta and then to the whole body. A similar process happens to the RV, except that it pumps the blood through the right atrium to the pulmonary artery and eventually to the lung. This completes a cardiac cycle, which is repeated throughout the life span of a person. In this cycle, the contraction phase is called the systolic phase, and the relaxation phase is called the diastolic phase. Hence, the most relaxed phase and the most contracted phase of the ventricles are called end-diastole (ED) and end-systole (ES) phases. Figure 1.1 illustrates a cut section of a normal human heart. Although not noted on the image, the papillary muscles of both ventricles are easily seen in the figure. They are rooted in the myocardium wall and connect to the mitral valves for left ventricle and the tricuspid valves for right ventricle. Superior vena cava Pulmonary valve Right atrium Tricuspid valve Right ventricle Inferior vena cava Aorta Pulmonary artery Septum Left ventricle Mitral valve Aortic valve Left atrium Myocardium Figure 1.1: Illustration of a cut section of a normal human heart 2 As shown in Fig. 1.1, there are papillary muscles inside both left and right ventricles. They are connected with the chordae tendineae, which attach to the tricuspid valve in the right ventricle and the mitral valve in the left ventricle. In diastole, the ventricles relax and the papillary muscles contracts to pull the valves open, which allows blood to flow into the ventricles from the left and right atriums. In systole, the ventricles contract and the valves are closed by the blood pressure, which blocks the blood from flowing back to the atria and forces it to flow to the aorta (for the left ventricle) or the pulmonary arteries (the right ventricle). 1.2 Cardiac Magnetic Resonance Imaging Magnetic resonance imaging (MRI) [40] is a non-invasive imaging technique. It is free of ionizing radiation and able to generate good-quality images of the imaged tissue. Cardiac MRI provides information on morphology and function of the cardiovascular system [71]. For clinical cardiac applications, the imaging task has been more difficult compared to imaging of other parts of the body due to the constant movement of the human heart and breathing. As technologies advance, modern cardiac MRI is able to generate a series of good-quality images from multiple imaging planes for the deforming heart in a time that can be endured by many patients. These images provide a topographical field of view of the heart and its internal structures. As a result, cardiac MRI is widely used in hospitals and research institutions to assist clinical diagnosis and assessment of various cardiac functions. With cardiac MRI, images of the heart are acquired at different phases throughout the cardiac cycle, from systole to diastole, triggered by electrocardiogram (ECG). The number of phases imaged in one cycle varies, typically from 10 to 30. Generally, three perpendicular imaging planes, including one short-axis plane and two long axis planes are used to acquire three-dimensional data of the heart at each phase of the cardiac cycle. In the short-axis plane, images of multiple slices (typically 10 to 14) are acquired, covering the whole range from the base to the apex of the heart. The long axis planes include the two-chamber view (cross-section of the left ventricle and left atrium) and the four-chamber view (cross section 3 of all 4 chambers). Usually only one slice is imaged for each long axis plane. Sometimes more long axis planes are imaged. For example, one can acquire 12 long axis planes 15 degrees apart centered at the center line of the left ventricle. The two most commonly used cardiac MR imaging protocols are cine MR imaging and tagged MR imaging. Cine MR imaging generates images of the myocardium with a high soft-tissue contrast against the flowing blood. But it provides little contrast inside a tissue, such as the myocardium. Tagged MR imaging, on the other hand, has a very low soft-tissue contrast. However, tagged MR images provides contrast inside a tissue, which allows more accurate evaluation of tissue deformation during the cardiac cycle. To give a concrete concept of the analysis process, we will describe a typical set of images (called a study) generated by cine MRI. A study generated by tagged MRI will have the same structure, only with tagged images. To simplify the description, all the methods explained later will be illustrated based on this exemplar study. In reality, studies can be different in certain aspects, but the structure will be similar. In our prototypical cine MRI study, a patient is scanned to generate three groups of images, four-chamber view (4CH) group, two-chamber view (2CH) group and short-axis view (SA) group. The three groups are perpendicular to each other. An SA image is shown in Fig. (1.2) on the left, the 4CH and 2CH view planes are illustrated on top of it by dotted lines. The actual 4CH and 2CH view images and slice projections from each one?s perspective are shown in the middle and on the right. The SA group contains 13 parallel slices covering the LV and RV from base to apex. The SA slices are prescribed parallel to the LV/RV base plane. The 4CH view group contains one slice that cuts through the LV, RV, left atrium and right atrium. The 2CH view containing one slice is perpendicular to the 4CH view and only slices the LV and left atrium. Notice that the 4CH view and the 2CH view slices intersect at the center line of the LV. For every slice, a total of 20 time frames are acquired throughout the cardiac cycle beginning at ED. In total, there are 260 SA images, 20 2CH images and 20 4CH images, giving a final sum of 300 images. 4 Figure 1.2: Slice prescriptions from SA, 4CH and 2CH view perspectives. Left: SA view image with 4CH (horizontal) and 2CH (vertical) slice projections; middle: 4CH view image with SA (parallel lines) and 2CH (single line) slice projection; right: 2CH view image with SA and 4CH slice projections In systole, as the ventricles contract, the myocardium shortens longitudinally, twists and contracts circumferentially, and stretches radially. In diastole, the LV and RV relaxes and the reverse happens In the 4CH view, the basal points of both the LV lateral wall and the RV lateral wall are called the mitral annulus. Since the apex of the ventricles are relatively fixed during cardiac cycle, the mitral annulus provides a good indication of how much the ventricles shorten longitudinally at every time frame. Since the volume of the ventricles are a 3D measurement, one has to take account of the longitudinal shortening when computing them. Figure 1.3 shows both cine and tagged MR images from short-axis, 4-chamber and 2-chamber slices at ED and ES. 1.3 Overview of Cardiac Functional Analysis From a clinical perspective, the functional analysis of the heart can be roughly divided into three categories: global volumetric measurements of ventricles, ventricular curvature measurements and local myocardial strain measurements. The first category is the global volumetric measurements of the ventricles. Typical pa- rameters include ventricular volumes, ventricular mass and ventricular wall thickness. These parameters are usually measured at end-diastole and end-systole (LVEDV: left-ventricular end-diastole volume; LVESV: left-ventricular end-systole volume). Furthermore, one can 5 (a) (b) Figure 1.3: ED and ES images from short-axis, 4-chamber and 2-chamber slices of cine (a) and tagged (b) MR imaging. also deduce parameters such as ventricular ejection fraction (EF) from the above measure- ments. Specifically, the left-ventricular ejection fraction is defined as the ratio between the ejection volume (LVEDV-LVESV) and the end-diastole volume (LVEDV) expressed in percentage. For example, a 60% LVEF means that 60% of the left-ventricular end-diastole volume is ejected out of the LV. More advanced volumetric rate analysis, such as peak ejection rate and peak filling rate, can also be achieved by measuring ventricular volumes throughout the cardiac cycle. The second category is local ventricular curvature measurements. The curvature mea- surements can be used to deduce hypertrophy of the ventricles [56]. For example, a higher 6 curvature will indicate eccentric hypertrophy, in which the ventricle is enlarged like a bas- ketball. On the other hand, a lower curvature measurement will indicate concentric hy- pertrophy, in which case the ventricular muscle thickens to assume a conical shape like an American football. The last category is local myocardial strain measurements. The strain measurement at a point is basically a measurement of the degree of contraction or stretching of the my- ocardium at that point. For example, the circumferential strain of the LV is a measurement of the contraction or stretching of the LV myocardium in the circumferential direction. For a normal heart, the left ventricle will usually have a higher circumferential strain (more contraction) at the free wall (on the opposite side of the right ventricle) than at the septum wall (shared with the right ventricle). Since the ventricles shorten both circumferentially and longitudinally, measurements in both directions are important. For the left ventricle, circumferential shortening is approximately twice as important compared to the longitudinal strain. For the right ventricle, the importance of the two are reversed because the dominant deformation of the right ventricle is along the longitudinal direction. The peak strains are reached when the ventricles deform the most, which is around end-systole. Similar to the global volumetric rate analysis, one can obtain local strain rate analysis by measuring the strains throughout the cardiac cycle. Other important ventricular parameters include torsion and rotation. The different cardiac functional parameters described above are measured with different imaging modalities. All of them can be measured using cardiac MR imaging with both cine and tagged MRI. Naturally, each MR imaging protocol has advantages and disadvantages in measuring different parameters. The global volumetric parameters as well as the curvature parameters are best measured with cine MR imaging due to its high soft-tissue contrast. To achieve these measurements, one will need to segment the ventricles in the cardiac cycle. Through various kinds of segmentation techniques, one achieves a set of contours delineat- ing the ventricles. For this reason, the segmentation is also called myocardial contouring. Manual myocardial contouring can be prohibitive due to the large number of images to be 7 contoured. For the prototypical study, at least 80% of the 300 images need to be contoured if the volumetric rates analysis is desired. For this reason, many automatic myocardial contouring algorithms have been proposed in the last two decades. However, while these automatic myocardial contouring algorithms can be successful for a certain group of studies, they usually fail to generate accurate results for other studies. Thus, manual correction is required to guarantee the accuracy of measurements, which can still be labor intensive. On the other hand, semi-automatic myocardial contouring algorithms have also been explored. Although these methods require user input to initiate the segmentations, the user input can be limited to a small amount that is very much practical. The benefit is that they are more robust in generating accurate results and require much less post-correction compared to automatic contouring algorithms. The main cause of inaccuracies facing both automatic and semi-automatic myocardial contouring with cine MR imaging is the presence of the papillary muscles. This is especially severe for left-ventricular segmentation. In the left ventricle, the papillary muscles are much more prominent than in the right ventricle. They are usually separated from the left- ventricular endo surface at end-diastole. As the LV contracts, they touch the myocardium and show up merged with the myocardium in the cine MR images with no contrast between the them. The local myocardial strain analysis is traditionally measured with tagged MR images. The advantage of tagged MR imaging is that it provides contrast inside the myocardium by altering the magnetic properties of different areas inside the myocardium. This shows up in the tagged MR images as tag lines, as shown in Fig. 1.3. The tag lines deform with the myocardium, thus recording the deformations. The measured deformations can then be used to compute strains. However, there are some disadvantages of tagged MR imaging. First, since the tags are generated by altered magnetic fields, they fade as the heart deforms. Normally, the tag lines are reasonably well defined for 300 milliseconds. A resting normal human heart beats about 70 times per second, which is equivalent to about 850 milliseconds per cardiac cycle. Hence the tag lines are good for less than half of the cardiac 8 cycle. Since the tag lines are initiated at end-diastole, this means that one can get fairly good measurements through systole, and maybe into early diastole. For most of diastole, tagged MR imaging will not produce good results. The other problem facing tagged MR imaging is the limited resolution of the tag lines. Due tophysical limitations, the spacingbetween taglinescannotbetoosmall, usually about 10 millimeters, which is comparable to the thickness of the left-ventricular wall for a normal heart. Hence the tag lines inside the myocardium will be sparse. The sparse measurements are prone to imaging noise, which limits the accuracy of strain measurements from tagged MR imaging. This problem gets worse for the right ventricle since the right ventricular wall is much thinner (about 4?1 millimeters for a normal heart) than the left-ventricular wall. So the strain measurement for the right ventricular wall is essentially an interpolation of the deformation of the blood pool and the surrounding tissues. This results in inaccurate measurements of right-ventricular myocardial strain. 1.4 Cardiac Functional Analysis: A Non-Rigid Registration Approach The above described cardiac functional analysis tasks have been addressed with various kinds of algorithms. For example, myocardial segmentation has been tackled with tradi- tional image segmentation techniques, deformable models (active contours) and statistical shape models, to name a few. The strain analysis has been carried out in different ways. One way is a two-step approach, where the tag lines are identified first, then a deformation model is fitted to the identified tag lines to generate the estimated deformation. Another way is through harmonic phase (HARP) analysis in the Fourier frequency domain. In this research, we have formed a consistent approach of measuring almost all the above cardiac functional parameters using non-rigid image registration. The case for using non-rigid image registration for cardiac functional analysis comes from the physiology of the heart. Unlike most other parts of the human body, the heart is constantly deforming. Thus the fundamental question regarding cardiac function assessment is how the heart deforms. A perfectly healthy looking heart at the relaxed state (end-diastole) can show its 9 dysfunction (cardiac disease) when it does not deform as it is supposed to. By definition, image registration is a technique to recover the deformations of certain parts inside the images. So it is well suited for cardiac functional analysis. For the global volumetric analysis, the work is concentrated on myocardial contouring. While image registration does not segment myocardium by itself, it can be an effective tool to achieve fast myocardial segmentation. As stated above, automatic contouring algorithms lacks robustness and can be inaccurate. So we have drawn the conclusion that a semi- automatic contouring algorithms is the best way to solve this problem. Thus, we propose a dual contour propagation algorithm for myocardial segmentation. This algorithm requires the user to draw myocardial contours at end-diastole and end-systole frames of the MR images, which are usually drawn manually anyway in most clinical institutions. Non-rigid image registration is then performed between all the frames to estimate the myocardial deformation throughout the cardiac cycle. The ED and ES contours are then propagated using the estimated deformations to all the other frames. This generates two sets of contours for each frame outside of ED and ES. These two sets of contours are then combined using a weighting scheme to generate the final myocardial contours. The curvature analysis of the ventricles is then easily done with the propagated con- tours. For the local myocardial strain analysis, which is a deformation analysis by definition, non-rigid image registration is again a natural fit. As is documented in Section 1.3, the tagged MR imaging is limited in both scope and accuracy for determining both left and right ventricular myocardial strain throughout the cardiac cycle. We propose to measure myocardial strain with cine MR imaging. One of the benefits is that with cine MR imaging, one has high-quality images throughout the cardiac cycle. So diastolic functions can be as accurately measured as systolic functions, which is not the case with tagged MR imaging. The high-quality cine MR images also help with right ventricular deformation analysis, since the myocardium can be easily identified, as it has high soft-tissue contrast with both 10 the blood pool and the surrounding tissues. With tagged MR imaging, even identifying the location of the right ventricular myocardium is a difficult task. The disadvantage of cine MR imaging for deformation analysis is that it provides no contrast inside the myocardium. This renders it difficult to measure deformation of the mid-wall myocardium. However, both the endo (inside) and epi (outside) surfaces of the ventricles assume a certain shape, which will help guide the image registration to recover a reasonably accurate deformation even inside the myocardium. On the other hand, tissues such as the papillary muscles that are close to the endo surface will misguide the image registration to generate false deformations. This calls for the need of proper regularization to correct for the misguidance. We have developed different kinds of regularization schemes to improve the performance of cardiac image registration. These include contour regularization and polar regularization for left-ventricular registration and custom shape regularization for right-ventricular registration. 1.5 Adaptive and Topology Preserving Consistent Image Registration As the main analysis tool for this dissertation, non-rigid image registration has in itself been a major research field and generated numerous published papers [51, 75, 27]. Image registration can be viewed as a matching process that recovers the deformations necessary to deform parts of one image to match parts of the other image. By the above definition, the first question that needs to be resolved is the matching criterion, which we call the similarity measure. The similarity measure can be quite different for different applications. A common similarity measure is based on the image intensity differences for intra-model registration, where the image intensities for similar objects remain similar. Another prominent similarity measure is the mutual information measure, which measures the structural similarity between images. Mutual information is best suited for inter-model registration, where the intensity values of different images do not match, such as when registering an MR image and a Positron Emission Tomography (PET) image. 11 However, the most challenging part of designing an image registration algorithm is not only the formulation of the similarity measure. There are two other major issues that have to be resolved to achieve a well-defined registration algorithm. The first issue is what to use as a deformation model. As with any other numerical problems, the image registration problem has to be solved numerically. This means that the final solution has to be solvable in the first place. The necessary condition for a numerically solvable solution is to have a finite degree of freedom. However, a deformation is naturally defined as a continuum. Without any restrictions, it has an infinite degree of freedom and belongs to an infinite- dimensional space. To make it solvable, a subspace approximation has to be imposed. Common deformation models include Fourier series [45], B-spline functions [89], thin-plate spline functions [68], etc. B-spline functions are a popular choice for approximating smooth functions, since they have several advantages, including local support (which facilitates fast computation) and small approximation error [8]. The other issue that has to be resolved for any specific image registration applications is regularization. More often than not, a particular image registration application is an ill-posed problem, meaning that it has more than one possible solution. Adding regular- ization terms will condition the problem to be well-posed so that it has a unique solution. For example, in this research, for left-ventricular short-axis image registration, we have developed a polar regularization term that enforces smoothness in polar directions of the estimated deformation. Properly designed, a regularization term can improve the accuracy of the registration considerably. The smoothness regularization term is a generic term used by many algorithms. There are two other categories of regularization terms that are more deeply related to the reg- istration problem. The first is inverse consistency. Traditionally, the registration is only performed in one direction from one image to the other. This gives a single direction es- timated deformation that minimizes the given cost function. If one performs a similar registration in the other direction, in an ideal world, one should expect to get a deforma- tion that is approximately the inverse of the first estimated deformation. However, this is 12 usually not the case largely due to the numerical nature of the algorithm. Hence there is a disparity between the two deformations, which indicates that there is error from one or both registrations. This disparity is called inverse inconsistency. Consistent image registration algorithms are proposed to solve this problem. In consistent image registration, the similarity measures as well as the deformation models in both directions are included in the cost function. In addition, an inverse consis- tency term is added to the cost function that penalizes the inverse inconsistency between the two models. The other category of regularization term that is more involved with the registration problems is topology-preserving regularization. For most applications, especially medical image registration, the topology of the tissue remains constant during deformation. In plain words, thereisnofoldingortearingoftissuestructuresduringdeformation. Mathematically, this can be formulated as a topology preservation term that enforces the positivity of the Jacobian determinant of the deformation matrix. Once the similarity measure and regularization terms are chosen, the image registra- tion problem turns into an optimization problem, where an objective cost is formulated containing the similarity measure and the regularization terms. The goal is then to com- pute a solution that minimizes the cost. The optimization problem can be solved with all kinds of algorithms, such as a gradient descent algorithm, Newton?s algorithm, or even linear programming if one can manage to formulate the cost function into one. Different optimization algorithms will lead to different performances in terms of computational load, speed and accuracy. The common goal is then to minimize computational load, increase speed and improve registration accuracy. The computational load for registering two standard-size images can turn out to be prohibitively high. For example, when registering two two-dimensional images of 256?256 with a deformation model of the resolution 64?64, the solution space would have a degree of freedom of 642 ?2?2 = 16384 for a consistent image registration setup. With Newton?s 13 method, this would involve solving a linear system of size 16384?16384, which is impractical with common computation devices. In this research, we propose an adaptive and consistent image registration algorithm that also preserves topology. The adaptive algorithm curtails wasteful computation in the registration process and improves registration speed. 1.6 Notations The notation in this dissertation is chosen to be easily understood. Some specific rules are as follows. A bold small letter is used to denote a vector. For example, in 2D, a point written as x is equivalent to (x,y). The control parameters vector of a B-spline deformation is written as ? and is equivalent to (?1,?2,...,?n), wherenis the number of control points. A matrix is usually represented by a regular capital letter. Other specific notation is defined as it appears in the context. 1.7 Overview of the Dissertation This dissertation is structured as follows. In Chapter 2, we first describe in detail the formulation of image registration algorithms. Then we introduce the proposed fast, adaptive and topology preserving consistent image registration algorithm and show the experimental results. Then in Chapter 3, we describe the proposed dual contour propagation algorithm and demonstrate the clinical results that we achieved with the algorithm on patient MR imaging data. In Chapter 4, we describe the applications of the non-rigid registration algorithm in both left and right ventricular strain analysis and present clinical results. Finally, in Chapter 5, we conclude the research presented in this dissertation. 14 Chapter 2 Adaptive and Topology-Preserving Consistent Image Registration 2.1 Introduction Non-rigid image registration has been an important tool in medical image analysis ap- plications. In brain MR imaging, it is used to register patient data with a brain atlas to identify pathology. In cardiac MR imaging, myocardial segmentation can be achieved using contour propagation based on non-rigid image registration. Except the fundamental goal of matching the registered image data, there are three aspects of the non-rigid registration problem that have to be treated carefully to achieve clinically and practically sound registra- tion algorithms. The first aspect is computational efficiency, which is critical in applications where a large amount of image data needs to be processed in a relatively short time period. The second aspect is inverse consistency of the estimated deformations. Traditional image registration algorithms only try to estimate a deformation in one direction between two images, which may generate biased deformations. The third aspect is topology preserva- tion, which is valid in almost all medical image applications. A deformation can be inverse consistent yet changes topology of the underlying tissue structure. On the other hand, al- though a pair of topology-preserving deformations are both guaranteed invertible, they may not be inverse consistent. These aspects of the non-rigid registration problem have been tackled in a number of ways. For example, adaptive registration algorithms [79, 68] have been developed to improve computation speed. Consistent image registration algorithms [16, 15, 43, 26, 11, 78, 3, 5] have been proposed to generate inverse-consistent deforma- tions. Topology preservation constraints are sometimes incorporated in various registration algorithms [60, 34]. 15 Ever since its introduction [18, 101], mutual information (MI) showed great promise and has been used in various applications in hundreds of papers [75]. It is best suited for multi-modality image registration, such as images from CT and MRI, where intensity values are not consistent between images. For intra-model registration, however, it does not have any advantage over traditional intensity-based registration, such as those based on the sum of squared differences (SSD) measure. In fact, it is more complicated than SSD to implement due to its more sophisticated formulation. In this research, since our focus is on cardiac MR imaging, only SSD measure is considered. For traditional registration, we denote the two images to be registered as template image It and source image Is. Here the template image is the image whose samples are taken at fixed regular grids in the process of registration. The source image, on the other hand, will be sampled arbitrarily according to the estimated deformation. Generally the cost function can be written as, E(Is,It;T) = Esim +Ereg (2.1) where Esim and Ereg are the similarity measure and the regularization measure, and T is the mapping to be solved for. Using the SSD measure, the similarity measure can be written as [38] Esim = bardblIt ?T ?Isbardbl2 (2.2) where T represents the mapping from the template to the source. The similarity measure for traditional registration function (2.2) is ?inconsistent? in that it does not treat the two images exactly the same under a numerical setting. In the registration process, the template (It) is sampled on a fixed grid all the time, while the source (Is) is sampled arbitrarily in every iteration. A small yet expanding region in the template image will be mapped to a larger region in the source image, resulting in a smaller number of samples taken at the template image (thus carrying less weight). 16 On the other hand, a big yet shrinking region will have more samples in the registration process. This means that the algorithm is numerically biased toward shrinking regions. This inconsistency is philosophically unfounded since under many conditions, there is no solid reason for one to favor either image as the template (or the source). We call this bias as data inconsistency. Another problem with the traditional registration formulation is that the estimated deformation is uni-directional, from the template to the source. The data inconsistency means that if one were to register the two images in the reverse direction, the estimated deformation could be much different from that estimated in the forward direction. In other words, the two estimated deformations are supposed to be inverse of each other, but actually they are far from it. This is called inverse inconsistency. Futhermore, even with a inverse consistent registration, the estimated solution may be physically invalid given the physical property of the objects in the image. For medical image registration applications, the tissue structures are usually assumed to preserve topology. A change in topology, such as ?folding? or ?tearing? of a tissue structure should be avoided. We call this topology preservation. In [23, 22], the two images are registered in both directions. The data inconsistency is alleviated by combining the two mappings in both directions. However, it does not help with the inverse inconsistency since the two registrations are totally independent of each other. To improve both the data consistency and the inverse consistency, consistent image registration algorithms have been proposed [16, 15, 43, 26, 11, 78, 3, 5]. In [16, 15, 43, 26, 11, 78], the formulation of the cost function is essentially an extension of traditional registration. The similarity term can be expressed as follows, Esim = bardblIt ?T ?Isbardbl2 +bardblIs ?R?Itbardbl2. (2.3) As seen in (2.3), the data inconsistency problem is eliminated since there is no numerical bias toward either image in the registration process. To minimize the inverse inconsistency, 17 the following constraint term is added to the cost function: Eic = ?ic parenleftBigvextenddoublevextenddouble T ?R?1vextenddoublevextenddouble2 +vextenddoublevextenddoubleR?T?1vextenddoublevextenddouble2 parenrightBig , (2.4) where ?ic is the weighting parameter for inverse consistency. The first half of Eq. (2.4), vextenddoublevextenddoubleT ?R?1vextenddoublevextenddouble2, is illustrated in Fig. 2.1. Figure 2.1: Illustration of the inverse consistency constraint scheme used in traditional CIR algorithms. With the above two terms, one can write the basic cost function for consistent image registration as E(Is,It;T,R) = Esim +Eic +Ereg, (2.5) where again Ereg is an optional regularization term. Note that the cost function is the same if the template and source images are swapped. In [15, 16], the mappings T and R are modeled with a 3D Fourier series representation. The optimization problem is solved iteratively and alternatively. In each iteration, the mapping in one direction is fixed and the mapping in the other direction is updated. In each iteration, the inverse of a mapping needs to be approximated in order to proceed, as seen in (2.4). Inverting a 2D or 3D deformation is a difficult problem and generally can only be approached with numerical approximation [16, 11, 2, 1]. The inversion is computationally expensive, and approximation errors will occur. 18 Following the lead in [15], a method for inverting a mapping using Taylor series ap- proximation is used in [43]. Under this approximation, the inversion consistency can be implemented implicitly, thus reducing computation complexity. In [11], a symmetrization technique is proposed to avoid inversion. However, numerical instability can occur with the proposed method. A more recent development of the work in [15, 16] is shown in [26], in which the registration is extended to a manifold of images and transitive consistency is enforced on the solution. Other methods formulate the consistent registration problem as a continuous process and constrain the temporal smoothness of the deformation. Examples are the diffeomorphic registration proposed in [3, 5]. In this dissertation, we propose an adaptive, topology preserving consistent image reg- istration algorithm with a B-spline deformation model. With the B-spline model, the cost function is conveniently formulated in a different way to eliminate the need to invert a deformation. Furthermore, both the gradient and Hessian of the cost function can be eas- ily derived analytically. This allows us to implement fast-convergent algorithms such as Levenberg algorithm. It is numerically stable, and the implementation is fairly simple. Fur- thermore, a customizable topology preservation term is used to constrain physical validity of estimated deformations. Finally, the efficiency of the optimization process is improved by elliminating redundant computations through an adaptive optimization scheme. This chapter is organized as follows. In Section 2.2, the B-spline deformation model that is used throughout the dissertation is introduced. In Section 2.3, we describe in de- tail a generic formulation of the traditional non-rigid image registration algorithm. We will derive both the gradient and Hessian of the cost function mathematically. We also describe the optimization algorithm for this formulation. Then in Section 2.4, the proposed topology-preserving consistent image registration algorithm is described. The mathematical derivation of both the gradient and Hessian of the new regularization terms is presented. Section 2.5 presents the adaptive optimization strategy that we developed to improve com- putational efficiency. Then in Section 2.6, experiments are shown to illustrate the superior 19 performance of the proposed registration algorithm compared to the traditional non-rigid image registration algorithm. 2.2 The B-spline Deformation Model It has been shown that B-spline functions are a good fit for medical image analysis applications in that they are smooth functions and facilitate fast computation [92, 93]. Theoretically, there are infinitely many B-spline basis functions available. In practice, only a few are commonly used. These include the B-spline basis functions from degree 0 to degree 3. In the language of interpolation, a degree 0 B-spline corresponds to nearest neighbor interpolation and degree 3 B-spline corresponds to cubic interpolation. If we let ?n(x) be the n-th degree B-spline basis function, then it can be expressed in the convolution form as ?n(x) = ?0(x)??0(x)?...??0(x)bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipupright n+1 . (2.6) That is, a degree n B-spline basis is the result of (n+1) times convolution of the degree 0 B-spline basis. Degree 0 through degree 3 B-spline basis functions are shown in Fig. 2.2. 0 1 ?2 0 2 0 1 ?2 0 2 Figure 2.2: B-spline bases of degree 0 (top left), degree 1 (top right), degree 2 (bottom left) and degree 3 (bottom right) A B-spline function f(x) is a weighted sum of uniformly shifted and scaled B-spline basis functions: f(x) = summationdisplay k?K ck?n(xh ?k), (2.7) 20 where h is the spacing between adjacent B-spline bases (usually called knot spacing), ck is the weight parameter for the k-th B-spline basis, and K is the set of knot locations that define the support of the B-spline function f(x). A degreenB-spline functionf(x) is a smooth function and is continuously differentiable upton?1times. Inreality, B-splinefunctionsareusuallyconsideredinfinitelydifferentiable, since it may not be differentiable only at a finite number of locations. The smoothness of the B-spline functions is a desirable property especially for medical image analysis [93]. For example, in image interpolation, a smooth image is usually assumed. Also it has been shown that the B-spline functions generate smaller approximation errors when fitted to general smooth functions [8]. In this dissertation, B-spline functions are used for two purposes. The first is image interpolation, for which the cubic B-spline basis (degree 3) is used. The other application of B-spline functions is to construct the deformation model. B-spline basis functions are especially fitted for constructing deformation models for cardiac applications because the deformation of the heart is inevitably smooth. Considering both speed and accuracy, a quadratic B-spline basis (degree 2) is chosen as the basis for the deformation model. In 2D, a deformation can be written as m(x;?) = summationdisplay k,l?K?L ?k,l?2( xh x ?k)?2( yh y ?l), (2.8) where ?k,l is the (k,l)-th B-spline weight parameter, which we call a control point. K?L is the index set of all the control points and the cardinality of K?L is called the Degrees Of Freedom (DOF) of the deformation model. (xk,yl) is the location of the (k,l)-th control point. Here a uniform B-spline with constant spacing between control points is used. hx andhy are the scaling parameters inxandy directions that map (x,y) in the image plane to their corresponding locations in the deformation plane. When embedded in the registration formulated as an optimization problem, it is usually convenient to view the variables, which are the control points of the B-spline deformation, as a vector instead of a 2D array. Also 21 when there is no confusion, the degree of the B-spline basis is left out for the sake of brevity. So in the following mathematical derivations, we express the B-spline deformation model as m(x;?) = summationdisplay k?K ?k?( xh x ?xk)?( yh y ?yk). (2.9) Note that with the image size fixed, as the scaling parameters hx and hy decreases, the number of B-spline control points increases. This could pose a problem for computation when the images are large. For a standard 2D image size of 256?256, with a 4 : 1 pixel-to- control-point ratio, one would have a B-spline deformation model of roughly 64?64 = 4096, which is a fairly large number of variables. 2.3 Traditional Non-Rigid Image Registration 2.3.1 The Cost Function The traditional non-rigid image registration is formulated in a single direction from the template image It to the source image Is. The similarity measure using sum of squared differences is shown in Eq. (2.2) and can be expanded as follows: Esim = 1N bardblIt ?T ?Isbardbl2 = 1N summationdisplay k?K bracketleftbigI t(xk)?Is (xk +m(xk;?)) bracketrightbig2. (2.10) 2.3.2 Gradient Differentiating the cost function (2.10) with respect to control point ?i, we get 22 ?Esim(?) ??i = 2 N summationdisplay xk?K bracketleftbigI t(xk)?Is(xk +m(xk;?)) bracketrightbig? bracketleftBigg ??Is(tk)?t k vextendsinglevextendsingle vextendsinglevextendsingle tk=xk+m(xk;?) bracketrightBiggT ?m(xk;?) ??i (2.11) Let ek = It(xk)?Is(xk +m(xk;?)). For 2D image registration, bracketleftbigg?I s(tk) ?tk bracketrightbiggT = bracketleftbigg?I s(tk) ?tk,x , ?Is(tk) ?tk,y bracketrightbigg = [I?t,x(tk),I?t,y(tk)] and ?m(xk;?) ??i = bracketleftbigg?m x(xk;?) ??i , ?my(xk;?) ??i bracketrightbiggT =bracketleftbigm?x,i(xk),m?y,i(xk)bracketrightbigT , where tk,x = xk +mx(xk;?) and tk,y = yk +my(xk;?) are the mapped coordinates of the uniform sampling point (xk,yk). Plugging the above into (2.11), we get ?Esim(?) ??i = ? 2 N summationdisplay xk?K ek[I?s,x(tk),I?s,y(tk)] ? ?? m?x,i(xk) m?y,i(xk) ? ?? = ? 2NGTi EdF where Gi =bracketleftbigm?x,i(x1),m?y,i(x1),...,m?x,i(xM),m?y,i(xM)bracketrightbigT , (2.12) F =bracketleftbigI?t,x(t1),I?t,y(t1),...,I?t,x(tM),I?t,y(tM)bracketrightbigT , (2.13) 23 and Ed = Diag[e1,e1,e2,e2,...,eM,eM]T. (2.14) Finally, the gradient of the cost function (2.10) can be written in matrix form as ?Esim = bracketleftbigg?E sim ??1 ?Esim ??2 ... ?Esim ??C bracketrightbiggT = ?2GTEdF, where G = [G1,G2,...,GC]. (2.15) 2.3.3 Hessian Starting from (2.11), differentiating again with respect to ?j, we get ?2Esim(?) ??i??j = 2 N summationdisplay xk?K braceleftBiggbracketleftbigg ?m(xk;?) ??j bracketrightbiggT parenleftBigg?I s(tk) ?tk bracketleftbigg?I s(tk) ?tk bracketrightbiggTparenrightBigg?m(x k;?) ??i ? ek bracketleftbigg?m(x k;?) ??j bracketrightbiggT parenleftbig ?2Is(tk)parenrightbig?m(xk;?)?? i bracerightBigg = 2N summationdisplay xk?K braceleftBiggbracketleftbigg ?m(xk;?) ??j bracketrightbiggT parenleftBigg?I s(tk) ?tk bracketleftbigg?I s(tk) ?tk bracketrightbiggT ?ekparenleftbig?2Is(tk)parenrightbig parenrightBigg ? ?m(xk;?) ??i bracerightbigg In matrix form, we can write the Hessian matrix of the cost function as ?Esim = 2GBGT, (2.16) 24 where G was defined in (2.15) and B = ? ?? ?? ? B(t1) ... B(tM) ? ?? ?? ? with B(tk) = ?Is(tk)?t k bracketleftbigg?I s(tk) ?tk bracketrightbiggT ?ekparenleftbig?2Is(tk)parenrightbig. 2.3.4 Multi-resolution Optimization To improve convergence speed and avoid local optima, a multi-resolution approach is used. Suppose there are n layers of multi-resolution for the registration ordered from top to bottom, from coarse to fine. At each layer except the last one on the bottom, both It and Is are usually downsampled to a smaller size. The higher the layer, the bigger the downsampling ratio, and the smaller the image size. The deformation model is designed in a similar way. At each multi-resolution layer, the relative relationship between the image size and the deformation size is represented by the pixel-to-control-point ratio. The optimization starts at the highest layer with the lowest resolutions for both the images and the deformations. The deformations are optimized with the images at current resolution until convergence. Then the optimization advances to the next layers, where at least one of the following two requirements compared to the current resolution has to be met: a) the size of the images is increased; b) the size of the deformations is increased. This guarantees the improved accuracy of the registration in the next layer. In this dissertation, an exponential of 2 is used for the downsampling ratio of both the images and the deformations for each multi-resolution layer. For the deformations, this results in a series of nested subspaces from small to large [92]. Since the linear space formed by the B-spline basis in a coarse layer is a subspace of that formed by the basis in a fine layer, the estimated coarse B-spline deformation can be mapped exactly into the fine layer 25 with no approximation error at all. Since a 2D or 3D deformation is a vector function formed with scalar functions of single directional deformations (x, y and/or z directions), it suffices to illustrate the linear relationship between different resolution layers in single directional deformations. Mathematically, let mj and mj+1 be the corresponding single directional deformations (one of x and y directions for 2D deformation) at the j-th and (j +1)-th multi-resolution layer. Let Vj, j = 1,2,...,J be the linear spaces formed by the j-th layer B-spline basis. Then, V1 ?V2 ?...?VJ, as illustrated in Fig. 2.3. Let ... ... Figure 2.3: Illustration of embeded subspaces formed with multi-resolution B-spline bases. mj(x;?j) = summationdisplay k?Kj ?jk?( xhj x ?xjk)?( yhj y ?yjk) (2.17) be the representation of mj in Vj. For the sake of simplicity, we use the B-spline control point vector?j to representmj. Let the mapping ofmj intoVj+1 be represented by?j?(j+1). After optimization convergence in the (j+1)-th layer, represent the estimated deformation by ?j+1. Then the refinement for the estimated deformation achieved in the (j+1)-th layer is ?j+1? = ?j+1 ??j?(j+1). (2.18) 26 Using ?j+1? for j = 1,2,...,J?1, one can represent the estimated deformation at the J-th layer with the multi-resolution B-spline control point vector braceleftbig?1,?2 ?,...,? J ? bracerightbig. (2.19) Mapping the initial deformation ?1 and the following deformation refinements ?2?,...,?J?1? into the J-th layer, one gets ?J = ?1?J +?2?J? +...+?(J?1)?J? +?J?. (2.20) From Eq. (2.20), it is clear that the estimated final deformation after the J-layer optimiza- tion is a linear combination of deformation updates for all the multi-resolution layers. This is similar to the approach used in [79]. However, since a radial basis function was used instead of B-spline functions, the subspace nesting property is not valid for the algorithm in [79]. One can not map coarse resolution deformations or updates into the linear space formed by the basis functions in a finer resolution layer without incurring approximation error. Hence all the deformation updates have to be stored and re-evaluated during opti- mization. With a B-spline deformation model, there is no need to store the deformation updates in coarse layers since they dissolve seamlessly into a finer layer. 2.3.5 Optimization Algorithm Given the gradient and the Hessian, several optimization schemes are possible. The simplest one would be steepest descent (SD) algorithm using only the gradient. In each iteration, the SD algorithm searches for a local optimal solution that minimizes the cost function along the negative gradient direction. Once an optimal solution is found, the gradient is recomputed and the line search is repeated. The SD algorithm will terminate when the maximum number of iterations is reached or when the cost reduction is too small. Since it is a first order optimization algorithm, the convergence speed of SD is slow. 27 Toimprovetheconvergencespeed, especiallysincewehavebothgradientandHessianof the cost function, we can exploit the benefits of second-order optimization algorithms, which assumes a quadratic approximation of the cost function in each iteration. The prototypical second-order optimization algorithm is Newton?s method. Let ?i denote the estimated deformation for iteration i, then the estimate of the next deformation ?i+1 is ?i+1 = ?i ?(?E)?1?E. (2.21) However, the update equation 2.21 requires the Hessian of the cost function to be not only invertible, but also well conditioned for numerical stability. This requirement can be violated, in which case Newton?s method will fail numerically. To solve the numerical instability problem, various adaptations of Newton?s method have been proposed. A popular adaptation is the Levenberg-Marquardt algorithm with the following update equation: ?i+1 = ?i ?(?E +?i Diag?E)?1?E, (2.22) where?i is a tuning factor to scale the contribution of the diagonal of the Hessian. However, this still requires the matrix sum ?E +?i Diag?E to be well conditioned, which may not be the case for image registration. For the above reasons, the Levenberg algorithm is used for optimization of image registration in this dissertation. The Levenberg algorithm has the following update equation: ?i+1 = ?i ?(?E +?iI)?1?E, (2.23) where ?i controls the balance between linear approximation (SD) and quadratic approxi- mation (Newton?s method) in each iteration. When ?i is small, the algorithm approaches Newton?s method. On the other hand, when ?i is large, it approaches the SD algorithm. Notice that the numerical instability problem is solved with the Levenberg optimization 28 algorithm since the sum ?E +?iI is guaranteed well conditioned with a probability 1 as long as ?i is not too small. In implementation, the best way to compute the update in each iteration is not to follow the update equation strictly as to compute the inverse of the matrix ?E+?iI. The reason is that this matrix of size DOF2 can be very large, as we have mentioned previously. Inverting a large matrix is computationally prohibitive, which renders the algorithm impractical. A common practice is to use a linear system solver to solve a linear system, which is well developed and much faster. However, it is still impractical when the size of the linear system gets too large. At a certain size, the benefit of fast convergence speed of the second-order optimization will be compromised by the increased computation demand that comes from solving a large linear system. We solve this problem by observing and exploiting the banded structure of the Hessian matrix. This will be described in detail in Section 2.5.2. 2.4 Topology-Preserving Consistent Image Registration Let the template image It(x) and the source image Is(x) be defined on the sample set ? = {xi : i = 1...N.}. By defining both images It and Is on the same sample set ?, it is assumed that the boundaries of both images correspond with the underlying deformations. The goal of consistent image registration is to find a forward deformation mf and a backward deformation mb that minimize a cost function formed with It and Is. This section is structured as follows. In Section 2.4.1, we will describe the formulation of topology-preserving consistent image registration and mathematically derive both the gradient and Hessian of the cost function. This includes the gradient and Hessian formula for the inverse consistency term in Section 2.4.2 and the gradient and Hessian for the topology preservation regularization in Section 2.4.3. 29 2.4.1 Consistent Image Registration with Topology Preservation The registration problem is formulated as an optimization problem where the following cost function is minimized: E = Esim +Eic +Etopo. (2.24) Here, Esim is the image data matching term defined with SSD as Esim = 12N Nsummationdisplay i=1 bracketleftbiggparenleftBig It(xi)?Is(xi +mf(xi)) parenrightBig2parenleftBig Is(xi)?It(xi +mb(xi)) parenrightBig2bracketrightbigg . (2.25) The second and third terms are the inverse consistency constraint and the topology preser- vation constraint. We describe them in detail in the following sections. Inverse Consistency The term Eic measures the inverse consistency cost and is defined also with SSD as Eic = ?ic2N Nsummationdisplay i=1 bracketleftbiggparenleftBig mf(xi)+mb(xi +mf(xi)) parenrightBig2 + parenleftBig mb(xi)+mf(xi +mb(xi)) parenrightBig2bracketrightbigg . (2.26) In contrast to other consistent image registration algorithms [16, 15, 43, 26, 11, 78], the inverse consistency term defined in (2.26) requires no inverse of any deformation. It avoids the inversion by concatenating the forward and backward deformations to form a loop. This is illustrated in Fig. 2.4. ?ic is the weighting parameter for the inverse consistency cost. A bigger ?ic will enforce a stronger inverse consistency. Topology Preservation In medical images, the topology represents the anatomical tissue structure and should be preserved during registration in most cases. Mathematically, topology variation can be measured using the Jacobian determinant of the deformation. In 2D, the Jacobian 30 Figure 2.4: Illustration of the inverse consistency constraint formulated by concatenating the deformations. determinant of a deformation is defined as Jm(x) = det vextendsinglevextendsingle vextendsinglevextendsingle vextendsinglevextendsingle vextendsingle m?xx(x) m?xy(x) m?yx(x) m?yy(x) vextendsinglevextendsingle vextendsinglevextendsingle vextendsinglevextendsingle vextendsingle , (2.27) where m?xx(x) = dmx(x)dx and m?xy, m?yx and m?yy are defined similarly. Topology is pre- served when Jm(x) > 0, for all x ? ?. Furthermore, the value of Jm(x) is an index for the contraction or expansion of the structure at x. For example, Jm(x) > 1 represents tissue expansion and 0