A faster algorithm for the discrete Fréchet distance under translation ∗ RinatBenAvraham HaimKaplan MichaSharir † ‡ § Abstract 5 ThediscreteFréchetdistanceisausefulsimilaritymeasureforcomparingtwosequencesofpoints 1 0 P = (p1,...,pm) and Q = (q1,...,qn). In many applications, the quality of the matching can be 2 improvedifweletQundergosometransformationrelativetoP. Inthispaperweconsidertheproblem n offindingatranslationofQthatbringsthediscreteFréchetdistancebetweenP andQtoaminimum. a We devise an algorithm that computes the minimum discrete Fréchet distance under translation in R2, J and runs in O(m3n2(1+log(n/m))log(m+n)) time, assuming m n. This improves a previous ≤ 5 algorithmofJiangetal.[10],whichrunsinO(m3n3log(m+n))time. 1 ] G C . s c [ 1 v 4 2 7 3 0 . 1 0 5 1 : v i X r a ∗WorkbyHaimKaplanhasbeensupportedbyIsraelScienceFoundationgrantno. 822/10,andGerman-IsraeliFoundationfor ScientificResearchandDevelopment(GIF)grantno. 1161/2011,andIsraeliCentersofResearchExcellence(I-CORE)program (CenterNo. 4/11). WorkbyMichaSharirhasbeensupportedbyGrant892/13fromtheIsraelScienceFoundation,bytheIsraeli Centers of Research Excellence (I-CORE) program (Center No. 4/11), and by the Hermann Minkowski–MINERVA Center for GeometryatTelAvivUniversity. WorkbyMichaSharirandRinatBenAvrahamhasbeensupportedbyGrant2012/229fromthe U.S.-IsraeliBinationalScienceFoundation. †SchoolofComputerScience,TelAvivUniversity,TelAviv69978,Israel;[email protected] ‡SchoolofComputerScience,TelAvivUniversity,TelAviv69978,Israel;[email protected] §SchoolofComputerScience,TelAvivUniversity,TelAviv69978,Israel;[email protected] 1 1 Introduction Shape matching is an important area of computational geometry, that has applications in computer vision, patternrecognition,andotherfieldsthatareconcernedwithmatchingobjectsbyshapesimilarity. Generally, inshapematchingwearegiventwogeometricobjectsAandB andwewanttomeasuretowhatextentthey are similar. Usually we may allow certain transformations, like translations, rotations and/or scalings, of oneobjectrelativetotheother,inordertoimprovethequalityofthematch. Inmanyapplications,theinputdataconsistsoffinitesetsofpointssampledfromtheactualobjects. To measuresimilaritybetweenthesampledpointsets,variousdistancefunctionshavebeenused. Onepopular functionistheHausdorffdistancethatequalstothemaximumdistancefromapointinonesettoitsnearest point in the other set. However, when the objects which we compare are curves, sequences, or contours of larger objects, and the sampled points are ordered along the compared contours, the discrete Fréchet distancemay bea moreappropriate similarity measure. This isbecausethe discreteFréchet distancetakes intoaccounttheorderingofthepointsalongthecontourswhichtheHausdorffdistanceignores. Comparing curves and sequences is a major task that arises in computer vision, image processing and bioinformatics (e.g.,inmatchingbackbonesequencesofproteins). Thediscrete Fréchetdistance between asequence of pointsP and anothersequence of pointsQis de- finedastheminimum,overallpossibleindependent(forward)traversalsofthesequences,ofthemaximum distance between the current point of P and the current point of Q during the traversals. See below and in Section2foramoreformaldefinition. In this work, we focus on the problem of computing the minimum discrete Fréchet distance under translation. Thatis,giventwosequencesP andQofmandnpoints,respectively,intheplane,wewishto translateQbyavectort R2 suchthatthediscreteFréchetdistancebetweenP andQ+tisminimized. ∈ Background. TheFréchetdistancehasbeenextensivelystudiedduringthepast20years. Themainvariant, the continuous Fréchet distance, where no transformation is allowed, measures similarity between (polyg- onal) curves. It is the smallest δ for which there exist forward simultaneous traversals of the two curves, from start to end, so that at all times the distance between the corresponding points on the curves is at most δ. The discrete Fréchet distance considers sequences P and Q of points instead of curves. It is de- fined analogously, where (a) the simultaneous traversals of the sequences are represented as a sequence of pairs (p(1),q(1)),...,(p(t),q(t)), where p(i) P, q(i) Q, for i = 1,...,t, (b) the first (resp., last) pair ∈ ∈ consists of the starting (resp., terminal) points of the two sequences, and (c) each (p(i),q(i)) is obtained from (p(i 1),q(i 1)) by moving one (or both) point(s) to the next position in the corresponding sequence. − − Moststudiesoftheproblemconsiderthesituationwherenotranslation(oranyothertransformation)isal- lowed. Inthis“stationary”case,thediscreteFréchetdistanceintheplanecanbecomputed,usingdynamic programming, in O(mn) time (Eiter and Mannila [9]). Agarwal et al. [1] slightly improve this bound, (cid:18) (cid:19) mnloglogn and show that the (stationary) discrete Fréchet distance can be computed in O time on logn a word RAM, and a very recent result of Bringmann [5] indicates that a substantially subquadratic solu- tion (one that runs in time O((mn)1 δ), for some δ > 0) is unlikely to exist. Alt and Godau [3] showed − that the (stationary) continuous Fréchet distance of two planar polygonal curves with m and n edges, re- spectively, can be computed, using dynamic programming, in O(mnlogmn) time. This has been slightly improved recently by Buchin et al. [6], who showed that the continuous Fréchet distance can be computed in O(N2(logN)1/2(loglogN)3/2) time on a pointer machine, and in O(N2(loglogN)2) time on a word RAM (here N = m = n denotes the number of edges in each curve). In short, the best known algorithms forthestationarycase,forbothdiscreteandcontinuousvariants,hoveraroundthequadratictimebound. 1 Notsurprisingly,theproblemsbecomemuchharder,andtheirsolutionsmuchlessefficient,whentrans- lations (or other transformations) are allowed. For the problem of computing the minimum continuous Fréchetdistanceundertranslation,Altetal.[4]giveanalgorithmwithO(m3n3(m+n)2log(m+n))run- ning time, where m and n are the number of edges in the curves. They also give a (1+ε)-approximation algorithm for the problem, that runs in O(ε 2mn) time. That is, they compute a translation of one of the − curvesrelativetotheother,suchthattheFréchetdistancebetweentheresultingcurvesisatmost(1+ε)times theminimumFréchetdistanceunderanytranslation. Inthreedimensions,Wenk[14]showedthat,giventwo polygonalchainswithmandnedgesrespectively,theminimumcontinuousFréchetdistancebetweenthem, underanyreasonablefamilyoftransformations,canbecomputedinO((m+n)3f+2log(m+n))time,where f isthenumberofdegreesoffreedomformovingonechainrelativetotheother. Sowithtranslationsalone (f = 3), the minimum continuous Fréchet distance in R3 can be computed in O((m+n)11log(m+n)) time,andwhenbothtranslationsandrotationsareallowed(f = 6),thecorrespondingminimumcontinuous FréchetdistancecanbecomputedinO((m+n)20log(m+n))time. The situation with the discrete Fréchet distance under translation is somewhat better, albeit still ineffi- cient. Jiangetal.[10]showthat,giventwosequencesofpointsintheplane,theminimumdiscreteFréchet distancebetweenthemundertranslationcanbecomputedinO(m3n3log(m+n))time. Forthecasewhere bothrotationsandtranslationsareallowed, theygiveanalgorithmthatrunsinO(m4n4log(m+n))time. They also design a heuristic method for aligning two sequences of points under translation and rotation in three dimensions. Mosig et al. [12] present an approximation algorithm that computes the discrete Fréchet distanceundertranslation,rotationandscalingintheplane,uptoafactorcloseto2,andrunsinO(m2n2) time. Our results. Our algorithm improves the bound of Jiang et al. [10] by a nearly linear factor, with running timeO(m3n2(1+log(n/m))log(m+n)),assumingm n. Itusesa0/1-matrixM(P,Q)ofsizem n, ≤ × whose rows (resp., columns) correspond to the points of P (resp., of Q). Assuming a stationary situation, or, rather, a fixed translation of Q, an entry in the matrix is equal to 1 if and only if the distance between thetwocorrespondingpointsisatmostδ,whereδ issomefixeddistancethreshold. Weuse(i,j)todenote an entry in the matrix that corresponds to the points p and q , and we use M to denote its value. The i j i,j discreteFréchetdistanceisatmostδ ifandonlyifthereisarow-andcolumn-monotonepathofonesinM thatstartsat(1,1)andendsat(m,n)(seeSection2foramoreprecisedefinition). We can partition the plane of translations into a subdivision with O(m2n2) regions, so that, for all δ A translationsinthesameregion, thematrixM isfixed(forthefixedδ). Wethentraversetheregionsof , δ A moving at each step from one region to a neighboring one. Assuming general position, in each step of our traversal exactly one entry of M changes from 1 to 0 or vice versa. We present a dynamic data structure Γ(M)thatsupportsanupdateofanentryofM,inO(m(1+log(n/m)))time,assumingm n,1 andthen ≤ re-determines whether there is a monotone path of ones from (1,1) to (m,n), in O(1) additional time. If wefindsuchamonotonepathinM,wehavefoundatranslationt(actuallyawholeregionoftranslations2) such that the discrete Fréchet distance between P and Q + t is at most δ. Otherwise, when we traverse the entire and fail after each update, we conclude that no such translation exists. Using this procedure, δ A combinedwiththeparametricsearchingtechnique[11],weobtainanalgorithmforcomputingtheminimum discreteFréchetdistanceundertranslation. We reduce the dynamic maintenance of M to dynamic maintenance of reachability in a planar graph, as edges are inserted and deleted to/from the graph. Specifically, we can think of (the 1-entries of) M as a representation of a planar directed graph with N mn nodes. Each 1-entry of M corresponds to a node ≤ 1ThisiswithoutlossofgeneralityaswecanchangetherolesofmandnbyflippingM. 2Foracriticalvalueofδ,theregioncandegeneratetoasinglevertexofA ;seeSections3and6fordetails. δ 2 in the graph, and each possible forward move in a joint traversal is represented by an edge (see Section 2 for details). Then, determining whether there is a row- and column-monotone path of ones from (1,1) to (m,n)correspondstoareachabilityqueryinthegraph(from(1,1)to(m,n)). A data structure for dynamic maintenance of reachability in directed planar graphs was given by Sub- ramanian[13]. ThisdatastructuresupportsupdatesandreachabilityqueriesinO(N2/3logN)time,where N is the number of nodes in the graph. Diks and Sankowski [8] improved this data structure, and gave a structurethatsupportsupdatesandreachabilityqueriesinO(N1/2log3N)time. WegiveasimplerandmoreefficientstructureformaintainingreachabilityinM thatexploitsitsspecial structure. OurstructurecanupdatereachabilityinformationinM inO(m(1+log(n/m)))time,assuming m n,andanswersreachabilityquery(from(1,1)to(m,n))inO(1)time. Incontrast,thedatastructure ≤ of [8] applied in our context performs an update and a query in O((mn)1/2log3(mn)) time. Using our structure, we obtain an algorithm for computing the minimum discrete Fréchet distance under translation thatrunsinO(m3n2(1+log(n/m))log(m+n))time(again,assumingm n). ≤ Tosummarizethecontributionsofthispaperaretwofold: (a)Thereductionoftheproblemofcomputing the minimum discrete Fréchet distance to a dynamic planar directed graph reachability problem. (b) An efficient data structure for this reachability problem. For m n our structure is faster than the general ≈ reachabilitystructureof[8]byapolylogarithmicfactor,andwhenm ntheimprovementisconsiderably (cid:112) (cid:28) moresignificant(roughlybyafactor n/m). Moreover,ourdatastructureissimplerthanthatofDiksand Sankowski. 2 Preliminaries We now define the (stationary) discrete Fréchet distance formally. Let P = (p ,...,p ) and Q = 1 m (q ,...,q )bethetwoplanarsequencesdefinedintheintroduction. 1 n Forsomefixeddistanceδ > 0wedefinea0/1-matrixM (P,Q)formallyasfollows. Therows(resp., δ columns) of M (P,Q) correspond to the points of P (resp., of Q) in their given order. An entry (i,j) of δ M (P,Q)is1ifthedistancebetweenp andq isatmostδ,andis0otherwise. wedenoteM (P,Q)byM δ i j δ whenP andQandδ areclearfromthecontext. ThedirectedgraphG (P,Q)associatedwithP,Qandδhasavertexforeachpair(p ,q ) P Qand δ i j ∈ × anedgeforeachpairofadjacentonesinM (P,Q). Specifically,wehaveanedgefrom(p ,q )to(p ,q ) δ i j i+1 j ifandonlyifboth(i,j)and(i+1,j)are1inM,anedgefrom(p ,q )to(p ,q )ifandonlyifboth(i,j) i j i j+1 and(i,j+1)are1inM,andanedgefrom(p ,q )to(p ,q )ifandonlyifboth(i,j)and(i+1,j+1) i j i+1 j+1 are1inM. wedenoteG (P,Q)byGwhenP andQandδ areclearfromthecontext. δ The(stationary)discreteFréchetdistancebetweenP andQ,denotedbyδ (P,Q),isthesmallestδ > 0 ∗ for which (p ,q ) is reachable from (p ,q ) in G . Informally, think of P and Q as two sequences of m n 1 1 δ steppingstonesandoftwofrogs,theP-frogandtheQ-frog,wheretheP-froghastovisitalltheP-stones in order and the Q-frog has to visit all the Q-stones in order. The frogs are connected by a rope of length δ, and are initially placed at p and q , respectively. At each move, either one of the frogs jumps from its 1 1 current stone to the next one and the other stays at its current stone, or both of them jump simultaneously fromtheircurrentstonestothenextones. Furthermore,suchajumpisallowedonlyifthedistancesbetween the two frogs before and after the jump are both at most δ. Then δ (P,Q) is the smallest δ > 0 for which ∗ thereexistsasequenceofjumpsthatgetsthefrogstop andq ,respectively. m n TheproblemofcomputingtheminimumdiscreteFréchetdistanceundertranslation,asreviewedinthe introduction,istofindatranslationtsuchthatδ (P,Q+t)isminimized. ∗ We say that an entry (i,j) of M is reachable from an entry (k,l), with k i,l j, if (p ,q ) is i j ≤ ≤ 3 reachablefrom(p ,q )inG. Apathfrom(p ,q )to(p ,q )inGcorrespondstoa(weakly)row-monotone k l k l i j and column-monotone sequence of ones in M connecting the one in entry (k,l) to the one in entry (i,j). This is sequence consists of three kinds of moves: 1) upward moves between entries of the form (r,s) to (r+1,s)inwhichtheP-frogmovesfromp top , 2)rightmovesbetweenentriesoftheform(r,s)to r r+1 (r,s+1) in which the P-frog moves from q to q , and 3) diagonal moves between entries of the form s s+1 (r,s) to (r+1,s+1) in which the P-frog moves from q to q both frogs move simultaneously — the s s+1 P-frogfromp top , andtheQ-frogfromq toq . SeeFigure1. Wecallsuchamonotonesequence r r+1 s s+1 of ones in M a path in M from (k,l) to (i,j). To determine whether δ (P,Q) δ, we need to determine ∗ ≤ whetherthereissuchapathinM thatstartsat(1,1)andendsat(m,n). Wesaythatanentry(i,j)ofM is reachableifthereisapathfrom(1,1)to(i,j). Wedenotetheconcatenationoftwopathsπ ,π byπ π ,assumingthatthelastentryofπ isthefirst 1 2 1 2 1 · entryofπ ;thisentryappearsonlyonceintheconcatenation. 2 p8 1 1 1 p8 1 1 1 M p7 1 1 1 1 1 1 p7 1 1 1 1 1 1 p6 1 1 1 1 p6 1 1 1 1 p5 1 1 1 1 1 1 1 p5 1 1 1 1 1 1 1 p4 1 1 1 1 1 1 1 p4 1 1 1 1 1 1 1 p3 1 1 1 1 1 1 1 1 p3 1 1 1 1 1 1 1 1 p2 1 1 1 p2 1 1 1 p p 1 1 1 1 1 1 1 1 1 1 1 1 1 1 q q q q q q q q q q q q q q q q q q q q q q q q 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 (a) (b) Figure1: (a)Areachabilitypathfrom(4,3)to(7,7). (b)Areachabilitypathfrom(1,1)to(8,12). WeuseadecompositionofM intoblocks. AblockisasubmatrixofM thatcorrespondstocontiguous subsequencesP andQ ofP andQ,respectively. WedenotebyM(P ,Q)theblockofM(P,Q)formed (cid:48) (cid:48) (cid:48) (cid:48) by P and Q. Consider a block M(P ,Q) of M(P,Q). Let p (resp., p+) denote the first (resp., last) (cid:48) (cid:48) (cid:48) (cid:48) − pointofP andletq (resp.,q+)denotethefirst(resp.,last)pointofQ. WecalltheentriesofM(P ,Q) (cid:48) − (cid:48) (cid:48) (cid:48) correspondingto p Q P q theinputboundaryofM(P ,Q),anddenoteitbyM(P ,Q) − (cid:48) (cid:48) − (cid:48) (cid:48) (cid:48) (cid:48) − { }× ∪ ×{ } (thecommonentrycorrespondingto p q appearsonlyonceintheboundary). Similarly,wecallthe − − { }×{ } entriesofM(P ,Q)correspondingto p+ Q P q+ theoutputboundaryofM(P ,Q),anddenote (cid:48) (cid:48) (cid:48) (cid:48) (cid:48) (cid:48) { }× ∪ ×{ } it by M(P ,Q)+ (with a similar suppression of the duplication of the common element p+ q+ ). (cid:48) (cid:48) { } × { } Note that there is a two-entry overlap between the input and output boundaries. We enumerate the entries of M(P ,Q) by first enumerating the entries of p Q from right to left (i.e., backwards) and then (cid:48) (cid:48) − − (cid:48) { }× theremainingentriesofP q frombottomtotop(forward). WeenumeratetheentriesofM(P ,Q)+ (cid:48) − (cid:48) (cid:48) ×{ } by first enumerating the entries of P q+ from bottom to top (forward) and then the remaining entries (cid:48) ×{ } of p + Q fromrighttoleft(backwards). Informally, M(P ,Q) isenumeratedin“clockwise”order, (cid:48) (cid:48) (cid:48) − { } × while M(P ,Q)+ is enumerated in “counterclockwise” order; see Figure 2. For two entries i,j of this (cid:48) (cid:48) enumerationofaninputoroutputboundaryB ofM(P ,Q),weuse[i,j]todenotethesequenceofentries (cid:48) (cid:48) (i,i+1,...,j)ofB. We also use the following definitions. We call the entries corresponding to P q the vertical (cid:48) − × { } input boundary of M(P ,Q), and denote it by M(P ,Q) . We call the entries corresponding to P (cid:48) (cid:48) (cid:48) (cid:48) − (cid:48) × q+ theverticaloutputboundaryofM(P ,Q),anddenoteitbyM(P ,Q)+. Thatis,M(P ,Q) and (cid:48) (cid:48) (cid:48) (cid:48) (cid:48) (cid:48) − { } M(P ,Q)+ aretheverticalpartsofM(P ,Q) andM(P ,Q)+,respectively. Weenumeratetheentries (cid:48) (cid:48) (cid:48) (cid:48) − (cid:48) (cid:48) ofeachverticalboundaryfrombottomtotop. 4 output-boundaryM(P0,Q0)+ p8 1 1 1 M p7 1 1 1 1 1 1 91 8 7 6 51 41 M(P0,Q0) p6 1 1 1 1 8 3 p5 1 1 1 1 1 1 1 1 1 1 1 p4 1 1 1 1 1 1 1 7 1 1 21 p3 1 1 1 1 1 1 1 1 p2 1 1 1 6 5 1 4 31 2 1 1 p 1 1 1 1 1 1 1 q q q q q q q q q q q q 1 2 3 4 5 6 7 8 9 10 11 12 input-boundaryM(P0,Q0)− (a) (b) Figure2: (a)A(highlighted)blockM(P ,Q)ofM(P,Q),whereP =(p ,p ,p ,p )andQ =(q ,q ,...,q ). (cid:48) (cid:48) (cid:48) 5 6 7 8 (cid:48) 7 8 12 (b)TheinputboundaryM(P ,Q) andtheoutputboundaryM(P ,Q)+ ofM(P ,Q)aremarked,withtheorder- (cid:48) (cid:48) − (cid:48) (cid:48) (cid:48) (cid:48) ingsoftheirelements. 3 The subdivision of the plane of translations We first consider the corresponding decision problem. That is, given a value δ > 0, we wish to decide whetherthereexistsatranslationt R2 suchthatδ (P,Q+t) δ. ∗ ∈ ≤ For a point x R2, let D (x) be the disk of radius δ centered at x. Given two points p P and δ i ∈ ∈ q Q,considerthediskD (p q ),andnoticethatt D (p q )ifandonlyif (p q ) t δ(or j δ i j δ i j i j ∈ − ∈ − (cid:107) − − (cid:107) ≤ p (q +t) δ). Thatis,D (p q )ispreciselythesetoftranslationstforwhichq +tisatdistance i j δ i j j (cid:107) − (cid:107) ≤ − atmostδ fromp . i Weconstructthearrangement = (P,Q)ofthedisksin = D (p q ) (p ,q ) P Q . δ δ δ i j i j A A D { − | ∈ × } Weassumegeneralpositionofthepoints. Thatis,weassumethat(a)nomorethantwoboundariesofthese disksintersectinacommonvertexof ,and(b)nopairofthedisksaretangenttoeachother. Nevertheless, δ A such a degeneracy can arise when δ is a critical value (see Section 6 for details about critical values of δ that arise during the optimization procedure), but we assume that at most one such degeneracy can happen foragivenδ. Sincethenumberofdisksismn,thecombinatorialcomplexityof isO(m2n2). Letf be δ A afaceof ofanydimension0,1or2(byconvention,f isassumedtoberelativelyopen),andlett f be δ A ∈ atranslation. Then,forpointsp P,q Q,q +tisatdistanceatmostδ fromp ifandonlyifthedisk i j j i ∈ ∈ D (p q )containsf (otherwise,thediskisdisjointfromf). Sincethisholdsforeveryt f,itfollows δ i j − ∈ thatf correspondstoauniquepairwise-distancesmatrixM(P,Q+t),foranyt f. Wedenotethismatrix ∈ byM(P,Q+f),forshort. The setup just described leads to the following naive solution for the decision problem. Construct the arrangement for the given distance δ, and traverse its faces. For each face f , form the δ δ A ∈ A correspondingpairwise-distancesmatrixM(P,Q+f),andsolvethe(stationary)discreteFréchetdistance decisionproblemforP andQ+f usingastraightforwarddynamicprogrammingonM(P,Q+f)(orthe moresophisticatedslightlysubquadraticalgorithmofAgarwaletal.[1]). Ifδ (P,Q+f) δforsomeface ∗ ≤ f,weconcludethatthereexistsatranslationtsuchthatδ (P,Q+t) δ (anytranslationt f woulddo). ∗ ≤ ∈ Iftheentirearrangement istraversedandnofacef of satisfiesδ (P,Q+f) δ,wedeterminethat δ δ ∗ A A ≤ δ (P,Q+t) > δ for all translations t R2. The complexity of is O(m2n2), and solving the discrete ∗ δ ∈ A Fréchetdistancedecisionproblemforeachfaceof takesO(mn)time(orslightlyless,asin[1]). Hence, δ A thesolutionjustdescribedforthedecisionproblemtakes(slightlylessthan)O(m3n3)time. Jiang et al. [10] used an equivalent solution for the decision problem, that takes the same asymptotic runningtime. Rephrasingtheirprocedureintermsof ,theytestwhetherδ (P,Q+t) δfortranslations δ ∗ A ≤ t corresponding to the vertices of , and over an additional set of mn translations, one chosen from the δ A 5 boundaryofeachdisk. Thecorrectnessofthisapproachfollowsbyobservingthatiff isafaceof andt δ A isanypointon∂f,thenallthe1-entriesofM(P,Q+f)arealso1-entriesofM(P,Q+t),soitsufficesto testtheverticesoff,or,iff hasnovertices,totestanarbitrarypointt ∂f. Wewillusethisobservation ∈ inourimplementationoftheoptimizationprocedure. Our naive solution is similar to the algorithm of Jiang et al. [10], in the sense that they both discretize the set of possible translations. However, our solution is more suitable for the improvement of this naive bound,thatwepresentinSection5,sinceitallowsustotraversethesetofpossibletranslationsinamanner that introduces only a single change in M(P,Q+f), when we move from one face f of translations to a neighboringone. To exploit this property we need a data structure that maintains reachability data for M, and updates it efficiently after each change. We present this structure in two stages. First, in Section 4, we present a compact reachability structure for blocks of M(P,Q+f), which is the main building block of the overall structure. Then, in Section 5, we present the overall data structure, and show how to use it to improve the naivesolutionsketchedabovebyanearlylinearfactor. 4 Compact representation of reachability in a block Let B be a block of M = M(P,Q + f) of size r c, and suppose that we have already computed the reachableentriesofB andwethenwishtocompute×thereachableentriesofB+. IftheσAe(1n)triesoftheblock − aregivenexplicitly,thiscanbedoneinO(rc)timeusingdynamicprogramming(orslightlyfasterusingthe algorithmof[1]). Ourgoalinthissectionistodesignadatastructure,thatwedenoteasΦ(B),thatallows ustocomputethereachableentriesofB+ fromthereachableentriesofB ,inO(r+c)time. Theoverall − data structure itself is constructed recursively from these block structures (see Section 5 for details), and implicitly accesses all the entries of B. The advantage of using this block decomposition is that updating thestructurecanbedonemoreefficiently. σZ(j) σZ(i) B 1 1 1 1 1 1 1 σ(i) σ(j) B+ 1 1 1 1 1 1 1 1 1 σA(j) π(e,σ(i)) 1 1 1 1 1 e 1 1 1 1 1 1 1 1 σA(i) π(j,e) π(e,σ(j)) 1 1 1 1 1 1 1 1 1 1 1 B+ j B 1 1 1 1 − π(i,e) 1 1 1 1 1 1 B− B 3··· j 1 1 1 1 1 1 1 1 11 13··· σZ(1) 1 1 1 1 1 12 ··· 2 i 3 2 1 1 1 1 13 211 σA(1) i (a) (b) Figure3: (a)Twoentriesi,jofB andtwoentriesσ(i)andσ(j)ofB+thatarereachablefromiandj,respectively. − Since i < j and σ(i) > σ(j), σ(j) is reachable also from i, and σ(i) is reachable also from j. (b) The intervals [σ (k),σ (k)], for any 1-entry k of B , are either disjoint or overlap in a common subinterval. Neither of these A Z − intervalscanstrictlycontainbothendpointsoftheother. Theseintervalsaredefined(andshowninthefigure)onlyfor 1-entriesofB . − 6 Observation4.1. LetBbeablockofM andleti,jbetwoentriesofB suchthatj > i(inthe“clockwise” − orderdefinedinSection2). Letσ(i)beanentryofB+ thatisreachablefromi,andletσ(j)beanentryof B+ thatisreachablefromj. Ifσ(j) < σ(i)(inthecorresponding“counterclockwise”order)thenσ(j)is alsoreachablefromi,andσ(i)isalsoreachablefromj. Proof. SeeFigure3(a). Sinceσ(i)isreachablefromi,thereisa(monotone)pathπ(i,σ(i))fromitoσ(i)in M. Similarly,sinceσ(j)isreachablefromj,thereisa(monotone)pathπ(j,σ(j))fromjtoσ(j). Sincei < j andσ(j) < σ(i),π(i,σ(i))mustcrossπ(j,σ(j))(i.e.,thereexistsa1-entrye π(i,σ(i)) π(j,σ(j)). ∈ ∩ Hence, π(i,σ(i)) can be decomposed into two subpaths π(i,e),π(e,σ(i)) such that π(i,σ(i)) = π(i,e) · π(e,σ(i)),andπ(j,σ(j))canbesimilarlydecomposedasπ(j,σ(j)) = π(j,e) π(e,σ(j)). Asaresult,the · pathsπ(i,σ(j)) = π(i,e) π(e,σ(j))andπ(j,σ(i)) = π(j,e) π(e,σ(i))arealso(monotone)paths, and · · theclaimfollows. Corollary 4.2. Let B be a block of M, let i be an entry of B and let σ (i),σ (i) be two entries in B+ − 1 2 thatarebothreachablefromi,withσ (i) < σ (i). Ifthereexistsanentryσ(j)thatisreachablefromsome 1 2 j B ,suchthatσ (i) < σ(j) < σ (i),thenσ(j)isalsoreachablefromi. − 1 2 ∈ Proof. ByObservation4.1, ifi < j, thenσ(j)isreachablefromisinceσ(j) < σ (i). Ifi > j, thenσ(j) 2 isreachablefromisinceσ (i) < σ(j). 1 The corollary is applied as follows. Let i be an entry of B , let σ (i) and σ (i) denote the first and − A Z lastentriesinB+ thatarereachablefromi. (Notethatfortheseentriestobedefined,thevalueoftheentry i must be 1. Symmetrically, the values of both σ (i) and σ (i), if defined, must be equal to 1.) Then the A Z interval[σ (i),σ (i)]canonlycontainentriesofthefollowingthreetypes. A Z 1. 1-entriesthatarereachablefromi. 2. 0-entries. 3. 1-entriesthatarenotreachablefromi,norfromanyotherentryofB . − In other words, [σ (i),σ (i)] cannot contain 1-entries that are reachable from some j in B and not A Z − fromi. Corollary 4.3. Let B be a block of M and let i and j be two entries of B such that j > i. Then − σ (j) σ (i)andσ (j) σ (i). A A Z Z ≥ ≥ Proof. Assumetothecontrarythatσ (i) > σ (j). Then,accordingtoObservation4.1,σ (j)isreachable A A A fromi. Hence,σ (i) σ (j),acontradiction. Similarly, ifσ (j) < σ (i),thenObservation4.1implies A A Z Z ≤ thatσ (i)isreachablefromj. Hence,σ (j) σ (i),contradiction. Z Z Z ≥ Inotherwords,[σ (i),σ (i)]and[σ (j),σ (j)]canbeeitherdisjointoroverlapinacommonsubinter- A Z A Z val,buttheycannotbeproperlynestedinsideoneanother(thatis,neitheroftheseintervalscancontainboth endpointsoftheotherinits“interior”). Note,however,thatoneintervalcanweaklycontaintheother. That is,ifoneintervalcontainstheother,theneitherσ (i) = σ (j)orσ (i) = σ (j),orboth. SeeFigure3(b). A A Z Z Let B be a block of M of size r c. We construct a data structure Φ(B) for B, which stores the × followinginformation. (Hereweonlyspecifythestructure;itsconstructionisdetailedinSection5.) 1. Foreach1-entryiofB westore − (a) thefirstentryσ (i)ofB+ thatisreachablefromi,and A 7 (b) thelastentryσ (i)ofB+ thatisreachablefromi. Z 2. Foreach1-entryj ofB+ westore (a) aflagf(j)indicatingwhetherj isreachablefromsomeentryofB . − (b) alistL (j)ofthe1-entriesi B suchthatσ (i) = j,and A − A ∈ (c) alistL (j)ofthe1-entriesi B suchthatσ (i) = j. Z − Z ∈ Lemma4.4. GiventhedatastructureΦ(B)forablockB, andgiventheentriesofB thatarereachable − from(1,1),wecandetermine,inO(r+c)time,theentriesofB+ thatarereachablefrom(1,1). Proof. Wegooverthereachable1-entriesofB inorder. Foreachsuchentryi,wegoovertheentriesinthe − intervalI(i) = [max σ (i ),σ (i) ,σ (i)]ofB+,wherei isthepreviousreachable1-entryofB (for Z − A Z − − { } the first reachable entry i of B , I(i) = [σ (i),σ (i)] and i is undefined). Note that, by Corollary 4.3, − A Z − max σ (i ),σ (i) [σ (i),σ (i)]soI(i) [σ (i),σ (i)]. (Theentriesof[σ (i),σ (i)]thatprecede Z − A A Z A Z A Z { } ∈ ⊆ max σ (i ),σ (i) were already processed when we went over I(i ) or over intervals associated with Z − A − { } earlierindices.) For each 1-entry j of I(i) that is reachable from some entry of B (according to the flag f(j)), we − determine that j is reachable also from (1,1). Since we traverse each interval [σ (i),σ (i)] starting from A Z max σ (i ),σ (i) ,theinternalportionsofthesubintervalsthatweinspectarepairwisedisjoint,implying Z − A { } thattherunningtimeislinearinr+c. Weomitthestraightforwardproofofcorrectnessofthisprocedure. 5 Dynamic maintenance of reachability in M(P,Q + f) We present a data structure, that uses the compact representation of reachability in a block of the previous section, to support an update of a single entry of M in O(m(1+log(n/m))) time, assuming m n. We ≤ presentthisdatastructureintwostages. First,inSection5.1,weshowhowtosupportanupdateofasingle entry in O(m) time, in the case where M is a square matrix of size m m. Then, in Section 5.2, we × generalize this data structure to support an update of a single entry in O(m(1 + log(n/m))) time, in the general case where M is an m n matrix with m n (the case m n is treated in a fully symmetric × ≤ ≥ manner). In Section 5.3, we describe the overall decision procedure that improves the naive solution sketched in Section3,usingthisdynamicdatastructure. 5.1 Adynamicdatastructureforreachabilitymaintenanceinasquarematrix We store the reachability data of M(P,Q+f) (of some arbitrary face f from which we start the traversal of the arrangement ) in a so-called decomposition tree Γ, by halving P and Q alternately. That is, the δ A root v of Γ corresponds to the entire matrix M(P,Q+f) and we store at v the reachability information Φ(M(P,Q+f)), as described in the previous section. (The actual construction of the reachability data, at all nodes of Γ, is done bottom-up, as described below.) In the next level of Γ we partition P into two subsequences P ,P , of at most m/2 + 1 points each, such that the last point of P is the first point 1 2 1 (cid:98) (cid:99) of P , and obtain a corresponding “horizontal” partition of M(P,Q+f) into two blocks M(P ,Q+f), 2 1 M(P ,Q+f), each of size at most ( m/2 +1) m, with a common “horizontal” boundary. We create 2 (cid:98) (cid:99) × twochildrenv ,v ofv andstoreateachv thereachabilityinformationΦ(M(P ,Q+f)),fori = 1,2. In 1 2 i i 8 1 By=M(Pv∪Pw,Qy) σy(i) B+ σZy(‘) σAy(‘) B+ w w Bw=M(Pw,Qy) σw(j) Bw=M(Pw,Qy) Bw− Bw− j σv(‘) σv(‘) Z A B+ B+ v v Bv=M(Pv,Qy) σv(k) ‘ i Bv− k Bv− (a) (b) Figure 4: A block B = M(P ,Q ), corresponding to a node y of Γ, is composed of the blocks of the children y y y v,wofy. TheblockB =M(P ,Q )correspondingtotheleftchildvliesbelowtheblockB =M(P ,Q )that v v v w w w correspondstotherightchildw,andwehaveP =P P andQ =Q =Q . y v w y v w ∪ (a)σy(i)andσy(k)areexamplesofreachableentriesofB+ andwehaveσy(i) = σw(σv(i))andσy(k) = σv(k). y σw(j)isanentryofB+thatisreachablefromB ,butitisnotareachableentryofB+(fromB )sinceallthepaths w w− y y− inB thatleadtoσw(j)gothroughentriesofB+thatarenotreachablefromB . y v v− (b)σy((cid:96))=σw(σv((cid:96)))andσy((cid:96))=σw(σv((cid:96))). A A A Z Z Z the next level of Γ, we partition Q into two subsequences Q ,Q , of at most m/2 +1 points each, such 1 2 (cid:98) (cid:99) thatthelastpointofQ isthefirstpointofQ ,andobtainacorresponding“vertical”partitionofeachblock 1 2 M(P ,Q+f),i 1,2 ,intotwoblocksM(P ,Q +f),j 1,2 ,eachofsizeatmost( m/2 +1) i i j ∈ { } ∈ { } (cid:98) (cid:99) × ( m/2 +1), withacommonverticalboundary. Weconstructfourrespectivegrandchildren, andstorethe (cid:98) (cid:99) correspondingreachabilitystructuresΦ(M(P ,Q +f))atthesenodes. Wecontinuerecursivelytopartition i j eachblockbyhalvingithorizontallyorvertically,alternately,inthesamemanner,untilwereachblocksof size 2 2. For each node v of Γ, let P and Q denote the subsequences of P and Q that form the block v v × M(P ,Q +f)thatisassociatedwithv. Tosimplifythenotation,wedenoteΦ(M(P ,Q +f))asΦ ,for v v v v v eachnodev. ThereachabilitydataΦ atthenodesvofΓiscomputedbyabottom-uptraversalofΓ,startingfromthe v leaves. TheconstructionofΦ(M(P ,Q +f))ataleafv istrivial,andtakesconstanttime. Thefollowing v v lemmaprovidesanefficientprocedureforconstructingthereachabilitydataatinnernodesofΓ. Lemma5.1. LetybeaninnernodeofΓwithleftandrightchildrenvandw,wheretheblocksstoredatv,w haveacommonhorizontalboundary. GiventhereachabilitydataΦ ,Φ ,thedataΦ canbecomputedin v w y O( P + Q ) time. An analogous statement holds when the common boundary of the children blocks is y y | | | | vertical. Proof. Notethatinthesetupofthelemma,wehaveQ = Q = Q andP = P P . Byconstruction, y v w y v w ∪ M(P ,Q )liesbelowM(P ,Q ). DenoteM(P ,Q )byB ,M(P ,Q )byB ,andM(P ,Q )byB . v y w y v y v w y w y y y For each entry i of B , denote by σy(i) (resp., σy(i)) the first (resp., last) entry of B+ that is reachable y− A Z y fromi. Wealsouseσy(i)todenoteanentryofB+ thatisreachablefromiinB . Analogousnotationsare y y usedforthechildrenblocksB ,B . SeeFigure4. v w We first copy the reachability information from the boundaries of B and B to the boundary of B v w y (except for the “interior” portion B of the common boundary B+ B of B and B , which is not a v∗w v ∩ w− v w 9