Higher-dimensional Voronoi diagrams in linear expected time

This work is the first to validate theoretically the suspicions of many researchers — that the “average” Voronoi diagram is combinatorially quite simple and can be constructed quickly. Specifically, assuming that dimension <italic>d</italic> is fixed, and that <italic>n</italic> input points are chosen independently from the uniform distribution on the unit <italic>d</italic>-ball, it is proved that<list><item>the expected number of simplices of the dual of the Voronoi diagram is &THgr;(<italic>n</italic>) (exact constants are derived for the high-order term), and
</item><item>a relatively simple algorithm exists for constructing the Voronoi diagram in &THgr;(<italic>n</italic>) time.
</item></list>
It is likely that the methods developed in the analysis will be applicable to other related quantities and other probability distributions.


Introduction
The Voronoi diagram is a natural and intuitively appealing structure.Named for the mathematician Voronoi [22], it has been reinvented by risearchers in several fields; in particular, meteorologists associate the tw+ dimensional version with the name Thiessen [21], and physicists honor Wigner and Seits [24] for the threedimensional version.It has been used by geologists, foresters, agriculturalists, medical researchers, geographers, crystallographers, and astronomers.Within the domain of the mathematical sciences, it is applied to simulate differential equations by finite element methods, to interpolate surfaces in geometric modeling systems, and to solve geometric problems such as finding Euclidean minimum spanning trees and largest empty circles.(Avis & Bhattacharya [l] present an extensive list of references for applications.) The Voronoi diagram of a set of points -called rites is a partition of Rd that assigns a surrounding region of "nearby" points to each of the sites.Each region is the dpolytope containing the points lying nearer to the site in its interior than to any other site.Formally, the (nearest-site) Voronoi  persphere passing through these sites contains no other sites.(We call such a hypersphere empty or site-free.)l (k + 1) sites are vertices of a k-face of some Voronoi dual simplex if and only if some empty hypetsphere passes through these sites.
l Every convex-hull face is a face of some Voronoi dual simplex.
l The Votonoi dual partitions the convex hull of K, into d-simplices.
Many researchers have considered Voronoi-diagram construction.Previous work on the two-dimensional case has been surveyed elsewhere [7].Browstow et aL [5], Finney [lo], and T anemura et al. [20] have addressed the construction of three-dimensional diagrams; none of these authors present detailed analyses of running time.
Bowyer [3] and Watson [23] describe algorithms for higher dimensions, but neither analyzes his algorithm rigorously.Bowyer gives some theoretical and empirical evidence suggesting that his algorithm requires O(nl+'ld) time on average for points uniform in a ddimensional hypercube.Watson claims O(ns-'ld) time in the worst case.Avis and Bhattacharya's algorithms for the Voronoi diagram and its dual [l] rely heavily on the simplex method for linear programming; since this method has exponential worst-case running time, they focus mainly on experimental studies of their algorithms' performance and of the expected complexity of the diagrams for points distributed uniformly in the unit hypercube.Their results suggest combinatorial complexity proportional to input size for this distribution.
Brown [4] was the first to observe a pleasing connection between Voronoi diagrams in d dimensions and convex hulls in d+ 1 dimensions that allows any convex-hull algorithm to be used to construct Voronoi diagrams.The gift-wrapping algorithm [6,2,19] may be used in O(nS,,) time, where S,, is the number of dual simplices in the result.Or Seidel's shelling algorithm [17] may be used in O(n' + S,, log n) time.Seidel [16,18] has showed that S,, can be extremely large -Q(nL(d+l)l'J)in the worst case.On the other hand, it is not difficult to construct families of problem instances for which S,, = O(n).Thus probabilistic estimates of the average value of the two quantities are useful.intensity in Rd; Meijering showed that the expected number of nearest-site Voronoi neighbors of a site depends only on d; in particular, it is 6 for d = 2 and NN 15.54 for d = 3.Such a set of sites may be thought of as an infinite set of sites drawn from a uniform distribution over all of Rd.In computational practice, however, one must deal with finite sets of sites drawn from a particular distribution, e.g., a uniform distribution on the interior of some convex body like a hypercube or hypersphere.Two sites are neighbors if and only if they lie on the surface of some ball that contains no other site.In the Poisson case, a pair of distant neighbors is always unlikely since it implies the existence of a large empty ball.In the case of a bounded set of sites, it can still be shown that sites far from the boundary of the body probably have only nearby neighbors, but some distant pairs of neighbors always occur near the boundary of the body, where most of the empty ball may he outside the support of the distribution.Thus results dealing only with the Poisson case are inadequate for the average-case analysis of algorithms.
The next section describes a new method for determining ES,,, the expected number of simplices in the dual of the Voronoi diagram.In $3 this method is ap plied to the analysis of the asymptotic behavior of ES, for sites drawn independently from the uniform distribution in the unit d-ball.94 presents a variation of the gift-wrapping algorithm that uses standard bucketing techniques.Finally, in $5, the methods of 52 are ap plied to show that the algorithm requires only linear time on average for i.u.d. points in the unit d-ball.This section describes a general method for bounding ES,, the expected number of simplices in the duals of nearest-and furthest-site Voronoi diagrams of random point sets.
The first d + 1 points zr , . . .,2&r define a &simplex with probability one.Let us first reckon the probabiiity P,, that they also define a simplex in the dual of the nearest-site Voronoi diagram.This is just the probability that the other n -d -1 points lie outside the hypersphere passing through the d + 1 points.Writing g(e) for the density function of the zi and I' for the probability content of interior of the hypersphere, we see that this probability is Meijering [13] and Gilbert [ll]  We divide the domain of integration of the integral of (2.2) into eight regions corresponding to the possible patterns of intersection of the two balls B and 24, and compute I', 8, and esimp for each case.A relatively simple case is Case 1 (Figure 1); the most difficult is perhaps Case 6 (Figure 2).
Rather lengthy calculations show that the contribution of Case 1 dominates, and that   These values are not obviously inconsistent with the values 6.31 and 25.6 found empirically by Avis & Bhattacharya for (rather smaII) samples of 1000 points chosen from the unit hypercube [l,Table 11; it is reasonable to conjecture that (3.1) in fact holds for point sets chosen from a uniform distribution on any convex body.

A Fast Algorithm for the Unit d-Ball
It is immediate from Theorem 1 and the connection to convex hulls mentioned in the introduction that the Voronoi diagram of random points from a d-ball can be constructed in O(3) time on average by either the shelling algorithm or the gift-wrapping aigorithm.In this section, we describe a faster algorithm requiring only O(n) time on average.
Our algorithm enumerates the Umplices of the Voronoi dual.It is similar to Maus' planar algorithm [12]: it employs standard bucketing techniques, and its operation in Rd corresponds to the operation of the gifiwrapping algorithm for convex hulla in Rd+'.It will be convenient to calI the d-simplices of the Voronoi dual celb and the (d -1)-simplices ficetr (since they are facets of the cells); likewise, we wiB caU an empty dsphere defined by the vertices of a cell a cell sphere and a (d -1)sphere defined by the vertices of a f&et a @et sphere.The algorithm proceeds by repeatedly finding a new cell adjacent to a known facet.Except for facets that are also facets of the (d-dimensional) convex hull, every facet belongs to exactly two cells.We maintain a dictionary of facets for which only one cell is known.At each step a facet is removed from the dictionary and its unknown cell (if it exists) is found by searching for the unknown (d+ 1)st vertex (the site search).The remaining facets of the new cell are searched for in the dictionary (facet searches).Each that is found is deleted, since both of its cells are already known.Each that is not found is inserted so that its unknown cell will be searched for in some later step.
The facet dictionary is organized as a linear array of n buckets; a random facet Us into a particular bucket with probability l/a.Within each bucket facets may be organized in a balanced search tree to insure good (logarithmic) worst-case performance, but a simple linear function Find,aiZe(f, X) -U is the unit &ball.
-SiteJuZ(z,F) is the signed distance of the center of the d-ball defined by z and 3 above the hyperplane of 3. -BozIvZ(B, 3) = mi+(qfjp~) Site,ZuZ(y, 3).-q is a priority queue of boxes ordered by BozJvl.Insert(Boz(center of facet sphere of 3), q); an8 := nil; while (q # 0) A(BotIuZ(Findmin(q), 3) < SiteJvZ(unr, 3)) do curBoz := DeZetemin(q); for zr E curBoz do if (25 E X) A (SileluZ(a, 3) < SiteJwZ(onr,3)) then mu := z; for newBoz adjoining curBoz do if ((newBoz n 3-1 I-I 2.4) # 0) A (newBoz $ q) then Insert(newBoz, q); return ans; end Find-site  zdn/pd hypercubic "boxes" of side L(h/n)'IdJ and volume about pa/n.Boxes intersecting the unit ball willcontain in expectation at most one site each.On a particular call to FindAte, let 3 and 3t be the actual arguments (a known facet and a ha&pace to search), let 2.4 be the unit &ball, and let 23 be the cell ball of the unknown cell ifit exists, If the cell exists, it is necessary and sufficient to examine those boxes intersecting (13 n 3c n 24) to determine it.If it does not exist, it is necessary and sufficient to examine the boxes intersecting (3c nu) to determine this (Figures 4, 5).The function Find-site examines almost exactly those boxes.
SiteAr(a, 3) is the center of the &sphere determined by the point z and the (d -1) vertices of the facet 3.If 3c is the halfspace (8,~) > 0, then the goal of the site search is to find the site z E 7-i for which SiteluZ(z,3) = (a,Site,ctr(z,F)) is minimized.The priority queue is primed by inserting the box containing the center of the facet sphere of 3. When a box is removed from the priority queue and its contents are examined, the 2d boxes that share a facet with it are entered into the queue for possible examination later.The boxes in the priority queue are ordered by the quantity vmw thus it guarantees that the correct site will be found, while perhaps causing a few boxes to be examined unnecessarily.The computation involves projecting the center of the box onto the hyperplane defined by 3, solving 8 quadratic equation, then choosing one of its roots.(See Figure 6, where K is the center of the box B, and dist(C, P) is the value of Boz~luZ.) The priority queue operations Insert, Find,min, and Deletemin can be implemented so that only O(logn) time is required for each, but a naive linked-list implementation in which each operation requires time proportional to the length of the list suffices for the purposes of our average-case analysis.Proof Sketch.To manage the facet dictionary, we maintain n buckets and hash facets into buckets by computing the exclusive-or of the binary-representations of the indices of the sites defining the t&et.Within each bucket we maintain a linked list of the facets iu the bucket; each search or insertion in the bucket takes time proportional to the length of this list.It is not hard to see that this scheme hashes an equal number of the d-subsets of {1,2,. . ., n) to each bucket, but we must show that those d-subsets actually defining facets are not correlated in a way that causes them to be hashed to a small number of buckets.We now turn to the search for the (d+ 1)st site completing a cell with a known facet.We call a site search "successful" if a site is found, and "unsuccessful" if no site is found because the facet lies on the boundary of the convex hull.(5.1) The use of a priority queue in Find-rite guarantees that no box is examined unless its circumsphere is contained diagram of the set x, = {WQ, . ...z.,) of n sites in Rd is the set of n convex regions Vi = { 2 1 Vj : dist(a, xi) 5 dist(2, zj) } for 1 5 i 5 n.The straight-line dual of the Voronoi diagram in the plane is called the Dehnay triungdation.In the planar case, sites zi and zj are joined by an edge in the Delaunay triangulation if and only if Y, and Yj share an edge.In d dimensions, sites zio,zil,, .e,Zir define a k-face of the Voronoi dual if and only if YiO,Yil,. ..,Yi, share a (d -k)-face in the Voronoi diagram.The Voronoi disc gram can be constructed easily from its dual and vice uersu.If, as is assumed in the sequel, no d + 2 sites fall on the same hypersphere, the following facts are easily verified.

Figure 4 :
Figure 4: What Find&e must search when successful.

Figure 5 :Figure 6 :
Figure 5: What Find-site must search when unsuccessful.

5
Analysis of the AlgorithmIn this section we show that all the facet and sites searches made by Algorithm A in the previous section can be completed in O(n) time on average.Lemma 2 The facet searches can be completed in O(n) ezpected time.
Let C(n, d) = c), let Dr,D,, . ..,Z)C(,,,~) be the dsubsets of &, and let Fi represent the condition 'Vi is a dual facet".The conditional probability Pr(Fj 1 Fi) depends only on ]Dj \ Di).If ]Dj \ Pi) = m > 0, then there are at most C(n -d,m -l)/C(m, m -1) = O(n*-l ) ways to choose the m elements of Dj.\ Vi so that 'Df and Dj fall into the same bucket.This holds because choosing the first m-1 elements fixes the mth.On the other hand, there are C(n-d, m) = O(n*) ways to choose the m elements of ZJj \ Vi overall.Thus only a fraction 0(1/n) of the facets fall into the same bucket as Dd.Roughly, this shows that any correlation is at the level of constant firctors only. q

Lemqa 3
The successfit site searches can be completed in O(n) ezpected time.Proof Sketch.If we define the distance between a point x and a set Y by dist(z, Y) = mi; dist(z, y), all boxes with circumsphcres intersecting B n X nZ.4 are completely contained by the set A = (z 1 dist(z, B flu) < ,&',&)"d}.