GEOMETRICMETHODSFOR IMAGEREGISTRATION ANDANALYSIS by AnandArvindJoshi A DissertationPresented tothe FACULTYOF THEGRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA InPartial Fulfillmentofthe RequirementsfortheDegree DOCTOROF PHILOSOPHY (ELECTRICAL ENGINEERING) August2008 Copyright 2008 AnandArvindJoshi Dedication Idedicatethisthesistomyparents Arvindand SuvarnaJoshi. ii Acknowledgments I would like to express my deepest respect and gratitude to my advisor, Prof. Richard Leahy. His contribution has been invaluable; he supported me financially throughout my graduate studies in USC, guided me in selecting and solving research problems, introduced me to distinguished researchers and scientists, provided me a friendly and welcome environment, gave me freedom to select my own research directions and let me enjoy long holidays. He is a great mentor and a source of inspiration for me, and it has really been an honor working with him. I also wish to thank the members of my guidance committee: Dr. Krishna Nayak and Dr. Francis Bonahon for their sugges- tions and valuable feedback. My appreciations go to Dr. David Shattuck and Dr. Paul Thompsonat theUniversityofCalifornia,Los Angelesformanyfruitful collaborations anddiscussions. IcouldnothaveselectedbettercolleaguesthanDr. AbhijitChaudhari, Sangeetha Somayajula, Sanghee Cho, Joyita Dutta, Dr. Sangtae Ahn, Dr. Quanzheng LiandDr. DimitriosPantazis. Wehavespendqualitytimetogether,discussingresearch problems as well as life experiences. I must also mention that Abhijit, Quanzheng and Dimitrios motivated me from time to time to look for research problems. I also would liketo expressmygratitudeto Dr. IlyaEcksteinforfruitfulcollaborations. In theUniversityofSouthern CaliforniaIhaveenjoyedthecozyand warmenviron- ment of an excellent research lab. I want to thank its members: Sangtae Ahn, Abhijit Chaudhari, Sanghee Cho, Belma Dogdas, Hua Hui, Zheng Li, Sangeetha Somayajula, iii Juan Luis Poletti Soto, Evren Asma, YuTeng Chang, David Wheland and Syed Ashra- fullaforaddingpositivelytomyresearch and academicexperience. Most importantly, I want to thank my family for providing unlimited support and helpingmerealize mydreams. iv Table of Contents Dedication ii Acknowledgments iii Abstract xii Chapter 1: Introduction 1 Chapter 2: CorticalSurface Parameterization 5 2.1 Parameterizationand theCoordinateSystem . . . . . . . . . . . . . . . 6 2.1.1 Validationofp-harmonicmappings . . . . . . . . . . . . . . . 9 Chapter 3: CorticalSurface Registration 13 3.1 ThinPlateSplinesRegistrationintheIntrinsicGeometryoftheCortical Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.1 MathematicalFormulation . . . . . . . . . . . . . . . . . . . . 16 3.1.2 DiscretizationAlgorithm . . . . . . . . . . . . . . . . . . . . . 20 3.1.3 BendingEnergy Minimization . . . . . . . . . . . . . . . . . . 21 3.1.4 ValidationTPS surfaceregistration . . . . . . . . . . . . . . . . 24 3.2 AFiniteElementMethodforSimultaneousRegistrationandParameter- ization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2.1 Surface Registration . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.2 MathematicalFormulation . . . . . . . . . . . . . . . . . . . . 28 3.2.3 FiniteElementFormulation . . . . . . . . . . . . . . . . . . . . 29 3.2.4 Resultsand Validation . . . . . . . . . . . . . . . . . . . . . . 33 3.3 OptimumChoiceofSulcal SubsetforRegistration . . . . . . . . . . . . 35 3.3.1 RegistrationError . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.3.2 ProbabilisticModeloftheSulcal Errors . . . . . . . . . . . . . 39 3.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.4 Discussion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Chapter 4: Processing ofDatainthe Surface Geometry 49 4.1 ImageFilteringonSurfaces . . . . . . . . . . . . . . . . . . . . . . . . 50 v 4.2 MathematicalFormulation . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2.1 Isotropicfiltering . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.2.2 Anisotropicfiltering . . . . . . . . . . . . . . . . . . . . . . . . 53 4.3 Discretizationand Numerical Method . . . . . . . . . . . . . . . . . . 55 4.3.1 DiscretizationAlgorithm . . . . . . . . . . . . . . . . . . . . . 55 4.3.2 TheHeat Equationin theIntrinsicGeometry . . . . . . . . . . . 61 4.4 TheHeat Kernelas aPDF . . . . . . . . . . . . . . . . . . . . . . . . . 62 Chapter 5: VolumetricRegistrationusing HarmonicMaps 66 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.2 ProblemStatement andFormulation . . . . . . . . . . . . . . . . . . . 68 5.3 IndirectMappingApproach . . . . . . . . . . . . . . . . . . . . . . . . 70 5.3.1 MathematicalFormulation . . . . . . . . . . . . . . . . . . . . 73 5.3.2 InitializationProcedure . . . . . . . . . . . . . . . . . . . . . . 75 5.3.3 MappingtotheUnitBall B(0,1) . . . . . . . . . . . . . . . . . 75 5.3.4 HarmonicMappingBetween theTwo Brains . . . . . . . . . . 77 5.3.5 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.4 Direct MappingApproach . . . . . . . . . . . . . . . . . . . . . . . . 79 5.4.1 MathematicalFormulation . . . . . . . . . . . . . . . . . . . . 80 5.4.2 HarmonicMapping . . . . . . . . . . . . . . . . . . . . . . . . 80 5.5 VolumetricIntensityRegistration . . . . . . . . . . . . . . . . . . . . . 81 5.5.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.5.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.6 Resultsand Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Chapter 6: Conclusions andFuture Work 94 6.1 GeometricFeatures andManual LandmarksbasedSurface Registration 96 6.2 RegistrationofDTIimages . . . . . . . . . . . . . . . . . . . . . . . . 97 Bibliography 98 vi List of Figures 1.1 ThecorticalsurfaceofthehumanbraindepictedonaMRdata(toprow) andrendered as asurface(bottomrow). . . . . . . . . . . . . . . . . . 2 2.1 Sulcal TracingTools. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Thefigureshowsthecorticalsurfaceanditsmaptoasquare. Thecorpus callosumis constrainedtolieontheboundaryofthesquare. . . . . . . 9 2.3 Thep-harmonicmapsofthelefthemisphereofan individualcortex. . . 10 2.4 The figure shows smoothed histograms for angle distortion and area distortion respectively. In the angle distortion plot, angle distortion increases with the value of p. In the area distortion plot, the distortions for p=4,6,8 are less than that for p=2 and most of the points have small angledistortiononly. Howeverthereisnoobservabletrendforthevalue ofpin eithercase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1 (a) A cortical surface with hand labeled sulci; (b) A flat map of the two cortical surface. The arrows show connectivity at points along the boundary of the square. Due to the spherical topology of the cortical surface,wecanassigntoitacoordinatesystemthatallowsustocompute partial derivatives across the interhemispheric fissure. (c) Chessboard texturemappedto thesurface usingthesquaremaps. . . . . . . . . . . 15 3.2 (upper) The figure shows the warping field computed on the surface. Thedeformationfieldissmoothlyvarying. Thisisachievedbecausethe bending energy regularization was performed in the intrinsic geometry of the surface. The color indicates the magnitude of the deformation. (lower) The thin-plate spline deformation field applied to a regular grid representingleft and righthemispheres. . . . . . . . . . . . . . . . . . 22 vii 3.3 Alignment of the sulcal landmarks: 6 brains are registered to a com- mon cortical surface using their p-harmonic maps in the plane. They are approximately aligned by the p-harmonic maps justifying our small deformation linear model (thin plate bending energy model) which is used for landmark alignment. After applying the covariant TPS defor- mation field to the surface parameterization, we can see that the sulci showbetteralignment. . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.4 (a),(b)Thetwocorticalsurfaceswithhandlabeledsulciascoloredcurves; (c),(d) flat maps of a single hemisphere for the two brains without the sulcalalignmentconstraint;(e)overlayofsulcalcurvesontheflat maps without alignment; (f),(g) flat maps with sulcal alignment; (h) overlay ofsulcalcurveson theflat mapswithalignment. . . . . . . . . . . . . . 30 3.5 RMSerrorand percentage overlapin theflattened map as afunctionofσ. 33 3.6 Mapping of sulcal landmarks from 5 subjects to the atlas brain (left) withoutand (right)withthesulcalalignmentconstraint. . . . . . . . . . 34 3.7 The complete set of candidate sulcal curves from which we select an optimalsubsetforconstrainedcortical registration . . . . . . . . . . . . 37 3.8 (a)Registrationoftwocorticalsurfacesbasedontheflatmappingmethod; (b) Parcellation of the cortex into regions surrounding the traced sulci; (c) Registration error for two corresponding sulci where e (s) are sam- n plesoftheregistrationerror. . . . . . . . . . . . . . . . . . . . . . . . 38 3.9 Samplecovariance matrices for the x, y, and z componentsof the regis- trationerror, represented as colorcoded images. . . . . . . . . . . . . . 42 3.10 Optimal subsets of sulci for cortical registration. Each row gives the indicesoftheoptimalsubsetofsulcithatminimizetheregistrationerror againstall other combinationswith equal numberof constrained curves (also see Fig. 3.7). The three right columns show that the estimated (est.) error is close to the calculated actual (act.) error when actual reg- istrations with the same constrained curves are performed. Our method predicts the registration error both for the training (trn) and the testing (tst)set ofbrains. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.11 Optimalsulcal setsfor5, 10,and 15curves. . . . . . . . . . . . . . . . 45 viii 3.12 Top row: subjectiveselection of 6 curves, with preference on long sulci distant from each other that are expected to minimize cortical registra- tion error; bottom row: optimal sulcal set with the 6 curves selected by ourmethod. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.1 Theimpulseresponseofthe isotropicsmoothingfilters are displayedin the parameter space and on the surface [JSTL05]. It can be seen that when the surface metric is used to compute the Laplace-Beltrami, the impulseresponsekernelisnotisotropicintheparameterspace,however itisisotropiconthesurface. . . . . . . . . . . . . . . . . . . . . . . . 58 4.2 left: The mean curvature of the cortical surface plotted on a smoothed representation (for improved visualization of curvature in sulcal folds; right: The mean curvature plotted in 2D parameter space for a single cortical hemisphere. Isotropic diffusion blurs the regions as well as the edges separating them while while anisotropic diffusion reduces noise whilepreservingedges. . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.3 The figures shows the heat kernels estimated to fit the two datasets for MEG somatosensory data. For each of the datasets the estimated pdf is displayedin theparameterspaceand onthecortical surface. . . . . . . 64 4.4 Theclassifier: Red andBlue regionsshowsthetwo decisionregions . . 64 5.1 Cortical surface alignment after using AIR software for intensity based volumetricalignmentusinga168parameter5th orderpolynomial. Note that although the overall morphology is similar between the brains, the twocortical surfaces do notalignwell. . . . . . . . . . . . . . . . . . . 70 5.2 Illustration of our general framework for surface-constrained volume registration. We first compute the map v from brain manifold (N,I) totheunitballtoformmanifold(N,h). Wethencomputeamapu˜from brain (M,I) to (N,h). The final harmonic map from (M,I) to (N,I) isthen givenbyu = v−1 u˜. . . . . . . . . . . . . . . . . . . . . . . . 74 ◦ ix 5.3 Initialization for harmonic mapping from M to N. First we gener- ate flat square maps of the two brains, one for each hemisphere, with pre-aligned sulci. The squares corresponding to each hemispheres are mapped to a disk and the disks are projected onto the unit sphere. We thengenerateavolumetricmapsfromeachofthebrainstotheunitball. Sinceallthesemapsarebijective,theresultingmapresultsinabijective pointcorrespondencebetween thetwobrains. However,thiscorrespon- dence is not optimal with respect to the harmonic energy of maps from thefirstbraintothesecond,butisusedasaninitializationforminimiza- tionof(5.6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.4 Illustration of the deformation induced with respect to the Euclidean coordinates by mapping to the unit ball. Shown are iso-surfaces corre- spondingtotheEuclideancoordinatesfordifferentradiiintheunitball. Distortions become increasingly pronounced towards the outer edge of thespherewhere theentireconvolutedcortical surface is mappedto the surfaceoftheball. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.5 Schematic of the intensity alignment procedure. Once harmonic maps uM and uN are computed, we refine these with intensity driven warps wM and wN while imposing constraints so that the final deformations areinverseconsistent. . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.6 Illustration of the effects of the two stages of volumetric matching is shownby applyingthe deformationstoa regularmeshrepresenting one slice. Since the deformation is in 3D, the third in-paper value is repre- sented by color. (a) Regular mesh representing one slice in the subject; (b) the regular mesh warped by the harmonic mapping which matches the subject cortical surface to the template cortical surface. Note that deformation is largest near the surface since the harmonic map is con- strained only by thecortical surface; (c) Regular mesh representing one slice in the harmonically warped subject; (d) the intensity-based refine- ment now refines the deformation of the templateto improvethe match between subcortical structures. In this case the deformation is con- strained to zero at the boundary and are confined to the interior of the volume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.7 Examples of direct mapping approach. (a) Original subject volume; (b) original template; (c) registration of subject to template using surface constrainedharmonicmapping,notethatthesurfacematchesthatofthe template;(d)intensity-basedrefinementoftheharmonicmapofsubject totemplateto completeregistrationprocedure . . . . . . . . . . . . . . 88 x