Fast Tetrahedral Meshing in the Wild
Yixin Hu, Teseo Schneider, Bolun Wang, Denis Zorin, Daniele Panozzo
FFast Tetrahedral Meshing in the Wild
YIXIN HU,
New York University, USA
TESEO SCHNEIDER,
New York University, USA
BOLUN WANG,
Beihang University, China
DENIS ZORIN,
New York University, USA
DANIELE PANOZZO,
New York University, USA >1m >2m >4m >8m >16m >32m0%10%20%30%40%50%
TetGen CGAL TetWild Ours
Ours
Input
Fig. 1. The bar charts show the percentage of models requiring more than the indicated time for the different approaches over 4 540 inputs (the subset ofThingi10k where all 4 compared algorithms succeed). Our algorithm successfully meshes 98.7% of the input models in less than 2 minutes, and processes allmodels within 32 minutes. The comparison has been done using the experimental setup of TetWild [Hu et al. 2018] and selecting a similar target resolution forall methods. The CGAL surface approximation parameter has been selected to be comparable to the envelope size used for TetWild and for our method. Theimages above the plot show a mouse skull model (from micro-CT) tetrahedralized with fTeWild (right) compared with other popular tetrahedral meshingalgorithms.
We propose a new tetrahedral meshing method, fTetWild, to convert trian-gle soups into high-quality tetrahedral meshes. Our method builds on theTetWild algorithm, replacing the rational triangle insertion with a new in-cremental approach to construct and optimize the output mesh, interleavingtriangle insertion and mesh optimization. Our approach makes it possible tomaintain a valid floating-point tetrahedral mesh at all algorithmic stages,eliminating the need for costly constructions with rational numbers usedby TetWild, while maintaining full robustness and similar output quality.This allows us to improve on TetWild in two ways. First, our algorithm issignificantly faster, with running time comparable to less robust Delaunay-based tetrahedralization algorithms. Second, our algorithm is guaranteedto produce a valid tetrahedral mesh with floating-point vertex coordinates,while TetWild produces a valid mesh with rational coordinates which is notguaranteed to be valid after floating-point conversion. As a trade-off, our
Authors’ addresses: Yixin Hu, New York University, USA, [email protected]; TeseoSchneider, New York University, USA, [email protected]; Bolun Wang, BeihangUniversity, China, [email protected]; Denis Zorin, New York University, USA,[email protected]; Daniele Panozzo, New York University, USA, [email protected] to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than ACMmust 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 afee. Request permissions from [email protected].© 2020 Association for Computing Machinery.XXXX-XXXX/2020/1-ART $15.00https://doi.org/10.1145/nnnnnnn.nnnnnnn algorithm no longer guarantees that all input triangles are present in theoutput mesh, but in practice, as confirmed by our tests on the Thingi10kdataset, the algorithm always succeeds in inserting all input triangles.
ACM Reference Format:
Yixin Hu, Teseo Schneider, Bolun Wang, Denis Zorin, and Daniele Panozzo.2020. Fast Tetrahedral Meshing in the Wild. 1, 1 (January 2020), 17 pages.https://doi.org/10.1145/nnnnnnn.nnnnnnn
Tetrahedral meshes are commonly used in graphics and engineeringapplications. Tetrahedral meshing algorithms usually take a 3Dsurface triangle mesh as input and output a volumetric tetrahedralmesh filling the volume bounded by the input mesh. Traditionaltetrahedral meshing algorithms have strong assumptions on theinput, requiring it to be a closed manifold, free of self-intersectionsand numerical unstably close elements, and so on. However, thoseassumptions often do not hold on imperfect 3D geometric data inthe wild.The recently proposed Tetrahedral Meshing in the Wild (TetWild)[Hu et al. 2018] algorithm makes it possible to reliably tetrahedralizetriangle soups by combining exact rational computations with ageometric tolerance to automatically address self-intersections, gapsand other imperfections in the input. The algorithm imposes no for-mal assumptions on the input mesh and is extremely robust, opening , Vol. 1, No. 1, Article . Publication date: January 2020. a r X i v : . [ c s . G R ] J a n • Hu, Y. et al the door to automatic processing and repair of large collections of3D models.However, TetWild has two downsides, one theoretical and onepractical. The theoretical downside is that it does not guarantee thegeneration of a floating point tetrahedral mesh: the algorithm inter-nally uses rational numbers, which are then converted to floatingpoint in the process of mesh optimization. While quite unlikely, itis possible that the mesh optimization stage will be unable to roundall coordinates of the output mesh to floating point. The practicaldownside is the long running time compared with Delaunay-basedtetrahedralization algorithms.We introduce fTetWild, a variant of the TetWild algorithm ad-dressing both these limitations while keeping the important proper-ties of TetWild: robustness to imperfect input and ability to batchprocess large collections of models without parameter tuning, whileproducing high-quality tetrahedral meshes. Differently from TetWild,which generates a polyhedral rational mesh inserting all trianglesat once, we start from a floating point tetrahedral mesh, insert oneinput triangle at a time and re-tetrahedralize locally, rejecting theoperations producing inverted or degenerate elements. We thenimprove the quality of the mesh iteratively, and attempt to insertthe rejected triangles into a higher quality mesh, which is less likelyto fail.Our algorithm always guarantees to generate a valid tetrahe-dral mesh with floating point vertex positions, independently fromthe stopping criteria or quality of the mesh. It might fail to insertfew input triangles leading to a “less accurate” boundary preserva-tion, however we never observe this behavior in our experiments.The new algorithm can be implemented using floating point con-structions, 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 therunning time on large models. Consequently, our new algorithm issignificantly faster than TetWild, with running times comparable toDelaunay-based algorithms (Figure 1), while providing the strongerguarantee of always producing a valid floating point output at thesame time.These improvements make fTetWild more practical than TetWildnot only for volumetric meshing problems, but also for mesh repairand approximate mesh arrangements. By combining fTetWild andsome elements of [Zhou et al. 2016], we obtain an approximatemesh arrangement algorithm for input triangle soups guaranteedto produce a valid floating point output. In comparison, the origi-nal algorithm presented in [Zhou et al. 2016] may fail to producea floating-point output due to impossibility of rounding after therational-arithmetic arrangement computation.We demonstrate the robustness and practical utility of our algo-rithm by computing tetrahedral meshes on the Thingi10k dataset(10 000 models) and computing approximate Booleans. We use thegenerated tetrahedral meshes to solve elasticity, fluid flow, and heatdiffusion equations on complex geometric domains. The completeimplementation of fTetWild is provided in the additional material,together with scripts to reproduce all results in the paper. 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 overviewof the topic. We also review mesh repair and mesh arrangementalgorithms (Section 2.2), since our technique can be also used inthese settings to enable processing of imperfect geometry.
Delaunay Meshing.
The most studied and most widely used algo-rithms to generate tetrahedral meshes are based on the Delaunaycondition [Alliez et al. 2005a; Aurenhammer 1991; Aurenhammeret al. 2013; Bishop 2016; Boissonnat et al. 2002; Boissonnat andOudot 2005; Busaryev et al. 2009; Chen and Xu 2004; Cheng et al.2008, 2012; Chew 1989, 1993; Cohen-Steiner et al. 2002; Dey andLevine 2008; Du and Wang 2003; George et al. 2003; Jamin et al.2015; Murphy et al. 2001; Remacle 2017; Ruppert 1995; Sheehy 2012;Shewchuk 1996, 1998, 1999, 2002a; Si 2015; Si and Gartner 2005; Siand Shewchuk 2014; Tournois et al. 2009; Weatherill and Hassan1994]. These methods are efficient and are widely used in commer-cial software. They can be applied to either point clouds inputs, or totessellate the interior of manifold non self-intersecting meshes withno degenerate faces, but are not designed to deal with imperfectinput, and, as a consequence, these techniques fail on a significantfraction of data in the wild [Hu et al. 2018]. We provide a directcomparison with TetGen [Si 2015], the most commonly used code inthis category, which builds on the techniques developed in [Georgeet 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 backgroundgrid as a starting point [Baker et al. 1988; Bern et al. 1994; Bridsonand Doran 2014; Bronson et al. 2013; Doran et al. 2013; Labelleand Shewchuk 2007; Molino et al. 2003; Yerry and Shephard 1983].These algorithms fill the entire bounding box of the input with aregular lattice or with a hierarchical space partitioning, optionallyintersect the background mesh with the input surface, and thendiscard the elements outside of the input. These methods are simplerand more robust than Delaunay methods, but still struggle withimperfect input geometry, and create high-quality elements only inthe interior of the mesh, where the background mesh is preservedexactly. However, placing badly shaped triangles on the boundary isproblematic for many applications. Our algorithm borrows the ideaof a background mesh from these methods, but inserts the elementsincrementally, interleaving mesh optimization stages to ensure thatthe final quality of the mesh is uniformly high.
Front-Advancing Methods.
Another family of methods starts fromthe boundary, and inserts one element at a time, growing the vol-umetric mesh (i.e. marching in space), until the entire volume isfilled [Alauzet and Marcum 2014; Cuilliere et al. 2013; George 1971;Haimes 2014; Peraire et al. 1987; Sadek 1980]. These methods createhigh quality elements close to the boundary, but introduce manycorner cases in the interior regions where the fronts meet, loweringthe quality of the elements and making a robust implementationchallenging. , Vol. 1, No. 1, Article . Publication date: January 2020. Envelope Meshing.
All methods discussed above assume a valid,manifold, non self-intersecting boundary input mesh, and are notdesigned to handle the imperfections which are common in real-world CAD and scanned data. This issue has been tackled for surfacemeshes in Mandad et al. [2015] by creating a surface approximationwithin a tolerance volume using a modified Delaunay refinementprocess. A similar idea has been exploited for volumetric meshingin TetWild [Hu et al. 2018], and its 2D counterpart TriWild [Huet al. 2019]. The main idea of this work is to combine exact compu-tation, using a hybrid kernel similar to [Attene 2017], and a surfaceenvelope [Hu et al. 2017], which allows the resulting mesh to ap-proximate the input instead of reproducing it exactly. Our methodclosely follows [Hu et al. 2018], but we design our algorithm toavoid the use of exact computation. We compare the two techniquesin Section 4.
Mesh Improvement.
Many algorithms have been proposed to im-prove the quality of an existing tetrahedral mesh by displacingvertices or changing the local connectivity [A. Freitag and Ollivier-Gooch 1998; Alexa 2019; Alliez et al. 2005b; Canann et al. 1996, 1993;Chen and Xu 2004; Faraj et al. 2016; Feng et al. 2018; Hu et al. 2018;Lipman 2012]. Our method relies on the algorithm proposed in Huet al. [2018], which uses a set of local operations to optimize theconformal AMIPS energy [Fu et al. 2015; Rabinovich et al. 2017]. Weparallelized some of the steps of that algorithm (Section 3), whichis easier in our case since we only have floating point coordinates.
Mesh Repair.
Since our algorithm can be used for mesh repair, wereview the most recent works on this topic, and we refer to [Atteneet al. 2013] for a complete overview.MeshFix [Attene 2010, 2014] detects problematic regions in tri-angle meshes, and uses a set of local operations to heal them. Thetool is very effective, but due to its use of a greedy algorithm itmight delete large parts on a mesh. The most recent mesh repairtechnique has been introduced in [Hu et al. 2018]: the algorithmgenerates a tetrahedral mesh and discards the generated tetrahedra,only keeping the boundary surface. While simple and effective, thistechniques is computationally expensive, and thus only usable inbatch processing applications. Our algorithm can be used in thesame way, but its higher efficiency makes it more practical. Wealso propose a simple modification to the surface mesh extractionprocedure to guarantee a manifold output.
Booleans and Mesh Arrangements.
Many approaches to perform-ing Boolean operations on meshes were proposed, with some meth-ods emphasizing robustness, other methods aiming to produce exactresults, 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-intersectionsare allowed, or most restrictively vertices may be assumed in generalposition.CGAL, one of the most robust implementations of Boolean oper-ations available [Granados et al. 2003], relies on exact arithmetic,and uses a very general structure of Nef polyhedra [Bieri and Nef1988] 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 latterconstraint could be eliminated).Another approach to achieve robustness at the expense of ac-curacy, is to convert input meshes to level sets e.g. by sampling asigned distance function for each object [Museth et al. 2002] andperform all operations on the level set functions. The obvious dis-advantage of these methods is that their accuracy is limited by theresolution of the grid; the original mesh geometry is lost, and it isnon-trivial to maintain even explicitly tagged features. These down-sides are partially addressed by adaptive [Varadhan et al. 2004] andhybrid [Pavic et al. 2010; Wang 2011; Zhao et al. 2011], the latterpreserving mesh geometry away from intersections. All these meth-ods rely on well-defined signed distance function, i.e., assume thatinput meshes are closed, and may still significantly alter the inputgeometry near intersections. [Schmidt and Singh 2010] does not usea signed distance function, but resembles these methods, in that itremoves existing geometry near intersections and replaces it by newmesh connecting the two objects and approximating the result of theBoolean. Binary Space Partitioning (BSP) based methods, startingfrom [Naylor et al. 1990; Thibault and Naylor 1987] are closest intheir approach to ours. Using BSP trees preserves the input moreaccurately, and, along the way, creates a volume partition. However,it is prone to errors due to numerical instability of intersection calcu-lations, and, due to global intersections of triangle planes, performsexcessive refinement. [Bernstein and Fussell 2009] addresses theissue of non-robustness by using exact predicates, and [Campen andKobbelt 2010] reduces refinement by creating localized BSP trees inan octree. Examples of highly efficient but non-robust software forcomputing Booleans are [Douze et al. 2015], [Barki et al. 2015], and[Bernstein 2013]. A general position assumption is often requiredexplicitly or implicitly. In [Zhou et al. 2016] a robust way to compute mesh arrangements is introduced, with Boolean operations as anapplication. Robustness is achieved by using rational numbers forcritical computations. To perform Booleans the mesh is required tobe Positive Winding Number (PWN), which does not always holdin 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 ba-sic operations needed for Boolean operations performed in floats.Their method achieves very high efficiency, at the expense of some-what lower robustness compared to the state of the art [Granadoset al. 2003; Zhou et al. 2016]. Their method assumes that the inputmeshes enclose solids and are free of self-intersections. [Magalhãeset al. 2017] is an efficient technique using simulation-of-simplicitytechniques to handle general intersections between objects, self-intersections or holes are not handled. [Paoluzzi et al. 2017] con-siders a general problem of arrangements of complexes in 2D and3D, presenting a theoretical general merge algorithm, but do notconsider the questions of robustness and handling imperfect inputs.Compared to existing methods, the application of fTetWild toBoolean operations is more conservative, in terms of mesh geome-try changes and refinement, compared to level set and BSP-basedmethods, while maintaining their level of robustness. At the sametime, thanks to the geometric tolerance, fTetWild is capable of , Vol. 1, No. 1, Article . Publication date: January 2020. • Hu, Y. et al
Input
Fig. 2. Example of an input surface mesh with self-intersections and a badtriangulation on the base. fTetWild converts this model into a high-qualitytetrahedral mesh. eliminating near-degenerate or overly refined triangles in the inputmodel, which [Zhou et al. 2016] cannot do. We also make fewerassumptions on the inputs, allowing gaps, self-intersections, anddegeneracies. fTetWild takes as input a 3D triangle soup (i.e., a set of arbitrarilyconnected, potentially intersecting triangles with vertices poten-tially duplicated) whose vertices are represented in floating-pointcoordinates, representing the surface of an object. The algorithm hastwo user-defined parameters: target edge length ℓ , and envelope size ϵ . The ϵ -envelope represents the maximal deviation from the inputsurface the user is willing to accept. For instance, in additive manu-facturing applications ϵ can be the machining precision. It outputsa volumetric tetrahedral mesh, with floating-point vertex coordi-nates, whose elements are (1) non-inverted (i.e., positive volume)and (2) with boundary faces conforming to the input soup within auser-defined ϵ -envelope. fTetWild makes no assumptions on theinput triangle soup and it is robust when handling imperfect inputwith self-intersections or small gaps. This robustness is achievedby allowing the faces of the tetrahedral mesh corresponding to theinput surface to move inside an ϵ -envelope (up to ϵ far from theinput): self-intersections, degenerate and near-degenerate faces andgaps contained in the envelope are automatically removed whencombined with proper mesh improvement operations (Figure 2). Similarities and Differences to Existing Face Insertion Algorithms.
The main challenge tackled in many existing tetrahedral meshing al-gorithm is the preservation of the input faces, which can be exact orapproximate. One of the best known algorithms exactly preservingthe input faces is [George et al. 2003], which proposed to subdividea background mesh by intersecting it with input faces. This proce-dure can, however, introduce inverted elements due to floating-pointrounding, which then need to be untangled, a difficult task for whichno robust algorithm currently exists. A robust solution is proposedin TetWild [Hu et al. 2018], that initially inserts the faces exactly using rational numbers to avoid numerical problems, but is thenforced to allow them to move to round the rational coordinates backto floating point. Although robust and conservative, this solutionrelies on expensive rational constructions, and it is not guaranteedto succeed in the rounding phase.Our method follows an approach similar to TetWild (see Appen-dix A for a brief description of the algorithm), enabling small andcontrolled deviations from the input surface, but sidesteps the needfor constructing a rational mesh, always using floating-point coor-dinates, while inheriting the robustness of TetWild. Algorithmically,there are three major differences:(1) fTetWild preserves the input faces by inserting one inputtriangle at a time into an existing background tetrahedralmesh. To facilitate the insertion it relaxes the insertion witha snapping tolerance (relatively larger than floating pointmachine precision) which is only possible thanks to the ϵ -envelope.(2) fTetWild always tetrahedralizes the region affected by thenewly inserted face by looking up a pre-computed table andalways maintains a valid inversion-free tetrahedral mesh(using exact predicates).(3) fTetWild represents the vertices using only floating pointcoordinates, reducing the running time and memory con-sumption.We note that inserting a triangle might fail due to limitations ofthe floating-point representation. For instance, the inserted face canbe arbitrarily close to one of the existing vertices and the insertionwill introduce a tetrahedron with a volume numerically equal tozero. In this scenario, we rollback the problematic operation, markthe problematic face as un-inserted, iteratively perform mesh im-provement operations on the whole mesh, and try to insert the faceagain when the mesh quality has increased. This procedure showsthe only possible failure of fTetWild: the impossibility of addingsome input faces. While this is indeed possible, it never manifestedin our experiments. Note that even if some faces could not be in-serted, fTetWild still outputs a valid mesh conforming to all otherfaces. Our algorithm consists of four phases (Figure 3): (1) the input tri-angle soup is simplified while ensuring it stays in the ϵ -envelopeof the input (Section 3.3), (2) a background mesh is generated andthe triangles are iteratively inserted into it, if the insertion does notintroduce inverted elements (Section 3.4), (3) the mesh is improvedusing local operations (Section 3.5) and at the end of every threeimprovement iteration we re-attempt the insertion of input trian-gles that could not be inserted at phase 2, (4) the mesh elements areoptionally filtered to remove the elements outside the surface or toperform Boolean operations (Section 3.6).During the whole procedure we ensure that the tetrahedral meshremains valid , that is, we ensure that (1) each element has posi-tive volume (checked using exact predicates [Lévy 2019; Shewchuk1997]) and (2) all successfully inserted triangles, from now on calledthe tracked surface , stay inside the ϵ -envelope of the input. , Vol. 1, No. 1, Article . Publication date: January 2020. Input Preprocessing Triangle Insertion Mesh Improvement Output
Fig. 3. Overview of our algorithm. From left to right, the input mesh is simplified, a background mesh is created and the input faces are inserted, the meshquality is optimized, and the final result is obtained by filtering the elements lying outside the input surface.Fig. 4. A two-dimensional example of edge coloring. From left to right: oneparallel-independent edge is selected (in red), its vertex-adjacent trianglesare colored in black. The algorithm proceeds until no edges can be marked(last image). In the end, all red edges can be safely collapsed in parallelwithout affecting other red edges.
Throughout the algorithm, we consider a distance between twopoints zero if it is below a numerical tolerance ϵ zero . Similarly, weuse ϵ , ϵ for areas and volumes respectively. We found thatthe performance of the algorithm are not heavily affected by thistolerance, as long as ϵ zero > − : in our experiments we used ϵ zero = − . We use the envelope definition and the algorithms introduced in[Hu et al. 2018] to build the envelope and check if a triangle iscontained in it. In particular, testing if a triangle is contained withinthe envelope is done by sampling the input triangle and checkingif the samples are all within a slightly smaller envelope with thesampling error conservatively compensated [Hu et al. 2018].
We use the same preprocessing procedure proposed in [Hu et al.2018] for simplifying the input: we merge vertices closer than ϵ zero and collapse an edge if: (1) it is a manifold edge (has no more thantwo incident triangles) and vertex-adjacent edges are also manifold,and (2) the collapsing operation does not move triangles outsidea smaller envelope of size ϵ prep < ϵ . At this stage, we use ϵ prep = . ϵ since this value gives space for snapping in triangle insertion(Section 3.4), and prevents vertices to be too close to the boundaryof the envelope, thus leaving more freedom for surface vertices tomove in the mesh improvement stage (Section 3.5). On our dataset,we observed that changing this parameter has a minor impact onthe running time and negligible effect on the output when in therange 0.7 to 0.999. Note that it cannot be set to 1 because it willprevent snapping (Section 3.4). We use the value 0.8 since it is farfrom the bounds of this range. Since the preprocessing step is computationally expensive, due tothe envelope containment checks, we propose a basic parallelizationstrategy which leads, on average, to a 4x speedup of the preprocess-ing step when using 8 cores. Our parallel edge collapsing procedureuses a serial 2-coloring pass over all input faces. We mark all inputtriangles white in the initial stage. Then, iteratively, we mark alledges as parallel-independent if all its vertex-adjacent triangles arewhite, and then mark these triangles black (Figure 4). At this point,we can safely collapse all marked parallel-independent edges inparallel. We iterate this procedure until we are unable to removemore than 0 .
01% of the original input vertices.
The triangle insertion stagerequires a background mesh (which does not necessarily conformto the input triangles) which we create using Delaunay tetrahedral-ization [Lévy 2019] on the vertices from the preprocessing stage.Since we allow the surface to move within an ϵ -envelope, we gen-erate the background mesh for a bounding box 2 ϵ larger than thebounding box of the input vertices. Similarly to TetWild, additionalpoints are added uniformly inside the box and at least ϵ away fromthe input faces before Delaunay tetrahedralization to obtain moreuniformly-shaped initial elements. More precisely, the additionalpoints are added in a regular grid with spacing of d /
20 (where d isthe diagonal of bounding box of the input mesh), skipping insertingthe additional points with distance to the input faces smaller than ϵ . The key component of our algo-rithm is three-stage procedure for inserting one triangle T into avalid tetrahedral mesh M , adding new vertices and tetrahedra, andadjusting mesh connectivity, to minimize the number of insertionfailures and number of badly shaped tetrahedra created by insertion.Our algorithm uses ideas from marching tetrahedra [Doi and Koide1991] and other tetrahedralization methods [George et al. 2003;Weatherill and Hassan 1994], as well as marching cubes [Lorensenand Cline 1987]. It consists of the following steps: (1) Find the set T I consisting of the tetrahedra of M that triangle T cuts, as defined be-low; (2) Compute the intersection points of the plane spanned by T (denoted as P ) and the edges of the tetrahedra in T I ; (3) Subdivide all , Vol. 1, No. 1, Article . Publication date: January 2020. • Hu, Y. et al M p q p p p p M p q v ( v is moved to line pq. ) p p p v ( v is not moved to line pq. ) p p p Traingle Insertion with SnappingTraingle Insertion without Snapping Case 1Case 2Case 3 . (a) Initial Mesh (b) Cutting Triangles (c) Snapping (d) Intersection (e) Triangles Subdivision (f) Final Mesh .
Fig. 5. Segment insertion into a triangle mesh (a 2D analog of triangle insertion) with and without snapping. (a) Insertion of segment pq into mesh M . (b)Identification of cut triangles T I (in red). (c) Snapping vertex v to line pq and updating T I , where v is δ -close to pq . v is moved to its closest points on line pq if this does not invert any elements of M (case 1). Vertex v is added in V δ both if v is moved (case 1) or not (case 2). The yellow triangle is added to T I . (d)Computing the intersection of line pq with the edges of T I (points p , p , p ). (e) Triangles requiring subdivision (shown in red). (f) The final mesh aftersubdivision. T T p p p p p p p p p T ∈ T I T (cid:60) T I Fig. 6. Examples of tetrahedra T included into, or excluded from, T I . Theintersections are shown in red. Left two: T intersects a face of T at a segment( [ p , p ] ) or a polygon ( [ p , p , p ] ) that contains interior points of both T and the intersected face of T . In this case, we put T into T I . Right two: T intersects a face of T (in light red) at a segment ( [ p , p ] ) that does notcontain any interior points of either T or the intersected face of T . In thiscase, we do not put T into T I . cut tetrahedra using a connectivity pattern from a pre-computed tet-subdivision table . These patterns guarantee that a valid tetrahedralmesh 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 T if (1) it is completelycontained inside T , or (2) it cuts through at least one face of T (Figure 6). We initialize 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 [Lévy 2019; Shewchuk 1997] for checkingif a triangle is contained inside a tetrahedron. To detect if one trian-gle cuts through another, we combine the exact predicates with thealgorithm [Guigue and Devillers 2003]. The use of predicates ensurestopological correctness when using floating-point coordinates. Plane-Tetrahedra Intersection.
To insert a triangle T defining aplane P into T I (Figure 5(c)(d)), we need to ensure that after theinsertion: (1) for every point p ∈ T there is a face F of the refined tetsin the set T I such that min p ∈ F ∥ p − q ∥ < δ ; (2) the projection of faces F in T I , that are within the distance δ from T to the plane P covers T . We call sets with these properties covers of T . We allow the coverof triangles to deviate from P . This is crucial for robustly insertingtriangles using floating point computations: without it, we observe p p p TP T v F (a) TP T v F (b) p p TP T v F (c) T p p P T v F (d) Fig. 7. Plane P intersects T I ( T ∈ T I ). (a) The faces of F (marked in yellow)are the covering of triangle T . (b) Snapping v to its closest point on P and expanding F to include red triangles makes F safely covering T . (c)Snapping a boundary vertex p of F to v changes the area of F and make itnot covering T . (d) Snapping an interior vertex p of F to v does not changethe boundary of F : F is still covering T . 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 passof triangle insertion (i.e., before any mesh improvement is done), weuse a larger δ = max ( ϵ zero , − ϵ ) , while for all subsequent passeswe reduce it to δ = ϵ zero .We start with the idealized case of infinite-precision arithmetics.In this case, we can easily realize the covering of T by intersectingall the faces of the tetrahedra in T I with plane P . This generates aplanar polygonal mesh F on P covering T and the maximal distancefrom T to F is zero (Figure 7 (a)). The vertices of F are intersectionpoints of P and edges of T I .When the floating point representation is used to represent thecoordinates of vertices, round-off error may result in degenerate , Vol. 1, No. 1, Article . Publication date: January 2020. p qv p qv p qvv v v ’ v v (1) (2) (3) Fig. 8. 2D illustration of step (1) to (3) of snapping when inserting segment [ p , q ] . The cut triangles in T I are marked in red. (1) Vertex v is within δ distance to segment [ p , q ] . (2) Check the effect of moving v to line [ p , q ] .(Vertex v cannot be moved to [ p , q ] in this case, because a triangle revertsits orientation.) (3) One-ring triangle [ v , v , v ] (marked in yellow) of v isintersecting with line [ p , q ] , and the projection of its edges [ v , v ] , [ v , v ] on line [ p , q ] , segment [ v ′ , v ′ ] , [ v ′ , v ′ ] , intersect with segment [ p , q ] . or inverted tetrahedra. Our approach is to reject insertion in thesecases. However, to minimize the number of triangles that have tobe 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 T I (thus changing M ) changes the cover F , because P intersects the edges of T I in different locations. As nointersection of P with an edge of T I disappears (at most, it may moveto an endpoint), and no new intersections appear (other than at theendpoints shared with already intersected edges), the connectivityof T I can be viewed as unchanged, possibly with some zero-lengthedges. We can view snapping vertices of T I to P as a deformation of F , keeping it on plane P . If the affected vertices are in the interiorof F , F still covers T since the boundary of F does not change.However, if moving v changes the boundary of F , the coveringmight be invalidated (Figure 7 (b)). In this case, before moving aboundary vertex, we extend T I by adding its 1-ring neighbourhood,intersect it with P , and extend F accordingly. We repeat this processuntil 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, insteadof moving v to P , we deform F by moving some vertices of F to v ,which is at δ distance from P by definition (Figure 7 (c)(d)). Similarlyto the previous case, this operation can only be applied on interiorvertices of F . We thus extend T I if this operation affects vertices onthe boundary of F .In practice, we never explicitly compute F on the plane P sinceit is uniquely defined by the intersection points (Appendix E), butinstead use the following 4 steps, that directly determine the verticesof F (the faces of F are obtained by table-based refinement of T I ).(1) Find all vertices of the tetrahedra in T I , with distance to P smaller than δ and put them in V δ (e.g., vertex v in Fig-ure 8(1)).(2) Move vertices in V δ to their closest points on P if it does notinvert any elements of M (Figure 8(2)).(3) For each vertex in V δ , add all of its vertex-adjacent tetrahedrato T I if these are cut by P and have faces covering T (i.e., theprojection of the face to P intersects with T ). For example, [ v , v , v ] in Figure 8(3) is added.(4) Repeat steps (1) to (3) until no more new tetrahedra are addedto T I . Table-based Tetrahedron Subdivision.
All tetrahedra sharing theedges of tetrahedra in T I cut by P are subdivided into sub-tetrahedra Table 1. A subset of the tet-subdivision table, the complete table is pro-vided in the additional material. The first row corresponds to the case of atetrahedron without cut edges, the second and third to the case of one cutedge, e and e respectively, and the last row to two cut edges e and e .All tetrahedra shown in the table have same edge label. I II · · · ( ) = ( ) e e e e e e · · · ( ) = ( ) · · · ( ) = ( ) · · · ( ) = ( ) · · · ... ... ... . . . according to the tet-subdivision table. Note that this set of tetra-hedra usually contains some tetrahedra from T I (red elements inFigure 5(e)) and some neighboring tetrahedra of T I (yellow elementsin Figure 5(e)).Since an edge can have at most one intersection point with P ,the decomposition of the subdivided tetrahedron is largely (butnot entirely) determined by which edges are cut. (An edge will becut if it intersects with P and neither endpoints are snapped.) Werecord all possible decompositions of a tetrahedron in a subdivisiontable, indexed by primary cut indices and secondary cut indices . Theprimary index (I) (Table 1), is a binary string, indicating which edgesare cut. If two edges on a face are cut (3 edges of a face cannot be cutat the same time), there are two possible triangulations of this faceand also multiple decompositions of the tetrahedron; the secondaryindex (II) is the number of a specific decomposition (Table 1). Aprimary index paired with a secondary index retrieves a uniquedecomposition of a tetrahedron.For an oriented tetrahedron T , there are 2 =
64 combinations ofpossible intersection points on its edges, but not all 64 combinationscan be practically realized. A direct enumeration shows that thefollowing edge-cut configurations are impossible: (1) T has six cutedges, (2) T has five cut edges, (3) T has 4 cut edges, and 3 of themare on the same face, and (4) T has 3 cut edges on the same face.In total, there are (cid:0) (cid:1) + (cid:0) (cid:1) + · + =
23 impossible edge-cutconfigurations.The remaining 41 realizable edge-cut configurations cover allsubdivision cases and we can categorize them into 7 symmetryclasses (Figure 9). Five of them were discussed in [Schweiger andArridge 2016] and used for a tetrahedron cut by a plane. We need2 extra configurations (Figure 9, (4) and (6)) for subdividing theneighboring tetrahedra of T I with only certain edges cut by P (yellowelements in Figure 5(e)). , Vol. 1, No. 1, Article . Publication date: January 2020. • Hu, Y. et al (1) No cut (2) 1 edge cut (3) 2 edge cut (4) 2 edge cut (5) 3 edge cut (6) 3 edge cut (7) 4 edge cut Fig. 9. 7 symmetry classes of edge-cut configurations. (4) and (6) can only happen on neighboring tetrahedra of T I with only certain edges cut by P . We retrieve a list of decompositions of T corresponding to aprimary index; we now need to select a secondary index corre-sponding to a decomposition that preserves the validity of mesh M after the subdivision, that is, we want M to have a valid topologyand no inverted tetrahedra. To ensure a valid topology, two adjacent v v v p p T subdivided tetrahedra must have the sametriangulation on the shared face. We set arule for choosing such triangulation usingonly the global ordering of the vertices of M .For a face [ v , v , v ] of tetrahedra T withtwo intersection points p , p on it (see in-set), we select the triangulation containing the edge [ p , v ] if theunique integer label of vertex v is larger than the one of v . Other-wise, we select the triangulation containing the edge [ p , v ] . Thissimple rule completely identifies a secondary index and preservesthe topology of the mesh. For completeness, we remark that someconfigurations might require additional vertices (Appendix C). How-ever, our rule automatically excludes them. We attach the visual-ization of the tet-subdivision table in the supplementary material.We then check if all sub-tetrahedra have volume larger than ϵ (since we observed that elements with positive but extremely smallvolume could delay later insertions in this local region) and rejectthe insertion if this is not the case. After triangle insertion,the input edges shared by two non-coplanar triangles are preservedthrough the insertion of adjacent triangles, as the plane of the nextinserted triangle will intersect the cover F of a previously insertedtriangle. This does not hold for boundary edges. An edge is an open-boundary edge if it has only one incident triangle or has multiplecoplanar incident triangles on the same side of the edge in theircommon plane.To preserve an open-boundary edge e of a triangle T , we project e and the cover F of T to the plane P ′ that best fits the faces of F . Then, we compute the intersection of the projection of e to P ′ with the projections of the faces of F . The intersection points ofthe projection of e and T are then computed on P ′ and are liftedback to the corresponding faces of F . Since we have a set of edgescut into two, we can subdivide the affected tetrahedra using theprevious table-based tetrahedron subdivision. An example can befound in Appendix D.Note that the open-boundary edge preservation might fail dueto numerical reasons, in this case we rollback the operation andpostpone the insertion of the open-boundary triangle to later stages. Input
Fig. 10. Example of model where the numerical instability of the AMIPSenergy causes over-refinement (middle). By evaluating the energy usingrational numbers (when it is above ) the issue disappears (right). We adapt the mesh improvement framework proposed in TetWild [Huet al. 2018] that optimizes the conformal AMIPS 3D energy [Ra-binovich et al. 2017] to increase the mesh quality, but avoid theoverhead introduced by the hybrid kernel by specializing the frame-work for floating point computation. Note that, as mentioned in [Huet al. 2018], we use the AMIPS energy since it is differentiable andscale-invariant. We made three changes to the original optimization:(1) We try to insert the uninserted input faces every three meshimprovement iterations until all input faces are inserted orthe mesh improvement terminates.(2) We parallelize the vertex smoothing step using a simple graphcoloring strategy.(3) We discovered an instability in the evaluation of the AMIPSenergy computation in floating points, which sometimes leadsto overrefinement in TetWild. We propose a fix using a hybridevaluation that uses rational numbers to compute intermedi-ate results.(1) is a change required by the new algorithm, since not all facescan be inserted when computations are done in floating point. (2) is aminor, yet effective, modification that slightly improves performance(Figure 14). (3) is a subtle problem, which we now explain in moredetail. The conformal AMIPS 3D energy is a Jacobian-based energydefined as: AMIPS = tr ( J T J ) det ( J ) / , (1)where J is the Jacobian of the transformation from a regular tetra-hedron to the tetrahedron T . The larger the energy is, the worsethe quality of T is. The minimal value is 3, the energy of a regulartetrahedron. The AMIPS energy is invariant under permutation of , Vol. 1, No. 1, Article . Publication date: January 2020. Input
Fig. 11. Input with open boundary on the bottom (left). The output tetrahe-dral mesh preserves the input geometry and closes the open side (middle).Users can choose to enable an additional smoothing process for smoothingthe open region (right). the vertices of T , however its numerical evaluation in floating-pointarithmetic is not, due to floating-point rounding. Usually the differ-ence is negligible, but when the energy of T is large (on the orderof 10 ), the floating point computation becomes unstable and theresulting energy could differ by two orders of magnitude, whichmeans that the descent direction that appears to be decreasingthe energy may be determined incorrectly (see Appendix B for aconcrete example). This numerical instability might prevent meshimprovement and thus lead to over-refinement, since the algorithmis trying to add degrees of freedom unnecessary to improve thequality (Figure 10).To address this issue, we raise the energy to the third degree mak-ing it completely rational, and evaluate it using rational computationfor elements with energy larger than 10 . We round the computedvalue of its third degree to the 64-bit floating point representation,and then compute the cubic root. The rational computation is moreaccurate (and permutation invariant) but significantly slower. How-ever, the cases of precision loss in the energy are rare, and theoverall computational overhead is negligible. This change has a ma-jor effect on the speed and effectiveness of our mesh optimization(Figure 10), avoiding unnecessary refinement and decreasing theoverall runtime.The mesh improvement terminates when either a user-specifiedmesh quality or a user-controlled number of iterations is reached.To ensure a fair comparison, for the large dataset testing and allexamples in the paper, we use the same stopping criteria (maxAMIPS energy is smaller than 10 or the number of optimizationiterations reaches 80) and input parameters (envelope size ϵ = − d ,targeted edge length ℓ = d /
20, where d is the diagonal’s length ofthe bounding box of the input mesh) as in [Hu et al. 2018].We use the same strategy as in TetWild for handling inputs withopen boundaries. We track the vertices on the open boundary andproject 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 ofinput with open boundary is shown in Figure 11. The output of the mesh improvement step is a volumetric tetrahedralmesh of the expanded bounding box of the input triangle soup, withthe preprocessed input triangles inserted. We provide two waysof 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 [Huet al. 2018]. The second application is a volumetric extension of themesh arrangement algorithm [Zhou et al. 2016]. In this case, theinput becomes a set of triangle soups, coupled with a set of Booleanoperations to perform on them. During the triangle insertion stage,we keep track of the provenance of each triangle, and use it atthe end to compute a set of generalized winding numbers (one foreach tracked input surface) for the centroids of all tetrahedra in thevolumetric mesh. We use the set of winding numbers to decide whichtetrahedron to keep by checking if it is supposed to be contained inthe result of the Boolean operation. For instance, when intersectingtwo triangle soups, we keep all tetrahedra that are inside both inputtriangle soups, according to the winding number definition.There are three major advantages of this method over [Zhouet al. 2016]: (1) Boolean operations can be performed on non-PWNsurfaces, (2) the output is equipped with a tetrahedral mesh, whichcould be useful in downstream applications, and (3) the surfacequality is high since the algorithm is allowed to remesh within the ϵ envelope. Our algorithm is implemented in C++ and uses Eigen [Guennebaudet al. 2010] for the linear algebra routines. We perform a large-scalecomparison of our method with other meshing methods on theThingi10k dataset [Zhou and Jacobson 2016], which contains 10 000real-world surface triangle meshes. We run our experiments oncluster nodes with a Xeon E5-2690v4 2.6GHz, allowing every modelto use up to 8 threads, 128GB memory, and 24 hours running time.The reference implementation and the scripts used to generate theresults are attached to the submission and will be released as anopen-source project.
With the above memory and time constraints, fTetWild success-fully tetrahedralizes 100% of the 10 000 input meshes (Figure 12).Most of the input models can be tetrahedralized with less than 1GBof RAM as detailed in Figure 13. Note that very complex modelsmight require more memory, for instance the one in Figure 21 usesaround 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 summarizethe results on the whole Thingi10k dataset in Table 2. Note thatonly our method and TetWild have high success rates: our averagetiming is however seven times faster than TetWild. , Vol. 1, No. 1, Article . Publication date: January 2020.
Fig. 12. 1000 random samples of fTetWild output on Thingi10k dataset.
Memory usage of fTetWild (MB)
Fig. 13. Histogram of memory usage of fTetWild over the Thingi10k dataset(data truncated at 2GB).Table 2. Comparison of code robustness and performance on the Thingi10kdataset.
Method Successrate Out ofmemory( > > CGAL
TetGen
TetWild *99.89% 0.05% 0.11% 0.00% 360.0
Ours **99.97% 0.02% 0.03% 0.00% 49.8
Note:
The maximum resources allowed for each model are 3 hours and 32GB ofmemory. The first 3 lines of data are from [Hu et al. 2018], Table 2. Note that theaverage time (last column) is computed over all the models for which each methodsucceeded, and it is thus not directly comparable between different methods. *:TetWild exceeds the 3h time on 11 models. If 27 hours of maximal running time areallowed, TetWild achieves 100% success rate. **: Our method exceeds the 3h time limiton 3 models. If 11 hours of maximal running time are allowed, fTetWild achieves100% success rate.
Thingi10k Dataset (10000 Models).
We compare the running timeof our method with TetWild. For a fair comparison, we disable ourcode optimizations that could be easily ported to TetWild, such asparallelization of the preprocessing and smoothing step, and usingthe recent fast winding number algorithm for the final filtering.Without these optimizations, our algorithm is 4 times faster thanTetWild on average (80.4s vs 360s). With code optimizations, wefurther improve our running time to 49.8s on average on a machinewith 8 cores, which is 7 times faster than the serial implementation >60s >120s >240s >480s >960s0%20%40%60%
TetWildTetWild Ours unoptimizedOurs unoptimized Ours optimizedOurs optimized
Fig. 14. Percentage of models requiring more than a certain time for ourparallel and serial algorithm compared with TetWild on Thingi10k dataset.
Input
Fig. 15. Example of a challenging model where fTetWild is 17 times fasterthan TetWild. of TetWild (Figure 14). On more complex examples, like the modelin Figure 15, our method is up to 17 times faster than TetWild.
Reduced Thingi10k Dataset (4540 Models).
We use a reduced datasetcontaining the intersection of the Thingi10k models that TetGen,CGAL, TetWild, and our method all succeed on. The dataset con-tains 4540 models, and allows us to fairly compare the performanceof the different methods. On average, our method is comparable(18.5s) to the widely used, Delaunay-based tetrahedral mesher Tet-Gen (22s), and is faster than CGAL (95s) and TetWild (107s), whilerobustly handling imperfect inputs. Figure 1 shows the number ofmodels requiring more than a given time. For example, within lessthan 2 minutes, our method successfully tetrahedralizes 98.7% ofthe inputs. It is interesting to note that the tail of the distribution ofour method is shorter than both TetGen and CGAL. For instance,there are only 4 models where our method requires more than 16minutes, differently from TetGen, CGAL, and TetWild which have20, 122, and 25 models, respectively. , Vol. 1, No. 1, Article . Publication date: January 2020. Input
Fig. 16. Our method (right) produces high-quality tet-meshes that aresimilar to TetWild (middle).
The geometric quality of meshes produced by our algorithm is simi-lar to the meshes produced by TetWild (Figure 16), since our methodimplements a similar mesh optimization strategy. We quantitativelyevaluate and compare the element quality of TetWild and our outputusing five different measures:(1) AMIPS energy (Equation (1)), range [ , + ∞) , optimal 3,(2) Minimal dihedral angle, range ( , . ] , optimal 1.23,(3) Volume-to-edge ratio 6 √ V / ℓ , range ( , ] , optimal 1,(4) Aspect ratio (cid:112) / h min / ℓ max , range ( , ] , optimal 1,(5) Radius-to-edge ratio 2 √ r in / ℓ max , range ( , ] , optimal 1,where V is the volume, ℓ max is the longest edge, h min the minimumheight, and r in the radius of the inscribed circle of a tetrahedron T . We use (3), (4) and (5) since these are standard measures fortetrahedral quality [Shewchuk 2002b].Figure 17 shows the histograms of worst and average elementquality of 10 000 output meshes of TetWild and our method. Thequality of our outputs are quite similar to TetWild’s output. Werefer to the study in [Hu et al. 2018, Figure 14] for the full qualitycomparison of TetWild and other tetrahedral meshing algorithms. Compared with TetWild, our method generates meshes of similardensity (Figure 18). Both TetWild and our method aim to generateas-coarse-as possible meshes while preserving the input surface.This choice is useful in downstream applications to reduce their com-putational cost. Optionally, the algorithm supports a user-specifiedsizing field to increase the density if desired.In contrast to our method, TetGen preserves the input surface ge-ometry exactly and thus generates a dense tetrahedral mesh aroundthe surface if the input surface mesh is dense, as visible in the modelshown in Figure 1. CGAL approximates the surface by means ofan implicit function, but sometimes over-refines sharp features and
TetWild OursAverage AMIPS energyMaximum AMIPS energy (truncated at 11)Average smallest dihedral angleMinimum smallest dihedral angleAverage volume-to-edge ratioMinimum volume-to-edge ratioAverage aspect ratioMinimum aspect ratioAverage radius-to-edge ratioMinimum radius-to-edge ratio
Fig. 17. Histogram for mesh quality comparison of TetWild (red) and ourmethod (green) in five different quality measures. The statistic is based onthe output of the whole Thingi10k dataset. , Vol. 1, No. 1, Article . Publication date: January 2020.
Comparison of mesh size between TetWild and fTetWild TetWild Ours
Fig. 18. Histograms of number of tetrahedra in log scale for the outputmeshes of Thingi10k dataset.
Input
Fig. 19. Example of repairing an invalid triangular mesh (left) with MeshFix(middle) and our algorithm (right). MeshFix is fast but loses details duringprocessing, while our method preserves them. The max AMIPS energy ofour intermediate tetrahedral mesh is 1975. Here we stop mesh improvementwhen maximum energy reaches 2000. tiny artifacts as illustrated in Figure 1, where the dark spots areover-refined regions.
Similarly to TetWild, our algorithm can be used to repair imper-fect triangle meshes by tetrahedralizing the volume and extractingthe surface of the generated tetrahedral mesh. However, the meshimprovement step of our method (Section 3.5) can be stopped atany time since we maintain an inversion-free floating point tetra-hedral mesh at all stages of our algorithm. High tetrahedral meshquality is not required for this application, and we can stop meshoptimization as soon as all input faces are inserted, further reducingthe running time. We compared our result with the state-of-the-artmesh repairing tool MeshFix [Attene 2010] in Figure 19. Our method,while slower, provides a higher-quality result with controllable geo-metric error. A minor, yet important, observation is that keepingonly the boundary of a valid tetrahedral mesh might generate anon-manifold surface mesh (Figure 20). To avoid this problem, weidentify the non-manifold edges and split them. Then we duplicateevery non-manifold vertex to ensure a global manifold output, us-ing the algorithm proposed in [Attene et al. 2009]. Note that thisprocedure ensures manifoldness, but introduces vertices in the samegeometric position. With this minor change, our algorithm can be
Input
Fig. 20. Example of a non-manifold surface mesh (left) which is automati-cally repaired by our algorithm (right second). used to repair triangle meshes, guaranteeing the extraction of anhigh-quality, manifold boundary surface mesh within the prescribeddistance from the input triangle soup.We also tested an extremely challenging model coming from anindustrial application in additive manufacturing (the part is copy-righted by Velo3D): the design of an exhaust pipe using a volumefilled with a structure based on the gyroid triply periodic minimalsurface. The model has a multitude of issues introduced during themodeling phase, but it can be cleaned up by our algorithm within 55minutes (or 122 minutes with the envelope size decreased by a factorof two), compared to around two weeks of manual labor requiredby Velo3D’s current processing pipeline. Our output mesh (Figure21) is directly usable for FEA, further editing, or fabrication. As areference point, the original implementation of TetWild takes 215minutes with a default envelope size. Another challenging modelwe tested contains complex thin structures coming from architec-ture (Figure 22). The method in [Masoud 2016; Tabatabaei Ghomiet al. 2018] optimizes for the layout of a graph, then replaces thegraph edges with cylinders of varying radii. To ensure solidity of thefinal structure, all cylinders are intersecting as shown in the closeup. Although the mesh contains many irregularities, fTetWildsuccessfully meshes the domain into an analysis-ready mesh.
Zhou et al. [2016] proposes to compute the arrangement betweenmultiple surfaces using an algorithm to map Boolean operationsinto simple algebraic expressions involving the winding numberof the input surfaces. Their method is robust, but only supportsclean PWN surfaces as input. We propose a simple extension of thisalgorithm (as explained in Section 3.6) to arbitrary triangle soups.The advantages of our method is evident when the input surfacescome from CAD models containing small gaps or self-intersections:both Mesh Arrangements [Zhou et al. 2016] and CGAL [Hachen-berger and Kettner 2019] are unable to perform the operation (sinceit is not well-defined for non-PWN surface), while fTetWild cancompute an approximate (since it allows for an ϵ -deviation from the , Vol. 1, No. 1, Article . Publication date: January 2020. Envelope size = 1e-3 d, d, Fig. 21. Meshing a complex model with 93 million vertices and 31 millionfaces with different envelope sizes (top). The input mesh contains degeneratetriangles and severe self-intersections. Our output tetrahedral meshes arein geometric high quality with either default envelope size (middle) or halfenvelope size (bottom).
Input
Fig. 22. Example of an architectural application with
80 999 self-intersectingfaces. The cylinders in the input are intersecting with each other as shownin the closeup. fTetWild successfully cleaned and tetrahedralized this input.Here we stop mesh optimization when maximum energy reaches 50. input surfaces) union, difference, and intersection between them(Figures 23, 24), providing robust (but slower) Boolean operationson imperfect geometries. The output is a tetrahedral mesh, whichcan be useful in downstream applications, and its boundary is ahigh quality surface triangular mesh.
The main application of tetrahedral meshing is physical simula-tions, and the high-quality of our results makes them ideal to bedirectly used in downstream finite element software (Figure 25).Additionally, the recently proposed a priori p -refinement [Schnei-der et al. 2018] is an ideal fit for our approach when targeting FEMapplications, since fTetWild always produces a valid floating-pointmesh. Schneider et al. [2018] provides a simple formula to determinethe order of each element to compensate for its, possibly bad, shape.We can use this criterion to terminate the mesh optimization earlyin our algorithm (thus reducing the meshing time) without affectingthe quality of the simulation, Figure 26.We use the Boolean difference (Section 5.2) to generate the back-ground mesh required for simulating the fluid flow on a cylindricaltube containing an obstacle (Figure 27). We introduced fTetWild, a novel robust tetrahedral meshing algo-rithm for triangle soups which combines the robustness of TetWildwith a running time comparable to Delaunay-based methods. Theimproved performance makes this algorithm suitable not only for ap-plications requiring a volumetric discretization, but also for surfacemesh repair and Boolean operations.Our current naive parallelization approach shows that our algo-rithm benefits from shared-memory parallelization; exploring moreadvanced parallelization techniques and extending it to distributedcomputation on HPC clusters are important directions for futurework. Our iterative triangle insertion algorithm could be used indynamic remeshing tasks, potentially allowing to reuse an existingmesh and insert new faces only in regions with high deformation.While conceptually trivial, extending our algorithm to 2D trianglemeshing could improve the performance of [Hu et al. 2019].Our algorithm uses the conformal AMIPS energy [Rabinovichet al. 2017] to measure and optimize the quality of the tetrahedra.An interesting alternative has been introduced concurrently to ourwork by [Alexa 2019]: they propose to optimize directly for theDirichlet energy of the tetrahedralization and show that this mea-sure is effective at removing slivers, while being computationallyefficient to evaluate. A comparative study of the two measures wouldbe interesting, and using the Dirichlet energy could lead to furtherreductions in the running time of our method.
REFERENCES
L. A. Freitag and C. Ollivier-Gooch. 1998. Tetrahedral Mesh Improvement UsingSwapping and Smoothing.
Internat. J. Numer. Methods Engrg.
40 (05 1998). https://doi.org/10.1002/(SICI)1097-0207(19971115)40:213.0.CO;2-9F. Alauzet and D. Marcum. 2014. A Closed Advancing-Layer Method With ChangingTopology Mesh Movement for Viscous Mesh Generation. In
Proceedings of the 22ndInternational Meshing Roundtable . Springer International Publishing, Cham, 241–261.https://doi.org/10.1007/978-3-319-02335-9_14, Vol. 1, No. 1, Article . Publication date: January 2020.
Two objects for Boolean Input
Fig. 23. Three Boolean operations computed on non-manifold, self-intersecting, and non-PWN input surface meshes. The left are two objects for Booleanoperation. The middle is the input surface mesh of fTetWild. The right are our output meshes after computing the union, difference, and intersection betweenthe two objects. The average max AMIPS energy of outputs and average time of different operations are with small variance.
U UU- ))(( =
Fig. 24. Four Boolean operations among 5 objects. fTetWild takes 34s andproducts output with
Input
Fig. 25. Example of non-linear elastic deformation of a body (right).
Input ≤ , 107s p ≤ , 69s Fig. 26. Two different stopping criteria of our algorithm. The full optimiza-tion (middle) improves the mesh to high quality, while using the criterionin [Schneider et al. 2018] (right) results in lower mesh quality but fastermeshing and smaller mesh size. The color shows the solution of the volu-metric Laplace equation.
M. Alexa. 2019. Harmonic Triangulations.
ACM Transactions on Graphics (Proceedingsof SIGGRAPH)
38, 4 (2019), 54. Input
Fig. 27. Streamlines of a fluid (right) moving in a cylindrical pipe (left top)with a complicated obstacle (left bottom) in the center. The backgroundmesh (middle) is obtained by subtracting the obstacle from a cylinder usingour method.
P. Alliez, D. Cohen-Steiner, M. Yvinec, and M. Desbrun. 2005a. Variational TetrahedralMeshing.
ACM Transactions on Graphics
24, 3 (07 2005), 617. https://doi.org/10.1145/1073204.1073238P. Alliez, D. Cohen-Steiner, M. Yvinec, and M. Desbrun. 2005b. Variational TetrahedralMeshing.
ACM Trans. Graph.
24, 3 (July 2005), 617–625. https://doi.org/10.1145/1073204.1073238M. Attene. 2010. A lightweight approach to repairing digitized polygon meshes.
The Visual Computer
26, 11 (01 Nov 2010), 1393–1406. https://doi.org/10.1007/s00371-010-0416-3M. Attene. 2014. Direct Repair of Self-intersecting Meshes.
Graph. Models
76, 6 (Nov.2014), 658–668. https://doi.org/10.1016/j.gmod.2014.09.002M. Attene. 2017.
ImatiSTL - Fast and Reliable Mesh Processing with a Hybrid Kernel .Springer Berlin Heidelberg, Berlin, Heidelberg, 86–96.M. Attene, M. Campen, and L. Kobbelt. 2013. Polygon Mesh Repairing: An ApplicationPerspective.
ACM Comput. Surv.
45, 2, Article 15 (March 2013), 33 pages. https://doi.org/10.1145/2431211.2431214M. Attene, D. Giorgi, M. Ferri, and B. Falcidieno. 2009. On converting sets of tetrahedrato combinatorial and PL manifolds.
Computer Aided Geometric Design
26, 8 (2009),850 – 864. https://doi.org/10.1016/j.cagd.2009.06.002F. Aurenhammer. 1991. Voronoi Diagrams&Mdash;a Survey of a Fundamental Geo-metric Data Structure.
ACM Comput. Surv.
23, 3 (Sept. 1991), 345–405. https://doi.org/10.1145/116873.116880F. Aurenhammer, R. Klein, and D.-T. Lee. 2013.
Voronoi Diagrams and Delaunay Trian-gulations
Discrete & Computational Geometry
3, 2 (01 Jun 1988), 147–168. https://doi.org/10.1007/BF02187904G. Barill, N. Dickson, R. Schmidt, D. I. Levin, and A. Jacobson. 2018. Fast WindingNumbers for Soups and Clouds.
ACM Transactions on Graphics
37, 4 (2018), 43:1–43:12.H. Barki, G. Guennebaud, and S. Foufou. 2015. Exact, robust, and efficient regularizedBooleans on general 3D meshes.
Computers and Mathematics with Applications J. Comput.System Sci.
48, 3 (1994), 384 – 409. https://doi.org/10.1016/S0022-0000(05)80059-5G. Bernstein. 2013. Cork Boolean Library . https://github.com/gilbo/cork.G. Bernstein and D. Fussell. 2009. Fast, Exact, Linear Booleans. In
SGP . EurographicsAssociation, Aire-la-Ville, Switzerland, Switzerland, 1269–1278.H. Bieri and W. Nef. 1988. Elementary Set Operations with D-dimensional Polyhedra.In
Proc. IWCGA . Springer-Verlag, Berlin, Heidelberg, 97–112.C. J. Bishop. 2016. Nonobtuse Triangulations of PSLGs.
Discrete & ComputationalGeometry
56, 1 (2016), 43–92.J.-D. Boissonnat, O. Devillers, S. Pion, M. Teillaud, and M. Yvinec. 2002. Triangulationsin CGAL.
Computational Geometry
22 (2002), 5–19.J.-D. Boissonnat and S. Oudot. 2005. Provably Good Sampling and Meshing of Surfaces.
Graphical Models
67, 5 (09 2005), 405–451. https://doi.org/10.1016/j.gmod.2005.01.004R. Bridson and C. Doran. 2014. Quartet: A tetrahedral mesh generator that doesisosurface stuffing with an acute tetrahedral tile. https://github.com/crawforddoran/quartet.J. R. Bronson, J. A. Levine, and R. T. Whitaker. 2013. Lattice Cleaving: ConformingTetrahedral Meshes of Multimaterial Domains With Bounded Quality. In
Proceedingsof the 21st International Meshing Roundtable . Springer Berlin Heidelberg, Berlin,Heidelberg, 191–209. https://doi.org/10.1007/978-3-642-33573-0_12O. Busaryev, T. K. Dey, and J. A. Levine. 2009. Repairing and Meshing Imperfect Shapeswith Delaunay Refinement. In . ACM, New York, NY, USA, 25–33. https://doi.org/10.1145/1629255.1629259M. Campen and L. Kobbelt. 2010. Exact and Robust (self-)intersections for PolygonalMeshes.
Comput. Graph. Forum
29, 2 (2010), 397–406.S. A. Canann, S. N. Muthukrishnan, and R. K. Phillips. 1996. Topological refinementprocedures for triangular finite element meshes.
Engineering with Computers
12, 3(01 Sep 1996), 243–255. https://doi.org/10.1007/BF01198738S. A. Canann, M. B. Stephenson, and T. Blacker. 1993. Optismoothing: An optimization-driven approach to mesh smoothing.
Finite Elements in Analysis and Design
13, 2(1993), 185 – 190. https://doi.org/10.1016/0168-874X(93)90056-VL. Chen and J.-c. Xu. 2004. Optimal Delaunay Triangulations.
Journal of ComputationalMathematics
22, 2 (2004), 299–308.S.-W. Cheng, T. K. Dey, and J. A. Levine. 2008. A Practical Delaunay Meshing Algorithmfor a Large Class of Domains. In
Proceedings of the 16th International MeshingRoundtable . Springer, Springer Berlin Heidelberg, Berlin, Heidelberg, 477–494.S.-W. Cheng, T. K. Dey, and J. Shewchuk. 2012.
Delaunay Mesh Generation . Chapmanand Hall/CRC, Boca Raton, Florida.L. P. Chew. 1989. Constrained delaunay triangulations.
Algorithmica
4, 1 (01 Jun 1989),97–108. https://doi.org/10.1007/BF01553881L. P. Chew. 1993. Guaranteed-Quality Mesh Generation for Curved Surfaces. In
Pro-ceedings of the ninth annual symposium on Computational geometry - SCG ’93 . ACMPress, New York, NY, USA, 274–280. https://doi.org/10.1145/160985.161150D. Cohen-Steiner, E. C. de Verdière, and M. Yvinec. 2002. Conforming DelaunayTriangulations in 3D. In
Proceedings of the eighteenth annual symposium on Com-putational geometry - SCG ’02 . ACM Press, New York, NY, USA, 217 – 233. https://doi.org/10.1145/513400.513425J.-C. Cuilliere, V. Francois, and J.-M. Drouet. 2013. Automatic 3D Mesh Generation ofMultiple Domains for Topology Optimization Methods. In
Proceedings of the 21stInternational Meshing Roundtable . Springer Berlin Heidelberg, Berlin, Heidelberg,243–259. https://doi.org/10.1007/978-3-642-33573-0_15T. K. Dey and J. A. Levine. 2008. Delpsc: A Delaunay Mesher for Piecewise SmoothComplexes. In
Proceedings of the twenty-fourth annual symposium on Computationalgeometry - SCG ’08 . ACM Press, New York, NY, USA, 220–221. https://doi.org/10.1145/1377676.1377712A. Doi and A. Koide. 1991. An efficient method of triangulating equi-valued surfacesby using tetrahedral cells.
IEICE TRANSACTIONS on Information and Systems
74, 1(1991), 214–224.C. Doran, A. Chang, and R. Bridson. 2013. Isosurface Stuffing Improved: Acute Latticesand Feature Matching. In
ACM SIGGRAPH 2013 Talks on - SIGGRAPH ’13 . ACMPress, New York, NY, USA, 38:1–38:1. https://doi.org/10.1145/2504459.2504507M. Douze, J.-S. Franco, and B. Raffin. 2015.
QuickCSG: Arbitrary and Faster BooleanCombinations of N Solids . Technical Report 01121419. Inria Research Centre Grenoble,Rhone-Alpes.Q. Du and D. Wang. 2003. Tetrahedral Mesh Generation and Optimization Based onCentroidal Voronoi Tessellations.
International journal for numerical methods inengineering
56, 9 (2003), 1355–1373.N. Faraj, J.-M. Thiery, and T. Boubekeur. 2016. Multi-Material Adaptive VolumeRemesher.
Compurer and Graphics Journal (proc. Shape Modeling International2016)
58 (2016), 150 – 160.L. Feng, P. Alliez, L. Busé, H. Delingette, and M. Desbrun. 2018. Curved OptimalDelaunay Triangulation.
ACM Trans. Graph.
37, 4 (2018), 61:1–61:16.X.-M. Fu, Y. Liu, and B. Guo. 2015. Computing Locally Injective Mappings by AdvancedMIPS.
ACM Trans. Graph.
34, 4, Article 71 (July 2015), 12 pages. https://doi.org/10. 1145/2766938J. A. George. 1971.
Computer Implementation of the Finite Element Method . Ph.D.Dissertation. Stanford University, Stanford, CA, USA. AAI7205916.P. L. George, H. Borouchaki, and E. Saltel. 2003. ’Ultimate’ robust-ness in meshing an arbitrary polyhedron.
Internat. J. Numer. Meth-ods Engrg.
58, 7 (2003), 1061–1089. https://doi.org/10.1002/nme.808arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.808M. Granados, P. Hachenberger, S. Hert, L. Kettner, K. Mehlhorn, and M. Seel. 2003.Boolean operations on 3D selective Nef complexes: Data structure, algorithms,and implementation. In
Proc. ESA . Springer Berlin Heidelberg, Berlin, Heidelberg,654–666.G. Guennebaud, B. Jacob, et al. 2010. Eigen v3.P. Guigue and O. Devillers. 2003. Fast and Robust Triangle-Triangle Overlap TestUsing Orientation Predicates.
Journal of graphics tools
8, 1 (2003), 39–52. https://doi.org/10.1080/10867651.2003.10487580P. Hachenberger and L. Kettner. 2019. 3D Boolean Operations on Nef Polyhedra.In
CGAL User and Reference Manual (4.14 ed.). CGAL Editorial Board. https://doc.cgal.org/4.14/Manual/packages.html
Proceedings of the 22ndInternational Meshing Roundtable . Springer International Publishing, Cham, 75–91.https://doi.org/10.1007/978-3-319-02335-9_5K. Hu, D. Yan, D. Bommes, P. Alliez, and B. Benes. 2017. Error-Bounded and FeaturePreserving Surface Remeshing with Minimal Angle Improvement.
IEEE Transactionson Visualization and Computer Graphics
23, 12 (Dec 2017), 2560–2573. https://doi.org/10.1109/TVCG.2016.2632720Y. Hu, T. Schneider, X. Gao, Q. Zhou, A. Jacobson, D. Zorin, and D. Panozzo. 2019.TriWild: Robust Triangulation with Curve Constraints.
ACM Trans. Graph. (2019).Y. Hu, Q. Zhou, X. Gao, A. Jacobson, D. Zorin, and D. Panozzo. 2018. TetrahedralMeshing in the Wild.
ACM Trans. Graph.
37, 4, Article 60 (July 2018), 14 pages.https://doi.org/10.1145/3197517.3201353C. Jamin, P. Alliez, M. Yvinec, and J.-D. Boissonnat. 2015. CGALmesh: A GenericFramework for Delaunay Mesh Generation.
ACM Trans. Math. Software
41, 4 (102015), 1–24. https://doi.org/10.1145/2699463F. Labelle and J. R. Shewchuk. 2007. Isosurface Stuffing: Fast Tetrahedral Meshes WithGood Dihedral Angles. In
ACM SIGGRAPH 2007 papers on - SIGGRAPH ’07 . ACMPress, New York, NY, USA, 57. https://doi.org/10.1145/1275808.1276448B. Lévy. 2019. Geogram. http://alice.loria.fr/index.php/software/4-library/75-geogram.html.Y. Lipman. 2012. Bounded Distortion Mapping Spaces for Triangular Meshes.
ACMTrans. Graph.
31, 4 (2012), 108.W. E. Lorensen and H. E. Cline. 1987. Marching Cubes: A High Resolution 3D SurfaceConstruction Algorithm.
SIGGRAPH Comput. Graph.
21, 4 (Aug. 1987), 163–169.https://doi.org/10.1145/37402.37422S. V. Magalhães, W. R. Franklin, and M. V. Andrade. 2017. Fast exact parallel 3D meshintersection algorithm using only orientation predicates. In
Proceedings of the 25thACM SIGSPATIAL International Conference on Advances in Geographic InformationSystems . ACM, ACM, New York, NY, USA, 44.M. Mandad, D. Cohen-Steiner, and P. Alliez. 2015. Isotopic Approximation Withina Tolerance Volume.
ACM Trans. Graph.
34, 4, Article 64 (July 2015), 12 pages.https://doi.org/10.1145/2766950A. Masoud. 2016.
3D Graphical Statics Using Reciprocal Polyhedral Diagrams . Ph.D.Dissertation. ETH Zruich, Stefano Franscini Platz 5, Zurich, CH, 8093.N. Molino, R. Bridson, and R. Fedkiw. 2003. Tetrahedral Mesh Generation for DeformableBodies. In
Proc. Symposium on Computer Animation .M. Murphy, D. M. Mount, and C. W. Gable. 2001. A Point-Placement Strategyfor Conforming Delaunay Tetrahedralization.
International Journal of Compu-tational Geometry & Applications
11, 06 (12 2001), 669–682. https://doi.org/10.1142/s0218195901000699K. Museth, D. E. Breen, R. T. Whitaker, and A. H. Barr. 2002. Level set surface editingoperators.
ACM Trans. Graph.
21, 3 (2002), 330–338.B. Naylor, J. Amanatides, and W. Thibault. 1990. Merging BSP trees yields polyhedralset operations. In
Proc. SIGGRAPH . ACM, New York, NY, USA, 115–124.A. Paoluzzi, V. Shapiro, and A. DiCarlo. 2017. Arrangements of cellular complexes.
CoRR abs/1704.00142 (2017). arXiv:1704.00142 http://arxiv.org/abs/1704.00142D. Pavic, M. Campen, and L. Kobbelt. 2010. Hybrid Booleans.
Comput. Graph. Forum
J. Comput. Phys.
72, 2 (Oct. 1987), 449–466.https://doi.org/10.1016/0021-9991(87)90093-3M. Rabinovich, R. Poranne, D. Panozzo, and O. Sorkine-Hornung. 2017. ScalableLocally Injective Mappings.
ACM Trans. Graph.
36, 2 (April 2017), 16. https://doi.org/10.1145/2983621J.-F. Remacle. 2017. A Two-Level Multithreaded Delaunay Kernel.
Computer-AidedDesign
85 (04 2017), 2–9. https://doi.org/10.1016/j.cad.2016.07.018J. Ruppert. 1995. A Delaunay Refinement Algorithm for Quality 2-Dimensional MeshGeneration.
Journal of Algorithms
18, 3 (05 1995), 548–585. https://doi.org/10.1006/, Vol. 1, No. 1, Article . Publication date: January 2020. jagm.1995.1021E. A. Sadek. 1980. A scheme for the automatic generation of tri-angular finite elements.
Internat. J. Numer. Methods Engrg.
ACM SIGGRAPH 2010 Talks . ACM, ACM, New York, NY, USA, 6.T. Schneider, Y. Hu, J. Dumas, X. Gao, D. Panozzo, and D. Zorin. 2018. Decouplingsimulation accuracy from mesh quality.
ACM Transactions on Graphics
37, 6 (dec2018), 1–14. https://doi.org/10.1145/3272127.3275067M. Schweiger and S. Arridge. 2016. Basis mapping methods for forward and inverseproblems: BASIS MAPPING METHODS.
Internat. J. Numer. Methods Engrg.
109 (052016). https://doi.org/10.1002/nme.5271D. R. Sheehy. 2012. New Bounds on the Size of Optimal Meshes.
Computer GraphicsForum
31, 5 (08 2012), 1627–1635. https://doi.org/10.1111/j.1467-8659.2012.03168.xB. Sheng, P. Li, H. Fu, L. Ma, and E. Wu. 2018a. Efficient non-incremental constructivesolid geometry evaluation for triangular meshes.
Graphical Models
97 (2018), 1–16.B. Sheng, B. Liu, P. Li, H. Fu, L. Ma, and E. Wu. 2018b. Accelerated robust Booleanoperations based on hybrid representations.
Computer Aided Geometric Design
Unstructured Mesh Generation . Chapman and Hall/CRC, Boca Raton,Florida, Chapter 10, 257 – 297.J. R. Shewchuk. 1996. Triangle: Engineering a 2D quality mesh generator and Delaunaytriangulator. In
Applied Computational Geometry Towards Geometric Engineering ,Ming C. Lin and Dinesh Manocha (Eds.). Springer Berlin Heidelberg, Berlin, Heidel-berg, 203–222.J. R. Shewchuk. 1997. Adaptive Precision Floating-Point Arithmetic and Fast RobustGeometric Predicates.
Discrete & Computational Geometry
18, 3 (Oct. 1997), 305–363.J. R. Shewchuk. 1998. Tetrahedral Mesh Generation by Delaunay Refinement. In
Proceedings of the fourteenth annual symposium on Computational geometry - SCG’98 . ACM Press, New York, NY, USA, 86–95. https://doi.org/10.1145/276884.276894J. R. Shewchuk. 1999. Lecture Notes on Delaunay Mesh Generation. (1999).J. R. Shewchuk. 2002a. Constrained Delaunay Tetrahedralizations and Provably GoodBoundary Recovery. In
Eleventh International Meshing Roundtable . Sandia NationalLaboratories, 193–204.J. R. Shewchuk. 2002b. What is a good linear element? interpolation, conditioning, andquality measures. In
In 11th International Meshing Roundtable . 115–126.H. Si. 2015. TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator.
ACM Trans.Math. Softw.
41, 2, Article 11 (Feb. 2015), 36 pages. https://doi.org/10.1145/2629697H. Si and K. Gartner. 2005. Meshing Piecewise Linear Complexes by Constrained Delau-nay Tetrahedralizations. In
Proceedings of the 14th international meshing roundtable .Springer, Springer Berlin Heidelberg, Berlin, Heidelberg, 147–163.H. Si and J. R. Shewchuk. 2014. Incrementally Constructing and Updating ConstrainedDelaunay Tetrahedralizations With Finite-Precision Coordinates.
Engineering withComputers
30, 2 (04 2014), 253–269. https://doi.org/10.1007/s00366-013-0331-0A. Tabatabaei Ghomi, M. Bolhassan, A. Nejur, and M. Akbarzadeh. 2018. Effect ofSubdivision of Force Diagrams on the Local Buckling, Load-Path and Material Useof Founded Forms. In
Proceedings of the IASS Symposium 2018, Creativity in StructuralDesign . MIT, Boston, USA.W. C. Thibault and B. F. Naylor. 1987. Set operations on polyhedra using binary spacepartitioning trees. In
Proc. SIGGRAPH . ACM, New York, NY, USA, 153–162.J. Tournois, C. Wormser, P. Alliez, and M. Desbrun. 2009. Interleaving DelaunayRefinement and Optimization for Practical Isotropic Tetrahedron Mesh Generation.
ACM Transactions on Graphics
28, 3 (07 2009), 1. https://doi.org/10.1145/1531326.1531381G. Varadhan, S. Krishnan, T. Sriram, and D. Manocha. 2004. Topology preservingsurface extraction using adaptive subdivision. In
SGP . ACM, New York, NY, USA,235–244.C. C. L. Wang. 2011. Approximate Boolean Operations on Large Polyhedral Solids withPartial Mesh Reconstruction.
IEEE Trans. Vis. Comput. Graph.
17, 6 (2011), 836–849.N. P. Weatherill and O. Hassan. 1994. Efficient three-dimensional Delaunay triangulationwith automatic point creation and imposed boundary constraints.
Internat. J. Numer.Methods Engrg.
37, 12 (1994), 2005–2039. https://doi.org/10.1002/nme.1620371203arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620371203R. Wein, E. Berberich, E. Fogel, D. Halperin, M. Hemmer, O. Salzman, and B. Zuk-erman. 2018. 2D Arrangements. In
CGAL User and Reference Manual (4.13ed.). CGAL Editorial Board. https://doc.cgal.org/4.13/Manual/packages.html
IEEE Computer Graphics and Applications
3, 1 (Jan 1983), 39–46.https://doi.org/10.1109/MCG.1983.262997H. Zhao, C. C. Wang, Y. Chen, and X. Jin. 2011. Parallel and efficient Boolean onpolygonal solids.
The Visual Computer
27, 6-8 (2011), 507–517.Q. Zhou, E. Grinspun, D. Zorin, and A. Jacobson. 2016. Mesh Arrangements for SolidGeometry.
ACM Transactions on Graphics (TOG)
35, 4 (2016), 39. Q. Zhou and A. Jacobson. 2016. Thingi10K: A Dataset of 10, 000 3D-Printing Models.
CoRR abs/1605.04797 (2016). arXiv:1605.04797
A A BRIEF DESCRIPTION OF THE TETWILDALGORITHM
The TetWild algorithm [Hu et al. 2018] takes a 3D triangle soupas input and generates a tetrahedral mesh that (1) has no invertedor degenerate tetrahedra and (2) contains an approximation of theinput surface within a user-defined ϵ -envelope .The method starts with an initial background mesh generatedusing an unconstrained Delaunay tetrahedralization on the inputpoints plus an additional set of evenly-spaced points sampled froma regular grid. These additional points are added to improve theshape of the tetrahedra in the background mesh.This step generates tetrahedra that might not represent the inputfaces of the triangle soup: to ensure that they are preserved TetWilduses a Binary Space Partitioning (BSP) subdivision step. Each inputtriangle is converted into a plane that cuts the tetrahedra of thebackground mesh. The output of this stage is a polyhedral mesh. Toavoid numerical issues and to guarantee that the sub-elements inthe polyhedral mesh are convex and non-inverted, TetWild convertsthe coordinates of the vertices of the initial background mesh intorational numbers and performs all computations using rationalnumbers.Since any convex polyhedron can be trivially subdivided intotetrahedra by adding an additional point in its barycenter, a tetra-hedral mesh that exactly preserves input triangles can be naturallyobtained after BSP subdivision. However, the vertices of this tetrahe-dral mesh are represented in rational coordinated. Rounding themto floating point is not simple, since the BSP subdivision introducesbadly-shaped tetrahedra which could flip after rounding. TetWildthus increases the quality of the elements using an hybrid optimiza-tion procedures that mixes floating point and rational representa-tion.During mesh improvement, the preserved input triangle’s facesare tracked and are allowed to move inside the ϵ -envelope. The ϵ -envelope limits the tracked surface from deviating from the inputfurther than ϵ .TetWild uses four local operations for mesh improvement: (1)edge splitting, (2) edge collapsing, (3) edge swapping, and (4) vertexsmoothing. Every operation is rolled back if the tracked surfaceleaves the envelope after the operation or if any tetrahedra areinverted, ensuring a valid output. Differently from other mesh im-provement methods, TetWild uses the 3D conformal AMIPS energyfor measuring the geometric quality of the tetrahedra. The AMIPSenergy is scaling-invariant and easily differentiable, which booststhese traditional local operations.As the quality of the mesh is improved, the rational coordinatescan be gradually rounded into floating points. Theoretically, theremight be some unroundable vertices, but it does not occur on theten thousand models that TetWild has been tested on.The final step is the removal of the tetrahedra outside of thetracked surface. To handle potentially noisy inputs, TetWild com-putes the winding number of the centroids of all tetrahedra with , Vol. 1, No. 1, Article . Publication date: January 2020. respect to the tracked surface, and filters out all tetrahedra withcentroid’s winding-number larger than 0.5. B EXAMPLE OF UNSTABLE AMIPS ENERGY
If we compute the 3D AMIPS energy for a tetrahedron with these 4vertices v = ( . , . , . ) v = ( . , . , . ) v = ( . , . , . ) v = ( . , . , . ) we obtain AMIPS = . = . = . = . , where the subscript indicates the vertex permutations. There are24 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 obtainfluctuations AMIPS = . = . = . = . . As reference the correct value computed with rational number isAMIPS = . . C UNUSED DECOMPOSITIONS OF A TETRAHEDRON
We enumerated all the possible decompositions of a tetrahedronand discovered two symmetry classes (Figure 28) of triangulation offaces whose decomposition requires an additional internal vertex.Note that these two cases are never selected by our algorithm (weinclude them here for completeness), as our rule (Section 3.4.2) neverselects these two cases.We show that our rule does not select case 1 (Figure 28 left). Bycontradiction: since the configuration is selected then the edges [ p , v ] , [ p , v ] and [ p , v ] are present, thus v > v , v > v , and v > v , according to our rule. Combining these inequalities, theindices of the vertices must satisfy v > v > v > v , which isimpossible. Case 2 (Figure 28 right) is also not selected following asimilar argument. D AN EXAMPLE FOR OPEN-BOUNDARY EDGEPRESERVATION
If triangle T is the only inserted triangle and is entirely containedinside a tetrahedron T (Figure 29(1)), the intersection of the plane P and T will be a larger polygon (marked in yellow) containing T .In this case, the edges of T , which are open-boundary edges, are p p T p v v v v p p T p v v v v p Case 1 Case 2
Fig. 28. Two unused configurations requiring an additional vertex. T T P p p p p e (1) (2) (3) Fig. 29. Example for preserving an open-boundary edge e of triangle T . (1)Insert T and T I = {T } in this case. (2) The sub-tetrahedra of T after subdi-vision. (Only sub-tetrahedra behind T are shown for better visualization.)(3) Inserting edge e and get the intersection points (in green).Table 3. Edge-cut configurations of a cutting tetrahedron before and aftersnapping. 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) not preserved. To preserve them, we subdivide the tetrahedra oncemore.In Figure 29(1), T first get decomposed into sub-tetrahedra (Fig-ure 29(2)). Then the faces covering T are F = {[ p , p , p ] , [ p , p , p ]} Figure 29(3). The open-boundary edge e and the faces in F are pro-jected to the best-fitting plane of p , p , p , and p . The intersectionpoints of the projection of e and T are then computed in 2D andare lifted to 3D (3 green points in Figure 29(3)). Now there are 3edges [ p , p ] , [ p , p ] , [ p , p ] cut into two. We thus subdivide allthe neighbouring tetrahedra with the table-based subdivision. E CHANGES OF EDGE-CUT CONFIGURATION AFTERSNAPPING
Table 3 shows all possible edge-cut configurations of a cutting tetra-hedron T after snapping. The final configurations have no morethan two vertices which makes the triangulation of F uniquelydefined by the points. The table includes only the 4 symmetryclasses where T is cut by plane P and contains a face in F (Figure 9(2)(3)(5)(7)), but excludes the remaining 3 classes where T is not cutor is just affected by their neighbors (Figure 9 (1)(4)(6)).A tetrahedron T can have at most 3 vertices snapped. If T hasall its 4 vertices within a δ distance to the P , we only snap the 3vertices closer to P ..