Delaunay Triangulation

Definition and the Empty Circumcircle Property

A triangulation of a point set $S$ in $\mathbb{R}^2$ is a maximal set of non-crossing edges connecting points of $S$, decomposing the convex hull into triangles.

The Delaunay triangulation $\text{DT}(S)$ is the unique triangulation (assuming general position: no four cocircular points) satisfying the empty circumcircle property:

The circumscribed circle of every triangle contains no other site in its interior.

Equivalently, $\text{DT}(S)$ maximizes the minimum angle over all triangulations of $S$. This makes it the "most well-shaped" triangulation, avoiding skinny triangles when possible.

Key Properties

• Uniqueness: unique when no four points are cocircular. With cocircular points, the Delaunay triangulation is not unique (can flip the shared diagonal), but the Delaunay complex (which may include non-triangular faces) is unique
• Planarity: $O(n)$ triangles, edges, and vertices (at most $2n - 5$ triangles, $3n - 6$ edges for $n \geq 3$ points)
• Max-min angle: among all triangulations, DT maximizes the sorted vector of angles (lexicographically). This is the angle optimality property
• Closest pair: the edge connecting the closest pair of points in $S$ is always a Delaunay edge
• EMST subset: the Euclidean minimum spanning tree is a subgraph of the Delaunay triangulation
• Nearest neighbor graph: every point's nearest neighbor is connected to it by a Delaunay edge

Duality with Voronoi Diagram

The Delaunay triangulation is the dual graph of the Voronoi diagram:

• Sites $s_i$ and $s_j$ are connected by a Delaunay edge iff their Voronoi cells $V(s_i)$ and $V(s_j)$ share an edge
• Three sites form a Delaunay triangle iff their Voronoi cells share a common vertex

Circumcenter = Voronoi vertex: the circumcenter of each Delaunay triangle is the corresponding Voronoi vertex.

Lifting interpretation: map each point $(x, y)$ to $(x, y, x^2 + y^2)$ on the paraboloid $z = x^2 + y^2$. The Delaunay triangulation is the projection of the lower convex hull of the lifted points. This connects Delaunay to convex hull algorithms in $\mathbb{R}^3$.

Construction Algorithms

Incremental Insertion

Algorithm (Bowyer-Watson, 1981): insert points one at a time.

For each new point $p$:

1. Locate: find the triangle containing $p$ (or the edge if $p$ is on it)
2. Identify violated triangles: find all triangles whose circumcircle contains $p$ (the "cavity"). These form a star-shaped polygon
3. Retriangulate: delete violated triangles, connect $p$ to the boundary of the cavity

Point location: walk from a known triangle toward $p$ using orientation tests. With randomized insertion order, expected walk length is $O(\sqrt{n})$ without optimization, or $O(\log n)$ with a point location structure (e.g., DAG of triangle history).

Expected complexity: $O(n \log n)$ with random insertion order (each insertion takes $O(\log n)$ amortized by backward analysis: expected number of triangles in the cavity is $O(1)$ for random order).

Divide and Conquer

Algorithm (Lee and Schachter, 1980; Guibas and Stolfi, 1985): $O(n \log n)$ worst case.

1. Sort points by $x$-coordinate
2. Recursively triangulate left and right halves
3. Merge by finding the base edge (lowest mutual tangent) and computing the merge chain upward, performing InCircle tests to decide left vs. right candidate

The merge step uses the quad-edge data structure and takes $O(n)$ time. The algorithm is cache-efficient and highly parallelizable.

Edge Flipping

Start with any triangulation (e.g., obtained by incremental insertion without the Delaunay criterion). Repeatedly find a non-Delaunay edge (shared by two triangles forming a convex quadrilateral where the InCircle test fails) and flip it (replace diagonal).

Lawson's algorithm: flip non-Delaunay edges in any order. Terminates in $O(n^2)$ flips (since flips increase the minimum angle vector, and there are finitely many triangulations). Practical but not asymptotically optimal.

InCircle Predicate

The key geometric predicate for Delaunay construction. Given four points $a, b, c, d$ (with $a, b, c$ in counterclockwise order):

$\text{InCircle}(a, b, c, d) = \begin{vmatrix} a_x - d_x & a_y - d_y & (a_x - d_x)^2 + (a_y - d_y)^2 \\ b_x - d_x & b_y - d_y & (b_x - d_x)^2 + (b_y - d_y)^2 \\ c_x - d_x & c_y - d_y & (c_x - d_x)^2 + (c_y - d_y)^2 \end{vmatrix}$

Positive if $d$ is inside the circumcircle of $\triangle abc$; zero if cocircular; negative if outside. This is a $4 \times 4$ determinant (in the non-reduced form) and must be computed with exact arithmetic for robustness (Shewchuk's incircle predicate).

Constrained Delaunay Triangulation (CDT)

A constrained Delaunay triangulation respects a given set of constraint edges (which must appear in the triangulation), while being "as Delaunay as possible."

Definition: triangle $t$ satisfies the constrained empty circumcircle property if no vertex visible from $t$'s interior (without crossing a constraint edge) lies inside $t$'s circumcircle.

Properties:

• The CDT always exists for a planar straight-line graph (PSLG) input
• It is unique (given the constraints and assuming general position)
• It maximizes the minimum angle subject to the constraints
• Reduces to the standard DT when there are no constraints

Construction: insert constraint edges into an existing DT by:

1. Find and remove all triangles intersecting the constraint edge
2. Retriangulate the two resulting polygonal cavities
3. Flip non-constrained Delaunay edges as needed

Applications: mesh generation (constraint edges represent domain boundaries), terrain modeling (breaklines), geographic information systems.

Delaunay Refinement and Mesh Generation

Ruppert's Algorithm

Goal: produce a triangulation of a planar domain with no angles smaller than a threshold (typically 20-30 degrees).

Algorithm (Ruppert, 1995):

1. Start with the CDT of the input PSLG
2. If a triangle has circumradius-to-shortest-edge ratio exceeding a bound $B$ (i.e., has a small angle), insert the circumcenter of that triangle as a new point ("Steiner point")
3. If a new point encroaches on a constraint segment (lies inside its diametral circle), instead split the segment at its midpoint
4. Repeat until all triangles meet the quality criterion

Guarantees:

• Terminates with all angles $\geq \arcsin(1/(2B))$, which is about $20.7°$ for $B = \sqrt{2}$
• Output size: $O(n)$ times the optimal mesh size (proportional to local feature size)
• Graded: triangles are small near small input features and large far away

Chew's Algorithm

Chew's second algorithm (1993): similar to Ruppert's but inserts the circumcenter of the worst triangle (largest circumradius). Produces minimum angle $\geq 30°$ for convex domains. Simpler but less flexible than Ruppert's.

Quality Guarantees

The local feature size $\text{lfs}(p)$ at point $p$ is the distance from $p$ to the second-nearest feature of the input PSLG. Refinement algorithms produce meshes where edge lengths are proportional to local feature size:

$\text{edge length at } p = \Theta(\text{lfs}(p))$

This yields size-optimal meshes: the number of triangles is within a constant factor of any mesh achieving the same minimum angle.

Mesh Generation for FEM

The Delaunay triangulation and its refinements are the workhorse of finite element mesh generation:

2D Mesh Generation Pipeline

1. Input: domain boundary (PSLG)
2. CDT: construct constrained Delaunay triangulation
3. Refinement: Ruppert/Chew to eliminate bad triangles
4. Smoothing: Laplacian or optimization-based smoothing to improve mesh quality
5. Adaptation: refine/coarsen based on solution error estimates

3D Mesh Generation

Tetrahedral meshing: far more challenging. Key algorithms:

• Delaunay-based: incremental insertion with Bowyer-Watson, constrained Delaunay tetrahedralization (not always exists for arbitrary PSLCs), refinement via circumcenter insertion
• TetGen (Si, 2015): state-of-the-art Delaunay-based tetrahedral mesh generator. Constrained DT + Steiner point insertion + quality refinement
• CGAL 3D mesh generation: handles complex domains via Delaunay refinement with restricted Delaunay triangulation for surface approximation

Mesh Quality Metrics

• Aspect ratio: circumradius / inradius ($= 1$ for equilateral, $\to \infty$ for degenerate)
• Minimum angle: Delaunay maximizes minimum angle among all triangulations
• Edge ratio: max edge / min edge
• Condition number: relates to FEM solution accuracy

Delaunay in Higher Dimensions

In $\mathbb{R}^d$, the Delaunay triangulation has complexity $O(n^{\lceil d/2 \rceil})$:

• $d = 2$: $O(n)$ (optimal for geometric applications)
• $d = 3$: $O(n^2)$ worst case (but $O(n)$ for well-distributed points)
• $d \geq 4$: can be exponential in $n$

Construction: generalize incremental or divide-and-conquer. The Bowyer-Watson algorithm extends naturally: insert point, find cavity (violated simplices), retriangulate.

Practical note: for $d \geq 5$, Delaunay triangulations become impractical due to combinatorial explosion. Alternative spatial data structures (kd-trees, ball trees) are preferred for nearest-neighbor queries in high dimensions.

Computational Considerations

Robustness: the InCircle and Orientation predicates must be exact. Floating-point errors cause:

• Missing or duplicate triangles
• Infinite loops in point location
• Invalid triangulations with crossing edges

Standard tools:

• Triangle (Shewchuk, 1996): the gold-standard 2D Delaunay mesh generator. Uses exact arithmetic, robust predicates, and Ruppert's refinement. Produces quality meshes for FEM
• CGAL: comprehensive computational geometry library with exact arithmetic and certified algorithms
• Qhull: general-dimension convex hull / Delaunay via the lifting transform
• scipy.spatial.Delaunay: Python interface to Qhull

Performance: for $n$ points in 2D, modern implementations triangulate $10^6$-$10^7$ points per second. Parallelization is possible but challenging due to the inherently sequential nature of incremental insertion (GPU Delaunay triangulation exists but is complex).