It is useful in pattern recognition particle physics computer graphics medical imaging of organs and statistical error analysis We describe an algorithm for 64257nding an equation for a multidimensional ellipsoid 64257t to a distribution of discrete ID: 82150
Download Pdf The PPT/PDF document "Preprint SMUHEP Multidimensional Ellipso..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
2fz:zTCz+bTz=1;z2dg;(4)whereC=C 1�bTCbandbT=�2bTC 1�bTCb:(5)Notethatthetermlinearinzdisappearsforanellipsoidcenteredattheorigin.Theequalitysignin(3)emphasizesthatweareonlydealingwiththesurfaceofanellipsoid.Ifitwerereplacedbyalessthanorequaltosign,wewouldhavetheequationforalledellipsoid.Further,weknowthatacenteredellipsoidalignedwiththecoordinateaxesobeystheformuladXi=1z2i R2i=1;(6)whereRiarethesemi-axislengthsoftheellipsoid.Byanalogytothis,wecanseethattheeigenvectorsofCaretheaxesoftheellipsoidandifciaretheeigenvaluesofC,thenthesemi-axislengthsoftheellipsoidaregivenby1=p ci.III.TWO-DIMENSIONALELLIPSEFITTINGOurgoalistotaketwo-dimensionalprojectionsofdatapointsonad-dimensionalellipsoidandthentoreconstructthefullellipsoidusingtheprojections,weneedamethodtondtheellipticbound-aryoftheprojectedpoints.Theboundaryofasetofdatapointscanbefoundbyusingwell-establishedconvexhullalgorithms,suchastheonein[4]orthoseintheFastGEOlibrary(availableathttp://www.partow.net/projects/fastgeo/index.html).Oncetheboundaryisfound,itcanbettoatwo-dimensionalellipse.ThealgorithmweuseistheleastsquaresminimizationmethoddescribedbyFitzgibbonetal.in[1].Theystartbydeningthe\algebraicdistance"ofapoint(x;y)tothesurfaceofageneralconicbyF(a;x)=ax=ax2+bxy+cy2+dx+ey+f;(7)wherea=(abcdef)T,andx=(x2xyy2xy1)T.TheconicisdenedbyF(a;x)=0,andsoaconicmaybettoasetofpointsthroughaleastsquaresminimizationofalgebraicdistances.Inordertoensurethatthebesttconicfoundisanellipse,oneonlyhastoimposetheconstraintthatb2�4ac0.Thiscanbedonebynoticingthat(7)allowsforarbitraryrescalingofallthecoecients.Henceonecaninsteadimposetheconstraint4ac�b2=1,whichcanbeexpressedasaTCa=1,whereC=0BBBBBB@0020000�100002000000000000000000000001CCCCCCA:(8)IntroducingamatrixDdenedbyD=(x1x2:::xN),whereNisthetotalnumberofpoints,theproblembecomesanoptimizationproblemofminimizingPNi=1F(a;xi)2=jjDajj2,subjecttotheconstraintaTCa=1. 3ThiscanbesolvedusingaLagrangemultiplierapproach.ThusintroducingaLagrangemultiplier,,weneedtominimizethequantityjjDajj2�aTCa:(9)Di erentiating(9)withrespecttoa,weseethatwejustneedtondthevectorathatsatisesSa=Ca;aTCa=1;(10)whereS=DTD.TheequationforthesolutionellipseisthengivenbyF(a;x)=0(wherexisnotapointontheellipse,butissimplythecoordinatevectorgivenbelow(7)).Thesolutionaisoneofsixgeneralizedeigenvectorsobtainedfromsolvingtherstequationin(10).In[1],Fitzgibbonetal.furthershowthatofthesixeigenvectors,onlyonewillbeassociatedwithapositiveeigenvalue,andthatthiseigenvectoristhedesiredsolution.HoweveritisimportanttonoticethatwhenalloftheNdatapointslieontheellipse,Sa=0andsothesolutionaisassociatedwithazeroeigenvalueandnotapositiveeigenvalue.Infact,aspointedoutin[5],duetonumericalerror,itispossiblethattheeigenvaluethatlabelsthesolutioneigenvectorisevenasmallnegativenumber.Appropriatenumericalprecautionsshouldbetaken.Onecanhandlethissituationbyseparatingitintothreecases:#1rank(S)=6,#2rank(S)=5,and#3rank(S)5.(Ofcourseheretheequalityandinequalitysymbolsarenotstrictandareusedlooselytomean`withinnumericalerror.')Therstcaseistheidealoneinwhichonemaysimplyndthegeneralizedeigenvaluesof(10)andpicktheeigenvectorassociatedwiththepositiveeigenvalue.Inthesecondcase,eitheronlyvepointsweregiven,orallinputpointslieontheellipse,andnumericalroundingmaymakeitdiculttochoosethecorrecteigenvectoreverytime.SinceSissymmetric,onecaninsteadcalculatethesingularvaluedecompostion(SVD)ofS,S=OTOwhereOisorthogonalandisadiagonalmatrixwiththesingularvaluesofSalongitsdiagonal.Sincerank(S)=5,oneofthediagonalentriesofiszero{withoutlossofgenerality,saythat66=0.AlsosinceOSa=OSOTOa=Oau=0;(11)u=(00000u6)T,andu6maybefoundexplicitlyfromaTCa=1.Also,a=OTOa=OTu;(12)andthesolutionvector,a,canbetakentoasequaltothesixthcolumnofOT(sincethesolutionmaybearbitrarilyscaled).ThismethodispreferredsincendingtheSVDofSisbothcomputationallylessexpensiveandlessambiguousthanndingtheeigenvaluesandeigenvectorsofS.Thethirdcaseofrank(S)5occurseitherwhentoofewpointshavebeenspeciedorwhenthedatapointsarearrangedinsuchawaythatnumerically,rank((S))5.Inthiscasethereisnouniquesolutionto(10).Notethatatleastvepointsareneededtouniquelydetermineatwo-dimensionalellipse.Ingeneral,(4)showsthatforanellipsoidinddimensions,d(d+3)=2pointsareneededtouniquelydeterminetheellipsoid.IV.FITTINGAMULTI-DIMENSIONALELLIPSOIDWITHAKNOWNCENTERForthispart,thealgorithmistheonedevelopedbyKarlin[2].ThegoalofKarl'salgorithmistotad-dimensionalellipsoidtod-dimensionaldatapointsusinglowerdimensionalprojectionsoftheellipsoid,providedthatthecenterisknown.Ifthecenteroftheellipsoidisknown,itiseasytotranslatetheellipsoidtotheoriginwherethemathiscleaner,soallcalculationsarecarriedoutasthoughtheellipsoidwerecenteredattheorigin.Inthiscase,theellipsoidcanbecompletelyrepresentedbyasymmetric,PSDmatrix. 4Thesetupoftheproblemisasfollows.LetXbeadxdsymmetric,PSDmatrixrepresentingthed-dimensionalellipsoidandletPbeadxmnon-degenerateprojectionmatrixwhosemcolumnsformanorthonormalbasisforanm-dimensionalsubspaceofthed-dimensionalspace.Themxmsymmetric,PSDmatrix,Y,representingthem-dimensionalprojectionoftheellipsoidisthengivenbyY=PTXP:(13)Hencethereisanice,linearrelationshipbetweenthePSDrepresentationofthed-dimensionalellipsoid,X,anditsm-dimensionalprojection,Y.Evenso,Karlintroducesamoreconvenientandnon-redundantwaytoexpressthisrelationship,exploitingthesymmetryoftheXandYmatrices.Karlintroducesasymmetricmatrixinnerproductspace,inwhichd-dimensionalsymmetricmatricesarerepresentedbyad(d+1)=2-dimensionalvectors,andtheinnerproductbetweentwosymmetricmatricesAandBisgivenbyhA;Bi=tr(ATB).LetxandybethesymmetricspacevectorrepresentationsoftheellipsoidsintermsofthesymmetricmatricesXandYintheorthonormalbasesM(d)jandM(m)i.Then(13)isexpressedsimplyasy=~Px;(14)wherethecomponentsofx,y,and~Paregivenbyxj=hX;M(d)jiyi=hY;M(m)ii~Pij=hM(m)i;PTM(d)jPi(15)Asimpleformtousefortheorthonormalbasisisthesetofmatriceswhichcontainallzerosexceptforeithera1inasinglelocationalongthediagonalora1=p 2intwolocations,symmetricaboutthediagonal(e.g.M(d)kl=M(d)lk=1=p 2for1l,kdandl6=k).Ifthereareqellipsoidprojections,thenthesolutionvectorxisobtainedbystackingindividualprojectionequationsasin(14)andsolvingthemsimultaneously:0BBBB@y1 y2 ... yq1CCCCA=0BBBB@~P1 ~P2 ... ~Pq1CCCCAx:(16)Forsimplicitywerewritethisasystack=~Pstackx:(17)Auniquesolutionforxexistsifandonlyifthematrix~Pstackhasfullcolumnrank.Inthiscase,thesolutionisgivenbyx=~P+stackystack=~PTstack~Pstack�1~PTstackystack;(18)where~P+stackisthefullcolumnrankMoore-Penroseinverseof~Pstack.TheMoore-Penroseinverseisapseu-doinversewhichprovidestheleastsquaressolutiontoanover-orunder-determinedproblem.NotethatthenalsolutionforxdoesnotnecessarilyrepresentaPSDmatrix.Hencethenalsolutionisnotguarenteedtobeanellipsoid,andadditionalconstraintsmustbeemployedtoattainthis.KarlsuggestsmethodsforaddingaPSDconstrainttothenalsolutionin[2].Ananalyticalformforcomputingthenearest(intheFrobeniusnorm)symmetricPSDmatrixtoanarbitraryrealmatrixisalsogivenin[6]. 5V.FINDINGTHECENTEROFAD-DIMENSIONALELLIPSOIDThecenterofad-dimensionalellipsoidmaybefoundinamanneralongthelinesofthatintroducedintheprevioussection.Recallingtheequationgivenforanellipsoidwithanarbitrarycentergivenby(4),considertheequivalentstackedmatrixequation z p z!T C 0 0 Diag[b]! z p z!=1;(19)whereextraminussignsobscuredbytakingthesquareroothavebeenabsorbedintoDiag[b].FindingCandbisthenpossiblebyusinganadditionallystackedversionofthealgorithmdescribedintheprevioussection.Letcdenotethevectordescribingthecenterofthed-dimensionalellipsoid.ThecenteroftheellipsoidcanthenbeobtainedfromCandbbysolvingthefollowingmatrixequation:c=(�2C)�1b:(20)VI.THECODEThecode,writtentotad-dimensionalellipsoidtod-dimensionaldatapoints,takestwo-dimensionalpro-jectionsofthedata,tsellipsestothetwo-dimensionalprojectionsusingthemethodgiveninSectionIII,andreconstructsthed-dimensionalellipsoidfromthetwo-dimensionalprojectionsusingthemethodgiveninSectionIV.Inparticular,allpossibletwo-dimensionalprojectionsoftheellipsoidontoplanesdenedbytwoofthecoordi-nateaxesareused.Alld(d�1)=2projectionsareneededtoensurethat~Phasfullcolumnrank.Additionally,whenusingtheseprojections,thecalculationofthe~Piissimpliedtothepointwheretheycanbereducedtoafewif-thenstatementsratherthanhavingthecodecreatethebasesanddoallofthematrixmultiplications.VII.USAGEINFORTRANTheellipsoidttingprogramusestheLAPACK(withBLAS)linearalgebralibrary(availableathttp://www.netlib.org/lapack/)whichneedstobeinstalledinorderforittowork.TheprogramcanbecalledasaFortransubroutine:subroutineEllipsoid t(inputdata,ndim,maxpoints,needcenter,printinfo,FinalMatrix,eigenvecs,eigenvals,center).Inputparameters :inputdataisadoubleprecisionarray;inputdata(ndim,maxpoints).Itcontainsall(maxpoints)ofthendim-dimensionalpoints.Therstindex(orrow)ofthearrayinputdatalabelstheaxis,i.e.theithcoordinateofthendim-dimensionalpoint.ndimisanintegerparameterwhichspeciesthenumberofdimensionsoftheto-be-reconstructedellipsoid.maxpointsisanintegerparameterwhichisequaltothenumberofdatapointsinputted.needcenterisaninteger.{needcenter=0forcestheellipsetohaveacenterattheorigin,and{needcenter=1tellsthesubroutinetondthecenteroftheellipse.printinfoisaninteger. 6{printinfo=0printsnothing,and{printinfo=1printsoutthesingularvaluesandrankofthematrixS,whichisusedtotthetwo-dimensionalellipseprojections.Italsoprintsouttheleastsquareserrorineachtwo-dimensionalt.{printinfo=2printseverythingthatprintinfo=1prints.Italsoprintsouttheellipseaxesandthecorrespondingaxislengths.Outputparameters :FinalMatrixisadoubleprecisionarray;FinalMatrix(ndim,ndim)whichcontainsthepositivesemi-denite(PSD)matrixrepresentationofthecenteredellipse.eigenvecsisadoubleprecisionarray;eigenves(ndim,ndim)whichcontainsthendim-dimensionalaxesofthereconstructedellipsesidebysideinitscolumns(secondindex).ThesearetheeigenvectorsofFinalMatrix.eigenvalsisadoubleprecisionarray;eigenvals(ndim)whichcontainsthendimaxislengthsoftheellipse.Theyarelistedinthesameorderastheaxesoftheellipseinthecolumnsofeigenvecs.ThesearetheeigenvectorsofFinalMatrix.centerisadoubleprecisionarray;center(ndim)whichisavectorrepresentingthecenteroftheellipse.VIII.CODEPERFORMANCEThecodeisfastandingeneral,itperformswell.SeveralexamplesareshownbelowinFig.1. FIG.1:(a)showsatwo-dimensionalellipsettoslightlynoisyellipticdata,(b)showsatwo-dimensionalellipsettonoisylineardata,and(c)showsatwo-dimensionalprojectionofathree-dimensionalellipsoidttonoisydata.Theellipsein(c)doesnotlooklikethebesttellipsesinceitisjustaprojectionofabesttellipsoidandinthiscase,theellipsoidtwasforcedtobecenteredattheorigin.Therearesomecaveatstothisapproach.Themethodofndingthecenteroftheellipsoidisunstableanddoesnotworkverywellwithnoisyinputdatapoints.Additionally,thoughthecodeguaranteesthatforeachtwo-dimensionalprojection,thedatawillbettoanellipse,thereis(currently)nopositivesemi-denite(PSD)constraintonthenalmatrixrepresentingthed-dimensionalellipsoid,hencethenalsolutionisnotguaranteedtobeanellipsoid.Yetusingsimulateddatawithaknowncenter,analellipsoidisalmostalwaysobtained.See[2]foraquantitativedescriptionoftheamountofnoiseintheinputdatawhichwillstillguaranteeasolutionwhichisPSD.AddinganalPSDconstraintandimprovingthemethodofndingtheequationofanon-centered,d-dimensionalellipsoidaretwopathsforfurtherworkonthisprogram. 7AcknowledgementsThisworkwasdoneasapartofa2010summerresearchprojectinthePhysicsDepartmentatSouthernMethodistUniversityinDallas,TX.ItwassupportedbyP.Nadolsky'sstart-upgrant.TheauthorthanksP.Nadolskyforhisguidanceinthiswork. [1]A.Fitzgibbon,M.Pilu,andR.Fisher,Directleastsquarettingofellipses.IEEETrans.PAMI,21(5),1999[2]W.C.Karl,Reconstructingobjectsfromprojections.PhDThesis,MIT,Dept.ofElectricalEngineeringandComputerScience,1991.[3]H.L.Lai,J.Huston,Z.Li,P.Nadolsky,J.Pumplin,D.Stump,andC.-P.Yuan,UncertaintyinducedbyQCDcouplinginCTEQ-TEAglobalanalysisofpartondistributions,arXiv:1004.4624v2[hep-ph].[4]C.Barber,D.Dobkin,andH.Huhdanpaa.TheQuickhullalgorithmforconvexhulls.ACMTrans.onMathematicalSoftware,22:469-483,1997.[5]R.Halr,J.Flusser,Numericallystabledirectleastsquaresttingofellipses.InSkala,V.,ed.:Proc.Int.Conf.inCentralEuropeonComputerGraphics,VisualizationandInteractiveDigitalMedia(WSCG98),1998.[6]N.J.Higham,Computinganearestsymmetricpositivesemi-denitematrix,LinearAlgebraAppl.,103,pp.103118,1988.