    # A Nearly Optimal Algorithm for the Geodesic Voronoi Diagram in a Simple Polygon

The geodesic Voronoi diagram of m point sites inside a simple polygon of n vertices is a subdivision of the polygon into m cells, one to each site, such that all points in a cell share the same nearest site under the geodesic distance. The best known lower bound for the construction time is Omega( n + m log m ), and a matching upper bound is a long-standing open question. The state-of-the-art construction algorithms achieve O( (n+m) log (n+m) ) and O( n+m log m log^2 n ) time, which are optimal for m=Omega(n) and m=O( n / log^3 n ), respectively. In this paper, we give a construction algorithm with O( n+m( log m+ log^2 n) ) time, and it is nearly optimal in the sense that if a single Voronoi vertex can be computed in O( log n ) time, then the construction time will become the optimal O( n+m log m ). In other words, we reduce the problem of constructing the diagram in the optimal time to the problem of computing a single Voronoi vertex in O( log n ) time.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The geodesic Voronoi diagram of point sites inside a simple polygon of vertices is a subdivision of the polygon into cells, one to each site, such that all points in a cell share the same nearest site where the distance between two points is the length of the shortest path between them inside the polygon. The common boundary between two cells is a Voronoi edge, and the endpoints of a Voronoi edge are Voronoi vertices. A cell can be augmented into subcells such that all points in a subcell share the same anchor, where the anchor of a point in the cell is the vertex of the shortest path from the point to the associated site that immediately precedes the point. An anchor is either a point site or a reflex polygon vertex. Figure 1(a) illustrates an augmented diagram.

The size of the (augmented) diagram is  . The best known construction time is   and  . They are optimal for and for , respectively, since the best known lower bound is . The existence of a matching upper bound is a long-standing open question by Mitchell .

Aronov  first proved fundamental properties: a bisector between two sites is a simple curve consisting of straight and hyperbolic arcs and ending on the polygon boundary; the diagram has vertices, of which are Voronoi vertices. Then, he developed a divide-and-conquer algorithm that recursively partitions the polygon into two roughly equal-size sub-polygons. Since each recursion level takes time to extend the diagrams between every pair of sub-polygons, the total time is .

Papadopoulou and Lee  combined the divide-and-conquer and plane-sweep paradigms to improve the construction time to . First, the polygon is triangulated and one resultant triangle is selected as the root such that the dual graph is a rooted binary tree and each pair of adjacent triangles have a parent-child relation; see Figure 1(b). For each triangle, its diagonal shared with its parent partitions the polygon into two sub-polygons: the “lower” one contains it, and the “upper” one contains its parent. Then, the triangles are swept by the post-order and pre-order traversals of the rooted tree to respectively build, inside each triangle, the two diagrams with respect to sites in its lower and upper sub-polygons. Finally, the two diagrams inside each triangle are merged into the final diagram.

Very recently, Oh and Ahn  generalized the notion of plane sweep to a simple polygon. To supplant the scan line, one point is fixed on the polygon boundary, and another point is moved from the fixed point along the polygon boundary counterclockwise, so that the shortest path between the two points will sweep the whole polygon. Moreover, Guibas and Hershberger’s data structure for shortest path queries [5, 7] is extended to compute a Voronoi vertex among three sites or between two sites in time. This technique enables handling an event in time, leading to a total time of .

Papadopoulou and Lee’s method  has two issues inducing the time-factor. First, while sweeping the polygon, an intermediate diagram is represented by a “wavefront” in which a “wavelet” is associated with a “subcell.” Although this representation enables computing the common vertex among three subcells in time, since it takes  time to update such a wavefront, the vertices lead to the factor. Second, when a wavefront enters a triangle from one diagonal and leaves from the other two diagonals, it will split into two. Since there are triangles, there are split events, and since a split event takes time, the factor arises again.

### 1.1 Our contribution

We devise a construction algorithm with time, which is slightly faster than Oh and Ahn’s method  and is optimal for . More importantly, our algorithm is, to some extent, nearly optimal since the factor solely comes from computing a single Voronoi vertex. If the computation time can be improved to , the total construction time will become , which equals the optimal since for and for . In other words, we reduce the problem of constructing the diagram in the optimal time by Mitchell  to the problem of computing a single Voronoi vertex in time.

At a high level, our algorithm is a new implementation of Papadopoulou and Lee’s concept  using a different data structure of a wavefront, symbolic maintenance of incomplete Voronoi edges, tailor-made wavefront operations, and appropriate amortized time analysis.

First, in our wavefront, each wavelet is directly associated with a cell rather than a subcell. This representation makes use of Oh and Ahn’s  -time technique of computing a Voronoi vertex. Each wavelet also stores the anchors of incomplete subcells in its associated cell in order to enable locating a point in a subcell along the wavefront.

Second, if each change of a wavefront needs to be updated immediately, a priority queue for events would be necessary, and since the diagram has vertices, an time-factor would be inevitable. To overcome this issue, we maintain incomplete Voronoi edges symbolically, and update them only when necessary. For example, during a binary search along a wavefront, each incomplete Voronoi edge of a tested wavelet will be updated.

Third, to avoid split operations, we design two tailor-made operations. If a wavefront will separate into two but one part will not be involved in the follow-up “sweep”, then, instead of using a binary search, we “divide” the wavefront in a seemingly brute-force way in which we traverse the wavefront from the uninvolved part until the “division” point, remove all visited subcells, and build another wavefront from those subcells. If a wavefront propagates into a sub-polygon that contains no point site, then we adopt a two-phase process to build the diagram inside the sub-polygon instead of splitting a wavefront many times.

Finally, when deleting or inserting a subcell (anchor), its position in a wavelet is known. Since re-balancing a red-black tree (RB-tree) after an insertion or a deletion takes amortized time [8, 12, 13], by augmenting each tree node with pointers to its predecessor and successor, an insertion or a deletion with a known position takes amortized  time.

This paper is organized as follows. Section 2 formulates the geodesic Voronoi diagram, defines a rooted partition tree, and introduces Papadopoulou and Lee’s two subdivisions ; Section 3 summarizes our algorithm; Section 4 designs the data structure of a wavefront; Section 5 presents wavefront operations; Section 6 implements the algorithm with those operations.

## 2 Preliminary

### 2.1 Geodesic Voronoi diagrams

Let be a simple polygon of vertices, let denote the boundary of , and let be a set of point sites inside . For any two points in , the geodesic distance between them, denoted by , is the length of the shortest path between them that fully lies in , the anchor of with respect to is the last vertex on the shortest path from to before , and the shortest path map (SPM) from in is a subdivision of such that all points in a region share the same anchor with respect to . Each edge in the SPM from is a line segment from a reflex polygon vertex of to along the direction from the anchor of (with respect to ) to , and this line segment is called the SPM edge of (from ).

The geodesic Voronoi diagram of in , denoted by , partitions into cells, one to each site, such that all points in a cell share the same nearest site in under the geodesic distance. The cell of a site can be augmented by partitioning the cell with the SPM from into subcells such that all points in a subcell share the same anchor with respect to . The augmented version of is denoted by . With a slight abuse of terminology, a cell in indicates a cell in together with its subcells in . Then, each cell is associated with a site, and each subcell is associated with an achor. As shown in Figure 1(a), is the anchor of the shaded subcell (in ’s cell), and the last vertex on the shortest path from to any point in the shaded subcell before is . Figure 1: (a) Augmented geodesic Voronoi diagram Vor∗P(S). (b) Rooted partition tree T.

A Voronoi edge is the common boundary between two adjacent cells, and the endpoints of a Voronoi edge are called Voronoi vertices. A Voronoi edge is a part of the geodesic bisector between the two associated sites, and consists of straight and hyperbolic arcs. Endpoints of these arcs except Voronoi vertices are called breakpoints, and a breakpoint is incident to an SPM edge in the SPM from one of the two associated sites, indicating a change of the corresponding anchor. There are Voronoi vertices and breakpoints .

In our algorithm, each anchor refers to either a reflex polygon vertex of or a point site in ; we store its associated site , its geodesic distance from , and its anchor with respect to . The weighted distance from to a point is .

Throughout the paper, we make a general position assumption that no polygon vertex is equidistant from two sites in and no point inside is equidistant from four sites in . The former avoids nontrivial overlapping among cells , and the latter ensures that the degree of each Voronoi vertex with respect to Voronoi edges is either 1 (on ) or 3 (among three cells).

The boundary of a cell except Voronoi edges are polygonal chains on . For convenience, these polygonal chains are referred to as polygonal edges of the cell, the incidence of an SPM edge onto a polygonal edge is also a breakpoint, and the polygonal edges including their polygon vertices and breakpoints also count for the size of the cell.

([10, Lemma 5 and 14]) It takes time to compute the degree-1 or degree-3 Voronoi vertex between two sites or among three sites after -time preprocessing.

### 2.2 A rooted partition tree

Following Papadopoulou and Lee , a rooted partition tree for and is built as follows: First, is triangulated using Chazelle’s algorithm  in time, and all sites in are located in the resulting triangles by Edelsbrunner et al’s approach  in time. The dual graph of the triangulation is a tree in which each node corresponds to a triangle and an edge connects two nodes if and only if their corresponding triangles share a diagonal. Then, an arbitrary triangle whose two diagonals are polygon sides, i.e., a node with degree 1, is selected as the root, so that there is a parent-child relation between each pair of adjacent triangles. Figure 1(b) illustrates a rooted partition tree .

For a diagonal , let and be the two triangles adjacent to such that is the parent of , and call the root diagonal of ; also see Figure 1(b). partitions into two sub-polygons: contains and contains . and are said to be “below” and “above” , respectively. Assume that ; let and , which indicate the two respective subsets of sites below and above .

In this paper, we adopt the following convention: for each triangle , let be its root diagonal, let be its parent triangle, let be the other two diagonals of , and let be the other two diagonals of . Assume that is the diagonal between and its parent, and that , , and are incident to , and , , and are incident to . Let be the vertex shared by and , and let be the vertex shared by and . Denote the set of sites in as . Figure 2(a) shows an illustration. Figure 2: (a) △ and △′. (b) SD∩△. (c) SD′∩△. (Borders on d are indicated by arrows.)

### 2.3 Subdivisions

Papadopoulou and Lee  introduced two subdivisions, and , of , which can be merged into . For each triangle with a root diagonal , and respectively contain and ; also contains . Since and (resp.  and ) may differ, a border forms along in (resp. ) to “remedy” the conflicting proximity information. Figure 2(b)–(c) illustrate and .

The incidence of a Voronoi edge or an SPM edge in onto a border in or is called a border vertex, and the border vertices partition a border into border edges. Both and have border vertices , of which are induced by Voronoi edges. Hereafter, a diagram vertex means a Voronoi vertex, a breakpoint, a border vertex, or a polygon vertex.

## 3 Overview of the algorithm

We compute in the following three steps:

1. Build the rooted partition tree for and in time. (Section 2.2)

2. Construct and in time by sweeping the polygon using the post-order and pre-order traversals of , respectively. (Section 6.1 and Section 6.2)

3. Merge and into in time using Papadopoulou and Lee’s method [11, Section 7].

By the above-mentioned running times, we conclude the total running time as follows. can be constructed in time.

## 4 Wavefront structure

A wavefront represents the “incomplete” boundary of “incomplete” Voronoi cells during the execution of our algorithm, and wavefronts will “sweep” the simple polygon triangle by triangle to construct and . To avoid excessive updates, each “incomplete” Voronoi edge, which is a part of a Voronoi edge and will be completed during the sweep, is maintained symbolically, preventing an extra time-factor. During the sweep, candidates for Voronoi vertices in and called potential vertices will be generated in the unswept part of .

### 4.1 Formal definition and data structure

Let be a diagonal or a pair of diagonals sharing a common polygon vertex, and let be a subset of lying on the same side of . A wavefront represents the sequence of Voronoi cells in appearing along , and each appearance of a cell induces a wavelet in . The unswept area of is the part of on the opposite side of from . Since in the unswept area has not yet been constructed, Voronoi and polygonal edges incident to are called incomplete. Each wavelet is bounded by two incomplete Voronoi or polygonal edges along , and its incomplete boundary comprises its two incomplete edges and the portion of between them. When the context is clear, a wavelet may indicate its associated cell.

is stored in an RB-tree in which each node refers to one wavelet and the ordering of nodes follows their appearances along . The RB-tree is augmented such that each node has pointers to its predecessor and successor, and the root has pointers to the first and last nodes, enabling efficiently traversing wavelets along and accessing the two ending wavelets.

The subcells of a wavelet are the subcells in its associated cell incident to its incomplete boundary. The list of their anchors is also maintained by an augmented RB-tree in which their ordering follows their appearances along the incomplete boundary. Due to the visibility of a subcell, each subcell appears exactly once along the incomplete boundary. Since the rebalancing after an insertion or a deletion takes amortized time [12, 13, 8], inserting or deleting an anchor at a known position in an RB-tree, i.e., without searching, takes amortized time.

### 4.2 Incomplete Voronoi and polygonal edges

When a wavefront moves into its unswept area, incomplete Voronoi edges will extend, generating new breakpoints. If each breakpoint needs to be created immediately, all candidates for breakpoints should be maintained in a priority queue, leading to an running time due to breakpoints. To avoid these excessive updates, we maintain each incomplete Voronoi edge symbolically, and update it only when necessary. For example, when searching along the wavefront or merging two wavefronts, the incomplete Voronoi edges of each involved wavelet (cell) will be updated until the diagonal or the pair of diagonals.

Since a breakpoint indicates the change of a corresponding anchor, for a Voronoi edge, if the anchors of its incident subcells on “each side” are stored in a sequence, the Voronoi edge can be computed in time proportional to the number of breakpoints by scanning the two sequences [10, Section 4]. Following this concept, for each incomplete Voronoi edge, we store its fixed Voronoi vertex (in the swept area), its last created breakpoint, and its last used anchors on its two sides, so that we can update an incomplete Voronoi edge by scanning the two lists of anchors from the last used anchors. When creating a breakpoint, we also build a corresponding SPM edge, and then remove the corresponding anchor from the anchor list.

Each polygonal edge is also maintained symbolically in a similar way; in particular, it will also be updated when a polygon vertex is inserted as an anchor. Meanwhile, the SPM edges incident to a polygonal edge will also be created using its corresponding anchor list.

### 4.3 Potential vertices

We process incomplete Voronoi edges to generate candidates for Voronoi vertices called potential vertices. For each incomplete Voronoi edge, since its two associated sites lie in the swept area, one endpoint of the corresponding bisector lies in the unswept area and is a degree-1 potential vertex. For each two adjacent Voronoi edges along the wavefront, their respective bisectors may intersect in the unswept area, and the intersection is a degree-3 potential vertex. By Lemma 1, a potential vertex can be computed in time.

Potential vertices are stored in their located triangles; each diagonal of a triangle is equipped with a priority queue to store the potential vertices in the triangle associated with sites on the opposite side of the diagonal, where the key is the distance to the diagonal.

## 5 Wavefront operations

We summarize the eight wavefront operations, where , and are the numbers of involved (i.e., processed and newly created) potential vertices, visited anchors, and created diagram vertices, respectively:

Initiate:

Compute and initialize in time.

Extend:

Extend one wavefront into a triangle from one diagonal to the opposite pair of diagonals to build the diagram inside the triangle in plus amortized time.

Merge:

Merge two wavefronts sharing the same diagonal or the same pair of diagonals together with merging the two corresponding diagrams in plus amortized time where is the underlying triangle.

Join:

Join two wavefronts sharing the same diagonal with building the border on the diagonal but without merging the two corresponding diagrams inside the underlying triangle in plus amortized time.

Split:

Split a wavefront using a binary search in time.

Divide:

Divide a wavefront by traversing one diagonal in amortized time.

Insert:

Insert into to be in time.

Propagate:

Propagate into , provided that , to build in time.

Merge and Join operations differ in that the former also merges the two corresponding Voronoi diagrams in the underlying triangle, and the latter does not.

Readers could directly read the algorithm in Section 6 without knowing the detailed implementation for wavefront operations. For the operation times, we have three remarks. During Extend and Propagate operations and at the end of the other operations, potential vertices will be generated according to new incomplete Voronoi edges. It will be clear in Section 6 that the total number of potential vertices in the construction of and is , so a priority queue takes time for an insertion or an extraction.

As stated in Section 4.2, we maintain incomplete Voronoi/polygonal edges symbolically. For the sake of simplicity, we charge the time to update an incomplete edge to the created breakpoints, and assign the corresponding costs to their located triangles.

Since a wavelet (resp. anchor) to remove from a wavefront (resp. wavelet) must be inserted beforehand, we charge its removal cost at its insertion time. For a wavelet, the cost is , and for an anchor, since the position is always known, the cost is amortized . Similarly, we charge the cost to delete a diagram vertex at its creation time.

### 5.1 Initiate operation

An Initiate operation processes a triangle to compute and initiate . Since no polygon vertex lies in a triangle, each resulting Voronoi cell has only one subcell. For the sake of simplicity, we assume that no diagonal of is a polygon side; the other cases are similar.

An Initiate operation consists of two simple steps: First, is computed by constructing the Voronoi diagram of without considering and trimming the diagram by . Second, by tracing along the pair of diagonals, is built as the aforementioned data structure (Section 4.1). During the tracing, if (resp. ) is a reflex polygon vertex, then it will be inserted into the anchor list of its located wavelet (cells) in . Note that if is a reflex polygon vertex, then will be inserted as an anchor after has been divided or split at .

The operation time is . The first step takes time to construct the Voronoi diagram using a standard divide-and-conquer or a plane-sweep algorithm [2, Section 3], and takes time to trim the diagram. The tracing in the second step simply takes time by first locating in , and then repeatedly traversing the boundary of the located Voronoi cell to find the next cell along until reaching . The data structure can also be built from scratch in time. Since a cell in has only one subcell, whose anchor is exactly the associated site, the position to insert or into a wavelet as an anchor can be computed in time, and although inserting an anchor at a known position takes amortized time, since the worst-case time to insert an anchor is , the time to insert and is dominated by the time-factor . Finally, since there are incomplete Voronoi edges, there are potential vertices, leading to the time factor . (Computing a potential vertex takes time (Lemma 1), locating it in a triangle of takes time, and inserting it in the priority queue takes time.)

### 5.2 Extend operation

An Extend operation extends a wavefront from one diagonal of a triangle to the opposite pair of diagonals to construct and , where , , and , and lies on the opposite side of from .

This operation is equivalent to sweeping the triangle with a scan line parallel to and processing each hit potential vertex. The next hit potential vertex is provided from the priority queue associated with (defined in Section 4.3), and will be processed in three phases: First, its validity is verified: for a degree-1 potential vertex, its incomplete Voronoi edge should be alive in the wavefront, and for a degree-3 one, its two incomplete Voronoi edges should be still adjacent in the wavefront. Second, the one or two incomplete Voronoi edges are updated up to the potential vertex. Since a degree-1 potential vertex is incident to a polygonal edge, the polygonal edge is also updated.

Finally, when a potential vertex becomes a Voronoi vertex, a wavelet will be removed from the wavefront. For a degree-1 potential vertex, since the wavelet lies at one end of the wavefront, a polygonal edge is added to its adjacent wavelet; for a degree-3 potential vertex, since the removal makes two wavelets adjacent, a new incomplete Voronoi edge is created, and one degree-1 and two degree-3 potential vertices are computed accordingly.

After the extension, if is a reflex polygon vertex, the wavefront is further processed by three cases. If neither nor is a polygon side, will later be inserted as an anchor while dividing or splitting at . If both and are polygon sides, all the subcells will be completed along and the wavefront will disappear. If only (resp. ) is a polygon side, all the subcells along (resp. ) excluding the one containing will be completed, and will be inserted into the corresponding wavelet as the last or first anchor.

The operation time is plus amortized , where is the number of involved (i.e., processed and newly created) potential vertices and is the number of created diagram vertices. First, extracting a potential vertex from the priority queue takes time, and verifying its validity takes time. Second, updating incomplete Voronoi and polygonal edges takes time linear in the number of created breakpoints (Section 4.2), i.e., time in total. Third, since computing a potential vertex takes time, locating it takes time, and inserting it into a priority queue takes time, creating a new potential vertex takes time. Finally, completing subcells takes time, and inserting takes amortized time since its position in the anchor list is known.

### 5.3 Merge operation

A Merge operation merges and into together with and into where is the underlying triangle. In our algorithm, either , , and or , , and ; in both cases, . The border will form on , i.e., for the former case and for the latter case. Although a wavefront only stores incomplete Voronoi cells, its associated diagram can still be accessed through the Voronoi edges of the stored cells. After the merge, new incomplete Voronoi edges will form, and their potential vertices will be created.

The Merge operation consists of two phases: (1) merge and into and (2) merge and into .

The first phase is to construct so-called merge curves. A merge curve is a connected component consisting of border edges along and Voronoi edges in associated with one site in and one site in ; the ordering of mergve curves is the ordering of their appearances along . This phase is almost identical to the merge process by Papadopoulou and Lee [11, Section 5], but since our data structure for a wavefront is different from theirs, a binary search along a wavefront to find a starting endpoint for a merge curve requires a different implementation. We first state this difference here, and will describe the details of tracing a merge curve in Section 5.3.1.

Assume to be oriented from to for and from to for . By [11, Lemma 4–6], a merge curve called initial starts from for the latter case (), but all other merge curves have both endpoints on , and those endpoints are associated with one site in . Let be the set of sites in that have a wavelet in . If , a site in can have two wavelets in , and with an abuse of terminology, such a site is imagined to have two copies, each to one wavelet. Since , finding a starting endpoint for each merge curve except the initial one is to test sites in following the ordering of their wavelets in along . After finding a starting endpoint, the corresponding merge curve will be traced; when the tracing reaches again, a stopping endpoint forms, and the first site in lying after the site inducing the stopping endpoint will be tested.

Let be the next starting endpoint, which is unknown, and let be the next site in to test. A two-level binary search on determines if induces , and if so, further determines the site that induces with as well as the corresponding anchor.

The first-level binary search executes on the RB-tree for the wavelets in , and each step determines for a site if its cell lies before or after ’s cell along or if . Let and denote the two intersection points between and the two incomplete edges of (in ), where lies before along , and let and be the two points defined similarly for (in ). Since lies in , the distance between and any point in can be computed in time. The two incomplete edges of will be updated until and , so that the distance from to (resp. to ) can be computed from the corresponding anchor. For example, if is the anchor of the subcell that contains , .

The determination considers four cases. (Assume and induce the “starting” endpoint.)

• lies before (resp.  lies after ): ’s cell lies after (resp. before) ’s cell.

• lies before and lies between and (resp.  lies after and lies between and ): if is closer to than to (resp.  is closer to than to ), then ’s cell lies after (resp. before) ’s cell; otherwise, is .

• Both and lie between and : is .

• Both and lie between and : let be the projection point of onto .

• If lies before , then ’s cell lies before ’s cell.

• If lies after : if is closer to than to , ’s cell lies after ’s cell; if both and are closer to than to , ’s cell lies before ’s cell; otherwise, .

• If lies between and : if is closer to than to , ’s cell lies before ’s cell; otherwise, .

If the first-level search does not find , then does not induce the next starting endpoint .

The second-level binary search executes on the RB-tree for ’s anchor list to either determine the next starting endpoint and ’s corresponding anchor or report that does not induce . Let be the current anchor of to test, and let be the projection point of onto . ’s “interval” on can be decided by checking ’s two neighboring anchors. If ’s interval lies after , lies before ’s interval; otherwise, if both endpoints of ’s interval are closer to (resp. farther from) than to (resp. from) , lies after (resp. before) ’s interval, and if one endpoint is closer to than to but the other is not, then lies in ’s interval and can be computed in time since . If the second-level binary search does not find such an interval, does not induce the next starting endpoint .

The second phase (i.e., merging and into ) splits and at the endpoints of merge curves, and concatenates active parts at these endpoints where a part is called active if it contributes to . In fact, the active parts along alternately come from and . At each merging endpoint, the two cells become adjacent, generating a new incomplete Voronoi edge. Potential vertices of these incomplete Voronoi edges will be computed and inserted into the corresponding priority queues. For each ending polygon vertex of , if it is reflex but has not yet been an anchor of , it will be inserted into its located wavelet as the first or the last anchor.

The total operation time is plus amortized , where is the number of created potential vertices and is the number of created diagram vertices while merging the two diagrams. First, since each two-level binary search takes time, finding starting points takes time. Second, by Section 5.3.1 or by [11, Section 5], tracing a merge curve takes time linear in the number of deleted and created vertices, but the time to delete vertices has been charged at their creation, implying that tracing all the merge curves takes time.

Third, an incomplete Voronoi edge generates at least one potential vertex, so the number of new incomplete Voronoi edges is . Since an endpoint of a merge curve corresponds to a new incomplete Voronoi edge, there are split and concatenation operations, and since each operation takes time, it takes time to merge the two wavefronts. By the same analysis in Section 5.2, creating potential vertices takes time. Finally, inserting an ending polygon vertex of as an anchor takes amortized time.

#### 5.3.1 Tracing merge Curves in Merge Operation

We make some further definitions. The bisector between and is the collection of points with the same minimum geodesic distance from both and , namely . In our algorithm, each anchor refers to either a reflex polygon vertex of or a point site in ; we store its associated site , its geodesic distance from , and its anchor with respect to . The weighted distance from a point to is , and the weighted bisector between two anchors and is .

The process to trace a merge curve from the starting endpoint is identical to that presented by Papadopolou and Lee. Note that when a wavelet (incomplete cell) is visited, its incomplete edges will be updated to , enabling the tracing between (incomplete) subcells. The merge curve begins with the weighted bisector between the two anchors that induce the starting endpoint. When the merge curve changes the underlying subcell in one of the two diagrams, it continues with the weighted bisector between the new anchor and the other original anchor. When the merge curve hits , the tracing continues along in the direction along which the next point is closer to the current site in than to the current site in ; until reaching an equidistant point, it turns into the interior of again following the corresponding weighted bisector. The merge curve may visit several times. Finally, it reaches at the stopping endpoint. Border vertices form on when the merge curve enters or leaves and when the merge curve changes the underlying subcell in one side of .

Since tracing a merge curve takes time proportional to the number of deleted and created vertices but the time to delete vertices has been charged at their creation, tracing all the merge curves takes time in total, where is the number of newly created diagram vertices.

### 5.4 Join operation

A Join operation joins and into with building the border on but without merging and into . In the construction of (Section 6.2), our algorithm will actually join and into , and join and into . To interpret a Join operation, for the former case, one would replace , , , and with , , , and , respectively, and for the latter case, one would replace , , , and with , , , and , respectively.

The reason for not merging and into is that contributes nothing to . The border on while merging and into is the collection of points in that are closer to than to , namely .

At a high level, a Join operation first builds the border on , and then joins the two wavefronts according to the border. This operation is conceptually identical to the process in [11, Section 6] except that our data structure of a wavefront is different.

With a slight abuse of terminology, let denote the collection of connected components of border edges on in .111In [11, Section 6], the authors define as the Voronoi edges in that are associated with one site in and one site in , but only compute , which is exactly . These components in are called joint segments and are ordered from to . By [11, Lemma 14–15], satisfies the following conditions:

• At most one joint segment starts at , and such a joint segment is called initial.

• Except the initial one, both endpoints of each joint segment are associated with a site in . The first endpoint is called starting, and the second endpoint is called stopping.

The second condition results from the fact that the Voronoi cells associated with and do not “cross” each other along , and this desirable condition enables locating a starting endpoint of a non-initial joint segment by searching with sites in .

Assume that is oriented from to . Joint segments will be constructed one by one from to . For each joint segment (except the initial one), its starting endpoint is first located on by searching and then it is traced from its starting endpoint along until reaching its stopping endpoint. We remark that the initial joint segment will be directly traced from .

Let be the set of sites in that have a wavelet in . By the second condition, finding the starting endpoint of a joint segment is to “locate” sites of in since each endpoint of a joint segment (except the initial one) is associated with a site in . Therefore, sites in will be tested following the ordering of their wavelets in from to ; after tracing a joint segment, the test will continue on the next site in , which lies after the site in associated the last traced stopping endpoint. Note that each site in has exactly one wavelet in .

Let be the next starting endpoint, which is unknown, and let be the site in to test. A two-level binary search on , which is identical to the one in Section 5.3, can determine if induces , and if so, determine the site that induces with together with the corresponding anchor. The total time to find starting points is trivially .

The process to trace a joint segment from the starting endpoint is simpler than tracing a merge curve in Section 5.3. The process walks along from the starting endpoint. Every time when one of the two corresponding anchors in and changes, a border vertex is created accordingly. The process ends when reaching a point equidistant from the site of its located wavelet in and the site of its located in ; this point is a border vertex and also the second endpoint of the joint segment. Note that the change of an anchor can be determined by the anchor lists, and the distance from a site to a point can be determined through the corresponding anchor. The time to trace all the joint segments is proportional to the number of created border vertices, i.e., time in total, where is the number of newly created diagram vertices.

After constructing joint segments, we split each of and at the endpoints of these joint segments, and concatenate the active parts at these endpoints into , where “active” means having a wavelet in . For each ending polygon vertex of , i.e., or , if it is reflex but has not been an anchor of , it will be inserted into its located wavelet as the first or the last anchor.

The time to update the wavefronts is plus amortized , where is the number of created potential vertices. First, an incomplete Voronoi edge generates at least one potential vertex, so the number of new incomplete Voronoi edges is . Since an endpoint of a joint segment corresponds to a new incomplete Voronoi edge, there are split and concatenation operations, and since each operation takes time, it takes time to join the two wavefronts. Moreover, by the same analysis in Section 5.2, creating potential vertices takes time. Finally, inserting or as an anchor takes amortized time.

To conclude, a Joint operation takes plus amortized time.

### 5.5 Split operation

A Split operation splits a wavefront associated with a pair of diagonals at the common polygon vertex using a binary search. First, the operation locates the wavelet that contains the common polygon vertex and the corresponding anchor. Second, the operation splits the wavelet, i.e., the corresponding list of anchors, at the located anchor, and duplicates the located anchor since it appears in both resulting wavelets. Third, the operation splits the wavefront between the two resulting wavelets. Finally, if the common polygon vertex is reflex, the operation inserts it into both of the duplicate wavelets as an anchor.

The operation time is . Since a wavefront has wavelets and a wavelet has anchors, both locating the common polygon vertex and spliting the wavelet and wavefront take time. Although inserting the common polygon vertex takes amortized time, since the worst-case time to insert an anchor is , the time to insert the common polygon vertex is dominated by the time-factor . Since there is no new site, there is no new incomplete Voronoi edge and no new potential vertex.

### 5.6 Divide operation

A Divide operation divides a wavefront associated with a pair of diagonals at the common polygon vertex by traversing one diagonal instead of using a binary search. Although a Divide operation seems a brute-force way compared to a Split operation, since a Split operation takes time and there are events to separate a wavefront, if only Split operations are adopted, the total construction time would be .

First, the wavefront is traversed from the end of the selected diagonal subcell by subcell, i.e., anchor by anchor, until reaching the common vertex. Then, the wavefront is separated at the common polygon vertex by removing all the visited anchors except the last one, duplicating the last one, and building a new wavefront for these “removed” anchors and the duplicate anchor from scratch. Finally, if the common polygon vertex is reflex, it is inserted into its located wavelets in both resultant wavefronts as an anchor without a binary search (since it is the first or last anchor of its located wavelets).

The total operation time is amortized , where is the number of visited anchors. Since each wavelet (resp. anchor) records its two neighboring wavelets (resp. anchors) in the augmented RB-tree, the time to locate the common polygon vertex is . Recall that a cell must have one subcell, and the time to remove a wavelet or an anchor has been charged when it was inserted. Finally, building the new wavefront from scratch takes amortized time, and inserting the common vertex takes amortized time.

### 5.7 Insert operation

For a triangle , an Insert operation inserts into to form . An Insert operation executes in the construction of , but its outcome will be used to construct . Let denote the intermediate site set during the Insert operation, so that at the beginning.

For each site , let be its vertical projection point onto , and conduct a binary search on the wavefront to locate the anchor whose subcell contains . During the binary search, the two incomplete Voronoi/polygonal edges of each visited wavelet will be updated up to . Recall that the weighted distance between a point and an anchor associated with a site is .

If is closer to than to the located anchor under the weighted distance, traverse anchors in from along each direction of until either touching a polygon vertex of or reaching a point equidistant from and the current visited anchor under the weighted distance. If all points in the “interval” of a visited anchor on are fully closer to than to the anchor, remove the anchor from its associated wavelet, and if all the anchors of a wavelet have been removed, remove the wavelet from . After the traversal, insert into , namely create and insert its wavelet into . Since and lie on the same side of <