Two Philosophies for Triangulating a Point Cloud — Bowyer–Watson and Advancing Front

In 1934 Boris Delaunay pulled a dual out of a diagram named for his teacher, Georgy Voronoi. Join the points into triangles, and keep every other point outside every circumcircle. Ninety years on, mesh generators still run on that rule.

Twelve points from seed 3 come back as 17 triangles. Before reading that count, check which way a diagonal goes when four points make a quad. Drag the top vertex up and down. The diagonal flips the moment it crosses the shared edge.

Watch where vertex T enters the circumcircle of the lower triangle. The diagonal switches from LR to BT on that single move. That move is a Lawson flip, the smallest edit that restores the empty circumcircle locally.

Two philosophies for covering a point cloud#

A point set has many triangulations. CFD dislikes the thin, sharp ones. A flat sliver inflates derivative error. It also worsens the condition number of the linear system.

Delaunay triangulation maximizes the minimum angle among the options. The sharpest corner gets as blunt as the points allow. Positions alone fix the connectivity, so reconnecting the same points later does not reshuffle the edges.

The advancing front runs the job backward. It does not scatter points first. The boundary is the start, and triangles are glued inward one at a time. New points appear only when a triangle needs them. Each step is local greed: take the best element available right now.

The empty circumcircle — the property that defines Delaunay#

One sentence fixes the mesh. No triangle's circumscribed circle contains another vertex in its interior. That is the empty-circumcircle condition. Whether a point $d$ lies inside the circumcircle of triangle $abc$ is one determinant. Order $abc$ counter-clockwise and compute:

Here $a,b,c$ are the triangle's three vertices and $d$ is the point under test. A positive determinant means $d$ is inside the circle, negative means outside, and zero means exactly on it. No division, no square root. Only the sign matters, so the test is relatively robust against floating-point error. How that sign breaks and recovers on a single insert is written up in the empty circumcircle and Bowyer–Watson insertion.

Bowyer–Watson: dig a cavity, then refill it#

Repeating Lawson flips also builds a Delaunay mesh. The straighter path is the incremental insert Adrian Bowyer and David Watson published side by side in the same journal in 1981.

Drop a new point $p$ into a mesh that is already Delaunay, and every triangle whose circumcircle contains $p$ goes invalid. Delete them and a polygonal hole remains, a cavity. Connect each boundary edge of that cavity to $p$ and the hole is filled.

Three steps.

1. Find all triangles (the bad triangles) whose circumcircle contains the new point $p$.
2. Extract the boundary edges of the cavity they form. A boundary edge belongs to exactly one bad triangle.
3. Delete the bad triangles and, for each boundary edge, form a new triangle with $p$.

Start from one giant super-triangle that encloses every input point. After the last insert, strip any triangle still touching a super-triangle vertex. Only the final Delaunay mesh remains.

Click the canvas to drop points. Each insert deletes the bad triangles and refills the cavity.

Turn the circumcircles on. None of them should hold another vertex. Pack the points as tightly as you like. The property still holds. That is what Delaunay means.

Python — building incremental Delaunay from scratch#

The circumcircle test and the insertion kernel, moved to pure Python with no libraries.

import numpy as np
 
def circumcircle(a, b, c):
    """Return the circumcircle center and squared radius. None if collinear."""
    ax, ay = a; bx, by = b; cx, cy = c
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if abs(d) < 1e-12:
        return None
    a2, b2, c2 = ax * ax + ay * ay, bx * bx + by * by, cx * cx + cy * cy
    ux = (a2 * (by - cy) + b2 * (cy - ay) + c2 * (ay - by)) / d
    uy = (a2 * (cx - bx) + b2 * (ax - cx) + c2 * (bx - ax)) / d
    return (ux, uy), (ax - ux) ** 2 + (ay - uy) ** 2
 
def in_circumcircle(p, tri, pts):
    """True if point p lies inside the circumcircle of triangle tri."""
    cc = circumcircle(pts[tri[0]], pts[tri[1]], pts[tri[2]])
    if cc is None:
        return False
    (ux, uy), r2 = cc
    return (p[0] - ux) ** 2 + (p[1] - uy) ** 2 < r2 - 1e-12
 
def bowyer_watson(points):
    """Build a Delaunay triangulation by incremental insertion."""
    pts = list(points)
    xs = [p[0] for p in pts]; ys = [p[1] for p in pts]
    dmax = 20 * max(max(xs) - min(xs), max(ys) - min(ys))
    cx, cy = (min(xs) + max(xs)) / 2, (min(ys) + max(ys)) / 2
    base = len(pts)                                    # first super-triangle vertex index
    pts += [(cx - dmax, cy - dmax), (cx + dmax, cy - dmax), (cx, cy + dmax)]
    tris = [(base, base + 1, base + 2)]
 
    for pi in range(base):                             # insert points one by one
        p = pts[pi]
        bad = [t for t in tris if in_circumcircle(p, t, pts)]
        edge_count = {}                                # cavity boundary = edges seen once
        for t in bad:
            for e in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]:
                key = tuple(sorted(e))
                edge_count[key] = edge_count.get(key, 0) + 1
        boundary = [e for e, n in edge_count.items() if n == 1]
        tris = [t for t in tris if t not in bad]
        tris += [(a, b, pi) for a, b in boundary]      # connect p to each boundary edge
 
    return [t for t in tris if all(i < base for i in t)]  # drop the super-triangle
 
if __name__ == "__main__":
    rng = np.random.default_rng(3)
    P = [tuple(xy) for xy in rng.random((12, 2))]
    T = bowyer_watson(P)
    print(f"{len(P)} points -> {len(T)} triangles")
    print("first triangle vertex indices:", T[0])

The run prints 12 points -> 17 triangles. The triangle count inside the convex hull is roughly $2n - 2 - h$ ($n$ is the number of points, $h$ the number of hull points). The pure insertion kernel fits in about 60 lines.

Put twelve points into that count and $2\times 12-2=22$. The printout is 17, so when the formula is exact the hull carries $22-17=5$ points. The super-triangle scale is 20 times the span of the cloud. The sample is drawn from the unit square, so the span does not exceed 1, and the three vertices sit that scale away from the box center, holding all 12 points inside. The 20 is slack for enclosure. It says nothing about how skinny a triangle is.

The super-triangle in the code has a base of twice the scale and a height of twice the scale, so its area is twice the scale squared. Set the cloud span to 1 and the scale is 20, the area 800. Against the unit square of area 1 that holds the points, the super-triangle is 800 times larger. Stripping every triangle that still touches a super-triangle vertex throws that 800 away. What remains is the 17 triangles built from the 12 points.

The determinant above only reads a sign. The code calls a point inside only when the squared distance is smaller than the squared radius minus 1e-12. A point sitting almost on the circle is pushed outside. Three nearly collinear points hit the same 1e-12 and the center is abandoned. On the cavity side, only the visit count of an edge matters. An edge seen twice is interior to the hole. An edge seen once becomes the base of a new triangle, and the mesh gains one triangle per such edge.

How do you honor the boundary — constraints and the pipe#

Up to here, points were enough. A real CFD shape is not. An airfoil surface or a cylinder wall carries boundary lines that must exist as mesh edges. Pure Delaunay is free to slice across them.

The fix is a constrained Delaunay triangulation. To recover a lost boundary line $AB$, first find the band of triangles that $AB$ crosses, the pipe. Start from a triangle attached to point $A$, walk each neighbor that $AB$ punches through, and stop at $B$. The pipe is complete.

Re-triangulate the interior so $AB$ returns as an edge, by swapping of diagonals or by divide-and-conquer. The boundary survives. Triangles next to it may no longer be perfectly Delaunay. Geometric fidelity is traded against element quality.

The advancing front — stacking elements by greed#

The advancing front technique (AFT) turns the boundary problem around. The boundary is the starting line.

Pick one reference segment $AB$ from the boundary list. Find the point $C$ that forms the best triangle with $AB$. Two kinds of candidate: a point already on the front, or a new interior point $I$. Take the higher quality, build the triangle, and update the front. New edges join the front. Edges already on it are treated as closed and removed.

Quality is a product of two factors.

Here $\alpha$ is the triangle's shape factor (1 for equilateral, toward 0 as it flattens) and $\delta$ is how well the actual edge length matches the requested element size (1 when they agree). The larger $\lambda$, the better the element.

The shape factor codes up like this.

import math
 
def shape_quality(a, b, c):
    """Shape factor alpha: 1 for equilateral, toward 0 for slivers."""
    l1 = math.dist(a, b); l2 = math.dist(b, c); l3 = math.dist(c, a)
    area = abs((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) / 2
    denom = l1 * l1 + l2 * l2 + l3 * l3
    return 4 * math.sqrt(3) * area / denom if denom > 0 else 0.0

When the front empties, the closed boundary has collapsed inward and the mesh is finished. Every step takes a local optimum, so an AFT initial mesh is usually cleaner than a Delaunay one. The gap is widest on strongly directional meshes, boundary layers in particular.

For an equilateral triangle the area is the side squared times $\sqrt{3}/4$, and the sum of the squared sides is three times the side squared. Multiply the area by the code's $4\sqrt{3}$ and divide by that sum, and $4\sqrt{3}\cdot(\sqrt{3}/4)/(3)=1$. Stretch one side and drop the height, and the numerator shrinks while the denominator grows, so $\alpha$ heads toward 0. Even with $\delta$ at 1, the product $\lambda$ shrinks with it. That is why a flat element does not stay on the front.

Speed is the other side of the ledger. Every new element must be tested against the other front edges for intersection. Written naively, that is $O(n)$ per element and $O(n^2)$ overall. Unless a background grid or a quadtree localizes the test, a longer front spends its time on intersection checks instead of elements.

When to choose which#

The two methods split the work more than they compete.

Delaunay returns a unique, fast answer once the points are fixed. It is the right tool for filling a convex hull, interpolating scattered data, and cleaning a mesh after the fact. It will not honor a boundary by itself. That pass is separate.

AFT is stronger on boundary fidelity, on the first mesh it produces, and on anisotropy. It is harder to implement, and the data structures decide the runtime.

Production meshers blend them. AFT chooses the order of placement. The Delaunay criterion chooses the connections. That is the front-Delaunay split.

A cut airfoil edge ends pure insertion#

Twelve points becoming 17 triangles is the story for points that already sit in the domain. If an airfoil surface is sliced by a diagonal, some of those 17 no longer represent the shape. Either re-triangulate that band as a pipe, or leave the cut edge on the front as the reference segment and let the next point $C$ be chosen against it.

Stacking triangles with a large $\lambda = \alpha\delta$ does not bring a cut boundary back. Where to add points so a constrained edge returns is covered in recovering a cut boundary with Steiner points. Get the 60-line insertion kernel right first. Move to that procedure only when the geometry actually loses an edge.