Bowyer-Watson Algorithm

Purpose: Build a Delaunay triangulation incrementally, in $O(n\log n)$ expected time with randomized insertion order ($O(n^2)$ worst case).

Delaunay triangulation

A triangulation in which no point lies inside the circumcircle of any triangle. Equivalently, it maximises the minimum angle over all triangulations — it avoids skinny triangles. Its dual is the Voronoi diagram.

Algorithm

1. Start with a super-triangle large enough to contain every input point.
2. Insert points one at a time (in random order). For each new point $p$:
   • find all triangles whose circumcircle contains $p$ — these are the “bad” triangles;
   • their union is a star-shaped cavity; collect its boundary edges (edges belonging to exactly one bad triangle);
   • delete the bad triangles, and create one new triangle from $p$ to each boundary edge.
3. At the end, delete every triangle touching a super-triangle vertex.

Code sketch

for (Point p : points) {
    vector<int> bad;
    for (int i = 0; i < (int)tri.size(); i++)
        if (inCircumcircle(p, tri[i])) bad.push_back(i);
 
    vector<Edge> boundary;
    for (int i : bad)
        for (Edge e : edgesOf(tri[i]))
            if (countIn(bad, e) == 1) boundary.push_back(e);   // not shared
 
    removeAll(tri, bad);
    for (Edge e : boundary) tri.push_back(Triangle{e.a, e.b, p});
}

The key predicate

Everything rests on the in-circle test: is $d$ inside the circumcircle of $a,b,c$ (given counter-clockwise orientation)?

$$\mathrm{inCircle}(a,b,c,d)=\mathrm{sign}\left|\begin{array}{ccc}{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{array}\right|>0.$$

With integer coordinates up to $10^9$ this determinant reaches $10^{45}$ — well beyond long long and beyond double’s 53-bit mantissa. Use __int128 with reduced coordinate ranges, exact rational arithmetic, or adaptive-precision predicates (Shewchuk). Getting this predicate wrong is the single most common cause of a broken Delaunay implementation.

Complexity

• Expected: $O(n\log n)$ with randomized insertion and a point-location structure
• Naive (scan all triangles per insertion): $O(n^2)$
• Worst case with adversarial order: $O(n^2)$
• Space: $O(n)$ — a planar triangulation has at most $2n-5$ triangles

Delaunay construction methods

Method	Time	Notes
Bowyer-Watson	$O(n\log n)$ expected	easiest to understand
Divide and conquer (Guibas-Stolfi)	$O(n\log n)$ worst case	quad-edge structure
Sweep line (Fortune)	$O(n\log n)$	see Fortune’s algorithm
Lift to 3D + convex hull	$O(n\log n)$	lift $(x,y)\to (x,y,x^2+y^2)$; the lower hull is the Delaunay triangulation
Flip algorithm	$O(n^2)$	start from any triangulation, flip illegal edges

The lifting method is the most elegant: Delaunay in 2D is exactly the lower convex hull of the points lifted onto a paraboloid in 3D. If you already have a 3D hull routine, you get Delaunay for free.

What Delaunay gives you

• Voronoi diagram — the dual graph, for free
• Euclidean MST — a subgraph of the Delaunay triangulation, so run Kruskal on $O(n)$ edges instead of $O(n^2)$
• Nearest neighbour graph — also a subgraph
• Closest pair — the shortest Delaunay edge
• Largest empty circle, mesh generation, terrain interpolation

Variants / Use Cases

• Delaunay Triangulation — the topic page
• Fortune’s algorithm — the sweep-line construction of the Voronoi dual
• Constrained Delaunay — force certain edges to appear (for polygon meshing)
• 3D Delaunay / tetrahedralisation — the same algorithm with circumspheres, $O(n^2)$ worst case