# The Geodesic Farthest-point Voronoi Diagram in a Simple Polygon^{1}^{1}1This work was supported by the NRF grant 2011-0030044 (SRC-GAIA) funded by the government of Korea and the MSIT(Ministry of Science and ICT), Korea, under the SW Starlab support program(IITP-2017-0-00905) supervised by the IITP(Institute for Information & communications Technology Promotion).

###### Abstract

Given a set of point sites in a simple polygon, the geodesic farthest-point Voronoi diagram partitions the polygon into cells, at most one cell per site, such that every point in a cell has the same farthest site with respect to the geodesic metric. We present an -time algorithm to compute the geodesic farthest-point Voronoi diagram of point sites in a simple -gon. This improves the previously best known algorithm by Aronov et al. [Discrete Comput. Geom. 9(3):217-255, 1993]. In the case that all point sites are on the boundary of the simple polygon, we can compute the geodesic farthest-point Voronoi diagram in time.

## 1 Introduction

Let be a simple polygon with vertices. For any two points and in , the geodesic path is the shortest path contained in connecting with . Note that if the line segment connecting with is contained in , then is a line segment. Otherwise, is a polygonal chain whose vertices (other than its endpoints) are reflex vertices of . The geodesic distance between and , denoted by , is the sum of the Euclidean lengths of the line segments in . Throughout this paper, when referring to the distance between two points in , we mean the geodesic distance between them unless otherwise stated. We refer the reader to the survey by Mitchell [12] in the handbook of computational geometry for more information on geodesic paths and distances.

Let be a set of point sites contained in . For a point , a (geodesic) -farthest neighbor of , is a site (or simply ) of that maximizes the geodesic distance to . To ease the description, we assume that every vertex of has a unique -farthest neighbor. This general position condition was also assumed by Aronov et al. [3] and Ahn et al. [2].

The geodesic farthest-point Voronoi diagram of in is a subdivision of into Voronoi cells. Imagine that we decompose into Voronoi cells (or simply ) for each site , where is the set of points in that are closer to than to any other site of . Note that some cells might be empty. The set defines the (farthest) Voronoi tree of with leaves on the boundary of . Each edge of this diagram is either a line segment or a hyperbolic arc [3]. The Voronoi tree together with the set of Voronoi cells defines the geodesic farthest-point Voronoi diagram of (in ), denoted by (or simply if is clear from context). We indistinctively refer to as a tree or as a set of Voronoi cells.

There are similarities between the Euclidean farthest-point Voronoi diagram and the geodesic farthest-point Voronoi diagram (see [3] for further references). In the Euclidean case, a site has a nonempty Voronoi cell if and only if it is extreme, i.e., it lies on the boundary of the convex hull of the set of sites. Moreover, the clockwise sequence of Voronoi cells (at infinity) is the same as the clockwise sequence of sites along the boundary of the convex hull. With these properties, the Euclidean farthest-point Voronoi diagram can be computed in linear time if the convex hull of the sites is known [1]. In the geodesic case, a site with nonempty Voronoi cell lies on the boundary of the geodesic convex hull of the sites. The order of sites along the boundary of the geodesic convex hull is the same as the order of their Voronoi cells along the boundary of . However, the cell of an extreme site may be empty, roughly because the polygon is not large enough for the cell to appear. In addition, the complexity of the bisector between two sites can be linear in the complexity of the polygon.

##### Previous Work.

Since the early 1980s many classical geometric problems have been studied in the geodesic setting. The problem of computing the geodesic diameter of a simple -gon (and its counterpart, the geodesic center) received a lot of attention from the computational geometry community. The geodesic diameter of is the largest possible geodesic distance between any two points in , and the geodesic center of is the point of that minimizes the maximum geodesic distance to the points in .

Chazelle [5] gave the first algorithm for computing the geodesic diameter of , which runs in time using linear space. Suri [17] reduced the time complexity to without increasing the space complexity. Later, Hershberger and Suri [9] presented a fast matrix search technique, one application of which is a linear-time algorithm for computing the diameter of .

The first algorithm for computing the geodesic center of was given by Asano and Toussaint [4], and runs in time. This algorithm computes a super set of the vertices of , where is the set of vertices of . In 1989, Pollack et al. [16] improved the running time to . In a recent paper, Ahn et al. [2] settled the complexity of this problem by presenting a -time algorithm to compute the geodesic center of .

The problem of computing the geodesic farthest-point Voronoi diagram is a generalization of the problems of computing the geodesic center and the geodesic diameter of a simple polygon. For a set of points in , Aronov et al. [3] presented an algorithm to compute in time. While the best known lower bound is , which is a lower bound known for computing the geodesic convex hulls of , it is not known whether or not the dependence on , the complexity of , is linear in the running time. In fact, this problem was explicitly posed by Mitchell [12, Chapter 27] in the Handbook of Computational Geometry.

##### Our Result.

In this paper, we present an -time algorithm for computing for a set of points in a simple -gon. To do this, we present an -time algorithm for the simpler case that all sites are on the boundary of the simple polygon and use it as a subprocedure for the general algorithm.

Our result is the first improvement on the computation of geodesic farthest-point Voronoi diagrams since 1993 [3]. It partially answers the question posed by Mitchell. Moreover, our result suggests that the computation time of Voronoi diagrams has only almost linear dependence in the complexity of the polygon. We believe our results could be used as a stepping stone to solve the question posed by Mitchell [12, Chapter 27]. Indeed, after the preliminary version [14] of this paper had been presented, Oh and Ahn [13] presented an -time algorithm for this problem. They observed that the adjacency graph of the Voronoi cells has complexity smaller than the complexity of the Voronoi diagram and presented an algorithm for the geodesic farthest-point Voronoi diagram based on a polygon-sweep paradigm, which is optimal for a moderate-sized point-set.

##### Outline.

We first assume that the site set is the vertex set of the input simple polygon. Then we present an algorithm for computing , which will be extended to handle the general cases in Section 6 and Section 7. This algorithm consists of three steps. Each section from Section 3 to Section 5 describes a step of the algorithm. In the first step, we compute the geodesic farthest-point Voronoi diagram restricted to the boundary of the polygon. In the second step, we decompose recursively the interior of the polygon into smaller cells, not necessarily Voronoi cells, until the complexity of each cell becomes constant. In the third step, we explicitly compute the geodesic farthest-point Voronoi diagram in each cell and merge them to complete the description of the Voronoi diagram.

In the first step, we compute the restriction of to in linear time, where denotes the boundary of . A similar approach was used by Aronov et al. [3]. However, their algorithm spends time and uses completely different techniques. The main tool used to speed up the algorithm is the matrix search technique introduced by Hershberger and Suri [9] which provides a “partial” description of (i.e., the restriction of to the vertices of .) To extend it to the entire boundary of , we borrow some tools used by Ahn et al. [2]. This reduces the problem to the computation of upper envelopes of distance functions which can be completed in linear time.

In the second step, we recursively subdivide the polygon into cells. To subdivide a cell whose boundary consists of geodesic paths, we construct a closed polygonal path that visits roughly endpoints of the geodesic paths at a regular interval. Intuitively, to choose these endpoints, we start at the endpoint of a geodesic path on the boundary of the cell. Then, we walk along the boundary, choose another endpoint after skipping of them, and repeat this. We consider the geodesic paths, each connecting two consecutive chosen endpoints. The union of all these geodesic paths can be computed in time linear in the complexity of the cell [15] and subdivides the cell into smaller simple polygons. By recursively applying this procedure on each resulting cell, we guarantee that after rounds the boundary of each cell consists of a constant number of geodesic paths. While decomposing the polygon, we also compute restricted to the boundary of each cell. However, the total complexity of restricted to the boundary of each cell might be in the worst case. To resolve this problem, we subdivide each cell further so that the total complexity of restricted to the boundary of each cell is for every iteration. Each round can be completed in linear time, which leads to an overall running time of . After the second step, we have cells in the simple polygon and we call them the base cells.

In the third step, we explicitly compute the geodesic farthest-point Voronoi diagram in each of the base cells by applying the linear-time algorithm of computing the abstract Voronoi diagram given by Klein and Lingas [11]. To apply the algorithm, we define a new distance function for each site whose Voronoi cell intersects the boundary of each cell such that the distance function is continuous on and the total complexity of the distance functions for all sites is . We show that the abstract Voronoi diagram restricted to a base cell is exactly the geodesic farthest-point Voronoi diagram restricted to . After computing the geodesic farthest-point Voronoi diagrams for every base cell, they merge them to complete the description of the Voronoi diagram.

For the case that the sites lie on the boundary of the simple polygon, we cannot apply the matrix searching technique directly although the other procedures still work. To handle this, we apply the matrix search technique with a new distance function to compute restricted to the vertices of . Then we consider the general case that the sites are allowed to lie in the interior of the simple polygon. We subdivide the input simple polygon in a constant number of subpolygons, and apply the previous algorithm for sites on the boundary to these subpolygons. The overall strategy is similar to the one for sites on the boundary, but there are a few nontrivial technical issues to be addressed.

## 2 Preliminaries

For any subset of , let and denote the boundary and the interior of , respectively. For any two points , we use to denote the line segment connecting and . Let be a simple -gon and be a set of point sites contained in . Let be the set of the vertices of . A vertex of a simple polygon is convex (or reflex) if the internal angle at with respect to the simple polygon is less than (or at least) .

For any two points and on , let denote the portion of from to in clockwise order. We say that three (nonempty) disjoint sets and contained in are in clockwise order if for any and any . To ease notation, we say that three points are in clockwise order if and are in clockwise order.

### 2.1 Ordering Lemma

Aronov et al. [3] gave the following lemma which they call
Ordering Lemma. We make use of this lemma to compute restricted to .
Before introducing the lemma, we need to define the geodesic convex hull of a set of points in .
We say a subset of is geodesically convex if for any two points .
The geodesic convex hull of is defined to be the intersection of all geodesic convex sets containing .
It can be computed in time [8]. ^{2}^{2}2The paper [8] shows that
their running time is . But it is . To see this,
observe that it is for . Also, it is for ..

###### Lemma 1 ([3, Corollary 2.7.4]).

The order of sites along is the same as the order of their Voronoi cells along , where is the geodesic convex hull of with respect to .

### 2.2 Apexed Triangles

An apexed triangle with apex is an Euclidean triangle contained in with an associated distance function such that (1) is a vertex of , (2) there is an edge of containing both and , and (3) there is a site of , called the definer of , such that

where denote the Euclidean distance between and .

Intuitively, bounds a constant complexity region where the geodesic distance function from can be obtained by looking only at the distance from . We call the side of an apexed triangle opposite to the apex the bottom side of . Note that the bottom side of is contained in an edge of .

The concept of the apexed triangle was introduced by Ahn et al. [2] and was a key to their linear-time algorithm to compute the geodesic center. After computing the -farthest neighbor of each vertex in linear time [9], they show how to compute apexed triangles in time with the following property: for each point , there exists an apexed triangle such that and . By the definition of the apexed triangle, we have . In other words, the distance from each point of to its -farthest neighbor is encoded in one of the distance functions associated with these apexed triangles.

More generally, we define a set of apexed triangles whose distance functions encode the distances from the points of to their -farthest neighbors. We say a weakly simple polygon is a funnel of a point if its boundary consists of three polygonal curves , and for some two points .

###### Definition 2.

A set of apexed triangles covers if for any site , the union of all apexed triangles with definer is a funnel of such that .

Ahn et al. [2] gave the following lemma. In Section 6 and Section 7, we show that we can extend this lemma to compute a set of apexed triangles covering for any set of points in a simple polygon.

###### Lemma 3 ([2]).

Given a simple -gon with vertex set , we can compute a set of apexed triangles covering in time.

While Lemma 3 is not explicitly stated by Ahn et al. [2], a closer look at the proofs of Lemmas 5.2 and 5.3, and Corollaries 6.1 and 6.2 reveals that this lemma holds. Lemma 3 states that for each vertex of , the set of apexed triangles with definer forms a connected component. In particular, the union of their bottom sides is a connected chain along . Moreover, these apexed triangles are interior disjoint by the definition of apexed triangles.

### 2.3 The Refined Geodesic Farthest-point Voronoi Diagram

Assume that we are given a set of apexed triangles covering . We consider a refined version of which we call the refined geodesic farthest-point Voronoi diagram defined as follows: for each site , the Voronoi cell of is subdivided by the apexed triangles with definer . That is, for each apexed triangle with definer , we define a refined cell , where is the union of and its bottom side (excluding the corners of ). Since any two apexed triangles and with the same definer are interior disjoint, and are also interior disjoint. We denote the set by . Then, forms a tree consisting of arcs and vertices. Notice that each arc of is a part of either the bisector of two sites or a side of an apexed triangle. Since we assume that the number of the apexed triangles is , the complexity of is still . (Lemma 2.8.3 in [3] shows that the complexity of is .)

###### Lemma 4.

For an apexed triangle and a point in , the line segment is contained in , where is the point on the bottom side of hit by the ray from towards .

###### Proof.

Let be a point on . We have . On the other hand, by the triangle inequality for any site . Since for any site other than , we have , which implies that lies in . ∎

###### Corollary 5.

For any site and any point , the line segment is contained in , where is the point on hit by the ray from the neighbor of along towards .

Throughout this paper, we use to denote the number of edges of for a simple polygon . For a curve , we use to denote the number of the refined cells intersecting . For ease of description, we abuse the term ray slightly such that the ray from in a direction denotes the line segment of the halfline from in the direction, where is the first point of encountered along the halfline from .

From Section 3 to Section 5, we will make the assumption that is the set of the vertices of . This assumption is general enough as we show how to extend the result to the case when is an arbitrary set of sites contained in (Section 6) and in (Section 7). The algorithm for computing consists of three steps. Each section from Section 3 to Section 5 describes each step.

## 3 Computing Restricted to

Using the algorithm in [2], we compute a set of apexed triangles covering in time. Recall that the apexed triangles with the same definer are interior disjoint and have their bottom sides on whose union forms a connected chain along . Thus, such apexed triangles can be sorted along with respect to their bottom sides.

###### Lemma 6.

Given a set of all apexed triangles of with definer for a site of , we can sort the apexed triangles in along with respect to their bottom sides in time.

###### Proof.

The bottom side of an apexed triangle is contained in an edge of , and the other two sides are chords of (possibly flush with ). Assume that these chords are oriented from its apex to its bottom side. Using a hash-table storing the chords of the apexed triangles in , we can link each of these chords to its neighboring triangles (and distinguish between left and right neighbors). In this way, we can retrieve a linked list with all the triangles in in sorted order along in time. ∎

### 3.1 Computing the -Farthest Neighbors of the Sites

The following lemma was used by Ahn et al. [2] and is based on the matrix search technique proposed by Hershberger and Suri [9].

###### Lemma 7 ([9]).

We can compute the -farthest neighbor of each vertex of in time.

Using Lemma 7, we mark the vertices of that are -farthest neighbors of at least one vertex of . Let denote the set of marked vertices of . Note that consists of the vertices of each of whose Voronoi region contains at least one vertex of .

We call an edge a transition edge if . Let be a transition edge of such that is the clockwise neighbor of along . Recall that we already have and and note that are in clockwise order by Lemma 1. Let be a vertex of such that are in clockwise order. By Lemma 1, if there is a point on whose -farthest neighbor is , then must lie on . In other words, the Voronoi cell restricted to is contained in and hence, there is no vertex of such that . For a nontransition edge such that , we know that for any point . Therefore, to complete the description of restricted to , it suffices to compute restricted to the transition edges.

### 3.2 Computing Restricted to a Transition Edge

Let be a transition edge of such that is the clockwise neighbor of . Without loss of generality, we assume that is horizontal and lies to the left of . Recall that if there is a site with , then lies in . Thus, to compute , it is sufficient to consider the apexed triangles of with definers in . Let be the set of apexed triangles of with definers in .

We give a procedure to compute in time using the sorted lists of the apexed triangles with definers in . Once it is done for all transition edges, we obtain the refined geodesic farthest-point Voronoi diagram restricted to in time. Let be the sites lying on in counterclockwise order along . See Figure 1.

#### 3.2.1 Upper Envelopes and

Consider any functions with for , where is a subset of . We define the upper envelope of ’s as the piecewise maximum of ’s. Moreover, we say that a function appears on the upper envelope if and at some point for any other functions .

Each apexed triangle has a distance function such that for a point and for a point . In this subsection, we restrict the domain of the distance functions to . By definition, the upper envelope of for all apexed triangles on coincides with in its projection on . We consider the sites one by one from to in order and compute the upper envelope of for all apexed triangles on .

While the upper envelope of for all apexed triangles is continuous, the upper envelope of of all apexed triangles with definers from up to (we simply say the upper envelope for sites from to ) might be discontinuous at some point on for . We let be the leftmost connected component of the upper envelope for sites from to along . By definition, is the upper envelope of the distance functions of all apexed triangles in . Note that for some apexed triangle . Thus the distance function of some apexed triangle might not appear on . Let be the list of the apexed triangles sorted in the order of their distance functions appearing on . If for any two consecutive apexed triangles and of , the bisector of and crosses the intersection of the bottom sides of and .

#### 3.2.2 Computing the Upper Envelope

Suppose that we have and for some index . We compute and from and as follows. We use two auxiliary lists and which are initially set to and . We update and until they finally become and , respectively. For simplicity, we use , and .

Let be the list of the apexed triangles of with definer sorted along with respect to their bottom sides. For any apexed triangle , we denote the list of the apexed triangles in overlapping with in their bottom sides by . Also, we denote the lists of the apexed triangles in lying left to and lying right to along with respect to their bottom sides by and , respectively. See Figure 1.

Let denote the the rightmost apexed triangle of along . With respect to , we partition into three disjoint sublists , and . We can compute these sublists in time.

##### Case 1 : Some apexed triangles in overlap with (i.e. ).

Let be the leftmost apexed triangle in along . We compare the distance functions and on . That is, we compare and for .

(1) If there is a point on that is equidistant from and , then appears on . Moreover, the distance functions of the apexed triangles in also appear on , but no distance function of the apexed triangles in appears on by Lemma 1. Thus we append the triangles in . We also update accordingly. Then, and are and , respectively.

(2) If for all points , then and its distance function do not appear on and , respectively, by Lemma 1. Thus we do nothing and scan the apexed triangles in , except , from left to right along until we find an apexed triangle such that there is a point on which is equidistant from and . Then we apply the procedure in (1) with instead of . If there is no such apexed triangle, then and are and , respectively.

(3) Otherwise, we have for all points . Then the distance function of does not appear on . Thus, we remove and its distance function from and , respectively. We consider the apexed triangles in from right to left along . For an apexed triangle , we do the following. Since is updated, we update to the rightmost element of along . We check whether for all points if overlaps with . If so, we remove from and update again. We do this until we find an apexed triangle such that this test fails. Then, there is a point on which is equidistant from and . After we reach such an apexed triangle , we apply the procedure in (1) with instead of .

##### Case 2 : No apexed triangle in overlaps with (i.e. ).

We cannot compare the distance function of any apexed triangle in with the distance function of directly, so we need a different method to handle this. There are two possible subcases: either or . Note that these are the only possible subcases since the union of the apexed triangles with the same definer is connected. For the former subcase, the upper envelope of sites from to is discontinuous at the right endpoint of the bottom side of along . Thus does not appear on for any apexed triangle . Thus and are and , respectively. For the latter subcase, at most one of and has a Voronoi cell in by Lemma 1. We can find a site ( or ) which does not have a Voronoi cell in in constant time once we maintain some geodesic paths. We describe this procedure at the end of this subsection.

If does not have a Voronoi cell in , then and are and , respectively. If does not have a Voronoi cell in , we remove all apexed triangles with definer from and their distance functions from . Since such apexed triangles lie at the end of consecutively, this removal process takes the time linear in the number of the apexed triangles. We repeat this until the rightmost element of and the rightmost element of overlap in their bottom sides along . When the two elements overlap, we apply the procedure of Case 1.

In total, the running time for computing is since each apexed triangle in is removed from at most once. Thus, we can compute is time in total.

##### Maintaining Geodesic Paths for Subcase of Case 2 : and .

We maintain and its geodesic distance during the whole procedure (for all cases), where is the site we consider and is the projection of the rightmost breakpoint of onto . That is, is the projection of the common endpoint of the two rightmost pieces of onto . Recall that changes from to . By definition, lies in the bottom side of the rightmost apexed triangle of . Thus we can evaluate in constant time. Note that the two points and the apexed triangle change during the procedure. Whenever they change, we update and its geodesic distance using the previous geodesic path. One of and does not have a Voronoi cell in in this subcase. But it is possible that neither nor has a Voronoi cell in . We can decide which site does not have a Voronoi cell in in constant time: if , then does not have a Voronoi cell. Otherwise, does not have a Voronoi cell.

We will show that the update of the geodesic path takes time in total for all transition edges. Let denote the region bounded by , , and . The sum of the complexities of for all transition edges is and they can be computed in time (Corollary 3.8 [2]). Moreover, is (Lemma 5.2 [2]). The total complexity of the shortest path trees rooted at and in is , and therefore we can compute them in time [7]. We compute them only one for each transition edge during the whole procedure.

The edges in , except the edge adjacent to , are also edges of the shortest path trees, and thus we can update them by traversing the shortest path trees in time linear in the amount of the changes on . Therefore, the following lemma implies that maintaining and its length takes time for each transition edge .

###### Lemma 8.

The amount of the changes on is during the whole procedure for .

###### Proof.

We claim that each edge of the shortest path trees is removed from at most times during the whole procedure for . Assume that we already have and we are to compute . There are three different cases: (1) lies to the left of ( is removed) along , (2) lies to the right of (a new apexed triangle is inserted to ) along , and (3) we consider a new site (that is, and .)

For the first and the second cases, it is possible that we remove more than one edge from . We prove the claim for the first case only. The claim for the second case can be proved analogously. See Figure 2(a). Let be an edge in which is not adjacent to and is not in with . Let and be the points on hit by the rays from towards and towards , respectively. The right endpoint of the bottom side of lies to the right of since contains . Moreover, lies in . Thus, lies in .

There are two possible subcases: is in , or in . There is at most one site in such that an apexed triangle with definer has its apex in by the construction of the set of apexed triangles in [2]. (In this case, the apex lies in .) When is deleted, all such apexed triangles are also deleted from . After is deleted, no apexed triangle with definer in and with apex in is inserted to again. Therefore, the number of deletions of due to the first subcase is only one. For the second subcase, notice that once is removed from , no apexed triangle with definer in is added to again. Thus, the number of deletions of due to the second subcase is also one.

For the third case, lies after from in counterclockwise order along . It occurs when we finish the procedure for handling . After we consider the site , we do not consider any site from to again. Consider an edge removed from due to this case. Let be the point on hit by the extension of . See Figure 2(b). If contains for some and some , we have . This means that once is removed due to the last case, does not appear on the geodesic path again in the remaining procedure. Thus, the number of deletions of each edge due to the last case is also one. ∎

Therefore, we can complete the first step in time and we have the following theorem.

###### Theorem 9.

The geodesic farthest-point Voronoi diagram of the vertices of a simple -gon restricted to the boundary of can be computed in time.

## 4 Decomposing the Polygon into Smaller Cells

Now we have of size . We add the points in (degree- vertices of ) to the vertex set of , and apply the algorithm to compute the apexed triangles with respect to the vertex set of again [2]. There is no transition edge because no additional vertex has a Voronoi cell and every degree-1 Voronoi vertex is a vertex of . Thus the bottom sides of all apexed triangles are interior disjoint. Moreover, we have the set of the apexed triangles sorted along with respect to their bottom sides.

###### Definition 10.

A simple polygon is called a -path-cell for some if it is geodesically convex and all its vertices are on among which at most are convex.

In this section, we subdivide into -path-cells recursively for some until each cell becomes a base cell. There are three types of base cells. The first type is a quadrilateral crossed by exactly one arc of through two opposite sides, which we call an arc-quadrilateral. The second type is a -path-cell. Note that a -path-cell is a pseudo-triangle. The third type is a region of whose boundary consists of one convex chain and one geodesic path (concave curve), which we call a lune-cell. Note that a convex polygon is a lune-cell whose concave chain is just a vertex of the polygon.

Let be the sequence such that and . Initially, itself is a -path-cell. Assume that the th iteration is completed. We show how to subdivide each -path-cell with into -path-cells and base cells in the th iteration in Section 4.1. A base cell is not subdivided further. While subdividing a cell into a number of smaller cells recursively, we compute the refined geodesic farthest-point Voronoi diagram restricted to the boundary of each smaller cell (of any kind) in time. In Section 5, we show how to compute the refined geodesic farthest-point Voronoi diagram restricted to a base cell in time once we have . Once we compute restricted to every base cell, we obtain .

### 4.1 Subdividing a -path-cell into Smaller Cells

In this subsection, we are to subdivide each -path-cell into -path-cells and base cells. If a -path-cell is a lune-cell or is at most three, the cell is a already base cell and we do not subdivide it further. Otherwise, we subdivide it using the algorithm described in this subsection.

The subdivision consists of three phases. In Phase 1, we subdivide each -path-cell into -path-cells by a curve connecting at most convex vertices of the -path-cell. In Phase 2, we subdivide each -path-cell further along the arcs of crossing the cell if there are such arcs. In Phase 3, we subdivide the cells that are created in Phase 2 and have vertices in into -path-cells and lune-cells.

#### 4.1.1 Phase 1. Subdivision by a Curve Connecting at Most Vertices

Let be a -path-cell computed in the th iteration. Recall that is a simple polygon which has at most convex vertices. Let be the largest integer satisfying that is less than the number of the convex vertices of . Then we have .

We choose vertices from the convex vertices of at a regular interval as follows. We choose an arbitrary convex vertex of and denote it by . Then we choose the th convex vertex of from in clockwise order and denote it by for all . We set . Then we construct the closed curve (or simply when is clear from context) consisting of the geodesic paths . See Figure 3(a). In other words, the closed curve is the boundary of the geodesic convex hull of . Note that does not cross itself. Moreover, is contained in since is geodesically convex.

We compute in time linear in the number of edges of using the algorithm in [15]. This algorithm takes source-destination pairs as input, where both sources and destinations are on the boundary of the polygon. It returns the geodesic path between the source and the destination for every input pair assuming that the shortest paths do not cross (but possibly overlap) one another. Computing the geodesic paths takes time in total, where is the complexity of the polygon. In our case, the pairs for are input source-destination pairs. Since the geodesic paths for all input pairs do not cross one another, can be computed in time. Then we compute in time using obtained from the th iteration. We will describe this procedure in Section 4.2.

The curve subdivides into -path-cells. To be specific, consists of at least connected components. Note that the closure of each connected component is a -path-cell. Moreover, the union of the closures of all connected components is exactly since is simple. These components define the subdivision of induced by .

#### 4.1.2 Phase 2. Subdivision along an Arc of

After subdividing into -path-cells by the curve , an arc of may cross for some . We say an arc of crosses a cell if intersects at least two distinct edges of . For example, in Figure 3(c), crosses while does not cross because crosses only one edge of . In Phase 2, for each arc crossing , we isolate the subarc . That is, we subdivide further into three subcells so that only one of them intersects . We call such a subcell an arc-quadrilateral. Moreover, for an arc-quadrilateral created by an arc crossing , we have .

###### Lemma 11.

For a geodesic convex polygon with convex vertices (), let be a simple closed curve connecting at most convex vertices of lying on such that every two consecutive vertices are connected by a geodesic path. Then, each arc of intersecting intersects at most three cells in the subdivision of induced by and at most two edges of .

###### Proof.

Consider an arc of intersecting . The arc is a part of either a side of some apexed triangle or the bisector of two sites. For the first case, the arc is a line segment. Thus intersects at most three cells in the subdivision of by and at most two edges of . For the second case, is part of a hyperbola. Let and be the two sites defining in . The combinatorial structure of the geodesic path from (or ) to any point in is the same. This means that is contained in the intersection of two apexed triangles and , one with definer and the other with definer . Observe that intersects at most twice and contains no vertex of in its interior. By construction, intersects at most two edges and of , and thus so does . For a cell in the subdivision of by , the arc intersects if and only if contains or on its boundary. Thus there exist at most three such cells in the subdivision of by and the lemma holds for the second case. See Figure 3(b-c). ∎

First, we find for every arc of crossing . If consists of at most two connected components (line segments) for every apexed triangle , we can do this by scanning all points in along . However, might consist of more than two connected components (line segments) for some apexed triangle . See Figure 4. Despite of this fact, we can compute all such arcs in time by the following lemma.

###### Lemma 12.

For every arc of crossing , we can find the part of contained in in time in total. Moreover, for each such arc , the pair of apexed triangles such that can be found in the same time.

###### Proof.

For each apexed triangle intersecting , we find all connected components of . Since we already have , this takes time for all apexed triangles in intersecting . There are at most two edges of that are intersected by due to Lemma 11. Let and be such edges, and we assume that contains the point in closest to without loss of generality. We insert all connected components of in the clockwise order along into a queue. Then, we consider the connected components of in the clockwise order along one by one.

To handle a connected component of , we do the following. Let be a point in the first element of the queue. If the line passing through and intersects , then we remove from the queue and check whether and are incident to the same refined cell . If so, we compute the part of the arc defined by and inside , and return the part of the arc and the pair . If and are not incident to the same refined cell, we remove from the queue since. We repeat this until the line passing through a point of the first element of the queue does not intersect . Then we handle the connected component of next to .

Every arc computed from this procedure is an arc of crossing . The remaining work is to show that we can find all arcs of crossing using this procedure. Consider an arc of crossing . There are two connected components