Delaunay Triangulations and Randomized Incremental Construction

Delaunay Triangulations and Randomized Incremental Construction

Part 16 partitioned the plane by nearest site. This chapter turns that partition inside out. Connect two sites whenever their Voronoi cells share an edge, and the resulting straight-line graph is a Delaunay triangulation—a mesh whose triangles have empty circumcircles. The dual view makes proximity structure combinatorial; the circle view makes it testable; randomized incremental construction makes it efficient.

The goal is not merely to draw triangles. We need to know which triangulations are legal, why local edge flips produce a global optimum, how a shuffled insertion order avoids expensive histories in expectation, and where exact predicates are essential.

What counts as a triangulation?

Let \(P\) be a finite set of planar points, not all collinear. A triangulation of \(P\) is a maximal noncrossing straight-line graph with vertex set \(P\). Its bounded faces are triangles, and its outer boundary is the convex hull of \(P\). Maximality matters: if another noncrossing edge can be inserted, the graph is not yet a triangulation.

If \(h\) points lie on the convex hull, Euler's formula gives every triangulation the same counts:

\[ T=2n-h-2, \qquad E=3n-h-3, \]

where \(T\) is the number of bounded triangular faces and \(E\) the number of edges. Different triangulations therefore differ in geometry, not size. A quadrilateral already demonstrates the choice: either diagonal triangulates it, but one may create a much thinner triangle.

Empty circles and the Voronoi dual

Three noncollinear sites \(a,b,c\) define a unique circumcircle. Triangle \(abc\) is Delaunay when the open disk bounded by that circle contains no point of \(P\). Its circumcenter is then a point equidistant from \(a,b,c\) and no closer to another site: precisely a Voronoi vertex. Conversely, a Voronoi vertex incident to three cells yields an empty-circle Delaunay triangle.

The duality needs one qualification. Unbounded Voronoi edges correspond to convex-hull edges, and a Voronoi vertex incident to more than three cells corresponds to more than three cocircular sites. The abstract dual then has a polygonal face. Choosing noncrossing diagonals inside that face produces a Delaunay triangulation, but the choice is not unique.

Why Delaunay triangles avoid small angles

Among all triangulations of a point set in general position, the Delaunay triangulation lexicographically maximizes the sorted list of triangle angles. In particular, it maximizes the minimum angle. This does not promise well-shaped triangles when the input itself forces skinny geometry, but it avoids unnecessarily skinny triangles caused only by a poor diagonal choice.

Consider two triangles sharing an interior edge \(ac\), with opposite vertices \(b\) and \(d\). If their union is a convex quadrilateral, the edge is locally Delaunay when \(d\) is not inside the circumcircle of \(abc\). Otherwise \(ac\) is illegal; flipping it to \(bd\) improves the local angle vector. Repeatedly flipping illegal edges terminates at a Delaunay triangulation. This is both an algorithm and a proof that a local circle test enforces the global empty-circle structure.

Not every edge can flip. The two incident triangles must form a convex quadrilateral; hull edges have only one incident triangle; constrained edges may be required even when they fail the Delaunay test. A constrained Delaunay triangulation changes visibility and legality definitions, so the unconstrained algorithm below cannot simply be told to “never flip” selected edges and expected to recover all constraints.

The in-circle predicate

For counterclockwise \(a,b,c\), the sign of

\[ \det\begin{pmatrix} a_x&a_y&a_x^2+a_y^2&1\\ b_x&b_y&b_x^2+b_y^2&1\\ c_x&c_y&c_x^2+c_y^2&1\\ d_x&d_y&d_x^2+d_y^2&1 \end{pmatrix} \]

is positive when \(d\) is inside the circumcircle, negative outside, and zero on it. If \(a,b,c\) are clockwise, reverse the sign. That orientation dependency is a common source of silently inverted legality tests.

Floating-point evaluation can fail when points are nearly cocircular or coordinates are large enough that squaring loses significant bits. Use adaptive exact predicates or exact arithmetic for sign decisions. Computing a circle center and comparing rounded distances is less reliable and introduces divisions. Predicates decide combinatorics; approximate coordinates may still be used for display.

def orient(a, b, c):
    return ((b[0] - a[0]) * (c[1] - a[1])
            - (b[1] - a[1]) * (c[0] - a[0]))

def incircle(a, b, c, p):
    # Translate p to the origin; use an adaptive-exact version in production.
    ax, ay = a[0] - p[0], a[1] - p[1]
    bx, by = b[0] - p[0], b[1] - p[1]
    cx, cy = c[0] - p[0], c[1] - p[1]
    det = ((ax*ax + ay*ay) * (bx*cy - by*cx)
         - (bx*bx + by*by) * (ax*cy - ay*cx)
         + (cx*cx + cy*cy) * (ax*by - ay*bx))
    return det if orient(a, b, c) > 0 else -det

def locally_delaunay(a, b, c, d):
    """Edge ac, with b and d on opposite sides, is legal at zero or below."""
    return incircle(a, b, c, d) <= 0

Randomized incremental construction

The Bowyer–Watson form of incremental construction inserts one point at a time. Begin with an artificial supertriangle large enough to contain every input point. For a new point \(p\), locate a triangle whose circumcircle conflicts with \(p\), find all conflicting triangles, delete them, and join \(p\) to every edge on the cavity boundary. After all insertions, delete triangles incident to a supertriangle vertex.

The conflicting triangles form a connected, star-shaped cavity for a valid Delaunay triangulation. An edge belongs to its boundary exactly when it occurs in one conflicting triangle but not two. Representing triangle adjacencies explicitly lets a flood fill discover the cavity from one conflicting seed without scanning the whole mesh.

shuffle(P)
T <- triangulation containing only the supertriangle
for p in P:
    seed <- locate a triangle whose circumcircle contains p
    bad <- flood across adjacent triangles while in_circle(triangle, p) > 0
    boundary <- edges incident to exactly one triangle in bad
    delete bad
    for directed edge (u, v) around boundary:
        add counterclockwise triangle (u, v, p), repairing adjacencies
remove every triangle incident to a supertriangle vertex

Point location is the difference between an elegant method and a quadratic scan. Walking from the previously inserted point works well in practice but has distribution-sensitive guarantees. A history DAG records how deleted triangles were replaced: start at the supertriangle and descend through replacement children containing the query. A conflict graph instead stores relationships between uninserted points and current triangles whose circumcircles contain them; insertion examines and updates only affected conflicts.

With a uniformly random insertion order, backward analysis asks how much work the last inserted point could have caused. In a triangulation of \(i\) points there are only \(O(i)\) edges and triangles, and the newest point has expected constant degree because it is equally likely to be any vertex. Thus the expected structural change per insertion is \(O(1)\). With an appropriate history or conflict structure, locating the affected region contributes expected \(O(\log i)\), giving expected \(O(n\log n)\) total time and linear final output size. The expectation is over the algorithm's permutation, not an assumption that the input coordinates are random. The randomized-algorithms deep dive develops this backward-analysis pattern and explains what “expected” does—and does not—guarantee.

Degeneracies and engineering boundaries

Deduplicate exact input coordinates before construction, or assign stable point identities and explicitly suppress zero-length edges. A repeated point has no new Voronoi cell and should not create a cavity. Collinear inputs have no planar triangulation; return the sorted chain or a documented special result. Nearly collinear triples stress orientation just as nearly cocircular quadruples stress in-circle.

For exact cocircularity, choose a policy. One option accepts either diagonal consistently, using a symbolic perturbation or stable identifier tie-break. Another preserves polygonal Delaunay faces rather than forcing triangles. Never perturb displayed coordinates merely to obtain predicate signs: simulated perturbation should resolve combinatorics without corrupting the geometry.

The supertriangle also needs care. “Very large” floating-point coordinates can overflow squared terms in the in-circle predicate. Construct a moderate enclosing triangle from the input bounding box, tag its vertices as artificial, and remove all incident faces at the end. Verify that every cavity boundary is a simple cycle, every output triangle is counterclockwise, adjacent faces agree on edge ownership, and no ordinary point lies strictly inside an output circumcircle.

Finally, ordinary Delaunay triangulation does not honor polygon boundaries, holes, prescribed segments, minimum-angle thresholds, or maximum-area constraints. Those requirements lead to constrained Delaunay triangulation, segment recovery, refinement with Steiner points, or domain-specific mesh generation. The max–min angle property compares triangulations of the same vertices; it is not a universal mesh-quality certificate.

The route from Part 16 is now complete: nearest-site cells and empty circles yield the same planar structure from dual viewpoints. Part 18 broadens that idea. Instead of pairing Voronoi cells with triangles, it will map points to lines, convert incidences into crossings, and use arrangements and levels to turn geometric questions into combinatorial ones.