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