Voronoi To The World
Voronoi diagrams arise naturally in many spatial analyses. More importantly for this author they also tend to be aesthetically attractive, making them useful for crafting compelling illustrations.
Efficiently constructing a Voronoi diagram, however, is not particularly straightforward. In this note we’ll work through the geometric structure of Voronoi diagrams and an algorithm that leverages that structure to implement Voronoi diagrams with the maximum possible performance.
1 Tessellating A Plane
A Voronoi diagram, also known as a Dirichlet tessellation, is an organization of a plane into subsets defined by the nearest proximity to given collection of points.
From what I can tell the standard reference for Voronoi diagrams, and other topics in computational geometry, is Berg et al. (2008). This text is sometimes referred to as “the 4Ms” which I can only infer refers to three out of the four authors having the first name “Mark” with the forth, non-Mark author being an honorary “M”. I humbly suggest that “Mark, Marky, Mark, and the Funky Bunch” is more appropriate, but to each their own.
To define a Voronoi diagram formally consider a two-dimensional Euclidean plane \mathbb{R}^{2} and a collection of sites ( \, s_{1} = (x_{1}, y_{1}), \ldots, s_{N} = (x_{N}, y_{N}) \, ) \subset \mathbb{R}^{2}. Give a particular site s_{n} we can construct a subset of points, also known as a cell, that are closer to that site than any other site, \mathsf{c}_{n} = \{ p \in \mathbb{R}^{2} \mid d( p, s_{n} ) < d( p, s_{n'} ) \; \mathrm{for all} \; n' \ne n \}, where d(p, p') is a distance function. The Voronoi diagram of a collections of sites is the corresponding collection of nearest-proximity subsets (Figure 1).
The standard Voronoi diagram considers the Euclidean distance function, d( p, p' ) = \sqrt{ (x - x')^{2} + (y - y')^{2} }. This construction, however, can be generalized to other distance functions. At the same time Voronoi diagrams can be defined on higher-dimensional real spaces and even more general metric spaces.
A Voronoi diagram almost forms partition of a Euclidean plane. The obstruction is that points that are equally distant from two or more sites don’t fall into any of the cells, but rather form a singular boundary between the cells. The combination of the Voronoi cells and the boundary between them, however, does form a proper partition.
Because each cell contains one, and only one, site, we can always label the cells by the associated site (Figure 2 (a)). Similarly we can label the linear segment that forms the boundary between two neighboring cells by the sites associated with those two cells, and the punctual boundaries separating the corners of three neighboring cells by the sites associated with those cells (Figure 2 (b), Figure 2 (b)).
2 The Geometry of Voronoi Diagrams
The boundary between the cells of a Voronoi diagram is useful in its own right. In particular it defines an undirected graph, with linear edges separating neighboring pairs of cells and punctual vertices separating neighboring triplets of cells.
Often it is easier in practice to construct this Voronoi graph and then derive the Voronoi cells from it rather than trying to construct the Voronoi cells directly. The vertices and edges of the Voronoi graph exhibit rich geometric structure that provides the foundation for powerful algorithms.
2.1 Voronoi Vertices
The vertices of the Voronoi graph are defined by points that are the same distance from three neighboring Voronoi cells (Figure 3). Consequently we can uniquely label each vertex with a triplet the sites contained by those cells, (s_{n_{1}}, s_{n_{2}}, w_{n_{3}} ). That said, not every triplet of Voronoi cells are neighbors and not every triplet of cells defines a valid vertex. Fortunately Euclidean geometry provides tools for identifying all of the valid vertices.
Each triplet of sites falls onto a unique circumcircle with center \begin{align*} x_{c} &= \frac{1}{2} \frac{ \left( x_{n_{1}}^{2} + y_{n_{1}}^{2} \right) \, (y_{n_{2}} - y_{n_{3}}) + \left( x_{n_{2}}^{2} + y_{n_{2}}^{2} \right) \, (y_{n_{3}} - y_{n_{1}}) + \left( x_{n_{3}}^{2} + y_{n_{3}}^{2} \right) \, (y_{n_{1}} - y_{n_{2}}) }{ x_{n_{1}} \, (y_{n_{2}} - y_{n_{3}} ) + x_{n_{2}} \, (y_{n_{3}} - y_{n_{1}} ) + x_{n_{3}} \, (y_{n_{1}} - y_{n_{2}} ) } \\ y_{c} &= \frac{1}{2} \frac{ \left( x_{n_{1}}^{2} + y_{n_{1}}^{2} \right) \, (x_{n_{3}} - x_{n_{2}}) + \left( x_{n_{2}}^{2} + y_{n_{2}}^{2} \right) \, (x_{n_{1}} - x_{n_{3}}) + \left( x_{n_{3}}^{2} + y_{n_{3}}^{2} \right) \, (x_{n_{2}} - x_{n_{1}}) }{ x_{n_{1}} \, (y_{n_{2}} - y_{n_{3}} ) + x_{n_{2}} \, (y_{n_{3}} - y_{n_{1}} ) + x_{n_{3}} \, (y_{n_{1}} - y_{n_{2}} ) } \end{align*} and radius \begin{align*} r &= \sqrt{ (x_{c} - x_{n_{1}})^{2} + (y_{c} - y_{n_{1}})^{2} } \\ &= \sqrt{ (x_{c} - x_{n_{2}})^{2} + (y_{c} - y_{n_{2}})^{2} } \\ &= \sqrt{ (x_{c} - x_{n_{3}})^{2} + (y_{c} - y_{n_{3}})^{2} }. \end{align*} By construction the center of the circumcircle is equidistant from all three sites, and hence defines a potential vertex in the Voronoi graph.
The three sites defining the circumcircle are the three closest sites to this potential vertex if and only if no other sites are inside of the circumcircle. In other words Voronoi vertices are defined by empty circumcircles.
2.2 Voronoi Edges
The edges in a Voronoi graph are comprised of points equidistant to two neighboring cells (Figure 5), allowing us to uniquely label each Voronoi edge with the pair of sites contained by those cells. Once, again, however, we have to take care because arbitrary pairs of sites do not always define valid edges.
Before worrying out isolating proper edges let’s first work out how to identify the points that define potential edges. Any pair of sites defines a line that connects them, and a mid point on that line that is equally distant to the sites. The perpendicular bisector of two sites is the line perpendicular to this connecting line that intersects with the midpoint (Figure 6 (a)). Every point along the perpendicular bisector between two sites is equally distant from those two sites and, consequently, could contribute to a Voronoi edge (Figure 6 (b)).
In order to determine which segment, if any, of a perpendicular bisector forms a Voronoi edge we come back to circles. Any point along a perpendicular bisector is the center of a unique circle that includes the two defining sites. If no other sites are on or inside of this circle then those two sites are the nearest neighbors to each other and that point is part of a Voronoi edge.
As we move along the perpendicular bisector these circles will eventually intersect with another site. When they do the point along the bisector is equally distant to three sites – the two initial sites and this new site – but the interior of the circle will remain empty. In other words the Voronoi edge terminates as soon as we hit a Voronoi vertex.
2.3 Voronoi Faces
The edges of a Voronoi graph trace around each site, defining faceted faces in the graph that correspond to the Voronoi cells (Figure 8). When working over the entire, unbounded, Euclidean plane not all of these faces will be closed. While some of the sites will be entirely fully enclosed by a path of edges, some will be bounded on only one side by a path of edges. I will refer to these as closed faces and open faces, respectively.
3 Fortune’s Algorithm
There are here are multiple ways to construct a Voronoi graph, most of them suffer from wasteful computation. For large enough numbers of sites these approaches quickly become too expensive to be practical. Fortunately there is one elegant approach that, while not at all initially obvious, is able to construct a Voronoi graph as efficiently as possible.
3.1 Inefficient Approaches
Before diving into efficient solutions let’s consider some of the more straightforward but ultimately inefficient approaches.
One way to identify all of the vertices in a Voronoi graph is to loop over each triplet of sites and then construct the corresponding circumcircles. For each of these circumcircles we then loop over the remaining sites to check for inclusion, with the centers of the empty circumcircles defining vertices.
We could then construct candidate Voronoi edges by looping over each pair of sites and scanning across the points between the two Voronoi vertices associated with those cites. Next we could take any point along each candidate edge, construct a circle centered on that point and intersecting with the two associated sites, and check if any other sites are inside. The candidates with empty circles define Voronoi edges.
Because there are { N \choose 3 } = \frac{ N \, (N - 1) \, (N - 2) }{ 6 } triplets of sites and { N \choose 2 } = \frac{ N \, (N - 1) }{ 2 } pairs of sites to consider, however, the number of operations we need to implement this brute-force approach scales cubically with the number of sites.
That said, much of this cost is wasted on the repeated querying of each site; we should be able to construct Voronoi graphs more efficiently if we can query each site fewer times. One more efficient approach is to loop over each site, construct the half-planes defined by the perpendicular bisector between the active site and each other site, and then construct the corresponding Voronoi face from the intersection of those half-planes. This cost of this approach requires only quadratic number of operations.
We can, however, do even better. By systematically scanning across a plane Fortune’s algorithm Berg et al. (2008) is able construct Voronoi graph with a number of operations proportional to N \log N. This performance turns out to be theoretically optimal. The main downside of Fortune’s algorithm is that is requires some care to understand and then implement in a way that achieves that optimal scaling.
3.2 The Sweep Line
Fortune’s algorithm processes the sites not by working through a list but rather by systematically scanning across the ambient plane with the use of a sweep line. Sites on one side of the sweep line are active and can contribute to the construction of the Voronoi graph while sites on the other side are inactive and cannot yet contribute. As the sweep line progresses more sites are processed and more of the Voronoi graph is built up.
In theory a sweep line can progress in any direction, but most references consider vertical lines that scan horizontally from left to right, horizontal lines that scan vertically from top to bottom, or horizontal lines that scan vertically from bottom to top. These axis-aligned orientations of the sweep line make geometric calculations substantially more straightforward.
Here I will consider a vertical sweep line that scans horizontally from left to right. I will denote the horizontal position of the scan at any time by S. I will also refer to any site to the left of the sweep line, x_{n} < S as active.
3.3 The Beach Line
A sweep line partitions the ambient plane into two pieces, allowing the construction of the Voronoi graph to focus on just a subset of the entire plane. That said, not all of the Voronoi graph to the left of the sweep line is completely determined by the active sites; some features of the Voronoi diagram can be influenced by sites just beyond the sweep line.
Consider, for example, a single active site s_{n} = (x_{n}, y_{n}) and the sweep line position S > x_{n}. If a point (x, y) on the left of the sweep line is closer to s_{n} than the sweep line, (x - x_{n})^{2} - (y - y_{n})^{2} < (S - x)^{2}, then we it will be closer to s_{n} than any other site that might be hiding beyond the sweep line. On the other any point to the left of the sweep line that is closer to the sweep line than it is to s_{n}, (x - x_{n})^{2} - (y - y_{n})^{2} > (S - x)^{2}, could be be closer to an inactive site on the other side of the sweep line.
When determining geometric features such as the nearest sites we cannot analyze the entire region to the left of the sweep line. Instead we can analyze only the subset of points that are closer to s_{n} than they are to the sweep line. The boundary of this analyzable subset is defined by the points equidistant to s_{n} and the sweep line, which traces out a rotated parabola (Figure 9) x = f(y) = \frac{1}{2} \left( S + x_{n} - \frac{( y - y_{n} )^{2}}{S - x_{n}} \right).
When there are multiple active sites then the region that we can safely analyze the points that are closer to any of the active sites than the sweep line. The boundary of this region is defined by the segmented envelope of the parabolas between each active site and the sweep line (Figure 10). Given its resemblance to the shape made by waves washing up on a beach (Figure 11) this boundary is known as the beach line.
