Measuring similarity between curves on 2-manifolds via homotopy area

Measuring the similarity of curves is a fundamental problem arising in many application fields. There has been considerable interest in several such measures, both in Euclidean space and in more general setting such as curves on Riemannian surfaces or curves in the plane minus a set of obstacles. However, so far, efficiently computable similarity measures for curves on general surfaces remain elusive. This paper aims at developing a natural curve similarity measure that can be easily extended and computed for curves on general orientable 2-manifolds. Specifically, we measure similarity between homotopic curves based on how hard it is to deform one curve into the other one continuously, and define this "hardness" as the minimum possible surface area swept by a homotopy between the curves. We consider cases where curves are embedded in the plane or on a triangulated orientable surface with genus $g$, and we present efficient algorithms (which are either quadratic or near linear time, depending on the setting) for both cases.


Introduction
Measuring curve similarity is a fundamental problem arising in many application fields, including graphics, computer vision, and geographic information systems.Traditionally, much research has been done on comparing curves embedded in the Euclidean space.However, in many cases it is natural to study curves embedded in a more general space, such as a terrain or a surface.
In this paper, we study the problem of measuring curve similarity on surfaces.Specifically, given two simple homotopic curves embedded on an orientable 2-manifold (including the plane), we measure their similarity by the minimum total area swept when deforming one curve to the other (the "area" of the homotopy between them), and present efficient algorithms to compute this new measure.

Related work.
From the perspective of computational geometry, the most widely studied similarity measures for curves is the Fréchet distance.Intuitively, imagine that a man and his dog are walking along two paths with a leash between them.The Fréchet distance between these two paths is the minimum leash length necessary for them to move from one end of the paths to the other end without back-tracking.Since the Fréchet distance takes the "flow" of the curves into account, in many settings it is a better similarity measure for curves than alternatives such as the Hausdorff distance [5,6].
Given two polygonal curves P and Q with n total edges in IR d , the Fréchet distance can be computed in O(n 2 log n) time [4].An Ω(n log n) lower bound for the decision problem in the algebraic computation tree model is known [11], and Alt has conjectured that the decision problem is 3SUM-Hard [2].Recently, Buchin et al. [12] show that there is a real algebraic decision tree to solve the Fréchet problem with subquadratic depth, suggesting that perhaps this is not the case.They also give an improved algorithm which runs in O(n 2 √ log n(log log n) 2 ) time.Very recently, Agarwal et al. present a novel approach to compute the discrete version of the Fréchet distance between two polygonal curves in sub-quadratic time [1].This is the first algorithm for any variant of the Fréchet distance to have a sub-quadratic running time for general curves.No previous algorithm, exact or approximate, with running time o(n 2 ) is known for general curves, although sub-quadratic approximation algorithms for special families of curves are known [6,7,23].
While the Fréchet distance is a natural curve similarity measure, it is sensitive to outliers.Variants of it, such as the summed-Fréchet distance, and the partial Fréchet similarity, have been proposed [14,15,24], usually at the cost of further increasing the time complexity.
The problem of extending and computing the Fréchet distance to more general metric space has also received much attention.Geodesic distance between points is usually considered when the underlying domain is not IR d .For example, Maheshwari and Yi [29] computed the geodesic Fréchet distance between two polygonal paths on a convex polytope in roughly O(n 3 K 4 log(Kn)) time, where n and K are the complexity of the input paths and of the convex polytope, respectively.Raichel and Har-Peled consider approximating the weak Fréchet distance between simplicial complexes in R d [27].Geodesic Fréchet distance between polygonal curves in the plane within a simple polygon has also been studied [8,19,25].
Rather than comparing distance between only two curves, Buchin et.al. [13] propose the concept of a median in a group of curves (or trajectories, in their setting).They give two algorithms to compute such a median.The first is based simply on the concept of remaining in the middle of the set of curves; this algorithm, while fast and simple, has a drawback in that the representative curve might not capture relevant features shared by a majority of the input curves.Their second algorithm addresses this issue by instead isolating a subset of relevant curves which share the same homotopy type with respect to obstacles that are placed in empty regions of the blame; it then computes a medial curve from this relevant subset.
One issue with generalizing Fréchet distance directly to surfaces is that the underlying topology is not taken into account; for example, in geodesic Fréchet distance, while the length of the leash varies continuously, the actual leash itself does not.As a result, several measures of similarity have been proposed which take the underlying topology into account.Chambers et al. [16] proposed the so-called homotopic Fréchet distance and gave a polynomial (although not efficient) algorithm for when the curves reside in a planar domain with a set of polygonal obstacles.The extra requirement for this homotopic Fréchet distance is that the leash itself and not just its length has to vary in a continuous manner, essentially restricting the homotopy class which the leash is in.A stronger variant called isotopic Fréchet distance has also been proposed and investigated, although no algorithms at all are known to even approximate this distance [17].
Orthogonal to homotopic Fréchet distance is the concept of the height of a homotopy; instead of minimizing the maximum leash length, this measure views the homotopy as tracing a way for the first curve to deform to the second curve, where the goal is to minimize the longest intermediate curve length.Introduced independently in two very different contexts [10,18], it is not even known if the problem is in NP.
Recent work on approximating the homotopy height and the homotopic Fréchet distance has yielded efficient O(log n) approximation algorithms for both of these problems [26].However, exact algorithms on surfaces for either problem are still unknown.

New work.
In this paper, we develop a natural similarity measure for curves on general surfaces that can be computed both quickly and exactly.Intuitively, we measure distances between homotopic curves based on how hard it is to deform one curve into the other one, and define this "hardness" as the minimal total surface area swept by a homotopy between them, which we call the optimal homotopy area.Our similarity measure is natural, and robust against noise (as the area in a sense captures average, instead of maximum, deviation from one curve to the other).To the best of our knowledge, this is the first similarity measure for curves on general surfaces with efficient polynomial-time algorithms to compute it exactly.
It is worth noting that this definition in a way combines homotopic Fréchet distance with homotopy height; those measures compute the "width" and "height" of the homotopy, while our measure calculates the total area.It is thus interesting that while no exact algorithms are known for either of those measures on surfaces, we are able to provide a polynomial running time for computing the area of a homotopy.
We consider both cases where curves are embedded in the plane, or on a closed, triangulated orientable surface with genus g.For the former case, our algorithm runs in O(n log n + I 2 log I) time, where n is the total complexity of input curves and I is the number of intersections between them.On a surface, if the input is a triangulation of complexity N , then our algorithm runs in time O(I 2 log I + ng log n + N ).While our similarity measure is more expensive to compute for the case of curves in the plane than the Fréchet distance when I = ω(n), one major advantage is that this measure can be computed on general orientable surfaces efficiently.In fact, the ideas and algorithms behind the planar case form the foundation for the handling of the case on general surfaces.
The main ideas behind our approach are developed by examining some properties of one natural class of homotopies, including a relation with the winding number of a closed curve.Specifically, the use of the winding number enables us to compute the optimal homotopy area efficiently in the plane, where the homotopy is restricted to be piecewise differential and regular.This forms the basis of our dynamic programming framework to compute similarity between curves in the plane.We also show how to build efficient data structures to keep the total cost of the dynamic program low.
For the case where the underlying surface is a topological sphere, we extend the winding number in a natural way and show how to adapt our planar algorithm without additional blow-up in the time complexity.For the case when the surface has non-zero genus, we must extend our algorithm to run efficiently in the universal cover (which is homeomorphic to the plane) by using only a small portion of it.
We remark that the idea of measuring deformation areas has been used before in practice [20,30].For example, similarity between two convex polygons can be measured by their symmetric difference [3,36]; we note that this is not equivalent to homotopy area, although it may be the same value in some situations.In another paper, the area sandwiched between an x-monotone curve and another curve is used to measure their similarity [9].However, computing the "area" between general curves has not been investigated prior to this work.

Definitions and Background
Paths and cycle.We will assume that we are working on an orientable 2-manifold M (which could be the plane).A curve (or a path) on a surface M is a map P : [0, 1] → M ; a cycle (or a loop) is a continuous map γ : S 1 → M where S 1 is the unit circle.A curve P or a cycle γ is simple if P (t 1 ) = P (t 2 ) (resp.γ(t 1 ) = γ(t 2 )) for any t 1 = t 2 .
Homotopy A homotopy between two paths P and Q (with the same endpoints) is a continuous map H : A homotopy describes a continuous deformation between the two paths or curves: for any value t ∈ [0, 1], we let H t = H(t, •) be the intermediate curve at time t, where H 0 = P and H 1 = Q.
We define the area of a homotopy H to be the total area covered by the image of the homotopy on the surface, where an area that is covered multiple times will be counted with multiplicity.More precisely, given a homotopy H whose image is piecewise differentiable, . The minimum homotopy area between P and Q is the infimum of the areas of all homotopies between P and Q, denoted by σ(P, Q).If such an infimum does exist and can be achieved by a homotopy, we call that homotopy an optimal homotopy.
We note that it is not immediately clear that this value exists, depending on the curves and underlying homotopy.Minimum area homotopies were considered by Douglas [22] and Rado [32] in the context of Plateau's problem; they noted that not only is the integral improper in general, but the infimum itself may not be continuous.The eventual proof that these exist in R n relies on a definition using Dirichlet integrals which ensure (almost) conformal parameterizations of the homotopy.See the book by Lawson [28] for an overview of this result as well as several extensions to minimal area submanifolds in more general settings.
However, beyond a proof of existence, we are interested in computing such homotopies, or at least measuring their actual area, in much simpler settings such as R 2 or a surface.To this end, we restrict the input curves to be simple curves which consist of a finite number of piecewise analytic components.We also need H to be continuous and piecewise differentiable, so that the integral can be defined.Finally, we will also require that at any time t, the intermediate curve H t is regular (see [37] for smooth curves and [31] for piecewise-linear curves).Intuitively, this means that the deformation is "kink"-free [31], and cannot create or destroy a local loop as shown in the right figure (the singular point in the right curve is a kink).Note that this is required for the minimum homotopy to even exist; again we refer the reader to the book by Lawson [28] for details.
Decomposing arrangements.Consider two simple piecewise analytic curves P and Q with the same endpoints.Their concatenation forms a (not necessarily simple) closed curve denoted by C = P • rev(Q), where rev(Q) is the reversal of Q.Let Arr(C) denote the arrangement formed by C, where vertices in Arr(C) are the intersection points between P and Q.An edge / arc in Arr(C) is a subcurve of either P or Q.
We give C (and thus P and Q) an arbitrary orientation.Hence we can talk about the sidedness with respect to C at a point p ∈ P .Specifically, a point x ∈ IR 2 is to the right of C at p if it is a counter-clockwise turn from the orientation of the vector px the orientation of (tangent of) C at p (see the right figure for an example).Given two oriented curves γ 1 and γ 2 , an intersection point p of them is positive if it is a counter-clockwise turn from the orientation of γ 1 to that of γ 2 at p.For a curve γ and a point x ∈ γ, the index of x is the parameter of x under the arc-length parameterization of γ.We sometimes use x to represent its index along γ when its meaning is clear from the context.Given two points x, y ∈ γ, we will use γ[x, y] to denote the unique sub-curve of γ between points x and y.
We say that a homotopy H from P to Q is right sense-preserving if for any t, s ∈ [0, 1], we have that either H t+dt (s) = H t (s) or H t+dt (s) is to the right of the oriented curve H t at H t (s).If it is the former case, then we say that p = H t (s) is a fixed point at time t.Similarly, we say that H is left sense-preserving if for any t, s ∈ [0, 1], H t (s) is either a fixed point or deforms to the left of the curve H t .Our homotopy H is sense-preserving if it is either right or left sense-preserving.The sense-preserving property means that we can continuously deform the curve P always in the same direction, without causing local folds in the regions swept.Intuitively, any optimal homotopy should have this property to some extent, which we will make more precise and prove later.

Structure of Optimal Homotopies
Given two simple curves P and Q (with the same end points) embedded on an orientable 2-manifold M , let X = {x 1 , . . ., x I } denote the set of I intersection points between them, sorted by their order along P .Given a homotopy H from P to Q, a point p ∈ M is called an anchor point with respect to H if it remains on H(t, •) = H t at all times t ∈ [0, 1].Of course not all intersection points are anchor points.However, if p is an anchor point, then it is necessarily an intersection point between P and Q, as p ∈ H 0 = P and p ∈ H 1 = Q.We exclude the beginning and ending end points of P and Q from the list of anchor points, as they remain fixed for all homotopies.In what follows, we show that any optimal homotopy can be decomposed by anchor points such that each of the resulting smaller homotopies has a simple structure.
Specifically, consider an arbitrary optimal homotopy H * .Let B = {b 1 , . . ., b k } be the set of anchor points with respect to H * , the minimum area homotopy.We order the b i 's by their indices along P .It turns out that the order of their indices along Q is the same, and the proof of this simple observation is in Appendix A. This observation implies that we can decompose H * into a list of sub-homotopies, where i is necessarily optimal, and it induces no anchor points.The following result states that an optimal homotopy without anchor points has a simple structure, which is sense-preserving.Intuitively, if any point changes its deformation direction at any moment, the deformation will sweep across some area redundantly and thus cannot be optimal.The detailed proof is in Appendix B. Lemma 3.2 If an optimal homotopy H from P to Q has no anchor points, then it is sense-preserving.

Minimum Area Homotopies In The Plane
In this section, we consider the case where the input consists of two simple polygonal curves in the plane.We develop an algorithm to compute the similarity between P and Q in O(I 2 log I + n log n) time, where n is the total complexity of input curves and I is the number of intersections.Note that I = Θ(n 2 ) in the worst case, although of course it may be much smaller in some cases.Although efficient algorithms for comparing curves in the plane exist (such as the Fréchet distance), our planar algorithm will be the fundamental component for comparing curves on general surfaces in the next section.It turns out that our approach can easily be extended to measure similarity between simple cycles in the plane; see Appendix D.

Relations to Winding Numbers
We are given two simple (open) curves in the plane which share common endpoints.Previously, we have shown that if an optimal homotopy does not induce anchor points, then it is sense-preserving.The implication of this result is manifested by using the winding number, defined for a loop γ in the plane at a base point x.
Intuitively, imagine starting from a point y on γ, and connecting x and y by a string.The winding number at x w.r.t.γ, denoted by wn(x; γ), is an integer measuring how many times this string winds, in a clockwise manner, around x as y traverses γ.More formally, consider an infinite ray f based at x which is generic (so it has a finite set of transversal intersections / crossings with γ).Consider a crossing γ(t) between the ray f and γ.This crossing is positive if the triangle x, γ(t), and γ(t + ǫ) is oriented counterclockwise, and is negative if oriented clockwise.Then wn(x; γ) is the number of positive crossings minus the number of negative crossings with respect to any generic ray from x.We say an oriented curve γ has consistent winding numbers if wn(x, γ) is either all non-negative, or all non-positive, for all x ∈ IR 2 .Note that for a curve with consistent winding numbers, we can always orient the curve appropriately so that wn(x, γ) is all non-negative.Two examples are shown in the figure on the right, where the second example has consistent winding numbers.Let Arr(γ) denote the arrangement formed by the curve γ.All points in the same cell of the arrangement Arr(γ) of γ have the same winding number, and the winding numbers of two neighboring cells differ by 1.The relation of consistent winding numbers and sense-preserving homotopies is given below, and the proof can be found in Appendix ??.Lemma 4.1 If there is a sense-preserving homotopy H from P to Q, then the closed curve P • rev(Q) has consistent winding numbers.
Proof: Without loss of generality, assume that the map H is left sense-preserving, always deforming an intermediate curve to its left.Consider the time-varying function F : [0, 1] × IR 2 → Z, where F (t, x) = wn(x; H t ) is the winding number at x ∈ IR 2 with respect to the curve parameterized by H t .Obviously, F (0, x) = wn(x; P • rev(Q)), and F (1, x) = 0.During the deformation, F (t, x) changes by either 1 or −1 whenever the intermediate curve sweep over it.Since the homotopy is left sense-preserving, when an intermediate curve sweeps x, x always moves from the left side of the intermediate curve to its right side.Hence the winding number x decreases monotonically.Since in the end, the winding number at each point is zero, wn(x; If the map H is left sense-preserving, then a symmetric argument shows that wn(x; Next, we describe two results to connect the above lemma to the computation of optimal homotopy.First, we define the total winding number Tw(γ) of a curve γ as where dν(x) is the area form 1 .The following observation is straightforward.Observation 4.2 For any P and Q in the plane, Proof: Take any regular homotopy H from P to Q.The area of a regular homotopy H in our setting can be reformulated as an integral on the image domain as where deg H (x) is defined as the number of connected components in the pre-image of x under H.In other words, deg H (x) is the number of times that any intermediate curve H t sweeps through x.Now consider the function = 0, and F is a continuous function.Furthermore, each time the winding number at a point x changes by 1 for some t ∈ [0, 1], it means that some intermediate curve H(t) sweeps through it.Hence |wn(x)| is a lower bound for deg H (x). We thus have that for any regular homotopy H, implying that Proof: We provide a sketch of the proof here to illustrate the main idea; see [?] for the full proof.We prove the claim by induction on the number of intersections between P and Q.The base case is when there is no intersection between P and Q.In this case, γ is a Jordan curve which decomposes IR 2 into two regions, one inside γ and one outside.The optimal homotopy area σ(P, Q) in this case is the area of the bounded cell.
All points in the bounded cell have winding number 1 (or −1) and the claim follows.Now assume that the claim holds for cases with at most k − 1 intersections.We next prove it for the case with k intersections.Let an X-arc denote a subcurve of curve X.Consider the arrangement Arr(γ) formed by γ = P •rev(Q).Since P and Q are simple, every cell in this arrangement has boundary edges alternating between P -arcs and Q-arcs.Assume without loss of generality that γ has all non-negative winding numbers.Consider a cell R ∈ Arr(γ) with largest (and thus positive) winding number.Since its winding number is greater than all its neighbors, it is necessary that all boundary arcs are oriented consistently as shown in Figure 1 (a), where the cell R (shaded region) lies to the left of its boundary arcs.
If R has only two boundary arcs, e from P and e ′ from Q, respectively, then we can morph P to another simple curve P ′ by deforming e through R to rev(e ′ ) as illustrated on the right.The area swept by this deformation is exactly the area of cell R. Furthermore, after the deformation, every point x ∈ R decreases their winding number by 1, and no other point changes its winding number.Since any point in this cell initially has strictly positive winding number, the resulting curve γ ′ = P ′ • Q still has all non-negative winding numbers.The number of intersections between Otherwise, the cell R has more than one P -arc.Take the P -arc e 1 with the smallest index along P , and let p be its ending point.Let e 2 be the next P -arc along the boundary of R, and q its starting point, and Q[q, p] the Q-arc between e 1 and e 2 , denoted by ē in Figure 1.P [p, q] and −Q[p, q] bound a simple polygon, which we denote by Ω.Since Ω does not intersect R, either Ω is on the opposite side of the Q-arc ē from the interior of R (Figure 1 (b)), or they are on the same side (Figure 1 (d)).It turns out that in both cases, we can deform P to another simple curve P ′ such that (i) the number of intersections is reduced by 2, and (ii) P ′ • Q still has consistent winding numbers.For example, in the case of Figure 1 (b), P is then deformed to sweep the region Ω as shown in Figure 1 (c).By applying the induction hypothesis to P ′ • Q, we are able to obtain that σ(P, Q) = Tw(γ).The details are in Appendix C.

The Algorithm
Lemma 3.2 and 4.1 imply that if the closed curve P • rev(Q) produces both positive and negative winding numbers, then any optimal homotopy from P to Q must have at least one anchor point.On the other hand, if it has consistent winding numbers, then by Lemma 4.3 we can compute the optimal cost to deform them by simply computing the total winding number.This leads to a simple dynamic-programming (DP) approach to compute σ(P, Q).
Specifically, let x 0 , x 1 , . . ., x I denote the intersection points between P and Q, ordered by their indices along P , where x 0 and x I are the beginning and ending points of P and Q, respectively.Let T (i) be the cost of the optimal homotopy between P [0, x i ] and Q[0, x i ], and C[i, j] the closed curve formed by We say that a pair of indices (i, j) is valid if (1) x i and x j have the same order along P and along Q; and (2) the closed curve C[i, j] has consistent winding numbers.We have the following recursion:

Time Complexity Analysis
The main components of the DP framework described above are (i) to compute Tw(C[i, j]) for all pairs of i, js, and (ii) to check whether each pair (i, j) is valid or not.These can be done in O(I 2 n) total time in a straightforward manner.We now show how to compute them in O(I 2 log I) time after O(I log I + n log n) pre-processing time.Specifically, we describe how to compute such information in O(I log I) time for all C[r, i]s for a fixed r ∈ [1, I] and all indices i > r.
To simplify the description of the algorithm, we extend Q on both sides until infinity, and obtain Q.Now collect all intersection points between P and Q, {x 1 , . . ., xI }, which is a super-set of previous intersection points, and sort them by their order along the curve P .The algorithm can be made to work with Q directly, but using Q makes the intuition behind our algorithm, as well as its description, much more clear.
Note that Q divides the plane into two half-planes.For illustration purpose, we will draw Q as a horizontal line, and use the upper and lower half-planes to refer to these two sides of Q.Another way to see that regarding Q as a horizontal line does not cause any loss of generality is that one can always find a homeomorphism from IR 2 → IR 2 such that the image of Q is a horizontal line under this homeomorphism.Now for a fixed integer r, we traverse P starting from xr .We aim to maintain appropriate data structures so that each time we pass through an intersection point xi with Q, we can, in O(log I) time, (1) check whether (r, i) is valid, and (2) obtain total winding number for C[r, i]. with Arr(C[r, u]), regardless of where r is, only points within R u will change their winding number, either all by +1 or all by −1, depending on whether R u is to the right side or the left side of the P -arc P [x u , xu+1 ], respectively.The winding numbers for points outside R u are not affected.Hence the change in the total winding number is simply α u Area(R u ), where α u is either +1 or −1.See Figure 2, where all points in R u will decrease their winding number by 1 as we move from C[r, u] to C[r, u + 1].
We can pre-compute the area of R u 's for all u in O(n log n + I log I) time, by observing that the set of R u s satisfy the parentheses property: Namely, either R u and R v are disjoint in their interior, or one contains the other.Specifically, first, we compute the arrangement of Arr(P + Q) and the area of all cells in it in O(n log n+ I) time.Each R u is the region bounded between a P -arc P [x i , xi+1 ] and a corresponding Q-segment Q[x i , xi+1 ].Since no two P -arcs intersect, the containment relationship between such P -arcs satisfies parentheses property.In particular, we can use a collection of trees to represent the containment relation among all regions R u s.See Figure 3 for an illustration.The difference between the region represented at a parent node and the union of regions represented by all its children is a cell in Arr(P + Q).For example, the shaded cell in the right figure is the difference between R 0 and its children R 2 and R 4 .We can thus compute the area of all R u s by a bottom-up traversal of these trees.Computing these trees take O(I log I) time by first sorting all intersection points with respect to their order along Q. Traversing these trees to compute all R u s takes O(I) time.Putting everything together, we need O(n log n + I log I) time.
With the area of all R u s known, updating the total winding number from C[r, u] to C[r, u + 1] takes only constant time.
Checking the validity of (r, i)s.To check whether (r, i) is valid or not, we need to check whether all cells in the arrangement Arr(C[r, i]) have consistent winding numbers.First observe that for any r and i, Arr(P + Q) is a refinement of the arrangement Arr(C[r, i]).That is, a cell in Arr(P + Q) is always contained within some cell in Arr(C[r, i]).Hence all points within the same cell of Arr(P + Q) always have the same winding number with respect to any C[r, i]), and we simply need one point from each cell in Arr(P + Q) to maintain the winding number for all cells in Arr(C[r, i]), for any r and i.We now describe how to maintain the winding number for cells of Arr(P + Q) (thus for Arr(C[r, u])s) as we pass each u > r, so that we can check whether C[r, u] has consistent winding numbers or not efficiently.The collection of such representative points hit all cells in Arr(P + Q).(It does not matter whether there may be more than one point taken from a cell of Arr(P + Q).) Hence Arr(C[r, i]) has consistent winding number if and only if all these representative points have consistent winding numbers.Next, we build a data structure to maintain the winding numbers for these points as i increases.Specifically, let U be the set of representatives that are to the right of Q, which are the stars above Q in the right figure.(Those to the left of Q will be handled in a symmetric manner).Each point has a key associated with it which is its index along Q.We build a standard balanced 1-D range tree on U based on such keys, where each leaf f stores a point from U .Every internal node v is associated with an interval where l v and r v are the smallest and largest keys stored in the subtree rooted at v. In other words, all representatives with an index along Q within [l v , r v ] are stored in the subtree rooted at v. At every node v, interior or not, we also store a value addW v .To compute the winding number for the representative point p f stored at a leaf node f , we identify the path {v 0 , v 1 , . . ., v a = f } from the root v 0 to f .The winding number for p f is simply a i=0 addW v i .Finally, each internal node v also stores the maximum and minimum winding numbers associated with all leaves in its subtree.At the beginning, all winding numbers are zero.The size of this tree is O(I) with height O(log I), and can be built in O(I log I) time once the arrangement Arr(P + Q) is known.
Let q i denote the index of point xi along Q (or can be considered as the x-coordinate of xi ).As we move from C[r, u] to C[r, u + 1], cells of Arr(P + Q) contained in R u should either all increase or all decrease their winding number by 1.Note that representatives of these cells are simply those contained in the interval [q u , q u+1 ] (or [q u+1 , q u ] if q u+1 < q u ).Hence updating the winding number is similar to an interval query of [q u , q u+1 ], and the O(log I) number of nodes in the canonical decomposition of [q u , q u+1 ] update their addW v values by either +1 or −1 depending on the sideness of R u with respect to the arc P [x u , xu+1 ].The minimum and maximum winding numbers can also be updated O(1) time per visited node.The entire process visits O(log I) nodes, and thus takes O(log I) time as i increases from u to u + 1.To see whether C[r, u + 1] has consistent winding numbers or not, we only need to check the minimum and maximum winding numbers stored at the root of the tree, denoted by w min and w max , respectively.If w min × w max equals to zero, then all winding numbers w.r.t.C[r, u + 1] are either all non-negative or all non-positive.Otherwise, (r, u + 1) is not valid.
Repeat the above process for every r ∈  The case where we have two simple cycles P and Q in the plane is discussed in Appendix D, and we obtain the following extension: Corollary 4.5 Given two polygonal cycles P and Q in the plane of n total complexity and with I intersection points, we can compute the optimal homotopy and its area in O(I(I 2 log I + n log n)) time if I > 0; and compute the optimal homotopy area in O(n log n) time if I = 0.

Minimum Area Homotopies on 2-Manifolds
In this section, we consider optimal homotopy between curves P and Q on an orientable and triangulated 2-manifold M without boundary.Our input is a triangulation K of M with complexity N , and two simple homotopic polygonal curves P and Q sharing endpoints.Edges in P and Q are necessarily edges from the triangulation K.The total complexity of P and Q is n, and there are I number of intersections between them.Note that in this setting, I = O(n).Below we discuss separately the cases when M has non-zero genus and when M is a topological sphere.

Surfaces with Positive Genus
Given an orientable 2-manifold M , let U (M ) be a universal covering space of M with φ : U (M ) → M the associated covering map.Note that φ is continuous, surjective, and a local homeomorphism.(For full details on covering spaces, we refer the reader to topology textbooks that address this area [33]; we will also build on existing algorithmic techniques developed for the computing and working in the universal cover [21,34].) For any path γ in M , if we fix the lift (pre-image) of its starting point, then it lifts to a unique path γ in U (M ), such that φ(γ) = γ.Since P and Q are homotopic with common endpoints, the closed curve formed by C = P • Q is contractible on M , and the lift of C, denoted by C, is a closed curve in U (M ).By the Homotopy Lifting Property of the universal cover [33], we have: Observation 5.1 Once we fix the lift of the starting point of P and Q in U (M ), there is a one-to-one correspondence between homotopies between P and Q in M and those between P and Q in U (M ).
We now impose an area measure in U (M ) by lifting the area measure in M ; this can be done via the map φ, which is a local homeomorphism.Now the area of a homotopy in M is the same as the area of its lift in U (M ).As such, we can convert the problem of finding an optimal homotopy in M to finding one in U (M ).Furthermore, for any orientable compact 2-manifold with genus g > 0, its universal cover is topologically equivalent to IR 2 .Intuitively, this means that we can then apply algorithms and results from previous section to the universal covering space.
More specifically, our algorithm proceeds as follows: Step 1: Compute relevant portion of U (M ).We will construct a portion of a universal covering space U (M ) made from polygonal schema of M [35,21].Specifically, we use the algorithm from [21] to construct a reduced polygonal schema T in O(N ) time.The universal covering space consists of an infinite number of copies of this polygonal schema glued together appropriately.We call each copy of the polygonal schema in the constructed universal covering space a tile.
Recall that the universal covering space U (M ) is homeomorphic to IR 2 .We fix a lift of the starting endpoint of P and Q in U (M ) and obtain a specific lift P and Q for P and Q respectively.Since P and Q are homotopic, P and Q form a closed curve, denoted by C = P • rev( Q).Note that the number of intersection points between P and Q is at most I, as every intersection point in the lift necessarily maps to an intersection point of P and Q under φ, but not vice versa.
Consider the arrangement formed by Arr( C) in the planar domain U (M ).We will construct the portion of the universal covering space U ⊆ U (M ) which is the union of tiles that intersect or are contained inside of Arr( C).
From [21], we know that the lifted curve C passes through O(n) tiles in U (K).However, while the total number of tiles in the interior of Arr( C) is O(n) for the case where g > 1, it can be Θ(n 2 ) for the case when g = 1.Hence we will separate the case for g = 1 and g > 1, since we wish to avoid the O(n 2 ) overhead in the genus 1 case.
For the case g > 1, we use the algorithm by Dey and Schipper [21] to compute the relevant portion U of the universal covering space in O(n log g + N ) time.The output contains all O(n) copies of the polygonal schema in U , where each tile is represented by a reduced 4g-gon without being explicitly filled with triangles from K.However, recall that P and Q are curves which follow edges of the triangulation; in this construction of the polygonal schema tiles, each edge of K can be broken into O(g) pieces.So in the worst case, we must break each edge in P or Q into O(g) pieces, giving a total complexity for P and Q is O(ng) in this representation of U .Once these are known, we can compute the combinatorial structure of the arrangement of C in U , as well as the description of the set of tiles each cell in Arr( C) intersects or contains, in O(ng + I log I) time.
For the case g = 1, the input manifold is a torus, and the canonical polygon schema for it is a rectangle with oriented boundary arcs aba −1 b −1 .Imagine now that we give the base polygonal schema T 0 (which is the tile that contains the lift of the starting point of P and Q) a coordinate (0, 0); we must now assign a (0, 0) (0, −1) coordinate for every other copy of the polygon (as shown in Figure 4(a)).Specifically, a copy of polygonal schema T has coordinate (i, j) if the closed loops whose lifts start in T 0 and end in T have the same homotopy type as a i b j .We can obtain the sequence of the rectangles (and their coordinates) that the curve C will pass through in O(n + N ) time [21].Once these coordinates are known, the combinatorial structure of the arrangement of C in U can also be computed in O(n + I log I) time.Note that in this case, we have avoided explicitly enumerating the set of tiles fully enclosed within Arr( C) (the shaded tiles in Figure 4 (b)), whose number can be Θ(n 2 ) instead of O(n) when g = 1.
Step 2: Area of cells in Arr( C).In order to perform our algorithm introduced in Section 4 to the lifted curves P and Q, in addition to the combinatorial structure of Arr( C), we also need the area of each cell in Arr( C).We first describe how to compute it for the case g = 1.
Take any cell X in Arr( C) and assume the boundary of X intersects m copies of polygonal schema.Even though that X may contain Θ(m 2 ) copies of (rectangular) tiles in its interior, we do not need to enumerate these interior tiles explicitly to compute their total area.Indeed, by a scanning algorithm from left to right, we can compute in O(m) time how many tiles are completely contained inside X (heavily-shaded regions in Figure 4 (b)) (note that the coordinates of each tile traversed by the boundary of X are known).Since the area of every polygonal schema is simply the total area of the input triangulation, we can compute the total area of tiles contained inside X in O(m) time.Now let R be the collection of tiles that intersect the boundary of X.It remains to compute the total area of R ∩ X.Call each region in T ∩ X a sub-cell, for any tile T ∈ R. Let G denote the boundary curves of the polygonal schema T .There are two types of sub-cells: the essential ones which contain at least one intersection point between P and Q as their vertices, and the non-essential ones which have no intersection; see Figure 4 (b) for examples.Note that a non-essential cell is bounded by arcs from G alternating with P/Q-arcs from P or Q, since there is no intersection of P and Q along the boundary of a non-essential cell.(Here, a P/Q-arc refers to either a P -arc or a Q-arc).
First let us consider the collection of non-essential sub-cells formed by alternating G-arcs (boundary arcs of a tile) and arcs from P and Q, and compute the area of each such non-essential sub-cells.If we plot all the P -arcs within a single tile T , no two P -arcs can intersect in this tile, since P is a simple curve.Imagine that we pick an arbitrary but fixed point on the boundary G of the polygonal schema T as the origin o.Each P -arc α subdivides T into two regions; we let T α denotes the canonical one excluding o.Note that since P is a simple curve, the set of canonical regions T α s for all P -arcs must satisfy the parenthesis property, and these regions and their areas, called canonical areas, can be computed in O(ng log n + N ) time using a data structure similar to one used in Section 4.3 to compute the area of R u s.See Figure 5  What remains is to compute the area of all essential sub-cells.Note that there are O(I) essential subcells since each contains an intersection between P and Q.Let a P Q-arc to refer to an arc that starts and ends with points on G (the boundary of the polygonal schema T ) and consists of alternating P -and Q-arcs.An essential sub-cell is either completely contained within a polygonal schema, or its boundary consists of P Q-arcs, G-arcs, P -arcs and Q-arcs where no two such arcs can be consecutive: they are separated by G-arcs.Now collect all P -arcs and Q-arcs that are involved in the boundary arcs of those essential sub-cells completely contained within a tile.Plot them within the same tile T and compute their arrangement A as well as the area for each cell in A. This can be done in O(ng log n + N ) time.Since A can have only O(I) vertices in the interior of the tile T , A contains O(I) cells.If an essential sub-cell X is completely contained within a polygonal schema, then it is a union of a set of cells from A. We can simply spend O(I) time to go through cells in A, identify those contained in X and return their total area.Hence it takes O(I 2 ) time to compute the area of all such O(I) essential sub-cells.If an essential sub-cell X has G-arcs on its boundary, then we need a slightly more complicated way to handle it.
Specifically, for all the remaining essential sub-cells, there can be O(I) number of P Q-arcs along their boundaries, denoted by L. We collect all P -arcs and Q-arcs involved in L and plot them in the same tile T and compute their arrangement Arr(L).Each P Q-arc α ∈ L divides the tile T into two regions, and we define T α to be the canonical one that excluding a specific origin o on G similar to before.T α consists of a union of cells from the arrangement Arr(L), and we can compute the area of T α in O(I) time since Arr(L) has O(I) cells.Overall, in O(I 2 ) time, we can compute the area of all T α s for all P Q-arcs α ∈ L Now take an essential sub-cell X that has s number of P -, Q-, or P Q-arcs along its boundary, denoted by α 1 , α 2 , . . ., α s .Let α 1 be the arc (which can be P -, Qor P Q-arc) whose endpoints along G spans the largest interval.Then, X can be represented as X = T α 1 − i∈[2,s] T α i , where T α i is the canonical region defined by an arc α i .Since the area of all canonical regions are known (for P -arcs or Q-arcs, we have computed their canonical areas before), X's are a can be computed in O(s) time.Computing the area of all remaining essential sub-cells thus takes O(I 2 + n) time.
Putting everything together, the total time needed to compute the area of all cells in Arr( C) is O(ng log n+ N + I 2 ) when g = 1.The case when g > 1 is similar but much simpler.Indeed, we now can afford to compute all the tiles contained within any cell of Arr( C) explicitly, as their total number is bounded by O(n) [21,34].The areas of essential and non-essential sub-cells are computed using the same algorithm as above.The total time complexity is O(ng log n + N + I 2 ).
Step 3: Putting everything together.With the combinatorial structure of Arr( C) and the area of each cell computed, we now apply the algorithm from Section 4.2 to compute the optimal homotopy in O(I 2 log I + ng log n) time in U (M ), which, by 5.1, gives the optimal homotopy between P and Q in M in the same time bound.The total time complexity for the entire algorithm is O(ng log n + I 2 log I + N ).

The Case of the Sphere
We now consider the remaining case where the input has g = 1, so M is a (topological) sphere S. All paths on S are homotopic.The universal cover of a sphere is itself, and hence is compact.However, the previous algorithm in Section 4.2 works for a domain homeomorphic to IR 2 and cannot be directly applied.We now sketch how we handle the sphere case.Missing details can be found in Appendix E. We observe that the results in Section 3 still hold.However, as the sphere is compact, the winding number is not well-defined.For example, see Figure 6, where there are two ways that the curve γ winds around the point p.In the first case, the winding number at p is 0, while in the second case, the winding number is −1.In order to use a dynamic programming framework as before to compute the optimal homotopy between P and Q, we need to develop analogs of Lemma 4.1 and 4.3 for curves on the sphere.
To this end, note that if we remove one point, say z ∈ S from the sphere S, then the resulting space S z = S − z is homeomorphic to IR 2 , and the concept of the winding number is well defined for S z .Specifically, z can be considered as the point of infinity in IR 2 .The winding number of x ∈ S z w.r.t.C and z, denoted by wn(x; z, C) (C omitted when its choice is clear), is simply the summation of signed crossing numbers for any path connecting x to z.As in the planar case, we say that C is consistent w.r.t.z if wn(x; z, C) is either non-negative, or non-positive for all x ∈ S z .Similar to before, we define the total winding number w.r.t. a base point z as Tw(C; z) = Sz wn(x; z, C)dx.Let σ(P, Q; Ω) denote the best cost to morph P to Q within domain Ω.Observation 5.2 If there is an optimal homotopy between P and Q that does not sweep through some point z ∈ S, then we have σ(P, Q; S) = σ(P, Q; S z ).
Observation 5.3 Suppose H * is an optimal homotopy between P and Q with no anchor points.For any cell R in Arr(P + Q), if H * sweeps through one point in its interior, then it sweeps through all points in R.
The simple proof for the above observation is in Appendix E.1.The key result is the following lemma, the proof of which can be found in Appendix E.2.

Lemma 5.4
If there is an optimal homotopy H * of P and Q with no anchor point, then the image of this optimal homotopy cannot cover all points in S.
Given two homotopic paths P and Q from S sharing common endpoints, Lemma 5.4 and Observation 5.2 imply that if P can be morphed to Q optimally without anchor points, then there exists some point z ∈ S such that σ(P, Q; S) = σ(P, Q; S z ).For this choice of z, it is necessary that the closed curve P • Q has consistent winding numbers in S z .Once this z is identified, σ(P, Q; S z ) is simply the total winding number of P • Q w.r.t.z, as suggested by Lemma 4.3, because S z is homeomorphic to the plane.Furthermore, by Observation 5.3, we only need to pick one point from each cell of Arr(P + Q) to check for the potential z.Specifically, let {z 1 , . . ., z l } be a set of such representatives, where l = O(I).The optimal homotopy area σ(P, Q) is simply the smallest of all Tw(P • rev(Q); z i ) for those z i s with respect to whom the curve P • rev(Q) has consistent winding numbers.Hence if we assume that if there is an optimal homotopy between P and Q with no anchor points, then we can compute σ(P, Q).

Overview of the algorithm for sphere case.
To compute the optimal homotopy between P and Q, we follow the same dynamic programming framework as before.If there is no anchor point in an optimal homotopy, then we use the discussion above to compute the optimal homotopy area.Otherwise, we identify the intersection point that serves as next anchor point, and recurse.The main difference lies in the component of computing σ(i, j) := σ(P [x i , x j ], Q[x i , x j ]), assuming that there is an optimal homotopy from with no anchor points.Previously, this is done by checking whether P ′ • Q ′ has consistent winding numbers.Now, we need to check the same condition but against l = O(I) number of potential representatives {z 1 , . . ., z l } as the potential point of infinity.This gives a linear-factor blow-up in the time complexity compared to the algorithm for the planar case.However, we show that this linear blowup can be tamed down and we can again compute all σ(r, j)s for all rs and all j > r in O(I 2 log I) time, after O(n log n + N ) pre-processing time.See Appendix E.3 for details.Overall, the total time complexity remains the same as before.
Putting both cases (g > 0 and g = 0) together, we conclude with the following main result.
Theorem 5.5 Given a triangulation K of an orientable compact 2-manifold M with genus g, let N be the complexity of K. Given two homotopic paths P and Q of total complexity n with I intersection points, we can compute an optimal homotopy and its area σ(P, Q; M ) in O(I 2 log I + ng log n + N ) time.

Conclusion
In this paper, we propose a new curve similarity measure which captures how hard it is to deform from one curve to the other based on the amount of total area swept.It is robust to noise (as it is area-based), and can be computed efficiently; to our knowledge, there is no other efficiently computable similarity measure for curves on surfaces.Our algorithm can be extended for cycles in the plane (see Appendix D).It appears that our algorithm can also be extended to cycles on the surfaces.Indeed, if the optimal free homotopy has an anchor point, then we can break cycles into curves that share a common start and end point, which then reduces to the problem of comparing curves on surfaces.However, on a surface the analog of Lemma D.1 no longer holds, so that two curves may intersect in M but not have an anchor point in the optimal homotopy; in this case, it is not immediately clear how to bound the size of the universal cover necessary for our algorithm.We leave the problem of working out these details, as well as improving the time complexity for comparing cycles, as an immediate future direction.
Currently, we assume that two input paths are simple paths which share starting and ending points, which makes it easier to define homotopy equivalence.This leads to two natural questions, namely how to handle curves which do not share endpoints and how to deal with non-simple curves.Another interesting problem is to compute optimal isotopy area where we require that any intermediate curve during the deformation is also simple.
Measuring similarity of curves on surfaces is an interesting problem, and many open areas remain.Geodesic Fréchet-based measures ignore the topological constraints of underlying surface, while the homotopic Fréchet distance, homotopy height, and our method require identification of a homotopy which optimizes some cost.As far as other measures of similarity which may be tractable, one interesting new idea would be to develop an area-based curve similarity measure that allows topological changes, such as allowing a region to be swept as long as it has trivial homology.Other directions include developing efficient curve simplification algorithms based on this measure, and studying similarity between curves from more general simplicial complexes than considered in this paper (such as a manifold with boundary or holes, or non-manifolds).
Figure 7: (a) and (b) p is fixed from γ to γ + , but not so in γ ++ .Sweeping γ to γ ++ directly through the shaded region in (c) has a smaller area than first to γ + then to γ ++ (see shaded region in (d)).The darker shaded region in (d) is swept twice.(e) If the deformation changes orientation at γ, then there is a local fold in the regions swept.twice.Hence this homotopy cannot have minimal area, since we can create a smaller one which stops at time t and moves directly to some intermediate curve at a time greater than t + without sweeping any portion twice.See Figure 7 (e).Now suppose that some portion of γ is deforming to one direction and another is morphing in the opposite direction.Since H is a homotopy and is therefore continuous, this means that there is at least one interval of fixed points between these two regions (which may possibly consist of a single point).Let p be an extremal point on this interval; see Figure 7 (a) for a picture when p is the only fixed point.In addition, since p is a fixed point but not an anchor point, where know there is some t + = t + dt where p is still on h t + = γ + and another t ++ = t + + dt where p is not on h t ++ = γ ++ .Now we have several cases to consider.First, consider if H has directly reversed the direction of either portion of the curve (before p or after p), we are in a similar situation to the one previously discussed, since the curve goes from locally forward to locally backward (or vice versa).In this case, we again know that some area of the domain has been swept twice, which means γ ++ has been swept over once and then was returned to, so we can reduce the area swept by H by reparameterizing H to move directly to γ ++ without passing it and then reversing.(See Figure 7 (b),(c), and (d) for an illustration.) Now if neither portion directly reverses, then γ ++ must also be deforming to two different directions.Also, we know that γ ++ must intersect γ at some point q = p, since we are essentially rotating around a central set of fixed points on these curves.In this case, we can again alter H to attain a smaller area swept by simply sweeping directly from γ to γ ++ ; this will reduce the area since the triangular region in the center bounded by γ, γ + , and γ ++ will be swept only one time instead of twice.
The claim thus follows, since any homotopy with no anchor points that is not sense preserving cannot be minimal.

C Proof for Lemma 4.3
We prove the claim by induction on the number of intersections between P and Q.The base case is when there is no intersection between P and Q.In this case, γ is a Jordan curve which decomposes IR 2 into two regions, one inside γ and one unbound.By orienting γ appropriately, every point in the bounded cell has winding number 1 and the claim follows.Now assume that the claim holds for cases with at most k − 1 intersections.We now prove it for the case with k intersections.Let an X-arc denote a subcurve of curve X.Consider the arrangement Arr(γ) formed by γ = P •rev(Q).Since P and Q are simple, every cell in this arrangement has boundary edges alternating between P -arcs and Q-arcs.Assume without loss of generality that γ has all non-negative winding numbers.Consider a cell R ∈ Arr(γ) with largest (and thus positive) winding number.Since its winding number is greater than all its neighbors, it is necessary that all boundary arcs are oriented consistently as shown in Figure 8 (a), where the cell R (shaded region) lies to the right of its boundary arcs.

Q
If R has only two boundary arcs, e from P and e ′ from Q, respectively, then we can morph P to another simple curve P ′ by deforming e through R to −e ′ (where '−' means reversing the orientation).See the right figure for an illustration.The area swept by this deformation is exactly the area of cell R. Furthermore, after the deformation, every point x ∈ R decreases their winding number by 1, and no other point changes its winding number.Since points in this cell initially has strictly positive winding number, the resulting curve γ ′ = P ′ • Q still has all non-negative winding number.The number of intersections between P ′ and Q is k − 2. By induction hypothesis, σ(P ′ , Q) = Tw(γ ′ ).Since Tw(γ) − Tw(γ ′ ) = Area(R), we have that Tw(γ) = σ(P ′ , Q) + Area(R).It then follows from Observation 4.2 and the fact σ(P, Q) ≤ σ(P ′ , Q) + Area(R) that σ(P, Q) = Tw(γ).
Otherwise, the cell R has more than one P -arc.Take the P -arc e 1 with the smallest index along P , and let p be the ending endpoint of it.Let e 2 be the next P -arc along the boundary of R, and q its starting endpoint, and Q[p, q] the Q-arc between e 1 and e 2 , denoted by ē in Figure 8. Obviously, the subcurve P [p, q] cannot intersect R, and P [p, q] and −Q[p, q] bound a simple polygon, which we denote by Ω.Either Ω is on the opposite side of the Q-arc ē from the interior of R (Figure 8 (b)), or they are on the same side (Figure 8 (d)).
Case (1): R and Ω are on the opposite side of ē.In this case, the region Ω is to the right of the oriented arc P [p, q].Note that P does not intersect the interior of Ω; as otherwise, P will either intersect itself or intersect ē, neither of which is possible.Hence only Q can intersect Ω.Since Q is also a simple curve, there is no vertices of Arr(γ) contained in the interior of Ω.As a result, every cell of Arr(γ) contained in Ω must have at least one boundary edge coming from P [p, q].This implies that each cell contained in Ω has strictly positive winding number; that is, wn(x; γ) > 0 for any x ∈ Ω.This is because if a cell ξ ⊆ Ω has winding number 0, then its neighbor across its boundary on the other side of P [p, q] will have winding number −1, as Ω is to the right of P [p, q].This violates the condition that γ has all non-negative winding numbers and thus cannot happen.
Case (2): R and Ω are both from the same side of ē.We now consider the remaining case as shown in Figure 8 (d).Take the unbounded region Ω := IR 2 \ Ω which is the complement of Ω.This unbounded region lies to the right of the oriented curve P [p, q].Since both P and Q are simple, only Q can intersect Ω and there is no vertices of Arr(γ) contained in the interior of Ω.First, observe that it is not possible that Ω ∩ Q = ∅.This is because otherwise, Ω is the unbounded face of Arr(γ) and thus the winding number for all points in Ω is 0. This however is not possible as this will imply that any point y to the immediate left of P [p, q] has winding number −1, violating our assumption that all cells in Arr(γ) have consistent (non-negative) winding numbers.
Hence Ω ∩ Q = ∅, and there are set of arcs from Q intersecting P [p, q].Then there must exist a bigon cell R ′ bounded by only two arcs, one P -arc P [a, b] ⊆ P [p, q] and a Q-arc β.See Figure 8 (d).Similar to the argument from the previous paragraph, we can show that points in R ′ must have strictly positive winding number.Now deform P to P ′ by sweeping P [a, b] through R ′ to β as shown in Figure 8 (e).P ′ is still simple, the number of intersection points between P ′ and Q is now k − 2. Only points in R ′ reduce their winding number by 1, and the resulting arrangement still has consistent winding numbers.As such, by induction hypothesis, we have that σ

D Cycles in the Plane
We now consider the case where we have two simple cycles P and Q in the plane.We have the following characterization: Lemma D.1 If the two simple cycles P and Q intersect, then there is an anchor point in the optimal homotopy between them.
Proof: Suppose that P and Q intersect but there is no anchor point in the optimal homotopy H * .By Lemma 3.2, we know that H * must be sense preserving.However, this means that H * continually moves one curve to the other in one local direction, which means that one curve must be entirely contained within the other, contradicting the assumption that they intersect.
At this point, the algorithm for cycles which intersect each other reduces to the one for curves: if we know which intersection point between P and Q is the anchor point, we can simply "break" the cycles into two curves at this point; this will become the start and end point for each of the two curves.Hence our algorithm for cycles will take a multiplicative factor of O(I) extra time than the algorithm for curves, since we need to try each possible intersection point as the required anchor point.
It then follows from Theorem 4.4 that: Corollary D.2 Given two polygonal cycles P and Q in the plane of n total complexity and with I > 0 intersection points, we can compute the optimal homotopy and its area in O(I(I 2 log I + n log n)) time.
The remaining case is that when the two polygonal cycles P and Q are disjoint.If one of the cycle contains the other, then the area of the optimal homotopy is simply the area sandwiched between these two simple cycles.This can be computed in O(n log n) time easily.
However, if the cycles are disjoint but neither contains the other, then in a sense the "optimal" free homotopy between them will simply be the area bounded by each curve, since the homotopy can collapse each curve separately to a point and then deform the points to each other.Indeed, that the sum of area bounded by the two Jordan cycles P and Q is the minimum possible homotopy area follows from a similar argument as the proof of Observation 4.2.However, in this case, the free homotopy described above is not regular since it collapses a curve to a single point.Nevertheless, one can argue that there exists a sequence of regular homotopies whose areas converge to this sum.In other words, the optimal area homotopy is still well-defined (as the infinum of the area of regular free homotopies between P and Q), although there does not exist a regular homotopy to achieve this optimal area.(This is analogous to similar issues that arise in the general case, and is the reason for introducing more restricted integrals in the mathematical literature [28].)E Missing Details for the Sphere Case Observation E.1 Given a closed curve Γ and any two points z, w ∈ S, we have that: wn(x; w) = wn(x; z) + wn(z; w) (all winding numbers are w.r.t the curve Γ).In particular, for any two points z 1 , z 2 from the same cell of Arr(Γ), we have that wn(x; z 1 ) = wn(x; z 2 ) for all x = z 1 , z 2 .
Proof: Let γ(x, y) be a path connecting point x to y.Note that the concatenation between γ(x, z) and γ(z, w) is a path from x to w.Since wn(x; w) is simply the summed signed crossing number of any path from x to w with respect to Γ, the claim follows immediately.

E.1 Proof for Observations 5.3
Proof of Observation 5.3 Suppose x and y are two points from the interior of R such that H * sweeps through x, but not y.Connect x with y by any path γ in the interior of R.This has to intersect the boundary of the region swept by H * , and let z be one such intersection point on γ.Obviously, there is a local fold in the optimal homotopy as it sweeps through z; namely, some intermediate curve will touch z and immediately trace back.Thus the input homotopy H * cannot be sense-preserving.Contradiction.Hence H * sweeps y as well.

E.2 Proof for Lemma 5.4
We prove the lemma by induction on the number of intersections between P and Q.When there is no intersection between P and Q (other than the common endpoints), the Jordan curve P • rev(Q) divides the sphere into two connected components, and the optimal homotopy is the smaller area of the two.The lemma holds for this base case.Now assume that the lemma holds for P and Q with at most k intersection points.We wish to show the result for the case where P and Q have k +1 intersection points.Since H * has no anchor points, this optimal homotopy is sense-preserving by Lemma 3.2.Assign an orientation to the closed curve C = P • Q so that locally, every point on the curve P will continuously deform to its right during the optimal homotopy.Now pick an arbitrary point z not on P and Q, and compute the winding number for each cell of Arr(P + Q) w.r.t.z.Take the cell R with the largest winding number.We assume that z / ∈ R. Suppose this is not the case and that z ∈ R. Then we show that we can change the choice of z to make this hold.
Specifically, if z ∈ R, then the cell R must have winding number 0. Now take the cell R ′ of Arr(P + Q) with the smallest winding number, and let w be a point from R ′ .Obviously, wn(w; z) ≤ wn(x; z) ≤ 0 for any x ∈ S. Now we consider the winding numbers w.r.t. to w instead of z.By Observation E.1 we have that wn(x; w) = wn(x; z) + wn(z; w).On the other hand, we have that wn(z; w) = −wn(w; z).Hence wn(z; w) ≥ wn(x; z) ≥ 0 for any x ∈ S. In other words, for this new choice of point w, we have that R still has the largest winding number and in this case, w / ∈ R. Similar to the proof of Lemma 4.3, the boundary of this cell consists of alternating arcs from P and from Q, and they necessarily have the orientation as shown in Figure 1 (a) (otherwise, one of the neighboring cell if R will have a larger winding number).Choose the P -arc e 1 that appears earliest along P , with p being its ending endpoint.Let P [q] be the next intersection between P and R. We have that P [p, q] and ē = Q[p, q] do not intersect each other.The Jordan curve P [p, q] • (rev(Q[p, q])) bounds two regions on the sphere (instead of a bounded one and an unbounded one in the planar case as shown in Figure 1 (b) and (d)).We consider the region that lies to the right of P [p, q] (thus left of ē), and denote it by Ω. See Figure 9 (a).Since through Ω ′ , and then optimally morph P ′ to Q. On the other hand, by the induction hypothesis, there is an optimal homotopy H ′ from P ′ to Q that does not sweep some point, say z 1 in S.There are now two cases: (i) If z 1 ∈ S − Ω ′ , then there is an optimal homotopy from P to Q that does not sweep z 1 as well.The induction step then holds and the claim follows.
(ii) Otherwise, z 1 ∈ Ω ′ .Consider the cell R ′ ∈ Arr(P ′ +Q) that contains z 1 .Note that R ′ ∩(S−Ω ′ ) = ∅, as there is no vertices of Arr(P ′ + Q) contained neither on nor inside Ω ′ .Hence R ′ must also contain some point, say z 2 , that is outside of Ω ′ .It then follows from Observation 5.3 that z 2 is not swept either.This leads us back to case (i), and the induction step again holds.
The claim then follows by induction.

E.3 Details of Algorithm for Sphere Case
Dynamic programming framework.Similar to the planar case, let x 0 , . . ., x I denote the intersection points between P and Q, ordered by their indices along P , with x 0 and x I being the beginning and ending points of P and Q.Let T (i) denote the optimal homotopy area between P [0, x i ] and Q[0, x i ], and C[i, j] the closed curve formed by P [x i , x j ] • Q[x j , x i ].However, now we say that a pair of indices (i, j) is valid as long as x i and x j have the same order along P and along Q.This is different from the definition of valid pairs of indices as in the planar case.Specifically, for any closed curve γ, it turns out that γ can always have consistent winding number for some choices of the point of infinity z: that is, there always exists z ∈ S such that wn(x, z; γ) is consistent for all x ∈ S z .We call such choices of z consistent representatives w.r.t.γ.Let Tw * (γ) := min z |Tw(γ; S z )| where z ranges over all possible choices of consistent representatives w.r.t.γ.Then, Lemma 4.3 and 5.4 imply that if there is an optimal homotopy between P [x i , x j ] and Q[x i , x j ] with no anchor points, then the optimal homotopy area is Tw * (γ).However, different from the planar case, Tw * (C[i, j]) is defined for all valid pairs of i, js, and it may not in general be the optimal homotopy area for P [x i , x j ] and Q[x i , x j ].We now have the following recurrence: 0, if i == 0 min j<i and (j, i) is valid { Tw * (C[j, i]) + T (j) }, otherwise intersection point x i .To this end, we use the same range tree data structure as in Section 4.2.Specifically, we use this data structure to maintain the winding number information w.r.t a fixed based point z 1 .The wn-min and wn-max indices can be easily maintained by storing at each internal node the minimum and maximum winding number within its subtree.We can also maintain the total winding number w.r.t. the base point z 1 .The time complexity for updates remains the same as before (i.e,O(log I) time per update).
The remaining question is to compute the valid total winding numbers as i increases.Let A denote the total area of topological sphere S. First, observe that for a fixec curve C, by Observation E.1, we have wn(x; z, C) = wn(x; z 1 , C) − wn(z; z 1 , C) with respect the fixed based point z 1 .Hence we have that: ( Assume that i 1 is the wn-min index and i 2 is the wn-max index.Hence we can compute Tw(C; z i 1 ) and Tw(C; z i 2 ) in O(1) time using Eqn (1), since Tw(C; z 1 ), wn(z i 1 ; z 1 , C) and wn(z i 2 ; z 2 , C) are all maintained as C changes from C[r, u] to C[r, u + 1].
Putting everything together, with O(n log n + N ) pre-processing time, we can compute all σ(r, j)s for all j > r for any fixed r, in O(I log I) and thus computing all Tw * (C[r, u])s for all r ∈ [1, I] and all u < r, in O(I 2 log I) total time.Putting everything together, the dynamic programming problem can be solved in O(n log n + I 2 log I + N ) total time.

Observation 3 . 1
The order of b i 's along P and along Q are the same.

Figure 1 :Lemma 4 . 3
Figure 1: (a) The cell R with highest positive winding number.Its boundary consists of alternating P -arcs (red) and Q-arcs (green).The two cases of relations between P [p, q] and R are shown in (b) and (d), respectively.For case (b), we can deform P to sweep through Ω as shown in (c), and reduce the number of intersections by 2.

Figure 2 :
Figure 2: Illustration of the regions R u s.

Figure 3 :
Figure 3: The containment relations of all R u regions can be represented as a forest structure on the right.
To this end, take four points around each intersection point xi of P and Q (shown as stars in the right figure).

[ 1 ,
I]. Overall, after O((n + I) log n) pre-processing, we can check whether (r, i) is valid or not and compute Tw(C[r, i]) for all r ∈ [1, I] and all i > r in O(I 2 log I) time.Putting everything together, we have the following result.

Theorem 4 . 4
Given two simple polygonal chains P and Q (with the same endpoints) in the plane of n total complexity, and with I intersection points between them, we can compute the optimal homotopy and its area in O(I 2 log I + n log n) time and O(I 2 + n) space.

Figure 4 :
Figure 4: (a) A combinatorial view of the universal covering space U(M ). a and b are the generators and we can give each cell a coordinate.(b)The lift of P (solid curve) and the lift of Q (dashed curve).The heavily shaded region are copies of polygonal schema contained inside cells of Arr( C), and their total number can be easily computed by a scanning algorithm.R 1 is an essential cell; R 2 and R 3 are two non-essential cells.

Figure 5 :
Figure 5: (a)  We overlay all non-essential sub-cells involving P -arcs into one copy of the polygonal schema.(b) An example of the canonical region T α is shown for arc α (shaded region in the top-right corner).The shaded region in the middle is a sub-cell X which can be computed asT β 0 − T β 1 − T β 2 − T β 3, where β i s are the boundary P -and Q-arcs for X.Among these P /Q-arcs, β 0 is the top-most arc in the containment relation.

Figure 6 :
Figure 6: Two ways of sweeping a curve on sphere from base point p.

Figure 8 :
Figure 8: (a) The cell R with highest positive winding number.It boundary consists of alternating P -arcs (red) andQ-arcs (green).The two cases of relations between P [p, q] and R are shown in (b) and (d), respectively.For case (b), we can deform P to sweep through Ω as shown in (c), and reduce the number of intersections by 2. Similarly, for case (d), we can identify any bigon R ′ and deform P to reduce the number of intersections by 2 as well.