pure-cdt

Constrained Delaunay triangulation and mesh refinement in 2D, written in Haskell with no bindings to a C or C++ library.

A gear with a keyed bore: constrainedTriangulate on the two outlines alone, then refine to a 28° minimum angle and a maximum triangle area. The only input points are the polygon vertices.

Usage

import Data.List.NonEmpty (NonEmpty (..))
import Triangulation

-- Delaunay triangulation of a point set (Nothing if it spans no triangle).
Just t = triangulate [Point 0 0, Point 4 0, Point 2 3, Point 2 1]
triangles t  -- :: [Triangle]
hull t       -- :: Polygon, the convex hull

-- A polygon with holes: the boundary and hole edges are forced into the
-- triangulation, triangles outside the boundary or inside a hole are dropped.
-- The third argument holds interior points, if there are any.
boundary = Polygon (Point 0 0 :| [Point 0 10, Point 10 10, Point 10 0])
hole     = Polygon (Point 4 4 :| [Point 4 6, Point 6 6, Point 6 4])
Just ts = constrainedTriangulate boundary [hole] []

-- Mesh refinement (Ruppert's algorithm): insert vertices until no angle is
-- below 25° and no triangle is larger than 0.5, keeping the polygon edges.
Just fine = refine defaultQuality {minAngle = 25, maxArea = Just 0.5} ts

-- The indexed form finite element code wants: nodes, and triangles as
-- counter-clockwise triples of node indices; plus the boundary as index pairs.
mesh = toMesh fine
nodes mesh         -- :: Vector Point
elements mesh      -- :: Vector (Int, Int, Int)
boundaryEdges mesh -- :: [(Int, Int)], the edges that belong to a single triangle

Two conventions to know before the first call: the y axis points up, and polygons are wound clockwise, so that walking along a boundary the interior is on the right. Edge and Triangle are abstract; build them with mkEdge/mkTriangle and take them apart with the read-only patterns Edge a b/Triangle a b c.

What this package is for

It does not compete on speed. A mature Delaunay implementation in C or C++ is an order of magnitude or two faster, and nothing here will close that gap. What it offers instead:

• A pure API. A handful of ordinary functions that can be mapped over a list or used inside a QuickCheck property. A point set that spans no triangle comes back as Nothing: there are no partial functions and no exceptions in the library.
• No toolchain beyond GHC. The dependencies are base, deepseq, hashable, parallel, random, unordered-containers and vector, so the package builds wherever GHC does, including under Nix, when cross-compiling, and on the JavaScript and WebAssembly backends.
• Degenerate input handled, not merely tolerated. The predicates are exact and the algorithm runs on symbolically perturbed points (Edelsbrunner and Mücke's Simulation of Simplicity), so lattices, regular polygons and concentric shapes are triangulated correctly instead of coming out with holes or overlaps. They cost about three times as long as points in general position, and that is a constant factor, not a worse complexity.
• A small readable core. The code reads as geometry rather than as bookkeeping of array indices.

Making it faster would mean giving up the last of those. Passing indices into an array instead of Point values through the algorithm is likely worth a factor of two or three, and that trade has been declined deliberately.

Performance

Wall-clock time, median of three runs on an Intel Core i7-14700KF (eight performance cores plus twelve efficiency cores), GHC 9.10.3:

Halves of a large point set are triangulated in parallel and then merged. The speedup plateaus between eight and twelve threads and falls away after that, because the merges along the spine of the recursion are sequential and the parallel collector has more capabilities to synchronise.

Building

stack build          # library and the example executable
stack test           # tasty: unit tests and QuickCheck properties
stack bench          # tasty-bench: triangulate, refine, toMesh
stack exec pure-cdt-examples   # renders the examples to assets/*.png

# wall-clock scaling, and the same on a regular lattice
stack bench pure-cdt:pure-cdt-scaling --ba '50000 500000 +RTS -N8'
stack bench pure-cdt:pure-cdt-scaling --ba 'lattice 100000 +RTS -N8'

tasty-bench reports CPU time, which grows with the number of cores and so says nothing about a parallel speedup; the table above comes from pure-cdt-scaling, which measures the clock. The example renderer sits behind the cabal flag examples, off by default, so that the library carries no image dependencies; stack.yaml turns it on for development, and plain cabal takes -f examples.

Formatting is fourmolu, linting is hlint; both run in CI.

Limitations

• Polygons must be simple and must not cross each other. A hole may touch the boundary at a vertex, in which case the corner is cut off.
• Input whose points all lie on one line has no triangulation, and triangulate returns Nothing for it.
• Ruppert's refinement is guaranteed to terminate for angle bounds up to about 20.7° when no two polygon edges meet at less than 60°, and in practice works up to about 30°. A demanding bound next to a small input angle can exhaust the insertion budget; refine then returns Nothing, and refineWithBudget returns the mesh it reached.
• refine works on the triangles of a region, so an interior constraint edge with triangles on both sides is not preserved. It splits segments at points computed in floating point, so a vertex it adds to the boundary lies on the original polygon edge only up to rounding.

License

BSD-3-Clause.