xx. Delaunay and Voronoi

Delaunay

A triangulation of a point set in which no point lies inside the circumcircle of any triangle. It is the dual of the Voronoi diagram, and it is unique when no four points are concyclic.

Why it is the “right” triangulation

Among all triangulations of the same points, Delaunay:

• maximises the minimum angle (avoids skinny triangles) — and in fact maximises the whole sorted angle vector lexicographically;
• minimises the maximum circumcircle radius;
• contains the closest pair as an edge;
• contains the Euclidean minimum spanning tree as a subgraph;
• contains the nearest-neighbour graph as a subgraph.

The last two are the reasons it matters algorithmically.

Proof: Delaunay maximises the minimum angle, and contains the closest pair

Minimum angle. Consider two triangles sharing an edge $AC$, with opposite vertices $B,D$ forming a (possibly non-convex) quadrilateral. The Delaunay flip lemma: if $ABCD$ is convex and $D$ is inside the circumcircle of $ABC$ (i.e. the diagonal $AC$ is the “wrong” one), flipping to diagonal $BD$ strictly increases the minimum of the six angles of the two triangles. This is a direct computation: the sum of the two angles at $B$ and $D$ (opposite the shared edge $AC$) is $>180°$ when $D$ is inside circle $ABC$ (by the inscribed angle theorem — the locus where $\angle ABC+\angle ADC=180°$ is exactly the circle through $A,B,C$, and being inside makes the angle at $D$ larger, pushing the sum over), and swapping the diagonal is known to replace the two largest-angle-sum triangles with the configuration whose minimum angle is provably no smaller (a classical inequality on the six angles). Since any non-Delaunay triangulation has some edge violating the in-circle test, repeatedly flipping (each flip strictly increases the lexicographic angle vector, so there’s no cycling) terminates only at the Delaunay triangulation, which must therefore be a global — not just local — maximum of the minimum angle, since every other triangulation is connected to it by a sequence of flips, each non-decreasing the minimum angle.
Contains the closest pair. Suppose the closest pair $(p,q)$ is not an edge of some triangulation $T$. Then the segment $pq$ crosses some edge of $T$, say $(r,s)$, with $r,s$ on opposite sides. By the triangle inequality in the two triangles formed, either $|pr|<|pq|$ or $|qs|<|pq|$ or similar (one of the four “other” distances among $p,q,r,s$ is shorter than $pq$) — contradicting that $(p,q)$ is the closest pair overall. So the closest pair must be an edge of every triangulation, Delaunay included. $■$

Construction

Method	Time	Note
Flip algorithm	$O(n^2)$	start anywhere, flip illegal edges until none remain
Bowyer-Watson	$O(n\log n)$ expected	incremental; conceptually the simplest
Divide and conquer (Guibas-Stolfi)	$O(n\log n)$ worst case	quad-edge structure
Sweep line (Fortune)	$O(n\log n)$	build Voronoi and dualise
Lift to 3D + convex hull	$O(n\log n)$	the most elegant, and reuses [[19-geometry/19-3d-geometry

The lifting map

Send $(x,y)\mapsto (x,y,x^2+y^2)$ onto the paraboloid. Then: the lower convex hull of the lifted points, projected back to the plane, is the Delaunay triangulation. A point is inside a circumcircle in 2D exactly when its lift is below the corresponding plane in 3D — so the in-circle predicate becomes an orientation test in 3D (proved in Primitives, via the same lifting argument). If you have a 3D hull routine, you get Delaunay for free, and it inherits the 3D hull’s exactness.

The in-circle predicate — the hard part

$\mathrm{inCircle}(a,b,c,d)=\mathrm{sign}\left|\begin{array}{ccc}{a}_x-d_x& a_y-d_y& |a-d{|}^2\\ b_x-d_x& b_y-d_y& |b-d{|}^2\\ c_x-d_x& c_y-d_y& |c-d{|}^2\end{array}\right|$

with $a,b,c$ counter-clockwise; positive means $d$ is strictly inside. (Proved via the paraboloid lifting in Primitives.)

This determinant is degree 4

With coordinates up to $10^9$ it reaches $10^{45}$ — far beyond long long and beyond double’s 53-bit mantissa. Getting the sign wrong produces a non-triangulation, not a slightly wrong one.

Options: restrict coordinates to $\le 10^4$ and use __int128; use exact rational or big-integer arithmetic; or use adaptive-precision predicates (Shewchuk). This is the single biggest obstacle to implementing Delaunay correctly — in practice it is also the single most common source of a subtly broken from-scratch implementation, which is why this note gives the algorithm and the predicate in full but does not publish a hand-rolled Bowyer-Watson or Fortune implementation as “tested code”: getting every degenerate case (near-cocircular points, collinear hull boundaries, duplicate points) exactly right needs the adaptive-precision machinery referenced above, and a version without it is a liability dressed up as a solution. If you need Delaunay in a contest, reach for a library implementation (e.g. a well-tested boost::polygon / CGAL-style routine, or a vetted competitive-programming library snippet with its own stress-tested history) rather than deriving one under time pressure.

Uses

Problem	How
Euclidean MST	Delaunay has $O(n)$ edges and contains the EMST $\to$ Kruskal on those
Closest pair	the shortest Delaunay edge
Nearest neighbour graph	a subgraph
Largest empty circle	centred at a Voronoi vertex = a Delaunay circumcentre
Terrain / mesh generation	maximising the minimum angle avoids numerical trouble
Interpolation (natural neighbour)	uses the Voronoi structure
Point location	walk the triangulation

Euclidean MST, concretely

1. Build the Delaunay triangulation, $O(n\log n)$.
2. Take its $O(n)$ edges with Euclidean weights.
3. Run Kruskal, $O(n\log n)$.

Versus $O(n^2\log n)$ on the complete graph. For $n=10^5$ that is the difference between feasible and not.

Proof: the Delaunay triangulation contains the Euclidean MST

Let $(p,q)$ be an edge of the EMST, and suppose it is not a Delaunay edge. Then there exists a circle through $p,q$ with no other point inside it only if $(p,q)$ is Delaunay (an edge is Delaunay iff some empty circle passes through its two endpoints) — so if $(p,q)$ is not Delaunay, every circle through $p,q$ contains some other point $r$ inside it. Taking the circle with $pq$ as diameter (the smallest such circle), $r$ inside it means $\angle prq>90°$, which implies $\max (|pr|,|qr|)<|pq|$ (the side opposite the largest angle in a triangle is the longest; since $\angle prq$ is obtuse, it’s the triangle’s largest angle, so $pq$ — opposite it — is the longest side, meaning both $pr$ and $qr$ are shorter than $pq$). But then replacing edge $(p,q)$ in the MST with $(p,r)$ and $(r,q)$ (removing $(p,q)$ disconnects the tree into two components; $r$ is in one of them, say with $p$, so adding $(r,q)$ reconnects with strictly smaller total weight than keeping $(p,q)$ and whichever of $r$‘s existing tree edges — more carefully, the standard cycle-property argument: $(p,q)$ being the maximum-weight edge on some cycle formed with a shorter path would contradict MST-optimality) contradicts minimality. So every EMST edge must be Delaunay. $■$

Constrained Delaunay

Forces certain edges to appear (for meshing a polygon with fixed boundaries). The result is no longer strictly Delaunay — triangles may violate the empty-circle property, but only because of the forced edges.

3D and beyond

Delaunay tetrahedralisation exists in 3D but can have $\Theta (n^2)$ tetrahedra, and its construction is much harder. The lifting map still works (to 4D), and the in-sphere predicate has degree 5.

The contest reality

Full Delaunay implementations are long and precision-critical. Before writing one, check whether you actually need it:

• Euclidean MST with $n\le 5000$ $\to$ the complete graph is $1.2\cdot 10^7$ edges; just run Kruskal.
• Nearest neighbour queries $\to$ a KD-tree.
• Closest pair $\to$ the sweep.
• Largest empty circle $\to$ enumerate candidate centres.

Delaunay earns its place at $n\ge 10^5$ with a genuine Euclidean MST requirement.

Bowyer-Watson, conceptually

Build incrementally, $O(n\log n)$ expected with randomised insertion order ($O(n^2)$ worst case):

1. Start with a super-triangle large enough to contain every input point.
2. Insert points one at a time (in random order). For each new point $p$:
   • find all triangles whose circumcircle contains $p$ (the in-circle predicate above) — the “bad” triangles;
   • their union is a star-shaped cavity; collect its boundary edges (edges belonging to exactly one bad triangle);
   • delete the bad triangles, and create one new triangle from $p$ to each boundary edge.
3. At the end, delete every triangle touching a super-triangle vertex.

Proof: Bowyer-Watson is correct

Invariant. After inserting a prefix of the points (plus the super-triangle), the current triangle set is exactly the Delaunay triangulation of those points (and the super-triangle’s vertices). Finding bad triangles. A triangle’s circumcircle containing the new point $p$ is exactly the condition that makes that triangle not Delaunay once $p$ is added (violating the empty-circle property) — so it must be removed. The cavity is star-shaped from $p$. Because the bad triangles’ union is the set of points whose Delaunay-defining circumcircle is invalidated by $p$, and circumcircles vary continuously across adjacent triangles sharing an edge, the bad region is connected and — crucially — every boundary edge of the cavity is visible from $p$ (no other bad-region geometry occludes it), which is what makes “connect $p$ directly to every boundary edge” produce a valid, non-self-intersecting retriangulation. Boundary edges via counting. An edge shared by two bad triangles is interior to the cavity and disappears; an edge belonging to exactly one bad triangle is on the cavity’s boundary — the standard “appears exactly once” trick (the same idea used for 3D hull’s horizon edges). $■$

Delaunay construction methods compared

Method	Time	Notes
Bowyer-Watson	$O(n\log n)$ expected	easiest to understand
Divide and conquer (Guibas-Stolfi)	$O(n\log n)$ worst case	quad-edge structure
Sweep line (Fortune)	$O(n\log n)$	see [[#Fortune, conceptually
Lift to 3D + convex hull	$O(n\log n)$	lift $(x,y)\to (x,y,x^2+y^2)$; the lower hull is the Delaunay triangulation
Flip algorithm	$O(n^2)$	start from any triangulation, flip illegal edges

Fortune, conceptually

Constructs the Voronoi diagram of $n$ points in $O(n\log n)$ using a sweep line and a “beach line” of parabolic arcs — one of the most intricate standard algorithms in computational geometry to implement fully (hundreds of lines in a production-quality version), presented here for the technique and its correctness, not as code to hand-roll under contest time pressure (same caveat as the in-circle predicate above).

Algorithm

1. Sweep a horizontal line downward across the plane.
2. Maintain a beach line: an ordered sequence of parabolic arcs, where each arc is the set of points equidistant from one site and the sweep line.
3. Maintain two event queues ordered by sweep position: site events (the line reaches a new point) and circle events (three consecutive arcs are about to converge to nothing).
4. On a site event, find the arc directly above the new site, split it in two around a new zero-width arc for the site, and check the newly formed arc triples for potential circle events.
5. On a circle event, the middle arc of the converging triple vanishes; its convergence point is a Voronoi vertex. Record the diagram edges meeting there, remove the arc, and check the newly adjacent arcs for fresh circle events.
6. Continue until both queues are empty. The final beach line traces the diagram’s unbounded edges, completing the Voronoi diagram.

Proof: Fortune's algorithm is correct

Invariant. At sweep position $y$, the beach line is exactly the lower envelope of parabolas focused at all already-processed sites with directrix at the sweep line — i.e. every point below the beach line is strictly closer to some processed site than to the sweep line, so its Voronoi cell assignment is already finalised and will never change (any site below the current sweep line, not yet processed, cannot be closer, since it’s farther from every point above the sweep line than the sweep line itself is).
Site events. When the sweep reaches a new site $s$, the arc directly above $s$ belongs to the closest already-processed site at that $x$-coordinate (by the invariant). Splitting that arc and inserting a zero-width arc for $s$ correctly extends the invariant to include $s$.
Circle events. A Voronoi vertex is a point equidistant from three (or more) sites with no site closer. This occurs exactly when three consecutive beach-line arcs converge to a single point as the sweep advances — precisely the condition a circle event detects, since the circle through those three sites has its bottommost point exactly on the sweep line at that moment, with no other site inside it (by the beach-line invariant).
Completeness and termination. Each of the $n$ sites generates exactly one site event, and each Voronoi vertex corresponds to exactly one valid circle event (the circle through its three defining, mutually adjacent sites, empty of any other site — at most $2n-5$ such vertices in a planar Voronoi diagram). Since every event either inserts, splits, or removes a bounded piece of beach-line structure, the queues hold $O(n)$ total events and the algorithm terminates having constructed every edge and vertex exactly once. $■$

Variants / Use Cases

• Delaunay triangulation — the dual graph of the Voronoi diagram, obtainable directly as a byproduct of Fortune’s sweep.
• Nearest-neighbour search / closest pair of points — Voronoi diagrams answer nearest-site queries directly.
• Weighted Voronoi diagrams — additively or multiplicatively weighted variants for uneven site “importance”.
• Higher-order Voronoi diagrams — generalisation to $k$-th nearest neighbour regions.

Voronoi

Partition the plane by “which site is nearest”. The Voronoi cell of a site $s$ is
$V(s)=\{p:|p-s|\le |p-t|\forall t\}.$
Each cell is a convex polygon (possibly unbounded) — the intersection of the half-planes defined by the perpendicular bisectors to every other site.

Structure

For $n$ sites in general position:

Quantity	Value
Cells	$n$
Edges	$\le 3n-6$
Vertices	$\le 2n-5$
Total complexity	$O(n)$

Voronoi vertices are circumcentres of triples of sites, and each is equidistant from three sites with no fourth inside that circle.

Proof: a Voronoi cell is a convex polygon, and the edge/vertex bounds

Convexity. $V(s)={\bigcap }_{t\ne s}\{p:|p-s|\le |p-t|\}$, an intersection of half-planes (each “closer to $s$ than to $t$” region, bounded by the perpendicular bisector of $st$) — an intersection of convex sets (half-planes) is convex.
Complexity bounds. The Voronoi diagram’s edges and vertices are the dual of the Delaunay triangulation’s, which (as a planar graph) satisfies Euler’s formula $V-E+F=2$; a triangulation of $n$ points has at most $2n-5$ interior triangles and $3n-6$ edges (standard planar-graph bounds for a maximal planar graph, since each face other than the outer one is a triangle), and dualising maps Delaunay triangles to Voronoi vertices, giving the matching $O(n)$ bounds on the Voronoi side. $■$

Construction

Method	Time	Note
Half-plane intersection per cell	$O(n^2\log n)$	simple, fine for $n\le 2000$
Fortune’s sweep	$O(n\log n)$	the classic
Dual of Delaunay	$O(n\log n)$	build Delaunay, then take the dual
Lift to 3D + convex hull	$O(n\log n)$	elegant; via the paraboloid
Divide and conquer	$O(n\log n)$

Build Delaunay instead

The Voronoi diagram and the Delaunay triangulation are duals: a Voronoi edge separates two sites exactly when those sites are joined by a Delaunay edge. Delaunay is the one to build if you need the explicit structure — but see the precision caveat above before committing to a from-scratch implementation.

The duality in detail

Voronoi	Delaunay
cell (a site)	vertex
edge (between 2 cells)	edge
vertex (meeting of 3 cells)	triangle
unbounded cell	a hull vertex
Voronoi vertex position	the triangle’s circumcentre

What it answers

Query	Via Voronoi/Delaunay
Nearest site to a query point	point location in the Voronoi diagram, $O(\log n)$
Closest pair	the shortest Delaunay edge
Euclidean MST	a subgraph of Delaunay — run Kruskal on $O(n)$ edges
Nearest neighbour graph	a subgraph of Delaunay
Largest empty circle	centred at a Voronoi vertex (or on the hull boundary)
Farthest-point Voronoi	gives the minimum enclosing circle
All nearest neighbours	one pass over Delaunay edges
Motion planning with maximum clearance	travel along Voronoi edges

The Euclidean MST reduction is the most valuable: it takes an $O(n^2)$ complete-graph problem down to $O(n\log n)$.

Variants

• Weighted (power) Voronoi — cells defined by $|p-s{|}^2-w_s$; still polygonal, dual to the regular triangulation.
• Multiplicatively weighted — cells bounded by Apollonius circles, no longer polygonal.
• Farthest-point Voronoi — “which site is farthest”; only hull vertices have non-empty cells.
• Order-$k$ Voronoi — regions with the same set of $k$ nearest sites.
• Voronoi under $L_1$ / $L_{\infty }$ — cells are still polygonal but with axis-parallel and diagonal edges.
• Voronoi of segments — cells bounded by parabolic arcs; substantially harder.

In competitive programming

Explicit Voronoi construction is rare — it is long, precision-sensitive, and usually avoidable. The three things worth knowing:

1. Nearest-site queries $\to$ a KD-tree is simpler and fast enough.
2. Euclidean MST $\to$ build Delaunay, then Kruskal.
3. Largest empty circle $\to$ its centre is a Voronoi vertex, a hull edge midpoint, or determined by the bounding region; enumerate the $O(n)$ candidates.

For small $n$, the $O(n^2\log n)$ per-cell half-plane intersection is entirely acceptable and far easier to get right.