Delaunay Triangulation

Delaunay triangulation is a core geometric structure in computational geometry and scientific computing, introduced by Soviet mathematician Boris Nikolaevich Delone (Delaunay) in 1934. Given a finite set P = {p₁, p₂, …, p_n} of points in the Euclidean plane ℝ² (not all collinear), a Delaunay triangulation DT(P) is a triangulation of the convex hull conv(P) having vertex set P such that the open circumcircle of every triangle in DT(P) contains no point of P in its interior. In d-dimensional space ℝ^d, this empty circumsphere property generalizes to decomposing conv(P) into d-simplices (such as tetrahedra in ℝ³) whose open circumspheres contain no points of P.

Under the planar general position assumption that no four points of P are cocircular, the Delaunay triangulation DT(P) is unique and forms the planar straight-line dual graph of the Voronoi diagram Vor(P). Because the empty-circumcircle condition maximizes the minimum interior angle across all triangulations of P—avoiding thin "sliver" triangles that cause large interpolation errors and ill-conditioned stiffness matrices—Delaunay triangulation serves as the foundation for finite-element mesh generation, terrain digital elevation modeling, natural-neighbor spatial interpolation, and neighborhood graph construction for DBSCAN and k-means clustering.

Two vertices p_i, p_j ∈ P are connected by a Delaunay edge if and only if there exists a closed disk passing through p_i and p_j that contains no other point of P. Equivalently, in terms of the Voronoi diagram Vor(P):

• Two sites p_i and p_j share a Delaunay edge in DT(P) if and only if their Voronoi cells V(p_i) and V(p_j) share a boundary segment of positive length.
• Three sites p_i, p_j, p_k form a Delaunay triangle in DT(P) if and only if their three Voronoi cells meet at a common Voronoi vertex—the circumcenter of triangle Δp_ip_jp_k.

By Euler's polyhedral formula V − E + F = 2 applied to the planar graph of n vertices whose outer boundary of conv(P) has h ≥ 3 hull vertices, every triangulation of P in ℝ² contains 2n − 2 − h triangles and 3n − 3 − h edges (both O(n)). In dimensions d ≥ 3, Klee's and Seidel's upper bound theorem shows that a d-dimensional Delaunay triangulation of n points can contain up to Θ(n^⌈d/2⌉) simplices in the worst case.

Paraboloid Lifting to a (d+1)-Dimensional Convex Hull

In 1979, Kevin Q. Brown established an equivalence between Delaunay triangulations in ℝ^d and convex hulls in ℝ^d+1. Mapping each point p = (x₁, …, x_d) ∈ ℝ^d onto the rotational paraboloid z = ‖p‖₂² in ℝ^d+1 sends every sphere in ℝ^d to the intersection of the paraboloid with a hyperplane in ℝ^d+1. Projecting the lower facets of the (d + 1)-dimensional convex hull of the lifted points {(p_i, ‖p_i‖₂²)} down onto ℝ^d yields DT(P).

The InCircle Determinant Predicate and Sibson's Max-Min Angle Theorem

Given a counterclockwise-oriented triangle (a, b, c) in ℝ² and a query point d, whether d lies inside the circumcircle of Δabc is decided by the sign of the 4 × 4 InCircle determinant:

InCircle(a, b, c, d) = det [ a_x−d_x, a_y−d_y, (a_x−d_x)²+(a_y−d_y)² ; b_x−d_x, b_y−d_y, (b_x−d_x)²+(b_y−d_y)² ; c_x−d_x, c_y−d_y, (c_x−d_x)²+(c_y−d_y)² ] > 0.

In floating-point geometry kernels (such as Shewchuk's Triangle library), adaptive-precision floating-point arithmetic evaluates InCircle to prevent topological inconsistencies. Robin Sibson (1978) proved that among all triangulations of P in ℝ², DT(P) lexicographically maximizes the sorted vector of all 3t interior angles (α₁ ≤ α₂ ≤ ⋯ ≤ α_3t).

Four main algorithmic paradigms compute DT(P) in the plane in optimal O(n log n) time:

Lawson's Edge-Flipping and the Bowyer–Watson Cavity Algorithm

In 1977, Charles L. Lawson proved that any triangulation of a planar point set P can be transformed into DT(P) by at most O(n²) local edge flips: whenever two adjacent triangles Δabc and Δcbd sharing an interior edge bc form a convex quadrilateral with ∠cab + ∠bdc > 180° (violating InCircle), replacing edge bc with diagonal ad increases the minimum angle of the six angles in the two triangles.

The Bowyer–Watson algorithm (1981) generalizes to arbitrary dimension d: starting from a bounding super-simplex, points p_k are inserted one by one along a Hilbert curve to achieve O(n log n) expected time. For each point p_k, a graph walk finds all simplices whose circumspheres contain p_k; deleting these leaves a star-shaped polygonal Delaunay cavity, whose boundary facets are connected to p_k.

Every Delaunay triangulation in ℝ^d contains a nested chain of proximity graphs on P:

Euclidean Minimum Spanning Tree (EMST) ⊆ Relative Neighborhood Graph (RNG) ⊆ Gabriel Graph (GG) ⊆ Delaunay Triangulation (DT).

Because the Euclidean Minimum Spanning Tree of n planar points is always a subgraph of DT(P) (which has at most 3n − 6 edges), constructing DT(P) in O(n log n) time and running Kruskal's algorithm on its O(n) edges solves the planar Euclidean MST problem in O(n log n) time.

In finite-element simulation (where solving PDEs uses Cholesky decomposition, LU decomposition, or the conjugate gradient method) and in terrain navigation for a particle filter or simplex optimization via the Nelder–Mead method, Constrained Delaunay Triangulation (CDT) forces prescribed domain boundary edges into the mesh, and Ruppert's Delaunay refinement inserts circumcenters of skinny triangles (Steiner points) until every triangle has all interior angles ≥ 20.7°.

What happens when four or more points lie on the same empty circle?

When k ≥ 4 points of P are cocircular on an empty circle, the Delaunay decomposition defines a convex k-gon, which can be subdivided into k − 2 triangles in multiple ways sharing the same circumcircle. Symbolic perturbation or deterministic vertex ordering selects a unique triangulation.

Does the 3D Delaunay tetrahedralization still maximize the minimum dihedral angle?

In 3D (d = 3), the empty-circumsphere criterion minimizes the maximum enclosing sphere radius (Rajan's theorem), but it does not maximize the minimum dihedral angle: four nearly coplanar points around the equator of a sphere can form a flat, near-zero-volume sliver tetrahedron whose circumsphere is empty. 3D meshers eliminate slivers using weighted Delaunay triangulations and sliver exudation.

1. Delaunay, B. (1934). Sur la sphère vide: A la mémoire de Georges Voronoï. Bulletin de l'Académie des Sciences de l'URSS, Classe des Sciences Mathématiques et Naturelles, 6, 793–800. https://doi.org/10.1007/978-3-662-04245-8_1
2. Guibas, L., & Stolfi, J. (1985). Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams. ACM Transactions on Graphics, 4(2), 74–123. https://doi.org/10.1145/282918.282923
3. Shewchuk, J. R. (1997). Adaptive precision floating-point arithmetic and fast robust geometric predicates. Discrete & Computational Geometry, 18(3), 305–363. https://doi.org/10.1007/PL00009321
4. Cheng, S.-W., Dey, T. K., & Shewchuk, J. R. (2012). Delaunay Mesh Generation. Chapman and Hall/CRC. https://doi.org/10.1201/b12987