I went through a period where I was continuously dreaming about triangles. Whenever I have a hard problem I find myself mulling it over at night, but around 2003 I was working on navigation meshes and seeing triangles in my sleep every night. I have since broken away from computational geometry and love my new work in lighting, but I would like to reminisce about Delaunay triangulation for old times’ sake.
Delaunay triangulations are more of a natural discovery than a mathematical construction. They are a ‘rich’ structure and giving a single definition seems to under sell them, but we must start somewhere:
The Delaunay triangulation of a set of given points in the plane is a set of triangles where the circumcircles defined by every triangle contain no points in their interior.
Such a triangulation always exists but is only unique if no more than 3 vertices exist on each circumcircle. For example, if 4 points form a square then all 4 points are on the same circumcircle and you can triangulate the square 2 ways. By our definition, both options would be valid Delaunay triangulations. The same principle applies to the vertices of a pentagon, hexagon and so on – they are “co-circular” points. In practice, triangulation algorithms simply impose some form of tie-breaker to ensure they output triangles. This arbitrary decision is a bit messy, so I prefer to think of there being a more natural structure you might call Delaunay Polygonisation which would leave co-circular regions as polygons.
Delaunay triangulations have several nice properties, the main one being that they produce the most maximally convex triangles. Avoiding thin long triangles is a universally good thing, although it should be noted that the Delaunay property doesn’t prevent thin triangles but it will avoid them where possible. From a practical viewpoint, Delaunay triangulations are a good basis for interpolation and helpful in most mesh structures while being fairly inexpensive to compute. But in my eyes, their beauty comes from the surprising relations that connection them to other structures.
A Voronoi diagram is the closest structure and the dual of a Delaunay triangulation. By “dual” we mean that for every Delaunay polygon we can find a corresponding vertex in our Voronoi diagram, and likewise for every Voronoi cell we can find a vertex in our Delaunay triangulation. This connection is so close you can think of them as different representations of the same information. A Voronoi diagram encodes nearest neighbour information explicitly and is composed of convex “cells” rather than triangles. (Cells are like polygons but can be open). The duality between Voronoi and Delaunay allows us to get much of the useful information we may like in a Voronoi diagram, but while working with a triangulation and avoiding the fiddly line-line intersection accuracy issues that can arise in constructing Voronoi diagrams in practice.
Voronoi diagrams have a somewhat simpler and more natural definition than Delaunay triangulation. The cells in a Voronoi diagram describe the region of the plane that is closest to a given point. The edges are therefore the boundaries in the plane where we are equidistant from two points. Vertices are locations where we are equidistant from more than 2 points. If we have a Voronoi diagram and a point of interest, by locating the cell that contains the point we can immediately state our closest point and our closest neighbours.
Before describing an algorithm to compute these diagrams, I’d like to show you how to draw them by hand.
Grab an envelope and try placing some random dots on a page, then draw a convex hull around them. Another remarkable property of Delaunay triangulations is that their boundary is always the convex hull of the set of points. So if you compute the Delaunay triangulation you get the convex hull ‘for free’.
Random points with convex boundary
You might now be able to see triples of vertices, which if you drew a circle through them would not enclose any other points. If you find one, draw in the triangle. There is always only 1 circle you can draw through 3 points but it can be hard to draw free hand for large circles or where more than 3 points are close to being co-circular. However you can improve your accuracy by first finding the centre of the circle. If you take two of the points and draw a separating line half way between them your circumcircle must have a centre on this line. With a bit of imagination and/or scribbling you should be able to see that for every point on this separating line you can draw a circle centred on the line that passes through your two points. The third point then constrains the choice down to a single circle.
Points with a circumcircle and separating line
If this approach becomes too inaccurate you can build on the separating lines to construct the Voronoi diagram locally, and use the duality to find the Delaunay triangulation. As you might expect, the separating lines are the basis for edges in the Voronoi diagram. If you draw these relatively accurately you can inspect where they intersect. At intersections the nearest neighbours change, and you can manually erase segments that are no longer relevant, giving you the Voronoi diagram. The duality tells us that for every edge in this Voronoi diagram we must have an edge in the Delaunay triangulation – even if the edge in the Voronoi diagram is not directly between the two points that created it, which is sometimes the case. So for every edge you find in your Voronoi diagram, join the two points that created it.
Delaunay triangulation with some Voronoi edges
When you have co-circular points you will see three or more separating lines intersecting at the same location. If you only joined points based on edges in the Voronoi diagram this would leave you with polygons around these intersection points. So vertices with degrees greater than 3 in the Voronoi diagram correspond to polygons in the Delaunay triangulation.
By far the most popular algorithm for computing Delaunay triangulations is Steve Fortune’s sweep line algorithm, which is O(n log n) in the number of points. This running time comes from the fact it must sort the points. The algorithm is conceptually elegant and worthy of its own article, although sadly robust implementations are necessarily full of obscure conditions for edge cases and far from easy to read. I actually find the disparity between the elegance of the theory of the sweep algorithm and its real world incarnation interesting. It appears we do not yet have a way to express the edge cases without enumerating them by hand. I would hope that some more general abstraction exists that could generate the implementation.
Conveniently, there is a simpler O(n^2) algorithm that exploits another fascinating property of these structures: The Delaunay triangulation of a set of 2D points can be found by computing the 3D convex hull of the same points with a generated z coordinate.
Intuitively, what links these structures is the notion of ‘nearest’. The shrink-wrapping effect produced by a convex hull turns out to be intimately related to the nearest neighbour property of the Voronoi diagram and similarly, the maximally convex property of the Delaunay triangulation. It is this tight inter-relationship and natural construction that makes structures them so interesting.
Given our points in the plane, we shall add a new z coordinate based on their squared distance from the origin. The location of the origin itself is unimportant, but it is common to locate it in the centre of the set of points where we can most clearly see that the points now form a bowl. Perhaps surprisingly, all the downward faces in the convex hull of this new set of points is a triangle (or polygon) in the Delaunay triangulation. So once we have the convex hull, we can simply read off the Delaunay triangulation.
Relation of Delaunay triangulations to convex hulls