Manifold Library

This page describes the theory and interesting implementation details of these algorithms. As I haven't had time to write any of this up for proper peer review yet, this is as close to a paper as you'll get for now. Before the technical write-up, here is a road map of the features implemented so far and what could come next:

Road Map

Manifoldness

What does it mean to be manifold and why should you care? Here we are dealing with triangle meshes, which are a simple form of B-rep, or boundary representation, which is to say a solid that is represented implicitly by defining only its boundary, or surface. The trouble is, this implicit definition of solidity is only valid if the boundary is manifold, meaning has no gaps, holes or edges. Think of a Mobius strip: it is a surface, but it does not represent the boundary of any solid.

This property is not important in 3D rendering, and is thus largely ignored in the animation world, and often even in the CAD world. It is important in engineering and manufacturing, because a non-manifold mesh does not define a solid and so any processing related to analyzing or building that object will either fail or become unreliable. Thus a large amount of middle-ware, both automated and manual, exists for fixing CAD models so they can be used downstream. Closing the design loop with powerful analysis tools like FEA and CFD in an automated optimization could vastly improve the efficiency of design, especially in conjunction with additive manufacturing's ability to easily build complex geometry. However, the lack of reliability of the computational geometry foundations, combined with the manual nature of CAD, means this level of round-trip automation is not currently feasible.

The goal of this library is to create a robust computational geometry foundation from which advanced use cases like the above can be realized.

Manifoldness definition

There are many definitions of manifoldness, most of them from the realm of mathematics. Common definitions will tell you things like the set of manifold objects is not closed under Boolean operations (joining solids) because if you join two cubes along an edge, that edge will join to four surfaces instead of two, so it is no longer manifold. This is not a terribly useful definition in the realm of computational geometry.

Instead, we choose here to separate the concepts of geometry and topology. We will define geometry as anything involving position in space. These will generally be represented in floating-point and mostly refers to vertex properties like positions, normals, etc (see Overlaps). Topology on the other hand, will be anything related to connection. These will generally be represented as integers and applies to faces and edges, which refer to each other and vertices by index.

The key thing to remember is that topology is exact, while geometry is inexact due to floating-point rounding. For this reason, we will choose a definition of manifoldness which relies solely on topology, so as to get consistent results.

As such we will use the same definition of manifoldness as is used in the 3D Manufacturing Format spec: 3MF Core Spec, section 4.1 (disclosure: I wrote most of the 3MF core spec while I was a member of their working group). Paraphrased:

Every edge of every triangle must contain the same two vertices (by index) as exactly one other triangle edge, and the start and end vertices must switch places between these two edges. The triangle vertices must appear in counter-clockwise order when viewed from the outside of the manifold.

This definition has the convenient property that the set of manifold meshes is closed under Boolean operations, since duplicate vertices are allowed as that is a geometric property.

Mesh Boolean

A mesh Boolean operation is the 3D equivalent of the logical Boolean operation, where every point in space is either true (inside the object) or false (outside). The result of a Boolean operation between two solids is another solid, where a Boolean OR is the union of the solids, AND NOT is the difference, while AND is the intersection.

The method used here to create a numerically stable, guaranteed manifold result comes from Julian Smith, whose dissertation includes not only several very powerful computational geometry algorithms, but also a very clear description of the numerical difficulties associated with computational geometry generally. I highly encourage reading it in its entirety if you want a solid grounding in this field. A very brief synopsis is that Euclidean geometry is based on real numbers, and does not hold in a finite-precision computer. Most existing methods either assume "general position" (none of the inputs are close to coincident or coplanar, which is rarely true in practice) or they use exact or arbitrary precision arithmetic, which is slow, complex, and often still suffers failures in edge cases.

Smith's approach was to instead accept the error inherent both in floating-point arithmetic and in float-point representations (after all, even the inputs are rounded so treating them as exact is somewhat meaningless) and focus instead on creating a result that is manifold by construction. One of the best things about his solution is that it does not rely on any kind of exact arithmetic kernel, but instead ensures that each geometric computation is based on those before it, in such a way that the same question is never asked in two different ways. This is key because while Euclidean geometry might say that one result will imply another, with floating-point error that is often no longer true, resulting in conflicting answers that derail an algorithm.

Novel contributions

My contributions, besides building this as an open source library, include parallelizing his algorithm and adding a parallel broad phase collider to get it down from O(n²) to O(nlogn) - more detail below. I also changed his symbolic perturbation to give more useful and intuitive results. The Boolean does not account for any tolerance, so symbolic perturbation is used to break ties for geometry that is exactly floating-point equal in any of the three axes. For axis-aligned, coincident faces (which are very common), it is convenient to avoid creating infinitely-thin shards and holes. This version perturbs the first mesh in the surface normal direction for a Union operation and the opposite for a Difference or Intersection operation, thus ensuring that e.g. touching cubes are merged, equal-height differences produce through-holes, and a mesh minus itself produces an empty mesh.

A couple of key points glossed over in Smith's dissertation were the difficulties of triangulation and degenerate removal, which I have now overcome after much effort. My other contributions are in developing what I hope is a clear and pleasant API to use, including the ability to keep track of arbitrary properties on these manifolds as they undergo operations so that things like textures can be naturally supported without putting any restrictions on which properties or how they should be handled. I also created a way to smoothly refine a manifold, allowing simple triangle meshes to approach the arbitrary precision of NURBS, which are the basis of most CAD software.

Performance

It's difficult to apply the big-O notion of performance to the mesh Boolean operation because of its very nature. Even if we assume that both input meshes have O(n) triangles, the resulting number of intersections (new vertices/triangles) could be anywhere from zero to O(n²) for sufficiently contrived geometry. However, to make the analysis somewhat meaningful for realistic cases, it is convenient to consider spheres, where the input's O(n) refers to partitioning the surface into n triangles of roughly equal edge lengths. If the two input spheres have a non-trivial overlap, then it is easy to see that there will be O(n^1/2) intersections, since n refers to a two-dimensional surface, while the intersections refer to a one-dimensional curve.

Smith's algorithm would still be O(n²), since each triangle of the first mesh would need to be checked for intersection with each triangle from the second mesh. By using a Bounding Volume Hierarchy (BVH) broad phase, this is reduced to O(nlogn) work, and better yet it is parallelized across a GPU. The rest of the computation ideally involves just the intersections, and so is only O(n^1/2), though in practice this often takes longer than the broad phase. Largely this is because the triangulation and degenerate removal steps have not been parallelized.

This code base makes extensive use of NVidia's Thrust library, which is convenient for the vector-based parallelization used here and also allows the code to be compiled for either CPU or GPU. This is especially important since many computers don't have an NVidia GPU, for instance the free CI built into GitHub Actions. However, optimizing code for GPU parallelization also helps performance on CPU-based machines, as these practices aid in efficient cache usage and in pipelining operations. Now that C++20 parallel algorithms are available, it's on the road map to switch, since that's what the Thrust library led to anyway and it will simplify and probably speed up the code.

Overlaps

While manifoldness is the exact, topological part of our system, this is a geometry library and the geometric aspects are the important and tricky parts to solve. The whole point of building a mesh Boolean algorithm is to take care of geometric overlaps in a manifold way. If that is not necessary for your application, you can simply use something like CSG (Constructive Solid Geometry), which tends to just take care of rendering the result correctly, rather than actually calculating all of the intersections.

Given that the whole point of the mesh Boolean is to remove overlaps, why not make that part of the guarantee? The reason is better described by Julian Smith, but to paraphrase: it's hard. Limited machine precision means it is nearly impossible to get consistent geometric results. This is a problem not just for a Boolean algorithm, but even for simple linear operations like rotation. If we loosely define a marginal mesh as one which is close to self-overlapping, but not quite (by some arbitrary metric), then the rounding involved with rotating it by some angle can easily make it cross over to actually self-overlapping.

Smith's algorithm is actually capable of handling self-overlapping input, but this creates self-overlapping polygons in the output, and triangulating these is so arbitrary that it often turns small overlaps into huge overlaps. Therefore, this 3D Boolean instead only guarantees a solution in the case of non-overlapping input. However, given the output of one operation should be valid input to the next, and given the above difficulties with marginal input, we need a different definition of valid input.

Definition of ε-valid

We define non-overlapping geometry in the usual fashion: taking the inputs as exact, every point in space that is not on the boundary has a winding number of either 0 or 1 (outside or inside, respectively). An ε-valid input is one for which there exists a perturbation of the vertices (each by a magnitude < ε) for which this perturbation is non-overlapping.

For convenience, when describing a manifold as self-overlapping, what we mean precisely is: not ε-valid. ε is known as the manifold's precision, which is generally a bound on its accrued numerical error rather than the precision of its floating-point representation. This definition applies equally to 3D manifolds and 2D polygons.

By this definition, the inputs and results are considered approximate, bounded by their respective precisions which are tracked internally. Therefore, any geometric differences within this precision are considered "in the noise" and are not considered to affect the geometric validity of the result. This is used extensively by the triangulation and degenerate removal steps to simplify the manifold and increase robustness of the algorithm.

Removing self-overlaps

Of course wouldn't it be great to remove the overlaps from a single mesh rather than just between two independent meshes? Yes it would, but Julian Smith describes why his Boolean algorithm cannot be applied to the single mesh case. A sweep plane algorithm should be capable of this, but it won't be as parallelizable or simple, but it's still on my road map.

Polygon triangulation

This mesh Boolean algorithm takes as input triangle meshes, but in its natural form, returns a polygonal mesh. Since the purpose here is to cycle back into further Boolean calculations, this necessitates the triangulation of these polygonal faces. This has turned out to be perhaps the trickiest part of this whole Boolean algorithm, in part because it was largely glossed over in the original paper.

Prior art

Why so tricky? A number of polygon triangulators exist, in both open and closed source. I had hoped to use one, but determined that in fact none of them met my needs. Probably the most promising was FIST, but it does not maintain manifoldness because instead of respecting the input edge orientation, it determines the orientation itself based on how it finds the contours to encapsulate each other.

Also promising was Earcut, based on Dave Eberly's code and writeup, but the problem here is in the handling of degenerate cases and self-intersections. The Boolean is robust to, and tends to generate, degenerate shapes (e.g. zero-length edges, coincident edges/faces, etc). Therefore it is critical that the triangulator take degenerate (ε-valid) polygons (and these are not necessarily simple polygons either; any depth of nested contours is possible) and output an ε-valid triangulation. This is a problem I have not seen any polygon triangulator attempt to solve. Usually they throw up their hands when self-intersection is detected, which is usually a toss-up in the case of degeneracy.

My solution

The triangulators linked above work by ear-clipping, which is a method of triangulating a simple polygon that is fairly straight-forward. However, in order to handle non-simple polygons, they must first keyhole the contours together to form a simple polygon. Unfortunately this key-holing step is where most of the code complexity comes in, and tends to be the source of most of the bugs. I opted to solve polygon triangulation using monotone subdivision instead, which is a sweep-line algorithm. It has the advantage of being O(nlogn) instead of O(n²) like ear-clipping. Some say it is more complicated than ear-clipping, but I would say perhaps not more complicated than key-holing, at least if one is going for a solution robust to degeneracies.

The trouble with ε-valid input is that the triangulator effectively has to find a valid perturbation, and this perturbation affects things like sweep-line ordering. Therefore my version of monotone subdivision actually runs the sweep-line forwards and backwards, where the first pass reorders the degenerate points based on non-degenerate geometry that may be distantly-connected, while the second pass does the subdivision based on the now valid perturbation. Unlike Smith's work on polygonal manifolds, I don't think I can prove this triangulator never fails, but it is exceptionally robust from having a lot of truly hideous degenerate polygons with lots of holes thrown at it by the Boolean op. Even so, I eventually realized that the degenerates would just keep piling up as more Boolean operations were strung together and I'd never get robustness or clean output unless I removed them.

Degenerate removal

Mesh simplification by edge collapse & edge swap is pretty well-trod territory, and considering I'm not even trying to operate on non-degenerate situations, it seemed that it should be fairly straightforward. Still, it was educational for me since there are a lot of tricky edge collapse situations that break manifoldness, generally by making four or more triangles meet at one edge. I handle these by duplicating the vertices and splitting the edge. An interesting effect of this strategy is that it only reduces the genus (number of handles) of the mesh, which means it helps clean up the topology of the mesh as well as the degenerate triangles. Note that it can also split the mesh into separate manifolds, as increasing the number of closed meshes also decreases the genus. The only way the genus increases is when a mesh collapses so far that it is completely removed, which I also consider topological cleanup.

Topological cleanup is an key step in this algorithm because the mesh Boolean operates exactly using symbolic perturbation, but on rounded floating-point input, which means it tends to create not just a lot of degenerate triangles, but also non-Euclidean topology in geometric features thinner than ε. To avoid confusing other computational geometry algorithms down the line, it's important to remove these invisible bits of cruft.

I stress-test degenerate removal (and triangulation and the Boolean itself) by creating a Menger Sponge through three difference operations (cutting the holes through a cube in the X, Y, and Z directions). At the fourth order, this object can be shown to have a genus of 26,433. And since nearly every intersection is degenerate, even though the degenerate removal step in this case consumes a significant fraction of the run time, it still runs much faster than without it because it decreases the number of triangles being processed by the other steps by nearly an order of magnitude.

It turns out I cannot guarantee that all degenerate triangles (height < ε) are removed, since there are cases where the only way to do so would involve perturbing the shape of the object by more than ε. Therefore you cannot expect the output to be completely free of degenerate triangles, but rather to have the vast majority removed so as to increase the reliability and speed of many chained operations.

Vertex properties

Handling vertex properties is tricky because I did not want to restrict what kinds of properties one could use, but that meant I also would not have enough information to update them internally. I settled instead on returning a MeshRelation object that for each output triangle, gives the meshID and triangle index it was cut from and for each of its vertices, either the original triangle vertex it corresponds to or its barycentric coordinates. This allows any type of property assigned to the original vertices to be interpolated as desired.

It's important to handle the relation at the triangle rather than vertex level, since along a Boolean intersection curve each vertex will have different properties on each side. Graphics libraries are generally not concerned with manifoldness, so they keep the properties 1:1 with the vertices and situations like this are handled by duplicating these vertices and effectively carving up the manifold into disjoint meshes. However, this operation cannot in general be undone because merging the nearby vertices would be a geometric operation and floating-point error means it would not necessarily come back to a manifold topological result.

One of the classes of degeneracy the previous algorithm removes is coplanar triangles, since the triangulation step adds many "flat" edges, which often become extraneous when intersected during a subsequent operation. However, if some of these vertices have different properties on different triangles, this restricts what edges can be collapsed without affecting the interpolated surface properties. To handle this, I allow an arbitrary list of n floating-point properties to be given for the input mesh triangles, along with precision for each property. Internally this is used to determine which triangles are coplanar, not just geometrically, but in property space as well. It's still up to the user to re-apply these properties to the output mesh as desired, in whatever format is appropriate.

Smooth mesh interpolation

Most commonly, the graphics world tends to work in faceted triangle meshes. However, both artists and CAD engineers tend to want to describe smooth shapes with arbitrary precision from relatively few control points. Most of the methods for this are based on quads of some kind: animators often use subdivision surfaces, while CAD tends to use NURBS patches. The trouble is that when intersected they are no longer quads, and while any polygon can be triangulated, not everything can be divided into quads again. Even if it could, it would not necessarily be able to represent the same surface as its parent. For this reason, CAD software tends to avoid calculating explicit geometric intersections and instead retains the original NURBS quads and keeps track of which ones connect to each other. This works fine for graphics purposes and pushes off the numerical difficulties of the Boolean operation to the analysis programs that take CAD data as input.

More to come...

Provide feedback

Saved searches

Use saved searches to filter your results more quickly