dtdtApathrisageodesicofasurfaceSifLsrislocallyminimized DateNovember162007 2JBAEKADEOPURKARANDKREDFIELD2MidpointMethodInthissectionourgoalisto ndanintuitivemethodforiterativelycomputing ID: 95041
Download Pdf The PPT/PDF document "FINDINGGEODESICSONSURFACESTEAM2:JONGMINB..." 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.
FINDINGGEODESICSONSURFACESTEAM2:JONGMINBAEK,ANANDDEOPURKAR,ANDKATHERINEREDFIELDAbstract.Wediscussdi erentmethodstocomputegeodesicpathsonsurfacesusingmethodsfromgraphtheoryandnumer-icalanalysis.1.IntroductionThispaperfocusesontheproblemofcomputinggeodesicsonsmoothsurfaces.Inthe rstfewsectionsweassumewearegivenanapprox-imatepathtostartfromwhenattemptingtocomputeageodesicbe-tweentwopoints.Insection2weattempttocomputethegeodesicbetweentwopointsiterativelyusingthemidpointsofanapproximatepathbetweenthem.Insection3weexploreasimilarmethod,gradi-entdescent,toiterativelyupdatethepathapproximatingthegeodesic.Insection4wecomputegeodesicsbynumericallysolvingthesystemofdi erentialequationsgoverningthem.Insection5wemodelthesurfaceasa nitegraphandapplygraphsearchmethodsto ndanapproximategeodesic.Insection6wetestthemethodsondi erentsurfacesandcomparetheresults.Thefollowingterminologywillbeusedregularlythroughoutthispaper.Inthesede nitions,Ireferstotheclosedinterval[0;1].AsurfaceisasmoothfunctionfromanopensetUR2toR3whereI2U.Apath\risapiecewisesmoothmapde nedonI.Wewillsometimerefertoasequenceofpointsasapath,thisreferstothepathconstructedbyconnectingconsecutivepoints.If\risapathfromItoI2andSisasurface,thenS\risapathfromItoR3.Thearclengthofthepathde nedbyS\risdenotedLs(\r),whereLs(\r)=10d(S\r)(t) dtdt.Apath\risageodesicofasurfaceSifLs(\r)islocallyminimized. Date:November16,2007. 2J.BAEK,A.DEOPURKAR,ANDK.REDFIELD2.MidpointMethodInthissectionourgoalisto ndanintuitivemethodforiterativelycomputingageodesicfromapointptoanotherpointqonasurfaceSwhenwearegivenanapproximatepath\rtostartfrom.WestartbyconsideringthisprobleminR2inhopesofdiscoveringamethodwhichcanbeeasilygeneralizedtotheproblemofcomputingageodesiconanysurface.2.1.GeodesicsonaPlane.Tobegin,weattempttocomputeageodesicfromapointptoanotherpointq,forp;q2R2.Weknowthatthispathshouldbethestraightlinebetweenpandq.ThemidpointmethodofcomputingageodesicbetweentwodistinctpointspandqinR2canbedescribedasfollows:(1)Startwithagivenapproximation\r1=p1;:::;pnmadeupofa nitesetofpointsinR2suchthatpiandpi+1areseperatedbyasmalldistanceandtheendpointsarep1=pandpn=q.(2)Computethemidpointsofeachsegment,andconnectthemtogetanewpath.Placep1atthebeginningofthispathandpnattheendofit.Wecalltheresultingpath\r2.(3)Iterativelyperformstep2toeachnew\r,untilyouhaveapathcloseenoughtothegeodesic.Lemma2.1.Let\rbeapathoflengthl1de nedbya niteseriesofpointsinR2thathasdistinctendpointspandq.Letlengthl2ofthepathproducedbyaniterationofthemidpointmethodappliedto\r.Thenl2l1.Thislemmafollowsdirectlyfromthetriangleinequality.Theorem2.2.Let\rbeapathde nedbya niteseriesofpointsinR2thathasendpointspandq.Performingthemidpointmethodonthispathwillproduceaseriesofpointsde ningthestraightlineconnectingpandq.Proof.Withoutlossofgeneralityconsiderthepathde nedbythe nitesetofpointsp1;:::;pnwherepi=(xi;yi)andinparticularp1=(0;0)andpn=(a;0).Thepathresultingfrommiterationsofthemidpointmethodismadeupofthepoints(xi;y0i)wherey0i= min(m;n )j=0y+j(mj) 2m.Becausethegeodesicconnectingp1andpnshouldbeasegmentofthex-axis,wesaythattheerrorofthisapproximationisthemagnitudeoftheyvaluefurthestfrom0.Sosince, FINDINGGEODESICSONSURFACES3limm !1Pmin(m;n i)j=0yi+j mj 2m=0weknowthatbyperformingthemidpointmethodonourpath,asthenumberofrecursivecallsincreasesthepathreturnedapproximatesthestraightlineconnecting(0;0)and(a;0)withanerrorthatgoesto0.Inparticular,wenotethat: Pmin(m;n i)j=0yi+j mj 2m Pmin(m;n i)j=0 yi+j mj 2mmin(m;n i)Xj=0yi+jj mm=2 2mnXj=1yjj mm=2 2mSosincePnj=1yjisjustaconstant,themidpointmethod'sx-thiterationproducesanapproximationofthegeodesicwhoseerrordi-minishesas( xx=2=2x),whichgoesto0. 0 5 10 15 20 25 0 5 10 15 20 25 (a)InitialPath 0 5 10 15 20 25 0 5 10 15 20 25 (b)Iteration100Figure1.MidpointMethodonplane.2.2.GeneralizingtheMidpointMethodtoR3.Whilethismid-pointmethodworkswellinR2,itdoesnotworkonacurvedsurfaceSbecausethemidpointmofanytwoconsecutivepointsonourpathisnotguaranteedtobeonS.Toremedythisproblem,wecalculatethepointm0onSclosesttominspaceandusethatpointinourpathapproximationinsteadofm.Forsimplesurfacessuchasspheresandtoruses,m0canbecalculatedeasily.Onarandomsurface,m0 4J.BAEK,A.DEOPURKAR,ANDK.REDFIELDcanbemuchmoredicultto ndsoweusethebuilt-inmalabfunc-tionfminsearchtocomputeit.Nowwecaniteratewiththemidpointmethodasbeforetocomptethegeodesiconanarbitrarysurface.Basedonourtestcases,thismethodcontinuestobee ective. (a)InitialPath (b)Iteration100 (c)Iteration200Figure2.MidpointMethodonsphere.Considerthecaseofasphere.Thegeodesicbetweentwopointsonasphereisalwaysasegmentofagreatcircle.Recallthatagreatcircleonasphereisonewhosediameterpassesthroughthecenterofthesphere.Inthescopeofourtestingonspheres,themidpointmethodconsistentlyproducedpathswhichapproachedthegeodesic;seeFigure2foranexample.After200iterations,thepathproducedbythemidpointmethodappearsmuchclosertothegeodesicthantheinitialpathprovided. (a)InitialPath (b)Iteration100Figure3.MidpointMethodontorus.Itisimportanttonotethatageodesicisnotnecessarilytheglobalshortestpathbetweenpandq.TakeasanexamplethetorusinFigure3.Whilethepathouralgorithmproducedapproachsageodesic,itisclearlynottheglobalshortestpathbetweenpandqoveroursurface.Becausetheoriginalpathwrappedaroundthetorus,sodoesour nal FINDINGGEODESICSONSURFACES5path.Thismethodhasthesamee ectasloopingastringaroundthetorusfollowingtheoriginalpathgiven,andthenpullingthestringtaut.Thestringwill ndtheshortestpathinaneighborhoodoftheoriginalpath,butthestructureofthesurfacepreventsitfrommovingtotheglobalshortestpathfromitsstartingpoint.Overall,themidpointmethodgeneralizesratherintuitivelytobeusedincalculatinggeodesicsovercurvedsurfaces,andwilloften ndagoodapproximationofageodesicafterarelativelysmallnumberofiterations.3.GradientDescentInthissectionweagainattempttoiterativelycomputeageodesicfromanapproximatepathusingamethodsimilartothatdescribedinSection2.2.Now,insteadofexaminingthemidpointofeachsegmentofthepath,weexaminethemiddlepointofeachconsecutivesequenceofthreepointsde nedbythepathandattempttoimprovethepath'sapproximationofthegeodesicbyalteringpointsonebyone.Letusconsiderasimplemodelofiterativere nementofapathde nedbythreeconsecutivepoints.Wecallthissequencepi;pi+1;pi+2inI2.Thissection'smethodinvolves rst xingpiandpi+2,andthenselectinganewmiddlepointufromaneighborhoodoftheoriginalpi+1whichallowsthepathtobetterapproximatethegeodesicconnectingpiandpi+2.ThelowerboundonthearclengthofthegeodesicconnectingS(pi)andS(pi+2)whereSisasurfaceistheEuclideandistancekS(pi) S(pi+2)k.Assumingthatthispathmustalsocontainpi+1,thelowerboundbecomesLB(pi+1)whereLB(x)=kS(pi) S(x)k+kS(x) S(pi+2)k:IfthesurfaceSislocallyplanar,andthepointsinthesequencearecloseenough,thislowerboundwillactuallyapproximatetheminimalgeodesiclength.Therefore,itwouldgiveusabetterapproximationoftheminimalgeodesiclengthifwedecreasedthislowerboundonourarclengthbychoosinguthatminimizesLB(u).TominimizeLB(u)withrespecttou,weemploygradientdescent;thatis,aniterativemethodinwhichuisupdatedinthedirectionopposingthegradientLBatthecurrentvalueofu,usingsomeweight �0:u u LB(u):Inourimplementation, ischosensuchthatitminimizesthefol-lowingsum: 6J.BAEK,A.DEOPURKAR,ANDK.REDFIELDkS(u LB(u)) S(pi)k+kS(u LB(u)) S(pi+2)k:Toextendthismethodfromourshortpathtoapathofarbitrarylength,weconsiderthreeconsecutivepointsatatimeandupdatethemiddlepointofeachofthesesubsequencespi;pi+1;pi+2ofthepathaswego. (a)InitialPathonttplane (b)InitialPathonttsphere (c)InitialPathontttorus (d)Iteration500 (e)Iteration500 (f)Iteration500Figure4.GradientDescentonplane,sphereandtorus.Weempiricallyobservedthatthegradientdescentise ectiveonlywhenkS(pi) S(pi+1)kandkS(pi+1) S(pi+2)karecomparable.Oth-erwise,theupdateisveryslow.Hence,itisbene cialtodroppi+1fromthesequenceifonelegbecomestoosmall,andinsteadaddapointtothepartofthesequencethatcontainsalongleg.InFigure4weseetheresultsofusingthegradientdescentmethodonarandompathoveraplane,sphereandtorusafter500iterations.Ingeneral,whileeachindividualiterationofthegradientdescentmethodismuchfasterthanthemidpointmethod,itoveralltakesmanymoreiterationsthanthemidpointmethodtocomputesimilarlyaccurateapproximationsofageodesiconanygivensurface. FINDINGGEODESICSONSURFACES74.DifferentialEquationApproachInthissectionweformulatetheproblemof ndinggeodesicsonasurfaceastheproblemofsolvingasystemofordinarydi erentialequations.AlthoughwetreatthecaseofsurfacesinR3inthispaper,wecaneasilygeneralizeittoarbitrarysubmanifoldsofRN.4.1.TheGeodesicEquation.ConsiderasmoothsurfaceMR3givenbyaparametricmapS:U!R3;whereUR2isanopenset.WeassumethatSisasmoothimmersion.Inotherwords,weassumethatdSxisinjectiveforallx2U.Apath\r:I!UgivesapathS\r:I!Monthesurface.ThelengthofthispathisgivenbyLS(\r)=10d(S\r) dt(t)dt:Bythechainrule,weobtainLS(\r)=10kdS\r(t)d\rtkdt=10q d\rTtdST\r(t)dS\r(t)d\rtdt(1)De neaninnerproductonthetangentspacetoMatp=S(u)bygu(v;w)=hdSuv;dSuwi=vTdSTudSuwUsingthisinnerproduct(1)canberewrittenasLS(\r)=10q g\r(t)(_\r(t);_\r(t))dt:Foru2U,weletubethe2by2matrixd Txd xcorrespondingtotheinnerproductgu.Observethat,uisinvertiblesinceitcorrespondstoapositivede nite(andhencenondegenerate)innerproduct.Denotetheentriesofbygijandtheentriesof 1bygij.For1i;j;k2,de netheChristo elsymbols kijasfollows(2) kij=1 2nXl=1gkl@gjl @xi+@gli @xj @gij @xl:Thesystemofdi erentialequationssatis edbygeodesicpathscanbegivenintermsoftheChristo elsymbols. 8J.BAEK,A.DEOPURKAR,ANDK.REDFIELDTheorem4.1.Let\r:I!Ubeasmoothmap.Denotethecoordi-natesof\rbyx1andx2.ThenthepathS\risageodesiconMif\risatisfythedi erentialequations(3)xk+Xi;j kij_xi_xj=0:Example4.2.LetuscalculatetheChristo elsymbolsandtheequa-tions(3)forasphere.RecallthattheunitsphereinR3canbegivenparatemricallyusingthesphericalcoordinatesby :(0;2)(0;)!R3 :(; )7!(cossin ;sinsin ;cos )ThentheJacobianmatrixd andtheinnerproductmatrixaregivenbyd =0@ sinsin coscos cossin sincos 0sin 1A=sin2 001TheChristo elsymbolsare 1=0cot cot 0 2= sin cos 000Hence,apath\r:I!(0;2)(0;)representsageodesiconthesphereifthecoordinates(; )of\rsatisfythedi erentialequationsgivenby+2__ cot =0 _2sin cos =0:Havingobtainedthedi erentialequationforgeodesics,weturntothenumericalsolutionofthedi erentialequation.Assumethesetupdescribedatthebeginningofthissection.GiventwopointspandqinU,thetaskisto ndageodesiconMjoiningS(p)toS(q).Equiva-lently,wewantto ndapath\r:I!Uwhosecoordinatessatisfythedi erentialequations(3).Namely,for1k2wemusthave,xk+Xi;j kij_xi_xj=0:Wesolvethesystemusingthemethodof nitedi erences.WechoosealargepositiveintegerNandsubdividetheintervalIintoNequalpartsinordertoreplacederivativesby nitedi erences.Fornotationalsimplicity,de ne=1=N.Furthermore,letxp=\r(p)andxpjbethe FINDINGGEODESICSONSURFACES9jthcoordinateofxpfor1j2and1pN.Forbrevity,weletxpk=xp+1k xp 1k.Observethat_xk(p)=1 2xpk+(2);xk(p)=1 2(xp+1k+xp 1k 2xpk)+(2)Weapproximatethe rstderivativeas_xk(p)xpk=2andthesecondderivativeasxk(p)(xp+1k+xp 1k 2xpk)=2.Substitutingin(3),weobtainthesystemxpk=xp+1k+xp 1k 2+1 4Xk;j ki;j(xp)xpixpj(4)subjecttotheboundaryconditionsx0=pandxN=q.Wediscusstwoiterativemethodstosolve(4):oneusingfunctionaliterationandoneusingNewton-Raphsonmethod.Weapproximateapath\r:I!UbytheN+1tuple(\r(0);:::;\r(i);:::;\r(1))ofpointsinU.Inbothmethods,weassumethatwearegivenaninitialpathfromptoqrepresentedbytheN+1tuples=(s0;:::;sN);whereeachspisapointinUfor1pN.Startingwiththisinitialpath,we ndasequenceofpathsinU(representedasN+1tuplesofpointsinU)thatconvergestoasolutionof(4).4.2.FunctionalIteration.Wephrasetheproblemofsolving(4)astheproblemof ndinga xedpointofamap:UN+1!R2(N+1).WedenoteelementsofR2(N+1)by(x0;:::;xN)wherexp=[xp1;xp2]Tfor0pN,wherexp1;xp22R.De nethemapby(x0;:::;xN)=(y0;:::;yN);wherey0=x0;yN=xNypk=xp+1+xp 1 2+1 4X1i;j2 ki;j(xp)(xp+1i xp 1i)(xp+1j xp 1j);for1k2and1pN 1(5)Wecan nda xedpointofbyiteratingstartingwiththegivenpointsofR2(N+1).Inotherwords,wede neasequencefx(i)gi2NofelementsofR2(N+1)byx(0)=sx(i+1)=(x(i))fori1. 10J.BAEK,A.DEOPURKAR,ANDK.REDFIELDWeconjecturethatthesequenceiswellde nedandconvergentifmaxi;psp+1i spiissucientlysmallandsissucientlyclosetoa xedpointof.Inotherwords,thesequenceconvergesifsucces-sivepointsintheinitialsequence(s0;:::;sN)arecloseenoughandsapproximatesapathwhichissucientlyclosetoageodesic.Ifthesequencex(i)converges,itisclearthatthelimitingvalueisa xedpointof.Furthermore,observethatx(i)0=s0andx(i)N=sNforalli2N.Hence,thelimitingvalueapproximatesapathinUfromptoq,anditsimageunderSapproximatesageodesicjoining (p)and (q).Observethatforaplane,alltheChristo elsymbols ki;jarezero.Hence,thefunctionaliterationmethodgivesthemidpointmethodof ndinggeodesics.AssumingthatapplyingthegivenparametrizationfunctionStakesconstanttime,letuscalculatetherunningtimeofeachiterationofthefunctionaliterationmethod.Usingthenotationof(5),weseethateachiterationinvolvescalculatingnewypkfor1k2and1pN.Foragivenkandp,calculatingnewypkinvolvescalculatingtheChristo elsymbolsatxp.Thiscanbedoneinconstanttimebycalculatingtheentriesgij(xp)ofxp,theentriesof 1xp,thederivatives@gij @xlandusing(2).Derivativescanbecomputednumerically,orsymbolicallyiftheparametrizationSadmitssymbolictreatment.Thus,eachypkcanbecomputedinconstanttime,leadingto(N)timeforeachiteration.Figure6andFigure5showthepathsobtainedbyapplyingthefunctioniterationmethodtoaninitialpathonaplaneandatorus.4.3.Newton-RaphsonMethod.Wecanrephrasethequestionof ndinga xedpointofasthequestionof ndingazeroof(x) x.Wecanthen ndazerobyNewton-Raphsonmethod,startingwiththeinitialguesss.However,thefunction(x) xhassingularJacobian;hencewehavetomakeaslightadjustmenttotousetheNewton-Raphsonalgorithm.Recallthatinthesequencex(i)=(x0(i);:::;xN(i))de nedabove,theextremevaluesx0(i)andxN(i)stayconstantats0andsNrespec-tively.Hence,wecanrestrictouriterationtotheintermediatevalues.Inotherwords, xx0=s0,xN=sNandde neH:UN 1!R2(N 1)byH(x1;:::;xN 1)=(y1;:::;yN 1); FINDINGGEODESICSONSURFACES11 (a)Initialpath (b)Iteration25 (c)Iteration100 (d)Iteration250Figure5.FunctionalIterationonaToruswhereypk=xp+1+xp 1 2+1 4X1i;j2 ki;j(xp)(xp+1i xp 1i)(xp+1j xp 1j);for1k2and1pN 1Thus,if(x1;:::;xN 1)isazeroofH(x) x,then(x0;x1;:::;xN 1;xN)isasolutionof(4).WeiterateusingNewton-Raphsonalgorithmto ndazeroofH(x) x,startingwiththeinitialguess(s1;:::;sN 1).Inotherwords,wede nethesequencefy(i)gofpointsofUbyy(0)=(s1;:::;sN 1)y(i+1)=y(i) (dHy(i) I) 1(H(y(i)) y(i));fori1:(6)Weconjecturethatifmaxi;psp+1i spiissucientlysmallandsisclosetoasolutionof(4)thenthesequencede nedby(6)iswellde nedandconvergestoasolutionof(4).Aswedidforthefunctioniterationmethod,letuscomputethenumberofstepsrequiredforeachiterationoftheNewton-Raphsonmethod.Observethateachiterationy(i)7!y(i+1)involvesthecomputationofH(y(i)),whichcanbedonein(n)timeasdescribed 12J.BAEK,A.DEOPURKAR,ANDK.REDFIELD (a)Initialpath (b)Iteration25 (c)Iteration100 (d)Iteration250Figure6.FunctionalIterationonaPlanebefore.However,wealsoneedtocomputetheJacobiandHy(i)andtheinverseof(dHy(i) I.SincedHy(i)isa(2N 4)(2N 4)matrix,thiswouldusuallytake(N3)time.However,weonlyneedtocalculate(dHy(i) I) 1(H(y(i)) y(i)),whichisequivalenttosolvingforZthesystem(7)(dHy(i) I)Z=H(y(i)) y(i):Furthermore,observethaty(i+1)pdependsonlyony(i)p1.HencedHy(i) Ihas(N)nonzeroentries,whicharelocatedonornearthediagonal.Inotherwords,dHy(i) Iisabandedmatrix,forwhich(7)canbesolvedin(N)steps.Thus,eachiterationoftheNewton-Raphsonmethodcanbedonein(N)steps.Figure7andFigure8showtheresultsofapplyingtheNewton-Raphsonmethodtoaninitialpathonatorusandaplanerespectively.Finally,weremarkthatalthoughcalculatingtheChristo elsymbols ki;jatapointisaconstanttimeoperation,numericallycalculatingthemcanbecomputationallyquiteintensiveinpractice.Boththemethodsoutlinedinthissectionrunmuchfasterifwecancompute ki;jsymbolicallyonceandforallandevaluatethesesymbolicexpressionsatparticularpoints.Forexample,inourimplementation,oneiteration FINDINGGEODESICSONSURFACES13 (a)Initialpath (b)Iteration1 (c)Iteration2 (d)Iteration3Figure7.Newton-RaphsononaTorus (a)Initialpath (b)Iteration1Figure8.Newton-RaphsononaPlaneofNewton-Raphsonforasphereusingsymbolic ki;jtook1second,whereasoneiterationusingnumericallycomputed ki;jtook69seconds(forN=26). 14J.BAEK,A.DEOPURKAR,ANDK.REDFIELD5.GeneratingAnInitialPathUsingGraphTheorySofarwehaveconsideredvariouswaysofiterativelyre ningapathtoageodesic.Inordertoemploythesemethods,however,onerequiresanapproximateshortestpathtobeginwith.Thissectionexploresapotentialmethodforgeneratingsuchpathusinggraphtheory.Justasapathcanbeapproximatedbyasequenceofpoints,asur-facecanalsobeapproximatedasameshofregularly-spacedpoints.Forinstance,computer-generatedvisualizationsofsurfacesnecessarilyinvolveadiscreteapproximation.Ifameshcansuccessfullyapprox-imateasurface,thenwecantreattheproblemof ndingageodesicastheproblemof ndingtheshortestpathbetweentwoverticesonagraph.Aswewillshow,theshortestpathonthemeshapproximationofthesurfaceisnottoofarfromtheminimalgeodesic.5.1.Induced-MeshApproximation.Firstwereviewthede nitionofaweightedgraphandconsiderhowtoconstructaweightedgraphthatapproximatesagivensurfaceS.De nition5.1.Aweightedgraphisanorderedpair(V;E)corre-spondingtothesetofverticesandthesetofedges,alongwithaweightfunctionw:E!Rspecifyingtheweightofeachedge.AnaturalwaytocreateagraphonagivensurfaceSistopickregularly-spacedpointsinI2asvertices,andtojoineachpointtoitsfourimmediateneighbors,wherethecostofeachedgee=fv1;v2gistheEuclideandistancebetweenS(v1)andS(v2).De nition5.2.Forn2Zgreaterthan1,aninducedmeshonasurfaceSistheweightedgraph=(V;E)whereV=0;1 n;2 n;:::;12I2;E=ffv1;v2gjv1andv2areadjacent.g;withw:fv1;v2g7!kS(v1) S(v2)k;iffv1;v2g2E.HeretwoverticesareadjacentiftheyareneighborlatticepointsinI2.WedenotetheinducedmeshbyMn(S),andcallnthesizeoftheinducedmesh.OnewaytocharacterizetheinducedmeshistorealizethatEistheapproximationofthemetricinducedonI2byS.Weareestimatingtheminimalgeodesicdistancebetweenv1andv2bytheEuclideandistance,whichisreasonableaslongasSislocallyplanarandv1;v2areclose.As FINDINGGEODESICSONSURFACES15anexample,Figure9comparesaparticularsurfaceSwithaninducedmeshonS. (a)ThesurfaceS. 0 0.5 1 0 0.5 1 -1 0 1 X Y Z (b)TheinducedmeshM10(S).Figure9.ComparisonofthesurfaceS(u;v)=(u;v2;sin((u v)))withaninducedmeshonit.De nition5.3.Let=(V;E)beaweightedgraphwithwasitsweightfunction.Thenapathbetweens;d2Visthesequenceofvertices(v1;v2;:::;vk)wherev1=sandvk=dandfvi;vi+1g2Eforalli=1;:::;k 1.Thesequenceisashortestpathincaseitminimizesthesumk 1Xi=1w(vi;vi+1):Hence,theshortestpathonaninducedmeshisequivalenttothephysicallyshortestpassageonthemeshalongthegridlines.Thereareseveralknownalgorithmsthat ndtheshortestpathonaweightedgraph.Dijkstra'sAlgorithmisaveryecientsuchalgorithm.Itspseudocodeisgivenintheappendix.Theorem5.4.[CL,Theorem24.6]Let=(V;E)beaweightedgraphwithw:E!Rasitsweightfunction.GiventwopointsinVsuchthatapathexistsbetweenthetwo,Dijkstra'sAlgorithm ndstheshortestpathbetweentheminasymptoticruntimeof(V2+E).Infact,Dijkstra'sAlgorithmcanbeexecutedin((E+V)logV)whenthegraphissparse,usingabinaryheapandanadjacencymatrix.IncaseofMn(S),wehaveE=2n(n+1)andV=(n+1)2.Remark5.5.GivenasurfaceS,Dijkstra'sAlgorithmcan ndtheshortestpathonMn(S)inasymptoticruntimeof(n2logn). 16J.BAEK,A.DEOPURKAR,ANDK.REDFIELDTherealsoexistsageneralizationofDijkstra'sAlgorithm,calledA-Star,incasealowerboundonthedistancetodcanbegivenateachvertexv.InthecaseofaninducedmeshonS,theEuclideandistancekS(d) S(v)kcanserveasthelowerbound.A-StarhasthesameasymptoticruntimeasDijkstra,butoftenrunsfaster.A-starreturnsasequenceofvertices(v1;:::;vk)thatresidesinI2.TotranslatethesequenceintoapathonI2,asimplelinearinterpola-tionwillsuce,asillustratedinthefollowingde nition.De nition5.6.Let(v1;:::;vk)beasequenceofpointsinI2.Thenthelinearinterpolationofthesequenceisthepiecewisefunction\r:I!I2where\r:t7!8���������&#x-0.1;蚅怀&#x-0.1;蚅怀&#x-0.1;蚅怀&#x-0.1;蚅怀&#x-0.1;蚅怀&#x-0.1;蚅怀&#x-0.1;蚅怀&#x-0.1;蚅怀&#x-0.1;蚅怀:v1(1 tk)+v2(tk)if0t1=kv2(2 tk)+v3(tk 1)if1=kt2=k......vi(i tk)+vi+1(tk i+1)if(i 1)=kti=k......vk 1(k 1 tk)+vk(tk k+1)if(k 1)=kt1Inotherwords,welinearlyinterpolatebetweeneachpairofconsec-utivevertices.Notethatthecomposition =S\rconstitutesapathonthesurfaceS.Figure10providesanexample. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 uv (a)Apathona10-by-10latticeonI2. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 x y z (b)Thesamepathonthein-ducedmeshonS. (c)CompositionofSwiththelinearinterpolationofthepath.Figure10.TranslationofasequenceofpointsinI2intoapathonS.Therefore,thefollowingrecipeyieldsanapproximateshortestpathonagivensurfaceS.(1)GenerateaninducedmeshMn(S). FINDINGGEODESICSONSURFACES17(2)ExecuteA-startogenerateasequenceofpointsbetweentwovertices.(3)Constructalinearinterpolation\rofthesequenceinI2.(4)Return =S\r.5.2.Convergence.Thepathgeneratedbytheaboverecipebegsthequestionofaccuracy.Itdi ersfromtheminimalgeodesicsincethepathmustalwaystravelalongthegridlinesofthemesh,andeachpiecewiselegofthepathisnotitselfaminimalgeodesic.Figure11illustratesthispoint:whereasthegeodesicjoining(0;0)and(1;1)inI2hasarclength 2,alinearinterpolationofthesequencereturnedbyA-starwillhavearclengthof2.Nonetheless,wemaystillprovethatthearclengthofthispathisnottoofarfromthatoftheminimalgeodesic. 0 0.5 1 0 0.5 1 XY (a)Thegeodesicjoining(0;0)and(1;1)inI2. 0 0.5 1 0 0.5 1 XY (b)Thelinearinterpolationofoneshortestpathontheinducedmeshofsize10.Figure11.Comparisonofthegeodesicjoining(0;0)and(1;1)inR2andtheapproximationviathemethodinSection5.1.Throughoutthissection,we xtwopointsp1;p22I2,letSbeasur-face,andlet~\r:I!I2bethepathbetweenp1andp2thatminimizesLS(),parametrizedbythearclengthofitscompositionwithS.TakeaninducedmeshMn(S).IdeallywewanttocompareS~\rwith =S\r,theapproximatepathconstructedbyourrecipeinSection5.1.However,sincewedonotknowapriorithepathreturnedbyA-star,weinsteadgenerateadi erentpath asdescribedbelow:ConsiderallsquaresonthegridofMn(S)thatS~\rpassesthrough.Thepathentersthroughoneofthefouredges,andleavesthroughanother.Therearetwocases:Thesetwoedgesarethesameoradjacent. 18J.BAEK,A.DEOPURKAR,ANDK.REDFIELDThesetwoedgesareparallel.Inthiscase,wekeepattachingsquaresthroughwhichS~\rentersandleavesthroughparalleledges,untilS~\rleavesthroughanedgeadjacenttotheentranceedge.SeeFigure12(b). ab (a)Adjacenten-tranceandexitedges. ab ab (b)Non-adjacententranceandexitedges.Inthiscase,weaddsquaresoneithersideuntilthepathexitsasquarethroughanon-paralleledge.Figure12.Constructionofanapproximatepathalongthemesh.Theredcurvesindicate~\r,andthebluelinesindicateourconstruction.Inthe rstcase,jointhetwoentraceandexitpointsalongthegrid,choosingthedirection(clockwiseorcounterclockwise)tominimizethelengthinI2.Iftheentranceandexitpointslieonthesameedge,wesimplyconnectthem;otherwise,thepathwillcontainthevertexsharedbythetwoedges,asseeninFigure12(a).Inthesecondcase,theentranceandexitpointsof~\rforthekconsec-utivesquaresmayormaynotlieonthesameside.SeeFigure12(b).Iftheydo,simplyconnectthem;ifnot,makeatransitionalongoneoftheparalleledges.Nowwejoinalltheselineswehavegeneratedforallsquares~\rpassesthrough.Theselinesformapiecewise-straight,continuouspathonI2.Furthermore,itisavalidpathonMn(S),sinceitentirelyconsistsofgridlines;eachedgecanneverbepartiallycovered,becauseifitwere,thenwemustbebacktracking,andwecansimplyremovetherepetition.NowwehavesuccessfullyobtainedapathinI2alongtheinducedmesh,anddenoteitby\r.Let beaparametrizationofS\rbyitsarclength. FINDINGGEODESICSONSURFACES19WecompareS~\rwith ineachofthesquare(s)weconsideredabove.Fornow,assumethe rstcase(Figure12(a).)TheentranceandexitpointsareaandclyinginI2,withbeingthecornervertex.Weexaminethefollowingarclengths:L1,thearclengthoftheminimalgeodesicbetweenS(a)andS(c)onthesurface.L2,theEuclideandistancebetweenS(a)andS(c).L3,thesumoftheEuclideandistancesbetweenS(a)andS(),andbetweenS()andS(c).L4,thearclengthofthecompositionofSwiththelinearinter-polationoffa;b;cg.L1correspondstothearclengthoftheportionofS~\rcontainedwithinthesquare,whereasL4correspondstothearclengthoftheportionof containinedwithinthesquare.Weremindthereaderthattheselengthsarethearclengthsofthepartcontainedwithinasinglesquare,notthelengthoftheentiregeodesic.Lemma5.7.L1andL4areatmost2 nq 2x+2y+2z;wherexisthemaximumofkrSxkinI2,andy;Gzareanalogouslyde ned.Proof.TakeasimpleinterpolationbetweentheendpointsofS~\rinsidethesquareofsidelength1=nbyconnectingthemlinearlyinsideI2andsendingthislineinI2throughS.Thenthearclengthisboundedaboveby 2 nq 2x+2y+2z:Toseethis,notethattherateofincreaseinarclengthrelativetoaunitstepinI2(inanarbitrarydirection)isatmostp 2x+2y+2z.Sincewemoveatotaldistanceofatmostp 2 ninsideI2(becauseweareinsideasmallsquare,andthemostwecanmoveinasingledirectionisalongthemaindiagonal),thetotalarclengthmustbeboundedaboveassuch.BecauseL1isshorterthanthisinterpolation,thelemmaforL1follows.AsforL4,sinceittakesaright-angledturninI2,wemoveatotaldistanceofatmost2 ninsideI2,whichindicatesthatthearclengthisatmost2 np 2x+2y+2zasdesired.Lemma5.8.Let:I!S(I2)beapathonthesurfaceSparametrizedbyarclength.Fix0t1t21.Thenthearclengthofbetween 20J.BAEK,A.DEOPURKAR,ANDK.REDFIELDt=t1andt=t2di ersfromthestraightlinedistancek(t2) (t1)kbyatmost(t2 t1)2q H2x+H2y+H2z;whereHxisthemaximumofkr2SxkinI2,andHy;Hzareanalogouslyde ned.Proof.Thearclengthofbetweent=t1andt=t2is,(8)L=t2t1s dx dt2+dy dt2+dz dt2dt:Lett;x;y;zbethedi erencest2 t1;x(t2) x(t1);y(t2) y(t1);z(t2) z(t1),respectively.Then,byMeanValueTheorem,thereexistssomet02[t1;t2]suchthatdx dt(t0)=x t:Thisyields,forr2[t1;t2], dx dt(r) = dx dt(t0)+rt0d2x dt2dt dx dt(t0) +rt0 d2x dt2 dt dx dt(t0) +rt0Hxdt dx dt(t0) +tHx:Therefore,forthesamevaluesofr,dx dt(r)2 x t +tHx:2:Similarinequalitiesholdfordy dt(r)2anddz dt(r)2.PluggingthesebackintoEquation(8),weobtainLt2t124s x t +tHx2+:::35dt:Sincetheintegranddoesnotdependont,wegetLq (x+t2Hx)2+(y+t2Hy)2+(z+t2Hz)2:Bythetriangleinequality,weobtain FINDINGGEODESICSONSURFACES21Lp x2+y2+z2+t2q H2x+H2y+H2z=k(t2) (t1)k+(t2 t1)2q H2x+H2y+H2z:SincethearclengthoftheminimalgeodesicisboundedbelowbytheEuclideandistancebetweentheendpoints,namelyk(t2) (t1)k,thelemmafollows.Theorem5.9.L2=L1+1 n2;andL4=L3+1 n2:Proof.ByLemma5.7,thelengthofL1is(1=n).Then,inapplicationofLemma5.8toL1,thedi erencet2 t1mustalsobe(1=n)becauseS~\risparametrizedbyarclength.Thenitfollowsthatthedi erencebetweenL2andL1is(1=n)2q H2x+H2y+H2z=(1=n)2;asdesired.ThestatementforL4 L3followssimilarly.Lemma5.10.L3 2L2:Proof.NotethatL2takesstraightlineswhereasL3travelsalongthegrid.Inotherwords,L2isthehypotenusesofatriangleofwhichthelegsformL3.Becausetheratioofthesumofthelegstothehypotenuseisatmost 2,thelemmafollows.Sofarwehaveassumedthecaseinwhich~\rentersandleavesasquareonMn(S)throughadjacentedges.Intheothercase,showninFigure12(b),analogousresultstoLemma5.10andTheorem5.9canalsobeobtained.Weomittheirproofforbrevity,butthereadershouldbeabletoconvincehimselfoftheomittedfact.Theorem5.11.Asnapproaches1,thearclengthofthepathgener-atedbytherecipeinSection5.1convergestoatmost 2LS(~\r).Proof.CombiningTheorem5.9withLemma5.10yieldstherelationL4 2L1=1 n2:Summingoverallsquares,weobtainthatLS(\r)islessthanthesumof 2LS(~\r)withan(1=n)term,sincethenumberofsquaresthrough 22J.BAEK,A.DEOPURKAR,ANDK.REDFIELDwhichthepathpassesis(n).Thisquantityconvergestozeroasnapproaches1,asdesired.Now,sincetherecipereturnsanoptimalpathwithsmallerarclengththan ,thesameinequalitymuststillhold.Theorem5.11tellsusthattherecipeinSection5.1isareasonableinitializationfortheaforementionediterativemethods,suchasmid-pointmethod,gradientdescent,andNewton-Raphson.FromFigure11,weseethattheboundinTheorem5.11istightinthesensethatreplacing 2withasmallerconstantwillinvalidatethetheorem.Theconstant 2seemstoarisefromthestructureofthemesh,however,ratherthanthemethodsintheproofofthetheorem.Weconjecturethataugmentingthemeshstructurebyaddingdiagonalsineachsquare,etcetera,shouldreducethisconstant.6.ComparisonStudyInthissection,weapplythemethodsdiscussedintheprevioussec-tionstoobtaingeodesicpathondi erentsurfaces.Wecomparetheresultsafterdi erentnumberofiterationsandsummarizeour ndings.6.1.Methodology.WeconstructedfoursurfacesusingMatlabandselectedtwopointsoneach.TheparametrizationsofthesurfacesaregiveninTable1oftheappendix.Foreachsurface,aninitialpathwasgeneratedbyA-staronaninducedmeshofsize40;anotherinitialpathwasgeneratedbyconnectingthetwopointswithasinusoidalwaveoftwofullperiodsinI2.Foreachsurface-pathpair,threemethodswereusedtoiterativelyre nethepathintoageodesic:1)midpointmethod,2)gradientdescent,and3)Newton-Raphson.TheNewton-RaphsonmethodwasusedonlyfortheinitialpathconstructedfromA-star,asthesinusoidalpathistoofarfromgeodesicandthemethoddoesnotconverge.Figures14through16intheAppendixdisplaytheresultsonthefoursurfaces.Part(a)istheapproximatepathgeneratedbyA-star;Part(b)-(d)aretheresultingpathsaftersomenumberofiterationsofeachofthethreemethods.ThethreemethodsareabbreviatedasMP(Midpointsearch),GD(GradientDescent)andNR(Newton-Raphson),andthenumberofiterationsusedisgiveninparentheses.Part(e)-(g)showthedecreaseinarclengthoveriterationsforeachmethod.Part(h)-(l)areanalogousresultsfrominitializingwithasinusoidalpath. FINDINGGEODESICSONSURFACES236.2.Results.Ingeneral,theNewton-Raphsonmethodusingdi er-enceequationsconvergedinthefewestiterations,stabilizinginacoupleiterationsatmost.WeemphasizethattogenerateFigures13through16,theNewton-Raphsonmethodranfor5iterations,whereasmidpointsearchandgradientdescentrequired200iterationstogetreasonablyclosetotheactualgeodesic.However,eachiterationofthesemeth-odsconsumedverylittletime,whereastheNewton-Raphsonmethodsu ersfromthebottleneckofcalculatingChristo elsymbols.Incasetheparametrizationofthesurfaceisknownandtractable,onecanpre-computetheChristo elsymbolssymbolically,ratherthannumerically,toconservetime.Inthatcase,Newton-Raphsonhasacomparableruntimeperiteration.SeethediscussionattheendofSection4.TheNewton-Raphsonmethodalsorequiredtheinitialpathtobere-liableandeachconsecutivepairofpointstobeveryclose.Forinstance,theNewton-Raphsonmethodwoulddivergewheninitializedwiththesinusoidalpathsgiveninpart(h)ofFigures13through16.Ontheotherhand,theothertwonumericalmethodsdidnotrequiretheinitialpathtobeveryreliable,andfunctionedsatisfactorilywheninitializedwithsuboptimalpaths.Betweenthemidpointsearchandgradientdescent,midpointsearchperformsslightlybetter,especiallywhenthesurfaceisnotveryplanarortheinitialpathiswavy.Weremarkthatonecanusethesetwomethodsinconjunction|bypickingtheEuclideanmidpointand nd-ingtheclosestpointthatliesonthesurface,andthenimprovingthispointbyrunninggradientdescent.7.PrimaryAuthorsSection1throughSection3wereprimarilyauthoredbyKatherineRed eld;Section4wasprimarilyauthoredbyAnandDeopurkar.Sec-tions5and6,alongwiththeappendix,wereprimarilyauthoredbyJongminBaek. 24J.BAEK,A.DEOPURKAR,ANDK.REDFIELD8.Appendix Program1MidpointSearch //Iterativelyrefinesagivenpath//S:asurfaceI^2---�R^3//pts:asequenceofpointsinI^2//n:thenumberofiterationsfunctionMidpoint_Search(S,pts,n)forniterationsfori=1tolengthofpts//MapthepointsontoSpts2(i)=S(pts(i))fori=1tolengthofpts//Takemidpointsmidpts(i)=(pts2(i)+pts2(i+1))*0.5fori=1tolengthofmidpts//findtheclosestpointonthesurfacetomidpts(i)pts(i)=argmin_p||S(p)-S(midpts(i))||//argmincanbedoneinmatlabusingfminsearch. Program2GradientDescent //Iterativelyrefinesagivenpath//S:asurfaceI^2---�R^3//pts:asequenceofpointsinI^2//n:thenumberofiterationsfunctionGradient_Descent(S,pts,n)forniterationsfori=1tolengthofpts-1//Takemidpointsmidpts(i)=(pts(i)+pts(i+1))*0.5fori=1tolengthsofmidpts//Performgradientdescent,startingfrommidpts(i)f(p):=||S(pts(i)-p)||+||S(pts(i+2)-p||grad=argmax_g(df(p)/dx,df(p)/dy).*gevaluatedatp=midpts(i)beta=argmin_bf(p-b*grad)//argmincanbedoneusingfminbndinmatlabpts(i)=midpts(i)-beta*grad FINDINGGEODESICSONSURFACES25 Program3FunctionIterationandNewtonRaphson \\Computesoneiterationofthefunctionaliteration\\method\\S=MapfromI^2--�I^3\\pts=SequenceofN-2pointsinI^2\\init=Initialpoint(inI^2)\\final=Finalpoint(inI^2)functionIter(S,init,final,pts)x(0):=initx(N-1):=finalx(i):=pts(i)for0iN-1forp=1toN-2fork=1,2y(p,k)=(x(p+1,k)+x(p-1,k))/2+1/4*Sum[Gamma(k,i,j)(x(p))*(x(p+1,i)-x(p-1,i))*(x(p+1,j)-x(p-1,j))]//Sumtakenoveralli,jin{1,2}//CalculatingGammaisstraightforward.Allthe//derivativesencounteredcanbecomputed//numerically.returny\\ComputesoneiterationoftheNewtonRaphsonmethod\\S,pts,init,finalareasinIterfunctionNRIter(S,init,final,pts)\\Lookatptsasa2x(N-2)arrayy:=Iter(S,init,final,pts)dH:=JacobianofH(X)-&#x-514;&#x.492;Iter(S,init,final,X)atX=pts\\Thiscanbecomputednumerically.\\Lookingatptsas2N-4tuple,dHwouldbea\\2N-4by2N-4matrix.Z:=Inverse(dH-I)*(y-pts)\\Canbeobtainedbysolving(dH-I)Z=y-pts\\UseamethodthatexploitsthatdH-Iis\\boundedforoptimalperformance.\\LookatZasa2byN-2array.return(pts-Z) 26J.BAEK,A.DEOPURKAR,ANDK.REDFIELD Program4Dijkstra'sAlgorithm //FindstheshortestpathbetweensanddinthegraphGfunctionDijkstra(G=(V,E,w),s,d)foreachvinVdist[v]=infinityprev[v]=nulldist[s]=0whileVnotemptyandunotequaltod//removethevertexuwiththeleastvalueofdist[u]u=extract_min(V)foreachneighborvofualt=dist[u]+w(u,v)ifaltdist[v]dist[v]=altprev[v]=u Name Parametrization plane S(u;v)=(u;v;0). sphere S(u;v)=(cos( )cos();cos( )sin();sin( ))where =(x 1=2)and=(y 1=2). torus S(u;v)=(3cos()+cos( )cos();3sin()+cos( )sin();sin( ))where=2uand =2v. saddle S(u;v)=(x;y;x2 y2)wherex=2u 1;y=2v 1. Table1.ParametrizationsofthesurfacesusedinSec-tion6. FINDINGGEODESICSONSURFACES27 (a)InitialPath (b)MP(200) (c)GD(200) (d)NR(5) 0 100 200 0.8 1 1.2 1.4 IterationsArc length (e)MP 0 100 200 0.8 1 1.2 1.4 IterationsArc length (f)GD 0 2 4 0.8 1 1.2 1.4 IterationsArc length (g)NR (h)InitialPath (i)MP(200) (j)GD(200) 0 100 200 0.8 1 1.2 1.4 IterationsArc length (k)MP 0 100 200 0.8 1 1.2 1.4 IterationsArc length (l)GDFigure13.Resultonplane. 28J.BAEK,A.DEOPURKAR,ANDK.REDFIELD (a)InitialPath (b)MP(200) (c)GD(200) (d)NR(5) 0 100 200 2.5 3 3.5 IterationsArc length (e)MP 0 100 200 2.5 3 3.5 IterationsArc length (f)GD 0 2 4 2.5 3 3.5 IterationsArc length (g)NR (h)InitialPath (i)MP(200) (j)GD(200) 0 100 200 2.5 3 3.5 IterationsArc length (k)MP 0 100 200 2.5 3 3.5 IterationsArc length (l)GDFigure14.Resultonsphere. FINDINGGEODESICSONSURFACES29 (a)InitialPath (b)MP(200) (c)GD(200) (d)NR(5) 0 100 200 6 6.5 7 7.5 8 IterationsArc length (e)MP 0 100 200 6 6.5 7 7.5 8 IterationsArc length (f)GD 0 2 4 6 6.5 7 7.5 8 IterationsArc length (g)NR (h)InitialPath (i)MP(200) (j)GD(200) 0 100 200 6 8 10 12 IterationsArc length (k)MP 0 100 200 6 8 10 12 IterationsArc length (l)GDFigure15.Resultontorus. 30J.BAEK,A.DEOPURKAR,ANDK.REDFIELD (a)InitialPath (b)MP(200) (c)GD(200) (d)NR(5) 0 100 200 1 2 3 4 IterationsArc length (e)MP 0 100 200 1 2 3 4 IterationsArc length (f)GD 0 2 4 1 2 3 4 IterationsArc length (g)NR (h)InitialPath (i)MP(200) (j)GD(200) 0 100 200 1 2 3 4 IterationsArc length (k)MP 0 100 200 1 2 3 4 IterationsArc length (l)GDFigure16.Resultonsaddle. FINDINGGEODESICSONSURFACES31References[CL]T.E.Cormen,C.E.Leiserson,R.L.Rivest,C.Stein:\IntroductiontoAlgorithms."MITPress,Cambridge1990.