A linear-time algorithm for the geodesic center of a simple polygon

Given two points in a simple polygon $P$ of $n$ vertices, its geodesic distance is the length of the shortest path that connects them among all paths that stay within $P$. The geodesic center of $P$ is the unique point in $P$ that minimizes the largest geodesic distance to all other points of $P$. In 1989, Pollack, Sharir and Rote [Disc. \&Comput. Geom. 89] showed an $O(n\log n)$-time algorithm that computes the geodesic center of $P$. Since then, a longstanding question has been whether this running time can be improved (explicitly posed by Mitchell [Handbook of Computational Geometry, 2000]). In this paper we affirmatively answer this question and present a linear time algorithm to solve this problem.

distance between them. Given a point x ∈ P , a (geodesic) farthest neighbor of x, is a point f P (x) (or simply f (x)) of P whose geodesic distance to x is maximized. To ease the description, we assume that each vertex of P has a unique farthest neighbor. We can make this general position assumption using simulation of simplicity [9].
Let F P (x) be the function that, for each x ∈ P , maps to the distance to a farthest neighbor of x (i.e., F P (x) = |π(x, f (x))|). A point c P ∈ P that minimizes F P (x) is called the geodesic center of P . Similarly, a point s ∈ P that maximizes F P (x) (together with its farthest neighbor) is called a geodesic diametral pair and their distance is known as the geodesic diameter. Asano and Toussaint [3] showed that the geodesic center is unique (whereas it is easy to see that several geodesic diametral pairs may exist).
In this paper, we show how to compute the geodesic center of P in O(n) time.

Previous Work
Since the early 1980s the problem of computing the geodesic center (and its counterpart, the geodesic diameter) has received a lot of attention from the computational geometry community. Chazelle [7] gave the first algorithm for computing the geodesic diameter (which runs in O(n 2 ) time using linear space). Afterwards, Suri [24] reduced it to O(n log n)-time without increasing the space constraints. Finally, Hershberger and Suri [13] presented a fast matrix search technique, one application of which is a linear-time algorithm for computing the diameter.
The first algorithm for computing the geodesic center was given by Asano and Toussaint [3], and runs in O(n 4 log n)-time. In 1989, Pollack, Sharir, and Rote [22] improved it to O(n log n) time. Since then, it has been an open problem whether the geodesic center can be computed in linear time (indeed, this problem was explicitly posed by Mitchell [19,Chapter 27]).
Several other variations of these two problems have been considered. Indeed, the same problem has been studied under different metrics. Namely, the L 1 geodesic distance [6], the link distance [23,14,8] (where we look for the path with the minimum possible number of bends or links), or even rectilinear link distance [20,21] (a variation of the link distance in which only isothetic segments are allowed). The diameter and center of a simple polygon for both the L 1 and rectilinear link metrics can be computed in linear time (whereas O(n log n) time is needed for the link distance).
Another natural extension is the computation of the diameter and center in polygonal domains (i.e., polygons with one or more holes). Polynomial time algorithms are known for both the diameter [4] and center [5], although the running times are significantly larger (i.e., O(n 7.73 ) and O(n 12+ε ), respectively).

Outline
In order to compute the geodesic center, Pollack et al. [22] introduce a linear time chord-oracle. Given a chord C that splits P into two sub-polygons, the oracle determines which sub-polygon contains c P . Combining this operation with an efficient search on a triangulation of P , Pollack et al. narrow the search of c P within a triangle (and find the center using optimization techniques). Their approach however, does not allow them to reduce the complexity of the problem in each iteration, and hence it runs in Θ(n log n) time.
The general approach of our algorithm described in Section 6 is similar: partition P into O(1) cells, use an oracle to determine which cell contains c P , and recurse within the cell. Our approach differs however in two important aspects that allows us to speed-up the algorithm. First, we do not use the chords of a triangulation of P to partition the problem into cells. We use instead a cutting of a suitable set of chords. Secondly, we compute a set Φ of O(n) functions, each defined in a triangular domain contained in P , such that their upper envelope, φ(x), coincides with F P (x). Thus, we can "ignore" the polygon P and focus only on finding the minimum of the function φ(x).
The search itself uses ε-nets and cutting techniques, which certify that both the size of the cell containing c P and the number of functions of Φ defined in it decrease by a constant fraction (and thus leads to an overall linear time algorithm). This search has however two stopping conditions, (1) reach a subproblem of constant size, or (2) find a triangle containing c P . In the latter case, we show that φ(x) is a convex function when restricted to this triangle. Thus, finding its minimum becomes an optimization problem that we solve in Section 7 using cuttings in R 3 .
The key of this approach lies in the computation of the functions of Φ and their triangular domains. Each function g(x) of Φ is defined in a triangular domain contained in P and is associated to a particular vertex w of P . Intuitively speaking, g(x) maps points in to their (geodesic) distance to w. We guarantee that, for each point x ∈ P , there is one function g defined in a triangle containing x, such that g(x) = F P (x). To compute these triangles and their corresponding functions, we proceed as follows.
In Section 3, we use the matrix search technique introduced by Hershberger and Suri [13] to decompose the boundary of P , denoted by ∂P , into connected edge disjoint chains. Each chain is defined by either (1) a consecutive list of vertices that have the same farthest neighbor v (we say that v is marked if it has such a chain associated to it), or (2) an edge whose endpoints have different farthest neighbors (such edge is called a transition edge).
In Section 4, we consider each transition edge ab of ∂P independently and compute its hourglass. Intuitively, the hourglass of ab, H ab , is the region of P between two chains, the edge ab and the chain of ∂P that contains the farthest neighbors of all points in ab. Inspired by a result of Suri [24], we show that the sum of the complexities of each hourglass defined on a transition edge is O(n). In addition, we provide a new technique to compute all these hourglasses in linear time. In Section 5 we show how to compute the functions in Φ and their respective triangles. We distinguish two cases: (1) Inside each hourglass H ab of a transition edge, we use a technique introduced by Aronov et al. [2] that uses the shortestpath trees of a and b in H ab to decompose H ab into O(|H ab |) triangles with their respective functions (for more information on shortest-path trees refer to [10]).
(2) For each marked vertex v we compute triangles that encode the distance from v. Moreover, we guarantee that these triangles cover every point of P whose farthest neighbor is v. Overall, we compute the O(n) functions of Φ in linear time.

Hourglasses and Funnels
In this section, we introduce the main tools that are going to be used by the algorithm. Some of the results presented in this section have been shown before in different papers.

Hourglasses
Given two points x and y on ∂P , let ∂P (x, y) be the polygonal chain that starts at x and follows the boundary of P clockwise until reaching y.
For any polygonal chain C = ∂P (p 0 , p 1 , . . . , p k ), the hourglass of C, denoted by H C , is the simple polygon contained in P bounded by C, π(p k , f (p 0 )), ∂P (f (p 0 ), f (p k )) and π(f (p k ), p 0 ); see Figure 1. We call C and ∂P (f (p 0 ), f (p k )) the top and bottom chains of H C , respectively, while π(p k , f (p 0 )) and π(f (p k ), p 0 ) are referred to as the walls of H C .
We say that the hourglass H C is open if its walls are vertex disjoint. We say C is a transition chain if f (p 0 ) = f (p k ) and neither f (p 0 ) nor f (p k ) are interior vertices of C. In particular, if an edge ab of ∂P is a transition chain, we say that it is a transition edge (see Figure 1). The following lemma is depicted in Figure 1 and is a direct consequence of the Ordering Lemma proved by Aronov et al. [2,Corollary 2.7.4].
Lemma 2.2. Let C 1 , C 2 , C 3 be three edge disjoint transition chains of ∂P that appear in this order when traversing clockwise the boundary of P . Then, the bottom chains of H C1 , H C2 and H C3 are also edge disjoint and appear in this order when traversing clockwise the boundary of P .
Let γ be a geodesic path joining two points on the boundary of P . We say that γ separates two points x 1 and x 2 of ∂P if the points of X = {x 1 , x 2 } and the endpoints of γ alternate along the boundary of P (x 1 and x 2 could coincide with the endpoints of γ in degenerate cases). We say that a geodesic path γ separates an hourglass H if it separates the points of its top chain from those of its bottom chain. Let Γ = {π(v i , v j ) : 1 ≤ i < j ≤ 4} and note that v 1 , . . . , v 4 split the boundary of P into at most four connected components. If a chain C i is completely contained in one of these components, then one path of Γ separates the top and bottom chain of H Ci . Otherwise, some vertex v j is an interior vertex of C i . However, because the chains C 1 , . . . , C r are edge disjoint, there are at most four chains in this situation. For each chain C i containing a vertex v j , we add the geodesic path connecting the endpoints of C i to Γ. Therefore, Γ consists of O(1) geodesic paths and each hourglass H Ci has its top and bottom chain separated by some path of Γ. Since only O(1) additional paths are computed, this can be done in linear time.
A chord of P is an edge joining two non-adjacent vertices a and b of P such that ab ⊆ P . Therefore, a chord splits P into two sub-polygons. Proof. Note that chords can only appear on walls of hourglasses. Because hourglasses are open, any chord must be an edge on exactly one wall of each of these hourglasses. Assume, for the sake of contradiction, that there exists two points s, t ∈ P whose chord st is in three hourglasses H Ci , H Cj and H C k (for some 1 ≤ i < j < k ≤ r) such that s visited before t when going from the top chain to the bottom one along the walls of the three hourglasses. Let s i and t i be the points in the in the top and bottom chains of H Ci , respectively, such that π(s i , t i ) is the wall of H Ci that contains st (analogously, we define s k and t k ) Because C j lies in between C i and C k , Lemma 2.2 implies that the bottom chain of C j appears between the bottom chains of C i and C k . Therefore, C j lies between s i and s k and the bottom chain of H Cj lies between t i and t k . That is, for each x ∈ C j and each y in the bottom chain of H Cj , the geodesic path π(x, y) is "sandwiched" by the paths π(s i , t i ) and π(s k , t k ). In particular, π(x, y) contains st for each pair of points in the top and bottom chain of H Cj . However, this implies that the hourglass H Cj is not open-a contradiction that comes from assuming that st lies in the wall of three open hourglasses, when this wall is traversed from the top chain to the bottom chain. Analogous arguments can be used to bound the total number of walls that contain the edge st (when traversed in any direction) to O(1). Lemma 2.5. Let x, u, y, v be four vertices of P that appear in this cyclic order in a clockwise traversal of ∂P . Given the shortest-path trees T x and T y of x and y in P , respectively, such that T x and T y can answer lowest common ancestor (LCA) queries in O(1) time, we can compute the path π(u, v) in O(|π(u, v)|) time. Moreover, all edges of π(u, v), except perhaps one, belong to T x ∪ T y .
Proof. Let X (resp. Y ) be the set containing the LCA in T x (resp. T y ) of u, y, and of v, y (resp. u, x and x, y). Note that the points of X ∪ Y lie on the path π(x, y) and can be computed in O(1) time by hypothesis. Moreover, using LCA queries, we can decide their order along the path π(x, y) when traversing it from x to y. (Both X and Y could consist of a single vertex in some degenerate situations). Two cases arise: Case 1. If there is a vertex x * ∈ X lying after a vertex y * ∈ Y along π(x, y), then the path π(u, v) contains the path π(y * , x * ). In this case, the path π(u, v) is the concatenation of the paths π(u, y * ), π(y * , x * ), and π(x * , v) and that the three paths are contained in T x ∪ T y . Moreover, π(u, v) can be computed in time proportional to its length by traversing along the corresponding tree; see Figure 2 (top).
Case 2. In this case the vertices of X appear before the vertices of Y along π(x, y). Let x (resp. y ) be the vertex of X (resp. Y ) closest to x (resp. y).
Let u be the last vertex of π(u, x) that is also in π(u, y). Note that u can be constructed by walking from u towards x y until the path towards y diverges. Thus, u can be computed in O(|π(u, u )|) time. Define v analogously and compute it in O(|π(v, v )|) time.
Let P be the polygon bounded by the geodesic paths π(x , u ), π(u , y ), π(y , v ) and π(v , x ). Because the vertices of X appear before those of Y along π(x, y), P is a simple polygon; see Figure 2 (bottom). In this case the path π(u, y) is the union of π(u, u ), π(u , v ) and π(v , v). Because π(u, u ) and π(v , v) can be computed in time proportional to their length, it suffices to compute π(u , v ) in O(|π(u , v )|) time.
Note that P is a simple polygon with only four convex vertices x , u , y and v , which are connected by chains of reflex vertices. Thus, the shortest path from x to y can have at most one diagonal edge connecting distinct reflex chains of P . Since the rest of the points in π(u , v ) lie on the boundary of P and from the fact that each edge of P is an edge of T x ∪ T y , we conclude all edges of π(u, v), except perhaps one, belong to T x ∪ T y .
We want to find the common tangent between the reflex paths π(u , x ) and π(v , y ), or the common tangent of π(u , y ) and π(v , x ) as one of them belongs to the shortest path π(u , v ). Assume that the desired tangent lies between the paths π(u , x ) and π(v , y ). Since these paths consist only of reflex vertices, the problem can be reduced to finding the common tangent of two convex polygons. By slightly modifying the linear time algorithm to compute this tangents, we can make it run in O(|π(u , v )|) time.
Since we do not know if the tangent lies between the paths π(u , x ) and π(v , y ), we process the chains π(u , y ) and π(v , x ) in parallel and stop when finding the desired tangent. Consequently, we can compute the path π(u, v) in time proportional to its length. Lemma 2.6. Let P be a simple polygon with n vertices. Given k disjoint Proof. Because the given transition chains are disjoint, Lemma 2.2 implies that the bottom chains of their respective hourglasses are also disjoint. Therefore, the sum of the complexities of all the top and bottom chains of these hourglasses is O(n). To bound the complexity of their walls we use Lemma 2.4. Since no chord is used more than a constant number of times, it suffices to show that the total number of chords used by all these hourglasses is O(n).
To prove this, we use Lemma 2.3 to construct O(1) split chains γ 1 , . . . , γ t such that for each 1 ≤ i ≤ k, there is a split chain γ j that separates the top and bottom chains of H Ci . For each 1 ≤ j ≤ t, let H j = {H Ci : the top and bottom chain of H Ci are separated by γ j }.
Since the complexity of the shortest-path trees of the endpoints of γ j is O(n) [10], and from the fact that the chains C 1 , . . . , C k are disjoint, Lemma 2.5 implies that the total number of edges in all the hourglasses of H j is O(n). Moreover, because each of these edges appears in O(1) hourglasses among C 1 , . . . , C k , we conclude that Since we have only O(1) split chains, our result follows.

Funnels
Let C = (p 0 , . . . , p k ) be a chain of ∂P and let v be a vertex of P not in C. The funnel of v to C, denoted by S v (C), is the simple polygon bounded by C, π(p k , v) and π(v, p 0 ); see Figure 3 (a). Note that the paths π(v, p k ) and π(v, p 0 ) may coincide for a while before splitting into disjoint chains. See Lee and Preparata [15] or Guibas et al. [10] for more details on funnels.
A subset R ⊂ P is geodesically convex if for every x, y ∈ R, the path π(x, y) is contained in R. This funnel S v (C) is also known as the geodesic convex hull of C and v, i.e., the minimum geodesically convex set that contains v and C.
Given two points x, y ∈ P , the (geodesic) bisector of x and y is the set of points contained in P that are equidistant from x and y. This bisector is a curve, contained in P , that consists of circular arcs and hyperbolic arcs. Moreover, this curve intersects ∂P only at its endpoints [1, Lemma 3.22].
The (farthest) Voronoi region of a vertex v of P is the set of points R(v) = {x ∈ P : F P (x) = |π(x, v)|} (including boundary points).
Lemma 2.7. Let v be a vertex of P and let C be a transition chain such Proof. Let a and b be the endpoints of C such that a, b, f (a) and f (b) appear in this order in a clockwise traversal of ∂P . Because R(v) ∩ ∂P ⊂ C, we know that v lies between f (a) and f (b).
Let α (resp. β) be the bisector of v and f (a) (resp. f (b)). Let h a (resp. h b ) be the set of points of P that are farther from v than from f (a) (resp. f (b)). Note that α is the boundary of h a while β bounds h b .
By definition, we know that . By continuity of the geodesic distance, the boundaries of h a ∩ h b and S v (C) must intersect. Because a / ∈ h a and b / ∈ h b , both bisectors α and β must have an endpoint on the edge ab. Since the boundaries of By the triangle inequality and since w cannot be a vertex of P as w intersects ∂P only at its endpoints, we get that

Decomposing the boundary
In this section, we decompose the boundary of P into consecutive vertices that share the same farthest neighbor and edges of P whose endpoints have distinct farthest neighbors.
Using a result from Hershberger and Suri [13], in O(n) time we can compute the farthest neighbor of each vertex of P . Recall that the farthest neighbor of each vertex of P is always a convex vertex of P [3] and is unique by our general position assumption.
We mark the vertices of P that are farthest neighbors of at least one vertex of P . Let M denote the set of marked vertices of P (clearly this set can be computed in O(n) time after applying the result of Hershberger and Suri). In other words, M contains all vertices of P whose Voronoi region contains at least one vertex of P .
Given a vertex v of P , the vertices of P whose farthest neighbor is v appear contiguously along ∂P [2]. Therefore, after computing all these farthest neighbors, we effectively split the boundary into subchains, each associated with a different vertex of M ; see Figure 3 (b).
Let a and b be the endpoints of a transition edge of ∂P such that a appears before b in the clockwise order along ∂P . Because ab is a transition edge, we know that f (a) = f (b). Recall that we have computed f (a) and f (b) in the previous step and note that f (a) appears also before f (b) along this clockwise order. For every vertex v that lies between f (a) and f (b) in the bottom chain  segment (a, b). In other words, the Voronoi region R(v) restricted to ∂P is contained in (a, b).

Building hourglasses
Let E be the set of transition edges of ∂P . Given a transition edge ab ∈ E, we say that H ab is a transition hourglass. In order to construct the triangle cover of P , we construct the transition hourglass of each transition edge of E. By Lemma 2.6, we know that ab∈E |H ab | = O(n). Therefore, our aim is to compute the cover in time proportional to the size of H ab .
By Lemma 2.3 we can compute a set of O(1) separating paths such that for each transition edge ab, the transition hourglass H ab is separated by one (or more) paths in this set. For each endpoint of the O(1) separating paths we compute its shortest-path tree [10]. In addition, we preprocess these trees in linear time to support LCA queries [12]. Both computations need linear time per endpoint and use O(n) space. Since we do this process for a constant number of endpoints, overall this preprocessing takes O(n) time.
Let γ be a separating path whose endpoints are x and y. Note that γ separates the boundary of P into two chains S and S such that S ∪S = ∂P . Let H(γ) be the set of each transition hourglass separated by γ whose transition edge is contained in S (whenever an hourglass is separated by more than one path, we pick one arbitrarily). Note that we can classify all transition hourglasses into the sets H(γ) in O(n) time (since O(1) separating paths are considered).
We claim that we can compute all transition hourglass of H(γ) in O(n) time. By construction, the wall of each of these hourglasses consists of a (geodesic) path that connects a point in S with a point in S . Let u ∈ S and v ∈ S be two vertices such that π(u, v) is the wall of a hourglass in H(γ). Because LCA queries can be answered in O (1)

Covering the polygon with apexed triangles
An apexed triangle = (a, b, c) with apex a is a triangle contained in P with an associated distance function g (x), called the apex function of , such that (1) a is a vertex of P , (2) b, c ∈ ∂P , and (3) there is a vertex w of P , called the definer of , such that In this section, we show how to find a set of O(n) apexed triangles of P such that the upper envelope of their apex functions coincides with F P (x). To this end, we first decompose the transition hourglasses into apexed triangles that encode all the geodesic distance information inside them. For each marked vertex v ∈ M we construct a funnel that contains the Voronoi region of v. We then decompose this funnel into apexed triangles that encode the distance from v.

Inside the transition hourglass
Let ab be a transition edge of P such that b is the clockwise neighbor of a along ∂P . Let B ab denote the bottom chain of H ab after removing its endpoints. As noticed above, a point on ∂P can be farthest from a vertex in B ab only if it lies in the open segment ab. That is, if v is a vertex of B ab such that R(v) = ∅, then R(v) ∩ ∂P ⊂ ab.
In fact, not only this Voronoi region is inside H ab when restricted to the boundary of P , but also R(v) ⊂ H ab . The next result follows trivially from Lemma 2.7.
Our objective is to compute O(|H ab |) apexed triangles that cover H ab , each with its distance function, such that the upper envelope of these apex functions coincides with F P (x) restricted to H ab where it "matters".
The same approach was already used by Pollack et al. in [22,Section 3]. Given a segment contained in the interior of P , they show how to compute a linear number of apexed triangles such that F P (x) coincides with the upper envelope of the corresponding apex functions in the given segment. While the construction we follow is analogous, we use it in the transition hourglass H ab instead of the full polygon P . Therefore, we have to specify what is the relation between the upper envelope of the computed functions and F P (x). We will show that the upper envelope of the apex functions computed in H ab coincides with F P (x) inside the Voronoi region R(v) of every vertex v ∈ B ab .
Let T a and T b be the shortest-path trees in H ab from a and b, respectively. Assume that T a and T b are rooted at a and b, respectively. We can compute these trees in O(|H ab |) time [10]. For each vertex v between f (a) and f (b), let v a and v b be the neighbors of v in the paths π(v, a) and π(v, b), respectively. We say that a vertex v is visible from ab if v a = v b . Note that if a vertex is visible, then the extension of these segments must intersect the top segment ab. Therefore, for each visible vertex v, we obtain a triangle v as shown in Figure 4.
We further split v into a series of triangles with apex at v as follows: Let u be a child of v in either T a or T b . As noted by Pollack et al., v can be of three types, either (1) u is not visible from ab (and is hence a child of v in both T a and T b ); or (2) u is visible from ab, is a child of v only in T b , and v b vu is a left turn; or (3) u is visible from ab, is a child of v only in T a , and v a vu is a right turn.
Let u 1 , . . . , u k−1 be the children of v of type (2)  For each 1 ≤ i ≤ k − 1, extend the segment u i v past v until it intersects ab at a point s i . Let s 0 and s k be the intersections of the extensions of vv a and vv b with the segment ab. We define then k triangles contained in v as follows. For each 0 ≤ i ≤ k − 1, consider the triangle (s i , v, s i+1 ) whose associated apexed (left) function is In a symmetric manner, we define a set of apexed triangles induced by the type (3) children of v and their respective apexed (right) functions. Let g 1 , . . . , g r and 1 , . . . , r respectively be an enumeration of all the generated apex functions and triangles such that g i is defined in the triangle i . Because each function is determined uniquely by a pair of adjacent vertices in T a or in T b , and since these trees have O(|H ab |) vertices, we conclude that r = O(|H ab |).
Note that for each 1 ≤ i ≤ r, the triangle i has two vertices on the segment ab and a third vertex, say a i , called its apex such that for each x ∈ i , g i (x) = |π(x, w i )| for some vertex w i of H ab . We refer to w i as the definer of i . Intuitively, i defines a portion of the geodesic distance function from w i in a constant complexity region.
Lemma 5.2. Given a transition edge ab of P , we can compute a set A ab of O(|H ab |) apexed triangles in O(|H ab |) time with the property that for any point p ∈ P such that f (p) ∈ B ab , there is an apexed triangle ∈ A ab with apex function g and definer equal to f (p) such that 1. p ∈ and 2. g(p) = F P (p).
Proof. Because p ∈ R(f (p)), Lemma 5.1 implies that p ∈ H ab . Consider the path π(p, f (p)) and let v be the neighbor of p along this path. By construction of A ab , there is a triangle ∈ A ab apexed at v with definer w that contains p. The apex function g(x) of encodes the geodesic distance from x to w. Because F P (x) is the upper envelope of all the geodesic functions, we know that g(p) ≤ F P (p).
To prove the other inequality, note that if v = f (p), then trivially g(p) = |pv| + |π(v, w)| ≥ |pv| = |π(p, f (p))| = F P (p). Otherwise, let z be the next vertex after v in the path π(p, f (p)). Three cases arise: (a) If z is invisible from ab, then so is f (p) and hence, (b) If z is a child of type (2), then z plays the role of some child u j of v in the notation used during the construction. In this case: (c) If z is a child of type (3), then analogous arguments hold using the (right) distance d r .
To bound the running time, note that the recursive functions d l , d r and c can be computed in O(|T a | + |T b |) time. Then, for each vertex visible from ab, we can process it in time proportional to its degree in T a and T b . Because the sum of the degrees of all vertices in T a and T b is O(|T a | + |T b |) and from the fact that both |T a | and |T b | are O(|H ab |), we conclude that the total running time to construct A ab is O(|H ab |).
In other words, Lemma 5.2 says that no information on farthest neighbors is lost if we only consider the functions in A ab within H ab .In the next section we use a similar approach to construct a set of apexed triangles (and their corresponding apex functions), so as to encode the distance from the vertices of M .

Inside the funnels of marked vertices
Recall that for each marked vertex v ∈ M , we know at least of one vertex on ∂P such that v is its farthest neighbor. For any marked vertex v, let u 1 , . . . , u k−1 be the vertices of P such that v = f (u i ) and assume that they appear in this order when traversing ∂P clockwise. Let u 0 and u k be the neighbors of u 1 and u k−1 other than u 2 and u k−2 , respectively. Note that both u 0 u 1 and u k−1 u k are transition edges of P . Thus, we can assume that their transition hourglasses have been computed.
Let C v = (u 0 , . . . , u k ) and consider the funnel S v (C v ). We call C v the main chain of S v (C v ) while π(u k , v) and π(v, u 0 ) are referred to as the walls of the funnel. Because v = f (u 1 ) = f (u k−1 ), we know that v is a vertex of both H u0u1 and H u k−1 u k . By definition, we have π(v, u 0 ) ⊂ H u0u1 and π(v, u k ) ⊂ H u k−1 u k . Thus, we can explicitly compute both paths π(v, u 0 ) and π(v, u k ) in O(|H u0u1 | + |H u k−1 u k |) time. So, overall, the funnel S v (C v ) can be constructed in O(k + |H u0u1 | + |H u k−1 u k |) time. Recall that, by Lemma 2.6, the total sum of the complexities of the transition hourglasses is O(n). In particular, we can bound the total time needed to construct the funnels of all marked vertices by O(n).
Since the complexity of the walls of these funnels is bounded by the complexity of the transition hourglasses used to compute them, we get that Since v = f (x), we know that x ∈ R(v) and hence that x ∈ S v (C v ).
We now proceed to split a given funnel into O(|S v (C v )|) apexed triangles that encode the distance function from v. To this end, we compute the shortest-path [11]. We consider the tree T v to be rooted at v and assume that for each node u of this tree we have stored the geodesic distance |π(u, v)|.
Start an Eulerian tour from v walking in a clockwise order of the edges. Let Let w 1 be the first leaf of T v found, and let w 2 and w 3 be the next two vertices visited in the traversal. Two cases arise: Case 1: w 1 , w 2 , w 3 makes a right turn. We define s as the first point hit by the ray apexed at w 2 that shoots in the direction opposite to w 3 .
We claim that w 1 and s lie on the same edge of the boundary of S v (C v ). Otherwise, there would be a vertex u visible from w 2 inside the wedge with apex w 2 spanned by w 1 and w 3 . Note that the first edge of the path π(u, v) is the edge uw 2 . Therefore, uw 2 belongs to the shortest-path T v contradicting the Eulerian order in which the vertices of this tree are visited as u should be visited before w 3 . Thus, s and w 1 lie on the same edge and s can be computed in O(1) time.
At this point, we construct the apexed triangle (w 2 , w 1 , s) apexed at w 2 with apex function We modify tree T v by removing the edge w 1 w 2 and replacing the edge w 3 w 2 by the edge w 3 s; see Figure 5. Case 2: w 1 , w 2 , w 3 makes a left turn and w 1 and w 3 are adjacent, then if w 1 and w 3 lie on the same edge of ∂P , we construct an apexed triangle (w 2 , w 1 , w 3 ) apexed at w 2 with apex function Otherwise, let s be the first point of the boundary of S v (C v ) hit by the ray shooting from w 3 in the direction opposite to w 2 . By the same argument as above, we can show that w 1 and s lie on the same edge of the boundary of S v (C v ) (and thus, we can compute s in O(1) time). We construct an apexed triangle (w 2 , w 1 , s) apexed at w 2 with apex function We modify the tree T v by removing the edge w 1 w 2 and adding the edge w 3 s; see Figure 5 for an illustration. Moreover, for each point x ∈ R(v), there is an apexed triangle with apex function g(x) such that (1) x ∈ and (2) g(x) = F P (x). Proof. The above procedure splits S v (C v ) into apexed triangles, such that their apex function in each of them is defined as the geodesic distance to v. By . Therefore, there is an apexed triangle with apex function g(x) such that x ∈ and g(x) = |π(x, v)| = F P (x). Consequently, we obtain properties (1) and (2).
We now bound the running time of the algorithm. The shortest-path tree [10]. For each leaf of T v we need a constant number of operations to determine in which of the cases we are in (and to treat it as well). Therefore, it suffices to bound the number of times these steps are performed. Note that a leaf is removed from the tree in each iteration. Since the number of leaves strictly decreases each time we are in Case 2, this step cannot happen more than O(|S v (C v )|) times. In Case 1 a new leaf is added if w 1 and w 3 do not lie on the same edge of ∂P . However, the number of leaves that can be added throughout is at most the number of edges of T v . Note that the edges added by either Case 1 or 2 are chords of the polygon and hence do not generate further leaves. Because |T v | = O(|S v (C v )|), we conclude that both Case 1 and 2 are only executed O(|S v (C v )|) times.

Prune and search
With the tools introduced in the previous sections, we can proceed to give the prune and search algorithm to compute the geodesic center. The idea of the algorithm is to partition P into O(1) cells, determine on which cell of P the center lies and recurse on that cell as a new subproblem with smaller complexity.
Naturally, we can discard all apexed triangles that do not intersect the new cell containing the center. Using the properties of the cutting, we can show that both the complexity of the cell containing the center, and the number of apexed triangles that intersect it decrease by a constant fraction in each iteration of the algorithm. This process is then repeated until either of the two objects has constant descriptive size.
Let τ be the set all apexed triangles computed in previous sections. Lemmas 2.6 and 5.4 directly provide a bound on the complexity of τ . Let φ(x) be the upper envelope of the apex functions of every triangle in τ (i.e., φ(x) = max{g i (x) : g i (x) ∈ τ }). The following result is a direct consequence of Lemmas 5.2 and 5.4, and shows that the O(n) apexed triangles of τ not only cover P , but their apex functions suffice to reconstruct the function F P (x).
Lemma 6.2. The functions φ(x) and F P (x) coincide in the domain of points of P , i.e., for each p ∈ P , φ(p) = F P (p).
Given a chord C of P a half-polygon of P is one of the two simple polygons in which C splits P . A 4-cell of P is a simple polygon obtained as the intersection of at most four half-polygons. Because a 4-cell is the intersection of geodesically convex sets, it is also geodesically convex.
Let R be a 4-cell of P and let τ R be the set of apexed triangles of τ that intersect R. Let m R = max{|R|, |τ R |}. Recall that, by construction of the apexed triangles, for each triangle of τ R at least one and at most two of its boundary segments is a chord of P Let C be the set containing all chords that belong to the boundary of a triangle of τ R . Therefore, |τ R | ≤ |C| ≤ 2|τ R |.
To construct an ε-net of C, we need some definitions (for more information on ε-nets refer to [17]). Let ϕ be the set of all open 4-cells of P . For each t ∈ ϕ, let C t = {C ∈ C : C ∩ t = ∅} be the set of chords of C induced by t. Finally, let ϕ C = {C t : t ∈ ϕ} be the family of subsets of C induced by ϕ.
Let ε > 0 (the exact value of ε will be specified later). Consider the range space (C, ϕ C ) defined by C and ϕ C . Because the VC-dimension of this range space is finite, we can compute an ε-net N of (C, ϕ C ) in O(n/ε) = O(n) time [17]. The size of N is O( 1 ε log 1 ε ) = O(1) and its main property is that any 4-cell that does not intersect a chord of N will intersect at most ε|C| chords of C.
Observe that N partitions R into O(1) sub-polygons (not necessarily 4cells). We further refine this partition by performing a 4-cell decomposition. That is, we shoot vertical rays up and down from each endpoint of N , and from the intersection point of any two segments of N , see Figure 6. Overall, this partitions R into O(1) 4-cells such that each either (i) is a convex polygon contained in P of at most four vertices, or otherwise (ii) contains some chain of ∂P . Since |N | = O(1), the whole decomposition can be computed in O(m R ) time (the intersections between segments of N are done in constant time, and for the ray shooting operations we walk along the boundary of R once).
In order to determine which 4-cell contains the geodesic center of P , we extend each edge of a 4-cell to a chord C. This can be done with two rayshooting queries (each of which takes O(m R ) time). We then use the chordoracle from Pollack et al. [22,Section 3] to decide which side of C contains c P . The only requirement of this technique is that the function F P (x) coincides with the upper envelope of the apex functions when restricted to C. Which is true by Lemma 6.2 and from the fact that τ R consists of all the apexed triangles of τ that intersect R.
Because the chord-oracle described by Pollack et al. [22,Section 3] runs in linear time on the number of functions defined on C, we can decide in total O(m R ) time on which side of C the geodesic center of P lies. Since our decomposition into 4-cells has constant complexity, we need to perform O(1) calls to the oracle before determining the 4-cell R that contains the geodesic center of P .
The chord-oracle computes the minimum of F P (x) restricted to the chord before determining the side containing the minimum. In particular, if c P lies on any chord bounding R , then the chord-oracle will find it. Therefore, we can assume that c P lies in the interior of R . Moreover, since N is a ε-net, we know that at most ε|C| chords of C will intersect R .
Using a similar argument, we can show that the complexity of R also decreases: since |C| ≤ 2|τ R | ≤ 2m R , we guarantee that at most 2εm R apexed triangles intersect R . Moreover, each vertex of R is in at least one apexed triangle of τ R by Lemma 6.2, and by construction, each apexed triangle can cover at most three vertices. Thus, by the pigeonhole principle we conclude that R can have at most 6εm R vertices. Thus, if we choose ε = 1/12, we guarantee that both the size of the 4-cell R and the number of apexed triangles in τ R are at most m R /2.
In order to proceed with the algorithm on R recursively, we need to compute the set τ R with the at most ε|C| apexed triangles of τ R that intersect R (i.e., prune the apexed triangles that do not intersect with R ). For each apexed triangle ∈ τ R , we can determine in constant time if it intersects R (either one of the endpoints is in R ∩∂P or the two boundaries have non-empty intersection in the interior of P ). Overall, we need O(m R ) time to compute the at most ε|C| triangles of τ R that intersect R .
By recursing on R , we guarantee that after O(log m R ) iterations, we reduce the size of either τ R or R to constant. In the former case, the minimum of F P (x) can be found by explicitly constructing function φ in O(1) time. In the latter case, we triangulate R and apply the chord-oracle to determine which triangle will contain c P . The details needed to find the minimum of φ(x) inside this triangle are giving the next section. Lemma 6.3. In O(n) time we can find either the geodesic center of P or a triangle containing the geodesic center.

Solving the problem restricted to a triangle
In order to complete the algorithm it remains to show how to find the geodesic center of P for the case in which R is a triangle. If this triangle is in the interior of P , it may happen that several apexed triangles of τ fully contain R . Thus, the pruning technique used in the previous section cannot be further applied. We solve this case with a different approach.
Recall that φ(x) denotes the upper envelope of the apex functions of the triangles in τ , and the geodesic center is the point that minimizes φ. The key observation is that, as it happened with chords, the function φ(x) restricted to R is convex.
Let 1 , 2 , . . . , m be the set of m = O(n) apexed triangles of τ that intersect R . Let g i (x) be the apex function of i such that where a i and w i are the apex and the definer of i , respectively, and κ i = |π(a i , w i )| is a constant. By Lemma 6.2, φ(x) = F P (x). Therefore, the problem of finding the center is equivalent to the following optimization problem in R 3 : (P1). Find a point (x, r) ∈ R 3 minimizing r subject to x ∈ R and Thus, we need only to find the solution to (P1) to find the geodesic center of P . We use some remarks described by Megiddo in order to simplify the description of (P1) [18].
To simplify the formulas, we square the equation |xa i | ≤ r − κ i : And finally for each 1 ≤ i ≤ m, we define the function h i (x, r) as follows: Therefore, our optimization problem can be reformulated as: (P2). Find a point (x, r) ∈ R 3 such that r is minimized subject to x ∈ R and h i (x, r) ≤ 0 and r > max{κ i }, for 1 ≤ i ≤ m.
Let h i (x, r) = x 2 −2x·a i + a i 2 −r 2 +2rκ i −κ 2 i be a function defined in the entire plane and let (P2 ) be an optimization problem analogous to (P2) where every instance of h i (x, r) is replaced by h i (x, r). The optimization (P2 ) was studied by Megiddo in [18]. We provide some of the intuition used by Megiddo to solve this problem.
Although the functions h i (x, r) are not linear inside i , they all have the same non-linear terms. Therefore, for i = j, we get that h i (x, r) = h j (x, r) defines a separating plane As noted by Megiddo [18], this separating plane has the following property: If the solution (x, r) to (P2 ) is known to lie to one side of γ i,j , then we know that one of the constraints is redundant.
Thus, to solve (P2 ) it sufficed to have a side-decision oracle to determine on which side of a plane γ i,j the solution lies. Megiddo showed how to implement this oracle in a way that the running time is proportional to the number of constraints [18].
Once we have such an oracle, Megiddo's problem can be solved using a prune and search similar to that introduced in Section 6: pair the functions arbitrarily, and consider the set of m/2 separating planes defined by these pairs. For some constant r, compute a 1/r-cutting in R 3 of the separating planes. A 1/r-cutting is a partition of the plane into O(r 2 ) convex regions each of which is of constant size and intersects at most m/2r separating planes. A cutting of planes can be computed in linear time in R 3 for any r = O(1) [16]. After computing the cutting, determine in which of the regions the minimum lies by performing O(1) calls to the side-decision oracle. Because at least (r − 1)m/2r separating planes do not intersect this constant size region, for each of them we can discard one of the constraints as it becomes redundant. Repeating this algorithm recursively we obtain a linear running time.
To solve (P2) we follow a similar approach, but our set of separating planes needs to be extended in order to handle apex functions as they are only defined in the same way as in (P2 ) in a triangular domain. Note that the vertices of each apexed triangle that intersect R have their endpoints either outside of R or on its boundary.

Optimization problem in a convex domain
In this section we describe our algorithm to solve the optimization problem (P2). To this end, we pair the apexed triangles arbitrarily to obtain m/2 pairs. By identifying the plane where P lies with the plane Z 0 = {(x, y, z) : z = 0}, we can embed each apexed triangle in R 3 . A plane-set is a set consisting of at most five planes in R 3 . For each pair of apexed triangles ( i , j ) we define a plane-set as follows: For each chord bounding either i or j , consider the line extending this chord and the vertical extrusion of this line in R 3 , i.e., the plane containing this chord orthogonal to Z 0 . Moreover, consider the separating plane γ i,j . The set containing these planes is the plane-set of the pair ( i , j ).
Let Γ be the union of all the plane-sets defined by the m/2 pairs of apexed triangles. Thus, Γ is a set that consists of O(m) planes. Compute an 1/r-cutting of Γ in O(m) time for some constant r to be specified later. Because r is constant, this 1/r-cutting splits the space into O(1) convex regions, each bounded by a constant number of planes [16]. Using a side-decision algorithm (to be specified later), we can determine the region Q of the cutting that contains the solution to (P2). Because Q is the region of a 1/r-cutting of Γ, we know that at most |Γ|/r planes of Γ intersect Q. In particular, at most |Γ|/r plane-sets intersect Q and hence, at least (r − 1)|Γ|/r plane-sets do not intersect Q.
Let ( i , j ) be a pair such that its plane-set does not intersect Q. Let Q be the projection of Q on the plane Z 0 . Because the plane-set of this pair does not intersect Q, we know that Q intersects neither the boundary of i nor that of j . Two cases arise: Case 1. If either i or j does not intersect Q , then we know that their apex function is redundant and we can drop the constraint associated with this apexed triangle.
Case 2. If Q ⊂ i ∩ j , then we need to decide which constrain to drop. To this end, we consider the separating plane γ i,j . Notice that inside the vertical extrusion of i ∩ j (and hence in Q), the plane γ i,j has the property that if we know its side containing the solution, then one of the constraints can be dropped. Since γ i,j does not intersect Q as γ i,j belongs to the plane-set of ( i , j ), we can decide which side of γ i,j contains the solution to (P2) and drop one of the constraints.
Regardless of the case if the plane-set of a pair ( i , j ) does not intersect Q, then we can drop one of its constraints. Since at least (r − 1)|Γ|/r planesets do not intersect Q, we can drop at least (r − 1)|Γ|/r constraints. Because |Γ| ≥ m/2 as each plane-set contains at least one plane, by choosing r = 2, we are able to drop at least |Γ|/2 ≥ m/4 constraints. Consequently, after O(m) time, we are able to drop m/4 apexed triangles. By repeating this process recursively, we end up with a constant size problem in which we can compute the upper envelope of the functions explicitly and find the solution to (P2) using exhaustive search. Thus, the running time of this algorithm is bounded by the recurrence T (m) = T (3m/4) + O(m) which solves to O(m). Because m = O(n), we can find the solution to (P2) in O(n) time.
The last detail is the implementation of the side-decision algorithm. Given a plane γ, we want to decide on which side lies the solution to (P2). To this end, we solve (P2) restricted to γ, i.e., with the additional constraint of (x, r) ∈ γ. This approach was used by Megiddo [18], the idea is to recurse by reducing the dimension of the problem. Another approach is to use a slight modification of the chord-oracle described by Pollack et al. [22,Section 3].
Once the solution to (P2) restricted to γ is known, we can follow the same idea used by Megiddo [18] to find the side of γ containing the global solution to (P2). Intuitively, we find the apex functions that define the minimum restricted to γ. Since φ(x) = F P (x) is locally defined by this functions, we can decide on which side the minimum lies using convexity. We obtain the following result.
Lemma 7.1. Let R be a convex trapezoid contained in P such that R contains the geodesic center of P . Given the set of all apexed triangles of τ that intersect R , we can compute the geodesic center of P in O(n) time.
The following theorem summarizes the result presented in this paper.
Theorem 7.2. We can compute the geodesic center of any simple polygon P of n vertices in O(n) time.