Improved Dynamic Geodesic Nearest Neighbor Searching in a Simple Polygon

03/15/2018
by   Pankaj K. Agarwal, et al.
0

We present an efficient dynamic data structure that supports geodesic nearest neighbor queries for a set S of point sites in a static simple polygon P. Our data structure allows us to insert a new site in S, delete a site from S, and ask for the site in S closest to an arbitrary query point q ∈ P. All distances are measured using the geodesic distance, that is, the length of the shortest path that is completely contained in P. Our data structure achieves polylogarithmic update and query times, and uses O(n^3n m + m) space, where n is the number of sites in S and m is the number of vertices in P. The crucial ingredient in our data structure is an implicit representation of a vertical shallow cutting of the geodesic distance functions. We show that such an implicit representation exists, and that we can compute it efficiently.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

07/10/2017

Dynamic Geodesic Nearest Neighbor Searching in a Simple Polygon

We present an efficient dynamic data structure that supports geodesic ne...
09/24/2021

Dynamic Data Structures for k-Nearest Neighbor Queries

Our aim is to develop dynamic data structures that support k-nearest nei...
04/17/2020

On the Approximate Nearest Neighbor Queries among Curves under the Fréchet Distance

Approximate nearest neighbor search (ANNS) is a long-studied problem in ...
12/12/2014

Representing Data by a Mixture of Activated Simplices

We present a new model which represents data as a mixture of simplices. ...
04/30/2022

Chromatic k-Nearest Neighbor Queries

Let P be a set of n colored points. We develop efficient data structures...
05/05/2018

On k Nearest Neighbor Queries in the Plane for General Distance Functions

We study k nearest neighbor queries in the plane for general (convex, pa...
10/21/2020

On Adaptive Distance Estimation

We provide a static data structure for distance estimation which support...
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

Nearest neighbor searching is a classic problem in computational geometry in which we are given a set of point sites , and we wish to preprocess these points such that for a query point , we can efficiently find the site closest to . We consider the case where is a dynamic set of points inside a simple polygon . That is, we may insert a new site into or delete an existing one. We measure the distance between two points and by the length of the geodesic , that is, the shortest path connecting and that is completely contained in . We refer to this distance as the geodesic distance .

Related work.

It is well known that if we have only a fixed set of sites, we can answer nearest neighbor queries efficiently by computing the Voronoi diagram of and preprocessing it for planar point location. This requires preprocessing time, the resulting data structure uses linear space, and we can answer queries in time. Voronoi diagrams have also been studied in case the set of sites is restricted to lie in a simple polygon , and we measure the distance between two points and by their geodesic distance  [3, 27, 18, 25]. The approach of Hershberger and Suri [18] computes the geodesic Voronoi diagram in time, where is the total number of vertices in the polygon , and is applicable even if has holes. Very recently, Oh and Ahn [25] presented an time algorithm. When this improves the previous results. All these approaches allow for time nearest neighbor queries. However, they are efficient only when the set of sites is fixed, as inserting or deleting even a single site may cause a linear number of changes in the Voronoi diagram.

To support nearest neighbor queries, it is, however, not necessary to explicitly maintain the (geodesic) Voronoi diagram. Bentley and Saxe [5] show that nearest neighbor searching is a decomposable search problem. That is, we can find the answer to a query by splitting into groups, computing the solution for each group individually, and taking the solution that is best over all groups. This observation has been used in several other approaches for nearest neighbor searching with the Euclidean distance [1, 11, 7]. However, even with this observation, it is hard to get both polylogarithmic update and query time. Chan [7] was the first to achieve this. His data structure can answer Euclidean nearest neighbor queries in time, and supports insertions and deletions in and amortized time, respectively. Recently, Kaplan et al. [19] extended the result of Chan to more general, constant complexity, distance functions.

Unfortunately, the above results do not directly lead to an efficient solution to our problem. The function describing the geodesic distance may have complexity , and thus the results of Kaplan et al. [19] do not apply. Moreover, even directly combining the decomposable search problem approach with the static geodesic Voronoi diagrams described above does not lead to an efficient solution, since every update incurs an cost corresponding to the complexity of the polygon. Only the very recent algorithm of Oh and Ahn [25] can be made amendable to such an approach. This results in an size data structure with query time and updates. Independently from Oh and Ahn, we developed a different data structure yielding similar results [2]. The core idea in both data structures is to represent the Voronoi diagram implicitly. Moreover, both approaches use similar primitives. In this manuscript, we build on the ideas from our earlier work, and significantly extend them to achieve polylogarithmic update and query times.

Our results.

We develop a fully dynamic data structure to support nearest neighbor queries for a set of sites inside a (static) simple polygon . Our data structure allows us to locate the site in closest to a query point , to insert a new site into , and to delete a site from . Our data structure supports queries in time, insertions in amortized expected time, and deletions in amortized expected time. The space usage is .

Furthermore, we show that using a subset of the tools and techniques that we develop, we can build an improved data structure for when there are no deletions. In this insertion-only setting, queries take worst-case time, and insertions take amortized time. We can also achieve these running times in case there are both insertions and deletions, but the order of these operations is known in advance. The space usage of this version is .

2 An overview of the approach

As in previous work on geodesic Voronoi diagrams [3, 27], we assume that and are in general position. That is, (i) no two sites and in (ever) have the same geodesic distance to a vertex of , and (ii) no three points (either sites or vertices) are colinear. Note that (i) implies that no bisector between sites and contains a vertex of .

Throughout the paper we will assume that the input polygon has been preprocessed for two-point shortest path queries using the data structure by Guibas and Hershberger [13] (see also the follow up note of Hershberger [17]). This takes time and allows us to compute the geodesic distance between any pair of query points in time.

Dynamic Euclidean nearest neighbor searching.

We briefly review the data structures for dynamic Euclidean nearest neighbor searching of Chan [7] and Kaplan et al. [19], and the concepts they use, as we will build on these results.

Let be a set of bivariate functions, and let denote the arrangement of the (graphs of the) functions in . In the remainder of the paper we will no longer distinguish between a function and its graph. A point has level if the number of functions in that pass strictly below is . The at most -level is the set of all points in for which the level is at most , and the -level is the boundary of that region.

Consider a collection of points in (e.g. a line segment), and let denote the set of functions from intersecting . We refer to as the conflict list of . Furthermore, let denote the vertical (downward) projection of onto the -plane.

A pseudo-prism is a constant complexity region in that is bounded from above by a function, unbounded from below, and whose sides are surfaces vertical with respect to the -direction. A a -shallow -cutting is a collection of such pseudo-prisms with pairwise disjoint interiors whose union covers and for which each pseudo-prism intersects at most functions in  [24]. Hence, for each region (pseudo-prism) , the conflict list contains functions. The number of regions in is the size of the cutting. Matoušek [24] originally defined a shallow cutting in terms of simplicies, however using pseudo-prisms is more convenient in our setting. In this case the parameter cannot become arbitrarily large, however we are mostly interested in -shallow -cuttings. In the remainder of the paper we simply refer to such a cutting as a -shallow cutting. Observe that each pseudo-prism in is intersected by functions.

Let be the distance from to . The data structures of Kaplan et al. [19] and Chan [7] actually maintain the lower envelope of the set of distance functions . To find the site closest to a query point they can simply query the data structure to find the function that realizes at . The critical ingredient in both data structures, as well as our own data structure, is an efficient algorithm to construct a shallow cutting of the functions in . Chan and Tsakalidis [8] show that if is a set of linear functions (i.e. planes in ), we can compute a -shallow cutting of size in time. Kaplan et al. [19] show how to compute a -shallow cutting of size , in time for a certain class of constant complexity algebraic functions . Since the geodesic distance function of a single site may already have complexity , any -shallow cutting of such functions may have size . To circumvent this issue, we allow the regions (pseudo-prisms) to have non-constant complexity. We do this in such a way that we can compactly represent each pseudo-prism, while still retaining some nice properties such as efficiently testing if a point lies inside it. Hence, in the remainder of the paper we will drop the requirement that the regions in a cutting need to have constant complexity.

The general approach.

The general idea in our approach is to recursively partition the polygon into two roughly equal size sub-polygons and that are separated by a diagonal. For the sites in the “left” subpolygon , we then consider their geodesic distance functions restricted to the “right” subpolygon , that is, , for . The crucial part, and our main contribution, is that for these functions we can represent a vertical shallow cutting implicitly. See Fig. 1 for a schematic illustration. More specifically, in expected time, we can build a representation of the shallow cutting of size . We can then use this algorithm for building implicitly represented shallow cuttings in the data structure of Chan [7] and Kaplan et al. [19]. That is, we build and maintain the lower envelope . Symmetrically, for the sites in , we maintain the lower envelope that their distance functions induce in . When we get a query point , we use to find the site in closest to in time. To find the site in closest to , we recursively query in sub-polygon . In total we query in levels, leading to an query time. When we add or remove a site we, insert or remove its distance function in lower envelope data structures (one at every level). Every insertion takes amortized expected time, and every deletion takes amortized expected time. Since every site is stored times, this leads to the following main result.

Let be a simple polygon with vertices. There is a fully dynamic data structure of size that maintains a set of point sites in and allows for geodesic nearest neighbor queries in worst case time. Inserting a site takes amortized expected time, and deleting a site takes amortized expected time.

Figure 1: A schematic drawing of the downward projection of an implicit -shallow cutting in . The faces are pseudo-trapezoids. Neighboring pseudo-trapezoids share a vertical segment, or part of a bisector. We store only the degree one and three vertices (fat), and their topology.

The main complexity is in developing our implicit representation of the -shallow cutting, and the algorithm to construct such a cutting. Once we have this algorithm we can directly plug it in into the data structure of Chan [7] and Kaplan et al. [19]. Our global strategy is similar to that of Kaplan et al. [19]: we compute an approximate -level of –in our case an implicit representation of this approximate -level– and then argue that, under certain conditions, this approximate -level is actually a -shallow cutting of . Our approximate -level will be a -level, for some appropriate , on a random sample of the functions in . So, that leaves us two problems: i) computing an implicit representation of a -level, and ii) computing the conflict lists for all pseudo-prism in our cutting .

For problem i), representing the -level implicitly, we use the connection between the -level and the -order Voronoi diagram. In Section 3 we describe a small, implicit representation of the -order Voronoi diagram that still allows us to answer point location queries efficiently. Initially, we use the recent algorithm of Oh and Ahn [25] to construct this representation. In Section 5 we then design an improved algorithm for our particular use case.

For problem ii), computing the conflict lists, we will use a data structure developed by Chan [6] together with the implicit Voronoi diagrams that we developed in Section 5. We describe these results in more detail in Section 6. In Section 7, we then show in detail how we can combine all the different parts into a fully dynamic data structure for nearest neighbor queries. Finally, we describe our simpler data structures for insertion-only and offline-updates in Section 8.

3 Implicit representations

Let denote the set of geodesic distance functions inside the entire polygon . Our implicit representation of a -shallow cutting is based on an implicit representation of the -level in . To this end, we first study higher order geodesic Voronoi diagrams.

Higher order Voronoi diagrams.

Consider a set of sites , a domain , and a subset of size . The -order Voronoi region is the region in in which the points are closer to (a site in) , with respect to some distance metric, than to any other subset of size . The -order Voronoi diagram is the partition of into such maximal regions over all subsets of size  [22]. Liu et al. [23] study the geodesic -order Voronoi diagram. They show that has complexity . In particular, it consists of degree one and degree three vertices, and degree two vertices (and by our general position assumption, there are no vertices of degree more than three).

Consider a Voronoi region . Let , be the edges bounding , let be the set of sites defining the other Voronoi region incident to , and observe that contains a single site  [22]. Let be the set of sites neighboring . Observe that these results imply that two adjacent regions and in are separated by a part of a bisector , where and .

By combining the above two observations we can represent implicitly. That is, we store only the locations of these degree one and degree three vertices, the adjacency relations between the regions, and the pair of labels corresponding to each pair of neighboring regions. See also the recent result of Oh and Ahn [25]. It follows that the size of this representation is linear in the number of degree one and degree three vertices of . We refer to this as the topological complexity of .

Representing the -level.

Consider the partition of into maximally connected regions in which all points in a region have the same nearest site in . Observe that this partition corresponds to the downward projection of the -level . As we argue next, this partition is closely related to the -order Voronoi diagram defined above.

Lee observes that there is a relation between the -order Voronoi diagram and the -order Voronoi diagram [22]. In particular, he shows that we can partition each -order Voronoi region into order Voronoi regions by intersecting with the (first order) Voronoi diagram of the set of sites neighboring . More specifically:

Observation .

Let be a geodesic -order Voronoi region, and let be the sites neighboring . For any point , the -closest site from is the site for which .

The topological complexity of is .

Proof.

By Observation 3, we can obtain from by partitioning every region using the Voronoi diagram of the sites neighboring , and merging regions that have the same nearest neighbor in the resulting subdivision. Thus, the vertices in appear either in or in one of the newly created Voronoi diagrams.

There are degree one and degree three vertices that are also vertices in . Since the (first order) geodesic Voronoi diagram has degree one and degree three vertices, a region creates new degree one and three vertices. Summing over all faces in this then sums to . ∎

As with we can represent implicitly by storing only the locations of the degree one and three vertices and the topology of the diagram. Note that if we can efficiently locate the region of that contains a query point , we can also easily compute the -coordinate of at , simply by computing the geodesic distance . Since we preprocessed for two-point shortest path queries this takes only time (in addition to locating the region of containing ).

Figure 2: A pseudo prism and its projection (darker red) onto .
Representing a vertical decomposition.

For every degree three and degree one vertex in we now extend a vertical segment up and down until it hits another edge of . Observe that the topological complexity of the resulting implicit vertical decomposition that we obtain still has topological complexity . Furthermore, each region in is a pseudo-trapezoid that is bounded on the left and right either by a vertical segment or a piece of polygon boundary, and on the top and bottom by pieces of bisectors or a piece of polygon boundary. We refer to the four degree one or degree three vertices on the boundary of as the corners of . See Fig. 2 for an illustration. In the remainder of the paper, we will no longer distinguish between and its implicit representation. In Section 6, we will use such an implicit vertical decomposition to obtain an implicit representation of a shallow cutting.

Computing implicit representations.

We can use the algorithm of Oh and Ahn [25] to compute the implicit representation of in time. To compute (the representation of) we first construct , and for each region in , we again use their algorithm to compute the Voronoi diagrams of the set neighbors . We then clip these diagrams to . This clipping can be done by a breadth first search in , starting with one of the vertices that is also in . This takes time in total. Summing over all faces in gives us again a running time of . Finally, to compute we need to insert two vertical extension segments at each vertex of . We can find the other endpoint of each extension segment using a point location query. Thus, we obtain the following result.

An implicit representation of the -level that uses space, can be computed in time. Using this representation, the pseudo-prism containing a query point (if it exists) can be determined in time.

In Section 5 we will show that if the sites defining the functions in lie in one half of the polygon and we restrict the functions to the other half we can improve these results. Moreover, we can then compute an implicit representation of a -shallow cutting of .

4 Approximating the -level

An -monotone surface is an -approximation of if and only if lies between and . Following the same idea as in Kaplan et al. [19] we construct an -approximation of as follows. We choose a random sample of of size and consider the -level of for some randomly chosen level in the range . Here and are some constants, , and is the desired approximation ratio. We now argue that is an -approximation of .

Consider the range space , where each range in is the subset of functions of intersected by a downward vertical ray in the -direction. See Har-Peled [15] for details on range spaces. An important concept is the VC-dimension of , defined as the size of the largest subset for which the number of sets in is .

[Lemma 2.3.5 of Aronov et al. [4]] Let , , and be three sites in . Their bisectors and intersect in at most a single point.

The VC-dimension of the range space is finite.

Proof.

The range space is equivalent to where is the family of subsets of that lie within some geodesic distance of some point . That is . We now observe that any three points , , and , define at most one geodesic disk (Lemma 4). Namely, the disk centered at the intersection point and radius . This means the same argument used by as Har-Peled [15, Lemma 5.15] now gives us that the shattering dimension of (see [15]) is constant (three). It then follows that the VC-dimension of is finite. ∎

Since has finite VC-dimension (Lemma 4), and has size , for and for some sufficiently large constant

, it follows that with high probability,

is a relative -approximation for  [16]. So, for every range , we have (whp.) that

(1)

Using exactly the same argument as Kaplan et al. [19] we then obtain the following result.

The level is an -approximation of the -level .

What remains is to show that the (expected) topological complexity of the -level is small, that is, that the expected number of degree three and degree one vertices is at most .

The expected topological complexity of level is .

Proof.

The lower envelope of has topological complexity , so by Clarkson and Shor [9] the total topological complexity of all levels in the range is . Hence, the expected topological complexity of a level , with randomly chosen from this range, is . Substituting and , we get as claimed. ∎

Again as in Kaplan et al. [19], if , we can skip the sampling of the set , and directly take a random level in in the range . This level has also an expected topological complexity of at most . Using Lemma 3 (and restarting the computation if the size of the cutting exceeds ) we then get:

An -approximation of the -level of that has topological complexity can be computed in expected time.

The main difference between our approach and that of Kaplan et al. [19] is the range space used. In our approach, the ranges are defined by downward vertical rays, whereas in Kaplan et al. the ranges are defined by more general objects. For example, their range space includes a range consisting of the functions intersected by some other constant complexity algebraic function. This allows them to directly turn their approximate level into a shallow cutting. Unfortunately, this idea does not directly extend to the geometric setting, as the VC-dimension of such a range space may again depend on the complexity of the polygon. Therefore, we will use a different approach in Section 6.

5 Computing implicit representations in subpolygon

Consider a diagonal that splits the polygon into two subpolygons and , and assume without loss of generality that is vertical and that lies right of . We consider only the sites in , and we restrict their functions to . We now present a more efficient algorithm to compute the implicit representation of in this setting. To this end, we show that the two-point shortest path query data structure of Guibas and Hershberger [13] essentially gives us an efficient way of accessing the bisector between a pair of sites without explicitly computing it. See Appendix A. In Section 5.1 we first use this to compute an implicit representation of the Voronoi diagram of in . Building on these results, we can compute an implicit representation of the -order Voronoi diagram in (Section 5.2), the -level of in (Section 5.3), and finally an implicit vertical decomposition of in (Section 5.3).

5.1 Computing an implicit Voronoi diagram

The Voronoi diagram in is a forest [3]. We now show that we can compute (the topology of) this forest efficiently by considering it as an abstract Voronoi diagram [21]. Our forest stores only the locations of the degree one and degree three vertices and their adjacencies. This turns out to be sufficient to still answer point location queries efficiently.

Assuming that certain geometric primitives like computing the intersections between “related” bisectors take time we can construct an abstract Voronoi diagram of sites in expected time [21]. We will show that is a actually a Hamiltonian abstract voronoi diagram, which means that it can be constructed in time [20]. We show this in Section 5.1.1. In Section 5.1.2 we discuss the geometric primitives used by the algorithm of Klein and Lingas [20]; essentially computing (a representation of) the concrete Voronoi diagram of five sites. We show that we can implement these primitives in time by computing the intersection point between two “related” bisectors and . This then gives us an time algorithm for constructing . Finally, in Section 5.1.3 we argue that having only the topological structure is sufficient to find the site in closest to a query point .

5.1.1 Hamiltonian abstract Voronoi diagrams

In this section we show that we can consider as a Hamiltonian abstract Voronoi diagram. A Voronoi diagram is Hamiltonian if there is a curve –in our case the diagonal – that intersects all regions exactly once, and furthermore this holds for all subsets of the sites [20]. Let be the set of sites in that we consider, and let be the subset of sites from whose Voronoi regions intersect , and thus occur in .

The Voronoi diagram in is a Hamiltonian abstract Voronoi diagram.

Proof.

By Lemma A any bisector intersects the diagonal at most once. This implies that for any subset of sites , so in particular for , the diagonal intersects all Voronoi regions in at most once. By definition, intersects all Voronoi regions of the sites in at least once. What remains is to show that this holds for any subset of . This follows since the Voronoi region of a site with respect to a set is contained in the voronoi region of with respect to . ∎

Computing the order along .

We will use the algorithm of Klein and Lingas [20] to construct . To this end, we need the set of sites whose Voronoi regions intersect , and the order in which they do so. Next, we show that we can maintain the sites in so that we can compute this information in time.

Let denote the sites in ordered by increasing distance from the bottom-endpoint of , and let be the subset of sites whose Voronoi regions intersect , ordered along from bottom to top. For any pair of sites and , with , we have that .

Proof.

Since and both contribute Voronoi regions intersecting , their bisector must intersect in some point in between these two regions. Since it then follows that all points on below , so in particular the bottom endpoint , are closer to than to . Thus, . ∎

Lemma 5.1.1 suggests a simple iterative algorithm for extracting from .

Given , can be computed from in time.

Proof.

We consider the sites in in increasing order, while maintaining as a stack. More specifically, we maintain the invariant that when we start to process , contains exactly those sites among whose Voronoi region intersects , in order along from bottom to top.

Let be the next site that we consider, and let , for some , be the site currently at the top of the stack. We now compute the distance between and the topmost endpoint of . If this distance is larger than , it follows that the Voronoi region of does not intersect : since the bottom endpoint of is also closer to than to , all points on are closer to than to .

If is at most then the Voronoi region of intersects (since was the site among that was closest to before). Furthermore, since is closer to than to the bisector between and must intersect in some point . If this point lies above the intersection point of with the bisector between and the second site on the stack, we have found a new additional site whose Voronoi region intersects . We push onto and continue with the next site . Note that the Voronoi region of every site intersects in a single segment, and thus correctly represents all sites intersecting . If lies below then the Voronoi region of , with respect to , does not intersect . We thus pop from , and repeat the above procedure, now with at the top of the stack.

Since every site is added to and deleted from at most once the algorithm takes a total of steps. Computing takes time, and finding the intersection between and the bisector of and takes time (Lemma A). The lemma follows. ∎

We now simply maintain the sites in in a balanced binary search tree on increasing distance to the bottom endpoint of . It is easy to maintain this order in time per upate. We then extract the set of sites that have a Voronoi region intersecting , and thus , ordered along using the algorithm from Lemma 5.1.1.

5.1.2 Implementing the required geometric primitives

Figure 3: The part of , the tree representing the Hamiltonian Voronoi diagram (green), that lies inside the Voronoi region (blue) of a new site is a subtree (fat). We can compute , by exploring from a point inside . In case is the first site in the ordering along we can start from the root of the “first” tree in .

In this section we discuss how to implement the geometric primitives needed by the algorithm of Klein and Lingas [20]. They describe their algorithm in terms of the following two basic operations: (i) compute the concrete Voronoi diagram of five sites, and (ii) insert a new site into the existing Voronoi diagram . In their analysis, this first operation takes constant time, and the second operation takes time proportional to the size of that lies inside the Voronoi region of in . We observe that that to implement these operations it is sufficient to be able to compute the intersection between two “related” bisectors and –essentially computing the Voronoi diagram of three sites– and to test if a given point lies on the -side of the bisector (i.e. testing if is “closer to” than to ). We then show that in our setting we can implement these operations in time, thus leading to an time algorithm to compute .

Inserting a new site.

Klein and Lingas [20] sketch the following algorithm to insert a new site into the Hamiltonian Voronoi diagram of a set of sites . We provide some missing details of this procedure, and briefly argue that we can use it to insert into a diagram of three sites. Let denote the domain in which we are interested in (in our application, is the subpolygon ) and let be the curve that intersects all regions in . Recall that is a forest. We root all trees such that the leaves correspond to the intersections of the bisectors with . The roots of the trees now corresponds to points along the boundary of . We connect them into one tree using curves along . Now consider the Voronoi region of with respect to , and observe that is a subtree of . Therefore, if we have a starting point on that is known to lie in (and thus in ), we can compute simply by exploring . To obtain we then simply remove , and connect up the tree appropriately. See Fig. 3 for an illustration. We can test if a vertex of is part of simply by testing if lies on the -side of the bisector between and one of the sites defining . We can find the exact point where an edge of , representing a piece of a bisector leaves by computing the intersection point of with and .

We can find the starting point by considering the order of the Voronoi regions along . Let and be the predecessor and successor of in this order. Then the intersection point of with must lie in . This point corresponds to a leaf in . In case is the first site in the ordering along we start from the root of the tree that contains the bisector between the first two sites in the ordering; if this point is not on the -side of the bisector between and one of the sites defining then forms its own tree (which we then connect to the global root). We do the same when is the last point in the ordering. This procedure requires time in total (excluding the time it takes to find in the ordering of ; we already have this information when the procedure is used in the algorithm of Klein and Lingas [20]).

We use the above procedure to compute the Voronoi diagram of five sites in constant time: simply pick three of the sites , , and , ordered along , compute their Voronoi diagram by computing the intersection of and (if it exists), and insert the remaining two sites. Since the intermediate Voronoi diagrams have constant size, this takes constant time.

Computing the intersection of bisectors and .

Since is a Hamiltonian Voronoi diagram, any any pair of bisectors and , with , intersect at most once (Lemma 4). Next, we show how to compute this intersection point (if it exists).

Given and , finding the intersection point of and (if it exists) takes time.

Figure 4: We find the intersection point of the two bisectors by binary searching along .

Proof. We will find the edge of containing the intersection point by binary searching along the vertices of . Analogously we find the edge of containing . It is then easy to compute the exact location of in constant time.

Let be the starting point of , i.e. the intersection of with , and assume that is closer to than , that is, (the other case is symmetric). In our binary search, we now simply find the last vertex for which . It then follows that lies on the edge of . See Fig. 4. Using Lemma A we can access any vertex of in time. Thus, this procedure takes time in total. ∎

Note that we can easily extend the algorithm from Lemma 5.1.2 to also return the actual edges of and that intersect. With this information we can construct the cyclic order of the edges incident to the vertex of representing this intersection point. It now follows that for every group of sites in , we can compute a representation of of size in time.

5.1.3 Planar point location in

In this section we show that we can efficiently answer point location queries, and thus nearest neighbor queries using .

For , the part of the bisector that lies in is -monotone.

Figure 5: A non -monotone bisector can occur only in degenerate inputs.

Proof. Assume, by contradiction, that is not -monotone in , and let be a point on such that is a local maximum. Since is not -monotone, it intersects the vertical line through also in another point further along . Let be a point in the region enclosed by the subcurve along from to and . See Fig. 5. This means that either or is non -monotone. Assume without loss of generality that it is . It is now easy to show that must pass through . However, that means that (and thus ) touches the polygon boundary in . By the general position assumption has no points in common with other than its end points. Contradiction. ∎

Since the (restriction of the) bisectors are -monotone (Lemma 5.1.3) we can preprocess for point location using the data structure of Edelsbrunner and Stolfi [12]. Given the combinatorial embedding of , this takes time. To decide if a query point lies above or below an edge we simply compute the distances and between and the sites and defining the bisector corresponding to edge . This takes time. Point lies on the side of the site that has the shorter distance. It follows that we can preprocess in time, and locate the Voronoi region containing a query point in time. We summarize our results in the following Lemma.

Given a set of sites in , ordered by increasing distance from the bottom-endpoint of , the forest representing the Voronoi diagram of in can be computed in time. Given , finding the site closest to a query point requires time.

5.2 Computing an implicit -order Voronoi diagram

Based on the relation between the -order Voronoi diagram and the -order Voronoi diagram (see Section 3) Lee developed an iterative algorithm to compute the Euclidean -order Voronoi diagram in time. His algorithm extends to any distance metric. Since the geodesic distance is a metric, this approach, together with our algorithm from the previous section gives us a way to compute in . We obtain a time algorithm. This then results in an time algorithm for computing a decomposition of (space below the) -level into pseudo-prisms.

An implicit representation of the -th order Voronoi diagram of in can be constructed in time.

Proof.

We use Lee’s algorithm to iteratively build the -th order Voronoi diagram in . Consider a cell in the -order Voronoi diagram. We collect the set of sites neighboring by traversing , and order them on increasing distance to the bottom endpoint of . Note that each edge on corresponds to a piece of bisector , and the site that we are looking for is either or , depending on which side of our cell lies. This means that we can collect the sites in , and order them on increasing distance to the bottom endpoint of in time. We then construct (an implicit representation) of its voronoi diagram in in time (Lemma 5.1.3). Finally, we clip to . Since all intersection points of with are vertices of on (see Lee [22]), all that remains is to find these points in . We use a point location query to find one of the vertices of in , and then find the remaining vertices in using a breadth first search in . It follows that in total we spend at most time. Summing over all cells in , and all rounds, gives us a running time of ) as claimed. ∎

Similar to in Theorem 5.2 we can compute the downward projection of the -level.

5.3 Computing an implicit vertical decomposition of

We now show how to turn our implicit representation of into an implicit vertical decomposition as follows. We note that this same procedure applies for computing a vertical decomposition of .

A representation of the -level consisting of pseudo-trapezoids can be computed in time. Given a query point , the -nearest site in can be reported in time.

Proof.

For every vertex , we know the faces of incident to directly above and below . The upward and downward extension segments will be contained in these faces, respectively. For each face , we collect the vertices whose upward extension segment will be contained in , and use a simple sweep line algorithm to compute which edge of (the implicit representation of) each such extension segment hits. For each vertex we then know that the upper endpoint of its upward extension segment lies on a bisector , for some . To find the exact location of this endpoint, we use a binary search along . We use the same approach for finding the bottom endpoint of the downward extension segment. Finding the edges of the implicit representation of hit, takes time per face, summing over all faces this solves to . The final binary search to find the exact location of the endpoint takes time per point (Theorem A). It follows that we spend time to compute all extension segments. This is dominated by the time it takes to compute itself. Note that in the resulting subdivision, all faces are again monotone (ignoring the boundary of ), so we can preprocess it for efficient point location as in Section 5.1. ∎

6 An implicit shallow cutting of the geodesic distance function

Let again denote the set of geodesic distance functions that the sites in induce in . We now argue that we can compute an implicit -shallow cutting for these functions.

As in Section 4, let be our random sample of size , and let be our approximate -level of . Let be the vertical decomposition of . We now raise every pseudo-trapezoid in to the -level. Denote the result by . Let denote the conflict list of , i.e., the functions intersecting the vertical downward half-line starting in .

Let be a pseudo prism in . The conflict list of is the union of the conflict lists of its corners , i.e. .

Figure 6: Since the bisectors restricted to are -monotone it follows that if a site conflicts with a prism , it must conflict with a corner of .

Proof. Let be the function defining the ceiling of . We have that by definition, so we focus on proving . Assume by contradiction that , but . So, there is a point for which , but for all corners . Hence, all four corners lie on the “-side” of , whereas lies on the “-side” of