Voronoi Diagram and Delaunay Triangulation in O(n log n) with Fortune's Algorithm

Revision en4, by Monogon, 2020-12-17 01:32:36

What's a Voronoi Diagram?

Given a set $$$S$$$ of $$$n$$$ points in the 2D plane, the Voronoi diagram is a partition of the plane into regions. The region associated with a point $$$p\in S$$$ contains precisely the points $$$q$$$ such that $$$q$$$ is closer to $$$p$$$ than any other point in $$$S$$$. In other words, a point $$$q$$$ belongs to the region of its nearest neighbor.

The Voronoi diagram can be represented by the planar graph of boundaries between regions. A vertex in this graph is where three or more segments meet (or a point at infinity), and an edge is a segment connecting two of these vertices. The way the regions were defined, we see that a segment contains the points equidistant to two points in $$$S$$$, so it is part of their perpendicular bisector. Similarly, a Voronoi vertex is equidistant to three or more vertices in $$$S$$$, so it is their circumcenter. Recall that the circumcenter is the intersection of perpendicular bisectors.

What's a Delaunay Triangulation?

The Delaunay triangulation is the dual graph of the Voronoi diagram. That is, a Voronoi vertex is a Delaunay face and a Delaunay vertex is a Voronoi vertex. A Delaunay edge connects two points if and only if their corresponding Voronoi regions share a border.

In the following image, the set $$$S$$$ corresponds to the black points. The red vertices and edges are the Voronoi diagram. The black points and edges are the Delaunay triangulation.

The Delaunay triangulation only involves edges between existing points, while the Voronoi diagram creates new vertices and needs a way to represent edges going infinitely in one direction. For this reason, it is often more convenient to compute the Delaunay triangulation and only convert to the Voronoi diagram if needed.

Why are these useful?

The definition of the Voronoi diagram immediately shows signs of applications. Given a set $$$S$$$ of $$$n$$$ points and $$$m$$$ query points $$$p_1,\ldots, p_m$$$, we can answer for each query point, its nearest neighbor in $$$S$$$. This can be done in $$$O((n + q)\log(n+q))$$$ offline by sweeping the Voronoi diagram and query points. Or it can be done online with persistent data structures.

Other than this, the Delaunay triangulation has many amazing properties of its own that can make it extremely useful for solving problems.

  • For each Delaunay triangle, its circumcircle does not strictly contain any points in $$$S$$$. (In fact, you can also consider this the defining property of Delaunay triangulation)
  • The number of Delaunay edges is at most $$$3n-6$$$, so there is hope for an efficient construction.
  • Each point $$$p\in S$$$ is adjacent to its nearest neighbor with a Delaunay edge.
  • The Delaunay triangulation maximizes the minimum angle in the triangles among all possible triangulations.
  • The Euclidean minimum spanning tree is a subset of Delaunay edges.

Degenerate Cases

Because this is geometry, there are some degenerate cases I have hidden from you up to this point. The first case is when all the points in $$$S$$$ are contained in a line. In this case, the Voronoi diagram has all edges connecting two points at infinity, and the Voronoi diagram is disconnected (in all other cases it's connected). Also in this case, the dual graph is no longer a triangulation. We may still decide to represent it by the path of edges connecting the points.

Another nasty case is co-circular points. Here, there are Voronoi vertices where more than three segments intersect. In this case, the dual graph has a face with more than three sides. However, we want it to be a triangulation, so we can also say that any triangulation of such faces is valid, and we will allow non-unique Delaunay triangulations.

Algorithms

Many algorithms exist to compute the Delaunay triangulation and/or Voronoi diagram of a set of points. Each have their advantages and disadvantages.

Flip Algorithms

In this approach, you start with an arbitrary triangulation of the points. Then you repeatedly fix adjacent triangle pairs until the triangulation has the Delaunay property. The main advantage of this method is its simplicity: an arbitrary triangulation can be constructed with an incremental convex hull algorithm, and the disadvantage is that it could require $$$\Omega(n^2)$$$ flips in the worst case.

Incremental Algorithms

In this approach, you always maintain the Delaunay triangulation of a subset of the points and insert the points one at a time. When a point is inserted, you should delete all triangles whose circumcircle contains the newly added point. After deleting this connected set of triangles, you should repair the face it exposes. For each edge in this exposed face, you should add a triangle with that edge and the new point. This gives the straightforward $$$O(n^2)$$$ Bowyer-Watson algorithm.

If we insert the points in a random order, we can significantly reduce the number of triangle insertions and deletions. But to get expected $$$O(n\log n)$$$ runtime, we would need to avoid processing all the unaffected triangles. And this makes the implementation far more challenging. We need to maintain a search structure of the triangles in order to locate which one contains the new point.

The advantage of this method is that it achieves expected $$$O(n\log n)$$$ runtime with only integer computations, and the construction is online. The disadvantage is that it has a bad constant, isn't deterministic, and the triangle search structure can make for a tricky implementation.

Divide and Conquer

The divide and conquer algorithm first sorts the points by $$$x$$$-coordinate. Then it splits the points into two sets, recursively computes their Delaunay triangulations, and merges them together. As long as you can perform the merge in $$$O(n)$$$ time, the overall algorithm is $$$O(n\log n)$$$. When merging the triangulations, it is actually quite hard to keep track of the important adjacency information. The most popular strategy is perhaps the Quad-Edge data structure.

The advantage is that this gives a deterministic $$$O(n\log n)$$$ algorithm with only integer computations. The disadvantage is the complexity of the quad-edge data structure and the costly memory usage that comes with it.

3D Convex Hull Reduction

An interesting fact is that the problem of Delaunay triangulation can be reduced to 3D convex hull. If for each point $$$(x,y)$$$ we also give it the $$$z$$$-coordinate $$$x^2+y^2$$$ and compute the 3D convex hull, then the downward-facing convex hull faces correspond exactly to the Delaunay triangles. The advantage is that if you happen to have a good 3D convex hull implementation, you get Delaunay triangulation for free. The disadvantages are that $$$O(n\log n)$$$ 3D convex hull is a whole beast of its own, and the large $$$z$$$-coordinates will exacerbate precision and overflow errors in the 3D convex hull implementation. The main strategies for 3D convex hull are incremental and divide-and-conquer algorithms, which are only more complicated than their counterparts in the special case of Delaunay triangulation.

My previous blog on 3D convex hull gave an $$$O(n\log n)$$$ randomized incremental implementation, but it turns out to have quite a bad constant for practical use.

Fortune's Algorithm

In this blog, I will focus on Fortune's Algorithm. It takes a sweep-line approach to the problem, and achieves a deterministic $$$O(n\log n)$$$ time complexity. The main disadvantage is that computations are done with floating points instead of integers, opening up the door to precision errors. However, keep in mind that the algorithms using only integers have a risk of integer overflow even on reasonably small constraints.

The main challenge in designing a sweep-line algorithm for Voronoi diagrams is that a point after the sweep-line can affect the diagram before the sweep-line. That is, if the sweep-line denotes the discovery time of a point, we cannot know the entire Voronoi diagram before the sweep-line. The way Fortune's algorithm overcomes this challenge is by maintaining the Voronoi diagram only for the area we can be certain about, regardless of what points lie after the sweep-line.

Now, which points can belong to the region of a point to the right of the sweep-line? Well, such a point must be closer to the sweep-line than to any point in $$$S$$$ to the left of the sweep-line. Recall that the set of points equidistant to a point and a line is given by a parabola, where the point is called the focus and the line is called the directrix. If the focus is a point $$$p\in S$$$ to the left of the sweep-line and the directrix is the sweep-line, then everything to the left of this parabola we can be certain of its region. Now, if we make such a parabola for every point $$$p\in S$$$ to the left of the sweep-line, everything to the right of all these parabolas is precisely what might belong to the region of a point to the right of the sweep-line. The boundary of the certainty and uncertainty regions is called the beach line, consisting of parabolic arcs.

Now we are ready to understand Fortune's algorithm abstractly. We scan a sweep-line from left to right, maintaining the combinatorial structure of the beach line. The parabolas are only for visual aid, and the ordered list of foci is what is actually stored. As points get added and parabolic arcs shrink to length $$$0$$$, these events tell us precisely the structure of the Voronoi diagram and Delaunay triangulation.

There are two types of events which change the beach line structure: point events and vertex events.

Point Events

In a point event, the sweep-line discovers a point $$$p\in S$$$ and inserts its parabolic arc into the beach line. In order to do this efficiently, we need to search for where it belongs in the beach line. Then it will split the arc it sees into two. An example of a point event is shown below (here the sweep-line goes from top to bottom instead).

Vertex Events

In a vertex event, a parabolic arc disappears from the beach line, corresponding to a Voronoi vertex. To handle a vertex event, we simply remove the point whose arc length shrinks to length $$$0$$$. Such an event is generated dynamically by a triple of points adjacent on the beach line. Every time we modify the beach line, we update the relevant points for their vertex events. To compute the time a vertex event should be handled, we take the circumcircle of the three points generating the event. The rightmost point of the circumcircle is the time we should process this event.

Summary

In summary, the algorithm maintains a priority queue of events. Initially, we push all point events into the priority queue. When processing a point event, we find the correct arc in the beach line, split it up, add the new point, and add new vertex events if needed. When processing a vertex event, we add the generating triple as a Delaunay triangle, remove the arc that should disappear, and add new vertex events if needed.

Fortune's Algorithm Implementation

My code doesn't work when two points have the same $$$x$$$ coordinate. This is handled simply by rotating all input points by $$$1$$$ radian and praying to the geometry gods.

The code is supposed to be in this spoiler but don't blame me: https://codeforces.net/blog/entry/85432

Example Problems

Tags geometry, voronoi diagrams, delaunay triangulation, fortunes algorithm

History

 
 
 
 
Revisions
 
 
  Rev. Lang. By When Δ Comment
en7 English Monogon 2020-12-17 06:57:06 8 Tiny change: 'a Voronoi vertex. A Delaun' -> 'a Voronoi region. A Delaun'
en6 English Monogon 2020-12-17 03:02:57 11
en5 English Monogon 2020-12-17 02:09:29 1790 Tiny change: '</spoiler>' -> '</spoiler>\n\nReferences\n------------------' (published)
en4 English Monogon 2020-12-17 01:32:36 7082 Tiny change: 'um=1520)\n\n<spoiler summary="Caravans Solution">\nugh\n</spoiler>' -> 'um=1520)\n'
en3 English Monogon 2020-12-17 01:05:21 3951 Tiny change: 'hm-slowed.gif"/>\n\nFor' -> 'hm-slowed.png"/>\n\nFor'
en2 English Monogon 2020-12-16 22:44:18 4091 Tiny change: 'rithms\n - asdf\n\nF' -> 'rithms\n * asdf\n\nF'
en1 English Monogon 2020-12-16 20:13:46 4155 Initial revision (saved to drafts)