Tag line tracking and Cardiac Motion Modeling from Tagged MRI Except where reference is made to the work of others, the work described in this dissertation is my own or was done in collaboration with my advisory committee. This dissertation does not include proprietary or classified information. Jin Li Certificate of Approval: Stanley J. Reeves Professor Electrical and Computer Engineering Thomas S. Denney, Jr., Chair Professor Electrical and Computer Engineering Soo-Young Lee Professor Electrical and Computer Engineering Stephen L. McFarland Dean Graduate School Tag line tracking and Cardiac Motion Modeling from Tagged MRI Jin Li 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 Dec. 15, 2006 Tag line tracking and Cardiac Motion Modeling from Tagged MRI Jin Li Permission is granted to Auburn University to make copies of this dissertation at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iii Vita Jin Li, daughter of Zhikuan Li and Cuiping Wang, was born in Pingdingshan, Henan Province, China. Before she started pursuing her M.S.E.E. and Ph.D. at Auburn University, she was with a joint venture company by Alcatel and China Telecom in Shanghai, China. Her research interests are general medical image processing and analysis, as well as software engineering. Jin married Haihua Yan in 2000. iv Dissertation Abstract Tag line tracking and Cardiac Motion Modeling from Tagged MRI Jin Li Doctor of Philosophy, Dec. 15, 2006 (Master of Science in Electrical Engineering, Auburn University, 2003) (Master of Science in BioMedical Engineering, Southeast University, 1999) 135 Typed Pages Directed by Thomas S. Denney Jr. Magnetic resonance (MR) tagging magnetically labels specified regions of the my- ocardium, which appear in the MR images with a spatially encoded pattern of dark stripes called tag lines. The deformation of these tag lines reflects the deformation of the underlying tissue, making it possible to quantitatively evaluate the regional myocardium deformation and strain. This is particularly valuable in the diagnosis of ischemia and infarction. In this dissertation, three new algorithms are presented to track the tag lines and re- construct the cardiac left ventricle (LV) motion from tagged cardiac MR images. In the new statistical tag point classification algorithm, the candidate tag point positions are not pre-smoothed during tracking, allowing smoothness constraints to be applied only as the de- formation model is fitted to the tag points. Two new algorithms based on three-dimensional (3-D) B-splines in prolate spheroidal coordinates are also presented in order to reconstruct the cardiac LV motion. One of these algorithms is used for 3-D LV motion reconstruction from tracked tag lines. The other combines tag tracking and motion reconstruction. These methods model the left ventricle with prolate spheroidal B-spline models, which o?er sev- eral advantages. First, spatially localized B-spline basis functions o?er better local control v than globally-defined spherical harmonics. Second, their domain more closely matches the shape of the LV wall than Cartesian or cylindrical models. Third, the models can enforce smoothing across the apex and compute the strain at that location. Human and animal MR studies and simulations were used to validate the e?ectiveness of the new methods. The experimental results verified the e?ectiveness of the new tag line tracking and 3-D LV motion reconstruction algorithms. vi Acknowledgments I would like to thank my research advisor, Professor Thomas S. Denney Jr., for his guidance throughout my doctoral program. I thank him for leading me into the research field. His passion for scientific research has set a model for me to follow. I would like to thank my committee members Professor Stanley J. Reeves and Professor Soo-Young Lee for their careful review and helpful suggestions on the dissertation. I would like to thank Professor Amnon J. Meir for serving as the outside reader. I would like to thank the National Institutes of Health, my research advisor, and the Electrical and Computer Engineering department for providing the financial support for my graduate studies. I would like to thank my husband Haihua and my parents for their patience and support. Thank you for sharing my pain and joy. My thanks also go to the classmates and friends in Auburn. Thanks for the precious time we spent together. vii Style manual or journal used Transactions of the American Mathematical Society (together with the style known as ?auphd?). Bibliograpy follows van Leunen?s A Handbook for Scholars. Computer software used The document preparation package T E X (specifically L A T E X) together with the departmental style-file auphd.sty. viii Table of Contents List of Tables xi List of Figures xiv 1 Introduction 1 1.1 CardiacImaging ................................. 2 1.1.1 CardiacAnatomy............................. 2 1.1.2 CardiacPhysiology............................ 4 1.1.3 CardiacImagingModalities....................... 6 1.1.4 Cardiac MR Tagging . . . . . . . .................... 10 1.2 Tagged Cardiac MR Image Analysis . . .................... 13 1.3 ContributionsofThisResearch ......................... 15 1.4 OrganizationoftheDissertation ........................ 17 2 Tag Point Classification in Tagged Cardiac MR Images 18 2.1 Introduction.................................... 18 2.2 AlgorithmDevelopment ............................. 18 2.3 ExperimentsandResults............................. 25 2.3.1 ParameterSensitivity .......................... 26 2.3.2 Validation................................. 31 2.4 Discussion..................................... 34 2.5 Appendix ..................................... 36 3 Left Ventricular Motion Reconstruction 37 3.1 Introduction.................................... 37 3.2 Background.................................... 38 3.2.1 The Prolate Spheroidal Coordinate System . ............. 38 3.2.2 ContinuumMechanics.......................... 40 3.3 AlgorithmDevelopment ............................. 41 3.3.1 SpatialDisplacementFieldFitting................... 42 3.3.2 MaterialDisplacementFieldFitting .................. 47 3.3.3 LagrangianStrainComputation .................... 48 3.3.4 Parameters ................................ 50 3.4 ExperimentalMethods.............................. 51 3.5 Results....................................... 52 3.5.1 CenterSensitivityTest.......................... 52 3.5.2 ValidationonInVivoHumanandCanineData............ 53 ix 3.5.3 StrainMapsonaNormalHeart..................... 54 3.5.4 NoiseImmunity.............................. 58 3.5.5 SpatialResolution ............................ 63 3.5.6 PathologyDatabase ........................... 65 3.6 Discussion..................................... 67 4 Combined Tag Tracking and Strain Reconstruction 72 4.1 Introduction.................................... 72 4.2 AlgorithmDevelopment ............................. 73 4.2.1 CoordinateSystemandInitialization.................. 73 4.2.2 ImagePre-processing........................... 75 4.2.3 DeformationModel............................ 75 4.2.4 ModelFitting............................... 79 4.3 ParameterSensitivity .............................. 83 4.4 ValidationMethods................................ 85 4.4.1 AnalysisofTagTrackingErrors .................... 91 4.4.2 3-DStrainReconstruction........................ 92 4.4.3 StatisticalMethods............................ 92 4.5 Validation Results on parallel-tagged MRI data . . . ............. 92 4.5.1 AnalysisofTagTrackingErrors .................... 92 4.5.2 3-DStrainReconstruction........................ 93 4.5.3 Statistics ................................. 97 4.6 Validation Results on Grid-tagged Human Imaging Studies . . . ...... 100 4.6.1 AnalysisofTagTrackingErrors .................... 100 4.6.2 3-DStrainReconstruction........................ 102 4.6.3 Statistics ................................. 102 4.7 Discussion..................................... 108 4.8 Appendix ..................................... 111 5 Conclusion and future work 112 Bibliography 114 x List of Tables 2.1 End-systolic tag point classification accuracies with various annealing rates (?) and empirical parameter values ? = 10, T init = 10, T final =0.5, T SAinit2 =3,T LAinit2 =2,nCP SA =6andnCP LA =8overthe10hu- mandevelopmentimagingstudies. ....................... 27 2.2 End-systolic tag point classification accuracies with various initial regulariza- tion weights (?), optimal annealing rate ? =0.75, and empirical parameter values T init = 10, T final =0.5, T SAinit2 =3,T LAinit2 =2,nCP SA =6and nCP LA =8overthe10humandevelopmentimagingstudies. ........ 27 2.3 End-systolic tag point classification accuracies with various initial temper- atures (T init ), optimal parameters ? =0.75 and ? =10 3 , and empirical parameter values T final =0.5, T SAinit2 =3,T LAinit2 =2,nCP SA =6and nCP LA =8overthe10humandevelopmentimagingstudies. ........ 28 2.4 End-systolic tag point classification accuracies with various final tempera- tures (T final ), optimal parameters ? =0.75, ? =10 3 and T init = 10, and empirical parameter values T SAinit2 =3,T LAinit2 =2,nCP SA =6and nCP LA =8overthe10humandevelopmentimagingstudies. ........ 28 2.5 End-systolic tag point classification accuracies with various initial temper- atures (T SAinit2 ) for the subsequent short-axis slices, optimal parameters ? =0.75, ? =10 3 , T init =10andT final =0.001, and empirical parameter values nCP SA = 6 over the short-axis slices of the 10 human development imagingstudies................................... 29 2.6 End-systolic tag point classification accuracies with various initial temper- atures (T LAinit2 ) for the subsequent long-axis slices, optimal parameters ? =0.75, ? =10 3 , T init =10andT final =0.001, and empirical parame- ter values nCP LA = 8 over the long-axis slices of the 10 human development imagingstudies................................... 29 2.7 End-systolic tag point classification accuracies with various numbers of con- trol points (nCP SA ) for short-axis slices, optimal parameters ? =0.75, ? =10 3 , T init = 10, T final =0.001, and T SAinit2 =4overtheshort-axis slicesofthe10humandevelopmentimagingstudies. ............. 30 xi 2.8 End-systolic tag point classification accuracies with various numbers of con- trol points (nCP LA ) for 2-chamber and 4-chamber slices, optimal parameters ? =0.75, ? =10 3 , T init = 10, T final =0.001, and T LAinit2 =4overthe2- chamber and 4-chamber slices of the 10 human development imaging studies. 30 2.9 End-systolic tag point classification accuracies over the 10 human develop- mentand10humanevaluationimagingstudies................. 31 2.10 The 20 grid-tagged human imaging studies for validating the tag classification algorithm...................................... 36 3.1 ParametersforthequadraticB-splinereconstructions............. 51 4.1 RMS tag distance error in pixels of the quadratic B-spline reconstructions with 3 radial ? 4 longitudinal ? 6 circumferential control points and various ? over the 10 parallel-tagged human development imaging studies. . . . . . 85 4.2 RMS tag distance error in pixels of the quadratic closed-apex PS-COTTER reconstruction on raw image data with various numbers of control points and regularization weight ? =10 ?5 over the 10 parallel-tagged human develop- mentimagingstudies. .............................. 86 4.3 RMS tag distance error in pixels of the quadratic open-apex PS-COTTER reconstruction on raw image data with various numbers of control points and regularization weight ? =10 ?7 over the 10 parallel-tagged human develop- mentimagingstudies. .............................. 86 4.4 RMS tag distance error in pixels of the quadratic closed-apex PS-COTTER reconstruction on detetced tag points with various numbers of control points and regularization weight ? =10 ?4 over the 10 parallel-tagged human devel- opmentimagingstudies.............................. 87 4.5 RMS tag distance error in pixels of the quadratic open-apex PS-COTTER reconstruction on detetced tag points with various numbers of control points and regularization weight ? =10 ?7 over the 10 parallel-tagged human devel- opmentimagingstudies.............................. 87 4.6 Optimal parameters and the corresponding RMS tag distance errors in pixels of the quadratic PS-COTTER reconstructions over the 10 parallel-tagged humandevelopmentimagingstudies....................... 88 4.7 RMS tag distance error in pixels of the quadratic B-spline reconstructions with 3 radial ? 4 longitudinal ? 6 circumferential control points and various ? over the 12 grid-tagged human development imaging studies. . ...... 88 xii 4.8 RMS tag distance error in pixels of the quadratic closed-apex PS-COTTER reconstruction on filtered image data with various numbers of control points and regularization weight ? =10 ?6 over the 12 grid-tagged human develop- mentimagingstudies. .............................. 88 4.9 RMS tag distance error in pixels of the quadratic open-apex PS-COTTER reconstruction on filtered image data with various numbers of control points and regularization weight ? =10 ?6 over the 12 grid-tagged human develop- mentimagingstudies. .............................. 89 4.10 RMS tag distance error in pixels of the quadratic closed-apex PS-COTTER reconstruction on detected tag points with various numbers of control points and regularization weight ? =10 ?4 over the 12 grid-tagged human develop- mentimagingstudies. .............................. 90 4.11 RMS tag distance error in pixels of the quadratic open-apex PS-COTTER reconstruction on detetced tag points with various numbers of control points and regularization weight ? =10 ?5 over the 12 grid-tagged human develop- mentimagingstudies. .............................. 90 4.12 Optimal parameters and the corresponding RMS tag distance errors in pixels of the quadratic PS-COTTER reconstructions over the 12 grid-tagged human developmentimagingstudies. .......................... 91 4.13 RMS tag distance errors in pixels of the quadratic PS-COTTER reconstruc- tions over the 8 parallel-tagged human evaluation imaging studies. . . . . . 93 4.14 Correlation Coe?cients (?), Di?erence (Mean ? SD) between the average midwall strains over the 8 parallel-tagged human evaluation imaging stud- ies computed using closed-apex PS-COTTER on raw image data and those computedwithuser-supervisedanalysis(TEAandDMF). ......... 100 4.15 RMS tag distance errors in pixels of the quadratic PS-COTTER reconstruc- tions with the optimal parameters over the 41 grid-tagged human evaluation imagingstudies................................... 102 4.16 Correlation Coe?cients (?), Di?erence (Mean ? SD) between the mid-wall average strains for the 41 grid-tagged human evaluation imaging studies com- puted using closed-apex PS-COTTER on detected tag points and those com- putedwiththeAPSBmodel. .......................... 104 4.17DevelopmentStudies............................... 111 4.18EvaluationStudies ................................ 111 xiii List of Figures 1.1 Cardiacanatomy................................. 3 1.2 Slice prescriptions in cardiac imaging. Left: illustration of the slice prescrip- tions. Right: cardiac MR images of the corresponding slices. The bold lines indiagramsontheleftdenotetheslicepositions................ 5 1.3 Cardiac MR tagging. Left: at end-diastole. Right: at end-systole. Top: illustrations of tag planes and image planes. Bottom: MR images of the LV. 11 1.4 Displacement measurements from planar tagged MR images. . . ...... 12 1.5 Domains of 3-D LV deformation reconstruction models. (a): Cartesian do- main; (b): cylindrical domain; (c): prolate spheroidal domain. The gray objectinsideeachdomainrepresentstheLVwall................ 14 2.1 Tag tracking results from a human heart. (a) Result of deterministic an- nealing algorithm. Synthetic tag lines are shown in magenta. Detected tag points are shown in yellow. (b) Final classification result. Each detected tag line is displayed with a color indicating the tag line to which it belongs. Tag pointsclassifiedasnotbelongingtoataglineareshowninblue. ...... 20 2.2 Tag line registration across multiple short-axis slices. Each detected tag line is displayed with a color indicating the tag line to which it belongs. Tag pointsclassifiedasfalsepositivesareshowninblue. ............. 33 3.1 The prolate spheroidal coordinate system and Cartesian coordinate system. n ? , n ? , n ? are the orthogonal unit vectors in ? (radial), ? (longitudinal), ? (circumferential)directionsrespectively. .................... 39 3.2 Strain error standard deviation (SD) versus center point perturbation SD with the quadratic B-spline reconstructions on the analytically defined math- ematicalLVdeformationmodel. ........................ 53 3.3 Circumferential shortening (E cc ) strain map at the mid-wall for the high resolution normal human study over t =[41,481] ms. ............. 55 3.4 Longitudinal shortening (E ll ) strain map at the mid-wall for the high reso- lution normal human study over t =[41,481] ms. . . ............. 56 xiv 3.5 Radial thickening (E rr ) strain map at the mid-wall for the high resolution normal human study over t =[41,481] ms. . . . . . . ............. 57 3.6 End-systolic strain maps of radial thickening (E rr , left), circumferential shortening (E cc , middle) and longitudinal shortening (E ll , right) computed with the APSB (top) and PSB (bottom) models on an open-apex mesh. The redballsintheplotsdenotetheLVmid-septum. ............... 59 3.7 End-systolic strain maps of radial thickening (E rr , left), circumferential shortening (E cc , middle) and longitudinal shortening (E ll , right) computed with the cylindrical B-spline model (bottom) and its regularized version (top) on an open-apex mesh. The red balls in the plots denote the LV mid-septum. 60 3.8 End-systolic strain maps of radial thickening (E rr , left), circumferential shortening (E cc , middle) and longitudinal shortening (E ll , right) computed with the APSB (top) and PSB (bottom) models on a closed-apex mesh. The redballsintheplotsdenotetheLVmid-septum. ............... 61 3.9 The mean and standard deviation on error on Lagrangian strain E in a local coordinate system (radial, circumferential or longitudinal). The center of each errorbar is on the mean value and the bar extends to one standard deviationaboveandbelowthemean. ..................... 62 3.10 Spatial resolution for each displacement direction (radial, circumferential or longitudinal). Top: for di?erent cloud widths, the per cent of displacement response for each method; bottom, the corresponding values in mm for a responseof50%. ................................. 64 3.11 RMS error on tag distances for di?erent types of motion: human normal, infarcted or dilated cardiomyopathy, canine debutamine, RA or RV paced. Each set of error bars represents one patient which is numbered from 1 to 5. 66 4.1 (a): Initial shape of the prolate spheroidal B-spline model. (b): Undeformed tag planes defined within the prolate spheroidal B-spline model as uniformly spacedmaterialplanes............................... 74 4.2 A grid-tagged MR image and its spectrum (top), and the spectrums of Gaus- sian filters (bottom) for obtaining the parallel-tagged images. . . ...... 76 4.3 Grid-tagged MR images (left) and the corresponding filtered images (middle, right)nearend-diastole(top)andnearend-systole(bottom). ........ 77 4.4 Tagpointsdetectedbyusingthepeakdetectionalgorithm[48]........ 77 4.5 Simulated tag points on the intersection lines of the imaging planes and tag planes........................................ 80 xv 4.6 Simulated tag lines of PS-COTTER for a mid-ventricular short-axis slice with a tag plane angle of 0 o (left) and 90 o (middle), and one long-axis slice (right). The top, middle and bottom rows correspond to early, mid- and late systole. 94 4.7 RMS tag di?erence in pixels between the simulated tags generated with closed-apex PS-COTTER on raw image data and a set of tag lines tracked with user-supervised processing for the 8 parallel-tagged human evaluation imagingstudies................................... 95 4.8 Maps of 3-D radial thickening (E rr ), circumferential shortening (E cc ), and longitudinal shortening (E ll ) reconstructed with closed-apex PS-COTTER on raw image data for a normal human (left) and a patient with an anterior myocardial infarction (right). The red balls denote the septum and the green ball denotes the apex. Each image is color-coded, with yellow indicating thickeningandblueindicatingshortening.................... 96 4.9 Correlation(left) and Bland-Altman (right) plots for comparing the aver- age midwall strains computed using closed-apex PS-COTTER on raw image data and those computed using TEA over the 8 parallel-tagged human eval- uation studies. In the correlation plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence line represent the mean ? 2SDofthedi?erence........................... 98 4.10 Correlation(left) and Bland-Altman (right) plots for comparing the aver- age midwall strains computed using closed-apex PS-COTTER on raw image data and those computed using DMF over the 8 parallel-tagged human eval- uation studies. In the correlation plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence line represent the mean ? 2SDofthedi?erence........................... 99 4.11 Simulated tag lines of PS-COTTER for a mid-ventricular short-axis slice (left), a 2-chamber slice (middle), and a 4-chamber slice (right). The top two andbottomtworowscorrespondtoearlyandlatesystole........... 101 4.12 Maps of 3-D radial thickening (E rr , left), circumferential shortening (E cc , middle), and longitudinal shortening (E ll , right) reconstructed with closed- apex PS-COTTER on detected tag points from grid-tagged MR data for a normal human, a hypertensive patient, a diabetic patient, and a patient with mitral regurgitation. The red ball denotes the anterior intersection, the green ball denotes mid-wall and the blue ball denotes the inferior intersection. Each image is color-coded, with yellow indicating thickening and blue indicating shortening. .................................... 103 xvi 4.13 Correlation(left) and Bland-Altman (right) plots over the 41 grid-tagged hu- man evaluation imaging studies for comparing average basal midwall strains computed using closed-apex PS-COTTER on detected tag points and those computed using the APSB algorithm. In the correlation plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence line represent mean ? 2SDofthedi?erence............. 105 4.14 Correlation(left) and Bland-Altman (right) plots over the 41 grid-tagged hu- man evaluation imaging studies for comparing average mid-ventricle midwall strains computed using closed-apex PS-COTTER on detected tag points and those computed using the APSB algorithm. In the correlation plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence line represent mean ? 2SDofthedi?erence. ........ 106 4.15 Correlation(left) and Bland-Altman (right) plots over the 41 grid-tagged hu- man evaluation imaging studies for comparing average apical midwall strains computed using closed-apex PS-COTTER on detected tag points and those computed using the APSB algorithm. In the correlation plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence line represent mean ? 2SDofthedi?erence............. 107 xvii Chapter 1 Introduction The estimation of three-dimensional (3-D) cardiac motion and deformation from two- dimensional (2-D) images is currently receiving a great deal of attention in the field of medical image analysis. One well-established technique for non-invasively imaging the car- diac motion is magnetic resonance (MR) imaging. This approach has the advantages of lack of radiation, adjustable soft tissue contrast, good spatial resolution and 3-D imag- ing with highly reproducible measurements. Cardiac MR images can provide a series of measurements that can be used to evaluate both global and local functions of the heart. Tagged MR imaging and analysis has received considerable attention in recent years due to its ability to track material points non-invasively. In tagged MR imaging, specified regions of the heart muscle are magnetically labeled, which appear in the MR images with a spatially encoded pattern of dark stripes called tag lines. Since the labeling temporarily modifies the material properties of the tissue, the deformation of the tag lines reflects the deformation of the underlying tissue. This allows quantitative evaluation of the regional heart wall deformation and strain. In this dissertation, we are concerned with tag line tracking and modeling the defor- mation of the cardiac LV from tagged MR images. This chapter focuses primarily on the background of cardiac imaging and tagged cardiac MR image analysis, and the contributions of this research. This chapter is organized as follows. Section 1.1 reviews cardiac anatomy and cardiac imaging modalities and Section 1.2 discusses existing techniques for tagged cardiac MR 1 image analysis. Then, in Section 1.3, the contributions of this dissertation are outlined. Finally, an overview of the organization of this dissertation is provided in Section 1.4. 1.1 Cardiac Imaging 1.1.1 Cardiac Anatomy The heart is an organ that pumps blood around the body and is located between the lungs, and behind and to the left of the chest bone. As shown in Figure 1.1, it has four chambers: the Right Atrium (RA), Left Atrium (LA), Right Ventricle (RV) and Left Ventricle (LV). The RV and LV share a common wall, named the interventricular septum. Between the chambers there are four one-way valves, which ensure the flow of blood in the proper direction. The top portion of the heart containing the atria and valves is known as the base, and the bottom portion containing the closed ends of the ventricles is called the apex. The heart muscle, or myocardium, contains three layers: the epicardium, midwall, and endocardium. The epicardium is the outermost layer of the heart and consists of a loose connective tissue of fibroblasts, collagen fibers and adipose tissue. The endocardium is the innermost layer of the heart and is composed of a continuous smooth endothelial layer that covers all the inner surfaces of the heart. Between the epicardium and the endocardium lies the main muscle mass of the heart, which is composed primarily of striated muscle cells and fibers forming a conduction system. Figure 1.2 shows both the short-axis and long-axis views of the heart. The long- axis views include 2-chamber, 4-chamber and radial slice prescriptions. In the short-axis prescription, the image slices are perpendicular to the longitudinal direction of the heart. Thus, the LV wall appears circular in the short-axis images. In long-axis prescriptions, the 2 Figure 1.1: Cardiac anatomy 3 image slices are in the longitudinal direction and the LV wall therefore appears elliptical, with a closed end. 1.1.2 Cardiac Physiology The heart works continuously throughout an organism?s entire life span. It performs the central function of the circulatory system, pumping blood through an intricate system of blood vessels to bring life-giving oxygen to every cell in the body. The function of the right side of the heart (the RA and RV) is to receive the oxygen- depleted blood from the superior and inferior vena cava, and to pump it to the lungs. Arriving in the lungs, the blood receives oxygen and disposes of the carbon-dioxide that the body has produced. Then, the oxygenated blood passes to the left side of the heart (the LA and LV), where the LV pumps the oxygenated blood to the rest of the body. As blood from the lungs and the rest of the body enters the atria, each atrium con- tracts, forcing the blood into its respective ventricle. The RV and LV then contract nearly simultaneously and force the blood into the lungs (the RV) and the rest of the body (the LV). Thus, the heart undergoes two motion phases in a cardiac cycle, namely the systolic and diastolic phases. During the systolic phase the oxygenated blood is pumped out; during the diastolic phase, the ventricles fill with blood. When the blood supply for the heart muscle is reduced for any reason, the heart fails to receive enough of the energy it needs to function; thus ischemic heart diseases exhibit reduced contractions of the myocardium. This restriction in blood supply is called ischemia. Cardiac imaging of the heart under conditions of both stress and rest o?ers great potential for assessment and diagnosis of the abnormal cardiac motion and/or underlying myocardial viability of an infarcted heart [1]. 4 Left ventricle Right ventricle Right ventricle Left ventricle Left ventricle Left atrium Left ventricle Right ventricle Left ventricle Left atrium Left ventricle Right ventricle (a) Short-axis prescription (b) 2-chamber prescription (c) 4-chamber prescription (d) Radial slice prescription Right ventricle Right atrium Figure 1.2: Slice prescriptions in cardiac imaging. Left: illustration of the slice prescriptions. Right: cardiac MR images of the corresponding slices. The bold lines in diagrams on the left denote the slice positions. 5 1.1.3 Cardiac Imaging Modalities The main purpose of cardiac imaging is to provide both qualitative and quantitative information about cardiac anatomy and function. There are various imaging modalities that can be used to view the cardiac anatomical structures and measure the heart wall geome- try. These imaging modalities can be split into four categories based on the type of imaging physics involved: X-ray techniques (including angiocardiography and computed tomogra- phy), nuclear imaging, ultrasound (echocardiography) and magnetic resonance imaging. These cardiac imaging modalities will be briefly reviewed in the following sections. Angiocardiography Angiocardiography consists of X-ray radiography with contrast enhancement and can be used to measure ventricular volume and mass. In this technique, a contrast medium is injected before imaging. X-ray images are then generated by projecting a wide beam of X-rays through the object onto a film. The injected substances in the tissue absorb and scatter the incoming energy of the X-rays, resulting in an inverse mapping of tissue opacity. Since X-ray images are composed of the superposition of the projections of di?erent tissues of the imaging volume in the X-ray projection path, the e?ective resolution of this method depends on the relative opacity of the neighboring objects, the dispersion of the X-ray beam, and the distance of objects from the X-ray tube focal spot and the film [2]. Although this technique is basically a 2-D approach, projections of orthogonal views can be used to depict the short-axis and long-axis dimensions of the heart. The main drawbacks of angiocardiography are its use of the ionizing radiation, and its invasive nature due to the need for a contrast injection. The advantage of angiocardiography is that X-ray projection systems are relatively inexpensive, easy to use and can be used to provide real-time feedback. 6 Computed Tomography Computed tomography(CT) is another imaging modality that utilizes X-rays. It is the first non-invasive radiological method that allows the generation of tomographic images. A single-slice CT image is formed by projecting many X-ray beams through the object using a fan-beam geometry. The X-ray source is moved in a circle around the object. At fixed angles of the circle, an X-ray fan-beam is emitted, and detectors positioned around the circumference of the scanner and at the opposite side of the object collect the attenuated signals. Then, di?erent views of the projections are combined to enable the tomographic image to be reconstructed by using a backprojection algorithm. The latest CT techniques include helical CT and multi-slice CT scans. In helical CT, the X-ray source and the detectors continuously rotate around a circle while the object is moved continuously through the circle?s center. In multi-slice CT, several rows of detectors are used to rapidly gather a cone of X-ray data, which is the 2-D projection of the 3-D object. These CT techniques are both important due to their ability to rapidly acquire 3-D data. Nuclear Imaging Nuclear imaging is used to selectively measure functional parameters of di?erent or- gans such as myocardial perfusion (blood flow into the cardiac muscle). It usually uses a gamma-ray emitting tracer bound to a chemical with known physiological properties. The concentration of the tracer in the body is measured as a function of position and time by detecting the emitted gamma-rays. Nuclear imaging modalities include conventional radionuclide imaging, single photon emission computed tomography (SPECT) and positron emission tomography (PET). Dif- ferent radio-active tracers are used in those modalities. Conventional radionuclide imaging 7 and SPECT use single photon tracers. In these scanners, a 2-D Anger scintillation camera is used to detect the gamma rays. PET utilizes the tracers that generate photon pairs after their radionuclide atoms decay. Conventional radionuclide imaging is comparable to conventional X-ray imaging. It produces planar projection images with all the depth in- formation lost. SPECT and PET produce images of slices within the body. In SPECT scanners, the Anger camera is rotated around the body. Since the Anger camera is a 2-D imager, SPECT is fundamentally 3-D. In PET scanners, opposing detectors in the detec- tion ring are combined to detect the pairs of photons that travel in opposite directions. Computed tomography techniques such as filtered projection are applied to reconstruct the tomographic images[3, 4]. Nuclear imaging is invasive since it may involve the injection of radio-active tracers before imaging. The advantage of nuclear imaging, however, is that it can provide functional measurements with a higher signal to noise ratio than is possible by using CT, MRI or ultrasound imaging [3]. Echocardiography(Cardiac Ultrasound) Echocardiography, which is often referred to as cardiac ultrasound, is an imaging tech- nique based on the partial reflection of acoustic waves at tissue boundaries. The di?erence in acoustic properties between adjacent tissues governs the image contrast. The commonly used cardiac ultrasound techniques include M-mode, B-Mode and Doppler scan. AnM-modescanconsistsofapulse-echoscaninwhichthebeampointsinasingle direction and a B-mode scan is a pulse-echo scan in which the beam sweeps across a tissue plane. Both M-scans and B-scans use the amplitude of the reflected echo and the elapsed time from the pulse emission to the reflected echo to reconstruct the image. In color Doppler imaging, the Doppler shift is measured to create a 2-D velocity image, which is color-coded 8 and displayed on the background of a 2-D B-mode tissue image. Transthoracic imaging allows short-axis, long-axis and oblique-image views, but it is limited due to the presence of attenuating tissues in the interface between the ultrasound transducer and the organ of interest. The use of transesophageal imaging overcomes these windowing constraints at the cost of being highly invasive as it places the transducer in close proximity to the heart [5, 6, 2]. Magnetic Resonance Imaging Magnetic resonance imaging (MRI) is based on nuclear magnetic resonance phenomenon. Standard clinical MRI looks at the protons in the nuclei of the hydrogen atoms in water molecules. Each proton has a spin, so the nucleus has a magnetic moment. Under normal circumstances, these spins align randomly in a macroscopic sample and their magnetic moments cancel each other out. If an external static magnetic field is applied, however, these spins precess about the field at an angular frequency known as the Larmor frequency and their magnetic moments line up, which results in a net magnetization of the object. In order to obtain an MR image, circularly polarized magnetic fields are applied to make the frequency of precession spatially dependent. Then, selected regions within the object are excited by applying radio-frequency (RF) pulses. The RF pulses interact with the spins to tip them away from their original orientations. After an RF pulse is turned o?, the spin-spin and spin-lattice interactions cause the magnetic moments to return to their equilibrium state. The precession of the magnetic moments generates a RF electromagnetic signal that can be sensed with an antenna [4]. The MR images are then reconstructed from the signal samples with Fourier reconstruction techniques. Compared with other imaging modalities such as conventional X-rays and CT, MRI has many advantages. First, it does not utilize ionizing radiation and thus avoids any negative 9 e?ects on the human body. Second, excellent contrast between soft tissues can be achieved by MRI and this contrast can be adjusted by applying di?erent pulse sequences. Third, by appropriately steering the magnetic fields, slices can be arbitrarily oriented. Thus, MR images can be obtained over a wide range of orientations without the need to change the position of the patient. 1.1.4 Cardiac MR Tagging MRI can also be used to image objects that are in motion. However, since motion can only be measured perpendicular to the image intensity gradient, and not along it, an aperture problem exists with standard MRI. Other image modalities, for example CT and ultrasound, also su?er from this problem. Tagged MRI is a major advance in this respect. It can be used to measure myocardial motion non-invasively by labeling specified regions of the myocardium which serve as markers during the cardiac cycle. Tagged cardiac MR images are acquired by magnetically labeling planar slices of the heart (called tag planes) near end-diastole and acquiring a series of images in slices perpen- dicular to the tag planes (see [8] for details of the imaging protocol). Figure 1.3 illustrates this process, in which the intersections of tag planes with image planes are tag lines, which appear as dark stripes in the images. Since the labeling temporarily modifies the material properties of the tissue, the deformation of the tag lines reflects the deformation of the underlying tissue. In tagged MR images, only motions perpendicular to the tag planes can be measured. As shown in Figure 1.4, a point identified on a deformed tag line, s, originates from some- where on the reference tag plane, but only the displacement component in the tag plane normal direction is known. Therefore, a 1-D measurement of myocardial motion can be 10 Tag Lines Myocardium Image Plane Reference Tag Planes Tag Lines Deformed Myocardium Deformed Image Plane Tag Planes Deformed Figure 1.3: Cardiac MR tagging. Left: at end-diastole. Right: at end-systole. Top: illustrations of tag planes and image planes. Bottom: MR images of the LV. 11 Image Plane Reference Tag Plane Deformed Tag Plane p s ? ? us n us n Figure 1.4: Displacement measurements from planar tagged MR images. 12 obtained from each point identified on a tag line. 3-D deformation is measured by acquiring images with tag planes in di?erent orientations and imaging multiple slices of the heart. 1.2 Tagged Cardiac MR Image Analysis In tagged MR images, the deformed tag lines can be analyzed in order to reconstruct the 3-D LV deformation and to compute the 3-D LV strain, which is an important measure of regional myocardial contractile performance. The 3-D deformation of the LV wall is typically reconstructed by first tracking the tag lines in each image in the data set and then fitting a deformation model to the tag line positions. It can also be reconstructed by fitting a deformation model directly to the image data, thus combining tag tracking and motion reconstruction. Several tag tracking algorithms have been proposed over the years. The approach typically used is to initialize an active contour [11, 12, 13] or a 3-D deformation model [14, 15] close to the tag lines in the image and then evolve the active contour or deformation model until it matches the tag lines in the image. However, there are some potential problems with this approach. Active contour algorithms use smoothness constraints in order to reduce the e?ects of noise and improve tracking accuracy, but ideally this smoothness should be applied in the deformation model that is fit to the tag line positions. Deformation-model algorithms do not su?er from this problem since they fit a deformation model directly to the image data, but these approaches are sensitive to the registration errors that occur between slices acquired during di?erent breath-holds. Once tag lines have been tracked, a deformation model is fit to the tag line points [16, 13, 17, 14, 18, 15, 19, 20, 21] to reconstruct the LV motion. Model-free methods such as HARP [24] can also be used to measure strain, but these methods are currently only capable of measuring 2-D strains. Many 3-D LV deformation models have been proposed, 13 (a) (b) (c) Figure 1.5: Domains of 3-D LV deformation reconstruction models. (a): Cartesian domain; (b): cylindrical domain; (c): prolate spheroidal domain. The gray object inside each domain represents the LV wall. including methods based on prolate spheroidal basis functions [16], finite elements [13], finite di?erences [17], free-form deformations [22, 23], and B-spline functions [14, 18, 15, 19, 20, 21]. B-spline methods, in particular, have received considerable attention because B-splines are spatially localized and B-spline models can be used to represent complex deformations with relatively few parameters. In the work of Chen and Amini [14], Huang et al [19], ? Ozt?urk and McVeigh [20], and Radeva et al [21], the LV deformation is modeled with either 3-D or 4-D (3-D+t) B-spline tensor products in Cartesian coordinates . The 4-D models employ a 1-D B-spline interpolation over time. The LV, however, is roughly elliptical in shape, and the domain of these Cartesian models is rectangular. This shape mismatch can potentially cause instabilities in the model fit because it is possible for a B-spline span to contain few or no data points [15]. This problem has been addressed by either explicit regularization [21] or adding derived measurements to the periphery of the LV [20]. In 14 an attempt to avoid this shape mismatch problem, Deng and Denney [18] used a 3-D B- spline deformation model defined in cylindrical coordinates. However, although cylindrical models match the LV?s shape near mid-ventricle, they do not match near the apex (see Figure 1.5(a)(b)). Several methods that combine tag tracking and motion reconstruction into a single step have been proposed recently for tagged cardiac MR image analysis. These techniques use either B-spline models [14, 15] or finite elements [25]. The 4-D B-spline model (3-D B-solid + time) [14] and the finite elements method [25] require user-defined myocardial contours in order to limit the computation of image forces to points inside the myocardium or to define the initial model geometry. However, identifying the myocardial contours is both labor intensive and time consuming. UNTETHER [26] is an unsupervised 3-D tagged MR image analysis technique that does not require user defined myocardial contours. In the UNTETHER algorithm, however, tags are tracked in each slice independently, and it is possible that tag line positions could be inconsistent from slice to slice. The COTTER algorithm [15] o?ers another method that does not require user-defined myocardial contours. It fits a 3-D cylindrical cardiac deformation model directly to the image data, which ensures that the tag line positions are consistent from slice to slice. However, since this model is defined in cylindrical coordinates, it incorporates a singularity at the apex because the radial coordinate is zero there, and thus it cannot compute the strain at the apex. 1.3 Contributions of This Research As discussed in Section 1.2, there are several problems with the existing cardiac LV motion reconstruction techniques. First, the existing tag tracking methods either require application of smoothness constraints or are sensitive to the registration errors between slices acquired during di?erent breath-holds. Second, the existing B-spline motion reconstruction 15 models are all defined by using either Cartesian or cylindrical coordinates. Their domains therefore fail to match the shape of the LV, which can lead to inaccuracies in the model fit. Third, the existing combined tag tracking and motion reconstruction approaches either require user-defined contours or include a singularity at the apex. In this dissertation, the objective is the quantitative estimation of the regional my- ocardium deformation, including tag line tracking and cardiac motion modeling. The con- tributions of our work include the development of new tag line tracking and 3-D cardiac motion reconstruction algorithms to address the problems in existing techniques. First, a tag point classification algorithm [47] is presented that is suitable for use in a new technique for tracking tag lines in tagged cardiac MR images. Instead of tracking tag lines from frame to frame in an image sequence with active contours or similar techniques, a set of candidate tag points are detected in each image which are then classified as either a false positive or belonging to a particular tag line. The advantages of this approach are that the tag point positions are not pre-smoothed during tracking, allowing the smoothness constraint to be applied only in the deformation model fit to the tag points, and it is not sensitive to the misregistered slices that are due to a patient breath-hold shift since it is a 2-D method. Second, a new 3-D myocardial displacement and strain reconstruction algorithm has been developed based on B-spline basis functions in prolate spheroidal coordinates. This method can be applied to the tag line points tracked by any tag tracking algorithm. We developed two deformation models ? one is based on the prolate spheroidal B-splines alone, named PSB, the other also incorporates a model for global a?ne deformation, called APSB. The APSB model has been published in [55, 56]. These models more closely match the shape of the LV wall than either of existing Cartesian or cylindrical coordinate models (see Figure 1.5(c)). The PSB and APSB deformation models also enforce smoothness across the apex 16 and can compute strain at that location. Validation experiments have shown that this method can be used to accurately reconstruct deformation and strain in the LV wall from tagged MR images. Third, a new combined tag tracking and strain reconstruction method called PS- COTTER is presented. In this algorithm, the PSB deformation model described above is fit directly to the detected tag points, raw or filtered image data. The tag points can be detected by using peak detection [48, 27], template matching [11], or other techniques. The myocardium motion field and strains are computed based on the results of the tag tracking procedure without the need for intermediate steps. In contrast to the existing techniques in [14, 25], this algorithm does not require user-defined myocardial contours. In contrast to COTTER [15], the domain of PS-COTTER matches the shape of the left ventricle more closely, which makes it perform better than COTTER. Validation experiments have shown that PS-COTTER can accurately reconstruct the deformation field for studies acquired with di?erent imaging protocols. 1.4 Organization of the Dissertation This dissertation is organized as follows. Chapter 2 presents the statistical tag point classification algorithm for tag line tracking and Chapter 3 presents a 3-D myocardial dis- placement and strain reconstruction method based on a?ne deformation and a prolate spheroidal B-spline model. Chapter 4 describes the new combined tag tracking and strain reconstruction method PS-COTTER based on both open-apex and closed-apex PSB mod- els. Finally, conclusions and future work are discussed in Chapter 5. 17 Chapter 2 Tag Point Classification in Tagged Cardiac MR Images 2.1 Introduction In this chapter, a new method is proposed for tracking tag lines in a sequence of tagged cardiac MR images. First, a set of candidate tag line points are identified in an image using peak detection [48, 27], template matching [11], or other techniques. A statistical classifica- tion technique is then used to group the tag points into tag lines, establish correspondence between tag lines at di?erent time frames and remove false positives. This approach allows tag points to be tracked without any smoothness constraints applied in the tracking process. In addition, this approach is robust to large frame-to-frame tag line displacements, which can cause active contour and model-based algorithms to lose tag line correspondence. This chapter is organized as follows. In Section 2.2, we derive the statistical tag point classification algorithm, which is based on a non-rigid feature registration method used in brain mapping [29]. We conduct validation experiments on human tagged cardiac MR image data in Section 2.3, and discuss the results in Section 2.4. 2.2 Algorithm Development In this section, we develop the tag point classification algorithm. We assume that a set of candidate tag line points (yellow-colored points in Figure 2.1(a)) have been identified in the image data using peak detection [48], template matching [11], or other techniques. The classification algorithm will then determine if a tag point belongs to a tag line or not, and, if it does belong to a tag line, which tag line it belongs to. The algorithm begins with a 18 set of ?synthetic? tag lines (magenta-colored lines in Figure 2.1(a)), which are estimates of the positions of the true tag lines. Each candidate tag point in the image is then assigned a set of membership probabilities that specify the probability that the point belongs to each tag line in the image. The positions of the synthetic tag lines are then updated based on the candidate tag points and their membership probabilities. Tag points with a relatively high probability of belonging to a tag line have greater influence over the updated position of the corresponding synthetic tag line. Note that the synthetic tag lines are not intended to exactly match the detected tag points. Rather, their purpose is to act as cluster centers for tag point classification. The two steps described above are iterated with a deterministic annealing schedule for a fixed number of steps. Tag points are then either classified as not belonging to a tag line or assigned to a tag line depending on the final membership probabilities. For simplicity, in the following development, we only present an algorithm for classifying tag points in a parallel-line tag pattern. In grid-tagged images, the classification algorithm would be applied twice ? once for tag points identified in one tag line direction and again to tag points identified in the other direction. The statistical classification algorithm requires two models: a model for updating the synthetic tag line positions based on tag points and membership probabilities, and a model for assigning membership probabilities between tag points and tag lines. To update the positions of the synthetic tag lines, we assume that deformed tag line positions can be described by displacements from a known, straight reference position. The reference position is the position of the tag line when it was first applied to the tissue, and this position is known from the imaging protocol. While tag lines can deform in 2- D in an image, tag lines only provide information about deformation in the undeformed tag line normal direction [58], so we only model deformation of the synthetic tag lines in this direction. We describe the tag line normal displacement with a scalar-valued, 2-D, 19 Short-axis 2-Chamber 4-Chamber End-diastole End-systole (a) Short-axis 2-Chamber 4-Chamber End-diastole End-systole (b) Figure 2.1: Tag tracking results from a human heart. (a) Result of deterministic annealing algorithm. Synthetic tag lines are shown in magenta. Detected tag points are shown in yellow. (b) Final classification result. Each detected tag line is displayed with a color indicating the tag line to which it belongs. Tag points classified as not belonging to a tag line are shown in blue. 20 tensor-product B-spline model in Cartesian coordinates: u ? (r)= N x summationdisplay i=1 N y summationdisplay j=1 B i (?x)B j (?y)C ij , where r is a point on a reference position tag line, N x and N y are the numbers of control points in x (undeformed tag line normal) and y (undeformed tag line parallel) directions, respectively, B i (?) are B-spline functions, C ij are the control points, and (?x, ?y)arethex and y-components of r in the B-spline parametric space. The deformed position of the i-th point on the j-th synthetic tag line (r ij )isgivenby f(r ij )=r ij + u ? (r ij )n, (2.1) where n is a unit vector normal to the undeformed tag lines. To ensure a well-conditioned problem when fitting the B-spline model to tag points, we impose a Laplacian smoothness constraint on the B-spline model, which can be expressed probabilistically as p(f)= 1 Z f exp parenleftBig ??bardblLfbardbl 2 parenrightBig . (2.2) To assign membership probabilities, we assume that the distance between a tag point and a tag line is specified by a Gaussian mixture model where each tag line is a cluster center. We assume that each cluster has the same variance (? 2 ) and that the cluster occupancy probabilities are uniform, i.e., p(j)= 1 N t ,wherej is the index of the j-th tag line and N t is the number of tag lines in the region of interest. Given that a tag point belongs to the j-th tag line, it position, v i , has the probability density p(v i |f,j,?)= 1 ? 2?? exp parenleftBigg ? bardblv i ? f(r ij )bardbl 2 2? 2 parenrightBigg , (2.3) 21 where f(r ij ) is the position of the point closest to v i on the j-th synthetic tag line. Assuming each tag point is independent, the joint distribution of the detected tag points is p(V |f,?)= N p productdisplay i=1 ? ? N t summationdisplay j=1 p(v i |f,j,?)p(j) ? ? , (2.4) where N p is the number of detected points, and V is a vector containing all the v i . With the models described above, the MAP estimate of f is to minimize the objective function E(f)=? N p summationdisplay i=1 log N t summationdisplay j=1 1 ? 2?? exp parenleftbigg ? bardblv i ?f(r ij )bardbl 2 2? 2 parenrightbigg + ?bardblLfbardbl 2 . (2.5) It has been shown in [46, 29] that the minimum of the objective function over f is equivalent to the minimum of E(f,M)= N p summationdisplay i=1 N t summationdisplay j=1 m ij bardblv i ? f(r ij )bardbl 2 + T N p summationdisplay i=1 N t summationdisplay j=1 m ij log m ij + ?bardblLfbardbl 2 (2.6) over f and M subject to the constraint that m ij > 0and summationtext N t j=1 m ij =1,whereT =2? 2 , M is a matrix containing all the m ij . Note that the variance of the tag point distance model, ? 2 has been replaced by the temperature, T, which will be varied in the deterministic annealing algorithm described below. This is a half-quadratic problem. The first term in Equation (2.6) is a data matching term, which forces the synthetic tag line deformation function, f, to match the detected tag points, v i . The influence of each detected tag point on f is weighted by the corresponding membership probability, m ij , which is the probability that the i-th detected tag point belongs to the j-th tag line. The second term is an entropy barrier. The third term regularizes the synthetic tag line deformation function. 22 The objective function E(f,M) is non-convex, so we minimize this function by alter- nately fixing one variable and optimizing with respect to the other and embedding this alternate minimization in a deterministic annealing algorithm similar to the one used in [29]. In a given annealing step, we first update the membership variables with f fixed. We then fix the m ij and update the synthetic tag deformation function with a weighted least-squares fit. These steps are described below: 1. Initialize the temperature, T = T init , and the control points of the B-spline deforma- tion model, C. 2. While T>T final (a) Update the membership variables m ij =exp parenleftBigg ? bardblv i ? f(r ij )bardbl 2 T parenrightBigg , (2.7) and normalize the membership variables first on all the detected points then on all the tag lines: m ij = m ij summationtext N p i=1 m ij , then m ij = m ij summationtext N t j=1 m ij . (2.8) (b) Update the synthetic tag line deformation function by fitting the control points, C, of the B-spline model f =argmin C N p summationdisplay i=1 N t summationdisplay j=1 m ij bardblv i ? f(r ij )bardbl 2 + ?bardblLfbardbl 2 (2.9) (c) Update the temperature and regularization parameter: T = T ? ?, ? = ?? ?. 23 where ? is the annealing rate. Note that the temperature controls the variance of the membership probability density in (2.3). Initially, this variance is relatively large, which allows many detected tag points to a?ect the position of a given synthetic tag line. As the annealing temperature cools, the variance gets smaller, and only detected tag points close to a synthetic tag line a?ect its position. A similar schedule is used with the regularization parameter, ?. Initially, when many detected tag points a?ect the position of a synthetic tag line, the deformation function is heavily smoothed to reduce the e?ect of detected tag points distal to the tag line. As the annealing temperature cools, the need for smoothing is reduced since the synthetic tag lines are only a?ected by detected points close to the tag line. After the deterministic annealing algorithm completes, the detected tag points are assigned a tag line index based on their maximum membership values. Since tag points are only detected in the neighborhood of pixel locations in the image, it is possible that a point on a tag line could have multiple detected tag points assigned to it. In this case, the detected tag point with the highest membership value is retained, and the remaining tag points are classified as false positives. The tag point classification algorithm can be applied to multiple slices in an imaging study. To ensure consistent tag line indices across multiple slices, the deformation function control points can be initialized to the control points from the previous slice, and a smaller initial temperature can be applied to those subsequent slices to make the classification and tracking process faster. If slices are misregistered due to a patient breath-hold shift, this shift will not cause a problem since the tag lines are applied in the same position relative to the slice prescription, which is typically consistent from slice to slice. If there is a breath- hold shift, the heart may shift relative to the slice, but the tags will be in the same position relative to the slice. 24 In many clinical studies, only the end-systolic strain in the LV wall is of interest, so only the end-diastolic (ED) and end-systolic (ES) time frames need to be analyzed. In short-axis slices, the classification algorithm is robust enough to inter-frame changes in tag line positions that it can classify detected tag points at ED and ES without tracking tags through intermediate time frames. The deformation control points at ED and ES are both initialized to zeros. In long-axis images, larger deformations can occur between ED and ES, and the increased size of the blood pool makes false positives more likely in the detected tag points. For these reasons, long-axis images are tracked in each time frame from ED to ES using the final control points from the previous time frame as the initial control points for the current frame, and a smaller initial temperature can be applied for the subsequent time frames. 2.3 Experiments and Results To perform a validation of the classification algorithm, 20 humans subjects were imaged with a grid-tagged imaging protocol with tag spacing of 7 pixels. The studies are identified in Section 2.5. Each study contained around 10 short-axis slices, one 2-chamber slice and one 4-chamber slice and 20 time frames were imaged during the entire cardiac cycle. Candidate tag points were detected in each image by pre-filtering the image with a composite Gabor filter and performing peak detection [48]. Tag points were detected and classified separately in each tag line direction in the grid. The 20 studies were randomly split into two groups of ten. One group, called development studies, were used to test and set parameters for the classification algorithm. The other group, called evaluation studies, were used to evaluate the performance of the algorithm. 25 2.3.1 Parameter Sensitivity The parameters of the classification algorithm were determined using the development studies. Considering the regions that the left ventricle occupies in long-axis images are usually wider than those in short-axis images, we tested the sensitivities to some parameters separately for short-axis and long-axis slices. The set of parameters includes the annealing rate ?, the initial temperatures T init , T SAinit2 (for subsequent short-axis slices) and T LAinit2 (for subsequent long-axis slices), the final temperature T final , the regularization weight ?, and the numbers of the B-spline control points N xSA , N ySA , N xLA and N yLA .Forthe grid tag pattern, since the initial tag lines are oriented perpendicularly, the numbers of the B-spline control points were assumed the same in the two orthogonal directions, i.e., N xSA = N ySA = nCP SA and N xLA = N yLA = nCP LA . Di?erent parameter values were applied for the parameter sensitivity test. The tracked tag lines were compared to tag lines tracked and edited by expert users. Only tag points inside the LV wall were used for classification accuracy analysis, where the classification accuracy was defined as the ratio of the number of correctly tracked tag lines / the total number of tag lines. The parameters were optimized one by one. The value that results in the maximum tag classification accuracy was chosen as the optimal value for a parameter. As the sensitivity of one parameter was tested with various values, the other parameters were set to empirically chosen values or the optimal values if they had already been obtained. The order of the parameters to be tested was ?, ?, T init , T final , T SAinit2 , T LAinit2 , and the numbers of control points. Tables 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7 and 2.8 show that the set of parameters ? =0.75, ? =10 3 , T init = 10, T final =0.001, T SAinit2 =4,T LAinit2 =4,nCP SA =6and nCP LA = 9 maximized the classification accuracy over the ten development studies. 26 Table 2.1: End-systolic tag point classification accuracies with various annealing rates (?) and empirical parameter values ? = 10, T init = 10, T final =0.5, T SAinit2 =3,T LAinit2 =2, nCP SA =6andnCP LA = 8 over the 10 human development imaging studies. ? classification accuracy 0.50 83.3% 0.55 82.4% 0.60 84.9% 0.65 88.8% 0.70 89.1% 0.75 91.1% 0.80 90.1% 0.85 88.4% 0.90 84.1% 0.95 77.0% Table 2.2: End-systolic tag point classification accuracies with various initial regularization weights (?), optimal annealing rate ? =0.75, and empirical parameter values T init = 10, T final =0.5, T SAinit2 =3,T LAinit2 =2,nCP SA =6andnCP LA =8overthe10human development imaging studies. ? classification accuracy 10 ?2 91.0% 10 ?1 91.0% 1 91.0% 10 91.1% 10 2 90.9% 10 3 91.9% 10 4 88.2% 10 5 86.2% 10 6 78.0% 27 Table 2.3: End-systolic tag point classification accuracies with various initial temperatures (T init ), optimal parameters ? =0.75 and ? =10 3 , and empirical parameter values T final = 0.5, T SAinit2 =3,T LAinit2 =2,nCP SA =6andnCP LA = 8 over the 10 human development imaging studies. T init classification accuracy 1 86.6% 5 89.1% 10 91.9% 10 2 42.9% 10 3 27.3% Table 2.4: End-systolic tag point classification accuracies with various final temperatures (T final ), optimal parameters ? =0.75, ? =10 3 and T init = 10, and empirical parameter val- ues T SAinit2 =3,T LAinit2 =2,nCP SA =6andnCP LA = 8 over the 10 human development imaging studies. T final classification accuracy 0.001 93.2% 0.005 90.9% 0.01 89.9% 0.05 89.8% 0.1 89.9% 0.5 91.9% 1 87.6% 28 Table 2.5: End-systolic tag point classification accuracies with various initial temperatures (T SAinit2 ) for the subsequent short-axis slices, optimal parameters ? =0.75, ? =10 3 , T init =10andT final =0.001, and empirical parameter values nCP SA =6overtheshort- axis slices of the 10 human development imaging studies. T SAinit2 classification accuracy 1 93.8% 2 93.7% 3 93.0% 495.9% 5 92.2% 6 89.8% 7 83.2% 8 80.3% 9 77.6% 10 67.5% Table 2.6: End-systolic tag point classification accuracies with various initial temperatures (T LAinit2 ) for the subsequent long-axis slices, optimal parameters ? =0.75, ? =10 3 , T init = 10 and T final =0.001, and empirical parameter values nCP LA = 8 over the long-axis slices of the 10 human development imaging studies. T LAinit2 classification accuracy 1 92.6% 2 93.8% 3 94.4% 4 94.7% 5 93.5% 6 94.1% 7 91.7% 8 93.2% 9 88.2% 10 87.9% 29 Table 2.7: End-systolic tag point classification accuracies with various numbers of control points (nCP SA ) for short-axis slices, optimal parameters ? =0.75, ? =10 3 , T init = 10, T final =0.001, and T SAinit2 = 4 over the short-axis slices of the 10 human development imaging studies. nCP SA classification accuracy 3 92.9% 4 92.8% 5 88.6% 6 95.9% 7 86.4% 8 85.6% 9 83.9% Table 2.8: End-systolic tag point classification accuracies with various numbers of control points (nCP LA ) for 2-chamber and 4-chamber slices, optimal parameters ? =0.75, ? =10 3 , T init = 10, T final =0.001, and T LAinit2 = 4 over the 2-chamber and 4-chamber slices of the 10 human development imaging studies. nCP LA classification accuracy 3 82.3% 4 87.0% 5 84.3% 6 92.3% 7 89.9% 8 94.7% 9 95.6% 30 2.3.2 Validation The optimal parameters obtained from the parameter sensitivity test were used in processing the evaluation studies. The synthetic tag line deformation model used quadratic B-splines with N xSA = N ySA = 6 for short-axis slices and N xLA = N yLA = 9 for long-axis slices. In the first processed slice and time frame, an initial temperature of T init =10was used. Since T determines the variance in the tag point position membership function in Equation (2.3), this value of T init is proportional to the tag spacing. In subsequent slices and time frames, a lower initial temperature, T SAinit2 =4andT LAinit2 =4forshortand long-axis slices, was used since the control points were initialized to the result of the previous slice or time frame. In all slices and time frames, T final =0.001, ? =10 3 and ? =0.75. Table 2.9: End-systolic tag point classification accuracies over the 10 human development and 10 human evaluation imaging studies. Development Evaluation classification accuracy 95.8% 98.2% Figure 2.1 shows the results of the deterministic annealing for a human study near end- diastole and near end-systole. The detected tag points are shown in yellow in Figure 2.1(a) and in blue in Figure 2.1(b), and the synthetic tag lines are shown in magenta. Contours of the left ventricle were displayed in each plot with white for reference purposes, but contours were not used at all in the detection or classification algorithms. Note that the synthetic tag lines are not intended to exactly match the detected tag points. Rather, their purpose is to act as cluster centers for tag point classification. Figure 2.2 shows the result of tag point classification in di?erent slices. Each detected tag point is displayed with a color indicating the tag line to which it belongs. Tag points classified as false positives are shown in blue. Note that tag points are classified according to their proximity to a given synthetic tag line. 31 Points between two synthetic tag lines have low maximum membership probabilities and are classified as false positives. In some cases, a connected set of detected points will change index in the middle of the line, but this primarily occurs in the blood pool, lung, or other regions of the image that are not tagged. The classification algorithm was implemented in MATLAB and took approximately 10 seconds for a short-axis slice, 50 seconds for a long-axis image, and 12 minutes of processing time for a single study on a 3.2GHz PC. To quantify the performance of the classification algorithm, the ten evaluation studies, which were not used for development or setting parameters, were processed. The results were compared to tag lines tracked and edited by an expert user. A tag point was considered classified correctly if its tag line index matched the index of the closest expert-edited tag line point. Only tag points inside the LV wall were included in the analysis. The classification accuracy (number of correctly classified tag lines/total number of tag lines) was 98.2%. The processing time for each study can be improved by applying di?erent annealing schedules. With the deformation model fixed, the computation time is determined by the initial and final temperatures and the cooling speed. An experiment with T init = 10, T final =0.5, T SAinit2 =3,T LAinit2 =2,? =10 3 and ? =0.75 over the ten evaluation studies showed the processing time was about 5 seconds for a short-axis image, 12 seconds for a long-axis image, and 4 minutes for a single imaging study on a 3.2GHz PC with a classification accuracy of 94.3%. In [56], we used T init = 10, T final =1,T SAinit2 =3, T LAinit2 =2,? =10 4 and ? =0.9 in the experimental validation, which resulted in a classification accuracy of 97.86%. For comparison purposes, we also implemented tag point classification by using the simulated annealing algorithm. Instead of optimizing the energy E(f,M) in Equation (2.6) with a deterministic schedule, we employed a stochastic process to update the control points 32 End-diastole End-systole Base Mid-ventricle Apex Figure 2.2: Tag line registration across multiple short-axis slices. Each detected tag line is displayed with a color indicating the tag line to which it belongs. Tag points classified as false positives are shown in blue. 33 of the B-spline deformation model. The same annealing schedule as in the deterministic annealing was used in the simulated annealing algorithm. Experiments on the 10 human evaluation imaging studies showed the classification accuracy was 69.0% and it took about 40 minutes to process a single study on a 3.8GHz PC. By using slower cooling rate, simulated annealing has the potential to get better classification accuracy. However, this will take a much longer time. Since the classification accuracy of the deterministic annealing algorithm is already 98.2%, there is not much room for improvement. It is more practical to keep the computation speed of deterministic annealing than sacrificing the speed for a less than 1.8% improvement. 2.4 Discussion In this chapter, we presented an algorithm for classifying a set of tag points detected in a tagged cardiac MR image into a set of tag lines. This algorithm allows tag lines to be tracked without smoothness constraints, which should be imposed when a 3-D deformation model is fitted to the tag data and not during tracking. A validation experiment on data not used for algorithm development resulted in a high classification rate when compared to expert-edited tag lines. The classification algorithm is robust enough to inter-frame changes in tag line position that, in short-axis images, end-systolic time frames can be processed without tracking tag lines through intermediate time frames, which can save the processing time. We explicitly tested the robustness of the parameter choices we made in the classifi- cation algorithm. A measure of robustness was demonstrated in our experimental results where high classification rates were obtained on both the development and evaluation data sets. 34 We also implemented the simulated annealing algorithm for comparison purposes. With the same annealing schedule used in the deterministic annealing, the simulated annealing resulted in a lower tag point classification accuracy. Additionally, it took a longer time to classify tag points with the simulated annealing schedule. We chose the deterministic annealing schedule in our algorithm because it is a controllable process and is faster. 35 2.5 Appendix Table 2.10: The 20 grid-tagged human imaging studies for validating the tag classification algorithm. Development Evaluation Study ID Visit Study ID Visit 05DAC041 BAS 05DAC005 M03 05DAC053 BAS 05DAC047 BAS 05DAC057 BAS 05DAC050 BAS 05P1MRR009 BAS 05P1VOL507 BAS 05P1MRR012 BAS 05P1VOL514 BAS 05P1VOL507 BAS 05P1VOL516 BAS 05P1VOL514 BAS 05P1VOL521 BAS 05P3CLN004 BAS 05P3CLN001 BAS 05P3CLN007 BAS 05P3CLN008 BAS 05P3CLN011 BAS 05P3CLN009 BAS 36 Chapter 3 Left Ventricular Motion Reconstruction 3.1 Introduction This chapter presents a new 3-D myocardial displacement and strain reconstruction method based on a?ne deformation and PSB deformation models. The domain more closely matches the LV?s shape than those of both Cartesian and cylindrical coordinate models. In this new deformation model, the spline basis functions are oriented radially, circumferentially, and longitudinally. This arrangement is consistent with the muscle fiber orientations in the LV wall, where the muscle fibers are oriented nominally in the circum- ferential direction with collagen struts [49] connecting the fibers in the cross-fiber direction. Also the deformation model includes the LV apex, which is di?cult to model in cylindrical coordinates. 3-D deformation and strain is computed independently for each imaged time frame using a three-step procedure. First, a backward fit is performed by fitting the deforma- tion model to measurements derived from tag lines identified in the tagged images. This deformation field maps each point in the deformed myocardium back to its reference po- sition, which we define as end-diastole. Next, a forward deformation field is computed by inverting the backward deformation field. This deformation maps each material point in the myocardium to its deformed position. In each of the first two steps, there is a two-stage fit where an a?ne deformation model is fit to the data, and a prolate spheroidal B-spline model is fit to the residuals. We call this deformation model the APSB model. Finally, Lagrangian strain is computed by analytically di?erentiating the forward deformation field. 37 This chapter is organized as follows. In Section 3.2, we briefly review the prolate spheroical coordinate system and continuum mechanics used in our algorithm. In Section 3.3, the closed-apex APSB deformation model is presented. In Sections 3.4 and 3.5, the closed-apex APSB model is validated and its performance is compared with those of a sim- plified model without the a?ne deformation called the PSB model and cylindrical models. The validation experiments were performed with an analytically defined LV deformation model and the tag displacement field tests presented in [10]. Finally, the results are dis- cussed in Section 3.6. 3.2 Background 3.2.1 The Prolate Spheroidal Coordinate System The prolate spheroidal coordinate system has been used in the cardiac LV motion modelling [16]. See Figure 3.1. The relationship between a point in Cartesian coordinates (x,y,z) and its corresponding prolate spheroidal coordinates (?,?,?)is: x = f sinh?sin?cos ? y = f sinh?sin?sin ? z = f cosh ?cos ?, (3.1) where ? ? [0,?), ? ? [0,?], ? ? [0,2?)andf is the focal radius. The relationship between a vector (u x ,u y ,u z ) in Cartesian coordinates and its counterpart (u ? ,u ? ,u ? )inprolate spheroidal coordinates is given by: u x = cosh ?sin ?cos ? radicalBig sinh 2 ? +sin 2 ? u ? + sinh?cos ?cos ? radicalBig sinh 2 ? +sin 2 ? u ? ? sin ? ? u ? 38 z x y ? ? n? n? n? foci Figure 3.1: The prolate spheroidal coordinate system and Cartesian coordinate system. n ? , n ? , n ? are the orthogonal unit vectors in ? (radial), ? (longitudinal), ? (circumferential) directions respectively. 39 u y = cosh ?sin ?sin ? radicalBig sinh 2 ? + sin 2 ? u ? + sinh?cos ?sin ? radicalBig sinh 2 ? +sin 2 ? u ? +cos? ? u ? u z = sinh?cos ? radicalBig sinh 2 ? +sin 2 ? u ? ? cosh ?sin ? radicalBig sinh 2 ? +sin 2 ? u ? , (3.2) where ?, ?,and? are the coordinates of the vector origin. Then the 3-D Laplacian of a function u(?,?,?) in prolate spheroidal coordinates is given by triangleinv 2 u = 1 sinh 2 ? +sin 2 ? [(csch 2 ? +csc 2 ?) ? 2 u ?? 2 +cot? ?u ?? + ? 2 u ?? 2 +coth? ?u ?? + ? 2 u ?? 2 ]. (3.3) 3.2.2 Continuum Mechanics The motion of a material point, which is fixed on the object, is described by a smooth, reversible function s = s(p,t) (3.4) that maps a material point p =[X,Y,Z] T to a spatial point s =[x,y,z] T , which is fixed in space. The inverse of this function is called a reference map. Material and spatial points are equal in the reference configuration, which we define as end-diastole. Let u p (p) denote the displacement of the material point p and u s (s) denote the dis- placement of the corresponding spatial point s. Then, s = p + u p (p) p = s + u s (s) (3.5) 40 and ds = Fdp =(I + triangleinvu p )dp (3.6) where F is the material coordinate deformation gradient tensor, F = ? ? ? ? ? ? ? ? ?x ?X ?x ?Y ?x ?Z ?y ?X ?y ?Y ?y ?Z ?z ?X ?z ?Y ?z ?Z ? ? ? ? ? ? ? ? , (3.7) triangleinvu p is the material coordinate displacement gradient tensor and I is the identity matrix. The Lagrangian strain tensor, E, is defined as E = 1 2 parenleftBig F T F?I parenrightBig , (3.8) where T denotes transpose. 3.3 Algorithm Development In this section, we develop the APSB algorithm for reconstructing 3-D deformation and strain from a set of deformed tag line points at a single phase in the cardiac cycle. This process can be repeated for each imaged cardiac phase to obtain a sequence of deformations and strains. We assume that the tag lines have been identified in each image in the study by one of the several published techniques [51, 12, 13, 48, 56]. The deformation reconstruction algorithm reconstructs 3-D displacement and strain from the 1-D displacement measurements obtained from the tag tracking procedure (see Figure 1.4). The algorithm is implemented in three steps. First, a 3-D spatial coordinate displacement field is computed to map points in the deformed myocardium back to their undeformed positions. Second, this spatial coordinate displacement field is converted to 41 a material displacement field, which maps material points in the undeformed myocardium to their deformed positions. Finally, strain is computed by analytically di?erentiating the material coordinate B-spline displacement field with respect to space and using Equation (3.8). Since cardiac motion includes not only contraction and relaxation but also bulk shear, bulk rotation and bulk displacement, the deformation reconstruction algorithm uses a two- stage fit similar to that used in [16] where an a?ne deformation model is first fit to the data, and a prolate spheroidal B-spline model is fit to the residuals. In each prolate spheroidal B-spline model fitting, the origin and focal radius of prolate spheroidal coordinate system are determined as follows. The z-axis of a Cartesian coordinate system is defined along the LV longitudinal direction from the apex to the base with the x-y plane parallel to the short-axis image planes. The origin is defined to be one-third of the distance from the base to the apex at the (x,y) centroid of the tag points [50]. The focal radius is chosen to make the ?-coordinate of the LV wall as close as possible to one [50] in the mean square sense, i.e., f =argmin f summationdisplay k ||1 ? ? k || 2 , (3.9) where ? k is the ?-coordinate of the kth tag point and the sum is over all tag points. 3.3.1 Spatial Displacement Field Fitting In this step, a spatial displacement field is computed that maps points in the deformed myocardium back to their undeformed positions. At a spatial point s, we define its motion as p =(A s + I)s + d s + R(s)u rs (s) (3.10) 42 where s is a point in Cartesian coordinates, and A s is a 3 ? 3 a?ne deformation matrix, which includes bulk rotation, linear shears, and linear contraction/expansion. d s is a 3?1 bulk displacement vector. u rs (s) is a non-a?ne displacement described by the PSB model. R(s) is the rotation matrix defined by Equation (3.2) that transforms a vector in prolate spheroidal coordinates into Cartesian coordinates. The 3-D displacement of a spatial point s,isthengivenby u s (s)=p?s = A s ?s + d s + R(s)u rs (s). (3.11) The 1-D displacement measurement, d(s), obtained from a single deformed tag line point is (see Figure 1.4) d(s)=n(s) ?u s (s) = n(s) ? [A s s + d s + R(s)u rs (s)] , (3.12) where n(s) is a vector normal to the undeformed tag plane from which the spatial point s originated. To fit the a?ne deformation model, we set the non-a?ne deformation to zero, in which case Equation (3.12) becomes d(s)=n(s) ? [A s ?s + d s ]. (3.13) 43 The displacement measurements for all tag line points in all slices and all tag line orienta- tions in a single time frame can be stacked into the matrix equation D = NSA ds , (3.14) where D is a Q ? 1 vector of 1-D displacement measurements, N is a Q ? 3Q matrix of normal vectors, S is a 3Q?12 sparse matrix of the spatial point positions, and A ds is a 12?1 vector containing the components of a?ne deformation matrix A s and bulk displacement vector d s . Equation (3.13) is solved for A s and d s via linear least squares. Note that both the a?ne and non-a?ne displacement fields could be fit at the same time. In experiments using mathematical models and in vivo MR image data, however, the two-stage fit was slightly faster and more accurate. In the second stage of spatial fitting, the prolate spheroidal B-spline model is fit to the displacement residuals, which are defined as d r (s)=n(s)[R(s)u s (s) ? (A s s + d s )] = n(s) ? [R(s)u rs (s)] =[R T (s)n(s)] ?u rs (s) = n ? (s)u ? (s)+n ? (s)u ? (s)+n ? (s)u ? (s), (3.15) or, in matrix form, D r = D?NSA ds = NU rs , (3.16) 44 where D r is a Q?1 vector of 1-D displacement measurement residuals, and U rs is a 3Q?1 vector of 3-D spatial coordinate displacement residuals. We model the 3-D non-a?ne displacements with a 3-D B-spline, which is the tensor product of three 1-D B-splines. Let u rs =[u ? (l,p,t),u ? (l,p,t),u ? (l,p,t)] T denote the non- a?ne displacement vector, where u ? (l,p,t)= L?1 summationdisplay i=0 P?1 summationdisplay j=0 T?1 summationdisplay k=0 B o i (l)B o j (p)B c k (t)c ? ijk u ? (l,p,t)= L?1 summationdisplay i=0 P?1 summationdisplay j=0 T?1 summationdisplay k=0 B o i (l)B o j (p)B c k (t)c ? ijk u ? (l,p,t)= L?1 summationdisplay i=0 P?1 summationdisplay j=0 T?1 summationdisplay k=0 B o i (l)B o j (p)B c k (t)c ? ijk . (3.17) B o i (?)isanopenB-spline,B c k (?) is an closed B-spline, c ijk is a control point, and L,P,T are the number of control points in the radial (?), longitudinal (?), and circumferential (?) directions respectively. Note that we use open B-splines in the radial and longitudinal directions and closed B-splines in the circumferential direction. In the radial direction, the knots are uniformly spaced between the minimum and maximum radial coordinates (? min and ? max )ofallde- formed tag points. In the longitudinal direction, the knots are uniformly spaced between the most basal tag line points and the apex (? min and ? max = ?). In the circumferen- tial direction, the knots are uniformly spaced between 0 and 2?. The prolate spheroidal coordinates ?, ? and ? are related to the B-spline parametric coordinates l, p and t by l = ?? ? min ? max ? ? min ? (L ? O +1) p = ?? ? min ? max ? ? min ?(P ? O +1) t = ? 2? ? T (3.18) 45 where O is the spline order. To relate the control points to the 1-D displacement measurements, we evaluate the B-spline displacement field at the same positions as the deformed tag line points and stack them in matrix form to get U rs = B s C s (3.19) where U rs is the vector of 3-D spatial coordinate displacement residuals in (3.16), C s is a 3LPT ? 1 vector containing the control points for all three displacement components, and B s is a 3Q ? 3LPT matrix of B-spline basis functions evaluated at the tag point positions using (3.18). Substituting (3.19) into (3.16) yields D r = NB s C s . (3.20) To ensure that the non-a?ne deformation is smooth, we impose Laplacian regulariza- tion on each of Cartesian displacement components u x (?,?,?), u y (?,?,?)andu z (?,?,?), which are obtained via rotating the prolate spheroidal displacements by (3.2). This regular- ization is imposed at the knot positions in the B-spline model with the exception of those at the apex. Regularization is not imposed at the apex because the Laplacian in prolate spheroidal coordinates (Equation (3.3)) has a singularity there. Since there are typically more tag line points than control points, (3.20) is over- determined. We solve it using least squares (3.21) while penalizing the Laplacian of each displacement component. We also enforce three equality constraints at the apex to ensure a smooth deformation across each point where the circumferential closed splines shrink to one dot and the longitudinal open splines meet. We constrain the circumferential control points to be equal at the apex, and we constrain the first and second order derivatives of each pair of antiparallel longitudinal displacement components to be equal at the apex. 46 These constraints can be written in matrix form as E s C s = 0,where0 is a vector of zeros. The optimal control points are therefore given by C s =arg min C s E s C s = 0 bardblD r ?NB s C s bardbl 2 + ?bardblL s C s bardbl 2 (3.21) where ? is a weighting parameter and L s C s is the Laplacian at each knot location (except the apex) in matrix form. To solve this constrained linear least squares fitting problem, we used the MATLAB command lsqlin, which uses an an active set quadratic programming method [52]. 3.3.2 Material Displacement Field Fitting While the spatial coordinate field can provide useful information about cardiac function, Lagrangian strain is usually desired for cardiac motion analysis. To compute Lagrangian strain, the spatial coordinate displacement field computed above must be converted to a material coordinate displacement field. To perform this conversion, we use the spatial displacement field to compute measurements of the material coordinate displacement field, u p (p). For a given spatial point, s, its corresponding material point is given by p = s+u s (s). So for a given set of spatial points, s, we compute a set of material coordinate displacement field measurements using Equation (3.5), i.e. u p (p)=s?p. (3.22) 47 Any set of spatial points could be used, but we used the deformed tag line points that were used to compute the spatial coordinate displacement field because we expect the B-spline deformation field to be most accurate at these points. The material coordinate a?ne and prolate spheroidal B-spline displacement fields are computed from these measurements in a manner similar to that used for the corresponding spatial displacement fields. The only di?erence is that the material displacement field measurements are 3-D, so the tag plane normals are not required. 3.3.3 Lagrangian Strain Computation Lagrangian strain is computed by analytically di?erentiating the material coordinate displacement field with respect to space to obtain the deformation gradient, F. Since none of the independent coordinates ?, ? and ? is a distance, it is easier to compute spatial derivatives in Cartesian coordinates. In Cartesian coordinates, a material point (X,Y,Z) T translating and deforming to a spatial position (x,y,z) T maintains the relationship ? ? ? ? ? ? ? ? x y z ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? X Y Z ? ? ? ? ? ? ? ? + A p ? ? ? ? ? ? ? ? X Y Z ? ? ? ? ? ? ? ? + d p + ? ? ? ? ? ? ? ? u X (X,Y,Z) u Y (X,Y,Z) u Z (X,Y,Z) ? ? ? ? ? ? ? ? (3.23) where u. denotes the material displacement field computed with the prolate spheroidal B- spline model (3.17) in Cartesian coordinates. The conversion between prolate spheroidal displacements and Cartesian displacements is given in (3.2). The Cartesian coordinate 48 deformation gradient tensor is then F = I + A + ? ? ? ? ? ? ? ? ?u X (?,?,?) ?X ?u X (?,?,?) ?Y ?u X (?,?,?) ?Z ?u Y (?,?,?) ?X ?u Y (?,?,?) ?Y ?u Y (?,?,?) ?Z ?u Z (?,?,?) ?X ?u Z (?,?,?) ?Y ?u Z (?,?,?) ?Z ? ? ? ? ? ? ? ? . (3.24) The partial derivatives with respect to material coordinates X,Y ,Z in (3.24) are computed by chain rule ? ?X = ?? ?X ? ?? + ?? ?X ? ?? + ?? ?X ? ?? ? ?Y = ?? ?Y ? ?? + ?? ?Y ? ?? + ?? ?Y ? ?? (3.25) ? ?Z = ?? ?Z ? ?? + ?? ?Z ? ?? + ?? ?Z ? ?? , where ?? ?X = 1 radicalbig (r 1 + r 2 ) 2 ? 4f 2 parenleftbigg X r 1 + X r 2 parenrightbigg ?? ?Y = 1 radicalbig (r 1 + r 2 ) 2 ? 4f 2 parenleftbigg Y r 1 + Y r 2 parenrightbigg ?? ?Z = 1 radicalbig (r 1 + r 2 ) 2 ? 4f 2 parenleftbigg Z + f r 1 + Z ? f r 2 parenrightbigg ?? ?X = 1 radicalbig 4f 2 ? (r 1 ? r 2 ) 2 parenleftbigg X r 2 ? X r 1 parenrightbigg ?? ?Y = 1 radicalbig 4f 2 ? (r 1 ? r 2 ) 2 parenleftbigg Y r 2 ? Y r 1 parenrightbigg ?? ?Z = 1 radicalbig 4f 2 ? (r 1 ? r 2 ) 2 parenleftbigg Z ? f r 2 ? Z + f r 1 parenrightbigg ?? ?X = ? Y X 2 + Y 2 ?? ?Y = ? X X 2 + Y 2 ?? ?Z =0, 49 and r 1 = radicalbig X 2 + Y 2 +(Z + f) 2 , r 2 = radicalbig X 2 + Y 2 +(Z ? f) 2 , ?= cosh ?1 r 1 +r 2 2f ? > 0, ?= cos ?1 r1?r2 2f 0 ? ? ? ?, ?= tan ?1 Y X 0 ? ? < 2?. The strain tensor is then computed using Equation (3.8) and rotated to a prolate spheroidal coordinate system. This procedure results in three normal strains ? radial thickening (E rr ), circumferential shortening (E cc ), and longitudinal shorting (E ll )?and three shear strains (E rc , E rl and E cl ). At the apex, the radial direction is uniquely defined in prolate spheroidal coordinates, but the circumferential and longitudinal directions are not. Typically strains are computed on an N r ? N c ? N l mesh of material points, so at the apex, for each radial point, N c prolate-spheroidal strain tensors can be computed. These tensors are averaged to obtain a single strain tensor for each radial point at the apex. 3.3.4 Parameters The user-defined parameters of the prolate spheroidal B-spline deformation model are the orders of the B-splines, numbers of control points in the three directions (L, P,andT in the radial, longitudinal, and circumferential directions, respectively), and the regularization weighting parameter ?. The same parameters are used for both the spatial and material coordinate displacement fields. The parameters L, P, T,and? as well as the orders of the B-spline used in the APSB model were chosen to minimize the reconstructed strain errors 50 when the APSB model was fit to simulated tag line data from the analytically-defined LV deformation model in [16]. In this chapter, quadratic B-splines in all three directions were used in all reconstruc- tions. For comparison purposes, we implemented a closed-apex PSB reconstruction where the prolate spheroidal B-spline deformation model directly fits to the data in each of the spatial and material estimations. The cylindrical B-spline deformation reconstruction [18] and its regularized version were also implemented. In the regularized cylindrical B-spline model, Laplacian regularization is imposed on each component of the 3-D displacement as in Section 3.3.1 to smooth the reconstructed displacement field. Table 3.1 shows the optimal parameters for each method. Table 3.1: Parameters for the quadratic B-spline reconstructions. regularization weight ?nCP ? nCP ? nCP ? APSB 0.019 3 5 8 PSB 0.010 3 5 8 regularization weight ?nCP r nCP z nCP ? Reg.Cyl.B. 0.008 3 4 10 Cyl.B. ? 3 4 10 3.4 Experimental Methods Two sets of experiments were performed to validate the APSB model. The first set of experiments studied the sensitivity of the model to perturbations of the coordinate system center. Details of this experiment are given below. The second set of experiments compared the APSB, PSB and cylindrical B-spline models using the validation framework described in [10]. This framework consists of a sequence of tests on in vivo human and canine data. 51 Center Sensitivity Test The centers of the prolate spheroidal coordinate systems used in the APSB and PSB algorithms are automatically computed from the data. We studied the sensitivity of the reconstructed strains to this center point with Monte Carlo simulations. First, a set of tag line data was simulated using the analytical LV deformation model of [16] with the following parameters: 7mm tag spacing; 6 short-axis imaging planes with 10-mm slice thickness; and 6 long-axis imaging planes with 30 o angular spacing. About 4000 tag points are generated with these parameters. The noise-free coordinate system origin was computed from this tag line data on a 3 radial ? 12 circumferential ? 20 longitudinal cylindrical mesh of material points using the procedure described in Section 3.3. In each trial, 1-D, zero-mean Gaussian noise was then added to each of the x, y and z coordinates of the noise-free origin and strain was reconstructed from the simulated data using the perturbed origin. The focal point of the prolate spheroidal coordinate systems was computed from the perturbed origin and tag line data for each trial. The strain error was defined as the root-mean-square (RMS) di?erence between the noise-free strain and noisy strain computed over all points in the material point mesh. 3.5 Results 3.5.1 Center Sensitivity Test Figure 3.2 shows the normal strain error standard deviation (SD) versus the center per- turbation SD for the APSB, PSB and cylindrical B-spline models. Statistics were computed over 400 trials. The di?erence between statistics computed over 300 trials and 400 trials was negligible. The noise-free RMS values of the E rr , E cc and E ll were 0.381, 0.131 and 0.295 respectively. The circumferential (E cc ) and longitudinal (E ll ) strains were relatively 52 0 1 2 3 4 5 0 0.05 0.1 0.15 0.2 0.25 0.3 error in E rr 0 1 2 3 4 5 0 0.05 0.1 0.15 0.2 0.25 0.3 standard deviation of noise in each dimension (mm) error in E cc APSB PSB Reg. Cyl. B. Cyl. B. 0 1 2 3 4 5 0 0.05 0.1 0.15 0.2 0.25 0.3 error in E ll Figure 3.2: Strain error standard deviation (SD) versus center point perturbation SD with the quadratic B-spline reconstructions on the analytically defined mathematical LV defor- mation model. robust to the perturbation of origin. The radial (E rr ) strain, however, was more sensitive. The APSB and PSB reconstructions were more sensitive to the center position than those from the cylindrical models. This increase in sensitivity is primarily because the APSB and PSB models must compute a focal point from the origin and tag line data, whereas this parameter is not required in the cylindrical model. 3.5.2 Validation on In Vivo Human and Canine Data The APSB, PSB and cylindrical B-spline models were validated and compared using the framework in [10]. First, strains were reconstructed from a high resolution tagged cardiac imaging study of a normal human volunteer. Next, the sensitivity of each algorithm to perturbations of the tag line data was studied using Monte Carlo simulation. Then, the spatial resolution of each algorithm was computed. Finally, all four deformation models were fit to tag line data from a pathology database. Details for each of these experiments can be found in [10]. The round robin test in [10] was not performed since, as mentioned in [10], this test is not very discriminating. 53 3.5.3 Strain Maps on a Normal Heart In this experiment, Lagrangian strain was reconstructed from a high-resolution normal human data set on an 18 circumferential ? 11 longitudinal cylindrical mesh located at the midwall. The computation time of APSB was about 12 seconds per time frame on a Pentium(R)4 3.2GHz PC. Figures 3.3, 3.4 and 3.5 show the strain maps computed with the APSB, PSB and cylindrical methods. For clarity, only one point out of two is displayed in both the circumferential and longitudinal directions. Each plot represents the evolution over time of the strain, starting at end-diastole, covering the systole and ending after the rapid filling of diastole. The arrangement of the plots represents a map of the myocardial mid-wall, from base to apex and around the septal, anterior, lateral and then posterior walls. The time evolution of strains computed with all B-spline algorithms were generally con- sistent with the results reported in [10]. The strains increase in magnitude through systole and then decrease in magnitude during diastole. The temporal evolution of circumferen- tial shortening (E cc ) and longitudinal shortening (E ll ) was smoother than radial thickening (E rr ). The higher variability in E rr was caused by the relatively sparse sampling of my- ocardial motion in the radial direction. In all of the APSB, PSB and cylindrical models, strain is reconstructed independently for each time frame, so the temporal strain evolution was not as smooth as the four-dimensional models in [10], which use temporal smoothing. Of the four models compared in [10], the temporal strain evolution with the APSB and PSB models was closest to the TEA method, which also uses a prolate spheroidal coordi- nate system. The strains computed from the APSB, PSB and cylindrical B-spline models were similar except near the apex, where the prolate spheroidal and cylindrical coordinate systems are most dissimilar. The spatial strain evolution from the APSB and PSB models, 54 septum anterio r latera l posterio r base ?0.3 ?0.2 ?0.1 0 apex APSB PSB Reg.Cyl.B. Cyl.B. Figure 3.3: Circumferential shortening (E cc ) strain map at the mid-wall for the high reso- lution normal human study over t =[41,481] ms. 55 septum anterio r latera l posterio r base ?0.2 ?0.1 0 0.1 apex APSB PSB Reg.Cyl.B. Cyl.B. Figure 3.4: Longitudinal shortening (E ll ) strain map at the mid-wall for the high resolution normal human study over t =[41,481] ms. 56 ?0.1 0 0.1 0.2 septum anterior lateral posterior base apex APSB PSB Reg.Cyl.B. Cyl.B. Figure 3.5: Radial thickening (E rr ) strain map at the mid-wall for the high resolution normal human study over t =[41,481] ms. 57 which used explicit regularization, was relatively smooth. However, although explicit regu- larization was also used in the regularized version of the cylindrical model, there was trivial di?erence in the strainmaps of the both cylindrical models. The strains in the cylindrical B-spline models had more artifacts. The 3-D end-systolic strains of the high resolution normal volunteer study are also shown in Figures 3.6 and 3.7 on the 3 radial ? 12 circumferential ? 11 longitudinal cylin- drical mesh. The normal volunteer exhibits fairly uniform radial thickening (yellow) and circumferential and longitudinal shortening (blue) in the 3-D strain maps reconstructed with the APSB and PSB. The cylindrical B-spline reconstructions of the normal volun- teer show nonphysiological radial thickening reduction in the normal sub-endocardium and sub-epicardium, expressed as negative artifacts in the radial thickening strain map. In the circumferential and longitudinal strain maps, the cylindrical B-spline reconstruction has inconsistently distributed strains close to the apex. These strain map artifacts are due to the mismatch between the cylindrical domain and the shape of the left ventricle. In contrast to the cylindrical B-spline models, the APSB and PSB models can compute strains at the apex. Figure 3.8 shows end-systolic strains on a 3 radial ? 24 circumferential ? 19 longitudinal mesh that includes the apex. The circumferential strain increases near the apex, whereas the radial and longitudinal strain stay relatively constant. 3.5.4 Noise Immunity Tag line detection and tracking algorithms have been shown to be accurate [53] with a precision inversely proportional to the tag line contrast-to-noise ratio (CNR) [54]. The purpose of this experiment was to access the e?ect of tag detection uncertainty on the reconstructed strains. We studied the sensitivity of spatial noise of each reconstruction with a Monte Carlo simulation as described in [17, 10]. 58 Figure 3.6: End-systolic strain maps of radial thickening (E rr , left), circumferential short- ening (E cc , middle) and longitudinal shortening (E ll , right) computed with the APSB (top) and PSB (bottom) models on an open-apex mesh. The red balls in the plots denote the LV mid-septum. 59 Figure 3.7: End-systolic strain maps of radial thickening (E rr , left), circumferential short- ening (E cc , middle) and longitudinal shortening (E ll , right) computed with the cylindrical B-spline model (bottom) and its regularized version (top) on an open-apex mesh. The red balls in the plots denote the LV mid-septum. 60 Figure 3.8: End-systolic strain maps of radial thickening (E rr , left), circumferential short- ening (E cc , middle) and longitudinal shortening (E ll , right) computed with the APSB (top) and PSB (bottom) models on a closed-apex mesh. The red balls in the plots denote the LV mid-septum. 61 The end-diastole and end-systole time frames were used in the test. The tag points at end-diastole and end-systole were perturbed with 1-D Gaussian noise with a chosen standard deviation ? perpendicular to the direction of the tag plane normal. 200 trials were generated for each noise ? of 0.1, 0.2, 0.3 and 0.4mm. Then, the displacement field was reconstructed for the 200 trials and the strain tensors were computed on the midwall mesh. The mean and standard deviation of the di?erence, which is the noisy strain tensor subtracted to the noise-free strain tensor, were computed over all the mesh points for all the trials to evaluate the sensitivity to spatial noise in the displacement field reconstruction. 0 0.1 0.2 0.3 0.4 0.5 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 0.02 0.03 0.04 Error in E rr 0 0.1 0.2 0.3 0.4 0.5 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 0.02 0.03 0.04 Error in E rc 0 0.1 0.2 0.3 0.4 0.5 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 0.02 0.03 0.04 Error in E rl 0 0.1 0.2 0.3 0.4 0.5 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 0.02 0.03 0.04 Error in E cc 0 0.1 0.2 0.3 0.4 0.5 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 0.02 0.03 0.04 Error in E cl Standard deviation of 1D noise(mm) 0 0.1 0.2 0.3 0.4 0.5 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 0.02 0.03 0.04 Error in E ll APSB PSB Reg. Cyl. B. Cyl. B. Figure 3.9: The mean and standard deviation on error on Lagrangian strain E in a local coordinate system (radial, circumferential or longitudinal). The center of each errorbar is on the mean value and the bar extends to one standard deviation above and below the mean. 62 Figure 3.9 shows the RMS error on the Lagrangian strain in local coordinates at dif- ferent noise levels. Except the error in E rr , the typical strain errors are below 0.02 at end-systole for these data. Cylindrical B-spline reconstructions show an E rl error of about 0.025. Only E rr shows a significant bias and a high sensitivity, which can be due to the tag data are more sparsely distributed in radial direction and the Lagrangian strain tensor is intrinsically biased when some noise is added to the data [17]. In this test, except for E rr , the biases for the other strains are subtle, and the errors are similar to these reported in [10]. The APSB and PSB reconstructions have smaller RMS strain errors for each noise level, so they are less sensitive to noise than the cylindrical B-spline reconstructions. The APSB reconstruction is the least sensitive. 3.5.5 Spatial Resolution A displacement cloud experiment [17, 10] was performed to measure the spatial resolu- tion of the 3-D motion reconstruction algorithms. A displacement cloud is a displacement field that is zero except for a localized, cloud-shaped motion abnormality. The smallest cloud that can be reconstructed with a certain accuracy is a measure of spatial resolution. A series of 3-D clouds were generated. Each cloud consisted of displacement vectors oriented in a particular direction and centered at a mesh point. The displacement vectors of the neighboring points of the cloud center were defined with a Gaussian roll-o? from the central point. One cloud experiment was generated for each point of the mid-wall mesh described in the previous section for each of the radial, circumferential and longitudinal directions and for a sequence of cloud widths from 4 to 20mm in increments of 1 or 2mm. For each cloud deformation, a displacement field was reconstructed and the displacement response of the cloud center was computed. The displacement response was defined as the magnitude of 63 the reconstructed cloud center displacement divided by the magnitude of the original center displacement and multiplying by 100%. 0 0.05 0.1 0.15 0.2 0.25 0.3 0 10 20 30 40 50 60 70 80 90 100 % displacement respons e Direction R 0 0.1 0.2 0.3 0 10 20 30 40 50 60 70 80 90 100 1/cloud width (1/mm) Direction C APSB PSB Reg. Cyl. B. Cyl. B. 0 0.05 0.1 0.15 0.2 0.25 0.3 0 10 20 30 40 50 60 70 80 90 100 Direction L cloud direction Methods R C L APSB 9.2 9.2 9.7 PSB 9.1 9.2 9.7 Reg.Cyl.B. 10.0 10.1 8.0 Cyl.B. 10.0 10.1 8.0 Figure 3.10: Spatial resolution for each displacement direction (radial, circumferential or longitudinal). Top: for di?erent cloud widths, the per cent of displacement response for each method; bottom, the corresponding values in mm for a response of 50%. Figure 3.10 shows the percentage displacement responses for di?erent cloud widths in the radial, circumferential and longitudinal directions, where the cloud width is the diameter of the cloud at one-half of its maximum magnitude. The displacement response becomes greater as the cloud width increases. If the spatial resolution is defined as the cloud width for which the response is 50%, it is ? 9mm for the APSB and PSB reconstruction, and is ? 10mm for the cylindrical B-spline reconstructions. 64 3.5.6 Pathology Database In this experiment, the performance of the four models was studied on di?erent types of cardiac motions. Each model was used to reconstruct a displacement field for six groups of subjects, including five normal humans volunteers, five patients with infarcted myocardium, four patients with dilated cardiomyopathy, five canines under debutamine infusion (10 ug Kg ?1 min ?1 ), five canines with right atrial pacing, and five canines with right ventricular pacing. For each method and each subject, the displacement field was reconstructed. The accuracy of the fit was evaluated by computing the signed distance error between actual tag line points and end-systolic simulated tag line points [10]. Simulated tag points were generated by deforming simulated tag planes by the reconstructed displacement field and computing points along the intersection between the deformed simulated tag planes and the image planes. Figure 3.11 shows the distance errors for di?erent types of cardiac motions. Each set of error bar represents one patient which is numbered from 1 to 5. As in the previous plots, the center of the error bar is on the mean value and the extension of the error bar is ?1 SD. The experimental results indicate that all methods can reconstruct the cardiac motion with similar tag distance errors. The mean error is generally around zero, showing the estimation of the motion is unbiased. The magnitude of the distance error for the normal human data is generally higher than those for the other human data and canine data. The reason is that the normal human heart has greater strains. The tagged distance errors are mostly in the 0.3mm to 0.6mm range, which is similar to those reported in [10]. In two of the human normal and two of the human infarct studies, however, the PSB model had 65 0 1 2 3 4 5 6 ?1 ?0.5 0 0.5 1 HUMAN NORMAL Tag distance error (mm) 0 1 2 3 4 5 6 ?1 ?0.5 0 0.5 1 HUMAN INFARCT Tag distance error (mm) 0 1 2 3 4 5 6 ?1 ?0.5 0 0.5 1 HUMAN DCM Tag distance error (mm) 0 1 2 3 4 5 6 ?1 ?0.5 0 0.5 1 CANINE DOBUTAMINE Tag distance error (mm) 0 1 2 3 4 5 6 ?1 ?0.5 0 0.5 1 CANINE RA PACED Tag distance error (mm) ?1 ?0.5 0 0.5 1 CANINE RV PACED Tag distance error (mm) 0 1 2 3 4 5 6 APSB PSB Reg. Cyl. B. Cyl. B. Figure 3.11: RMS error on tag distances for di?erent types of motion: human normal, infarcted or dilated cardiomyopathy, canine debutamine, RA or RV paced. Each set of error bars represents one patient which is numbered from 1 to 5. 66 larger errors (0.7mm to 1.0mm). The cylindrical models both had a higher than normal error (1.1mm) for a single normal study. 3.6 Discussion In this chapter, we presented a new method for reconstructing 3-D myocardial defor- mation and strain from planar tagged MR images based on a prolate spheroidal B-spline model and a?ne deformation. In each of spatial and material estimation, the algorithm uses a two-stage fit where the a?ne deformation model is fit to the data, and the prolate spheroidal B-spline model is fit to the residuals (APSB). The strains are then computed analytically from the reconstructed material displacement field for any material point in the myocardium including the apex. The PSB model we implemented for comparison purpose is a simplified version of the APSB model. It fits to the data directly with the 3-D prolate spheroidal B-splines. The APSB model, however, could be extended to 4-D by adding a B-spline in the time direction in the tensor-product model (3.17). An open B-spline would be used if tag line data is only available from end-diastole to end-systole, whereas a closed B-spline could be used if tag line data is available throughout the cardiac cycle. Compared to previously published methods for reconstructing myocardial deformation from tagged MR images, our APSB and PSB models construct 3-D myocardial deformation at each imaged cardiac phase similar to the finite element method of Young et al [13], the prolate spheroidal model of O?Dell et al [16], the discrete model of Denney and McVeigh [17], and the cylindrical B-spline model of Deng and Denney [18]. Unlike the 4-D (3-D+t) models of Huang et al[19] and ? Ozt?urk and McVeigh [20], the prolate spheroidal B-spline model can not reconstruct myocardial strain between imaged cardiac phases, but it can reconstruct the deformation at a single time frame without having to process the entire data set. This is particularly useful in studies where only end-systolic strain is measured. 67 We validated the APSB and PSB methods using the same data and tests used in [10], so comparisons can be made with the four techniques in that paper. The APSB model had less noise sensitivity and accurately reconstructed all the end-systolic deformations in the pathology database. It performed better than PSB. When comparing the time-evolution of strains in a normal human LV, the APSB and PSB methods are closest to the TEA technique [16], which also uses a deformation model defined in prolate-spheroidal coordinates. The noise immunities of the APSB and PSB algorithms are similar to the TTT method [20], but the APSB and PSB reconstructed radial strains with better accuracy and precision. This improvement in radial strain sensitivity is probably due to the use of a prolate-spheroidal coordinates in the APSB and PSB versus the Cartesian coordinate system used in the TTT. The APSB and PSB methods were also more immune to noise than the TEA method. This di?erence is probably due to the use of B-spline basis functions in the APSB and PSB models, which are spatially localized, versus the globally-defined spherical harmonics used in TEA. Spatially localized basis functions limit the e?ect of reconstruction artifacts whereas globally-defined basis functions can spread this e?ect over the entire model. The spatial resolution of the APSB method was similar to that of the TEA and TTT methods, except in the longitudinal direction where the prolate spheroidal methods (APSB, PSB and TEA) have a significantly higher spatial resolution. When reconstructing deformation from the pathology database, the tag distance errors in APSB model were closest to the TEA method and demonstrated the same general variation in performance with pathology as observed in [10]. The APSB method was also compared with the cylindrical B-spline method and its regularized version using the data and tests from [10]. The di?erence between the two cylindrical models were subtle in each test. In the time-evolution of strains in a normal hu- man LV, the APSB, PSB and cylindrical strains showed the same overall temporal trends, 68 but the APSB and PSB strains were spatially smooth compared to the cylindrical method, which showed more temporal and spatial variation. All of the APSB, PSB and cylindrical methods had a noise sensitivity with acceptable limits (below 10% of the average end- systolic strain [10], but the APSB and PSB methods were noticeably less sensitive to tag line noise than the cylindrical methods. The APSB and PSB methods had a higher spatial resolution in the radial and circumferential dimensions than the cylindrical methods. The cylindrical methods, however, had a higher spatial resolution in the longitudinal direction. There are two primary reasons for this increased resolution. First, in the cylindrical model, the longitudinal deformation is reconstructed almost entirely from the long-axis tag line points, whereas, in the APSB and PSB models, each deformation component is a function of the tag lines in all three directions. Second, in the cloud experiment data set, there are only six, radially-prescribed long-axis image planes [10], so the long-axis data is relatively sparse (?2000 tag points out of a total of ?10,000 in all three directions). Consequently, a longitudinally-displaced cloud of tag points will have more e?ect on the longitudinal displacement reconstructed using the cylindrical model because there will be a larger per- centage of the total number of points. If, for example, a longitudinal displacement cloud contains 250 points (typical for an 8mm cloud), this would comprise 12.5% of the points used to compute the longitudinal displacement in the cylindrical models but only 2.5% of the points in the APSB and PSB models. As a result, the cylindrical model is slightly more responsive to displacement clouds in the longitudinal direction than the APSB and PSB models. When reconstructing deformation from the pathology database data, the cylin- drical methods had slightly lower tag distance errors than the APSB and PSB methods in most cases. We performed the tests described above on a regularized version of the cylindrical method and the di?erence between the results with and without regularization was negligible. From this we conclude that the improved performance of the APSB relative 69 to the cylindrical models is primarily due to the prolate-spheroidal coordinate system and other di?erences such as removal of a?ne deformation modes before B-spline model fitting, and enforcing spatial smoothness across the apex. Another advantage of the APSB and PSB models over the cylindrical models is that the APSB and PSB models can enforce smoothing across and compute strain at the apex. The cylindrical B-spline models have a singularity at the apex because the radial coordi- nate is zero there. The Cartesian-coordinate models (DMF [17] and TTT) do not have this problem, but they do not enforce deformation smoothness in physiologically-relevant direc- tions. The APSB and PSB methods, however, use equality constraints to stitch together the deformation model at the apex. While the smoothness directions do not correspond to the exact fiber directions in the apical myocardium, they are significantly closer than the smoothness directions in Cartesian-coordinate models. One disadvantage of the APSB and PSB models over cylindrical and Cartesian models is that, in addition to a central axis, the focal point of the prolate spheroidal coordinate system must be computed for each fit. However, our center sensitivity test demonstrated that while the APSB and PSB methods were more sensitive to perturbations in the center than the cylindrical methods, the APSB and PSB sensitivities were still acceptable for circumferential and longitudinal strains. Radial strains were quite sensitive to center perturbations at relatively high noise levels, but radial strains are di?cult to measure with any method from parallel or grid-tagged data [10]. In conclusion, the APSB method can accurately reconstruct deformation and strain in the LV wall from tagged MR images and has several advantages relative to existing tech- niques. Relative to Cartesian methods such as the DMF and TTT, the APSB enforces smoothness in directions that are nominally aligned with the fiber and cross-fiber directions 70 in the myocardium. Relative to cylindrical coordinate techniques, the APSB method re- constructed strain with fewer artifacts, was less sensitive to tag line noise, and had higher radial and circumferential spatial resolution. Relative to prolate spheroidal methods with global basis functions such as TEA, the locally-supported basis functions used by the APSB o?er improved sensitivity to tag line noise while preserving spatial resolution. 71 Chapter 4 Combined Tag Tracking and Strain Reconstruction 4.1 Introduction Three-dimensional tagged cardiac MR image analysis typically consists of two steps. First, tag lines are tracked in each image in the data set. Then, a deformation model is fitted to the tag line positions. Methods that combine tag tracking and motion reconstruction into a single step have been proposed recently. However, as discussed in Section 1.2, there are some limitations in existing techniques. This chapter presents a new combined tag tracking and motion reconstruction method for computing 3-D myocardial strain from tagged cardiac MR images. We call this al- gorithm the Prolate Spheroidal COmbined Tag Tracking and Strain (E) Reconstruction (PS-COTTER) algorithm. In PS-COTTER, a 3-D prolate spheroidal B-spline (PSB) de- formation model presented in Chapter 3 is fitted to detected tag points, raw or filtered image data, thus combining tag tracking and motion reconstruction. The set of candidate tag line points can be identified in the images using peak detection [48, 27], template matching [11], or other techniques. The deformation model ensures that tag line positions are consistent from slice to slice. In contrast to some existing techniques in [14, 25], this algorithm does not rely on user-defined myocardial contours. In contrast to the COTTER algorithm [15], PS-COTTER uses a prolate spheroidal domain. This domain matches the shape of the LV better than the cylindrical domain of COTTER, and with PS-COTTER strain can be computed at the apex. 72 The initial geometry of the deformation model is determined from user-selected end- diastolic circular (for short-axis images) and elliptical (for long-axis images) regions of inter- est (ROIs). The algorithm then optimizes the control points of a prolate spheroidal B-spline model so that tag planes deformed by the model match the tag lines in the images. The myocardium motion field and strain can be reconstructed directly from the results of the combined tag tracking procedure without intermediate steps such as backward and forward displacement field reconstructions. The PS-COTTER algorithm was validated on the imaging studies from both normal human volunteers and patients with myocardial infarction, hypertension, diabetics, or mi- tral regurgitation. The displacement reconstruction accuracy was evaluated by computing the root-mean-square (RMS) tag location di?erence between user-defined tag points and simulated tag points with our algorithm. This chapter is organized as follows. In Section 4.2, the PS-COTTER algorithm is developed based on the PSB model. In Section 4.3, we conduct a parameter sensitivity test for PS-COTTER. Section 4.4 presents the methods used for the algorithm validation. Then, the algorithm is validated on parallel-tagged human MRI data in Section 4.5 and on grid-tagged human MRI studies in Section 4.6. Finally, Section 4.7 contains a discussion of results and conclusion. 4.2 Algorithm Development 4.2.1 Coordinate System and Initialization A prolate spheroidal B-spline deformable model as shown in Figure 4.1(a) is used to describe the deformation of the LV. The initial shape of the model is determined by user- selected ROIs defined in the end-diastolic images. The LV apex and base positions are 73 z x y ? ? foci Undeformed Tag Planes (a) (b) Figure 4.1: (a): Initial shape of the prolate spheroidal B-spline model. (b): Undeformed tag planes defined within the prolate spheroidal B-spline model as uniformly spaced material planes. approximated from the elliptical ROIs for the long-axis images. Simulated tag planes in orthogonal orientations are then defined within the model as uniformly spaced material planes as shown in Figure 4.1(b). The positions of the undeformed planes are specified by the user. The prolate spheroidal (PS) coordinate system is tied to the image prescription. For parallel-tagged MRI, the short-axis and radial image prescriptions were used to construct the PS coordinate system. For grid-tagged MRI, the short-axis, 2-chamber and 4-chamber image prescriptions were used to construct the PS coordinate system. The x and y axes of the prolate spheroidal model are parallel to the row and column directions of short-axis images, and the z?axis is perpendicular to the short-axis image planes. The origin is set to be one-third of the distance from the base to the apex on the common intersection line 74 of the long-axis image planes for the parallel-tagged group, and on the line across the (x,y) mass center of the ROI circles for the grid-tagged group; the focal radius is set to make the radial (?) coordinate of the LV apex equal to 1. 4.2.2 Image Pre-processing PS-COTTER can be applied to raw or filtered image data. In either case, the 3-D PSB deformation model is fit to parallel tag lines. For parallel-tagged MRI, the 3-D PSB model is applied directly to the raw images. However, to e?ectively track tag lines in grid-tagged MR images, filtering is used to eliminate the tag lines in one direction. Gaussian filters as shown in Figure 4.2 are applied to obtain the estimated parallel-tagged images from the grid-tagged ones. We choose Gaussian filters because they have optimal localization properties in both spatial and frequency domain. Figure 4.3 shows an example of the grid-tagged image and the corresponding filtered images with tag lines oriented in a single direction. PS-COTTER can also be applied to tag points detected in a set of images. In this case, the peak detection algorithm [48] is implemented to detect tag points. In [48], a banana- shaped composite Gabor filter is first applied to enhance the tag lines in one direction and to eliminate the tag lines in the other direction. Then tag points are detected by detecting local negatives or ?valleys? in the filtered image. This algorithm has high detection rate. However, it can also introduce false positives. Figure 4.4 shows some tag detection results of this algorithm. Detected points of the tag lines in di?erent directions are displayed in white. 4.2.3 Deformation Model The deformation of a tag plane is described by specifying the spatial position, s,ofa uniformly spaced mesh of material points, p,wherep =[?,?,?] T . As in Section 3.2.2, the 75 Figure 4.2: A grid-tagged MR image and its spectrum (top), and the spectrums of Gaussian filters (bottom) for obtaining the parallel-tagged images. 76 Figure 4.3: Grid-tagged MR images (left) and the corresponding filtered images (middle, right) near end-diastole (top) and near end-systole (bottom). Figure 4.4: Tag points detected by using the peak detection algorithm [48]. 77 spatial and material points are related by s = p + u, (4.1) where u is the 3-D displacement vector at p. We model the 3-D displacements with a 3-D B-spline, which is the tensor product of three 1-D B-splines. Let u =[u ? (l,p,t),u ? (l,p,t),u ? (l,p,t)] T denote the displacement vector, then u ? (l,p,t)= L?1 summationdisplay i=0 P?1 summationdisplay j=0 T?1 summationdisplay k=0 B o i (l)B o j (p)B c k (t)c ? ijk u ? (l,p,t)= L?1 summationdisplay i=0 P?1 summationdisplay j=0 T?1 summationdisplay k=0 B o i (l)B o j (p)B c k (t)c ? ijk (4.2) u ? (l,p,t)= L?1 summationdisplay i=0 P?1 summationdisplay j=0 T?1 summationdisplay k=0 B o i (l)B o j (p)B c k (t)c ? ijk , where B o i (?) is an open B-spline, B c k (?) is a closed B-spline and c ijk is a control point in the dimension denoted by its superscript, and L, P, T are the numbers of control points in radial (?), longitudinal (?), and circumferential (?) directions, respectively. We use open B-splines in the radial and longitudinal directions and closed B-splines in the circumferential direction. In the radial direction, the knots are uniformly spaced between the minimum and maximum radial coordinates (? min and ? max ) of the uniformly spaced mesh of material points. In the longitudinal direction, the knots are uniformly spaced between the most basal points (? min )andthemostapicalpoints(? max )forthe open-apex model or the apex (? max = ?) for the closed-apex model. In the circumferential direction, the knots are uniformly spaced between 0 and 360 degrees. The prolate spheroidal 78 coordinates ?, ? and ? are related to the B-spline parametric coordinates l, p and t by l = ?? ? min ? max ? ? min ? (L? O +1) p = ? ? ? min ? max ? ? min ? (P ? O +1) t = ? 2? ? T, where O is the spline order. The PSB model can be either closed-apex or open-apex. For closed apex, we penalize the di?erence between the circumferential control points at the apex, and the di?erence between the first and that between the second order derivatives of each pair of antiparallel longitudinal displacement components at the apex. This is to ensure a smooth deformation across each point where the circumferential closed splines shrink to one point and the lon- gitudinal open splines meet. Details of the apex penalty function are given in Section 4.2.4 (Equation (4.9)). For open apex, no penalty is imposed at the apex. 4.2.4 Model Fitting To track the tag features in the image, we define a set of simulated tag line points on the intersections of the tag planes and image planes, as shown in Figure 4.5. Given an image plane and a deformed tag plane, the locations of the simulated deformed tag line points are computed by linearly interpolating the spatial positions s similar to the technique in [14, 15]. The control points in Equation (4.2) are optimized such that the simulated tag lines move towards the real tag lines in the tagged MR images. The control points of the first time frames are initialized to zeros. The initial control points of a later time frame are based on those of its neighboring previous time frames as in [15]. The optimal control 79 Tag Planes Image Plane Left Ventricle Simulated Tag Points Figure 4.5: Simulated tag points on the intersection lines of the imaging planes and tag planes. points are computed by minimizing a weighted sum of the external and internal energy for the open-apex PS-COTTER: [C ? ? ,C ? ? ,C ? ? ]= arg min E external + ?E internal . C ? ? ,C ? ? ,C ? ? (4.3) For the closed-apex PS-COTTER, an additional apical penalty is included in the energy equation: [C ? ? ,C ? ? ,C ? ? ]= arg min E external + ?E internal + E apex . C ? ? ,C ? ? ,C ? ? (4.4) The MATLAB command lsqnonlin was used for the optimization. The external energy E external measures the degree of match between simulated de- formed tag line points and the real tag lines in the images. 80 For PS-COTTER on detected tag points, the closeness of a deformed simulated tag point to a real tag point is computed by its distance to the nearest candidate tag line point. The external energy E external is defined as the sum of the squared distance over all simulated tag points ? t k (i,j): E external = summationdisplay i summationdisplay j summationdisplay k dist 2 ( ? t k (i,j)), (4.5) where i is the index of an image plane, j is the index of a tag plane, and k is the index of a simulated tag point on the tag line where image plane i and tag plane j intersect. For PS-COTTER on raw or filtered image data, the formula used in [26, 15] is applied to compute the external energy. For a deformed simulated tag point ? t k (i,j), its closeness to a real tag point is computed by the squared di?erence between an ideal tag profile centered at ? t k (i,j) and the local image data [26], which is L( ? t k (i,j)) = [ ? tag ( ? t k (i,j))? ? w] T [ ? tag ( ? t k (i,j))? ? w] (4.6) where ? w is a 1 ? (2N s + 1) neighborhood of image pixels along a line perpendicular to the tag plane and centered at the pixel nearest to ? t k (i,j). N s is the next integer greater than the FWHM (Full Width at Half Maximum) of the tag line in pixels. ? tag is a 1?(2N s +1) ideal local tag profile centered at ? t k (i,j). The ideal tag profile is defined as the element-wise multiplication of the inverted Gaussian tag profile [11] centered at ? t k (i,j) and an estimate of the untagged magnitude image data in the same neighborhood as ? w. Then the external energy is defined as E external = summationdisplay i summationdisplay j summationdisplay k L( ? t k (i,j)), (4.7) where the notations of i, j and k are the same as in (4.5). 81 The internal energy E internal measures the degree of bending of the 3-D displacement field. It smoothes the reconstructed displacement field. We define E internal as E internal = summationdisplay i summationdisplay j summationdisplay k parenleftBig triangleinv 2 u ? (? ijk ) parenrightBig 2 + parenleftBig triangleinv 2 u ? (? ijk ) parenrightBig 2 + parenleftBig triangleinv 2 u ? (? ijk ) parenrightBig 2 , (4.8) where ? ijk is a knot located in the B-spline model except the apex, u ? ,u ? ,u ? denote the displacement field defined in (4.2), and triangleinv 2 denotes the Laplacian operator in prolate spheroidal coordinates. Because Laplacian operator has a singularity at ? = ?,wedonot impose the regularization at the apex. The apical energy E apex measures the consistency of the apical displacement field. E apex is computed only at the apex, and is given by E apex = L?1 summationdisplay i=0 T?1 summationdisplay k=1 bracketleftbigg parenleftBig c ? iA0 ? c ? iAk parenrightBig 2 + parenleftBig c ? iA0 ? c ? iAk parenrightBig 2 + parenleftBig c ? iA0 ? c ? iAk parenrightBig 2 bracketrightbigg + L?1 summationdisplay i=0 T 2 ?1 summationdisplay k=0 ? ? ?u ? ?? vextendsingle vextendsingle vextendsingle vextendsingle ? iAk ? ?u ? ?? vextendsingle vextendsingle vextendsingle vextendsingle ? iA(k+ T 2 ) ? ? 2 + L?1 summationdisplay i=0 T 2 ?1 summationdisplay k=0 ? ? ? ? 2 u ? ?? 2 vextendsingle vextendsingle vextendsingle vextendsingle vextendsingle ? iAk ? ? 2 u ? ?? 2 vextendsingle vextendsingle vextendsingle vextendsingle vextendsingle ? iA(k+ T 2 ) ? ? ? 2 , (4.9) where c iAk is an apical control point, ? iAk is a B-spline knot located at the apex and the number of circumferential control points T is an even integer. The first term is the summation of the squared di?erences between the circumferential control points. The second term is the summation of the squared di?erences between the first order derivatives of each pair of opposing longitudinal displacement components. And the third term is the summation of the squared di?erences between the second order derivatives of each pair of opposing longitudinal displacement components. The apical penalty helps make a consistent deformation field across the apex. 82 Lagrangian strain can be computed in Cartesian coordinates directly from the recon- structed displacement field as in Section 3.3.3. For consistency with previously published results [10, 18, 15, 17, 26, 8, 16, 20], the Lagrangian strain tensor is rotated to the prolate spheroidal coordinate system used in these papers. This results in three normal strains (E rr : radial thickening, E cc : circumferential shortening and E ll : longitudinal shorting) and three shear strains (E rc , E rl and E cl ). 4.3 Parameter Sensitivity For parallel-tagged MR data, the parameters of the PS-COTTER algorithm were op- timized on a set of 10 human imaging studies, composed of six normal volunteers and four patients with myocardial infarction. Another 8 human imaging studies were used for evalua- tion. Each study had 6 short-axis image planes with no interslice gap, 6 radially prescribed long-axis image planes with 30 o angular spacing. In each study, the slice thickness was 10mm, the tag spacing was 7mm, the tag line FWHM was about 1 pixel, and the temporal resolution was 32.5msec. For grid-tagged MR data, the parameters of PS-COTTER were optimized on a set of 12 human imaging studies, composed of 3 normal human and 3 groups of patients, each group consisting of 3 patients of hypertension, diabetics or mitral regurgitation. Another 41 human imaging studies, 12 from each of the normal, hypertensive and diabetic groups, and 5 from the mitral regurgitation group, were used for evaluation. Each imaging study had around 10 short-axis image planes with no interslice gap, one 2-chamber slice and one 4-chamber slice. Some studies had an additional long-axis slice. In each study, the slice thickness was around 6.36mm, the tag spacing was around 7 pixels with pixel spacing 1.5625mm?1.5625mm, and the temporal resolution was around 56msec. 83 The user-defined parameters of PS-COTTER are the numbers of control points in the three dimensions (radial, longitudinal, and circumferential), and the regularization weight ?. Quadratic B-splines were used in the reconstructions for good local control and smoothness of the reconstructed displacement field. First, we determined the regularization weight ? in (4.3) and (4.4). We empirically set the numbers of control points for the PSB model and compute the RMS tag distance error over a range of ? values over the development studies. The ? values were chosen on a logarithmic scale from a number greater than or equal to 1 down to 10 ?7 . The di?erence between the e?ect of a regularization weight less than 10 ?7 and that of 10 ?7 was not expected to be significant. The RMS tag distance error was defined as RMS tag error = radicalBigg summationtext (tag user ? tag sim ) 2 N , (4.10) where N is the total number of tag points, tag user are the tag points tracked with a user- supervised process, and tag sim are the simulated tag points with the PS-COTTER algo- rithm. Tables 4.1 and 4.7 show the RMS tag distance errors in pixels of the quadratic B-spline reconstructions with 3 radial ? 4 longitudinal ? 6 circumferential control points and various ? over the 10 parallel-tagged development studies (Table 4.1) and over the 12 grid-tagged development studies (Table 4.7). For each reconstruction, the ? value that resulted in the minimum RMS tag tracking error was selected for the remaining experiments. Then, we studied the e?ect of the numbers of B-spline control points on the tag tracking accuracy at end-systole over the development studies for each of closed-apex/open-apex PS-COTTER reconstructions on detected tag points, raw or filtered image data, as shown in Tables 4.2, 4.3, 4.4 and 4.5 for the parallel-tagged image data and in Tables 4.8, 4.9, 4.10 and 4.11 for the grid-tagged image data. For quadratic B-spline reconstructions, the 84 number of control points in each dimension should be at least 3. For the closed-apex B-spline model, however, due to the apical penalty, the number of control points in the circumferential direction should be an even number and thus it is at least 4. The set of control point numbers that resulted in minimum RMS tag distance error was used in the following validations. Table 4.6 shows the optimal parameters and the corresponding RMS tag errors for each of the quadratic PS-COTTER reconstructions on the 10 parallel-tagged human development imaging studies. Table 4.12 shows the optimal parameters and the corresponding RMS tag errors over the 12 grid-tagged development studies. Table 4.1: RMS tag distance error in pixels of the quadratic B-spline reconstructions with 3radial? 4 longitudinal ? 6 circumferential control points and various ? over the 10 parallel-tagged human development imaging studies. On raw image data On detected tag points ? Closed-apex Open-apex Closed-apex Open-apex 10 ?7 0.6667 0.65 0.84 0.78 10 ?6 0.6663 0.68 0.85 0.80 10 ?5 0.6657 0.68 0.79 0.82 10 ?4 0.67 0.72 0.76 0.85 10 ?3 0.6664 0.73 0.86 1.30 10 ?2 0.68 0.82 0.87 1.59 10 ?1 0.73 1.24 1.36 1.74 1 0.80 1.51 1.69 2.13 10 1.22 1.82 2.11 2.25 10 2 1.57 2.20 2.44 2.32 The PS-COTTER algorithm was implemented in MATLAB. It took about 10 minutes to process a single time frame of a human MR imaging study on a 3.0GHz PC. 4.4 Validation Methods PS-COTTER was validated on a set of 8 parallel-tagged human imaging studies and another set of 41 grid-tagged human imaging studies that were not used in optimizing 85 Table 4.2: RMS tag distance error in pixels of the quadratic closed-apex PS-COTTER reconstruction on raw image data with various numbers of control points and regularization weight ? =10 ?5 over the 10 parallel-tagged human development imaging studies. L PT 345 4 0.88 0.80 0.80 360.72 0.68 0.71 8 0.68 0.69 0.69 4 0.78 0.72 0.75 460.67 0.656 0.70 8 0.68 0.68 0.662 4 0.75 0.79 0.78 560.664 0.71 0.69 8 0.69 0.74 0.71 Table 4.3: RMS tag distance error in pixels of the quadratic open-apex PS-COTTER re- construction on raw image data with various numbers of control points and regularization weight ? =10 ?7 over the 10 parallel-tagged human development imaging studies. L PT 34 3 0.94 0.94 4 0.81 0.82 350.71 0.73 6 0.73 0.69 7 0.71 0.69 8 0.74 0.69 3 0.94 0.93 4 0.79 0.80 450.68 0.68 6 0.65 0.68 7 0.66 0.69 8 0.66 0.66 3 1.03 0.97 4 0.81 0.89 550.84 0.77 6 0.75 0.76 7 0.76 0.77 8 0.77 0.79 86 Table 4.4: RMS tag distance error in pixels of the quadratic closed-apex PS-COTTER reconstruction on detetced tag points with various numbers of control points and regular- ization weight ? =10 ?4 over the 10 parallel-tagged human development imaging studies. L PT 3456 4 0.82 0.79 0.80 0.84 360.86 0.76 0.79 0.89 8 0.91 0.78 0.86 0.93 4 0.91 0.86 0.79 0.78 460.76 0.75 0.72 0.78 8 0.85 0.79 0.83 0.75 4 0.91 0.88 0.82 0.75 560.88 0.80 0.80 0.76 8 0.88 0.82 0.77 0.76 Table 4.5: RMS tag distance error in pixels of the quadratic open-apex PS-COTTER recon- struction on detetced tag points with various numbers of control points and regularization weight ? =10 ?7 over the 10 parallel-tagged human development imaging studies. L PT 345 3 0.88 0.84 0.84 4 0.83 0.81 0.80 350.75 0.74 0.76 6 0.78 0.80 0.76 7 0.77 0.80 0.81 8 0.83 0.82 0.78 3 0.90 0.88 0.84 4 0.79 0.83 0.77 450.77 0.74 0.74 6 0.78 0.69 0.73 7 0.81 0.76 0.73 8 0.81 0.77 0.78 3 0.99 0.94 0.93 4 0.98 0.87 0.86 550.87 0.87 0.90 6 0.97 0.94 0.84 7 1.03 0.93 0.86 8 1.02 0.95 0.91 87 Table 4.6: Optimal parameters and the corresponding RMS tag distance errors in pixels of the quadratic PS-COTTER reconstructions over the 10 parallel-tagged human development imaging studies. On raw image data On detected tag points Closed-apex Open-apex Closed-apex Open-apex L 4354 P 4444 T 6666 ? 10 ?5 10 ?7 10 ?4 10 ?7 RMS tag di?erence 0.66 0.65 0.72 0.69 Table 4.7: RMS tag distance error in pixels of the quadratic B-spline reconstructions with 3radial? 4 longitudinal ? 6 circumferential control points and various ? over the 12 grid-tagged human development imaging studies. On filtered image data On detected tag points ? Closed-apex Open-apex Closed-apex Open-apex 10 ?7 1.069 1.09 1.11 1.090 10 ?6 1.067 1.07 1.13 1.16 10 ?5 1.11 1.12 1.08 1.087 10 ?4 1.12 1.15 1.04 1.17 10 ?3 1.15 1.14 1.08 1.16 10 ?2 1.28 1.19 1.16 1.15 10 ?1 1.32 1.24 1.38 1.35 1 1.45 1.34 1.58 1.67 Table 4.8: RMS tag distance error in pixels of the quadratic closed-apex PS-COTTER reconstruction on filtered image data with various numbers of control points and regular- ization weight ? =10 ?6 over the 12 grid-tagged human development imaging studies. L PT 345 4 1.09 1.08 1.06 361.07 1.04 1.13 8 1.11 1.12 1.10 4 1.08 1.05 1.05 461.07 1.08 1.10 8 1.10 1.14 1.12 88 Table 4.9: RMS tag distance error in pixels of the quadratic open-apex PS-COTTER recon- struction on filtered image data with various numbers of control points and regularization weight ? =10 ?6 over the 12 grid-tagged human development imaging studies. L PT 345 3 1.15 1.10 1.09 4 1.11 1.07 1.10 351.14 1.08 1.06 6 1.10 1.06 1.06 7 1.16 1.09 1.10 8 1.22 1.09 1.18 3 1.13 1.11 1.09 4 1.09 1.05 1.09 451.06 1.04 1.03 6 1.07 1.11 1.10 7 1.09 1.10 1.14 8 1.10 1.11 1.15 3 1.07 1.03 1.05 4 1.04 1.00 1.01 551.02 0.98 1.01 6 1.04 1.02 1.04 7 1.10 1.07 1.08 8 1.12 1.10 1.12 3 1.10 1.06 1.06 4 1.03 1.01 1.07 651.04 1.05 1.05 6 1.06 1.07 1.06 7 1.05 1.08 1.11 8 1.15 1.11 1.20 89 Table 4.10: RMS tag distance error in pixels of the quadratic closed-apex PS-COTTER reconstruction on detected tag points with various numbers of control points and regular- ization weight ? =10 ?4 over the 12 grid-tagged human development imaging studies. L PT 345 4 0.89 0.87 0.94 360.96 0.91 0.90 8 1.00 0.89 0.98 4 0.95 0.95 0.92 461.04 0.91 0.93 8 1.21 1.01 0.93 Table 4.11: RMS tag distance error in pixels of the quadratic open-apex PS-COTTER reconstruction on detetced tag points with various numbers of control points and regular- ization weight ? =10 ?5 over the 12 grid-tagged human development imaging studies. L PT 34 56 3 1.01 1.00 0.99 0.99 4 0.95 0.92 0.901 0.91 350.93 0.904 0.93 0.97 6 0.99 1.15 0.99 1.03 7 1.03 1.03 1.04 1.01 8 1.12 1.02 1.05 1.21 3 1.03 0.98 0.94 0.94 4 0.93 1.01 0.92 0.93 450.99 0.97 1.02 1.05 6 1.09 1.08 1.21 1.18 7 1.15 1.15 1.18 1.26 8 1.24 1.25 1.12 1.22 90 Table 4.12: Optimal parameters and the corresponding RMS tag distance errors in pixels of the quadratic PS-COTTER reconstructions over the 12 grid-tagged human development imaging studies. On filtered image data On detected tag points Closed-apex Open-apex Closed-apex Open-apex L 4445 P 3533 T 6544 ? 10 ?6 10 ?6 10 ?4 10 ?5 RMS tag di?erence 1.04 0.98 0.87 0.90 parameters. The studies are identified in the appendix in Section 4.8. These imaging studies, composed of normal volunteers and patients, were used to evaluate the tag tracking accuracy and compare strains reconstructed with PS-COTTER to those reconstructed with user-defined analysis. 4.4.1 Analysis of Tag Tracking Errors The tag tracking error was analyzed by computing the RMS tag position di?erence between the simulated tags generated with PS-COTTER and the tag lines tracked with user-supervised processing [12]. For the parallel-tagged group, the tag tracking error was analyzed for each time frame in the sequence for comparing the performance of PS-COTTER and COTTER [15]. In many clinical studies, only the end-systolic strain in the LV wall is of interest, so only the end-systolic tag tracking error needs to be analyzed. Thus, for the grid-tagged group, only the end-systolic tag errors were studied. The RMS tag position di?erence was computed with Equation (4.10). 91 4.4.2 3-D Strain Reconstruction The myocardium strain was computed directly from the reconstructed displacement field of the combined tag tracking procedure as in Section 3.3.3 without intermediate steps, such as backward and forward deformation model fittings. For consistency with previously published results [10, 18, 15, 17, 26, 8, 16, 20], the Lagrangian strain tensor was rotated to the prolate spheroidal coordinate system used in these papers. 4.4.3 Statistical Methods The end-systolic strains measured by PS-COTTER were compared with those com- puted with other methods for each of the human evaluation imaging studies by correlation analysis and Bland-Altman plots. For the parallel-tagged MRI data, we compared the end-systolic average midwall strains measured by PS-COTTER to those computed by user- supervised methods TEA [16] and DMF [17]. For the grid-tagged MRI data, we compared the end-systolic average midwall strains computed with PS-COTTER and the closed-apex APSB model presented in Chapter 3. 4.5 Validation Results on parallel-tagged MRI data 4.5.1 Analysis of Tag Tracking Errors The final tag tracking results are shown on short-axis and long-axis images from a parallel-tagged human imaging study in Figure 4.6. A mid-ventricular short-axis slice with a tag plane angle of 0 o and 90 o , and a long-axis slice are displayed. Only the images within ROI are shown. The simulated tag lines are consistent with the dark stripes in the images. 92 Table 4.13: RMS tag distance errors in pixels of the quadratic PS-COTTER reconstructions over the 8 parallel-tagged human evaluation imaging studies. On raw image data On detected tag points Closed-apex Open-apex Closed-apex Open-apex RMS tag di?erence 0.55 0.64 0.65 0.69 The end-systolic RMS tag position di?erences over the eight evaluation studies are shown in Table 4.13. The closed-apex PS-COTTER on raw image data resulted in the minimum tag distance error over the eight human evaluation studies, which was 0.55 pixels. The RMS di?erence between tags simulated by closed-apex PS-COTTER on raw image data and tags tracked with a user-supervised procedure over the eight human evaluation imaging studies is plotted vs. time after detection of the electrocardiograph QRS complex in Figure 4.7. The RMS tag di?erence increases with time because the algorithm uses the tags of previous time frame as the initial condition of the current time frame, thus, there is an accumulation of error. However, all the RMS tag errors are less than 0.8 pixels. 4.5.2 3-D Strain Reconstruction Three-dimensional end-systolic strains of a normal volunteer and a patient with a myocardial infarction reconstructed with PS-COTTER are shown in Figure 4.8. The normal volunteer exhibits fairly uniform radial thickening (positive strain) and circumferential and longitudinal shortening (negative strains). In contrast, the patient shows negative radial thickening or radial thinning in the anterior wall, which was infarcted, as well as positive circumferential and longitudinal shortening. 93 Figure 4.6: Simulated tag lines of PS-COTTER for a mid-ventricular short-axis slice with a tag plane angle of 0 o (left) and 90 o (middle), and one long-axis slice (right). The top, middle and bottom rows correspond to early, mid- and late systole. 94 0 33 65 98 130 163 195 228 260 293 325 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Time from QRS in msec RMS difference in puxels Figure 4.7: RMS tag di?erence in pixels between the simulated tags generated with closed- apex PS-COTTER on raw image data and a set of tag lines tracked with user-supervised processing for the 8 parallel-tagged human evaluation imaging studies. 95 Figure 4.8: Maps of 3-D radial thickening (E rr ), circumferential shortening (E cc ), and longitudinal shortening (E ll ) reconstructed with closed-apex PS-COTTER on raw image data for a normal human (left) and a patient with an anterior myocardial infarction (right). The red balls denote the septum and the green ball denotes the apex. Each image is color-coded, with yellow indicating thickening and blue indicating shortening. 96 4.5.3 Statistics Table 4.14 shows the correlation coe?cients between end-systolic radial thickening E rr , circumferential shortening E cc , longitudinal shortening E ll and the minimum principle strain E min computed with closed-apex PS-COTTER on raw image data and the user- supervised reconstruction methods TEA and DMF in the midwall. A P-value of 0.05 or less is considered statistically significant. The average midwall strain data of the eight human evaluation imaging studies show that PS-COTTER showed good correlation with TEA and DMF in E cc , E ll and E min . There was no statistically significant correlation between the E rr strains. E rr correlated weaker because the parallel tag pattern more sparsely samples the radial motion. Comparisons between the three sets of strains are also shown with scatter plots and Bland-Altman plots in Figures 4.9 and 4.10. Most of the points are between the mean ? 2SD lines, indicating good agreement between PS-COTTER and TEA/DMF strains. 97 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.1 0.2 0.3 0.4 TEA PS?COTTER Err ( ?=0.65) 0.2 0.25 0.3 0.35 0.4 ?0.2 ?0.1 0 0.1 0.2 (PS?COTTER+TEA)/2 PS?COTTER?TEA Err ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 ?0.3 ?0.25 ?0.2 ?0.15 ?0.1 TEA PS?COTTER Ecc ( ?=0.99) ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 ?0.01 ?0.005 0 0.005 0.01 (PS?COTTER+TEA)/2 PS?COTTER?TEA Ecc ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.2 ?0.15 ?0.1 ?0.05 0 TEA PS?COTTER Ell ( ?=0.82) ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.05 0 0.05 0.1 0.15 (PS?COTTER+TEA)/2 PS?COTTER?TEA Ell ?0.25 ?0.2 ?0.15 ?0.25 ?0.2 ?0.15 TEA PS?COTTER Emin ( ?=0.98) ?0.25 ?0.2 ?0.15 ?0.02 ?0.01 0 0.01 (PS?COTTER+TEA)/2 PS?COTTER?TEA Emin Figure 4.9: Correlation(left) and Bland-Altman (right) plots for comparing the average midwall strains computed using closed-apex PS-COTTER on raw image data and those computed using TEA over the 8 parallel-tagged human evaluation studies. In the correlation plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence line represent the mean ? 2SD of the di?erence. 98 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.1 0.2 0.3 0.4 DMF PS?COTTER Err ( ?=0.55) 0.2 0.25 0.3 0.35 0.4 ?0.2 0 0.2 0.4 0.6 (PS?COTTER+DMF)/2 PS?COTTER?DMF Err ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 ?0.3 ?0.25 ?0.2 ?0.15 ?0.1 DMF PS?COTTER Ecc ( ?=0.99) ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 ?15 ?10 ?5 0 5 x 10 ?3 (PS?COTTER+DMF)/2 PS?COTTER?DMF Ecc ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.2 ?0.15 ?0.1 ?0.05 0 DMF PS?COTTER Ell ( ?=0.92) ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.05 0 0.05 (PS?COTTER+DMF)/2 PS?COTTER?DMF Ell ?0.25 ?0.2 ?0.15 ?0.25 ?0.2 ?0.15 DMF PS?COTTER Emin ( ?=0.96) ?0.25 ?0.2 ?0.15 ?0.04 ?0.02 0 0.02 (PS?COTTER+DMF)/2 PS?COTTER?DMF Emin Figure 4.10: Correlation(left) and Bland-Altman (right) plots for comparing the average midwall strains computed using closed-apex PS-COTTER on raw image data and those computed using DMF over the 8 parallel-tagged human evaluation studies. In the correla- tion plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence linerepresentthemean? 2SD of the di?erence. 99 Table 4.14: Correlation Coe?cients (?), Di?erence (Mean ? SD) between the average midwall strains over the 8 parallel-tagged human evaluation imaging studies computed using closed-apex PS-COTTER on raw image data and those computed with user-supervised analysis (TEA and DMF). TEA ?P Mean ? SD Mid-wall Err 0.6460 0.0835 ?0.0231 ? 0.0568 Ecc 0.9935 < 0.0001 ?0.0000 ? 0.0048 Ell 0.8244 0.0118 0.0071 ? 0.0258 Emin 0.9816 < 0.0001 ?0.0068 ? 0.0061 DMF ?P Mean ? SD Mid-wall Err 0.5496 0.1582 0.0471 ? 0.0778 Ecc 0.9927 < 0.0001 ?0.0046 ? 0.0047 Ell 0.9158 0.0014 0.0028 ? 0.0192 Emin 0.9594 0.0002 ?0.0051 ? 0.0090 4.6 Validation Results on Grid-tagged Human Imaging Studies 4.6.1 Analysis of Tag Tracking Errors The final tracking result for a grid-tagged human imaging study is shown in Figure 4.11. A mid-ventricular short-axis slice, a 2-chamber and a 4-chamber slice are displayed. Only the images within ROI are shown. For clarity, tag lines oriented in a single direction are displayed in each image. The simulated tag lines are consistent with the true tag lines in the images. The end-systolic RMS di?erences between tags simulated with PS-COTTER and tags tracked with a user-supervised procedure over the 41 grid-tagged human evaluation imaging studies are shown in Table 4.15. The closed-apex PS-COTTER on detected tag points had the minimum tag distance error, which was 0.71 pixels. 100 Figure 4.11: Simulated tag lines of PS-COTTER for a mid-ventricular short-axis slice (left), a 2-chamber slice (middle), and a 4-chamber slice (right). The top two and bottom two rows correspond to early and late systole. 101 Table 4.15: RMS tag distance errors in pixels of the quadratic PS-COTTER reconstructions with the optimal parameters over the 41 grid-tagged human evaluation imaging studies. On filtered image data On detected tag points Closed-apex Open-apex Closed-apex Open-apex RMS tag di?erence 0.77 0.73 0.71 0.77 4.6.2 3-D Strain Reconstruction Three-dimensional end-systolic strains of a normal volunteer and three patients recon- structed with PS-COTTER are shown in Figure 4.8. One patient had hypertension and was resistant to medical therapy. One patient had diabetics and a myocardial infarction. One patient had mitral regurgitation. The normal volunteer exhibited fairly uniform radial thickening and circumferential and longitudinal shortening. In contrast, the hypertensive patient showed reduced circumferential shortening and longitudinal shortening and the di- abetic patient showed reduced strains in all of E rr , E cc and E ll . The patient with mitral regurgitation showed slightly reduced circumferential shortening. 4.6.3 Statistics Table 4.16 shows the correlation coe?cients between end-systolic E rr , E cc , E ll and E min strains computed with closed-apex PS-COTTER on detected tag points and the strains computed with expert-edited tag points and the closed-apex APSB model in the midwall. The average midwall strain data from the 41 human evaluation imaging studies show that PS-COTTER showed good correlation with APSB in E cc and E min . E ll showed reasonable correlation, but not as strong as E cc and E min .TheE rr correlation was not statistically significant. The basal and mid-ventricular strains correlated stronger than the apex strains. 102 Normal Resistant hypertensive Diabetics with mitral regurgitation Mitral regurgitation Figure 4.12: Maps of 3-D radial thickening (E rr , left), circumferential shortening (E cc ,mid- dle), and longitudinal shortening (E ll , right) reconstructed with closed-apex PS-COTTER on detected tag points from grid-tagged MR data for a normal human, a hypertensive pa- tient, a diabetic patient, and a patient with mitral regurgitation. The red ball denotes the anterior intersection, the green ball denotes mid-wall and the blue ball denotes the in- ferior intersection. Each image is color-coded, with yellow indicating thickening and blue indicating shortening. 103 Table 4.16: Correlation Coe?cients (?), Di?erence (Mean ? SD) between the mid-wall average strains for the 41 grid-tagged human evaluation imaging studies computed using closed-apex PS-COTTER on detected tag points and those computed with the APSB model. ?PMean ? SD Base Err 0.0017 0.9916 0.0201 ? 0.1040 Ecc 0.8159 < 0.0001 ?0.0201 ? 0.0227 Ell 0.5772 0.0001 ?0.0203 ? 0.0400 Emin 0.8252 < 0.0001 ?0.0273 ? 0.0180 Mid-Ventricle Err 0.3125 0.0466 ?0.0578 ? 0.0896 Ecc 0.8326 < 0.0001 ?0.0265 ? 0.0198 Ell 0.7796 < 0.0001 ?0.0125 ? 0.0240 Emin 0.8211 < 0.0001 ?0.0217 ? 0.0212 Apex Err -0.0366 0.8248 ?0.0075 ? 0.1899 Ecc 0.4835 0.0018 ?0.0237 ? 0.0407 Ell 0.5016 0.0011 0.0058 ? 0.0525 Emin 0.6794 < 0.0001 ?0.0588 ? 0.0345 Comparisons between the three sets of strains are also shown with scatter plots and Bland-Altman plots in Figures 4.13, 4.14 and 4.15. Most of the points are between the mean ? 2SD lines, indicating good agreement between PS-COTTER and APSB. 104 ?0.1 ?0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 ?0.5 0 0.5 APSB PS?COTTER Err ( ?=0.00) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 ?0.2 0 0.2 0.4 0.6 (PS?COTTER+APSB)/2 PS?COTTER?APSB Err ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.4 ?0.3 ?0.2 ?0.1 0 APSB PS?COTTER Ecc ( ?=0.82) ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 ?0.1 ?0.05 0 0.05 0.1 (PS?COTTER+APSB)/2 PS?COTTER?APSB Ecc ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.4 ?0.3 ?0.2 ?0.1 0 APSB PS?COTTER Ell ( ?=0.58) ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.2 ?0.1 0 0.1 (PS?COTTER+APSB)/2 PS?COTTER?APSB Ell ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 ?0.3 ?0.25 ?0.2 ?0.15 ?0.1 APSB PS?COTTER Emin ( ?=0.83) ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 ?0.1 ?0.05 0 0.05 0.1 (PS?COTTER+APSB)/2 PS?COTTER?APSB Emin Figure 4.13: Correlation(left) and Bland-Altman (right) plots over the 41 grid-tagged human evaluation imaging studies for comparing average basal midwall strains computed using closed-apex PS-COTTER on detected tag points and those computed using the APSB algorithm. In the correlation plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence line represent mean ? 2SD of the di?erence. 105 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 0.1 0.2 0.3 0.4 APSB PS?COTTER Err ( ?=0.31) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 ?0.4 ?0.2 0 0.2 0.4 (PS?COTTER+APSB)/2 PS?COTTER?APSB Err ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.4 ?0.3 ?0.2 ?0.1 0 APSB PS?COTTER Ecc ( ?=0.83) ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 ?0.1 ?0.05 0 0.05 0.1 (PS?COTTER+APSB)/2 PS?COTTER?APSB Ecc ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.2 ?0.15 ?0.1 ?0.05 0 APSB PS?COTTER Ell ( ?=0.78) ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.1 ?0.05 0 0.05 0.1 (PS?COTTER+APSB)/2 PS?COTTER?APSB Ell ?0.4 ?0.35 ?0.3 ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.3 ?0.25 ?0.2 ?0.15 ?0.1 APSB PS?COTTER Emin ( ?=0.82) ?0.3 ?0.25 ?0.2 ?0.15 ?0.1 ?0.1 ?0.05 0 0.05 0.1 (PS?COTTER+APSB)/2 PS?COTTER?APSB Emin Figure 4.14: Correlation(left) and Bland-Altman (right) plots over the 41 grid-tagged human evaluation imaging studies for comparing average mid-ventricle midwall strains computed using closed-apex PS-COTTER on detected tag points and those computed using the APSB algorithm. In the correlation plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence line represent mean ? 2SD of the di?erence. 106 ?0.2 ?0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 ?1 ?0.5 0 0.5 1 APSB PS?COTTER Err ( ?=?0.04) ?0.2 ?0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 ?0.5 0 0.5 (PS?COTTER+APSB)/2 PS?COTTER?APSB Err ?0.4 ?0.35 ?0.3 ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.4 ?0.3 ?0.2 ?0.1 0 APSB PS?COTTER Ecc ( ?=0.48) ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 ?0.2 ?0.1 0 0.1 0.2 (PS?COTTER+APSB)/2 PS?COTTER?APSB Ecc ?0.2 ?0.15 ?0.1 ?0.05 0 0.05 0.1 0.15 0.2 ?0.5 0 0.5 APSB PS?COTTER Ell ( ?=0.50) ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.2 ?0.1 0 0.1 0.2 (PS?COTTER+APSB)/2 PS?COTTER?APSB Ell ?0.4 ?0.35 ?0.3 ?0.25 ?0.2 ?0.15 ?0.1 ?0.05 0 ?0.4 ?0.3 ?0.2 ?0.1 0 APSB PS?COTTER Emin ( ?=0.68) ?0.35 ?0.3 ?0.25 ?0.2 ?0.15 ?0.15 ?0.1 ?0.05 0 0.05 (PS?COTTER+APSB)/2 PS?COTTER?APSB Emin Figure 4.15: Correlation(left) and Bland-Altman (right) plots over the 41 grid-tagged human evaluation imaging studies for comparing average apical midwall strains computed using closed-apex PS-COTTER on detected tag points and those computed using the APSB algorithm. In the correlation plots, the dashed line represents perfect correlation. In the Bland-Altman plots, the dashed line represents the mean di?erence. The lines above and below the mean di?erence line represent mean ? 2SD of the di?erence. 107 4.7 Discussion This chapter presented the PS-COTTER algorithm for combined tag tracking and motion reconstruction without user-defined myocardial contours. The algorithm fits the 3-D PSB deformation model as in Chapter 3 directly to the detected candidate tag line points, raw or filtered image data to track tags. PS-COTTER optimizes the control points of the PSB model so that tag planes deformed by the model match the real tag lines. The myocardium motion field and strain are computed from the results of tag tracking procedure without the need of intermediate steps. The PS-COTTER algorithm is implemented in a prolate spheroidal domain enclosed by the ROIs, which are user-selected end-diastolic circular and elliptical regions. This prolate spheroidal domain matches the actual shape of the LV more closely than the cylindrical domain in [15]. At end-systole, the tag lines in the myocardium have significant deformation, while the tags around the LV are relatively static. Those static tags can prevent the PSB model from deforming. By using the prolate spheroidal domain, the e?ect of those static tags on motion reconstruction is reduced. Another advantage of the PS-COTTER relative to COTTER [15] is that strain can be computed anywhere in the myocardium. In COTTER, the radial coordinate is zero at the apex and strain cannot be computed there. However, in the closed-apex PS-COTTER, strain at the apex can be computed as described in Chapter 3. The same parallel-tagged MR data set as used in COTTER [15] was employed in the experiments for PS-COTTER. Compared with COTTER (end-systolic RMS tag distance error of 0.66 pixels), open-apex PS-COTTER on raw image data (0.64 pixels) and closed- apex PS-COTTER on detected points (0.65 pixels) had marginal increase in tag tracking accuracy. Closed-apex PS-COTTER on raw image data had a smaller end-systolic RMS 108 tag distance error (0.55 pixels) on the 8 human evaluation imaging studies, resulting in a 17% increase in tag tracking accuracy. On detected tag points, the closed-apex PS-COTTER outperformed the open-apex PS-COTTER in tag tracking accuracy. The closed-apex PS-COTTER uses apical penalty to stitch together the deformation model at the apex. At end-systole, the apex twists more significantly and experiences larger displacement than the other part of the left ventricle. The apical penalty helps make the reconstructed displacement field consistent across the apex, and thus the tag tracking was more accurate. On raw or filtered image data, PS-COTTER had smaller RMS tag distance errors for the parallel-tagged MR images than for the grid-tagged MR data. There are several reasons. First, the tag lines in the parallel-tagged images were much thinner than those in the grid-tagged images. The wider the tag lines, the more significant the di?erence between the automatic tracked results and those with a user-supervised process. Second, di?erent temporal resolutions were used to acquire the two groups of images. The parallel-tagged image data were acquired with a smaller temporal interval (32.5msec) than the grid-tagged image data (56msec). Thus the algorithm is more likely to lose track of the tag lines in the grid-tagged studies than in the parallel-tagged studies. PS-COTTER on detected tag points has the advantage in this aspect since the algorithm is less a?ected by the fading of the tag lines during the tag tracking procedure compared to PS-COTTER on raw or filtered image data. Third, di?erent imaging protocols were applied for image acquisition. For the parallel-tagged group, six long-axis slices were acquired at each time slot. However, usually a 2-chamber slice and a 4-chamber slice were obtained for each grid-tagged MR study. Some of the grid-tagged studies had an additional long-axis slice. This means the parallel-tagged studies had more longitudinal constraints on the tag tracking procedure than the grid-tagged studies. Therefore, PS-COTTER on raw or filtered image data had better 109 tag tracking accuracy over the parallel-tagged studies than over the grid-tagged studies in the validation experiments. For both parallel-tagged and grid-tagged MRI data, the E cc and E min strains com- puted with PS-COTTER and the user-supervised reconstruction methods TEA/DMF or APSB showed better correlation than E rr and E ll in the midwall. E ll showed reasonable correlation, but not as strong as E cc and E min . E rr correlated weaker because the parallel or grid tag pattern more sparsely samples the radial motion. Due to the large deformation of the apex, the strains correlated stronger at the base and mid-ventricle than at the apex. Based on the experimental results, for parallel-tagged data, closed-apex PS-COTTER on raw image data is recommended. For grid-tagged data, closed-apex PS-COTTER on detected tag points is recommended. In conclusion, the PS-COTTER method can accurately reconstruct deformation and strain in the LV wall from tagged MR images and has several advantages relative to existing techniques. Relative to APSB and the cylindrical B-spline model, the PS-COTTER com- bines tag tracking and motion reconstruction into a single step and strains can be computed directly from the reconstructed displacement field. Relative to COTTER, the PS-COTTER method reconstructed the deformation with less RMS tag distance error. 110 4.8 Appendix The grid-tagged development and evaluation human imaging studies used in this chap- ter: Table 4.17: Development Studies Study ID Visit Study ID Visit Study ID Visit Study ID Visit 05DAC011 BAS 05P3CLN001 BAS 05P1VOL503 BAS 05P1MRR001 BAS 05DAC016 BAS 05P3CLN004 BAS 05P1VOL507 BAS 05P1MRR002 BAS 05DAC017 BAS 05P3CLN010 BAS 05P1VOL510 BAS 05P1MRR006 BAS Table 4.18: Evaluation Studies Study ID Visit Study ID Visit Study ID Visit Study ID Visit 05DAC005 BAS 05P1VOL504 BAS 05P3CLN002 BAS 05P1MRR003 BAS 05DAC022 BAS 05P1VOL506 BAS 05P3CLN005 BAS 05P1MRR007 BAS 05DAC023 BAS 05P1VOL511 BAS 05P3CLN007 BAS 05P1MRR008 BAS 05DAC030 BAS 05P1VOL514 BAS 05P3CLN008 BAS 05P1MRR009 BAS 05DAC031 BAS 05P1VOL516 BAS 05P3CLN009 BAS 05P1MRR002 M06 05DAC033 BAS 05P1VOL517 BAS 05P3CLN011 BAS 05DAC037 BAS 05P1VOL518 BAS 05P3CLN012 BAS 05DAC038 BAS 05P1VOL521 BAS 05P3CLN014 BAS 05DAC039 BAS 05P1VOL524 BAS 05P3CLN015 BAS 05DAC052 BAS 05P1VOL526 BAS 05P3CLN501 BAS 05DAC057 BAS 05P1VOL002 BAS 05P3CLN502 BAS 05DAC060 BAS 05P1VOL002 M06 05P3CLN504 BAS 111 Chapter 5 Conclusion and future work In this dissertation, three new algorithms for tag line tracking and 3-D cardiac LV motion reconstruction were proposed to quantitatively estimate regional myocardium defor- mation from tagged MR images. They are the statistical tag point classification algorithm, the 3-D deformation reconstruction algorithm and the PS-COTTER algorithm. In the tag point classification algorithm, rather than tracking tag lines from frame to frame in an image sequence with active contours or similar techniques, a set of candidate tag points are detected in each image which are then classified as either false positives or belonging to a particular tag line. The advantage of this approach is that the tag point positions are not pre-smoothed during tracking, allowing smoothness constraints to be applied only to the deformation model fitted to the tag points. The 3-D deformation reconstruction algorithm described in Chapter 3 is based on B- spline basis functions in prolate spheroidal coordinates. Two prolate spheroidal B-spline models, PSB and APSB, were developed for this method. This algorithm o?ers several advantages. First, the use of spatially localized B-spline basis functions limits the e?ects of reconstruction artifacts. Second, their domain more closely matches the actual shape of the LV wall. This increases the accuracy of the displacement field reconstruction. Third, the PSB and APSB models can enforce smoothing across and compute strain at the apex. The PS-COTTER algorithm employs the PSB deformation model for combined tag tracking and motion reconstruction. PS-COTTER uses a prolate spheroidal domain en- closed by the ROIs, which are user-selected end-diastolic circular and elliptical regions. 112 The PS-COTTER algorithm conserves the advantages of B-spline basis functions and the prolate spheroidal domain. In PS-COTTER, the myocardium motion model is fit directly to the detected tag points, raw or filtered image data without an intermediate tag tracking step and backward and forward motion field fitting. Future work on the tag point classification algorithm can include post-processing the tracking results to ensure that connected sets of tag line points have the same tag line index. Future work on the 3-D deformation reconstruction algorithm could extend the APSB model to 4-D by adding a B-spline in the time direction in the tensor-product model (Equa- tion (3.17)). An open B-spline can be used if tag line data is only available from end-diastole to end-systole, whereas a closed B-spline can be used if tag line data is available throughout the cardiac cycle. Future work on PS-COTTER could employ di?erent tag point detection or image filtering so as to improve the tag tracking performance. The image registration technique can also be incorporated in PS-COTTER to make the image slices consistent. 113 Bibliography [1] N. Reichek, MRI myocardial tagging, J. Magn. Reson. Imaging 1999, 10:609-616. [2] W. G. O?Dell and A. D. McCulloch, Imaging three-dimensional cardiac function, Annu. Rev. Biomed. Eng. 2000, 02:431-456. [3] P. Suetens, Fundamentals of medical imaging, Cambridge University Press 2002. [4] J. L. Prince and J. M. Links, Medical imaging signals and systems, Printice Hall. [5] T. Taxt, A. Lundervold, J. Strand and S. Holm, Advances in medical imaging, ICPR?98 1998. [6] A. F. Frangi, W. J. Niessen and M. A. Viergever, Three-dimensional modeling for functional analysis of cardiac images: a review, IEEE Trans. Med. Imaging 2001, 20(1):2-25. [7] C. C. Moore, C. H. Logo-Olivieri, E. R. McVeigh and E. A. Zerhouni, Three- dimensional systolic strain patterns in the normal human left ventricle: characteri- zation with tagged MR imaging, Radiology 2000, 214: 453-466. [8] E. R. McVeigh, MRI of myocardial function: motion tracking techniques, J. Magn. Reson. Imaging 1996, 14:137-150 [9] R. J. van der Geest and J. H. C. Reiber, Quantification in cardiac MRI, J. Magn. Reson. Imaging 1999, 10:602-608. [10] J. Declerk, T. S. Denney Jr, C. ? Ozt?urk, W. O?Dell and E. R. McVeigh, Left ventricular motion reconstruction from planar tagged MR images: a comparison, Phys. Med. Biol. 2000, 45:1611-32. [11] T. S. Denney, Estimation and detection of myocardial tags in MR image without user-defined myocardial contours, IEEE Trans. Med. Imag. 1999, 18(4): 330-44. [12] M. A. Guttman, E. A. Zerhouni, and E. R. McVeigh, Analysis and visualization of cardiac function from MR imgages, IEEE Computer Graphics and Applications 1997, 17(1): 30-38. [13] A. A. Young, D. L. Kraitchman, L. Dougherty, and L. Axel, Tracking and finite element analysis of stripe deformation in magnetic resonance tagging, IEEE Trans. Med. Imag. 1995, 14(3): 413-421. 114 [14] Y. Chen and A. Amini, A MAP framework for tag line detection in SPAMM data using Markov randon fields on the B-spline solid, IEEE Trans. Med. Imag. 2002, 21 1110-22. [15] X. Deng, and T. S. Denney Jr., Combined tag tracking and strain reconstruction from tagged cardiac MR images without user-defined myocardial contours, J. Magn. Reson. Imaging, 2005 21(1): 12-22. [16] 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 1995, 195: 829-835. [17] T. S. Denney and E. R. McVeigh, Model-free reconstruction of three-dimensional my- ocardial strain from planar tagged MR images, J. Magn. Reson. Imaging 1997, 7: 799-810. [18] X. Deng and T. S. Denney Jr., Three-dimensional myocardial strain reconstruction from tagged MRI using a cylindrical B-spline model, IEEE Trans. Med. Imag. 2004, 23: 861-7. [19] J. Huang, D. Abendschein, V. G. D?avila-Rom?an and A. A. Amini, Spatio-temporal tracking of myocardial deformations with a 4-D B-spline model from tagged MRI, IEEE Trans. Med. Imag. 1999, 18: 957-972. [20] C. ? Ozt?urk and E. R. McVeigh, Four-dimensional B-spline based motion analysis of tagged MR images: introduction and in vivo validation, Phys. Med. Biol. 2000, 45: 1683-1702. [21] P. Radeva, A. Amini and J. Huang, Deformable B-solids and implicit snakes for 3D localization and tracking of SPAMM MRI data, Comput. Vis. Image Understanding 1997, 66: 163-178. [22] R. Chandrashekara, R. H. Mohiaddin and D. Rueckert, Analysis of 3-D myocardial motion in tagged MR images using nonrigid image registration, IEEE Trans. Med. Imag. 2004, 23: 1245-50. [23] D. Perperidis, R. H. Mohiaddin, D. Rueckert, Spatio-temporal free-form registration of cardiac MR image sequences, Med. Image Anal. 2005, 9:441-56. [24] N. F. Osman, W. S. Kerwin, E. R. McVeigh and J. L. Prince, Cardiac motion tracking using CINE harmonic phase (HARP) magnetic resonance imaging, J. Magn. Reson. Med. 1999, 42(6):1048-60. [25] A. A. Young, Model Tags: Direct 3D tracking of heart wall motion from tagged mag- netic resonance images, Med. Image Anal. 1999, 3:361-372. [26] T. S. Denney Jr., B. L. Gerber, and L. Yan, Unsupervised reconstruction of a three- dimensional left ventricular strain from parallel tagged cardiac images, J. Magn. Reson. Med. 2003, 49(4):743-54. 115 [27] Z. Qian, A. Montillo, D. Metaxas and L. Axel, Segmenting cardiac MRI tagging lines using Gabor filter banks, in Proc. IEEE EMBS Int. Conf. 2003, pp 630-3. [28] H. Frank and S. Globits, Magnetic Resonance imaging evaluation of myocardial and pericardial disease, J. Magn. Reson. Imaging 1999, 10:617-626. [29] H. Chui, L. Win, R. Schultz, J. S. Duncan, and A. Rangarajan, A unified non-rigid feature registration method for brain mapping, Med. Image Anal. 2003, 7(2):113-130. [30] H. Chui and A. Rangarajan, A new point matching algorithm for non-rigid registration, Computer Vision and Image Understanding 2003, 89(2):114-141. [31] J. Montagnat,H. Delingette and N. Ayache, A review of deformable surfaces: topology, geometry and deformation, Image and Vision Computing 2001, 19(14):1023-40. [32] F. L. Bookstein, Principal Warps: Thin-plate splines and the decomposition of defor- mations, IEEE Trans. Pattern Analysis and Machine Intelligence 1989, 11(6):567-85. [33] R. Chandrashekara, R. H. Mohiaddin and D. Ruechert, Analysis of 3-D myocardial motion in tagged MR images using nonrigid image registration, IEEE Trans. Med. Imaging 2004, 23(10):1245-50. [34] A. Ross, S. Dass and A. Jain, A deformable model for fingerprint matching, Pattern Recognition 2005, 38(1):95-103. [35] N. Rougon, C. Petitjean, F. Preteux, P. Cluzel and P. Grenier, A non-rigid registra- tion approach for quantifying myocardial contraction in tagged MRI using generalized information measures, Med. Image Anal. 2005, 9(4):353-75. [36] S. Kaneko, T. Kondo and A. Miyamoto, Robust matching of 3D contours using it- erative closest point algorithm improved by M-estimation, Pattern Recognition 2003, 36(9):2041-7. [37] A. Rangarajan, H. Chui and J. S. Duncan, Rigid point feature registration using mutual information, Med. Image Anal. 1999, 3(4):425-40. [38] S. Gold, A. Rangarajan, C.-P. Lu , S. Pappu and E. Mjolsness, New algorithms for 2D and 3D point matching: pose estimation and correspondence, Pattern Recognition 1998, 31(8):1019-31. [39] N. Paragios, M Rousson. and V. Ramesh, Non-rigid registration using distance func- tions, Computer Vision and Image Understanding 2003, 89(2):142-65. [40] T. McInerney and D. Terzopoulos, Deformable models in medical image analysis, Pro- ceedings of MMBIA?96 171-180. [41] A. Rangarajan, Chui H. and E. Mjolsness, A relationship between spline-based de- formable models and weighted graphs in non-rigid matching, IEEE Conference Com- puter Vision and Pattern Recognition (CVPR), 2001. 116 [42] J. Zhang and A. Rangarajan, A?ne image registration using a new information metric, IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2004. [43] H. Guo, A. Rangarajan, S. Joshi and L. Younes, Non-rigid registration of shapes via di?eomorphic point matching, IEE International Symposium on Biomedical Imaging (ISBI), 2004. [44] J. Rexilius, S. K. Warfield, C. R. G. Guttmann, X. Wei, R. Benson, L. Wolfson, M. Shenton, H. Handels and R. Kikinis, A novel nonrigid registration algorithm and appli- cations, Proceedings of the 4th International Conference on Medical Image Computing and Computer-Assisted Intervention 2001: 923-31. [45] J. Amores and P. Radeva, Registration and retrieval of highly elastic bodies using contextual information, Pattern Recognition Letters 2005, 26(11):1720-31. [46] A.L. Yuille, P. Stolorz and J. Utans, Statistical physics, mixtures of distributions, and the EM algorithm, Neural Comput. 1994, 6(2):334-40. [47] J. Li, C. Davis, T. S. Denney, Jr., H. Gupta, S. Lloyd and L. Dell?Italia, Tag point classification in tagged cardiac MR images, International Symposium on Biomedical Imaging, 2006. [48] C. Davis, J. Li, and T. S. Denney Jr., Analysis of spectral changes and filter design in tagged cardiac MRI, International Symposium on Biomedical Imaging, 2006. [49] J. Ohayon and R. S. Chadwick, E?ects of collagen microstructure on the mechanics of the left ventricle, Biophysical J. 1988, 54:1077-1078. [50] A. A. Young, R. Orr, B. H. Smaill and L. J. Dell?Italia, Three-dimensional changes in left and right ventricular geometry in chronic mitral regurgitation, AM J Physiol. Heart Circ. 1996, 271 H2689-H2700. [51] A. A. Amini, R. Curwen, R. T. Constable and J. C. Gore, Spring Symposium series. Applications of computer vision in medical image processing. 1994, (Menlo Park : AAAI press) p 126 [52] P. E. Gill, W. Murray and M. H. Wright, Practical optimization (London, UK: Aca- demic Press) 1981. [53] C. C. Moore, S. B. Reeder and E. R. McVeigh, Tagged MR imaging in a deforming phantom: photographic validation, Radiology 1994, 190:765-9. [54] E. Atalar and E. R. McVeigh, Optimization of tag thickness for measuring position with magnetic resonance imaging, IEEE Trans. Med. Imag. 1994, 13:152-60. [55] J. Li and T. S. Denney, Evaluation of B-spline cardiac deformation models, ISMRM?05 p1653 [56] J. Li and T. S. Denney, Left ventricular motion reconstruction with a prolate spheroidal B-spline model, Phys. Med. Biol. 2006, 51(3):517-37. 117 [57] D. J. Hawkes, D. Barratt, J. M. Blackall, C. Chan, P. J. Edwards, K. Rhode, G. P. Penney, J. McClelland and D. L. G. Hill, Tissue deformation and shape models in image-guided interventions: a discussion paper, Med. Image Anal. 2005, 9(2):163-175. [58] T. S. Denney and J. L. Prince, Optimal brightness functions for optical flow estimation of deformable motion, IEEE Trans. Imag. Proc. 1994, 3(2):178-91. 118