Fast and Exact Constrained Delaunay Triangulation in C++

A constrained Delaunay triangulation (CDT) triangulates a set of 2D points and keeps chosen edges in the result. It turns outlines into meshes, cuts holes, and holds feature lines. trueform builds it on exact predicates. Constraints may cross; the crossings are resolved. It is header-only C++17. This page covers plain Delaunay, outlines with holes, crossing constraints and feature lines, region labels, quality refinement, and index maps.

Install

Header-only. It needs C++17 and oneTBB. Add it to your project via CMake FetchContent:

FetchContent_Declare(
  trueform
  GIT_REPOSITORY https://github.com/polydera/trueform
)
FetchContent_MakeAvailable(trueform)
target_link_libraries(my_app PRIVATE tf::trueform)

Also available via conan, nuget, and pip. See the installation guide for all options.

Loading the points

The input is a flat array of 2D coordinates and, optionally, a flat array of constraint edges that index into it. Wrap them as views. No copy is made:

#include <trueform/trueform.hpp>

std::vector<float> xy;  // x0, y0, x1, y1, ...
for (int i = 0; i < 1000; ++i) {
  xy.push_back(tf::random(-1.f, 1.f));
  xy.push_back(tf::random(-1.f, 1.f));
}
auto points = tf::make_points<2>(xy);

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, tf::make_cdt returns the Delaunay triangulation of the convex hull:

auto tri = tf::make_cdt(points);
auto faces = tri.faces();          // triangles, int indices
auto out_points = tri.points();    // 2D points, float

The result is a tf::polygons_buffer. The faces index its own points, not your input array. The order can differ. To map between them, ask for the index map (below).

An outline with a hole

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

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) {                             // interior 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);

float area = 0;
for (auto t : polys.polygons()) area += tf::area(t);
// area == 3.36 = 4 - 0.64

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

Crossing constraints and feature lines

To keep a constraint as an edge without making it a wall, pass a per-edge mask. true marks a wall, false a feature line that splits nothing. Constraints may cross. The crossing point is added to the output:

std::vector<float> square = {-1, -1, 1, -1, 1, 1, -1, 1};
std::vector<int> lines = {0, 1, 1, 2, 2, 3, 3, 0,  // outline
                          0, 2, 1, 3};             // diagonals
std::vector<bool> is_boundary = {true, true, true, true, false, false};

auto sq = tf::make_cdt(tf::make_points<2>(square),
                       tf::make_edges(tf::make_range(lines)),
                       tf::make_range(is_boundary));
// 4 triangles, 5 points: the corners and the crossing at (0, 0)

If crossings are not allowed in your data, pass false as the last argument. The call then refuses on a crossing and returns an empty buffer. The empty result is the answer to "do these constraints cross", not an error:

auto strict = tf::make_cdt(tf::make_points<2>(square),
                           tf::make_edges(tf::make_range(lines)),
                           tf::make_range(is_boundary), false);
bool crossed = strict.polygons().size() == 0;  // true: the diagonals cross

Region labels

tf::make_cdt keeps only the interior. For every triangle with a label, use the tf::constrained_delaunay_triangulator class:

tf::constrained_delaunay_triangulator<int, float> cdt;
cdt.build(points, edges);

auto all_faces = cdt.make_faces();     // every triangle of the hull
auto labels = cdt.region_labels();     // per triangle, parity of walls crossed
auto all_points = cdt.converted_points();  // float, indexed by make_faces()
// labels[t] % 2 == 1: interior. Otherwise outside or in the hole.

cdt.points() gives the same points on the exact integer lattice the predicates run on. Use converted_points() for coordinates in your own type.

For an id per region cut out by the walls, with 0 for the hull exterior, build in components mode. An island inside a hole then gets its own id:

cdt.build(points, edges, true, tf::cdt_region_mode::components);
// labels: 1 for the ring, 2 for the hole

The true is the split_constraints flag. It must be spelled out before the mode.

Quality refinement

A CDT keeps the triangles the points give it, slivers included. tf::cdt_refiner inserts points until every triangle meets a quality floor (Ruppert refinement):

tf::cdt_refine_config quality;
quality.min_quality = 0.3f;        // 1 is equilateral; values above 0.45 are clamped
quality.split_encroached = true;   // allow midpoint splits of constraints

tf::cdt_refiner<int, float> refiner;
if (refiner.build(points, edges, quality)) {
  auto rfaces = refiner.make_faces();
  auto rlabels = refiner.region_labels();
  auto rpoints = refiner.converted_points();
}

Every split of a constraint is recorded as a parameter along the original input edge, so several triangulations that share a polyline can agree on the same split points.

Working with results

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

auto [mapped, im] = tf::make_cdt(points, edges, tf::return_index_map);
auto sentinel = im.f().size();
int dropped = 0;
for (auto out_id : im.f()) dropped += out_id == sentinel;

im.f()[i] is the output index of input point i, or the sentinel if the point fell outside the kept region. im.kept_ids()[k] is the input index of output point k, or the sentinel for a point made at a constraint crossing. Coincident input points are merged.

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

auto p2 = polys.polygons();
tf::points_buffer<float, 3> lifted;
for (auto p : p2.points()) lifted.emplace_back(p[0], p[1], 0.f);
auto mesh3 = tf::make_polygons(p2.faces(), lifted.points());
tf::write_obj(mesh3, "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 C++ — plane sections, layers and slabs of a mesh.
• Fast Mesh Booleans in C++ — exact union, intersection and difference of two meshes.
• Fast and Robust Isocontours on a Triangle Mesh in C++ — curves and bands of equal value on a scalar field.
• Fast STL and OBJ Reading in C++ — load and write meshes, the first step of every pipeline.

Performance

Arrays in, triangle arrays out, at 1M random points, in C++. Region labels are excluded for every library. Every row is validated.

Apple M4 Max 16 threads Clang -O3

Of the seven libraries measured, trueform is the fastest on both tasks at 1M points. Delaunay32 is 1.4× slower on plain Delaunay and has no constrained mode. The constrained case uses non-crossing constraints, one per about ten points. trueform runs on float32 coordinates with exact int32 predicates. The other libraries use their own native kernels.

Full benchmarks · Documentation

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