Incremental Delaunay Triangulation

The incremental algorithm builds the Delaunay triangulation of a point set $P$ by adding the points one at a time. After each insertion of a new point $p$, every triangle whose circumcircle contains $p$ is illegal and must be removed; the cavity left behind is re-triangulated by connecting $p$ to each edge on the cavity boundary. When $p$ lies outside the current hull, the convex hull grows and new "phantom" triangles are added between $p$ and the visible hull edges.

A textbook clean-room trick to avoid special cases at the boundary is to start with one huge super-triangle that contains every input point, insert the real points into it, and at the end discard every triangle that still references a super-triangle vertex. The code below follows that recipe and keeps the implementation small by scanning the triangle list linearly on each insertion — the more sophisticated conflict graph of de Berg et al. brings the expected running time down to $O(n \log n)$.

Codes

Python

def solve(xs):
    points = xs[0]
    n = len(points)
    if n < 3:
        return []

    M = 0
    for x, y in points:
        M = max(M, abs(x), abs(y))
    big = 10 * (M + 1)
    p_super = [(-big, -big), (3 * big, -big), (-big, 3 * big)]
    all_pts = list(points) + p_super
    s0, s1, s2 = n, n + 1, n + 2
    triangles = [(s0, s1, s2)]

    def in_circle(a, b, c, d):
        ax, ay = a[0]-d[0], a[1]-d[1]
        bx, by = b[0]-d[0], b[1]-d[1]
        cx, cy = c[0]-d[0], c[1]-d[1]
        det = (ax*(by*(cx*cx+cy*cy) - cy*(bx*bx+by*by))
             - ay*(bx*(cx*cx+cy*cy) - cx*(bx*bx+by*by))
             + (ax*ax+ay*ay)*(bx*cy - by*cx))
        return det > 0

    def orient(a, b, c):
        return ((b[0]-a[0])*(c[1]-a[1]) - (b[1]-a[1])*(c[0]-a[0]))

    def ccw(i, j, k):
        if orient(all_pts[i], all_pts[j], all_pts[k]) >= 0:
            return (i, j, k)
        return (i, k, j)

    for idx in range(n):
        p = all_pts[idx]
        bad = []
        keep = []
        for t in triangles:
            a, b, c = all_pts[t[0]], all_pts[t[1]], all_pts[t[2]]
            if orient(a, b, c) < 0:
                a, b, c = all_pts[t[0]], all_pts[t[2]], all_pts[t[1]]
            if in_circle(a, b, c, p):
                bad.append(t)
            else:
                keep.append(t)

        edges = {}
        for t in bad:
            for e in ((t[0], t[1]), (t[1], t[2]), (t[2], t[0])):
                key = tuple(sorted(e))
                edges[key] = edges.get(key, 0) + 1

        triangles = keep
        for (u, v), count in edges.items():
            if count == 1:
                triangles.append(ccw(u, v, idx))

    result = []
    for t in triangles:
        if t[0] < n and t[1] < n and t[2] < n:
            result.append(sorted(t))
    result.sort()
    return result

C++

#include <array>
#include <map>
#include <vector>
#include <cmath>

using Pt = std::array<double, 2>;

static long double orient(Pt a, Pt b, Pt c) {
    return (long double)(b[0]-a[0])*(c[1]-a[1])
         - (long double)(b[1]-a[1])*(c[0]-a[0]);
}

static bool in_circle(Pt a, Pt b, Pt c, Pt d) {
    if (orient(a, b, c) < 0) std::swap(b, c);
    long double ax = a[0]-d[0], ay = a[1]-d[1];
    long double bx = b[0]-d[0], by = b[1]-d[1];
    long double cx = c[0]-d[0], cy = c[1]-d[1];
    long double det = ax*(by*(cx*cx+cy*cy) - cy*(bx*bx+by*by))
                    - ay*(bx*(cx*cx+cy*cy) - cx*(bx*bx+by*by))
                    + (ax*ax+ay*ay)*(bx*cy - by*cx);
    return det > 0;
}

std::vector<std::array<int,3>> incremental_delaunay(std::vector<Pt> pts) {
    int n = pts.size();
    double M = 1.0;
    for (auto& p : pts) M = std::max({M, std::fabs(p[0]), std::fabs(p[1])});
    double big = 10.0 * (M + 1.0);
    pts.push_back({-big, -big});
    pts.push_back({3*big, -big});
    pts.push_back({-big, 3*big});

    std::vector<std::array<int,3>> tris = {{n, n+1, n+2}};
    for (int idx = 0; idx < n; ++idx) {
        std::vector<std::array<int,3>> keep;
        std::map<std::pair<int,int>, int> edge_count;
        for (auto& t : tris) {
            if (in_circle(pts[t[0]], pts[t[1]], pts[t[2]], pts[idx])) {
                for (int k = 0; k < 3; ++k) {
                    int u = t[k], v = t[(k+1)%3];
                    if (u > v) std::swap(u, v);
                    edge_count[{u, v}]++;
                }
            } else {
                keep.push_back(t);
            }
        }
        tris = keep;
        for (auto& [e, c] : edge_count) {
            if (c == 1) tris.push_back({e.first, e.second, idx});
        }
    }
    std::vector<std::array<int,3>> out;
    for (auto& t : tris)
        if (t[0] < n && t[1] < n && t[2] < n) out.push_back(t);
    return out;
}

Example

Loading Python runner...

Description

Run time analysis

$O(n^2)$ in the worst case for this naive variant: each of the $n$ insertions scans the current triangle list, which has $O(n)$ triangles. The randomized-insertion variant of de Berg et al. with a conflict graph improves the expected running time to $O(n \log n)$.

Space analysis

$O(n)$. The triangulation of $n$ points in the plane has $O(n)$ triangles, and the super-triangle adds three phantom vertices.

Proof of correctness

The super-triangle is chosen large enough that every real input point lies strictly inside it, so the initial triangulation is a valid Delaunay triangulation of the three super-triangle vertices. The insertion step preserves the empty-circle property by Lawson's theorem: removing all triangles whose circumcircles contain the new point $p$ produces a star-shaped cavity, and reconnecting $p$ to each boundary edge of the cavity yields a triangulation whose new triangles are all Delaunay with respect to the current set. Finally, discarding the triangles that touch super-triangle vertices leaves exactly the Delaunay triangulation of the real points. See de Berg et al., Chapter 9, for the full argument.

Extensions

Applications

References

• de Berg, Cheong, van Kreveld, Overmars. Computational Geometry: Algorithms and Applications, 3rd ed., Chapter 9 (Delaunay Triangulations).
• Guibas, Knuth, Sharir. "Randomized incremental construction of Delaunay and Voronoi diagrams." Algorithmica, 1992.
• Bowyer, A. "Computing Dirichlet tessellations." The Computer Journal, 1981.
• Watson, D. F. "Computing the $n$-dimensional Delaunay tessellation with application to Voronoi polytopes." The Computer Journal, 1981.
• cp-algorithms — Delaunay triangulation