Mesh Arrangements for Solid Geometry
Qingnan Zhou1
1

Eitan Grinspun2

New York University

Denis Zorin1
2

Alec Jacobson2

Columbia University

Figure 1: Our method takes as input any number of meshes (three shown in this 2D illustration). We resolve intersections and assign a
winding number vector to every delineated cell. Different boolean results are extracted according to these winding number vectors.

Abstract
Many high-level geometry processing tasks rely on low-level constructive solid geometry operations. Though trivial for implicit representations, boolean operations are notoriously difficult to executerobustly for explicit boundary representations. Existing methods
for 3D triangle meshes fall short in one way or another. Some methods are fast but fail to produce closed, self-intersection free output.
Other methods are robust but place prohibitively strict assumptions
on the input, e.g., no hollow cavities, non-manifold edges or selfintersections. We propose a systematic recipe for conducting a family of exact constructive solid geometry operations. The two-stage
method makes no general position assumptions and does not resort
to numerical perturbation. The method is variadic, operating on
any number of input meshes. This generalizes unary mesh-repair
operations, classic binary boolean differencing, and n-ary operations such as finding all regions inside at least k out of n inputs.
We demonstrate the superior effectiveness and robustness of our
method on a dataset of 10,000 “real-world” meshes from a popular online repository. To encourage development, validation, andcomparison, we release both our code and dataset to the public.
Keywords: arrangements, constructive solid geometry, booleans
Concepts: •Computing methodologies → Computer graphics;

1

Introduction

Geometric modeling tools are more accessible than ever, scanningtechnologies are available at the commodity level, additive fabrication technologies rapidly grow in popularity, and platforms haveemerged for sharing 3D models online. Unfortunately, the resulting
wealth of 3D models comes with a catch. The diversity of thesemodels coincides with unpredictable mesh quality and structure.
For example, manually sculpted models are created using a broad
range of software by designers with varying skill and intention.

available algorithms impose strict requirements on their inputs. As
the complexity and amount of 3D data grow, manual one-by-one
preprocessing is no longer acceptable. For example, cumbersome
mesh cleanup of solid mesh may take more time than a subsequentphysical simulation of it.
An important class of mesh repair and solid operations aims topreserve input geometry as much as possible when recombining
it in new ways or converting it to suitable form for a downstream
application. These operations are elegantly viewed as operations
on space partitions defined by mesh arrangements. A mesh arrangement, built from a collection of (possibly non-manifold, openboundary, self-intersecting, with degenerate triangles, etc.) meshes,partitions space into a number of cells. This view unifies tasks often
viewed as distinct, such as mesh repair and boolean operations.
We introduce mesh arrangements constructed from a restrictedclass of meshes: those with piecewise-constant winding number, or
PWN. By surveying 10,000 popular meshes, we show that PWN
meshes cover a large fraction of practically relevant situations.
An arrangements of a PWN mesh is a lightweight yet powerfulrepresentation. Casually, a PWN mesh arrangement is composedfrom possibly (self-)overlapping components bounding a numberof solids. These arrangements enable a set of highly robust andconservative algorithms. By robust, we mean that algorithms successfully produce output for all PWN meshes. By conservative, wemean that arrangements preserve original mesh geometry exactly.
Our general approach can be separated into two stages: adding
meshes to an arrangement—independent of the desired operation—and extracting the boundary of the result according to an extraction
function describing the desired boundary in terms of the winding
numbers of the region it bounds (see Figure 1).
The first stage, in turn, consists of resolving intersections, partitioning space into cells, and labeling each cell with winding numberswith respect to each input mesh.

Meanwhile, geometry processing operations increase in sophistication: from remeshing to physical simulation. Yet, many, if not most,In the second stage, all classic boolean operations (e.g., union, intersection, difference) have trivial extraction function definitions.
We also explore other interesting functions, such as extracting thePermission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are notmade or distributed for profit or commercial advantage and that copies bear
this notice and the full citation on the first page. Copyrights for componentsof this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on
servers or to redistribute to lists, requires prior specific permission and/or
a fee. Request permissions from permissions@acm.org. © 2016 Copyrightheld by the owner/author(s). Publication rights licensed to ACM.
SIGGRAPH ’16 Technical Paper, July 24 - 28, 2016, Anaheim, CA
ISBN: 978-1-4503-4279-7/16/07
DOI: http://dx.doi.org/10.1145/2897824.2925901

Figure 2: Self-intersections are not just “artifacts.” They also occur intentionally during modeling. A coffee mug handle is extruded
then curled inward on itself (blue). The self-intersections (orange)
do not prohibit constructing the shape’s underlying torus-topology
self-union with our algorithm (green).

Figure 3: Simple and complex variadic operations cost the sameusing our mesh arrangements. Converting variadic operations to acascade of binary operations is worst-case exponential in time.
“self-union” of a single, overlapping input mesh or extracting all
regions of space inside at least two of many input meshes.
Our method guarantees as output a solid mesh. Informally, a solidmesh is the non-degenerate boundary of a solid subregion of R3 .
Our method also guarantees that the output is exact, i.e., interpreting the input positions as exact rationals, all intersections result
from exact construction. The output coordinates may optionallybe converted to floating point in a post-process.
We validate our method on the 10,000 triangle meshes from theonline 3D printing repository Thingiverse, as well as benchmarksof previous work. We present extensive comparisons with state-ofthe art methods, all of which fail with significant frequency, either
rejecting the input, failing to produce any output, or producing an
output that is not a solid mesh. We also evaluate the conversion of
our exact results to floating-point positions; in this case, we outperform existing floating point methods along the same criteria.

2

Related work

Most previous work separate into two broad groups: boolean operations on different classes of objects and “mesh repair”, in particular,
elimination of self-intersections, computing outer hulls and similar
operations. While both require intersecting meshes or surfaces, inmost cases the problems have been treated disjointly.
Previous boolean methods define a restrictive class of input 3D pointsets that are closed with respect to set
operations. The methods output a restrictive class of boundary representation or spatial partition. In almost all cases, it is assumed
that creating a valid boundary representation from a broader classof inputs is a separate task, delegated to the user or preprocessing.
Inputs not meeting the strict requirements are not handled directly.

Boolean operations

Previous works differ by generality of input/output representations
and tradeoffs between performance and robustness. We compare
directly to the state of the art in Sections 5 and 7. For now, we categorize approaches, highlighting salient similarities and distinctions.
The current standard in robustness is CGAL’s exact-arithmetic implementation [Granados et al. 2003] of Nef polyhedra [Bieri & Nef
1988]. CGAL’s implementation requires a valid Nef polyhedron
data structure (Sphere Maps and Selective Nef Complexes) as input. Current tools in CGAL and O PEN SCAD (a modeling tool
bootstrapping CGAL) only construct this form via embedded polyhedra, excluding inputs with self-intersections, non-manifold features, and inner cavities, although the latter two could be represented by the Nef representation. More fundamentally, the generality of the representation requires a complex heavy-weight data
structure and has a significant impact on performance.
On the other extreme, Douze et al. [2015] restrict inputs to embed-

Figure 4: Combinatorially, the multi-component shapes on the left
and right are the same, though their outer hulls (thick) differ.
ded polyhedra with vertices in general position, but is extremely
efficient and is capable of handling very large meshes. Douze et
al. introduce the concept of variadic boolean operations: immediately evaluating boolean operations involving many inputs rather
than decomposing into a tree of binary operations (see Figure 3).
Recently, Barki et al. [2015] use a more general yet lightweight
representation and exact rational arithmetic to handle a variety ofnear degeneracies. Their efficient algorithm maintains robustness
on a 26-model benchmark.
All methods so far are similar to ours in that they add intersectionsof boundary representations, and proceed to classify elements of theboundary to construct the result. However, these algorithms assume
that input surfaces are free of self-intersections. Self-intersections
are not only a ubiquitous meshing artifact, but also a common way
to model interesting topologies Figure 2.
Many early methods based on intersection and classification suffer from robustness problems, before exact-arithmetic based methods became practical. This leads to development of more robust
approaches, initially based on conversion to volumetric representations, starting from [Museth et al. 2002]. However, the approximation of the original meshes depends on the chosen resolution level
for the volume grid, and high accuracy requires high tessellation.
Different approaches accelerate this approach, reduce complexityof the output (e.g., using adaptivity [Varadhan et al. 2004]), andattempt to preserve the original mesh as much as possible [Pavic
et al. 2010; Wang 2011; Zhao et al. 2011]. The fundamental issue
with such techniques is the approximate and grid-dependent natureof volumetric calculations: while increasing robustness, these maylead to unwieldy topology changes and geometric deviations.
The space-partition view of boolean operations has appeared mostclearly in binary space partitioning (BSP) methods, starting with[Thibault & Naylor 1987; Naylor et al. 1990]. Bernstein & Fussell[2009] combined this with robust predicates to develop an efficientand robust way to compute booleans on surfaces in BSP representation. As in most other works on booleans, conversion to this representation is viewed as a preprocess, with the range of inputs thispreprocess can handle not precisely defined. While compared tovolumetric-grid approaches, BSP methods increase mesh complexity more moderately and input geometry is better preserved, yetsignificant refinement is still needed. Campen & Kobbelt [2010a]
localized their BSP-based method using an octree and perform refinement only locally near intersections. Importantly, this workFigure 5: The outer-hull (e.g., [Attene 2014; Campen & Kobbelt
2010a] is simpler to extract, but not always appropriate. Hollow
cavities should be maintained for 3D printing not only to reducematerial costs, but to maintain functional properties like balancing.

points out that BSP trees provide a general representation of spacepartition that flexibly performs both boolean operations, outer hullcomputations, and other operations. We expand this idea, but use ahigher-level and more compact space partition.

Figure 6: By outputting an exact solid mesh, our method ensures
there are no inversions (illustrated as correct signed distance onleft). Previous methods may create inversions, and even small inversions catastrophically affect signed distance to the output (right).

3
Mesh repair techniques historically deviate fromboolean operations by focusing on converting a maximally broad
range of input meshes to a normalized representation (e.g., closed
manifold meshes without self-intersections). These methods often
rely on volumetric approximation and for certain problems (e.g.,hole-filling) this may be unavoidable. If possible, preserving original geometry is desirable. A common example of mesh repair of
this type is computation of the outer hull, though both state-of-theart methods [Campen & Kobbelt 2010a; Attene 2014] do not appear
to disambiguate nested and non-nested components (see Figure 4).
For example, when preparing models for 3D printing, the outer hull
may be inappropriate as it removes inner cavities (see Figure 5).
Mesh repair

Approaches to robustness, and sources of non-robustness

Many implementations (e.g., [Bernstein 2013; Mei & Tipper 2013;
Douze et al. 2015]) explicitly assume general position of inputs(no four points on a circle, no co-planar intersections, etc.) anddo not attempt to handle numerical non-robustness. The development of new boolean and mesh repair techniques was driven, to alarge extent, by robustness considerations. We consider more explicitly how robustness was addressed in different contexts, and theunresolved problems we are addressing in our approach.
The most comprehensive approach is to represent all objects using exact arithmetic. With advances in filtered predicates for efficiency, this approach is increasingly preferred. We (like others[Granados et al. 2003; Barki et al. 2015]) largely follow it. Earlier BSP-based work used the observation [Sugihara & Iri 1992]
that representing points as the intersection of original planes eliminates the need for exact computations (only exact predicates). This,in principle makes it possible to do most computations robustly infloating point, but some constructions or rounding still inevitably
appear in all methods, and lead to non-robustness, often subtly.
For example, Banerjee & Rossignac [1996] and later Xu & Keyser[2013] build exact topologies but fixed-precision floating point vertex positions, leading to self-intersections, inversions (see Figure 6)
and degeneracies in the output. Campen & Kobbelt [2010a] round
all input vertices aggressively to ensure exact plane intersections
for a BSP-representation. We discuss several other problems in existing techniques in greater detail as we describe our method. For
our approach (and all exact methods), conversion of the output to afloating point (if such a conversion is desired) is a potential sourceof problems, although we have observed it in a very small proportion of cases and attack it with an additional heuristic utilizing ourcore method in Section 6.1.

Concepts

The crux of our method is construction of the mesh arrangement
data structure, consisting of cells annotated with winding numbers,
patches and their adjacency graph, that allows us to extract resultsof a variety of operations from the arrangement.
This structure is a relatively lightweight representation of a spacepartition (cf. BSP trees). Based on compound surface objects(patches), it allows for complex cells (does not require them to beconvex or even topological balls). Yet, it allows us to perform all
operations robustly and efficiently.
The inputs to our arrangement creation algorithm are n piecewiseconstant winding number triangle meshes A1 , . . . , An . For extraction of the results of specific operations from the arrangement, we use a variadic extraction function f to determine thesolid mesh boundary of which region(s) of space carved out by
A1 , . . . , An to output.

3.1

Piecewise-constant winding number meshes

A triangle mesh is a set of 3D vertices (some of which may be geometrically coinciding) and a set of triangles connecting these vertices, each triangle represented by a triplet of vertices, with orientation implied by the vertex order for non-degenerate triangles. We
may view triangles combinatorially as triplets of vertices as well as
geometrically as pointsets in 3D.
Effectively, any valid triangular mesh in Wavefront OBJ-like format, is a valid input, subject to one general condition: we requirethat triangle meshes Ai induce a piecewise-constant integer generalized winding number (PWN) field wi [Jacobson et al. 2013]:
wi (p) ∈ Z∀ p ∈ R \ |Ai |,(1)

where |Ai | denotes the union of all triangles of Ai viewed as pointsets. For a triangle mesh, this is simply the sum of the signed solid
angles Ωt (p) of each oriented triangle t:wi (p) =1 XΩt (p).
4π t∈A(2)i

We call meshes with this property piecewise-constant winding number meshes or PWN meshes.
A PWN mesh Ai can be interpreted as dividing all of R3 into regions that are outside (wi = 0) or inside (wi 6= 0) of the “solid
implied by Ai .” This allows multiplicity (|wi | > 0) for parts ofMethod
CGAL[Campen... 2010a]
Q UICK CSG
C ARVEC ORK[Bernstein... 2009][Barki et al. 2015]

Our method

C AVITIESN ON - MANIFOLDE XACT•

D EGENERACIESS ELF - INTER .

|w| > 1X

M ULTI - COMP.
•S EAMS

X

X

X

X

X

X

X

X

X

X

X

•

X

•
•
•
•
•
•

X

X

X

X

X

X

X

•
•
•
•
•

•
•

X

X

X

X

X

X

X

X

X

X

X

•
•
•

•
•
•

X

X

X

•
•

X

X

•

•

•
•

Table 1: Mesh boolean algorithm input preconditions feature chart: Previous techniques have severe input restrictions.
space considered to be twice, thrice, etc.
inside and also allows for negative insideness (wi < 0) for parts of space inside
an inverted part of Ai (see inset).
PWN meshes may exhibit commonly witnessed “artifacts” making
them unsuitable for previous algorithms (see Table 1 and Figure 7).
N ON - MANIFOLD: PWN mesh connectivity may be non-manifold
at vertices or edges.

Figure 7: Our assumptions allow a wide class of inputs with selfintersections, multiple components and non-manifold connections,but does not include open boundaries or non-manifold flaps.

C OPLANAR, D UPLICATE: Co-planar and duplicate facets (i.e., geometrically identical, but logically distinct: either belonging to different meshes, or using vertices with the same positions) do notnecessarily invalidate a PWN mesh. For example, the conjoining of
two cubes abutting along a triangle is a PWN mesh. The conjoining of two entirely identical cubes is also a PWN mesh. However,a cube with a single duplicated triangle is not a PWN mesh.

3.2

C AVITIES, M ULTI - COMP., |w| > 1: The winding number elegantly handles correctly oriented boundaries of multiple connectedcomponents or nested shells, such as a hollowed-out sphere with anouter boundary and inversely oriented boundary of the inner cavity.
If the inner boundary were not inversely oriented then the core has
w = 2 > 1 and is considered twice inside.

In general, the extraction function f takes as input a winding number vector w = [w1 , . . . , wn ] corresponding to the winding numberof each input mesh at the points of a given cell of the space partition
defined by the mesh arrangement. The function f returns “true” if
a region with this winding number vector is to be included in the
output solid, and “false” otherwise.

S EAMS: A mesh with a combinatorially open boundary does notnecessarily imply that it is not PWN. Open boundaries are permissible so long as they meet up geometrically along seams.

For example, to implement n-way union, one would provide:
inside any


}|
{
z
funion (w) = true if ∃ i | wi 6= 0,
false otherwise.

E XACT: As discussed in Section 5, the vertices of PWN mesh may
be defined with rational coordinates, not just floating-point.
S ELF - INTER .:A PWN mesh may have structured selfintersections. For example, two overlapping spheres constitute avalid PWN mesh. Similarly, a vertex-displacement of a PWN meshwithout seams is also a PWN mesh regardless of any incurred selfintersections [Sacht et al. 2013]. We will say that a mesh is freeof self-intersections if any two geometric triangles of the mesh intersect only over a (combinatorially) shared edge or vertex, or are
combinatorially identical. We allow combinatorially duplicate triangles and exploit these to cope with degenerate configurations.
D EGENERACIES: Geometrically degenerate triangles (zero area)
do not affect the winding number. This implies that onedimensional “needles” will be ignored entirely. One can formalize
ignoring degenerate triangles as reconstructing a discretization of
the discontinuity sheets of the winding number field. For boolean
set operations, we operate on open-set interiors without boundary.
Verification We may verify whether a mesh Ai is PWN by resolving self-intersections (see Section 5.1) and then checking thatthe total signed incidence of every edge in the result is zero. For anyedge e = {i, j} an oriented triangle f = {i, j, k} contributes +1

to the total signed incidence of e. An oppositely oriented triangle
g = {j, i, `} contributes −1.

Variadic extraction function(3)

When n = 1, this function will identify a mesh’s self-union. Similarly for n-way intersection:
inside allz }| {
fintersect (w) = true if wi 6= 0 ∀ i,(4)
false otherwise.
Some extractions are asymmetric, e.g., subtraction (A1 \ A2 ):fminus (w) = w1 6= 0 and w2 = 0| {z }| {z }
inside of A1(5)

outside of A2

One can also design more esoteric functions, such as extracting allparts of space inside at least two of the inputs:(
true if ∃ i and j 6= i | wi , wj 6= 0,
fmin-2 (w) =(6)false otherwise.
Changing the two-sided inequalities above (e.g., wi 6= 0) to singlesided inequalities (e.g., wi > 0) results in orientation-sensitive op-erations [Campen & Kobbelt 2010a; Jacobson et al. 2013]. Orientation sensitivity is useful in some cases, where inversion moreintuitively denotes the exterior or void space of an input shape.

3.3

Solid meshes

Our algorithm’s output meshes belong to a special subclass of PWN
meshes that we call solid meshes. Solid meshes are free of selfintersections, degenerate triangles or duplicate triangles, and theirgeneralized winding number field is either zero or one.
Note that even if the input meshes A1
and A2 are manifold polyhedra, the output of C = A1 ∪ A2 may be a nonmanifold solid mesh (e.g., if A1 bounds
the unit cube and A2 bounds the unit
cube offset by (1, 1, 0) then C will contain a non-manifold edge where A1 andA2 “kiss”, see inset).

4

Overview

Before worrying about details of themethod, we review the key aspects ofeach stage. For now, we consider theusual binary boolean operations on two
meshes A and B. The insets in this section illustrate the stages of the computation of the asymmetric difference A \ B.
In thefirst stage, we resolve all intersections
between input meshes using exact arithmetic. We add new triangles by subdividing the inputs so that all intersections
occur exactly at edges and vertices. All
refined triangles retain references to theoriginal triangles of A and B.

Arrangement construction

In the second stage, we determine adjacency information between cells defining a space partition. We organize themesh resulting from resolving intersections in the first stage into patches of triangles connected by manifold edges. By
definition, patches are incident to eachother along non-manifold mesh edges. Two cells are adjacent via
a shared oriented boundary patch. Two patches incident on thesame non-manifold edge may bound the same cell. We determine
the patch-cell relations by sorting facets from all incident patches
around this edge. In this way, we determine the cell adjacency foreach connected component of the adjacency graph of patches. To
ensure correct cell adjacency of nested components, we identify
a boundary facet of the ambient cell surrounding each component
and determine if it is contained in an interior (non-ambient) cell ofanother component, via point location (see Section 5.5.1). Aftermerging the cell adjacency graphs across connected components,
there remains a single ambient cell outside of all components.
In the third stage, we assign winding
numbers with respect to A and B toeach cell, w = [wA , wB ]. Having constructed the cell-patch adjacency information in the previous stage, this step ispurely combinatorial. The ambient “0”cell is defined to [0, 0].

Figure 8: Alternative extractions (A ∪ B, A ∩ B, B \ A) share thesame input mesh arrangement.
Remaining cells are labeled via a traversal of the cells: we add +1 or −1 tothe winding number of the originatingmesh of the patch crossed between cells
depending on the patch orientation. For
example, consider an unlabeled cell adjacent to a cell labeled [a, b] via a patch
originating from input A. That unlabeled cell will receive [a + 1, b] if crossing into the oriented patch or [a − 1, b]
if crossing out.
We identify desired output cells purely by their assigned winding numbers. For the difference A \ B, we collect the cells labeled [1, 0]. We return the facets of thepatches separating these desired cells from undesired cells, reversing orientations if necessary. Different extractions reuse the sameintersection resolution, cell adjacency, and winding-number labeling of the first stage (see Figure 8).
Extracting the result

5

Algorithms

In this section, we consider in detail the algorithms for each of thesteps overviewed in Section 4. We will break each stage into core
subroutines. For each subroutine, we will provide preconditions on
its input and postconditions on its output. We start with specifying
preconditions and postconditions of our method as a whole.
The method accepts as input a sequence of PWN
meshes, and an extraction function whose arguments correspond tothe mesh sequence. The mesh vertex coordinates are assumed to berational coordinates, a property we call (E XACT ). We review exactness in the Appendix. In accepting as input the broader class ofexact coordinates rather than floating point values, we accommodate upstream operations, whose output is exact.
Preconditions

The output of our algorithm is guaranteed to be
an exact solid mesh. As such it is a valid input to a downstream
application of our own algorithm or another module in CGAL. Observe that while our input preconditions permit self-intersections,and co-planar/degenerate/duplicate facets—by design, these will
never occur in our output.
Postconditions

5.1

Intersection resolution

The first stage of arrangement construction resolves all triangletriangle intersections, enriching the mesh combinatorics so that allintersections are exactly represented by shared vertices and edges.
We consider all input meshes as a single mesh, i.e., we make no
distinction between intersections and self-intersections.

Preconditions

The input is an exact PWN mesh A.

The output is an exact PWN mesh free of selfintersections, co-incident vertices, and degenerate triangles, inducing exactly the same winding number field as the input mesh.
Postconditions

Self-intersection resolution consists of four steps: discard exactly zero area input triangles as they do not affect the winding number, compute the intersection between every pair of triangles, conduct a constrained Delaunay triangulation for every coplanar cluster of intersections, and extract and replicate subtriangles from each triangle’s cluster’s triangulation.

Algorithm

For now we set aside conservative culling for performance acceleration. We consider all pairs of triangles a and b in A. The intersection intersect(a, b) between these triangles can be one of thefollowing four cases: empty, a single point, a line segment, or aconvex polygon (see Figure 9). This intersection must be computed
exactly, therefore intersect(a, b) commutes.

Figure 11: We retain the relationship between the output triangles
and the inputs. Because of our exactness, attributes like texture
coordinates are losslessly maintained.
the CDT cluster those sub-triangles {t1 , t2 , . . .} whose three vertices, according to exact 2D predicates, are not strictly outside t.
We clone each sub-triangle ti and orient it to match t, again using exact 2D predicates. In order to label winding number vectors(Section 5.3), the cloned oriented subtriangle stores a reference to t.
This reference is also useful for applications requiring interpolationof texture coordinates, colors, or other attributes onto the boolean
output mesh (see Figure 11).
We clean up by purging geometrically duplicate vertices. As we areusing exact vertex representation, this can be done efficiently using
lexicographical sorting and unique entry extraction from the list ofvertices. The result is a possibly non-manifold mesh with possible
duplicate triangles, but no self-intersections.

Figure 9: Blue triangles intersect the green at a point, segments
and a polygon (left, 3D). Subdivided constraints reduce to coplanar
points and segments (middle, 2D). The original green triangle isreplaced with replications of the green triangles of a constrained
Delaunay triangulation (CDT) of the coplanar cluster (right, 2D).
Next we replace each input triangle with a triangulation containing
those elements resulting from the intersections.
Previous methods construct this triangulation independently foreach triangle [Jacobson et al. 2013; Attene 2014; Barki et al. 2015],
but this approach may introduce inconsistencies between overlapping triangles due to non-general position configurations (see Figure 10). Inconsistent triangulations of coplanar intersections resultin violating the precondition of the following stage that the mesh is
free of intersections as defined in Section 3.1.

Duplicate triangles need to be retained at this stage, as their removal
requires knowledge of the extraction function.
The output mesh has exactly the same
winding number field as the input. This
immediately follows from the fact thatin the result all non-zero area trianglesof the original meshes are retained, possibly in the subdivided form as a resultof intersection resolution. One aspect of
ensuring this is cloning sub-triangles at
co-planar intersections. In the inset figure, the orange and blue shapes share aside. If only one set of faces is kept in
the output the result is not a PWN mesh.

5.2

Partitioning space into cells

The second stage explicitly constructs a
set of cells. We define a cell as a region bound by the union oforiented patches forming a closed manifold mesh with no selfintersections. The set of cells forms a space partition.
Figure 10: Two overlapping, co-planar right triangles (left) admit
multiple constrained Delaunay triangulations (CDTs). Independent triangulation could lead to inconsistency (middle and right).
Instead, we gather clusters of triangles connected via non-trivial
co-planar intersections (i.e., intersections resulting in convex polygons). By construction all triangles in a cluster share the same
supporting plane. We compute a 2D constrained Delaunay triangulation (CDT) of the convex hull of each cluster. The original
constraints collected from triangles in the cluster are the points,
segments and polygons resulting from intersections with all othertriangles in the input mesh A, as well as the vertex points and edgesof the cluster triangles themselves. We further subdivide segment
constraints so that all intersections are resolved as constraint Steiner
points. Finally, we compute the CDT of the convex hull of thesepoints and segments; no additional Steiner vertices are required.
With CDTs constructed for each cluster, we iterate over each original triangle t to collect its respective subdivisions. We select among

Each patch is a subset of triangles of the input mesh and inherits
their orientation; the patch is a maximal connected set of faces withall edges shared by two faces from the set being manifold. The
condition implies that a boundary edge of a patch (if it exists) is a
non-manifold junction with neighboring patches.
Boundary patches of a cell may be geometrically coplanar, producing zero-volume cells; by the absence of self-intersections, such
patches must consist of single triangles sharing the same vertices.
The input mesh is
PWN without degenerate triangles, selfintersections, and co-incident vertices.
Preconditions

The output is a bipartite directed graph encoding of cellpatch incidences. Each patch node has
one incoming and one outgoing edge tocell nodes, representing the volumetric
Postconditions

regions on the positive and negative sides of the (oriented) patch, respectively, which we call above and below cells. In the following,
we refer to “patch” (the combinatorial and geometric data structure) and “patch node” (the bipartite graph node) interchangeably,
and likewise for cells.

C − (ri+1 ) are unassigned, we assign them both to a newly createdcell node. If only one is unassigned, we set it to the other cell node.
If both are assigned to distinct cell nodes, we merge the two nodes
in the bipartite graph. After visiting all patch-patch incidences, thegraph for one connected component is complete.

The output also includes mutual references between patches and input mesh triangles, i.e., each patch node contains a list of triangles,
and each triangle has a pointer to the patch it is contained in.

After completing each connected component, it remains to merge the bipartite
graphs. In particular, the ambient cell ofa nested component must be equated tothe corresponding internal cell of the enclosing component (see inset). As a special case, the ambient cells of all nonnested components must be equated.

Geometrically, the cells cover all R3 .
Some cells will have zero geometric volume. These cells are always bound byexactly two clones of the same geometric triangle: i.e., two patches, each with
one triangle. There may be many suchdegenerate cells stacked on the same
multiply cloned triangle.

As detailed in Section 5.5.2, for each component Ki , we identify its
ambient cell A(Ki ).

Our cell partitioning algorithm first separates the input
mesh into connected components of triangles. Two triangles are
considered connected if and only if they share an edge.

We now iterate over each component, determine whether it isnested, and find the cell to which its ambient cell should be equal.
Visiting component Ki , we select an arbitrary point p ∈ Ki . We
build a set of candidates E = {Kj | j 6= i, p 6∈ A(Kj )} consistingof every component Kj 6= Ki whose ambient cell does not containp. To determine whether a cell contains p, we use point location(Section 5.5.1): given a point p and component Kj find the containing cell L(Kj , p) and distance from p to Kj .

For each such connected component, we construct a cell-patch
graph independently of the other components. First, we cluster triangles into patches. Starting with any unassigned triangle we growa new patch traversing across manifold edges until the boundary of
the patch is either empty or consists only of non-manifold edges.

If the set E is not empty, then Ki is a nested component. Among thecandidates E, we select the one (and only) component Kj closest
to p. We merge the bipartite graph nodes A(Ki ) and L(Kj , p),equating the ambient cell of the nested component to the interiorcell of its enclosing component, respectively.

During clustering, we record, for each non-manifold edge, its incident patches.

If the set E is empty, then Ki is not a nested component. In thiscase, we equate its ambient cell A(Ki ) with the universal ambient
cell C0 , defined as the cell containing “all points at infinity”, bymerging these two nodes in the bipartite graph.

When the input to this stage is the resolved intersection (see Section 5.1) of n piecewise-constant winding number inducing meshes,
then the output is denoted a valid cell-patch data structure.
Algorithm

The adjacency between patches is encoded as a matrix A, settingA(p, q) = e. which means that patch p is incident to patch q sharing with it a representative non-manifold edge e. Incident patches
p and q may share multiple non-manifold edges, and the choice ofrepresentative A(p, q) is arbitrary.
We now construct the bipartite graph of cells and patches encoding
the volumetric partition: while the patches have already been established above, it remains to construct the cells and to add, for eachpatch p, one outgoing edge to C ↑ (p) and one incoming edge fromC ↓ (p), the cells above and below p, respectively.
We will traverse all patch-patch incidences in arbitrary order.
When visiting an incident pair (p, q),
we retrieve the representative edge e =A(p, q). As detailed in Section 5.5.3,
we sort the e-incident patches cyclicallyaround e, with respect to p’s orientation.
Suppose the sort results in the ordering[p, r1 , r2 , . . . , rk ], so that r1 and rk are
immediately “above” and “below” p, respectively (see inset).

After we have processed all components, the ambient cell of eachnested component has been equated with the interior cell of its enclosing component, and the ambient cell of all non-nested components is C0 . The bipartite graph is now connected; each patch isincident to two cells in a consistent manner.

5.3

Winding number labelingIn the third stage, we compute the winding number of each cell withrespect to each input mesh.
Preconditions The input is a valid cell-patch graph. The universal ambient cell is seeded with a known winding number vector;
by default w = [0, . . . , 0], signifying that infinity lies outside all
shapes. Only the combinatorial (not geometric) aspects of the input are considered by this algorithm: instead of computing winding
numbers geometrically, we use property that the winding number
changes by 1 or −1, whenever a surface is crossed, and thus can becomputed by propagation along the cell-patch bipartite graph.

If we think of this sorting order as “upward,” then each patch’s
own orientation is either consistent or inconsistent with the sortingorder (e.g. in inset, r2 is inconsistent). Let C + (ri ) ≡ C ↑ (ri ) andC − (ri ) ≡ C ↓ (ri ) if patch ri is oriented consistently with the sort,otherwise let C + (ri ) ≡ C ↓ (ri ) and C − (ri ) ≡ C ↑ (ri ).

Postconditions

We now propagate cell assignments by iterating over each consecutive pair of e-incident patches in order, beginning with (p, r1 ) andending with (rk , p). Visiting the pair (ri , ri+1 ), our task is to identify the cell references C + (ri ) ≡ C − (ri+1 ). If both C + (ri ) andThis process assigns winding numbers to zero-volume cells, formedby duplicate triangles; although no points have this winding number, this records what interior points would have if the cell boundaries were separated consistently with patch ordering.

The output is a valid cell-patch data structurewith consistently labeled winding number vector for each cell.
Neighboring cells will differ in winding number vector by exactly+1 or −1 in a single entry corresponding to the originating meshof the patch between them, signed according to its orientation.

Algorithm Given cell C with a known winding number vector
wc , we assign the winding number vector of neighboring cells via
breadth first traversal. For each oriented patch p separating cell
C from neighbor cell N , if the winding number vector wn is still
unknown we set it to wc adjusted to account for crossing patch p,
originating from mesh Ai :wn ← wc + sp [δi1 . . . δin ],(7)

where sp is +1 if cell C lies above p and N below and −1 if viceversa, and δij is Kronecker’s delta. We then add cell N to the queueof cells to process later. When the queue is empty, all cells havebeen labeled, and the algorithm has completed.
Complements By default, all input meshes Ai are assumed torepresent bounded solids. Under this convention, the winding number at infinity is zero, wj (∞) = 0. Winding numbers elegantly
handle complements by subtraction from 1. If Aj is the complement of Ai = Aci , thenwj = 1 − |wi |or

wj = 1 − wi ,(8)depending on whether the complement operator is orientationinsensitive or -sensitive, respectively (see inset).1

As a consequence, if Aj represents the unbounded complement ofsome bounded solid, then the seeded winding number vector at infinity should be wj (∞) = 1.
In this way the winding number elegantly captures set identities. In
particular, we produce exactly the same result for A\B and A∩Bc .

5.4

Operation result extraction

With the arrangement data structure constructed, we may perform
arbitrary extraction operations. We extract the triangulated boundary of all cells for which an extraction function f is true.
Inputs are a valid cell-patch data structure withconsistently labeled winding number vectors and a function f (w)returning true or false for a given winding number vector w. As in
the previous stage, this stage makes use only of combinatorial (notgeometric) aspects of the input.
Preconditions

Postconditions

The output is a solid mesh.

We flag all cells that pass f , collecting all patches separating a flagged cell from an unflagged cell, and then collecting all
triangles of those patches, flipping the orientation of triangles frompatches with a flagged cell above and unflagged cell below. We
then purge possible boundaries arising from zero-volume symbolic
cells. Since these always occur as perfect combinatorial duplicatesof single triangles, we need only remove all triangles with zero total
signed occurrence: sum of +1 if oriented i, j, k and −1 if k, j, i.

Algorithm

5.5Core low-level subroutines

The robustness of our method rests on the correctness of several
core low-level subroutines.
1 In the orientation-insensitive case, the complement of the complement

results in taking the absolute value of the winding number w, preserving
orientation-insensitive insideness defined as w 6= 0.

5.5.1

Point location

We are asked to locate in which cell a given point lies. This is a
special case of the fundamental point location problem in computational geometry. We take special care to solve this problem robustly
and in the presence of zero-volume cells (e.g., due to resolved coplanar intersections).
Inputs are a query point q ∈ R3 and a valid cellpatch data structure. The query must not lie exactly on any patch.
Preconditions

Postcondition

The output is the unique cell containing the query.

Searching over all triangles, we find a triangle t containing the point c on the input mesh closest to the query point q.
Point c lies either exactly at a vertex v of t, or else along an edgee, or else within the interior of t (but not on its boundary). Noneof these cases is trivial. The vertex v or edge e could be a nonmanifold junction of many cells, and the triangle t could lie deep ina “stack” of zero-volume cells due to duplicated faces.

Algorithm

In fact, only at edges can we robustly determine the symbolic andgeometric cell arrangement. The cyclic ordering of cells incident
on an edge e is consistent, and since q does not lie on the input
mesh it must lie in one of the incident cells. We insert a dummy
facet connecting e and q into the sorted list of facets incident on e.
The next facet after the dummy (or previous facet before) must bepart of a patch bounding the cell containing q.
The ambiguity in the case where c lies within the triangle t arisesin the presence of duplicates of the triangle t. Since all duplicatesshare the same three edges, we choose one arbitrarily as the sortingedge e, and insert a dummy as in the edge case above.
If the closest point c lies at a vertex v, we identify a good sortingedge e (i.e., on the convex hull of v and its vertex neighbors) andagain insert a dummy as above. Identifying the containing cell ofa query whose point of closest approach is a vertex will also arise
when identifying the ambient cell of a component. In this case, we
project edges incident on v onto the plane formed by q − c and any
orthogonal vector (rather than the xy plane) and then follow the restof the ambient cell identification algorithm in Section 5.5.2.
5.5.2

Ambient cell identification

In this subroutine, we identify the ambient cell (containing all
points at infinity) of a given mesh.
Inputs are a (non-manifold) self-intersection-freetriangle mesh and corresponding cell-patch data structure.
Preconditions

The output is a facet guaranteed to contain anouter vertex (a vertex on the convex hull of the vertices) and participate in a patch forming part of the boundary of the ambient cell of
this mesh containing all points at infinity.
Postconditions

We could solve this problem by identifying the cell
containing some arbitrary far away point using the point location
algorithm of Section 5.5.1. However, we enjoy the performance
benefits of avoiding closest point computation by choosing a query
point with a vertex as its known closest point.
Algorithm

We locate a vertex v with the maximum x-coordinate magnitude,
breaking ties arbitrarily. It follows immediately that q = v +(1, 0, 0) lies in the desired ambient cell and that v is the point ofclosest approach of q to the input mesh.

Figure 13: Consider five co-planar facets (lines) incident on twooriented edges (dots). Using the edge orientation to sign these indices during sorting ensures that orderings from either edge areconsistent (left). Otherwise the sort is the same regardless of theedge orientation leading to inconsistent ordering.
Figure 12: Consider a tetrahedron (orange) inside a wedge (blue),connected at a non-manifold vertex. One face normal of the interior
tetrahedron (orange) has a larger x-component than the normals of
exterior triangles (blue) incident on this non-manifold outer vertex.
We will identify the ambient cell by finding a facet incident on v
that is part of a patch on the ambient cell’s boundary. To find an
incident outer facet, we first select an edge incident on this vertex
that also lies on the convex hull. Then we sort facets cyclicallyaround this edge and select one of the two facets that are part ofambient-cell boundary patches.
Sorting facets around an edge is discussed in Section 5.5.3, so it
remains to identify an edge of the convex hull edge on the vertex
with maximal x-coordinate. We rely on our exact representationof the input mesh and the ability to determine predicates exactly(e.g., is a point below, on, or above a plane?). We sort incident
edges with respect to their projection on the xy-plane. We selectthe edge whose projected edge-vector e = (ex , ey ) is most orthogonal to the x-axis. Our particular exact representation kernel allows
construction of quotients (but not square roots), so we identify theedge with maximum slope as a line function of x: that is, according
to |ey /ex |. We may break ties arbitrarily because all edges withmaximum slope must lie on the convex hull.

tions). It is assumed that if two triangles are coplanar, then either
they intersect only along the edge {i, j} (their dihedral angle is180◦ ) or they are geometrically identical (same third vertex position and their dihedral angle is 0◦ ).
This subroutine outputs a sorted (clockwise) ordering of the triangles, looking down the edge {i, j}. Geometrically distinct triangles are sorted cyclically according to their dihedral angle with the first triangle t1 . Duplicate triangles—without
loss of generality all are {i, j, k}—are sorted consistently in thesense that their relative ordering is maintained when sorting around{i, j}, {j, k}, or {k, i} and their ordering is reversed when sorting
around {j, i}, {k, j}, or {i, k}.
Postconditions

Let each triangle t be positively incident on {i, j} if
t = {i, j, k} and otherwise negatively (i.e., if t = {k, j, i}).
Algorithm

Let pk refer to the vertex position of triangle tk ’s third “flap” vertex not lying on the shared edge {i, j}. Our recursive divide-andconquer algorithm begins by selecting a starting triangle t = t0
and sorting each other triangle tk into one of four groups, based onwhether pk lies (1) co-planar with t0 and on the same side of {i, j}
as p0 , (2) co-planar with t0 , and on the opposite side of {i, j} as
p0 , (3) below the plane of t0 , (4) above the plane of t0 .

Remark Attene [2014] proceeds in a similar way by finding a maximal x-coordinate vertex and then chooses the incident triangle
whose normal has the largest magnitude x-component. This criterion cannot be applied if the vertex is non-manifold: an inner“flap” might have a more outward-pointing normal than the trueouter facets. For a concrete counterexample, consider extruding
the triangle {(0, 0), (1, 1), (0, 2)} two units in the z-direction, thenmove the lower-right corner to (2, 1, 0) and add a inner tetrahedron
connecting that vertex to the top-left corners and any interior vertex
floating below (see Figure 12).

We sort within groups (1) and (2) by simulating simplicity à la[Edelsbrunner & Mücke 1990]. Duplicate triangles are sorted according to their uniquely assigned index εk . To ensure that this ordering is consistent and not erroneously reversed when viewed froma different edge incident on the same replicated triangles, we signthese indices based on the signed incidence of each triangle withrespect to {i, j} (see Figure 13). This symbolic perturbation will
differ depending on input indices, but is always consistent. Only
the ordering of zero-volume cells are effected, so different orderings will always produce the same geometric result.

5.5.3

Triangles in groups (3) and (4) are sorted by recursive calls. The
complete output is then simply the merger of the four sorted groups.

Cyclical sort triangles about a common edge

The cell partitioning, point location, and ambient cell identification
subroutines depend on the ability to sort triangles about a commonedge robustly. We sort only at one representative edge between
incident patches, rather than at every edge of the triangulation (cf.
[Attene 2014; Barki et al. 2015]).
Sorting triangles around a common edge is misleadingly innocuous. This subroutine must (and will) ensure consistent ordering ofexactly duplicate triangles (e.g., resulting from resolved co-planar
input triangles) and geometrically correct ordering of numericallynearly co-planar triangles.
The input is a set of m non-degenerate triangles
t1 , . . . , tm incident on a mutual (non-degenerate) edge {i, j}. Eachtriangle tk is endowed with a globally assigned index εk (e.g., its
index in the non-manifold output mesh after resolving all intersecPreconditions

6

Implementation

We implemented the algorithm in C++ utilizing the exact arithmetic kernel of the popular CGAL library. We specifically use
its subroutines for: exact testing and construction triangle-triangle
intersections; 2D constrained Delaunay tessellation (CDT); pointtriangle closest point queries and point-plane predicates. We found
CGAL’s CDT implementation to be robust on all examples if constraints are subdivided at intersections as a preprocess.
We also use CGAL’s built-in bounding-box-based spatial acceleration for collecting a list of candidate triangle-triangle intersections.
We further accelerate the exact triangle-triangle intersection detection and construction by processing candidates in parallel. Due tothe reference counting employed by CGAL’s deferred evaluation

Figure 14: We conduct extensive evaluation of the robustness of our method on a dataset of 10,000 popular real-world models.
exact number type (CGAL::Lazy_exact_nt), seemingly readonly simultaneous access of triangle data is unsafe. Fortunately intersection detection and construction is compute-bound, so despite
placing mutex locks around every mesh vertex leads to parallelism
performance gains.
We also use CGAL’s axis-aligned bounding-box hierarchy forpoint to triangle-soup closest point querying. We further accelerate
the point location in Section 5.5.1 by culling points entirely outsideof the bounding box of a component (the query point must then lie
in that component’s ambient cell).

6.1

Converting to floating-point

While input meshes with floating-point vertex position coordinates
losslessly convert to our exact representation, the reverse is not true
about our output exact meshes. Naively rounding a solid exact meshto floating point may result in a non-solid mesh due to newly introduced self-intersections. This occurs in 2.19% of our output meshes
in the Thingi10K dataset.
This problem is known as vertex rounding. Without allowing subdivision of facets and insertion of new vertices, this problem is NPhard [Milenkovic & Nackman 1990]. Allowing for re-triangulation,a robust—albeit slow and complicated—solution to this problem
exists in theory [Fortune 1997].
To fit into floating-point pipelines and fairly compare to previous
methods producing floating-point output meshes (e.g., [Bernstein
2013; Attene 2014; Douze et al. 2015]), we propose a heuristic forrounding our exact output meshes to floating-point. Our heuristic is
related to the method proposed in [Sacks & Milenkovic 2014].
Preconditions We assume the input to be a solid triangle-meshwith exact coordinates.

Though we can make no guarantees of convergence, our exact method equipped with this rounding heuristic
successfully finds self-intersection free floating-point meshes for
99.95% of the dataset. Otherwise, we can only claim the output toremain a PWN mesh.
Postconditions

Given a solid mesh with exact vertices, we iteratively
apply the following steps: (1) round all vertices to double precision floating-point coordinates, (2) find all triangles participating
in self-intersections (if none, then return), (3) round all vertices ofthese triangles to single precision floating-point coordinates, and(4) compute the self-union of the resulting mesh.

Heuristic

Figure 15: Real-world meshes with real-world problems.

7

Experiments and results

Constructed or procedurally generated examples may help investigate corner cases, but do not necessarily report how robustly analgorithm will perform in practice. To this end, we gather a datasetof 10,000 meshes from “the wild,” and test our method and previous works against it. Considering these meshes as a representative
sampling of a general population of meshes encountered in practice,
we evaluate the restrictiveness of preconditions and the robustnessof claimed postconditions across methods.

7.1

Thingi10K dataset

The Thingi10K dataset contains
the first 10,000 meshes of “Featured” models on thingiverse.com, a
popular shape repository. These models are heavily biased towardmodels designed by amateurs or semi-professionals for 3D printing(though there is no official restrictive policy). We therefore interpret these models as a representative sampling of the population ofmeshes intended to model a solid 3D object.

Contents and methodology

Each “thing” featured on thingiverse.com may contain several
distinct mesh files. We collected the first 2011 things, totalling 10,000 meshes (see Figure 14). All things have free licenses (GPL, LGPL, CC, BSD, or public domain). The original meshes came in a biased variety of file formats: 9956
.stl, 42 .obj, one .off, and one .ply. The vast majority of meshes have single-precision vertex-coordinates. Since.stl files store triangle streams rather than meshes, we
immediately merge exactly duplicate
corners. The number of faces in eachmesh follows a log-normal distribution
with geometric mean and standard deviation 5077.6 ± 8.5 (see inset).

Figure 16: Previous methods frequently failed to produce an output without self-intersections, without open boundaries, and withpiecewise-constant winding number.

Figure 17: The performance of our method is competitive with existing floating-point methods and faster than the state-of-the-art exact method [CGAL 2015]. Geometric means given for each method.
Unequal histogram areas correspond to success rate.

Of the10,000 meshes, 8616 meet our PWN
precondition. Of these, 5113 are solid;of those 4963 are manifold polyhedra.
Comparing preconditions

The 10,000 meshes exhibit a variety of
typical problematic cases: open boundaries, self-intersections, non-manifold
elements, multiple components, etc. (see Figure 15). Among the4524 meshes containing self-intersections, 3082 contain coplanar
self-intersections. This quantifies an approximation of the fractionof models deviating from the general positioning assumption.
Many “problematic” meshes seem to result from modeling with
self-intersections (see Figure 2) and overlapping, independently
modeled components or from previous failed boolean operations.

7.2

Testing self-union

Assuming each mesh in the Thingi10K dataset to represent an intended solid, we compare extracting a valid boundary of this solidwith available implementations of five previous works: “CGAL”[CGAL 2015], “Carve” [CARVE 2014], “Cork” [Bernstein 2013],“QuickCSG” [Douze et al. 2015], and “Attene” [Attene 2014].
We emphasize that this experimental comparison reflects both algorithmic limitations and implementation deficiencies, so a different
implementation of any given method could potentially and performbetter. The no-longer-maintained implementation of [Bernstein &Fussell 2009] failed on most examples. Similarly, the web-service
implementation of [Campen & Kobbelt 2010a] failed to produce a
result roughly 40% of the time. We are unable to obtain implementations or outputs for other methods (e.g., [Barki et al. 2015]).
We limit our comparison to the 8616 PWN meshes. Of these, only
3413 contain self-intersections. Nonetheless, we consider all 8616
PWN meshes as implementations relying on internal rounding (e.g.
[Attene 2014; Bernstein 2013]) often also stumble on nearly selfintersecting meshes.
Attene’s mesh repair method computes the outer hull, rather than
the self-union [2014]. The other implementations do not provide
an explicit API for conducting self-union, so we intersect the input
model with its conservative bounding box.

Figure 18: The main bottleneck of our algorithm is triangletriangle intersection resolution. Within this major stage, the intersection detection and exact construction dominate.
Self-union should output a solid
mesh. We report whether a method successfully produced any output and if it met certain necessary postconditions. Strictly testingsolidity requires a correct implementation of cyclic facet ordering
around a non-manifold edge to determine that all incident cells arealternating zero/one winding number. Absent trusted third-party
code, we test necessary (but not sufficient) conditions: lack of selfintersections, open boundaries, and non-zero total signed incidenceedges. Such meshes form a strict subclass of PWN meshes, but asuperclass of solids. Our exact method succeeds with 100% success across all criteria (see Figure 16). Previous methods fall shortin at least one criteria. This unique success places our exact method
robustly into the exact geometry pipeline.
Comparing postconditions

In a floating-point context, our method also out-performs all others. Our heuristic for converting our exact outputs to floating-point
meshes in Section 6.1 succeeds in removing new self-intersections
all but five cases out of the 8616. These meshes fail to converge
after 20 iterations. Rates of closedness and total signed edgeincidence are—by construction—maintained at 100%.
The specific causes of failure of the previous methods are difficult todetermine. We can identify robustness flaws associated with characteristics of the input meshes. Methods assuming general positioning or resorting to numerical perturbation [Bernstein 2013; Douzeet al. 2015] will struggle in the presence of coplanar intersections.

Figure 19: Self-intersections confuse per-vertex ambient occlusionand sharp-line detection. Rendering the outerhull ameliorates this.

Figure 21: Numerical perturbation can produce spurious artifacts.

Figure 20: David emerges from a block by repeatedly subtracting
the Minkowski sum of a drill bit along piecewise-linear paths.
Figure 22: We reproduce and exhaustively expand the pairwise
testing in [Barki et al. 2015] (see supplemental material).
Attene assumes accurate floating-point normals during self-intersection culling and outerhull extraction [2014], but inputs may containdegenerate or nearly degenerate triangles with
untrustworthy normals. The inset highlights
self-intersections (orange) and an open boundary (red) on a problematic output of Attene’s.

The Minkowski sum of a solid mesh A along
a segment {s, d} can be computed as the union of A at s, A at d,and the union of all prisms formed by triangles of A along {s, d}:!
[[
A + {s, d} =A + s, A + d,t + {s, d} . (9)Minkowski sumst∈A

We collected timing information across the 8616 self-unions for our
method and four others (CGAL, Carve,
QuickCSG, Cork) locally on a machine with
an 8-core Intel Xeon 3GHz processor with
16GB of memory.2 The violin histograms of
running timings in Figure 17 show that whileours is not the fastest, it is competitive.
Performance

The Thingi10K dataset also provides means to further examine theperformance of our individual subroutines. The profile in Figure 18
reveals that resolving intersections is the dominating bottleneck.

7.3

General discussion

Outer hull The outer hull of an input triangle mesh is definedas those triangles reachable from infinity by some (possibly nonstraight) path that does not intersect the mesh [Campen & Kobbelt
2010b; Attene 2014]. In general, the outer hull cannot be categorized in terms of the winding number: boundaries with inner hollow
cavities with zero winding number are not part of the outer hull. For
some applications, retaining these inner cavities is crucial (see Figure 5). For other applications, such as rendering, the outer hull may
be appropriate and desired (see Figure 19). We can easily adapt
our algorithm to compute outer hulls. We construct cell partition
according to Section 5.2, find the ambient cell according to Section 5.5.2, and simply extract its boundary as per Section 5.4.
2 Attene provided results independently.

Explicitly computing the union of all triangular prisms via our meshboolean algorithm would produce the correct result but after too
much unnecessary computation: most neighboring prisms are exact duplicates. We cull the union of prisms with a pre-process, removing all facets with zero total signed occurrence, replacing all
instances of facets with ±2k total signed occurrences with ±k positive/negative clones. This proof-of-concept inherits the robustnessof our method, but is likely suboptimal in terms of performancecompared to specialized methods [Campen & Kobbelt 2010b]. In
Figure 20, we simulate a CNC-milling tool.
A traditional test for a boolean
algorithm is to select two meshes, randomly rotate them, then conduct a binary operation (union, intersection, difference, etc.) andinvestigate the result for artifacts or errors. This type of testing encourages the general positioning assumption and may give a false
sense of robustness in cases with coplanar intersections and exact
co-incidences. An extreme case is taking the intersection of an object with a clone of itself (see Figure 21). Methods based on numerical perturbation, such as [Bernstein 2013], panic in the presence ofso many co-planar intersections.
Traditional binary boolean tests

For completeness, we reproduce and expand the testing in [Barkiet al. 2015]. Barki et al. compute the union and intersection for 22
pairs of meshes from a collection of 26 standard computer graphics
meshes (Armadillo, Dino, etc.). We exhaustively compute union,intersection, and both asymmetric differences for all pairs (see Figure 22). All 4(26(26 + 1))/2 = 1404 tests result in valid solidmeshes. For two of these meshes, Barki et al. also compute theunion and intersection of the mesh and a clone rotated by randomFigure 26: Disconnected stars and overlapping spike components
correctly unites with a nested turtle (slice view above).
Figure 23: We reproduce and expand upon the A-union-rotated-A
style test in [Barki et al. 2015].

Figure 24: The letters are separated overlapping components on
the cryptex. Exact union reveals disjoint rings, shown in different
colors, but small tolerances between parts cause inexact methods
to merge over zealously.

Figure 27: We reproduce the carving example of [Bernstein &Fussell 2009] by subtracting 10,000 dodecahedra from a box.
An interesting and challenging application of booleans is to “undo”
boolean subtractions given only the argument and the result, produced by an unknown boolean implementation (see Figure 28).

Figure 25: A popular commercial app produces three different results (presumably due to randomization), yet all are incorrect.
rotation. For all 26 meshes, we compute the union of the mesh and
10 clones rotated by π/10, 2π/10, . . . , π about the same axis (seeFigure 23). All our results are valid solid meshes.
The robustness of several previous works rely on the assumptionthat input vertices lie on a regular grid [Campen & Kobbelt 2010a]
or at general positions [Bernstein 2013; Douze et al. 2015]. However, input rounding or perturbation—no matter how subtle—may
introduce unnecessary intersections that merge disjoint components(see Figure 24) or cause numerical problems (see Figure 21). In
contrast, the exact nature of our approach allows us to only resolveintersections already present in the inputs.
Others (e.g., [Barki et al. 2015]) have demonstrated robustness issues with boolean implementations in commercial software suchas M AYA. We add to this by comparing to Trimble’s S KETCH U P
P RO. S KETCH U P P RO consistently fails to intersect four randomly
rotated icosahedra (see Figure 25). We also attempted to union eachof the 26 models of [Barki et al. 2015] with a clone rotated by π/10
(a simplified version of Figure 23). After a day of computation,
only 12 produced an output, and none were without flaws (all werecombinatorially open, only two were without self-intersections).
Although very common, inputs with multiple, possibly nested,components are overlooked in previous works. The implementation of [Campen & Kobbelt 2010a] assumes single component input, and [CGAL 2015] does not detect nested voids automatically.
Our algorithm correctly handles both cases (see Figure 26).
In addition to standard tests on common computer
graphics models, we also stress test our algorithm on challenging
examples. Our algorithm is robust for carrying out consecutive
boolean operations because the output solid mesh is trivially a valid
PWN input for the following operations (see Figure 27).

Stress tests

Besides conventional boolean operations, the spacepartition defined by mesh arrangement is useful for many important geometry processing applications. In Figure 29, the outer hullcomputation is a necessary preprocessing step for generating a volumetric discretization for structural analysis [Zhou et al. 2013]. The
outer hull is also useful for culling extra internal complexity (seeFigure 30).
Generality

Our variadic formulation also allows us to compute regions insideat least k of the input meshes without the combinatorial explosion
associated with binary boolean operations (see Figure 31). Figure 3 considers ten intersecting tetrahedra. Our variadic union of
all ten tets is roughly twice as fast as decomposing the union intoa cascading tree of binary union operations (and 6.5× faster than alinear chain of binary unions). Intuitively, this is because our intersection resolution is the most economical for this arrangement. In
contrast, repeated unions will require resolving intersections withprevious results, aggregating unnecessary complexity: though geometrically identical, our result has 384 triangles, compared to thecascading tree’s 1000. We also construct extraction of the region
inside at least five input tetrahedra. Decomposing this task into abinary tree of cascading operations leads to an exponential numberof operations in the number of tets (5 10− 1 = 1259 binary oper5
ations for ten tets). The aggregation of complexity is catastrophic,
leading to performance measured in weeks. Instead, extracting thisresult from our arrangement requires the same cost as extracting theunion: just a few seconds.
Lastly, the cell data structure used by our algorithm can be easily
extended for customized applications. For example, it is easy toeliminate small cells immersed inside a shape (see Figure 32).

8

Limitations and future work

A limitation of this method is the requirement that the input meshhave no open boundaries or non-manifold “flaps”. Of the 10,000
meshes in our Thingi10K dataset, 18% did not meet our precon-Figure 28: The groved, yellow frog is the result of subtracting thestripy, blue flog from an unknown (presumably solid) frog using
some boolean implementation (not ours). Our robust union recovers the original frog (blue and yellow).

Figure 30: The coin mesh has been modeled with many overlapping components (gold). Constructing the outer hull adds many
new vertices, but removes hidden interior geometry. After decimation, the few triangles are well spent on the visible surface (silver).
Naive decimation wastes precious triangles on the hidden interior(bronze).

Figure 29: Tetrahedralization of the input foot mesh fails due tooverlapping input components, but succeeds after self-union. The
volume mesh helps analyze the shapes structure.

ditions. While repairing invalid input meshes is beyond the scopeof this method, many previous works are available [Attene 2010;
Campen et al. 2012; Jacobson et al. 2013]. On the other hand, thehigh-level structure of our approach would clearly extend to mesheswith boundaries and unstructured non-manifoldness via the generalized winding number [Jacobson et al. 2013]. However, our lowlevel combinatorial and sorting based subroutines would need to be
replaced with robust or exact evaluations of the generalized winding number. While Barki et al. [2015] provide a partial solution forsimple cases (e.g., clipping a closed model with a plane), a robust
solution for arbitrary geometrically open models is elusive.
Our method is variadic, but does not optimize operations based on
the requested extraction and inputs. For example, consider conducting the 1000-way union of 999 overlapping spheres enclosed
and their conservative bounding box. Clearly resolving the intersections between the 999 spheres is overkill. It would be interesting
to explore conservative optimization.
In the hopes of fostering continued work and more exhaustive testing in geometry processing at large, we release our code (now in
LIBIGL [Jacobson et al. 2016]) and our Thingi10K dataset.

Acknowledgments
We thank G. Bernstein for sharing code and M. Attene for testing on
the Thingi10K dataset. We thank M. Campen, A. Fleming H. Maia,
J. Panetta, R. Sawhney, O. Stein, P. Thamjaroenporn, O. Winn, andE. Yao for early feedback and proofreading. Funded in part by NSF
grants CMMI-11-29917, IIS-14-09286, and IIS-17257.

References
ATTENE , M. 2010. A lightweight approach to repairing digitized
polygon meshes. The Visual Computer.
ATTENE , M. 2014. Direct repair of self-intersecting meshes.
Graphical Models.

Figure 31: Our method supports n-ary operations. For example,
extracting all regions inside at least k of the input spheres centeredat each corner of the unit cube.

Figure 32: The self-union of a wheel of cheese retains its internal
bubbles after slicing (intersecting) with a knife (blue wedge). Slicing the outer hull reveals no bubbles. Eliminating small volumes
cells in the self-union before extraction, produces a few bubbles.
BANERJEE , R. P., AND ROSSIGNAC , J. R. 1996. Topologicallyexact evaluation of polyhedra defined in CSG with loose primitives. Comput. Graph. Forum.
BARKI , H., G UENNEBAUD , G., AND F OUFOU , S. 2015. Exact,robust, and efficient regularized booleans on general 3d meshes.
Computers and Mathematics with Applications.
B ERNSTEIN , G., AND F USSELL , D. 2009. Fast, exact, linear
booleans. In Proc. SGP.
B ERNSTEIN , G., 2013. Cork boolean libraryhttps://github.com/gilbo/cork.

.

B IERI , H., AND N EF, W. 1988. Elementary set operations withd-dimensional polyhedra. In Proc. IWCGA.
C AMPEN , M., AND KOBBELT, L. 2010. Exact and robust(self-)intersections for polygonal meshes. Comput. Graph. Forum.
C AMPEN , M., AND KOBBELT, L. 2010. Polygonal Boundary
Evaluation of Minkowski Sums and Swept Volumes. Comput.
Graph. Forum.

C AMPEN , M., ATTENE , M., AND KOBBELT, L., 2012. A practicalguide to polygon mesh repairing. Eurographics Tutorial.

Z HAO , H., WANG , C. C., C HEN , Y., AND J IN , X. 2011. Paralleland efficient boolean on polygonal solids. The Visual Computer.

CARVE, 2014. C ARVE: An efficient and robust library for boolean
operations on polyhedra. http://carve-csg.com/.

Z HOU , Q., PANETTA , J., AND Z ORIN , D. 2013. Worst-case structural analysis. ACM Trans. Graph..

CGAL, 2015. C GAL, Computational Geometry Algorithms Library. http://www.cgal.org.

Appendix: Exact representation

D OUZE , M., F RANCO , J.-S., AND R AFFIN , B. 2015. QuickCSG:
Arbitrary and faster boolean combinations of n solids. Tech. Rep.
01121419, Inria Research Centre Grenoble, Rhone-Alpes.
E DELSBRUNNER , H., AND M ÜCKE , E. P. 1990. Simulation of
simplicity: A technique to cope with degenerate cases in geometric algorithms. ACM Trans. Graph..
F ORTUNE , S. 1997. Vertex-rounding a three-dimensional polyhedral subdivision. Discrete Comput. Geom.
G RANADOS , M., H ACHENBERGER , P., H ERT, S., K ETTNER , L.,
M EHLHORN , K., AND S EEL , M. 2003. Boolean operations
on 3d selective nef complexes: Data structure, algorithms, andimplementation. In Proc. ESA.
JACOBSON , A., K AVAN , L., AND S ORKINE -H ORNUNG , O. 2013.
Robust inside-outside segmentation using generalized winding
numbers. ACM Trans. Graph..
JACOBSON , A., PANOZZO , D., ET AL ., 2016. libigl: A simple
C++ geometry processing library. http://libigl.github.io/libigl/.
M EI , G., AND T IPPER , J. C. 2013. Simple and robust
boolean operations for triangulated surfaces. arXiv preprint
arXiv:1308.4434.
M ILENKOVIC , V. J., AND NACKMAN , L. R. 1990. Finding compact coordinate representations for polygons and polyhedra.
M USETH , K., B REEN , D. E., W HITAKER , R. T., AND BARR ,A. H. 2002. Level set surface editing operators. ACM Trans.
Graph..
NAYLOR , B., A MANATIDES , J., AND T HIBAULT, W. 1990. Merging bsp trees yields polyhedral set operations. In Proc. SIGGRAPH.
PAVIC , D., C AMPEN , M., AND KOBBELT, L. 2010. Hybrid
Booleans. Comput. Graph. Forum.
S ACHT, L., JACOBSON , A., PANOZZO , D., S CHÜLLER , C., ANDS ORKINE -H ORNUNG , O. 2013. Consistent volumetric discretizations inside self-intersecting surfaces. Proc. SGP.
S ACKS , E., AND M ILENKOVIC , V. 2014. Robust cascading ofoperations on polyhedra. Computer-Aided Design (Tech. Note).
S UGIHARA , K., AND I RI , M. 1992. Construction of the Voronoi
diagram for one million generators in single-precision arithmetic. Proc. of the IEEE.
T HIBAULT, W. C., AND NAYLOR , B. F. 1987. Set operations onpolyhedra using binary space partitioning trees. In Proc. SIGGRAPH.
VARADHAN , G., K RISHNAN , S., S RIRAM , T., AND M ANOCHA ,D. 2004. Topology preserving surface extraction using adaptive
subdivision. In Proc. SGP.
WANG , C. C. L. 2011. Approximate boolean operations on large
polyhedral solids with partial mesh reconstruction. IEEE TVCG.
X U , S., AND K EYSER , J. 2013. Fast and robust booleans on polyhedra. Computer-Aided Design.

Since a computer cannot represent all
points in R3 , we assume all verticesof our input and output meshes are
given in the “exact” rational coordinate
space Q3 . Exactness means that points,
line segments and convex polygons withendpoints in Q3 form a group closed under spatial intersection.
The space of the double precision floating-point coordinate space
F3 does not satisfy these criteria: e.g., the point of intersection
between line segments [(0, 0, 0), (2, 1, 0)] and [(1, 0, 0), (0, 1, 0)]is the non-floating-point position (2/3, 1/3, 0), see inset. The exact rational space Q3 contains the floating-point space as a subset:
F3 ⊂ Q3 . So, if given a mesh Ai with floating-point vertex positions, we can losslessly cast them to our exact space.