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.


Introduction
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 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 Θ(n + m) [1]. The best known construction time is O (n + m) log(n + m) [11] and O(n + m log m log 2 n) [10]. They are optimal for m = Ω(n) and for m = O( n log 3 n ), respectively, since the best known lower bound is Ω(n + m log m). The existence of a matching upper bound is a long-standing open question by Mitchell [9].
Aronov [1] first proved fundamental properties: a bisector between two sites is a simple curve consisting of Θ(n) straight and hyperbolic arcs and ending on the polygon boundary; the diagram has Θ(n+m) vertices, Θ(m) of which are Voronoi vertices. Then, he developed a divide-and-conquer algorithm that recursively partitions the polygon into two roughly equalsize sub-polygons. Since each recursion level takes O (n + m) log(n + m) time to extend the diagrams between every pair of sub-polygons, the total time is O (n + m) log(n + m) log n .
Papadopoulou and Lee [11] combined the divide-and-conquer and plane-sweep paradigms to improve the construction time to O (n+m) log(n+m) . First, the polygon is triangulated and one resultant triangle is selected as the root such that the dual graph is a rooted binary

Our contribution
We devise a construction algorithm with O n + m(log m + log 2 n) time, which is slightly faster than Oh and Ahn's method [10] and is optimal for m = O( n log 2 n ). More importantly, our algorithm is, to some extent, nearly optimal since the log 2 n factor solely comes from computing a single Voronoi vertex. If the computation time can be improved to O(log n), the total construction time will become O n + m(log m + log n) , which equals the optimal O(n+m log m) since m log n = O(n) for m = O( n log n ) and log n = O(log m) for m = Ω( n log n ). In other words, we reduce the problem of constructing the diagram in the optimal time by Mitchell [9] to the problem of computing a single Voronoi vertex in O(log n) time.
At a high level, our algorithm is a new implementation of Papadopoulou and Lee's concept [11] 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 [10] O(log 2 n)-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 Θ(m+n) vertices, an (m+n) log(m+ n) 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 Ω(n) 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 propa- gates 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 O(1) 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 O(1) 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 [11]; 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.

Geodesic Voronoi diagrams
Let P be a simple polygon of n vertices, let ∂P denote the boundary of P , and let S be a set of m point sites inside P . For any two points p, q in P , the geodesic distance between them, denoted by d(p, q), is the length of the shortest path between them that fully lies in P , the anchor of q with respect to p is the last vertex on the shortest path from p to q before q, and the shortest path map (SPM) from p in P is a subdivision of P such that all points in a region share the same anchor with respect to p. Each edge in the SPM from p is a line segment from a reflex polygon vertex v of P to ∂P along the direction from the anchor of v (with respect to p) to v, and this line segment is called the SPM edge of v (from p).
The geodesic Voronoi diagram of S in P , denoted by Vor P (S), partitions P into m cells, one to each site, such that all points in a cell share the same nearest site in S under the geodesic distance. The cell of a site s can be augmented by partitioning the cell with the SPM from s into subcells such that all points in a subcell share the same anchor with respect to s. The augmented version of Vor P (S) is denoted by Vor * P (S). With a slight abuse of terminology, a cell in Vor * P (S) indicates a cell in Vor P (S) together with its subcells in Vor * P (S). Then, each cell is associated with a site, and each subcell is associated with an achor. As shown in Figure 1(a), v is the anchor of the shaded subcell (in s 1 's cell), and the last vertex on the shortest path from s 1 to any point x in the shaded subcell before x is v.
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 Θ(m) Voronoi vertices and Θ(n) breakpoints [1].
In our algorithm, each anchor u refers to either a reflex polygon vertex of P or a point site in S; we store its associated site s, its geodesic distance from s, and its anchor with respect to s. The weighted distance from u to a point x is d(s, u) + |ux|.
Throughout the paper, we make a general position assumption that no polygon vertex is equidistant from two sites in S and no point inside P is equidistant from four sites in S. The former avoids nontrivial overlapping among cells [1], and the latter ensures that the degree of each Voronoi vertex with respect to Voronoi edges is either 1 (on ∂P ) or 3 (among three cells).
The boundary of a cell except Voronoi edges are polygonal chains on ∂P . 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.

A rooted partition tree
Following Papadopoulou and Lee [11], a rooted partition tree T for P and S is built as follows: First, P is triangulated using Chazelle's algorithm [3] in O(n) time, and all sites in S are located in the resulting triangles by Edelsbrunner et al's approach [4] in O(n+m log n) 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 T .
For a diagonal d, let (d) and (d) be the two triangles adjacent to d such that (d) is the parent of (d), and call d the root diagonal of (d); also see Figure 1(b). d partitions P into two sub-polygons: P (d) contains (d) and P (d) contains (d). P (d) and P (d) are said to be "below" and "above" d, respectively. Assume that d ⊆ P (d); let S(d) = S ∩ P (d) and S (d) = S \ S(d), which indicate the two respective subsets of sites below and above d.
In this paper, we adopt the following convention: for each triangle , let d = v 1 v 2 be its root diagonal, let be its parent triangle, let d 1 , d 2 be the other two diagonals of , and let d 3 , d 4 be the other two diagonals of . Assume that d 4 is the diagonal between and its parent, and that

Subdivisions
Papadopoulou and Lee [11] introduced two subdivisions, SD and SD , of P , which can be merged into Vor * P (S). For each triangle with a root diagonal d, SD and SD respectively contain Vor * P S(d) ∩ and Vor * P S (d) ∩ ; SD also contains Vor * P (S) ∩ . Since S(d) and S(d 4 ) (resp. S (d) and S (d 4 )) may differ, a border forms along d in SD (resp. SD ) to "remedy" the conflicting proximity information. Figure 2(b)-(c) illustrate SD and SD . The incidence of a Voronoi edge or an SPM edge in Vor * P (S) onto a border in SD or SD is called a border vertex, and the border vertices partition a border into border edges. Both SD and SD have O(n + m) border vertices [11], O(m) of which are induced by Voronoi edges. Hereafter, a diagram vertex means a Voronoi vertex, a breakpoint, a border vertex, or a polygon vertex.

Overview of the algorithm
We compute Vor * P (S) in the following three steps: 1. Build the rooted partition tree T for P and S in O(n + m log n) time. (Section 2.2) 2. Construct SD and SD in O n + m(log m + log 2 n) time by sweeping the polygon using the post-order and pre-order traversals of T , respectively. (Section 6.1 and Section 6.2) 3. Merge SD and SD into Vor * P (S) in O(n + m) time using Papadopoulou and Lee's method [11,Section 7].
By the above-mentioned running times, we conclude the total running time as follows.

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 P triangle by triangle to construct SD and SD . 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 log n time-factor. During the sweep, candidates for Voronoi vertices in SD and SD called potential vertices will be generated in the unswept part of P .

Formal definition and data structure
Let η be a diagonal or a pair of diagonals sharing a common polygon vertex, and let S be a subset of S lying on the same side of η. A wavefront W η (S ) represents the sequence of Voronoi cells in Vor * P (S ) appearing along η, and each appearance of a cell induces a wavelet in W η (S ). The unswept area of W η (S ) is the part of P on the opposite side of η from S . Since Vor * P (S ) 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. W η (S ) 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 O(1) time [12,13,8], inserting or deleting an anchor at a known position in an RB-tree, i.e., without searching, takes amortized O(1) time.

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 Ω(n log n) running time due to Ω(n) 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.

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 O(log 2 n) 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.

Wavefront operations
We summarize the eight wavefront operations, where K inv , A vis and I new are the numbers of involved (i.e., processed and newly created) potential vertices, visited anchors, and created diagram vertices, respectively: 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.
Remark. 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 SD and SD is O(m), so a priority queue takes O(log m) time for an insertion or an extraction.
Remark. 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.
Remark. 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 O(log m), and for an anchor, since the position is always known, the cost is amortized O (1). Similarly, we charge the cost to delete a diagram vertex at its creation time.

Initiate operation
An Initiate operation processes a triangle to compute Vor (S ) and initiate W (d,d2) (S ). 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, Vor (S ) is computed by constructing the Voronoi diagram of S without considering P and trimming the diagram by . Second, by tracing Vor (S ) along the pair (d, d 2 ) of diagonals, W (d,d2) (S ) is built as the aforementioned data structure (Section 4.1). During the tracing, if v 1,2 (resp. v 1 ) is a reflex polygon vertex, then it will be inserted into the anchor list of its located wavelet (cells) in W (d,d2) (S ). Note that if v 2 is a reflex polygon vertex, then v 2 will be inserted as an anchor after W (d,d2) (S ) has been divided or split at v 2 .
The operation time is O |S | · (log m + log 2 n) . The first step takes O(|S | log |S |) time to construct the Voronoi diagram using a standard divide-and-conquer or a plane-sweep algorithm [
This operation is equivalent to sweeping the triangle with a scan line parallel tod and processing each hit potential vertex. The next hit potential vertex is provided from the priority queue associated withd (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ṽ 1,2 is a reflex polygon vertex, the wavefront is further processed by three cases. If neitherd 1 nord 2 is a polygon side,ṽ 1,2 will later be inserted as an anchor while dividing or splitting W (d1,d2) (Q) atṽ 1,2 . If bothd 1 andd 2 are polygon sides, all the subcells will be completed along (d 1 ,d 2 ) and the wavefront will disappear. If onlyd 1 (resp.d 2 ) is a polygon side, all the subcells alongd 1 (resp.d 2 ) excluding the one containing v 1,2 will be completed, andṽ 1,2 will be inserted into the corresponding wavelet as the last or first anchor.

Merge operation
A Merge operation merges W η (Q) and W η (Q ) into W η (Q ∪ Q ) together with Vor * P (Q) ∩ and Vor * P (Q )∩ into Vor * P (Q∪Q )∩ where is the underlying triangle. In our algorithm, either Q = S , Q = S(d 1 ), and η = (d, d 2 ) or Q = S ∪ S(d 1 ), Q = S(d 2 ), and η = d; in both cases, S ⊆ Q. The border will form on ∂ \ η, i.e., d 1 for the former case and (d 1 , d 2 ) 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 Vor * P (Q) ∩ and Vor * P (Q ) ∩ into Vor * P (Q ∪ Q ) ∩ and (2) merge W η (Q) and W η (Q ) into W η (Q ∪ Q ). 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 Vor * P (Q ∪ Q ) ∩ associated with one site in Q and one site in Q ; 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 v 1 to v 1,2 for η = (d, d 2 ) and from v 1 to v 2 for η = d. By [11, Lemma 4-6], a merge curve called initial starts from v 1,2 for the latter case (η = d), but all other merge curves have both endpoints on η, and those endpoints are associated with one site in S . Let Q be the set of sites in S that have a wavelet in W η (Q). If η = (d, d 2 ), a site in Q can have two wavelets in W η (Q), and with an abuse of terminology, such a site is imagined to have two copies, each to one wavelet. Since Q ⊆ Q, finding a starting endpoint for each merge curve except the initial one is to test sites in Q following the ordering of their wavelets in W η (Q) 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 Q lying after the site inducing the stopping endpoint will be tested.
Let x be the next starting endpoint, which is unknown, and let s be the next site in Q to test. A two-level binary search on W η (Q ) determines if s induces x, and if so, further determines the site t ∈ Q that induces x with s as well as the corresponding anchor.
The first-level binary search executes on the RB-tree for the wavelets in W η (Q ), and each step determines for a site q ∈ Q if its cell lies before or after t's cell along η or if q = t. Let y 1 and y 2 denote the two intersection points between η and the two incomplete edges of s (in W η (Q)), where y 1 lies before y 2 along η, and let z 1 and z 2 be the two points defined similarly for q (in W η (Q )). Since s lies in , the distance between s and any point in η can be computed in O(1) time. The two incomplete edges of q will be updated until z 1 and z 2 , so that the distance from q to z 1 (resp. to z 2 ) can be computed from the corresponding anchor. For example, if u is the anchor of the subcell that contains z 1 , d(z 1 , q) = |z 1 u| + d(u, q).
The determination considers four cases. (Assume s and t induce the "starting" endpoint.) z 2 lies before y 1 (resp. z 1 lies after y 2 ): t's cell lies after (resp. before) q's cell. z 1 lies before y 1 and z 2 lies between y 1 and y 2 (resp. z 2 lies after y 2 and z 1 lies between y 1 and y 2 ): if z 2 is closer to q than to s (resp. z 1 is closer to q than to s), then t's cell lies after (resp. before) q's cell; otherwise, t is q. Both y 1 and y 2 lie between z 1 and z 2 : t is q. Both z 1 and z 2 lie between y 1 and y 2 : let x s be the projection point of s onto η. If x s lies before z 1 , then t's cell lies before q's cell. If x s lies after z 2 : if z 2 is closer to q than to s, t's cell lies after q's cell; if both z 1 and z 2 are closer to s than to q, t's cell lies before q's cell; otherwise, t = q. If x s lies between z 1 and z 2 : if z 1 is closer to s than to q, t's cell lies before q's cell; otherwise, t = q. If the first-level search does not find t, then s does not induce the next starting endpoint x.
The second-level binary search executes on the RB-tree for t's anchor list to either determine the next starting endpoint x and t's corresponding anchor or report that s does not induce x. Let u be the current anchor of t to test, and let x s be the projection point of s onto η. u's "interval" on η can be decided by checking u's two neighboring anchors. If u's interval lies after x s , x lies before u's interval; otherwise, if both endpoints of u's interval are closer to (resp. farther from) t than to (resp. from) s, x lies after (resp. before) u's interval, and if one endpoint is closer to t than to s but the other is not, then x lies in u's interval and can be computed in O(1) time since d(x, s) = |xu| + d(u, t). If the second-level binary search does not find such an interval, s does not induce the next starting endpoint x.
The second phase (i.e., merging W η (Q) and W η (Q ) into W η (Q ∪ Q )) splits W η (Q) and W η (Q ) at the endpoints of merge curves, and concatenates active parts at these endpoints where a part is called active if it contributes to W η (Q ∪ Q ). In fact, the active parts along η alternately come from W η (Q) and W η (Q ). 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 W η (Q ∪ Q ), it will be inserted into its located wavelet as the first or the last anchor.

Tracing merge Curves in Merge Operation
We make some further definitions. The bisector B(Q, Q ) between Q and Q is the collection of points with the same minimum geodesic distance from both Q and Q , namely B(Q, Q ) = {x ∈ P | min q∈Q d(x, q) = min q ∈Q d(x, q )}. In our algorithm, each anchor u refers to either a reflex polygon vertex of P or a point site in S; we store its associated site s(u), its geodesic distance w(u) from s(u), and its anchor a(u) with respect to s(u). The weighted distance d w (x, u) from a point x to u is |xu| + w(u), and the weighted bisector between two anchors u 1 and u 2 is {x ∈ P | d w (x, u 1 ) = d w (x, u 2 )}.
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 Q than to the current site in Q ; 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 O(I new ) time in total, where I new is the number of newly created diagram vertices.

Join operation
In the construction of SD (Section 6.2), our algorithm will actually join W d1 S(d 2 ) ∪ S and W d1 (S (d)) into W d1 (S (d 1 )), and join W d2 S(d 1 ) ∪ S and W d2 (S (d)) into W d2 (S (d 2 )). To interpret a Join operation, for the former case, one would replace d, d 3 At a high level, a Join operation first builds the border on d, 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 b (d) denote the collection of connected components of border edges on d in SD . 1 These components in b (d) are called joint segments and are ordered from v 2 to v 1 . By [11,, b (d) satisfies the following conditions: At most one joint segment starts at v 2 , and such a joint segment is called initial. Except the initial one, both endpoints of each joint segment are associated with a site in S . 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 S(d 3 ) and S(d 4 ) do not "cross" each other along d, and this desirable condition enables locating a starting endpoint of a non-initial joint segment by searching W d S (d 4 ) with sites in S .
Assume that d is oriented from v 2 to v 1 . Joint segments will be constructed one by one from v 2 to v 1 . For each joint segment (except the initial one), its starting endpoint is first located on d by searching W d S (d 4 ) and then it is traced from its starting endpoint along d until reaching its stopping endpoint. We remark that the initial joint segment will be directly traced from v 2 .
Let S |d be the set of sites in S that have a wavelet in W d S(d 3 )∪S . By the second condition, finding the starting endpoint of a joint segment is to "locate" sites of S |d in W d S (d 4 ) since each endpoint of a joint segment (except the initial one) is associated with a site in S |d . Therefore, sites in S |d will be tested following the ordering of their wavelets from v 2 to v 1 ; after tracing a joint segment, the test will continue on the next site in S |d , which lies after the site in S |d associated the last traced stopping endpoint. Note that each site in S |d has exactly one wavelet in Let x be the next starting endpoint, which is unknown, and let s be the site in S |d to test. A two-level binary search on W d S (d 4 ) , which is identical to the one in Section 5.3, can determine if s induces x, and if so, determine the site t ∈ S (d 4 ) that induces x with s together with the corresponding anchor. The total time to find starting points is trivially O |S |(log n + log m) .
The process to trace a joint segment from the starting endpoint is simpler than tracing a merge curve in Section 5. After constructing joint segments, we split each of W d S(d 3 ) ∪ S and W d S (d 4 ) at the endpoints of these joint segments, and concatenate the active parts at these endpoints into W d S (d) , where "active" means having a wavelet in W d S (d) . For each ending polygon vertex of d, i.e., v 1 or v 2 , if it is reflex but has not been an anchor of W d S (d) , it will be inserted into its located wavelet as the first or the last anchor. The

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 O(log n + log m). Since a wavefront has O(m) wavelets and a wavelet has O(n) anchors, both locating the common polygon vertex and spliting the wavelet and wavefront take O(log n + log m) time. Although inserting the common polygon vertex takes amortized O(1) time, since the worst-case time to insert an anchor is O(log n), the time to insert the common polygon vertex is dominated by the time-factor O(log n + log m). Since there is no new site, there is no new incomplete Voronoi edge and no new potential vertex.

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 Ω(log n + log m) time and there are Ω(n) events to separate a wavefront, if only Split operations are adopted, the total construction time would be Ω n(log n + log m) .
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 O(A vis + 1), where A vis 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 O(A vis ). 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 O(A vis ) time, and inserting the common vertex takes amortized O(1) time.

Insert operation
For a triangle , an Insert operation inserts S into W d1 S(d 2 ) to form W d1 S(d 2 ) ∪ S . An Insert operation executes in the construction of SD, but its outcome will be used to construct SD . Let T denote the intermediate site set during the Insert operation, so that T = S(d 2 ) at the beginning.
For each site s ∈ S , let x s be its vertical projection point onto d 1 , and conduct a binary search on the wavefront W d1 (T ) to locate the anchor whose subcell contains x s . During the binary search, the two incomplete Voronoi/polygonal edges of each visited wavelet will be updated up to d 1 . Recall that the weighted distance between a point y and an anchor u associated with a site t is |yu| + d(u, t).
If x s is closer to s than to the located anchor under the weighted distance, traverse anchors in W d1 (T ) from x s along each direction of d 1 until either touching a polygon vertex of d 1 or reaching a point equidistant from s and the current visited anchor under the weighted distance. If all points in the "interval" of a visited anchor on d 1 are fully closer to s 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 W d1 (T ). After the traversal, insert s into W d1 (T ), namely create and insert its wavelet into W d1 (T ). Since S(d 2 ) and S lie on the same side of d 1 and all sites in S(d 2 ) lie outside , each cell of a site in S(d 2 ) will not be separated by a cell of a site in S "along d 1 ", supporting the correctness of the above process.
After processing all the sites in S , check the two polygon vertices, v 1 and v 1,2 , of d 1 , and if v 1 (resp. v 1,2 ) should belong to a wavelet associated a site in S , insert v 1 (resp. v 1,2 ) into the wavelet as the first or last anchor. For each inserted wavelet, generate its incomplete Voronoi or polygonal edges, and compute the corresponding potential vertices. The

Propagate operation
A Propagate operation propagates a wavefront W d S (d) into P (d), i.e., from the upside to the downside of d, To some extent, a Propagate operation is a generalized version of an Extend operation underlying a sub-polygon instead of a triangle.
The operation consists of two phases. The first phase constructs Vor P (S) ∩ P (d), and the second phase refines each cell into subcells to obtain Vor * P (S) ∩ P (d). The first phase "sweeps" the triangles in P (d) by a preorder traversal of the subtree of T rooted at (d). Similar to the Extend operation, this sweep processes potential vertices inside each triangle to construct Voronoi edges and to update the wavefront accordingly. However, this sweep will not process the polygon vertices of P (d), so that the anchor lists will not be updated, preventing constructing a Voronoi edge using the two corresponding anchor lists. Fortunately, Oh and Ahn [10] gave another technique that obtains in O(log n) time the two anchor lists of a Voronoi edge provided that the two Voronoi vertices are given, so a Voronoi edge can still be built in time proportional to log n plus the number of its breakpoints.
Therefore, the first phase takes O K d (log m+log 2 n)+|Vor P (S)∩P (d)| time, where K d is the number of involved potential vertices. Note that for the Voronoi edges that intersect d, their breakpoints outside P (d) will be counted in the respective triangles in P (d).
The second phase triangulates each cell in Vor P (S) ∩ P (d) and constructs the SPM from the associated site in each triangulated cell. Since Chazelle's algorithm [3] takes linear time to triangulate a simple polygon and Guibas et al's algorithm [6] takes linear time to build the SPM in a triangulated polygon, the second phase takes O(|SD ∩ P (d)|) time.
Remark. Although all the sites lie outside P (d), the information stored in W d S (d) allows us not to conduct Guibas et al's algorithm from scratch. For example, for each anchor u, its anchor a(u) is also stored, and the SPM edge of u is the line segment from u to the polygon boundary along the direction from a(u) to u.

58:15
To sum up, a Propagate operation takes O K d (log m + log 2 n) + |SD ∩ P (d)| time.

6
Subdivision construction

Construction of SD
To construct SD, we process each triangle by the postorder traversal of the rooted partition tree T and build SD ∩ . We first assume that no diagonal of is a polygon side, and we will discuss the other cases later. Let d be the root diagonal of and adopt the convention in Section 2.2 and Section 2.3. When processing , since its two children, (d 1 ) and (d 2 ), have been processed, W d1 S(d 1 ) and W d2 S(d 2 ) are available.
The processing of each triangle consists of 8 steps: and W d2 S(d 1 ) ∪ S will be used to construct SD . If exactly one diagonal d of is not a polygon side, it is either the root triangle or a leaf triangle. For the former, compute Vor (S ), extend W d (S(d)) into to build Vor * P (S(d)) ∩ , and merge Vor (S ) and Vor * P (S(d)) ∩ into Vor * P (S) ∩ = SD ∩ ; for the latter, compute Vor (S ) and initiate W d (S  To apply Lemma 3, we bound K , I and A by the following two lemmas. Proof. We consider Step 4 (divide W (d,d2) S(d 1 ) ∪ S along d 2 ), which is similar to Step 6. Since there are O(n + m) anchors, it is sufficient to bound the number of anchors that are visited by Step 4 but still involved in the future construction of SD, namely SD ∩ P (d).
Since the subcell of each visited anchor intersects d 2 , if the subcell does not intersect d, its anchor will not be involved in constructing SD ∩ P (d). A subcell of a visited anchor intersects d in two cases. In the first case, the subcell contains v 2 and thus intersects both d and d 2 . Since there are O(n) triangles, the total number for the first case is O(n). In the second case, the subcell intersects both d and d 2 but does not contain v 2 . By the definition of W (d,d2) S(d 1 ) ∪ S , its associated "site" belongs to either S or S(d 1 ). For the former, since its anchor must be the site itself, the total number is S = m. For the latter, since all the sites in S(d 1 ) lie outside , only one site in S(d 1 ) can own a cell intersecting both d and d 2 . Moreover, due to the visibility of a subcell, only one subcell in such a cell can intersect both d and d 2 . Since there are O(n) triangles, the total number is O(n).
By Lemma 3, 4, and 5, we conclude the construction time of SD as follows.
We first analyze Split and Propagate operations. Step 2), although each of the two resultant wavefronts could be joined with another wavefront in Step 3, since the wavefront along the traversed diagonal will propagate into the corresponding sub-polygon in Step 4 and will not separate anymore, the reasoning of Lemma 5 implies that the total time for all the Divide operations is O(n + m). Finally, by Lemma 8, the total time for all the Propagate operations (Step 4) is O n + m(log m + log 2 n) , leading to the statement.