Delaunay Triangulation

How Delaunay triangulation works

Give me a set of points scattered on a plane. There are many ways to connect them into triangles — combinatorially, often more than you'd want to enumerate. The Delaunay triangulation is the one specific choice that satisfies the empty circumcircle property: for every triangle in the triangulation, draw its circumscribed circle (the unique circle passing through its three vertices); no other point in the set lies strictly inside that circle.

This single property has remarkable consequences:

• Maximum minimum angle. Among all possible triangulations of the same point set, Delaunay maximizes the smallest interior angle. No other triangulation can produce thinner slivers than Delaunay.
• Closest pair guaranteed. The two closest points in the set are always connected by a Delaunay edge.
• Minimum spanning tree subset. Every edge of the Euclidean minimum spanning tree appears in the Delaunay triangulation.
• Voronoi duality. Two points share a Delaunay edge if and only if their Voronoi cells share a boundary. Either structure gives you the other in linear time.

That last property is the deepest. The Voronoi diagram partitions the plane into cells — each cell is "everywhere closer to point P than to any other point." The Delaunay triangulation connects every pair of points whose cells touch. So Delaunay and Voronoi are two views of the same combinatorial structure.

Construction — randomized incremental

The classic algorithm builds the triangulation by adding points one at a time:

1. Initialize with a "super triangle" — a single huge triangle big enough to contain every input point.
2. For each input point p:
   1. Find the triangle T currently containing p (point location — O(log n) with a history structure, O(√n) with a walk).
   2. Split T into three new triangles by connecting p to T's three vertices.
   3. Edge flip: for each new edge, check whether the opposite vertex of its neighboring triangle lies inside the circumcircle. If so, the edge is illegal — flip it (replace it with the perpendicular edge of the quadrilateral) and recurse on the newly exposed edges. The flips propagate until every triangle's circumcircle is again empty.
3. Remove the super triangle and any triangles that reference its vertices.

Expected runtime O(n log n) with randomized insertion order (Guibas-Knuth-Sharir, 1992). Worst case Θ(n²) on pathological inputs — for instance, all points on a convex curve. Real-world implementations either accept the average case or use Bowyer-Watson with a spatial index for guaranteed log-time point location.

When to use Delaunay triangulation

• Finite-element meshing. FEA simulations need triangles or tetrahedra with bounded aspect ratios — Delaunay's maximum-minimum-angle property is exactly that. Ansys, COMSOL, deal.II all use Delaunay or Delaunay-refinement (Ruppert, Chew).
• Terrain models (TINs). Triangulated Irregular Networks store elevation samples as Delaunay triangles. ArcGIS, QGIS, and every GIS package since the 80s ships Delaunay TIN tools.
• Computer graphics surface reconstruction. Given a 3D point cloud from LIDAR or photogrammetry, ball-pivoting and alpha-shape algorithms use Delaunay tetrahedralization to reconstruct surfaces.
• Pathfinding nav meshes. Game AI navigation meshes are usually Delaunay triangulations of walkable terrain — even angles give better pathfinding heuristics.
• Nearest-neighbor preprocessing. Build the Delaunay triangulation once; nearest-neighbor queries on the Voronoi dual are O(log n) per point.
• Natural-neighbor interpolation. Scientific data interpolation that respects local neighbors uses the Delaunay structure to weight nearby samples.

Delaunay vs Voronoi vs convex hull vs Gabriel graph

	Delaunay triangulation	Voronoi diagram	Convex hull	Gabriel graph
Elements	Triangles (edges between points)	Cells (regions of plane)	Outer polygon	Edges between points
Defining property	Empty circumcircle	Closer to one point than any other	All inputs inside	Empty diameter circle on edge
Construction	O(n log n) expected	O(n log n)	O(n log n)	Subgraph of Delaunay
Relation	Dual of Voronoi	Dual of Delaunay	Subset of Delaunay edges	Subgraph of Delaunay
Key use	Meshing, FEA	NN queries, region maps	Outer envelope	Sparse spanning graph
Generalizes to 3D	Tetrahedralization Θ(n²) worst	3D Voronoi cells	O(n log n)	Direct extension
Bounded angle	Yes — max-min angle property	N/A	N/A	Some slivers possible
Famous algorithm	Bowyer-Watson, Fortune	Fortune's sweep	Graham scan, Quickhull	Derived from Delaunay

Practical rule: Delaunay for meshing (FEA, terrain, nav meshes). Voronoi for nearest-neighbor regions (service areas, biology cell layouts). Convex hull for envelope. Gabriel graph for sparse spanners in wireless network topology.

Pseudo-code — Bowyer-Watson incremental

// Insert points one at a time, maintaining the Delaunay invariant.

triangulate(points):
    T = { SuperTriangle big enough to contain all points }
    for p in random_order(points):
        // Find every triangle whose circumcircle contains p — the "bad" triangles.
        bad = { t in T : circumcircle(t).contains(p) }

        // The boundary of those bad triangles is the polygon we need to retriangulate.
        polygon = boundary_edges(bad)

        // Remove bad triangles; reconnect each boundary edge to p.
        T = T - bad
        for edge in polygon:
            T = T + Triangle(edge.a, edge.b, p)

    // Drop any triangle that uses a super-triangle vertex.
    return [t in T : t doesn't reference super-triangle]

JavaScript implementation

Production work uses Delaunator (Mapbox) which builds Delaunay over 100k points in 100 ms. The teaching implementation below shows the algorithm explicitly:

// Delaunay triangulation via Bowyer-Watson incremental insertion.
// Points are {x, y}; triangles are { a, b, c } (vertex indices).

function inCircumcircle(a, b, c, p) {
  // Returns true if p is strictly inside the circumcircle of triangle abc.
  // Uses the determinant test — exact for floats when |coords| < 2^26 or so.
  const ax = a.x - p.x, ay = a.y - p.y;
  const bx = b.x - p.x, by = b.y - p.y;
  const cx = c.x - p.x, cy = c.y - p.y;
  const det = (ax * ax + ay * ay) * (bx * cy - cx * by)
            - (bx * bx + by * by) * (ax * cy - cx * ay)
            + (cx * cx + cy * cy) * (ax * by - bx * ay);
  return det > 0;
}

function delaunay(points) {
  // 1. Super triangle that encloses every input point.
  let minX = Infinity, minY = Infinity, maxX = -Infinity, maxY = -Infinity;
  for (const p of points) {
    if (p.x < minX) minX = p.x; if (p.y < minY) minY = p.y;
    if (p.x > maxX) maxX = p.x; if (p.y > maxY) maxY = p.y;
  }
  const dx = maxX - minX, dy = maxY - minY;
  const deltaMax = Math.max(dx, dy) * 10;
  const midx = (minX + maxX) / 2, midy = (minY + maxY) / 2;
  const sp = [
    { x: midx - 20 * deltaMax, y: midy - deltaMax },
    { x: midx,                  y: midy + 20 * deltaMax },
    { x: midx + 20 * deltaMax, y: midy - deltaMax },
  ];
  const verts = [...points, ...sp];
  const N = points.length;
  let triangles = [{ a: N, b: N + 1, c: N + 2 }];

  // 2. Insert each input point.
  for (let i = 0; i < N; i++) {
    const p = verts[i];
    const bad = [], badEdges = [];
    for (const t of triangles) {
      if (inCircumcircle(verts[t.a], verts[t.b], verts[t.c], p)) {
        bad.push(t);
        badEdges.push([t.a, t.b], [t.b, t.c], [t.c, t.a]);
      }
    }
    // Find boundary edges — appearing exactly once.
    const seen = new Map();
    for (const [u, v] of badEdges) {
      const key = u < v ? `${u}|${v}` : `${v}|${u}`;
      seen.set(key, (seen.get(key) || 0) + 1);
    }
    const boundary = badEdges.filter(([u, v]) => {
      const key = u < v ? `${u}|${v}` : `${v}|${u}`;
      return seen.get(key) === 1;
    });

    triangles = triangles.filter(t => !bad.includes(t));
    for (const [u, v] of boundary) triangles.push({ a: u, b: v, c: i });
  }

  // 3. Drop triangles that reference super-triangle vertices.
  return triangles.filter(t => t.a < N && t.b < N && t.c < N);
}

// Usage
const pts = Array.from({ length: 100 }, () => ({ x: Math.random() * 100, y: Math.random() * 100 }));
const tris = delaunay(pts);

This sketch is O(n²) because each insert scans all current triangles. Production-grade implementations like Delaunator use a half-edge structure with an spatial walk that achieves O(n log n) expected.

Python implementation

# Production Python: use scipy.spatial.Delaunay (Qhull) — robust and fast.
from scipy.spatial import Delaunay
import numpy as np

points = np.random.rand(100, 2)
tri = Delaunay(points)
# tri.simplices is the (n_triangles, 3) array of point indices.

# For Voronoi (the dual):
from scipy.spatial import Voronoi
vor = Voronoi(points)
# vor.regions, vor.vertices, vor.ridge_points


# Teaching implementation — Bowyer-Watson, O(n^2) but illustrative.
import math
from dataclasses import dataclass
from typing import List

@dataclass
class Pt:
    x: float; y: float

def in_circumcircle(a, b, c, p):
    ax, ay = a.x - p.x, a.y - p.y
    bx, by = b.x - p.x, b.y - p.y
    cx, cy = c.x - p.x, c.y - p.y
    det = ((ax*ax + ay*ay) * (bx*cy - cx*by)
         - (bx*bx + by*by) * (ax*cy - cx*ay)
         + (cx*cx + cy*cy) * (ax*by - bx*ay))
    return det > 0

def delaunay(pts: List[Pt]):
    if not pts:
        return []
    minx = min(p.x for p in pts); maxx = max(p.x for p in pts)
    miny = min(p.y for p in pts); maxy = max(p.y for p in pts)
    dmax = max(maxx - minx, maxy - miny) * 10
    midx = (minx + maxx) / 2; midy = (miny + maxy) / 2
    sp = [Pt(midx - 20*dmax, midy - dmax),
          Pt(midx, midy + 20*dmax),
          Pt(midx + 20*dmax, midy - dmax)]
    verts = list(pts) + sp
    N = len(pts)
    triangles = [(N, N+1, N+2)]

    for i, p in enumerate(pts):
        bad = [t for t in triangles
               if in_circumcircle(verts[t[0]], verts[t[1]], verts[t[2]], p)]
        edges = []
        for t in bad:
            edges += [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]
        edge_count = {}
        for (u, v) in edges:
            k = (min(u,v), max(u,v))
            edge_count[k] = edge_count.get(k, 0) + 1
        boundary = [e for e in edges
                    if edge_count[(min(e), max(e))] == 1]
        triangles = [t for t in triangles if t not in bad]
        for (u, v) in boundary:
            triangles.append((u, v, i))

    return [t for t in triangles if all(v < N for v in t)]

Variants

• Bowyer-Watson (1981). The incremental algorithm above. Conceptually clean; most teaching implementations use it.
• Divide-and-conquer (Lee-Schachter 1980). Sort points, recursively triangulate halves, merge along the seam. Worst-case O(n log n).
• Fortune's sweep (1986). Originally for Voronoi but adaptable to Delaunay. Sweeps a horizontal line down the plane; events trigger structural changes. O(n log n) worst-case.
• Constrained Delaunay. Forces specific edges (e.g., the outline of a polygon) into the triangulation even if they violate the empty-circumcircle property. Used in CAD, GIS for polygon meshing.
• Delaunay refinement (Ruppert 1995, Chew 1989). Iteratively inserts points to enforce minimum angle and maximum area bounds. Standard for FEM mesh generation.
• Weighted (regular) triangulation. Each point has a weight; the empty-circumcircle property is replaced by an empty power-circle test. Generalizes Delaunay; dual of the power diagram.
• 3D Delaunay (tetrahedralization). Same idea, one dimension higher. Worst case is Θ(n²) — n points can produce quadratic tetrahedra. CGAL and TetGen handle million-point datasets in seconds.

Performance — costed claims

• Expected runtime. Randomized incremental construction gives O(n log n) expected — proven by Guibas-Knuth-Sharir 1992. Worst case is Θ(n²) on adversarial inputs like collinear points or convex curves.
• JavaScript benchmarks. Delaunator triangulates 100k random 2D points in ~100 ms on modern hardware. d3-delaunay does the same in ~200 ms with a slightly richer API.
• C++ libraries. CGAL, Triangle, and Geogram triangulate one million points in one to two seconds. CGAL adds robust predicates (exact arithmetic) which protect against floating-point degeneracies at the cost of ~2× slowdown.
• 3D Delaunay. TetGen tetrahedralizes one million random 3D points in 5–10 seconds. Worst-case cubic-grid inputs can push this to 30+ seconds due to the Θ(n²) tetrahedra count.
• Memory. Half-edge structure uses ~12 references per triangle. For one million 2D points (~2 million triangles), memory footprint is ~200 MB including coordinates and topology.
• Numerical robustness. The in-circumcircle test is a 4×4 determinant; with double-precision floats it can return wrong sign on near-cocircular points. Production libraries use either exact arithmetic (CGAL) or adaptive precision (Shewchuk's predicates).

Common bugs and edge cases

• Numerical degeneracies. Four points exactly on a circle make the in-circumcircle predicate ambiguous. Production code uses Shewchuk's robust predicates or exact rational arithmetic.
• Collinear points. Three points on a line don't form a triangle; the algorithm needs to handle this either by perturbation (move a coordinate by epsilon) or by extending the triangulation to a degenerate polyline.
• Duplicate points. Inserting the same point twice triggers the in-circumcircle test on a zero-area triangle — undefined behavior. Always de-duplicate inputs first.
• Super-triangle too small. If the super-triangle doesn't enclose every input point, some inputs fall outside and produce wrong results. Make it generous — 20× the bounding box diagonal is standard.
• Forgetting to drop super-triangle vertices. The final triangulation includes triangles touching the super-triangle; remove them or you'll render imaginary triangles extending to infinity.
• Insertion order matters. Worst-case quadratic behavior on adversarial orders (sorted inputs). Always shuffle inputs before insertion; this is what makes "randomized incremental" the standard.