Delaunay Triangulation and the Empty Circumcircle Test

A triangulation of a point set is a maximal set of non-crossing straight-line edges among the points, partitioning their convex hull into triangles. For n points with h of them on the hull, every triangulation has the same number of triangles, 2n − 2 − h, and the same number of edges, 3n − 3 − h. The counts follow from Euler’s formula and do not depend on which diagonals are chosen. What does depend on the choice is the shape of the triangles, and the Delaunay triangulation is the one that avoids slivers: among all triangulations, it maximises the smallest angle, and it is the unique triangulation in which no point lies strictly inside any triangle’s circumcircle, whenever no four of the points are cocircular.

That empty-circumcircle property is the operational definition, and it converts a global optimisation over triangulations into a local test. A triangulation is Delaunay exactly when every interior edge is locally Delaunay: for the two triangles sharing the edge, neither opposite vertex lies inside the other’s circumcircle. Local optimality implies global optimality here, which is what makes incremental repair work.

The in-circle determinant

For a triangle with vertices a, b, c in counterclockwise order and a query point d, d lies inside the circumcircle when the following determinant is positive:

| ax-dx  ay-dy  (ax-dx)² + (ay-dy)² |
| bx-dx  by-dy  (bx-dx)² + (by-dy)² |
| cx-dx  cy-dy  (cx-dx)² + (cy-dy)² |

The translation by d is not cosmetic; it reduces the degree of the terms and improves conditioning. Expanded, the determinant is a sum of products of degree four in the coordinate differences, which makes it the most numerically delicate of the standard predicates. In double precision, coordinates around 10^4 produce intermediate magnitudes around 10^16, and the cancellation in the final sum can leave a result whose sign is decided by rounding.

The predicate is used in the exact or adaptive form described in exact geometric predicates, and it is not optional here. A single incorrect in-circle result flips one edge, and a flipped edge can violate the empty-circle property for neighbouring triangles, so the errors propagate through the mesh and produce overlapping triangles that are individually well formed.

Incremental construction by cavity retriangulation

The Bowyer–Watson formulation builds the triangulation one point at a time. Start with a triangle large enough to contain every input point, inserted as scaffolding and removed at the end. For each new point p, find every triangle whose circumcircle contains p. Those triangles form a connected region, the cavity, whose boundary is a star-shaped polygon visible from p. Delete the cavity’s triangles and connect p to each boundary edge. The result is Delaunay, because every new triangle has p as a vertex and its circumcircle is empty by construction.

def insert(p, triangles, in_circle):
    cavity = [t for t in triangles if in_circle(t, p) > 0]
    boundary = edges_on_boundary(cavity)     # edges shared by one cavity triangle
    triangles = [t for t in triangles if t not in cavity]
    triangles += [make_triangle(e, p) for e in boundary]
    return triangles

The cost of one insertion is proportional to the number of triangles in the cavity, which is O(1) on average for points inserted in random order and O(n) in the worst case. Scanning all triangles to find the cavity is O(n) per insertion, giving O(n²) overall, which is the version most implementations start with. Adding a point-location structure, such as a walk from a known triangle toward p across edges or a hierarchy over the triangles, brings the expected total to O(n log n).

The random insertion order is load-bearing for the analysis. Inserting points in a spatially sorted order, which is what happens when input comes from a scan or a file written in coordinate order, makes each new point adjacent to the previous cavity and the total work quadratic. Shuffling the input, or generating a random permutation of indices, restores the expected bound at the cost of one array.

Flipping instead of rebuilding

The alternative to cavity retriangulation is edge flipping. Insert p into the triangle containing it, splitting it into three, then examine the edges opposite p. If the neighbouring triangle’s circumcircle contains p, flip the shared edge, replacing the two triangles’ common diagonal with the other diagonal of the quadrilateral, and continue from the new edges. Each flip strictly increases the minimum angle of the affected quadrilateral, so the process terminates, and the total number of flips over the whole construction is linear in expectation.

Flipping has an advantage that matters in practice: it is local, so it can be applied to an existing triangulation to restore the Delaunay property after a vertex moves. That is the basis of mesh smoothing and of remeshing in finite-element codes, where the point set changes continuously and rebuilding from scratch is wasteful. The price is that the flip loop must be careful about non-convex quadrilaterals, where the alternative diagonal lies outside the union of the two triangles and the flip would create an invalid mesh.

Cocircular points and non-uniqueness

Four points on a common circle make the in-circle determinant exactly zero, and the empty-circle property then holds for both possible diagonals. The Delaunay triangulation is not unique; it is a family, and any member of the family satisfies the definition. An implementation that treats zero as “inside” produces one member, one that treats zero as “outside” produces another, and both are correct.

The practical consequence is that the output is implementation-dependent on grid data, symmetric configurations, and any input with cocircular quadruples. Tests that compare a triangulation against a stored reference will fail on such inputs even when both are valid. Comparing invariants such as triangle count, the empty-circle property, and the absence of crossing edges, rather than exact edge lists, avoids that class of false failure.

Collinear triples are a milder version of the same issue: three points on a line produce a degenerate triangle with zero area. Most implementations remove them before triangulation, since a zero-area triangle has an undefined circumcircle and breaks the predicate’s assumptions.

When required edges override the property

The empty-circumcircle property concerns the point set only. When the input includes edges that must appear in the output, such as polygon boundaries, breaklines in terrain models, or material interfaces in a simulation, the constrained Delaunay triangulation relaxes the property to apply only where it does not conflict with a required edge. The result is still well shaped in the unconstrained regions, and the standard algorithm inserts the constraint edges and then restores local Delaunayhood by flips that do not cross them.

What the empty-circle test buys and what it does not

The mesh quality argument deserves its own accounting, because it is often overstated. Maximising the minimum angle means the smallest angle in the Delaunay triangulation is at least as large as in any other triangulation of the same points. It does not mean the smallest angle is large. Four points forming a thin quadrilateral still produce a thin triangle in the Delaunay triangulation, because no diagonal choice makes both triangles well shaped. The property is comparative, and comparisons between bad options remain bad.

The formal statement, due to Sibson, is that among all triangulations of a point set the Delaunay triangulation maximises the minimum angle. The practical value is that slivers appear only where the point set forces them, so mesh generation can concentrate on placing points rather than on choosing diagonals. The two concerns separate cleanly: Delaunay decides connectivity, and point placement decides shape.

Two limits are worth stating. First, Delaunay maximises the minimum angle but does not bound the maximum angle or guarantee a minimum edge length; a Delaunay mesh of badly distributed points still contains thin triangles, and mesh generation adds Steiner points to fix that. Second, the triangulation is of the convex hull, so a non-convex domain needs either the constrained variant or a post-pass that removes triangles outside the domain, and the removal must be done by flood fill from the exterior rather than by testing triangle centroids, which misclassifies slivers that straddle the boundary.

The dual graph of the triangulation is the Voronoi diagram: each input point’s Voronoi cell is the region closer to it than to any other point, and the cell boundaries are perpendicular bisectors of triangulation edges. Constructing either structure gives the other in linear time, which is why implementations usually build the triangulation and derive the diagram rather than computing both.