Meshes and the Half-Edge Data Structure

Store a triangle mesh so that "what is across this edge?" and "what surrounds this vertex?" take constant time, check it with Euler's formula, and flip edges until the mesh is Delaunay.

A list of triangles cannot tell you who is next door

The usual way to store a triangle mesh is two arrays: the vertex positions, and three vertex indices per triangle. That is what OBJ and PLY files contain and what a GPU draws. It is compact, and it answers exactly one question quickly: "what are the corners of triangle t?"

Mesh algorithms ask other questions. Smoothing needs the neighbors of each vertex. Lesson 7's Delaunay test needs the triangle on the other side of an edge. Hole filling needs to walk along the boundary. With a plain triangle list, the only way to find the triangle across edge a-b is to look at every triangle until one contains b followed by a. This program counts how many triangles that search looks at, for grids of increasing size:

main.cppc++23

#include <array>
#include <print>
#include <vector>

using Triangle = std::array<int, 3>;  // three vertex indices, counter-clockwise

// A k-by-k grid of squares, each cut into two triangles.
std::vector<Triangle> gridMesh(int k) {
    std::vector<Triangle> tris;
    auto id = [k](int i, int j) { return j * (k + 1) + i; };
    for (int j = 0; j < k; ++j)
        for (int i = 0; i < k; ++i) {
            tris.push_back({id(i, j), id(i + 1, j), id(i + 1, j + 1)});
            tris.push_back({id(i, j), id(i + 1, j + 1), id(i, j + 1)});
        }
    return tris;
}

// The triangle on the other side of edge a->b: the one that contains b->a.
// With only a list of triangles, the sole option is to look at every one.
int neighborByScan(const std::vector<Triangle>& tris, int a, int b, long& looked) {
    for (int t = 0; t < static_cast<int>(tris.size()); ++t) {
        ++looked;
        for (int k = 0; k < 3; ++k)
            if (tris[t][k] == b && tris[t][(k + 1) % 3] == a) return t;
    }
    return -1;  // a->b is on the boundary
}

int main() {
    std::println("{:>6} {:>12} {:>22}", "tris", "looked at", "per neighbor query");
    for (int k : {4, 8, 16, 32}) {
        auto tris = gridMesh(k);
        long looked = 0, queries = 0;
        for (const Triangle& t : tris)
            for (int e = 0; e < 3; ++e, ++queries)
                neighborByScan(tris, t[e], t[(e + 1) % 3], looked);
        std::println("{:>6} {:>12} {:>22.1f}", tris.size(), looked,
                     double(looked) / double(queries));
    }
}

Output

Press Run to compile and run this example.

  tris    looked at     per neighbor query
    32         1832                   19.1
   128        26800                   69.8
   512       410336                  267.1
  2048      6425536                 1045.8

The search costs about half the mesh per query, so asking it once for every edge costs time proportional to the square of the mesh size. A 2,048-triangle grid already needs more than six million triangle visits. Real meshes have millions of triangles.

This lesson builds the half-edge data structure, which answers every local question (next edge around a face, the triangle across an edge, the next edge around a vertex) by following one stored index. We build it from a triangle list, check it with a validity checker and Euler's formula, and then use its one essential editing operation, the edge flip, to turn any triangulation into a Delaunay triangulation. That uses the in-circle predicate from lesson 7 and the orientation test from lesson 1.

Every edge twice, with a direction

The idea is to split each undirected edge a-b into two directed half-edges, a->b and b->a. Each half-edge belongs to the face on its left. Because faces are counter-clockwise, as everywhere in this course, the three half-edges of a triangle run around it counter-clockwise, and the two half-edges of an interior edge point in opposite directions, one for each neighboring triangle.

Each half-edge stores four integers:

origin: the vertex it leaves. Its destination is the origin of its twin.
twin: the half-edge of the same edge in the opposite direction.
next: the following half-edge around the same face.
face: the triangle on its left.

Take a square with corners 0 (0,0), 1 (2,0), 2 (2,2), 3 (0,2), cut along the diagonal 0-2 into the triangles 0 1 2 and 0 2 3. Triangle 0 1 2 owns the half-edges 0->1, 1->2, 2->0; triangle 0 2 3 owns 0->2, 2->3, 3->0. The diagonal appears twice, as 2->0 and 0->2, and these two are twins. To get from triangle 0 1 2 to its neighbor, take any of its half-edges, jump to the twin, and read the twin's face. One step, independent of the mesh size.

The four outer edges have a triangle on one side only. We give them twins anyway: an outside half-edge with face = -1, running clockwise around the mesh. That costs a few extra entries and removes every special case from the traversals, because twin always exists and the boundary becomes a loop you can walk with next like any face.

We store all half-edges in one std::vector and refer to them by index, not by pointer. Indices stay valid when the vector grows and reallocates, they can be printed and compared in tests, and they take 4 bytes instead of 8.

Building half-edges from a triangle list

Building the structure has two passes. The first creates three half-edges per triangle, with next linking them around the triangle, and records each directed edge (from, to) in a map. The second finds each half-edge's twin by looking up the reversed pair; a half-edge without a partner is on the boundary and gets an outside twin. A final loop links the outside half-edges into boundary loops: the one that ends at vertex a continues with the one that leaves a.

The program uses std::map, which costs O(log E) per lookup and makes the build O(T log T) for T triangles. A hash map would make it expected linear; that is fine too, as long as nothing iterates the hash map to produce output, since its order is not specified.

Building is also where bad input is caught. If the same directed edge appears twice, then either two triangles with opposite orientation share an edge, or three or more triangles share one edge. Neither can be represented, and the builder reports it instead of producing a structure that quietly lies. It returns std::expected, so the caller has to look.

main.cppc++23

#include <array>
#include <expected>
#include <format>
#include <map>
#include <print>
#include <string>
#include <utility>
#include <vector>

using Triangle = std::array<int, 3>;  // vertex indices, counter-clockwise

struct HalfEdge {
    int origin;  // vertex this half-edge leaves
    int twin;    // the opposite half-edge of the same edge
    int next;    // the following half-edge around the same face
    int face;    // triangle on the left, or -1 for the outside
};

struct Mesh {
    std::vector<HalfEdge> he;
    std::vector<int> out;  // per vertex: one outgoing half-edge, -1 if unused
    int dest(int h) const { return he[he[h].twin].origin; }
};

std::expected<Mesh, std::string> build(int vertexCount,
                                       const std::vector<Triangle>& tris) {
    Mesh m;
    m.out.assign(vertexCount, -1);
    std::map<std::pair<int, int>, int> directed;  // (from, to) -> half-edge
    for (int f = 0; f < static_cast<int>(tris.size()); ++f)
        for (int k = 0; k < 3; ++k) {
            int a = tris[f][k], b = tris[f][(k + 1) % 3];
            int h = 3 * f + k;
            if (a == b)
                return std::unexpected(std::format("triangle {} is degenerate", f));
            if (!directed.emplace(std::pair{a, b}, h).second)
                return std::unexpected(std::format("edge {}->{} is used twice", a, b));
            m.he.push_back({a, -1, 3 * f + (k + 1) % 3, f});
            m.out[a] = h;
        }
    // Pair each half-edge with its twin; an edge with no partner gets an
    // outside half-edge (face -1), so that every twin exists.
    std::map<int, int> boundaryOut;  // vertex -> outside half-edge leaving it
    const int interior = static_cast<int>(m.he.size());
    for (int h = 0; h < interior; ++h) {
        int a = m.he[h].origin, b = m.he[m.he[h].next].origin;
        if (auto it = directed.find({b, a}); it != directed.end()) {
            m.he[h].twin = it->second;
            continue;
        }
        int g = static_cast<int>(m.he.size());
        m.he.push_back({b, h, -1, -1});
        m.he[h].twin = g;
        if (!boundaryOut.emplace(b, g).second)
            return std::unexpected(std::format("vertex {} touches the outside twice", b));
        m.out[b] = g;  // boundary vertices start their ring on the outside
    }
    // Chain the outside half-edges into boundary loops.
    for (int g = interior; g < static_cast<int>(m.he.size()); ++g)
        m.he[g].next = boundaryOut.at(m.dest(g));
    return m;
}

int main() {
    // A 2-by-2 square, cut along the diagonal 0-2.
    //  3 --- 2
    //  |   / |
    //  | /   |
    //  0 --- 1
    auto mesh = build(4, {{0, 1, 2}, {0, 2, 3}});
    if (!mesh) {
        std::println("error: {}", mesh.error());
        return 1;
    }
    std::println(" h  edge  twin  next  face");
    for (int h = 0; h < static_cast<int>(mesh->he.size()); ++h) {
        const HalfEdge& e = mesh->he[h];
        std::println("{:>2}  {}->{}  {:>4}  {:>4}  {:>4}", h, e.origin, mesh->dest(h),
                     e.twin, e.next, e.face);
    }
}

Output

Press Run to compile and run this example.

 h  edge  twin  next  face
 0  0->1     6     1     0
 1  1->2     7     2     0
 2  2->0     3     0     0
 3  0->2     2     4     1
 4  2->3     8     5     1
 5  3->0     9     3     1
 6  1->0     0     9    -1
 7  2->1     1     6    -1
 8  3->2     4     7    -1
 9  0->3     5     8    -1

Read the table against the square. Half-edges 0 to 2 are triangle 0, half-edges 3 to 5 are triangle 1, and 2 and 3 are twins (the diagonal). Half-edges 6 to 9 are the outside ones with face -1; following next from 6 gives 9, 8, 7 and back to 6, the boundary walked clockwise: 1->0, 0->3, 3->2, 2->1.

Walking a mesh: faces, one-rings and boundary loops

Three traversals cover most of what mesh algorithms need, and each is a short loop.

Around a face: follow next until you are back at the start.
Around a vertex (its one-ring, the vertices joined to it by an edge): from a half-edge h leaving v, next(twin(h)) is the next half-edge leaving v in clockwise order. Repeat until you are back. For this the mesh keeps, per vertex, one outgoing half-edge out[v]. On a boundary vertex it is the outside half-edge, so the walk starts at the gap and covers the whole fan.
Along the boundary: start at any half-edge with face == -1 and follow next.

The next program builds a hexagon fan: vertex 0 in the middle, vertices 1 to 6 around it.

main.cppc++23

#include <array>
#include <expected>
#include <format>
#include <map>
#include <print>
#include <string>
#include <utility>
#include <vector>

using Triangle = std::array<int, 3>;  // vertex indices, counter-clockwise

struct HalfEdge {
    int origin;  // vertex this half-edge leaves
    int twin;    // the opposite half-edge of the same edge
    int next;    // the following half-edge around the same face
    int face;    // triangle on the left, or -1 for the outside
};

struct Mesh {
    std::vector<HalfEdge> he;
    std::vector<int> out;  // per vertex: one outgoing half-edge, -1 if unused
    int dest(int h) const { return he[he[h].twin].origin; }
};

std::expected<Mesh, std::string> build(int vertexCount,
                                       const std::vector<Triangle>& tris) {
    Mesh m;
    m.out.assign(vertexCount, -1);
    std::map<std::pair<int, int>, int> directed;  // (from, to) -> half-edge
    for (int f = 0; f < static_cast<int>(tris.size()); ++f)
        for (int k = 0; k < 3; ++k) {
            int a = tris[f][k], b = tris[f][(k + 1) % 3];
            int h = 3 * f + k;
            if (a == b)
                return std::unexpected(std::format("triangle {} is degenerate", f));
            if (!directed.emplace(std::pair{a, b}, h).second)
                return std::unexpected(std::format("edge {}->{} is used twice", a, b));
            m.he.push_back({a, -1, 3 * f + (k + 1) % 3, f});
            m.out[a] = h;
        }
    // Pair each half-edge with its twin; an edge with no partner gets an
    // outside half-edge (face -1), so that every twin exists.
    std::map<int, int> boundaryOut;  // vertex -> outside half-edge leaving it
    const int interior = static_cast<int>(m.he.size());
    for (int h = 0; h < interior; ++h) {
        int a = m.he[h].origin, b = m.he[m.he[h].next].origin;
        if (auto it = directed.find({b, a}); it != directed.end()) {
            m.he[h].twin = it->second;
            continue;
        }
        int g = static_cast<int>(m.he.size());
        m.he.push_back({b, h, -1, -1});
        m.he[h].twin = g;
        if (!boundaryOut.emplace(b, g).second)
            return std::unexpected(std::format("vertex {} touches the outside twice", b));
        m.out[b] = g;  // boundary vertices start their ring on the outside
    }
    // Chain the outside half-edges into boundary loops.
    for (int g = interior; g < static_cast<int>(m.he.size()); ++g)
        m.he[g].next = boundaryOut.at(m.dest(g));
    return m;
}

// The corners of one face, by following next until we are back.
std::vector<int> faceLoop(const Mesh& m, int start) {
    std::vector<int> v;
    int h = start;
    do { v.push_back(m.he[h].origin); h = m.he[h].next; } while (h != start);
    return v;
}

// The neighbors of v, clockwise: next(twin(h)) is the next half-edge
// leaving v after h when you turn clockwise around v.
std::vector<int> oneRing(const Mesh& m, int v) {
    std::vector<int> ring;
    int h = m.out[v];
    do { ring.push_back(m.dest(h)); h = m.he[m.he[h].twin].next; } while (h != m.out[v]);
    return ring;
}

void show(const char* label, const std::vector<int>& v) {
    std::print("{:<28}", label);
    for (int x : v) std::print(" {}", x);
    std::println("");
}

int main() {
    // A hexagon fan: vertex 0 in the middle, 1..6 around it counter-clockwise.
    std::vector<Triangle> tris;
    for (int i = 1; i <= 6; ++i) tris.push_back({0, i, i % 6 + 1});
    auto mesh = build(7, tris);
    if (!mesh) {
        std::println("error: {}", mesh.error());
        return 1;
    }
    show("corners of face 2:", faceLoop(*mesh, 3 * 2));
    show("one-ring of 0 (interior):", oneRing(*mesh, 0));
    show("one-ring of 1 (boundary):", oneRing(*mesh, 1));

    std::vector<bool> seen(mesh->he.size());
    for (int h = 0; h < static_cast<int>(mesh->he.size()); ++h) {
        if (mesh->he[h].face != -1 || seen[h]) continue;
        for (int g = h; !seen[g]; g = mesh->he[g].next) seen[g] = true;
        show("boundary loop:", faceLoop(*mesh, h));
    }
}

Output

Press Run to compile and run this example.

corners of face 2:           0 3 4
one-ring of 0 (interior):    6 5 4 3 2 1
one-ring of 1 (boundary):    6 0 2
boundary loop:               2 1 6 5 4 3

The interior vertex 0 has all six neighbors, listed clockwise. The boundary vertex 1 has three: its two boundary neighbors 6 and 2 and the center. Each traversal step costs O(1), so a one-ring costs O(degree) and a boundary loop O(length), regardless of the mesh size.

Why clockwise? next(twin(h)) turns clockwise; turning counter-clockwise needs the previous half-edge in the face, twin(prev(h)). In a triangle, prev(h) is next(next(h)), but the outside loops are longer than three. Libraries that want both directions cheaply store a prev field as well.

A validity checker and Euler's formula

Every operation on a half-edge mesh relies on invariants. If one breaks, traversals return wrong answers or never stop. So we write a function that checks them all, and call it in tests after every edit:

twin(twin(h)) == h and twin(h) != h.
The destination of h is the origin of next(h), and next(h) has the same face as h.
Every interior face is a cycle of exactly three half-edges.
out[v] leaves v, and one walk around v meets every half-edge that leaves v.

The last check is how the structure detects a non-manifold vertex. A mesh is a manifold (with boundary) when every edge has one or two triangles and the triangles around every vertex form one fan, closed or open. Two cones touching at their tips are not: walking around the shared vertex covers one cone and misses the other.

Euler's formula is a second, independent check. For a mesh with V vertices, E edges and F triangles, the number V - E + F is the Euler characteristic. For a connected, orientable surface with g handles (its genus) and b boundary loops, it equals 2 - 2g - b. A closed surface without handles, such as a tetrahedron or a sphere mesh, gives 2; a flat disk gives 1 (2 if you count the outside as a face, the form of the formula for planar graphs); a cylinder gives 0, and so does a torus. If your mesh code reports something else, the connectivity is wrong.

main.cppc++23

#include <array>
#include <expected>
#include <format>
#include <map>
#include <print>
#include <string>
#include <utility>
#include <vector>

using Triangle = std::array<int, 3>;  // vertex indices, counter-clockwise

struct HalfEdge {
    int origin;  // vertex this half-edge leaves
    int twin;    // the opposite half-edge of the same edge
    int next;    // the following half-edge around the same face
    int face;    // triangle on the left, or -1 for the outside
};

struct Mesh {
    std::vector<HalfEdge> he;
    std::vector<int> out;  // per vertex: one outgoing half-edge, -1 if unused
    int dest(int h) const { return he[he[h].twin].origin; }
};

std::expected<Mesh, std::string> build(int vertexCount,
                                       const std::vector<Triangle>& tris) {
    Mesh m;
    m.out.assign(vertexCount, -1);
    std::map<std::pair<int, int>, int> directed;  // (from, to) -> half-edge
    for (int f = 0; f < static_cast<int>(tris.size()); ++f)
        for (int k = 0; k < 3; ++k) {
            int a = tris[f][k], b = tris[f][(k + 1) % 3];
            int h = 3 * f + k;
            if (a == b)
                return std::unexpected(std::format("triangle {} is degenerate", f));
            if (!directed.emplace(std::pair{a, b}, h).second)
                return std::unexpected(std::format("edge {}->{} is used twice", a, b));
            m.he.push_back({a, -1, 3 * f + (k + 1) % 3, f});
            m.out[a] = h;
        }
    // Pair each half-edge with its twin; an edge with no partner gets an
    // outside half-edge (face -1), so that every twin exists.
    std::map<int, int> boundaryOut;  // vertex -> outside half-edge leaving it
    const int interior = static_cast<int>(m.he.size());
    for (int h = 0; h < interior; ++h) {
        int a = m.he[h].origin, b = m.he[m.he[h].next].origin;
        if (auto it = directed.find({b, a}); it != directed.end()) {
            m.he[h].twin = it->second;
            continue;
        }
        int g = static_cast<int>(m.he.size());
        m.he.push_back({b, h, -1, -1});
        m.he[h].twin = g;
        if (!boundaryOut.emplace(b, g).second)
            return std::unexpected(std::format("vertex {} touches the outside twice", b));
        m.out[b] = g;  // boundary vertices start their ring on the outside
    }
    // Chain the outside half-edges into boundary loops.
    for (int g = interior; g < static_cast<int>(m.he.size()); ++g)
        m.he[g].next = boundaryOut.at(m.dest(g));
    return m;
}

// Every invariant the rest of the code relies on. Returns "" if all hold.
std::string validate(const Mesh& m) {
    const int n = static_cast<int>(m.he.size());
    for (int h = 0; h < n; ++h) {
        const HalfEdge& e = m.he[h];
        if (e.twin == h || m.he[e.twin].twin != h) return std::format("twin of {}", h);
        if (m.he[e.next].origin != m.dest(h)) return std::format("next of {}", h);
        if (m.he[e.next].face != e.face) return std::format("face of {}", h);
        if (e.face >= 0 && m.he[m.he[e.next].next].next != h)
            return std::format("face {} is not a triangle", e.face);
    }
    std::vector<int> degree(m.out.size()), ring(m.out.size());
    for (int h = 0; h < n; ++h) ++degree[m.he[h].origin];
    for (int v = 0; v < static_cast<int>(m.out.size()); ++v) {
        if (m.out[v] < 0) continue;
        if (m.he[m.out[v]].origin != v) return std::format("out[{}] is stale", v);
        int h = m.out[v];
        do {
            ++ring[v];
            h = m.he[m.he[h].twin].next;
        } while (h != m.out[v] && ring[v] <= n);  // the bound stops a broken loop
        // One walk around v must meet every half-edge that leaves v.
        if (ring[v] != degree[v])
            return std::format("walk around vertex {} misses half-edges", v);
    }
    return "";
}

void report(const char* name, int vertexCount, const std::vector<Triangle>& tris) {
    auto m = build(vertexCount, tris);
    std::string error = m ? validate(*m) : m.error();
    if (!error.empty()) return std::println("{:<14} rejected: {}", name, error);
    int V = 0, loops = 0;
    for (int h : m->out) V += h >= 0;
    std::vector<bool> seen(m->he.size());
    for (int h = 0; h < static_cast<int>(m->he.size()); ++h) {
        if (m->he[h].face != -1 || seen[h]) continue;
        ++loops;
        for (int g = h; !seen[g]; g = m->he[g].next) seen[g] = true;
    }
    int E = static_cast<int>(m->he.size()) / 2, F = static_cast<int>(tris.size());
    std::println("{:<14} V={:<3} E={:<3} F={:<3} loops={}  V-E+F={}", name, V, E, F,
                 loops, V - E + F);
}

// A k-by-r grid of squares cut into triangles; columns and rows can wrap around.
std::vector<Triangle> grid(int k, int r, bool wrapX, bool wrapY) {
    int cols = wrapX ? k : k + 1, rows = wrapY ? r : r + 1;
    auto id = [&](int i, int j) { return (j % rows) * cols + i % cols; };
    std::vector<Triangle> t;
    for (int j = 0; j < r; ++j)
        for (int i = 0; i < k; ++i) {
            t.push_back({id(i, j), id(i + 1, j), id(i + 1, j + 1)});
            t.push_back({id(i, j), id(i + 1, j + 1), id(i, j + 1)});
        }
    return t;
}

int main() {
    report("tetrahedron", 4, {{0, 1, 2}, {0, 3, 1}, {1, 3, 2}, {0, 2, 3}});
    report("disk (grid)", 16, grid(3, 3, false, false));
    report("cylinder", 16, grid(4, 3, true, false));
    report("torus", 16, grid(4, 4, true, true));
    report("fin", 5, {{0, 1, 2}, {1, 0, 3}, {0, 1, 4}});
    report("flipped face", 4, {{0, 1, 2}, {0, 3, 2}});
    report("bow-tie", 5, {{0, 1, 2}, {0, 3, 4}});
    report("two tetrahedra", 7, {{0, 1, 2}, {0, 3, 1}, {1, 3, 2}, {0, 2, 3},
                                 {0, 4, 5}, {0, 6, 4}, {4, 6, 5}, {0, 5, 6}});
}

Output

Press Run to compile and run this example.

tetrahedron    V=4   E=6   F=4   loops=0  V-E+F=2
disk (grid)    V=16  E=33  F=18  loops=1  V-E+F=1
cylinder       V=16  E=40  F=24  loops=2  V-E+F=0
torus          V=16  E=48  F=32  loops=0  V-E+F=0
fin            rejected: edge 0->1 is used twice
flipped face   rejected: edge 2->0 is used twice
bow-tie        rejected: vertex 0 touches the outside twice
two tetrahedra rejected: walk around vertex 0 misses half-edges

Each count matches the formula: the tetrahedron is a sphere (2), the disk has one boundary loop (1), the cylinder two (0), the torus has one handle and no boundary (0). Four bad inputs are rejected, three by the builder and one only by the checker. The "fin" has three triangles on one edge. The "flipped face" lists one triangle clockwise. The "bow-tie" is two triangles that share only a vertex, which gives that vertex two separate gaps in its fan. The two tetrahedra glued at a vertex have no boundary at all, so only the walk around vertex 0 finds the problem.

Flipping an edge in constant time

An edge flip takes two triangles a b c and b a d that share the edge a-b, removes that edge, and inserts the other diagonal c-d of the quadrilateral. The result is the triangles d c a and c d b. The vertices, the four outer edges and the region covered do not change; only the diagonal does.

In a half-edge mesh a flip reuses all six half-edges and rewires them. Name them as in the code: h = a->b with h1 = b->c and h2 = c->a, and its twin t = b->a with t1 = a->d and t2 = d->b. After the flip, h becomes d->c and t becomes c->d. The new triangle d c a is the cycle h, h2, t1, and c d b is t, t2, h1. Two half-edges, t1 and h1, change face. Finally, if a or b used h or t as their stored outgoing half-edge, that entry must change, because h and t no longer leave those vertices.

Here is a flip on a flat quadrilateral whose long diagonal is a bad choice: the fourth point lies inside the circumcircle of the other triangle.

main.cppc++23

#include <array>
#include <map>
#include <print>
#include <utility>
#include <vector>

struct Point { long long x, y; };
using Triangle = std::array<int, 3>;

struct HalfEdge { int origin, twin, next, face; };

struct Mesh {
    std::vector<HalfEdge> he;
    std::vector<int> out;  // per vertex: one outgoing half-edge
    int dest(int h) const { return he[he[h].twin].origin; }
    int prev(int h) const { return he[he[h].next].next; }  // triangles only
};

// The builder from earlier, without its error checks (the input is known to
// be a valid triangulation here).
Mesh build(int vertexCount, const std::vector<Triangle>& tris) {
    Mesh m;
    m.out.assign(vertexCount, -1);
    std::map<std::pair<int, int>, int> directed;
    for (int f = 0; f < static_cast<int>(tris.size()); ++f)
        for (int k = 0; k < 3; ++k) {
            int a = tris[f][k], b = tris[f][(k + 1) % 3];
            directed[{a, b}] = 3 * f + k;
            m.he.push_back({a, -1, 3 * f + (k + 1) % 3, f});
            m.out[a] = 3 * f + k;
        }
    std::map<int, int> boundaryOut;
    const int interior = static_cast<int>(m.he.size());
    for (int h = 0; h < interior; ++h) {
        int a = m.he[h].origin, b = m.he[m.he[h].next].origin;
        if (auto it = directed.find({b, a}); it != directed.end()) {
            m.he[h].twin = it->second;
            continue;
        }
        int g = static_cast<int>(m.he.size());
        m.he.push_back({b, h, -1, -1});
        m.he[h].twin = g;
        m.out[b] = boundaryOut[b] = g;
    }
    for (int g = interior; g < static_cast<int>(m.he.size()); ++g)
        m.he[g].next = boundaryOut.at(m.dest(g));
    return m;
}

// Replace edge a-b (half-edge h, triangles a b c and b a d) by edge c-d.
// Every half-edge keeps its index; six of them get new links.
void flip(Mesh& m, int h) {
    int t = m.he[h].twin;
    int h1 = m.he[h].next, h2 = m.he[h1].next;  // b->c, c->a
    int t1 = m.he[t].next, t2 = m.he[t1].next;  // a->d, d->b
    int a = m.he[h].origin, b = m.he[t].origin;
    int c = m.he[h2].origin, d = m.he[t2].origin;
    int f0 = m.he[h].face, f1 = m.he[t].face;
    m.he[h] = {d, t, h2, f0};  // triangle d c a:  d->c, c->a, a->d
    m.he[h2].next = t1;
    m.he[t1].next = h;
    m.he[t1].face = f0;
    m.he[t] = {c, h, t2, f1};  // triangle c d b:  c->d, d->b, b->c
    m.he[t2].next = h1;
    m.he[h1].next = t;
    m.he[h1].face = f1;
    if (m.out[a] == h) m.out[a] = t1;  // h and t no longer leave a and b
    if (m.out[b] == t) m.out[b] = h1;
}

// 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");

// > 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);
}

void printFaces(const Mesh& m) {
    for (int h = 0; h < static_cast<int>(m.he.size()); ++h)
        if (m.he[h].face >= 0 && h < m.he[h].next && h < m.prev(h))  // once per face
            std::print("  face {}: {} {} {}", m.he[h].face, m.he[h].origin,
                       m.he[m.he[h].next].origin, m.he[m.prev(h)].origin);
    std::println("");
}

int main() {
    // A flat quadrilateral, cut along its long diagonal 0-2.
    //        3
    //  0           2
    //        1
    std::vector<Point> p{{0, 2}, {5, 0}, {10, 2}, {5, 4}};
    Mesh m = build(4, {{0, 1, 2}, {0, 2, 3}});
    int h = 2;  // the half-edge 2->0 of face 0
    std::println("edge {}-{}: is 3 inside the circle of 2, 0, 1? {}", m.he[h].origin,
                 m.dest(h), inCircle(p[2], p[0], p[1], p[3]) > 0 ? "yes" : "no");
    printFaces(m);
    flip(m, h);
    std::println("after the flip, half-edge {} is {}->{}", h, m.he[h].origin, m.dest(h));
    printFaces(m);
    std::println("is 0 inside the circle of 3, 1, 2? {}",
                 inCircle(p[3], p[1], p[2], p[0]) > 0 ? "yes" : "no");
}

Output

Press Run to compile and run this example.

edge 2-0: is 3 inside the circle of 2, 0, 1? yes
  face 0: 0 1 2  face 1: 0 2 3
after the flip, half-edge 2 is 3->1
  face 1: 0 1 3  face 0: 1 2 3
is 0 inside the circle of 3, 1, 2? no

Before the flip, vertex 3 lies inside the circle through 2, 0 and 1, which makes edge 0-2 illegal in Delaunay terms. After it, the diagonal is 1-3 and the opposite check fails, so the new edge is legal. Note the third line: half-edge 2 still exists, but now it is the edge 3->1. Flips keep indices and change what they mean, which matters for any code that holds on to an index across an edit.

Lawson's algorithm: flip until Delaunay

A triangulation is Delaunay when no point lies strictly inside the circumcircle of any triangle (lesson 7). Checking that directly costs O(n) per triangle. The key fact that makes flips useful is local: an edge a-b with neighboring triangles a b c and b a d is locally Delaunay if d is not strictly inside the circle through a, b, c. A triangulation in which every interior edge is locally Delaunay is Delaunay everywhere. This is the Delaunay lemma; de Berg et al. prove it in chapter 9, in the form that every legal triangulation is a Delaunay triangulation.

That gives Lawson's flip algorithm (Lawson, 1977): start from any triangulation of the points, push every edge on a stack, and repeatedly pop an edge, flip it if it is illegal, and push the four outer edges of its quadrilateral, the only ones whose status can change. It stops when the stack is empty, and then every edge is legal.

Why it stops, and how soon, is clearest with the lifting map from lesson 7. Lift each point (x, y) to (x, y, x^2 + y^2) on a paraboloid. A point is inside the circle through three others exactly when its lift lies below the plane through their lifts. So an illegal edge is a fold where the lifted surface bends upward, and flipping it replaces two lifted triangles by two lower ones. The surface only ever moves down. A flipped-out edge lies strictly above the surface from then on, so it can never come back. There are at most n(n-1)/2 possible edges, so Lawson's algorithm performs at most that many flips: O(n^2) in the worst case. The order of that bound cannot be improved: Hurtado, Noy and Urrutia (1999) constructed pairs of triangulations of the same points that are Ω(n^2) flips apart, so for at least one of the two, any sequence of flips to the Delaunay triangulation is quadratic. In practice the count is much smaller, as the program shows.

An illegal edge always has a strictly convex quadrilateral (the far vertex is inside the circle and on the other side of the edge), so Lawson's algorithm never makes the bad flip that the pitfall below shows.

The starting triangulation. Any triangulation of the points works. The cheapest to write is a sweep that reuses Andrew's monotone chain from lesson 4: sort the points by x then y, maintain the lower and upper hull chains, and every time a new point pops a vertex off a chain, the popped vertex, its predecessor and the new point form a triangle. Popping only on strict turns keeps collinear hull points connected. The result is valid but full of slivers, a good test for the flips.

Checking the result. Three properties, tested on 3,000 seeded point sets, many with coordinates between 0 and 2 or 0 and 4 so that duplicates, collinear points and cocircular points are common (duplicates are removed first; a triangulation cannot use the same location twice):

Every triangle turns left.
No point is strictly inside any triangle's circumcircle, checked by brute force over all points.
The triangle count is 2n - 2 - h, where h is the number of points on the hull boundary. This follows from Euler's formula: with n - E + T = 1 and 3T + h = 2E (each triangle has three edges, interior edges are counted twice, hull edges once), solve for T.

This program is the longest in the lesson because it carries everything it needs. The new parts are sweepTriangulate, illegal, lawson and isDelaunay. The builder and flip repeat the code above, and inCircle is lesson 7's, exact in __int128 for coordinates up to 5 * 10^8.

main.cppc++23

#include <algorithm>
#include <array>
#include <map>
#include <numeric>
#include <print>
#include <utility>
#include <vector>

struct Point { long long x, y; };
using Triangle = std::array<int, 3>;

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);
}
int orientation(Point a, Point b, Point c) {
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}

// 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");

// > 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);
}

struct Lcg {
    unsigned long long s;
    unsigned next() {
        s = s * 6364136223846793005ULL + 1442695040888963407ULL;
        return unsigned(s >> 33);
    }
    int range(int lo, int hi) {  // inclusive
        return lo + int(next() % unsigned(hi - lo + 1));
    }
};

// Any triangulation will do as a start. This is the monotone chain from
// lesson 4: every point it pops off a hull chain becomes one triangle.
std::vector<Triangle> sweepTriangulate(const std::vector<Point>& p) {
    std::vector<int> order(p.size());
    std::iota(order.begin(), order.end(), 0);
    std::ranges::sort(order, {}, [&](int i) { return std::pair{p[i].x, p[i].y}; });
    std::vector<Triangle> tris;
    std::vector<int> lo, up;  // lower and upper hull chains
    for (int q : order) {
        // Pop strict turns only: collinear hull points must stay connected.
        while (lo.size() > 1 && orientation(p[lo.end()[-2]], p[lo.back()], p[q]) < 0) {
            tris.push_back({lo.end()[-2], q, lo.back()});
            lo.pop_back();
        }
        while (up.size() > 1 && orientation(p[up.end()[-2]], p[up.back()], p[q]) > 0) {
            tris.push_back({up.end()[-2], up.back(), q});
            up.pop_back();
        }
        lo.push_back(q);
        up.push_back(q);
    }
    return tris;
}

struct HalfEdge { int origin, twin, next, face; };

struct Mesh {
    std::vector<HalfEdge> he;
    std::vector<int> out;  // per vertex: one outgoing half-edge
    int dest(int h) const { return he[he[h].twin].origin; }
    int prev(int h) const { return he[he[h].next].next; }  // triangles only
};

// The builder from earlier, without its error checks (the input is known to
// be a valid triangulation here).
Mesh build(int vertexCount, const std::vector<Triangle>& tris) {
    Mesh m;
    m.out.assign(vertexCount, -1);
    std::map<std::pair<int, int>, int> directed;
    for (int f = 0; f < static_cast<int>(tris.size()); ++f)
        for (int k = 0; k < 3; ++k) {
            int a = tris[f][k], b = tris[f][(k + 1) % 3];
            directed[{a, b}] = 3 * f + k;
            m.he.push_back({a, -1, 3 * f + (k + 1) % 3, f});
            m.out[a] = 3 * f + k;
        }
    std::map<int, int> boundaryOut;
    const int interior = static_cast<int>(m.he.size());
    for (int h = 0; h < interior; ++h) {
        int a = m.he[h].origin, b = m.he[m.he[h].next].origin;
        if (auto it = directed.find({b, a}); it != directed.end()) {
            m.he[h].twin = it->second;
            continue;
        }
        int g = static_cast<int>(m.he.size());
        m.he.push_back({b, h, -1, -1});
        m.he[h].twin = g;
        m.out[b] = boundaryOut[b] = g;
    }
    for (int g = interior; g < static_cast<int>(m.he.size()); ++g)
        m.he[g].next = boundaryOut.at(m.dest(g));
    return m;
}

// Replace edge a-b (half-edge h, triangles a b c and b a d) by edge c-d.
// Every half-edge keeps its index; six of them get new links.
void flip(Mesh& m, int h) {
    int t = m.he[h].twin;
    int h1 = m.he[h].next, h2 = m.he[h1].next;  // b->c, c->a
    int t1 = m.he[t].next, t2 = m.he[t1].next;  // a->d, d->b
    int a = m.he[h].origin, b = m.he[t].origin;
    int c = m.he[h2].origin, d = m.he[t2].origin;
    int f0 = m.he[h].face, f1 = m.he[t].face;
    m.he[h] = {d, t, h2, f0};  // triangle d c a:  d->c, c->a, a->d
    m.he[h2].next = t1;
    m.he[t1].next = h;
    m.he[t1].face = f0;
    m.he[t] = {c, h, t2, f1};  // triangle c d b:  c->d, d->b, b->c
    m.he[t2].next = h1;
    m.he[h1].next = t;
    m.he[h1].face = f1;
    if (m.out[a] == h) m.out[a] = t1;  // h and t no longer leave a and b
    if (m.out[b] == t) m.out[b] = h1;
}

// Edge h is illegal if the far corner of the other triangle lies strictly
// inside the circumcircle of h's triangle.
bool illegal(const Mesh& m, const std::vector<Point>& p, int h) {
    int t = m.he[h].twin;
    if (m.he[h].face < 0 || m.he[t].face < 0) return false;  // hull edges stay
    int a = m.he[h].origin, b = m.he[t].origin;
    int c = m.he[m.prev(h)].origin, d = m.he[m.prev(t)].origin;
    return inCircle(p[a], p[b], p[c], p[d]) > 0;
}

// Lawson: flip illegal edges until none is left. Returns the number of flips.
long lawson(Mesh& m, const std::vector<Point>& p) {
    std::vector<int> stack(m.he.size());
    std::iota(stack.begin(), stack.end(), 0);
    long flips = 0;
    while (!stack.empty()) {
        int h = stack.back();
        stack.pop_back();
        if (!illegal(m, p, h)) continue;
        flip(m, h);
        ++flips;
        int t = m.he[h].twin;  // the four outer edges of the quad may now be illegal
        for (int g : {m.he[h].next, m.prev(h), m.he[t].next, m.prev(t)})
            stack.push_back(g);
    }
    return flips;
}

std::vector<Point> randomPoints(Lcg& rng, int n, int range) {
    std::vector<Point> p;
    for (int i = 0; i < n; ++i) p.push_back({rng.range(0, range), rng.range(0, range)});
    auto key = [](Point q) { return std::pair{q.x, q.y}; };
    std::ranges::sort(p, {}, key);
    auto [first, last] = std::ranges::unique(p, {}, key);
    p.erase(first, last);  // a duplicate point would need a zero-area triangle
    return p;
}

// Brute force: every triangle turns left and has an empty circumcircle.
bool isDelaunay(const Mesh& m, const std::vector<Point>& p) {
    for (int h = 0; h < static_cast<int>(m.he.size()); ++h) {
        if (m.he[h].face < 0) continue;
        Point a = p[m.he[h].origin], b = p[m.he[m.he[h].next].origin];
        Point c = p[m.he[m.prev(h)].origin];
        if (orientation(a, b, c) <= 0) return false;
        for (Point q : p)
            if (inCircle(a, b, c, q) > 0) return false;
    }
    return true;
}

int main() {
    Lcg rng{2024};
    std::println("{:>5} {:>7} {:>12}", "n", "flips", "flips per n");
    for (int n : {100, 200, 400, 800}) {
        auto p = randomPoints(rng, n, 1'000'000);
        Mesh m = build(static_cast<int>(p.size()), sweepTriangulate(p));
        long flips = lawson(m, p);
        std::println("{:>5} {:>7} {:>12.2f}", p.size(), flips,
                     double(flips) / double(p.size()));
    }
    int passed = 0;
    const int cases = 3000;
    for (int i = 0; i < cases; ++i) {
        // Tiny coordinate ranges force duplicates, collinear and cocircular points.
        auto p = randomPoints(rng, rng.range(1, 30), std::array{2, 4, 10, 1000}[i % 4]);
        const int n = static_cast<int>(p.size());
        auto tris = sweepTriangulate(p);
        Mesh m = build(n, tris);
        lawson(m, p);
        int t = static_cast<int>(tris.size());
        int hull = static_cast<int>(m.he.size()) - 3 * t;  // outside half-edges
        passed += isDelaunay(m, p) && (t == 0 || t == 2 * n - 2 - hull);
    }
    std::println("random point sets: {}/{} Delaunay", passed, cases);
}

Output

Press Run to compile and run this example.

    n   flips  flips per n
  100     260         2.60
  200     664         3.32
  400    1683         4.21
  800    4010         5.01
random point sets: 3000/3000 Delaunay

All 3,000 cases pass, including the ones where all points are collinear (no triangles) or many points lie on one circle (the Delaunay triangulation is not unique there, and the strict > 0 test accepts any of them). The number of flips per point grows slowly on this input, from 2.6 at 100 points to 5.0 at 800, and stays far below the quadratic worst case.

Polygon Clipping: Sutherland-Hodgman and Boolean Operations

Project: A Small, Robust 2D Geometry Toolkit