Voronoi Diagrams

An ordinary Voronoi diagram is formed by a set of points in the plane called the generators or generating points. Every point in the plane is identified with the generator which is closest to it by some metric. The common choice is to use the Euclidean $L_2$ distance metric

$$\vert\mathbf{x}_1 - \mathbf{x}_2\vert = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}$$

where $\mathbf{x}_1=(x_1,y_1)$ and $\mathbf{x}_2=(x_2,y_2)$ are any two points in the plane. The set of points in the plane identified with a particular generator form that generator's Voronoi region, and the set of Voronoi regions covers the entire plane. Figure 1 illustrates a set of generating points and their associated Voronoi regions.

Figure 1: Voronoi diagram with generators (large dots) and centroids (small dots)

We implemented the the fast 3D graphics hardware-based algorithm in Hoff [5] and originally in [11] to compute our Voronoi diagrams. The algorithm draws a set of right cones with their apexes at each generator. The cones all have the same height and are viewed from above the apexes with an orthogonal projection. In addition, each cone is given a unique colour which acts as the generator's identity. Since the cones must intersect if there is more than one generator, the z-buffer determines for each pixel which cone is closer to the viewer and assigns that pixel the appropriate colour value. We can then scan the resulting image and determine which generator is closest to each pixel by using the unique colours. This technique allows us to compute discrete Voronoi diagrams extremely quickly and perform computations on the resulting regions.

Centroidal Voronoi Diagrams

A centroidal Voronoi diagram has the odd property that each generating point lies exactly on the centroid of their Voronoi region. The centroid of a region is defined as

$$\mathbf{C}_i = \frac{\int_A{\mathbf{x}\rho(\mathbf{x}){\textstyle dA}}} {\int_A{\rho(\mathbf{x}){\textstyle dA}}}$$ (1)

where $A$ is the region, $\mathbf{x}$ is the position and $\rho(\mathbf{x})$ is the density function. For a region of constant density $\rho$, the centroid can be considered as the centre of mass. Figure 1 has the centroids of each region marked with small circles.

A centroidal Voronoi diagram is a minimum-energy configuration in the sense that it minimises $\int_A{\rho(\mathbf{x})\vert\mathbf{C}_i - \mathbf{x}\vert^2}$ [2]. Practically speaking, a centroidal distribution of points is useful because the points are well-spaced in a definite sense. Figure 2 shows a centroidal Voronoi diagram.

Figure 2: Centroidal Voronoi diagram with generators

Generating Centroidal Voronoi Diagrams

Lloyd's method [8] is an iterative algorithm to generate a centroidal Voronoi diagram from any set of generating points. The algorithm can simply be stated:

Figure 1 relaxes under Lloyd's algorithm to become Figure 2. The convergence of Lloyd's algorithm to a centroidal Voronoi diagram has been proven for the one-dimensional case and higher dimensions seem to act similarly [2]. There are several different convergence criteria which should be equivalent in the limiting case as the algorithm runs forever. The obvious criteria to use would be that the computed centroids are numerically equal to the generating points. However, for most applications this criteria is far too stringent and it would be perhaps better to look at the average distance moved by all generating points. Since we are interested in generating well-spaced sets of points, we look at the average change in inter-point distance, or equivalently, the average change in Voronoi region area.

Efficient Computation of Centroids

Calculating the centroids requires efficiently evaluating the integrals in equation (1). Since the integrals are over arbitrary Voronoi regions, we convert to iterated integrals and integrate the region row by row. In this manner we can precompute much of the integral.

The denominator of the centroid is transformed as follows:

$$\begin{eqnarray*} \int_A{\rho(\mathbf{x}){\textstyle dA}} & = & \int_{y_1}^{y_2... ...xdy} \\ & = & \int_{y_1}^{y_2}\left[P\right]_{x_1}^{x_2}dy \\ \end{eqnarray*}$$

where $P\equiv P(x,y)\equiv\int_0^x\rho(s,y)ds$ can be precomputed from the density function². Note that we cannot precompute the entire integral because we do not know the boundaries of the Voronoi regions beforehand.

The numerator of the y-coordinate of the centroid is transformed similarly:

$$\begin{eqnarray*} \int_A{y\rho(x,y){\textstyle dA}} & = & \int_{y_1}^{y_2}\int_... ...dy} \\ & = & \int_{y_1}^{y_2}y\left[P\right]_{x_1}^{x_2}dy \\ \end{eqnarray*}$$

The numerator of the x-coordinate of the centroid involves integration by parts:

$$\begin{eqnarray*} \int_A{x\rho(x,y){\textstyle dA}} & = & \int_{y_1}^{y_2}\int_... ...\\ & = & \int_{y_1}^{y_2}\left[xP - Q\right]_{x_1}^{x_2}dy \\ \end{eqnarray*}$$

where $Q\equiv Q(x,y)\equiv\int_0^xP(s,y)ds$ can also be precomputed from the density function.

Note that the final expressions require numerical integration only in the y-direction and otherwise involve expressions only at the region boundaries $x_1$ and $x_2$. $P$ and $Q$ are precomputed once from the density function and then evaluated at the horizontal endpoints $x_1$ and $x_2$ as needed. This allows us to compute the integrands only at region boundaries and not at every pixel. Otherwise we would have to compute the integrands $x\rho$ and $y\rho$ for every span of pixels across a region and numerically integrated. The above integrand computation is particularly simple - at worst two lookups for $P$ and $Q$, a multiplication and a subtraction.

Resolution of Voronoi Calculation

One disadvantage of using a discrete calculation of the Voronoi regions is the centroids calculation is effected by the resolution of the diagram. The relative error of the calculated centroid location will increase as the number of pixels in the Voronoi region decreases. A related problem is that if the resolution is low enough, two generating points can effectively overlap and one of the regions will disappear. The solution, as described in [5], is to split the diagram into tiles and compute each tile at the full resolution available and then stitch the full diagram back together at a higher virtual resolution. The virtual resolution can be increased arbitrarily to meet a lower bound on Voronoi region pixel area.

Footnotes

... function²

Recall that $\left[\int_0^xf(s)ds\right]_a^b = \int_0^bf(s)ds - \int_0^af(s)ds = \int_a^bf(s)ds$