Delaunay Triangulation and Voronoi Diagrams
Build the Delaunay triangulation of a point set with Bowyer-Watson on an exact in-circle predicate, verify it by brute force, and read the Voronoi diagram off its dual.
Good triangles and bad ones
A surveyor measures the height of the ground at a few hundred points, and you need a surface through them: connect the points into triangles and interpolate linearly on each one. A finite element solver needs a mesh of the region it works on. In both cases a point set has many possible triangulations, and they are not equally good. Long thin triangles (slivers) connect points that are far apart while skipping a closer neighbor. They make the interpolated terrain show ridges that are not there, and they make the numerical error of a finite element solution larger.
The ear clipping of the previous lesson does nothing to avoid slivers: it takes the first ear it finds. The Delaunay triangulation is the standard answer. It is defined by one geometric rule, it is unique except in degenerate configurations, and among all triangulations of a point set it has the largest smallest angle.
Building it needs a second predicate next to orientation: is a point inside the circle through three others? That predicate involves squares of coordinates, so it is the place where the exactness questions of lesson 2 get sharper. This lesson derives it, makes it exact under a stated bound, uses it to build the triangulation with the Bowyer-Watson algorithm, and checks the result by brute force on thousands of degenerate inputs. At the end, the same structure turns into the Voronoi diagram.
The empty circle rule
Every triangle has a circumcircle, the unique circle through its three corners. A triangulation of a point set is Delaunay if no point of the set lies strictly inside the circumcircle of any triangle. Points exactly on a circumcircle are allowed.
The definition, due to Boris Delaunay (1934), has consequences that make it useful. The ones this lesson relies on, for a set of n points in the plane, not all on one line:
- It exists, for every point set.
- It is unique when no four points lie on a common circle. When four or more points are cocircular, any triangulation of them is equally Delaunay, and a tie-breaking rule has to choose.
- It maximizes the minimum angle: among all triangulations of the point set, a Delaunay triangulation has the largest smallest angle (de Berg et al., chapter 9). This is the formal version of "no unnecessary slivers".
- It covers the convex hull, and like every triangulation of the set that uses all its points, it has
2n - 2 - htriangles, where h is the number of points on the boundary of the hull (including points in the middle of hull edges). - Its dual is the Voronoi diagram, the subject of the last section.
Take four points that form a convex quadrilateral a, b, c, d. There are two ways to triangulate it: cut along ac or along bd. The empty circle rule decides between them. If d is inside the circle through a, b, c, the diagonal ac is not Delaunay and bd is; if it is outside, the reverse; if all four points are on one circle, both are fine. A whole triangulation is Delaunay exactly when this test passes for every interior edge and its two triangles, which is what makes local repairs possible.
The in-circle predicate, from the lifting map
The question "is d inside the circle through a, b, c?" looks like it needs the circle's center and radius: a construction with a division. It does not. Lift every point onto the paraboloid z = x^2 + y^2, so that (x, y) becomes (x, y, x^2 + y^2). A circle in the plane lifts to the intersection of the paraboloid with a plane, and the points inside the circle lift to points below that plane. So "inside the circle" becomes "below the plane through three lifted points": an orientation test in three dimensions, whose sign is a determinant.
Translate everything so that d is at the origin first. With adx = a.x - d.x, ady = a.y - d.y and so on, and alift = adx^2 + ady^2, the test is the sign of
| adx ady alift || bdx bdy blift || cdx cdy clift |
which is positive when d is strictly inside the circle through the counter-clockwise triangle abc, negative when outside, and zero when on it. (Translating changes each lift by a linear function of x and y, which column operations cancel. The determinant keeps its value, with smaller numbers in it.) Expanding along the third column gives three products of a lift and a 2 by 2 minor.
The bound. With every coordinate in [-M, M]:
- each difference is at most
2Min absolute value; - each lift is at most
(2M)^2 + (2M)^2 = 8M^2; - each minor, like
bdx * cdy - cdx * bdy, is at most2 * (2M)^2 = 8M^2; - the determinant, three products of a lift and a minor, is at most
3 * 8M^2 * 8M^2 = 192 M^4.
Lifts and minors fit in long long for any M up to 10^9 (8 * 10^18 is below 9.22 * 10^18). The products do not, so they are computed in __int128, a GCC and Clang extension whose largest value is 2^127 - 1, about 1.70 * 10^38. With M = 10^9 the bound is 1.92 * 10^38, which is too large: the orientation bound from lesson 1 does not carry over. With M = 5 * 10^8 it is 1.2 * 10^37, which fits with room to spare. That is this lesson's coordinate bound, and the static_assert repeats the arithmetic at compile time.
#include <print>
struct Point { long long x, y; };
// In-circle needs a smaller bound than orientation. With |x|, |y| <= M, every
// difference is at most 2M, every "lift" dx^2 + dy^2 at most 8M^2, every 2x2 minor
// at most 8M^2, and the determinant (three lift * minor terms) at most 192 M^4.
constexpr long long kMaxCoord = 500'000'000;
using Wide = __int128; // a GCC/Clang extension, not standard C++
constexpr Wide kWideMax = Wide((static_cast<unsigned __int128>(1) << 127) - 1);
static_assert(192 * Wide(kMaxCoord) * kMaxCoord * kMaxCoord * kMaxCoord <= kWideMax,
"in-circle determinant must fit in __int128");
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
// > 0 if d is strictly inside the circle through a, b, c (counter-clockwise),
// < 0 if strictly outside, 0 if on it. Exact for coordinates within kMaxCoord.
int inCircle(Point a, Point b, Point c, Point d) {
long long adx = a.x - d.x, ady = a.y - d.y;
long long bdx = b.x - d.x, bdy = b.y - d.y;
long long cdx = c.x - d.x, cdy = c.y - d.y;
long long aLift = adx * adx + ady * ady;
long long bLift = bdx * bdx + bdy * bdy;
long long cLift = cdx * cdx + cdy * cdy;
Wide det = Wide(aLift) * (bdx * cdy - cdx * bdy)
+ Wide(bLift) * (cdx * ady - adx * cdy)
+ Wide(cLift) * (adx * bdy - bdx * ady);
return (det > 0) - (det < 0);
}
const char* name(int s) { return s > 0 ? "inside" : s < 0 ? "outside" : "on the circle"; }
int main() {
// The circle through these three points has center (0, 0) and radius 5.
Point a{5, 0}, b{0, 5}, c{-5, 0};
std::println("abc turns {}", cross(a, b, c) > 0 ? "left (CCW)" : "right");
for (Point d : {Point{0, 0}, Point{3, 4}, Point{4, 4}, Point{0, -5}, Point{6, 0}})
std::println("({:>2}, {:>2}): {}", d.x, d.y, name(inCircle(a, b, c, d)));
// The same circle moved near the corner of the allowed box: still exact.
auto shift = [](Point p) {
return Point{p.x + kMaxCoord - 5, p.y - kMaxCoord + 5}; };
std::println("shifted (3, 4): {}", name(inCircle(shift(a), shift(b), shift(c),
shift({3, 4}))));
std::println("shifted (4, 4): {}", name(inCircle(shift(a), shift(b), shift(c),
shift({4, 4}))));
}Press Run to compile and run this example.
abc turns left (CCW)
( 0, 0): inside
( 3, 4): on the circle
( 4, 4): outside
( 0, -5): on the circle
( 6, 0): outside
shifted (3, 4): on the circle
shifted (4, 4): outside(3, 4) is exactly on the circle of radius 5 (it is the 3-4-5 triangle), and the predicate says so, also after moving the whole picture to the corner of the allowed box. The predicate returns a sign, never the determinant, for the same reasons as orientation in lesson 1.
Testing the predicate against a different construction
A predicate this important deserves a test that does not share its formula. The reference below takes the long way round: it constructs the circumcenter as an exact fraction, then compares the squared distance from d to the center with the squared radius, everything multiplied by the common denominator so it stays in integers. That is only safe for small coordinates, so the test draws its points from a 13 by 13 grid, which makes cocircular quadruples common. Each case is also translated by a random offset up to the coordinate bound, and the fast predicate must give the same answer there: a circle does not care where it is.
#include <print>
#include <utility>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
// In-circle needs a smaller bound than orientation. With |x|, |y| <= M, every
// difference is at most 2M, every "lift" dx^2 + dy^2 at most 8M^2, every 2x2 minor
// at most 8M^2, and the determinant (three lift * minor terms) at most 192 M^4.
constexpr long long kMaxCoord = 500'000'000;
using Wide = __int128; // a GCC/Clang extension, not standard C++
constexpr Wide kWideMax = Wide((static_cast<unsigned __int128>(1) << 127) - 1);
static_assert(192 * Wide(kMaxCoord) * kMaxCoord * kMaxCoord * kMaxCoord <= kWideMax,
"in-circle determinant must fit in __int128");
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
// > 0 if d is strictly inside the circle through a, b, c (counter-clockwise),
// < 0 if strictly outside, 0 if on it. Exact for coordinates within kMaxCoord.
int inCircle(Point a, Point b, Point c, Point d) {
long long adx = a.x - d.x, ady = a.y - d.y;
long long bdx = b.x - d.x, bdy = b.y - d.y;
long long cdx = c.x - d.x, cdy = c.y - d.y;
long long aLift = adx * adx + ady * ady;
long long bLift = bdx * bdx + bdy * bdy;
long long cLift = cdx * cdx + cdy * cdy;
Wide det = Wide(aLift) * (bdx * cdy - cdx * bdy)
+ Wide(bLift) * (cdx * ady - adx * cdy)
+ Wide(cLift) * (adx * bdy - bdx * ady);
return (det > 0) - (det < 0);
}
// Reference for small coordinates, by a different route: construct the circumcenter
// as a fraction u / D (relative to a) and compare squared distances, all scaled by D.
int referenceInCircle(Point a, Point b, Point c, Point d) {
long long bx = b.x - a.x, by = b.y - a.y, cx = c.x - a.x, cy = c.y - a.y;
long long D = 2 * (bx * cy - by * cx); // > 0 for counter-clockwise abc
long long ux = cy * (bx * bx + by * by) - by * (cx * cx + cy * cy);
long long uy = bx * (cx * cx + cy * cy) - cx * (bx * bx + by * by);
long long px = D * (d.x - a.x) - ux, py = D * (d.y - a.y) - uy;
long long dist2 = px * px + py * py, radius2 = ux * ux + uy * uy; // both times D^2
return (dist2 < radius2) - (dist2 > radius2);
}
int main() {
Lcg rng{42};
const int R = 6; // a tiny range, so that four points are often cocircular
long cases = 0, wrong = 0, tally[3] = {};
while (cases < 1'000'000) {
Point p[4];
for (Point& q : p) q = {rng.range(-R, R), rng.range(-R, R)};
long long turn = cross(p[0], p[1], p[2]);
if (turn == 0) continue; // no circle through collinear points
if (turn < 0) std::swap(p[1], p[2]); // make abc counter-clockwise
int want = referenceInCircle(p[0], p[1], p[2], p[3]);
// Translate all four points far out: the answer must not change.
long long ox = rng.range(-int(kMaxCoord) + R, int(kMaxCoord) - R);
long long oy = rng.range(-int(kMaxCoord) + R, int(kMaxCoord) - R);
Point t[4];
for (int i = 0; i < 4; ++i) t[i] = {p[i].x + ox, p[i].y + oy};
wrong += inCircle(p[0], p[1], p[2], p[3]) != want;
wrong += inCircle(t[0], t[1], t[2], t[3]) != want;
++tally[want + 1];
++cases;
}
std::println("cases {}: outside {}, on the circle {}, inside {}", cases, tally[0],
tally[1], tally[2]);
std::println("disagreements with the reference: {}", wrong);
}Press Run to compile and run this example.
cases 1000000: outside 563913, on the circle 24492, inside 411595
disagreements with the reference: 0Almost 25,000 of the cases are exactly on the circle, and both versions agree on every one. This test covers small configurations anywhere in the box. For large configurations, where every difference is near 2M, a separate check outside the playground compared the predicate with a 4 by 4 determinant evaluated in Python's unbounded integers, on 300,000 quadruples within the bound (uniformly random, crowded into the corners of the box, and taken from a large circle near the corners), 60,422 of them exactly cocircular. It found no disagreement.
Flipping an edge: which diagonal is Delaunay
Before building whole triangulations, look at the smallest case again. Four points form a long flat quadrilateral, and the program checks which diagonal passes the empty circle test and what that does to the smallest angle.
#include <algorithm>
#include <cmath>
#include <numbers>
#include <print>
struct Point { long long x, y; };
// In-circle needs a smaller bound than orientation. With |x|, |y| <= M, every
// difference is at most 2M, every "lift" dx^2 + dy^2 at most 8M^2, every 2x2 minor
// at most 8M^2, and the determinant (three lift * minor terms) at most 192 M^4.
constexpr long long kMaxCoord = 500'000'000;
using Wide = __int128; // a GCC/Clang extension, not standard C++
constexpr Wide kWideMax = Wide((static_cast<unsigned __int128>(1) << 127) - 1);
static_assert(192 * Wide(kMaxCoord) * kMaxCoord * kMaxCoord * kMaxCoord <= kWideMax,
"in-circle determinant must fit in __int128");
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
// > 0 if d is strictly inside the circle through a, b, c (counter-clockwise),
// < 0 if strictly outside, 0 if on it. Exact for coordinates within kMaxCoord.
int inCircle(Point a, Point b, Point c, Point d) {
long long adx = a.x - d.x, ady = a.y - d.y;
long long bdx = b.x - d.x, bdy = b.y - d.y;
long long cdx = c.x - d.x, cdy = c.y - d.y;
long long aLift = adx * adx + ady * ady;
long long bLift = bdx * bdx + bdy * bdy;
long long cLift = cdx * cdx + cdy * cdy;
Wide det = Wide(aLift) * (bdx * cdy - cdx * bdy)
+ Wide(bLift) * (cdx * ady - adx * cdy)
+ Wide(cLift) * (adx * bdy - bdx * ady);
return (det > 0) - (det < 0);
}
// Smallest angle of triangle abc in degrees: a construction, for display only.
double minAngle(Point a, Point b, Point c) {
auto angle = [](Point p, Point q, Point r) { // at q, between q->p and q->r
double ux = double(p.x - q.x), uy = double(p.y - q.y);
double vx = double(r.x - q.x), vy = double(r.y - q.y);
return std::abs(std::atan2(ux * vy - uy * vx, ux * vx + uy * vy));
};
double m = std::min({angle(c, a, b), angle(a, b, c), angle(b, c, a)});
return m * 180 / std::numbers::pi;
}
int main() {
// A convex quadrilateral a, b, c, d (counter-clockwise), long and flat.
Point a{0, 0}, b{6, -1}, c{12, 0}, d{6, 1};
std::println("d in circle(a, b, c): {}", inCircle(a, b, c, d));
std::println("b in circle(c, d, a): {}", inCircle(c, d, a, b));
std::println("diagonal ac: smallest angle {:.1f} degrees",
std::min(minAngle(a, b, c), minAngle(a, c, d)));
std::println("diagonal bd: smallest angle {:.1f} degrees",
std::min(minAngle(a, b, d), minAngle(b, c, d)));
std::println("a in circle(b, c, d): {}", inCircle(b, c, d, a));
}Press Run to compile and run this example.
d in circle(a, b, c): 1
b in circle(c, d, a): 1
diagonal ac: smallest angle 9.5 degrees
diagonal bd: smallest angle 18.9 degrees
a in circle(b, c, d): -1The diagonal ac fails: d is inside the circle through a, b, c, and symmetrically b is inside the circle through c, d, a. Replacing ac by bd is an edge flip, and after it the test passes (a is outside the circle through b, c, d). The smallest angle doubles, as the max-min angle property predicts.
This gives a second algorithm for Delaunay triangulations, due to Charles Lawson (1977): start from any triangulation, and flip every edge that fails the test until none does. Each flip makes the triangulation strictly better in a well-defined sense, so the process terminates, and it needs at most O(n^2) flips. Flips need a data structure that can find the triangle on the other side of an edge in constant time; the half-edge lesson builds one and runs Lawson's algorithm on it.
Bowyer-Watson: dig a hole, fill it with a fan
The algorithm this lesson implements was published independently by Adrian Bowyer and by David Watson, in the same issue of The Computer Journal in 1981. It inserts the points one at a time and keeps the triangulation Delaunay after every insertion:
- Find every triangle whose circumcircle strictly contains the new point
p. These triangles are in conflict withp: they cannot stay. - Remove them. Together they form a hole (the cavity) that contains
p, and the cavity is star-shaped as seen fromp. - Connect
pto every edge on the boundary of the cavity.
The boundary edges are easy to find without any geometry: an edge of a removed triangle is on the boundary if the triangle on its other side was not removed, which in code means its reverse direction does not appear among the removed triangles' edges.
Where to start, and the super-triangle. A triangulation needs a first triangle. The textbook trick is a super-triangle: a huge triangle containing all points, whose three corners are deleted at the end with every triangle that touches them. It is a trap. A finite corner is not infinitely far away, so it can lie inside the circumcircle of a flat triangle along the hull, and then the triangles that should cover that part of the hull are never created. On 1,000 random sets of 5 to 40 integer points in a 100 by 100 box, a super-triangle about ten times the size of the box produced the wrong number of triangles for 192 sets. A super-triangle a thousand times the box size had no failures on those sets, but the circumcircle of three nearly collinear integer points can have a radius of the order of the cube of the box size, so no fixed factor is safe without a proof, and a large factor eats into the coordinate bound.
Ghost triangles avoid the problem instead. Add one symbolic vertex at infinity, and for each edge of the hull a ghost triangle made of that edge and the infinite vertex, standing for the region outside the edge. CGAL's triangulations use the same idea, which they call the infinite vertex. A ghost's "circumcircle" is the open half-plane beyond its hull edge, plus the open hull edge itself (the limit of a circle through two points whose center moves to infinity). A new point outside the hull then conflicts with the ghosts of exactly the hull edges it can see, and the same three steps create both the new real triangles and the new ghosts. No coordinate is invented, so the predicates stay exact and within the bound.
Two degenerate inputs need a decision. If all points are collinear there is no triangle, and the constructor returns an empty triangulation. A duplicate point conflicts with nothing (it is on every circle through it), so insert reports it and skips it.
#include <algorithm>
#include <array>
#include <print>
#include <utility>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
// In-circle needs a smaller bound than orientation. With |x|, |y| <= M, every
// difference is at most 2M, every "lift" dx^2 + dy^2 at most 8M^2, every 2x2 minor
// at most 8M^2, and the determinant (three lift * minor terms) at most 192 M^4.
constexpr long long kMaxCoord = 500'000'000;
using Wide = __int128; // a GCC/Clang extension, not standard C++
constexpr Wide kWideMax = Wide((static_cast<unsigned __int128>(1) << 127) - 1);
static_assert(192 * Wide(kMaxCoord) * kMaxCoord * kMaxCoord * kMaxCoord <= kWideMax,
"in-circle determinant must fit in __int128");
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
// > 0 if d is strictly inside the circle through a, b, c (counter-clockwise),
// < 0 if strictly outside, 0 if on it. Exact for coordinates within kMaxCoord.
int inCircle(Point a, Point b, Point c, Point d) {
long long adx = a.x - d.x, ady = a.y - d.y;
long long bdx = b.x - d.x, bdy = b.y - d.y;
long long cdx = c.x - d.x, cdy = c.y - d.y;
long long aLift = adx * adx + ady * ady;
long long bLift = bdx * bdx + bdy * bdy;
long long cLift = cdx * cdx + cdy * cdy;
Wide det = Wide(aLift) * (bdx * cdy - cdx * bdy)
+ Wide(bLift) * (cdx * ady - adx * cdy)
+ Wide(cLift) * (adx * bdy - bdx * ady);
return (det > 0) - (det < 0);
}
// A triangle (a, b, c), counter-clockwise. c == kInfinite marks a ghost triangle:
// it stands for the outside of the hull edge from b to a.
constexpr int kInfinite = -1;
using Tri = std::array<int, 3>;
struct Delaunay {
const std::vector<Point>& pts;
std::vector<Tri> tris;
long inCircleTests = 0;
// Is p inside the circumcircle of t? For a ghost, "inside" means the open
// half-plane beyond its hull edge, plus the open hull edge itself.
bool conflicts(const Tri& t, Point p) {
++inCircleTests;
Point a = pts[t[0]], b = pts[t[1]];
if (t[2] != kInfinite) return inCircle(a, b, pts[t[2]], p) > 0;
long long side = cross(a, b, p);
if (side != 0) return side > 0;
long long dot = (p.x - a.x) * (b.x - a.x) + (p.y - a.y) * (b.y - a.y);
long long len2 = (b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y);
return dot > 0 && dot < len2; // strictly between a and b
}
// Bowyer-Watson step: remove every triangle whose circle contains p, and
// connect p to the boundary of the hole. Returns false for a duplicate point.
bool insert(int ip) {
std::vector<Tri> keep, cavity;
for (const Tri& t : tris) (conflicts(t, pts[ip]) ? cavity : keep).push_back(t);
if (cavity.empty()) return false;
std::vector<std::pair<int, int>> edges;
for (const Tri& t : cavity)
for (int e = 0; e < 3; ++e) edges.push_back({t[e], t[(e + 1) % 3]});
for (auto [u, v] : edges) {
if (std::ranges::find(edges, std::pair{v, u}) != edges.end()) continue;
// (u, v) is on the hole's boundary: make (u, v, p), infinity last.
if (u == kInfinite) keep.push_back({v, ip, kInfinite});
else if (v == kInfinite) keep.push_back({ip, u, kInfinite});
else keep.push_back({u, v, ip});
}
tris = std::move(keep);
return true;
}
explicit Delaunay(const std::vector<Point>& points) : pts(points) {
const int n = int(pts.size());
int i1 = 1, i2 = 2; // find a first non-degenerate triangle
while (i1 < n && pts[i1].x == pts[0].x && pts[i1].y == pts[0].y) ++i1;
for (i2 = i1 + 1; i2 < n && cross(pts[0], pts[i1], pts[i2]) == 0;) ++i2;
if (i2 >= n) return; // all points on one line: no triangles
if (cross(pts[0], pts[i1], pts[i2]) < 0) std::swap(i1, i2);
tris = {{0, i1, i2}, {i1, 0, kInfinite}, {i2, i1, kInfinite}, {0, i2, kInfinite}};
for (int i = 1; i < n; ++i)
if (i != i1 && i != i2) insert(i);
}
std::vector<Tri> finite() const {
std::vector<Tri> out;
for (const Tri& t : tris)
if (t[2] != kInfinite) out.push_back(t);
return out;
}
};
int main() {
// Six points; the four corners of the square are cocircular.
std::vector<Point> pts{{0, 0}, {6, 0}, {6, 6}, {0, 6}, {3, 1}, {9, 3}};
Delaunay dt(pts);
for (auto [a, b, c] : dt.finite())
std::println("({}, {}) ({}, {}) ({}, {})", pts[a].x, pts[a].y, pts[b].x,
pts[b].y, pts[c].x, pts[c].y);
std::println("{} triangles, {} ghost triangles (hull edges)", dt.finite().size(),
dt.tris.size() - dt.finite().size());
Lcg rng{11};
std::println("{:>6} {:>16} {:>8}", "n", "in-circle tests", "/ n^2");
for (int n : {100, 200, 400, 800}) {
std::vector<Point> cloud;
for (int i = 0; i < n; ++i)
cloud.push_back({rng.range(0, 1'000'000), rng.range(0, 1'000'000)});
Delaunay big(cloud);
std::println("{:>6} {:>16} {:>8.2f}", n, big.inCircleTests,
double(big.inCircleTests) / n / n);
}
}Press Run to compile and run this example.
(0, 0) (6, 0) (3, 1)
(6, 0) (6, 6) (3, 1)
(6, 6) (0, 6) (3, 1)
(0, 6) (0, 0) (3, 1)
(6, 6) (6, 0) (9, 3)
5 triangles, 5 ghost triangles (hull edges)
n in-circle tests / n^2
100 9700 0.97
200 39400 0.98
400 158800 0.99
800 637600 1.00The six points have five hull points, and 2 * 6 - 2 - 5 = 5 triangles. The square's four corners are cocircular, so its two diagonal triangulations would both be Delaunay; inserting (3, 1) inside the square decides the question by connecting it to all four corners.
The counting table shows the cost of this implementation: n^2 in-circle tests, because finding the conflicting triangles scans all of the roughly 2n triangles for every insertion. The cavity itself is small on average, so the scan is the waste. Practical implementations find one conflicting triangle by point location (walking through the triangulation towards p, or keeping a history of destroyed triangles) and grow the cavity from there through neighbors. With the points inserted in random order, Guibas, Knuth and Sharir (Algorithmica, 1992) showed that this takes O(n log n) expected time. The worst case of incremental insertion stays O(n^2), since a single insertion can destroy a linear number of triangles. Algorithms with an O(n log n) worst case exist: divide and conquer (Guibas and Stolfi, 1985) and Fortune's sweep line for the Voronoi diagram (1987).
Checking the triangulation by brute force
Every property from the definition can be checked directly, independently of how the triangulation was built:
- every triangle is counter-clockwise, with positive area;
- empty circles: no input point is strictly inside any triangle's circumcircle, tested for every triangle against every point,
O(n * t); - the number of triangles is
2n - 2 - h, with n counting distinct points and h counted by brute force (a point is on the hull boundary when some line through it and another point has every point on one side); - the triangle areas add up to the area of the convex hull, from the monotone chain of lesson 4.
Together they say that the triangles cover the hull without gaps or overlaps and satisfy the Delaunay rule. The inputs are deliberately hostile: up to 40 points on grids as small as 2 by 2, so that duplicates, collinear triples and cocircular quadruples are everywhere.
#include <algorithm>
#include <array>
#include <print>
#include <set>
#include <utility>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
// In-circle needs a smaller bound than orientation. With |x|, |y| <= M, every
// difference is at most 2M, every "lift" dx^2 + dy^2 at most 8M^2, every 2x2 minor
// at most 8M^2, and the determinant (three lift * minor terms) at most 192 M^4.
constexpr long long kMaxCoord = 500'000'000;
using Wide = __int128; // a GCC/Clang extension, not standard C++
constexpr Wide kWideMax = Wide((static_cast<unsigned __int128>(1) << 127) - 1);
static_assert(192 * Wide(kMaxCoord) * kMaxCoord * kMaxCoord * kMaxCoord <= kWideMax,
"in-circle determinant must fit in __int128");
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
// > 0 if d is strictly inside the circle through a, b, c (counter-clockwise),
// < 0 if strictly outside, 0 if on it. Exact for coordinates within kMaxCoord.
int inCircle(Point a, Point b, Point c, Point d) {
long long adx = a.x - d.x, ady = a.y - d.y;
long long bdx = b.x - d.x, bdy = b.y - d.y;
long long cdx = c.x - d.x, cdy = c.y - d.y;
long long aLift = adx * adx + ady * ady;
long long bLift = bdx * bdx + bdy * bdy;
long long cLift = cdx * cdx + cdy * cdy;
Wide det = Wide(aLift) * (bdx * cdy - cdx * bdy)
+ Wide(bLift) * (cdx * ady - adx * cdy)
+ Wide(cLift) * (adx * bdy - bdx * ady);
return (det > 0) - (det < 0);
}
// A triangle (a, b, c), counter-clockwise. c == kInfinite marks a ghost triangle:
// it stands for the outside of the hull edge from b to a.
constexpr int kInfinite = -1;
using Tri = std::array<int, 3>;
struct Delaunay {
const std::vector<Point>& pts;
std::vector<Tri> tris;
long inCircleTests = 0;
// Is p inside the circumcircle of t? For a ghost, "inside" means the open
// half-plane beyond its hull edge, plus the open hull edge itself.
bool conflicts(const Tri& t, Point p) {
++inCircleTests;
Point a = pts[t[0]], b = pts[t[1]];
if (t[2] != kInfinite) return inCircle(a, b, pts[t[2]], p) > 0;
long long side = cross(a, b, p);
if (side != 0) return side > 0;
long long dot = (p.x - a.x) * (b.x - a.x) + (p.y - a.y) * (b.y - a.y);
long long len2 = (b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y);
return dot > 0 && dot < len2; // strictly between a and b
}
// Bowyer-Watson step: remove every triangle whose circle contains p, and
// connect p to the boundary of the hole. Returns false for a duplicate point.
bool insert(int ip) {
std::vector<Tri> keep, cavity;
for (const Tri& t : tris) (conflicts(t, pts[ip]) ? cavity : keep).push_back(t);
if (cavity.empty()) return false;
std::vector<std::pair<int, int>> edges;
for (const Tri& t : cavity)
for (int e = 0; e < 3; ++e) edges.push_back({t[e], t[(e + 1) % 3]});
for (auto [u, v] : edges) {
if (std::ranges::find(edges, std::pair{v, u}) != edges.end()) continue;
// (u, v) is on the hole's boundary: make (u, v, p), infinity last.
if (u == kInfinite) keep.push_back({v, ip, kInfinite});
else if (v == kInfinite) keep.push_back({ip, u, kInfinite});
else keep.push_back({u, v, ip});
}
tris = std::move(keep);
return true;
}
explicit Delaunay(const std::vector<Point>& points) : pts(points) {
const int n = int(pts.size());
int i1 = 1, i2 = 2; // find a first non-degenerate triangle
while (i1 < n && pts[i1].x == pts[0].x && pts[i1].y == pts[0].y) ++i1;
for (i2 = i1 + 1; i2 < n && cross(pts[0], pts[i1], pts[i2]) == 0;) ++i2;
if (i2 >= n) return; // all points on one line: no triangles
if (cross(pts[0], pts[i1], pts[i2]) < 0) std::swap(i1, i2);
tris = {{0, i1, i2}, {i1, 0, kInfinite}, {i2, i1, kInfinite}, {0, i2, kInfinite}};
for (int i = 1; i < n; ++i)
if (i != i1 && i != i2) insert(i);
}
std::vector<Tri> finite() const {
std::vector<Tri> out;
for (const Tri& t : tris)
if (t[2] != kInfinite) out.push_back(t);
return out;
}
};
// Monotone chain (lesson 4), dropping collinear points: hull vertices, CCW.
std::vector<Point> hull(std::vector<Point> p) {
std::ranges::sort(p, {}, [](Point q) { return std::pair{q.x, q.y}; });
std::vector<Point> h(2 * p.size());
std::size_t k = 0;
for (int pass = 0; pass < 2; ++pass, std::ranges::reverse(p))
for (std::size_t i = 0, start = k; i < p.size(); h[k++] = p[i++])
while (k >= start + 2 && cross(h[k - 2], h[k - 1], p[i]) <= 0) --k;
h.resize(k > 1 ? k - 1 : k);
return h;
}
// Independent checks: orientation, empty circles, triangle count, area.
bool valid(const std::vector<Point>& pts, const std::vector<Tri>& tris) {
std::set<std::pair<long long, long long>> distinct;
for (Point p : pts) distinct.insert({p.x, p.y});
long long n = std::ssize(distinct), onHull = 0; // h counts mid-edge points too
for (auto [x, y] : distinct)
onHull += std::ranges::any_of(pts, [&](Point q) { // a supporting line p-q
return (q.x != x || q.y != y) && std::ranges::all_of(pts, [&](Point r) {
return cross({x, y}, q, r) >= 0; });
});
auto h = hull(pts);
long long hullArea = 0, sum = 0;
for (std::size_t i = 0; i < h.size(); ++i)
hullArea += cross({0, 0}, h[i], h[(i + 1) % h.size()]);
for (auto [a, b, c] : tris) {
if (cross(pts[a], pts[b], pts[c]) <= 0) return false;
sum += cross(pts[a], pts[b], pts[c]);
for (Point p : pts)
if (inCircle(pts[a], pts[b], pts[c], p) > 0) return false;
}
if (hullArea == 0) return tris.empty(); // all points on one line
return std::ssize(tris) == 2 * n - 2 - onHull && sum == hullArea;
}
int main() {
Lcg rng{11};
int tests = 0, passed = 0;
for (; tests < 3000; ++tests) {
int n = rng.range(1, 40), side = rng.range(1, 12);
std::vector<Point> pts;
for (int i = 0; i < n; ++i)
pts.push_back({rng.range(0, side), rng.range(0, side)});
Delaunay dt(pts);
passed += valid(pts, dt.finite());
}
std::println("random point sets: {} of {} passed", passed, tests);
}Press Run to compile and run this example.
random point sets: 3000 of 3000 passedIn a longer background run, 40,000 more point sets, a third of them drawn from the 12 integer points on a circle of radius 5 (where every quadruple is cocircular) mixed with random points, also passed every check, and a pairwise overlap test of the triangles found no overlaps.