4D Strain from Tagged Magnetic Resonance Images by Unwrapping the
Harmonic Phase Images
by
Bharath Ambale Venkatesh
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
May 9th 2011
Keywords: magnetic resonance imaging, tagged MRI, harmonic phase, phase unwrapping,
ventricles, strain, torsion, hysteresis
Copyright 2010 by Bharath Ambale Venkatesh
Approved by
Thomas S. Denney Jr., Chair, Professor of Electrical and Computer Engineering
Stanley J. Reeves, Professor of Electrical and Computer Engineering
Jitendra Tugnait, Professor of Electrical and Computer Engineering
Abstract
Accurate assessment of ventricular function is clinically important. Tagged magnetic
resonance imaging (MRI) is the modality of choice for calculating parameters such as torsion,
rotation and strain that characterize regional myocardial function. In this dissertation, a
method for reconstructing three-dimensional strain over time from unwrapped harmonic
phase measurements obtained from tagged MRI is presented. The procedure involved the
placement of branch cuts before unwrapping the phase images to remove inconsistencies
and obtaining displacement measurements. An affine prolate spheroidal B-spline (APSB)
model was then used to reconstruct strain from the estimated displacement measurements.
Initially, this method was implemented to obtain strain and torsion in the left ventricle (LV)
by the manual placement of branch cuts through a graphical user interface (GUI) before
unwrapping. This procedure is known as manual Strain from Unwrapped Phase (mSUP).
The strain obtained was compared to 3D strain obtained from a feature-based method and
2D strain obtained from harmonic phase strain measurements.
A method called the biventricular strain from unwrapped phase (BiSUP) for reconstruct-
ing three-dimensional plus time (3D+t) biventricular strain maps from unwrapped harmonic
phase (HARP) images was then introduced. This is an extension of mSUP to obtain 3D+t
strain in the right ventricle (RV). Accurate assessment of RV function is clinically important.
Compared to LV, however, analysis of RV function is relatively difficult because of a lack of
geometric symmetry and the comparatively thinner myocardium. Displacement estimates
were obtained from unwrapped harmonic phase measurements, in a procedure similar to
the mSUP, and were then used to reconstruct 3D strain using a discrete model free (DMF)
approach. The DMF method of reconstruction does not require the model to have geometric
symmetry and so performed efficiently to produce accurate strains. The BiSUP strain and
ii
displacements were compared to those estimated by a 3D feature-based (FB) technique and
a 2D+t HARP technique.
A computer-assisted Strain from Unwrapped Phase (caSUP) procedure to automatically
place branch cuts was then devised. A combination of simulated annealing and exhaustive
search methods were used to obtain the optimal branch cut configuration. The energy
function that was minimized in the search methods aimed to maintain spatial and temporal
continuity. This process reduced the manual intervention required to obtain 3D+t strain to
nearly one-third of that of the GUI-based approach. The strain obtained was compared to
3D strain obtained from mSUP.
To summarize, these works included validation and detailed comparisons of strains,
torsion and rotation parameters of both the LV and the RV in human subjects of different
pathologies and morphologies.
iii
Acknowledgments
I would like to take this opportunity to thank my advisor Dr. Thomas S. Denney Jr. I
feel priveleged to have had the chance to work with him for the past few years. He has been
both inspirational and immensely supportive to my work. He has given me ample freedom
to explore my thoughts in my research, given ideas/suggestions when necessary and guided
me through the various roadblocks I faced in my research. I will forever be grateful to him.
I would like to thank my advisory committee members, Dr. Stanley J. Reeves and Dr.
Jitendra Tugnait for providing valuable suggestions at various times. I would also like to
thank Dr. Amnon J. Meir for being the outside reader for my dissertation.
I would like to express my sincere gratitude to my labmates. The stimulating discussions
and seminars we had certainly sparked plenty of ideas in my research. They have been both
helpful and supportive, providing me with a great workplace atmosphere.
My friends and roommates have been great. They made moving to a new place, a
daunting task in my mind, fun and easy. I will never forget the laughs, the parties, and the
random sports and science discussions. They have been a part of some of the best memories
of my life. Thank you. I would like to thank the cricket team for a lot of great moments, I
really enjoyed the cameraderie and it was a privelege to be a part of the team.
Finally, I would like to thank my family for showering me with unconditional love and
affection. This would not have been possible without your constant encouragment and moral
support. I dedicate this dissertation to you.
iv
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Structure and Functioning of the Heart . . . . . . . . . . . . . . . . . . . . . 2
1.1.1 Structure of the Heart . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Functioning of the Heart . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.3 Cardiovascular Diseases . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Cardiac Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.1 MR Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.2 Cardiac Cine MR Imaging . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.3 Cardiac MR Tagging . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 Applications of Tagged MRI . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.1 Strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.2 Torsion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.3.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2 Tagged Image Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1 Feature-based Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.1 Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.1.2 Tag Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.1.3 3D Motion Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2 Harmonic Phase (HARP) MRI . . . . . . . . . . . . . . . . . . . . . . . . . . 31
v
2.2.1 HARP Refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.2 Extended HARP Tracking . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2.3 3D HARP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3 Optical-Flow Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.4 Non-rigid Registration Techniques . . . . . . . . . . . . . . . . . . . . . . . . 40
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3 Other MR Imaging Methods to Compute Myocardial Function . . . . . . . . . . 43
3.1 3D Myocardial Tagging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.2 Slice Following HARP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.3 zHARP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.4 Combined Harmonic Phase and Strain Encoded MRI . . . . . . . . . . . . . 46
3.5 DENSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5 Phase Unwrapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.2 Two-dimensional Phase Unwrapping Terminology . . . . . . . . . . . . . . . 57
5.2.1 Residues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.2 Dipoles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.3 Branch Cuts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.4 Quality Maps and Masks . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.3 Two-dimensional Phase-Unwrapping Algorithms . . . . . . . . . . . . . . . . 59
5.3.1 Path-Following Algorithms . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3.2 Least-squares Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6 3D Left Ventricular Strain from Unwrapped Phase: A GUI-based approach . . . 67
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
vi
6.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6.2.1 Human Subjects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6.2.2 Data Acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6.2.3 Unwrapping HARP Phase . . . . . . . . . . . . . . . . . . . . . . . . 69
6.2.4 Quality-Guided Phase Unwrapping . . . . . . . . . . . . . . . . . . . 71
6.2.5 Residues and Branch Cuts . . . . . . . . . . . . . . . . . . . . . . . . 72
6.2.6 Demodulated HARP Phase . . . . . . . . . . . . . . . . . . . . . . . 76
6.2.7 Inter-frame Phase Consistency . . . . . . . . . . . . . . . . . . . . . . 77
6.2.8 Motion and Strain Estimation . . . . . . . . . . . . . . . . . . . . . . 78
6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.3.1 Computation Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.3.2 Quantitative Comparisons . . . . . . . . . . . . . . . . . . . . . . . . 81
6.3.3 mSUP Strains in Normals and Patients . . . . . . . . . . . . . . . . . 84
6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7 Comparison Of 2D And 3D Torsion Measured From Tagged MRI . . . . . . . . 88
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
7.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
7.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
7.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
8 3D+t Bi-Ventricular Strain using Phase Unwrapped HARP . . . . . . . . . . . . 93
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
8.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
8.2.1 Human Subjects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
8.2.2 Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
8.2.3 Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
8.2.4 Displacement Estimation . . . . . . . . . . . . . . . . . . . . . . . . . 95
8.2.5 Strain Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
vii
8.2.6 Validation Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 98
8.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
8.3.1 Computation time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
8.3.2 Comparison of Displacement Estimates with FB Technique . . . . . . 100
8.3.3 Comparison of End-Systolic Strains with Feature-Based Technique . . 100
8.3.4 Comparison with HARP Strains . . . . . . . . . . . . . . . . . . . . . 101
8.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
9 Computer-Assisted Strain from Unwrapped Phase . . . . . . . . . . . . . . . . . 111
9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
9.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
9.2.1 Human Subjects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
9.2.2 Data Acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
9.2.3 Phase Unwrapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
9.2.4 Simulated Annealing . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
9.2.5 caSUP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
9.2.6 Motion and Strain Estimation . . . . . . . . . . . . . . . . . . . . . . 118
9.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
9.3.1 Computation Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
9.3.2 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
9.3.3 caSUP in Normals and Patients . . . . . . . . . . . . . . . . . . . . . 124
9.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
10 Torsion Hysteresis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
10.2 Magnetic Resonance Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . 127
10.3 Geometric analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
10.4 Torsional Deformation Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 128
10.5 Torsional Hysteresis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
viii
10.6 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
11 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
ix
List of Figures
1.1 The anterior view of the human heart. . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 The cardiac cycle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 LV Hypertrophy illustration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Slice prescriptions in cardiac MR imaging. . . . . . . . . . . . . . . . . . . . . . 10
1.5 Tag patterns illustration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6 Illustration of tag planes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.7 Illustration of displacement measurement using tags. . . . . . . . . . . . . . . . 12
1.8 Strain characterization in different coordinate systems. . . . . . . . . . . . . . . 14
1.9 AHA 17-segment model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.10 Depiction of basal and apical rotation in the LV. . . . . . . . . . . . . . . . . . 16
2.1 Contour propagation illustration. . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2 ML/MAP tag tracking results. . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3 Coupled B-snake tracker. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4 Deformable B-solid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.5 Thin-plate spline 3D deformation grid. . . . . . . . . . . . . . . . . . . . . . . . 26
x
2.6 Model tags. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.7 HARP filtering illustration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.8 Refined HARP tracking. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.9 Extended HARP tracking. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.10 3D HARP - material mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.11 3D HARP - In-plane motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1 SF HARP illustration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.2 SF HARP Planes Illustration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.3 DENSE: flood-fill unwrapping. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.4 DENSE: procedure illustration. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.1 1D phase unwrapping. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.2 The original MRI image and corresponding residues. . . . . . . . . . . . . . . . 58
5.3 The effect of branch cuts on unwrapping. . . . . . . . . . . . . . . . . . . . . . . 62
5.4 Unwrapping using least-squares methods. . . . . . . . . . . . . . . . . . . . . . 64
6.1 LV HARP illustration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.2 Placement of branch cuts in mSUP. . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.3 Phase unwrapping result from mSUP. . . . . . . . . . . . . . . . . . . . . . . . 74
6.4 Demodulated HARP and mSUP. . . . . . . . . . . . . . . . . . . . . . . . . . . 76
xi
6.5 Displacement vectors from mSUP. . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.6 Strain maps: mSUP vs. FB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
6.7 Torsion angle from mSUP in normals and patients. . . . . . . . . . . . . . . . . 85
7.1 3D vs. 2D rotation and torsion plots. . . . . . . . . . . . . . . . . . . . . . . . . 90
7.2 Peak mitral annulus displacement. . . . . . . . . . . . . . . . . . . . . . . . . . 90
8.1 PHTN heart illustration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
8.2 Segmentation of heart for BiSUP. . . . . . . . . . . . . . . . . . . . . . . . . . . 96
8.3 Branch cut placement for BiSUP. . . . . . . . . . . . . . . . . . . . . . . . . . . 97
8.4 Strain maps: BiSUP vs. FB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
8.5 Bland-Altman and correlation plots in LV: BiSUP vs. FB. . . . . . . . . . . . . 103
8.6 Bland-Altman and correlation plots in RV: BiSUP vs. FB. . . . . . . . . . . . . 104
8.7 Strain over time plots: BiSUP vs. HARP. . . . . . . . . . . . . . . . . . . . . . 106
8.8 Emin strain over time from BiSUP. . . . . . . . . . . . . . . . . . . . . . . . . . 108
9.1 Possible branch cut configurations given the no. of residues. . . . . . . . . . . . 114
9.2 Branch cut perturbation illustration. . . . . . . . . . . . . . . . . . . . . . . . . 117
9.3 Time taken for simulated annealing and exhaustive searches. . . . . . . . . . . . 118
9.4 Bland-Altman and correlation plots: mSUP vs caSUP. . . . . . . . . . . . . . . 121
9.5 Strain maps: mSUP vs. caSUP. . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
xii
9.6 Emin over time using caSUP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
10.1 Volume time curve for the 3 groups. . . . . . . . . . . . . . . . . . . . . . . . . 128
10.2 Torsion time curve for the 3 groups. . . . . . . . . . . . . . . . . . . . . . . . . 129
10.3 Schematic diagram for calculating torsional hysteresis. . . . . . . . . . . . . . . 130
10.4 Torsion Volume loops. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
10.5 Torsion hysteresis data in 3 groups. . . . . . . . . . . . . . . . . . . . . . . . . . 132
xiii
List of Tables
6.1 Displacement measurements: mSUP vs FB. . . . . . . . . . . . . . . . . . . . . 82
6.2 Strain comparison: mSUP vs. FB. . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.3 Strain comparison: mSUP vs. HARP. . . . . . . . . . . . . . . . . . . . . . . . 84
7.1 Comparison of 3D and 2D rotation and torsion. . . . . . . . . . . . . . . . . . . 91
8.1 Displacement measurements: BiSUP vs. FB. . . . . . . . . . . . . . . . . . . . . 100
8.2 Strain Comparison: BiSUP vs. FB. . . . . . . . . . . . . . . . . . . . . . . . . . 101
8.3 Strain comparison: BiSUP vs. HARP. . . . . . . . . . . . . . . . . . . . . . . . 105
9.1 Strain comparison: mSUP vs caSUP. . . . . . . . . . . . . . . . . . . . . . . . . 124
xiv
Chapter 1
Introduction
Image processing has found a number of applications in medical imaging. It has been an
essential tool for diagnostic purposes. It has been observed that geometric and mechanical
properties of the heart are key parameters in the prediction of the condition of the heart
and the efficacy of therapies. Magnetic resonance imaging (MRI) in particular has a lot of
advantages over other modalities. First, it uses non-ionizing radiation, unlike CT or PET,
and therefore has no known harmful effect on the body. Second, it has good image quality
and provides for 3D imaging. Cine MRI captures slices of the heart at various angles that
can then be put together to form a 3D geometric model. It is capable of providing a series of
such images over the entire cycle of the heart. Most importantly, tagged MR imaging helps
in calculating key mechanical properties of the heart. In addition, non-invasive analysis of
cardiac motion is possible through MR imaging. It is, however, necessary that the image
analysis be both efficient and accurate for the purpose of clinical use. In this dissertation,
we are concerned with tracking cardiac tissue and modeling the deformation of the cardiac
LV from tagged MR images.
This chapter starts off by introducing the structure and functioning of the heart. It is
followed by an overview of cine MR imaging and tagged MR imaging. Chapter 2 reviews the
current techniques in analysis of tagged MR images and in particular the harmonic phase
method (HARP). Chapter 3 gives an overview of the latest techniques being developed for 3D
imaging. Chapter 4 lists the objectives and contributions of the current research. Chapter
5 lists some of the common phase unwrapping algorithms that were implemented to test
their performance on phase images from tagged MRI. Chapter 6 explains the formulation of
a procedure to obtain 3D strain from tagged MR images by using a graphical user interface
1
(GUI) approach to the phase unwrapping problem (mSUP). Chapter 7 compares 2D (HARP
strains) and 3D (mSUP) strains and explores the advantages of the obtained 3D strains.
Chapter 8 extends the mSUP method to obtain strain in the right ventricle (BiSUP). Chapter
9 details a method (caSUP) to reduce the manual intervention involved in the mSUP method.
In chapter 10, an explanation of the energy mechanics of the heart using commonly measured
parameters is explored. Chapter 11 lists suggestions for future work.
1.1 Structure and Functioning of the Heart
1.1.1 Structure of the Heart
The heart is a hollow muscular organ situated between the lungs and the ribcage, and
is surrounded by a fibrous sac called the pericardium [1], [2]. It is about the size of a fist
and is located underneath the sternum.
The human heart consists of four chambers (Figure 1.1): two atria and two ventricles.
The right side of the heart (the right atrium and the right ventricle) is in charge of circulation
of blood through the lungs or pulmonary circulation. The left side (the left atrium and the
left ventricle) is in charge of pumping blood to the rest of the body or systemic circulation.
There are two sets of vessels connected to the heart: arteries and veins. There are two
arteries connected to the heart: the aorta that takes oxygenated blood away from the left
ventricle; and the pulmonary artery that takes deoxygenated blood from the right ventricle
to the lungs. There are two sets of veins connected to the heart: the superior and inferior
vena cavae return deoxygenated blood to the right atrium; and the pulmonary veins return
oxygenated blood from the lungs to the left atrium.
The septum divides the right and left sides of the heart, preventing passage of blood from
one side to the other. The heart has four valves: pulmonary, aortic, mitral and tricuspid,
that allow for blood flow within the heart and throughout the body. Cusps are segments
in the valve that regulate the blood flow. Valves are either bicuspid or tricuspid. The
pulmonary valve guards the opening of the right ventricle to the pulmonary artery. The
2
aortic valve helps blood flow to the aorta from the left ventricle and prevents it from flowing
back. The mitral and tricuspid valves regulate blood flow from the atria to the respective
ventricles. They are connected to the ventricular walls through fibrous strands called chordae
tendinae.The papillary muscles limit the movement of the two atrio-ventricular valves and
prevent them from inverting.
The wall of the heart consists of three layers. The myocardium is the thick muscular
layer of the heart responsible for the contraction and relaxation of the heart. To the outside
of the myocardium is the epicardium, which consists mostly of tissue and fat. The inner
lining of the heart is made of smooth membrane and is called the endocardium.
The cells that comprise cardiac muscle are called cardiomyocytes [3]. The cardiac my-
ocyte is a specialized muscle cell that is approximately 25? in diameter and about 100?
in length. The myocyte is composed of bundles of myofibrils that contain myofilaments.
The myofibrils have distinct, repeating microanatomical units, termed sarcomeres, which
represent the basic contractile units of the myocyte.
1.1.2 Functioning of the Heart
The heart contracts rhythmically to move blood throughout the body through the circu-
latory system. Each beat of the heart encompasses a cardiac cycle. A cardiac cycle consists
of three phases: the atrial systole, the ventricular systole and the complete cardiac diastole.
Figure 1.2 explains the complete cardiac cycle. Atrial systole involves contraction of atrial
walls to pump blood to the ventricles. Ventricular systole involves contraction of ventricu-
lar walls to pump blood out to the aorta and the pulmonary arteries. The atrioventricular
valves close, preventing blood flowing back to the atria, while the pulmonary and aortic
valves open. The beating sound of the heart is from the opening and closing of these valves.
The complete cardiac diastole involves relaxation of the entire heart. Blood fills up from the
vena cavae and the pulmonary veins into the right and left atria respectively.
3
Figure 1.1: The anterior view of the human heart [4].
4
Figure 1.2: The cardiac cycle [2]. The arrows indicate the direction of blood flow.
The systematic and rhythmic beating of the heart is controlled by electrical impulses.
These electrical impulses are provided by two electrical nodes: the sinoatrial node and the
atrioventricular node. During diastole, there are no electrical signals traveling through the
heart. The sinoatrial node sends out an electrical signal, which starts the contraction of
atria leading to atrial systole. At the end of the atrial systole these signals reach the atri-
oventricular valve, which then triggers ventricular systole.
1.1.3 Cardiovascular Diseases
There are many diseases that affect the heart, but we shall now take a look at a few we
presently study [5].
5
Ischemia is a restriction in blood supply, generally due to factors in the blood vessels,
with resultant damage or dysfunction of tissue. Myocardial infarction is a disease state
that occurs when blood supply to a part of the heart is interrupted or severely reduced.
This happens when one of the coronary arteries is blocked. The blockage is caused by the
building up of plaque inside the arteries; this is called atherosclerosis. When the plaque
deposits disintegrate at once, a blood clot forms and blocks the flow of blood. This leads to
a heart attack. By looking at strain maps, the area of infarction can be identified, and the
progress observed.
Mitral regurgitation, also known as mitral insufficiency, is the abnormal leaking of blood
through the mitral valve, from left ventricle into the left atrium. The mitral valve is com-
posed of valve leaflets, the mitral valve annulus (which forms a ring around the valve leaflets),
papillary muscles and chordae tendinae. A dysfunction of any of these can cause mitral re-
gurgitation. Long-standing mitral regurgitation may show evidence of left atrial enlargement
and left ventricular hypertrophy. Different phases of the disease show different symptoms.
Hypertension or High Blood Pressure is defined in an adult as a systolic pressure in
excess of 140 mm Hg and/or a diastolic pressure in excess of 90 mm Hg. The heart has to
work harder to maintain blood pressure, and thus, would be more prone to injuries. This
leads to ventricular remodeling and other heart diseases. By studying the increase in volume
of the heart, progress of the hypertensive heart can be monitored.
Ventricular remodeling refers to changes in size, shape and function of the heart after
injury to the left ventricle. It might be because of any number of causes, acute myocardial
infarction, chronic hypertension, and valvular heart diseases included. Left ventricular hy-
pertrophy is the thickening of the myocardium of the left ventricle. There are two types of
hypertrophy. Concentric hypertrophy is caused by increased systolic loads (pressure over-
load), leading to an increase in the mass of the chamber. Sarcomeres are added in parallel
and the myofibers get fatter. Eccentric hypertrophy is caused by increased diastolic volume
(volume overload), leading to an increase in myocardial mass. Sarcomeres increase in series
6
Figure 1.3: Illustration of the change in size and shape of the left ventricular wall for con-
centric and eccentric hypertrophy compared to normal heart [2].
and the myofibers get longer. The thickness of the affected chamber walls may be the same
as in a healthy heart or thinner. An illustration of hypertrophy is shown in Figure 1.3.
Diabetes is a disease where the body does not produce enough insulin or when the
insulin is not utilized properly (insulin resistance). Insulin resistance is associated with
fatty buildups in arteries (atherosclerosis) and blood vessel disease, even before diabetes is
diagnosed. Therefore, diabetes raises the risk of heart attacks and strokes.
1.2 Cardiac Imaging
Many imaging modalities like Computed Tomography (CT), Magnetic Resonance Imag-
ing (MRI), X-rays, Echocardiography (ultrasound) can be used to image the heart, but we
shall limit our discussion to cardiac MR imaging. This report concentrates on images ob-
tained by the method described in Prince and Links [6]. As mentioned earlier, MRI provides
many advantages over other imaging modalities and is especially useful for some applica-
tions. It is non-ionizing, non-invasive, allows for 3-D imaging and tagged MR imaging has
the ability to measure tissue motion. We shall take a look at the basis for MR imaging and
how cardiac MR imaging is performed in this section.
1.2.1 MR Imaging
Standard MRI is based on exciting the protons in hydrogen atoms. The tissues in
the human body consist mostly of molecules of water and fat, which are the source of
7
hydrogen atoms. The protons possess a spin, and when a magnetic field B0 is introduced,
these spins align themselves with the main field. These spin at a frequency known as the
Larmor frequency, ? = ?B0, which is proportional to the magnetic field applied. ? is the
gyromagnetic ratio, which is 42.58 MHz/T for protons.
This magnetization is then excited by the application of a radio-frequency (RF) pulse
such that it introduces a second magnetic field B1 perpendicular to the main field. This
field introduces a transverse magnetization vector, Mxy = M0 sin(?), orthogonal to the main
field and a vector, Mz = M0 cos(?), parallel to it. M0 is the initial equilibrium longitudinal
magnetization, and ? = null ?B1dt is the flip angle. The proton precesses around the main
field at the Larmor frequency and the rotating transverse signal produces a current in an
RF coil placed near the signal. This is the signal measured by an MR scanner.
After the excitation, Mz recovers exponentially to M0 with a rate constant of T1, and
Mxy decays exponentially with a rate constant of T2. T1 and T2 are a result of spin-lattice
and spin-spin interactions respectively. These parameters depend on the tissue and provide
a source of contrast for imaging. If each region experiences a unique magnetic field, then we
can image the position of the region. This is accomplished by spatially varying the magnetic
field.
The second step of spatial localization is called phase encoding. A magnetic gradient
field is applied briefly in one of the directions (y). As the change in frequency is very brief,
when the gradient is switched off, it causes a change in phase that is proportional to the
distance. The protons in the x direction have similar phase. The protons in the y direction
have different phases. In a standard spin echo, the number of phase-encoding steps is equal
to the number of lines in the matrix. Each step is performed with an incremental change in
the strength of the phase-encoding gradient. The last step of spatial localization is frequency
encoding, applied in the direction perpendicular to the phase-encoding direction (x). The
precessional frequency of columns of nuclear spins varies in the direction of the frequency-
encoding gradient. This gradient is applied during data acquisition.
8
The MR signal is a mix of all these frequencies (encoding in the frequency-encoding di-
rection) and phase shifts (encoding in the phase-encoding direction). The Fourier Transform
(FT) of the obtained signal gives describes the frequency spectrum. This is sampled for each
row of the slice to obtain the entire Fourier space for the image. An inverse FT then yields
the image. Back-projection techniques can be used when polar schemes are used to obtain
the signals.
1.2.2 Cardiac Cine MR Imaging
The technique in Section 1.2.1 describes how to obtain cardiac images known as cine
images. Generally, the images are obtained for an entire cardiac cycle. Images are obtained
as a series of short-axis (SA) images that cut across the length of the heart, from the base to
the apex. The left ventricle appears circular in these images. A long-axis (LA) slice would
be orthogonal to the SA slices. It can be prescribed either as a set of radial slices, or as
2-chambered and 4-chambered LA images. Figure 1.4 illustrates how slices are prescribed.
Cine images are used to measure parameters such as wall thickening, stroke volume, ejection
fraction and wall thickness. In this research, we are concerned with the measurement of
strain, for which tagged MR images are used. The next section discusses how MR tagging
is done.
1.2.3 Cardiac MR Tagging
One of the major advantages of using MRI for cardiac imaging is that it can measure
heart motion [8?10]. There are many ways for estimating heart motion using MRI, like phase
contrast MRI (PCMRI), tagged MRI (TMRI) and displacement encoding with stimulated
echoes (DENSE). In this dissertation, we are concerned only about extracting information
from tagged MR images.
Tags or temporary fiducial markers that change with the underlying tissue can be placed
in the image. The tag patterns can be of many types, but the two most popular ones are the
9
Figure 1.4: Slice prescriptions in cardiac MR imaging. Right: The slice prescriptions. The
bold lines indicate the plane of the slices. Left: The corresponding cardiac MR images [7].
10
(a) (b)
Figure 1.5: A midventricular slice for two different tag patterns: parallel taglines (a) and
grid tags (b).
parallel stripe pattern (Figure 1.5(a) ), and the grid pattern with perpendicular tag planes
(Figure 1.5(b)).
These tags are generated by a procedure known as SPAtial Modulation of Magnetization
(SPAMM) tagging procedure [11]. Tags are created by dephasing the magnetizing gradient
by the application of a spoiler gradient perpendicular to the tags. The tags are introduced
at the initial time frame, usually at end-diastole, and are tracked for an entire cycle. The
tags fade over time due to T1 relaxation. The tags cut across as planes across the entire
cross-section of the slice, and the planes deform with time as shown in Figure 1.6.
The tags are encoded onto the tissue and can be considered as material points. They
are a measure of the displacement of the underlying tissue [12]. Tags, though, can measure
displacement only in the direction perpendicular to them as shown in Figure 1.7. To obtain
a 3D displacement field, tags in three orthogonal directions are required [13]. In the next
chapter, a discussion of the present methods for extracting information from tagged images
is presented.
11
Figure 1.6: Illustration of tag planes [7]. Left: Tag planes introduced at end-diastole. Right:
Deformation of the planes observed at end-systole.
Figure 1.7: Illustration of displacement measurement using tags [14]. As observed in the
figure, the point r originates from point p but the only displacement in the direction per-
pendicular to the tag plane can be measured.
12
1.3 Applications of Tagged MRI
Tagged MRI facilitates the capture of regional function of the heart [15,16]. Measures of
function can be obtained transmurally, at several levels between the base and the apex, and
along the different segments of the heart wall. These measures of function include strain,
strain rate, shear, rotation and torsion.
1.3.1 Strain
Strain defines the amount of stretch or compression along a material line or elements of
fibers, i.e. normal strain, and the amount of distortion associated with the sliding of plane
layers over each other, i.e. shear strain [17]. Strain, in general depends on the particular
direction being considered. During systole, typically, the heart wall may be thickening in
the radial direction, leading to increased radial strain (Err), and simultaneously shortening
in the circumferential and longitudinal directions, leading to decreased strains (Ecc and Ell
respectively). Also, the angle between two material directions in the wall will generally
change due to shear strains (Ecl,Erl,Ecl). Figure 1.8 shows the three normal strains and
three shear strains usually measured for a given 3D displacement field [18]. Strain can also
be characterized in a direction based on the orientation of fibers in the heart. Figure 1.8
illustrates both the above characterizations. Alternatively, strain can also be characterized
in a coordinated system-independent way through the corresponding eigenvectors and eigen-
values. This way of characterization gives maximal stretching (Emax) and shortening (Emin)
strains. The time derivative of the strain gives the strain rate.
The LV is segmented to represent the obtained regional strains in an arrangement that
makes it easy to compare the results in different conditions or between patients. Figure 1.9
illustrate the segmentation in more detail.
13
Figure 1.8: Strain characterization in two different coordinate systems [18].
1.3.2 Torsion
Torsion measures the twisting of an object due to an applied torque. The LV undergoes
a wringing motion to pump blood during systole, and subsequently, has torsion associated
with it. Figure 1.10 illustrates the measurement of torsion in LV. The apex and the base
both undergo rotation with respect to their position at end-diastole; the difference between
these rotational components gives the torsion. Twist is the ratio of the measured torsion
and the distance between the apex and base.
Torsion has been traditionally computed in a 2D sense. A basal and an apical slice
are captured using tagged MRI. The points are then tracked through the time frames using
either one of the mesh tracking methods [20?22] or the HARP method [23]. The rotation
is calculated in each of the two slices through the cycle. The difference in the rotational
components between the apex and the base gives the desired torsion measurement.
2D torsion measurements can be erroneous because there usually is a longitudinal short-
ening of the LV which results in a change in the position of the basal slice through the cycle.
As a result, 3D methods to compute the true torsion have been proposed [24?27].
14
Figure 1.9: The muscle and cavity of the left ventricle can be divided into a variable number
of segments [19]. Based on autopsy data the AHA recommends a division into 17 segments
for the regional analysis of left ventricular function or myocardial perfusion. The left ventricle
is divided into equal thirds perpendicular to the long axis of the heart. This generates three
circular sections of the left ventricle: basal, mid-cavity, and apical. Only slices containing
myocardium in all 360o are included. The basal part is divided into six segments of 60o each.
Similarly the mid-cavity part is divided into six 60osegments. Only four segments of 90o
each are used for the apex because of the myocardial tapering. The apical cap represents
the true muscle at the extreme tip of the ventricle where there is no longer cavity present.
15
Figure 1.10: Schematic depiction of basal and apical rotation during systole in the normal
LV [23]. Arrowheads indicate the direction and magnitude of rotation.
1.3.3 Applications
The regional measures of function help us better understand the normal functioning of
the heart [15,28?30]. The normal motion of the heart is complex and encompasses regional
non-uniformities. The differences in the measures of function are not only time-dependent,
but also transmural. They can also be used as a means to characterize different cardiovas-
cular diseases such as valve disease [31], hypertension [32], myocardial infarction [33] and
cardiomyopathy [34?37].
It has been observed that in patients with mitral regurgitation, the rate of untwisting
in early-diastole reduces when compared to normal volunteers [31,38]. Patients with hyper-
tension tend to have hearts that experience higher amount of torsion [39] and nonuniform
circumferential shortening [32] to compensate for the increased blood pressure. Infarcted ar-
eas of the myocardium tend to experience lower strains than normal myocardial tissue [33].
The effect of age, sex and size on the functioning of the heart has also been studied us-
ing tagged MRI [18,40]. As a clinical tool, tagged MRI provides a detailed functional and
quantitative assessment of the function of the heart.
16
Chapter 2
Tagged Image Analysis
Tagged MR images provide the ability to track tissue in the heart. Tracking the heart
tissue helps us derive regional parameters such as strain, and torsion. Tagged image analysis
is the process of extraction of this information from tagged MR images.
Tag tracking involves extracting displacement estimates from the tag lines in the image.
Various methods have been proposed for this procedure. We shall take a look at a few of
the most widely used methods ? optical flow-based methods [41?43], feature-based data
extraction methods [14,44?47] and the harmonic phase (HARP) based methods [27,48?52].
This procedure is complicated because of the complex and nonuniform motion of the
heart, the through-plane motion (i.e. the movement of the base of the heart towards the apex
as the heart contracts), the fading of the tag lines through the cardiac cycle, and various
artifacts found in images due to incorrect imaging. It is important that the process be quick
and accurate and involve the least amount of manual intervention, in spite of the difficulties,
to be clinically viable. In this chapter, each procedure is further discussed.
2.1 Feature-based Methods
Tagged image analysis using feature-based methods involves three steps:
1. Segmentation of the myocardium.
2. Tag tracking.
3. 3D motion reconstruction.
17
2.1.1 Segmentation
Delineating the endocardial and epicardial boundaries of the heart is important for a
number of reasons. Firstly, from the myocardial boundary contours, global functional pa-
rameters such as ejection fraction, volume, mass, etc, can be computed. Secondly, regional
myocardial functional parameters are computed in segments within the myocardial bound-
aries. The 17-segment AHA model is normally used for this purpose. Thirdly, data outside
the myocardium is generally noisy and may interfere in tracking of points within the heart
wall in most techniques. Finally, the contour data helps in registration of the SA and LA
images, and eventually, in the development of the 3D deformation model.
Automatic methods for myocardial segmentation in tagged images have to overcome
problems of high noise, intensity inhomogeneities and complex backgrounds. Morphological
image processing [44,53] was first used for automated tagged MR myocardial segmentation.
It involves a morphological closing operation to remove the tag lines found in the tissue,
before using intensity-based thresholding. Markov random fields [45,54] have been used
too. These methods involve the modeling of the background tissue, tagged tissue and the
myocardial tissue as random variables and then using a MAP decision rule to identify the
myocardial boundaries. Gabor filter [55] based methods involve the removal of the tag lines
using Gabor filters and using deformable models to obtain the final segmentation. These
methods though are not adequately accurate.
Manual input is normally used to segment the myocardium in at least one image of
the sequence. This is normally the end-diastolic time-frame image. This contour is then
propagated to the rest of the image sequence using deformable models [56,57] and non-rigid
registration [58]. Once the myocardial boundaries have been segmented, the next procedure
is tag tracking.
18
Figure 2.1: A result of myocardial segmentation in a SA slice through the cycle. The contours
at ED and ES were manually drawn and then propagated using non-rigid registration.
19
2.1.2 Tag Tracking
This method involves the use of tag lines as dark intensity lines. To reconstruct motion,
the tag lines are to be tracked over time. In most methods, tag lines are detected either
manually or semi-automatically [59]. The tags for image frames are obtained by matching a
template of the tags obtained in the previous time or a user-defined template for the most
probable tags. The tag lines themselves can be detected using various ways ? template
matching [44] and filters [46,54,60,61]. Methods that do not require user-defined myocardial
contours [45] have also been suggested. Once the displacement vectors are obtained for
the tag lines, data is interpolated to obtain the full-field displacement data. The major
disadvantage in these methods is that interpolation is required to obtain a dense displacement
field. Manual intervention to correct the tracked tag lines is often required.
Template Matching
Guttman et al. used morphological image processing [44] to detect tag lines in radial
tagged images. The first step is the detection of myocardial boundaries detected after pro-
cessing a morphologically closed image to remove the tags. The tag lines are next detected
in the initial image of the sequence. To detect the tags in the next image of the sequence, an
a priori position is first determined based on the expected rotation and the previous position
of the tags. A search window is then devised around the determined a priori point to find a
point in the image with a tag profile closest to the tag profile of the applied tag lines. This
method of tag tracking is a very basic method; more complicated and accurate algorithms
have since been developed.
ML/MAP Estimation of Tag Lines
Denney explained a procedure for automated tag identification and tracking in [45].
This technique does not require user-defined segmentation of the myocardium. The user
needs to specify a region of interest (ROI) enclosing the myocardium at the first time frame.
20
The tags at the first image are obtained from the imaging parameters. This method, called
the ML/MAP method, consists of three stages.
First, a set of candidate tag line centers are estimated across the entire region of interest
(ROI) with a snake algorithm. The initial snake is determined based on tags from the
previous time frame. This snake is then deformed to minimize an energy function consisting
of four terms:
E(?j) =
Ncnull1null
i=0
[?ijL(?ij,wij)
+?ijEsep(?ij+1 null?ij) +?ijEsep(?ij null?ijnull1)]
+Estretch(vj) +Ebend(vj).
(2.1)
The first term consists of image forces applied to each snake point determined from a maxi-
mum likelihood (ML) estimate of the tag center?. ?ij is theith tag center on thejth tag line.
Nc is the number of centers in each tag line. w is a set of image pixels along a line perpen-
dicular to the tag line. L(?ij,wij) is the log-likelihood function that is to be minimized to
match the template of the estimated tag center with the expected tag template. This forces
the snake to converge to the dark intensity lines in the image. ?ij is a normalizing function
that prevents the weight of the first term in the overall function from being too large.
The second term enforces a constraint to maintain separation between adjacent tag
lines.
Esep(x) =
null
nullnullnull
nullnullnull
0, for xnulldsep
1
2(dsep nullx)
2, for x<dsep
(2.2)
where dsep is the minimum distance between tag points. ?ij is the weght of the separation
constraint and is set to zero if ?ij is outside the myocardium.
The third and fourth terms are constraints on the snake to control the stretching and
bending of the snake. vj is the displacement of the tag line ?j from its initial undeformed
position pj such that ?j = pj +vj.
21
After optimization, the tag centers for the present time frame are obtained. Once the tag
centers in the ROI are obtained, next it is determined if the tag center is in the myocardium
or not. The tag center estimation algorithm above is not completely accurate, which may
lead to points in the myocardium that are not part of the tag tissue to be a candidate tag
center. Two MAP tests are applied to test if the tag centers are part of a tag line. Three
mutually exclusive hypotheses are considered:
null Hypothesis T: The neighborhood contains tagged myocardium. The signal model under
hypothesis T depends on the tag pattern, proton density, T1, T2, the imaging protocol
and noise.
null Hypothesis M: The neighborhood contains untagged bright tissue (presumably, but
not necessarily myocardium). The signal model under hypothesis B depends on signal
samples of the myocardium and noise.
null Hypothesis B: The tissue contains untagged dark tissue such as saturated blood (in
black blood images) and air. The signal model under hypothesis B depends on signal
samples of tissue outside the myocardium and noise.
In the first MAP test, a candidate tag centre is considered to be part of the myocardium
if it is more likely to be part of a tag line (HT) or in untagged myocardium (HM) than in
the background tissue (HB). In the second MAP test, a candidate tag centre is considered
to be part of a tag line if it is more likely to be part of a tag line (HT) than in untagged
myocardium (HM) or in the background tissue (HB). In the final step, a pruning algorithm
is used to remove small, isolated clusters of points that appear and disappear with time.
This pruning step removes points that are not persistently tagged and isolated tag points
created as a result of noise in the image.
A pre-processing step of computing image statistics necessary for the ML/MAP method
is performed. Parameters such as noise variance and mean and noise of the background pixels
22
(a) (b)
Figure 2.2: Results of the ML/MAP tag tracking algorithm for a midventricular slice for
two different tag angles (0 and 90 degrees) [45].
and myocardium pixels were calculated. Figures 2.2(a) and 2.2(b) show the results of the
ML/MAP estimator for two different tag angles for a mid-ventricular slice.
Combined Tag Tracking and Strain Reconstruction
In [47], a method for computing three-dimensional (3D) myocardial strain from tagged
cardiac magnetic resonance (MR) images called the COmbined Tag Tracking and strain (E)
Reconstruction (COTTER) algorithm is presented. This algorithm fits a deformation model
directly to the image data, thus combining the tag tracking and deformation reconstruction
into one step. The deformation model ensures that the tag lines identified in the image
data are consistent from slice to slice. The initial geometry of the deformation model is
determined from user-selected, end-diastolic elliptical regions of interest (ROI). The COT-
TER algorithm then optimizes the control points of the cylindrical B-spline model so that
tag planes deformed by the model match the tag lines in the image data. The myocardium
motion field and strains can be reconstructed directly from the results of the combined tag
tracking procedure without intermediate steps, such as inverse and forward motion field
fitting.
23
Figure 2.3: Results of couple B-snake tracker on a simulated image sequence [63].
Analysis of Tagged MR Cardiac Images with B-splines
In [62], tag lines are tracked with dynamically programmed B-snakes. Oriented tag
templates were first constructed by simulating the MR imaging process. The templates are
subsequently correlated with the image data, giving rise to an energy landscape, which is
optimized over with dynamic programming. Tracking the valleys of the energy landscape
with snakes allows for measurement of MR tag deformations due to tissue motion.
In [63,64], coupled B-snake grids and constrained thin-plate splines are used for anal-
ysis of 2-D tissue deformations. Coupled snake grids are a sequence of spatially ordered
deformable curves represented by B-spline bases, which respond to image forces, and track
non-rigid tissue deformations from SPAMM data. The spline grids are constructed by hav-
ing the horizontal and vertical grid lines share control points. By moving a spline control
point, the corresponding vertical and horizontal snakes deform. This representation involves
shared control points. B-spline bases provide local control of shape, compact representation,
and parametric continuity. Efficient spline warps are proposed which warp an area in the
plane such that two embedded snake grids obtained from two tagged frames are brought
into registration, interpolating a dense displacement vector field. The reconstructed vector
field adheres to the known displacement information at the intersections, forces correspond-
ing snakes to be warped into one another, and for all other points in the plane, where no
information is available, a C1 continuous vector field is interpolated (thin-plate splines).
24
(a) (b)
Figure 2.4: (a) Tag grids formed from implicit snakes determine a B-solid. (b) Heart volu-
metric model located in a B-solid [65].
In[65], anapproachfor3DtaglocalizationandtrackingofSPAMMdatabyadeformable
B-solid is described. The solid is designed in terms of a 3D tensor product B-spline. The iso-
parametric curves of the B-spline solid (termed implicit snakes) deform under image forces
from tag lines in different image slices. The localization and tracking of tag lines is performed
under constraints of continuity and smoothness of the B-solid. The framework unifies the
problems of localization, and displacement fitting and interpolation into the same procedure
utilizing B-spline bases for interpolation. To track motion from boundaries and restrict image
forces to the myocardium, a volumetric model is employed as a pair of coupled endocardial
and epicardial B-spline surfaces. To recover deformations in the LV, an energy-minimization
strategy is used where both tag and LV boundary data are used.
In [66], tags are tracked using surfaces represented by thin-plate splines. These surfaces
represent tags in one direction and are obtained using a stochastic optimality criterion. These
surfaces are estimated using tags in both past and present images maintaining temporal con-
tinuity (temporal filter). A continuous representation of each tag surface is also maintained
25
Figure 2.5: (a) 3 sets of tag surfaces are assembled into a deforming 3D grid; (b) the
intersections of the grid within the LV wall mark material points of tissue (c) these markers
move with the tissue, depicting 3D contraction [66].
(spatial filter). The method generates a continuous three-dimensional reconstruction of de-
formed tag planes. Tag plane movement is modeled as a time update process, where each
update consists of a quadratic component that describes the bulk motion of the tag planes
plus a stochastic component that characterizes the fine-scale variation in the motion.
Model Tags
A method is described for the reconstruction of 3-D heart wall motion directly from
tagged magnetic resonance images, without prior identification of ventricular boundaries or
tag stripe locations in [67,68]. The method utilizes a finite-element model to describe the
shape and motion of the heart. Initially, the model geometry is determined at the time of
26
(a) (b)
Figure 2.6: (a) Model tags at the time of creation [67]. The two lighter tags are from short-
axis images while the dark tag is from a long-axis image. (b) Model fit to guide points
at ED [67]. Endocardial surface shaded, guided data as +. Corresponding model locations
(joined to each guide point with a line) are shown with cross-hairs.
tag creation by fitting a small number of guide points which are placed interactively on the
images. Model tags are then created within the model as material surfaces which define
the location of the image tags. An objective function is derived to measure the degree
of match between the model tags and the image stripes. The objective is minimized by
allowing the model to deform directly under the influence of the images, utilizing an efficient
method for calculating image-derived motion constraints. The model deformation can also
be manipulated interactively by guide points.
2.1.3 3D Motion Reconstruction
Once the myocardium has been segmented and the points tracked, the 2D motion vectors
are combined to obtain a 3D deformation field. The most common method involves fitting a
model to the calculated 2D displacement estimates. This method also involves optimization
and smoothing of the deformation field, for both interpolating missing data points and to
compensate for erroneously tracked points. Many 3-D LV deformation models have been
27
proposed, including methods based on prolate-spheroidal basis functions [69], finite elements
[67], finite differences , free-form deformations [70,71], and B-spline functions [47,65,72?74].
Combined tag tracking and deformation field fitting algorithms have been proposed that
fit the deformation field directly to the image data [47,54]. This alleviates the need for
smoothing, but this method is sensitive to registration errors between slices.
TEA: Tag Strain(E) Analysis
For the 3D reconstruction of the displacement field, a series expansion is expressed in a
dimensionless prolate spheroidal coordinate system [?null,?,?]. ?null represents wall depth and
has values of zero at mid-wall, -0.5 at the endocardial surface, and 0.5 at the epicardial
surface. ?,? represent azimuthal and circumferential angles, respectively. Displacement
?? in the radial coordinate direction can be represented using associated Legendre poly-
nomial functions and ar,l,m, the unknown coefficients. The number of terms included in
the series expression is given by the parameters nr and nl. Similar expressions are used to
model the angular displacement in the longitudinal and circumferential directions, ?? and
??. The measured displacements in X, Y and Z directions are projected onto the local
prolate spheroidal coordinate unit vectors, and the ??, ?? and ?? series coefficients are de-
termined simultaneously. The backward displacement field is inverted at the points where
the strain is computed using a gradient descent approach. Iterations are stopped when an
arbitrary precision is reached. The fit to the prolate spheroidal coordinate-based model is
preceded by a correction for the bulk motions of the ventricle. This is necessary because the
origin and long axis of the prolate spheroidal coordinate system are singularity points which
should be kept out of the myocardium. This 3D kinematic model functions to interpolate
the myocardial motion between samples in space, to filter out noise in the sampled data and
to facilitate the computation of the local displacement gradients used for strain calculation.
28
DMF: Discrete Model-Free
The discrete model-free (DMF) algorithm [14] uses finite-difference analysis to recon-
struct displacement and strain from tag lines and user-defined myocardial contours. For each
imaged time frame, the DMF first constructs a 3D segmentation of the LV wall on a regularly
spaced discrete grid of points from the user-defined endocardial and epicardial contours. The
grid spacing is equal to the pixel size in the MR image. A displacement vector is estimated
for each grid point in the LV segmentation that maps the grid point back to its end-diastolic
position. The displacement vectors to be estimated at a certain time are obtained from
the one-dimensional constraints obtained from the tag line data. Since the grid spacing
is small compared with the tag line spacing, there are many more unknown displacement
vectors than tag line constraints, and a smoothness constraint is required. The smoothness
constraint minimizes the spatial variation of the displacement field. This is imposed by set-
ting the differences between pairs of displacement vectors on the spatial point grid equal to
zero-mean random variables. The linear minimum mean-square estimate from the tag line
and smoothness constraints is then obtained by solving a system of equations containing a
sum of the two constraints. A weighting parameter controls the resulting smoothness and
the spatial resolution of the reconstructed displacement field. Strain is computed by first
converting the reconstructed displacement field to material coordinates using linear inter-
polation. The displacement gradient is then computed using a finite-difference derivative
approximation averaged over a 3null3null3 point neighborhood. Finally, the Lagrangian strain
tensor is computed from the displacement gradient.
TTT: Tag Tissue Tracking
Initially, for each short axis image plane, time frame and short-axis tag stack, a two-
dimensional B-spline field is formulated to approximate the one-dimensional displacement
of a tag surface. Once displacement fields of both tag stacks at a given short-axis image
plane and time frame are calculated, two components of the past total trajectory of any tag
29
point can be determined. The final inverse motion field gives for each point the 3D vector
pointing to its reference material point at the undeformed state. The missing component,
the SA through-plane motion, is computed from LA images as follows: (a) a 2D B0-spline
field is fit for the displacement on each long-axis imaging plane; (b) this parametric field
is then evaluated at each SA image plane intersection line; (c) through-plane displacement
information coming from all LA images is used to fit a second field in each SA plane. The
same region of interest in both SA tag displacement field fitting and the last step of long-axis
analysis is used, and a single formula for a B-spline inverse deformation field is defined. The
forward deformation at the corresponding point at tagging time should be the inverse. The
forward deformation is also formulated as a three-dimensional B-spline tensor product for
a given time frame. After the 3D forward deformation field is constructed for each time
frame, a 1D time smoothing is performed on control points. The motion of the heart is then
described as a four-dimensional B-spline forward motion field.
FTEA: Four-D Tag Strain(E) Analysis
The backward transformations deform a point x onto a point for which the coordinates
are defined as a tensor product of B-spline curves, as for the B-solid defined in [65] or for the
forward transformation in TTT. The transformation is estimated by minimizing a criterion
which is the weighted sum of a data term, which is the sum of squared distances of each
deformed tag points from its tagging plane, and a smoothing term, which is the integral of the
square of the second derivatives of the backward transformations. The criterion is minimized
with a conjugate gradient method. The transformations are estimated sequentially, using
the deformation field at the previous time as the initial transformation for the conjugate
gradient. With all the backward transformations obtained, a 4D planispheric transformation
is computed. Displacements in the planispheric coordinates can be defined using a B-spline
tensor product in time and space. To compute the forward transformation at a particular
time, a least-squares criterion is estimated that minimizes the difference in the distance
30
between a tag point position and its position after performing the corresponding backward
and forward transformations at that point.
Summary
These methods were found to yield noticeable differences in the displacement and strain
[13]. DMF uses the least smoothing and was found to have the highest spatial resolution,
but it was also the most noise sensitive. FTEA uses the most smoothing and hence was
least noise sensitive but had the lowest spatial resolution. The TEA and TTT methods were
found to fall somewhere in the middle. It was also found that all the methods produced
strain maps that were useful for characterizing myocardial function with high precision. The
method chosen thus is a trade off between smoothing and noise sensitivity.
2.2 Harmonic Phase (HARP) MRI
HARP analysis concentrates on the Fourier transform of tagged MR images. Figure
2.7(a) shows a grid-tagged MR image and its Fourier transform (Figure 2.7(b)). One DC
peak at the center of Fourier space, and the harmonic spectral peaks can be seen in Figure
2.7(b). A bandpass filter is applied to isolate a spectral peak, which can be identified by the
index k. The inverse Fourier transform of the kth peak produces the complex image
Ik(y,t) = Dk(y,t)ej?k(y,t) (2.3)
where y is the location of a pixel on the image plane, t is the time sampled in discrete
steps and corresponding to one of the time frames, Dk is called the harmonic magnitude
image, and ?k is called the harmonic phase image. Since phase angles can only be measured
modulo 2pi, however, only a ??wrapped?? value of harmonic phase in the range of [nullpi,pi] can
31
be obtained, which can be denoted as
?k(y,t) = W(?k(y,t)) (2.4)
where W is a nonlinear wrapping function and ? is the HARP angle image. The HARP
angle image is shown in Figure 2.7(c). The main properties of the HARP image are:
null The locations of the wrapping transitions on the HARP angle image are closely coin-
cident with the locations of the tag lines during myocardial motion.
null The phase values on the HARP image are dense, providing more information than just
the locations of tags lines.
null The phase value of a material point is fixed from the tagging instant and does not
change, i.e., the HARP value is a material property of the underlying tagged tissue.
This phase invariance property is the basis of HARP tracking [27,49,75]. The slope of
the phase is used to compute the local tissue strain [50,76,77].
Due to the fact that the phase is wrapped, an assumption of small motion between two
adjacent time frames is needed to guarantee that the points can be uniquely associated. In
particular, if a point has a displacement of greater than one-half tag spacing in any of the tag
directions, HARP analysis does not guarantee correct displacement estimates. To be clear,
the HARP tracking algorithm will converge in all cases, but in the case of large motion, it
will converge to a point having the same pair of HARP values but is one or more tag periods
away from the true corresponding point [48].
2.2.1 HARP Refinement
In [48], a method called the HARP refinement method based on seeded region-growing is
introduced. Starting from a seed point, defined by the user, to be tracked correctly through
the entire sequence, this method tracks the motion of the entire tissue. Once the seed point
32
(a) (b) (c)
Figure 2.7: (a) a tagged short-axis image acquired near end-systole in a normal human; (b)
the Fourier transform log-magnitude of (a); (c) HARP image obtained by bandpass filtering
around one of the peaks in (b).
is tracked, its neighbors are added to a sequentially sorted list (SSL). These points are known
as boundary points. The first point on the SSL is then tracked with the initial guess for
HARP tracking being provided by the displacement of the seed point. After tracking, this
point is added to the list of seed points and its neighbors added to the SSL. The points in the
SSL are sorted based on a cost function. The cost function is defined as the wrapped phase
difference between the boundary point in the reference time frame and the initial guess
provided for HARP tracking in the target time frame. Illustration of the region-growing
method is shown in Figures 2.8(a)-2.8(f).This method was applied on tagged MR images of
the tongue. It was seen that some points that were mistracked using the traditional HARP
tracking method were corrected using the HARP refinement procedure. This method thus
makes it possible for more accurate HARP tracking.
2.2.2 Extended HARP Tracking
Due to the combination of the conical shape of the heart and its longitudinal motion,
new tissues appear during the systolic phase near the epicardial contour and other tissues
disappear near the endocardial contour. Since these new tissues are not present during
the time frame with the undeformed tag pattern, their angle values are not available and,
33
(a) (b) (c)
(d) (e) (f)
Figure 2.8: Illustration of the region-growing process [48]. (a) is the checkerboard image at
the first time frame by overlaying two perpendicularly oriented tagged images. (b) is the
checker board image at the second time frame. The red dot in (a) is the manually placed seed
point. (c-f) shows how the region grows. The green color means tracked points, brown means
boundary points, and blue means points that are not tracked and not boundary points.
34
consequently, these points are not tracked. However, disregarding these points can lead to
systematic errors in the strain results, due to the transmural gradient in the strain. In [77],
an extended HARP tracking method is described, which tracks the new points that enter
into the image plane during the systolic phase and recovers the points that reappear at the
endocardial contour during the diastolic phase. With this extended method, the authors aim
to track all parts of the myocardium completely during all phases of the cardiac cycle.
First, the myocardial contours for the target and reference time frames are computed.
All points inside a square region encompassing the myocardial contours and points within
an additional distance of two tag lines from the contours in the reference time are tracked.
The resulting tracked points are then classified into active and inactive points. Active points
are those that fall inside the myocardial contours in the target time and also satisfy a set
of conditions relating to the topology. HARP tracking can be erroneous in areas of large
displacements; the topology-based conditions remove these points off the active list. After
tracking all active points, the absolute displacement field is computed and used to estimate
the displacement of the inactive neighbors. If a point is currently inactive but one or more of
its near neighbors is active, the displacement of this inactive point is set equal to the averaged
displacement of its active neighbors. The resulting estimated displacement is applied to the
inactive point. If the estimated position falls inside the current contours, the point is made
active and it is tracked. This procedure is continued for all the time frames. Figure 2.9(b)
demonstrates the increased myocardial coverage as a result of the extended version of HARP
tracking.
Comparing the strain results obtained with both HARP tracking versions, it was ob-
served that, on average, the extended version leads to less circumferential shortening. The
decrease is expected, since the myocardium points near the endocardium have higher cir-
cumferential shortening and the original HARP tracking version is biased toward this region.
Therefore, in comparison with the original HARP tracking version, this extended version
leads to strain results that represent more of the myocardium.
35
(a) (b)
Figure 2.9: Visualization of the rectangular region of tracking encompassing the left ventric-
ular short-axis image at end-systole [77]. The red curves represent the myocardial contours,
the blue marks the inactive points and the yellow marks the active points. Since the en-
docardial contour is difficult to delineate, a point that is successfully tracked and falls less
than two pixels in distance outside the endocardial contour is labeled as active. (a) Original
HARP tracking version. (b) Extended HARP tracking version.
2.2.3 3D HARP
3D-HARP [27] is a fast, semiautomatic method for tracking three-dimensional (3-D)
cardiac motion from a temporal sequence of short-axis and long-axis tagged MR images.
This method extends the HARP approach, previously described for two-dimensional (2-D)
tag analysis, to 3-D. A 3-D material mesh model is built to represent a collection of material
points inside the left ventricle (LV) wall at a reference time. Harmonic phase, a material
property that is time-invariant, is used to track the motion of the mesh and deform the mesh
along with an additional constraint to control the internal properties of the mesh.
The basic algorithm is as follows. The user initializes the 3D material mesh (Figure
2.10) at the initial time frame. In the next time frame, intersection points of the mesh with
the image are determined. If the phase-time invariance is not satisfied, the in-plane motion is
determined using 2D HARP tracking. Based on the tracking, the material mesh is deformed.
This is illustrated in Figure 2.11. In (a), the material mesh intersects with an SA image
plane S. The intersection point y(?) is located on both the image plane and the material
mesh. ynull is the closest point to the intersection point y(?) on the image plane S that has
36
Figure 2.10: A material mesh built at the reference time [27]. S is an SA image plane and
L is an LA image plane.M is the material mesh. The material mesh is initialized as latitude
circles and longitude lines.
37
Figure 2.11: Diagram of the intersection of the material mesh and image planes, as well as a
description of the in-plane motion of intersection points in (a) SA case and (b) LA case [27].
38
the same 2-D phase value as the initial phase value, of ?. The displacement vector between
y and ynull describes the in-plane 2-D motion of the intersection point y(?). Similarly, in (b),
the material mesh intersects with an LA image plane L at the intersection point y(m), and
ynull is the closest point to the intersection point y(m) on the image plane L that has the same
1-D phase value as the initial phase value of m. The displacement vector between y and ynull
describes the in-plane 1D motion of the intersection point. This procedure is repeated till
the phase-time invariance is satisfied to obtain the material mesh at that time. This material
mesh is then used as the initial mesh for the next time to be tracked.
Various motion-related functional properties of the myocardium, such as circumferential
strain and left ventricular twist, are computed from the tracked mesh. This method though
is less robust to large interframe deformations and can track motion only in the middle third
of the myocardium.
2.3 Optical-Flow Methods
Optical-flow methods have been developed for the analysis of tagged cardiac MR images
[41?43]. These methods use low-frequency tag patterns, modeled as sinusoidal brightness
patterns. Dense motion estimates can be obtained using gradient-based optical flow motion
estimation techniques. The brightness constraint equation (Equation (2.5)) is the basis for
this method. The equation is
?xu+?yv+?t = 0, (2.5)
where ? denotes the pixel intensity and (u,v) are the velocity components to be determined.
This method involves estimating the spatial gradient of the intensity image, using finite
differences of the intensities at neighboring pixels, as well as a temporal intensity derivative
which is estimated by subtracting the intensity at a pixel from one time frame to the next in
a cine acquisition. The resulting formulation requires the solution of a pair of coupled partial
39
differential equations (PDEs) to find the velocity field estimate. This method assumes that
a material point possesses the same brightness throughout the entire cycle, but tags fade
over time in MR images, and this assumption fails. Variable-brightness optical-flow methods
have been suggested that model the tag fading over the time sequence and produce more
accurate displacement fields [43]. One problem with these methods is that, to measure large
motions, like the ones between end-diastole and end-systole, the motion between intermediate
time-frame images has to be measured, and tracking is necessary to obtain total motion.
2.4 Non-rigid Registration Techniques
Non-rigid registration methods use a similarity criterion and transformation models to
estimate the deformation in a sequence of images. Chandrashekara et al. [70] first developed
a method for computing myocardial motion using non-rigid registration. Recently, robust
methods for computing displacement fields in tagged MRI [78] as well as from cine MR
images [58] have been developed.
Chandrashekhara et al. [70] used normalized mutual information as the similarity mea-
sure to obtain the displacement field. This procedure involves the combined nonrigid regis-
tration of SA and LA images. Initially one volume is built from the SA stack and another
from the LA stack of images. Mutual information is then used to obtain the local affine
deformation field between two volumes from adjacent time frames. The similarity measure
used is a summation of the mutual information from the registration of the SA stack vol-
ume and the LA stack volume. The affine transformation obtained between successive time
frames is then aggregated to obtain the displacement field back to end-diastole.
Ledesma-Carbayo et al. [78] use a multi-source nonrigid registration algorithm to com-
pute the displacement field from vertical and horizontal tagged images. A series of trans-
formations between successive pairs of images is defined. The transformation field was then
optimized based on an energy term comprising two terms. The first term (data term) com-
putes the sum of squared differences between the reference image and the deformed image
40
from the two sources. The second term (regularization term) was chosen as the vector Lapla-
cian of the resulting transformation. This scheme also used multiresolution to improve the
speed and robustness of the process.
2.5 Conclusion
All of the above methods are ill posed: thermal noise and image artifacts can cause
ambiguities in the tracking of tag, optical flow or phase features. The above tracking methods
tend to be problematic near the boundary, where the epicardial and endocardial surfaces
often appear as dark lines and thus can be confused as tag features. The myocardium
near the boundary is also particularly difficult to characterize because tag lines can appear
or disappear due to through-plane motion of the myocardium. These conditions near the
boundaries can cause failures for all the algorithms in several ways:
1. The intensity template used by the feature-based methods may lock onto the boundary
instead of a tag.
2. The boundaries can occlude tags and can act as spurious matches that confound the
optical-flow solution and the non-rigid registration techniques.
3. The boundaries can introduce apparent shifts of the tagging pattern in subendocardial
and subepicardial tissue, which can cause inaccuracies in the motion recovered from
harmonic phase method.
As a result, researchers have used regularization methods to find a motion field which bal-
ances fidelity to tag features with a smoothness constraint imposed on the overall motion
field.
Tag tracking based on feature extraction has an advantage over the optical-flow, HARP
and non-rigid registration techniques, in that the tracked tags can be manually corrected.
The tag points that are detected using feature-based techniques can be moved around by
the user to the correct location after manual inspection. HARP and non-rigid registration
41
though produce denser displacement estimates and hence alleviate the need for interpolation.
2D HARP strain (without tracking) is the preferred method for computing 2D strain maps
as there is no need for tracking the material points. This method though can be erroneous;
moreover, 3D strain maps can be more beneficial for diagnostic purposes.
We have studied the procedures generally followed in the analysis of tagged MR images.
Many different techniques can be applied in the segmentation of the myocardium and the
tracking of tissue. All of these methods are plagued by inaccuracies as a result of various
imaging factors. The obtained displacement field is smoothed to alleviate the inconsistency
in the obtained deformation field. The scope of improvement in analysis techniques is limited
to the quality of images obtained. In the next chapter, we look at a few MR imaging methods
that seek to improve upon the MR tagging technique both in terms of ease of analysis and
the accuracy of measurements.
42
Chapter 3
Other MR Imaging Methods to Compute Myocardial Function
MR imaging methods based on MR tagging [11,79], phase-contrast velocity imaging
[80,81] and stimulated-echo imaging [82?84] have been used to study the regional functional
characteristics of the heart. Recent advances in imaging methods [27,85?87] have improved
the clinical viability of MR imaging techniques. Most of these advances are nascent and have
not yet found their way into regular clinical use.
3.1 3D Myocardial Tagging
The 3D-CSPAMM [79] tagging preparation is an extension of the 2D-CSPAMM se-
quence. CSPAMM is based on the subtraction of two images with complementary signed
tagging modulation. The subtraction technique reduces tagline intensity fading and con-
sequently allows observation of the heart motion throughout the heart cycle. The 3D-
CSPAMM tagging preparation generates a 3D tagging grid on the myocardium. Two 90o
block pulses, interspersed by a dephasing gradient, produce a sinusoidal modulation of z-
magnetization and thus a line-shaped tag pattern. To generate a 3D tagging grid, the
modulation is repeated in all three spatial directions. 3D-CSPAMM is a valuable technique
to register the 3D motion of the myocardium. It overcomes the matching problems that
may arise when combining short-axis and long-axis tagging images to generate 3D tagging
images. The major disadvantage is that the scan duration is 40 minutes which is almost
double that of SPAMM tagged MR images.
43
Figure 3.1: SF HARP illustration. A thin slice is selectively tagged while imaging a larger
slab that encompasses the moving tagged slice [87].
3.2 Slice Following HARP
The slice following (SF) CSPAMM technique [87], selectively tags a thin slice at an early
time frame A while imaging a large slab that always encompasses the moving tagged slice
even at a later time frame B (Figure 3.1). After subtracting the complementarily signed
tag-modulated images, the reconstructed images represent the same material slice at all
time points irrespective of where it has moved in the through-plane direction. Therefore,
performing 2D HARP point tracking on these SF images yields a true 2D displacement of
the material point. Now consider two orthogonal SF tagged slices, both of which have two
tag orientations, and consider a point lying on the intersection of the original tag planes of
these slices. Such points have two projected pathlines that can be computed using 2D HARP
tracking on orthogonal planes (Figure 3.2), and these projected pathlines can be combined
into a single 3D pathline representing the 3D motion of an MR marker. This can be done
on any such intersection point and can be computed automatically.
3.3 zHARP
zHARP [86] is a novel method for calculating cardiac 3D strain using a stack of two or
more images acquired in only one orientation. The zHARP pulse sequence encodes in-plane
44
Figure 3.2: SF HARP Planes Illustration. Points common to SA and LA slices would have
two pathlines tracked using 2D HARP. These are combined to give a 3D pathline [87].
motion using MR tagging and out-of-plane motion using phase encoding, and is capable of
computing 3D displacement within a single image plane. Here, data from two adjacent image
planes are combined to yield a 3D strain tensor at each pixel; stacks of zHARP images can
be used to derive stacked arrays of 3D strain tensors without imaging multiple orientations
and without numerical interpolation.
The proposed technique, zHARP, is similar to the standard slice-following CSPAMM
(SF-CSPAMM), except that during the imaging sequence, a small z-encoding gradient is
applied immediately before the readout and again with opposite polarity to the second or-
thogonal CSPAMM acquisition. This gradient adds a z-position dependent phase?z to every
material point in the acquired slice. This additional phase is linearly related to the distance
of the point from the iso-center. Unfortunately, susceptibility and general field inhomo-
geneities lead to an additional and artifactual phase accumulation. This erroneous phase
is identical in both the horizontally and vertically tagged images, which can be exploited
to separately compute the in-plane and through-plane displacements. The algorithm runs
on the two orthogonally tagged data. First, k-space data are processed using HARP, upon
extracting both the positive and negative harmonic peaks. Second, the harmonic phase maps
?x,?y and ?z, in the x,y and z directions are computed. Third, material points are tracked
45
in 3-D. Given an arbitrary material point, a tracking procedure is introduced to track the
in-plane motion using HARP on the ?x and ?y maps, from which the z-position can be
extracted using the ?z map.
The zHARP technique is a modification of SF-CSPAMM, where encoding gradients
add a through-plane motion-dependent phase to every material point in the acquired slice.
Data acquisition for zHARP requires no extra time over a conventional CSPAMM acquisi-
tion; however, the spatial resolution and accuracy of this method have yet to be carefully
evaluated.
3.4 Combined Harmonic Phase and Strain Encoded MRI
HARP-SENC MRI [85] is a method to measure the dynamic evolution of three normal
components of myocardial strain within a six-heartbeat acquisition per imaging slice using a
HARP-SENC MRI pulse sequence. The calculation of the true 3-D strain tensor cannot be
obtained using this technique described, as longitudinal shear components are not measured.
The post-processing of the datasets acquired from a six-heartbeat breath-hold provides full
quantification and visualization of three-component myocardial strains with good temporal
resolution in a single SA slice. For a stack of n contiguous SA slices prescribed in the
left ventricle, these images are acquired in a single breath-hold lasting 6n heartbeats. This
method is still under development, as some artifacts have been observed in the SENC images.
3.5 DENSE
Displacement encoding with stimulated echoes (DENSE) is a quantitative MRI tech-
nique that encodes tissue displacement into the phase of the complex MRI images [83,84,88].
Cine DENSE provides a time series of these images, thus facilitating the non-invasive study
of myocardial kinematics. DENSE images are used for the quantification of regional dis-
placement and strain as a function of time.
46
Spottiswoode et al. [82], developed a method for analyzing DENSE image sequences
using a 3D flood-fill phase unwrapping method. To achieve enhanced motion sensitivity,
relatively high displacement encoding frequencies are used to capture DENSE images and
phase wrapping typically occurs. Algorithms to phase unwrap a time series of cine-DENSE
images, thus enabling the computation of absolute myocardial motion trajectories are devel-
oped. The unwrapping was implemented by extending a 2D quality-guided path following
phase unwrapping method [89] to unwrap phase through time as well as space.
Let ?(i,j) be the wrapped phase and ?(i,j) be the true, unwrapped phase. The coor-
dinates i and j define the pixel position within the image. Phase wrapping can be defined
as ?(i,j) = ?(i,j)+2pik(i,j), where k(i,j) are integers that force ?(i,j) to lie between nullpi
and pi radians. Phase unwrapping entails determining all values ofk(i,j), and it can be done
along a particular path by integrating the wrapped phase differences of the initially wrapped
phases.
The 2D quality-guided phase unwrapping method uses a measure of phase ?quality?
to guide the path of unwrapping. A good measure of phase quality is a root-mean-square
measure of the variances of the spatial partial derivatives within an nnulln pixel region. The
unwrapping path is guided by a region-growing algorithm, which floods into regions of high-
quality phase before those with lower quality. Figure 3.3 captures six steps of the 2D quality-
guided flood-fill procedure commencing from an arbitrary seed point in the myocardium, and
using the phase quality map shown in Figure 3.4(c).
To obtain absolute phase measurements, it is necessary that the unwrapping path be
initiated on a portion of tissue with known absolute phase. Very little cardiac motion takes
place between the displacement encoding and the acquisition of the first frame of cine-DENSE
data. It is therefore characteristic that no phase wrapping will have occurred on the first
frame, and that all phase measurements in the myocardium will be absolute. Unwrapping the
phase through time will thus result in measurements of absolute phase for all frames. This
implicitly assumes that the true motion magnitude between frames for each pixel separated
47
in time amounts to less thannullpinullradians. This will always be the case given practical values of
pixel size and displacement encoding frequencies. If the phase quality maps for a cine series
are stacked on top of each other, then the quality-guided algorithm can easily be adapted to
unwrap in 3D (one temporal and two spatial dimensions).
The following algorithm for unwrapping cine DENSE phase images was implemented:
1. Calculate phase quality maps for all frames.
2. Manually place a seed point anywhere in the myocardium on the first frame.
3. Perform 3D quality-guided phase unwrapping, limiting the spatial unwrapping to a
specified radius per frame. As soon as an unwrapped pixel appears on a particular
frame, disallow all unwrapping for all previous frames.
4. Once unwrapped pixels are present on all frames, perform 2D quality-guided phase
unwrappingforeachframeindependently, startingatthecorrectlyunwrappedreference
pixels.
Step 2 serves to force the unwrapping through time using only high phase quality. If the
3D quality-guided flood-fill is used for the entire unwrapping process, phase inconsistencies
are introduced at the later stages of the procedure due to temporal unwrapping between
isolated regions of high-quality phase in the noisy areas (lungs and ventricular cavities) and
regions within the myocardium. The radius size is set based on the frame rate and the
maximum expected amount of cardiac displacement per frame. It should be large enough
to allow the 3D unwrapping to migrate spatially with the ventricle as it moves, but small
enough to restrict the 3D unwrapping to relatively high values of phase. A radius of 2 pixels
was used here.
Once the images were unwrapped and the trajectories tracked, a fifth-order Fourier series
curve is fitted to the tracked trajectories. This improves the accuracy of measurements for
noisy image data.
48
Figure 3.3: Progression of the guided flood-fill operation using the quality map in Figure
3.4(c). The calculation intervals are arbitrary.
(a) (b) (c) (d)
Figure 3.4: (a) Short-axis DENSE magnitude image. (b) Wrapped DENSE phase image for
horizontal motion. (c) Variance of the derivatives of the locally unwrapped phase for the
four neighboring pixels. (d) Unwrapped image with LV endocardial and epicardial contours
superimposed. [82]
A useful method for estimating where phase-unwrapping errors may have occurred in-
volves combining a map of unwrapped phase discontinuities with the phase quality map.
The discontinuity map specifies the locations of pixels that differ from a neighbor by more
than pi radians. Unwrapping errors are likely to have occurred within regions isolated by the
combined map. The combined maps were used to aid a visual inspection of the left ventricle
(LV) for unwrapping errors. If a single error was encountered, then the image was deemed
incorrectly unwrapped. It was shown that this method correctly unwrapped 786 out of the
800 images analyzed (98.3%).
The only manual interaction required was the placement of a seed point and the my-
ocardial contours. The rest of the algorithm is automated. It was shown that the strain
measurements obtained were comparable to those obtained from tagged image data.
49
3.6 Conclusion
MR tagging of the myocardium is a powerful method for assessment of regional my-
ocardial deformation. It has permitted new insights into normal physiology and myocardial
diseases. Applications of MR tagging in a routine clinical setting are yet to be commonplace
due to the long imaging times and expensive and time-consuming post-processing proce-
dures. Automated tagged image analysis procedures are an active research area. Many of
the concepts developed in the analysis of tagged images can be useful in the analysis of images
obtained from newer and faster imaging methods. With effective solutions in both analy-
sis procedures and imaging methods, tagged MR imaging is likely to become an important
technique for clinical assessment of cardiovascular function.
50
Chapter 4
Objectives
This chapter discusses the objectives and contributions of this dissertation. As discussed
in Chapter 2, there are several problems with the existing cardiac motion tracking methods
in tagged MRI. The feature-based methods involve tracking of tags as dark intensity lines.
Mistracking of tags is a big problem and there is a need for extensive manual intervention to
correct for the errors in tracking. This procedure is tedious and cumbersome when a series of
images have to be processed. Interpolation is also necessary to produce dense displacement
fields. The phase-based methods overcome the need for interpolation. They do not require
any manual intervention but, this is a disadvantage as it does not allow for correction of
errors after tracking. HARP tracking produces incorrect results when the displacement of
the tracked tissue between consecutive time frames is more than one-half the tag-spacing
(because of phase ambiguity). The 3D imaging methods discussed in Chapter 3, require
lesser analysis time but need twice the amount of time to acquire a slice.
The objective of this research was to develop a method to compute 4D strain in the
human heart using tagged MR images. The method:
1. must require minimal user intervention.
2. should be capable of measuring strain through to mid/end diastole.
3. should include a way to recognize and fix the errors from the automatic algorithm.
4. should produce dense displacement measurements.
5. should be robust to large inter-frame motion.
51
First, the phase ambiguity problem in HARP tracking was solved by unwrapping the
HARP phase image. Different methods of phase unwrapping were explored and a suit-
able procedure for unwrapping a sequence of HARP images was identified. The residue-
compensation procedure (see Chapter 5), where inconsistencies in the phase image are re-
moved using branch cuts was chosen. A graphical user interface to help the users place
branch cuts was developed and a procedure to efficiently and accurately compute displace-
ments by unwrapping a sequence of images was devised. The procedure developed was faster
to compute 3D displacement over time and produced denser displacement fields compared
to the feature-based methods. It was more accurate compared to the HARP tracking pro-
cedure. Since it used SPAMM tagged MR images, the acquistion of a slice was twice as fast
compared to other recent 3D MR imaging methods.
Secondly, the phase unwrapping method was extended to obtain 4D strain in the RV.
Obtaining 4D strain in the RV has been problematic because of the thin myocardium and
the lack of geometric symmetry in the RV. The thin myocardium causes problems in tissue-
tracking in feature-based (too few tag points to reconstruct deformation field) and phase-
based methods (interference from myocardial borders is greater). Using the phase unwrap-
ping procedure developed for the LV (Chapter6), a dense deformation field with no tracking
errors was obtained. The DMF method (see Chapter 2) is then used to compute 3D strain
from the unwrapped-phase-derived displacement measurements. The DMF method has the
advantage of not requiring the object being modeled to have geometric symmetry and hence
produces accurate biventricular strains.
Thirdly, a computer-assisted approach was developed for the placement of the branch
cuts and unwrap a sequence of images. Here, branch cuts are automatically placed using ei-
ther a simulated annealing procedure or an exhaustive search procedure. An energy function
was proposed to impose temporal continuity of the unwrapped phases between time frames.
If errors were detected after the automated placement of branch cuts, the user corrected the
52
branch cut configuration. This procedure reduced the time taken for manual intervention to
about one-third of that of the procedure where all branch cuts were manually placed.
Fourthly, An explanation of the energy mechanics in the heart using commonly measured
parameters is explored. Parameters that are measured easily such as volume, area, strain and
torsion are used to develop a framework to classify hearts based on their energy consumption
through a cycle.
Chapter 5 reviews the problems and current techniques associated with phase unwrap-
ping. Chapter 6 presents a method to manually correct the problems and then unwrap HARP
images to obtain strain (this work was published in JMRI 2010 [90] and ISBI 2008 [91]).
Chapter 7 compares 2D and 3D torsion measured from tagged MR images and explains
the effect of through-plane motion on obtained measurements (this work was published in
ISMRM 2009 [92]). Chapter 8 extends the concept of phase unwrapping to the right ventricle
(RV) to obtain 4D biventricular strain (this work was published in ISMRM 2010 [93] and
ISMRM 2008 [94]). Chapter 9 describes a computer-assisted technique for phase unwrapping
HARP images (this work was published in ISBI 2009 [95]). Chapter 10 explains the energy
mechanics in the heart using commonly measured parameters.
53
Chapter 5
Phase Unwrapping
5.1 Introduction
This chapter starts out with an introduction to phase unwrapping followed by a com-
parison of the results to determine the best method. Other implementation issues specific
to the geometry of the heart are also discussed.
Phase unwrapping is a well established field in other forms of imaging such as SAR
interferometry [96], optics [97] and brain MRI imaging [98]. It is a well-researched problem.
Our aim here is to unwrap the HARP images [49]. As we shall see through this chapter, it
is a difficult procedure as it is sensitive to noise and phase aliasing.
When continuous phase is wrapped, it lies between nullpi and +pi. This is normally so
because phase is obtained using the arctangent function. The procedure of making the dis-
crete wrapped phase continuous to give the continuous unwrapped phase that lies from nullnull
to +null is called phase unwrapping. It is important to understand the aspects involved in
one-dimensional phase unwrapping before we proceed to two-dimensional (images) phase un-
wrapping. Let us look at a simple example of a one-dimensional phase unwrapping problem.
Consider a complex signal,
s(t) = ej5pit, 0 nulltnull 1 (5.1)
We wish to reconstruct the continuous phase ?(t) = 5pit from s(t). We can only obtain
the wrapped phase estimate ?(t) by using the arctangent function of the imaginary phase
(null(s(t))) over the real phase (null(s(t))) to extract phase.
54
0 0.2 0.4 0.6 0.8 1?4
?2
0
2
4
6
8
10
12
14
16
t:Relative time
Phase in radians
Unwrapped Phase
Wrapped Phase
wraps
Figure 5.1: Illustration of one-dimensional phase unwrapping
?(t) = arctan null(s(t))null(s(t)) (5.2)
An illustration of one-dimensional phase unwrapping is given in Figure 5.1. The phase
wraps at pi and 3pi can be clearly seen. It is also clear that by adding 2pi at the first phase
wrap and to every subsequent point thereafter until the next wrap occurs, the wrapped and
the unwrapped phase would match. By adding 2pi every time the phase wraps (drops by)
2pi, we can reconstruct the continuous phase. The phase wraps themselves are detected by
observing the phase differences. A simple unwrapping operation can be defined as follows:
?(n)null?(nnull1) nullpi null ?(n) = ?(n)null2pi
?(n)null?(nnull1) nullpi null ?(n) = ?(n) + 2pi (5.3)
55
There are many problems associated with one-dimensional phase unwrapping that are
of interest to us. First is the problem of phase aliasing [99]. If the data is sampled at a
lower rate such that one of the phase wraps is completely missed, then the phase calculated
thereafter is all wrong. This is a particular problem that we experienced in tagged cardiac
MRI. For example, in Figure 5.3(a), we can see two tag lines merge (at 4 ?o clock). This
would amount to missing of a phase wrap while unwrapping. Another problem is the noise
associated with the signal. If the signal is noisy and a phase wrap is not detected in the
unwrapping procedure, the unwrapped phase would be faulty. These are two of the most
important factors that ultimately decide the applicability of phase unwrapping to a given
problem.
The easiest way one can think of implementing phase unwrapping for images is to extend
the one-dimensional algorithm to the two-dimensional case, but there are some problems that
are introduced when we want to use phase unwrapping in images that make it ineffective.
Shearing of the phase plane is one of the problems. This gives rise to both problems of phase
aliasing and introduction of singularities in the wrapped phase. Singularities are undesired
points of discontinuity. Rotation of the phase plane too is a problem as it introduces phase
aliasing. These are two of the more difficult ones to deal with because of the shearing and
rotation the heart undergoes at and near end-systole. In the regions of blood flow and air
in the image obtained, the MR signal is corrupted leading to a noisy phase. This noise in
the phase obtained at points outside the myocardium makes the algorithm highly sensitive
to the myocardial segmentation. It is important to realize that any algorithm designed
should deal with all these problems in a satisfactory way. With this introduction, we shall
look at the two-dimensional unwrapping terminology and the various methods that were
implemented [89]. Phase unwrapping has already been used for motion estimation using
CINE DENSE [82]; here we try to use it for SPAMM-tagged images.
56
5.2 Two-dimensional Phase Unwrapping Terminology
5.2.1 Residues
Residues arise out of singularities in the wrapped phase or zeros in the complex function
used to derive the wrapped phase. They are the cause for inconsistencies in the image. A
simple method for identifying the residues is by unwrapping around a small 2 null 2 block
around a pixel.
res =
null
2null2
null? (5.4)
This is the smallest loop one can obtain in an image. The sum of the phase differences
around the 2 null 2 block should be zero. If the result of the unwrapping, res, is not zero,
then a residue is said to be present. Residues in general add up to either +2pi or null2pi
and are respectively known as positive or negative residues. Residues are a major problem
in phase unwrapping and the presence of residues implies that the unwrapping is no more
path-independent. Removal of residues is essential to obtain unwrapped phase without any
discontinuities.
5.2.2 Dipoles
It is common to see positive and negative residues next to each other in the image.
These are known as dipoles and are removed from the image. The basic concept of residue
removal is the neutralization of the charges on the residues; connecting a positive residue to
a negative residue and vice versa. A dipole is already balanced.
5.2.3 Branch Cuts
The connection of a positive and a negative residue in an image is known as a branch cut.
A decision must be made as to which positive and negative residues the branch cut should
57
(a) (b)
Figure 5.2: (a) The original MRI image; (b) The residues are marked black and white for
negative and positive residues respectively. The number of residues is clearly larger in the
outer areas where there is air and areas of the image with large motion compared to the
areas with tissues as seen in the original MRI image.
connect. In general, every image has an equal number of positive and negative residues. The
two most popular and widely used methods for balancing the residues are:
null minimizing the aggregate length of the branch cuts
null choosing the shortest branch cut available every time
In our algorithms, the objectives are:
null to obtain the unwrapped phase only inside the myocardium
null not to isolate any portion of the myocardium
null keep the myocardium on a whole connected
These objectives are important when using path-following methods to do the unwrapping.
5.2.4 Quality Maps and Masks
In the path-following methods of phase unwrapping, it is necessary to base the inte-
gration on a certain quantity that will improve accuracy. This is done by the use of phase
58
quality maps [100]. The two most common quantities used to describe the quality are the
phase derivative and the phase derivative variance. In the phase derivative quality map, the
wrapped phase derivative at each point represents the quality map. The phase derivative
variance map z was obtained as follows:
zm,n =
nullnull
(?xi,j null?xm,n)2 +
nullnull
(?yi,j null?ym,n)2
k2 (5.5)
where for each sum i,j varies over the knullk window centered at the pixel (m,n). The terms
?xi,j and ?yi,j are the wrapped phase derivatives, and ?xm,n and ?ym,n are the average of the
derivatives in the knullk window. Normally, 3null 3 window is used. Both the methods were
tested and it was found that the phase variance performed better but had a longer execution
time by a factor of 20 times. The phase derivative quality map was less accurate but faster.
Quality maps can be thresholded and the resulting binary image used as a mask in both
the path-following methods and the least-squares methods. Thresholded quality maps are
also used to segment the myocardium when only the myocardial phase is of interest.
5.3 Two-dimensional Phase-Unwrapping Algorithms
Phase unwrapping algorithms can be classified into two main categories, the path fol-
lowing methods and the least-squares methods. Several different algorithms from each type
were implemented and their performance studied, [96,101?105]. We shall only take a look
at the basic ones and those that were implemented here.
5.3.1 Path-Following Algorithms
The path-following algorithms are a set of algorithms that perform phase integration so
as to cover all the pixels in the image. They can be classified into three groups:
null path-dependent algorithms
59
null quality-guided algorithms
null residue-compensation algorithms
Path-Dependent Algorithms
These are the simplest algorithms and perform integration over a fixed path. They use
a predefined strategy to do the unwrapping such as linear scanning, spiral scanning and
multi-directional scanning. These algorithms follow a fixed order where unwrapping starts
from a seed point and follows to a pixel in the predefined path. This continues till all the
pixels are unwrapped. This strategy implies that if unwrapping of one of the initial pixels
is wrong, all the pixels unwrapped henceforth are wrong. Moreover, the fixed order implies
that it does not have the ability to avoid residues present in the path.
Quality-guided Algorithms
The quality-guided algorithms are a very robust set of algorithms as far as dealing
with noisy phase data is concerned, but, since they are path-following methods, they fail
to unwrap correctly when there are residues in the image. The underlying logic in every
quality-guided method is simple, and they are generally fast methods. These characteristics
make them popular methods to use once it has been verified that there are no residues in the
image. For all the algorithms, a seed point is chosen initially. The choice of the seed point
is random in most applications. From the seed point, the algorithm travels to a pixel of
highest quality among its neighbors ( a 4-pixel neighborhood is ususally considered) and so
on till all the pixels in the image have been traversed. The quality being assessed is derived
from the quality maps. It can be the pixel of lowest phase-derivative variance or the lowest
phase derivative with respect to the present pixel or any other such measure. In our case,
the variance map is chosen, as it is the most accurate and fairly robust in its unwrapping. It
has to be used in conjunction with some other method to remove any residues or unbalanced
charges that remain in the unwrapping path. This is explored in detail in the next chapter.
60
Residue-Compensation Algorithms
These algorithms place branch cuts to balance the residual charges. The branch cuts
are placed, the unwrapping is path independent and any path can be chosen to unwrap the
images. In reality, the criterion to choose a branch cut is not so easy, and one cannot be
totally sure that the resulting phase is correct [106]. In our case it is even more difficult, as we
want all the areas inside the myocardium to be connected and therefore avoid all branch cuts
from passing through the myocardium. We prefer to unwrap the pixels in the myocardium
only because the pixels outside are generally too noisy. The difficulty is choosing the right
branch cut with respect to the wrapped phase observed in the heart. This is explained in
detail in future chapters.
We can look at an example where we look to unwrap the wrapped image from Figure
5.3(a), which has already been segmented. The point of residues clearly are the places where
the two phase (wraps) lines meet (4 o? clock) and the premature ending of the phase (wrap)
line near the epicardium (10 o? clock). It can be clearly noticed that, without branch cuts,
the phase is discontinuously unwrapped as in Figure 5.3(b). Figure 5.3(c) shows the result
that is unwrapped with branch cuts. The unwrapped phase appears smooth and would seem
to be the right phase when in actuality it is not. By rearranging the branch cuts as in
Figure 5.3(d), we again obtain a smooth phase but this time it is correct as compared to the
expected phase.
It is a fact that branch cuts may lead to smoothly unwrapped phase, but the unwrapped
phase is not necessarily correct. Therefore the residue-compensation methods must be used
with care. While they make the process on the whole slightly more accurate, they also make
it unpredictable. A procedure for placing correct branch cuts to obtain correctly unwrapped
phase is explained in Chapter 6.
61
(a) (b)
(c) (d)
Figure 5.3: (a) The wrapped HARP phase to be unwrapped; (b) Unwrapping done without
branch cuts and using the Quality-guided method; (c) Unwrapping done by introducing
wrong branch cuts; (d) Correctly unwrapped image by introducing the correct branch cuts.
62
5.3.2 Least-squares Methods
The second major approach to phase unwrapping is the application of a mathematical
formalism to reduce the problem to a least-squares method [97]. There are both weighted and
unweighted procedures for achieving the same. In the least-squares methods, the gradient
field is determined from the wrapped phase differences of the wrapped phase data and is
then integrated subject to the constraint of a smooth field. This calls for the minimization
of
nullnull?nullW[null?]null22, (5.6)
where ? is the unwrapped phase. It can be shown that the minimizer in Equation (5.6) is
the solution to a Poisson equation
null2? = nullnullW[null?] on ? (5.7)
with appropriate boundary conditions ?. Hunt?s matrix formulation [107] is used to convert
Equation (5.7) to a discretized Poisson equation of the form
(?i+1,j null2?i,j +?inull1,j) + (?i,j+1 null2?i,j +?i,jnull1) = ?i,j, (5.8)
where
?i,j = (?xi,j null?xinull1,j) + (?yi,j null?yi,jnull1)
and ?xi,j and ?yi,j are wrapped phase differences i. e. discretized versiion of the x and y com-
ponents of the gradients at pixel (i,j). The discrete cosine transform (DCT) is the most
63
(a) (b) (c)
Figure 5.4: (a) The wrapped HARP phase to be unwrapped; (b) The phase after unwrapping
using the DCT method; (c) The unwrapped phase image using the Preconditioned Conju-
gate Gradient method of unwrapping. The PCG algorithm uses weights to mask out the
myocardium unlike the DCT algorithm.
popular method used to solve the discrete Poisson equation of Equation (5.8). It is computa-
tionally efficient since it is based on the DCT algorithm which is both memory efficient and
very fast. It is also easy to implement as it does not require any prior information (quality
map, identify residues etc.). The DCT algorithm can be described as follows:
1. Compute the array of values for ?i,j as in Equation (5.8)
2. Compute the two-dimensional DCT of ?i,j to get ??i,j
3. Replace the values in ??i,j to obtain ??i,j as in Equation (5.9). This makes it memory
efficient.
??i,j = ??i,j
2cos(pii/M) + 2cos(pij/N)null4 (5.9)
4. Compute the two-dimensional inverse DCT of ??i,j to get ?i,j, the unwrapped phase
The preconditioned conjugate gradient (PCG) method is an improvement over the DCT
method in that it allows the user to give a weight matrix so that the user can supply the
quality map that describes the areas with high signal-to-nose-ratio. The algorithm works
iteratively by use of the full-lattice solution [89] as the preconditioner in the conjugate
64
gradient algorithm. The full-lattice problem is solved every iteration using the DCT method.
In our case, the segmented myocardial pixels were designated the value of highest quality of
one, and the others were designated to zero.
The results for the two implementations of the least-squares methods using the DCT
and PCG methods are as shown in Figure 5.4(b) and Figure 5.4(c) for a wrapped phase
given in Figure 5.4(a). The data points only inside the myocardium have been shown for
the sake of easier analysis. It is clear that the PCG method performs better than the DCT
method. The PCG method is fast and is a good choice compared to the DCT method, as
it gives a more accurate solution. The problem is that this algorithm is sensitive to large
inconsistencies in the data (large amount of residues) and does not always provide a solution
that allows for pixel-accurate measurements. It was tested on tagged MR images, and it
showed artifacts in the unwrapped phase images that reduce its usefulness.
A significant problem experienced with the least-squares methods is that an integral
multiple of 2pi is not necessarily added to ensure smoothness and this may lead to unwanted
smoothing and increased differences between the original and obtained results. Least-squares
methods show better performance than path-following methods when the data to be used
is very noisy and shows little or no aliasing. The data used in this research generally has
high SNR inside the area of interest (the myocardium) and the major problem as opposed to
noise is aliasing. The unwrapped phase data that we require also needs to be highly precise
to get pixel-accurate displacement measurements. Least-squares approaches are therefore
not pursued.
5.4 Conclusion
The chapter reviewed the basics of phase unwrapping and how the different methods
work on HARP images [108]. It has been recognized that residues caused mostly by phase
aliasing are the major hurdle facing unwrapping the HARP images. It is important to realize
that the residues occur when one can notice two wrapped phase lines (those that occur when
65
there is nullpi to pi shift) merge together or when a wrapped phase line ends prematurely. This
is an important concept that shall be utilized later.
It has been shown that the path-following methods are a better choice than the least-
squares methods at least when the input phase has a fairly high SNR. The quality-guided
method did not account for any residues present in the image. The residue-compensation
method, while being more unpredictable, also produces more accurate unwrapped phase
maps in the presence of residues when the branch cuts are placed correctly. In Chapter 6,
we explore the placement of branch cuts in harmonic phase images and in Chapter 9 we
automate this procedure.
66
Chapter 6
3D Left Ventricular Strain from Unwrapped Phase: A GUI-based approach
6.1 Introduction
Parameters of cardiac left-ventricular (LV) mechanical function are important for di-
agnosing and managing patients with heart disease and assessing the efficacy of therapies
over time. Imaging methods based on magnetic resonance (MR) tagging [11,79] and dis-
placement encoding with stimulated-echoes (DENSE) [82,88] have been used to study the
regional functional characteristics of the heart. While recent advances in imaging meth-
ods [27,85?87] and post-processing techniques [49,77] have improved the clinical viability
of these techniques, the ability to provide dense 3-D strain maps and torsion measurements
throughout the cardiac cycle in a reasonable amount of time is still an active area of research.
Tagged MRI is an established method for measuring parameters of LV deformation and
strain. Tagged MRI spatially modulates the longitudinal magnetization of the underlying
tissue before image acquisition. The result is a periodic tag pattern that deforms with the
tissue as shown in Figure 6.1a. Several techniques have been developed to measure LV de-
formation and strain from the deformation of the tag pattern. These techniques include
feature-based (FB) methods [14,44?47] , optical-flow and non-rigid registration based meth-
ods [43,70,78,109] , and HARP-based methods [27,49,51,52].
In HARP analysis, the tag pattern deformation is measured by the local change in phase
of the tag pattern. The tag pattern modulates the signal intensity of the image, producing
a series of peaks in the Fourier domain as shown in Figure 6.1b. One of these peaks is
filtered out and inverse Fourier transformed to produce a HARP phase image as shown in
Figures 6.1c and 6.1d. The HARP phase at a point in the image is a material property of
the underlying tissue and can be tracked through the image sequence or used to compute
67
2-D strain [49,50] or 3D strain in the middle third of the LV wall [27]. The HARP phase,
however, is wrapped because it can only be measured modulo 2pi. This wrapping can cause
tracking problems if a region of the tissue deforms more than one-half tag spacing i.e. a phase
shift of more than pi between time frames. Deformations of this magnitude are possible in
both healthy and diseased hearts.
In DENSE imaging [83,110], myocardial tissue displacement is encoded into the phase
of the MR image during acquisition. To reconstruct myocardial displacement from the phase
data, the phase must be unwrapped [82,83,110], which involves adding integer multiples of
2pi to the wrapped phase so that the unwrapped phase is continuous. Noise and imaging
artifacts, however, can cause phase inconsistencies where no integer multiple of 2pi can yield
a continuous unwrapped phase map. While DENSE phase unwrapping continues to improve
[88,110], it is still possible for phase errors to occur, which can result in errors in deformation
and strain estimates.
In this chapter, a new method is presented for measuring 3D deformation and strain from
tagged MRI data based on unwrapping HARP images called Strain from Unwrapped Phase
(mSUP). In contrast to a recently-proposed 3D wrapped-phase based method in [27], the
mSUP algorithm can measure strains in the entire LV. The mSUP algorithm can also easily
measure 3D strains in the entire LV wall over multiple time frames. The mSUP approach is
more robust to inter-frame motion than techniques such as [27] that are based on tracking
wrapped phase. With unwrapped phase, points can be tracked through an image sequence
as long as the average inter-frame deformation over the entire myocardium is less than one-
half tag spacing. Compared to methods based on tracking tag lines [14,111], displacement
measurements from unwrapped HARP phase images are dense, potentially resulting in more
accurate estimates of strain.
68
6.2 Materials and Methods
6.2.1 Human Subjects
The mSUP algorithm was validated on a cohort of 40 human subjects (10 normal volun-
teers, 10 patients with myocardial infarction, 10 patients with mitral regurgitation, and 10
patients with hypertension). All human studies were approved by the Institutional Review
Boards of both institutions (Auburn University and University of Alabama at Birmingham)
and informed consent was obtained from all participants.
6.2.2 Data Acquisition
All participants underwent MRI on a 1.5T MRI scanner (GE, Milwaukee, WI) optimized
for cardiac application. Tagged images were acquired in standard views (2-chamber and 4-
chamber long-axis and short-axis) with a fast gradient-echo cine sequence with the following
parameters: FOV = 300 mm, image matrix = 224x256, flip angle = 45, TE = 1.82ms,
TR = 5.2ms, number of cardiac phases = 20, slice thickness = 8 mm. A 2D fast gradient
recalled spatial modulation of magnetization (FGR-SPAMM) tagging preparation was done
with a tag spacing of 7 pixels. This protocol resulted in 2 long-axis slices and on average,
12 short-axis slices per study.
6.2.3 Unwrapping HARP Phase
HARP images [49] were obtained by filtering out a single spectral peak in a tagged image
and computing the phase in each pixel of the resulting complex-valued image. These phase
images are wrapped because phase can only be computed modulo 2pi. Phase unwrapping
involves adding integer multiples of 2pi to each pixel in the wrapped image so that:
1. The unwrapped phase is continuous.
2. If the unwrapped phase is re-wrapped, the result is equivalent to the original wrapped
image.
69
Figure 6.1: (a) A tagged short-axis image acquired near end-systole in a normal human; (b)
the Fourier transform log-magnitude of (a); (c) HARP image obtained by bandpass filtering
the lower circle in (b); (d) HARP image obtained by bandpass filtering the upper circle in
(b).
70
The relationship between the wrapped phase and the unwrapped phase may be stated as
?(x,y) = mod null[pi+?(x,y)],2pinullnullpi, nullpi<?(x,y) <pi (6.1)
where W is the wrapping operator. In this research, phase was unwrapped only inside the
LV wall. Pixels outside the myocardium were not processed.
Several phase-unwrapping algorithms have been proposed for MRI and other applica-
tions [89] (see Chapter 5). In particular, DENSE post-processing involves unwrapping the
DENSE phase maps [82,84,88]. In this research, a path-following method was used that
is similar to the one used for DENSE in [82,89], where phase was unwrapped along a path
in the image from one pixel to another. If the phase is unwrapped correctly, the phase of
the destination pixel is independent of the path taken. Two main issues arise with this ap-
proach: how to choose the path and how to handle regions of inconsistent phase. The path
was chosen by the quality-guided method described below. Phase inconsistencies are called
residues. These were handled by introducing branch cuts. These methods are described in
detail in the following sections.
6.2.4 Quality-Guided Phase Unwrapping
In path-following methods, it is necessary to determine a path that visits each pixel of
interest. The phase is unwrapped along this path. Quality-guided algorithms [100] use a
quality map, which indicates the reliability of the wrapped phase map at each pixel, to guide
the integration path. A starting point was arbitrarily chosen and the next pixel visited was
the neighbor with the highest quality. This process was repeated until all pixels of interest
were processed. Similar to DENSE [82] algorithms, quality in this research was defined as
the reciprocal of the phase difference variance [112], which is given by
zm,n =
nullnull
(?xi,j null?xm,n)2 +
nullnull
(?yi,j null?ym,n)2
9 (6.2)
71
where each sum is over a 3null3 window, centered at the pixel (m,n). The terms ?xi,j and ?yi,j
are the wrapped phase derivatives, and ?xm,n and ?ym,n are the average of the derivatives in
the 3null3 window.
6.2.5 Residues and Branch Cuts
A fundamental assumption in phase unwrapping algorithms is that the unwrapped phase
is smooth and that any large changes in phase between pixels are due to phase wrapping.
Noise and artifacts, however, can cause inconsistencies in the phase map where the phase
difference between pixels is large and cannot be corrected by adding an integer multiple of
2pi. Tag fading, myocardial motion, and the presence of boundaries can also cause phase
inconsistencies. Examples of phase inconsistencies are shown in Figure 6.2. +45 degree tag
lines disappear before they reach the endocardium at 12 o?clock and 5 o?clock in Figure 6.2c.
Also the tag-like phase wrap near the epicardium at 4 o?clock in Figure 6.2c is caused by the
epicardial boundary. Figure 6.2d shows a ?double wrap? artifact where two tag lines merge
(7 o?clock).
Phase inconsistencies can cause any phase-unwrapping algorithm to fail. Figure 6.3a
shows an example of a phase image with inconsistencies. Figure 6.3b shows the unwrapped
phase image obtained with the technique described above. The phase inconsistencies result
in discontinuities in the unwrapped phase. In this section we present a method for detecting
and correcting for phase inconsistencies using residues and branch cuts.
Phase inconsistencies were detected using the principle that the wrapped phase dif-
ference around any closed path is equal to zero if the wrapped phase is consistent [89].
Consequently, these phase inconsistencies were localized by summing the wrapped phase dif-
ferences around a 2null2 neighborhood of each pixel in the myocardium. It can be shown [89]
that this sum is either 0,+2pi, or null2pi. If the sum is zero, then no phase inconsistency exists.
Otherwise, either a positive residue (+2pi) or a negative residue (null2pi) is present at that
pixel. Examples of residues can be seen in Figures 6.2e and 6.2f.
72
Figure 6.2: (a-b) Tagged short-axis images acquired near end-systole in two different human
subjects; (c-d) HARP images obtained from a and b. Arrows indicate points of phase
discontinuity. Cyan contours delineate the myocardium; (e-f) negative (black) and positive
(white) residues in the HARP images in (c-d); (g-h) the corresponding HARP phase image
with branch cuts drawn (yellow lines).
73
Figure 6.3: Phase Unwrapping example in a short-axis slice of a normal human left ventricle:
(a) wrapped HARP phase from Figure 2c; (b) unwrapped HARP image without branch cuts;
(c) HARP image with branch cuts shown in red (d) unwrapped HARP image with branch
cuts.
74
In this research, phase inconsistencies were corrected using the residue compensation
method [89], where residues were removed by connecting residues with branch cuts. Branch
cuts are lines in the image where an unwrapping path is not allowed to cross. In our
implementation, pixels that touch the branch cut line were removed from the myocardium
segmentation. Note that the number of pixels removed due to branch cuts is negligible
compared to the total number of pixels used in the 3D deformation model fit described
below.
Branch cuts either connect a positive to a negative residue or, if only a region of the
image is unwrapped, a residue is connected to a point outside the myocardium segmentation.
The problem, of course, is which residues to connect. It is common to see positive and
negative residues next to each other in the image (see Figure 6.2e). These residues are
known as dipoles and can be automatically connected by a branch cut.
Connecting non-dipole residues (called monopoles), however, is a difficult problem. In
cases where there are an equal number of positive and negative residues, automated methods
have been proposed that, for example, minimize the average branch cut length [89,106].
Minimizing branch-cut length, however, does not always yield the best unwrapped phase.
As a result, the user must place the branch cuts interactively. This was done quickly with a
graphical user interface (GUI).
Figures 6.2g and 6.2h show branch cut placement examples. The phase wraps in the
second row of Figure 6.2 correspond to the +45 direction tag lines. If the phase inconsistency
was caused by a prematurely ending tag line (11 o?clock and 5 o?clock in Figure 6.2c) or tag
line merging (7 o?clock in Figure 6.2d), the branch cuts were placed along the missing or
merged tag line. If a phase wrap occurred where there was no tag line (4 o?clock in Figure
6.2c) the branch cut was placed to remove the spurious phase wrap from the myocardium
segmentation.
Once branch cuts were placed, the phase unwrapping technique described above resulted
in a smooth unwrapped phase map as shown in Figure 6.3d. Note that after branch cuts
75
Figure 6.4: (a) HARP image; (b) negative (black) and positive (white) residues in the HARP
image in (a); (c) the corresponding demodulated HARP phase image; (d) residues in the
demodulated HARP image in (c).
are placed, the unwrapped phase is independent of the path taken during the unwrapping
process.
6.2.6 Demodulated HARP Phase
Instead of unwrapping the HARP phase image, one could unwrap the demodulated
HARP phase image. The demodulated HARP phase image is computed by shifting the
filtered HARP peak to the center of k-space before computing the inverse Fourier transform.
Figures 6.4a and 6.4c show the HARP phase and demodulated HARP phase images of the
tagged image in Figure 6.2a. The demodulated HARP phase is similar to the phase obtained
with DENSE imaging [113]. At first, it might appear that the demodulated phase is easier
76
to unwrap than the modulated HARP phase because there are fewer phase wraps. Phase
wraps, however, are usually not a problem for phase unwrapping algorithms. The challenge
in phase unwrapping is how to handle phase inconsistencies. Since phase inconsistencies are
mainly caused by tag lines disappearing, tag lines merging, myocardial boundaries and other
inherent properties of the tissue being imaged, they appear in both the HARP phase and
demodulated HARP phase images. For example, Figures 6.4b and 6.4d show the residues for
both the HARP phase and demodulated HARP phase images. The residues in both images
are almost identical. While either the HARP phase or demodulated HARP phase image
could be unwrapped, the HARP phase image was used in this research because the phase
wraps correspond to tag lines and provide visual cues that make it easier for the user to
place branch cuts as seen in Figure 6.3c.
6.2.7 Inter-frame Phase Consistency
Each image was unwrapped independently, and the starting point for phase unwrapping
may correspond to a different material point in each image. As a result, the unwrapped
phases in two adjacent time frames may differ by an integer multiple of 2pi. This difference
was corrected byadding the integer multiple of 2pito all unwrappedphases in the currenttime
frame. In this research, the multiple was chosen to minimize theL1 norm of all displacements
between the current and previous time frames. Note that this correction assumes that the
average deformation of the heart between consecutive frames is less than one-half of the tag
spacing. Inter-frame deformation of more than one-half tag spacing in a localized region or
regions will not affect the unwrapped phase as long as the average deformation is less than
one-half tag spacing.
77
6.2.8 Motion and Strain Estimation
Calculating 1-D Displacement Measurements
Once the phases were unwrapped, 1-D displacements were measured from each pixel
in the myocardium similar to those measured from tag-line tracking data [14]. Each 1-
D displacement measurement is the displacement of a material point back to its reference
position in the direction normal to the tag plane. As in [13,14,114], we define the reference
time to be the time the tag pattern was applied. Note that this time is several milliseconds
before the first image is acquired.
To compute the 1-D displacements, an estimate of the reference phase plane must be
computed. The reference phase plane was estimated by fitting a plane to the unwrapped
phase plane in the first time frame in a given slice. This computation is analogous to the
one performed in [13,14,114] to determine the tag line center and spacing. There can be
some motion between the time the tags are applied and the acquisition of the first image,
but this is usually small. Since tag lines are parallel when they are applied, the reference
phase plane was specified by the 1-D linear equation
?r = ms+b (6.3)
where s is a parameter that sweeps out a line perpendicular to the tag line and m and b are
parameters of the reference plane fit to the first time frame unwrapped phase.
To compute the 1-D displacement measurements, let sd correspond to a point in a
deformed image with phase ?d(sd). From Equation (6.3), the point in the reference image
with the same phase is sr = (?d(sd) nullb)/m. The displacement measurement is then d =
sd nullsr. Figure 6.5 shows an example of displacement estimates measured from a short-axis
slice of a normal human volunteer at end-systole. Note that due to through-plane motion
a 1-D displacement does not establish a one-to-one correspondence between points at the
78
Figure 6.5: 1-D displacement vectors superimposed for (a) the +45 and (b) -45 tag directions.
The arrows point to the position of the points at the time of tag line imposition.
reference and deformed times. Rather, it is a 1-D constraint on the 3-D displacement of a
material point back to its position in the reference time [13,14,114].
Segmentation
Segmentation of the myocardium is the first step in most cardiac image analysis proce-
dures [44,50,82]. Both the phase unwrapping and displacement measurement steps require
that the LV wall be segmented in each slice and time of interest. In this algorithm, LV
contours were drawn semi-automatically by the user at end-diastole (ED) and end-systole
(ES). Non-rigid registration [115] was then used to propagate ED and ES contours to all the
time frames. ED contours are used in almost all tagged MRI analysis algorithms to define
a material-point mesh or other set of points of interest. ES contours were included because
they improve the ability of the propagated contours to distinguish between papillary muscles
and the LV wall.
Strain Calculation
The 1-D displacement measurements and a material-point mesh automatically con-
structed from the ED contours were used to compute 3-D LV deformation and strain in
79
each time frame. The deformation and strain were reconstructed using the Affine Prolate
Spheroidal B-Spline (APSB) method in [116], which fits a B-spline deformation model de-
fined in prolate-spheroidal coordinates to the displacement measurements. First, a backward
fit was performed in each time frame that maps each point back to its undeformed position.
A series of forward fits were then performed to compute the deformed position of each point
in the undeformed mesh of the initial time frame. The trajectory of each point in the
material-point mesh through the cardiac cycle is thus obtained. 3-D strain was computed
by spatially differentiating the deformation model.
Validation Experiments
The mSUP technique was validated on the cohort of 40 human subjects described above.
In each subject, 3D strain was computed in all imaged time frames using the mSUP tech-
nique. 3D strain was computed at end-systole (ES) with the feature-based (FB) technique
in [116]. Both the FB and mSUP methods use the APSB deformation model described
in [116]. The tags were tracked using the FB technique and then manually corrected by an
expert. 2D strain in each time frame was computed using HARP analysis [117].
To assess the accuracy of 1D displacement measurements obtained from unwrapped
phase images, 1D displacement measurements were computed from tag lines tracked in the
FB method and compared to 1D measurements computed from the unwrapped phase images
at the tracked tag line points. The accuracy of 3D strains at end-systole was assessed by
comparing mSUP strains and FB strains with a paired t-test. The accuracy of strains versus
time was assessed by comparing peak strains and strain rates computed using mSUP and
HARP with paired t-tests. In all statistical comparisons, a P value of 0.05 or less was
considered statistically significant.
80
6.3 Results
6.3.1 Computation Time
All studies were processed using the mSUP technique described above, which was im-
plemented in MATLAB (The Mathworks Inc, Natick, MA). Unwrapping all images in a
study took less than a minute on a 2.6GHz Core2 Duo processor with 4 GB of memory.
Approximately 30 minutes per study of user interaction were required to resolve residues.
Strain reconstruction took approximately 15 minutes per study. The total time required to
analyze 20 time frames of a typical study was 45 minutes.
6.3.2 Quantitative Comparisons
Comparison of Displacement Estimates with Feature-Based Technique
Table 6.1 shows statistics of the differences between 1D displacement measurements
computed from unwrapped phase images and tracked tag lines at end-systole. The mSUP
and tag line methods showed excellent agreement on long-axis images and short-axis images
in the proximal and middle thirds of the LV (relative to the base). The mean differences are
a few hundredths of a pixel, which means there is no bias toward under- or overestimation
of displacement. The difference standard deviation is around 1/3 of a pixel, which is close
to the tag line tracking accuracy [114,118]. In three slices at or next to the apex, (out of a
total of 429 processed) the maximum difference was large (null 8 pixels) due to partial volume
effects that severely diminished the contrast-to-noise ratio (CNR). These slices were removed
from the differences in Table 6.1, but left in for the comparison of strains in subsequent
comparisons.
Comparison of End-Systolic Strains with Feature-Based Technique
Figure 6.6 shows the maps of 3D circumferential strain (Ecc) in a representative normal
human volunteer computed using mSUP and the FB method. The same material-point mesh
81
Table 6.1: Comparison of displacement measurements computed from mSUP at tracked tag
lines at end-systole. Proximal, middle, and distal refer to the LV base. Difference = mSUP
- tag line. All measurements are in pixels. RMSE = root mean squared error.
Slices Mean Diff. Std. Dev. Max. Diff.
All slices 0.0012 0.3963
4-Chamber -0.0263 0.3364 2.0479
2-Chamber 0.0052 0.3349 2.7517
Short-axis Proximal 3rd -0.017 0.3572 2.5694
Short-axis - Middle 3rd -0.0217 0.3299 2.4348
Short-axis Distal 3rd 0.0787 0.5625 2.9898
and deformation model were used to reconstruct strain in both methods. Note that strain
is computed in the entire LV except the apex, which is typically not included in this type of
analysis. The strain maps are similar.
Table 6.2 shows statistics of the difference between the mSUP and FB methods. In all
strains in Table 6.2, the standard deviation of the difference is less than 5% of the average
of the two methods. The correlation between the methods is strongest in circumferential
strain. Only two long-axis images were acquired, so the correlation in longitudinal strain
(Ell) and torsion is lower but still good. In fact, the mSUP method may provide more
accurate measurements of Ell, torsion and other strains than the FB method, but further
investigation is needed to support this claim.
Table 6.2: Comparison of 3D, end-systolic strains computed from strain from unwrapped
phase (mSUP) and feature-based (FB) methods. Differences = mSUP-FB. Differences =
Mean null Standard Error. ? = Correlation coefficient. CV = coefficient of variation. Ecc =
Circumferential shortening strain. Ell = Longitudinal shortening strain. Emin = Maximal
shortening strain.
Strain Differences p ? p CV
Ecc 0.0009 null 0.0017 0.9105 0.9613 <0.001 2.04%
Ell -0.0022 null 0.0019 0.7752 0.9351 <0.001 2.88%
Emin -0.0028 null 0.0019 0.6717 0.9201 <0.001 1.59%
Torsion -0.3498 null 0.2251 0.5757 0.8697 <0.001 3.66%
82
Figure 6.6: Maps of circumferential strain (Ecc) using feature-based (left) and unwrapped
phase (right) methods. Strains are mapped from blue = -30% to yellow = +30%. The green
balls denote the basal mid septum.
Comparison with HARP Strains
Table 6.3 shows statistics of the difference between the 3D strains calculated from mSUP
and the 2D strains from the HARP method [117]. 2D LV torsion using HARP was measured
using the procedure described in [51]. One slice each at the basal and apical levels was
used to compute torsion. The basal slice chosen was the slice closest to and below the
mitral valve at the end-systolic (ES) phase. The apical slice chosen was the slice closest
to the apex where the blood pool could be identified. A mesh consisting of 3 concentric
rings and 24 circumferential points was defined in each slice at end-diastole (ED) from semi-
automatically-drawn contours. Each mesh was tracked through all time frames using the
improved harmonic phase (HARP) method for motion tracking [52].
The peak strain and systolic strain rate of the measured circumferential strain show
a high correlation. Correlation in longitudinal strain is lower because the mSUP method
uses 3D deformation model fitting to compute longitudinal strain from the motion of all tag
points, whereas HARP computes strain from two long-axis images. Also, the tag line CNR
83
decreases through the cycle due to T1 decay of the tag pattern, so early diastolic correlations
are lower than correlations for peak strain and systolic strain rate. It was seen that HARP
in general underestimated the strains as compared to the mSUP method. This could be
attributed to the improved accuracy of tracking using the mSUP method as a result of
removal of inconsistencies in the HARP image using branch cuts.
A few degrees of difference was seen in comparison of peak torsion. The 3D method
corrects for the through-plane motion, but the 2D method tracks points in a given slice
through time. In the 2D method, the tracked points in the basal slice start out near mid-
ventricle and become more basal through systole. Since the base typically rotates more than
the mid-ventricle, though-plane motion can distort the 2D rotation and torsion curves.
Table 6.3: Comparison of strains and torsion using the strain from unwrapped phase (mSUP)
and HARP methods. Difference = mSUP - HARP. Differences = Mean null Standard Error.
? = Correlation coefficient. CV = coefficient of variation. Ecc = Circumferential shortening
strain. Ell = Longitudinal shortening strain.
Strain Differences p ? p CV
Ecc
Peak Strain -0.0223 null 0.0005 0.004 0.8361 <0.001 3.91%
Systolic Rate -0.1004 null 0.0024 0.011 0.8387 <0.001 3.49%
Early-Diastolic Rate -0.0494 null 0.0068 0.502 0.6615 <0.001 8.10%
Ell
Peak Strain -0.0217 null 0.0005 <0.001 0.7023 <0.001 4.93%
Systolic Rate 0.0657 null 0.0042 0.132 0.624 <0.001 6.41%
Early-Diastolic Rate -0.2091 null 0.0111 0.052 0.5739 <0.001 12.77%
Torsion
Peak Torsion -0.0036 null 0.1092 0.998 0.6567 <0.001 7.00%
Systolic Rate -9.4941 null 0.5079 0.069 0.649 <0.001 7.053
Early-Diastolic Rate 3.9996 null 0.8812 0.533 0.2578 0.108 10.41%
6.3.3 mSUP Strains in Normals and Patients
Figure 6.7 shows plots of torsion versus time averaged over 10 normal human volunteers,
10 patients with myocardial infarction, 10 patients with hypertension and 10 patients with
mitral regurgitation. The hypertensive patients showed an increase in torsion, which has
84
Figure 6.7: Torsion angle (rotation of the apex relative to the base) versus % systolic inter-
val averaged over 10 normal human volunteers, 10 patients with myocardial infarction, 10
patients with hypertension and 10 patients with mitral regurgitation. Error bars represent
null one standard error.
been reported in the clinical literature [39]. The untwisting in early-diastole was slower
for patients with mitral regurgitation as compared to normal human volunteers, which is
consistent with previously reported results [119].
6.4 Discussion
The mSUP method was presented for measuring 3D strain from tagged MRI. Strains
computed using mSUP demonstrate excellent agreement with a previously published feature
based method. Some user interaction is also required to resolve phase-unwrapping residues
(especially later in the cycle), but the mSUP method can compute strains in the entire LV
instead of just the middle third.
85
The time required for the mSUP-measured 3D strains in the entire LV over 20 time
frames is approximately 45 minutes per study. The same analysis using the feature-based
method would require approximately 3 hours.
While the tagged imaging protocol used in this research is commercially available and is
widely used clinically, advanced tagged imaging techniques such as Complementary SPAMM,
can yield a higher tag CNR throughout the cardiac cycle. Higher CNR images will result in
fewer phase inconsistencies, and less user interaction will be required to correct them.
The use of residues to automatically detect phase inconsistencies in wrapped phase
images and the use of branch cuts to correct them can also be used in unwrapping the phase
in DENSE imaging.
Spottiswoode et al. [82] used a spatio-temporal method of phase unwrapping to maintain
inter-frame phase consistency in DENSE images, where a group of pixels are selected in the
first time frame and these pixels are sequentially unwrapped in each successive time frame.
This technique works well in DENSE images because there are only a few phase wraps across
the myocardium. The mSUP method, however, is based on HARP images, which make it
easier to place branch cuts. HARP images have 6-10 phase wraps across the myocardium.
Temporal phase unwrapping in HARP images is sensitive to large (greater than one-half tag
spacing) interframe displacements. Therefore, we used a method based on the phase of the
entire myocardium over time to maintain phase consistency. This allows for local tissue to
have displacements greater than one-half tag-spacing between time frames as long as the
average deformation of tissue between time frames is less than one-half tag-spacing.
The mSUP technique reconstructs full 3-D strain maps from tagged MR images through
the cardiac cycle and takes into account the through-plane motion of the heart. Slice-
following DENSE [83], zHARP [86] and HARP-SENC [85] can account for through-plane
motion, but these techniques require two breath-holds per slice. In contrast, tagged images
can be acquired with multiple slices per breath-hold, which allows the entire LV to be imaged
in significantly less time.
86
In conclusion, the strain from unwrapped phase (mSUP) method can compute 3-D
strains in the LV through the cardiac cycle in a reasonable amount of time and user inter-
action compared to other 3D analysis methods.
87
Chapter 7
Comparison Of 2D And 3D Torsion Measured From Tagged MRI
7.1 Introduction
Cardiac left ventricular (LV) torsion is the difference between rotation of the apex and
base (see Chapter 1). LV torsion can provide important information about both systolic
and diastolic myocardial function [120]. In MRI data, torsion is typically measured by
tracking an annular mesh in an apical and a basal slice with two-dimensional (2D) HARP
analysis [51,52]. 2D-techniques may not take into account the through-plane myocardial
motion that may significantly affect the torsion measurements. In Chapter 6, a method for
computing 3D torsion was presented [90]. The purpose of this research is to quantitatively
compare 2D and 3D LV rotation and torsion measurements computed from tagged MRI.
7.2 Methods
A group of 40 subjects were imaged consisting of 10 normal volunteers (NRM), 10
diabetics with myocardial infarction (DMI), 9 patients with mitral regurgitation (MRR),
and 10 patients with hypertension (HTN). All participants underwent MRI on a 1.5T MRI
scanner (GE, Milwaukee, WI) optimized for cardiac application. Both cine and tagged
images were acquired in standard views (2 and 4 chamber long axis and short axis) with
a fast gradient-echo cine sequence with the following parameters: FOV = 300 mm, image
matrix = 224x256, flip angle = 45, TE = 1.82ms, TR = 5.2ms, number of cardiac phases =
20, slice thickness = 10 mm. A 2D spatial modulation of magnetization tagging preparation
was done with a tag spacing of 7 pixels.
88
2D LV torsion was measured using the procedure described in [51]. One slice each at
the basal and apical levels was used to compute torsion. The basal slice chosen was the
slice closest to and below the mitral valve at the end-systolic (ES) phase. The apical slice
chosen was the slice closest to the apex where the blood pool could be identified. A mesh
consisting of 3 concentric rings and 24 circumferential points was defined in each slice at end-
diastole (ED) from semi-automatically-drawn contours. Each mesh was tracked through all
time frames using the improved harmonic phase (HARP) method for motion tracking [52].
Drawing contours and mesh tracking took approximately 5min per study.
3D LV torsion was measured using the method described in chapter 6 [90]. Myocardial
contours were semi-automatically drawn at ED and end-systolic (ES) time frames for all
slices and then propagated to all time frames [115]. 3D LV deformation in each time frame
was computed by fitting the 1D displacement measurements derived from unwrapping the
HARP phase in each image to the affine prolate spheroidal B-spline method in [116]. Drawing
contours and other processing took approximately 30min per study.
Displacement of the mitral annulus toward the apex during systole was measured by
tracking a user-specified point on the mitral annulus through each frame of a 4-chamber view
using non-rigid registration [115].
7.3 Results and Discussion
Figure 7.1 shows rotation and torsion versus time curves averaged over the studies in
each of the four patient groups. Differences between 2D and 3D curves can be attributed to
two factors. First, a few degrees of rotation occur between tag pattern application and the
first images. The 3D method corrects for this initial deformation, but the 2D method does
not because it is based on tracking points identified in the first image. This factor shifts the
3D curves by a constant angle through time. The second factor is through-plane motion ?
primarily motion of the base toward the apex during systole and back again during diastole.
The 3D method corrects for this motion, but the 2D method tracks points in a given slice
89
Figure 7.1: Plots of average base and apex rotation, and torsion obtained from 2D (black)
and 3D (red) methods through 70% of the cardiac cycle. Error bars represent null one standard
error. All angles are in degrees
Figure 7.2: Peak mitral annulus displacement.
90
Table 7.1: Comparison of differences between 3D and 2D rotation and torsion (3D-2D).
nullp< 0.05 relative to 3D.
Peak (deg) Sys Rate (deg/s) Early Dia Rate (deg/s)
Base Rotation
NRM 4.3 null 0.8* 7.4 null 4.8 -47.1 null 18.8*
HTN 3.1 null 1.3 7.4 null 6.3 -23.1 null 11.5
DMI 3.1 null 1.1 6.1 null 6 -1.1 null 14.5
MRR 3.9 null 1.4 11.3 null 7.8 -34.6 null 7.4*
Apex Rotation
NRM 1.9 null 0.4* 6.6 null 2.4 -6.8 null 2.8
HTN 1.3 null 0.8 7.6 null 2.5 -14.4 null 3.8*
DMI -0.9 null 1 -3.1 null 4.4 2.7 null 3.4
MRR -2.1 null 0.9 -7.9 null 5.1 6.2 null 4.8
Torsion
NRM 1.8 null 1 -3.6 null 6.2 -36.3 null 18.2
HTN 1.2 null 1.8 -0.9 null 6.8 -8 null 15.2
DMI 4 null 0.9 9 null 6.3 -7.6 null 13.5
MRR 5.8 null 1.6* 21 null 5.5* -43.5 null 6.1*
through time. In the 2D method, the tracked points in the basal slice start out near mid-
ventricle and become more basal through systole. Since the base typically rotates more than
the mid-ventricle, through-plane motion can distort the rotation and torsion curves.
The effect of through-plane motion can be seen by comparing the measurements of
mitral annulus displacement (MAD) in Figure 7.2 with the curves in Figure 7.1. The DMI
patients have the lowest amount of MAD, and the 3D rotation and torsion curves are shifted
versions of the 2D curves. Normals and MR patients have thr largest amount of MAD, and
2D curves are more distorted relative to 3D. Table 7.1 shows the difference between 2D and
3D peak torsion/rotation and peak torsion/rotation rates. The significant differences are
primarily in the normal and MR patient groups because of the relatively high through-plane
motion.
91
7.4 Conclusion
While 2D and 3D rotation and torsion measurements are similar in subjects with low
base-to-apex motion such as patients with MI, the 3D method measures larger and probably
more accurate rotation and torsion in subjects with normal or elevated base-to-apex motion.
This may allow improved assessment of cardiac mechanics in patients with a variety of
myocardial diseases.
92
Chapter 8
3D+t Bi-Ventricular Strain using Phase Unwrapped HARP
8.1 Introduction
Accurate measurement of cardiac ventricular mechanical function is important for di-
agnosing and managing patients with heart disease and assessing the efficacy of therapies
over time. The right ventricle (RV) plays a vital role in cardiac function and is adversely
affected by many cardiac and pulmonary diseases. Compared to the left ventricle (LV), RV
functional analysis is relatively difficult due to its lack of geometric symmetry and relatively
thin wall.
Accurate assessment of right ventricular (RV) function is particularly important in pa-
tients with pulmonary hypertension (PHTN) and congestive heart failure. Also, in PHTN,
relatively higher systolic pressure in the RV can cause excursion of the inter-ventricular sep-
tum into the LV cavity (see Figure 8.1a). As a result, the LV cavity can also lose its geometric
symmetry, which is assumed in many 3D LV deformation models [7,29,65,69,121].
In Chapter 6, the strain from unwrapped phase (mSUP) method was presented for
computing three-dimensional plus time (3D+t) strain in the LV [90]. The mSUP method
is based on unwrapping 2D harmonic phase (HARP) images [49] from tagged MRI and
a prolate-spheroidal deformation model [7]. The mSUP strains compared well with strains
from 3D feature-based methods [7] and (wrapped) 2D HARP [49], but the prolate-spheroidal
deformation model does not fit the RV geometry and is difficult to extend to a biventricular
model. The discrete model free (DMF) model proposed by Denney and McVeigh [114],
however, is geometry independent and can be easily adapted to biventricular morphology.
Also, the unwrapped HARP produces displacement measurements at each pixel, which is
useful when estimating the deformation of the thin-walled RV.
93
Figure 8.1: (a) Basal short-axis slice of a PHTN patient near end-systole showing RV excur-
sion. (b) Same image with the radially-oriented long axis prescriptions shown.
In this chapter, myocardial displacements from unwrapped HARP in both the RV and
LV were used with DMF deformation to reconstruct 3D+t biventricular strain in each im-
aged time frame through systole and early diastole. We refer to this method as the biven-
tricular strain unwrapped phase (BiSUP) technique. The BiSUP technique produces 3D+t
biventricular displacements and strains with minimal manual intervention and with high spa-
tial resolution. Results are presented in normal volunteers, patients with LV hypertension
(HTN), patients with PHTN, and diabetic patients with myocardial infarction (DMI).
8.2 Materials and Methods
8.2.1 Human Subjects
The BiSUP algorithm was validated on a cohort of 30 human subjects (10 normal volun-
teers, 7 patients with pulmonary hypertension, 8 patients with hypertension, and 5 diabetics
with myocardial infarction). All human studies were approved by the Institutional Review
Board of both institutions (Auburn University and University of Alabama at Birmingham)
and informed consent was obtained from all participants.
94
8.2.2 Materials
All participants underwent MRI on a 1.5T MRI scanner (GE Healthcare, Milwaukee,
WI) optimized for cardiac application. Standard short-axis cardiac views were obtained,
yielding on average 12 short-axis slices. Six equally-spaced, radially-oriented long-axis views
(see Figure 8.1b) were obtained, which resulted in 2-3 long-axis slices through the RV free
wall, with a fast gradient-echo cine sequence with the following parameters: FOV = 300 mm,
image matrix = 224x256, flip angle = 45, TE = 1.82ms, TR = 5.2ms, number of cardiac
phases = 20, slice thickness = 8 mm. A 2D fast gradient-recalled spatial modulation of
magnetization (FGR-SPAMM) tagging preparation was done with a tag spacing of 7 pixels.
8.2.3 Segmentation
Myocardial contours modeled using active contours were semi-automatically drawn at
end-diastole (ED) and end-systole (ES) time frames for all slices (see Figure 8.2). In the
short-axis images, two closed contours (LV endo-contour and LV epi-contour) were drawn
for the left ventricle, and two open contours (RV endo-contour and RV epi-contour) were
drawn for the right ventricle. In the long-axis images, all four were open contours. The RV
contours connected to the LV epi-contour at the insertion points. These contours were then
propagated to all time frames [115].
8.2.4 Displacement Estimation
3D displacement was measured from the unwrapped HARP phase method described
in [90] and in Chapter 6. Phase unwrapping involves adding integer multiples of 2pi to the
wrapped phase so that the unwrapped phase is continuous. A quality-guided path-following
method, also used by Spottiswoode et al. [82] to unwrap wrapped DENSE images, is used
to unwrap the wrapped HARP phase images.
95
Figure 8.2: A short-axis slice in a normal volunteer at end-systole showing user-input RV
and LV contours for the segmentation of the myocardium.
The unwrapped-phase method described uses residues [89] to detect inconsistencies in
the wrapped HARP phase. Residues were computed by locally unwrapping a 2 null 2 neigh-
borhood of each pixel in the region of interest. Phase inconsistencies were corrected using
the residue compensation method [89], where residues were removed by connecting residues
with branch cuts. Branch cuts are lines in the image where an unwrapping path is not
allowed to cross. Branch cuts are specified by the user with a graphical user interface (GUI).
Branch-cut placement in the RV is an extension of the method used for the LV in [90]. Since
the LV and RV contours are connected at the RV insertion points, the two ventricles are
unwrapped as a single entity.
Figure 8.3 shows the effect of branch cuts and the resulting unwrapped phase maps in
a short-axis and a long-axis image. Once branch cuts were placed, the phase unwrapping
technique described above resulted in a smooth unwrapped phase map as shown in Figure 8.3
(8.3g and 8.3h). Note that after branch cuts are placed, the unwrapped phase is independent
of the path taken during the unwrapping process.
Once the phases were unwrapped, 1-D displacements were measured from each pixel in
the myocardium similar to those measured from tag-line tracking data [14].
96
Figure 8.3: Phase Unwrapping example in a short-axis slice (left column) and a long-axis slice
(right column) of a patient with pulmonary hypertension: (a-b) Tagged short-axis images
acquired near end-systole; (c-d) HARP images obtained from (a) and (b) with branch cuts
drawn(red lines). White contours delineate the myocardium. ; (e-f) HARP negative (black)
and positive (white) residues in the HARP images in (c-d); (g-h) unwrapped HARP image
with branch cuts.
97
8.2.5 Strain Computation
The DMF technique [114] was used to reconstruct 3D biventricular strains at each
imaged time frame from 1D displacement measurements. The DMF technique makes no
assumptions about the geometry of the object being analyzed, so both LV and RV strains
can be reconstructed as a single object. In addition, because the DMF technique uses finite
difference analysis, it does not require a specific coordinate system for the reconstruction [13].
The method proceeds as follows to reconstruct displacement and strain from tag lines
and user-defined myocardial contours. For each imaged time frame, the DMF first constructs
a 3D segmentation of the LV wall on a regularly spaced discrete grid of points from the user-
defined endocardial and epicardial contours. The grid spacing is equal to the pixel size in the
MR image. A displacement vector is estimated for each grid point in the LV segmentation
that maps the grid point back to its end-diastolic position using finite-difference analysis. A
smoothness constraint minimizes the spatial variation of the displacement gradient across the
myocardium. The displacement gradient is then computed using a finite-difference derivative
approximation averaged over a 3null3null3 point neighborhood. Finally, the Lagrangian strain
tensor is computed from the displacement gradient.
The strains were computed in the circumferential direction (Ecc), longitudinal direction
(Ell), and in the direction of maximum principal shortening strain (Emin) in the LV. In the
RV, strains were computed in the tangential direction along the short axis of the LV (EttRV),
the longitudinal direction (EllRV) and the direction of maximum principal shortening strain
(EminRV). In addition, 3D torsion in degrees was computed for the LV.
8.2.6 Validation Experiments
The BiSUP technique was validated on the cohort of 30 human subjects described
above. In each subject, 3D strain was computed in all imaged time frames using the BiSUP
technique. To validate 3D strain reconstruction, 3D strain was computed at end-systole (ES)
withthefeature-based(FB)techniquein[114], whichalsousestheDMFstrainreconstruction
98
technique. The tags were tracked using the FB technique and then manually corrected by
an expert.
To assess the accuracy of 1D displacement measurements obtained from unwrapped
phase images, 1D displacement measurements were computed from tag lines tracked in the
FB method and compared to 1D measurements computed from the unwrapped phase images
at the tracked tag line points. The accuracy of 3D strains at end-systole was assessed by
comparing BiSUP strains and FB strains with paired t-tests, measuring the coefficient of
variation and Bland-Altman plots.
To validate the evolution of strains over time, BiSUP strains were compared to 2D strain
in each time frame computed using HARP analysis [49]. 2D LV torsion using HARP was
measured using the procedure described in [51]. One slice each at the basal and apical levels
was used to compute torsion. The basal slice chosen was the slice closest to and below the
mitral valve at the end-systolic (ES) phase. The apical slice chosen was the slice closest
to the apex where the blood pool could be identified. A mesh consisting of 3 concentric
rings and 24 circumferential points was defined in each slice at end-diastole (ED) from semi-
automatically-drawn contours. Each mesh was tracked through all time frames using the
improved harmonic phase (HARP) method for motion tracking [52]. The accuracy of strains
versus time was assessed by comparing peak strains and strain rates computed using BiSUP
and HARP with paired t-tests. In all statistical comparisons, a p-value of 0.05 or less was
considered statistically significant.
8.3 Results
8.3.1 Computation time
All studies were processed using the BiSUP technique described above, which was im-
plemented in MATLAB (The Mathworks Inc, Natick, MA). Unwrapping all images in a
study took less than a minute on a 2.6GHz Core2 Duo processor with 4 GB of memory.
Approximately 60 minutes per study of user interaction were required to resolve residues.
99
Strain reconstruction took approximately 30 minutes per study of automatic computation.
The total time required to analyze 20 time frames of a typical study was 90 minutes.
8.3.2 Comparison of Displacement Estimates with FB Technique
Table 8.1 shows statistics of the differences between 1D displacement measurements
computed from unwrapped phase images and tracked tag lines at end-systole. The BiSUP
and tag line methods showed excellent agreement on long-axis images and short-axis images
in the proximal and middle thirds of the LV (relative to the base). The mean differences are
a few hundredths of pixel, which means there is no bias toward under or over estimation of
displacement. The difference standard deviation is around 1/2 of a pixel, which is close to
the tag line tracking accuracy [114,118]. A small number of differences (< 1%) were in the
range of 2-5 pixels. Most of these were near the apex in both the LV and RV where partial
volume effects diminished the contrast-to-noise ratio (CNR). However, these measurements
were small in number and isolated and had little effect on the resulting strain map due to
the smoothness constraints used in the reconstruction.
Table 8.1: Comparison of displacement measurements computed from BiSUP at tracked tag
lines at end-systole. Proximal, middle, and distal refer to the LV base. Difference = BiSUP
- Tagline
Slices Mean Diff. Std. Dev. Diff >2 px
All slices 0.03 0.51 0.73%
Long-axis 0.03 0.44 0.75%
Short-axis Proximal 3rd 0.00 0.44 0.38%
Short-axis - Middle 3rd 0.06 0.55 0.79%
Short-axis Distal 3rd 0.10 0.69 1.46%
8.3.3 Comparison of End-Systolic Strains with Feature-Based Technique
Figure 8.4 shows the maps of 3D maximal shortening strain (Emin) in a representative
normal human volunteer computed using BiSUP method and the FB method. The same
DMF reconstruction method was used to reconstruct strain in both methods. Note that
100
strain is computed in the entire LV and RV, except the apex which is typically not included
in this type of analysis. The strains maps are similar. Table 8.2 shows statistics of the
difference between the BiSUP and FB methods. In all strains in Table 8.2, the standard
deviation of the difference is less than 4% of the average of the two methods. The correlation
between the methods is strongest in longitudinal strain. Six long-axis images were acquired,
allowing more comparisons so the correlation in longitudinal strain (Ell) and torsion is higher,
but the other strains measured showed a high correlation as well. Comparisons between the
two sets of strains are displayed in Bland-Altman and correlation plots shown in Figure 8.5
and Figure 8.6. The differences between the computed strains are small with no significant
difference found in the LV as shown in Table 8.2. The tangential strain in the RV was found
to be significantly different between the BiSUP and FB methods, but this could be because
of the increased number of short-axis measurements with BiSUP, which would have a larger
effect on the relatively sparsely sampled RV.
Table 8.2: Comparison of 3D, end-systolic strains computed from strain from unwrapped
phase (BiSUP) and feature-based (FB) methods. Differences = Mean Standard Error. ? =
Correlation coefficient. CV = coefficient of variation. Ecc = LV Circumferential shortening
strain. Ell = LV Longitudinal shortening strain. Emin = LV Maximal shortening strain.
EttRV = RV Tangential shortening strain. EllRV = RV Longitudinal shortening strain.
EminRV = RV Maximal shortening strain.
Strain Differences p ? p CV
Ecc 0.0095 null 0.0019 0.18 0.93 <0.001 2.01%
Ell -0.0046 null 0.0018 0.63 0.97 <0.001 1.62%
Emin -0.0048 null 0.0022 0.56 0.93 <0.001 1.59%
Torsion -1.8019 null 0.2461 0.10 0.97 <0.001 3.01%
EttRV 0.0212 null 0.0026 0.03 0.93 <0.001 3.11%
EllRV 0.0151 null 0.0025 0.15 0.95 <0.001 2.31%
EminRV 0.0154 null 0.0021 0.13 0.96 <0.001 1.51%
8.3.4 Comparison with HARP Strains
Table 8.3 shows statistics of the difference between the 3D strains calculated from BiSUP
and the 2D strains from the HARP method [49]. Significant differences were found between
101
Figure 8.4: Maps of maximum shortening strain (Emin) using feature-based (left) and BiSUP
(right) methods for a normal volunteer (NL), and patients with pulmonary hypertension
(PHTN), LV hypertension (HTN) and diabetics with myocardial infarction (DMI). Strains
are mapped from blue = -25% to yellow = +25%. The anterior LV wall is at the bottom of
each strain map.
102
Figure 8.5: Bland-Altman (left) and correlation (right) plots for comparing BiSUP and FB
strain measurements in the LV. In the Bland-Altman plots, lines above and below the zero-
difference line represent 2 SD of the difference. In the correlation plots, the dashed line
represents perfect agreement.
103
Figure 8.6: Bland-Altman (left) and correlation (right) plots for comparing BiSUP and FB
strain measurements in the RV. In the Bland-Altman plots, lines above and below the zero-
difference line represent 2 SD of the difference. In the correlation plots, the dashed line
represents perfect agreement.
104
the two techniques in maximum shortening strains (Emin), peak systolic torsion rate and
EminRV rate, but the coefficients of variation for these parameters were all less than 10%.
Peak strains showed the strongest correlations between the two techniques. Weakest correla-
tions and largest coefficients of variation were found in early diastolic rates due to tag fading
in late diastole. Tag fading causes phase inconsistencies, which are corrected in BiSUP with
branch cuts. HARP does not provide a mechanism for correcting phase inconsistencies.
Figure 8.7 shows strain and torsion curves for both BiSUP and HARP averaged over
the studies for each of the four patient groups showing good agreement between BiSUP and
HARP in the temporal evolution of strain.
Table 8.3: Comparison of strains and torsion using the strain from unwrapped phase (BiSUP)
and HARP methods. Differences = Mean Standard Error. ? = Correlation coefficient. CV
= coefficient of variation. Ecc = LV Circumferential shortening strain. Emin = LV Maximal
shortening strain. EttRV = RV Tangential shortening strain. EminRV = RV Maximal
shortening strain.
Strain Differences p ? p CV
Ecc
Peak Strain -0.0066 null 0.0009 0.39 0.61 <0.001 5.35%
Systolic Rate -0.0659 null 0.0045 0.10 0.61 <0.001 5.06%
Early-Diastolic Rate -0.0604 null 0.0127 0.52 0.53 <0.001 13.27%
Emin
Peak Strain -0.0224 null 0.0006 0.01 0.83 <0.001 2.64%
Systolic Rate -0.2604 null 0.0057 0.00 0.66 <0.001 5.99%
Early-Diastolic Rate 0.1328 null 0.0063 0.01 0.58 <0.001 8.81%
Torsion
Peak Torsion -1.5122 null 0.1322 0.23 0.70 <0.001 7.35%
Systolic Rate -16.1142 null 0.7465 0.01 0.58 <0.001 8.49
Early-Diastolic Rate 13.0258 null 0.8720 0.07 0.55 0.0018 9.48%
EttRV
Peak Strain 0.0011 null 0.0010 0.90 0.61 <0.001 6.86%
Systolic Rate 0.0934 null 0.0105 0.16 0.24 0.2031 10.81%
Early-Diastolic Rate -0.3137 null 0.0248 0.07 0.39 0.0315 20.08%
EminRV
Peak Strain -0.0119 null 0.0008 0.21 0.80 <0.001 3.13%
Systolic Rate -0.1659 null 0.0076 0.01 0.60 <0.001 7.47%
Early-Diastolic Rate 0.0150 null 0.0094 0.83 0.48 0.0069 11.44%
105
Figure 8.7: Plots of average mid-ventricular strain and torsion obtained from 2D HARP
(black) and 3D BiSUP (red) methods for normal volunteers (NL) and, patients with pul-
monary hypertension (PHTN), LV hypertension (HTN) and diabetics with myocardial in-
farction (DMI). X-axis represents in percentage systolic interval. Error bars represent one
standard error.
106
Figure 8.8 shows maximal shortening strain for a normal volunteer and patients with
pulmonary hypertension, left ventricular hypertension and a diabetic with myocardial in-
farction at different phases of the cardiac cycle. It is clear that the shortening in the RV is
depressed in patients with pulmonary hypertension relative to the normal volunteers (also
see Figure 8.7). Also, note the excursion of the RV into the LV at the septum.
8.4 Discussion
The application of tagged MRI to the RV has various difficulties; the thin wall of
the normal RV limits the number of tag lines, which may hinder accurate quantification
[122]. Several myocardial tagging studies have been performed on the RV to investigate and
characterize mechanical deformation. Klein et al. [123] used measures of percent segmental
shortening of the RV free wall. Zahi et al. [122] designed a specific breath hold imaging
sequence with 1-D tag lines for RV tagging acquisition and thus were able to characterize the
RV regional deformation. Young et al. [124] used a finite-element model for 3D reconstruction
of the in-plane deformations of the RV free wall as a surface. Haber et al. [29] also used
a finite-element model to obtain 3D strains in the RV. Using this model, they were able
to quantify components of the in-plane strain tensor in the RV free wall. Tustison et al.
[74] computed biventricular strains using volumetric deformable models with a non-uniform
rational B-splines (NURBS) basis. All the 3D methods [29,74,124] feature tag tracking to
obtain displacements.
In this chapter, an unwrapped phase technique to compute 3D biventricular strain in
about three-fourths of the cardiac cycle using tagged cardiac MR images has been described
and validated. The BiSUP technique consists of an unwrapped phase technique to compute
the 1D displacements followed by a discrete model-free reconstruction method to obtain 3D
strain measurements.
We demonstrated on human subjects of different pathologies that the circumferential,
longitudinal and maximal shortening strains measured using BiSUP are similar to those
107
Figure 8.8: Maps of maximum shortening strain (Emin) using the BiSUP method for a
normal volunteer (NL), and patients with pulmonary hypertension (PHTN), LV hypertension
(HTN) and diabetics with myocardial infarction (DMI) from early-diastole through mid-
diastole of a cardiac cycle. Strains are mapped from blue = -25% to yellow = +25%. The
anterior LV wall is at the bottom of each strain map.
108
obtained from a feature-based method that uses tag tracking to obtain the displacement
measurements. The time required for the BiSUP measured strains in the entire LV and
RV combined over 20 time frames is less than 1.5 hours per study with less than an hour
of manual intervention. The same analysis using the feature-based method would require
approximately 4 hours with significantly higher manual intervention. Compared to methods
based on tracking tag lines [29,122,124], displacement measurements from the BiSUP method
are dense, potentially resulting in more accurate estimates of strain.
The strains from the BiSUP method were also compared to the strains from HARP.
BiSUP-HARP agreement was not as strong as BiSUP-FB agreement because the BiSUP
computes 3D strain, while HARP computes 2D strain. Also, the tag line CNR decreases
through the cycle due to T1 decay of the tag pattern, so early diastolic correlations are
lower than correlations for peak strain and systolic strain rate. It was seen that HARP in
general underestimated the maximal shortening strains as compared to the BiSUP method.
This could be attributed to the improved accuracy of tracking using the BiSUP method as
a result of removal of inconsistencies in the HARP image using branch cuts. Also, HARP
analysis does not correct for motion through the stationary basal and apical image planes,
whereas BiSUP reconstructs 3D deformation and strain in each time frame, which corrects
for through-plane motion. With HARP, torsion is computed with HARP tracking, and large
deformations between time frames can cause errors. The unwrapped phase technique used in
BiSUP is more robust than HARP to large interframe deformations [90]. With unwrapped
phase, points can be tracked through an image sequence as long as the average inter-frame
deformation over the entire myocardium is less than one-half tag spacing.
The DMF method has advantages that are utilized for the reconstruction of strain in
the RV and in the LV in patients with pulmonary hypertension. The DMF method does
not need a specific coordinate system to do the LV reconstruction and uses less smoothing
than other reconstruction methods [13]. Therefore, the lack of a consistent geometry does
109
not affect the performance of the DMF method as opposed to finite-element methods and
prolate-spheroidal methods.
In conclusion, the BiSUP method can compute accurate 3D biventricular strains through
the cardiac cycle in a reasonable amount of time and user interaction compared to other 3D
analysis methods.
110
Chapter 9
Computer-Assisted Strain from Unwrapped Phase
9.1 Introduction
The unwrapped-phase (mSUP) [90] method previously described in Chapter 6 uses
residues [89] to detect inconsistencies in the wrapped HARP phase. Residues were computed
by locally unwrapping a 2 null 2 neighborhood of each pixel in the region of interest. Phase
inconsistencies were corrected by a technique called branch cuts, which removed certain
pixels from the region of interest. Branch cuts are specified by the user with a graphical user
interface (GUI). With unwrapped phase, points can be tracked through an image sequence
as long as the deformation between time frames is less than half a tag spacing. Compared
to other feature-based methods, displacement measurements from unwrapped HARP phase
images are denser than those obtained from tracking tag lines or tag intersections, potentially
resulting in more accurate estimates of strain. It can compute strains in the entire LV and is
also more robust to interframe motion compared to methods based on tracking the wrapped
phase. This method still requires user interaction to place the branch cuts.
In this chapter, a computer-assisted branch-cut placement in short-axis HARP images
is presented to compute the strain from unwrapped phase (caSUP). This method is based
on branch-cut placement using simulated annealing, which has been proposed for general-
purpose phase unwrapping [125]. The method in [125], however, minimizes the branch cut
length, which is a reasonable criterion in single-image unwrapping. In cardiac imaging, we
wish to unwrap a sequence of images. We propose a new energy function for simulated
annealing-based phase unwrapping that is based on temporal continuity of the unwrapped
phase. An energy term is proposed that enables the application of the phase unwrapping
111
method to a sequence of images. The caSUP method described here requires user interac-
tion other than segmentation of the myocardium, to correct for any detected errors in the
unwrapping. It saves about one-third of the time involved in manual branch-cut placement.
To compute 3D LV strain, the caSUP method is used to unwrap short-axis images and the
manual branch-cut method [90] is used to unwrap long-axis images.
This chapter is organized as follows. The caSUP method and strain reconstruction is
presented in Section 2. Experimental validation results on 40 human studies and conclusions
are presented in Sections 3 and 4.
9.2 Materials and Methods
9.2.1 Human Subjects
A cohort of 40 studies consisting of 10 normal volunteers, 10 diabetics with myocardial
infarction, 10 patients with mitral regurgitation, and 10 patients with hypertension was used
to validate the caSUP method. All human studies were approved by the Institutional Review
Board of both institutions (Auburn University and University of Alabama at Birmingham)
and informed consent was obtained from all participants.
9.2.2 Data Acquisition
All participants underwent MRI on a 1.5T MRI scanner (GE, Milwaukee, WI) optimized
for cardiac application. Tagged images were acquired in standard views (2-chamber and 4-
chamber long-axis and short-axis) with a fast gradient-echo cine sequence with the following
parameters: FOV = 300 mm, image matrix = 224x256, flip angle = 45, TE = 1.82ms, TR
= 5.2ms, number of cardiac phases = 20, slice thickness = 8 mm. A 2D fast gradient-
recalled spatial modulation of magnetization (FGR-SPAMM) tagging preparation was done
with a tag spacing of 7 pixels. This protocol resulted in 2 long-axis slices and on average 12
short-axis slices per study.
112
9.2.3 Phase Unwrapping
In the unwrapped-phase (mSUP) method previously described in [90], branch cuts are
specified by the user with a graphical user interface (GUI). In this chapter we explore the
possibility of an automated branch-cut placement method. Ideally, given a set of positive
and negative residues, the surest way of finding the correct branch-cut configuration is to
explore every possible combination (exhaustive search). The configuration that minimizes a
defined energy function is chosen as the best fit.
A new energy function is proposed for phase unwrapping that minimizes the difference
in unwrapped phase between successive time frames. The energy function Et is defined as,
Et =
null
null?t null?tnull1null, (9.1)
where ?t is the unwrapped phase at time frame t. In the phase unwrapping literature,
algorithms seek to minimize the aggregate branch-cut length, which is a reasonable criterion
in single-image unwrapping. In cardiac imaging, we wish to unwrap a sequence of images.
The proposed energy function allows us to maintain temporal continuity of the unwrapped
phase, enabling its application to unwrap a sequence of images. The unwrapped phase ?
is obtained by unwrapping the HARP phase (wrapped) using flood-fill quality-guided phase
unwrapping [89,90] after the application of the present branch-cut configuration.
The exhaustive search method is adequate when the number of residues is low, but as the
number of resiudes increases, exploration of all combinations is impractical as the number of
possible configurations becomes prohibitively large (see Figure 9.1). In such cases, it is more
economical to use a method such as the simulated-annealing-based branch-cut (SABrCut)
placement method described below.
113
Number of positive residues
Nu
mbe
r of nega
tiv
e r
esidues
 
 
0 1 2 3 4 50
1
2
3
4
5
2
4
6
8
10
12
14
Figure 9.1: A logarithmic (log2) plot showing the number of possible branch cut configura-
tions given the number of residues in the image.
9.2.4 Simulated Annealing
This method is based on branch-cut placement using simulated annealing, which has
been proposed for general-purpose phase unwrapping [125]. The method in [125], however,
minimizes the branch-cut length, while we minimize energy defined in Equation (9.1).
Simulated annealing based branch cut placement [125] consists of:
1. An initial random configuration of branch cuts, which can be addressed as a proposed
solution to the branch-cut placement problem. In this research, the initial branch-cut
configuration is chosen by nearest-neighbor method, but any automatic algorithm can be
used. In the nearest-neighbor method, each residue is connected to the closest residue of
opposite polarity or the border, whichever is the nearest.
2. An energy function E that describes the quality of a proposed solution. In this research,
we minimize energy defined in Equation (9.1).
3. A generation mechanism that applies small changes to produce a new branch-cut config-
uration. For each iteration k, a small change to the branch-cut configuration is applied.
114
A negative and a positive residue are randomly selected. Existing branch cuts involving
these residues are deleted, and new branch cuts are applied according to the scheme de-
scribed in Figure 9.2. This scheme is similar to the one used in [125] with the exception
of the endocardial/epicardial configuration, which is unique to our algorithm. If there
is more than one possible new branch-cut configuration, one is chosen at random. The
energy Ek+1t of the new configuration is then calculated. The energy difference is
?Ekt = Ek+1t nullEkt (9.2)
4. An acceptance criterion to decide whether the new proposed solution is accepted. The
new configuration is always accepted if ?E < 0. Accepting a new configuration only if
?E < 0, though, may mean that the algorithm is trapped in a local minimum. Therefore,
increasing values of E are also accepted at a Boltzmann probability. The acceptance
probability p, therefore is
p =
nullnull
nullnull
nullnullnull
1 , if ?E null 0
enull?E/T , if ?E > 0
(9.3)
5. A control parameter T to control the acceptance probability of proposed solutions and
a cooling schedule for T. The Boltzmann probability is controlled by the temperature
T. For high values of T, even high increases in energy could be accepted, while for low
values of T the acceptance probability of these moves is close to zero. When a number
of transitions are applied to the system at a constant temperature, the system reaches a
thermal equilibrium. The series of configurations resulting from these transitions is called
a homogeneous Markov chain. The SABrCut algorithm consists of several Markov chains
that start at high values of T. Reducing the temperature T slowly reduces the internal
energy of the system until it reaches a minimum, which is identified as the solution of our
115
optimization problem. The values of T are determined by a cooling schedule
Tk+1 = T
k
1 + [Tk ln(1 +?)/(3?k)], (9.4)
with Equation (9.4) controlling the cooling speed and ? the standard deviation of the
energy during the final Markov chain. The cooling speeds are kept low to allow the
algorithm to reach the global minimum. The values Tinitial = 10000 and ? = 0.25 were
chosen empirically.
6. A stopping criterion to decide whether the algorithm has found a final solution. When
there is no change in the accepted energy for a number of consecutive configuration
changes, thealgorithmisterminatedandthefinalconfigurationisacceptedasthesolution.
9.2.5 caSUP
Figure 9.3 shows the time taken using the exhaustive search and simulated annealing
methods to find a suitable solution for a given number of residues in the image. In this
research, we utilized the fact that the exhaustive search to determine the best branch-cut
configuration was faster than the SABrCut when the number of residues is small, but the
SABrCut method is faster for a greater number of residues. Therefore, the exhaustive search
method was used when the number of residues is less than 6, and the SABrCut method was
used otherwise.
As explained in [90], branch cuts are required because of prematurely ending tag lines,
tag-line merging or spurious phase wraps. The branch cuts can be seen as either the contin-
uation or deletion of a tag line. Therefore, branch cuts are drawn parallel to the tag lines
(i.e. along the tag line direction) when connecting to the myocardial borders.
A discontinuity map, similar to the one used in [82], was used to detect images with
incorrect unwrapping. A discontinuity map measures the gradient of the unwrapped phase,
and if the unwrapping is correct, the gradient magnitude of the unwrapped phase does not
116
OR
OR
OR
OR
OR
possible changesinitial configuration
border
positive 
source
negative 
source
pos. / neg. 
source
epicardium
endocardium
Figure 9.2: Generation mechanism to apply small changes to the existing set of branch
cuts. Two sources of opposite sign are selected randomly. Depending on the initial situation,
branch cuts are deleted and new branch cuts are set. A border connection is a connection
to either the epicardial or the endocardial boundary, each of which is chosen with equal
probability.
117
2 4 6 8 10
0
1
2
3
4
Number of residues
Logarithm of time in seconds  
 
Exhaustive
Simulated Annealing
Figure 9.3: The logarithm (log10) of the time taken in seconds using the exhaustive search
and simulated annealing methods to find a suitable solution depending on the number of
residues.
exceed pi. When discontinuity is detected in the unwrapping of a slice, user interaction (as
in [90]) is used to manually place the branch cuts for the slice. The entire procedure involved
is called computer-assisted strain from unwrapped phase (caSUP).
9.2.6 Motion and Strain Estimation
Each image is unwrapped independently, and the starting point for phase unwrapping
may be different in each image. As a result, the unwrapped phases in two adjacent time
frames may differ by an integer multiple of 2pi. This difference can be corrected by adding
an integer multiple of 2pi to all unwrapped phases in the current time frame. The multiple
is chosen to minimize the L1 norm of all displacements between the current and previous
frames. Note that this correction assumes that the L1 norm of the actual unwrapped phase
difference between consecutive frames is less than pi, which corresponds to one-half of the
tag spacing. Also, interframe deformation of more than one-half tag spacing in a localized
region or regions will not affect the unwrapped phase as long as the average deformation is
less than one-half tag spacing.
118
The caSUP method described above was used to unwrap all short-axis slices. Long-axis
slices were unwrapped using the manually placed branch-cut method described in [90]. Once
the phases were unwrapped, 1-D displacements were measured from each pixel in the my-
ocardium using the technique described in [90]. Each 1-D displacement measurement is the
displacement of the point back to its undeformed position along a line perpendicular to the
tag plane (see [90] for details). Both the phase alignment and displacement measurement
steps required that the LV wall be segmented in each slice and time of interest. This segmen-
tation was done using the automated dual-contour propagation method described in [126].
The 1-D displacement measurements and a material-point mesh automatically constructed
from the end-diastolic (ED) contours were used to compute 3-D LV deformation and strain
in each time frame. The deformation and strain are reconstructed using the Affine Prolate
Spheroidal B-Spline (APSB) method described in [90].
9.3 Experiments
9.3.1 Computation Time
All studies were processed using the caSUP technique described above, which was im-
plemented in MATLAB (The Mathworks Inc, Natick, MA). Unwrapping all images in a
study took less than a minute on a 2.6GHz Core2 Duo processor with 4 GB of memory.
Approximately 30 minutes per study of automated processing followed by 15 minutes of user
interaction were required to resolve residues. Strain reconstruction took approximately 5
minutes per study. The total time required to analyze 20 time frames of a typical study was
50 minutes. This saves 15-20 minutes of user interaction for the user.
9.3.2 Validation
The caSUP technique was validated on the cohort of 40 human subjects. In each subject,
3D strain was computed in all imaged time frames using the caSUP and mSUP techniques.
Both the caSUP and mSUP methods use the APSB deformation model described in [116].
119
The mSUP method was earlier validated with and compared to HARP strains [50] and
3D strain from a feature-based method [45,116]. The accuracy of 3D strains was assessed
by comparing caSUP strains and mSUP strains at end-systole with a paired t-test. In all
statistical comparisons, a P value of 0.05 or less was considered statistically significant.
A total of 349 short-axis slices were processed. In 338 of these slices, the caSUP algo-
rithm placed branch cuts correctly, and no user interaction was required. In eleven slices,
errors were detected automatically using the discontinuity map and manually corrected.
The errors occur predominantly at the apical slices, caused by partial volume effects and low
contrast-to-noise ratio (CNR).
Table 9.1 shows statistics of the difference in averaged mid-ventricular strains between
the unwrapped phase methods caSUP and mSUP. In all strains in Table 9.1, the caSUP
strains are highly correlated with mSUP strains, and the difference standard deviation is
less than 3% of the average of the two methods. The mean difference between caSUP
and mSUP methods is low across all measured quantities, showing no bias toward under- or
overestimation of strain. The longitudinal strain and torsion show a slightly lower correlation
compared to the circumferential and maximal shortening strains, because there are 8-12
short-axis slices per study but only two long-axis images. Also, the torsion measurements
show a lower correlation without corrections as most of the errors occur at the apex.
Figure 9.5 shows maps of 3D circumferential strain (Ecc) computed from the caSUP
and the mSUP methods for a representative from all the patient groups. Note that the same
material-point mesh and deformation model was used to reconstruct strain in all methods.
The strain maps are mostly similar. Differences, if any, are noted at the apex where the
increased through-plane motion and rotation affects the ability of the computer-assisted
branch cut placement method to obtain the correct branch-cut configuration.
120
ba
g h
fe
dc
?0.2 ?0.15 ?0.1 ?0.05?0.1
?0.05
0
0.05
(Ecc mSUP +Ecc caSUP )/2
(Ecc
mSUP
?Ecc
caSUP
)/2
?0.15 ?0.1 ?0.05
?0.15
?0.1
?0.05
Ecc mSUP
Ecc
caSUP
?0.2 ?0.15 ?0.1 ?0.05
?0.1
?0.05
0
0.05
0.1
(Ell mSUP +Ell caSUP )/2
(Ell
mSUP
?Ell
caSUP
)/2
?0.15 ?0.1 ?0.05
?0.2
?0.15
?0.1
?0.05
Ell mSUP
Ell
caSUP
?0.3 ?0.2 ?0.1
?0.1
?0.05
0
0.05
0.1
(Emin mSUP +Emin caSUP )/2
(Emin
mSUP
?Emin
caSUP
)/2
?0.2 ?0.15 ?0.1
?0.2
?0.15
?0.1
Emin mSUP
Emin
caSUP
10 20 30
?10
?5
0
5
10
(Torsion mSUP +Torsion caSUP )/2
(To
rsi
on
mSUP
?Tor
sion
caSUP
)/2
5 10 15 20 25
5
10
15
20
25
Torsion mSUP
Tor
sion
caSUP
Figure 9.4: Bland-Altman (left) and correlation (right) plots for comparing caSUP and
mSUP strain measurements. In the Bland-Altman plots, lines above and below the zero-
difference line represent 2 SD of the difference. In the correlation plots, the dashed line
represents perfect correlation.
121
NL
HTN
MRR
mSUP caSUP
DMI
Figure 9.5: Maps of circumferential (Ecc) strain using computer-assisted strain from un-
wrapped phase (left) and manual strain from unwrapped phase for normal (NL), hyperten-
sive (HTN), mitral regurgitant (MRR) and myocardial infarct (DMI) hearts at end-systole.
Strains are mapped from blue = -25% to yellow = 25%. The mid-septum is marked by the
green landmark.
122
Early-Systole Mid-Systole End-Systole Early-Diastole
NL
HTN
MRR
DMI
Figure 9.6: Maps of maximal shortening (Emin) strain using computer-assisted strain from
unwrapped phase for normal (NL), hypertensive (HTN), mitral regurgitant (MRR) and
myocardial infarct (DMI) hearts through the cycle. Strains are mapped from blue = -25%
to yellow = 25%. The mid-septum is marked by the green landmark.
123
Table 9.1: Comparison of strains and torsion computed using caSUP and mSUP meth-
ods. Differences are caSUP-mSUP. Differences = Mean null Standard Error. ?= Correlation
coefficient. For all correlation coefficients, p< 0.001. CV = coefficient of variation.
Strain Differences p ? CV
After corrections
Ecc -0.002 null 0.001 0.79 0.99 0.65%
Ell 0.001 null 0.002 0.87 0.96 2.78%
Emin 0.001 null 0.001 0.94 0.99 0.66%
Torsion 0.321 null 0.227 0.72 0.94 2.33%
Without corrections
Ecc -0.003 null 0.001 0.71 0.99 0.77%
Ell -0.001 null 0.001 0.97 0.97 1.72%
Emin -0.003 null 0.002 0.72 0.95 1.38%
Torsion 0.354 null 0.296 0.71 0.91 3.03%
9.3.3 caSUP in Normals and Patients
Figure 9.6 shows maps of 3D maximal shortening strain (Emin) computed from the
caSUP method for a representative from all the patient groups. The figure demonstrates
the ability of the caSUP method to produce strains in various pathologies and with different
cardiac morphologies.
9.4 Discussion
A computer-assisted strain from unwrapped phase (caSUP) procedure was presented
for measuring 3D deformation from tagged MRI. This method used automated branch cut
placement to resolve residues, before unwrapping HARP images. Manual intervention is used
to modify the branch-cut configuration when errors are discovered in the proposed branch
cuts. Displacement and 3D strain can then be calculated from the unwrapped phase images.
40 studies from different groups were used to validate the algorithm. caSUP strains demon-
strated excellent agreement with the previously presented manual strain from unwrapped
phase (mSUP).
The time required for the caSUP-measured 3D strains in the entire LV over 20 time
frames is approximately 45 minutes per study with about 10 minutes of user interaction
124
required. The same analysis would require around 30-40 minutes of user interaction using
the manual strain from unwrapped phase (mSUP). The feature-based method [45,116] would
require approximately 3 hours of user interaction. Most of the user interaction required is
usually later in the cycle (mid to late diastole) when the tag CNR reduces rapidly. Advanced
tagged imaging techniques such as Complementary SPAMM, can yield a higher tag CNR
throughout the cardiac cycle. Higher CNR images will result in fewer phase inconsistencies,
and therefore, automated resolving of branch cuts is faster and lesser user interaction will be
required to correct them. Consequently, if strains are to be computed only through systole
as is so often the case, required user interaction would drastically reduce.
The HARP method [50] is sensitive to large (greater than one-half tag spacing) in-
terframe displacements. Therefore, we used a method based on the phase of the entire
myocardium over time to maintain phase consistency. This allows for local tissue to have
displacements greater than one-half tag-spacing between time frames as long as the average
deformation of tissue between time frames is less than one-half tag-spacing. Moreover, it
allows the users to check for errors and correct them, unlike HARP.
In conclusion, the computer assisted strain from unwrapped phase (caSUP) method can
compute parameters such as 3-D strain, strain rate and torsion in the LV through the cardiac
cycle in a reasonable amount of time and with minimal user interaction compared to other
3D analysis methods.
125
Chapter 10
Torsion Hysteresis
10.1 Introduction
Left ventricular (LV) diastolic dysfunction is characterized by abnormal myocardial me-
chanical properties that include impaired diastolic distensibility, impaired LV filling and slow
or delayed myocardial relaxation [127]. Up to 50% of patients with heart failure have pre-
dominant diastolic dysfunction in the presence of preserved LV ejection fraction (EF) [128].
Measurement of LV pressure-volume relations is perhaps the reference standard for diagnos-
ing diastolic dysfunction. This is based on the premise that the LV diastolic pressure rises
excessively during ventricular filling due to impaired LV distensibility and/or compliance.
Quantifying the exponential fall of pressure during isovolumic relaxation phase by the time
constant of relaxation, ?, is considered a reliable, relatively load independent measure of
diastolic relaxation of the LV.
LV torsion/twist is an important mechanical property of the myocardium that promotes
more efficient cardiac function. It represents the myocardial rotation gradient from the base
to apex along a longitudinal axis [120]. Torsion is due to the criss-crossing helical orientation
of the myofibrils from the subendocardium to the subepicardium. Contraction of myocardial
bundles and their interaction with extracellular matrix during systole results in storage of
torsional potential energy. Torsional recoil during isovolumic relaxation and early diastole
releases the potential energy stored in the deformed matrix during systole [129,130]. Previous
work with cardiac magnetic resonance imaging (MRI) in an animal model [131] and recently
inhumans[132] hasshownthatLVearly untwistrate correlatesclosely withthetime constant
of LV relaxation.
126
The term ?hysteresis? is used in mechanical engineering, to describe history-dependent
stress-strain relationships in materials undergoing cyclic loading. In these applications, the
area between stress-strain loops is related to energy loss due to friction [133]. Research
groups have used hysteresis in stress-strain loops to study viscoelastic losses in isolated
myocytes [134,135], chick embryo heart [136?138] and other animal models [139].
In the current investigation, we apply the hysteresis concept to torsion-volume loops
obtained in normal subjects and in patients with hypertension (HTN) and mitral regurgi-
tation (MR) patients with preserved LVEF and no evidence of coronary artery disease. We
hypothesized that the area enveloped by these curves (torsional hysteresis) represents loss
of energy from work against viscoelastic properties of myocardium and hence will represent
global diastolic function.
10.2 Magnetic Resonance Imaging
Magnetic resonance imaging was performed on a 1.5T MRI scanner (Signa GE, Milwau-
kee, Wisconsin) optimized for cardiac application. Electrocardiographically gated breath-
holdsteady-statefreeprecision wasusedtoobtainstandard(2-, 3-, and4-chambershort-axis)
views using the following parameters: slice thickness of the imaging planes 8 mm, field of
view 44 null 44, scan matrix 256 null 128, flip angle 45, repetition/echo times 3.8/1.6 ms. The
temporal resolution of the cine MRI was < 50 ms. Tagged magnetic resonance images were
acquired on the same scanner with repetition/echo times 8/44 ms, and tag spacing 7 mm.
Typical temporal resolution of the tagged MRI scan was 80 ms. Imaged population included
37 normal volunteers (NL), 78 patients with hypertension (HTN) and patients with mitral
valve disease (MR).
10.3 Geometric analysis
LV geometric parameters were measured from endocardial and epicardial contours man-
ually traced on cine images acquired near end-diastole and end-systole. These contours were
127
Figure 10.1: Volume time curve for the 3 groups. Data points are expressed as mean null
standard error. Normal is represented in black, hypertension is represented in red and mitral
regurgitation is represented in blue. Early filling rate is significantly lower in Hypertension
group compared to the other two groups.
propagated throughout the cardiac cycle using in-house software [115,126] to create LV
time-volume curves 10.1.
10.4 Torsional Deformation Analysis
Two-dimensional (2D) strain rates were measured using harmonic phase (HARP) anal-
ysis [49?51]. 2D torsion and torsion rates were measured by tracking a circular mesh of
points in a basal and apical slice. The mesh was identified in the first time based on user-
defined contours and tracked through the remaining imaged phases using improved HARP
tracking [50,52]. Because the tag lines faded with time due to T1 relaxation, torsion was
only computed through systole and the first 70% of diastole (see Figure 10.2).
10.5 Torsional Hysteresis
The torsion hysteresis area was computed with the following procedure for each subject.
First, each point on the torsion-versus-time curve was divided by peak torsion, and each point
on the LV volume-versus-time curve was divided by EDV to give normalized values. The
128
Figure 10.2: Torsion time curve for the 3 groups. Data points are expressed as mean null
standard error. Normal is represented in black, hypertension is represented in red and
mitral regurgitation is represented in blue. Peak systolic torsion is highest in Hypertension
group.
areas under the systolic and diastolic sections of the torsion-volume curve were computed
numerically using the trapezoidal rule. Both areas were computed over the volume interval
between the ESV and the volume at 70% diastole as shown in Figure 10.3. The normalized
torsion hysteresis area was computed by subtracting the diastolic area from the systolic area.
Similarly non-normalized torsion hysteresis was also computed.
10.6 Results and Discussion
The LV torsional hysteresis as measured by the net area difference was highest in HTN
compared to MR or NL (p < 0.001) (Figures 10.4 and 10.5).
In this study we demonstrate for the first time that under distinct physiological condi-
tions, one of eccentric remodeling (MR) and the other of concentric remodeling (HTN), the
area within the torsion-volume loop (torsional hysteresis) is increased in HTN compared to
NL or MR. We believe that the underlying mechanism for torsional hysteresis is likely due
129
Figure 10.3: Schematic diagram for calculating torsional hysteresis. Torsional hysteresis
represents the difference in area within the systolic and diastolic arms of the normalized
torsion volume (TV) curves.
to loss of stored torsional energy through inefficient myocardial relaxation. Diastolic un-
twist/torsion is characterized by early rapid recoil with little change in LV volume, followed
by more gradual untwisting when the bulk of diastolic filling occurs.
Hysteresis in the stress-strain relationship for a material is due to frictional losses during
material motion; these stress-strain relationships may be linear (stretching); shear (sliding),
or torsional (twisting), reflecting degrees of freedom in the motion of the body. Measurement
of the torsional stress (i.e., the torque) is difficult in the heart in vivo. Therefore, we have
used the LV volume as a surrogate for torque, since a twisting force acting on a shape with
a larger radius will result in increased torque, compared to a shape with a smaller radius.
In the case of the LV, apart from the intrinsic passive myocardial stiffness, the torsion-
volume loop may also be influenced by the rate of LV filling and the active ATP dependent
relaxation process occurring in the myofibrils. In HTN due to a stiffer ventricle, a greater
amount of work is expended for both active and passive relaxation and therefore greater
torsional hysteresis.
In our study we included MR as another group with a distinctly different remodeling
compared to HTN. This group demonstrated eccentric remodeling with a dilated LV cavity
130
Figure 10.4: a) Normalized b) non-normalized torsion volume curves for the 3 groups. Data
points are expressed as mean null standard error. Normal is represented in black, hypertension
is represented in red and mitral regurgitation is represented in blue. Torsional hysteresis as
measured by the area within the systolic and diastolic arm is greater in Hypertension group
compared to the other two groups.
131
Figure 10.5: a) Normalized and b) non-normalized torsional hysteresis measured as described
in Figure 10.3. Torsional hysteresis is significantly higher in Hypertension group compared
to Mitral Regurgitation or Normal control groups.
132
but LVEF > 60%. We found no difference in the torsional hysteresis in this group compared
to NL. Previous investigators have demonstrated that there is a loss or at the most no net
gain of collagen in MR [140]. The present study and previous studies have also demonstrated
relatively preserved diastolic function as measured by conventional parameters in MR and
preserved LVEF [141].
In conclusion, torsional hysteresis may be an important global parameter to assess both
systolic and diastolic function. To our knowledge, this is the first study to propose the
concept of torsional hysteresis as a measure of diastolic dysfunction, and its application may
provide greater insight into multiple forms of cardiac hypertrophy and heart failure.
133
Chapter 11
Conclusion and Future Work
In this dissertation, a method to compute displacement measurements from tagged MR
images is presented that depends on unwrapping the harmonic phase images. A branch-cut
placement method was used to remove inconsistencies and unwrap the images. Temporal
consistency is imposed while unwrapping the images. Initially, a GUI-based approach, called
manual Strain from Unwrapped Phase (mSUP) was designed where the user placed branch-
cuts. This was later replaced by a semi-automated method to place branch cuts in short-
axis images, called computer-assisted Strain from Unwrapped Phase (caSUP). Using caSUP
reduced the amount of manual intervention considerably.
The viability of applying the same methods in the right ventricle was also explored.
Computing strains in the RV has proven especially difficult, but the application of unwrapped
phase HARP to the RV (BiSUP) yielded accurate and quality 3D strain maps. BiSUP also
demonstrated the ability to compute strains over time in hearts of different pathologies and
morphologies.
An area of future work could be filtering techniques to obtain better phase images
(with fewer residues). This would reduce the inconsistencies in the image, and consequently
make the caSUP method faster while also requiring lesser manual intervention. Conversely,
techniques such as CSPAMM could be used to obtain higher-quality tagged images, that
would lead to fewer inconsistencies in the image and also reduce tag fading later in the
cycle. Exploration of faster and more accurate algorithms to place branch cuts is another
prospective area of research. It would especially be useful if this could be extended to 3D
(in either time or spatial domain) in an efficient manner.
134
Finally, extensive data has been collected for this research. 110 patients and volunteers
with 3D+t and 2D+t strains have been used. This could be useful in studying relationships
between 3D and 2D strains and strain rates and specifically studying the advantages of 3D
strain.
135
Bibliography
[1] ?http : //en.wikipedia.org/wiki/cardiology.?
[2] ?http : //www.invisionguide.com/heart/.?
[3] ?http://www.cvphysiology.com/index.html.?
[4] ?http : //www.columbiasurgery.org/pat/cardiac/anatomy.html.?
[5] ?http : //www.americanheart.org/.?
[6] J. L. Prince and J. M. Links, Medical Imaging Signals and Systems. Prentice-Hall,
2001.
[7] J. Li, ?Tag line tracking and cardiac motion modeling from tagged MRI,? Ph.D. dis-
sertation, Auburn University, 2006.
[8] C. Ozturk, J. A. Derbyshire, and E. R. McVeigh, ?Estimating motion from MRI data,?
Proceedings of the IEEE, vol. 91, no. 10, pp. 1627?1648, 2003.
[9] S. Masood, G.-Z. Yang, D. J. Pennell, and D. N. Firmin, ?Investigating intrinsic my-
ocardial mechanics: The role of MR tagging, velocity phase mapping and diffusion
imaging,? Journal of Magnetic Resonance Imaging, vol. 12, pp. 873?878, 2000.
[10] J. Duncan, P. Shi, T. Constble, and A. Sinusas, ?Physical and geometrical modeling
for image-based recovery of left ventricle deformation,? Progress in Biophysics and
Molecular Biology, vol. 69, no. 2-3, pp. 333?351, 1998.
[11] E. A. Zerhouni, D. M. Parish, W. J. Rogers, A. Yang, and E. P. Shapiro, ?Human
heart: tagging with MR imaging?a method for noninvasive assessment of myocardial
motion,? Radiology, vol. 169, no. 1, pp. 59?63, 1988.
[12] E. Atalar and E. R. McVeigh, ?Optimization of tag thickness for measuring position
with magnetic resonance imaging,? IEEE Transactions on Medical Imaging, vol. 13,
no. 1, pp. 152?160, 1994.
[13] J. Declerck, T. S. Denney, C. Ozturk, W. O?Dell, and E. R. McVeigh, ?Left ventricular
motion reconstruction from planar tagged MR images: a comparison,? Physics in
Medicine and Biology, vol. 45, no. 6, pp. 1611?1632, 2000.
[14] T. S. Denney and J. L. Prince, ?Reconstruction of 3-D left ventricular motion from
planartaggedcardiac MRimages: anestimationtheoreticapproach,? Medical Imaging,
IEEE Transactions on, vol. 14, no. 4, pp. 625?635, 1995.
136
[15] J. Bogaert and F. E. Rademakers, ?Regional nonuniformity of normal adult human
left ventricle,? Am J Physiol Heart Circ Physiol, vol. 280, no. 2, pp. H610?620, 2001.
[16] C. C. Moore, C. H. Lugo-Olivieri, E. R. McVeigh, and E. A. Zerhouni, ?Three-
dimensional systolic strain patterns in the normal human left ventricle: Character-
ization with tagged MR imaging,? Radiology, vol. 214, no. 2, pp. 453?466, 2000.
[17] L. Axel, ?Biomechanical dynamics of the heart with MRI,? Annual Review of Biomed-
ical Engineering, vol. 4, pp. 321?347, 2002.
[18] C. Petitjean, N. Rougon, and P. Cluzel, ?Assessment of myocardial function: A review
of quantification methods and results using tagged MRI,? Journal of Cardiovascular
Magnetic Resonance, vol. 7, no. 2, pp. 501?516, 2005.
[19] M. D. Cerqueira, N. J. Weissman, V. Dilsizian, A. K. Jacobs, S. Kaul, W. K.
Laskey, D. J. Pennell, J. A. Rumberger, T. Ryan, and M. S. Verani, ?Standardized
myocardial segmentation and nomenclature for tomographic imaging of the heart: A
statement for healthcare professionals from the cardiac imaging committee of the
council on clinical cardiology of the American Heart Association,? Circulation, vol.
105, no. 4, pp. 539?542, 2002. [Online]. Available: http://circ.ahajournals.org
[20] M. B. Buchalter, J. L. Weiss, W. J. Rogers, E. A. Zerhouni, M. L. Weisfeldt, R. Be-
yar, and E. P. Shapiro, ?Noninvasive quantification of left ventricular rotational de-
formation in normal humans using magnetic resonance imaging myocardial tagging,?
Circulation, vol. 81, no. 4, pp. 1236?1244, 1990.
[21] R. E. Henson, S. K. Song, J. S. Pastorek, J. J. H. Ackerman, and C. H. Lorenz, ?Left
ventricular torsion is equal in mice and humans,? Am J Physiol Heart Circ Physiol,
vol. 278, no. 4, pp. H1117?1123, 2000.
[22] M. A. Fogel, K. B. Gupta, P. M. Weinberg, and E. A. Hoffman, ?Regional wall mo-
tion and strain analysis across stages of Fontan reconstruction by magnetic resonance
tagging,? Am J Physiol Heart Circ Physiol, vol. 269, no. 3, pp. H1132?1152, 1995.
[23] E. Pettersen, T. Helle-Valle, T. Edvardsen, H. Lindberg, H.-J. Smith, B. Smevik, O. A.
Smiseth, and K. Andersen, ?Contraction pattern of the systemic right ventricle: Shift
from longitudinal to circumferential shortening and absent global ventricular torsion,?
Journal of the American College of Cardiology, vol. 49, no. 25, pp. 2450?2456, 2007.
[24] E. W. Remme, K. F. Augenstein, A. A. Young, and P. J. Hunter, ?Parameter distri-
bution models for estimation of population based left ventricular deformation using
sparse fiducial markers,? Medical Imaging, IEEE Transactions on, vol. 24, no. 3, pp.
381?388, 2005.
[25] A. A. Young, C. M. Kramer, V. A. Ferrari, L. Axel, and N. Reichek, ?Three dimensional
left ventricular deformation in hypertrophic cardiomyopathy,? Circulation, vol. 90,
no. 2, pp. 854?867, 1994.
137
[26] I. Russel, M. Gotte, J. Kuijer, and J. T. Marcus, ?Regional assessment of left ventricu-
lar torsion by CMR tagging,? Journal of Cardiovascular Magnetic Resonance, vol. 10,
no. 1, p. 26, 2008.
[27] L. Pan, J. Prince, J. Lima, and N. Osman, ?Fast tracking of cardiac motion using 3D-
HARP,? Biomedical Engineering, IEEE Transactions on, vol. 52, no. 8, pp. 1425?1435,
2005.
[28] C. C. Moore, E. R. McVeigh, and E. A. Zerhouni, ?Quantitative tagged magnetic
resonance imaging of the normal human left ventricle,? Topics in Magnetic Resonance
Imaging, vol. 11, no. 6, pp. 359?371, 2000.
[29] I. Haber, D. N. Metaxas, and L. Axel, ?Three-dimensional motion reconstruction and
analysis of the right ventricle using tagged MRI,? Medical Image Analysis, vol. 4, no. 4,
pp. 335?355, 2000.
[30] J. Garcia-Barnes, D. Gil, J. Barajas, F. Carreras, S. Pujadas, and P. Radeva, ?Charac-
terization of ventricular torsion in healthy subjects using gabor filters and a variational
framework,? in Computers in Cardiology, 2006, 2006, pp. 877?880.
[31] R. Mankad, C. J. McCreery, W. J. Rogers, R. J. Weichmann, E. B. Savage, N. Reichek,
and C. M. Kramer, ?Regional myocardial strain before and after mitral valve repair
for severe mitral regurgitation,? Journal of Cardiovascular Magnetic Resonance, vol. 3,
no. 3, pp. 257 ? 266, 2001.
[32] S.-J. Dong, A. P. Crawley, J. H. MacGregor, Y. F. Petrank, D. W. Bergman, I. Be-
lenkie, E. R. Smith, J. V. Tyberg, and R. Beyar, ?Regional left ventricular systolic
function in relation to the cavity geometry in patients with chronic right ventricular
pressure overload : A three-dimensional tagged magnetic resonance imaging study,?
Circulation, vol. 91, no. 9, pp. 2359?2370, 1995.
[33] C. F. Azevedo, L. C. Amado, D. L. Kraitchman, B. L. Gerber, N. F. Osman, C. E. Ro-
chitte, T. Edvardsen, and J. A. C. Lima, ?Persistent diastolic dysfunction despite com-
plete systolic functional recovery after reperfused acute myocardial infarction demon-
strated by tagged magnetic resonance imaging,? Eur Heart J, vol. 25, no. 16, pp.
1419?1427, 2004.
[34] A. A. Young, H. Imai, C. N. Chang, and L. Axel, ?Two-dimensional left ventricular
deformation during systole using magnetic resonance imaging with spatial modulation
of magnetization,? Circulation, vol. 89, no. 2, pp. 740 ? 752, 1994.
[35] D. Bree, J. R. Wollmuth, B. P.Cupps, M. D. Krock, A. Howells, J. Rogers, N. Moazami,
and M. K. Pasque, ?Low-dose dobutamine tissue-tagged magnetic resonance imaging
with 3-dimensional strain analysis allows assessment of myocardial viability in patients
with ischemic cardiomyopathy,? Circulation, vol. 114, no. 1suppl, pp. I?33?36, 2006.
[36] P. Moustakidis, B. P. Cupps, B. J. Pomerantz, R. P. Scheri, H. S. Maniar, A. M. Kates,
R. J. Gropler, M. K. Pasque, and T. M. Sundt, ?Noninvasive, quantitative assessment
138
of left ventricular function in ischemic cardiomyopathy,? Journal of Surgical Research,
vol. 116, no. 2, pp. 187?196, 2004.
[37] H. Kanzaki, S. Nakatani, N. Yamada, S.-i. Urayama, K. Miyatake, and M. Kitakaze,
?Impaired systolic torsion in dilated cardiomyopathy: Reversal of apical rotation at
mid-systole characterized with magnetic resonance tagging method,? Basic Research
in Cardiology, vol. 101, no. 6, pp. 465?470, 2006, 10.1007/s00395-006-0603-6.
[38] E. Nagel, M. Stuber, M. Lakatos, M. B. Scheidegger, P. Boesiger, and O. M. Hess, ?Car-
diac rotation and relaxation after anterolateral myocardial infarction,? Coron Artery
Dis, vol. 11, pp. 261 ? 267, 2000.
[39] M. Stuber, M. B. Scheidegger, S. E. Fischer, E. Nagel, F. Steinemann, O. M. Hess, and
P. Boesiger, ?Alterations in the local myocardial motion pattern in patients suffering
from pressure overload due to aortic stenosis,? Circulation, vol. 100, no. 4, pp. 361?8,
1999.
[40] H. C. Oxenham, A. A. Young, B. R. Cowan, T. L. Gentles, C. J. Occleshaw, C. G.
Fonseca, R. N. Doughty, and N. Sharpe, ?Age-related changes in myocardial relaxation
using three-dimensional tagged magnetic resonance imaging,? Journal of Cardiovascu-
lar Magnetic Resonance, vol. 5, no. 3, pp. 421 ? 430, 2003.
[41] J. L. Prince and E. R. McVeigh, ?Motion estimation from tagged MR image sequences,?
IEEE Transactions on Medical Imaging, vol. 11, no. 2, pp. 238?249, 1992.
[42] J. L. Prince and E. R. McVeigh, ?Optical flow for tagged MR images,? IEEE Trans-
actions on Medical Imaging, vol. 7, pp. 2441?2444, 1991.
[43] J. L. Prince, S. N. Gupta, and N. F. Osman, ?Bandpass optical flow for tagged MRI.?
Med Phys, vol. 27, no. 1, pp. 108?118, 2000.
[44] M. A. Guttman, J. L. Prince, and E. R. Mcveigh, ?Tag and contour detection in tagged
MR images of the left ventricle,? IEEE Transactions on Medical Imaging, vol. 13, no. 1,
pp. 74?88, 1994.
[45] T. S. Denney, ?Estimation and detection of myocardial tags in MR image without
user-defined myocardial contours,? IEEE Transactions on Medical Imaging, vol. 18,
no. 4, pp. 330?334, 1999.
[46] Z. Qian, A. Montillo, D. Metaxas, and L. Axel, ?Segmenting MRI tagging lines us-
ing gabor filter banks,? International Conference on the Engineering in Medical and
Biology Society, pp. 630?633, 2003.
[47] D. Xiang and T. S. Denney, ?Combined tag tracking and myocardium motion recon-
struction from planar tagged MR image data without user-defined myocardial con-
tours,? in Biomedical Imaging: Nano to Macro, 2004. IEEE International Symposium
on, 2004, pp. 888?891 Vol. 1.
139
[48] L. Xiaofeng, E. Murano, M. Stone, and J. L. Prince, ?HARP tracking refinement
using seeded region growing,? in Biomedical Imaging: From Nano to Macro, 2007.
ISBI 2007. 4th IEEE International Symposium on, 2007, pp. 372?375.
[49] N. F. Osman and E. R. McVeigh, ?Imaging heart motion using harmonic phase MRI,?
IEEE Transactions on Medical Imaging, vol. 19, no. 3, pp. 186?202, 2000.
[50] N. F. Osman, W. S. Kerwin, E. R. McVeigh, and J. L. Prince, ?Cardiac motion tracking
using cine harmonic phase (HARP) magnetic resonance imaging,? Magnetic Resonance
in Medicine, vol. 42, no. 6, pp. 1048?1060, 1999.
[51] W. Liu, J. Chen, S. Ji, J. S. Allen, P. V. Bayly, S. A. Wickline, and X. Yu, ?Harmonic
phase MR tagging for direct quantification of Lagrangian strain in rat hearts after
myocardial infarction.? Magn Reson Med, vol. 52, no. 6, pp. 1282?1290, 2004.
[52] A. Khalifa, A. B. M. Youssef, and N. Osman, ?Improved Harmonic Phase (HARP)
method for motion tracking a tagged cardiac MR images.? Conf Proc IEEE Eng Med
Biol Soc, vol. 4, pp. 4298?4301, 2005.
[53] J. Goutsias and S. Batman, ?Morphological methods for biomedical image analysis.?
in Handbook of Medical Imaging: Medical Imaging Processing and Analysis, M. Sonka
and J. Fitzpatrick, Eds. SPIE, 2000, pp. 255?263.
[54] Y. Chen and A. Amini, ?A MAP framework for tagline detetction in SPAMM data
using Markov random fields on the B-spline solid,? IEEE Workshop on Mathematical
Models in Biomedical Image Analysis, pp. 131?138, 2001.
[55] D. N. Metaxas, L. Axel, H. Zhenhua, A. Montillo, K. Park, and Q. Zhen, ?Segmen-
tation and analysis of 3D cardiac motion from tagged MRI images,? in Engineering
in Medicine and Biology Society, 2003. Proceedings of the 25th Annual International
Conference of the IEEE, vol. 1, 2003, pp. 122?125 Vol.1.
[56] D. N. Metaxas, L. Axel, Z. Qian, and X. Huang, ?A segmentation and tracking system
for 4D cardiac tagged MR images,? in Engineering in Medicine and Biology Society,
2006. EMBS ?06. 28th Annual International Conference of the IEEE, 2006, pp. 1541?
1544.
[57] Z. Qian, D. Metaxas, and L. Axel, ?A learning framework for the automatic and accu-
rate segmentation of cardiac tagged MRI images,? in Computer Vision for Biomedical
Image Applications, 2005, pp. 93?102, 10.1007/11569541 11.
[58] W. Feng, T. S. Denney, S. Lloyd, L. Dell?Italia, and H. Gupta, ?Contour regularized
left ventricular strain analysis from cine MRI,? in Biomedical Imaging: From Nano to
Macro, 2008. ISBI 2008. 5th IEEE International Symposium on, 2008, pp. 520?523.
[59] L. Axel, A. Montillo, and D. Kim, ?Tagged magnetic resonance imaging of the heart:
A survey,? Medical Image Analysis, vol. 9, pp. 376?393, 2005.
140
[60] C. A. Davis, L. Jin, and T. S. Denney, ?Analysis of spectral changes and filter design
in tagged cardiac MRI,? in Biomedical Imaging: Nano to Macro, 2006. 3rd IEEE
International Symposium on, 2006, pp. 137?140.
[61] M. Marinelli, V. Positand, N. F. Osman, F. A. Recchia, M. Lombardi, and L. Landini,
?Automatic filter design in HARP analysis of tagged magnetic resonance images,? in
Biomedical Imaging: From Nano to Macro, 2008. ISBI 2008. 5th IEEE International
Symposium on, 2008, pp. 1429?1432.
[62] A. A. Amir, W. C. Rupert, and C. G. John, ?Snakes and splines for tracking non-rigid
heart motion,? 1996, 649040 251-261.
[63] A. A. Amini, Y. Chen, J. Sun, and V. Mani, ?Constrained thin-plate spline recon-
struction of 2D deformations from tagged MRI,? 1998, 836884 36.
[64] A. A. Amini, C. Yasheng, R. W. Curwen, V. Mani, and J. Sun, ?Coupled B-snake grids
and constrained thin-plate splines for analysis of 2-D tissue deformations from tagged
MRI,? Medical Imaging, IEEE Transactions on, vol. 17, no. 3, pp. 344?356, 1998.
[65] P. Radeva, A. Amini, H. Jiantao, and E. Marti, ?Deformable B-solids and implicit
snakes for localization and tracking of SPAMM MRI-data,? in Mathematical Methods
in Biomedical Image Analysis, 1996., Proceedings of the Workshop on, 1996, pp. 192?
201.
[66] W. S. Kerwin and J. L. Prince, ?Cardiac material markers from tagged MR images,?
Medical Image Analysis, vol. 2, no. 4, pp. 339?353, 1998.
[67] A. A. Young, D. L. Kraitchman, L. Dougherty, and L. Axel, ?Tracking and finite ele-
ment analysis of stripe deformation in magnetic resonance tagging,? Medical Imaging,
IEEE Transactions on, vol. 14, no. 3, pp. 413?421, 1995.
[68] A. A. Young, ?Model tags: direct three-dimensional tracking of heart wall motion from
tagged magnetic resonance images,? Medical Image Analysis, vol. 3, no. 4, pp. 361?372,
1999.
[69] W. G. O?Dell, C. C. Moore, W. C. Hunter, E. A. Zerhouni, and E. R. McVeigh, ?Three-
dimensional myocardial deformations: calculation with displacement field fitting to
tagged MR images,? Radiology, vol. 195, no. 3, pp. 829?835, 1995.
[70] R. Chandrashekara, R. H. Mohiaddin, and D. Rueckert, ?Analysis of 3-D myocardial
motion in tagged MR images using nonrigid image registration,? Medical Imaging,
IEEE Transactions on, vol. 23, no. 10, pp. 1245?1250, 2004.
[71] D. Perperidis, R. H. Mohiaddin, and D. Rueckert, ?Spatio-temporal free-form regis-
tration of cardiac MR image sequences,? Medical Image Analysis, vol. 9, no. 5, pp.
441?456, 2005.
141
[72] C. Ozturk and E. R. McVeigh, ?Four-dimensional B-spline based motion analysis of
tagged MR images: introduction and in vivo validation,? Physics in Medicine and
Biology, vol. 45, pp. 1683?1702, 2000.
[73] J. Huang, D. Abendschein, V. G. Davila-Roman, and A. A. Amini, ?Spatio-temporal
tracking of myocardial deformations with a 4-D B-spline model from tagged MRI,?
Medical Imaging, IEEE Transactions on, vol. 18, no. 10, pp. 957?972, 1999.
[74] N. J. Tustison and A. A. Amini, ?Biventricular myocardial strains via nonrigid regis-
tration of AnFigatomical NURBS models,? Medical Imaging, IEEE Transactions on,
vol. 25, no. 1, pp. 94?112, 2006.
[75] N. F. Osman and J. L. Prince, ?Visualizing myocardial function using HARP MRI,?
Phys. Med. Biol., vol. 45, pp. 1665?1682, 2000.
[76] Z. A.-E. Khaled, P. Vijay, and L. P. Jerry, ?Artifact reduction in HARP strain maps
using anisotropic smoothing,? M. Armando and A. A. Amir, Eds., vol. 6143. SPIE,
2006, p. 61432P, medical Imaging 2006: Physiology, Function, and Structure from
Medical Images 1.
[77] S. R. Tecelao, J. J. M. Zwanenburg, J. P. A. Kuijer, and J. T. Marcus, ?Extended
harmonic phase tracking of myocardial motion: improved coverage of myocardium and
its effects on strain results,? J Magn Res Im, vol. 23, pp. 682 ? 690, 2006.
[78] M. J. Ledesma-Carbayo, J. A. Derbyshire, S. Sampath, A. Santos, M. Desco, and
E. R. McVeigh, ?Unsupervised estimation of myocardial displacement from tagged
MR sequences using nonrigid registration,? Magnetic Resonance in Medicine, vol. 59,
no. 1, pp. 181?189, 2008, 10.1002/mrm.21444.
[79] S. Ryf, M. A. Spiegel, M. Gerber, and P. Boesiger, ?Myocardial tagging with 3D-
CSPAMM,? Journal of Magnetic Resonance Imaging, vol. 16, no. 3, pp. 320?325,
2002, 10.1002/jmri.10145.
[80] N. Pelc, R. Herfkens, A. Shimakawa, and D. Enzmann, ?Phase contrast cine magnetic
resonance imaging,? Magn Reson Q., vol. 7, no. 4, pp. 229?54, 1991.
[81] E. Bergvall, P. Cain, H. Arheden, and G. Sparr, ?A fast and highly automated ap-
proach to myocardial motion analysis using phase contrast magnetic resonance imag-
ing,? Journal of Magnetic Resonance Imaging, vol. 23, no. 5, pp. 652?661, 2006,
10.1002/jmri.20565.
[82] B. S. Spottiswoode, X. Zhong, A. T. Hess, C. M. Kramer, E. M. Meintjes, B. M.
Mayosi, and F. H. Epstein, ?Tracking myocardial motion from cine DENSE images
using spatiotemporal phase unwrapping and temporal fitting,? Medical Imaging, IEEE
Transactions on, vol. 26, no. 1, pp. 15?30, 2007.
[83] A. H. Aletras, S. Ding, R. S. Balaban, and H. Wen, ?DENSE: Displacement encoding
with stimulated echoes in cardiac functional MRI,? Journal of Magnetic Resonance,
vol. 137, no. 1, pp. 247?252, 1999.
142
[84] D. Kim, W. D. Gilson, C. M. Kramer, and F. H. Epstein, ?Myocardial tissue tracking
with two-dimensional cine displacement-encodedMR imaging: Development and initial
evaluation,? Radiology, p. 2303021213, 2004.
[85] S. Sampath, N. F. Osman, and J. L. Prince, ?A combined harmonic phase and strain-
encoded pulse sequence for measuring three-dimensional strain,? Magnetic Resonance
Imaging, vol. 27, no. 1, pp. 55?61, 2009.
[86] K. Z. Abd-Elmoniem, M. Stuber, and J. L. Prince, ?Direct three-dimensional myocar-
dial strain tensor quantification and tracking using zharp,? pp. 778?786, 2008.
[87] S. Sampath and J. L. Prince, ?Automatic 3d tracking of cardiac material markers using
slice-following and harmonic-phase MRI,? Magnetic Resonance Imaging, vol. 25, no. 2,
pp. 197?208, 2007.
[88] A. H. Aletras, R. S. Balaban, and H. Wen, ?High-resolution strain analysis of the
human heart with Fast-DENSE,? Journal of Magnetic Resonance, vol. 140, no. 1, pp.
41?57, 1999.
[89] D. C. Ghiglia and M. D. Pritt, Two-dimensional phase unwrapping : theory, algo-
rithms, and software. New York: Wiley, 1998.
[90] B. Ambale, H. Gupta, S. G. Lloyd, L. D. ?Italia, and T. S. Denney, ?3D
left ventricular strain from unwrapped harmonic phase measurements,? Journal of
Magnetic Resonance Imaging, vol. 31, no. 4, pp. 854?862, 2010. [Online]. Available:
http://dx.doi.org/10.1002/jmri.22099
[91] B. Ambale, T. S. Denney, H. Gupta, S. Lloyd, and L. Dell?Italia, ?Measuring 3D left
ventricular strain from unwrapped harmonic phase,? in Biomedical Imaging: From
Nano to Macro, 2008. ISBI 2008. 5th IEEE International Symposium on, 2008, pp.
1433?1436.
[92] B. Ambale, S. Lloyd, T. S. Denney, L. Dell?Italia, R. Benza, and H. Gupta, ?Com-
parison of 2D and 3D torsion measured from tagged MRI,? in 17th Meeting of the
International Society for Magnetic Resonance in Medicine, April 2009.
[93] B. Ambale, S. Lloyd, T. S. Denney, L. Dell?Italia, R. Benza, and H. Gupta, ?4D right
ventricular strain and geometry in pulmonary hypertension and normals,? in 18th
Meeting of the International Society for Magnetic Resonance in Medicine, April 2010.
[94] B. Ambale, S. Lloyd, T. S. Denney, L. Dell?Italia, R. Benza, and H. Gupta, ?3D right
ventricular strain and geometry in pulmonary hypertension and normals,? in 16th
Meeting of the International Society for Magnetic Resonance in Medicine, May 2008.
[95] B. Ambale, T. S. Denney, H. Gupta, S. Lloyd, and L. Dell?Italia, ?3D left ventricu-
lar strain by phase unwrapping: a simulated annealing based branch-cut placement
method,? in Proceedings of the Sixth IEEE international conference on Symposium
on Biomedical Imaging: From Nano to Macro. Boston, Massachusetts, USA: IEEE
Press, 2009, pp. 466?469.
143
[96] W. Xu and I. Cumming, ?A region-growing algorithm for InSAR phase unwrapping,?
Geoscience and Remote Sensing, IEEE Transactions on, vol. 37, no. 1, pp. 124?134,
1999.
[97] D. C. Ghiglia and L. A. Romero, ?Robust two-dimensional weighted and unweighted
phase unwrapping that uses fast transforms and iterative methods,? Journal- Optical
Society OF America, vol. 11, pp. 1?107, 1994.
[98] S. Moon-Ho Song, S. Napel, N. Pelc, and G. Glover, ?Phase unwrapping of MR phase
images using poisson equation,? Image Processing, IEEE Transactions on, vol. 4, no. 5,
pp. 667?676, 1995.
[99] U. Spagnolini, ?2-d phase unwrapping and instantaneous frequency estimation,? Geo-
science and Remote Sensing, IEEE Transactions on, vol. 33, no. 3, pp. 579?589, 1995.
[100] T. Flynn, ?Consistent 2-D phase unwrapping guided by a quality map,? Geoscience and
Remote Sensing Symposium, 1996. IGARSS ?96. ?Remote Sensing for a Sustainable
Future.?, International, vol. 4, pp. 2057?2059vol.4, 1996.
[101] S. Kim and Y.-S. Kim, ?Least squares phase unwrapping in wavelet domain,? Vision,
Image and Signal Processing, IEE Proceedings-, vol. 152, no. 3, pp. 261?267, 2005.
[102] B. Friedlander and J. Francos, ?Model based phase unwrapping of 2-D signals,? Signal
Processing, IEEE Transactions on [see also Acoustics, Speech, and Signal Processing,
IEEE Transactions on], vol. 44, no. 12, pp. 2999?3007, 1996.
[103] N. Ching, D. Rosenfeld, and M. Braun, ?Two-dimensional phase unwrapping using a
minimum spanning tree algorithm,? Image Processing, IEEE Transactions on, vol. 1,
no. 3, pp. 355?365, 1992.
[104] R. Sivley and J. Havlicek, ?Multidimensional phase unwrapping for consistent APF
estimation,? Image Processing, 2005. ICIP 2005. IEEE International Conference on,
vol. 2, pp. II?458?61, 2005.
[105] C. Chen and H. Zebker, ?Phase unwrapping for large SAR interferograms: statistical
segmentation and generalized network models,? Geoscience and Remote Sensing, IEEE
Transactions on, vol. 40, no. 8, pp. 1709?1719, 2002.
[106] H. A. Zebker and Y. Lu, ?Phase unwrapping algorithms for radar in-
terferometry: residue-cut, least-squares, and synthesis algorithms,? J. Opt.
Soc. Am. A, vol. 15, no. 3, pp. 586?598, 1998. [Online]. Available:
http://josaa.osa.org/abstract.cfm?URI=josaa-15-3-586
[107] B. R. Hunt, ?Matrix formulation of the reconstruction of phase values from phase
differences,? Journal of the Optical Society of America (1917-1983), vol. 69, p. 393,
1979.
144
[108] M. Pritt, ?Comparison of path-following and least-squares phase unwrapping algo-
rithms,? Geoscience and Remote Sensing, 1997. IGARSS ?97. ?Remote Sensing - A
Scientific Vision for Sustainable Development?., 1997 IEEE International, vol. 2, pp.
872?874vol.2, 1997.
[109] T. S. Denney and J. L. Prince, ?Optimal brightness functions for optical flow estimation
of deformable motion,? Image Processing, IEEE Transactions on, vol. 3, no. 2, pp.
178?191, 1994.
[110] H. Wen, K. A. Marsolo, E. E. Bennett, K. S. Kutten, R. P. Lewis,
D. B. Lipps, N. D. Epstein, J. F. Plehn, and P. Croisille, ?Adaptive
postprocessing techniques for myocardial tissue tracking with displacement-encoded
MR imaging,? Radiology, vol. 246, no. 1, pp. 229?240, 2008. [Online]. Available:
http://radiology.rsnajnls.org/cgi/content/abstract/246/1/229
[111] D. Xiang and T. S. Denney, ?Three-dimensional myocardial strain reconstruction from
tagged MRI using a cylindrical B-spline model,? Medical Imaging, IEEE Transactions
on, vol. 23, no. 7, pp. 861?867, 2004.
[112] M. D. Pritt, ?Phase unwrapping by means of multigrid techniques for interferometric
SAR,? Geoscience and Remote Sensing, IEEE Transactions on, vol. 34, no. 3, pp.
728?738, 1996.
[113] J. P. Kuijer, M. B. M. Hofman, J. J. M. Zwanenburg, J. T. Marcus, A. C. v. Rossum,
and R. M. Heethaar, ?DENSE and HARP: Two views on the same technique of phase-
based strain imaging,? Journal of Magnetic Resonance Imaging, vol. 24, no. 6, pp.
1432?1438, 2006, 10.1002/jmri.20749.
[114] T. S. Denney and E. R. McVeigh, ?Model-Free reconstruction of three-dimensional
myocardial strain from planar tagged MR images,? J Magn Reson Imaging., vol. 7,
no. 5, pp. 799?810, 1997.
[115] W. Feng, H. Nagaraj, H. Gupta, S. Lloyd, I. Aban, G. Perry, D. Calhoun,
L. Dell?Italia, and T. S. Denney, ?A dual propagation contours technique for
semi-automated assessment of systolic and diastolic cardiac function by CMR,?
Journal of Cardiovascular Magnetic Resonance, vol. 11, no. 1, p. 30, 2009. [Online].
Available: http://www.jcmr-online.com/content/11/1/30
[116] J. Li and T. S. Denney, ?Left ventricular motion reconstruction with a prolate
spheroidal B-spline model.? Phys Med Biol, vol. 51, no. 3, pp. 517?537, 2006.
[117] N. Osman and J. Prince, ?On the design of the bandpass filters in harmonic phase
MRI,? Image Processing, 2000. Proceedings. 2000 International Conference on, vol. 1,
pp. 625?628vol.1, 2000.
[118] T. S. Denney, L. Yan, and B. L. Gerber, ?Unsupervised reconstruction of a
three-dimensional left ventricular strain from parallel tagged cardiac images,?
Magnetic Resonance in Medicine, vol. 49, no. 4, pp. 743?754?, 2003. [Online].
Available: http://dx.doi.org/10.1002/mrm.10434
145
[119] A. N. Borg, J. L. Harrison, R. A. Argyle, and S. G. Ray, ?Left ventricular torsion in
primary chronic mitral regurgitation,? Heart, vol. 94, no. 5, pp. 597?603, 2008.
[120] P. P. Sengupta, A. J. Tajik, K. Chandrasekaran, and B. K. Khandheria, ?Twist me-
chanics of the left ventricle: Principles and application,? JACC: Cardiovascular Imag-
ing, vol. 1, no. 3, pp. 366 ? 376, 2008.
[121] D. Jerome, A. Nicholas, and R. M. Elliot, ?Use of a 4D planispheric transformation for
the tracking and analysis of LV motion with tagged MR images,? C. Chin-Tu and V. C.
Anne, Eds., vol. 3660. SPIE, 1999, pp. 69?80, medical Imaging 1999: Physiology and
Function from Multidimensional Images 1.
[122] Z. A. Fayad, V. A. Ferrari, D. L. Kraitchman, A. A. Young, H. I. Palevsky, D. C.
Bloomgarden, and L. Axel, ?Right ventricular regional function using MR tagging:
Normals versus chronic pulmonary hypertension,? Magn. Reson. Med., vol. 39, no. 1,
pp. 116?123, 1998. [Online]. Available: http://dx.doi.org/10.1002/mrm.1910390118
[123] S. S. Klein, T. P. Graham, and C. H. Lorenz, ?Noninvasive delineation of normal right
ventricular contractile motion with magnetic resonance imaging myocardial tagging,?
Annals of Biomedical Engineering, vol. 26, pp. 756?763, 1998, 10.1114/1.75. [Online].
Available: http://dx.doi.org/10.1114/1.75
[124] A. A. Young, Z. A. Fayad, and L. Axel, ?Right ventricular midwall surface
motion and deformation using magnetic resonance tagging,? Am J Physiol Heart
Circ Physiol, vol. 271, no. 6, pp. H2677?2688, Dec. 1996. [Online]. Available:
http://ajpheart.physiology.org/cgi/content/abstract/271/6/H2677
[125] B. Gutmann and H. Weber, ?Phase unwrapping with the branch-cut method: Clus-
tering of discontinuity sources and reverse simulated annealing,? Appl. Opt., vol. 38,
no. 26, pp. 5577?5593, 1999.
[126] W. Feng, H. Gupta, S. Lloyd, l. D. Italia, and T. S. Denney, ?Myocardial contour
propagation in cine cardiac MRI.? ISMRM, May 2007, p. 766.
[127] M. R. Zile, C. F. Baicu, and W. H. Gaasch, ?Diastolic heart failure?abnormalities in
active relaxation and passive stiffness of the left ventricle,? N Engl J Med, vol. 350,
no. 19, pp. 1953?9, 2004.
[128] M. M. Redfield, S. J. Jacobsen, B. A. Borlaug, R. J. Rodeheffer, and D. A. Kass,
?Age- and gender-related ventricular-vascular stiffening: a community-based study,?
Circulation, vol. 112, no. 15, pp. 2254?62, 2005.
[129] S. P. Bell, L. Nyland, M. D. Tischler, M. McNabb, H. Granzier, and M. M. LeWinter,
?Alterations in the determinants of diastolic suction during pacing tachycardia,? Car-
diology Unit, University of Vermont College of Medicine, Burlington, VT, USA., pp.
235?40, 2000.
146
[130] F. E. Rademakers, M. B. Buchalter, W. J. Rogers, E. A. Zerhouni, M. L. Weisfeldt,
J. L. Weiss, and E. P. Shapiro, ?Dissociation between left ventricular untwisting and
filling - accentuation by cathecholamines,? Circulation, vol. 85, no. 4, pp. 1572 ? 1581,
1992.
[131] S.-J. Dong, P. S. Hees, C. O. Siu, J. L. Weiss, and E. P. Shapiro, ?Mri assessment
of lv relaxation by untwisting rate: a new isovolumic phase measure of tau,? Am J
Physiol Heart Circ Physiol, vol. 281, no. 5, pp. H2002?2009, 2001. [Online]. Available:
http://ajpheart.physiology.org/cgi/content/abstract/281/5/H2002
[132] A. T. Burns, A. La Gerche, D. L. Prior, and A. I. MacIsaac, ?Left
ventricular untwisting is an important determinant of early diastolic function,?
J Am Coll Cardiol Img, vol. 2, no. 6, pp. 709?716, 2009. [Online]. Available:
http://imaging.onlinejacc.org/cgi/content/abstract/2/6/709
[133] W. Hosford. New York, NY: Cambridge University Press, 2005.
[134] J. M. Huyghe, D. H. van Campen, T. Arts, and R. M. Heethaar, ?The constitutive
behaviour of passive heart muscle tissue: a quasi-linear viscoelastic formulation,? J
Biomech, vol. 24, no. 9, pp. 841?9, 1991.
[135] J. Layland, I. S. Young, and J. D. Altringham, ?The length dependence of work pro-
duction in rat papillary muscles in vitro,? J Exp Biol, vol. 198, no. Pt 12, pp. 2491?9,
1995.
[136] C. E. Miller, M. A. Vanni, L. A. Taber, and B. B. Keller, ?Passive stress-strain mea-
surements in the stage-16 and stage-18 embryonic chick heart,? J Biomech Eng, vol.
119, no. 4, pp. 445?51, 1997.
[137] C. E. Miller and C. L. Wong, ?Trabeculated embryonic myocardium shows rapid stress
relaxation and non-quasi-linear viscoelastic behavior,? J Biomech, vol. 33, no. 5, pp.
615?22, 2000.
[138] C. E. Miller, C. L. Wong, and D. Sedmera, ?Pressure overload alters stress-strain
properties of the developing chick heart,? Am J Physiol Heart Circ Physiol, vol. 285,
no. 5, pp. H1849?56, 2003.
[139] J. H. Omens, D. E. Milkes, and J. W. Covell, ?Effects of pressure overload on the
passive mechanics of the rat left ventricle,? Ann Biomed Eng, vol. 23, no. 2, pp. 152?
63, 1995.
[140] T. D. Ryan, E. C. Rothstein, I. Aban, J. A. Tallaj, A. Husain, P. A. Lucchesi, and
L. J. Dell?Italia, ?Left ventricular eccentric remodeling and matrix loss are mediated
by bradykinin and precede cardiomyocyte elongation in rats with volume overload,?
Journal of the American College of Cardiology, vol. 49, no. 7, pp. 811?821, 2007.
[141] W. J. Corin, T. Murakami, E. S. Monrad, O. M. Hess, and H. P.
Krayenbuehl, ?Left ventricular passive diastolic properties in chronic mitral
147
regurgitation,? Circulation, vol. 83, no. 3, pp. 797?807, 1991. [Online]. Available:
http://circ.ahajournals.org/cgi/content/abstract/83/3/797
148