Erin Wolf Chambers
email: echambe5@slu.edu

Colin De Verdière
email: eric.colin.de.verdiere@ens.fr

Jeff Erickson
email: jeffe@cs.uiuc.edu

Sylvain Lazard
email: lazard@loria.fr

Francis Lazarus
email: francis.lazarus@gipsa-lab.inpg.fr

Shripad Thite
email: shripad@caltech.edu

Homotopic Fréchet Distance Between Curves or, Walking Your Dog in the Woods in Polynomial Time

Keywords: Homotopy, Similarity of curves, Metric space, Homotopic Fréchet distance, Geodesic leash map, Punctured plane

The Fréchet distance between two curves in the plane is the minimum length of a leash that allows a dog and its owner to walk along their respective curves, from one end to the other, without backtracking. We propose a natural extension of Fréchet distance to more general metric spaces, which requires the leash itself to move continuously over time. For example, for curves in the punctured plane, the leash cannot pass through or jump over the obstacles ("trees"). We describe a polynomial-time algorithm to compute the homotopic Fréchet distance between two given polygonal curves in the plane minus a given set of polygonal obstacles.

Introduction

Given two input curves, it is natural to ask how similar they are to each other. One common measure of curve similarity is the Hausdorff distance, which is the maximum distance between a point on one curve and its nearest neighbor on the other curve. While the Hausdorff metric does measure closeness in space, it does not take into account the flow of the curves, which is important for many applications, such as morphing in computer graphics.

The Fréchet distance between two curves, sometimes also called the dog-leash distance, is defined as the minimum length of a leash required to connect a dog and its owner as they walk without backtracking along their respective curves from one endpoint to the other. The Fréchet metric takes the flow of the two curves into account; the pairs of points whose distance contributes to the Fréchet distance sweep continuously along their respective curves. This property makes the Fréchet distance a better measure of similarity for curves than alternatives for arbitrary point sets such as Hausdorff distance. It is possible for two curves to have small Hausdorff distance but large Fréchet distance. Fréchet distance is used in many different applications; see [AG95, AB05, AHPK + 06, SKB07] and the references therein.

When the two curves are embedded in a more complex metric space, such as a polyhedral terrain or some Euclidean space with obstacles, the distance between two points on the curves is most naturally defined as the length of the shortest path between them. Variations on the resulting geodesic Fréchet distance have been studied by Efrat et al. [EGHP + 02], Maheshwari and Yi [START_REF] Maheshwari | On computing Fréchet distance of two paths on a convex polyhedron[END_REF], and more recently Cook and Wenk [START_REF] Atlas | Geodesic Fréchet and Hausdorff distance inside a simple polygon[END_REF][START_REF] Atlas | Geodesic Fréchet distance inside a simple polygon[END_REF][START_REF] Atlas | Geodesic Fréchet distance with polygonal obstacles[END_REF]. The definition of geodesic Fréchet distance allows the leash to switch discontinuously, without penalty, from one side of an obstacle or a mountain to another.

In this paper, we introduce a continuity requirement on the motion of the leash. We require that the leash cannot switch discontinuously from one position to another; in particular, the leash cannot jump over obstacles, and can sweep over a mountain only if it is long enough. We define the homotopic Fréchet distance between two curves as the Fréchet distance with this additional continuity requirement. Our continuity requirement is satisfied automatically for curves inside a simple polygon [CW07, CW08a, EGHP + 02], but not in more general environments like convex polyhedra [START_REF] Maheshwari | On computing Fréchet distance of two paths on a convex polyhedron[END_REF] or the plane with obstacles [START_REF] Atlas | Geodesic Fréchet distance with polygonal obstacles[END_REF].

The motion of the leash defines a correspondence between the two curves that can be used to morph between the two curves-two points joined by a leash morph into each other [EGHP + 02]. Thus, the homotopic Fréchet distance can be thought of as the minimal amount of deformation needed to transform one curve into the other.

Efficiently computing the homotopic Fréchet distance in general metric spaces is a new open problem. We present a polynomial-time algorithm for a special case of this problem, which is to compute the homotopic Fréchet distance between two polygonal curves in the plane minus a set of polygonal obstacles.

The current paper is structured as follows. In Section 2, we give formal definitions of leash maps, homotopic Fréchet distance, relative homotopy classes and related notions, and then describe some relevant preliminary results in Section 3. In Section 4, we present an algorithm that enumerates a finite set of relative homotopy classes of leashes, such that the homotopic Fréchet distance is realized by a leash within one of these classes. We describe an algorithm to compute the homotopic Fréchet distance between two curves in Section 5. In Section 6, we describe extensions of our algorithm to closed curves and to generalizations of homotopic Fréchet distance. Finally, we conclude by suggesting several open problems.

Definitions

Let S be a fixed Hausdorff metric space. A curve or path in S is a continuous function from the unit interval [0, 1] to S. We will sometimes abuse notation by using the same symbol to denote a curve A : [0, 1] → S and its image in S. A reparameterization of [0, 1] is a continuous, non-decreasing, surjection α : [0, 1] → [0, 1]. A reparameterization of a curve A : [0, 1] → S is any curve A • α, where α is a reparameterization of [0, 1]. The length of any curve A, denoted len(A), is defined by the metric of S; in particular, two reparameterizations of the same curve are considered to have the same length.

Given two parameters s and t, an (s, t)-leash between two curves A and B is another curve λ : [0, 1] → S such that λ(0) = A(s) and λ(1) = B(t). A leash is an (s, t)-leash for some parameters s and t. If either A or B intersects itself, two distinct leashes may be equal as curves while corresponding to different parameters s and t.

A leash map is a continuous function : [0, 1] 2 → S such that (•, 0) is a reparameterization of A, and (•, 1) is a reparameterization of B. A leash map describes the continuous motion of a leash between a dog walking along A and its owner walking along B; the curve (t, •) is the leash at time t. The length of a leash map , denoted len( ), is the maximum length of any leash (t, •). Finally, the homotopic Fréchet distance between two curves A and B, denoted F(A, B), is the infimum, over all leash maps between A and B, of the length of :

F(A, B) := inf leash map : [0,1] 2 →S max 0 t 1 len( (t, •)) .
In contrast, the classical [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF] and geodesic [START_REF] Atlas | Geodesic Fréchet distance with polygonal obstacles[END_REF] Fréchet distances (for which the leash is not required to move continuously) are defined only in terms of reparameterizations and distances:

F(A, B) := inf reparameterizations α,β : [0,1]→[0,1] max 0 t 1 dist(A(α(t)), B(β(t))) ,
where dist(a, b) denotes the distance between points a and b in the ambient metric space.

In spaces where shortest paths vary continuously as their endpoints move, such as the Euclidean plane or the interior of a simple polygon, the two notions F and F are equivalent. In general, however, the homotopic Fréchet distance F(A, B) between two curves A and B could be larger (but never smaller) than the classical Fréchet distance F(A, B) between them.

Leash maps are closely related to the standard topological notion of homotopy. Two curves λ and λ with the same endpoints are homotopic if λ can be continuously deformed into λ without moving the endpoints. More formally, λ and λ are homotopic if there is a continuous function h : [0, 1] 2 → S such that h(u, 0) = λ(u) and h(u, 1) = λ (u) for all u ∈ [0, 1], and h(0, v) = λ(0) = λ (0) and h(1, v) = λ (1) = λ (1) for all v ∈ [0, 1]. It is easy to prove that being homotopic is an equivalence relation over the set of curves with any fixed pair of endpoints, and thus determines homotopy classes.

An (s, t)-leash λ and an (s , t )-leash λ are homotopic relative to A and B, or simply relatively homotopic, if λ can be continuously deformed into λ while keeping each endpoint of the leash on its respective curve. More formally, λ is relatively homotopic to λ if there are three continuous functions α, β :

[0, 1] → [0, 1] and h : [0, 1] 2 → S, such that α(0) = s, α(1) = s , β(0) = t, β(1) = t , and such that h(u, 0) = λ(u) and h(u, 1) = λ (u) for all u ∈ [0, 1], and h(0, v) = A(α(v)) and h(1, v) = B(β(v)) for all v ∈ [0, 1].
Again, for any fixed curves A and B, being relatively homotopic is an equivalence relation, which defines relative homotopy classes of leashes. Clearly, all leashes (t, •) determined by a single leash map are relatively homotopic.

Preliminaries

In this paper, we develop a polynomial-time algorithm to compute the homotopic Fréchet distance between two polygonal curves A and B in the Euclidean plane E 2 minus a set P of polygonal obstacles. In most of the paper, in order to avoid some technicalities, we assume that the obstacles in P are open sets; however, we also consider the special case of point obstacles in Section 4.4. To simplify our exposition, we assume that no three vertices of the input (vertices of polygons in P or vertices of A and B) are collinear; this assumption can be enforced algorithmically using standard perturbation techniques [START_REF] Seidel | The nature and meaning of perturbations in geometric computing[END_REF]. Figure 1 illustrates leash maps for a few sample inputs where P is a set of very small polygonal obstacles.

Let E denote the space E 2 \ P with the metric defined by shortest path distances. In any leash map between curves A and B in E, the moving leash can neither intersect nor jump over any obstacle in P . Curves A and B may self-intersect and intersect each other. For ease of exposition, we will assume that the closures of the obstacle polygons in P are disjoint from each other and from the curves A and B; however, our algorithms can be easily adapted to avoid this restriction.

Let a 0 , a 1 , . . . , a m denote the ordered sequence of vertices of A; these points define a unique parameterization A : [0, 1] → E whose restriction to any range of the form [(i -1)/m, i/m] is an affine map onto the corresponding edge a i-1 a i . Similarly, the vertices b 0 , b 1 , . . . , b n of B define a unique piecewise-affine parameterization B : [0, 1] → E. Finally, let k denote the total number of vertices in all obstacle polygons.

Universal Cover

Given a topological space S, its universal cover S is a simply connected topological space that locally resembles S, but is (usually) infinitely larger. Each point x in S corresponds (in general) to an infinite number of points in the universal cover S, one for each homotopy class of curves with both endpoints equal to x. Universal covers are a fundamental concept in topology and using this standard tool considerably simplifies many of our proofs, especially in Section 4.3.

More formally, a continuous function p : S → S is a covering map if every point x ∈ S has an open neighborhood U such that p -1 (U ) is the union of disjoint open sets i V i , and the restriction of p to each open set V i is a homeomorphism from V i to U [START_REF] Munkres | Topology[END_REF]. If there is a covering map from S to S, then S is called a covering space of S. A point x in S is called a lift of its image p( x) in S; similarly, a path α in S is a lift of its image p( α) in S. Unless the covering map p is a homeomorphism, each point and path in S has several lifts in S. The universal cover S is the unique simply connected covering space of S. Two paths α and β in S are homotopic (with fixed endpoints) if and only if they have lifts α and β with the same endpoints in the universal cover S. The universal cover S naturally inherits the metric properties of S; for any path π in S, its length is defined as the length of its projection p • π in S, where p : S → S is the covering map. In particular, a path in S is as short as possible in its homotopy class (with fixed endpoints) if and only if it lifts to a globally shortest path in S. For further details, see Munkres [START_REF] Munkres | Topology[END_REF]. We sketch an equivalent constructive definition of the universal cover of a large bounded subset of E, due to Hershberger and Snoeyink [START_REF] Hershberger | Computing minimum length paths of a given homotopy class[END_REF]; see Figure 2. Let R be a large bounding rectangle strictly containing the obstacles in P and let S = R \ P . (The homotopic Fréchet distance within S is equal to the homotopic Fréchet distance in E.) Let ∆ be any triangulation of S whose vertices are exactly those of P and R and whose set of edges contains the edges of P and R. For example, we can take ∆ to be the constrained Delaunay triangulation [START_REF] Chew | Constrained Delaunay triangulations[END_REF] of the edges of S, minus the triangles inside the obstacles. A triangulation ∆ of the universal cover S of S is obtained by repeatedly gluing distinct copies of (closed) triangles in ∆ along their common edges ad infinitum. Thus, the dual graph of the resulting triangulation ∆ is an infinite tree. For example, in Figure 2, the dual graph of the lifted triangulation ∆ is an infinite path.

Although universal covers are a convenient tool for proving our results, we emphasize that our algorithm never explicitly constructs the universal cover.

Geodesics

A geodesic in a metric space S is a path that is locally as short as possible. More formally, a geodesic is a path π : [0, 1] → S such that for every parameter t ∈ [0, 1], the restriction of π to a sufficiently small neighborhood of t is a globally shortest path. A path in E is a geodesic if and only if every lift in E is also a geodesic.

The fact that the dual graph of the triangulation of E is a tree has important consequences for geodesics and shortest paths in E, as noted by Hershberger and Snoeyink [START_REF] Hershberger | Computing minimum length paths of a given homotopy class[END_REF] and other authors [GS98, CLMS04, EKL06, Bes03, Bes04]. Specifically, because S is simply connected and locally Euclidean, shortest paths between points in E are unique; indeed, every geodesic in E is a globally shortest path. Shortest paths in E are piecewise linear curves whose internal vertices are lifts of obstacle vertices. Furthermore, shortest paths in E vary continuously as the endpoints move continuously.

Hershberger and Snoeyink [START_REF] Hershberger | Computing minimum length paths of a given homotopy class[END_REF] describe how to algorithmically maintain a shortest path σ in E as the endpoints move continuously, by storing the sequence of lifted obstacle vertices that lie on σ in a double-ended queue or deque. These obstacle vertices partition σ into a sequence of straight line segments. If the first or last segment of σ collides with a lifted obstacle vertex, that vertex is pushed onto the appropriate end of the deque. Conversely, if the first two or last two segments of σ become collinear, the obstacle vertex joining those two segments is removed from the appropriate end of the deque. Except at these critical events, only the first and last segments of the geodesic change as the endpoints of the geodesic move.

Geodesic Leash Maps

A geodesic leash map is a leash map : [0, 1] × [0, 1] → E in which every leash (t, •) is a geodesic. We next prove that for any leash map , there is a geodesic leash map in the same relative homotopy class that is no longer than . Lemma 3.1. Suppose is a leash map between two curves A and B. There is a geodesic leash map between A and B such that, for all t ∈ [0, 1], the leash (t, •) is the shortest path homotopic to (t, •) with fixed endpoints. Additionally, the length of is at most the length of .

Proof. We lift to the universal cover E of E, obtaining a leash map between the lifts A and B of A and B respectively. For each t ∈ [0, 1], let (t, •) be the globally shortest path between the endpoints of (t, •). Because shortest paths in E vary continuously as their endpoints move continuously, is a continuous function in both arguments, and therefore a (geodesic) leash map in E. The projection of back to E is a (geodesic) leash map between A and B. For each t, the leash (t, •) is the shortest path in E that is homotopic with fixed endpoints to (t, •), so len( (t, •)) len( (t, •)). It follows that len( ) len( ). This lemma implies that the homotopic Fréchet distance F(A, B) is the infimum, over all relative homotopy classes h, of the classical Fréchet distance, where distances are defined by shortest paths in h:

F h (A, B) := inf reparameterizations α,β : [0,1]→[0,1] max t∈[0,1] dist h (A(α(t)), B(β(t))) F(A, B) := inf relative homotopy class h F h (A, B).
Here, dist h (u, v) denotes the length of the shortest path from u to v in relative homotopy class h.

For the rest of the paper, we restrict our attention to geodesic leashes and geodesic leash maps. We call a relative homotopy class h optimal if F(A, B) = F h (A, B). In Section 4, we show that there is at least one optimal relative homotopy class. We also prove a structural result about optimal relative homotopy classes, which leads to a polynomial-time algorithm to enumerate a subset of relative homotopy classes, at least one of which is optimal. Section 5 describes our polynomial-time algorithm to compute the Fréchet distance within a single homotopy class. Combining these two subroutines gives us a polynomial-time algorithm to compute the homotopic Fréchet distance.

Structural Properties of Optimal Relative Homotopy Classes

Minimality

For any relative homotopy class h and any parameters s, t ∈ [0, 1], let σ h (s, t) denote the shortest path in h between points A(s) and B(t). We define a partial order on relative homotopy classes as follows: For any two relative homotopy classes h and h , we write h h if and only if len(σ h (s, t)) ≤ len(σ h (s, t)) for all parameters s and t. We write h ≺ h whenever h h but h h.

Lemma 4.1. For any relative homotopy classes h and h , if h h , then

F h (A, B) F h (A, B).
Proof. Let be any leash map in relative homotopy class h : for some reparameterizations α and β of [0, 1], we have (

•, 0) = A(α(•)) and (•, 1) = B(β(•)).
Let be the geodesic leash map in relative homotopy class h defined by the same reparameterizations: (t, •) = σ h (α(t), β(t)) for all t (Lemma 3.1). The definition of implies that len( (t, •)) len( (t, •)) for all t; hence F h (A, B) len( ) len( ). But F h (A, B) is the infimum of all such len( ); this concludes the proof.

A relative homotopy class h is minimal if h

h implies h h . In other words, h is not minimal if there is another relative homotopy class h such that h ≺ h.

Existence of Minimal Relative Homotopy Classes

Lemma 4.2. For any relative homotopy class h, there is a minimal relative homotopy class h such that h h.

Proof. Assume, for the sake of contradiction, that there is no minimal relative homotopy class h such that h h. In other words, for any h h (including h = h), h is not minimal, so there is another relative homotopy class h such that h ≺ h h. Then, by induction, we can define an infinite descending chain of relative homotopy classes

h = h 0 h 1 h 2 • • • . To simplify notation, let σ n = σ hn (0, 0).
Consider the ordered list of obstacle vertices on each path σ n . There are finitely many such ordered lists, because len(σ n ) len(σ 0 ) for each n. Thus, for some pair of indices i < j, the paths σ i and σ j have the same endpoints (A(0) and B(0)) and visit the same ordered list of obstacle vertices. Thus, the paths σ i and σ j are identical, which implies that their relative homotopy classes h i and h j are equal. This is a contradiction. Lemmas 4.1 and 4.2 together imply that the homotopic Fréchet distance F(A, B) is the infimum of F h (A, B) over all minimal relative homotopy classes h:

F(A, B) = inf minimal relative homotopy class h F h (A, B).
In the remainder of this section, we prove that all minimal relative homotopy classes have a special form, which implies that the number of minimal relative homotopy classes is finite. (Thus, we can finally replace the infimum in the expression above with a minimum.) We also describe how to enumerate, in polynomial time, a finite set of relative homotopy classes that contains an optimal one. Our overall strategy is to compute F h (A, B) for each such candidate homotopy class h, and to return the smallest value obtained.

Structure of Minimal Homotopy Classes

We define a direct geodesic to be a geodesic in E that is either (1) a line segment from A to B, or (2) a geodesic that consists of a line segment from A to some obstacle vertex p, a globally shortest path from p to an obstacle vertex q, and a line segment from q to B. We will prove: Proposition 4.3. Every minimal relative homotopy class contains a direct geodesic.

Let h be an arbitrary minimal relative homotopy class. Let A and B be lifts of A and B in the universal cover E, such that for all s and t, the shortest path σ h (s, t) between A(s) and B(t) is a lift of σ h (s, t). Let P denote the set of all lifts of the vertices of obstacles in P ; every point in P lies on the boundary of E. Let π h denote the intersection of all shortest paths σ h (s, t). Proposition 4.3 follows directly from the following pair of lemmas.

Lemma 4.4. If π h = ∅, then h contains a direct geodesic of type (1): a line segment.

Proof. If the shortest path σ h (0, 0) is a line segment, then the geodesic σ h (0, 0) is also a line segment, and the proof is complete. Thus, we assume that σ h (0, 0) passes through at least one vertex in P .

Let p 1 , . . . , p κ be the sequence of lifted obstacle vertices on the shortest path σ h (0, 0). (The vertices p i are distinct, although their projections back into the plane might not be.) Because π h = ∅, there is, for each i, a pair of parameters (s i , t i ) such that σ h (s i , t i ) does not pass through p i .

We consider a continuous motion of the parameter point (s, t), starting at (s, t) = (0, 0) and then moving successively to each point (s i , t i ). Specifically, we define two continuous functions s : [0, κ] → [0, 1] and t : [0, κ] → [0, 1] such that s(0) = t(0) = 0, and for any integer i, we have s(i) = s i and t(i) = t i . To simplify notation, let σ(τ ) denote the shortest path σ h (s(τ ), t(τ )).

As the parameter τ ('time') increases, vertices in P are inserted into and deleted from the deque of obstacle vertices on σ(τ ). If the deque is empty at any time τ , then the shortest path σ(τ ) is a line segment, which implies that the projected path σ(τ ) is a line segment in E, concluding the proof. Thus, we assume to the contrary that the deque is never empty. Each vertex p 1 , . . . , p κ must be deleted from the deque at least once during the motion (but may be reinserted later).

Suppose p is the last vertex among p 1 , . . . , p κ to be removed from the deque for the first time. Without loss of generality, we assume p is first removed from the front of the deque at time τ 1 . Let q denote the second vertex in the deque just before p is removed; this vertex must exist, because the deque is never empty. The vertex p lies on the first segment a q of σ(τ 1 ), where a = A(s(τ 1 )).

By definition of p, vertex q must have been pushed onto the back of in the deque at some earlier time τ 2 < τ 1 . Just before q is inserted, the last vertex in the deque must be p. Moreover, q lies on the last segment p b of σ(τ 2 ), where b = B(t(τ 2 )). Thus, there is a line segment a b between a point in A and a point in B. The projection ab of this segment into E is a line segment in homotopy class h. Lemma 4.5. If π h = ∅, then h contains a direct geodesic of type (2): the concatenation of a line segment from A to some obstacle vertex p, a globally shortest path from p to an obstacle vertex q, and a line segment from q to B.

Proof. The path π h is a shortest path between some pair of lifted obstacle vertices p and q. (In the special case where π h is a single point, we have p = q = π h .) Now p and q are lifts of obstacle vertices p and q (which may be the same point, even if p and q are not), and π h is similarly a lift of some path π h with endpoints p and q.

Let σ(p, q) denote a globally shortest path from p to q, and suppose that it is strictly shorter than π h . For any parameters s and t, let τ (s, t) denote the curve obtained from σ h (s, t) by replacing the subpath π h with σ(p, q). All paths τ (s, t) belong to the same relative homotopy class, which we denote h . We now easily confirm that h ≺ h, contradicting our assumption that h is minimal. We conclude that π h is the shortest path from p to q.

It remains to show that there is a line segment from some point on A to p. (A similar argument implies that there is a line segment from q to some point on B.) For all s and t, the geodesic σ h (s, t) is the concatenation of a geodesic α(s) from A to p, the shortest path π h , and a geodesic β(t) from q to B. If α(0) is a line segment, our claim is proved. Thus, we assume that α(0) is not a line segment, which implies that the lifted path α(0) passes through at least one lifted obstacle vertex other than its endpoint p. Let p -be the last lifted obstacle vertex on α(0) before p. Let s 0 be the largest value such that α(s) contains p -for all 0 s s 0 . Because p -is not on the common subpath π h , it is not on every geodesic α(s), which implies that s 0 < 1. The geodesic α(s 0 ) is a line segment.

Corollary 4.6. We can enumerate a set of O(mnk 4 ) relative homotopy classes that contains at least one optimal relative homotopy class, in O(mnk 4 ) time.

Proof. For any points a ∈ A and b ∈ B, we call the line segment ab extremal if it satisfies one of the following conditions:

(i) The endpoints are vertices of A and B.

(ii) One endpoint is a vertex of A or B and the segment contains one vertex of P .

(iii) The segment contains two vertices of P .

Every line segment in E is relatively homotopic to at least one extremal line segment in E. Thus, to enumerate the relative homotopy classes that contain a line segment, it suffices to enumerate the extremal line segments in E.

There are O(mn) extremal segments of type (i), which we can easily enumerate in O(mn) time by brute force. Each vertex a ∈ A and vertex p ∈ P determine at most n extremal segments of type (ii), one for each intersection between the ray from a through p and B.

Similarly, each vertex b ∈ B and vertex p ∈ P determine at most m extremal segments of type (ii). Thus, there are O(mnk) extremal segments of type (ii); again, we can easily enumerate these in O(mnk) time. Finally, any two vertices p, q ∈ P determine O(mn) extremal segments of type (iii), distinguished by the intersection points of the line through p and q with A and B, so there are O(mnk 2 ) type-(iii) extremal segments in total. For any obstacle vertices p and q, we can compute the intersection points between the line through p and q and A or B in O(m + n) time, and then enumerate the extremal segments that contain both p and q in O(mn) time, again by brute force.

Altogether, we enumerate O(mnk 2 ) extremal line segments in O(mnk 2 ) time. To build all extremal line segments in E, we discard any line segment that intersects any obstacle polygon; this takes O(mnk 3 ) time in total.

To enumerate all other direct geodesics (of type (2)), we begin by computing shortest paths between every pair of obstacle vertices [START_REF] Hershberger | An optimal algorithm for euclidean shortest paths in the plane[END_REF]. (If there is more than one shortest path between any pair of obstacles vertices, we can break ties arbitrarily.) Next, for every obstacle vertex p, we want to find all (relative homotopy classes of) line segments starting at p and ending at a point on A or on B. We compute them as follows: for every obstacle vertex q = p, we shoot a ray from p in the direction of q until it reaches the interior of an obstacle (or infinity), and then compute all O(m + n) intersections between the resulting line segment (or ray) and the curves A and B. This gives us endpoints of line segments starting at p. To this list of line segments, we also add every segment in E from p to a vertex of A or B. We now have the complete list of potential initial and final segments of direct geodesics. Finally, we concatenate all O(mk) initial segments, O(k 2 ) shortest paths, and O(nk) final segments to obtain O(mnk 4 ) paths in O(mnk 4 ) time.

Proposition 4.3 is not a complete characterization of minimal relative homotopy classes. It is easy to find direct geodesics of type (2) whose relative homotopy classes are not minimal. However, the next lemma shows that every type (1) direct geodesic determines a minimal homotopy class.

Lemma 4.7. The relative homotopy class of any line segment is minimal.

Proof. Let σ be a line segment from A(s) to B(t), and let h be the relative homotopy class of σ. For any relative homotopy class h = h, the shortest path σ h (s, t) must be longer than σ = σ h (s, t), which implies that h h. We conclude that h is minimal.

There are input instances that admit Ω(mnk 2 ) distinct minimal relative homotopy classes. For example, Figure 3 shows such an example with k/3 triangular obstacles. If the triangles are sufficiently small, the line through any two obstacle vertices intersects a constant fraction of the edges of both A and B, defining Ω(mnk 2 ) type-(iii) extremal line segments. Lemma 4.7 implies that the homotopy classes of these extremal line segments are minimal. A constant fraction of these minimal homotopy classes contain at most four extremal segments.

Moreover, this example admits Ω(mnk 4 ) relative homotopy classes of type-(2) direct geodesics. Consider the direct geodesics whose first and last obstacle vertices lie on the convex hull of the obstacles. There are Ω(k 2 ) choices for the first and last obstacle vertices; for each such choice, there are Ω(mk) choices for the initial line segment and Ω(nk) choices for the final line segment. Thus, any improvement in this portion of the algorithm will require a finer characterization of minimal relative homotopy classes.

Point Obstacles

Our previous structural results also apply to degenerate obstacles, such as points or line segments, with little modification, by replacing them with sufficiently small or thin triangles. Corollary 4.6 then implies a bound of O(mnk 4 ) on the number of minimal homotopy classes for such inputs. The goal of this section is to provide a complete characterization of the minimal homotopy classes when all the obstacles are points, which yields a better bound of O(mnk 2 ) on their number.

Thus, in this section, let P be a set of obstacle points in general position. Because the obstacles are now closed sets, there are pairs of points in E = E 2 \ P that have no shortest path between them; more generally, there are homotopy classes of paths in E that contain no geodesics. In this setting, the distance between two points a and b (within any homotopy class) is properly defined as the infimum of the lengths of all paths (in that homotopy class) from a to b. For simplicity of exposition (and computation), we extend the definition of 'geodesic' to include any path in E 2 that arises as the limit of a converging sequence of paths in E in the same homotopy class (with fixed endpoints), whose lengths converge to the distance between the endpoints within that homotopy class. Geometrically, geodesics in E are now polygonal paths in E 2 whose internal vertices are obstacle points. (This extension is implicit in the works of Efrat et al. [START_REF] Efrat | Computing homotopic shortest paths efficiently[END_REF] and Bespamyatnikh [START_REF] Bespamyatnikh | Computing homotopic shortest paths in the plane[END_REF], who describe algorithms to compute 'shortest' paths homotopic to a given path, in the plane minus a set of points.)

However, in order to uniquely identify the relative homotopy class of a geodesic, some additional information is now required in addition to its geometry. Specifically, we associate a turning angle with each obstacle point that the geodesic touches. Consider a geodesic γ that passes through an obstacle point p. Let C ε be a circle centered at p with radius ε > 0, small enough to exclude every other obstacle in P . A turning angle of θ at an obstacle point p indicates that replacing the portion of γ inside C ε with an arc of length ε|θ| around C ε , which goes counterclockwise around C ε if θ > 0 and clockwise if θ < 0, yields a new path homotopic to γ. See Figure 4. A path could meet the same obstacle point more than once; we associate a different turning angle with each incidence. If γ is a geodesic, none of its turning angles is in the range (-π, π), since otherwise γ could be locally shortened.

We emphasize that Proposition 4.8 does not imply that the optimal leash map contains a line segment. At first glance, it may seem natural to conjecture that the optimal leash map must also contain a line segment; surprisingly, this conjecture is actually false. Lemma 4.10. There is a pair of polygonal curves and a set of point obstacles such that no optimal leash map contains a line segment.

Proof. Consider the instance shown in Figure 5(a). The vertices of A have coordinates (-2, 2), (-2, 4), (2, 4), and (2, 2), in that order; the vertices of B have coordinates (-2, -2), (-2, -4), (2, -4), and (2, -2), in that order; and the obstacle points have coordinates (1, 2), (-1, 2), (-1, -2), and (1, -2). (This instance is highly degenerate, but it can easily be perturbed into general position without affecting the result.) (a) Up to rotations and reflections, there are only four relative homotopy classes of segments with one endpoint on each curve. Figure 5(b) shows one line segment (dashed) and the initial and final leashes (solid) in each relative homotopy class. As suggested by the figure, we call these four classes I, Z, L, and C. We claim that class I is the only optimal relative homotopy class, and that the optimal leash map does not contain a line segment.

I Z L C (b) Õ 3 Õ 2 Õ 1 Õ 0 Õ 4 (c)
Figure 5(c) shows the first half of a leash map in class I, in which one endpoint of the leash traverses A completely before the other endpoint moves at all. The figure shows five critical leashes λ 0 , λ 1 , λ 2 , λ 3 , λ 4 ; between any two critical leashes, the length of the leash is either monotonically increasing or monotonically decreasing. The longest critical leash is λ 3 , which has length 1 + 3 √ 5 ≈ 7.708; this is also the length of the leash map. The final leashes in classes Z, L, and C have lengths 8, 4 + 2 √ 5 ≈ 8.472, and 10, respectively. Within each relative homotopy class, the length of the final leash is a lower bound on the length of any leash map. Thus, I is the unique optimal homotopy class. On the other hand, the shortest line segment in class I has length 8, which is longer than λ 3 . This completes the proof.

Non-polygonal Obstacles

Our proof of Proposition 4.3 can be extended to non-polygonal obstacles with only minor modifications; the obstacles need not be convex or have smooth boundaries. In this more general setting, the initial and final segments of a direct geodesic must be tangent to the obstacles at their endpoints; that is, these segments can be made slightly longer without intersecting any obstacle. The algorithmic results in Sections 5 and 6 similarly extend to non-polygonal objects, provided one can efficiently compute the visibility graph of the obstacles [PV95, PV96]; the running time of the resulting algorithm obviously depends on the exact representation of the objects. Further details of this extension are described by Chambers [START_REF] Chambers | Finding Interesting Topological Features[END_REF].

Computing Homotopic Fréchet Distance

Finally, we describe our algorithm to compute the homotopic Fréchet distance between two polygonal curves in the plane with polygonal obstacles. Our approach is to compute a set of relative homotopy classes that includes at least one optimal class, as described by Corollary 4.6, and then compute the Fréchet distance F h (A, B) within each homotopy class h in this set. Our algorithm to compute F h (A, B) is a direct adaptation of Alt and Godau's algorithm for computing the classical Fréchet distance between polygonal paths in the plane [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF].

Henceforth, to simplify notation, we consider that the polygonal chain A, whose ordered sequence of vertices is a 0 , a 1 , . . . , a m , is parameterized over the interval [0, m], instead of [0, 1] as in previous sections, so that A(i) = a i for each integer i between 0 and m. Similarly, we parameterize B over the interval [0, n]. As in the previous section, for any s ∈ [0, m] and t ∈ [0, n], let σ h (s, t) denote the shortest path from A(s) to B(t) in homotopy class h, and let dist h (s, t) = len(σ h (s, t)). For any ε > 0, let F ε ⊆ [0, m] × [0, n] denote the free space {(s, t) | dist h (s, t) ε}. Our goal is to compute the smallest value of ε such that F ε contains a monotone path from (0, 0) to (m, n); this is precisely the Fréchet distance F h (A, B).

The parameter space [0, m] × [0, n] decomposes naturally into an m × n grid; let i,j = [i -1, i] × [j -1, j] denote the grid cell representing paths from the ith edge of A to the jth edge of B.

Geodesic Distance Is Convex

In this section, we prove the following proposition, required by our generalization of Alt and Godau's algorithm:

Proposition 5.1. The restriction of the function dist h to any grid cell i,j is convex.

We first recall the following elementary classical properties of the Euclidean norm (denoted • ). Let α, β : [0, 1] → E be affine functions with (constant) derivatives a and b respectively, and let h be a relative homotopy class. For each t ∈ [0, 1], let σ(t) be the shortest path from α(t) to β(t) in relative homotopy class h, and let d(t) denote the length of σ(t).

Fix t ∈ [0, 1]. The shortest path σ(t) is a polygonal curve. Let u(t) be the unit vector representing the direction of the first line segment of σ(t) (at its initial point α(t)). Similarly, we denote by v(t) the unit representing the direction of the last line segment of σ(t).

Recall that, as t increases, the shortest path σ(t) encounters a finite number of events. Between every two consecutive events, the sequence of obstacle vertices at which σ(t) bends is the same.

Lemma 5.3. Between any two consecutive events, d is convex and of class C 1 . In particular,

d (t) = b • v(t) -a • u(t)
, where • denotes the inner product.

Proof. Fix two consecutive events t 0 and t 1 .

Assume first that for all t between t 0 and t 1 , the path σ(t) is not a line segment. Then for every t ∈ [t 0 , t 1 ], σ(t) is the concatenation of a line segment from α(t) to a fixed obstacle vertex p, a geodesic from p to another fixed obstacle vertex q, and the line segment from q to β(t). It follows that d(t) equals a constant plus ---→ pα(t) + ---→ qβ(t) . Our result is now a consequence of Lemma 5.2. Specifically, d is the sum of two convex functions, and is therefore convex. Since α and β do not meet obstacle vertices, the function d is C 1 in the interval [t 0 , t 1 ]. The chain rule implies the claimed expression for d . Specifically,

d dt ---→ pα(t) = d dt ϕ p (α(t)) = -→ ∇ϕ p (α(t)) • d dt α(t) = ---→ pα(t) ---→ pα(t) • a = -u(t) • a.
A similar derivation implies that

d dt ---→ qβ(t) = v(t) • b.
If σ(t) is a line segment whenever t 0 t t 1 , then d(t) = -----→ α(t)β(t) . Since the function t → -----→ α(t)β(t) is affine, Lemma 5.2 also implies that d is convex and of class C 1 , and that

d (t) = ( b -a) • -----→ α(t)β(t) -----→ α(t)β(t)
.

Finally, we observe that

u(t) = v(t) = -----→ α(t)β(t) -----→ α(t)β(t)
, which completes the proof.

Lemma 5.4. The function d is convex.

Proof. Lemma 5.3 implies that between consecutive events the function d is continuous and non-decreasing; indeed, a • u(t) = a cos θ(t) where θ(t) is the angle between a and u(t). The angle θ(t) is constant if point p and segment α([0, 1]) are collinear; otherwise, |θ(t)| is strictly increasing (from 0 to π if α([0, 1]) was an infinite line). Thus, -a • u(t) is non-decreasing; a similar argument implies that b • v(t) is non-increasing. Let t 0 be an arbitrary event. Since the functions t → u(t) and t → v(t) are continuous at t 0 , Lemma 5.3 implies that d is also continuous at t 0 . Thus, d is non-decreasing over the entire interval [0, 1], which implies that d is convex.

Decision Procedure

Like Alt and Godau, we first consider the following decision problem: Is F h (A, B) at least some given value ε? Equivalently, is there a monotone path in the free space F ε from (0, 0) to (m, n)? Our algorithm to solve this decision problem is identical to Alt and Godau's, except for the O(log k)-factor penalty for distance queries; we briefly sketch it here for completeness.

For any integers i and j, let h i,j denote the intersection of the free space F ε with the horizontal edge ([i -1, i], j), and let v i,j denote the intersection of F ε with the vertical edge (i, [j -1, j]). In the first phase of the decision procedure, we compute h i,j and v i,j for all i and j, using one distance query (and O(log k) time) for each edge of the parameter grid.

In the second phase of the decision procedure, we propagate in lexicographic order from 1,1 to m,n and determine which h i,j and v i,j are reachable via a monotone path from 1,1 . Since the free space in each i,j is convex, we can propagate through each cell in constant time.

Our decision algorithm returns true if and only if there is a monotone path that reaches (m, n). The total running time of our decision procedure is O(mn log k).

Optimization

Finally, we describe how to use our decision procedure to compute the minimum value ε * of ε such that the free space F ε contains a monotone path from (0, 0) to (m, n); this is the Fréchet distance F h (A, B).

We start by computing critical distances and the distance-query data structure in time O(mnk log k), as described in Lemma 5.6. We then sort the O(mn) critical distances. Using the decision procedure, we can compare the optimal distance ε * with any critical distance ε in O(mn log k) time. By binary search, we can, repeating this step O(log mn) times, compute an interval [ε -, ε + ] that contains ε * but no critical distances.

We then apply Megiddo's parametric search technique [START_REF] Megiddo | Applying parallel computation algorithms in the design of serial algorithms[END_REF]; see also [START_REF] Cole | Slowing down sorting networks to obtain faster sorting algorithms[END_REF][START_REF] Van Oostrum | Parametric search made practical[END_REF]. Parametric search combines our decision procedure with a 'generic' parallel algorithm whose combinatorial behavior changes at the optimal value ε * . Alt and Godau observe that one of two events occurs when ε = ε * :

• For some integers i, i , j, the bottom endpoint of v i,j and the top endpoint of v i ,j lie on the same horizontal line.

• For some integers i, j, j , the left endpoint of h i,j and the right endpoint of h i,j lie on the same vertical line.

Thus, it suffices to use a 'generic' algorithm that sorts the O(mn) endpoint values of all non-empty segments h i,j and v i,j , where the value of an endpoint (s, j) of h i,j is s, and the value of an endpoint (i, t) of v i,j is t.

We use Cole's parallel sorting algorithm [START_REF] Cole | Parallel merge sort[END_REF], which runs in O(log N ) parallel steps on O(N ) processors, as our generic algorithm. Each parallel step of Cole's sorting algorithm needs to compare O(mn) endpoints. The graph of an endpoint, considered as a function of ε, is monotone and made of O(k) hyperbolic arcs; see the proof of Lemma 5.5). It follows that the sign of a comparison between two endpoints may change at O(k) different values of ε, which can be computed in O(k) time. Applying the parametric search paradigm requires the following operations for each parallel step of the sorting algorithm:

• Compute the O(mnk) values of ε corresponding to the changes of sign of the O(mn) comparisons. This can be done in O(mnk) time and O(mnk) space.

• Apply binary search to these values by median finding, calling the decision procedure to discard half of them at each step of the search. This takes O(mnk + T d log(mnk)) time, where T d = O(mn log k) is the running time of our decision procedure. We obtain this way an interval for ε where each of the O(mn) comparisons has a determined sign.

• Deduce in O(mn log k) time the sign of each of the O(mn) comparisons within the previously computed interval.

Since the underlying sorting algorithm requires O(log mn) parallel steps, the resulting parametric search algorithm runs in time O(mn log(mn)(k + log k log(mnk)).

The distance query data structure requires O(mnk) space. We require O(mnk) additional space to simulate sequentially each parallel step of the sorting algorithm; we can re-use this space for subsequent parallel steps. Therefore, the total space complexity of our algorithm is O(mnk).

Lemma 5.7. Given a direct geodesic in a minimal relative homotopy class h, the Fréchet distance F h (A, B) can be computed in O(N 3 log N ) time and O(N 3 ) space, where N = m + n + k + 2 is the total input size.

Cook and Wenk propose a more efficient and practical randomized alternative to parametric search in their algorithm to compute geodesic Fréchet distance [START_REF] Atlas | Geodesic Fréchet and Hausdorff distance inside a simple polygon[END_REF][START_REF] Atlas | Geodesic Fréchet distance inside a simple polygon[END_REF]. A direct application of their technique reduces the time to search for the optimal distance ε * to O(N 2 log 2 N ) with high probability. Unfortunately, because our decision procedure relies on data structures built in the preprocessing stage, which requires O(N 3 log N ) time and O(N 3 ) space, Cook and Wenk's randomized technique does not improve the overall running time of our optimization algorithm.

Putting Everything Together

Finally, to compute the homotopic Fréchet distance F(A, B) in the plane minus a set of point obstacles, we enumerate the O(N 4 ) minimal homotopy classes and compute F h (A, B) for each minimal homotopy class h. Similarly, for polygonal obstacles, we construct a set of O(N 6 ) relative homotopy classes that includes at least one optimal class, and then for each class h in that set, we compute F h (A, B).

Theorem 5.8. The homotopic Fréchet distance between two polygonal curves in the plane minus a set of points can be computed in O(N 7 log N ) time and O(N 3 ) space, where N = n + m + k + 2 is the total input complexity.

Theorem 5.9. The homotopic Fréchet distance between two polygonal curves in the plane minus a set of polygonal obstacles can be computed in O(N 9 log N ) time and O(N 3 ) space, where N = n + m + k + 2 is the total input complexity.

Reduction to Geodesic Fréchet Distance?

Finally, we sketch a promising approach to a randomized algorithm which is faster with high probability than that of Section 5.5, using the recent algorithm of Cook and Wenk [START_REF] Atlas | Geodesic Fréchet and Hausdorff distance inside a simple polygon[END_REF][START_REF] Atlas | Geodesic Fréchet distance inside a simple polygon[END_REF] for computing the geodesic Fréchet distance between two polygonal curves inside a simple polygon. In this setting, the geodesic and homotopic Fréchet distances are identicalbecause the interior of a simple polygon is simply connected, there is only one relative homotopy class of leashes. Cook and Wenk's algorithm runs in O(p + M 2 log(pM ) log M ) time with high probability, and in O(p+M 3 log(pM )) time in the worst case, where M = m+n is the total complexity of the curves and p is the complexity of the enclosing polygon.

Recall that F h (A, B) is the equivalent to the classical Fréchet distance between A and B, except that distances are measured by shortest paths in relative homotopy class h. Let A and B be lifts of A and B to the universal cover E, chosen so that the shortest path σ h (s, t) between any two points A(s) and B(t) is a lift of a geodesic in homotopy class h. Then F h (A, B) is also equal to the geodesic Fréchet distance between A and B in E. However, the universal cover E is infinite, so we cannot pass it as input to Cook and Wenk's algorithm.

Recall that we are given a direct geodesic σ h (s, t) in class h. Let ∆ denote the infinite triangulation of E defined by Hershberger and Snoeyink [START_REF] Hershberger | Computing minimum length paths of a given homotopy class[END_REF]. Finally, let Π denote the union of triangles in ∆ that intersect either the lifted curves A or B or the shortest path σ h (s, t). Because the dual graph of ∆ is a tree, the shortest path from any point on A to any point on B must be contained in Π. It is not hard to prove that the region Π has complexity O(mk + nk), and it can be computed in O(mk + nk) time.

At this point, it is tempting to invoke Cook and Wenk's algorithm with the curves A and B and the region Π as input. However, we face a major technical hurdle: Π is not a simple polygon. Although Π is simply connected, has a locally Euclidean metric, and is bounded by straight line segments, it cannot be isometrically embedded in the plane without self-intersection.

We conjecture that Cook and Wenk's algorithm can be generalized to this setting with no loss of performance; this would require also generalizing the shortest-path query data structures of Guibas and Hershberger [START_REF] Guibas | Optimal shortest path queries in a simple polygon[END_REF][START_REF] Hershberger | A new data structure for shortest path queries in a simple polygon[END_REF] on which Cook and Wenk's algorithm relies. Specifically, we believe these algorithms can be modified to accept arbitrary simply connected boundary-triangulated 2-manifolds, as defined by Hershberger and Snoeyink [START_REF] Hershberger | Computing minimum length paths of a given homotopy class[END_REF], instead of simple polygons. If our conjecture is correct, this approach would imply an algorithm to compute F h (A, B) in O(N 2 log 2 N ) time with high probability, and still in O(N 3 log N ) in the worst case, using only O(N 2 ) space.

Extensions

Closed Curves

Formally, a closed curve in a topological space S is a continuous function from the circle S 1 = R/Z to S. Let A and B be two closed curves in S. A free homotopy between A and B is a continuous function h : S 1 × [0, 1] → S, such that h(•, 0) = A and h(•, 1) = B; unlike a homotopy between paths, a free homotopy between closed curves does not keep any point fixed. If such a function exists, the two closed curves are freely homotopic. A closed curve A is contractible if it is freely homotopic to a constant function (that is, a single point).

Equivalently, A is contractible if there is a continuous function from the unit disk into S whose restriction to the disk boundary is A. A reparameterization of S 1 is a continuous monotone surjection α : S 1 → S 1 of index 1. A reparameterization of a closed curve is the composition of this closed curve with a reparameterization of S 1 . Now homotopic Fréchet distance can be defined almost exactly as it is for paths. Specifically, a leash map between A and B is a free homotopy between some reparameterization of A and some reparameterization of B; the length of a leash map is the maximum length of any leash (t, •); and the homotopic Fréchet distance is the infimal length of any leash map:

F(A, B) := inf leash map : S 1 ×[0,1]→S max t∈S 1 len( (t, •)) .
If there is no leash map between A and B-that is, if A and B are not homotopic-then we define F(A, B) = ∞. Our algorithm will automatically detect this situation.

In this section, we show how to adapt our algorithm to compute homotopic Fréchet distance between paths to compute the homotopic Fréchet distance between closed curves. Our derivation is complicated by the fact that different geodesics with the same endpoints can be homotopic relative to the closed curves A and B.

For notational convenience, we define S 1 as R/Z, or equivalently, as the unit interval [0, 1] with its endpoints identified. A closed curve A shifted by s is the closed curve denoted by A + s where (A + s)(t) = A((s + t) mod 1). For any closed curve A, let A / : [0, 1] → E be the corresponding ordinary curve, defined by setting A / (t) = A(t) for all t ∈ [0, 1].

Consider a leash map between A and B, that is a free homotopy between some reparameterizations A • α and B • β. For every parameter u in S 1 , can be 'cut' along the leash (u, •) to form a leash map u/ between the paths (A + α(u)) / and (B + β(u)) / . Formally, ∀(s, t) ∈ [0, 1] × [0, 1], u/ (s, t) := ((u + s) mod 1, t). The cut leash map u/ and the original leash map have equal length. Also note that the initial and final leashes in the cut leash map coincide: u/ (0, .) = u/ (1, .). Conversely, for any parameters s and t in S 1 , a leash map / between the paths (A + s) / and (B + t) / can be 'glued' to form a leash map between the closed curves A and B, if and only if / (0, .) = / (1, .).

The following lemma further characterizes leashes that can appear in leash maps between closed curves. Suppose λ is a leash with endpoints A(s) and B(t). Let C(λ) denote the closed curve obtained by concatenating (A + s) / , followed by λ, followed by the reversal of (B + t) / , followed by the reversal of λ. We call a leash λ valid if C(λ) is contractible; see Figure 6. Lemma 6.1. Let A and B be closed curves in some topological space S. In every leash map between A and B, every leash is valid. Conversely, any valid leash belongs to at least one leash map between A and B.

Proof. Let be a leash map between A and B, and let (u, .) be a leash in . The cut leash map u/ is a continuous function defined on the topological disk [0, 1] × [0, 1]. By definition, the restriction C of u/ to the boundary of this disk (which we can identify to S 1 ) is contractible. Thus, since C is just a reparameterization of C( (u, .)), C( (u, .)) is freely homotopic to C and contractible. Conversely, consider any valid leash λ, with endpoints A(s) and B(t). Because λ is valid, the cycle C(λ) is contractible. Thus, there is a continuous map h from the unit disk into S whose restriction to the boundary of the disk is C(λ). By identifying the unit disk to [0, 1] × [0, 1], it appears that h is actually a leash map between (A + s) / and (B + t) / with the property that h(0, .) = h(1, .) = λ. It follows that h can be glued to form a leash map between A and B.

If two leashes are homotopic relative to (A + s) / and (B + t) / , then either both leashes are valid or neither leash is valid; we say that a leash map or a relative homotopy class is valid if its elements are valid leashes. Let F v ((A + s) / , (B + t) / ) denote the infimum length of a valid leash map between the ordinary curves (A + s) / and (B + t) / . Cutting every leash map along its leash between some A(s) and B(0) and using the previous lemma, we can express the homotopic Fréchet distance between A and B as:

F(A, B) = inf s∈S 1 F v ((A + s) / , B / ) Lemma 3.1 implies that F v ((A + s) / , B / )
is the length of a geodesic leash map in some valid homotopy class (relative to (A + s) / and B / ). The proofs of Lemmas 4.1 and 4.2 extend to valid relative homotopy classes without modification; thus,

F v ((A + s) / , B / ) = inf valid minimal homotopy class h F h ((A + s) / , B / ).
Finally, Propositions 4.8 and 4.3 also extend to valid relative homotopy classes. The extension uses the fact that replacing a common non-empty subpath of all the leashes of a valid leash map yields another valid leash map. The two propositions imply that any valid minimal homotopy class of leashes between (A + s) / and B / contains a direct geodesic.

For any valid direct geodesic γ and any shift value s, let [γ] s denote the homotopy class of γ relative to (A + s) / and B / . The preceding discussion implies that

F(A, B) = min valid direct geodesic γ F γ (A, B),
where F γ (A, B) := min s∈S 1 F [γ]s ((A + s) / , B / ). Alt and Godau extend their algorithm to compute the classical Fréchet distance to closed curves [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF]. Their algorithm solves the decision problem, whether the Fréchet distance is at most a given ε, by concatenating two copies of the free space diagram from Section 5 and therefore has 2m × n cells. The decision problem can be rephrased as follows: Is there a shift s, such that the free space diagram contains a monotone path from (s, 0) to (m + s, n)? Alt and Godau augment their representation of the free space diagram with additional pointers to answer this question efficiently.

We build an analogous free space diagram to determine, for any valid direct geodesic γ and threshold ε > 0, whether F γ (A, B) ε. Aside from the preprocessing described in Section 5.2, our approach is identical to that of Alt and Godau.

First, we modify our decision algorithm from Section 5.3 to determine whether there is a valid leash map in homotopy class h whose length is at most ε. The modified decision algorithm is more expensive by a factor of O(log mn) because it constructs a data structure analogous to that of Alt and Godau which helps to determine whether there exists a shift s, such that the 2m × n free space diagram contains a monotone path from (s, 0) to (m + s, n). The running time of the new decision procedure is O(mn log k log(mn)) = O(N 2 log 2 N ). The space complexity remains O(mnk) = O(N 3 ) because the space required for storing the pointers in the data structures is O(mnk).

Then, we modify our parametric search which determines the smallest ε for which there exists a valid leash map of length at most ε. As observed by Alt and Godau, we need to consider additional critical values of ε, but those also are the critical distances between a vertex of one curve and an edge of the other curve. Therefore, there are O(mnk) such critical values of ε, as argued in Section 5.2. Hence, the running time of the optimization procedure is now O(mn log(mn)(k + log k log(mn) log(mnk))) = O(N 3 log N ). We emphasize that the increased cost of the decision procedure does not lead to a similar increase in the overall running time of our optimization procedure.

Because each direct geodesic has complexity O(k) = O(N ), we can test whether a direct geodesic λ is valid in O(N 3/2 log N ) time using an algorithm of Cabello et al. [START_REF] Cabello | Testing homotopy for paths in the plane[END_REF] to decide whether C(λ) is contractible. Thus, we can compute a set of valid minimal homotopy classes that contains an optimal one in O(N 11/2 log N ) time for point obstacles, or in O(N 15/2 log N ) for polygonal obstacles. Finally, to compute the homotopic Fréchet distance, we compute the optimum leash map in the relative homotopy class of each valid direct geodesic. Theorem 6.2. The homotopic Fréchet distance between two closed polygonal curves in the plane minus a set of obstacles can be computed in O(N 7 log N ) time and O(N 3 ) space if the obstacles are points, or in O(N 9 log N ) time and O(N 3 ) space if the obstacles are polygons.

Variants of Fréchet Distance

Finally, we briefly consider two natural variants of homotopic Fréchet distance that can also be computed using our techniques.

The weak Fréchet distance is a variant of the classical Fréchet distance without the requirement that the endpoints move monotonically along their respective curves-the dog and its owner are allowed to backtrack to keep the leash between them short. Alt and Godau [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF] gave a simpler algorithm for computing the weak Fréchet distance, using a graph shortest-path algorithm instead of parametric search. A similar simplification of our algorithm computes the weak homotopic Fréchet distance between curves in the plane minus polygonal obstacles in polynomial time. The proofs of Lemmas 4.1 and 4.2, as well as Propositions 4.8 and 4.3, extend to the weak variant of homotopic Fréchet distance without modification. We reiterate that Propositions 4.8 and 4.2 do not imply that the optimal (weak) leash map contains a direct geodesic, only the optimal relative homotopy class. Thus, we can define and compute the weak homotopic Fréchet distance as the minimum, over all relative homotopy classes h that contain a direct geodesic, of the weak Fréchet distance with respect to shortest path lengths in h.

The discrete homotopic Fréchet distance, also called the coupling distance, is an approximation of the Fréchet metric for polygonal curves defined by Eiter and Mannila [START_REF] Eiter | Computing discrete Fréchet distance[END_REF]. The discrete Fréchet distance considers only positions of the leash where its endpoints are located at vertices of the two polygonal curves and never in the interior of an edge. This special structure allows the discrete Fréchet distance to be computed by an easy dynamic programming algorithm. As usual, Propositions 4.8 and 4.3 imply that we can define and compute the discrete homotopic Fréchet distance as the minimum, over all relative homotopy classes h that contain a direct geodesic, of the discrete Fréchet distance with respect to shortest path lengths in h.

Conclusion

In this paper, we introduced a natural generalization of the Fréchet distance between curves to more general metric spaces, called the homotopic Fréchet distance. We described a polynomial-time algorithm to compute the homotopic Fréchet distance between polygonal curves in the plane with point or polygon obstacles.

Improving the running time of our algorithms is the most immediate outstanding open problem. We described one promising approach in Section 5.6, which would require generalizing Cook and Wenk's algorithm for geodesic Fréchet distance to more general simply-connected spaces. We also conjecture that the running time can be improved by optimizing leash maps in every minimal homotopy class simultaneously. Since shortest paths between the same endpoints but belonging to different homotopy classes are related, we expect to (partially) reuse the results of shortest path computations in one homotopy class when we consider other homotopy classes. Finally, for polygonal obstacles, an exact characterization of minimal homotopy classes would almost certainly lead to a significantly faster algorithm.

It would be interesting to compute homotopic Fréchet distance in spaces more general than those considered in the current paper. In particular, we are interested in computing the homotopic Fréchet distance between two curves on a convex polyhedron, generalizing the algorithm of Maheshwari and Yi for classical Fréchet distance [START_REF] Maheshwari | On computing Fréchet distance of two paths on a convex polyhedron[END_REF]. The vertices of the polyhedron are 'mountains' over which the leash can pass only if it is long enough. Shortest paths on the surface of a convex polyhedron do not vary continuously as the endpoints move, because of the positive curvature at the vertices, so we cannot consider only geodesic leash maps.

Finally, it would also be interesting to consider the homotopic Fréchet distance between higher-dimensional manifolds; such problems arise with respect to surfaces in configuration spaces of robot systems. Ordinary Fréchet distance is difficult to compute in higher dimensions, although the weak Fréchet distance between two triangulated surfaces can be computed in polynomial time [START_REF] Alt | Semi-computability of the Fréchet distance between surfaces[END_REF].

Figure 1 :

 1 Figure 1: Leash maps for three inputs. Dashed curves between matching numbers represent intermediate leashes.

Figure 2 :

 2 Figure 2: (a) Triangulation ∆ of a region S, which is a rectangle minus a triangle. (b) Triangulation ∆ of the universal cover S of S.

Figure 3 :

 3 Figure3: Curves and obstacles that admit Ω(mnk 2 ) minimal relative homotopy classes and Ω(mnk 4 ) relative homotopy classes of direct geodesics.

Figure 5 :

 5 Figure 5: (a) An instance where no optimal leash map contains a line segment. (b) Up to symmetry, the only four relative homotopy classes of line segments. (c) Half of a symmetric leash map in homotopy class I.

Lemma 5. 2 .

 2 Let o be a fixed point in the plane, and let ϕ o : R 2 → R be the function p → -→ op . The gradient of ϕ o at any point p = o is -→ op/ -→ op . The function ϕ o is convex everywhere, and of class C 1 everywhere except at o.

Figure 6 :

 6 Figure 6: (a) Two invalid leashes and a valid leash between two homotopic closed curves. (b) The cycles C(λ) corresponding to each leash λ; only the last cycle is contractible.

inria-00438463, version 1 -3 2009

Acknowledgments

This research was initiated during a visit to INRIA Lorraine in Nancy, made possible by a UIUC-CNRS-INRIA travel grant. Research by Erin Chambers and Jeff Erickson was also partially supported by NSF grant DMS-0528086; Erin Chambers was additionally supported by an NSF graduate research fellowship. Research by Shripad Thite was partially supported by the Netherlands Organisation for Scientific Research (NWO) under project number 639.023.301 and travel by INRIA Lorraine.

We would like to thank the anonymous referees for their careful reading of the paper and their numerous suggestions for improvement. Finally, we thank Hazel Everett and Sylvain Petitjean for useful discussions, and Kira and Nori for great company and several walks in the woods.

Proposition 4.8. In the case of point obstacles, a relative homotopy class is minimal if and only if it contains a line segment.

Proof. One direction of the proof is straightforward: Let h be the relative homotopy class of the line segment σ from A(s) to B(t). By our non-degeneracy assumption, up to slightly moving σ, we may assume that it touches no obstacle point and apply Lemma 4.7.

To prove the opposite implication, we consider a minimal relative homotopy class h. As in the proof of Proposition 4.3, we define π h to be the intersection of the shortest paths between A(s) and B(t), for all s, t ∈ [0, 1]. If π h = ∅, then the proof of Lemma 4.4 already implies that h contains a line segment, so we assume that π h = ∅. In particular, there is a lifted obstacle point p such that for any s, t ∈ [0, 1], the shortest path σ h (s, t) passes through p. Let θ(s, t) denote the turning angle of σ h (s, t) at p.

For all s and t, the path σ h (s, t) is a shortest path, so θ(s, t) must lie outside the open interval (-π, π). This turning angle is a continuous function of s and t, so we can assume without loss of generality that it is always at least π. In other words, we assume that every path σ h (s, t) winds counterclockwise around p. Recall that no three vertices of the input are collinear by our non-degeneracy assumption, so the minimum of θ(s, t) is not a multiple of π, and can therefore be written as 2πx + y for some integer x and some angle y ∈ (-π, π). Now p is a lift of some obstacle point p, and σ h (s, t) similarly projects to a geodesic σ h (s, t). For each s and t, let τ (s, t) denote the path meeting the same obstacles in E in the same order and with the same turning angles as σ h (s, t), except that the turning angle at p is reduced by 2πx. All paths τ (s, t) belong to a single relative homotopy class, which we denote by h .

For every s and t, the paths τ (s, t) and σ h (s, t) have precisely the same length. This implies that σ h (s, t) is never longer than σ h (s, t); thus, h h. Now let s and t be parameters such that θ(s, t) is minimized, and write θ(s, t) = 2πx + y for some integer x and some y ∈ (-π, π). By construction, the turning angle of τ (s, t) at p equals y. Thus τ (s, t) is not a geodesic, so σ h (s, t) is strictly shorter than τ (s, t), which has the same length as σ h (s, t). With the previous paragraph, this proves h ≺ h, contradicting our initial assumption that h is minimal.

Corollary 4.9. In the case of point obstacles, we can enumerate a superset of the minimal relative homotopy classes of size O(mnk 2 ) in O(mnk 2 ) time.

Proof. In the proof of Corollary 4.6, we saw how to enumerate the O(mnk 2 ) relative homotopy classes of line segments in O(mnk 2 ) time; every minimal homotopy class contains a line segment by Proposition 4.8.

We now conclude:

Proof of Proposition 5.1. Let i,j be an arbitrary grid cell. Choose parameters s, s ∈ [i-1, i] and t, t ∈ [j -1, j]. We claim that the function ψ : [0, 1] → E defined by setting

respectively. Then ψ is exactly the function d that is proved convex in Lemma 5.4.

We conclude that the restriction of dist h to any line segment in i,j , which is a map of the form ψ above, is convex. This completes the proof that the restriction of dist h to any grid cell i,j is convex.

Preprocessing for Distance Queries

The only significant difference between our algorithm and Alt and Godau's is that we require additional preprocessing to compute several critical distances and an auxiliary data structure to answer certain distance queries. (If there are no obstacles, each critical distance can be computed, and each distance query can be answered, in constant time.)

There are three types of critical distances:

• endpoint distances dist h (0, 0) and dist h (m, n),

and

Given integers i and j and any real value ε, a horizontal distance query asks for all values of t ∈ [j -1, j] such that dist h (i, t) = ε, and a vertical distance query asks for all values of s ∈ [i -1, i] such that dist h (s, j) = ε. The convexity of dist h within any grid cell implies that any distance query returns at most two values.

We first describe how to preprocess a single vertical edge in the parameter grid to answer distance queries; critical values are automatically computed during the preprocessing. Obviously a similar result applies to horizontal grid edges.

Lemma 5.5. Suppose we are given a point p and a line segment = xy, parameterized over [0, 1], as well as the geodesic σ h (p, x) and its length dist h (p, x). In O(k log k) time, we can build a data structure of size O(k) such that for any ε, all values t ∈ [0, 1] such that dist h (p, (t)) = ε can be computed in O(log k) time. We also report the critical vertex-edge distance dist h (p, ), the path σ h (p, y), and its length dist h (p, y).

Proof. We first compute a constrained Delaunay triangulation [START_REF] Chew | Constrained Delaunay triangulations[END_REF] of the obstacles P , the segment , and point p in time O(k log k). This triangulation includes and the edges of polygons in P as edges.

We apply the following observations used in the funnel algorithm for computing shortest homotopic paths [START_REF] Chazelle | A theorem on polygon cutting with applications[END_REF][START_REF] Lee | Euclidean shortest paths in the presence of rectilinear barriers[END_REF][START_REF] Hershberger | Computing minimum length paths of a given homotopy class[END_REF]. The shortest homotopic paths σ h (p, x) and σ h (p, y) may share a common subpath and then split at some vertex v; this vertex is then the apex of two concave chains that form a funnel with base xy. Each concave chain has complexity at most k and intersects a given edge of the triangulation at most twice.

The geodesic from p to x may have complexity greater than O(k), but (as observed above) the concave chain from v to x has at most O(k) segments. Our goal is to find a vertex w on σ h (p, x) such that the path from w to x contains v. In other words, the chain from w to x along σ h (p, x) has complexity O(k) and contains the concave funnel path.

To find w, walk along the geodesic from x to p. If we find a vertex where the chain is not concave, we must have passed v, so we mark the non-concave vertex as w. If we ever re-cross a segment of the triangulation a second time, we again must have passed the funnel apex v so we can mark the second crossing as w. (We walked along O(k) edges of the chain to find w.) Let π h be the portion of σ h (p, x) between p and w, and τ 1 be the portion of σ h (p, x) between w and x.

We know that π h is contained in σ h (p, y), since w is before the apex of the funnel v. Let τ 2 be the portion of σ h (p, y) between w and y; this can be computed in O(k) time using the funnel algorithm. Given τ 2 , we can then find the apex of the funnel v in O(k) time.

Imagine extending each line segment on the concave chains until it intersects , the line connecting x and y. Between the two concave chains, the combinatorial description of the distance function changes only at points where the extended lines meet . To answer distance queries, we record the O(k) intersections of the extended lines with . For each of the resulting intervals, record the (fixed) length of the geodesic up to the first vertex in the extended line, as well as the equations of the two lines that bracket the interval. In constant time per interval, we can also compute and store the value t * ∈ [0, 1] such that dist h (p, (t * )) is minimized; this gives the desired value dist h (p, ).

The funnel data structure requires O(k) space to store the O(k) combinatorial changes to the leash as its endpoint sweeps = xy. Now given this data structure, we answer distance queries as follows. If the distance queried is smaller than dist h (p, ), we return the empty set. If it is equal to dist h (p, ), we return (t * ). If it is larger than dist h (p, ), we do two binary searches, one on the intervals between x and (t * ) and the other on the intervals between (t * ) and y.

Lemma 5.6. Given any shortest path in the minimal relative homotopy class h, we can compute all critical distances, in O(mnk log k) time and using O(mnk) space, as well as build a data structure of size O(mnk) that can answer any horizontal or vertical distance query in O(log k) time.

Proof. We preprocess each edge of the parameter grid as described in Lemma 5.5. We start from the vertex (i, j) that is our given input, either a straight line segment or a direct geodesic. We then walk on the edges of the grid, visiting each edge at least once and at most twice. During this walk, at each current vertex (i, j), we maintain the shortest homotopic path σ h (i, j) and its length dist h (i, j). Each time we walk along an edge, we apply Lemma 5.5 to preprocess it and to compute the shortest homotopic path corresponding to the target vertex of that edge. Each step takes O(k log k) time, and there are O(mn) edges, whence the running time. As we walk along an edge of the parameter grid, we use a deque to push and pop the obstacle vertices along the leash in constant time per operation. Since at most k vertices are pushed onto the deque for each grid edge, the total size of the deque is O(mnk).


Alexandre Reiffers-Masson

Yezekael Hayel

Eitan Altman

Guillaume Marrel

A Generalized Fractional Program for Maximizing Content Popularity in Online Social Networks

In this paper, we consider a "generalized" fractional program in order to solve a popularity optimization problem in which a source of contents controls the topics of her contents and the rate with which posts are sent to a time line. The objective of the source is to maximize its overall popularity in an Online Social Network (OSN). We propose an efficient algorithm that converges to the optimal solution of the Popularity maximization problem.

I. INTRODUCTION

As noticed in [START_REF] Heidemann | Online social networks: A survey of a global phenomenon[END_REF], one key challenge on Online Social Networks (OSNs) is to find an answer to the following question: how to become popular on an OSN? In the context of marketing and sales, the importance for a company to be popular on OSNs is of main interest. Indeed, in a recent report (http ://www.socialmediaexaminer.com), the numberone challenge marketers want to see resolved (92%) is: "What social tactics are most effective?". Many empirical studies have been developed to answer this issue [START_REF] Bandari | The pulse of news in social media: Forecasting popularity[END_REF], [START_REF] Lakkaraju | What's in a name? understanding the interplay between titles, content, and communities in social media[END_REF], [START_REF] Susarla | Social networks and the diffusion of user-generated content: Evidence from youtube[END_REF], [START_REF] Szabo | Predicting the popularity of online content[END_REF]. An optimization approach has been proposed in the case where a marketer, also called a source in our paper, will pay the OSNs to get his content promoted [START_REF] Huang | Recent developments in dynamic advertising research[END_REF], [START_REF] Pellegrini | Differential games of competition in online content diffusion[END_REF], [START_REF] Helmes | Optimal advertising and pricing in a class of general new-product adoption models[END_REF], [START_REF] Mesak | On the superiority of pulsing under a concave advertising market potential function[END_REF]. In this paper, we propose a novel approach, where a source controls the flow of posted contents and the topics associated to each content.

In order to provide an answer to this challenge of popularity in OSNs, our methodology pursues the following steps. We start by understanding how content popularity behaves on OSNs. A vast amount of literature about this question exists (for instance [START_REF] Szabo | Predicting the popularity of online content[END_REF], [START_REF] Lakkaraju | What's in a name? understanding the interplay between titles, content, and communities in social media[END_REF] and references therein). In this paper, we assume that the popularity of a content is highly related to its position in a News Feed. This assumption is satisfied in Facebook/Twitter-style content posting system. In Facebook for example, a News Feed is a feed where the user's friends and subscribed pages' posts are published. These posts can be displayed in a top-down anti-chronological order from the newest to the oldest if you use the most recent option. More explanations about the News Feed dynamics and models are provided in [START_REF] Altman | Competition over timeline in social networks[END_REF]. Then, we assume that the maximization of content popularity is similar to the oil producer problem [START_REF] Sethi | Optimal Control Theory-Applications to Management Science and Economics[END_REF]. This problem assumes that an oil producer would pump oil from a single well and after some time, when the return on the well decreases, he would abandon that source and start exploiting another one. In the case of OSNs, we assume that a well pumping is equivalent to a content staying in the first position of the News Feed. At some point, the popularity return will decrease and it is recommended to send a new content, 1 Indian Institute of Science, India 2 University of Avignon, France 3 INRIA Sophia Antipolis, France meaning to exploit a new source. Our model of popularity is expressed as a system of discrete differential equations. Finally, we study the maximization of content popularity. We assume that a company (or source) that tries to be popular will focus on two types of control: the number of posted contents and the topics of its contents. A continuous formulation of the popularity maximization problem is proposed and we prove its equivalence with a "generalized" fractional program.

The main contributions of this paper are as follows: First we formally describe the popularity maximization problem, an advertising optimization problem on online social networks. We also define pertinent popularity functions. We use the data to estimate the parameters of these popularity functions and assesses the quality of their correlation level with real data. Then we prove that for standard popularity functions, we can use convex optimization to design efficient algorithms. Finally, we list some of the constraints that generalized the popularity maximization problem and that can be expressed as convex inequalities.

Related works: Several publications have appeared in recent years documenting the maximization of popularity in online social networks. Usually there are two approaches. The first one is advertising allocation where a source decides the amount of money he is willing to pay to the OSN in order to get its content promoted. Papers about this subject are often linked to optimal control [START_REF] Helmes | Optimal advertising and pricing in a class of general new-product adoption models[END_REF], [START_REF] Huang | Recent developments in dynamic advertising research[END_REF]. The second approach adopted by researchers is to model the propagation of information through a cascade model in a social network. Then the question would be to find which node a source would need to target in order to maximize the spread of its contents [START_REF] Kempe | Maximizing the spread of influence through a social network[END_REF], [START_REF] Chen | Scalable influence maximization for prevalent viral marketing in large-scale social networks[END_REF].

We propose in this paper to define a new approach based on flow control. In this framework the amount of contents is sent by the source affects their popularity. To the authors' best knowledge, very few publications can be available in the literature that address this issue. In [START_REF] Altman | Competition over timeline in social networks[END_REF], [START_REF] Reiffers-Masson | A time and space routing game model applied to visibility competition on online social networks[END_REF] the authors have used a Poisson process to model the flow of contents on a unique news feed and define a visibility measure of a content. The major drawback of this approach is that the authors are concerned with visibility and not popularity, so the measure is not realistic. Our visibility measure is more general. Moreover, they do not list all the possible constraints that can occur in a popularity maximization problem and only deals with a simple optimization problem.

Organization of the paper. In section III, we propose a mathematical model of the popularity. The section IV is about the Popularity maximization problem. First, we provide its definition, and then we prove, for standard popularity Fig. 1: Histogram of first position fraction of popularity functions, its equivalence with a pseudoconcave optimization problem. Finally, an algorithm that converges to the optimal solution is provided.

II. DATA EXPERIMENT: IMPORTANCE OF THE FIRST

POSITION

In this section, we develop an empirical experiment to illustrate how the popularity of a content is related to the position of a content in the News Feed. We inspect specifically the relationship between the popularity (in terms of number of comments) and the duration a content spends in the first position in the News Feed. Our dataset is composed of contents posted on a specific Facebook page associated to a french political personality (Anne Hidalgo https://www.facebook.com/HidalgoAnne Paris mayor). We have extracted 892 contents and 39083 associated comments with timestamps using the Netvizz app [START_REF] Rieder | Studying facebook via data extraction: the netvizz application[END_REF]. The dataset was collected from June 15 2013 to December 01 2014. The histogram depicted on fig. 1 represents the proportion (or density) of contents and the relative proportion of popularity obtained in the first position in the News Feed. For example, this histogram shows that close to 50% of the contents get 90% of their popularity while being in the first position in the News Feed. Also, it shows that more than 70% of the posted contents reached 80% of their total popularity while the content is in the first position in the News Feed. Therefore, it makes sense to consider a popularity model, described in detail in the next section, by assuming that almost all the popularity of a content is obtained while the content is in the first position in the News Feed.

III. POPULARITY MODEL

The level of popularity associated with a topic can be defined as the total number of comments, given by OSN subscribers, of all contents related to this topic. We model the dynamic of the popularity of each topic through a discrete time dynamical system.

Let c ∈ C := {1, . . .C} be a topic index. Topics can be "Sport", "Arts", "Entertainment", etc. Note that our model is scalable and can be applied to any number of topics studied. The random state variable X c (t) ∈ N + is the total number of comments that gets all content related to topic c in a News Feed, at a given time t ∈ R + . For each topic c ∈ C , the arrival process of a type c content in the News Feed is a Poisson Point Process with intensity 0 < λ c < ∞. All arrival processes are independent. We also consider that there exists a flow, with intensity σ of posts which are not related to any topics. Let t m ∈ R + be the random variable denoting the arrival instant of the m th content, with m ∈ N * + . As the arrival processes are independent Poisson Point Processes, the inter-arrival time denoted by ∆ m between consecutive contents arrived at t m and t m-1 , follows an exponential distribution with rate Λ := ∑ C d=1 λ d . Therefore, the probability associated with the inter arrival time ∆ m of two consecutive contents is given by:

∀m, δ ∈ R + , P(∆ m < δ ) = δ 0 Λe -Λt dt. (1)
A content can also be called a post on an OSNs. We define by η(m) ∈ C the topic of the m th content. According to [START_REF] Durrett | Essentials of stochastic processes[END_REF], for each content m, η(m) is a random variable with the following distribution:

∀m, P(η(m) = c) = λ c Λ . ( 2
)
The evolution of the cumulative number of comments associated with topic c in a News Feed, is modeled by the following dynamical system:

X c (t m ) = X c (t m-1 ) + 1 {η(m-1)=c} Y c (∆ m ), (3)
X c (t 1 ) = 0, (4)
where the random variable Y c (δ ) ∈ N + , with δ ∈ R + , describes the cumulative popularity obtained by a content of type c that stays δ units of time on the first position of the News Feed.

The indicator 1 {η(m-1)=c} = 1 if the (m -1) th content is of type c, it is equal to 0 otherwise. We assume that for each c, the random variable Y c (δ ) is described by a non-homogeneous Poisson process with intensity E[Y c (δ )] := δ 0 f c (τ)dτ ≥ 0. The density function f c (.) are called popularity functions. In point of fact, in recent studies on analyzing social networks, nonhomogeneous Poisson processes have been used to model and predict the popularity of contents. For instance, in [START_REF] Ding | Video popularity prediction by sentiment propagation via implicit network[END_REF], [START_REF] Rizoiu | Can this video be promoted?-endogenous and exogenous popularity processes in social media[END_REF], the authors suggest to use Hawks processes in order to capture the evolution of the number of likes, comments of contents posted in online social networks. Hawks processes are composed of an exogenous and an endogenous part. The exogenous part corresponds to E[Y c (δ )] and the usual assumption about the shape of popularity functions is a uniform, an exponential or a power law distribution. Therefore, our model can be seen as the worst-case of the hawks model since the endogenous part is equal to zero.

From data experiments, we have been able to show that three types of popularity functions f c (.) have significant correlation levels with data. The first function is the constant function A closed form expression of the expected cumulative popularity for each topic c at any given time, is given in the following proposition.

f c (τ) = M c ,
Proposition 1: The expected cumulative popularity for each topic c at the time t m when the m th content is posted is given by:

∀c ∈ C ,t m ∈ R * + , E[X c (t m )] = m ∞ 0 λ c e -Λt t 0 f c (τ)dτdt. (5)
In table I we present mathematical closed form expression of the expected cumulative popularity E[X c (t m )] for different shapes of popularity functions of the non-homogeneous Poisson process.

f c (τ) M c a c τ + b c acbc (bc+τ) 2 E[X c (t m )] Mcλc Λ λ c -2ac Λ 3 -bc Λ a c λ c 1 Λ -b c e bcΛ Γ(0, b c Λ)
TABLE I: Examples of expected cumulative popularity and their closed form for the three different popularity functions.

IV. POPULARITY OPTIMIZATION

We now study the case where a single source maximizes the popularity of his content. The source determines the topics of his content and also the rate of sending contents per topic. Therefore, the source parameters are described by a vector (λ 1 , . . . , λ C ) where λ c ∈ [0, λ ] is the rate of contents related with topic c and λ is the maximum rate for each topic. Based on proposition 1, the average expected cumulative popularity X ∞ c associated with the content related to topic c posted by the source is equal to:

∀c ∈ C , X ∞ c := lim m→∞ 1 m E[X c (t m )] = ∞ 0 λ c e -Λt t 0 f c (τ)dτdt. (6)
Let λ cost := ∑ C c=1 γ c λ c be the cost perceived by the source for sending all his contents. This cost is considered to be in monetary units and we assume that the source wants to maximize the sum of its total average expected cumulative popularity minus the cost. Then, the net utility of the source is defined as follows:

U(λ 1 , . . . , λ C ) := C ∑ c=1 ∞ 0 λ c e -Λt t 0 f c (τ)dτdt - C ∑ c=1 γ c λ c . (7)
The Popularity Optimization Problem (PoP) is the following optimization problem defined below:

max λ c ∈[0,λ ], ∀c U(λ 1 , . . . , λ C ) (8)
A solution of this problem is denoted by a vector λ OA := (λ OA 1 , . . . , λ OA C ) ∈ argmax λ c ∈[0,λ ], ∀c U(λ 1 , . . . , λ C ). This optimization problem is note straightforward, and the purpose of the rest of this section is to find an equivalent form of the Popularity optimization problem that can be solved efficiently.

The next theorem proves that the optimal vector rate λ OA has a particular form. In fact, optimally, the source has to send content only related to a single and unique topic. Therefore, the Popularity optimization problem can be reduced to a simpler optimization problem.

Proposition 2: If the f c (•) is one of the standard functions used to model the popularity of a content (see table I), the Popularity optimization problem is equivalent to the following generalized fractional problem:

max c max λ c ∈[0,λ ] U c (λ c ) . (9)
According to [START_REF] Avriel | Generalized concavity[END_REF], optimization problems like the one defined in equation ( 9) is a pseudo-concave optimization problem. Then, for each c ∈ C , a projected gradient algorithm, converges to the optimal solution of (9). We have observed that under verified conditions, the Popularity optimization problem is a generalized fractional problem. Considering a gradient scheme, we are able to propose a global algorithm that first looks for the optimal rate for each topic c and then determines, by comparing all the outputs of the C gradient algorithms, on which topic the source has to send his contents.

We illustrate the convergence of the gradient scheme for the Popularity optimization problem on a numerical example, with f c (τ) = e -α c τ , where the parameters are C = 3 and (α 1 , α 2 , α 3 ) = (10, 5, 50). On Figure 4, each curve draws the iteration of the gradient algorithm for each type c. We observe the convergence of the gradient in each case as stated previously. Each dashed line shows the optimal flow λ OA c for each type c. The last iteration of the algorithm implies that the source has to create contents about the second topic, in order to become popular. Convex constraints: Some constraints can be introduced into the generalized the popularity optimization problem and that can be expressed as convex inequalities. For example, a minimum time T can be imposed on the interval between two contents in order to avoid spamming. This type of constraint is captured by 1 λ ≥ T . To go further in this direction, we can also assume that we are interested to bound the quantity of contents related to a specific topic sending per unit of time. In this case, the following constraint is modeled by

( λ c Λ ) k c ≤ d c
where k c is this bound for topic c. Both of these constraints are convex and can be easily integrated into our algorithm. Finally We can impose upper and lower bounds on the intensity ratio between two topics l c ≤ λ c λ d ≤ u c , as a measure of fairness between topics.

V. DISCUSSION AND CONCLUSION

In this paper, we define a popularity optimization problem in online social networks. The model is based on the observation that most of the popularity of a content is obtained during the time a content spent in the first position in a News Feed. We propose an optimization problem where a source determines optimally the rate of his content and also the topics. An equivalent fractional problem is obtained and an efficient gradient-based algorithm is proposed to determine the optimal solution. Finally, we provide to generalize our problem by considering constraints that do not alter the algorithm proposed.

  the second function is the linear function f c (τ) = a c τ + b c and the third function is the Michaelis-Menten model which is given by f c (τ) = a c b c (b c +τ) 2 . M c is the parameter of the first function and a c and b c are the parameters of the second and third functions. Examples of estimations are depicted in

Fig. 2 :

 2 Fig. 2: Estimation of popularity functions for three posts.

Fig. 2 .

 2 Fig. 2. The levels of correlation are described in Fig. 3 for three different posts. We observe that for the type 3 function the average level of correlation is 0.95 which is very significant. And even for the simplest constant function 1 f c (τ) = M c , the level of correlation has an average of 0.7. Moreover the median of the p-value associated to each curve is below 1e-05 (Function 1: 4.25e-06, Function 2: 1.09e-10, Function 3: 1.07e-17).

Fig. 3 :

 3 Fig. 3: Level of correlation for the three functions: function 1 f c (τ) = M c , function 2 f c (τ) = a c τ +b c and function 3 f c (τ) = a c b c (b c +τ) 2 .

Fig. 4 :

 4 Fig. 4: Convergence of the gradient algorithm.