Constrained Delaunay triangulation

Most triangulators want clean input: closed rings that never touch and constraint edges that never cross. Real outlines are not clean. Slices of a scanned part overlap. GIS layers share edges and cross at their seams. A drawing has a hatch line running through a hole. trueform takes that input as it comes, splits every crossing with exact arithmetic, and returns a constrained Delaunay triangulation (CDT) of what you drew. Move holes through a shape in the live example.

What a CDT is

A Delaunay triangulation joins 2D points into triangles so that no point lies inside any triangle's circumcircle. That rule avoids thin slivers as far as the points allow. It covers the convex hull and picks its own edges.

So a Delaunay triangulation ignores your lines. Draw a breakline through a point set, and the edges it picks run straight across it.

A constrained Delaunay triangulation keeps the edges you choose: an outline, the rim of a hole, a breakline, a fold. Every other edge is as Delaunay as those edges allow. The rule becomes: no point visible from inside a triangle lies in its circumcircle, and a constraint blocks the view. On the right, the shaded triangle's circle holds three points. All three sit behind the pink line, so the triangle stays. With outlines and holes, a CDT also decides which triangles are inside.

Readers come to it from many sides: polygon triangulation with holes, PSLG meshing, a base mesh for finite elements, a TIN with breaklines, caps for cut surfaces, letters and logos. Inside trueform, the same CDT cuts the surfaces of mesh booleans and re-triangulates the faces cut by isobands.

What happens to messy input

The input above is synthetic and built to be awkward. It has 1,186 points and 483 constraint edges. Here is what one call does with each part of it.

Crossing constraints are split. Two hole rims cross each other. Three feature lines cross each other, both rims, the third hole and the outline. Before the triangulation, every crossing between constraints is found: edge through edge, a vertex on an edge, two vertices at one spot. A point is added at each one, orange in the figure, and the constraints are split there. You do not pre-process anything.

Walls decide what is inside. By default every constraint is a wall. Walking in from outside, each wall you cross flips inside to outside. Triangles behind an odd number of walls are kept. A ring inside a ring is a hole. Where the two holes overlap, the walk crosses three walls, so the lens is solid again, in sand above. That is the even-odd rule. Each copy of a wall counts, so a rim given twice cancels; give a wall once to keep it. If you mean the union of the two holes, merge them before you triangulate.

Feature lines split nothing. The white lines are constraints with their wall flag off. They appear as edges of the output but leave inside and outside alone. Where they run outside the outline, they are dropped with the outside.

Loose points outside are dropped. Points in the holes and beyond the outline fall away with their triangles. An index map gives the output index of every input point, says which input points fell away, and marks the output points that were made at crossings.

You can refuse crossings instead. When your data must not cross, say so. The call then keeps every constraint whole and returns empty arrays on the first crossing. The empty result is the answer to "do these constraints cross", not an error. It makes the CDT a validity check for outlines.

faces, points = tf.cdt(points2d, edges, split_constraints=False)
valid = len(faces) > 0      # empty: two constraints cross

Regions get labels. Instead of the interior alone, ask for every triangle with a label: its nesting parity, or the id of the region the walls cut it into. Each hole, the lens and the solid then has its own id.

In short:

Exact, and what exact means here

A triangulator rests on a few yes-or-no tests. Is this point left of that line? Is it inside that circle? Do these two segments cross? Near-collinear and near-cocircular points make floating-point answers inconsistent, and inconsistent answers give crossed edges, holes in the mesh or a loop that never ends.

trueform welds the input onto an integer grid over its bounding box: 32-bit for float32 input, 64-bit for float64. On that grid, the orientation, in-circle and crossing tests are exact integer arithmetic. Whether a crossing lies inside a segment is decided on the exact ratio, not on a rounded value. A crossing point is then placed on the same grid and written out in your float type.

Two consequences follow. Input points that weld to the same grid point become one output point. And a crossing point in the output is rounded, because the output is a float array; no float array can hold an exact rational point. The decisions that build the mesh are exact, so the mesh is valid.

In a pipeline

• A quality floor. The CDT adds points only at crossings, so thin triangles along a long constraint stay thin. In C++, refine the result with tf::cdt_refiner, Ruppert's method, which splits encroached constraints at their midpoints. It takes constraints that do not cross.
• A planar cut in 3D. Points are (x, y). Project the outline onto its plane, triangulate, and lift the result back.

# outline: (N, 3) points on a plane through o, with unit in-plane axes u and v
uv = np.c_[(outline - o) @ u, (outline - o) @ v].astype(np.float32)
faces, points2d = tf.cdt(uv, edges)
points = o + points2d[:, :1] * u + points2d[:, 1:] * v     # back in 3D

How it compares

Libraries split into three groups on crossing constraints: those that resolve them, those that refuse them, and those that take points only. One more group is the dangerous one: libraries that accept crossing input and quietly lose a constraint. This table is from trueform's own benchmark harness, which declares each library's capabilities before it runs.

Refusing is correct behaviour; a library that refuses tells you so. The cdt2d case is the one to watch: the unit square with both diagonals came back as two triangles, one diagonal gone.

Benchmarks

Points in, triangles out. Apple M4 Max, 16 threads, trueform 0.10.10.

In Python, trueform 0.10.10 triangulates 1M points in 18.1 ms, against 427 ms for Triangle (24×) and 3,052 ms for SciPy's Qhull (169×). On a 100,000-vertex closed boundary it takes 23.9 ms, against Triangle's 35.8 ms (1.5×); at 200,000, 47.0 ms against 79.0 ms (1.7×); at 10,000 the margin is 1.3×. In Node.js, trueform npm 0.10.10 triangulates 1M points in 30.2 ms, against 456 ms for Delaunator (15×); at 10,000 points the margin is 1.1×. cdt2d tops out at 10,000 points.

Language

C++PythonTypeScript

C++1M points in 18.0 ms

Python1M points in 18.1 ms

TypeScript1M points in 30.2 ms

Every trueform benchmark, by module and language, is on the trueform benchmarks page.

Constrained Delaunay triangulation in code

Each example triangulates a square outline with a square hole and 500 loose points. The interior points inside the hole are dropped, and the area is the square minus the hole. The tutorials cover plain Delaunay, crossings, feature lines, region labels and index maps.

#include <trueform/trueform.hpp>

std::vector<float> xy = {
    -1, -1,  1, -1,  1,  1,  -1,  1,                        // outer square
    -0.4f, -0.4f, 0.4f, -0.4f, 0.4f, 0.4f, -0.4f, 0.4f};    // inner square
for (int i = 0; i < 500; ++i) {                             // loose points
  xy.push_back(tf::random(-1.f, 1.f));
  xy.push_back(tf::random(-1.f, 1.f));
}
std::vector<int> ring = {0, 1, 1, 2, 2, 3, 3, 0,   // outer square
                         4, 5, 5, 6, 6, 7, 7, 4};  // inner square: the hole

auto points = tf::make_points<2>(xy);
auto edges = tf::make_edges(tf::make_range(ring));
auto polys = tf::make_cdt(points, edges);

import numpy as np
import trueform as tf

outer = np.array([[-1, -1], [1, -1], [1, 1], [-1, 1]], dtype=np.float32)
inner = 0.4 * outer
loose = np.random.default_rng(0).uniform(-1, 1, (500, 2)).astype(np.float32)
points = np.vstack([outer, inner, loose])
edges = np.array([[0, 1], [1, 2], [2, 3], [3, 0],     # outer square
                  [4, 5], [5, 6], [6, 7], [7, 4]],    # inner square: the hole
                 dtype=np.int32)

faces, out_points = tf.cdt(points, edges)
print(tf.area((faces, out_points)))   # 3.36 = 4 - 0.64

faces, _ = tf.cdt(points, edges, split_constraints=False)
print(len(faces) > 0)                 # True: no two constraints cross

import * as tf from "@polydera/trueform";

const n = 500;
const xy = new Float32Array(2 * (8 + n));
xy.set([-1, -1, 1, -1, 1, 1, -1, 1,                      // outer square
        -0.4, -0.4, 0.4, -0.4, 0.4, 0.4, -0.4, 0.4]);    // inner square
for (let i = 16; i < xy.length; i++) xy[i] = Math.random() * 2 - 1;
const points = tf.ndarray(xy, [8 + n, 2]);

const edges = tf.ndarray(new Int32Array([
  0, 1, 1, 2, 2, 3, 3, 0,   // outer square
  4, 5, 5, 6, 6, 7, 7, 4,   // inner square: the hole
]), [8, 2]);

const { faces, points: outPoints } = tf.cdt(points, { edges });

The faces index the output points, not your input array, and the output order can differ. Float64 input gives float64 output on the 64-bit grid. In JavaScript, tf.async.cdt runs the same call off the main thread.

Further

Tutorials. C++, Python and JavaScript: plain Delaunay, outlines with holes, crossings and feature lines, region labels and index maps.

Documentation. Constrained Delaunay triangulation in the topology module, in C++, Python and TypeScript. The C++ page also documents tf::cdt_refiner.

Related. Outlines from 3D meshes: mesh slicing and cross-sections. The same CDT inside other operations: mesh booleans, CSG for engineering and isocontours.

Triangulation is one part of trueform. The same library does booleans, spatial queries, curvature and remeshing, in C++, Python and the browser.

@article{polydera:constrained-delaunay-triangulation,
  title={Constrained Delaunay triangulation},
  author={Sajovic, {\v{Z}}iga, Polydera},
  year={2026},
  url={https://polydera.com/algorithms/constrained-delaunay-triangulation},
  organization={Polydera}
}