Computational Geometry,Fundamental Structures:Voronoi Diagrams

Voronoi Diagrams

Recall the post-office problem mentioned in the introduction. Here we want to preprocess a set S of n points, referred to as sites, in the plane such that we can answer the query: which site in S is closest to a given query point q? In order to solve this problem we can divide the plane into regions according to the nearest-neighbor rule: each site s gets assigned the region which is closest to s. This subdivision, which compactly captures the distance information inherent in a given configuration, is called the Voronoi diagram of S—see Fig. 62.7(a).

More formally, the Voronoi diagram of a set of sites S = {s1,..., sn} in Rd, which we refer to as Vor(S), partitions space into n regions—one for each site—such that the region for a site si consists of all points that are closer to si than to any other site sj ∈ S. The set of points that are closest to a particular site si forms the so-called Voronoi cell of si, and is denoted by V (si). Thus, when S is a set of sites in the plane we have

where dist(., .) denotes the Euclidean distance.

Now consider the dual graph of the Voronoi diagram, that is, the graph that has a node for every Voronoi cell and an arc between any two Voronoi cells that share a common edge—see Fig. 62.7(b). (Observe that the concept of dual graph used here has nothing to do with the duality transform discussed in Section 62.2.3.) Suppose we embed this graph in the plane, by using the site si to represent the node corresponding to the cell V (si) and by drawing the edges as straight line segments, as in Fig. 62.7(c). Somewhat surprisingly perhaps, this graph is always planar. Moreover, it is actually a triangulation of the point set S, assuming that no four points in S are co-circular. More details on this special triangulation, which is called the Delaunay triangulation, will be given in Section 62.5.1.

There exists a fascinating connection between Voronoi diagrams in Rd and half-space intersections in Rd+1. Assume for simplicity that d = 2, and consider the transformation that maps a site s = (sx, sy ) in R2 to the non-vertical plane h(s) : z = 2sxx+2sy y −(s2 +s2 ) in . Geometrically, h(s) is the plane tangent to the unit paraboloid z = x +y at the point vertically above (sx, sy , 0). Let H(S) be the set of planes that are the image of a set of point sites S in the plane. Let S denote the convex polyhedron that is formed by the intersection of the positive half-spaces defined by the planes in ( ), that is, = n h∈H(S) h+, where h+ denotes the half-space above h. Surprisingly, the projection of the edges and vertices of S vertically downward on the xy-plane is exactly the Voronoi diagram of S.

The Voronoi diagram can be defined for various sets of sites, for example points, line segments, circles, or circular arcs, and for various metrics. Sections 62.4.1 and 62.4.2 discuss the complexity and construction algorithms of Voronoi diagrams for the usual case of point sites in Rd, while Section 62.4.3 describes some of the possible variations. Additional details and proofs of the material presented in this section can be found in [4, 5, 22].

Complexity

Let S be a set of n points in Rd. The Voronoi diagram of S is a cell complex in Rd. If d = 2 then the Voronoi cell of a site is the interior of a convex, possibly infinite polygon. Its boundary consists of Voronoi edges, which are equidistant from two sites, and Voronoi vertices, which are equidistant from at least three sites. The Voronoi diagram of n ≥ 3 sites has at most 2n − 5 vertices and at most 3n − 6 edges, which implies that the Delaunay triangulation has at most 2n − 5 triangles and 3n − 6 edges. What is the complexity of the Voronoi diagram in d ≥ 3? Here the connection between Voronoi diagrams in Rd and intersections of half-spaces in Rd+1 comes in handy: we know from the Upper Bound Theorem that the intersection of n half-spaces in Rd+1 has complexity O(nl(d+1)/2J ) = O(nld/21). This bound is tight, so the maximum complexity of the Voronoi diagram—and of the Delaunay triangulation, for that matter—in Rd is Θ(nld/21).

Construction

In this section we present several algorithms to compute the Voronoi diagram or its dual, the Delaunay triangulation, for point sites in Rd. Several of these algorithms can be generalized to work with metrics other than the Euclidean, and to sites other than points. Since the Voronoi diagram in the plane has only linear complexity one might be tempted to search for a linear time construction algorithm. However the problem of sorting n real numbers is reducible to the problem of computing Voronoi diagrams and, hence, any algorithm for computing the Voronoi diagram must take Ω(n log n) time in the worst case.

Two data structures that are particularly well suited to working with planar subdivi- sions like the Voronoi diagram are the doubly-connected edge list (DCEL) by Muller and Preparata [38] and the quad-edge structure by Guibas and Stolfi [29]. Both structures re- quire O(n) space, where n is the complexity of the subdivision, and allow to efficiently traverse the edges adjacent to a vertex and the edges bounding a face. In both structures, one can easily obtain in O(n) time a structure for the Voronoi diagram from a structure for the Delaunay triangulation, and vice versa. In fact, the quad-edge structure, as well as a variant of the DCEL [7], are representations of a planar subdivision and its dual at the same time, so there is nothing to convert (except, perhaps, that one may wish to explicitly compute the coordinates of the vertices in the Voronoi diagram, if the quad-edge structure or DCEL describes the Delaunay triangulation).

Divide-and-conquer. The first deterministic worst-case optimal algorithm to compute the Voronoi diagram was given by Shamos and Hoey [52]. Their divide-and-conquer algorithm splits the set S of point sites by a dividing line into subsets L and R of approximately the same size. The Voronoi diagrams Vor(L) and Vor(R) are computed recursively and then merged into Vor(S) in linear time. This algorithm constructs the Voronoi diagram of a set of n points in the plane in O(n log n) time and linear space.

Plane Sweep. The strategy of a sweep line algorithm is to move a line, called the sweep line, from top to bottom over the plane. While the sweep is performed, information is maintained regarding the structure one wants to compute. It is tempting to apply the same approach to Voronoi diagrams, by keeping track of the Voronoi edges that are currently intersected by the sweep line. It is problematic, however, to discover a new Voronoi region in time: when the sweep line reaches a new site, then the line has actually been intersecting the Voronoi edges of its region for a while. Fortune [23] was the first to find a way around this problem. Fortune’s algorithm applies the plane sweep paradigm in a slightly different fashion: instead of maintaining the intersection of the sweep line with the Voronoi diagram, it maintains information of the part of the Voronoi diagram of the sites above the line that can not be changed by sites below it. This algorithm provides an alternative way of computing the Voronoi diagram of n points in the plane in O(n log n) time and linear space.

Randomized incremental construction. A natural idea is to construct the Voronoi diagram by incremental insertion, that is, to obtain Vor(S) from Vor(S\{s}) by inserting the site s. Insertion of a site means integrating its Voronoi region into the diagram constructed so far. Unfortunately the region of s can have up to n − 1 edges, for |S| = n, which may lead to a running time of O(n2). The insertion process is probably better described and implemented in the dual environment, for the Delaunay triangulation DT : construct

DTi = DT ({s1,..., si}) by inserting si into DTi−1. The advantage of this approach over a direct construction of Vor(S) is that Voronoi vertices that appear only in intermediate diagrams but not in the final one need not be computed or stored. DTi is constructed by exchanging edges, using edge flips [36], until all edges invalidated by si have been removed.

Still, the worst-case running time of this algorithm can be quadratic. However, if we insert the sites in random order, and the algorithm is implemented carefully, then one can prove that the expected running time is O(n log n), and that the expected amount of storage is O(n).

Other approaches. Finally, recall the connection between Delaunay triangulations and convex hulls. Since there exist O(n log n) algorithms to compute the convex hull of points in R3 (see Section 62.3) we therefore have yet another optimal algorithm for the computation of Voronoi diagrams.

Variations

In this section we present some of the common variations on Voronoi diagrams. The first is the order-k Voronoi diagram of a set S of n sites, which partitions Rd on the basis of the first k closest point sites. In other words, each cell in the order-k Voronoi diagram of a set S of sites in the plane corresponds to a k-tuple of sites from S and consists of those points in the plane for which that k-tuple are the k closest sites. One might fear that the order-k Voronoi diagram has Θ(nk) cells, but this is not the case. In two dimensions, for example, its complexity is O(k(n−k)), and it can be computed in O(k(n−k) log n+n log3 n) expected time [2].

The furthest-site Voronoi diagram partitions Rd according to the furthest site, or equivalently, according to the closest n − 1 of n sites. The furthest-site Voronoi diagram can be computed in O(n log n) time in two dimensions, and in O(nld/21 ) in dimension d ≥ 3.

One can also consider different distance functions than the Euclidean distance. For ex- ample, one can alter the distance function by the addition of additive or multiplicative weights. In this case every point site si is associated with a weight wi and the distance function d(si, x) between a point x and a site si becomes d(si, x) = wi + dist(si, x) (additive weights) or d(si, x) = wi · dist(si, x) where dist(si, x) denotes the Euclidean distance between si and x. The Voronoi diagram for point sites in 2 dimensions with additive weights can be computed in O(n log n) time, for multiplicative weights the time increases to O(n2) time.

Finally the power diagram, or Laguerre diagram, is another Voronoi diagram for point sites si that are associated with weights wi. Here the distance function is the power distance intro- duced by Aurenhammer [3], where the distance from a point x to a site si is measured alon ga line tangent to the sphere of radius √wi centered at si, i.e., d(si, x) = jdist(si, x)2 − wi.

The power diagram can be computed in O(n log n) time in two dimensions.