Fast and Exact Constrained Delaunay Triangulation in Python

A constrained Delaunay triangulation (CDT) triangulates a set of 2D points and keeps chosen edges in the output. Those edges can be an outline, the rim of a hole, or a feature line. trueform does it with tf.cdt. Constraint edges may cross each other; the crossings are split for you. One pip install. NumPy arrays in, NumPy arrays out. This tutorial covers plain Delaunay, outlines with holes, crossing constraints, region labels and index maps, then shows measured times against Triangle and SciPy.

Install

pip install trueform

See the installation guide for all options.

Loading the points

tf.cdt takes an (N, 2) point array and, optionally, an (E, 2) int32 array of constraint edges that index into it. No mesh object is needed.

import numpy as np
import trueform as tf

rng = np.random.default_rng(0)
points = rng.uniform(-1, 1, (1000, 2)).astype(np.float32)

If your points come from a 3D file, read it with tf.read_stl or tf.read_obj and project to two coordinates first.

Delaunay triangulation

Without edges you get the Delaunay triangulation of the convex hull:

faces, out_points = tf.cdt(points)
print(faces.shape, faces.dtype)     # (1980, 3) int32
print(out_points.shape)             # (1000, 2)

The faces index out_points, not your input array. The output order can differ from the input order. To map between them, use return_index_map=True (below).

An outline with a hole

Pass the outline and the hole as constraint edges. Every constraint edge is a wall by default. Triangles inside an odd number of walls are kept, so a ring of walls inside another ring is a hole.

outer = np.array([[-1, -1], [1, -1], [1, 1], [-1, 1]], dtype=np.float32)
inner = 0.4 * outer
steiner = rng.uniform(-1, 1, (500, 2)).astype(np.float32)
ring_points = np.vstack([outer, inner, steiner])

ring_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(ring_points, ring_edges)
print(faces.shape)                  # (858, 3)
print(tf.area((faces, out_points))) # 3.36... = 4 - 0.64

The interior points inside the hole are dropped, and the area is the square minus the hole.

Crossing constraints and feature lines

Constraints may cross. The crossing point is added to the output:

diagonals = np.array([[0, 2], [1, 3]], dtype=np.int32)
faces, out_points = tf.cdt(outer, diagonals)
print(out_points.shape)             # (5, 2): the four corners and the crossing at (0, 0)

To keep a constraint as an edge without making it a wall, pass an edge_mask. True marks a wall, False a feature line that splits nothing:

edges = np.array([
    [0, 1], [1, 2], [2, 3], [3, 0],     # outline: walls
    [0, 2], [1, 3],                     # diagonals: feature lines
], dtype=np.int32)
edge_mask = np.array([True, True, True, True, False, False])

faces, out_points = tf.cdt(outer, edges, edge_mask=edge_mask)
print(faces.shape)                  # (4, 3): four quarter triangles, square stays solid

If crossings are not allowed in your data, pass split_constraints=False. The call then returns empty arrays on the first crossing. The empty result is the answer to "do these constraints cross", not an error.

faces, out_points = tf.cdt(outer, diagonals, split_constraints=False)
print(faces.shape[0] > 0)           # False: the diagonals cross

Region labels

By default you get only the interior. Pass region_labels to get every triangle with a label per triangle:

faces, out_points, labels = tf.cdt(ring_points, ring_edges, region_labels="nesting")
print(np.unique(labels))            # [0 1]: 0 outside or in the hole, 1 interior
interior = faces[labels == 1]       # the 858 triangles of the default call

"nesting" gives the parity of walls crossed from outside. "components" gives an id per region cut out by the walls, with 0 for the hull exterior, so an island inside a hole gets its own id:

faces, out_points, labels = tf.cdt(ring_points, ring_edges, region_labels="components")
print(np.unique(labels))            # [1 2]: the ring and the hole

Labels are int32, one per output triangle. region_labels needs edges.

Working with results

faces is (K, 3) int32 and out_points is (M, 2) in the input dtype. Pass float64 points for float64 output. The pair is a 2D mesh; tf.Mesh(faces, out_points) wraps it for further queries.

To trace input points through the triangulation, ask for the index map:

(faces, out_points), (f, kept_ids) = tf.cdt(ring_points, ring_edges, return_index_map=True)

sentinel = f.size
dropped = f == sentinel             # input points that fell outside the kept region
created = kept_ids == sentinel      # output points made at constraint crossings
print(dropped.sum(), created.sum()) # 75 0

f[i] is the output index of input point i. kept_ids[k] is the input index of output point k. With region_labels the return shape is ((faces, out_points), labels, (f, kept_ids)).

To lift a 2D triangulation into 3D, add a third coordinate and write it out:

points3 = np.column_stack([out_points, np.zeros(len(out_points), dtype=np.float32)])
tf.write_obj((faces, points3), "cdt.obj")

The same CDT triangulates the cut surfaces in mesh booleans. For slicing a 3D mesh into 2D outlines, see mesh slicing and cross-sections.

Related tutorials

• Fast and Robust Mesh Slicing and Cross-Sections in Python — plane sections, layers and slabs of a mesh.
• Fast Mesh Booleans in Python — exact union, intersection and difference of two meshes.
• Fast and Robust Isocontours on a Triangle Mesh in Python — curves and bands of equal value on a scalar field.
• Fast STL and OBJ Reading in Python — load and write meshes, the first step of every pipeline.

Performance

Delaunay triangulation of 1M random points, no constraints, for all three libraries. SciPy wraps Qhull, which computes an unconstrained Delaunay triangulation only.

Apple M4 Max 16 threads

Median of three fresh-process medians. Benchmark run with trueform 0.10.3, Triangle 20250106 and SciPy 1.17.1.

With constraints, on a wiggly closed contour of 100,000 vertices, the constrained triangulation takes 23.7 ms against Triangle's 39.7 ms (1.7×). SciPy has no constrained mode.

Full benchmarks · Documentation

@article{polydera:constrained-delaunay-triangulation-in-python,
  title={Fast and Exact Constrained Delaunay Triangulation in Python},
  author={Sajovic, {\v{Z}}iga, Polydera},
  year={2026},
  url={https://polydera.com/tutorials/constrained-delaunay-triangulation-in-python},
  organization={Polydera}
}