<%BANNER%>

Geometric Approaches to Group-Wise Image Analysis

Permanent Link: http://ufdc.ufl.edu/UFE0044361/00001

Material Information

Title: Geometric Approaches to Group-Wise Image Analysis
Physical Description: 1 online resource (100 p.)
Language: english
Creator: Xie, Yuchen
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2012

Subjects

Subjects / Keywords: image
Computer and Information Science and Engineering -- Dissertations, Academic -- UF
Genre: Computer Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: Group-wise analysis of a set of images is a fundamental problem in computer vision and medical image analysis. A structured representation of the image collection is extremely useful in many important applications such as image registration, segmentation, classification and recognition. In this dissertation, we present a set of techniques which perform group-wise image analysis based on the geometric structure of the data, and we investigate the problems and solutions in two different contexts: low-dimensional subspaces embedded in an Euclidean and non-Euclidean high-dimensional ambient spaces.  In particular, image atlas construction problem can be formulated as computing the mean template by using the intrinsic geometric structure of the image dataset.  In contrast, statistical analysis of tensor images such as diffusion tensor images (DTI) and conductivity tensor images (CTI) requires the ambient space to be non-Euclidean. We will start by presenting a novel method for constructing multiple image templates for a set of 2D or 3D images via an intrinsic geometric approach. Next, we will present a statistical framework for analyzing a group of tensor-valued images. Furthermore, we will focus on finding a small number of bases to represent the tensor-valued images, and we will present one approach based on nonnegative factorization in this dissertation. Finally, we will propose a dictionary learning framework for Riemannian manifolds. All of the above mentioned methods have been experimentally validated and tested on both synthetic and real image data sets.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Yuchen Xie.
Thesis: Thesis (Ph.D.)--University of Florida, 2012.
Local: Adviser: Vemuri, Baba C.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2012
System ID: UFE0044361:00001

Permanent Link: http://ufdc.ufl.edu/UFE0044361/00001

Material Information

Title: Geometric Approaches to Group-Wise Image Analysis
Physical Description: 1 online resource (100 p.)
Language: english
Creator: Xie, Yuchen
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2012

Subjects

Subjects / Keywords: image
Computer and Information Science and Engineering -- Dissertations, Academic -- UF
Genre: Computer Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: Group-wise analysis of a set of images is a fundamental problem in computer vision and medical image analysis. A structured representation of the image collection is extremely useful in many important applications such as image registration, segmentation, classification and recognition. In this dissertation, we present a set of techniques which perform group-wise image analysis based on the geometric structure of the data, and we investigate the problems and solutions in two different contexts: low-dimensional subspaces embedded in an Euclidean and non-Euclidean high-dimensional ambient spaces.  In particular, image atlas construction problem can be formulated as computing the mean template by using the intrinsic geometric structure of the image dataset.  In contrast, statistical analysis of tensor images such as diffusion tensor images (DTI) and conductivity tensor images (CTI) requires the ambient space to be non-Euclidean. We will start by presenting a novel method for constructing multiple image templates for a set of 2D or 3D images via an intrinsic geometric approach. Next, we will present a statistical framework for analyzing a group of tensor-valued images. Furthermore, we will focus on finding a small number of bases to represent the tensor-valued images, and we will present one approach based on nonnegative factorization in this dissertation. Finally, we will propose a dictionary learning framework for Riemannian manifolds. All of the above mentioned methods have been experimentally validated and tested on both synthetic and real image data sets.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Yuchen Xie.
Thesis: Thesis (Ph.D.)--University of Florida, 2012.
Local: Adviser: Vemuri, Baba C.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2012
System ID: UFE0044361:00001


This item has the following downloads:


Full Text

PAGE 1

GEOMETRICAPPROACHESTOGROUP-WISEIMAGEANALYSISByYUCHENXIEADISSERTATIONPRESENTEDTOTHEGRADUATESCHOOLOFTHEUNIVERSITYOFFLORIDAINPARTIALFULFILLMENTOFTHEREQUIREMENTSFORTHEDEGREEOFDOCTOROFPHILOSOPHYUNIVERSITYOFFLORIDA2012

PAGE 2

c2012YuchenXie 2

PAGE 3

Tomyparents 3

PAGE 4

ACKNOWLEDGMENTS Thisdissertationisnotpossiblewithoutthehelp,supportandguidancefrommanypeople.First,Iwouldliketothankmyadvisor,Dr.BabaC.Vemuri,whointroducedmetotheeldofcomputervisionandmedicalimageanalysis.Hisinsight,encouragementsandguidancehavehelpedmethroughmanyobstaclesduringmygraduatestudy.Ifeelsofortunatetohavetheopportunitytoworkwithhim.MydeepgratitudealsogoestoDr.JeffreyHo,whocloselyworkedwithmeinmostofmyprojectsandprovidednumerousinvaluablesuggestions.IamalsothankfultoDr.AnandRangarajan,Dr.ArunavaBanerjeeandDr.SergeiShabanovforbeingonmythesiscommitteeandbeingverygenerousoftheirtimeandsupport.IfeelveryprivilegedtospendtimewithmanywonderfulpeopleintheLabofComputerVision,GraphicsandMedicalImaging(CVGMI).IamheartilygratefultoGuangChengwhohasbeenextremelyhelpfulandpatient.IdeeplyappreciateTingChenforherfriendshipandsupportthroughoutmygraduatestudy.IwanttotakethisopportunitytothankMeizhuLiu,WenxingYe,BingJian,SanthoshKodipaka,AngelosBarmpoutis,RitwikKumar,DohyungSeo,SileHu,QiDeng,YanDeng,YuanxiangWangandotherCVGMIlabmembers.Wehavesharedsomuchjoytogether.IamalsothankfultomywonderfulfriendsinCISEandUF:YingXuan,TaoLi,LuChen,YanLi,MingZhang,ZheWang,LeiYang,XiaokeQin,KunLiandYuChenfortheirhelpandsupport.IamgratefultoDr.XiaodongTaoforgivingmetheopportunitytoworkwithhimasaninterninGEGlobalResearchandforhiscontinuedhelp.IalsoappreciatetheopportunitygiventomebyDr.DebarghaMukherjeetoworkatGoogleinsummer2011.Finallyandmostimportantly,Ithankmyparentsfortheirunconditionalloveandsupport.MygratitudetothemisbeyondwhatIcanexpressinwords.Thisdissertationisdedicatedtothem. 4

PAGE 5

TheresearchforthisdissertationwasnanciallysupportedbyUFAlumniFellowship,NIHgrantsNS046812,EB007082andNS066340toDr.BabaC.Vemuri,andNSFgrantIIS0916001toDr.JeffreyHoandDr.BabaC.Vemuri. 5

PAGE 6

TABLEOFCONTENTS page ACKNOWLEDGMENTS .................................. 4 LISTOFTABLES ...................................... 8 LISTOFFIGURES ..................................... 9 ABSTRACT ......................................... 11 CHAPTER 1INTRODUCTION ................................... 13 2MULTIPLEIMAGEATLASCONSTRUCTION ................... 16 2.1Motivation .................................... 16 2.2RelatedWork .................................. 20 2.3TheoreticalFrameworkforAtlasConstruction ................ 22 2.3.1RotationalInvarianceandtheGeometryofthequotientspaceQ 23 2.3.2GraphrepresentationofthemanifoldM ............... 24 2.3.3Estimatingmeanfrompairwisegeodesicdistances ......... 27 2.3.4Computingtheimageatlases ..................... 29 2.4Experiments .................................. 30 2.4.1ExperimentalResultsonCOIL-20 ................... 31 2.4.2ExperimentalResultsonSyntheticImages .............. 33 2.4.3MultipleatlasesforOASISdata .................... 35 2.4.4Atlasesformaleandfemalesubjects ................. 41 2.4.5Atlasesfordementiapatientsandhealthyelderlysubjects ..... 43 3AGEOMETRICFRAMEWORKFORTENSORFIELDSANALYSIS ....... 45 3.1MotivationandRelatedWork ......................... 45 3.2StatisticalAnalysisintheSpaceofTensorFields .............. 48 3.2.1GeometryofTensorFields ....................... 49 3.2.2StatisticsontheSpaceofTensorFields ............... 50 3.3TensorFieldsClassication .......................... 53 3.4ExperimentalResults ............................. 55 3.4.1Deformationtensorelds ....................... 55 3.4.2Diffusiontensorelds ......................... 57 4NONNEGATIVEFACTORIZATIONOFDIFFUSIONTENSORIMAGES ..... 61 4.1Motivation .................................... 61 4.2Preliminaries .................................. 63 4.3Algorithmandimplementationdetails .................... 65 4.3.1OptimizationwithrespecttothebasismatrixW ........... 65 6

PAGE 7

4.3.2OptimizationwithrespecttocoefcientmatrixH ........... 68 4.3.3Diffusiontensorimagesegmentation ................. 68 4.4Experiments .................................. 70 4.4.1Diffusiontensorimagesegmentation ................. 70 4.4.2Nonnegativefactorizationwithmultipletensorelds ......... 72 5DICTIONARYLEARNINGONRIEMANNIANMANIFOLDS ........... 75 5.1Motivation .................................... 75 5.2RelatedWork .................................. 78 5.3ABriefonRiemannianManifold ........................ 78 5.4DictionaryLearningonRiemannianManifolds ................ 79 5.5SymmetricPositiveDeniteMatrixSpace .................. 83 5.6Square-rootDensityFunctionSpace ..................... 85 5.7Experiments .................................. 87 5.7.1Brodatzdataset ............................. 88 5.7.2OASISdataset ............................. 89 6CONCLUSIONS ................................... 91 REFERENCES ....................................... 93 BIOGRAPHICALSKETCH ................................ 100 7

PAGE 8

LISTOFTABLES Table page 2-1Meansandstandarddeviationsofsubjectagesfortwoclassespartition .... 37 2-2Meansandstandarddeviationsofsubjectagesforthreeclassespartition ... 39 3-1TensorFieldsClassicationonOASIS ....................... 57 5-1ThetextureclassicationaccuracyfordifferentmethodsontheBrodatzdataset. 88 5-2TheclassicationaccuracyfordifferentmethodsontheOASISdataset. .... 89 8

PAGE 9

LISTOFFIGURES Figure page 2-1Intrinsicmeantemplateconstruction. ........................ 18 2-2Semi-supervisedgraphpartition. .......................... 24 2-3Samplesofinputduckimages. ........................... 31 2-4Experimentalresultsforduckimages. ....................... 32 2-5ExperimentalresultsforveobjectsfromCOIL-20. ................ 32 2-6Experimentalresultsonsyntheticimages. ..................... 34 2-7CamerapositionsofinputBunnyimages. ..................... 34 2-8ThehistogramofageforallsubjectsinOASISdataset.40binsareused. ... 35 2-9Sampleimagesfromyouth,middle-agedandelderlysubjects. ......... 36 2-10Twoclassespartition. ................................ 37 2-11TwoatlasesforOASISdata. ............................. 38 2-12Threeclassespartition. ............................... 39 2-13ThreeatlasesforOASISdata. ........................... 40 2-14AtlasesacrossdifferentagesforsubjectsinOASISdataset. .......... 41 2-15Femaleatlasesgeneratedbyoursemi-supervisedalgorithm. .......... 41 2-16Maleatlasesgeneratedbyoursemi-supervisedalgorithm. ............ 42 2-17VolumeratiosofCSF,GMandWMinatlasesforfemalesubjectsandmalesubjects. ....................................... 42 2-18Axialviewofatlasesforoldsubjectsanddementiapatients. ........... 43 2-19Coronalviewofatlasesforoldsubjectsanddementiapatients. ......... 43 3-1Exampleonsynthetictensorelds. ......................... 46 3-2Thespaceoftensorelds. ............................. 48 3-3Thelogmapontheunitsphere. .......................... 53 3-4StatisticalAnalysisfordeformationtensoreldsfromoldagegroup. ...... 55 3-5Comparisonswithvoxel-basedmethod. ...................... 56 3-6Fivesamplesfrominputratbraindiffusiontensorimages. ............ 58 9

PAGE 10

3-7Statisticalanalysisfordiffusiontensorimages. .................. 58 3-8Diffusiontensorimagereconstruction. ....................... 59 3-9Reconstructionerrorsbetweenthetesttensoreldanditsprojectionsonthegeodesicsubmanifoldsofdifferentdimensions. .................. 59 4-1Segmentationexperimentonsynthetictensoreld. ................ 70 4-2Diffusiontensorimagesegmentation. ....................... 71 4-3Nonnegativefactorizationonmultiplesynthetictensorelds. .......... 72 4-4Nonnegativefactorizationonrealdiffusiontensorimages. ............ 74 5-1DictionarylearningonRiemannianManifold. ................... 80 5-2SamplesofBrodatztexturesusedinourexperiments. .............. 88 10

PAGE 11

AbstractofDissertationPresentedtotheGraduateSchooloftheUniversityofFloridainPartialFulllmentoftheRequirementsfortheDegreeofDoctorofPhilosophyGEOMETRICAPPROACHESTOGROUP-WISEIMAGEANALYSISByYuchenXieAugust2012Chair:BabaC.VemuriMajor:ComputerEngineering Group-wiseanalysisofasetofimagesisafundamentalproblemincomputervisionandmedicalimageanalysis.Astructuredrepresentationoftheimagecollectionisextremelyusefulinmanyimportantapplicationssuchasimageregistration,segmentation,classicationandrecognition.Inthisdissertation,wepresentasetoftechniqueswhichperformgroup-wiseimageanalysisbasedonthegeometricstructureofthedata,andweinvestigatetheproblemsandsolutionsintwodifferentcontexts:low-dimensionalsubspacesembeddedinanEuclideanandnon-Euclideanhigh-dimensionalambientspaces.Inparticular,imageatlasconstructionproblemcanbeformulatedascomputingthemeantemplatebyusingtheintrinsicgeometricstructureoftheimagedataset.Incontrast,statisticalanalysisoftensorimagessuchasdiffusiontensorimages(DTI)andconductivitytensorimages(CTI)requirestheambientspacetobenon-Euclidean. Wewillstartbypresentinganovelmethodforconstructingmultipleimagetemplatesforasetof2Dor3Dimagesviaanintrinsicgeometricapproach.Next,wewillpresentastatisticalframeworkforanalyzingagroupoftensor-valuedimages.Furthermore,wewillfocusonndingasmallnumberofbasestorepresentthetensor-valuedimages,andwewillpresentoneapproachbasedonnonnegativefactorizationinthisdissertation.Finally,wewillproposeadictionarylearningframeworkforRiemannianmanifolds.Allof 11

PAGE 12

theabovementionedmethodshavebeenexperimentallyvalidatedandtestedonbothsyntheticandrealimagedatasets. 12

PAGE 13

CHAPTER1INTRODUCTION Findinganeffectiverepresentationtomodelacollectionofdataisoneofthecentraltopicsinstatisticalanalysis.Ithasalsobecomeimmenselyimportantinvariouscomputervisionandmedicalimagingapplications.Forinstance,Eigenfaceperformstheprincipalcomponentanalysis(PCA)[ 68 ]tondalow-dimensionallinearsubspaceformodelingalargesetoffaceimagesandusesitinfacerecognition.Activeshapemodel(ASM)[ 14 ]andactiveappearancemodel(AAM)[ 13 ]buildastatisticalmodelofobjectshapeandappearanceduringthetrainingphase.Themodelcanbeusedasapriorinlandmarksalignment,objectsegmentationandtracking.Inmedicalimageanalysis,therstorderstatistics(mean)ofimagesfromdifferentsubjectsisusuallyconsideredasanunbiasedtemplate(atlas)thatcanbeusefulinmedicaldiagnosis,particularlywithassociatedclinicalinformation. Theworkpresentedinthisdissertationaimstoconsidertheunderlyinggeometricstructureofthedatawhenmodelingasetofimages.Incomputervisionandmedicalimaging,imagesareoftenmodeledaspointsinahigh-dimensionalvectorspace(imagespace).Forexample,asmallimagewiththesize100100isconsideredasapointina104dimensionalimagespaceaftertheusualrasterization.Inrecentyears,manifoldlearningtechniquessuchasIsomap[ 34 ],LLE[ 53 ]andLaplacianEigenmaps[ 6 ]havegottenalotofattentionincomputervisionandmachinelearningcommunities.Thesemethodsassumethattheimagesofinterestlieonanembeddednon-linearspacesituatedinahigh-dimensionalambientspace.Giventhisassumption,it'smoreappropriatetoperformstatisticalanalysisontheembeddednon-linearspaceinsteadoftheambientimagespace.Forinstance,themeanofasetoffaceimagesshouldalsobeafaceimageinsteadofsomeimagethatremotelyresemblesaface,i.e.,itfallsoutsidethespaceoffaceimages.Extendingthisideatoimagetemplateconstruction,we 13

PAGE 14

proposeanovelapproachtoconstructmultipleimageatlasesviatheintrinsicaveragingontheembeddedspaceofthedataset. Furtherextensionscanbemadebyassumingthattheambientimagespaceitselfisanonlinearmanifold,andadiffusiontensorimage(DTI)offerstherstexample.Itusesasymmetricpositive-denitematrix(SPD)withineachvoxeltomodelthediffusionofwatermoleculesintissues,andeachdiffusiontensorimagecanbeconsideredasanSPDmatrixeld.However,thespaceofSPDmatricesisnotavectorspace,andthestatisticalmethodsdevelopedforvectorspacecannotbedirectlyusedfortensorimages.Toremedythis,wegeneralizethreepopularmethodsinthisthesis:principalcomponentanalysis(PCA),nonnegativematrixfactorization(NMF)[ 36 ]anddictionarylearningtoanalyzeagroupoftensorvalueddata.Specically,weconsideratensorimageasapointintheproductofSPDmatrixspaces.Thewell-knownPCAparadigmcanbeextendedtotheproductspacebyrstcomputingtheKarchermeanoftensorimagesusinganiterativemethod.Theclassesoftensorimagesarethenrepresentedbysubmanifoldsspannedviaasetofprincipalgeodesics,andthisapproachallowsustodirectlyextendthewell-knownclassicationframeworksuchasEigenfacetotensoreldsclassication.SimilartoPCA,NMFalsorepresentsthedataasalinearcombinationofbases.InsteadoftheorthonormalconstraintsplacedonthebasisvectorsinPCA,NMFrequireselementsinbothcoefcientsandbasestobenonnegative.Astheintensityvaluesinimagesarenonnegative,NMFhasbeenconsideredasmorenaturalforimagedatathanPCA.Inourgeneralizationofnonnegativefactorization,eachbasisisavectorofpositivesemi-denitematrices,wherenonnegativityessentiallycomesfromtheconestructureunderlyingthespaceofSPDmatrices.Thisnovelmethodoffersaneffectiveapproachtosimultaneouslyanalyzealargenumberofdiffusiontensorimages,anessentialtaskingroupstudiesandotherapplications.Dictionarylearningconsidersamodeltorepresentthedataasseverallowdimensionalsubspaces.Thedictionarytrainedfromacollectionofimagesisaset 14

PAGE 15

ofbasisimages(generators)andeachimagesamplecanbeeffectivelyapproximatedbychoosingonlyasmallnumberofgeneratorsfromthedictionary.Thistechniquehasbeensuccessfullyusedinimagerestoration,reconstructionandobjectrecognition.Inthisthesis,weproposeageneraldictionarylearningframeworkfordataonRiemannianmanifoldssuchasSPDmatrixspaceandhypersphere. Therestofthisdissertationisorganizedasfollows.Chapter2describesamultipleimageatlasconstructionmethod.InChapter3andChapter4,wegeneralizethePCAandNMFfortensorimagesanalysis.InChapter5,wepresentadictionarylearningapproachonRiemannianmanifolds.Finally,wesummarizethemaincontributionsofthisdissertationanddiscussafewthoughtsforfutureworkinChapter6.Partoftheworkpresentedinthisdissertationhasbeenpublishedin[ 80 82 ]. 15

PAGE 16

CHAPTER2MULTIPLEIMAGEATLASCONSTRUCTION 2.1Motivation Atlasesastheinformativerepresentativesofapopulationofimageshavebeenwidelyusedinmanycomputervisionandmedicalimagingapplicationssuchastemplate-basedimagealignment[ 64 ],atlas-basedimagesegmentation[ 5 21 65 73 ]andstatisticalanalysisacrosssubjects[ 3 30 ].Inaddition,forheterogeneousorlongitudinaldatasets,multipleatlasesareusuallyrequiredtoprovideaninformativeandcompleterepresentationoftheimagedata.Forexample,theagesofthesubjectsinOpenAccessSeriesofImagingStudies(OASIS)database[ 43 ]rangefrom18to96.Sincethestructuraldifferenceinthebrainacrossdifferentagegroupscanbesignicant[ 44 ],it'snotappropriatetoconstructasinglerepresentativeatlasimagefortheOASISdata.Notsurprisingly,algorithmsandmethodsforcomputingmultipleatlasesfromimagecollectionshaveattractedconsiderableamountofattentionrecently,andinparticular,theunbiaseddiffeomorphicatlasconstructionalgorithmproposedin[ 32 ]hasbeeninuentialforseveralnotablerecentdevelopmentssuchastheiClusteralgorithmproposedin[ 54 ]. TheiClusteralgorithmcomputesmultipleatlasesbyttingaGaussianmixturemodeltotheinputimages.Theresultingexpectation-maximizationalgorithm(EM)iterativelygroupstheimagesintodifferentclusters(theE-step)andcomputesasingleatlasforimagesineachclusterbyaveragingthemultiplealignedimages,anapproachthatissimilarinspirittotheunbiaseddiffeomorphicatlasconstructionalgorithm.Figure 2-13 displaysseveralatlasimagesproducedbytheiClustermethod.Acommonandvisiblefeatureoftheseatlasimagesistheircharacteristicblurrinessandtheapparentabsenceofclearandsharpanatomicaldetails.Inshort,withtheircharacteristiclackofimagedelity,thecomputedatlasesdonotcomeacrossasrealspecimenofMRIimagesofhumanbrains,ashortcomingthatcouldbecriticalinsome 16

PAGE 17

applications.Forinstance,in[ 78 ],ithasbeenarguedthataclearandsharpatlasisimportantfortheaccuracyofatlas-basedmultipleimageregistration.Intheclinicalsetting,acharacteristicallyclearandsharpatlasimageasaninformativerepresentativeofasubpopulationofpatientsiscertainlymoremeaningfulandaccessibletomedicalpractitionersthanafuzzyatlasimagewithoutclearlyidentiableanatomicaldetails. Thislossofimagedetailsandstructurescanbetracedbacktotheatlasconstructionstep(M-step)oftheiClusteralgorithminwhichtheatlasiscomputedasthearithmeticaverageofgroup-wiseregisteredimages,apopularatlas-constructionparadigmrstintroducedin[ 32 ].Itisourobservationthattheunbiaseddiffeomorphicatlasconstructionalgorithmoftenproducesblurryandfuzzyatlasimagesunlesstheinputimagesarehighlyhomogeneousintheirstructures.Unfortunately,thisisusuallynotthecaseandthecomputedatlasesalmostalwayssufferfromthelossofimagedetailsinvariousdegrees.Thisundesirableshortcomingcanbeexplainedinbothpracticalandtheoreticalterms.Inpractice,theexistingnumericalalgorithmsusedforcomputinggroup-wiseregistrationcannotguaranteeagloballyoptimalsolution,andunlesstheimagesareunusuallysimilar,non-negligibleregistrationerrorscannotbeavoidedingeneral.Coupledwiththesuboptimalchoicesofregularizationconstantsthatareoftendifculttodetermine,itisnotdifculttoimagineablurryandfuzzyatlascomputedasthepixel-wisearithmeticmeanofmanyandinaccuratelygroup-wiseregisteredimages.ThesecondandmoregeometricexplanationcanbeofferedasinFigure 2-1 .Thesmoothvariationoftheinputimagesallowsustoassume(e.g.,[ 16 ])thattheseimagesbelongtosomesmoothsubmanifoldMintheambientimagespaceI.Unfortunately,thismanifoldstructureisusuallydifculttomodeldirectlyandperhapsmoreimportantly,itisdifculttoutilizeinthegroup-wiseregistrationframework.Therefore,inessence,theunbiaseddiffeomorphicatlasconstructionalgorithmcomputesthemeanoftheimagesintheambientimagespaceI,andasthegureshows,thisdoesnotguaranteethatthe 17

PAGE 18

Figure2-1. Tenpointsonanone-dimensionalmanifold.Left:Meanforthetenpointscomputedusingthemetricoftheambientspace(iCluster).Right:MeancomputedusingthemetricofthesubmanifoldM(Theproposedmethod).Noticethegreaterpreservationofimagefeaturesandstructuresoftherightatlas. resultwillbeonornearthesubmanifoldM.Thisresultsinthelossofcommonimagefeaturesandrenderstheatlasimagesoftenblurryandfuzzy. Inthecontextofmultipleatlasconstructionfromalargeimagecollection,weproposeanovelatlasconstructionalgorithmbasedontwonovelconceptualimprovementsthataimtocorrecttheshortcomingsdescribedabove.First,unliketheEMalgorithmproposedin[ 54 ],wewillexplicitlydecoupletheclusteringandatlasconstructionsteps.Thisdecouplingprovidesuswithgreaterexibilityinapplyingmanifoldlearningmethods(e.g.,[ 6 34 ])tomodelthemanifoldM,anobjectivethatisdifculttoaccomplishusingthetraditionalEM-basedgaussianmixturemodelasin[ 54 ].Inparticular,wewillrequirethattheatlasimagesarerotationinvariantinthesensethatsupposeIistheatlasfortheimagecollectionfI1,,Ingandifndifferentrotationsgiareappliedtotheseimagestoformanewcollectionfg1(I1),,gn(In)g,thentheatlas~IforthelattercollectionshouldberelatedtoIbysomerotation.ThistechnicalrequirementmeansthatwecannolongerassumethattheinputimagesbelongtosomesubmanifoldMintheimagespaceI.Instead,thecorrectassumptionshouldbethatmodulorotations,theimagesbelongtoasubmanifoldMinan(abstract)quotientspaceQofI[ 45 46 ],andmanifoldlearningprovidesasuitablecontextandnaturalsolutionformodelingMinthisabstractsetting.Furthermore,itsexibilityalsoallowsustoreadilyincorporatelabeledimagesinasemi-supervisedcontext,againanobjectivethatwouldbeawkward 18

PAGE 19

toformulateusingtheframeworksproposedearlier(e.g.,[ 32 54 ]).Second,wewillabandontheusualparadigmofcomputingtheatlasimageasthearithmeticmeanofallgroup-wiseregisteredimages.Instead,followingthemanifoldassumption,wewillproposeanewparadigmforatlasimageconstructionbycomputingtheatlasastheweightedarithmeticmeanofafewgroup-wiseregisteredimages.Figure 2-1 illustratesthemainideathatonthemanifoldM,themean(atlas)shouldbeentirelydeterminedonlybyitsneighbors,anditscomputationcouldbeformulatedinatwo-stepprocessthatrstdeterminestheneighborsoftheabstractmeanonthemanifoldM(localization)andthenrendersthemeanasanatlasimageusingtheseneighbors(realization).Inparticular,thesetwostepsutilizetheinputimagesdifferently:inthelocalizationstep,allimagesareusedtodeterminetheneighborsoftheabstractmeanonMwhileinthefollowingrealizationstep,onlytheseneighborsareusedtorendertheatlasimagethroughweightedgroup-wiseregistration.Tothebestofourknowledge,theexplicitidenticationandinvestigationofthesetwoseparateandcriticalprocessesinatlasconstructionhavenotbeenreportedorproposedintheliterature. Theproposedmultipleatlasconstructionalgorithmincorporatestheideasdiscussedabove:Ak-nearestneighborgraph(k-NN)isusedtomodelthemanifoldstructureofM[ 6 34 ],andthegraphpartitionalgorithmisappliedtothek-NNgraphtocomputeclustersoftheinputimages.EachatlasiscomputedasthemeanmoftheimagesbelongtoeachclusteronM.Ingeneral,themetriconMisunknown;however,thepairwisegeodesicdistancesdij=dM(Ii,Ij)betweenimagescanbeestimatedusingthemetricinthequotientspaceQ.AnimportantelementintheproposedalgorithmistoestimatethegeodesicdistancesdM(m,Ii)betweentheabstractmeanmonMandimagesIiusingthepairwisegeodesicdistancesdijastheinputs.Weformulatethisasaconvexoptimizationproblem,andthesolutiontothisoptimizationproblemwillprovideanapproximationtodM(m,Ii).ThiswillallowustolocatethemeanmonthemanifoldMinthelocalizationstepbydeterminingitsneighbors(givenathreshold).In 19

PAGE 20

therealizationstep,weusetheweightedgroup-wiseregistrationmethodtoformtheatlasimagesineachclustergiventheestimatedneighborsoftheatlas.SeveralatlasescomputedbytheproposedmethodareshowninFigure 2-13 ,andincomparison,theyretainconsiderableamountofimportantimagefeaturesandmuchbetterpreservationofimagestructures.AnextensiveevaluationandvalidationoftheproposedalgorithmusingOASISdatawillbereportedlater,andtheexperimentalresultsonOASIShaveshownthat,withgreaterpreservationofimportantimagefeaturesanddetails,theimageatlasesproducedbytheproposedalgorithmenjoybetterimagesharpnesswhencomparedwithatlasesobtainedusingexitingmethods. 2.2RelatedWork ImageatlasasastandardtemplatehasbeenwidelyusedforthediscoveryofdensecorrespondencesbetweenMRIimagesofdifferentsubjects.Thusitisnaturaltoconstructtheimageatlasbasedongroup-wiseimageregistration(e.g.,[ 32 ]).Inthistypeofapproach,theatlasformationiscastintoaniterativeenergyminimizationframework:inputimagesarenon-rigidlyregisteredtoatemplateinthecommondomainandthetemplateisupdatedbyaveragingthealignedimages.Unfortunately,becauseofthelargenumberofvariablesandthenon-convexcostfunctional,thisalgorithm,similartomanyothers,cannotavoidthepricklytrapsoflocalminima.Moreover,therigidityoftheframeworkmakesitdifculttoincorporateotherusefulmodelsandinformationinthecomputation.Asalludedtointheprevioussection,theatlasgeneratedbythismethodessentiallycomputestheatlaswithrespecttothetheambientimagespacewhiledisregardingthepotentiallyusefulassumptionofanimpliedmanifoldstructurebytheinputimages. Inrecentyears,manifoldlearningtechniqueshavepercolatedintothemedicalimageanalysiscommunity,andtheyhavefoundsuccessesinseveralrecentworks.[ 26 ]performedbrainpopulationanalysisbyprojectinginputbrainimagesontoalow-dimensionalspaceusingIsomap[ 34 ].[ 29 ]proposedaframeworkforregistering 20

PAGE 21

imagestotheatlasonanatomicalmanifolds.Theirempiricalmanifoldisconstructedfrominputimagesasak-NNgraph.Asamplefromtheinputpopulationwiththesmallestsumofsquaresofgeodesicdistancestoallinputimagesisselectedastheatlas.[ 31 ]introducedahierarchicalgroupwiseregistrationframeworkusingthek-NNIsomapthatprovidestheintrinsicstructureoftheinputimagedatasetfortheiralgorithm.However,allthemethodsdiscussedaboveconstructasingleatlasfortheentireinputimageset,whichimplicitlyassumesahomogeneouspopulation,anassumptionthatisinappropriateforstudyingheterogeneousimagedatasetsthatarefarmorecommonandimportantinmedicalimaginganalysis.iClusterisanEM-basedalgorithmpresentedby[ 55 ]forcomputingmultipletemplateimagesforanimagepopulation.ThealgorithmtsaGaussianmixturemodeltotheinputimagesandtheinputimagesarenotassumedtopossessanintrinsicmanifoldstructure.Furthermore,itiswell-knownthattheoptimizationintheEM-algorithmcannotguaranteethegloballyoptimalsolution,andtherefore,goodinitializationsoftheEM-algorithm,whichmaynotbeavailableorevenpossible,arecriticalforthesuccessofthemethod. ThreepertinentfeaturesdistinguishingourmethodfromtheiCluster[ 55 ]andmanypreviouslypublishedmethods.First,theproposedmethodexplicitlymodelstheintrinsicmanifoldstructureoftheinputimagesusingmanifoldlearning(k-NNgraph)andthemanifoldmodelingtakesplaceintheabstractquotientspaceinsteadoftheusualimagespace.Second,ourmethoddecouplestheclusteringstepfromtheatlasconstructionstepthatallowstheapplicationofmorepowerfulclusteringalgorithmtoproducemorehomogeneousclusters.Finally,asmentionedinthemotivation,weexplicitlyidentifytwoimportantandseparatestepsintheatlasformationprocess:localizationandrealization.Weformulateaconvexcostfunctioninthelocalizationstepthatbypassesthedifcultyofgloballyminimizingthecomplicatednon-convexfunctional,andintherealizationstep,wegroup-wiseregisteronlyafewimages,ataskthatcanbeaccomplishedefcientlyandaccuratelyduetoboththenumberofimagesaswellastheirhomogeneity. 21

PAGE 22

2.3TheoreticalFrameworkforAtlasConstruction Thissectionpresentstheproposedalgorithmforcomputingmultipleatlasesfromacollectionofimages.Theprimaryaimofouralgorithmistoproducesharpatlasimagesthatretainclearimportantanatomicalstructurescommontothesubpopulationofimagestheyrepresent.Furthermore,wealsorequirethatthecomputedatlasesareinvariantwithrespecttoimagerotations,animportantaspectoftheproposedmethodthatwillbeelaboratedlater.Thestrategyforachievingthesetwosomewhatdisparategoalsistoutilizethenotionofmanifoldstructureimpliedbytheimages.Specically,wewillassume,asinmanymanifold-basedlearningmethods(e.g.,[ 7 88 ]),thattheinputimagesaresamplesfromanunknownmanifoldM,andthemanifoldstructureisthenmodeledusingagraphGcomputedfromtheinputimagesandsomemetric(e.g.,similarity)informationamongtheimages.TherotationalinvariancerequirementimpliesthatthemanifoldofinterestMshouldbeconsideredasasubmanifoldofaquotientspace(quotientbytherotations),andcomputationally,thisrequiresrotationally-invariantmetricforcomputingthegraphG.ThegraphGallowsustopartitiontheimagecollectionintosubcollections,inbothsupervisedandsemi-supervisedfashion,suchthatanatlascanbecomputedfromeachsubcollectionofimagesusingonlymetricalinformationbetweentheimages. Beforedelvingintothedetails,wewillxthenotationsthatwillbeusedthroughoutthefollowingdiscussion.LetIdenotethespaceofimages,andwedeneimagesinIformallyasL2-functionsonaniteimagedomainIRd(d=2for2Dimagesandd=3for3Dimages).LetCfI1,,IngdenoteaninputcollectionofnimagesandAfIa1,,IatgthecollectionoftatlasesfortheimagecollectionC.Arotationg2SO(d)intheimagedomaintransformsanimageIintoarotatedimageg(I)accordingtotheformulag(I)(x)=I(g)]TJ /F5 7.97 Tf 6.58 0 Td[(1(x))forx2andI2I. 22

PAGE 23

2.3.1RotationalInvarianceandtheGeometryofthequotientspaceQ Therotationinvariancepropertydescribedabovecanbeformulatedpreciselyasfollows.LetCdenotetheinputimagecollectionandC0fI01,,I0ngareimagesobtainedbyapplyingnrotationstoimagesinC:I0i=gi(Ii)forgi2SO(d),1in,theatlasesA0fI0a1,,I0atgcomputedforC0aredifferentfromtheatlasesinAuptorotations,i.e.,I0aj=rj(Iaj)forrj2SO(d),1jt.TherotationalinvariancepropertyrequiresustoworknotintheimagespaceIbutinitsquotientspaceQ,thequotientspaceofIbytherotationgroupSO(d)[ 45 46 ].ThequotientspaceQisaspacethatparameterizestheSO(d)-orbitsinIandinthiswork,wewillassumethatQhasamanifoldstructureinducedfromthemanifoldstructureofI.Let:I!QdenotethecanonicalprojectionmapthatsendseachimageI2ItotheuniqueSO(d)-orbit[I]2QcontainingI:(I)=[I]. Thereisawell-knownbijectivecorrespondence[ 45 ]betweenmetricsonQandSO(d)-invariantmetriconI.RecallthataSO(d)-invariantmetriconIsatisesthefollowingconditiondI(I1,I2)=dI(g(I1),g(I2)), foranytwoimagesI1,I22Iandg2SO(d).ThecorrespondencebetweenmetricsdQonQandSO(d)-invariantmetricsdI(x,y)onIisprovidedbytheformula,[ 45 ] dQ([I1],[I2])=minx2[I1],y2[I2]dI(x,y),(2) foranytwopoints[I1],[I2]2Q.Notethat[I1],[I2]arerealizedinIasSO(d)-orbits,anddQsimplycomputesthedistancebetweenthetwoorbitsinIasmeasuredbydIinI. ThetworelatedmetricsdI,dQallowsustogobetweenIandQ.Inparticular,anycomputationrelatingtoametriconQcanbeequivalentlyformulatedusingitscorrespondingSO(d)-invariantmetriconI.Fortunately,manymetricsinIcanbeeasilyshowntobeSO(d)-invariant,suchastheusualL2-metricandthefollowingmetric 23

PAGE 24

Figure2-2. Semi-supervisedgraphpartition.Mistheunderlyingmanifold.Redandbluenodesarelabeledverticesfortwodifferentclasses.Greennodesareunlabeledvertices.Thedashedcurveisthecutthatgeneratesthepartition(subgraphs). proposedin[ 16 ]: dI(I1,I2)=minhiZjI1(hi(x)))]TJ /F3 11.955 Tf 11.95 0 Td[(I2j2dx+Z10ZkLvik2dxds1=2(2) whereviisthetime-dependentvectoreldthatdenesthediffeomorphicowfromtheidentitytothediffeomorphismhi,andLisasecond-orderellipticoperator.GiventheimagesIiinC,themanifoldMthatwillbemodeledpertainstothepointsf[I1],,[In]ginthequotientspaceQ.AsasubmanifoldofQ,MisnaturallyequippedwiththeinducedRiemannianmetric,andinthefollowing,wewilldenotedM([Ii],[Ij])thegeodesicdistancefunctiononM. 2.3.2GraphrepresentationofthemanifoldM BecausethemanifoldMdoesnothaveexplicitrepresentation,wefollowthecommonapproachthatusesagraphG=(V,E)tocharacterizeitsmanifoldstructure,whereV,Earethenodeandedgesets,respectively.Inthisconstruction,thenodesetV=f[I1],,[In]gisprovidedbythepoints[Ii]inthequotientspaceQ.Theedgeset 24

PAGE 25

EthatdenestheconnectivityofGisdeterminedbythepair-wisedistancesamongthepoints[Ii]viathestandardk-nearestneighbor(k-NN)construction.Thatis,anedgeisformedbetweeneverynodeanditsknearestneighbors.AsimpleexampleisshowninFigure 2-1 ,andinthisexample,eachnodeisconnectedtoitstwonearestneighbors(k=2)andtheresultinggraphgivesagoodapproximationtotheunderlyingcurve. Becauseoftheheterogeneityofsomeimagesets,wepartitionthegraphGintotsubgraphsinordertocomputemultipleatlasestorepresenttheimagesets.SpectralclusteringandrelatedmethodsprovideareadilyavailablealgorithmsforpartitioningthegraphG.PopularandWell-knownalgorithmsincludeRatioCut[ 28 ]andNormalizedCut[ 57 ],andourunsupervisedgraphpartitionalgorithmbasedonNormalizedCutisoutlinedinAlgorithm1.ThealgorithmpartitionsthegraphusingthegraphLaplacianmatrixLthatiscomputedfromthesimilaritymatrixWthatencodesallthemetricalinformationprovidedbythenodes[Ii]: Wij=8>><>>:exp()]TJ /F5 7.97 Tf 10.49 6.11 Td[(dQ([Ii],[Ij])2 2)ifj2Niori2Nj,0otherwise,(2) whereNiandNjareknearestneighborhoods.Theparameterisempiricallyestimatedby=1 nPni=1dQ([Ii],[Iik])where[Iik]isthek-thnearestneighborof[Ii].ThegraphLaplacianmatrixisdenedas L=D)]TJ /F3 11.955 Tf 11.95 0 Td[(W(2) whereDisanndiagonalmatrixwithDii=Pnj=1Wij. Afterthepartition,weobtaintsubgraphsGi=(Vi,Ei),whereV1,...,VtisthepartitionofthesetVandEiEonlyincludesedgeswhosetwoverticesareinVi. Theaboveunsupervisedframeworkcanbeeasilyextendedtoasemi-supervisedframeworkwithpartiallylabeleddata.Thisextensionisrelevantbecauseinmanymedicalimagingapplications,asmallnumberofdatalabeledbyexpertsareusuallyavailable,andtheyshouldbeutilizedtoimprovetheanalysisandclassicationofthe 25

PAGE 26

Algorithm1Unsupervisedgraphpartitionbasedonspectralclustering 1: ComputesimilaritymatrixWandgraphLaplacianmatrixL. 2: Solvethegeneralizedeigenvectorproblem,Lu=Du. Letu1,...,utbeteigenvectorscorrespondingtothesmallesteigenvaluesandformthematrixU=[u1,u2,...,ut]2IRnt. 3: Letyi2IRtbethei-throwofU,wherei=1,...,nandclusterthepointsy1,,ynintotclustersusingK-meansalgorithm. remainingunlabeleddata.Inspiredbytheideasfromsemi-supervisedlearning[ 7 88 ],weconsiderthegraphpartitionproblemwithpartiallylabelednodesasaclassicationproblemthattransferstheknownlabelstotheunlabeleddata.Atwo-classexampleisshowninFigure 2-2 .Redandbluenodesarelabeledverticesfortwodifferentclasses.Greennodesareunlabeledverticeswhichwillbeclassiedusingtransductivelearningalgorithmsthatgivesthetwosubgraphswithredandbluenodesandtheiredges,respectively. Wewillconsiderthetwo-classpartitionproblemrstwiththelabelsetY=f)]TJ /F4 11.955 Tf 15.28 0 Td[(1,1g.Theclassierisdenedasafunction:f:M!Y.Becauseanyfunctionf2L2(M)canbewrittenintermsofeigenfunctionsoftheLaplacian,f(x)=P1i=1aiui(x),whereuiareeigenfunctions:ui=iuiandistheLaplace-Beltramioperator.Thuswecouldapproximatetheclassieronthebasisoflabeleddataanduseittoclassifyunlabeleddata.LetthesamplesetbeX=fx1,...,xs,xs+1,...,xngM.Withoutlossofgenerality,weassumetherstsdataarelabeledandothersareunlabeled.Thelabelsetisfc1,...,csg,whereci2Y.GiventhegraphLaplacianL,weneedtosolvetheeigenvectorproblemLu=u.Letu1,...,upbeeigenvectorscorrespondingtothepsmallesteigenvalues.ui=[ui1,...,uin]T,wherei=1,...,p.Thenwecomputetheapproximatedclassierfromthelabeleddatabyminimizingthefollowingcostfunction E(a)=sXi=1 ci)]TJ /F7 7.97 Tf 18.27 15.2 Td[(pXj=1ajuji!2=kUa)]TJ /F3 11.955 Tf 11.95 0 Td[(ck2,(2) 26

PAGE 27

wherethecoefcientsa=[a1,...,ap]T,thelabelsc=[c1,...,cs]TandthebasismatrixU=0BBBB@u11up1.........u1sups1CCCCA.Theclosed-formsolutionofthisoptimizationproblemisa=(UTU))]TJ /F5 7.97 Tf 6.59 0 Td[(1UTc.Fortheunlabeleddataxi,i>s,wesetthelabeltobe1ifPpj=1ajuji0,otherwisethelabelissetto-1.Fort-classesgraphpartition,t>2,wefollowtheone-versus-the-restapproachwhichiscommonlyusedinmulticlasssupportvectormachines[ 72 ].Webuildtseparatetwo-classclassiersfk,wherek=1,...,t.Thek-thclassierisestimatedbyusingthedatafromclass-kasthepositivesamplesandthedatafromremainingt)]TJ /F4 11.955 Tf 12.47 0 Td[(1classesasthenegativesamples.Thentcoefcientsa1,...,atarecomputedandwillbeusedtoclassifyunlabeleddata,whereak=[ak1,...,akp]T.Forxi,i>s,thelabelissettobeargmaxkfk(xi)=Ppj=1akjuji.ThistransductivegraphpartitionalgorithmisshowninAlgorithm2. Algorithm2semi-supervisedgraphpartition 1: ComputesimilaritymatrixWandgraphLaplacianmatrixLusingbothlabeledandunlabeleddata. 2: Solvetheeigenvectorproblem,Lu=u. Letu1,...,upbepeigenvectorscorrespondingtothepsmallesteigenvalues. 3: Fort-classlabeleddata,weestimatetseparatetwo-classclassiersbyminimizingthefollowingcostfunctions,E(ak)=sXi=1 ci)]TJ /F7 7.97 Tf 18.27 15.21 Td[(pXj=1akjuji!2 wherek=1,...,t. 4: Weusethecoefcientsa1,...,atcomputedfromstep3toclassifytheunlabeleddata.Forxi,i>s,thelabelissettoargmaxkPpj=1akjuji. 2.3.3Estimatingmeanfrompairwisegeodesicdistances Afterthegraphpartition,weneedtoconstructanatlasforeachsubgraph.Withoutlossofgenerality,wesupposefx1,,xkgarenodesofonesubgraphGaonamanifoldM.Letm2Mdenotethemeanofpointsfx1,,xkgonManddMdenotes 27

PAGE 28

theRiemanniangeodesicdistanceonM.AlthoughwedonothavetheanalyticrepresentationfordM,wecouldapproximatedM(xi,xj),1i,jkwithdQandthegraphstructureofGa.Foreachedge(xp,xq)2Ea,letthedistancedQ(xp,xq)beitsweight.Sinceageodesicistheshortestpathbetweentwopointsonthemanifold,wecoulduseDijkstra'salgorithmorFloyd-Warshallalgorithm[ 15 ]tondtheshortestpathbetweentwonodesxi,xjonthegraphanduseittoapproximatethegeodesiconthemanifoldM.Onepossiblewaytodeterminethemeanmisthefollowing.Letai=dM(xi,m),1ik.DeterminingaiisofcourseequivalenttolocatingmonM,andaicanbedeterminedasthesolutiontoanoptimizationproblemgivenby: mina21+a22++a2ks.t.ai0,i=1,kai+ajmij,ai+mijaj,aj+mijai,i,j=1,kHDCM(ai,mij), (2) wheremij=dM(xi,xj)andHDCM(ai,mij)denotetheremaininghigher-degreeconstraintsamongai,mij.ThecollectionofconstraintsinHDCMdependsonMasametricspaceandingeneral,fortwodifferentmanifolds,theircorrespondingHDCwillbedifferent.ThelinearinequalityconstraintsbetweenaiandmijareimposedbecausedMsatisesthetriangleinequality. Forourproblem,themanifoldManditsmetricdMareunknown,andconsequently,theconstraintsinHDCMcannotbeknown.However,asanapproximation,wesolvethefollowingoptimizationprobleminIRkandusethesolutiona1,,akasanapproximation 28

PAGE 29

tothesolutionoftheaboveproblem: mina21+a22++a2ks.t.ai0,i=1,kai+ajmij,ai+mijaj,aj+mijai,i,j=1,k (2) Theoptimizationproblemisclearlyconvexsincetheobjectivefunctionisstrictlyconvex(theHessianispositive-deniteeverywhere)andthedomain,whichistheintersectionofhalf-spaces,isalsoconvex.Thisparticulartypeofoptimizationproblem(withquadraticcostfunctionandlinearinequalityconstraints)canbesolvedefcientlyevenwithalargenumberofvariablesandconstraints.Andbecausetheobjectivefunctionisstrictlyconvex,thesolutionisunique[ 8 ].Furthermore,thesolutionisstablewithrespecttotheinputparametersmijinthesensethatsmallperturbationsofmijwillnotsignicantlyalterthesolutiona1,,ak[ 8 ].ForcomputingrobustL1-median[ 22 ]insteadofmean,wesimplyneedtochangethecostfunctiontoa1+a2++ak. 2.3.4Computingtheimageatlases TheresultofthepreviouslocalizationstepprovidesuswiththepairwisedistancesdM(m,[Ii]),andatthispoint,theatlashasbeenlocatedonlyintheabstractmanifoldMintermsofthesedistances.Intherealizationstep,wewillrendertheatlasimageusingthedistancedatadM(m,[Ii])andtheinputimages.Specically,giventhevertexsetf[I1],,[Ik]gofsubgraphGaandthecorrespondingnon-negativenumbersaiasestimatesonthegeodesicdistancesdM(m,[Ii]),wherei=1,...,k,werealizetheimageatlasusingaweightedgroup-wiseregistrationapproach.BecausethecloserasampleimagetotheatlasonthemanifoldM,themorecontributionitwillgiveinatlas 29

PAGE 30

realization,theweightisintheformofanexponentialfunctionexp()]TJ /F8 11.955 Tf 9.29 0 Td[(a2i=2).Sincetheweightwillbecloseto0ifaiislarge,wedetermineasmallnumberofpointsinf[I1],,[Ik]gthatareclosetomasmeasuredbydM,andmisthenapproximatedfromthesepointswithrespecttothemetricdQinQ.Inpractice,apositiveintegerKisspeciedandKpoints[Ii1],,[IiK]withshortestdistancestomarechosen.WerequirethatK1+dimensionofthemanifold. WecouldapproximatemwithaweightedmeanoftheKpoints[Iij],withweightswj,1jK:[I]=argmin[I0]2QKXj=1wjd2Q([Iij],[I0]) Theweightswjcanbeconstructedfromtheestimatedgeodesicdistancesaij:wj=bj= KXk=1bk! wherebj=exp()]TJ /F8 11.955 Tf 9.3 0 Td[(a2ij=2),1jK,forsome>0.UsingthemetricdenedbyEquation 2 ,thecorrespondingatlasIiscomputedbysolvingthevariationalproblem min~I,hj,sjKXj=1wjZjIij(hjsj(x)))]TJ /F4 11.955 Tf 10.8 2.66 Td[(~Ij2dx+Z10ZkLvjk2dxds.(2) Overall,ourmultipleatlasesconstructionapproachisoutlinedinAlgorithm3. Algorithm3Multipleatlasesconstructiononimagemanifolds 1: Constructk-NNgraphforimagedataonthemanifoldManduseAlgorithm1topartitionthegraphintotsubgraphs.Iftherearelabeleddata,applythesemi-supervisedgraphpartitiondiscussedinAlgorithm2. 2: Foreachoftsubgraphs, locatethemeanbysolvingtheoptimizationprobleminequation 2 realizetheimageatlasbytheweightedgroup-wiseregistration. 2.4Experiments Inthissection,wepresenttheexperimentalresultsofourmethodagainsttraditionalatlasconstructiontechniquesonsyntheticandrealimagedata.Thetraditionalmethods 30

PAGE 31

Figure2-3. Samplesofinputduckimages. usedascomparisonsinfollowingexperimentsarenon-rigidgroup-wiseregistrations[ 22 32 54 ]. 2.4.1ExperimentalResultsonCOIL-20 TherstexperimentisperformedonColumbiaObjectImageLibrary(COIL-20)[ 47 ]dataset.Theseimagesweretakenbyrotatingeachobjectthrough360degreeswithrespecttoaxedcamera.SinceinCOIL-20theunderlyingmanifoldof72imagesforeachobjectisone-dimensional,thisdatabaseiscommonlyusedtotestmanifoldlearningalgorithmsandisappropriatetoshowtheadvantagesofouratlas(mean)imageconstructionmethod.Inthisexperiment,thenearestneighborparameterKissetto4. 31

PAGE 32

Figure2-4. FirstRow:Karchermeanobtainedbythecompetingmethod(left)andouralgorithm(middle).Therightimagerepresentsthegroundtruth.SecondRow:Medianimageobtainedbythecompetingmethod(left)andouralgorithm(middle).Therightimagerepresentsourresultforpre-processeddatabyrotatingimageswithdifferentangles. Figure2-5. ExperimentalresultsforveobjectsfromCOIL-20.Theimagesinthetoprowarethemeansgeneratedbytraditionalmethod.Ourresultsareshowninthebottomrow. 32

PAGE 33

Wechoosethe`duck'imageasanexampletostartourdiscussion.36imagesatposesfrom5degreesto90degreesandfrom-5degreesto-90degreesarechosentobetheinputofourexperiment.SomesmaplesofinputimagesareshowninFigure 2-3 .TheKarchermeansobtainedbythecompetingnon-rigidgroup-wiseregistrationmethodandouralgorithmareshownintherstrowofFigure 2-4 .Sinceinputimagesaregeneratedbyrotatingtheobjectwithrespecttooneaxis,itisreasonabletousetheduckimageat0degreeposeasgroundtruthatlasthatisshowninFigure 2-4 .Itisobviousthatourresultisclosertothegroundtruththantheatlas(mean)imagegeneratedbythecompetingmethod.Becausethecompetingmethoddoesnotrespectthemanifoldonwhichthedataisresident,itisnotsurprisingthatitsresultsareblurryandlosesubstantialdetailsoftheinputduck.Incomparisontotherobustmedianimages(secondrowofFigure 2-4 ),ourmethodalsoyieldsabetterresult.FivemoreresultsforCOIL-20objectsarepresentedinFigure 2-5 Becauseourmethodisbasedonorbitspacewhichquotientsouttherotation,iftheinputimagesarerotatedbysomedegrees,theirprojectionsonorbitspaceareinvariantandthenewmeanwouldbeonthesameorbitoftheoldone,whichmeansthetwomeanswouldbedifferentonlybyarotation.Weconstructanewimagesetbyrotatingeachinputimageadifferentanglefrom20degreesto30degrees.AsshowninFigure 2-4 ,themeanforthisnewimagesetcomputedbyourmethodsatisesthepropertymentionedabove. 2.4.2ExperimentalResultsonSyntheticImages Inthesecondexperiment,wesyntheticallygenerate56imagesfrombunnymeshmodel.Themodeliscenteredattheoriginandwesample56locationsonaspheresurroundingthemodel.These56locationsare(approximately)uniformlysampledwithrespecttoacental(apex)location([001]),withmaximumangulardeviationfromthecentrallocationof30-degree.Inalltherenderings,weassumesimpleLambertianreectancemodelwiththedistantlightingdirectionparalleltotheviewpoint(camera 33

PAGE 34

Figure2-6. Lefttoright:StraightforwardAveraging,UnbiasedDiffeomorphicAtlasandOurAtlas. Figure2-7. Camerapositionsofinputimages.Redpointsrepresentmeanimage'sneighborscomputedbyourmethod 34

PAGE 35

Figure2-8. ThehistogramofageforallsubjectsinOASISdataset.40binsareused. location).Wechooseavalueof5fortheparameterKinthisexperiment.Figure 2-6 showstheatlasconstructedbyourmethodandthetwocomparisons(straightforwardaveragingandtheunbiaseddiffeomorphicatlas).Notethatoutputsfromthetwocompetingmethodsarenotsatisfactoryasmanydetailsandfeaturesarelostorsmoothedout. The56cameralocationsonthesphereareshownasbluepointsintheFigure 2-7 .Theredpointsrepresenttheveneighborsofmeancomputedbyouralgorithm.Fromtheresultsdepictedhere,itisevidentthattheouratlasisclosetotheexpectedmeanforthiscollectionofpointsonthesphere. 2.4.3MultipleatlasesforOASISdata Inthissection,weapplytheproposedalgorithmstoconstructatlasesfromcollectionsofMRimages.Specically,theimagedatausedinourexperimentaretheMRimagesfromthefreelyavailableOpenAccessSeriesofImagingStudies(OASIS) 35

PAGE 36

Figure2-9. Sampleimagesfromyouth,middle-agedandelderlysubjects. dataset[ 43 ].OASIScontainsT1weightedMRbrainimagesfromacross-sectionalpopulationof416subjects.EachMRIscanhasaresolutionof176208176voxels.Theagesofthesubjectsrangefrom18to96.ThehistogramofageforallOASISsubjectsisshowninFigure 2-8 Becausethestructuraldifferenceinthebrainacrossdifferentagegroupscanbesignicant[ 44 ],theOASISdatasetisheterogeneous.TheaxialviewsofthreesamplesubjectsareshowninFigure 2-9 .Therefore,it'sbettertoconstructmultipleatlasesinsteadofasinglemeantemplateimagetorepresenttheOASISdata. Weatrstconstructtwoatlasesforthecollectionof416subjectsinOASISset.Bothunsupervisedandsemi-supervisedmethodsareused.Inalltheexperiments,wesetk=10whenconstructingk-NNgraphs.Fortheunsupervisedgraphpartition,weusetwoclusters.Forthesemi-supervisedgraphpartition,werandomlychoose20subjectsyoungerthan50yearsoldand20subjectsolderthanorequalto50asthelabeledsamples.Othersubjectsareconsideredasunlabeledsamples.ThehistogramsofthepartitionresultsusingbothmethodsaredisplayedinFigure 2-10 .Table 2-1 comparesthemeansandstandarddeviationsofagesofsubjectsinbothclassescomputedbyEM-basediClusteralgorithm[ 55 ],ourunsupervisedandsemi-supervisedmethods.TheclusteringresultofourunsupervisedmethodissimilartoiCluster,whileoursemi-supervisedmethodgivesslightlysmallermeanagesofbothyoungandold 36

PAGE 37

(A)(B) Figure2-10. Twoclassespartition.A)unsupervisedmethod.B)semi-supervisedmethod. groups.Thepartiallabelsinsemi-supervisedlearningprovideusefulinformationtoclassifymoresubjectsagedbetween50and70intotheoldergroup.AsshowninFigure 2-11 ,theatlasesgeneratedbyiClusterareblurry,whiletheatlasescomputedbyourmethodsaresubstantiallysharperwithmorestructuraldetails.Thisisnotsurprisingbecauseouratlasesarecomputedasmeansonthemanifolds,insteadofthemeansintheambientspace.Thus,ouratlasesretainimportantbrainstructuresandtheygenuinelylooklikerealbrainswhencomparedwiththerealsamplesinFigure 2-9 Table2-1. Meansandstandarddeviationsofsubjectagesfortwoclassespartition iClusterunsupervisedsemi-supervised Class139.119.937.318.733.717.0Class277.89.377.09.975.39.6 Figure 2-8 alsotellsusthatthereisaclusterformiddleagedsubjectsinOASIScross-sectionalset.Thusinoursecondexperiment,wechoosethreeclustersforourunsupervisedatlasesconstructionalgorithm.Forthesemi-supervisedalgorithm,wedividetheOASISpopulationintothreegroupsbasedonthehistogramofagesshowninFigure 2-8 :young(subjectsyoungerthan40),middle-aged(between40and60)andolderadults(olderthan60)atrst.Thenwerandomlysample20subjectsfromeachgroupasthelabeleddataandconsidertheremainderunlabeled.The 37

PAGE 38

Figure2-11. TwoatlasesforOASISdata.Firstrow:atlasforyoungsubjects.Secondrow:atlasforoldersubjects.A)atlasesusingiCluster.B)atlasesusingourunsupervisedalgorithm.C)atlasesusingoursemi-supervisedalgorithm. histogramsofbothunsupervisedandsemi-supervisedpartitionresultsareshowninFigure 2-12 .Forourunsupervisedmethod,themiddleagedclusterhasasignicantoverlapinagewiththeoldercluster.Asimilarresultwasreportedin[ 55 ].Theseclusteringmethodscannotsimplyndthemiddleagemodeintheagedistribution,becausetheycanonlydiscoverthedominantstructuralmodesandthereisanobviousstructuraldifferencebetweenold(around60)andelderly(around75)subjects.Sincethesemi-supervisedmethodusesmoreinformationfrompartiallylabeleddata,itgivesabetterclassicationfordifferentmodesinagedistribution.Table 2-2 comparesthemeansandstandarddeviationsofagesofsubjectsinallthreeclassescomputedby 38

PAGE 39

(A)(B) Figure2-12. Threeclassespartition.A)unsupervisedmethod.B)semi-supervisedmethod. iCluster,ourunsupervisedandsemi-supervisedmethods,andFigure 2-13 showsthecorrespondingatlasesconstructedbythesethreealgorithms.TheclusteringresultofourunsupervisedmethodissimilartoiCluster,butourmethodbuildsmuchsharperatlases,whichareveryusefulinapplicationssuchasatlasbasedsegmentationandtensorbasedmorphometry.Themeanageofclass2returnedbyoursemi-supervisedpartitionmethodis54.6,whichismoreappropriateforthemiddleagedmodeintheagedistribution,comparedwithcorrespondingmeanagesreturnedbyiClusterandourunsupervisedmethods. Table2-2. Meansandstandarddeviationsofsubjectagesforthreeclassespartition iClusterunsupervisedsemi-supervised Class131.214.530.313.930.016.9Class268.913.667.115.254.614.9Class379.67.579.67.377.18.5 Withtheageinformationforsomesubjects,oursemi-supervisedmethodhasabilitytoconstructmorethanthreeatlasesacrossages,whichisakintothepopulationshaperegression[ 16 ].Theregressionmethodrequiresagesofallsubjectsinthedataset,whileourmethodjustneedsclasslabelsforsmallamountofimages.TheOASISdatacanbepartitionedintosevenclasseswith10yearsinterval.Werandomlychoose10 39

PAGE 40

Figure2-13. ThreeatlasesforOASISdata.Firstrow:atlasforyoungsubjects.Secondrow:atlasformiddleagedsubjects.Thirdrow:atlasforelderlysubjects.A)atlasesbyiCluster.B)atlasesbyourunsupervisedalgorithm.C)atlasesbyoursemi-supervisedalgorithm. 40

PAGE 41

Figure2-14. AtlasesacrossdifferentagesforsubjectsinOASISdataset. Figure2-15. Atlasesgeneratedbyoursemi-supervisedalgorithmforfemalesubjectsinOASISdataset.A)atlasforyoungsubjects.B)atlasformiddleagedsubjects.C)atlasforelderlysubjects. subjectsfromeachclassasthelabeleddataandconsidertheremainderunlabeled.Notethatonly16.8%ofimagesarelabeled.OuratlasesacrossdifferentagesforOASISdataareshowninFigure 2-14 .It'sobviousthattheventriclesenlargewithagingandtheexpansionisacceleratedafter60.Thisagreeswiththendingsreportedin[ 44 ]and[ 16 ]. 2.4.4Atlasesformaleandfemalesubjects TheOASISdatasetincludes160malesubjectsand256femalesubjects.Inthisexperiment,webuildthreeatlasesforbothsexgroupsusingoursemi-supervisedalgorithm.15subjectsfromeachageclassareusedaslabeleddata.OurbrainatlasesforfemaleareshowninFigure 2-15 andbrainatlasesformaleareshowninFigure 2-16 .SinceOASIScross-sectionaldatasetprovidesthesegmentationsofcerebrospinal 41

PAGE 42

Figure2-16. Atlasesgeneratedbyoursemi-supervisedalgorithmformalesubjectsinOASISdataset.A)atlasforyoungsubjects.B)atlasformiddleagedsubjects.C)atlasforelderlysubjects. (A)(B) Figure2-17. VolumeratiosofCSF,GMandWMinatlasesforA)femalesubjectsandB)malesubjects. uid(CSF),graymatter(GM)andwhitematter(WM)foreachimage,abraintissueclassicationmapcouldbegeneratedforeachatlasduringtheatlasconstructionprocess.ThevolumepercentageofCSF,GMandWMforfemaleandmaleatlasesareplottedinFigure 2-17 .Malesandfemalesappeartohavethesimilarvolumeratiochanges.ThepercentagesofCSFincreasewithageandthepercentagesofGMdecreasewithage.TherearenosignicantchangesfortheratioofWM. 42

PAGE 43

Figure2-18. AxialviewofatlasesforA)oldsubjectsandB)dementiapatients. Figure2-19. CoronalviewofatlasesforA)oldsubjectsandB)dementiapatients.C)thezoominofhippocampusregions. 2.4.5Atlasesfordementiapatientsandhealthyelderlysubjects TheOASISdatasetalsoincludesonehundredsubjectswhohavebeenclinicallydiagnosedwithverymildtomoderateAlzheimer'sdisease(AD)usingtheClinicalDementiaRating(CDR)scale.SincealltheADsubjectsinOASISdatasetareolderthan60,tomakeacomparisonofthenormalagedbrainandtheADpatient'sbrain,wechoose98healthysubjects(CDRscoreiszero)aged60orolderasthenormalgroup.ThenwebuildatlasesseparatelyforADsubjectsandhealthyeldersubjects.Theaxial 43

PAGE 44

viewsofouratlasesareshowninFigure 2-18 andthecorrespondingcoronalviewsareshowninFigure 2-19 .WecanndthattheADatlashasseverelyenlargedventriclesandshrinkageofhippocampus,comparedwiththeatlasforthenon-dementedoldsubjects. 44

PAGE 45

CHAPTER3AGEOMETRICFRAMEWORKFORTENSORFIELDSANALYSIS 3.1MotivationandRelatedWork InChapter2,weconstructafewmeantemplatestorepresentacollectionofimages.Methodsconsideringhigherorderstatisticssuchasvarianceareusuallyusedtoimprovethegroup-wiseimageanalysis.DiffusionTensorImaging(DTI)introducedbyBasseretal.[ 4 ]allowsonetomeasurethediffusionofwatermoleculesintissue,whichismodeledusingasecondordertensorwithineachvoxel.Thelocalpatternsofdiffusionareusefultorevealmicroscopicstructuralpropertiesofthemeasuredunderlyingtissue,whichmaynotbevisibleinotherimagingmodalitiessuchasanatomicalMRIandCT.BecauseDTimages,commonlyrepresentedaseldsofsymmetricpositive-denitematrices(SPDs),havespecialgeometricstructures,manytraditionalmethodsinmedicalimageanalysishavebeenextendedtothisnewtensorvalueddataform.Examplesincludediffusiontensorelddenoising[ 51 ],DTIsegmentation[ 37 75 ]andDTIregistration[ 85 87 ].TensorelddatasetsarealsocommonlyencounteredinTensor-BasedMorphometry(TBM)[ 38 39 ].Deformationeldisobtainedbywarpingeachsubjecttothetemplateusingnon-rigidimageregistrationtechniques.ThentheJacobianisusedtoconstructadeformationtensorineachvoxel.Voxel-basedstatisticalmethodsarecommonlyusedforlongitudinalstudyandcomparisonbetweenhealthyanddiseasedsubjects. Becauseoftheconestructureofthesymmetricpositive-denitematrixspace,theRiemanniangeometryoftheSPDtensorsanditsapplicationstomedicalimageanalysisproblemsthatrequirestatisticalanalysisofensemblesofSPDmatrices,havebeenthefocusofintensivestudyinthepastseveralyears,e.g.[ 23 51 56 ].In[ 51 ],Pennecetal.developedaRiemannianframeworkforcomputingstatisticsonSPDtensors.ARiemannianmetricwithsuperiorpropertiesoverEuclideanmetricwasintroduced.In[ 56 ],Schwartzmandiscussedthegeometryofpositive-denitematricesandstudied 45

PAGE 46

Figure3-1. Leftcolumn:Eightinputtensorelds.A)Mean,rstmodeandsecondmodeoftensoreldscomputedusingapixel-basedapproach:meansandmodesaredeterminedateachpixelindependentlyusingeighttensors.Notethatthetwomodesareconstantelds.B)Mean,rstmodeandsecondmodeoftensoreldscomputedusingtheproposedmethod. probabilitydistributionsdenedonSPDmatrices.PrincipalGeodesicAnalysis(PGA),asageneralizationofPrincipalComponentAnalysis(PCA)fordataonaRiemannianmanifold,wasintroducedin[ 24 ].In[ 23 ],FletcherandJoshidescribedPGA-basedmethodsforstatisticalanalysisofdiffusiontensors,inthecontextofcomputingintrinsicmeanofacollectionofSPDmatricesandcharacterizingtheirvariance. Weremarkthatalltheseearlierworkshavefocusedonthestatisticalanalysisofsymmetricpositive-denitematrices,whilethedataobtainedinmedicalimagingistensoreldinitsentirety.Thereisaneedtodevelopstatisticalanalysismethodsonacollectionoftensoreldswhereeacheldistreatedasasingleentity,similartothegroupanalysisofintensityimages.AlthoughthemethodsforSPDmatricescanbeappliedtostudythestatisticsoftensoreldsusingavoxel-basedapproachwheretensoreldsarealignedandthestatisticsaregatheredindependentlyforeveryvoxel, 46

PAGE 47

itisclearlyinsufcientandinadequateasitfailstocaptureglobalinteractionpatternsinthetensorelds.ThispointcanbebestillustratedusingasimpleexampleshowninFigure 3-1 .Here,wegeneratefourpairsof1010tensorelds.Eachtensoreldcontainstwoconstantsubeldsoccupyingthetopandbottomregionsofthedomain.Tensoreldsineachpairdifferbyareectionthatswapsthetopandbottomregions.Forthiscollectionofeighttensorelds,pixel-wisestatisticsfailstocapturetheglobalpatternsinthetensoreldsastherstandsecondmodesareallconstantelds.Thisresultisnotsurprisingasapixel-basedapproachonlyconsiderstensorsateachlocationindependently,anditcompletelyignoresthepossiblecorrelationbetweentensorsatdifferentlocations. WeproposeaRiemannianframeworkforstatisticalanalysisofasetoftensoreldsthatiscapableofcapturingtheglobalcorrelationwithintheelds.Specically,eachtensoreldcanberepresentedbyapointintheRiemannianproductspace.TheRiemannianexponentialmapofthetensoreldandthecorrespondinglogarithmicmaparedened.WealsointroduceaninvariantRiemannianmetricinthespaceoftensorelds.WeextendintrinsicmeancomputationalgorithmandPrincipalGeodesicsAnalysis(PGA)toRiemanniansymmetricproducts.ThisprovidesaprincipledmethodforlinearizingtheproblembymappingtensoreldstoaEuclideanspacethatisthetangentspaceattheirintrinsicmeantensoreld.TheglobalcorrelationsofthetensoreldarethencapturedusingPrincipleComponentAnalysis(PCA)inthetangentspace,andthemodesofvariationforthetensoreldscanbedeterminedrstinthetangentspaceandfollowedbytheexponentialmap.Moreover,dimensionalityreductionoftensorvaluedimagesisparticularlyimportantasonemajordifcultyofworkingwiththespaceoftensoreldsisitsdimension.Forexample,thespaceof100100100tensor-valuedimageshas6106dimension,anddimensionalityreductionisnecessaryformostapplicationsandanalysis.InourPGAframework,eachtensoreldisrecordedasjustalistofrealcoefcientsinsteadofamatrixvaluedimage.Thusmuchless 47

PAGE 48

Figure3-2. Thespaceoftensorelds.MandXaretwotensorelds.ThegrayplanerepresentsthetangentspaceatM.Y,Zaretwotangentvectors.ThelightbluecurveisthegeodesicbetweenMandX. spaceistakenforeachtensorimage.Wealsopresentanoveltensoreldclassicationapproachusingthegeodesicdistancestosubmanifoldsasthefeatureforclassication.TheproposedmethodisevaluatedusingtheOpenAccessSeriesofImagingStudies(OASIS)datasetandacollectionofratbraindiffusiontensorimages.Experimentalresultshavedemonstratedthesuperiorityofourmethodcomparedtothevoxel-basedapproach. 3.2StatisticalAnalysisintheSpaceofTensorFields Inthissection,wepresentthedetailsofstatisticalanalysisinthespaceoftensorelds,andwewillconsideronlyeldsofsymmetricpositive-denitematrices.LetP(n)denotethespaceofnnsymmetricpositive-denitematricesandSym(n)denotethespaceofnnsymmetricmatrices.AtensorelddenedonadomaininRKistreatedasafunctionf:!P(n).Sincebothdiffusiontensoreldsanddeformationtensoreldsusedinmedicalimagingarealmostalwaysdenedoveragrid(pixelsorvoxels),willbeacollectionofmpointsinRK,andweidentifythespaceoftensoreldson 48

PAGE 49

withtheproductP(n)m=P(n)P(n)P(n)| {z }m,whichisadifferentiablemanifoldofdimensionn(n+1)m=2.ThusatensoreldXinP(n)misrepresentedasanm-tuple(X1,X2,...,Xm),whereeachXiisasymmetricpositive-denitematrix,thevalueofXatonepointin. 3.2.1GeometryofTensorFields ThespaceP(n)isasymmetricRiemannianmanifold[ 23 ]withGL(n)asthesymmetrygroup.ThiscanbegeneralizeddirectlytoproductspacesP(n)musingproductRiemannianstructure,andinparticular,theRiemanniangeodesicdistances,logandexponentialmapshaveclosed-formexpressions.Figure 3-2 illustratesthespaceofSPDmatrixelds.Specically,thegroupGL(n)mactstransitivelyonP(n)mwiththeaction:GL(n)mP(n)m!P(n)mgivenbyG(X)=(G1X1GT1,...,GmXmGTm),whereeachGi2GL(n)isanninvertiblematrixandXiisannpositive-denitematrix.Intuitively,thegroupactionofGL(n)mcanbeinterpretedasachangeofcoordinatesintheproductspaceoftensoreldsP(n)m.ThetangentspaceofP(n)matanypointcanbeidentiedwithSym(n)msincethetangentspaceofaproductmanifoldistheproductoftangentspaces.LetY,Z2TMP(n)mbetwotangentvectorsatM2P(n)m.TheproductRiemannianmetricgivestheinnerproductbetweenthetwovectorsas hY,ZiM=mXi=1tr(YiM)]TJ /F5 7.97 Tf 6.59 0 Td[(1iZiM)]TJ /F5 7.97 Tf 6.59 0 Td[(1i).(3) ThisscaledinnerproductisinvariantunderthegroupactionofGL(n)m.Specically,giventhegroupactiononP(n)m,thedifferentialmapofisdG(Y)=(G1Y1GT1,...,GmYmGTm). ForanyG2GL(n)m, hdG(Y),dG(Z)iG(M)=hY,ZiM.(3) 49

PAGE 50

Usingthismetric,thegeodesicpassingthroughMinthedirectionofYisuniquelygivenby M(t;Y)=)]TJ /F8 11.955 Tf 5.48 -9.68 Td[(G1exp(G)]TJ /F5 7.97 Tf 6.59 0 Td[(11Y1G)]TJ /F7 7.97 Tf 6.58 0 Td[(T1t)GT1,...,Gmexp(G)]TJ /F5 7.97 Tf 6.59 0 Td[(1mYmG)]TJ /F7 7.97 Tf 6.59 0 Td[(Tmt)GTm(3) whereGi2GL(n)suchthatM=)]TJ /F8 11.955 Tf 5.48 -9.68 Td[(G1GT1,...,GmGTm.TheRiemannianexponentialmapatMwhichmapsthetangentvectorYtoapointinP(n)misdenedasthevalueofthegeodesicM(t;Y)att=1.Inourcase,theclosed-formrepresentationoftheRiemannianexponentialmapis ExpM(Y)=)]TJ /F8 11.955 Tf 5.48 -9.68 Td[(G1exp(G)]TJ /F5 7.97 Tf 6.58 0 Td[(11Y1G)]TJ /F7 7.97 Tf 6.59 0 Td[(T1)GT1,...,Gmexp(G)]TJ /F5 7.97 Tf 6.59 0 Td[(1mYmG)]TJ /F7 7.97 Tf 6.59 0 Td[(Tm)GTm(3) whereGi2GL(n)suchthatM=)]TJ /F8 11.955 Tf 5.48 -9.68 Td[(G1GT1,...,GmGTm. Itsinversemap,Riemmannianlogarithmicmapcanbeeasilyobtained.GivenX2P(n)m,thelogmapatMis LogM(X)=)]TJ /F8 11.955 Tf 5.48 -9.68 Td[(G1log(G)]TJ /F5 7.97 Tf 6.59 0 Td[(11X1G)]TJ /F7 7.97 Tf 6.58 0 Td[(T1)GT1,...,Gmlog(G)]TJ /F5 7.97 Tf 6.58 0 Td[(1mXmG)]TJ /F7 7.97 Tf 6.58 0 Td[(Tm)GTm.(3) UsingthisdenitionoflogmapinP(n)m,thegeodesicdistancebetweentwotensoreldsMandXiscomputedby d(M,X)=kLogM(X)k=vuut mXi=1tr)]TJ /F4 11.955 Tf 5.48 -9.68 Td[(log2(G)]TJ /F5 7.97 Tf 6.59 0 Td[(1iXiG)]TJ /F7 7.97 Tf 6.59 0 Td[(Ti).(3) 3.2.2StatisticsontheSpaceofTensorFields Usingtheformulaaboveforthegeodesicdistance,wedenetheintrinsicmeanofacollectionofNtensoreldsX1,...,XNasthetensoreldthatminimizesthesumofsquaredgeodesicdistances: M=argminM2P(n)m1 NNXi=1d(M,Xi)2.(3) ThesumofsquaresontheRHSabovecanbere-writtenasasumoverallpointsin.ThisimpliesthatthevalueofM(p)ofMatonepointp2istheusualKarchermean 50

PAGE 51

inP(n)ofX1(p),,XN(p).Inparticular,becausethespaceofSPDmatricesP(n)iscompleteandhasnon-positivesectionalcurvature,theKarchermeanisuniqueonP(n)[ 23 ].ThisshowsthatMwillbeuniqueaswell,anditcanbecomputedusingthefollowinggradientdescentalgorithm. Algorithm4Intrinsicmeantensoreld 1: SettheinitialmeantensoreldM0asoneofNinputtensoreldsXi.Settheinitialstepset=1anditerationk=0. 2: ComputetheaverageonthetangentspaceatMk:Yk=1 NNXi=1LogMkXi. 3: Ifk>0andkYkk>kYk)]TJ /F5 7.97 Tf 6.59 0 Td[(1k,adjuststepsize==2,setYk=Yk)]TJ /F5 7.97 Tf 6.59 0 Td[(1,Mk=Mk)]TJ /F5 7.97 Tf 6.59 0 Td[(1. 4: Projectbacktothespaceoftensorelds:Mk+1=ExpMk(Yk). 5: Stoppingcondition:ifkYkk<,outputMk+1astheintrinsicmeantensoreld,otherwisesetk=k+1andgotostep2. AfterobtainingtheintrinsicmeanMoftheinputtensoreldsX1,...,XN,wewilldeterminethemodesofvariationusingPGA.Specically,weusethelogmaptomapallthetensoreldstothetangentspaceatM,xi=logM(Xi).ThisisaEuclideanspaceinwhichwecananalyzethedatapointsx1,,xNusingprincipalcomponentanalysis.WedenetheprincipalvectorsV1,,VkinTMP(n)maccordingtothefollowingequations: V1=argmaxkVk=1NXi=1hV,LogM(Xi)i2M,Vk=argmaxkVk=1NXi=1k)]TJ /F5 7.97 Tf 6.58 0 Td[(1Xj=1hVj,LogM(Xi)i2M+hV,LogM(Xi)i2M. 51

PAGE 52

TheorthonormalvectorsVispannedaK-dimensionalsubspaceSKthatbestapproximatesx1,,xNintheleast-squaressense,andtheycanbecomputedusingPCA.ByexponentiatingvectorsinSK,ExpMPdk=1kVk,wherektellsthevariationofkthmode,weobtainthegeodesicsubmanifoldSKP(n)mthatcanserveasagoodapproximationoftheinputtensorelds. TherearetwoimportantdetailsthatdifferfromtheusualapplicationofPGA[ 23 ].First,exceptattheidentity,theinnerproductdenedinEquation 3 doesnotcorrespondtothestandardEuclideaninnerproduct,whichisrequiredforthefamiliarPCAalgorithm.Therefore,wersttransformthedatatothetangentspaceattheidentity,whichisaccomplishedviathefollowingtransform,X=(X1,Xm)2TMP(n)m dG)]TJ /F9 5.978 Tf 5.76 0 Td[(1(X)=)]TJ /F8 11.955 Tf 5.48 -9.69 Td[(G)]TJ /F5 7.97 Tf 6.59 0 Td[(11X1G)]TJ /F7 7.97 Tf 6.58 0 Td[(T1,...,G)]TJ /F5 7.97 Tf 6.59 0 Td[(1mXmG)]TJ /F7 7.97 Tf 6.59 0 Td[(Tm,(3) whereG=(G1,...,Gm)issuchthatM=)]TJ /F8 11.955 Tf 5.48 -9.68 Td[(G1GT1,...,GmGTm.OncethedatahavebeenmappedtoTIP(n)m,wecanapplytheusualPCAalgorithmtoobtainprincipalvectorsUi,i=1,,K.TheyarethentransformedbacktoTMP(n)musingVi=dG(Ui).Duetothehigh-dimensionalityofP(n)m,weusetheGrammatrixinsteadoftheusualcovariancematrixwhencomputingtheprincipalvectorsinthetangentspace,whichistheapproachusedinmanycomputervisionapplicationssuchastheEigenfaces[ 68 ].ThecompletealgorithmissummarizedinAlgorithm5. Algorithm5PGAforTensorFields 1: InputNtensoreldsX1,...,XN2P(n)m. 2: ComputeintrinsicmeanMofinputtensoreldsusingthemethodlistedAlgorithm4. 3: ComputeYi=LogM(Xi)fori=1,...,N. 4: TranslateYitothetangentspaceofidentityI. 5: PerformPCAinTIP(n)mandgeteigenvectorsUi. 6: TranslateUibacktogetViinthetangentspaceofM. 52

PAGE 53

Figure3-3. Thelogmapontheunitsphere.AandBaretwopointsontheunitsphere.A0andB0aretheirlogprojectionstothetangentplaneatM. 3.3TensorFieldsClassication Wecanformulateatensoreldclassicationalgorithmusingtheprincipaldirectionsandgeodesicsubmanifolds.OnecommonmethodforsolvingclassicationproblemsonRiemannianmanifoldistomapinputdatatothetangentspaceanddotheclassicationinthetangentspace[ 79 ].However,thisapproachdoesnotrespectthegeometryofthemanifoldasthegeodesicdistancebetweentwopointsonthemanifoldareusuallynotthesameorevencommensuratewiththedistancebetweentheirimagesinthetangentspace.OneexampleisshowninFigure 3-3 .LetA=(1,0,0),B=(0,1,0),M=(0,0,1)bethreepointsontheunitsphere.A0andB0arethelogprojectionsofAandBtothetangentplaneatM.ThespatialrelationofthreepointsonthesphereisnotpreservedonthetangentplaneatM,sincethedistancebetweenA0andB0(p 2 2)isnotthesameasthegeodesicdistancebetweenAandBontheunitsphere( 2).Amore 53

PAGE 54

principledapproachistousethedistancestogeodesicsubmanifoldsasthefeatureforclassication. Assumeabinaryclassicationproblem,andthetrainingtensoreldsarelabelledasoneofthetwoclasses.Foreachlabelk(k=1,2),wecomputealow-dimensionalgeodesicsubmanifoldSkusingtrainingtensoreldswithlabelk.ForatesttensoreldX,wecandetermineitsclassbycomparingthegeodesicdistancesdk=minY2Skd(X,Y).Atensoreldisclassiedasbelongingtoclasskifdkissmallerthantheothergeodesicdistance.ThekeystepinthisalgorithmistondtheminimizerinSkthatgivesthegeodesicdistancedk.SinceanypointinthegeodesicsubmanifoldSkcanbewrittenasExpMPdi=1iVi,whereMisthemeanand1,...,darerealcoefcients,dkcanbesolvedviathefollowingoptimizationprobleminRd, min1,...,dd X,ExpM dXi=1iVi!!2.(3) Unfortunately,minimizingEquation 3 canbetime-consumingforlargetensorelds.Therefore,weapproximatedkbythegeodesicdistanced(X,Z)betweenXandZdenedby Z=ExpM dXi=1VihVi,LogM(X)iM!.(3) Thatis,weobtainZbyrstmapXtothetangentspaceatMusingLogmapandprojectitontotheprincipalsubspaceSd.TheprojectionisthenmappeddowntothemanifoldusingtheexponentialmaptogetZ.Inourexperimentsdiscussedinthenextsection,thedifferencesbetweentheapproximateddistanceandtheonecomputedbytheoptimizationarelessthan1%andhavenomajorinuenceontheclassicationresults.WesummarizethetensoreldsclassicationalgorithminAlgorithm6.Ouralgorithmcanbeeasilyextendedtomulti-classcasesbycomputinggeodesicsubmanifoldsforeachclassinthetrainingstepandassigningthetesttensoreldtotheclasswiththesmallestgeodesicdistance. 54

PAGE 55

Algorithm6TensorFieldsClassication 1: TrainingcomputegeodesicsubmanifoldsS1andS2byusingPGAontrainingtensoreldsofdifferentclassesseparately. 2: Testingforeachtesttensoreld,computetheirgeodesicdistancesd1andd2tosubmanifoldsS1andS2respectively.Ifd1
PAGE 56

Figure3-5. Comparisonswithvoxel-basedmethod.(a)and(c)aretherstmode(1)andsecondmode(2)computedusingAlgorithmFive.(b)and(d)aretherstmodeandthesecondmodecomputedusingvoxel-wisePGA.FAisusedasthebackgroundinthedisplay. voxels.Becausethestructuraldifferenceinthebrainacrossdifferentagegroupscanbesignicant[ 44 ],theOASISdatasetisheterogeneous.Wedivided416subjectsthreegroups:youngsubjects(40oryounger),middle-agedsubjects(between40and60yearsofage)andoldsubjects(60orolder).ThereisasubsetofoldsubjectswhohavebeenclinicallydiagnosedwithverymildtomoderateAlzheimersdisease(AD)usingtheClinicalDementiaRating(CDR)scale.WecomputetheatlasforalltheMRimagesintheOASISdatasetusingagroup-wisenonrigidregistration[ 32 ],andthisalsogivesthethedeformationeldfromeachimagetotheatlas.Foreachvoxel,wecomputetheJacobianmatrixJofthedeformationeldandbuildthedeformationstraintensorS=(JTJ)1=2.Thisgivesadeformationtensoreldforeverysubject. Intherstexperiment,wecharacterizethemeanandvariationinthetensoreldswithinanagegroupbycomputingthemodesofvariationusingAlgorithm4andAlgorithm5.Thedimensionofthegeodesicsubmanifoldissetat20afterexaminingtheeigenvaluedistribution.Figure 3-4 displaystheintrinsicmeantensoreldandthevariationsalongthersttwoprincipaldirectionsfortheoldgroup.Thecomparison 56

PAGE 57

withthemeanandmodescomputedusingthevoxel-basedapproachinshowninFigure 3-5 .Wecanclearlyseethatthemodescomputedusingvoxel-basedmethodarefragmentaryandtheydonotreecttheglobalstructuresofthetensorelds.Thisisnotsurprisingbecausevoxel-basedmethoddoesnotconsidercorrelationsbetweendifferentvoxels. Forthesecondexperiment,wetestourtensor-eldclassicationalgorithm.Werandomlydividethebrainimagesforeachagegroupintofoursubsets.Imagesfromoneofthesubsetsarethetestimages,whileotherthreesubsetsgivethetrainingimages.ThetrainingimagesareusedtocomputethegeodesicsubmanifoldsSd,andweclassifythetestimagesforeverypairofagegroupsusingAlgorithm6.Weuseafour-foldcross-validationintheexperimentstofullyevaluatethealgorithmonOASISdataset.WecomparedtheperformanceofouralgorithmwiththeLog-Euclideannearestneighbormethodthatmapseachtensoreldtothetangentspaceofthemean.Low-dimensionalfeaturevectorsaregeneratedusingPCAprojection,andtheclassicationisdoneusingnearestneighborofthesefeaturevectors.Wehavealsotestedtheproposedalgorithmonclassicationforhealthyanddiseased(Alzheimer)brainimages.AllresultsaretabulatedinTable 3-1 .Theexperimentalresultsindicate Table3-1. TensorFieldsClassicationonOASIS Ovs.YOvs.MMvs.YADvs.Control NearestNeighbor92.43%87.74%78.42%84.29%SubmanifoldProjection96.43%90.23%84.32%88.57% thatthedeformationstraintensoreldsdocapturethesubtlestructuralchangesinthebrainimagesacrossdifferentagegroups(andalsofordiseasedsamples),andthegoodclassicationresultsshowtheimportanceandvalueofdoingstatisticalanalysisonthespaceoftensorelds. 3.4.2Diffusiontensorelds Inthethirdexperiment,weuseasetofratbraindiffusiontensorimages.Inthepreprocessingstep,allimagesarealignedusingsimilaritytransforms,andineach 57

PAGE 58

Figure3-6. Fivesamplesfrominputratbraindiffusiontensorimages.S0issetasthebackgroundinthedisplay.Wehaveenhancedtheanisotropyforbettervisualization. (A)(B)(C)(D)(E) Figure3-7. Statisticalanalysisfordiffusiontensorimages.A)meantensoreld.B-E)tensoreldvariationalongtherstfourprincipaldirections. image,aregionofinterestofsize46277ismanuallyselected.Figure 4.4.2 displaysvesampleslicesfromtheinput3Ddiffusiontensorimages.Weenhancethetensoranisotropyinthevisualizationtobetterdisplaytheratbrainstructures.Thediffusiontensoreldsfrom50subjectsarechosentoconstructtheprincipalgeodesicsubmanifoldapproximationoftheimageset.TheintrinsicmeantensoreldandthevariationsalongtherstfourprincipaldirectionsareshowninFigure 3-7 .Ourprincipalgeodesicsubmanifoldrepresentationclearlycapturesimportantanatomicalstructuressuchaswhitematter,putamenandnucleusaccumbens.Wealsotestthereconstruction 58

PAGE 59

(A)(B)(C)(D)(E) Figure3-8. Diffusiontensorimagereconstruction.A)onesliceoforiginaldiffusiontensorimage.B-E)reconstructedtensorimagesinprincipalsubspaceswithdimensions5,15,25,35respectively. Figure3-9. Reconstructionerrorsbetweenthetesttensoreldanditsprojectionsonthegeodesicsubmanifoldsofdifferentdimensions. 59

PAGE 60

performanceofourmodelbycomputingtheprojectionsofanewtensoreldtothegeodesicsubmanifoldsgeneratedbydifferentnumbersofprincipalcomponents.Figure 3-8 displaystheoriginalinputdiffusiontensorimageI0anditsprojectionsItothelowdimensionalgeodesicsubmanifoldscomputedusingAlgorithm5.Normalizedgeodesicdistance1 md(I0,I)isusedtoquantitativelymeasurethereconstructionerror,whichisplottedinFigure 3-9 .Asthenumberoftheprincipalcomponentsincreases,theprojectionsinthegeodesicsubmanifoldmoreaccuratelyapproximatedtheoriginaldiffusiontensorimage. 60

PAGE 61

CHAPTER4NONNEGATIVEFACTORIZATIONOFDIFFUSIONTENSORIMAGES 4.1Motivation Inthischapter,weintroducethenovelnotionofnonnegativefactorizationoftensoreldsanditsapplicationtothesegmentationofdiffusiontensorimages(DTI).Tensorimages(ortensorelds)aboundinmedicalimagingapplications,withwell-knownexamplessuchasDTIs,conductivitytensorimages(CTI)andelasticitytensorsinElastography.Thereisanincreasingdemandforaprincipledandversatilemethodthatcanautomaticallyextractimportantfeaturesfromsingleaswellasmultipletensorimages,andourworktakesastepforwardinthisdirectionbyinvestigatingthegeneralfactorizationproblemfortensorelds.Specically,weformulatethefactorizationasanoptimizationproblem,andasthemaintechnicalcontribution,wepresentaniterativealgorithmthatcanefcientlyandreliablyoptimizetheobjectivefunctionandcomputegood-qualityfactorizations. MatrixfactorizationintheformoffactoringadatamatrixVintoaproductofthebasismatrixWandthecoefcientmatrixHappearsfrequentlyincomputervisionandimageprocessing.Examplesincludetheprincipalcomponentanalysis(PCA)andtheTomasi-Kanadefactorizationfrommulti-view3Dstructurerecovery.Algebraically,thefactorizationattemptstodiscoverthelinearsubspacespannedbythecolumnsofWthatcanbeusedtoapproximatethedatagivenasthecolumnsofV.Insteadofthefullsubspace,thelinearmodelcanbefurtherconstrainedbyrequiringthatthedatavectorsbelongto(orcanbeapproximatedby)theconegeneratedbythecolumnsofW.Sinceapointinaconecanbewrittenasanonnegativelinearcombinationofthecone'sgenerators,thisgivesrisetothenonnegativematrixfactorization(NMF)rstintroducedin[ 36 ].Specically,usingimagesasdata,thematricesappearedinanNMFareallrequiredtohavenonnegativecomponents.Whiletheunconstrainedmatrixfactorizationadmitsacomparativelystraightforwardsolutionusinglinearalgebra, 61

PAGE 62

nonnegativematrixfactorizationisgenerallymoredifculttosolveasitsobjectivefunctionisnon-convexandthereexistsnoknownlinearalgebraicmethodsthatcancomputeitssolutioninclosedform.Nevertheless,theoptimizationisnotparticularlyinvolvedasthenonnegativeconstraintcanbeeasilyenforced,andNMFisalmostalwayssolvedusinganalternatingsequenceofnonnegativelinearleastsquareswithvariablesinWandHseparately.Moreimportantly,thenonnegativityconstrainthasbeenarguedintheoriginalpaper[ 36 ]tobemorecapable,comparedwiththeunconstrainedfactorization,ofcapturingfeaturesthatarecommonamongthedata.Inparticular,thecolumnsofthebasismatrixWcanbeconsideredasthediscoveredcommonpartsamongtheinputdata,andthecoefcientmatrixHprovidestheweightsforreconstructingeachinputdatausingtheseparts.ThepowerandversatilityofNMFhavebeendemonstratedtovariousdegreesbyawiderangeofsuccessfulapplicationsthatincludedataclustering[ 17 ],partsdiscovery[ 41 ]andMRimageanalysis[ 33 ]. Theproposednonnegativefactorizationoftensor-valuedimagesisanextensionofthenonnegativematrixfactorizationmethodabove.Inourcontext,componentsinthematricesVandWarerepresentedassymmetricpositivesemi-denite(PSD)tensorsofranktwoandHisthecoefcientmatrixwithnonnegativerealcomponents.Inparticular,thenonnegativityconstraintinourmoregeneralsettinghasnowacquiredtwodifferentrealizations:thenonnegativeconstraintonthecomponentsofHasintheoriginalNMFandthegeneralizednonnegativeconstraintonthecomponentsofWthattheybelongtothespaceofPSD(n)ofnnsymmetricpositivesemi-denitematrices.Geometrically,thiscorrespondstoreplacingthenonnegativerealline(intensityvalues)withthespacePSD(n)(diffusiontensors),anditistheirrespectiveconestructuresthatpermitustoformulatenonnegativefactorizationusingthesespaces.Whileourextensioniseasytounderstandconceptually,theresultingoptimizationprobleminourcaseisconsiderablymoredifculttosolvebecauseofthegeneralizednonnegativeconstraintonWandthedimensionofPSD(3).Theformerrequiresgeneralizedmatrix 62

PAGE 63

inequalitiestoenforcetheconstraint[ 8 ]andthelatterintroducesalargenumberofvariablesintheobjectivefunction.Therefore,amajorportionofthischapterisdevotedtoanoptimizationalgorithmthatcanreliablyandefcientlycomputethefactorization. Havingovercomethiscomputationalhurdle,wewillnextshowthattheproposedtensorimagesfactorizationmethodcanbesuccessfullyappliedtosegmentationproblemsforsingleandmultipleDTIimages.Forasingleimage,thedatamatrixVisan1narrayofPSDtensors,andadirectclusteringonthecoefcientmatrixHgivesthesegmentation.Formultipletensorimages,thebasismatrixWgivesasbeforeapart-decompositionofthecollection,andthecommonpartsgivenasthecolumnsofWcanusuallyberealizedastensoreldswithlocalsupport.Inthecurrentmedicalimagingliterature,algorithmsthatsegmentDTIimagescanbebroadlyclassiedintotwocategories,thelevel-setbasedmethods(e.g.,[ 37 75 ])andthemethodsbasedoncombinatorialoptimizationsuchasgraphcuts[ 76 ].Intermsofitsunderlyingmotivationandnumerics,ourfactorization-basedmethodoffersacompletelydifferentapproachtothesegmentationproblem.Wehavevalidatedtheproposedmethodusingbothsyntheticandrealtensorimages.PreliminarysegmentationresultsfromsingleandmultipleDTIimageshaveindicatedthattheproposedfactorization-basedapproachisaviableandcompetitivealternativemethodforsegmentingtensorimages. 4.2Preliminaries Givenannmnon-negativematrixV,eachofwhosecolumnvirepresentsaninputimage,non-negativematrixfactorization(NMF)[ 36 ]attemptstofactoritintotwomatriceswithnon-negativecomponents,VWH,whereWdenotesthenrbasismatrixandHdenotesthermcoefcientmatrix.Accordingtothegeometricinterpretationelucidatedin[ 18 ],NMFdeterminesaconeW=fxjx=Prj=1hjwj,hj0gthatapproximatestheinputimageswiththebasiselements(orgeneratorsofW)wjgivenbythej-thcolumnofW. 63

PAGE 64

Fortensor-valuedimages,thereisaddsymmetric,positivesemi-denitematrixassociatedtoeachpixel(orvoxel).d=3fordiffusiontensorimages.ThemathematicsfornonnegativefactorizationoftensorimagesisstraightforwardasitonlyrequiresaconestructurethatisprovidedbyPSD(d),thespaceofddsymmetricpositivesemi-denitematrices.AcollectionofmtensorimagesofsizencanbearrangedintoablockmatrixV=0BBBB@V11V1m.........Vn1Vnm1CCCCA,whereVki2PSD(d),k=1,...,nandi=1,...,m.Eachofthemcolumnsrepresentsonetensorimage.Asbefore,nonnegativefactorizationfortensorimagesattemptstofactorVintoaproductofthebasismatrixWthecoefcientmatrixHwhoseelementsarenon-negativerealnumbers: VWH=0BBBB@W11W1r.........Wn1Wnr1CCCCA0BBBB@h11h1m.........hr1hrm1CCCCA,(4) wheretheblocksWijinWarematricesinPSD(d),andtheblockwiseproductisdenedas WH=0BBBB@Prj=1W1jhj1Prj=1W1jhjm.........Prj=1Wnjhj1Prj=1Wnjhjm1CCCCA.(4)rintheaboveequationisthenumberofbasiselements(generators)usedinthefactorization,anditisclearthatournonnegativefactorizationreducestotheusualNMFwhend=1.Weremarkthatthenonnegativefactorizationcanbeconsideredasagenerativemodelforthecollectionofinputtensorimagesaseachinputtensorimageisapproximatedbyanon-negativelinearcombinationofrcolumnsofW,eachofwhichcanbeconsideredasatensorimage.TodetermineW,HfromthedatamatrixV,we 64

PAGE 65

formulateaconstrainedoptimizationproblemthatminimizesthecostfunction E(W,H)=1 2mXi=1nXk=1kVki)]TJ /F7 7.97 Tf 18.67 14.95 Td[(rXj=1Wkjhjik2F(4) withtheconstraintsWkj0andhji0,i=1,...,m,j=1,...,randk=1,...,n.denotesmatrixinequalityandkkFdenotestheFrobeniusnorm. 4.3Algorithmandimplementationdetails Inthissection,wepresentanefcientalgorithmthatsolvestheconstrainedoptimizationdenedabove.WhiletheobjectivefunctionE(W,H)isnotconvex,itisconvexwithrespecttothetwoblockvariablesWandH.Acommonapproachtosolvethistypeofconstrainedoptimizationproblemistheblockcoordinatedescentmethod[ 48 ],andinparticular,wecanalternativelyxoneblockvariableandimprovetheotherasshowninAlgorithm7.Grippoetal[ 27 ]haveshownthateverylimitpointofthe Algorithm7AlternatingNon-negativeFactorization InitializeH10. Fort=1,2,... Wt+1=argminWE(W,Ht),s.t.Wkj0,8k,j. Ht+1=argminHE(Wt+1,H),s.t.H0. sequencefWt,HtggeneratedbyAlgorithm7isastationarypointoftheoptimizationproblemdenedin(3).Sincebothsub-problemsareconvex,ouralgorithmcaneasilybeshowntobeprovablyconvergent.Thedetailsforsolvingthetwosub-problemsefcientlywillbediscussedbelow. 4.3.1OptimizationwithrespecttothebasismatrixW WhenHisxed,theoptimizationproblem(3)reducestoaquadraticsemi-deniteprogrammingproblemwithalargenumberofpositivesemi-denitematricesastheconstrainedvariables.Thisisachallengingoptimizationproblemwithoutreadilyavailablesolversasmostavailablesemi-deniteprogrammingsolverssuchasSeDuMi 65

PAGE 66

[ 62 ]andSDPT3[ 69 ]requirethelinearobjectivefunctions.Tohetal.[ 67 ]proposedaninexactprimal-dualalgorithmtosolveaspecialclassofconvexquadraticsemi-deniteprogrammingproblem.However,theiralgorithmonlydealswithasinglepositivesemi-denitematrixandcannotbedirectlyappliedtoouroptimizationproblem.Instead,wewillexploitthespecialfeatureinourproblemthatd=3isarelativelysmallnumberanddesignaspecicalgorithmbasedonprimal-dualpath-followinginterior-pointmethodtosolvethissubproblem. Theprimalproblemisgivenby minW1 2mXi=1nXk=1kVki)]TJ /F7 7.97 Tf 18.67 14.94 Td[(rXj=1Wkjhjik2Fs.t.Wkj0,k=1,...,n,j=1,...,r.(4) WeintroducetheddsymmetricmatricesZkjassociatedwiththematrixinequalitiesWkj0,wherek=1,...,n,j=1,...,r.TheLagrangianofproblem(4)isthen L(W,Z)=1 2mXi=1nXk=1kVki)]TJ /F7 7.97 Tf 18.67 14.94 Td[(rXj=1Wkjhjik2F)]TJ /F7 7.97 Tf 18.3 14.94 Td[(nXk=1rXj=1Tr(ZkjWkj).(4) IfwetakethederivativeofL(W,Z)withrespecttoWkjandsetittozero,wehave Zkj=)]TJ /F7 7.97 Tf 16.41 14.94 Td[(mXi=1(Vki)]TJ /F7 7.97 Tf 18.66 14.94 Td[(rXl=1Wklhli)hji.(4) SubstitutingthisexpressionbackintotheLagrangiangivesthedualproblem maxZmXi=1nXk=1 )]TJ /F4 11.955 Tf 10.5 8.08 Td[(1 2kVki)]TJ /F7 7.97 Tf 18.67 14.95 Td[(rXj=1Wkjhjik2F+!s.t.mXi=1 Vkihji)]TJ /F7 7.97 Tf 18.66 14.95 Td[(rXl=1Wklhjihli!+Zkj=0Zkj0,k=1,...,n,j=1,...,r.(4) 66

PAGE 67

Forprimal-dualinteriorpointmethod,weusetheperturbedKarush-Kuhn-Tucker(KKT)conditions: Zkj+mXi=1(Vki)]TJ /F7 7.97 Tf 18.67 14.94 Td[(rXl=1Wklhli)hji=0WkjZkj=kjIWkj0,Zkj0k=1,...,n,j=1,...,r(4) wherekjarepositiveparameters.Giventhecurrentiterate(W,Z),thesearchdirection(W,Z)ateachinterior-pointiterationisthesolutionofthefollowingNewtonsystem rXl=1Wkl(mXi=1hlihji))]TJ /F4 11.955 Tf 11.95 0 Td[(Zkj=Rdkj:=Zkj+mXi=1(Vki)]TJ /F7 7.97 Tf 18.66 14.95 Td[(rXl=1Wklhli)hji(4) HPkj(WkjZkj+WkjZkj)=kjI)]TJ /F8 11.955 Tf 11.96 0 Td[(HPkj(WkjZkj+WkjZkj)(4) wherekj==dand2(0,1)isthecenteringparameter.ThesymmetrizationschemeHPdenedasHP(M)=1 2[PMP)]TJ /F5 7.97 Tf 6.59 0 Td[(1+P)]TJ /F7 7.97 Tf 6.59 0 Td[(TMTPT]isrequiredheretogeneratesymmetricWkjandZkj.Inspiredby[ 66 ],wechoosePkj=T)]TJ /F5 7.97 Tf 6.58 0 Td[(1=2kj,whereTkjistheNesterov-ToddscalingmatrixsatisfyingTkjZkjTkj=Wkj.Thelinearizationofequation(10)gives Ekj(Wkj)+Fkj(Zkj)=Rckj:=kjI)]TJ /F8 11.955 Tf 11.95 0 Td[(HPkj(WkjZkj)(4) whereEkjandFkjarelinearoperators.ByeliminatingZkjinequations(9)and(11),weobtain rXl=1(mXi=1hlihji)Wkl+F)]TJ /F5 7.97 Tf 6.59 0 Td[(1kjEkj(Wkj)=Rdkj+F)]TJ /F5 7.97 Tf 6.59 0 Td[(1kjRckj.(4) Therefore,wejustneedtosolvenlinearsystemstogetW.Eachlinearsystemincludesd(d+1)r 2equations.Becausebothdandraresmallforourproblem,theselinearsystemscanbesolvedefciently,andZcanbecomputedeasilyusingequation(9).ThedetailstepsofthealgorithmisshowninAlgorithm8. 67

PAGE 68

4.3.2OptimizationwithrespecttocoefcientmatrixH WithaxedW,(3)becomesaconvexoptimizationproblem.Since @E @hji=)]TJ /F7 7.97 Tf 17.64 14.94 Td[(nXk=1Tr(VkiWkj)+nXk=1rXl=1Tr(WkjWkl)hli,(4) setting@E hji=0gives rXl=1(nXk=1)hli=nXk=1.(4) Thuswejustneedtosolveanon-negativeleastsquaresproblemArrHrm=Brmtogethji,whereAjl=Pnk=1andBji=Pnk=1,foralli,j,l.Inourimplementation,weusethefastactivesetmethodproposedbyVanBenthemandKeenan[ 71 ]tosolvethislarge-scalenonnegativeleastsquaresproblem. OncewehaveobtainedthebasismatrixW,wecaneasilycomputetheprojectionofanewdiffusiontensorimageX=(X1,...,Xn)Tbysolvingthefollowingoptimizationproblem miny0nXk=1kXk)]TJ /F7 7.97 Tf 18.67 14.94 Td[(rXj=1Wkjyjk2F(4) whereeachXkisadiffusiontensorandy=(y1,...,yr)Tisthenonnegativecoefcientvector.Thisproblem,similartothesubproblemwithrespecttoHabove,canalsobeefcientlysolvedusingthesamefastactivesetmethod. 4.3.3Diffusiontensorimagesegmentation BecauseofitsrelationstoK-means[ 17 ]andprobabilisticlatentsemanticanalysis[ 25 ],nonnegativematrixfactorizationhasbeenwidelyusedindataclustering,e.g.,[ 83 ].Inthissection,weformulatethediffusiontensorimagesegmentationproblemasaspecialcaseofourmoregeneralnonnegativefactorizationproblemwithspatialconstraints. Specically,givenadiffusiontensorimageofsizem,wearrangethemtensorsinarowtoformthedatamatrixV=(V1,...,Vm).Theoriginalnonnegativefactorization 68

PAGE 69

Algorithm8Quadraticsemi-deniteprogrammingfornonnegativefactorization 1: InitializationWkj0,Zkj0,8k,jand=0.9. 2: ConvergencetestStoptheiterationiftheaccuracymeasureissufcientlysmall.=Pnk=1Prj=1 1+jpobjj+jdobjj wherepobjanddobjarevaluesofprimalanddualobjectivefunctions. 3: PredictorstepComputethepredictorsearchdirection(W,Z)bychoosing=0. 4: Predictorstep-lengthComputep=min(1,).isthemaximumsteplengththatcanbetakensothatfork=1,...,nandj=1,...,r,Wkj+WkjandZkj+Zkjremainpositivesemidenite. 5: CenteringruleSet=Pk,j Pk,j. 6: CorrectorstepComputethesearchdirection(W,Z)usingRckj=kjI)]TJ /F8 11.955 Tf 11.96 0 Td[(HPkj(WkjZkj+WkjZkj). 7: Correctorstep-lengthComputecsimilartostep4with(W,Z)replacedby(W,Z). 8: Update(W,Z)tothenextiterate(W+,Z+).W+kj=Wkj+cWkj,Z+kj=Zkj+cZkj. 9: Updatethestep-lengthparameterby+=0.9+0.08c. problemismodiedas minW,H1 2mXi=1kVi)]TJ /F7 7.97 Tf 18.66 14.94 Td[(rXj=1Wjhjik2F+ 2X(k,l)2rXj=1(hjk)]TJ /F8 11.955 Tf 11.96 0 Td[(hjl)2(4) withnonnegativeconstraintsWj0andhji0,foralli,j.Thersttermissimplytheobjectivefunctionin(3)givenVasarowvector.Thesecondtermisthespatialsmoothness(soft)constraintthatrequiresneighboringpixelstohavesimilarcoefcients,andintheequationabove,denotestheedgesetofthediscreteimagegraph,andthecouplingparameter.Theoptimizationproblem(16)canbesolvedusingthesame 69

PAGE 70

Figure4-1. Segmentationexperimentonsynthetictensoreld.Left:Inputtensoreldwithoutnoise.Right:Segmentationaccuracyvs.differentlevelsofaddedGaussiannoiseN(0,)withcovariance. alternatingmethoddiscussedaboveasthesecondtermisalsoquadraticinH.OncethecoefcientmatrixHhasbeendetermined,weclusterthecolumnsofHtoproducethediffusiontensorimagesegmentation.Inourimplementation,weuseK-meansforthislastclusteringstep. 4.4Experiments Inthissection,wepresenttwosetsofexperimentalresults.Therstexperimentisondiffusiontensorimagesegmentationusingthesegmentationalgorithmoutlinedintheprevioussection.Forthesecondsetofexperiments,weworkwithmultipleimages,andtheresultsdemonstratedthat,asforscalar-valuedimages,meaningfulpartscanbediscoveredordetectedusingnonnegativefactorization. 4.4.1Diffusiontensorimagesegmentation SyntheticTensorImagesInthisexperiment,wetesttheaccuracyofthesegmentationalgorithmusingsynthetictensorimageswithvariouslevelsofaddednoise.Werstgeneratethe3232tensorimageTshowninFigure 4-1 thatcontainsonlytwotensors:diagonalmatricesdiag(0.5,0.25,0.25)anddiag(0.25,0.5,0.25),andthisdenesthegroundtruthofthesegmentation.DifferentlevelsofGaussiannoiseN(0,)areaddedtothetensorimageTtogeneratenoisytensorimages,andwe 70

PAGE 71

Figure4-2. Diffusiontensorimagesegmentation.FirstRow:Spinalcordofarat.SecondRow:Corpuscallosumofarat.ThirdRow:Hippocampusofarat.ColumnsA)Diffusiontensorimages.B)SegmentationresultsusingK-means.C)Segmentationresultsusingourmethod.Segmentsarecolor-coded. comparethesegmentationaccuracyofourmethodwiththesegmentationalgorithmbasedonclusteringthepixelsusingK-means(onthetensors).Inthisexperiment,thesegmentationaccuracyisdenedbythepercentageofcorrectlysegmentedpixels,andthecomparisonacrossdifferentnoiselevelsisplottedinFigure 4-1 .TheresultclearlyshowstherobustnessofourmethodwhencomparedwithK-means,especiallyinthepresenceofsubstantialamountofnoise. DiffusionTensorImagesWenextpresentsegmentationresultsonthreerealdiffusiontensorimagesshowninFigure 4-2 .TheseareDTIimagesofthespinalcord,corpuscallosumandanisolatedhippocampusofarat.ThedatawereacquiredusingaPGSEwithTR=1.5s,TE=28.3ms,bandwidth=35Khz.Twenty-onediffusionweighted 71

PAGE 72

Figure4-3. Nonnegativefactorizationonmultiplesynthetictensorelds.A)VisualizationofthedatamatrixV,eachofwhosecolumnsrepresentsa53tensorimage.B)VisualizationofthebasismatrixW,eachofwhosecolumnsrepresentsabasistensorimage. imageswithab-valueof1250s/mm2werecollected.Thesizesoftheregionsofinterestforratspinalcord,corpuscallosumandhippocampusare7161,10174and7139,respectively.Inthisexperiment,thenumberofclustersfortheseimagesarefour,eightandseven,respectively,andthenumberofbasiselements(columnsinW)issettoveforallthreeimages.Again,wecomparedourmethodwiththeK-meansbasedsegmentationandtheresultsconsistentdemonstratethatourmethodcanproduceanatomicallymoreaccurateandmeaningfulsegmentationthanthestraightforwardK-meansclustering.Finally,wenotethatbothalgorithmsuseK-meansforclusteringpixels.However,inourmethod,K-meansisappliedonlytothecoefcientvectorsinH,whileinthecomparisonmethod,K-meansisapplieddirectlytothetensors.Whiletheclusteringpowerofthenonnegativefactorizationformatricesarewell-knownforscalar-valueddata(e.g.,[ 17 ]),ourexperimentalresultsprovidetherstconvincingevidenceofitsclusteringpowerfortensor-valueddataasthetwosetsofexperimentshaveshownthatthereisnoreasonstoexpectthatadirectclusteringontensorswouldproducedesiredsegmentationresults.However,directclusteringoncoefcientvectorsdoesyieldsatisfactoryresults. 4.4.2Nonnegativefactorizationwithmultipletensorelds Forasuitablecollectionofimages,itiswell-knownthatNMFhastheabilitytoautomaticallydeterminesimilardecompositionsoftheimagesintotheircommonparts 72

PAGE 73

andprovidepart-basedrepresentationsfortheimagesinthecollection[ 18 ].Inthisexperimentusingsyntheticdata,weshowthatnonnegativefactorizationalsohassimilarcapabilityfortensorimages.Twenty-seventensoreldswithsize53(15pixels)aregeneratedandformthecolumnsofthedatamatrixVasshowninFigure 4-3 (a).WefactorVintothecoefcientmatrixHandthebasismatrixWshowninFigure 4-3 (b).Ourfactorizationalgorithmcorrectlydeterminestheninebasiselements(columnsofW)requiredtoformthedatamatrixV,andthecoefcientmatrixreturnedbyouralgorithmisasparsematrix.Furthermore,theL2factorizationerrorislessthan8.5710)]TJ /F5 7.97 Tf 6.59 0 Td[(10. Finally,wepresentapreliminaryresultonapplyingournonnegativefactorizationmethodtoautomaticallydiscoverandsegmentanatomicallyimportantregionsfromacollectionof53ratbraindiffusiontensorimages.Inthepreprocessingstep,allimagesarealignedusingsimilaritytransforms,andineachimage,aregionofinterestofsize46277ismanuallyselected.TheleftcolumnofFigure4displaysvesampleslicesfromtheinput3Ddiffusiontensorimages.Weapplythefactorizationalgorithmwithr=5(vebasiselements)tothiscollectionofDTIs,andonesampleslicefromeachofthevebasisimagesfoundbythealgorithmisshownontherightcolumninFigure4-4.Importantanatomicalregionssuchaswhitematter,putamenandnucleusaccumbensareclearlyrepresentedinthevebasisimages. 73

PAGE 74

Figure4-4. Leftcolumn:Fivesampleslicesfrom53inputratbraindiffusiontensorimages.Rightcolumn:Sampleslicesfromthevebasistensorimagesproducedbytheproposednonnegativefactorizationalgorithm. 74

PAGE 75

CHAPTER5DICTIONARYLEARNINGONRIEMANNIANMANIFOLDS 5.1Motivation Dictionarylearning,whichseekstondacollectionofatomsforsparserepresentationoftheinputdata,hasbeenwidelyusedincomputervisionapplicationssuchasimagerecognition,classicationandrestoration(e.g.,[ 2 ]).Underthismodel,eachdatapointisassumedtobegeneratedlinearlyusingonlyasmallnumberofatoms,andthislinearsparsityassumptionisresponsibleformuchofitsgeneralizationpowerandsuccess.However,theunderlyinglinearprocessrequiresthatthedatapointsaswellastheatomsaretreatedasvectorsinsomevectorspaceRd,andthedictionaryislearnedfromtheinputdatausingonlythevectorspacestructure(anditsassociatedinnerproduct).FormanyapplicationsincomputevisionandmedicalimageanalysisthatinvolvedatapointsbelongingtosomeknownRiemannianmanifoldssuchasthespaceofsymmetricpositive-denitematrices[ 23 ],hyper-spheresforparameterizingsquare-rootdensities[ 61 ],StiefelandGrassmannmanifolds[ 10 ],etc.,theexistingextrinsicapproachesthatcompletelyignorethepotentiallyimportantintrinsicstructureimpliedbythedataisclearlyinadequateandunsound(iftheycanbeappliedatall).Toremedythisdeciencyandinadequacy,wetakeasmallstepinthisdirectionbyinvestigatingtheproblemofextendingdictionarylearningframeworktoincorporateintrinsicgeometryimpliedbytheinputdata. Theapplicabilityandsuitabilityofapplyingexistingdictionarylearningmethodstosolvevisionproblemsthathavetodealwithmanifold-valueddatacanbetwothornyissuestoconsider.First,asaprerequisite,thedatamanifoldmustadmitanembeddingintosomeRdinordertobeabletoapplytheexistingdictionarylearningmethods.However,formostmanifolds,suchasGrassmannandStiefelmanifolds,theresimplydoesnotexistknowncanonicalembeddingintoRd(orsuchembeddingisdifculttocompute).Second,eveninthecasewhentheexistingmethodcanbeapplied, 75

PAGE 76

duetotheirextrinsicviewpoint,importantintrinsicpropertiesofthedatamaynotberepresentedinthedictionary.Thiscanbeillustratedbyasimpleexamplethatitispossiblethattwopointsx,yonthemanifoldMhavealargegeodesicdistanceseparatingthembutundertheembeddingi:M!Rd,i(x),i(y)hasasmalldistanceinRd.Therefore,sparsecodingusingdictionarylearnedinRdislikelytocodei(x),i(y)(andhencex,y)usingthesamesetofatomswithsimilarcoefcients.Clearly,thiswillbeundesirableandunsatisfactoryiftheapplicationsrequiretaskssuchasclassicationandclustering,forwhichonewouldpreferthesparsecodingtoreectsomedegreeofactualsimilarity(i.e.,geodesicdistance)betweenthetwosamplesx,y. WhiletheaboveexampleprovidesthemotivationforseekinganextensiontotheexistingdictionarylearningframeworktothemoregeneralRiemannainsetting,itisbynomeansobvioushowtheextensionshouldbecorrectlyformulated.LetMdenotetheRiemannianmanifoldonwhichacollectionofdatapointsx1,,xnaregiven.Attheminimum,thegoalofdictionarylearningonMistocomputeacollectionofatomsfa1,,amgM,alsopointsonM,suchthateachdatapointxicanbegeneratedusingonlyasmallnumberofatoms(sparsity).IntheEuclideansetting,thisisusuallyformulatedas minD,w1,,wnnXi=1kxi)]TJ /F8 11.955 Tf 11.96 0 Td[(Dwik2+Sp(wi),(5) whereDisthematrixwithcolumnscomposedoftheatomsai,withesparsecodingcoefcientsandSp(wi)thesparsitypromotingterm.Oneimmediatetechnicalhurdlethatanysatisfactorygeneralizationneedstoovercomeisthelackofagloballinearstructurethatwillallowthedatatobegeneratedfromtheatoms.Instead,theRiemanniangeometryprovidesonlylocallinearstructuresthroughtheRiemannianexponentialandlogarithmmaps[ 59 ],andbymovingtothemoregeneralRiemanniansetting,weessentiallytradetheuniquegloballinearstructurewithinnitelymanylocallinearstructures,whichisthemainsourceofthevarioustechnicaldifcultiespresentin 76

PAGE 77

ourgeneralization.However,thisdiversityoflinearstructuresalsoprovidesuswithanopportunitytoformulatethedictionarylearningusingadataspecicapproach. Specically,wewillformallymodifyeachsummandinEquation 5 sothatthesparsecodingofadataxiwithrespecttotheatomsfa1,,amgMisnowobtainedbyminimizing minwikmXj=1wijlogxiajk2xi+Sp(wi),(5) withtheimportantafneconstraintthatPmj=1wij=1,wherewi=(wi1,...,wim)T.Thatis,weareusingtheRiemannianexponentialandlogarithmmapsateachdatapointxtodenethegenerativeprocess,andthesparseapproximationofagivendatapointisrstcomputedinitstangentspaceTxMandthenrealizedonMbyapplyingtheexponentialmap.Weremarkthatthisformulationisentirelycoordinate-independentsinceeachlogxiajiscoordinate-independent,andEquation 5 canbeminimizedusinganylocalchartanditsassociatedbasisforTxiM(witharesultthatwillbeindependentofthesechoices).Furthermore,althoughthisafneconstraintisnotpresentinmostexistingdictionarylearningalgorithms,itistheconsequenceofourdata-specicapproach:asubspaceSinagivencoordinatessystemisrepresentedasanafnesubspaceinadifferentcoordinatessystemwithadifferentorigin.Inshort,Equation 5 isthedirectgeneralizationoflinearsparsityconditionwiththeexceptionnowtheoriginhasbeenmovedtothedatapointxi.Computationally,theresultingoptimizationproblemusingEquation 5 canbeeffectivelyminimized.Inparticular,thegradientofthecostfunctioninsomecasesadmitsclosed-formformulasandingeneral,theycanbeevaluatednumericallytoprovidetherequiredinputsforthegradient-basedoptimizationalgorithmonmanifolds[ 19 ]. Tovalidatetheproposedmethod,wehaveappliedthedictionarylearningframeworktoseveralclassicationproblemsincomputervisionandmedicalimaginganalysis.Thefeaturevectorsspecicfortheseproblemsbelongtoseveralwell-knownRiemannianmanifolds,andtheexperimentalresultsshowthatthedictionarylearned 77

PAGE 78

usingtheproposedframeworkcananddoesindeedproviderealimprovementswhencomparedwithotherdirectapproaches.Finally,beforeembarkingonpresentingthedetailsoftheproposedmethod,wementionthatfortechnicalreasons,wewillassumethattheRiemannianmanifoldMunderconsiderationiscompleteandforeachpointx2M,theRiemannianlogmaplogxatxcanbedenedonMoutsideofasubsetofcodimensionatleastone.Thistechnicalconditionistoensurethatforalmosteveryarbitrarypointy2M,logxycanbeuniquelydened,andfortheRiemannianmanifoldsconsideredinthischapter,thisconditionwillalwaysbesatised. 5.2RelatedWork Thesparsecodingframeworkwhichapproximatestheobservedsignalswiththelinearcombinationofasmallnumberofatomsfromaxeddictionaryhasbeensuccessfullyusedinmanycomputervisionandpatternrecognitionapplications[ 77 ].DictionarylearningproposedintheseminalworkbyOlshausenandField[ 50 ]seekstolearnacodebookleadingtothesparsedecompositionoftheobservations.Inrecentresearch,ithasbeenshownthatthelearneddictionarygivessuperiorperformanceinmanycomputervisiontaskssuchasimagedenoising[ 2 ]andimageclassication[ 84 ].In[ 42 ],Mairaletal.introducedanonlinedictionarylearningalgorithmwhichtsverylargedatasetsbetterthanthebatchprocedures.Thecombinationofgroup-structuredsparsityregularizationandonlinedictionarylearninghasbeendiscussedin[ 63 ].Sraetal.proposedageneralizeddictionarylearningmethodforsymmetricpositivedenite(SPD)matrices[ 60 ].TheydidnotconsidertheRiemannianstructureoftheSPDmatrixspaceandtheSPDmatrixwasapproximatedbynonnegativeweightedsumofrank-1matrices. 5.3ABriefonRiemannianManifold WenowpresentaverybriefintroductiontosomewellknownmaterialfromRiemmanianmanifoldsandreferthereaderto[ 59 ]fordetails.AmanifoldMofdimensiondisatopologicalspacethatislocallyhomeomorphictoopensubsetsof 78

PAGE 79

theEuclideanspaceRd.Withagloballydeneddifferentialstructure,manifoldMbecomesadifferentiablemanifold.Thetangentspaceatx2MindicatedbyTxMisavectorspacethatcontainsallthetangentvectorstoMatx.ARiemannianmetriconMassociatestoeachpointx2Madifferentiablevaryinginnerproducth,ixinthetangentspaceTxM.Letxi,xjbetwopointsonthemanifoldM.Thegeodesiccurve:[0,1]!Misasmoothcurvewithminimumlengthconnectingxiandxj.Letv2TxMbeatangentvectortothemanifoldatx.Thereexistsauniquegeodesicvsatisfyingv(0)=xwithinitialtangentvectorv.Thentheexponentialmapisdenedasexpx(v)=v(1).Theinverseoftheexponentialmapexpxiscalledthelogmapanddenotedbylogx:M!TxM.Giventwopointsxi,xj2M,logximapsthepointxjtotheuniquetangentvectoratxithatistheinitialvelocityofthegeodesicwith(0)=xiand(1)=xj.Thegeodesicdistancebetweenxiandxjiscomputedbydist(xi,xj)=klogxi(xj)kxi. Givenasetofpointsx1,...,xnonaRiemannianmanifoldM,wecancalculatetheintrinsicmeanonMbysolvingtheoptimizationproblem x=argminx2M1 nnXi=1dist(x,xi)2.(5) ThepropertiesoftheRiemannianintrinsicmeanhavebeenstudiedbyKarcher[ 35 ].Ingeneral,theintrinsicmeanxmaynothaveaclosedform.Thusoptimizationalgorithmsonmanifoldsareusuallyusedtocomputex[ 1 51 ]. 5.4DictionaryLearningonRiemannianManifolds Givenacollectionofsignalsx1,...,xn2Rd,classicaldictionarylearningmethods[ 2 40 49 ]trytondadictionaryD2RdmwhichincludesmatomssuchthateachsignalxicanberepresentedasasparselinearcombinationoftheseatomsxiDwi,wherewi2Rm.Asothers[ 42 84 ],weformulatethedictionarylearningproblemusingl1regularizationonwi minD,winXi=1)]TJ /F2 11.955 Tf 5.48 -9.69 Td[(kxi)]TJ /F8 11.955 Tf 11.96 0 Td[(Dwik22+kwik1(5) 79

PAGE 80

Figure5-1. Thebluecurvedenotesaone-dimensionalmanifoldM.xisadatapointonMandtheredpointsrepresentthedictionaryatoms.Left:thelinearapproximation^xwiththedictionaryatomsmaynotbeonthemanifoldM.Right:ourapproachprojectstheatomstothetangentspaceatxandperformsthelinearcombinationonthevectorspaceTxM.Theexponentialmapguaranteestheapproximation^xisalwaysonM. whereisaregularizationparameter.Recentextensionsoftheclassicaldictionarylearningincludeonlinedictionarylearning[ 42 ]andusingdifferentregularizationtermssuchasgroup-structuredsparsity[ 63 ]andlocalcoordinateconstraint[ 74 ]. WegeneralizetheclassicaldictionarylearningtechniquesinEuclideanspacetoRiemannianmanifolds.SupposeMrepresentsaRiemannianmanifold.Letx1,...,xn2MbeacollectionofndatapointsonthemanifoldM.Leta1,...,am2MbeatomsofthelearneddictionaryD=fa1,...,amg.BecauseofthegeometricstructureofM,wecannotusethelinearcombinationofatoms^xi=Pmj=1wijajtorepresentthedataxi,since^ximaynotevenbeonthemanifoldM.ThusweusethegeodesiclinearinterpolationonRiemannianmanifold,andxicanberepresentedby ^xi=expxi mXj=1wijlogxi(aj)!(5) whereexpxiandlogxiareexponentialandlogarithmicmapsatxi,andwij2Rareweights.Notethatinordertoapproximatedatapointxi,weprojectalltheatomsinthedictionarytothetangentspaceatxi.BecauseTxiMisavectorspace,wecanperformlinearcombinationvi=Pmj=1wijlogxi(aj)onTxiM.Thentheapproximation^xicanbe 80

PAGE 81

representedastheexponentialmapofviatxi.A1DcaseisillustratedinFigure 5-1 .Wehopetobuildadictionarythatminimizesthesumofreconstructionerrorforeachdatapoint.Herewedene Edata=nXi=1dist(xi,^xi)2=nXi=1klogxi(^xi)k2xi=nXi=1kmXj=1wijlogxi(aj)k2xi.(5) Byusingthel1sparsityregularization,thedictionarylearningonthemanifoldMcanbewrittenasthefollowingoptimizationproblem minW,DnXi=1kmXj=1wijlogxi(aj)k2xi+kWk1s.t.mXj=1wij=1,i=1,...,n(5) whereW2Rnmandthe(i,j)entryofWiswrittenaswij.TheafneconstraintPmj=1wij=1usedinourformulationhasalsoappearedinseveralrecentpapersondictionarylearningsuchassparsesubspaceclustering[ 20 ]andlocalcoordinatecoding[ 86 ].Theafneconstraintmeansthatweareusingafnesubspacestoapproximatethedatainsteadoftheusualsubspaces,whicharesimplyafnesubspacesbasedattheorigin.MovingawayfromvectorspacestogeneralRiemannianmanifolds,thereisnocorrespondingnotionoftheoriginthatcanbeusedtodenesubspaces,andthissimplegeometricfactrequirestheabandonmentoftheusualsubspacesinfavorofgeneralafnesubspaces.Wecanalsointroduceotherregularizationstoourdictionarylearningframeworkinsteadofthel1norm.Forexample,thelocalizationregularizationusedinlocalcoordinatecoding[ 86 ]canbegeneralizedtoEreg=Pni=1Pmj=1jwijjklogxi(aj)k2xi.Similartoclassicaldictionarylearningmethods,wecouldusetheiterativemethodtosolvethisoptimizationproblem: 1. Sparsecodingstep:xthedictionaryDandoptimizewithrespecttothecoefcientsW. 81

PAGE 82

2. Codebookoptimizationstep:xWandoptimizewithrespecttothedictionaryD. Therststepbecomesaregularsparsecodingproblemandmanyfastalgorithmscanbeusedtosolveit.Thesecondsubproblemismuchmorechallenging,becausethedictionaryatomsareonthemanifoldMandtheoptimizationmethodsinEuclideanspacearenotappropriateanymore. Algorithm9K-meansonmanifoldM Input:AsetofdataX=fx1,...,xngonthemanifoldMandthenumberofclustersK. Output:Clustersc1,...,cKandapartitionoftheinputdataL=fl1,...,lng,whereeachlabelli2f1,...,Kg 1.RandomlyselectKelementsfromXastheinitializationofclusters. 2.Foreachdataxi,setthelabelasli=argminkdist(xi,ck). 3.Foreachclusterk,letXk=fxijli=kg,themeancanbecomputedbyck=argminc2MXx2Xkdist(c,x)2. 4.StopifLdoesnotchange,otherwisegotostep2. WepresentalinesearchbasedalgorithmonRiemannianmanifoldtoupdatethedictionaryD.Letf(a1,...,am)bethecostfunctiontobeminimized.Atrst,weneedtoinitializetheatomsindictionary.Onepossibleinitializationisusingmclustersofthedatax1,...,xngeneratedbyK-meansalgorithm.Becauseofthemanifoldstructure,K-meansonMshowninAlgorithm9isalittledifferentfromthetraditionalK-meansalgorithminEuclideanspace.Thenweusethelinesearchonmanifoldtooptimizef(a1,...,am).Thebasicideaistondadescentdirectionvonthetangentspace,andthenwalkastepalongthegeodesicwhoseinitialvelocityisv.ThedetailsarelistedinAlgorithm10.Theconvergenceanalysisofthelinesearchmethodonmanifoldisdiscussedin[ 1 ]. 82

PAGE 83

Inthenexttwosections,wegivetwoconcreteexamplesanduseourdictionarylearningframeworkonsymmetricpositivedenite(SPD)matrixspaceandsquare-rootdensityfunctionspace. Algorithm10LinesearchonRiemannianmanifold Input:AsetofdataX=fx1,...,xngonthemanifoldM,coefcientsW2Rnmandinitialdictionaryatomsa01,...,a0m. Output:(a1,...,am)thatminimizethecostfunctionf(a1,...,am). 1.Setscalars>0,,2(0,1)andinitializek=0. 2.Computegradf(ak1,...,akm)=(@f(ak1) @a1,...,@f(akm) @am) 3.Pickk=(k1,...,km)=)]TJ /F4 11.955 Tf 9.3 0 Td[(gradf,whereki2TakiM. 4.Findthesmallestssuchthatf(expak1(sk1),...,expakm(skm))f(ak1,...,akm))]TJ /F7 7.97 Tf 17.07 14.94 Td[(mXi=1skkikaki. 5.Setak+1i=expaki(ski),i=1,...,m. 6.Stopiffdoesnotchangemuch,otherwisesetk=k+1andgobacktostep2. 5.5SymmetricPositiveDeniteMatrixSpace LetP(d)representtheddsymmetricpositivedenite(SPD)matrixspacewhichformsaconnectedRiemannianmanifold[ 23 51 ].TheRiemanniansymmetricspaceP(d)arisesnaturallyfromthegroupactionofthegenerallineargroupGL(d).ThegroupactionofGL(d)onP(d)isthetransformation:GL(d)P(d)!P(d)givenbyG(X)=GXGT,whereG2GL(d)isaddinvertiblematrixandX2P(d)isaddSPDmatrix.ThetangentspaceTXP(d)ateverypointX2P(d)canbeidentiedwithSym(d),whichrepresentstheddsymmetricmatrixspace. LetY,Z2TMP(d)betwotangentvectorsatM2P(d).TheRiemannianmetricgivestheinnerproductas hY,ZiM=tr(YM)]TJ /F5 7.97 Tf 6.58 0 Td[(1ZM)]TJ /F5 7.97 Tf 6.59 0 Td[(1)(5) 83

PAGE 84

wheretrrepresentsthematrixtrace.ThegeodesicpassingthroughM2P(d)inthedirectionofY2TMP(d)withrespecttotheinnerproductdenedaboveisuniquelygivenby (t)=Gexp)]TJ /F8 11.955 Tf 5.48 -9.69 Td[(G)]TJ /F5 7.97 Tf 6.59 0 Td[(1YG)]TJ /F7 7.97 Tf 6.59 0 Td[(TtGT(5) whereexprepresentsthematrixexponentialandG2GL(d)isanysquarerootofMsuchthatM=GGT.ThereforetheRiemannianexponentialmapatMwhichmapsYtoapointinP(d)canberepresentedas expM(Y)=Gexp)]TJ /F8 11.955 Tf 5.48 -9.68 Td[(G)]TJ /F5 7.97 Tf 6.59 0 Td[(1YG)]TJ /F7 7.97 Tf 6.58 0 Td[(TGT.(5) GiventwopositivedenitematricesX,M2P(d),theRiemannianlogarithmicmaplogM:P(d)!TMP(d)isdenedas logM(X)=Glog(G)]TJ /F5 7.97 Tf 6.59 0 Td[(1XG)]TJ /F7 7.97 Tf 6.58 0 Td[(T)GT(5) wherelogrepresentsthematrixlogarithm.UsingthedenitionoflogmapinP(d),thegeodesicdistancebetweenMandXcanbecomputedby dist(M,X)=klogM(X)k=q tr(log2(G)]TJ /F5 7.97 Tf 6.58 0 Td[(1XG)]TJ /F7 7.97 Tf 6.59 0 Td[(T)).(5) Inrecentyears,thereisagrowinginterestonsparsecoding,dictionarylearningandnonnegativefactorizationofSPDmatrixdata[ 58 80 ].However,intheseworksanSPDmatrixisrepresentedasanonnegativelinearcombinationofthedictionaryatoms,whichdoesnotrespecttheintrinsicRiemannianstructureoftheSPDmatrixspacediscussedabove.Inthissection,wepresentadictionarylearningmethodforsymmetricpositivedenitematrices,whichconsiderstheintrinsicgeometryofP(d). LetX1,...,Xn2P(d)beacollectionofddSPDmatrices.A1,...,Am2P(d)arematomsinthelearneddictionaryD.Wisanmmatrix.DifferentfromthelinearcombinationapproximationofapointXiinP(d),ourrepresentationwithintrinsic 84

PAGE 85

Riemannianstructurecanbewrittenas ^Xi=Giexp mXj=1wijlog(G)]TJ /F5 7.97 Tf 6.58 0 Td[(1iAjG)]TJ /F7 7.97 Tf 6.59 0 Td[(Ti)!GTi(5) whereGi2GL(d)suchthatGiGTi=Xi. Ifwejustusel1regularization,thedictionarylearningforSPDmatricesbecomes minW,DnXi=1mXj=1mXk=1wijwiktr(LijLik)+kWk1s.t.mXj=1wij=1,i=1,...,n(5) whereLij=log(G)]TJ /F5 7.97 Tf 6.59 0 Td[(1iAjG)]TJ /F7 7.97 Tf 6.58 0 Td[(Ti),i=1,...,nandj=1,...,m.Thisoptimizationproblemcanbesolvedusingthealgorithmdiscussedinsection4. 5.6Square-rootDensityFunctionSpace Probabilitydensityfunctions(pdfs)asaclassofconstrainednon-negativefunctionshavebeenwidelyusedinmanycomputervisionapplicationssuchastextureanalysisandshapeanalysis.Withoutlossofgenerality,werestricttothepdfsdenedontheinterval[0,T]inthissection: P=fp:[0,T]!Rj8s,p(s)0,ZT0p(s)ds=1g.(5) TheFisher-RaometrichasbeenintroducedtostudytheRiemannianstructureformedbythestatisticalmanifoldin[ 52 ].Forapdfpi2P,theFisher-Raometricisdenedas hvj,vki=ZT0vj(s)vk(s)1 pi(s)ds(5) wherevj,vk2Tpi(P).TheFisher-Raometricisinvarianttoreparameterizationsofthefunctions.InordertomaketheresultingmanifoldeasytocomputeonwithRiemannianoperations,thesquarerootdensityrepresentation =p pwasused[ 61 ].Thespaceof 85

PAGE 86

squarerootdensityfunctionsisdenedas =f :[0,T]!Rj8s, (s)0,ZT0 2(s)ds=1g(5) Aswecansee,formsaconvexsubsetoftheunitsphereinaHilbertspace.ThentheFisher-Raometriccanbeobtainedas hvj,vki=ZT0vj(s)vk(s)ds(5) wherevj,vk2T iaretangentvectors.Givenanytwofunctions i, j2,thegeodesicdistancebetweenthesetwopointsis dist( i, j)=cos)]TJ /F5 7.97 Tf 6.59 0 Td[(1(h i, ji)(5) whichisjusttheanglebetween iand jontheunithypersphere.Thegeodesicat iwithadirectionv2T i()isdenedas (t)=cos(t) i+sin(t)v jvj.(5) Thentheexponentialmapcanberepresentedas exp i(v)=cos(jvj) i+sin(jvj)v jvj.(5) Toensuretheexponentialmapisabijection,werestrictjvj2[0,).Thelogmapisthengivenby log i( j)=ucos)]TJ /F5 7.97 Tf 6.59 0 Td[(1(h i, ji)=p hu,ui(5) whereu= j)-222(h i, ji i. Usingthegeometricexpressionsofdiscussedabove,wecanperformthedictionarylearningonsquare-rootdensityfunctions.Letx1,...,xn2beacollectionofsquare-rootdensityfunctions,andai,...,am2beatomsinthedictionaryD.Wisa 86

PAGE 87

nmmatrix.Ifweusel1regularization,ourdictionarylearningframeworkbecomes minW,DnXi=1kmXj=1wijcos)]TJ /F5 7.97 Tf 6.59 0 Td[(1(hxi,aji)uij juijjk2xi+kWk1s.t.mXj=1wij=1,i=1,...,n.(5) whereuij=aj)-250(hxi,ajixi.Thisoptimizationproblemcanbeefcientlysolvedusingthealgorithmpresentedinsection5.4. 5.7Experiments Weevaluatedourdictionarylearningalgorithmsonclassicationexperiments.Letx1,...,xn2Mbeacollectionofntrainingdata.AtrstacodebookD=fa1,...,amgislearnedfromthetrainingdatabyusingthedictionarylearningtechniqueproposedinthischapter,whereeachaidenotesanatominthedictionaryD.Thenforeachxi,wegeneratethesparsecodewi2RmbasedonthelearneddictionaryD.NotethatthesparsecodingonRiemannianmanifoldisjustonestepinthedictionarylearningalgorithmwiththecodebookDxed.Similarto[ 84 ],alinearSVMistrainedforclassication.Inthetestingstage,lety1,...,yk2Mbeasetofktestingdata.Thesparsecodingcanbeperformedonthetestingdatainasimilarmanner.Giveneachsparsecodeasafeaturevector,thetrainedSVMmodelcanbeusedtoclassifythetestingdata. Inourexperiments,threemethodsareusedforcomparisonpurposes. GeodesicK-nearestneighbor(GKNN) SupportVectorMachine(SVM)onvectorizeddata SVMonsparsecodeswiththecodebooktrainedbyKSVD GKNNistheK-nearestneighborclassicationusingtheRiemanniangeodesicdistanceinsteadofthetraditionalEuclideandistance.It'scommonlyusedinliteratureforclassicationonmanifolds.Kissetto5inourexperiments.SincethefeaturevectorforSVMshouldbeinEuclideanspace,wevectorizethedataonmanifoldsandfeed 87

PAGE 88

Figure5-2. SamplesofBrodatztexturesusedinourexperiments. themtotheSVMclassier.TheLIBSVM[ 11 ]packageisusedinourexperiments.WealsoapplythepopularKSVDmethod[ 2 ]totrainadictionaryforsparsecoding,andthenalinearSVMonsparsecodesisemployedinclassication.Bycomparingtheproposedmethodwiththedictionarylearningmethodinvectorspace,wetrytodiscoverwhetherkeepingtheRiemannianstructureisimportantintheclassication.Theexperimentalcomparisonsareperformedontwopublicavailabledatasets:BrodatztexturedatasetandOASISbrainMRIdataset. 5.7.1Brodatzdataset Table5-1. ThetextureclassicationaccuracyfordifferentmethodsontheBrodatzdataset. SVMKSVD+SVMGKNNProposed 16classes93.3694.9295.7099.0232classes88.6790.8291.1195.70 Inthisexperiment,weevaluatethedictionarylearningalgorithmforSPDmatrices.Weuseregioncovariancedescriptorintroducedin[ 70 ]fortextureclassicationon 88

PAGE 89

Brodatzdataset[ 9 ].Withthesimilarexperimentalsettingin[ 58 ],weconstruct16-textureand32-texturesetsusingtheimagesfromBrodatzdataset.SomesampletexturesareshowninFigure 5-2 .Each256256textureimageisseparatedinto64blocksofsize3232.Insideeachblock,wecanconstructa55regioncovariancematrixusingthefeature F=I,@I @x,@I @y,@2I @x2,@2I @y2T.(5) Fortheexistenceofthederivativeateachpixel,thetextureimagesarepre-processedbygaussianlters.Tomakesureeachcovariancematrixstrictlypositivedenite,weaddE,whereisaverysmallpositiveconstantandEdenotestheidentitymatrix.Fork-classproblem,wherek=16or32,weobtain64kSPDmatrices,whichareequallypartitionedtothetrainingsetandthetestingset.ThesizeofthedictionaryDissetto5k.ThetextureclassicationresultsarereportedinTable 5-1 .Ourdictionarylearningbasedmethodperformsbetterthanotheralgorithmsincomparison.Aswecansee,althoughtheproposedmethodalsousesalinearSVMtoclassifythefeaturevectorinEuclideanspace,thesparsecodefromlearneddictionarypreservesthegeometricinformationbetterthanthevectorrepresentationofSPDmatrix. 5.7.2OASISdataset Table5-2. TheclassicationaccuracyfordifferentmethodsontheOASISdataset. SVMKSVD+SVMGKNNProposed Yvs.M90.0491.2991.8498.34Mvs.O97.3298.08100100Ovs.Y98.1299.46100100Yvs.Mvs.O91.9794.0493.1898.62 Inthisexperiment,weevaluatethedictionarylearningalgorithmforsquarerootdensitiesonthefreelyavailableOpenAccessSeriesofImagingStudies(OASIS)database[ 43 ].OASIScontainsT1weightedMRbrainimagesfromacross-sectionalpopulationof416subjects.EachMRIscanhasaresolutionof176208176voxels.Theagesofthesubjectsrangefrom18to96.WedividetheOASISpopulationinto 89

PAGE 90

threegroups.175subjectsyoungerthan40yearsoldaredesignatedasyoung(Y),and66subjectsbetween40and60aredesignatedasmiddle-aged(M).Weconsiderother195subjectsolderthan60asoldadults(O).Becauseofthestructuraldifferenceinthebrainacrossdifferentagegroups,wecanusetheMRimagestoperformagegroupclassication.Atrst,wealignalltheMRimagesinOASISdatasetusingthenonrigidgroup-wiseregistrationmethoddescribedin[ 32 ].Foreachimage,weobtainadisplacementeld.Thenthehistogramofthedisplacementvectorsisconstructedineachimageasthefeatureforclassication[ 12 ].Inourexperiment,thenumberofbinsineachdirectionissetto444.Thusthe64dimensionalhistogramisusedasthefeaturevectorfortheKSVDandSVMclassication,whilethesquarerootofthehistogramisusedinGKNNandourdictionarylearningbasedmethod.Thelearneddictionaryincludes100atomsinourexperiment.Weuse5-foldcrossvalidationandreporttheclassicationresultsinTable 5-2 .NotethatallfourmethodsgiveprettygoodclassicationaccuracyforYoungvs.OldandMiddle-agedvs.Old,butourmethodoutperformsthecompetingmethodsintheYoungvs.Middle-agedexperiment.Becausetheregionalbrainvolumeandcorticalthicknessofadultsarerelativelystablepriortoreaching60years[ 44 ],theclassicationofyoungandmiddle-agedsubjectsismorechallenging. 90

PAGE 91

CHAPTER6CONCLUSIONS Effectivegroup-wiseimageanalysisusingthegeometricstructureofthedataisabigtopic.Inthisdissertation,wehavepresentedseveraltechniquesalongthisdirectionforbothEuclideanimagespaceandnon-Euclideanhigh-dimensionalambientspace.Below,wewillsummarizethecontributions. Oneofthesimplebuteffectivewaystorepresentacollectionofdataistoconstructafewtemplates.Forinstance,populark-meansclusteringconstructsktemplatesandusesthislowdimensionalrepresentationtomodelthewholedataset.Withthesimilaridea,wehaveproposedanovelapproachtomultipleimageatlasconstructionproblem.Ak-NNgraphisemployedtocapturethegeometricstructureoftheimagedataset.Unsupervisedandsemi-supervisedgraphpartitionmethodsenableustobuildmultipleatlasesforheterogeneousimagedatasets.Theatlasescomputedusingtheproposedalgorithmpreserveimportantimagefeaturesandgenerallyenjoybetterimagesharpness. Somedatasuchasdiffusiontensorimageanddeformationtensoreldlieinanon-Euclideanambientspace.Thustraditionalgroup-wiseimageanalysismethodsinvectorspaceareunabletobedirectlyusedtoprocessthesetensorvalueddata.WehavepresentedaRimanniangeometricframeworkfortensorvaluedimageanalysis,whichextendsPrincipalComponentAnalysis(PCA)tothespaceoftensorelds.Thegeodesicsubmanifoldgeneratedbyourmethodprovidesaneffectivedimensionalityreductionofthetensorvaluedimagedatainthehigh-dimensionalproductspace. SimilartoPCA,NonnegativeMatrixFactorization(NMF)alsousesthelinearcombinationofasetofbasestoapproximateeachimagedata.ThedifferenceisthatPCArequiresthebasestobeorthonormal,whileNMFusesnonnegativeconstraintsonbothbasesandcoefcients.WehaveproposedanovelmethodasageneralizationofNMFtoapproximateacollectionofdiffusiontensorimagesusingnonnegativelinear 91

PAGE 92

combinationsofbasistensorimages.Anefcientiterativeoptimizationalgorithmisproposedtosolvethisfactorizationproblem. Inrecentyears,learneddictionarycombinedwithsparsecodingbecomesamajortoolforgroup-wisedataanalysis.Usingthesparseregularization,itcanautomaticallyselectlinearsubspacesspannedbyasmallnumberofbasesfromthedictionarytoapproximatethedata.ExistingdictionarylearningalgorithmsrelyheavilyontheassumptionthatthedatapointsarevectorsinsomeEuclideanspace.However,inmanyvisionapplications,featuresanddatapointsoftenbelongtosomeRiemannianmanifoldwithitsintrinsicmetricstructurethatispotentiallyimportantandcriticaltotheapplication.WedevelopedaverygeneraldictionarylearningframeworkfordatalyingonsomeknownRiemannnianmanifolds.UsingthelocallinearstructuresfurnishedbytheRiemanniangeometry,wehaveproposedanoveldictionarylearningalgorithmthatcanbeconsideredasdata-specic,afeaturethatisnotpresentintheexistingmethods. Forthefuturework,weplantoexploremore,innovativerepresentationsforgroup-wiseimageanalysiswithsolidcomputationalmodelsandrobustoptimizationmethods. 92

PAGE 93

REFERENCES [1] Absil,P.A.,Mahony,R.,andSepulchre,R.Optimizationalgorithmsonmatrixmanifolds.PrincetonUnivPr,2008. [2] Aharon,M.,Elad,M.,andBruckstein,A.K-SVD:AnAlgorithmforDesigningOvercompleteDictionariesforSparseRepresentation.SignalProcessing,IEEETransactionson54(2006).11:4311. [3] Ashburner,J.,Hutton,C.,Frackowiak,R.,Johnsrude,I.,Price,C.,andFriston,K.Identifyingglobalanatomicaldifferences:deformation-basedmorphometry.HumanBrainMapping6(1998).5-6:348. [4] Basser,P.J.,Mattiello,J.,andLeBihan,D.MRdiffusiontensorspectroscopyandimaging.Biophysicaljournal66(1994).1:259. [5] Bazin,P.L.andPham,D.L.Statisticalandtopologicalatlasbasedbrainimagesegmentation.Proceedingsofthe10thinternationalconferenceonMedicalimagecomputingandcomputer-assistedintervention.2007,94. [6] Belkin,M.andNiyogi,P.Laplacianeigenmapsandspectraltechniquesforembeddingandclustering.Advancesinneuralinformationprocessingsystems1(2002):585. [7] .Semi-supervisedlearningonRiemannianmanifolds.MachineLearning56(2004).1:209. [8] Boyd,A.andVandenberghe,L.ConvexOptimization.CambridgeUniversityPress,2004. [9] Brodatz,P.Textures:aphotographicalbumforartistsanddesigners,vol.66.DoverNewYork,1966. [10] Cetingul,H.E.andVidal,R.IntrinsicmeanshiftforclusteringonStiefelandGrassmannmanifolds.IEEEConferenceonComputerVisionandPatternRecognition(CVPR).2009,1896. [11] Chang,Chih-ChungandLin,Chih-Jen.LIBSVM:Alibraryforsupportvectormachines.ACMTransactionsonIntelligentSystemsandTechnology2(2011):27:1:27. [12] Chen,T.,Rangarajan,A.,andVemuri,B.C.Caviar:Classicationviaaggregatedregressionanditsapplicationinclassifyingoasisbraindatabase.IEEEInter-nationalSymposiumonBiomedicalImaging:FromNanotoMacro(ISBI).2010,1337. [13] Cootes,T.F.,Edwards,G.J.,andTaylor,C.J.Activeappearancemodels.IEEETransactionsonPatternAnalysisandMachineIntelligence23(2001).6:681. 93

PAGE 94

[14] Cootes,T.F.,Taylor,C.J.,Cooper,D.H.,Graham,J.,etal.Activeshapemodels-theirtrainingandapplication.Computervisionandimageunderstand-ing61(1995).1:38. [15] Cormen,T.H.,Leiserson,C.E.,Rivest,R.L.,andStein,C.Introductiontoalgorithm-s.TheMITpress,2001. [16] Davis,B.,Fletcher,P.,Bullitt,E.,andJoshi,S.PopulationShapeRegressionFromRandomDesignData.IEEEInternationalConferenceonComputerVision(ICCV).2007. [17] Ding,C.,He,X.,andSimon,H.D.Ontheequivalenceofnonnegativematrixfactorizationandspectralclustering.Proc.SIAMDataMiningConf.2005. [18] Donoho,D.L.andStodden,V.Whendoesnon-negativematrixfactorizationgiveacorrectdecompositionintoparts?NIPS.2004. [19] Edelman,A.,Arias,T.A.,andSmith,S.T.Thegeometryofalgorithmswithorthogonalityconstraints.Arxivpreprintphysics/9806030(1998). [20] Elhamifar,E.andVidal,R.Sparsesubspaceclustering.CVPR.2009,2790. [21] Fischl,B.,Salat,D.H.,Busa,E.,Albert,M.,Dieterich,M.,Haselgrove,C.,vanderKouwe,A.,Killiany,R.,Kennedy,D.,Klaveness,S.,etal.WholeBrainSegmentation::AutomatedLabelingofNeuroanatomicalStructuresintheHumanBrain.Neuron33(2002).3:341. [22] Fletcher,P.,Venkatasubramanian,S.,andJoshi,S.RobuststatisticsonRiemannianmanifoldviathegeometricmedian.IEEEConferenceonComputerVisionandPatternRecognition(CVPR).2008. [23] Fletcher,P.T.andJoshi,S.Riemanniangeometryforthestatisticalanalysisofdiffusiontensordata.SignalProcessing87(2007).2:250. [24] Fletcher,P.T.,Lu,C.,Pizer,S.M.,andJoshi,S.Principalgeodesicanalysisforthestudyofnonlinearstatisticsofshape.IEEETransactionsonMedicalImaging23(2004).8:995. [25] Gaussier,E.andGoutte,C.RelationbetweenPLSAandNMFandimplications.ACMSIGIR.2005. [26] Gerber,S.,Tasdizen,T.,ThomasFletcher,P.,Joshi,S.,andWhitaker,R.Manifoldmodelingforbrainpopulationanalysis.Medicalimageanalysis14(2010).5:643. [27] Grippo,L.andSciandrone,M.Ontheconvergenceoftheblocknonlineargauss-seidelmethodunderconvexconstraints.OperationsResearchLetters26(2000).3:127. 94

PAGE 95

[28] Hagen,L.andKahng,A.B.Newspectralmethodsforratiocutpartitioningandclustering.IEEETransactionsonComputer-AidedDesignofIntegratedCircuitsandSystems11(1992).9:1074. [29] Hamm,J.,Ye,D.H.,Verma,R.,andDavatzikos,C.GRAM:Aframeworkforgeodesicregistrationonanatomicalmanifolds.MedicalImageAnalysis14(2010).5:633. [30] Hua,X.,Leow,A.D.,Parikshak,N.,Lee,S.,Chiang,M.C.,Toga,A.W.,JackJr,C.R.,Weiner,M.W.,andThompson,P.M.Tensor-basedmorphometryasaneuroimagingbiomarkerforAlzheimer'sdisease:anMRIstudyof676AD,MCI,andnormalsubjects.Neuroimage43(2008).3:458. [31] Jia,H.,Wu,G.,Wang,Q.,andShen,D.ABSORB:Atlasbuildingbyself-organizedregistrationandbundling.NeuroImage51(2010).3:1057. [32] Joshi,S.,Davis,B.,Jomier,M.,andGerig,G.UnbiasedDiffeomorphicAtlasConstructionforComputationalAnatomy.NeuroImage23(2004):S151S160. [33] Joshi,S.,Karthikeyan,S.,Manjunath,BS,Grafton,S.,andKiehl,K.A.Anatomicalparts-basedregressionusingnon-negativematrixfactorization.CVPR.2010. [34] J.Tenenbaum,deSilva,V.,andLangford,J.Aglobalgeometricframeworkfornonlineardimensionalityreduction.Science290(5500)(2000):2319. [35] Karcher,H.Riemanniancenterofmassandmolliersmoothing.Communicationsonpureandappliedmathematics30(1977).5:509. [36] Lee,D.D.andSeung,H.S.Learningthepartsofobjectsbynon-negativematrixfactorization.Nature401(1999):788. [37] Lenglet,C.,Rousson,M.,andDeriche,R.DTIsegmentationbystatisticalsurfaceevolution.IEEETrans.onMedicalImaging25(2006).6:685. [38] Leow,A.D.,Klunder,A.D.,JackJr,C.R.,Toga,A.W.,Dale,A.M.,Bernstein,M.A.,Britson,P.J.,Gunter,J.L.,Ward,C.P.,Whitwell,J.L.,etal.LongitudinalstabilityofMRIformappingbrainchangeusingtensor-basedmorphometry.Neuroimage31(2006).2:627. [39] Lepore,N.,Brun,C.,Chou,Y.Y.,Chiang,M.C.,Dutton,R.A.,Hayashi,K.M.,Luders,E.,Lopez,O.L.,Aizenstein,H.J.,Toga,A.W.,etal.Generalizedtensor-basedmorphometryofHIV/AIDSusingmultivariatestatisticsondeformationtensors.IEEETransactionsonMedicalImaging27(2008).1:129. [40] Lewicki,M.S.andSejnowski,T.J.Learningovercompleterepresentations.Neuralcomputation12(2000).2:337. [41] Li,S.Z.,Hou,X.W.,Zhang,H.J.,andCheng,Q.S.Learningspatiallylocalized,parts-basedrepresentation.CVPR.2001. 95

PAGE 96

[42] Mairal,J.,Bach,F.,Ponce,J.,andSapiro,G.Onlinelearningformatrixfactorizationandsparsecoding.TheJournalofMachineLearningResearch11(2010):19. [43] Marcus,D.,Wang,T.,Parker,J.,Gsernansky,J.G.,Morris,J.,andBuckner,R.OpenAccessSeriesofImagingStudies(OASIS):Cross-SectionMRIDatainYoung,MiddleAged,Nondemented,andDementedOlderAdults.JournalofCognitiveNeuroscience19(2007):1498. [44] Mortamet,B.,Zeng,D.,Gerig,G.,Prastawa,M.,andBullitt,E.EffectsofhealthyagingmeasuredbyintracranialcompartmentvolumesusingadesignedMRbraindatabase.MICCAI.2005. [45] Mumford,D.,Fogarty,J.,andKirwan,F.GeometricInvariantTheory.Springer,1994. [46] Munkres,J.Topology.PrenticeHall,2000. [47] Nene,S.A.,Nayar,S.K.,andMurase,H.Columbiaobjectimagelibrary(coil-20).TechnicalReport005-96,CUCS(1996). [48] Nocedal,J.andWright,S.J.Numericaloptimization.Springer,2000. [49] Olshausen,B.A.andField,D.J.Sparsecodingwithanovercompletebasisset:AstrategyemployedbyV1?Visionresearch37(1997).23:3311. [50] Olshausen,B.A.etal.Emergenceofsimple-cellreceptiveeldpropertiesbylearningasparsecodefornaturalimages.Nature381(1996).6583:607. [51] Pennec,X.,Fillard,P.,andN.Ayache.ARiemannianFrameworkforTensorComputing.Int.J.ComputerVision66(2006).1:41. [52] Rao,C.R.Informationandaccuracyattainableintheestimationofstatisticalparameters.Bull.CalcuttaMath.Soc37(1945).3:81. [53] Roweis,S.T.andSaul,L.K.Nonlineardimensionalityreductionbylocallylinearembedding.Science290(2000).5500:2323. [54] Sabuncu,M.,Balci,S.,andGolland,P.DiscoveringModesofanImagePopulationthroughMixtureModeling.MICCAI.2008. [55] Sabuncu,M.R.,Balci,S.K.,Shenton,M.E.,andGolland,P.Image-drivenpopulationanalysisthroughmixturemodeling.IEEETransactionsonMedicalImaging28(2009).9:1473. [56] Schwartzman,A.Randomellipsoidsandfalsediscoveryrates:Statisticsfordiffusiontensorimagingdata.Ph.D.thesis,StanfordUniversity,2006. [57] Shi,J.andMalik,J.Normalizedcutsandimagesegmentation.PatternAnalysisandMachineIntelligence,IEEETransactionson22(2000).8:888. 96

PAGE 97

[58] Sivalingam,R.,Boley,D.,Morellas,V.,andPapanikolopoulos,N.Tensorsparsecodingforregioncovariances.ECCV.2010,722. [59] Spivak,M.Acomprehensiveintroductiontodifferentialgeometry.PublishorperishBerkeley,1979. [60] Sra,S.andCherian,A.Generalizeddictionarylearningforsymmetricpositivedenitematriceswithapplicationtonearestneighborretrieval.MachineLearningandKnowledgeDiscoveryinDatabases.2011,318. [61] Srivastava,A.,Jermyn,I.,andJoshi,S.Riemanniananalysisofprobabilitydensityfunctionswithapplicationsinvision.IEEEConferenceonComputerVisionandPatternRecognition(CVPR).2007. [62] Sturm,J.F.UsingSeDuMi1.02,aMATLABtoolboxforoptimizationoversymmetriccones.Optimizationmethodsandsoftware11(1999).1:625. [63] Szabo,Z.,Poczos,B.,andLorincz,A.Onlinegroup-structureddictionarylearning.IEEEConferenceonComputerVisionandPatternRecognition(CVPR).2011,2865. [64] Talairach,J.andTournoux,P.Co-planarstereotaxicatlasofthehumanbrain,vol.147.ThiemeNewYork,1988. [65] ThomasYeo,BT,Sabuncu,M.R.,Desikan,R.,Fischl,B.,andGolland,P.Effectsofregistrationregularizationandatlassharpnessonsegmentationaccuracy.MedicalImageAnalysis12(2008).5:603. [66] Todd,MJ,Toh,KC,andTutuncu,RH.OntheNesterov-Todddirectioninsemideniteprogramming.SIAMJournalonOptimization8(1998).3:769. [67] Toh,KC,Tutuncu,RH,andTodd,MJ.Inexactprimal-dualpath-followingalgorithmsforaspecialclassofconvexquadraticSDPandrelatedproblems.Pac.J.Optim3(2007):135. [68] Turk,M.A.andPentland,A.P.Facerecognitionusingeigenfaces.IEEEComputerSocietyConferenceonComputerVisionandPatternRecognition(CVPR).1991,586. [69] Tutuncu,R.H.,Toh,K.C.,andTodd,M.J.Solvingsemidenite-quadratic-linearprogramsusingSDPT3.MathematicalProgramming95(2003).2:189. [70] Tuzel,O.,Porikli,F.,andMeer,P.Regioncovariance:Afastdescriptorfordetectionandclassication.EuropeanConferenceonComputerVision(ECCV).2006,589. [71] VanBenthem,M.H.andKeenan,M.R.Fastalgorithmforthesolutionoflarge-scalenon-negativity-constrainedleastsquaresproblems.Journalofchemometrics18(2004).10:441. 97

PAGE 98

[72] Vapnik,V.N.Statisticallearningtheory.Wiley,NewYork,1998. [73] Vemuri,BC,Ye,J.,Chen,Y.,andLeonard,CM.Imageregistrationvialevel-setmotion:Applicationstoatlas-basedsegmentation.Medicalimageanalysis7(2003).1:1. [74] Wang,J.,Yang,J.,Yu,K.,Lv,F.,Huang,T.,andGong,Y.Locality-constrainedlinearcodingforimageclassication.IEEEConferenceonComputerVisionandPatternRecognition(CVPR).2010,3360. [75] Wang,Z.andVemuri,B.C.DTIsegmentationusinganinformationtheoretictensordissimilaritymeasure.IEEETrans.onMedicalImaging24(2005).10:1267. [76] Weldeselassie,Y.T.andHamarneh,G.DT-MRIsegmentationusinggraphcuts.Proc.ofSPIEMedicalImaging:ImageProcessing.2007. [77] Wright,J.,Ma,Y.,Mairal,J.,Sapiro,G.,Huang,T.S.,andYan,S.Sparserepresentationforcomputervisionandpatternrecognition.ProceedingsoftheIEEE98(2010).6:1031. [78] Wu,G.,Jia,H.,Wang,Q.,andShen,D.SharpMean:Groupwiseregistrationguidedbysharpmeanimageandtree-basedregistration.NeuroImage(2011). [79] Wu,J.,Smith,WAP,andHancock,ER.Genderclassicationusingshapefromshading.Britishmachinevisionconference.2007. [80] Xie,Y.,Ho,J.,andVemuri,B.NonnegativeFactorizationofDiffusionTensorImagesandItsApplications.InformationProcessinginMedicalImaging(IPMI).Springer,2011,550. [81] Xie,Y.,Ho,J.,andVemuri,B.C.Imageatlasconstructionviaintrinsicaveragingonthemanifoldofimages.IEEEConferenceonComputerVisionandPatternRecognition(CVPR).2010. [82] Xie,Y.,Vemuri,B.,andHo,J.Statisticalanalysisoftensorelds.MedicalImageComputingandComputer-AssistedIntervention(MICCAI)(2010):682. [83] Xu,W.,Liu,X.,andGong,Y.Documentclusteringbasedonnon-negativematrixfactorization.ACMSIGIR.2003. [84] Yang,J.,Yu,K.,Gong,Y.,andHuang,T.Linearspatialpyramidmatchingusingsparsecodingforimageclassication.IEEEConferenceonComputerVisionandPatternRecognition(CVPR).2009,1794. [85] Yeo,B.T.T.,Vercauteren,T.,Fillard,P.,Peyrat,J.M.,Pennec,X.,Golland,P.,Ayache,N.,andClatz,O.DT-REFinD:Diffusiontensorregistrationwithexactnite-straindifferential.IEEETransactionsonMedicalImaging28(2009).12:1914. 98

PAGE 99

[86] Yu,K.,Zhang,T.,andGong,Y.Nonlinearlearningusinglocalcoordinatecoding.AdvancesinNeuralInformationProcessingSystems(NIPS)22(2009):2223. [87] Zhang,H.,Yushkevich,P.A.,andGee,J.C.Registrationofdiffusiontensorimages.IEEEComputerSocietyConferenceonComputerVisionandPatternRecognition(CVPR).vol.1.2004. [88] Zhu,X.Semi-supervisedlearningliteraturesurvey.TechnicalReport1530DepartmentofComputerScience,UniversityofWisconsin,Madison.(2005). 99

PAGE 100

BIOGRAPHICALSKETCH YuchenXiewasborninDatong,Shanxi,P.R.China.HereceivedhisBachelorofSciencedegreefromSchoolofMathematicalSciencesatPekingUniversityin2006.HethenjoinedthePh.D.programoftheDepartmentofComputerandInformationScienceandEngineeringattheUniversityofFlorida.Hisresearchinterestsincludecomputervision,machinelearning,medicalimagingandnumericaloptimization. 100