Voronoi Diagrams for a Moderate-Sized Point-Set in a Simple Polygon

Given a set of sites in a simple polygon, a geodesic Voronoi diagram of the sites partitions the polygon into regions based on distances to sites under the geodesic metric. We present algorithms for computing the geodesic nearest-point, higher-order and farthest-point Voronoi diagrams of m point sites in a simple n-gon, which improve the best known ones for m<= n/ polylog n. Moreover, the algorithms for the geodesic nearest-point and farthest-point Voronoi diagrams are optimal for m<= n/ polylog n. This partially answers a question posed by Mitchell in the Handbook of Computational Geometry.


Introduction
The geodesic distance between any two points x and y contained in a simple polygon is the length of the shortest path in the polygon connecting x and y. A geodesic Voronoi diagram of a set S of m sites contained in a simple polygon P partitions P into regions based on distances to sites of S under the geodesic metric. The geodesic nearest-point Voronoi diagram of S partitions P into cells, exactly one cell per site, such that every point in a cell has the same nearest site of S under the geodesic metric. The higher-order Voronoi diagram, also known as the order-k Voronoi diagram, is a generalization of the nearest-point Voronoi diagram. For an integer k with 1 ≤ k ≤ m − 1, the geodesic order-k Voronoi diagram of S partitions P into cells, at most one cell per k-tuple of sites, such that every point in a cell has the same k nearest sites under the geodesic metric. Thus, the geodesic order-1 Voronoi diagram is the geodesic nearest-point Voronoi diagram. The geodesic order-(m − 1) Voronoi diagram of m sites is also called the geodesic farthest-point Voronoi diagram. The geodesic farthest-point Voronoi diagram of S partitions P into cells, at most one cell per site, such that every point in a cell has the same farthest site under the geodesic metric.
In this paper, we study the problem of computing the geodesic nearest-point, higher-order and farthest-point Voronoi diagrams of a set S of m point sites contained in a simple n-gon P . Each edge of a geodesic Voronoi diagram is either a hyperbolic arc or a line segment consisting of points equidistant from two sites under the geodesic metric [2,3,12]. The boundary between any two neighboring cells of a geodesic Voronoi diagram is a chain of O(n) edges. Each end vertex of the boundary is of degree 1 or 3 under the assumption that no point in the plane is equidistant from four distinct sites while every other vertex is of degree 2. There are Θ(k(m − k)) degree-3 vertices in the geodesic order-k Voronoi diagram of S [12]. Every degree-3 vertex is equidistant from three sites and is a point where three Voronoi cells meet. The number of degree-2 vertices is Θ(n) for both the geodesic nearest-point Voronoi diagram and the geodesic farthest-point Voronoi diagram [2,3]. For the geodesic order-k Voronoi diagram, the number of degree-2 vertices is O(kn) [12], but this bound is not tight.
The first nontrivial algorithm for computing the geodesic nearest-point Voronoi diagram was given by Aronov [2] in 1989, which takes O((n + m) log 2 (n + m)) time. Later, Papadopoulou and Lee [16] improved the running time to O((n + m) log(n + m)). However, there has been no progress since then while the best known lower bound of the running time remains to be Ω(n + m log m). In fact, Mitchell posed a question whether this gap can be resolved in the Handbook of Computational Geometry [14,Chapter 27].
For the geodesic order-k Voronoi diagram, the first nontrivial algorithm was given by Liu and Lee [12] in 2013 for polygonal domains with holes. Their algorithm works for m point sites in a polygonal domain with a total of n vertices and takes O(k 2 (n + m) log(n + m)) time. Thus, this algorithm also works for a simple polygon. They presented an asymptotically tight combinatorial complexity of the geodesic order-k Voronoi diagram of m points in a polygonal domain with a total of n vertices, which is Θ(k(m − k) + kn). However, it is not tight for a simple polygon: the geodesic order-(m − 1) Voronoi diagram of m points in a simple n-gon has complexity Θ(n + m) [3]. There is no bound better than the one by Liu and Lee known for the complexity of the geodesic order-k Voronoi diagram in a simple polygon.
For the geodesic farthest-point Voronoi diagram, the first nontrivial algorithm was given by Aronov et al. [3] in 1993, which takes O((n + m) log(n + m)) time. While the best known lower bound is Ω(n + m log m), there has been no progress until Oh et al. [15] presented an O((n + m) log log n)-time algorithm for the special case that all sites are on the boundary of the polygon in 2016. They also claimed that their algorithm can be extended to compute the geodesic farthest-point Voronoi diagram for any m points contained in a simple n-gon in O(n log log n + m log(n + m)) time.
Our results. Our main contributions are the algorithms for computing the nearest-point, higher-order and farthest-point Voronoi diagrams of m point sites in a simple n-gon, which improve the best known ones for m ≤ n/ polylog n. To be specific, we present • an O(n + m log m log 2 n)-time algorithm for the geodesic nearest-point Voronoi diagram, • an O(k 2 m log m log 2 n+min{nk, n(m−k)})-time algorithm for the geodesic order-k Voronoi diagram, and • an O(n+m log m+m log 2 n)-time algorithm for the geodesic farthest-point Voronoi diagram.
Moreover, our algorithms close the gaps of the running times towards the lower bounds. Our algorithm for the geodesic nearest-point Voronoi diagram is optimal for m ≤ n/ log 3 n. Since the algorithm by Papadopoulou and Lee is optimal for m ≥ n, our algorithm together with the one by Papadopoulou and Lee gives the optimal running time for computing the diagram, except for the case that n/ log 3 n < m < n. Similarly, our algorithm for the geodesic farthest-point Voronoi diagram is optimal for m ≤ n/ log 2 n. Since the algorithm by Aronov et al. [3] is optimal for m ≥ n, our algorithm together with the one by Aronov et al. gives the optimal running time for computing the diagram, except for the case that n/ log 2 n < m < n. This answers the question posed by Mitchell on the geodesic nearest-point and farthest-point Voronoi diagrams, except for the short intervals of n/ polylog n < m < n stated above.
For the geodesic order-k Voronoi diagram, we analyze an asymptotically tight combinatorial complexity of the diagram of m points in a simple n-gon, which is Θ(k(m−k)+min{nk, n(m−k)}).
Other contributions of this paper are the algorithms for computing the topological structures of the geodesic nearest-point, order-k and farthest-point Voronoi diagrams which take O(m log m log 2 n), O(k 2 m log m log 2 n) and O(m log m log 2 n) time, respectively. These algorithms allow us to obtain a dynamic data structure for answering nearest or farthest point queries efficiently. In this problem, we are given a static simple n-gon P and a dynamic point set S ⊆ P . We are allowed to insert points to S and delete points from S. After processing updates, we are to find the point of S nearest (or farthest) from a query point efficiently. This data structure supports each query in O( √ m log(n + m)) time and each update in O( √ m log m log 2 n) time, where m is the number of points in S at the moment.

Outline
Our algorithms for computing the geodesic nearest-point, higher-order and farthest-point Voronoi diagrams are based on a polygon-sweep paradigm. For the geodesic nearest-point and higher-order Voronoi diagrams, we fix a point o on the boundary of the polygon and move another point x from o in clockwise order along the boundary of the polygon. While x moves along the boundary, we compute the Voronoi diagram of sites contained in the subpolygon bounded by the shortest path between o and x and the part of the boundary of P from o to x in clockwise order. For the geodesic farthest-point Voronoi diagram, we sweep the polygon with a curve consisting of points equidistant from the geodesic center of the sites. The curve moves from the boundary towards the geodesic center. During the sweep, we gradually compute the diagram restricted to the region the curve has swept. To achieve algorithms running faster than the best known ones for m ≤ n/ polylog n, we first compute the topological structure of a diagram instead of computing the diagram itself directly. The topological structure, which will be defined later, represents the adjacency of the Voronoi cells and has complexity smaller than the complexity of the Voronoi diagram. Once we have the topological structure of a Voronoi diagram, we can compute the Voronoi diagram in O(T 1 + T 2 log n) time, where T 1 denotes the complexity of the Voronoi diagram and T 2 denotes the complexity of the topological structure of the diagram.
We define four types of events where the topological structure of the diagram changes. To handle each event, we compute a point equidistant from three points under the geodesic metric. There is no algorithm known for computing a point equidistant from three points efficiently, except an O(n)-time trivial algorithm. We present an O(log 2 n)-time algorithm assuming that the data structure by Guibas and Hershberger [9] is constructed for P . To obtain this algorithm, we apply two-level binary search on the regions of a subdivision of the polygon. This algorithm allows us to handle each event in O(polylog{n, m}) time.
One application of the topological structure of a diagram is a data structure for nearest (or farthest) point queries for a dynamic point set. To obtain this data structure, we apply the framework given by Bentley and Saxe [4] using the algorithm for computing the topological structure of the geodesic nearest-point (or farthest-point) Voronoi diagram. We subdivide the dynamic point set into √ m almost equal-sized subsets, where m is the number of the input point. Then compute the topological structure of the diagram for each subset. We observe that we can find the Voronoi cell of each diagram containing a query point in O(log(n + m)) time once we have the topological structure of the diagram, which leads to the query time of O( √ m log(n + m)).

Preliminaries
Let P be a simple polygon with n vertices and S be a set of m points contained in P . For ease of description, we use VD[S], k-VD[S] and FVD[S] (or simply VD, k-VD and FVD if they are understood in the context) to denote the geodesic nearest-point, order-k and farthest-point Voronoi diagrams of S in P , respectively. We assume the general position condition that no vertex of P is equidistant from two distinct sites of S and no point of P is equidistant from four distinct sites of S. This was also assumed by in previous work [2,3,12,16] on geodesic Voronoi diagrams. Consider any three points x, y and z in P . We use π(x, y) to denote the shortest path (geodesic path) between x and y contained in P , and d(x, y) to denote the geodesic distance between x and y. Two geodesic paths π(x, y) and π(x, z) do not cross each other, but may overlap with each other. We call a point x the junction of π(x, y) and π(x, z) if π(x, x ) is the maximal common path of π(x, y) and π(x, z). Refer to Figure 1(a).
Consider the set B of points q of P satisfying d(x, q) = d(y, q) for any two points x and y in P . Since x and y are not necessarily contained in S, the set B may contain a two-dimensional region if there is a vertex v of P satisfying d(x, v) = d(y, v) [3]. However, there are at most two such two-dimensional regions (including their boundaries) in B and the other points of B form a simple curve that is incident to the regions by the general position assumption. We call the curve the bisecting curve of x and y and denote it by b(x, y).
Given a point p ∈ P and a closed set A ⊆ P , we slightly abuse the notation π(p, A) to denote the shortest path contained in P connecting p and a point in A. Similarly, we abuse the notation d(p, A) to denote the length of π(p, A). It holds that d(p, A) ≤ d(p, q) + d(q, A) for any two points p, q ∈ P and any closed set A ⊆ P .
We say a set A ⊆ P is geodesically convex if π(x, y) ⊆ A for any two points x and y in A. The geodesic convex hull of S is the intersection of all geodesic convex sets in P that contain S. The geodesic convex hull of a set of m points in P can be computed in O(n + m log(n + m)) time [9]. The geodesic center of a simple polygon P is the point c ∈ P that minimizes max p∈P d(c, p). The center is unique [17] and can be computed in O(n) time [1]. Similarly, the geodesic center of S can be defined as the point c ∈ P that minimizes max s∈S d(c, s). It is known that the geodesic center of a set S of m points in P coincides with the geodesic center of the geodesic convex hull of S [3]. Therefore, we can compute the geodesic center of S by computing the geodesic convex hull of S and its geodesic center. This takes O(n + m log(n + m)) time in total.

Computing the Geodesic Center of Points in a Simple Polygon
We first present an O(log n)-time algorithm for computing the geodesic center of three points contained in P , assuming that we have the data structure of Guibas and Hershberger [9,10]. Using ideas from this algorithm, we present an O(m log m log 2 n)-time algorithm for computing the geodesic center of m points in Section 3.2. These algorithms will be used as subprocedures for computing the Voronoi diagrams of points in P .

Computing the Geodesic Center of Three Points
Let p 1 , p 2 and p 3 be three points in P , and let c be the geodesic center of them. The geodesic convex hull of p 1 , p 2 , p 3 is bounded by π(p 1 , p 2 ), π(p 2 , p 3 ), and π(p 3 , p 1 ). The geodesic convex hull may have complexity Ω(n), but its interior is bounded by at most three concave chains. This allows us to compute the geodesic center of it efficiently.
We first construct the data structure of Guibas and Hershberger [9,10] for P that supports the geodesic distance query between any two points in O(log n) time. To compute c, we compute the shortest paths π(p 1 , p 2 ), π(p 2 , p 3 ), and π(p 3 , p 1 ). Each shortest path has a linear size, but we can compute them in O(log n) time using the data structure of Guibas and Hershberger. Then we find a convex t-gon with t ≤ 6 containing c such that the geodesic path π(x, p i ) has the same combinatorial structure for any point x in the t-gon for each i = 1, 2, 3. To find such a convex t-gon, we apply two-level binary search. Then we can compute c directly in constant time inside the t-gon.
The data structure given by Guibas and Hershberger. Guibas and Hershberger [9,10] gave a data structure of linear size that enables us to compute the geodesic distance between any two query points lying inside P in O(log n) time. We call this structure the shortest path data structure. This data structure can be constructed in O(n) time.
In the preprocessing, they compute a number of shortest paths such that for any two points p and q in P , the shortest path π(p, q) consists of O(log n) subchains of precomputed shortest paths and O(log n) additional edges that connect the subchains into one. In the query algorithm, they find such subchains and edges connecting them in O(log n) time. Then the query algorithm returns the shortest path between two query points represented as a binary tree of height O(log n) [10]. Therefore, we can apply binary search on the vertices of the shortest path between any two points.
Computing the geodesic center of three points: two-level binary search. Let be the geodesic convex hull of p 1 , p 2 and p 3 . The geodesic center c of the three points is the geodesic center of [3], thus is contained in . If the center lies on the boundary of , we can compute it in O(log n) time since it is the midpoint of two points from p 1 , p 2 and p 3 . So, we assume that the center lies in the interior of . Let p i be the junction of π(p i , p j ) and π(p i , p k ) for three distinct indices i, j and k in {1, 2, 3}. See Figure 1(a).
We use the following lemmas to apply two-level binary search. Recall that we have already constructed the shortest path data structure for P .

Lemma 2 ([6]). Given a point p ∈
and a direction, we can find the first intersection point of the boundary of with the ray from p in the direction in O(log n) time.
Proof. Chazelle et al. [6] showed that the first intersection point can be found O(log n) time if we have balanced binary search trees representing the maximal concave curves lying on the boundary of . We can obtain such balanced binary search trees from the shortest path data structure in O(log n) time.
The first level. Imagine that we subdivide into O(n) regions with respect to p 1 by extending the edges of π(p 1 , p 2 ) ∪ π(p 1 , p 3 ) towards π(p 2 , p 3 ). See Figure 1(b). The extensions of the edges can be sorted in the order of their endpoints appearing along π(p 2 , p 3 ). Consider the subdivision of by the extensions, and assume that we can determine which side of a given extension in contains c in T (n) time. Then we can compute the region of the subdivision containing c in O(T (n) log n) time by applying binary search on the extensions. Note that any point x in the same region has the same combinatorial structure of π(x, p 1 ) (and π(x, p 1 )).
We also do this for p 2 and p 3 . Then we have three regions whose intersection contains c. Let D be the intersection of these three regions. We can find D in constant time by the following lemma.
Lemma 3. The intersection D is a convex polygon with at most six edges from extensions of the regions.
Proof. Let C i be the region of the subdivision with respect to p i containing c for each i = 1, 2, 3.
The boundary of C i consists of two extensions and a part of π(p j , p k ), where j and k are distinct indices in {1, 2, 3} \ {i}. This means that C i does not contain any concave curve which comes from π(p i , p j ) or π(p i , p k ). Therefore, the intersection of the three regions is a convex polygon with at most six edges.
We do not subdivide explicitly. Because we have π(p 1 , p 2 ) and π(p 1 , p 3 ) in binary trees of height O(log n), we can apply binary search on the extensions of the edges of the geodesic paths without subdividing explicitly. In this case, during the binary search, we compute the extension of a given edge of π(p 1 , p 2 ) ∪ π(p 1 , p 3 ) using Lemma 2, which takes O(log n) time.
There is a vertex p on the boundary of such that for any point x contained in D we have d(p 1 , x) = d(p 1 , p) + p − x , where p − x is the Euclidean distance between p and x. Moreover, we already have p from the computation of the region containing c in the subdivision with respect to p 1 . The same holds for p 2 and p 3 . Therefore, we can compute the point c that minimizes the maximum of d(c, p 1 ), d(c, p 2 ) and d(c, p 3 ) in constant time inside D.
Therefore, we have the following lemma. The second level. In the second level binary search, we determine which side of an extension e in contains c. Without loss of generality, we assume that e comes from the subdivision with respect to p 1 . Then π(p 1 , x) has the same combinatorial structure for any point x ∈ e.
This subproblem was also considered in a few previous works on computing the geodesic center of a simple polygon [1,17]. They first compute the point c e in e that minimizes max p∈P d(p, c e ), that is, the geodesic center of the polygon restricted to e. By using c e and its farthest point, Pollack et al. [17] presented a way to decide which side of e contains the geodesic center of the polygon in constant time. However, to compute c e , they spend O(n) time.
In our problem, we can do this in O(log n) time using the fact that the interior of is bounded by at most three concave chains. By this fact, there are two possible cases: c e is an endpoint of e, or c e is equidistant from p 1 and p i for i = 2 or 3. We compute the point on e equidistant from p 1 and p 2 , and the point on e equidistant from p 1 and p 3 . Then we find the point c e among the two points and the two endpoints of e. In the following, we show how to compute the point on e equidistant from p 1 and p 2 if it exists. The point on e equidistant from p 1 and p 3 can be computed analogously.
Observe that e can be subdivided into O(n) disjoint line segments by the extensions of the edges of π(p 2 , v 1 ) ∪ π(p 2 , v 2 ) towards e, where v 1 and v 2 are endpoints of e. See Figure 1(c). For any point x in the same line segment, π(p 2 , x) has the same combinatorial structure.
We first claim that there is at most one point on e equidistant from p 1 and p 2 . Assume to the contrary that there are two such points x and y. Without loss of generality, we assume that d(x, p 1 ) < d(y, p 1 ). By definition, d(x, p 1 ) = d(x, p 2 ) and d(y, p 1 ) = d(y, p 2 ). By the construction of e, d(y, p 1 ) = d(x, p 1 ) + d(x, y) = d(x, p 2 ) + d(x, y), but d(y, p 2 ) < d(x, p 2 ) + d(x, y) by triangle inequality and the construction of e. This contradicts that d(y, p 1 ) = d(y, p 2 ).
Thus we can apply binary search on the line segments in the subdivision of e. As we did before, we do not subdivide e explicitly. Instead, we use the binary trees of height O(log n) representing π(p 2 , v 1 ) and π(p 2 , v 2 ). For a point x in e, by comparing d(p 1 , x) and d(p 2 , x), we can determine which part of x on e contains the point equidistant from p 1 and p 2 in constant time. In this case, we can compute the extension from an edge towards e in constant time since e is a line segment. Thus, we complete the binary search in O(log n) time.
Therefore, we can compute c e in O(log n) time and determine which side of e in contains c in the same time using the method of Pollack et al [17]. The following lemma summarizes this section.
Lemma 5. Given any three points p 1 , p 2 and p 3 contained in a simple n-gon P , the geodesic center of p 1 , p 2 and p 3 can be computed in O(log 2 n) time after the shortest path data structure for P is constructed in linear time.
Remark. The observations in this section together with the tentative prune and search technique [11] yield an O(log n)-time algorithm for computing the geodesic center of any three points. Refer to [11,Section 3.4]. But the running time for computing the geodesic center of three points is subsumed by the overall running times for computing the Voronoi diagrams since the algorithms in Sections 3.1.1 and 3.1.2 take O(log 2 n) time. It seems unclear whether these algorithms can be improved to O(log n) time by applying this technique. Thus we do not provide details of the O(log n)-time algorithm for this problem here.

The Point Equidistant from Three Points
The geodesic center of three points in a simple polygon may not be equidistant from all of them. Moreover, three points in a simple polygon may have no point equidistant from them in the polygon. For example, any three points whose geodesic convex hull is an obtuse triangle contained in a simple polygon have their geodesic center at the midpoint of the longest side of the obtuse triangle, but it is not equidistant from the three points. If the three points are almost aligned along a line, they may have no equidistant point in the polygon.
Under the general position condition on the sites, there is at most one point equidistant from three sites. However, for three points which are not necessarily in S, there may be an infinite number of points equidistant from the three points. In this case, we compute the one closest to the three points.
We can compute the closest equidistant point from any three points efficiently using the algorithm for computing the center of them if any equidistant point exists. Proof. Let p 1 , p 2 and p 3 be three input points and c * be the closest equidistant point from them. We first compute the geodesic center c of the three points in O(log 2 n) time. If c is equidistant from the three points, we are done. Otherwise, it is equidistant from only two of them, say Figure 2: (a) The region (the gray region) in the subdivision with respect to p 1 containing c * is subdivided with respect to p 2 . (b) c * is the intersection point between b(p 1 , p 2 ) and b(p 1 , p 3 ).
p 2 and p 3 , and it lies on π(p 2 , p 3 ). Then we have d(p 1 , c) < d(p 2 , c) = d(p 3 , c). Recall that p i denote the junction of π(p i , p j ) and π(p i , p k ) for three distinct indices i, j and k in {1, 2, 3}. If c lies on π(p 2 , p 2 ), P has no point equidistant from the three points because d(p 1 , p 2 ) < d(p 3 , p 2 ), and therefore d(p 1 , x) < d(p 2 , x) = d(p 3 , x) for any point x on the bisecting curve of p 2 and p 3 . Similarly, if c lies on π(p 3 , p 3 ), P has no point equidistant from the three points. So we assume that c lies in an edge Then v 2 v 3 subdivides P into two subpolygons. See Figure 2(a). Assume that c * exists in P . Then it lies in the subpolygon P of P bounded by v 2 v 3 and not containing p 1 . We claim that the part of π(p 1 , c * ) contained in P is just a line segment. Assume to the contrary that π(p 1 , c * ) ∩ P is a chain of at least two line segments. Then each point q where the chain makes a turn is a vertex of P other than v 2 and v 3 . Then π(p 2 , c * ) or π(p 3 , c * ) also passes through q because q is a vertex of P and π(p i , c * ) do not cross each other for i = 1, 2, 3. This is a contradiction. To see this, observe that c * is an intersection point of b(p 1 , p 2 ), b(p 2 , p 3 ) and b(p 3 , p 1 ). Thus, for any point x in b(p 1 , p 2 ) (or b(p 3 , p 1 )), there is no vertex of P where both π(p 1 , x) and π(p 2 , x) (or π(p 3 , x)) make a turn. Therefore, the part of π(p 1 , c * ) contained in P is a line segment.
Consider the subdivision of P by the extensions of the edges of π(p 1 , v 2 ) ∪ π(p 1 , v 3 ) with respect to p 1 . Figure 2(a) shows a region (gray) in such a subdivision. We find the region C containing c * by applying binary search on the extensions. For an extension e, we determine which side of e in P contains c * as follows. There is at most one point on e equidistant from p 1 and p i for i = 2 or 3. Let c i be the intersection point of b(p 1 , p i ) with e. Figure 2(b) shows b(p 1 , p 2 ) and b(p 1 , p 3 ) contained in P . Each bisecting curve is a simple connected curve and it intersects e at most once by the construction of e.
If c 2 comes after c 3 along e from p 1 , the point c * lies in the side of e containing v 2 . Otherwise, c * lies in the side of e containing v 3 . Thus we can determine which part of e inside P contains c * in constant time after computing c 2 and c 3 in O(log n). To compute the extension from an edge, we can use the ray-shooting algorithm which takes O(log n) time [6] since the endpoints of the extensions lie on the boundary of P . Therefore, we can find the region C containing c * in the subdivision of P with respect to p 1 in O(log 2 n) time.
We let vv 2 and vv 3 be two extensions bounding C. See the gray region in Figure 2(a). By applying binary search on π(v 2 , v 2 ), we find the junction u 2 of π(p 2 , v 2 ) and π(p 2 , c * ) in O(log 2 n) time as follows. We find the point c 2 on the extension of an edge on π(v 2 , v 2 ) that is equidistant The subdivision of ch with respect to t 1 t 2 . (c) If c lies on π(p 2 , t 2 ), we first consider the subdivision P with respect to p 1 and find the region containing c * (the gray region). Then consider the subdivisions of the region with respect to t 1 t 2 and p 2 , find the regions containing c, and compute c directly inside their intersection.
from p 1 and p 2 in O(log n) time. Then we compare d(c 2 , p 1 ) and d(c 2 , p 3 ), which determines whether the junction u 2 lies before the edge from p 2 along π(p 2 , v 2 ) or not. We also do this for p 3 and find the junction We observe that the geodesic path π(p i , c * ) is the concatenation of π(p i , u i ) and the line segment u i c * . Therefore d(p j , ·) for j = 1, 2, 3 is a hyperbolic function with some domain D containing c * . We compute the three hyperbolic functions and find points where the three hyperbolic functions have the same value without considering D. Since we do not consider D, some point that we compute may not be equidistant from p 1 , p 2 and p 3 . We check additionally if each such point is equidistant from the three points. In this way, we can compute c * in O(log 2 n) time in total.

The Point Equidistant from Two Points and a Line Segment
Using a way similar to the one in Section 3.1.1, we can compute the point equidistant from two points in P and a line segment contained in P under the geodesic metric. Given two points and a line segment contained in P , if more than one point of P are equidistant from them, we choose the one closest to them. This will be also used as a subprocedure for computing Voronoi diagrams. Proof. Let p 1 and p 2 be any two points in P , and t 1 t 2 be any line segment contained in P . We compute two points t 1 and t 2 on t 1 t 2 that are closest to p 1 and p 2 in O(log n) time, respectively. We compute the junction p 1 of π(p 1 , p 2 ) and π(p 1 , t 1 ), and the junction p 2 of π(p 2 , p 1 ) and π(p 2 , t 2 ) in O(log n) time. Without loss of generality, we assume that t 1 , t 2 , p 2 and p 1 appear in order along the boundary of the geodesic convex hull ch of them as shown in Figure 3. Note that the interior of ch is connected.
Let c * be the closest equidistant point from p 1 , p 2 and t 1 t 2 in P . We first compute the geodesic center c of them in O(log 2 n) time as follows. If c lies on the boundary of ch, we can compute it in O(log n) time since it is the midpoint of two points from p 1 , p 2 , t 1 and t 2 . Thus we assume that c lies in the interior of ch.
Consider the subdivision of ch with respect to p 1 by the extensions of edges in π(p 1 , t 1 ) ∪ π(p 1 , p 2 ) in direction opposite to p 1 . We can determine which side of a given extension contains c in O(log n) time. Therefore, we can compute the region of the subdivision containing c in O(log 2 n) time without constructing the subdivision as we do for the first level binary search in Section 3.1. Similarly, we can compute the region of the subdivision of ch containing c with respect to p 2 by the extensions of edges in π(p 2 , t 2 ) ∪ π(p 2 , p 1 ) in direction opposite to p 2 in the same time. Now consider the subdivision of ch by the extensions of edges in π(t 1 , p 1 ) ∪ π(t 2 , p 2 ) in direction opposite to t 1 and t 2 . See Figure 3(b). We can compute the region of the subdivision containing c in the same time. Then the intersection of the three regions is a convex polygon with at most six vertices and we can find the geodesic center c in the intersection in constant time.
If c is equidistant from p 1 , p 2 and t 1 t 2 , we are done. Otherwise, it is equidistant from only two of them and it lies in the shortest path connecting the two in P . Let xy be the edge of the shortest path that contains c. If xy is not on the boundary of the interior of ch, then P has no point equidistant from p 1 , p 2 and t 1 t 2 by an argument similar to the one in the first paragraph of the proof of Lemma 6.
Assume that xy is on the boundary of the interior of ch. Then xy subdivides P into two subpolygons and c * , if it exists in P , lies in the subpolygon P not containing the interior of ch. In the case that xy lies on π(p 1 , p 2 ), consider the subdivision of P with respect to t 1 t 2 by the extensions of edges in π(t 1 , x) ∪ π(t 2 , y) towards P . For the other case, that is, xy lies on π(p i , t i ) for i = 1 or 2, consider the subdivision of P with respect to p j by the extension of edges in π(p j , x) ∪ π(p j , y) for j ∈ {1, 2} \ {i}. For an extension e, we determine which side of e in P contains c * as follows. We compute the intersection point q of b(p 1 , p 2 ) with e by performing binary search on the intersections of e by the extensions of edges on the shortest path from p 2 to an endpoint of e. Then we determine which side of e contains c by comparing the geodesic distances d(p 1 , q) and d(t 1 t 2 , q). This can be done in O(log n) time as the intersection point q and the distances can be computed in O(log n) time. See Figure 3(c). In any case, as we do in the proof of Lemma 6, we can find c * in P in O(log 2 n) time if it exists.

The Geodesic Center of Points in a Simple Polygon
Combining the result in the previous subsection with the algorithms for computing the center of points in the plane [13] and for computing the geodesic center of a simple polygon [17], we can compute the geodesic center of a set of m points contained in a simple polygon P with n vertices in O(m log m log 2 n) time after computing the shortest path data structure for P . This algorithm will be used as a subprocedure for computing the topological structure of the geodesic farthest-point Voronoi diagram of points in P .
To compute the center of a simple polygon, Pollack et al. [17] first triangulate the input polygon and construct a balanced binary search tree on the chords of the triangulation using the algorithm by Guibas et al [8]. Then, they find the triangle t of the triangulation containing the geodesic center by applying binary search with a chord-oracle to the balanced binary search tree. A chord-oracle is the procedure for determining which side of a given chord contains the geodesic center. They subdivide t further and locate a smaller triangle t such that the geodesic path from a vertex of P to any point in t has the same combinatorial structure. Finally, they find the center inside t which is the lowest point of the upper envelope of some distance functions within domain t .
To remove the linear dependency of n in the time complexity of our algorithm, we again use the shortest path data structure. We may assume that we already have the triangulation of P and the balanced binary search tree on the chords because the algorithm for constructing the shortest path data structure constructs them.
Computing the triangle of the triangulation containing the geodesic center. Let c be the geodesic center of S. We apply a chord-oracle described in Lemma 8 to O(log n) chords of the triangulation and find the triangle t containing c in O(m log 2 n) time.
Lemma 8. Given a chord of P , we can determine which side of the chord contains the geodesic center in O(m log n) time.
Proof. Let e be a given chord. We will find the point c e on e that minimizes the geodesic distance from its farthest site in S. Based on c e and its farthest sites in S, we can determine which side of e contains c in constant time using the method by Pollack et al [17].
To compute c e , we do the followings. The chord subdivides P into two regions. Let S 1 be the set of sites in S contained in one region of P and S 2 be the set of sites in S contained in the other region of P . We pair the sites in the same set S i for i = 1, 2. Then we have m/2 + 1 pairs. For any pair (s, s ) we have, the bisecting curve of s and s intersects e at most once because both s and s lie in the same side of e. We can prune and search on S with this property.
We compute the intersection p of the bisecting curve of s and s with e for each pair (s, s ) in O(log n) time as follows. Observe that e is subdivided into O(n) disjoint line segments by the extensions with respect to s and by the extensions with respect to s . We observe that p lies between x and x , where x and x are the points on e that are closest to s and s , respectively. We can compute x and x in O(log n) time. We apply binary search on the line segments in the subdivision lying between x and x to find p as follows. Initially, the search space for s is the set of the extensions of edges of π(s, x) ∪ π(s, x ) towards xx , and the search space for s is the set of the extensions of edges of π(s , x) ∪ π(s , x ) towards xx . For each iteration of the binary search, we choose the median of each search space, and let q and q be the intersection of xx with the two medians. We can choose one line segment which does not contain p among xq, xq , x q and x q in constant time based on distances d(s, q) and d(s , q ). This is because one of d(s, ·) and d(s , ·) is increasing and the other is decreasing in the domain of xx . Since at least one of the two search spaces is reduced by a constant fraction in every iteration, we can find p in O(log n) time.
We find the median of the m/2 + 1 intersections by the bisecting curves on e and determine which side of the median contains c e in m · O(log n) time using the convexity of d(s, ·) on e for each site s in S. Let e be the side (line segment) of the median in e that does not contain c e . There are at least m/4 pairs of sites whose intersection points lie on e . For such a pair (s, s ), we have either d(s, z) > d(s , z) or d(s, z) < d(s , z) for every point z ∈ e \ e . Thus we can discard at least one site for each such pair: discard s if d(s, z) < d(s , z), and discard s otherwise. Therefore, we can discard at least m/4 sites in each iteration.
After discarding at least m/4 sites, we update S 1 and S 2 accordingly. Then we pair the remaining sites in the same set and discard some sites repeatedly until there remain only a constant number of candidates of the sites in S farthest from c e . This takes O(m log n) time because the number of sites in S we consider decreases by a constant factor.
For a constant number of candidates, we compute c e directly in O(log n) time as we do in Section 3.1.  Applying binary search, we can find the region in the subdivision of t with respect to s containing c in O(m log 2 n) time. However, in this case, we have m sites. If we do this for every site, it takes O(m 2 log 2 n) time. To do this more efficiently, we apply additional prune and search.
For each site s, we choose the median extension in the subdivision of t with respect to s. Let L be the set of all median extensions from all sites. See Figure 4(b). For a site s, the search space (regions in the subdivision of t with respect to s) can be halved by determining which side of the median extension for s in L contains c. Once we have the region containing c in the arrangement of L, we can determine which side of the median extension for s for every site s. We find the region containing c in the arrangement of L using a (1/r)-cutting for L as follows.
Consider the range space Then X has finite VC-dimension, which can be shown by using a way similar to the one for lines in the plane [5]. Let r be a sufficiently large constant. We compute an (1/r)-net N ⊆ L for Finding the geodesic center inside the smaller triangle. For each site s, the geodesic path π(s, x) between s and any point x in τ has the same combinatorial structure, which means that π(s, x) is a hyperbolic function. Thus, our problem reduces to computing the lowest point of the upper envelope of O(m) functions. Pollack et al. [17] present a procedure for this problem that takes O(m) time.
Therefore, we have the following theorem.

Topological Structures of Voronoi Diagrams
In this section, we define the topological structure of a Voronoi diagram and show how to compute the Voronoi diagram from its topological structure. The topological structure of a Voronoi diagram represents the adjacency of the Voronoi cells.
The common boundary of any two adjacent Voronoi cells is connected for the nearest-point and farthest-point Voronoi diagrams of point sites in a simple polygon [2,3]. Similarly, this also holds for the higher-order Voronoi diagram of point sites in a simple polygon.

Lemma 10.
A Voronoi cell of k-VD is connected for any k with 1 ≤ k ≤ m − 1. Moreover, the common boundary of any two adjacent Voronoi cells in k-VD is connected.
Proof. Assume to the contrary that the Voronoi cell of a k-tuple is not connected. Consider two connected components C a and C b of the Voronoi cell, and let a and b be the two points with a ∈ C a and b ∈ C b that minimize d(a, b). Note that the line tangent to C a at a is orthogonal to the edge of π(a, b) incident to a. See Figure 5(a). Since a is on the boundary of the Voronoi cell, it lies on the bisecting curve of two sites, say s and s . We claim that C a and C b lie in different sides of b(s, s ), which contradicts that C a and C b are contained in the same Voronoi cell of k-VD. For any point x in π(a, s), we have d(x, s) ≤ d(x, s ). Similarly, for any point x in π(a, s ), we have d(x, s) ≥ d(x, s ). Also, consider the extension of the edge of π(a, s) incident to a towards a until it reaches the boundary of P . For any point x in this extension, we have d(x, s) ≥ d(x, s ). Similarly, we have d(x, s) ≤ d(x, s ) for any point x on the extension of the edge of π(a, s ) incident to a towards a until it reaches the boundary of P . This implies that for any point x in π(a, b), we have d(x, s) ≥ d(x, s ). Therefore the claim holds.
For the second part of the lemma, observe that the common boundary of the Voronoi cells of any two distinct k-tuples S k and S k of S is a part of the bisecting curve of two sites s and s with s ∈ S k \ S k and s ∈ S k \ S k . Let p be an endpoint of a connected component Γ of the common boundary. We show that b(s, s ) \ Γ is not contained in the common boundary. Note that p is a degree-3 vertex of k-VD, which is equidistant from s, s and another site, say s . Consider the nearest-point Voronoi diagram of s, s and s . Since the common boundary of any two Voronoi cells is connected for this diagram [2], we have d(s , x) < d(s, x) = d(s , x) for any point x lying on the connected component of b(s, s ) \ Γ with endpoint p. Therefore, the connected component of b(s, s ) \ Γ with endpoint p is not contained in the common boundary. Similarly, we can prove that the other component of b(s, s ) \ Γ is not contained in the common boundary, and thus the lemma holds.
The topological structure of k-VD is defined as follows for 1 ≤ k ≤ m − 1. Imagine that we apply vertex suppression for every degree-2 vertex of the Voronoi diagram while preserving the topology of the Voronoi diagram. Vertex suppression of a vertex v of degree 2 is the operation of removing v (and the edges incident to v) and adding an edge connecting the two neighboring vertices of v. Then the resulting graph consists of vertices of degree-1 and degree-3 and edges connecting the vertices. We call the dual of the resulting graph the adjacency graph of the Voronoi diagram. The adjacency graph is a planar graph with complexity O(k(m − k)), because the number of degree-1 and degree-3 vertices of k-VD is O(k(m − k)) [12]. The adjacency graph of a Voronoi diagram represents the topological structure of the Voronoi diagram.
Assume that we have the adjacency graph of the Voronoi diagram together with the exact positions of the degree-1 and degree-3 vertices of the Voronoi diagram. Each Voronoi cell is defined by k sites, but any two adjacent Voronoi cells share k − 1 sites. Consider two adjacent Voronoi cells V 1 and V 2 . Let s 1 and s 2 be the two sites defining V 1 and V 2 , respectively, which are not shared by them. Then the common boundary of V 1 and V 2 is a simple curve connecting two Voronoi vertices v and u of degree-1 or degree-3 such that each Voronoi edge in the common boundary is a part of b(s 1 , s 2 ) lying between v and u. See Figure 5(b).
To compute the Voronoi edges in the common boundary of V 1 and V 2 , we consider the geodesic paths π(s i , v) and π(s i , u) for i = 1, 2, where s i is the junction of π(s i , v) and π(s i , u). Then for any vertex x in π(s 1 , v) ∪ π(s 1 , u), there exists a point q in b(s 1 , s 2 ) lying between v and u such that π(s 1 , q) and π(s 1 , x) have the same combinatorial structure. The same holds for s 2 . This implies that the number of edges in the common boundary is bounded by the total complexity of π(s 1 , v) ∪ π(s 1 , u) and π(s 2 , v) ∪ π(s 2 , u). Thus, we compute the geodesic paths explicitly and consider every edge of the geodesic paths.
Therefore, we can compute the common boundary of two adjacent Voronoi cells in time linear to its complexity plus O(log n). This leads to O(T 1 + T 2 log n) time for computing the Voronoi diagram from the topological structure of the diagram, where T 1 is the combinatorial complexity of the Voronoi diagram and T 2 is the combinatorial complexity of the adjacency graph.
Lemma 11. We can compute the Voronoi diagram of m points in a simple polygon with n vertices in O(T 1 + T 2 log n) time once the adjacency graph and the exact positions of the degree-1 and degree-3 vertices of the Voronoi diagram are given, where T 1 is the combinatorial complexity of the Voronoi diagram and T 2 is the combinatorial complexity of the adjacency graph.
Therefore, in the following, we focus on computing the adjacency graphs and the exact positions of the degree-1 and degree-3 vertices of VD, k-VD and FVD. In our case, we sweep the polygon with a geodesic path π(o, x) for a fixed point o on the boundary of P and a point x moving along the boundary of P from o in clockwise order. The point x is called the sweep point. If we compute all degree-1, degree-2 and degree-3 vertices of the Voronoi diagram during the sweep, we may not achieve the running time better than O((n + m) log(n + m)) as there are O(n + m) such vertices. The key to improve the running time is to compute the topological structure of the Voronoi diagram first which consists of the degree-1 and degree-3 vertices of the Voronoi diagram and the adjacency of the Voronoi cells. Then we construct the complete Voronoi diagram, including degree-2 vertices, from its topological structure using Lemma 11.

The Geodesic Nearest-Point Voronoi Diagram
Let o be an arbitrary point on ∂P , where ∂P denotes the boundary of P . Consider the sweep point x that moves from o along ∂P in clockwise order. We use P (x) to denote the subpolygon of P bounded by π(o, x) and the part of ∂P from o to x in clockwise order. In other words, P (x) is the region swept by π(o, x). Note that P (x) is weakly simple. See Figure 6. Clearly, P (x 1 ) ⊆ P (x 2 ) for any two points x 1 and x 2 on ∂P such that x 1 comes before x 2 from o in clockwise order along ∂P .
For a site s ∈ P (x), let R s (x) be the set {p ∈ P (x) | d(p, s) ≤ d(p, π(o, x))}. By definition, R s (x 1 ) ⊆ R s (x 2 ) for any two points x 1 and x 2 on ∂P such that x 1 comes before x 2 from o in clockwise order along ∂P .
Lemma 12. If s lies on an edge of π(o, x), R s (x) is a line segment that is incident to s and orthogonal to the edge.
Proof. Let γ s be the boundary of R s (x). By definition, d(p, s) = d(p, π(o, x)) for any point p in γ s . Note that d(p, s) = d(p, π(o, x)) if and only if s is the point of π(o, x) closest to p under the geodesic metric. This is because s lies on π(o, x). Therefore, for any point p in γ s , the edge of π(s, p) incident to s is orthogonal to the edge of π(o, x) containing s.
Moreover, we claim that γ s consists of a single line segment. If π(p, s) consists of more than one line segment for some point p on γ s , we can choose a sufficiently small neighborhood N p of p contained in P such that the point on π(o, x) closest to p is s for any point p ∈ N p . This implies that N p is contained in R(x), which contradicts that γ s is a part of the boundary of R(x).
Lemma 13. For a site s ∈ P (x), R s (x) is connected.
Proof. If s ∈ π(o, x), R s (x) is a line segment by Lemma 12, and therefore it is connected. In the following, we assume that s is not on π(o, x) and show that π(p, s) ⊆ R s (x) for any point p ∈ R s (x). This implies that R s (x) is connected because s is contained in R s (x) if s ∈ P (x) by definition. Let p be a point in R s (x). Consider a point r ∈ π(p, s). We have d(r, s) = d(p, s) − d(p, r). Moreover, we have d(p, π(o, x)) − d(p, r) ≤ d(r, π(o, x)). Since p ∈ R s (x), it holds that d(p, s) ≤ d(p, π(o, x)). Thus, d(r, s) ≤ d(r, π(o, x)), and r is also in R s (x) by definition. Therefore R s (x) is connected.
We say that a subset A of P is weakly monotone with respect to a geodesic path γ if the intersection of π(p, γ) with A is connected for any point p ∈ A. We define a shaft from a point x ∈ P towards a direction to be the line segment connecting x and y, where y is the first intersection point of the boundary of P with the ray from x towards the direction.
Lemma 14. For a site s ∈ P (x) \ π(o, x), the boundary of R s (x) consists of one polygonal chain of ∂P (x) and a simple curve whose both endpoints lie on ∂P (x). Moreover, the simple curve is weakly monotone with respect to π(o, x).
Proof. By Lemma 13, R s (x) is connected. Thus to prove the first part of the lemma, it suffices to show that the boundary of R s (x) intersects ∂P (x). To do this, consider the shaft from a point p ∈ R s (x) in direction opposite to the edge of π(p, π(o, x)) incident to p. We claim that the shaft is contained in R s (x). This is because d(r, π(o, x)) = d(p, π(o, x))+d(p, r) ≥ d(p, s)+d(p, r) ≥ d(r, s) for any point r in the shaft by triangle inequality and the fact that R s (x) contains no point of π(o, x). Note that an endpoint of the shaft is on ∂P (x). Therefore, R s (x) intersects ∂P (x), and the first part of the lemma holds.
For the second part of the lemma, we claim that π(p, π(o, x)) \ {p} does not intersect R s (x) for any point p on ∂R s (x) \ ∂P (x). Assume to the contrary that a point p in the path is in R s (x). We already showed that the shaft from p in direction opposite to the edge of π(p , π(o, x)) incident to p is contained in R s (x). Note that the shaft contains p. This contradicts that p lies on the boundary of R s (x). Therefore, the second part of the lemma also holds. Now we consider the union of R s (x) for every site s in P (x) and denote it by R(x). Figure 6 shows three sites contained in P (x) for a simple polygon P . The dashed region in Figure 6(a) is R(x). Note that for any point p ∈ P (x) and any site s lying outside of P (x), it holds that d(p, π(o, x)) < d(p, s). This implies that for any point p ∈ R(x), the nearest site of p is in P (x). Therefore, for computing VD restricted to R(x), we do not need to consider any sites lying outside of P (x).
The following corollaries follow from the properties of R s (x).
Corollary 15. The closure of P (x) \ R(x) is weakly simple.

Corollary 16. Each connected component of R(x)
consists of a polygonal chain of ∂P and a simple curve with endpoints on ∂P unless π(o, x) contains a site. Moreover, the union of such curves is weakly monotone with respect to π(o, x) at any time.
We maintain the geodesic nearest-point Voronoi diagram restricted to R(x) of the sites contained in P (x) while x moves along ∂P . However, after the sweep, R(x) is still a proper subpolygon P , and therefore VD is not completed yet. To resolve this, we attach a long and very thin triangle to P to make a bit larger simple polygon P and set o to lie at the tip of the triangle, as illustrated in Figure 6, so that once we finish the sweep (when the sweep point returns back to o) in P we have the complete Voronoi diagram in P . (The triangle has height of the diameter of P .) The beach line and the breakpoints. There are O(m) connected components of R(x) each of whose boundary consists of a polygonal chain of ∂P and a connected simple curve. The beach line is defined to be the union of the O(m) simple curves of R(x). The beach line has properties similar to those of the beach line of the Euclidean Voronoi diagram. It consists of O(n + m) hyperbolic or linear arcs. We do not maintain them explicitly because the complexity of the sequence is too large for our purpose. Instead, we maintain the combinatorial structure of the beach line using O(m) space as follows.
A point on the beach line is called a breakpoint if it is equidistant from two distinct sites. Each breakpoint moves and traces out a Voronoi edge as the sweep point moves along ∂P . We represent each breakpoint symbolically. That is, a breakpoint is represented as the pair of sites equidistant from it. We say that such a pair defines the breakpoint. Given the pair defining a breakpoint and the position of the sweep point, we can find the exact position of the breakpoint using Lemma 7 in O(log 2 n) time.
Additionally, we consider the endpoints of the O(m) simple curves comprising the beach line. We simply call them the endpoints of the beach line. For an endpoint p of the beach line, there is a unique site s satisfying d(p, s) = d(p, π(o, x)). We say that s defines p. While maintaining B(x), we compute the adjacency graph of VD together with the degree-1 and degree-3 vertices of VD restricted to R(x). Specifically, when a new breakpoint defined by a pair (s, s ) of sites is added to B(x), we add an edge connecting the node for s and the node for s into the adjacency graph.
Events. We have four types of events: site events, circle events, vanishing events, and merging events. Every event corresponds to a key, which is a point on the boundary of P . We maintain the events with respect to their keys sorted in clockwise order from o along the boundary. The event occurs when the sweep point x passes through the corresponding key. The sequence B(x) changes only when x passes through the key of an event. Given a sorted sequence of ν events for ν ∈ N, we can insert a new event in O(log ν) time since we are given each point on the boundary of P together with the edge of P where it lies. We can delete an event from ν in O(log |ν|) time, and find the first event in ν and delete it in O(1) time.
The definitions of the first two event types are similar to the ones in Fortune's algorithm [7]. Each site s in S defines a site event. The key k of the site event defined by s is the point on ∂P that comes first from o in clockwise order along ∂P among the points x with s ∈ π(o, x ). An event of the other three event types is defined by a pair of consecutive points in B(x) or a single breakpoint in B(x). An event is said to be valid if the two points defining the event are consecutive, or the breakpoint defining the event is in B(·). Before the sweep point reaches the key of an event, the two points defining the event may become non-consecutive, or the breakpoint defining the event may disappear from B(x) due to the changes of B(·). In this case, we say that the event become invalid.
A pair (β 1 , β 2 ) of consecutive breakpoints in B(x) defines a circle event if there is a point c equidistant from s 1 , s 2 and s 3 under the geodesic metric, where (s 1 , s 2 ) and (s 2 , s 3 ) are two pairs of sites defining β 1 and β 2 , respectively. The key k of this circle event is the point on ∂P that comes first from o in clockwise order along ∂P among the points x with d(c, π(o, x )) = d(c, s 1 )(= d(c, s 2 ) = d(c, s 3 )). If this event is valid when x passes through k, c appears on the beach line at the time. Moreover, β 1 and β 2 disappear from the beach line, and a new breakpoint defined by (s 1 , s 3 ) appears on the beach line. (One may regard this as merging β 1 and β 2 into the breakpoint defined by (s 1 , s 3 ).) Each breakpoint β in B(x) defines a vanishing event. Let (s 1 , s 2 ) be the pair of sites defining β. Consider the two points on ∂P equidistant from s 1 and s 2 . We observe that exactly one of them lies outside of R(x). We denote it by p. See Figure 6(c). The key k of the vanishing event is the point on ∂P that comes first from o in clockwise order along ∂P among the points x with d(p, π(o, x )) = d(p, s 1 ). Assume that this event is valid when x passes through k. Then, β traces out a Voronoi edge and reaches p, which is a degree-1 Voronoi vertex. Moreover, B(x) changes accordingly as follows. Right before x reaches k, an endpoint defined by s 1 or s 2 is a neighbor of β in B(x). (Otherwise, the beach line is not weakly monotone.) Without loss of generality, we assume that an endpoint defined by s 1 is a neighbor of β. When x reaches the key, this endpoint and β disappear from B(x), and an endpoint defined by s 2 appears in B(x). (One may regard this as merging the endpoint defined by s 1 and β into the endpoint defined by s 2 .) A pair (β 1 , β 2 ) of consecutive endpoints in B(x) defines a merging event if β 1 and β 2 are endpoints of two distinct connected curves of the beach line. Let s 1 and s 2 be the sites defining β 1 and β 2 , respectively. Without loss of generality, assume that β 1 comes before β 2 from o in clockwise order along ∂P . Let p be the (unique) point on ∂P equidistant from s 1 , s 2 that comes after β 1 and before β 2 from o in clockwise order along ∂P . The key k of the merging event is the point on ∂P that comes first from o in clockwise order along ∂P among the points x with d(p, π(o, x )) = d(p, s 1 ). Assume that this event is valid when x passes through k. Then, as the sweep point x moves along ∂P , β 1 and β 2 become closer and finally meet each other at p. This means that the two connected curves, one containing β 1 and the other containing β 2 , merge into one connected curve. At this time, β 1 and β 2 merge into a breakpoint defined by (s 1 , s 2 ).

An Algorithm
Initially, B(x) = B(o) is empty, so there is no circle, vanishing, or merging event. We compute the keys of all site events in advance. For each site, we compute the key corresponding to its event in O(log n) time [8]. As B(x) changes, we obtain new pairs of consecutive breakpoints or endpoints, and new breakpoints. Then we compute the key of each event in O(log 2 n) time using the following lemmas and Lemma 6. In addition, as B(x) changes, some event becomes invalid. Once an event becomes invalid, it does not become valid again. Thus we discard an event if it becomes invalid. Therefore, the number of events we have is O(|B(x)|) at any time.
Lemma 17. Given a point p in P and a distance r ∈ R, the point that comes first from o in clockwise order along ∂P among the points x ∈ ∂P with d(p, π(o, x )) = r can be found in O(log 2 n) time.
Proof. We compute a point x p on ∂P such that π(o, x p ) contains p in O(log n) time [8]. Then d(p, π(o, x p )) = 0 by definition. As the sweep point x moves along ∂P from x p in clockwise order, d(p, π(o, x)) does not decrease. This is because the point p ∈ π(p, π(o, x 2 )) closest to p does not lie in the interior of P (x 1 ) for two points x 1 and x 2 such that x 1 comes before x 2 from o in clockwise order along ∂P . Let x * be the point that comes first from o in clockwise order along ∂P among the points x with d(p, π(o, x )) = r. By the monotonicity of d(p, π(o, ·)), we can apply binary search on the vertices of ∂P to compute the edge of P containing x * . Since we can compute the geodesic distance from p to π(o, x ) for any point x in ∂P in O(log n) time, we can find the edge e of P containing x * in O(log 2 n) time by applying binary search on the vertices of ∂P .
Consider the subdivision of e by the extensions of the edges of the geodesic paths between o and each endpoint of e. The subdivision consists of O(n) disjoint line segments on e. We apply binary search to find the line segment on e containing x * in O(log 2 n) time without computing the subdivision explicitly as we did in Section 3.1. For any point x in the line segment, π(o, x ) has the same combinatorial structure. Thus we can compute x * directly in constant time.
Therefore, we can compute the point on ∂P that comes first from o in clockwise order along ∂P among the points x with d(p, π(o, x )) = r in O(log 2 n) time in total.
Lemma 18. For any two sites s 1 and s 2 in P , we can compute the two points equidistant from s 1 and s 2 under the geodesic metric that lie on the boundary of P in O(log 2 n) time.
Proof. We first compute the balanced binary search tree of π(s 1 , s 2 ) in O(log n) time. Then we extend the edge of π(s 1 , s 2 ) incident to s i towards s i for i = 1, 2 until it reaches a point x i on ∂P . Then x 1 and x 2 partition ∂P into two parts, one from x 1 to x 2 and one from x 2 to x 1 in clockwise order. Note that each part contains one point equidistant from s 1 and s 2 under the geodesic metric. We show how to compute the equidistant point p on the part from x 1 to x 2 in clockwise order in O(log 2 n) time.
To compute the edge of P containing p, we apply binary search on the vertices of ∂P lying from x 1 to x 2 in clockwise order using the property that d(x, s 1 ) ≤ d(x, s 2 ) for any point x on ∂P from x 1 to p in clockwise order, and d(x, s 1 ) ≥ d(x, s 2 ) for any point x on ∂P from p to x 2 in clockwise order. We find the edge e of P containing p in O(log 2 n) time.
Then we consider the subdivision of the edge containing p by the extensions of the edges of π(s 1 , b) ∪ π(s 1 , b ) ∪ π(s 2 , b) ∪ π(s 2 , b ) for two endpoints b and b of the edge e. The subdivision consists of O(n) disjoint line segments on e. By applying binary search without computing the subdivision explicitly, we find the line segment containing p in O(log 2 n) time. Inside the line segment, we find p directly in constant time. Therefore, we can find p and p in total O(log 2 n) time.
Handling site events. To handle a site event defined by a site s with key k, we do the following. By definition, s appears on the beach line when the sweep point x passes through k. By Lemma 12, the part of the beach line defined by a site s is a single line segment when the sweep point passes through k. Therefore, B(k) contains two (degenerate) breakpoints defined by s and some other sites or two (degenerate) endpoints defined by s, which lie at the endpoint of the single line segment other than s. See Figure 6(a) and (b). We find the positions of them in B(x) and update B(x) by adding them using Lemma 19.
Lemma 19. We can obtain B(k) from B(k ) in O(log 2 n log |B(k )|) time, where k is the key previous to k.
Proof. By applying binary search on the breakpoints and the endpoints on B(k ), we compute the part (a single line segment) γ s of the beach line defined by s when the sweep point x passes through k as follows. Let β be a breakpoint or an endpoint of the beach line on B(k ). We determine the position of γ s on the beach line γ with respect to β in O(log 2 n) time. Recall that we maintain β symbolically. We have two sites (or one site if β is an endpoint) defining β, but not the exact position of β. If β is a breakpoint, we compute the exact position of β using Lemma 7. For the case that β is an endpoint, we can compute the exact position of β in O(log 2 n) time in a way similar to the one in Lemma 7. The order of γ s and β along γ is the same as the order of s and the point on π(o, k) closest to β under the geodesic metric by Lemma 16. Thus we can compute the order of γ s and β along γ in O(log 2 n) time in total. Using this property, we apply binary search on the breakpoints and the endpoints on B(k ), and find the position for γ s on the beach line in O(log 2 n log |B(k )|) time.
Based on this result, we compute B(k) by adding two (degenerate) breakpoints or two (degenerate) endpoints to B(k ). In any case, we can compute B(k) in O(log |B(k )|) time. Thus, the total running time is O(log 2 n log |B(k )|).
Adding two breakpoints or two endpoints to B(x) makes a constant number of new pairs of consecutive breakpoints or endpoints, which define circle or merging events. This also makes a constant number of new vanishing events. We compute the key for each such event in O(log 2 n) time and add the events to the event sequence in O(log |B(x)|) time. Recall that the number of events we have is O(|B(x)|) at any time. Therefore, each site event can be handled in O(log 2 n log |B(x)|) time.
Handling the other events. Let k be the key of a circle, vanishing, or merging event. For a circle event, two breakpoints defining the event disappear from the beach line, and a new breakpoint appears. For a vanishing event, the breakpoint defining the event and a neighboring endpoint disappear from the beach line, and a new endpoint appears. For a merging event, two endpoints defining the event are replaced with a new breakpoint. In any case, we can update After updating B(x), we have a constant number of new pairs of consecutive breakpoints or endpoints, and a constant number of new breakpoints. We compute the keys of events defined by them in O(log 2 n) time.
Analysis. For analysis of the correctness, we show that the combinatorial structure of the beach line changes only when the sweep point x passes through an event key.
Lemma 20. The combinatorial structure of the beach line changes only when the sweep point x passes through an event key.
Proof. Consider a breakpoint β on B(x) defined by two sites s 1 and s 2 with the sweep point at a point x on ∂P . If s 1 or s 2 is on π(o, x), x is the key of a site event defined by s 1 or s 2 . In the following, we assume that s 1 , s 2 ∈ P (x) but not on π(o, x). Then β lies on the common boundary of the Voronoi cells of s 1 and s 2 or a degree-3 Voronoi vertex of VD defined by three sites including s 1 and s 2 .
Consider the case that β lies on the common boundary of the Voronoi cells other than its endpoints. We claim that the combinatorial structure does not change in this case. There are two cases: β lies on a Voronoi edge or a degree-2 Voronoi vertex. For the first case, we let e be the edge containing β. For the second case, we let e be the Voronoi edge incident to the degree-2 Voronoi vertex such that d(β, s 1 ) < d(x, s 1 ) for any point x ∈ e. Let U p denote the union of the shafts from p in direction opposite to the edges of π(p, π(o, x)) incident to p for every point p ∈ π(s 1 , β) ∪ π(s 2 , β). Then U p ⊆ R s 1 (x) ∪ R s 2 (x) and U p intersects e at a point other than β unless s 1 ∈ π(β, π(o, x)) or s 2 ∈ π(β, π(o, x)). (But, it is not possible that s 1 , s 2 ∈ π(β, π(o, x)) since d(β, s 1 ) = d(β, s 2 ) = d(β, π(o, x)) and x is not the key of a site event.) By continuity of R s (·) for any site s ∈ P (·), this implies that there is a breakpoint defined by s 1 and s 2 before the sweep point reaches x. Thus, the combinatorial structure does not change in this case.
Consider the case that β lies on a degree-3 Voronoi vertex. We claim that x is the key of a circle event. Let s 3 be the site such that β is defined by s 1 , s 2 and s 3 . We again consider the union U p of the shafts from p in direction opposite to the edges of π(p, π(o, x)) incident to p for every point p ∈ π(s 1 , β) ∪ π(s 2 , β) ∪ π(s 3 , β). Then U p ⊆ R s 1 (x) ∪ R s 2 (x) ∪ R s 3 (x) and U p intersects two Voronoi edges incident to β. This implies that there are two consecutive breakpoints defined by (s 1 , s 2 ) and (s 2 , s 3 ) before the sweep point reaches x. So, x is the key of the circle event defined by this pair of breakpoints.
We can show that the same holds for the case of an endpoint β in a similar way. If β lies on ∂P which is not a Voronoi vertex, the combinatorial structure does not change. If β lies on a degree-1 Voronoi vertex, x is the key of a vanishing event or merging event.
For analysis of the running time, we give upper bounds on the length of B(x) and the total number of valid events we have handled.

The Complexity of the Diagram inside a Simple Polygon
Liu and Lee [12] presented an asymptotically tight complexity of the geodesic higher-order Voronoi diagram of points in a polygonal domain with holes. However, their bound is not tight for a simple polygon. We present an improved upper bound and prove that it is asymptotically tight.
A lower bound. Figure 7 shows an example of which the order-k Voronoi diagram has complexity of Ω(k(m − k) + min{nk, n(m − k)}). Let o be an arbitrary point in the plane. We construct a simple polygon P and a set of sites with respect to the point o. Let 1 be a sufficiently long horizontal line segment whose right endpoint is o and 2 be a sufficiently long line segment with a positive slope whose left endpoint is o. As a part of the boundary of P , we put a concave and y-monotone polygonal curve with n − 4 vertices below 2 such that the highest point of the curve is sufficiently close to u 1 . Then we put another four vertices of P sufficiently far from the sites such that P contains all sites and there is no simple polygon containing more Voronoi vertices than P contains.
Lemma 25. For any κ consecutive sites on 1 with 1 ≤ κ ≤ min{k, m − k}, there exist k − κ sites on 2 such that the κ sites on 1 and the k − κ sites on 2 define a non-empty Voronoi cell.
Proof. Consider κ consecutive sites v i , . . . , v i+κ−1 on 1 . We show that there is a point whose nearest k sites are v i , . . . , v i+κ−1 and u 1 , . . . , u k−κ . To show this, consider the point x equidistant This implies that the nearest k sites from x is v i , . . . , v i+κ−1 and u 1 , . . . , u k−κ . Therefore, the lemma holds.
Lemma 26. For every integer κ with 1 ≤ κ ≤ min{k, m − k}, there are Θ(n) Voronoi edges which come from the bisecting curve of v κ and u k−κ .
Proof. Consider the bisecting curve of v κ and u k−κ . Due to the concave curve with n − 4 vertices, the bisecting curve consists of O(n) hyperbolic arcs and one line segment. We claim that each hyperbolic arc is a Voronoi edge of k-VD. To see this, for any point p in a hyperbolic arc of the bisecting curve, observe that d(v κ , p) ≤ d(v i , p) for any i ≥ κ and d(u k−κ , p) ≤ d(u j , p) for any j ≥ k − κ. This is because p lies below 1 ∪ 2 . Therefore, there are O(n) Voronoi edges from the bisecting curve of v κ and u k−κ .
An upper bound. As Liu and Lee [12] shows, each Voronoi vertex of k-VD has degree 2 or 3, except for the vertices lying on the boundary of P . In the following, we show that k-VD has O(k(m − k) + min{nk, n(m − k)}) Voronoi edges, which implies that the complexity of k-VD is O(k(m − k) + min{nk, n(m − k)}) by the Euler's formula for planar graphs.
A Voronoi edge is a part of the bisecting curve of two sites. Consider Voronoi edges that come from the bisecting curve of two sites s 1 and s 2 in S. Such Voronoi edges form a curve that connects two degree-3 (or degree-1) vertices v and u. Let V (s 1 , s 2 ) is the set of vertices of P in π(s 1 , v) ∪ π(s 1 , u) ∪ π(s 2 , v) ∪ π(s 2 , u) excluding s 1 and s 2 , where s i is the junction of π(s i , u) and π(s i , v) for i = 1, 2. Recall that the number of Voronoi edges on the common boundary of V 1 and V 2 is O(1 + |V (s 1 , s 2 )|). (Refer to Figure 5 For each vertex w of P , we focus on the number of pairs (s 1 , s 2 ) of sites such that w is contained in V (s 1 , s 2 ). Observe that either s 1 or s 2 (but not both) is one of the k nearest sites from w. Without loss of generality, we assume that s 1 is one of the k nearest sites from w. We claim that there is no site s 3 distinct from s 1 such that w ∈ V (s 2 , s 3 ). This is because π(s 2 , x) intersects at most one Voronoi edge defined by s 2 and some other site for any point x in P . We also claim that there is no site s 3 distinct from s 2 such that w ∈ V (s 1 , s 3 ). Assume to the contrary that there is such a site s 3 . We observe that s 1 is one of the k nearest sites of any point in the region R bounded by π(s 1 , v) ∪ π(s 1 , u) and b(s 1 , s 2 ). A Voronoi edge defined by s 1 and s 3 does not intersect R. This means that V (s 1 , s 2 ) and V (s 1 , s 3 ) are interior disjoint. Thus w is the junction of π(s 1 , u ) and π(s 2 , v ), where u and v are the endpoints of the Voronoi edge defined by s 1 and s 2 . Thus it is not contained in V (s 1 , s 3 ), which is a contradiction. By the two claims, the number of pairs (s, s ) with w ∈ V (s, s ) is min{k, m − k} for each vertex w of P .
Therefore, the number of Voronoi edges is O(k(m − k) + min{nk, n(m − k)}), and the total complexity of k-VD is O(k(m − k) + min{nk, n(m − k)}).
Combining the lower bound example, we have the following lemma.
Lemma 27. The geodesic order-k Voronoi diagram of m points contained in a simple polygon with n vertices has complexity of Θ(k(m − k) + min{nk, n(m − k)}).

Computing the Topological Structure of the Diagram
We use the notations defined in Section 5. For VD, we maintain the combinatorial structure of the beach line as x moves along ∂P . For k-VD, we maintain the combinatorial structures of k curves each of which represents a level of the arrangement of some curves. For a site s ∈ P (x), recall that R s (x) = {p ∈ P (x) : d(p, s) ≤ d(p, π(o, x))} is connected. Moreover, the boundary of R s (x) consists of one polygonal chain of ∂P and a simple curve with endpoints on ∂P . We call the simple curve the wave-curve for s. A wave-curve is weakly monotone with respect to π(o, x) by Lemma 14.
We say that a point p ∈ P (x) lies above a curve if π(p, π(o, x)) intersects the curve. We use L i (x) to denote the region of P (x) consisting of all points lying above at least i wave-curves for sites contained in P (x). The ith-level of the arrangement of wave-curves of sites contained in P (x) is defined to be the boundary of L i (x) excluding ∂P . The ith-level for each i is weakly monotone with respect to π(o, x) and consists of at most m simple curves with endpoints lying on ∂P . Note that the beach line for VD coincides with the 1st-level of the arrangement of the wave-curves for all sites in P (x) and R(x) coincides with L 1 (x).
Lemma 28. k-VD[S] restricted to L k (x) coincides with k-VD[S ∩ P (x)] restricted to L k (x).
Proof. Consider a Voronoi cell of k-VD[S] containing a point p ∈ L k (x). We show that this cell is associated with k sites in S ∩ P (x). To see this, observe that d(p, s) > d(p, π(o, x)) for any site s / ∈ P (x). By the definition of L k (x), there are k sites in P (x) whose geodesic distance from p is at most d(p, π(o, x)). Thus the Voronoi cell of k-VD[S] containing p is associated with k sites in S ∩ P (x). Therefore, the lemma holds.
Therefore, we can obtain the topological structure of k-VD[S] by maintaining the topological structure of the kth-level of the arrangement of wave-curves. In addition to this, we maintain the topological structures of the ith-level of the arrangement for all i < k to detect the changes to the topological structure of the kth-level. A breakpoint of the ith-level is an intersection of two wave-curves for each i = 1, . . . , k, and an endpoint of the ith-level is an endpoint of a maximal connected curve in the ith-level. Let B i be the sequence of breakpoints and endpoints of the ith-level, which represents the combinatorial structure of the ith-level.
To avoid the case that R(x) does not contain P after the sweep, we consider a simple polygon P obtained from P by attaching a long and very thin triangle on ∂P as we do in Section 5.
Events. There are four types of events: site events, circle events, vanishing events and merging events. Every event corresponds to a key which is a point on ∂P . An event occurs when the sweep point x passes through its corresponding key. The definitions of the event types are analogous to the ones in Section 5.
Each site s in S defines a site event.
The key of the site event defined by s is the point on ∂P that comes first from o in clockwise order along ∂P among the points x with s ∈ π(o, x ). A pair (β 1 , β 2 ) of consecutive breakpoints in B i defines a circle event if there exists the point c equidistant from three sites defining β 1 or β 2 . The key of this circle event is the point on ∂P that comes first from o in clockwise order along ∂P among the points x with d(c, π(o, x )) = d(c, p), where p is one of the sites defining β 1 or β 2 . Each breakpoint β in B i defines a vanishing event.
Let (s 1 , s 2 ) be the pair defining β. Let p be the point on ∂P equidistant from s 1 and s 2 that lies outside of L i (x). The key k of the vanishing event is the point on ∂P that comes first from o in clockwise order along ∂P among the points x with d(p, π(o, x )) = d(p, s 1 ). A pair (β 1 , β 2 ) of consecutive endpoints in B i defines a merging event if β 1 and β 2 are endpoints of different connected curves of the ith-level. Let s 1 and s 2 be the sites defining β 1 and β 2 , respectively. Let p be the (unique) point on ∂P equidistant from s 1 , s 2 that comes after β 1 and before β 2 from o along ∂P in clockwise order. The key k of the merging event defined by (β 1 , β 2 ) is the point on ∂P that comes first from o in clockwise order along ∂P among the points x with d(p, π(o, x )) = d(p, s 1 ).
The changes of B i due to events are analogous to the ones in Section 5, except for the circle events. The sequence B i may change when the sweep point passes through the key of some circle event defined by a pair of consecutive breakpoints in B i−1 or B i−2 as illustrated in Figure 8.
An algorithm. Initially, B i is empty for all i = 1, . . . , k, so there is no event. We compute all site events in advance. For each site, we can compute its key in O(log n) time [8]. As we handle events, the sequence B i changes accordingly and we reflect these changes to B i . At the same time, we compute the topological structure of k-VD using B k .
To handle a site event defined by a site s, we find the position for s in B i and add s to B i for each i ≤ k. This takes O(log 2 n log |B i |) time for each i by Lemma 19. Adding s to B i makes new pairs of consecutive breakpoints or endpoints in B i , or new breakpoints, which define some events. We compute the key for each event in O(log 2 n) time using Lemma 17.
To handle a circle event defined by a pair (β, β ) of consecutive breakpoints in B i , we consider the following observations, which are analogous to the ones for Euclidean order-k Voronoi diagram given by Zavershynskyi and Papadopoulou [18]. Let (s 1 , s 2 ) and (s 2 , s 3 ) be the pairs of sites defining β and β , respectively. Let c be the point equidistant from s 1 , s 2 and s 3 , and k be the key of the event. When the sweep point x passes through k, two breakpoints merge into a breakpoint defined by (s 1 , s 3 ) in B i . For illustration, see Figure 8. Additionally, B i+1 and B i+2 change. The breakpoints β and β are also in B i+1 before x reaches k. Moreover, a breakpoint β defined by (s 1 , s 2 ) lies between β and β in B i+1 . When x reaches k, the three breakpoints merge into one (degenerate) breakpoint equidistant from s 1 , s 2 and s 3 . After x passes through k, the order of the breakpoints changes: β , β and β. For B i+2 , there is β before x reaches k. After x passes through k, β is replaced with β and β. The remaining levels do not change as x passes through x . We update B i , B i+1 and B i+2 accordingly After updating B i , B i+1 and B i+2 , we have new pairs of consecutive breakpoints or endpoints, or new breakpoints. We compute the events defined by them and add the events to the event queue in O(log 2 n + log |B i | + log |B i+1 | + log |B i+2 |) time. The other events can be handled in a way similar to the one in Section 5 in O(log 2 n + log |B i |) time.
An analysis. We have m site events and spend O( 1≤i≤k log 2 n log |B i |) time to handle each site event, where |B i | is the maximum length of B i during the sweep. We can handle all site events in O(km log 2 n log m) time in total by the following lemma.
Lemma 29. The length of B i is O(im) at any time.
Proof. For i = 1, the lemma holds by Lemma 21. Consider i ≥ 2. Initially, B i is empty. When the sweep point passes through the key of a site event, the length of B i increases by at most two. When the sweep point passes through the key of a circle, vanishing, or merging event caused by B i , the length of B i decreases. In contrast to the nearest-point Voronoi diagram, there is one more case that the length of B i increases in the higher-order Voronoi diagram. A circle event due to B i−1 or B i−2 (or, B i−1 for i = 2) increases the length of B i by at most one. Note that a circle event due to B j that is valid when the sweep point passes through its key corresponds to a degree-3 Voronoi vertex of j-VD. Thus there are O(im) circle events due to B i−1 or B i that are valid when the sweep point passes through their keys. Therefore, the length of B i is O(im), and the lemma holds.
There are O(k 2 m) circle events in total that are valid when the sweep point passes through their keys. This is because each such circle event corresponds to a degree-3 Voronoi vertex of i-VD for some i. There are O(k 2 m) vanishing or merging events in total that are valid when the sweep point passes through their keys. This is because each such event corresponds to a degree-1 Voronoi vertex. Recall that i-VD has O(i(n − i)) degree-1 or degree-3 vertices. Since we can handle each event in O(log m log 2 n) time, the running time for this algorithm is O(k 2 m log m log 2 n) time.
Therefore, we have the following lemma and theorem.
Lemma 30. After constructing the shortest path data structure for P and the shortest path map for P from any fixed point, we can compute k-VD in O(k 2 m log m log 2 n) time using O(km) space (excluding the two data structures).
Theorem 31. Given a set of m points contained in a simple polygon with n vertices, we can compute the geodesic order-k Voronoi of the points in O(k 2 m log m log 2 n + min{nk, n(m − k)}) time using O(n + km) space.

The Geodesic Farthest-Point Voronoi Diagram
In this section, we present two algorithms for computing the topological structure of FVD: an O(n + m log m + m log 2 n)-time algorithm and an O(m log m log 2 n)-time algorithm. The latter algorithm assumes that we have the shortest path data structure for P . The former algorithm is faster than the latter because computing the shortest path data structure for P takes O(n) time. However, the latter algorithm has an advantage if we compute FVD several times with different point sets. For an application, see Section 8.
The algorithm given by Aronov et al. Aronov et al. [3] showed that FVD forms a tree whose root is the geodesic center c of the sites. They first compute c and the geodesic convex hull of the sites. Then they compute FVD restricted to the boundary of the polygon in O((n + m) log(n + m)) time. That is, they compute all degree-1 vertices of FVD.
Then they compute the edges of FVD towards the geodesic center by a reverse geodesic sweep method. They showed that d(u, c) ≤ d(v, c) for any ancestor u of a Voronoi vertex v in FVD (a tree). Thus they compute the Voronoi edges of FVD from the boundary of P toward c one by one. This takes O((n + m) log(n + m)) time.
Computing the geodesic center of the sites and the geodesic convex hull of the sites. We follow the framework of the algorithm by Aronov et al. [3], but we can achieve a faster algorithm for m ≤ n/ log 2 n by applying their approach to compute the topological structure of FVD. We first compute the geodesic center of the sites in S. This takes O(n + m log m) time [1,9], or O(m log m log 2 n) time with the shortest path data structure for P by Theorem 9. The geodesic convex hull of the sites can be computed in O(n + m log m) time [9].
Computing the diagram restricted to the boundary of the polygon We can compute FVD restricted to the boundary of P in O(n + m) time using the algorithm by Oh et al. [15] once we have the sequence of the sites along the boundary of the geodesic convex hull. The algorithm uses the property that the order of the sites along the boundary of the geodesic convex hull of them is the same as the order for their Voronoi cells along the boundary of the polygon, which is shown by Aronov et al [3]. They consider the sites on the boundary of the geodesic convex hull one by one and compute the FVD restricted to the boundary of P along the boundary. Combining their approach with Lemma 18, we can compute FVD restricted to the boundary of P in O(m log 2 n) time.
Extending the diagram inside the polygon We apply the reverse geodesic sweep from the boundary of P towards the center c of the sites. In this problem, a sweep line is a simple curve consisting of points equidistant from c. Let B be the (circular) sequence of the sites that have their Voronoi cells on the sweep line in clockwise order. As an exception, in the initial state, we set B to be the sequence of the sites whose Voronoi cells appear on ∂P .
No site, vanishing, or merging event occurs during the sweep because FVD forms a tree. So we handle the circle event only. Every triple of consecutive sites in B defines a circle event. The key defined by a triple of sites is d(c, c ) for the point c equidistant from the three sites. Thus, we can compute the key of each triple of consecutive sites in B in O(log 2 n) time.
The analysis is similar to the one in Section 5. There are O(m) events in total, and each event can be handled in O(log m + log 2 n) time. Thus extending the diagram inside the polygon takes O(m log m + m log 2 n) time.
Lemma 32. We can compute the topological structure of FVD in O(m log m log 2 n) or O(n + m log m + m log 2 n) time once the shortest path data structure for P is constructed.

Dynamic Data Structures for Nearest or Farthest Point Queries
We showed that the topological structure of a Voronoi diagram can be computed without considering the whole polygon once we have the shortest path data structure for P . The topological structure can also be used for data structures for answering nearest or farthest point queries for a dynamic point set.
We apply the framework given by Bentley and Saxe [4]. At all times, we maintain O( To answer a query, we find the Voronoi cell containing the query point using the topological structure of a Voronoi diagram. To do this, when computing the adjacency graph, we also construct a data structure that supports a point location query as follows.
Point location. The key idea is to approximate the Voronoi diagram into a polygonal subdivision using its adjacency graph and the exact positions of degree-1 and degree-3 vertices. Consider an edge of the adjacency graph corresponding to two adjacent Voronoi cells V 1 and V 2 . We say that a polygonal curve approximates the common boundary of V 1 and V 2 if it connects the two endpoints of the common boundary, consists of at most three line segments, and is contained in the closure of V 1 ∪ V 2 .
Lemma 34. We can compute a polygonal curve approximating the common boundary of two adjacent Voronoi cells in O(log n) time.
Proof. Let u and v be the endpoints of the common boundary of two adjacent Voronoi cells V 1 and V 2 . Let s i be the site associated with V i and s i be the junction of π(s i , u) and π(s i , v) for i = 1, 2. The junction of π(u, s 1 ) and π(u, s 2 ) is u itself, and the junction of π(v, s 1 ) and π(v, s 2 ) is v itself under the general position condition. See figure 9(a) and (b). Therefore, the geodesic convex hull, denoted by ch , of s 1 , s 2 , u and v consists of at most four maximal concave curves. Using this fact, we first find two points t 1 , t 2 such that the geodesic convex hull ch of t 1 , t 2 , u, v is contained in the closure of V 1 ∪ V 2 and the boundary of ch consists of at most four maximal concave polygonal curves. Then we compute a polygonal curve approximating the common boundary in ch.
We first show how to find such points t 1 and t 2 . Consider the case that V 1 and V 2 are two adjacent Voronoi cells in the nearest-point Voronoi diagram. Observe that π(s i , u) ∪ π(s i , v) ⊆ V i for i = 1, 2. Thus, ch is contained in the closure of V 1 ∪ V 2 . Since the boundary of ch consists of at most four concave polygonal curves, π(s i , u), π(s i , v) for i = 1, 2, the two points s 1 and s 2 satisfy the condition for t 1 and t 2 , respectively.
For the case that V 1 and V 2 are two adjacent Voronoi cells in the farthest-point Voronoi diagram, consider two line segments contained in P with endpoint u such that each of them is collinear to the edge of π(s i , u) incident to u for i = 1, 2. Similarly, we consider two line segments contained in P with endpoint v such that each of them is collinear to the edge of π(s i , v) incident to v for i = 1, 2. If u or v is on ∂P , the line segments are u or v itself. It is known that the line segments incident to u are contained in V 1 and the other two line segments are contained in V 2 [3]. The four line segments subdivide P into subpolygons. Let Q denote the subpolygon that contains the common boundary of V 1 and V 2 . If u and v are on ∂P , the subpolygon Q is the whole polygon P . Then Q is contained in the closure of V 1 ∪ V 2 . See Figure 9(b).
Since each polygonal curve approximating the common boundary of two adjacent Voronoi cells described in Lemma 34 is contained in the closure of the union of the Voronoi cells, the following property is satisfied.
Corollary 35. No two polygonal curves constructed by Lemma 34 cross each other.
The following lemma allows us to construct a data structure on a adjacency graph that supports a point location query efficiently. In the lemma, we assume that we have a ray-shooting data structure for P .
Lemma 36. We can construct a data structure for the topological structure of the nearest-point We show how to find the Voronoi cell containing a query point q. Consider the subdivision P A of P with respect to A. Note that the complexity of P A is O( √ m + n). We first find the region in P A containing q as follows without constructing P A explicitly. Let r be the ray from q going upwards. We find the line segment in A that r hits first using the ray-shooting data structure constructed on A in O(log m) time. Then we find the edge of P that r hits first in O(log n) time. By comparing the vertical distances from q to the edge of P and the line segment in A, we can determine the edge of P A that r hits first. Thus we can find the region in P A containing q in O(log(n + m)) time. (If r hits an edge of P before hitting any approximate polygonal curve, we determine which region in P A contains q in O(log(n + m)) time by sorting the degree-1 vertices along ∂P in advance.) Let s be the site associated with the region we just found. Since we consider only the subdivision P A , s might not be the nearest-point (or the farthest-point) from q. In this case, s is contained in a region bounded by a chain γ of Voronoi edges defined by s and s , and a chainγ approximating γ. Note that s is the nearest-point (or the farthest-point) of q. This is because a polygonal curve approximating the common boundary of two Voronoi cells is contained in the union of them, and no two curves constructed from Lemma 34 cross each other.
We show how to find s . For a nearest-point query, the ray from q in direction opposite to the edge of π(s, q) incident to q does not intersect γ. Since q is contained in the region bounded by γ andγ, but the ray does not intersect γ, the ray intersectsγ. We computeγ by finding the point where the ray hits first in P A as we did before, and find the site s definingγ other than s. By comparing d(s, q) and d(s , q), we get the answer in O(log(n + m)) time.
For a farthest-point query, we do in a similar way. In this case, the extension of the edge of π(s, q) incident to q in direction opposite to q does not intersect γ, and thus it intersectsγ. We find the edge of P A that intersects the extension in O(log(n + m)) time, and compute s in a way similar to the one for a nearest-point query. Then by comparing d(s, q) and d(s , q), we get the answer in O(log(n + m)) time.