Fast Tetrahedral Meshing in the Wild

TL;DR

The paper introduces fTetWild, a fast and robust tetrahedral meshing algorithm that converts triangle soups into high-quality meshes using floating-point arithmetic, improving speed and validity guarantees over previous methods.

Contribution

fTetWild replaces rational triangle insertion with an incremental approach, enabling faster processing and guaranteed validity of tetrahedral meshes with floating-point coordinates.

Findings

Significantly faster than previous methods

Guarantees valid tetrahedral meshes with floating-point coordinates

Successfully inserts all input triangles in practical tests

Abstract

We propose a new tetrahedral meshing method, fTetWild, to convert triangle soups into high-quality tetrahedral meshes. Our method builds on the TetWild algorithm, replacing the rational triangle insertion with a new incremental approach to construct and optimize the output mesh, interleaving triangle insertion and mesh optimization. Our approach makes it possible to maintain a valid floating-point tetrahedral mesh at all algorithmic stages, eliminating the need for costly constructions with rational numbers used by TetWild, while maintaining full robustness and similar output quality. This allows us to improve on TetWild in two ways. First, our algorithm is significantly faster, with running time comparable to less robust Delaunay-based tetrahedralization algorithms. Second, our algorithm is guaranteed to produce a valid tetrahedral mesh with floating-point vertex coordinates, while TetWild produces a valid mesh with rational coordinates which is not guaranteed to be valid after floating-point conversion. As a trade-off, our algorithm no longer guarantees that all input triangles are present in the output mesh, but in practice, as confirmed by our tests on the Thingi10k dataset, the algorithm always succeeds in inserting all input triangles.

Tables

Table 1. Table 1. A subset of the tet-subdivision table, the complete table is provided in the additional material. The first row corresponds to the case of a tetrahedron without cut edges, the second and third to the case of one cut edge, e 0 subscript 𝑒 0 e_{0} and e 1 subscript 𝑒 1 e_{1} respectively, and the last row to two cut edges e 0 subscript 𝑒 0 e_{0} and e 1 subscript 𝑒 1 e_{1} . All tetrahedra shown in the table have same edge label.

	0	1	$\cdots$
$0_{\left(10\right)}={000000}_{\left(2\right)}$			$\cdots$
$1_{\left(10\right)}={000001}_{\left(2\right)}$			$\cdots$
$2_{\left(10\right)}={000010}_{\left(2\right)}$			$\cdots$
$3_{\left(10\right)}={000011}_{\left(2\right)}$			$\cdots$
$⋮$	$⋮$	$⋮$	$\ddots$

Table 2. Table 2. Comparison of code robustness and performance on the Thingi10k dataset.

Method	Success rate	Out of memory ($>$32GB)	Time exceeded ($>$3h)	Algorithm limitation	Average time(s)
CGAL	79.00%	0.00%	0.00%	21.00%	11.7
TetGen	49.50%	0.10%	1.70%	48.70%	32.3
TetWild	*99.89%	0.05%	0.11%	0.00%	360.0
Ours	**99.97%	0.02%	0.03%	0.00%	49.8

Table 3. Table 3. Edge-cut configurations of a cutting tetrahedron before and after snapping. Numbers corresponds to the configurations in Figure  9 .

Before	1 vertex snapped	2 vertices snapped	3 vertices snapped
(2)	(1)	(2)	(1)
(3)	(1)(2)	(1)(2)	(1)
(5)	(1)(3)	(1)(2)	(1)(2)
(7)	(3)	(1)(3)	(1)

Equations

$$\mathrm{AMIPS}=\frac{\mathrm{tr}({\mathbf{J}}^T\mathbf{J})}{\det (\mathbf{J}{)}^{2/3}},$$

$$v_1$$

$$v_2$$

$$v_3$$

$$v_4$$

$${\mathrm{AMIPS}}_{1234}$$

$${\mathrm{AMIPS}}_{2341}$$

$${\mathrm{AMIPS}}_{3412}$$

$${\mathrm{AMIPS}}_{4123}$$

$${\mathrm{AMIPS}}_{1234}^3$$

$${\mathrm{AMIPS}}_{2341}^3$$

$${\mathrm{AMIPS}}_{3412}^3$$

$${\mathrm{AMIPS}}_{4123}^3$$

$$\mathrm{AMIPS}=1.127562687503913\mathrm{e}11.$$

Full text

††copyright: acmcopyright

1. Introduction

Tetrahedral meshes are commonly used in graphics and engineering applications. Tetrahedral meshing algorithms usually take a 3D surface triangle mesh as input and output a volumetric tetrahedral mesh filling the volume bounded by the input mesh. Traditional tetrahedral meshing algorithms have strong assumptions on the input, requiring it to be a closed manifold, free of self-intersections and numerical unstably close elements, and so on. However, those assumptions often do not hold on imperfect 3D geometric data in the wild.

The recently proposed Tetrahedral Meshing in the Wild (TetWild) (Hu et al., 2018) algorithm makes it possible to reliably tetrahedralize triangle soups by combining exact rational computations with a geometric tolerance to automatically address self-intersections, gaps and other imperfections in the input. The algorithm imposes no formal assumptions on the input mesh and is extremely robust, opening the door to automatic processing and repair of large collections of 3D models.

However, TetWild has two downsides, one theoretical and one practical. The theoretical downside is that it does not guarantee the generation of a floating point tetrahedral mesh: the algorithm internally uses rational numbers, which are then converted to floating point in the process of mesh optimization. While quite unlikely, it is possible that the mesh optimization stage will be unable to round all coordinates of the output mesh to floating point. The practical downside is the long running time compared with Delaunay-based tetrahedralization algorithms.

We introduce fTetWild, a variant of the TetWild algorithm addressing both these limitations while keeping the important properties of TetWild: robustness to imperfect input and ability to batch process large collections of models without parameter tuning, while producing high-quality tetrahedral meshes. Differently from TetWild, which generates a polyhedral rational mesh inserting all triangles at once, we start from a floating point tetrahedral mesh, insert one input triangle at a time and re-tetrahedralize locally, rejecting the operations producing inverted or degenerate elements. We then improve the quality of the mesh iteratively, and attempt to insert the rejected triangles into a higher quality mesh, which is less likely to fail.

Our algorithm always guarantees to generate a valid tetrahedral mesh with floating point vertex positions, independently from the stopping criteria or quality of the mesh. It might fail to insert few input triangles leading to a “less accurate” boundary preservation, however we never observe this behavior in our experiments. The new algorithm can be implemented using floating point constructions, avoiding the overhead associated with rational numbers. The use of floating point numbers also simplifies parallelization, which we use during mesh optimization to further improve the running time on large models. Consequently, our new algorithm is significantly faster than TetWild, with running times comparable to Delaunay-based algorithms (Figure 1), while providing the stronger guarantee of always producing a valid floating point output at the same time.

These improvements make fTetWild more practical than TetWild not only for volumetric meshing problems, but also for mesh repair and approximate mesh arrangements. By combining fTetWild and some elements of (Zhou et al., 2016), we obtain an approximate mesh arrangement algorithm for input triangle soups guaranteed to produce a valid floating point output. In comparison, the original algorithm presented in (Zhou et al., 2016) may fail to produce a floating-point output due to impossibility of rounding after the rational-arithmetic arrangement computation.

We demonstrate the robustness and practical utility of our algorithm by computing tetrahedral meshes on the Thingi10k dataset (10 000 models) and computing approximate Booleans. We use the generated tetrahedral meshes to solve elasticity, fluid flow, and heat diffusion equations on complex geometric domains. The complete implementation of fTetWild is provided in the additional material, together with scripts to reproduce all results in the paper.

2. Related Work

We briefly review the literature on tetrahedral meshing (Section 2.1), with an emphasis on envelope-based techniques, and we refer to (Cheng et al., 2012; Shewchuk, 2012) for a more detailed overview of the topic. We also review mesh repair and mesh arrangement algorithms (Section 2.2), since our technique can be also used in these settings to enable processing of imperfect geometry.

2.1. Tetrahedral Meshing

Delaunay Meshing.

The most studied and most widely used algorithms to generate tetrahedral meshes are based on the Delaunay condition (Chew, 1993; Shewchuk, 1998; Ruppert, 1995; Sheehy, 2012; Remacle, 2017; Du and Wang, 2003; Alliez et al., 2005a; Tournois et al., 2009; Murphy et al., 2001; Cohen-Steiner et al., 2002; Chew, 1989; Si and Gartner, 2005; Shewchuk, 2002a; Si and Shewchuk, 2014; Cheng et al., 2008; Boissonnat and Oudot, 2005; Jamin et al., 2015; Dey and Levine, 2008; Chen and Xu, 2004; Shewchuk, 1996; Aurenhammer, 1991; Shewchuk, 1999; Cheng et al., 2012; Aurenhammer et al., 2013; Bishop, 2016; Busaryev et al., 2009; Boissonnat et al., 2002; Si, 2015; George et al., 2003; Weatherill and Hassan, 1994). These methods are efficient and are widely used in commercial software. They can be applied to either point clouds inputs, or to tessellate the interior of manifold non self-intersecting meshes with no degenerate faces, but are not designed to deal with imperfect input, and, as a consequence, these techniques fail on a significant fraction of data in the wild (Hu et al., 2018). We provide a direct comparison with TetGen (Si, 2015), the most commonly used code in this category, which builds on the techniques developed in (George et al., 2003; Weatherill and Hassan, 1994), and CGAL (Wein et al., 2018) in Section 4.

Grid Methods.

An alternative approach is the use of a background grid as a starting point (Yerry and Shephard, 1983; Baker et al., 1988; Bern et al., 1994; Molino et al., 2003; Bronson et al., 2013; Labelle and Shewchuk, 2007; Doran et al., 2013; Bridson and Doran, 2014). These algorithms fill the entire bounding box of the input with a regular lattice or with a hierarchical space partitioning, optionally intersect the background mesh with the input surface, and then discard the elements outside of the input. These methods are simpler and more robust than Delaunay methods, but still struggle with imperfect input geometry, and create high-quality elements only in the interior of the mesh, where the background mesh is preserved exactly. However, placing badly shaped triangles on the boundary is problematic for many applications. Our algorithm borrows the idea of a background mesh from these methods, but inserts the elements incrementally, interleaving mesh optimization stages to ensure that the final quality of the mesh is uniformly high.

Front-Advancing Methods.

Another family of methods starts from the boundary, and inserts one element at a time, growing the volumetric mesh (i.e. marching in space), until the entire volume is filled (George, 1971; Sadek, 1980; Peraire et al., 1987; Cuilliere et al., 2013; Alauzet and Marcum, 2014; Haimes, 2014). These methods create high quality elements close to the boundary, but introduce many corner cases in the interior regions where the fronts meet, lowering the quality of the elements and making a robust implementation challenging.

Envelope Meshing.

All methods discussed above assume a valid, manifold, non self-intersecting boundary input mesh, and are not designed to handle the imperfections which are common in real-world CAD and scanned data. This issue has been tackled for surface meshes in Mandad et al. (2015) by creating a surface approximation within a tolerance volume using a modified Delaunay refinement process. A similar idea has been exploited for volumetric meshing in TetWild (Hu et al., 2018), and its 2D counterpart TriWild (Hu et al., 2019). The main idea of this work is to combine exact computation, using a hybrid kernel similar to (Attene, 2017), and a surface envelope (Hu et al., 2017), which allows the resulting mesh to approximate the input instead of reproducing it exactly. Our method closely follows (Hu et al., 2018), but we design our algorithm to avoid the use of exact computation. We compare the two techniques in Section 4.

Mesh Improvement.

Many algorithms have been proposed to improve the quality of an existing tetrahedral mesh by displacing vertices or changing the local connectivity (Canann et al., 1996, 1993; A. Freitag and Ollivier-Gooch, 1998; Lipman, 2012; Chen and Xu, 2004; Alliez et al., 2005b; Faraj et al., 2016; Feng et al., 2018; Hu et al., 2018; Alexa, 2019). Our method relies on the algorithm proposed in Hu et al. (2018), which uses a set of local operations to optimize the conformal AMIPS energy (Fu et al., 2015; Rabinovich et al., 2017). We parallelized some of the steps of that algorithm (Section 3), which is easier in our case since we only have floating point coordinates.

2.2. Applications

Mesh Repair.

Since our algorithm can be used for mesh repair, we review the most recent works on this topic, and we refer to (Attene et al., 2013) for a complete overview.

MeshFix (Attene, 2010, 2014) detects problematic regions in triangle meshes, and uses a set of local operations to heal them. The tool is very effective, but due to its use of a greedy algorithm it might delete large parts on a mesh. The most recent mesh repair technique has been introduced in (Hu et al., 2018): the algorithm generates a tetrahedral mesh and discards the generated tetrahedra, only keeping the boundary surface. While simple and effective, this techniques is computationally expensive, and thus only usable in batch processing applications. Our algorithm can be used in the same way, but its higher efficiency makes it more practical. We also propose a simple modification to the surface mesh extraction procedure to guarantee a manifold output.

Booleans and Mesh Arrangements.

Many approaches to performing Boolean operations on meshes were proposed, with some methods emphasizing robustness, other methods aiming to produce exact results, and another set prioritizing performance. In most cases, non-trivial assumptions are made on the input meshes: most commonly, these are required to be closed; in other cases, no self-intersections are allowed, or most restrictively vertices may be assumed in general position.

CGAL, one of the most robust implementations of Boolean operations available (Granados et al., 2003), relies on exact arithmetic, and uses a very general structure of Nef polyhedra (Bieri and Nef, 1988) to represent shapes. This allows one to obtain exact Boolean results in degenerate cases (e.g., when the result is a line or a point). At the same time, the assumptions on the input are quite restrictive: the surfaces need to be closed and manifold (although the latter constraint could be eliminated).

Another approach to achieve robustness at the expense of accuracy, is to convert input meshes to level sets e.g. by sampling a signed distance function for each object (Museth et al., 2002) and perform all operations on the level set functions. The obvious disadvantage of these methods is that their accuracy is limited by the resolution of the grid; the original mesh geometry is lost, and it is non-trivial to maintain even explicitly tagged features. These downsides are partially addressed by adaptive (Varadhan et al., 2004) and hybrid (Pavic et al., 2010; Wang, 2011; Zhao et al., 2011), the latter preserving mesh geometry away from intersections. All these methods rely on well-defined signed distance function, i.e., assume that input meshes are closed, and may still significantly alter the input geometry near intersections. (Schmidt and Singh, 2010) does not use a signed distance function, but resembles these methods, in that it removes existing geometry near intersections and replaces it by new mesh connecting the two objects and approximating the result of the Boolean. Binary Space Partitioning (BSP) based methods, starting from (Thibault and Naylor, 1987; Naylor et al., 1990) are closest in their approach to ours. Using BSP trees preserves the input more accurately, and, along the way, creates a volume partition. However, it is prone to errors due to numerical instability of intersection calculations, and, due to global intersections of triangle planes, performs excessive refinement. (Bernstein and Fussell, 2009) addresses the issue of non-robustness by using exact predicates, and (Campen and Kobbelt, 2010) reduces refinement by creating localized BSP trees in an octree. Examples of highly efficient but non-robust software for computing Booleans are (Douze et al., 2015), (Barki et al., 2015), and (Bernstein, 2013). A general position assumption is often required explicitly or implicitly. In (Zhou et al., 2016) a robust way to compute mesh arrangements is introduced, with Boolean operations as an application. Robustness is achieved by using rational numbers for critical computations. To perform Booleans the mesh is required to be Positive Winding Number (PWN), which does not always hold in meshes in the wild (Zhou and Jacobson, 2016).

Sheng et al. (2018a, b) use a combination of plane-based and vertex-based representations of mesh faces to improve robustness of basic operations needed for Boolean operations performed in floats. Their method achieves very high efficiency, at the expense of somewhat lower robustness compared to the state of the art (Zhou et al., 2016; Granados et al., 2003). Their method assumes that the input meshes enclose solids and are free of self-intersections. (Magalhães et al., 2017) is an efficient technique using simulation-of-simplicity techniques to handle general intersections between objects, self-intersections or holes are not handled. (Paoluzzi et al., 2017) considers a general problem of arrangements of complexes in 2D and 3D, presenting a theoretical general merge algorithm, but do not consider the questions of robustness and handling imperfect inputs.

Compared to existing methods, the application of fTetWild to Boolean operations is more conservative, in terms of mesh geometry changes and refinement, compared to level set and BSP-based methods, while maintaining their level of robustness. At the same time, thanks to the geometric tolerance, fTetWild is capable of eliminating near-degenerate or overly refined triangles in the input model, which (Zhou et al., 2016) cannot do. We also make fewer assumptions on the inputs, allowing gaps, self-intersections, and degeneracies.

3. Method

fTetWild takes as input a 3D triangle soup (i.e., a set of arbitrarily connected, potentially intersecting triangles with vertices potentially duplicated) whose vertices are represented in floating-point coordinates, representing the surface of an object. The algorithm has two user-defined parameters: target edge length $\ell$, and envelope size $\epsilon$. The $\epsilon$-envelope represents the maximal deviation from the input surface the user is willing to accept. For instance, in additive manufacturing applications $\epsilon$ can be the machining precision. It outputs a volumetric tetrahedral mesh, with floating-point vertex coordinates, whose elements are (1) non-inverted (i.e., positive volume) and (2) with boundary faces conforming to the input soup within a user-defined $\epsilon$-envelope. fTetWild makes no assumptions on the input triangle soup and it is robust when handling imperfect input with self-intersections or small gaps. This robustness is achieved by allowing the faces of the tetrahedral mesh corresponding to the input surface to move inside an $\epsilon$-envelope (up to $\epsilon$ far from the input): self-intersections, degenerate and near-degenerate faces and gaps contained in the envelope are automatically removed when combined with proper mesh improvement operations (Figure 2).

Similarities and Differences to Existing Face Insertion Algorithms.

The main challenge tackled in many existing tetrahedral meshing algorithm is the preservation of the input faces, which can be exact or approximate. One of the best known algorithms exactly preserving the input faces is (George et al., 2003), which proposed to subdivide a background mesh by intersecting it with input faces. This procedure can, however, introduce inverted elements due to floating-point rounding, which then need to be untangled, a difficult task for which no robust algorithm currently exists. A robust solution is proposed in TetWild (Hu et al., 2018), that initially inserts the faces exactly using rational numbers to avoid numerical problems, but is then forced to allow them to move to round the rational coordinates back to floating point. Although robust and conservative, this solution relies on expensive rational constructions, and it is not guaranteed to succeed in the rounding phase.

Our method follows an approach similar to TetWild (see Appendix A for a brief description of the algorithm), enabling small and controlled deviations from the input surface, but sidesteps the need for constructing a rational mesh, always using floating-point coordinates, while inheriting the robustness of TetWild. Algorithmically, there are three major differences:

1. (1)

fTetWild preserves the input faces by inserting one input triangle at a time into an existing background tetrahedral mesh. To facilitate the insertion it relaxes the insertion with a snapping tolerance (relatively larger than floating point machine precision) which is only possible thanks to the $\epsilon$-envelope. 2. (2)

fTetWild always tetrahedralizes the region affected by the newly inserted face by looking up a pre-computed table and always maintains a valid inversion-free tetrahedral mesh (using exact predicates). 3. (3)

fTetWild represents the vertices using only floating point coordinates, reducing the running time and memory consumption.

We note that inserting a triangle might fail due to limitations of the floating-point representation. For instance, the inserted face can be arbitrarily close to one of the existing vertices and the insertion will introduce a tetrahedron with a volume numerically equal to zero. In this scenario, we rollback the problematic operation, mark the problematic face as un-inserted, iteratively perform mesh improvement operations on the whole mesh, and try to insert the face again when the mesh quality has increased. This procedure shows the only possible failure of fTetWild: the impossibility of adding some input faces. While this is indeed possible, it never manifested in our experiments. Note that even if some faces could not be inserted, fTetWild still outputs a valid mesh conforming to all other faces.

3.1. Algorithm Overview

Our algorithm consists of four phases (Figure 3): (1) the input triangle soup is simplified while ensuring it stays in the $\epsilon$-envelope of the input (Section 3.3), (2) a background mesh is generated and the triangles are iteratively inserted into it, if the insertion does not introduce inverted elements (Section 3.4), (3) the mesh is improved using local operations (Section 3.5) and at the end of every three improvement iteration we re-attempt the insertion of input triangles that could not be inserted at phase 2, (4) the mesh elements are optionally filtered to remove the elements outside the surface or to perform Boolean operations (Section 3.6).

During the whole procedure we ensure that the tetrahedral mesh remains valid, that is, we ensure that (1) each element has positive volume (checked using exact predicates (Shewchuk, 1997; Lévy, 2019)) and (2) all successfully inserted triangles, from now on called the tracked surface, stay inside the $\epsilon$-envelope of the input.

Throughout the algorithm, we consider a distance between two points zero if it is below a numerical tolerance $\epsilon_{\mathrm{zero}}$. Similarly, we use $\epsilon_{\mathrm{zero}}^{2}$, $\epsilon_{\mathrm{zero}}^{3}$ for areas and volumes respectively. We found that the performance of the algorithm are not heavily affected by this tolerance, as long as $\epsilon_{\mathrm{zero}}>10^{-20}$: in our experiments we used $\epsilon_{\mathrm{zero}}=10^{-8}$.

3.2. Envelope

We use the envelope definition and the algorithms introduced in (Hu et al., 2018) to build the envelope and check if a triangle is contained in it. In particular, testing if a triangle is contained within the envelope is done by sampling the input triangle and checking if the samples are all within a slightly smaller envelope with the sampling error conservatively compensated (Hu et al., 2018).

3.3. Preprocessing

We use the same preprocessing procedure proposed in (Hu et al., 2018) for simplifying the input: we merge vertices closer than $\epsilon_{\mathrm{zero}}$ and collapse an edge if: (1) it is a manifold edge (has no more than two incident triangles) and vertex-adjacent edges are also manifold, and (2) the collapsing operation does not move triangles outside a smaller envelope of size $\epsilon_{\mathrm{prep}}<{\epsilon}$. At this stage, we use $\epsilon_{\mathrm{prep}}=0.8{\epsilon}$ since this value gives space for snapping in triangle insertion (Section 3.4), and prevents vertices to be too close to the boundary of the envelope, thus leaving more freedom for surface vertices to move in the mesh improvement stage (Section 3.5). On our dataset, we observed that changing this parameter has a minor impact on the running time and negligible effect on the output when in the range 0.7 to 0.999. Note that it cannot be set to 1 because it will prevent snapping (Section 3.4). We use the value 0.8 since it is far from the bounds of this range.

Since the preprocessing step is computationally expensive, due to the envelope containment checks, we propose a basic parallelization strategy which leads, on average, to a 4x speedup of the preprocessing step when using 8 cores. Our parallel edge collapsing procedure uses a serial 2-coloring pass over all input faces. We mark all input triangles white in the initial stage. Then, iteratively, we mark all edges as parallel-independent if all its vertex-adjacent triangles are white, and then mark these triangles black (Figure 4). At this point, we can safely collapse all marked parallel-independent edges in parallel. We iterate this procedure until we are unable to remove more than $0.01\%$ of the original input vertices.

3.4. Incremental Triangle Insertion

3.4.1. Background Mesh Generation

The triangle insertion stage requires a background mesh (which does not necessarily conform to the input triangles) which we create using Delaunay tetrahedralization (Lévy, 2019) on the vertices from the preprocessing stage. Since we allow the surface to move within an $\epsilon$-envelope, we generate the background mesh for a bounding box $2\epsilon$ larger than the bounding box of the input vertices. Similarly to TetWild, additional points are added uniformly inside the box and at least $\epsilon$ away from the input faces before Delaunay tetrahedralization to obtain more uniformly-shaped initial elements. More precisely, the additional points are added in a regular grid with spacing of $d/20$ (where $d$ is the diagonal of bounding box of the input mesh), skipping inserting the additional points with distance to the input faces smaller than $\epsilon$.

3.4.2. Single Triangle Insertion

The key component of our algorithm is three-stage procedure for inserting one triangle $T$ into a valid tetrahedral mesh $M$, adding new vertices and tetrahedra, and adjusting mesh connectivity, to minimize the number of insertion failures and number of badly shaped tetrahedra created by insertion. Our algorithm uses ideas from marching tetrahedra (Doi and Koide, 1991) and other tetrahedralization methods (George et al., 2003; Weatherill and Hassan, 1994), as well as marching cubes (Lorensen and Cline, 1987). It consists of the following steps:

(1) Find the set $\mathcal{T}_{I}$ consisting of the tetrahedra of $M$ that triangle $T$ cuts, as defined below; (2) Compute the intersection points of the plane spanned by $T$ (denoted as $P$) and the edges of the tetrahedra in $\mathcal{T}_{I}$; (3) Subdivide all cut tetrahedra using a connectivity pattern from a pre-computed tet-subdivision table. These patterns guarantee that a valid tetrahedral mesh connectivity is maintained.

Finding Cut Tetrahedra.

We first define that object $A$ cuts though object $B$ if their intersection contains interior points of both $A$ and $B$. We say that triangle $T$ cuts tetrahedron $\mathcal{T}$ if (1) it is completely contained inside $\mathcal{T}$, or (2) it cuts through at least one face of $\mathcal{T}$ (Figure 6). We initialize $\mathcal{T}_{I}$ to be the set of the tetrahedra of $M$ that $T$ cuts. Note that this set will be iteratively expanded by the algorithm.

We use exact predicates (Shewchuk, 1997; Lévy, 2019) for checking if a triangle is contained inside a tetrahedron. To detect if one triangle cuts through another, we combine the exact predicates with the algorithm (Guigue and Devillers, 2003). The use of predicates ensures topological correctness when using floating-point coordinates.

Plane-Tetrahedra Intersection.

To insert a triangle $T$ defining a plane $P$ into $\mathcal{T}_{I}$ (Figure 5(c)(d)), we need to ensure that after the insertion: (1) for every point $p\in T$ there is a face $F$ of the refined tets in the set $\mathcal{T}_{I}$ such that $\min_{p\in F}\|p-q\|<\delta$; (2) the projection of faces $F$ in $\mathcal{T}_{I}$, that are within the distance $\delta$ from $T$ to the plane $P$ covers $T$. We call sets with these properties covers of $T$. We allow the cover of triangles to deviate from $P$. This is crucial for robustly inserting triangles using floating point computations: without it, we observe a significantly higher running time due to more insertion failures, which leads to additional iterations of mesh optimization. Also, more faces remain uninserted in the final output. For the first pass of triangle insertion (i.e., before any mesh improvement is done), we use a larger $\delta=\max(\epsilon_{\mathrm{zero}},10^{-3}{\epsilon})$, while for all subsequent passes we reduce it to $\delta=\epsilon_{\mathrm{zero}}$.

We start with the idealized case of infinite-precision arithmetics. In this case, we can easily realize the covering of $T$ by intersecting all the faces of the tetrahedra in $\mathcal{T}_{I}$ with plane $P$. This generates a planar polygonal mesh $\mathcal{F}$ on $P$ covering $T$ and the maximal distance from $T$ to $\mathcal{F}$ is zero (Figure 7 (a)). The vertices of $\mathcal{F}$ are intersection points of $P$ and edges of $\mathcal{T}_{I}$.

When the floating point representation is used to represent the coordinates of vertices, round-off error may result in degenerate or inverted tetrahedra. Our approach is to reject insertion in these cases. However, to minimize the number of triangles that have to be rejected, we either snap vertices of the tetrahedral mesh $M$ to $P$ (Figure 7 (b)) or snap intersection points to vertices of $M$ (Figure 7 (c)(d)). Moving a vertex $v$ of $\mathcal{T}_{I}$ (thus changing $M$) changes the cover $\mathcal{F}$, because $P$ intersects the edges of $\mathcal{T}_{I}$ in different locations. As no intersection of $P$ with an edge of $\mathcal{T}_{I}$ disappears (at most, it may move to an endpoint), and no new intersections appear (other than at the endpoints shared with already intersected edges), the connectivity of $\mathcal{T}_{I}$ can be viewed as unchanged, possibly with some zero-length edges. We can view snapping vertices of $\mathcal{T}_{I}$ to $P$ as a deformation of $\mathcal{F}$, keeping it on plane $P$. If the affected vertices are in the interior of $\mathcal{F}$, $\mathcal{F}$ still covers $T$ since the boundary of $\mathcal{F}$ does not change. However, if moving $v$ changes the boundary of $\mathcal{F}$, the covering might be invalidated (Figure 7 (b)). In this case, before moving a boundary vertex, we extend $\mathcal{T}_{I}$ by adding its 1-ring neighbourhood, intersect it with $P$, and extend $\mathcal{F}$ accordingly. We repeat this process until all affected vertices are in the interior.

Moving the point $v$ to the plane $P$ might not always be possible, since it could invert some tetrahedra in $M$. In these cases, instead of moving $v$ to $P$, we deform $\mathcal{F}$ by moving some vertices of $\mathcal{F}$ to $v$, which is at $\delta$ distance from $P$ by definition (Figure 7 (c)(d)). Similarly to the previous case, this operation can only be applied on interior vertices of $\mathcal{F}$. We thus extend $\mathcal{T}_{I}$ if this operation affects vertices on the boundary of $\mathcal{F}$.

In practice, we never explicitly compute $\mathcal{F}$ on the plane $P$ since it is uniquely defined by the intersection points (Appendix E), but instead use the following 4 steps, that directly determine the vertices of $\mathcal{F}$ (the faces of $\mathcal{F}$ are obtained by table-based refinement of $\mathcal{T}_{I}$).

1. (1)

Find all vertices of the tetrahedra in $\mathcal{T}_{I}$, with distance to $P$ smaller than $\delta$ and put them in $\mathcal{V}_{\delta}$ (e.g., vertex $v$ in Figure 8(1)). 2. (2)

Move vertices in $\mathcal{V}_{\delta}$ to their closest points on $P$ if it does not invert any elements of $M$ (Figure 8(2)). 3. (3)

For each vertex in $\mathcal{V}_{\delta}$, add all of its vertex-adjacent tetrahedra to $\mathcal{T}_{I}$ if these are cut by $P$ and have faces covering $T$ (i.e., the projection of the face to $P$ intersects with $T$). For example, $[v,v_{0},v_{1}]$ in Figure 8(3) is added. 4. (4)

Repeat steps (1) to (3) until no more new tetrahedra are added to $\mathcal{T}_{I}$.

Table-based Tetrahedron Subdivision.

All tetrahedra sharing the edges of tetrahedra in $\mathcal{T}_{I}$ cut by $P$ are subdivided into sub-tetrahedra according to the tet-subdivision table. Note that this set of tetrahedra usually contains some tetrahedra from $\mathcal{T}_{I}$ (red elements in Figure 5(e)) and some neighboring tetrahedra of $\mathcal{T}_{I}$ (yellow elements in Figure 5(e)).

Since an edge can have at most one intersection point with $P$, the decomposition of the subdivided tetrahedron is largely (but not entirely) determined by which edges are cut. (An edge will be cut if it intersects with $P$ and neither endpoints are snapped.) We record all possible decompositions of a tetrahedron in a subdivision table, indexed by primary cut indices and secondary cut indices. The primary index (I) (Table 1), is a binary string, indicating which edges are cut. If two edges on a face are cut (3 edges of a face cannot be cut at the same time), there are two possible triangulations of this face and also multiple decompositions of the tetrahedron; the secondary index (II) is the number of a specific decomposition (Table 1). A primary index paired with a secondary index retrieves a unique decomposition of a tetrahedron.

For an oriented tetrahedron $\mathcal{T}$, there are $2^{6}=64$ combinations of possible intersection points on its edges, but not all 64 combinations can be practically realized. A direct enumeration shows that the following edge-cut configurations are impossible: (1) $\mathcal{T}$ has six cut edges, (2) $\mathcal{T}$ has five cut edges, (3) $\mathcal{T}$ has 4 cut edges, and 3 of them are on the same face, and (4) $\mathcal{T}$ has 3 cut edges on the same face. In total, there are ${6\choose 6}+{6\choose 5}+3\cdot 4+4=23$ impossible edge-cut configurations.

The remaining 41 realizable edge-cut configurations cover all subdivision cases and we can categorize them into 7 symmetry classes (Figure 9). Five of them were discussed in (Schweiger and Arridge, 2016) and used for a tetrahedron cut by a plane. We need 2 extra configurations (Figure 9, (4) and (6)) for subdividing the neighboring tetrahedra of $\mathcal{T}_{I}$ with only certain edges cut by $P$ (yellow elements in Figure 5(e)).

We retrieve a list of decompositions of $\mathcal{T}$ corresponding to a primary index; we now need to select a secondary index corresponding to a decomposition that preserves the validity of mesh $M$ after the subdivision, that is, we want $M$ to have a valid topology and no inverted tetrahedra. To ensure a valid topology, two adjacent

subdivided tetrahedra must have the same triangulation on the shared face. We set a rule for choosing such triangulation using only the global ordering of the vertices of $M$. For a face $[v_{0},v_{1},v_{2}]$ of tetrahedra $\mathcal{T}$ with two intersection points $p_{1},p_{2}$ on it (see inset), we select the triangulation containing the edge $[p_{2},v_{1}]$ if the unique integer label of vertex $v_{1}$ is larger than the one of $v_{2}$. Otherwise, we select the triangulation containing the edge $[p_{1},v_{2}]$. This simple rule completely identifies a secondary index and preserves the topology of the mesh. For completeness, we remark that some configurations might require additional vertices (Appendix C). However, our rule automatically excludes them. We attach the visualization of the tet-subdivision table in the supplementary material. We then check if all sub-tetrahedra have volume larger than $\epsilon_{\mathrm{zero}}^{3}$ (since we observed that elements with positive but extremely small volume could delay later insertions in this local region) and reject the insertion if this is not the case.

3.4.3. Open-boundary Edge Preservation.

After triangle insertion, the input edges shared by two non-coplanar triangles are preserved through the insertion of adjacent triangles, as the plane of the next inserted triangle will intersect the cover $\mathcal{F}$ of a previously inserted triangle. This does not hold for boundary edges. An edge is an open-boundary edge if it has only one incident triangle or has multiple coplanar incident triangles on the same side of the edge in their common plane.

To preserve an open-boundary edge $e$ of a triangle $T$, we project $e$ and the cover $\mathcal{F}$ of $T$ to the plane $P^{\prime}$ that best fits the faces of $\mathcal{F}$. Then, we compute the intersection of the projection of $e$ to $P^{\prime}$ with the projections of the faces of $\mathcal{F}$. The intersection points of the projection of $e$ and $\mathcal{T}$ are then computed on $P^{\prime}$ and are lifted back to the corresponding faces of $\mathcal{F}$. Since we have a set of edges cut into two, we can subdivide the affected tetrahedra using the previous table-based tetrahedron subdivision. An example can be found in Appendix D.

Note that the open-boundary edge preservation might fail due to numerical reasons, in this case we rollback the operation and postpone the insertion of the open-boundary triangle to later stages.

3.5. Mesh Improvement

We adapt the mesh improvement framework proposed in TetWild (Hu et al., 2018) that optimizes the conformal AMIPS 3D energy (Rabinovich et al., 2017) to increase the mesh quality, but avoid the overhead introduced by the hybrid kernel by specializing the framework for floating point computation. Note that, as mentioned in (Hu et al., 2018), we use the AMIPS energy since it is differentiable and scale-invariant. We made three changes to the original optimization:

1. (1)

We try to insert the uninserted input faces every three mesh improvement iterations until all input faces are inserted or the mesh improvement terminates. 2. (2)

We parallelize the vertex smoothing step using a simple graph coloring strategy. 3. (3)

We discovered an instability in the evaluation of the AMIPS energy computation in floating points, which sometimes leads to overrefinement in TetWild. We propose a fix using a hybrid evaluation that uses rational numbers to compute intermediate results.

(1) is a change required by the new algorithm, since not all faces can be inserted when computations are done in floating point. (2) is a minor, yet effective, modification that slightly improves performance (Figure 14). (3) is a subtle problem, which we now explain in more detail. The conformal AMIPS 3D energy is a Jacobian-based energy defined as:

[TABLE]

where $\mathbf{J}$ is the Jacobian of the transformation from a regular tetrahedron to the tetrahedron $\mathcal{T}$. The larger the energy is, the worse the quality of $\mathcal{T}$ is. The minimal value is 3, the energy of a regular tetrahedron. The AMIPS energy is invariant under permutation of the vertices of $\mathcal{T}$, however its numerical evaluation in floating-point arithmetic is not, due to floating-point rounding. Usually the difference is negligible, but when the energy of $\mathcal{T}$ is large (on the order of $10^{8}$), the floating point computation becomes unstable and the resulting energy could differ by two orders of magnitude, which means that the descent direction that appears to be decreasing the energy may be determined incorrectly (see Appendix B for a concrete example). This numerical instability might prevent mesh improvement and thus lead to over-refinement, since the algorithm is trying to add degrees of freedom unnecessary to improve the quality (Figure 10).

To address this issue, we raise the energy to the third degree making it completely rational, and evaluate it using rational computation for elements with energy larger than $10^{8}$. We round the computed value of its third degree to the 64-bit floating point representation, and then compute the cubic root. The rational computation is more accurate (and permutation invariant) but significantly slower. However, the cases of precision loss in the energy are rare, and the overall computational overhead is negligible. This change has a major effect on the speed and effectiveness of our mesh optimization (Figure 10), avoiding unnecessary refinement and decreasing the overall runtime.

The mesh improvement terminates when either a user-specified mesh quality or a user-controlled number of iterations is reached. To ensure a fair comparison, for the large dataset testing and all examples in the paper, we use the same stopping criteria (max AMIPS energy is smaller than 10 or the number of optimization iterations reaches 80) and input parameters (envelope size $\epsilon=10^{-3}d$, targeted edge length $\ell=d/20$, where $d$ is the diagonal’s length of the bounding box of the input mesh) as in (Hu et al., 2018).

We use the same strategy as in TetWild for handling inputs with open boundaries. We track the vertices on the open boundary and project them back to the open boundary during mesh improvement (Section 3.5). Elements are classified as inside or outside the surface using the generalized winding number (Section 3.6). An example of input with open boundary is shown in Figure 11.

3.6. Filtering

The output of the mesh improvement step is a volumetric tetrahedral mesh of the expanded bounding box of the input triangle soup, with the preprocessed input triangles inserted. We provide two ways of optionally filter the result, targeting two different applications. The first strategy uses the fast winding number (Barill et al., 2018) to filter the tetrahedra outside of the preserved/tracked input (Hu et al., 2018). The second application is a volumetric extension of the mesh arrangement algorithm (Zhou et al., 2016). In this case, the input becomes a set of triangle soups, coupled with a set of Boolean operations to perform on them. During the triangle insertion stage, we keep track of the provenance of each triangle, and use it at the end to compute a set of generalized winding numbers (one for each tracked input surface) for the centroids of all tetrahedra in the volumetric mesh. We use the set of winding numbers to decide which tetrahedron to keep by checking if it is supposed to be contained in the result of the Boolean operation. For instance, when intersecting two triangle soups, we keep all tetrahedra that are inside both input triangle soups, according to the winding number definition.

There are three major advantages of this method over (Zhou et al., 2016): (1) Boolean operations can be performed on non-PWN surfaces, (2) the output is equipped with a tetrahedral mesh, which could be useful in downstream applications, and (3) the surface quality is high since the algorithm is allowed to remesh within the $\epsilon$ envelope.

4. Results

Our algorithm is implemented in C++ and uses Eigen (Guennebaud et al., 2010) for the linear algebra routines. We perform a large-scale comparison of our method with other meshing methods on the Thingi10k dataset (Zhou and Jacobson, 2016), which contains 10 000 real-world surface triangle meshes. We run our experiments on cluster nodes with a Xeon E5-2690v4 2.6GHz, allowing every model to use up to 8 threads, 128GB memory, and 24 hours running time. The reference implementation and the scripts used to generate the results are attached to the submission and will be released as an open-source project.

4.1. Success Rate

With the above memory and time constraints, fTetWild successfully tetrahedralizes 100% of the 10 000 input meshes (Figure 12). Most of the input models can be tetrahedralized with less than 1GB of RAM as detailed in Figure 13. Note that very complex models might require more memory, for instance the one in Figure 21 uses around 17GB of memory.

As observed in (Hu et al., 2018), most of the state-of-the-art tet-meshers have low success rate on in-the-wild data. We summarize the results on the whole Thingi10k dataset in Table 2. Note that only our method and TetWild have high success rates: our average timing is however seven times faster than TetWild.

4.2. Running Time

Thingi10k Dataset (10000 Models)

We compare the running time of our method with TetWild. For a fair comparison, we disable our code optimizations that could be easily ported to TetWild, such as parallelization of the preprocessing and smoothing step, and using the recent fast winding number algorithm for the final filtering. Without these optimizations, our algorithm is 4 times faster than TetWild on average (80.4s vs 360s). With code optimizations, we further improve our running time to 49.8s on average on a machine with 8 cores, which is 7 times faster than the serial implementation of TetWild (Figure 14). On more complex examples, like the model in Figure 15, our method is up to 17 times faster than TetWild.

Reduced Thingi10k Dataset (4540 Models)

We use a reduced dataset containing the intersection of the Thingi10k models that TetGen, CGAL, TetWild, and our method all succeed on. The dataset contains 4540 models, and allows us to fairly compare the performance of the different methods. On average, our method is comparable (18.5s) to the widely used, Delaunay-based tetrahedral mesher TetGen (22s), and is faster than CGAL (95s) and TetWild (107s), while robustly handling imperfect inputs. Figure 1 shows the number of models requiring more than a given time. For example, within less than 2 minutes, our method successfully tetrahedralizes 98.7% of the inputs. It is interesting to note that the tail of the distribution of our method is shorter than both TetGen and CGAL. For instance, there are only 4 models where our method requires more than 16 minutes, differently from TetGen, CGAL, and TetWild which have 20, 122, and 25 models, respectively.

4.3. Mesh Quality

The geometric quality of meshes produced by our algorithm is similar to the meshes produced by TetWild (Figure 16), since our method implements a similar mesh optimization strategy. We quantitatively evaluate and compare the element quality of TetWild and our output using five different measures:

1. (1)

AMIPS energy (Equation (1)), range $[3,+\infty)$, optimal 3, 2. (2)

Minimal dihedral angle, range $(0,1.23]$, optimal 1.23, 3. (3)

Volume-to-edge ratio $6\sqrt{2}\,V/\ell_{\mathrm{max}}^{3}$, range $(0,1]$, optimal 1, 4. (4)

Aspect ratio $\sqrt{3/2}\,h_{\mathrm{min}}/\ell_{\mathrm{max}}$, range $(0,1]$, optimal 1, 5. (5)

Radius-to-edge ratio $2\sqrt{6}\,r_{\mathrm{in}}/\ell_{\mathrm{max}}$, range $(0,1]$, optimal 1,

where $V$ is the volume, $\ell_{\mathrm{max}}$ is the longest edge, $h_{\mathrm{min}}$ the minimum height, and $r_{\mathrm{in}}$ the radius of the inscribed circle of a tetrahedron $\mathcal{T}$. We use (3), (4) and (5) since these are standard measures for tetrahedral quality (Shewchuk, 2002b).

Figure 17 shows the histograms of worst and average element quality of 10 000 output meshes of TetWild and our method. The quality of our outputs are quite similar to TetWild’s output. We refer to the study in  (Hu et al., 2018, Figure 14) for the full quality comparison of TetWild and other tetrahedral meshing algorithms.

4.4. Mesh Density

Compared with TetWild, our method generates meshes of similar density (Figure 18). Both TetWild and our method aim to generate as-coarse-as possible meshes while preserving the input surface. This choice is useful in downstream applications to reduce their computational cost. Optionally, the algorithm supports a user-specified sizing field to increase the density if desired.

In contrast to our method, TetGen preserves the input surface geometry exactly and thus generates a dense tetrahedral mesh around the surface if the input surface mesh is dense, as visible in the model shown in Figure 1. CGAL approximates the surface by means of an implicit function, but sometimes over-refines sharp features and tiny artifacts as illustrated in Figure 1, where the dark spots are over-refined regions.

5. Applications

5.1. Mesh Repair

Similarly to TetWild, our algorithm can be used to repair imperfect triangle meshes by tetrahedralizing the volume and extracting the surface of the generated tetrahedral mesh. However, the mesh improvement step of our method (Section 3.5) can be stopped at any time since we maintain an inversion-free floating point tetrahedral mesh at all stages of our algorithm. High tetrahedral mesh quality is not required for this application, and we can stop mesh optimization as soon as all input faces are inserted, further reducing the running time. We compared our result with the state-of-the-art mesh repairing tool MeshFix (Attene, 2010) in Figure 19. Our method, while slower, provides a higher-quality result with controllable geometric error. A minor, yet important, observation is that keeping only the boundary of a valid tetrahedral mesh might generate a non-manifold surface mesh (Figure 20). To avoid this problem, we identify the non-manifold edges and split them. Then we duplicate every non-manifold vertex to ensure a global manifold output, using the algorithm proposed in (Attene et al., 2009). Note that this procedure ensures manifoldness, but introduces vertices in the same geometric position. With this minor change, our algorithm can be used to repair triangle meshes, guaranteeing the extraction of an high-quality, manifold boundary surface mesh within the prescribed distance from the input triangle soup.

We also tested an extremely challenging model coming from an industrial application in additive manufacturing (the part is copyrighted by Velo3D): the design of an exhaust pipe using a volume filled with a structure based on the gyroid triply periodic minimal surface. The model has a multitude of issues introduced during the modeling phase, but it can be cleaned up by our algorithm within 55 minutes (or 122 minutes with the envelope size decreased by a factor of two), compared to around two weeks of manual labor required by Velo3D’s current processing pipeline. Our output mesh (Figure 21) is directly usable for FEA, further editing, or fabrication. As a reference point, the original implementation of TetWild takes 215 minutes with a default envelope size. Another challenging model we tested contains complex thin structures coming from architecture (Figure 22). The method in (Masoud, 2016; Tabatabaei Ghomi et al., 2018) optimizes for the layout of a graph, then replaces the graph edges with cylinders of varying radii. To ensure solidity of the final structure, all cylinders are intersecting as shown in the close up. Although the mesh contains many irregularities, fTetWild successfully meshes the domain into an analysis-ready mesh.

5.2. Mesh Arrangements

Zhou et al. (2016) proposes to compute the arrangement between multiple surfaces using an algorithm to map Boolean operations into simple algebraic expressions involving the winding number of the input surfaces. Their method is robust, but only supports clean PWN surfaces as input. We propose a simple extension of this algorithm (as explained in Section 3.6) to arbitrary triangle soups. The advantages of our method is evident when the input surfaces come from CAD models containing small gaps or self-intersections: both Mesh Arrangements (Zhou et al., 2016) and CGAL (Hachenberger and Kettner, 2019) are unable to perform the operation (since it is not well-defined for non-PWN surface), while fTetWild can compute an approximate (since it allows for an $\epsilon$-deviation from the input surfaces) union, difference, and intersection between them (Figures 23, 24), providing robust (but slower) Boolean operations on imperfect geometries. The output is a tetrahedral mesh, which can be useful in downstream applications, and its boundary is a high quality surface triangular mesh.

5.3. Simulation

The main application of tetrahedral meshing is physical simulations, and the high-quality of our results makes them ideal to be directly used in downstream finite element software (Figure 25).

Additionally, the recently proposed a priori $p$-refinement (Schneider et al., 2018) is an ideal fit for our approach when targeting FEM applications, since fTetWild always produces a valid floating-point mesh. Schneider et al. (2018) provides a simple formula to determine the order of each element to compensate for its, possibly bad, shape. We can use this criterion to terminate the mesh optimization early in our algorithm (thus reducing the meshing time) without affecting the quality of the simulation, Figure 26.

We use the Boolean difference (Section 5.2) to generate the background mesh required for simulating the fluid flow on a cylindrical tube containing an obstacle (Figure 27).

6. Concluding Remarks

We introduced fTetWild, a novel robust tetrahedral meshing algorithm for triangle soups which combines the robustness of TetWild with a running time comparable to Delaunay-based methods. The improved performance makes this algorithm suitable not only for applications requiring a volumetric discretization, but also for surface mesh repair and Boolean operations.

Our current naive parallelization approach shows that our algorithm benefits from shared-memory parallelization; exploring more advanced parallelization techniques and extending it to distributed computation on HPC clusters are important directions for future work. Our iterative triangle insertion algorithm could be used in dynamic remeshing tasks, potentially allowing to reuse an existing mesh and insert new faces only in regions with high deformation. While conceptually trivial, extending our algorithm to 2D triangle meshing could improve the performance of (Hu et al., 2019).

Our algorithm uses the conformal AMIPS energy (Rabinovich et al., 2017) to measure and optimize the quality of the tetrahedra. An interesting alternative has been introduced concurrently to our work by (Alexa, 2019): they propose to optimize directly for the Dirichlet energy of the tetrahedralization and show that this measure is effective at removing slivers, while being computationally efficient to evaluate. A comparative study of the two measures would be interesting, and using the Dirichlet energy could lead to further reductions in the running time of our method.

Appendix A A Brief Description of the TetWild Algorithm

The TetWild algorithm (Hu et al., 2018) takes a 3D triangle soup as input and generates a tetrahedral mesh that (1) has no inverted or degenerate tetrahedra and (2) contains an approximation of the input surface within a user-defined $\epsilon$-envelope.

The method starts with an initial background mesh generated using an unconstrained Delaunay tetrahedralization on the input points plus an additional set of evenly-spaced points sampled from a regular grid. These additional points are added to improve the shape of the tetrahedra in the background mesh.

This step generates tetrahedra that might not represent the input faces of the triangle soup: to ensure that they are preserved TetWild uses a Binary Space Partitioning (BSP) subdivision step. Each input triangle is converted into a plane that cuts the tetrahedra of the background mesh. The output of this stage is a polyhedral mesh. To avoid numerical issues and to guarantee that the sub-elements in the polyhedral mesh are convex and non-inverted, TetWild converts the coordinates of the vertices of the initial background mesh into rational numbers and performs all computations using rational numbers.

Since any convex polyhedron can be trivially subdivided into tetrahedra by adding an additional point in its barycenter, a tetrahedral mesh that exactly preserves input triangles can be naturally obtained after BSP subdivision. However, the vertices of this tetrahedral mesh are represented in rational coordinated. Rounding them to floating point is not simple, since the BSP subdivision introduces badly-shaped tetrahedra which could flip after rounding. TetWild thus increases the quality of the elements using an hybrid optimization procedures that mixes floating point and rational representation.

During mesh improvement, the preserved input triangle’s faces are tracked and are allowed to move inside the $\epsilon$-envelope. The $\epsilon$-envelope limits the tracked surface from deviating from the input further than $\epsilon$.

TetWild uses four local operations for mesh improvement: (1) edge splitting, (2) edge collapsing, (3) edge swapping, and (4) vertex smoothing. Every operation is rolled back if the tracked surface leaves the envelope after the operation or if any tetrahedra are inverted, ensuring a valid output. Differently from other mesh improvement methods, TetWild uses the 3D conformal AMIPS energy for measuring the geometric quality of the tetrahedra. The AMIPS energy is scaling-invariant and easily differentiable, which boosts these traditional local operations.

As the quality of the mesh is improved, the rational coordinates can be gradually rounded into floating points. Theoretically, there might be some unroundable vertices, but it does not occur on the ten thousand models that TetWild has been tested on.

The final step is the removal of the tetrahedra outside of the tracked surface. To handle potentially noisy inputs, TetWild computes the winding number of the centroids of all tetrahedra with respect to the tracked surface, and filters out all tetrahedra with centroid’s winding-number larger than 0.5.

Appendix B Example of Unstable AMIPS Energy

If we compute the 3D AMIPS energy for a tetrahedron with these 4 vertices

[TABLE]

we obtain

[TABLE]

where the subscript indicates the vertex permutations. There are 24 premutations in total and here we pick 4 of them as an example. Even if we use the cube of the energy without rational we obtain fluctuations

[TABLE]

As reference the correct value computed with rational number is

[TABLE]

Appendix C Unused Decompositions of a Tetrahedron

We enumerated all the possible decompositions of a tetrahedron and discovered two symmetry classes (Figure 28) of triangulation of faces whose decomposition requires an additional internal vertex. Note that these two cases are never selected by our algorithm (we include them here for completeness), as our rule (Section 3.4.2) never selects these two cases.

We show that our rule does not select case 1 (Figure 28 left). By contradiction: since the configuration is selected then the edges $[p_{1},v_{2}]$, $[p_{2},v_{3}]$ and $[p_{3},v_{1}]$ are present, thus $v_{2}>v_{1}$, $v_{3}>v_{2}$, and $v_{1}>v_{3}$, according to our rule. Combining these inequalities, the indices of the vertices must satisfy $v_{3}>v_{2}>v_{1}>v_{3}$, which is impossible. Case 2 (Figure 28 right) is also not selected following a similar argument.

Appendix D An Example for Open-Boundary Edge Preservation

If triangle $T$ is the only inserted triangle and is entirely contained inside a tetrahedron $\mathcal{T}$ (Figure 29(1)), the intersection of the plane $P$ and $\mathcal{T}$ will be a larger polygon (marked in yellow) containing $T$. In this case, the edges of $T$, which are open-boundary edges, are not preserved. To preserve them, we subdivide the tetrahedra once more.

In Figure 29(1), $\mathcal{T}$ first get decomposed into sub-tetrahedra (Figure 29(2)). Then the faces covering $T$ are $\mathcal{F}=\{[p_{1},p_{2},p_{3}]$, $[p_{4},p_{2},p_{3}]\}$ Figure 29(3). The open-boundary edge $e$ and the faces in $\mathcal{F}$ are projected to the best-fitting plane of $p_{1},p_{2},p_{3}$, and $p_{4}$. The intersection points of the projection of $e$ and $\mathcal{T}$ are then computed in 2D and are lifted to 3D (3 green points in Figure 29(3)). Now there are 3 edges $[p_{1},p_{2}],[p_{2},p_{3}],[p_{3},p_{4}]$ cut into two. We thus subdivide all the neighbouring tetrahedra with the table-based subdivision.

Appendix E Changes of Edge-cut Configuration After Snapping

Table 3 shows all possible edge-cut configurations of a cutting tetrahedron $\mathcal{T}$ after snapping. The final configurations have no more than two vertices which makes the triangulation of $\mathcal{F}$ uniquely defined by the points. The table includes only the 4 symmetry classes where $\mathcal{T}$ is cut by plane $P$ and contains a face in $\mathcal{F}$ (Figure 9 (2)(3)(5)(7)), but excludes the remaining 3 classes where $\mathcal{T}$ is not cut or is just affected by their neighbors (Figure 9 (1)(4)(6)).

A tetrahedron $\mathcal{T}$ can have at most 3 vertices snapped. If $\mathcal{T}$ has all its 4 vertices within a $\delta$ distance to the $P$, we only snap the 3 vertices closer to $P$.