Delaunay triangulation and Voronoi diagram
Also known as: Delaunay triangulation · Voronoi diagram · Dirichlet tessellation · Thiessen polygons · Bowyer–Watson algorithm · Lawson flips · empty circle · largest empty circle
The zone of each fire station in Almaty, the points for which it is the nearest, is a Voronoi diagram, and neighbours by zone join into the Delaunay triangulation, where no point lies inside the circumcircle of another triangle. How they are built, why empty circles give the “roundest” triangles, and what the paraboloid has to do with it.
OpenStreetMap marks 17 fire stations inside the boundary of Almaty. Divide the city into zones: each station gets the points to which it is closer than all others in a straight line. This gives a Voronoi diagram, polygons fitted together without gaps. Join the stations whose zones touch, and you get the Delaunay triangulation: the circumcircle of any of its triangles contains no other station. It is one structure seen from two sides, and it answers questions that would otherwise need brute force: which station is nearest, where the nearest one is furthest away, how to connect everyone most cheaply.
The point of the city furthest from the stations is on the southern boundary, in the mountains: 16.7 km to the nearest station. The question was taken literally, and the answer is literal: the city boundary takes in slopes where nobody lives. Ask differently, which of the 7909 street crossings is furthest from the stations, and you get 6.5 km. The algorithm is the same in both cases; the question differs.
The idea
Let points be given on the plane, called sites. The Voronoi cell of a site is all the points of the plane for which it is closer than the others. Points equidistant from two sites lie on the perpendicular bisector between them, and closer to site p than to q is a half-plane on one side of it. So a cell is an intersection of half-planes, a convex polygon: this is how the scene above builds the 17 zones. The boundary between two cells is a piece of a perpendicular bisector, and a vertex where three cells meet is a point equidistant from three sites, that is, the centre of the circle through them.
Join two sites with an edge if their cells share a side. This gives triangles, each with a Voronoi vertex at the centre of its circumcircle. That circle cannot contain another site: its centre is equally far from its three sites and closer to them than to any other, otherwise it would lie in another cell. This is the definition seen from the other side.
The Delaunay triangulation of a set of points is a triangulation in which the circumcircle of every triangle is empty: there is no point of the set inside it. If no four points lie on one circle, such a triangulation exists and is unique (this follows from the lift onto the paraboloid below), and its edges are exactly the pairs of sites whose Voronoi cells share a side. If four points are on one circle, like the corners of a rectangle, either diagonal can be drawn: there is one Voronoi diagram but several Delaunay triangulations.
Any triangulation of n points, h of which lie on the boundary of the convex hull (vertices and points on its sides), has exactly 2n − 2 − h triangles and 3n − 3 − h edges.
A triangulation is a planar graph: n vertices, E edges and F faces, of which T = F − 1 are triangles and one is the outer face. By Euler's formula n − E + F = 2. Count the sides of the faces: each triangle has three, the outer face h, and each edge is a side of exactly two faces. So 3T + h = 2E. Substitute F = T + 1 into Euler's formula: E = n + T − 1. Then 3T + h = 2n + 2T − 2, whence T = 2n − 2 − h and E = 3n − 3 − h.
For 17 stations, 5 of them on the hull, that is 27 triangles and 43 edges, for any triangulation, not only Delaunay's. An important consequence: both the triangulation and the Voronoi diagram have linear size, fewer than 3n edges, although there are n(n − 1)/2 pairs of points. Everything one needs to know about the closeness of points fits in O(n) numbers.
What lies inside
An edge pq belongs to the Delaunay triangulation if a circle can be drawn through p and q with no points of the set inside: if you inflate it while keeping it through p and q, it will hit a third point and become the circumcircle of a Delaunay triangle. Two useful facts follow at once from this test.
The nearest neighbour of any point is its neighbour in the Delaunay triangulation. The minimum spanning tree of points of the plane consists entirely of edges of the Delaunay triangulation.
Let q be the point nearest to p. The circle with diameter pq lies within the disc of radius |pq| centred at p, and there are no points inside that disc, otherwise there would be a point closer than q. So the circle on the diameter pq is empty, and pq is a Delaunay edge.
Let the edge pq belong to the minimum spanning tree, but let there be a point r inside the circle with diameter pq. Then the angle prq is obtuse, and in the triangle prq the side pq is the longest: |pr| and |rq| are less than |pq|. Remove pq from the tree: it falls into two parts, p in one, q in the other, and r in one of them, say together with p. The edge rq joins the parts, and the tree with it instead of pq is shorter, a contradiction. So the circle on the diameter pq is empty, and pq is a Delaunay edge.
The second fact turns the minimum spanning tree problem for points from “consider all n(n − 1)/2 pairs” into “consider fewer than 3n edges”. The listing checks this on the stations: the tree over 136 pairs and over 43 Delaunay edges is the same, 77.3 km. In the article on the minimum spanning tree the towns of a whole network are connected this way, and the saving is visible there too.
The third fact is about the largest empty circle. Within one Voronoi cell, the distance from a point of the plane to the nearest site grows with distance from the site, so in a bounded region the maximum is reached at a Voronoi vertex or on the edge of the region: where the edge passes from one cell to another, or at a vertex of the edge. This is a finite list of candidates: in the first scene there are 7252 of them, Voronoi vertices inside the city and points on district boundaries where the nearest station changes.
How to build it: Bowyer–Watson
The simplest way is to insert the points one at a time, maintaining the Delaunay triangulation of those already inserted. Start with a huge frame triangle that contains all the points. A new point p falls inside the circumcircles of several triangles, and these triangles have stopped being Delaunay triangles. Remove them: what remains is a hole, a polygon every vertex of which is visible from p. Join p to each edge of the hole's boundary. At the end throw away the triangles that touch the vertices of the frame.
Why the new triangles are Delaunay again: the circle of any new triangle pab is empty, because the old triangle beyond the edge ab was not knocked out, that is, p did not fall into its circle, and both circles pass through a and b. A careful proof that the hole is connected and entirely visible from p is in the book by de Berg and co-authors, “Computational Geometry”; on our data it is confirmed by the check in the listing: all circles are empty. For 17 stations 55 triangles are knocked out along the way, a little over three per point on average.
What is expensive here is something else: finding the triangles to knock out. The simple version checks all triangles, and for n points that is of the order of n2 checks. The fast one finds the triangle the point fell into by walking from the previous one, stepping across the edge the point lies to the right of, and looks for the other triangles to knock out only among neighbours: the hole is connected. If the points are inserted in random order, the expected number of knocked-out triangles per point is bounded by a constant, and all the work is O(n log n); Leonidas Guibas, Donald Knuth and Micha Sharir showed this in 1992. The scene below builds the triangulation of all the crossings of the central part of the city this way.
7909 points give 15,800 triangles, exactly 2n − 2 − h: this cloud has 16 points on its hull. Each insertion knocks out fewer than six triangles on average, and there are about a hundred thousand “inside the circle” checks instead of more than sixty million for the simple version. Switch to Voronoi cells: streets show as chains of elongated cells, blocks as large ones.
Precision is a separate topic. The “inside the circle” check is the sign of a determinant, and on nearly collinear points, of which a street grid has many, floating-point error can flip it, after which the triangulation breaks: triangles overlap or holes appear. The scene computes the sign first quickly, in floating point, with an error bound after Jonathan Shewchuk, and when the bound does not vouch for the answer, recomputes it in long arithmetic. The Python listing computes in integers from the start, with coordinates in metres, so there are no rounding errors there.
Flips: from any triangulation
There is a second way, and it explains what makes the Delaunay triangulation good. Take any triangulation, for example a “fan”: points by increasing x, each new one joined to all the hull edges it can see. Look at the edge ac between the triangles abc and acd. If d is inside the circle of abc, the edge is bad: replace it with the other diagonal of the quadrilateral, bd. This is a flip. Repeat while bad edges exist. Charles Lawson described this method in 1977.
For the stations the fan comes out with 15 angles under 10°. After 13 flips there are six such angles, and the result is exactly the Delaunay triangulation. This is no coincidence: among all triangulations of given points, Delaunay's has the largest smallest angle (Lawson, 1977), and if the sorted angles are compared in order, it is the best over the whole list (Robin Sibson, 1978). That is why it is used for meshes: in the finite element method elongated triangles spoil accuracy, and Delaunay avoids them as far as the points allow. The smallest angles, 0.5° and 0.6°, remain after the flips too. They belong to one triangle: three northern stations, No. 19, SPCh-14 and PCh-4, stand almost on one line near the edge of the hull, and since Delaunay has the largest smallest angle, any other triangulation has one no larger.
The paraboloid
Why flips eventually end, and why empty circles are enough for the triangulation to be unique, shows in three dimensions. Lift each point (x, y) onto the paraboloid, to the point (x, y, x2 + y2).
The point d lies inside the circle through a, b, c if and only if the lifted point d′ lies below the plane through the lifted a′, b′, c′.
The plane through a′, b′, c′ is z = αx + βy + γ for some α, β, γ (it is not vertical, since a, b, c are not on one line). The lifted point (x, y) lies below it when x2 + y2 < αx + βy + γ, that is, (x − α/2)2 + (y − β/2)2 < γ + (α2 + β2)/4. This is the interior of a disc. Its boundary is the points where the lifted point lies in the plane itself, and a′, b′, c′ lie there. So it is the disc bounded by the circle through a, b, c.
Everything else follows. The triangle abc is Delaunay when all lifted points are above its plane, that is, when the lifted triangle is a face of the lower side of the convex hull of the lifted points. The Delaunay triangulation is the projection of the lower hull, so it exists, and if no four lifted points lie in one plane, that is, no four points are on one circle, it is unique. Flipping a bad edge replaces two lifted triangles with two lying lower: the surface of triangles keeps going down, cannot repeat, and there are finitely many triangulations, so the flips end.
If a triangulation has no bad edges, it is the Delaunay triangulation.
Lift all the triangles onto the paraboloid: this gives a surface over the convex hull of the points, a function f(x, y) linear on each triangle. An edge is not bad, so by the claim above the lifted opposite vertex of the neighbouring triangle lies no lower than the plane of this triangle: at each edge the surface bends upwards or does not bend. A function that is piecewise linear and convex at every joint between pieces is convex on the whole convex region. The lower hull of the lifted points is the largest convex function not exceeding the lifted values at the points themselves. Our f is convex and coincides with them at the points, so f is no higher than the hull. And on each triangle the hull is no higher than the linear interpolation of the values at its vertices, by convexity, that is, no higher than f. So f is the lower hull, and the triangulation is its projection.
The same link between empty circles and convex hulls was noticed by Kevin Brown in 1979, and it gives one more consequence: building a Delaunay triangulation is no easier than sorting numbers. The numbers x1, …, xn are turned into points on a parabola, as in the article on the convex hull, and the order can be read from the triangulation in O(n). So with comparisons nothing faster than O(n log n) is possible.
Code
The listing builds the triangulation of the 17 stations by Bowyer–Watson in integers, checks that all circles are empty and Euler's formula, compares the minimum spanning tree over all pairs and over the Delaunay edges, looks for the place furthest from the stations in the whole city and among the crossings, and at the end builds a fan and flips its edges until the same result comes out.
In Go the “inside the circle” check with a huge frame does not fit in 64 bits, and the listing computes it with math/big; in Python integers are long by themselves. The search for the furthest place goes over all sides of the district boundaries: on each side it starts from the station nearest to its start and finds where along the side another station becomes the nearest; these are the points where the side crosses Voronoi edges.
What it costs
The output is linear: fewer than 2n triangles and 3n edges. The time depends on the method.
| Method | Time | Remark |
|---|---|---|
| Bowyer–Watson, checking all triangles | O(n2) | simple, fine for hundreds of points |
| Bowyer–Watson, walking and random order | O(n log n) on average | worst case O(n2); the most common in practice |
| Lawson flips from an arbitrary triangulation | O(n2) flips in the worst case | handy for repairing a mesh after edits |
| Fortune's sweep line | O(n log n) | builds the Voronoi diagram directly |
| divide and conquer | O(n log n) | split by x and stitch two triangulations together |
After building, the query “nearest site to a point” is answered by descending a point location structure in O(log n); in practice a k-d tree is often simpler, and for a dozen sites a scan, as in the scene. Memory is O(n): a triangle has three vertices and three neighbours.
Variants
By road, not in a straight line. The real zone of a fire station is determined by travel time. This is a Voronoi diagram on a graph: Dijkstra from all stations at once, with one shared queue, colours each crossing in the colour of the station it is closer to. Rivers, mountains and one-way streets make such zones look nothing like polygons.
With given edges. If some segments must be in the triangulation, such as a coastline, a road or the edge of a part, a constrained Delaunay triangulation is built: circles are empty, but points hidden behind a required edge do not count. By adding new points in the right places, meshes are made from it in which all triangles have angles no smaller than a given one; this is how mesh generators for simulations work.
With weights. If the sites have different “strengths”, such as transmitter power or the size of a shop, the boundary between cells shifts or curves: this gives power diagrams and multiplicatively weighted diagrams.
On a sphere. For points on a globe no lift is needed: the Delaunay triangulation on a sphere is the convex hull of the points themselves in space, for the same reason that a circle on the plane is a section of the paraboloid by a plane.
Where it is used
Maps and logistics. The nearest facility, a service zone, a place for a new facility: the problems of this article. Hydrologists compute average rainfall over an area by dividing it into station zones; they call these zones Thiessen polygons, after the meteorologist who proposed them in 1911.
Terrain and interpolation. From scattered points with heights a surface of triangles is built, an irregular network, and the Delaunay triangulation is used for it: long thin triangles would link distant points. In 1981 Robin Sibson proposed natural neighbour interpolation, which takes its weights from how much area a new point takes away from the neighbouring Voronoi cells.
Simulations. The finite element method divides a part into triangles or tetrahedra, and accuracy depends on the shape of the cells. Mesh generators build a constrained Delaunay triangulation and refine it.
Graphics and games. The Voronoi cells of random points give a “cellular” texture: stone, scales, cracks; the function that returns the distance from a point to the nearest site was described by Steven Worley in 1996, and it became a standard noise in graphics. Map generators divide the world into Voronoi cells, and a path that keeps away from obstacles runs along the edges of the Voronoi diagram.
History
Cells of the nearest point appeared in number theory. In 1850 Gustav Lejeune Dirichlet considered them for lattice points while studying quadratic forms, and in part of the literature they are called the Dirichlet tessellation. Georgy Voronoy in 1908 developed the theory for lattices in space of any dimension, and it is Voronoy whose name stuck in work on computation. In 1934 Boris Delaunay published the paper “Sur la sphère vide” (“On the empty sphere”), dedicated to Voronoy's memory: it is where the empty circle condition and the triangulation dual to the cells appeared. The same cells were found independently in meteorology (Alfred Thiessen, 1911) and in crystal physics (Eugene Wigner and Frederick Seitz, 1933).
Voronoi diagrams were brought into computer science by Michael Shamos and Dan Hoey in 1975: they showed that many problems about the closeness of points are solved through them in O(n log n). Lawson's flips date from 1977, Sibson's largest smallest angle property from 1978, Brown's link with the paraboloid from 1979. Adrian Bowyer and David Watson independently described insertion with knocking out triangles in the same issue of The Computer Journal in 1981, and Steven Fortune in 1987 the sweep line that builds the Voronoi diagram directly. The O(n log n) bound for random insertion order was proved by Guibas, Knuth and Sharir in 1992.
Problems
How many triangles Where to put a new stationData: Almaty fire stations © OpenStreetMap contributors, ODbL (export of 17 Oct 2026) Almaty roads © OpenStreetMap contributors, ODbL