Citation
Geometric Approaches to Group-Wise Image Analysis

Material Information

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

Thesis/Dissertation Information

Degree:
Doctorate ( Ph.D.)
Degree Grantor:
University of Florida
Degree Disciplines:
Computer Engineering
Computer and Information Science and Engineering
Committee Chair:
Vemuri, Baba C
Committee Members:
Rangarajan, Anand
Ho, Jeffrey
Banerjee, Arunava
Shabanov, Sergei
Graduation Date:
8/11/2012

Subjects

Subjects / Keywords:
Algorithms ( jstor )
Computer vision ( jstor )
Datasets ( jstor )
Distance functions ( jstor )
Factorization ( jstor )
Image analysis ( jstor )
Image rotation ( jstor )
Matrices ( jstor )
Riemann manifold ( jstor )
Tensors ( jstor )
Computer and Information Science and Engineering -- Dissertations, Academic -- UF
image
Genre:
bibliography ( marcgt )
theses ( marcgt )
government publication (state, provincial, terriorial, dependent) ( marcgt )
born-digital ( sobekcm )
Electronic Thesis or Dissertation
Computer Engineering thesis, Ph.D.

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. ( en )
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.
Thesis:
Thesis (Ph.D.)--University of Florida, 2012.
Local:
Adviser: Vemuri, Baba C.
Statement of Responsibility:
by Yuchen Xie.

Record Information

Source Institution:
UFRGP
Rights Management:
Copyright Xie, Yuchen. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
857912866 ( OCLC )
Classification:
LD1780 2012 ( lcc )

Downloads

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