Robust Mesh Booleans and the Exact-Arithmetic Escape

Boolean union, intersection, and difference on triangle meshes are specified by a set-theoretic statement: the result is the boundary of the region that is inside A and inside B, or inside either. Every real implementation of that specification is a sequence of geometric predicates, and every one of those predicates becomes ambiguous when its inputs are nearly or exactly degenerate. The difficulty is not the algorithm. It is that the algorithm’s branch structure depends on predicates whose answers are not representable.

Why the naive pipeline fails

The textbook procedure has three phases. Intersect every triangle of A with every triangle of B, producing a set of segments. Retriangulate each face along the segments lying in its plane. Classify each resulting face as inside, outside, or on the boundary of the other solid, and keep the ones that belong to the result.

Each phase has a failure mode. Pairwise triangle intersection needs a decision for two triangles that share a plane and overlap, and for two that touch along a single point. The segment produced by a transversal intersection has endpoints computed by solving a linear system, and those endpoints are irrational in general, so they cannot be stored exactly as floating-point coordinates. Retriangulation must handle segments that cross each other and segments that end in the interior of another segment. Classification casts a ray from a face and counts crossings, which is ambiguous whenever the ray hits an edge exactly.

The practical consequence is a class of bug that appears only on certain models: a coplanar top face touching another solid’s top face, two components meeting at an edge exactly, a rotationally symmetric part whose vertices land on lattice points. These are common in real datasets, not rare, because parts are frequently modelled on axis-aligned grids and assembled by snapping.

Exact predicates, approximate constructions

The split that makes robust implementations possible is between predicates and constructions. A predicate is a sign decision: is point $p$ above, below, or on the plane of $a, b, c$. It has a finite answer, and it can be computed exactly with enough precision. A construction is a coordinate: the point where a segment meets a plane. It generally has no exact representation, and its value is rounded to whatever the output format allows.

The classic example is the orientation predicate, the signed volume of the tetrahedron formed by four points:

$$ \text{orient}(a,b,c,d) = \det \begin{pmatrix} a_x - d_x & a_y - d_y & a_z - d_z \ b_x - d_x & b_y - d_y & b_z - d_z \ c_x - d_x & c_y - d_y & c_z - d_z \end{pmatrix}. $$

Computed in double precision, the sign of this determinant is reliable when the magnitude is far from zero and unreliable when it is near zero, because the subtraction of nearly equal coordinates loses bits and the products amplify the loss. The determinant of a $3\times3$ matrix of doubles is computed with a few hundred floating-point operations, and the error is bounded by a small multiple of the operation count times the machine epsilon times a bound on the intermediate magnitudes. That bound is computable, which makes the sign computable: compute the determinant in double precision along with its error bound; if the magnitude exceeds the bound, the sign is certain; otherwise recompute in higher precision or exactly.

Shewchuk’s adaptive predicates do exactly this, escalating from double precision to a sequence of increasingly exact expansions, and the escalation happens only for near-degenerate inputs, so the average cost stays near the fast path. This is the reason exactness is affordable at all: the expensive path is taken by a small fraction of predicates.

Exact constructions with rational coordinates

To place a new vertex where a segment crosses a plane exactly, the coordinates must be extended from floating-point numbers to a field closed under the operations used. If input coordinates are integers, all constructions are rational, and the rationals are closed under addition, subtraction, multiplication, and division. An intersection point computed as a ratio of integer determinants is therefore exactly representable as a ratio of two arbitrary-precision integers.

Working entirely in rationals is correct and slow: the bit lengths of the numerators and denominators grow with every chained construction, and a chain of dependent intersections produces numbers with hundreds of bits. The Cherchi–Pellacini–Attene approach avoids this by keeping the arrangement combinatorial. Predicates are evaluated exactly, but the new vertices are stored as the geometric primitives that produced them — a plane and a line, say — and coordinates are computed only when needed, in floating point, with the exact rational available as a tiebreaker.

The output therefore has rounding error, but no topological error. That is the crucial distinction: the combinatorial structure of the result — which faces exist, which edges bound which faces — is decided by exact predicates and is therefore a genuine arrangement of the input. The coordinates are approximations of the exact points, but the connectivity describes the correct surface. Without this separation, a rounded coordinate can land on the wrong side of a plane and flip a predicate, and the result’s connectivity contradicts its own geometry.

Coplanar faces and the arrangement view

The hardest common case is the coplanar one. Two solids whose top faces lie in the same plane have a two-dimensional region of contact, and the intersection is not a set of segments but a set of polygonal areas. Deciding which parts of A’s top face are inside B requires a two-dimensional arrangement of the overlap polygon, which is a different algorithm from the three-dimensional triangle intersection.

The robust formulation treats the whole problem as a single arrangement: overlay the two surfaces, compute all intersections in three dimensions including coplanar ones, split faces along every intersection curve, and then classify each cell by a winding number or a ray cast from a point known to be strictly interior to the cell. Coplanar contacts become zero-thickness cells, and the classification of those cells is the pivot: a cell on the boundary between inside and outside must be kept and duplicated with both orientations, or discarded, depending on whether the boolean is a union, intersection, or difference.

A winding-number classification has an advantage over ray casting: it is an integral, so it degrades gracefully. Computing the winding number of the other solid around a sample point requires summing signed solid angles of all its triangles, and the result is an integer for non-degenerate points. Near a boundary the sum is close to an integer but not equal; the difference indicates that the sample point was too close to the surface, and a new sample is chosen. A ray cast has no such self-diagnosis: a ray that grazes an edge gives a count that is wrong by one with no indication.

import math

def winding_number(point, tris):
    total = 0.0
    for a, b, c in tris:
        total += solid_angle(point, a, b, c) / (4.0 * math.pi)
    return total          # ~integer; a large fraction flags a bad sample point

The solid angle of a triangle as seen from a point is computed from the vectors to its vertices with the Van Oosterom–Strackee formula, which uses an arctangent of a determinant over a denominator of products of lengths. Its cost is constant per triangle, so a naive classification is linear in the other solid’s face count per query; acceleration structures reduce that to the faces near the query.

Where the cost lands and what the API hides

A robust boolean’s running time is dominated by the arrangement, which in the worst case has complexity proportional to the product of the two surface complexities: $n$ triangles of A against $m$ of B can produce $O(nm)$ intersections. In practice the surfaces intersect along curves and the actual count is far below the bound, which is why a broad-phase structure that rejects non-overlapping triangle pairs matters more than the predicate cost.

The memory cost is the exact-rational arithmetic. Numbers with a few hundred bits each, held for every constructed vertex, add up quickly; a model with $10^6$ output vertices is not going to be comfortable in a representation that stores a rational per coordinate. This is why production pipelines run the boolean exactly and then immediately round to floating point, keeping the exact data only for the duration of the operation.

The failure that remains is not numerical but semantic. A boolean of two solids that touch only along a shared face has no well-defined answer for a union: the result is the two solids joined, and whether the shared face should exist in the output is a convention. Some kernels remove it as an internal face, some keep it as a zero-area shell. Downstream code that expects one behaviour and receives the other will compute a wrong signed volume or fail a watertightness test, and neither the input nor the output is malformed in any way a checker would notice.

There is also a contrast with the field-based formulation. A signed distance field performs a union with a min and has none of these problems, because it never constructs an arrangement. What it loses is exactness and sharp edges, and what it gains is that the operation cannot fail. Choosing between the two is a choice about which kind of incorrectness is acceptable.