Delaunay Triangulation, in depth: the empty-circle rule, robust predicates, Bowyer-Watson and edge flips, SciPy in practice and the super-triangle trap

What makes a triangulation Delaunay

A triangulation of a point set P connects the points into triangles that cover the convex hull of P without overlapping. It is Delaunay when the circumcircle of every triangle, the unique circle through its three corners, contains no point of P in its interior. Points on the circle are allowed, which is where the degenerate cases come from.

Three facts explain why this definition matters. First, among all triangulations of P, the Delaunay one maximises the smallest angle, so it is the natural choice for meshes and interpolation, where slivers cause numerical trouble. Second, it is the dual of the Voronoi diagram: connect two points whenever their Voronoi cells share an edge and you get the Delaunay triangulation, with each Voronoi vertex at a triangle's circumcentre. Third, lift every point (x, y) to (x, y, x^2 + y^2) on a paraboloid; the downward-facing faces of the 3-D convex hull of the lifted points project back to exactly the Delaunay triangles. That lifting trick is how Qhull, and therefore SciPy, computes it.

The two predicates

Every Delaunay algorithm reduces to two sign tests. The orientation test says whether three points turn left, turn right or are collinear. The incircle test says whether a fourth point lies inside, on or outside the circle through the first three. Both are determinants, and the sign convention must be stated: the incircle formula below is positive when d is inside the circle only if a, b, c are in counter-clockwise order. Feed it a clockwise triangle and the sign flips, which is the most common bug in hand-written triangulators.

def orient(a, b, c):
    """> 0 if a, b, c turn counter-clockwise, < 0 if clockwise, 0 if collinear."""
    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])

def incircle(a, b, c, d):
    """> 0 if d is strictly inside the circumcircle of CCW triangle a, b, c."""
    adx, ady = a[0] - d[0], a[1] - d[1]
    bdx, bdy = b[0] - d[0], b[1] - d[1]
    cdx, cdy = c[0] - d[0], c[1] - d[1]
    ad, bd, cd = adx*adx + ady*ady, bdx*bdx + bdy*bdy, cdx*cdx + cdy*cdy
    return (adx * (bdy * cd - bd * cdy)
            - ady * (bdx * cd - bd * cdx)
            + ad * (bdx * cdy - bdy * cdx))

In floating point, both tests can return the wrong sign when the true value is near zero, and a triangulator that receives inconsistent answers can loop forever or produce overlapping triangles. Production code uses adaptive exact predicates, described in the robust orientation predicate article; the same technique extends to incircle.

Worked example: five points by hand

Take five points: p0 = (0, 0), p1 = (4, 0), p2 = (4, 3), p3 = (0, 3) at the corners of a 4 by 3 rectangle, and p4 = (2, 1) inside it. The Delaunay triangulation is the four triangles that join p4 to each side: (p0, p1, p4), (p1, p2, p4), (p2, p3, p4) and (p3, p0, p4).

Check one. The circle through p0, p1 and p4 has its centre on x = 2, the perpendicular bisector of p0 p1. Equating the distance to p0 and to p4 gives 4 + y^2 = (y - 1)^2, so y = -1.5 and the radius is 2.5. Point p2 is at distance sqrt(4 + 20.25), about 4.9, from the centre, and p3 is the same distance away, so the circle is empty and the triangle is Delaunay. The counts also check: n = 5 points with h = 4 on the hull should give 2n - h - 2 = 4 triangles and 3n - h - 3 = 8 edges, the four sides plus four spokes.

Now delete p4. The four rectangle corners lie on one circle, so both diagonals give a valid Delaunay triangulation and the answer is no longer unique. That is the cocircular degeneracy, and grids, pixel centres and CAD inputs are full of it. Libraries break the tie with a consistent rule; your code must not assume a particular diagonal.

Building it: Bowyer-Watson

Bowyer-Watson is the easiest correct algorithm to write. Start with a huge super-triangle that contains every point. Insert points one at a time. For each new point, find every triangle whose circumcircle contains it, the bad triangles; together they form a star-shaped cavity around the point. Delete them, then connect the new point to each edge on the cavity boundary. Finally discard any triangle that uses a super-triangle vertex.

def bowyer_watson(points):
    pts = list(points)
    xs, ys = [p[0] for p in pts], [p[1] for p in pts]
    cx, cy = (min(xs) + max(xs)) / 2, (min(ys) + max(ys)) / 2
    r = max(max(xs) - min(xs), max(ys) - min(ys)) * 1000 + 1   # see failure modes
    n = len(pts)
    pts += [(cx - r, cy - r), (cx + r, cy - r), (cx, cy + r)]   # CCW super-triangle
    tris = {(n, n + 1, n + 2)}
    for i in range(n):
        p = pts[i]
        bad = [t for t in tris if incircle(pts[t[0]], pts[t[1]], pts[t[2]], p) > 0]
        edges = {}
        for t in bad:
            for e in ((t[0], t[1]), (t[1], t[2]), (t[2], t[0])):
                key = frozenset(e)
                edges[key] = None if key in edges else e     # shared edges cancel
        tris.difference_update(bad)
        for e in edges.values():
            if e is not None:                                 # cavity boundary, CCW
                tris.add((e[0], e[1], i))
    return [t for t in tris if max(t) < n]                    # drop super vertices

Because every triangle is stored counter-clockwise and boundary edges keep their direction, the new triangles are counter-clockwise too, which keeps the incircle sign valid. On the five-point example this returns the same four triangles as SciPy. As written it scans every triangle per insertion, so it is O(n^2); it is a reference implementation for tests, not for a million points.

Edge flips and fast incremental insertion

The other classical view is local. Call an interior edge legal if the point opposite it in one adjacent triangle lies outside the circumcircle of the other. Lawson's theorem says a triangulation is Delaunay exactly when every edge is legal, and that repeatedly flipping illegal edges, replacing the diagonal of their convex quadrilateral with the other diagonal, always terminates at the Delaunay triangulation. Each flip increases the sorted vector of angles, so it cannot cycle.

Fast incremental algorithms combine both views. Locate the triangle containing the new point, split it into three, then push the three outer edges onto a stack and flip until the stack is empty. With random insertion order and a point-location structure the expected time is O(n log n); in practice libraries sort points along a space-filling curve and locate each one by walking from the previous triangle, which is close to linear on real data.

Counts and complexity

For n points with h of them on the convex hull and no three collinear on the boundary, every triangulation has exactly 2n - h - 2 triangles and 3n - h - 3 edges, by Euler's formula. The average vertex degree is therefore just under six. On 1,000 uniform random points SciPy produced 1,980 triangles with h = 18, and a mean degree of 5.96, which matches. These identities make cheap invariants for tests and for monitoring a pipeline: a triangle count that does not match tells you points were dropped or duplicated.

Using SciPy in practice

In Python, use scipy.spatial.Delaunay, a wrapper around Qhull. It returns triangle vertex indices, neighbours, a point-location routine, barycentric transforms and the vertex adjacency graph, which covers most uses: piecewise linear interpolation of scattered data, nearest-neighbour graphs for clustering or mesh processing, and terrain models in GIS.

import numpy as np
from scipy.spatial import Delaunay
from scipy.interpolate import LinearNDInterpolator

rng = np.random.default_rng(0)
pts = rng.random((1000, 2))
tri = Delaunay(pts)
h = len(tri.convex_hull)                        # hull edges = hull vertices in 2-D
assert len(tri.simplices) == 2 * len(pts) - h - 2

q = np.array([[0.5, 0.5], [1.5, 1.5]])
s = tri.find_simplex(q)                         # -1 means outside the hull

z = np.sin(pts[:, 0] * 6) + pts[:, 1]           # values measured at the points
T = tri.transform[s[0]]
b = T[:2] @ (q[0] - T[2])                       # barycentric coordinates
bary = np.append(b, 1 - b.sum())
print(bary @ z[tri.simplices[s[0]]])            # same as LinearNDInterpolator(tri, z)(q)[0]

indptr, nbrs = tri.vertex_neighbor_vertices     # CSR adjacency of the triangulation
print(nbrs[indptr[0]:indptr[1]])                # neighbours of point 0

Two details bite in practice. Points outside the hull get index -1 and interpolation returns NaN, so decide whether to extrapolate or reject. Duplicate points are dropped and reported in tri.coplanar rather than triangulated, so a downstream array indexed by triangle vertex can silently miss rows. For large inputs, incremental=True allows adding points later, and kd-trees as described in the k-d tree deep dive are the better tool when you only need nearest neighbours.

Constrained and 3-D variants

Real meshes have boundaries and holes. A constrained Delaunay triangulation forces given segments, such as a coastline or a part outline, to appear as edges and relaxes the empty-circle rule only where a constraint blocks visibility. Quality mesh generators such as Shewchuk's Triangle then insert extra Steiner points at circumcentres of poor triangles until every angle exceeds a bound, which is what finite element solvers need. In three dimensions the definition carries over with empty circumspheres, but the guarantees weaken: Delaunay tetrahedralisations can contain flat sliver tetrahedra, and the number of tetrahedra can grow quadratically in the worst case, so 3-D meshing pipelines add sliver removal.

Failure modes

• Super-triangle too small. With the scale factor at 20 instead of 1000, the code above lost hull triangles on 10 of 20 random point sets, because a super vertex sat inside the circumcircle of a thin hull triangle. Larger factors help but cost floating-point precision; robust implementations treat the extra vertices symbolically as points at infinity. Always assert the 2n - h - 2 count.
• Clockwise triangles. A single clockwise triangle flips the incircle sign and corrupts everything inserted after it.
• Floating-point predicates. Near-cocircular or near-collinear inputs give inconsistent signs; use exact or adaptive predicates.
• Duplicates and collinear inputs. All points on one line have no triangulation at all; duplicates are dropped. Deduplicate and check rank first.
• Assuming a unique answer. Cocircular points admit several Delaunay triangulations; tests that compare exact triangle lists across libraries will flake.
• Bad coordinate scaling. Geographic coordinates in degrees with metre-level detail, or large offsets like UTM eastings, lose precision in the squared terms of incircle. Translate to a local origin first.

Trade-offs

Use SciPy or Qhull for analysis and interpolation in Python; CGAL when you need exact predicates, constraints and 3-D in C++; Triangle for 2-D quality meshes; and your own Bowyer-Watson only as a test oracle or teaching tool. If what you actually need is the Voronoi cells, compute the Delaunay triangulation and take the dual. If you only need neighbour queries, a spatial index is simpler and faster to update.

What to do next

1. Run the predicates and Bowyer-Watson code above on the five-point example and confirm the four triangles.
2. Add the count invariant 2n - h - 2 as an assertion wherever you triangulate.
3. Triangulate your real data with SciPy and inspect tri.coplanar for dropped duplicates.
4. Translate coordinates to a local origin before triangulating large-offset data.
5. Build a linear interpolator from the triangulation and compare it with held-out measurements.
6. If you need fixed boundaries or quality guarantees, move to a constrained triangulator such as Triangle or CGAL.
7. Read the Voronoi article next to see the dual structure and Fortune's sweep.