Tetrahedral Mesh Generation by Delaunay Renement

Jonathan Richard Shewchuk

School of Computer Science Carnegie Mellon University Pittsburgh, Pennsylvania 15213 jrs@cs.cmu.edu


Given a complex of vertices, constraining segments, and planar straight-line constraining facets in  , with no input angle less than   , an algorithm presented herein can generate a conforming mesh of Delaunay tetrahedra whose circumradius-to-shortest edge ratios are no greater than two. The sizes of the tetrahedra can provably grade from small to large over a relatively short distance. An implementation demonstrates that the algorithm generates excellent meshes, generally surpassing the theoretical bounds, and is effective in eliminating tetrahedra with small or large dihedral angles, although they are not all covered by the theoretical guarantee.

!#"%$ '&)(10243 & $

Meshes of triangles or tetrahedra have many applications, including interpolation, rendering, and numerical methods such as the nite element method. Most such applications demand more than just a triangulation of the object or domain being rendered or simulated. To ensure accurate results, the triangles or tetrahedra must be wellshaped, having small aspect ratios or bounds on their smallest and largest angles. Mesh generation algorithms based on Delaunay renement are effective both in theory and in practice. Delaunay renement algorithms operate by maintaining a Delaunay or constrained Delaunay triangulation, which is rened by inserting carefully placed vertices until the mesh meets constraints on triangle quality and size. It is difcult to trace who rst used Delaunay triangulations for nite element meshing, and equally difcult to tell where the suggestion arose to use the triangulation to guide vertex creation. These ideas have been intensively studied in the engineering community since the mid-1980s; for instance, Frey [7] eliminates poorly shaped triangles from a triangulation by inserting new vertices at their circumcenters (dened in Section 2), whereas Weatherill [17] inserts new vertices at their centroids. These ideas bore vital theoretical fruit when the problem of mesh generation began to attract interest from the computational
Supported in part by the Advanced Research Projects Agency and Rome Laboratory, Air Force Materiel Command, USAF under agreement number F30602-96-1-0287, in part by the National Science Foundation under Grant CMS-9318163, and in part by a grant from the Intel Corporation.

geometry community in the early 1990s. The rst provably good Delaunay renement algorithm is due to Paul Chew [2], and takes as its input a set of vertices and segments that dene the region to be meshed. By inserting additional vertices, Chews algorithm generates a two-dimensional constrained Delaunay triangulation whose   and 687   . Dey, Bajaj, and Sugangles are bounded between 5 ihara [5] and Chew [4] generalize Chews algorithm to three dimensions, but only for unconstrained point set inputs. All three algorithms produce uniform meshes, whose triangles or tetrahedra are of roughly the same size. Uniform meshes sometimes have many more triangles or tetrahedra than are necessary, and thus impose an excessive computational load upon the applications that make use of them. Jim Ruppert [13] and Paul Chew [3] have each proposed two-dimensional Delaunay renement algorithms that produce meshes of well-shaped triangles whose sizes are graded, and Ruppert has furthermore proven that his algorithm produces meshes that are nicely graded in a theoretical sense described in Section 7. These algorithms are successful because they exploit several favorable characteristics of Delaunay triangulations. Delaunay triangulations have been extensively studied, and good algorithms are available. Inserting a vertex is a local operation, and is inexpensive except in unusual cases. In two dimensions, Delaunay triangulations maximize the minimum angle, compared with all other triangulations of the same vertex set [8]. The greatest advantage of Delaunay triangulations is less obvious. The central question of any Delaunay renement algorithm is where should the next vertex be inserted? As Section 3 will demonstrate, a reasonable answer is as far from other vertices as possible. If a new vertex is inserted too close to another vertex, the resulting short edge will engender thin triangles or tetrahedra. Because a Delaunay triangle has no vertices in its circumcircle (and a Delaunay tetrahedron has no vertices in its circumsphere), a Delaunay triangulation is an ideal search structure for nding points that are far from other vertices. Herein, I build upon the algorithmic and analytical framework of Ruppert to design a new tetrahedral Delaunay renement algorithm. This algorithm generates meshes whose tetrahedra have circumradius-to-shortest edge ratios (dened shortly) no greater than the bound 9A@B7 . Upon relaxing 9 to be greater than two, one can also guarantee good grading. My algorithm is distinguished from those of Dey et al. and Chew by this guarantee of good grading (as opposed to uniform meshes), and by its ability to handle the rather general constrained inputs described in Section 4 (as opposed to unconstrained sets of vertices). The main theoretical deciency of the algorithm is its requirement that incident input segments and C  . (The requirements are facets be separated by angles of at least less onerous in practice.) This deciency can be partly ameliorated; Section 9 includes a brief summary of the steps involved.

t v v

Needle / Wedge

Cap

Sliver

Figure 1: Tetrahedra with poor angles. Needles and wedges have edges of greatly disparate length; caps have a large solid angle; slivers have neither, but can have good circumradius-to-shortest edge ratios. Two other tetrahedral mesh generation algorithms (not based on Delaunay renement) have provable bounds. The octree-based algorithm of Mitchell and Vavasis [11, 12] is a theoretical tour de force, obtaining provable bounds on dihedral angles and grading, but its bounds are too weak to offer any practical reassurance (and have not been explicitly stated). An algorithm by Miller, Talmor, Teng, Walkington, and Wang [10] generates its nal vertex set before triangulating it. Their algorithm has provable bounds similar to those of the algorithm described herein, but the algorithm herein has the opportunity to stop early if fewer vertices are needed than the theory suggests. As a result, the present algorithm typically uses fewer vertices by several orders of magnitude; details are provided elsewhere [14].

Figure 2: Any triangle whose circumradius-to-shortest edge ratio is larger than some bound is split by inserting a vertex at its circumcenter. The Delaunay property is maintained, and the triangle is thus eliminated. Every new edge has length at least times that of shortest edge of the poor triangle.

0 3  01'&  3  3  

radius-to-shortest edge ratios typically arise in small numbers in practice. As Section 8 will demonstrate, the worst slivers can often be removed by Delaunay renement, even if there is no theoretical guarantee. Meshes with bounds on the circumradius-toshortest edge ratios of their tetrahedra are an excellent starting point for mesh smoothing and optimization methods that remove slivers and otherwise improve the quality of an existing mesh [6]. Even if slivers are not removed, the Voronoi dual of a tetrahedralization with bounded circumradius-to-shortest edge ratios has nicely rounded cells, and is sometimes ideal for use in the control volume method [9].

PRQS UTVB " (EUW0 S 3 $ (YX`a 0 $ BUbcBd $ ae $

Miller, Talmor, Teng, and Walkington [9] have pointed out that the most natural and elegant measure for analyzing Delaunay renement algorithms is the circumradius-to-shortest edge ratio of a triangle or tetrahedron. The circumsphere of a simplex is the unique circle or sphere that passes through all its vertices. The circumcenter and circumradius of a simplex are the center and radius of its circumsphere, respectively. The quotient of a simplexs circumradius and the length of its shortest edge is the metric that is naturally optimized by Delaunay renement algorithms. One would like this ratio to be as small as possible. But does optimizing this metric aid practical applications? In two dimensions, the answer is yes. A triangles circumradius-tois related to its smallest angle by the shortest edge ratio @ 6 7 . The smaller a triangles ratio, the formula larger its smallest angle. Chews two-dimensional algorithms produce meshes whose triangles circumradius-to-shortest edge ratios are bounded below  one, and hence their angles are bounded be  and tween 5 687 . Rupperts algorithm can produce meshes 7 , and hence their anwhose triangles ratios bounded below  and  are 65  . gles range between 7 In three dimensions, however, a mesh of tetrahedra whose circumradius-to-shortest edge ratios are bounded is not entirely adequate for the needs of interpolation. As Dey et al. illustrate, most tetrahedra with poor angles have circumcircles much larger than their shortest edges, including the needle, wedge, and cap illustrated in Figure 1, but there is one type called a sliver or kite tetrahedron that does not. The canonical sliver is formed by arranging four vertices, equally spaced, around the equator of a sphere, then perturbing one of the vertices slightly off the equator. A sliver can 7 , yet be have a circumradius-to-shortest edge ratio as low as 6 considered awful by most other measures, because its volume and its shortest altitude can be arbitrarily close to zero, and its dihedral   and angles can be arbitrarily close to 6 . Despite slivers, Delaunay renement methods are valuable for generating three-dimensional meshes. Slivers having good circum-

!#"%$ 2"3 !#5"%47$ 698@& '0) 1BA ED F

&('0) 1

(G D H

%"C

%G

The central operation of Chews, Rupperts, Deys, and my own Delaunay renement algorithms is the insertion of a vertex at the circumcenter of a triangle or tetrahedron of poor quality. The Delaunay property is maintained, perhaps by the Bowyer/Watson algorithm for the incremental update of Delaunay triangulations [1, 16]. The poor simplex cannot survive, because its circumsphere is no longer empty. For brevity, the act of inserting a vertex at a simplexs circumcenter is called splitting a simplex. If poor simplices are split one by one, either all will eventually be eliminated, or the algorithm will run forever. The main insight of all these Delaunay renement algorithms is that Delaunay renement is guaranteed to terminate if the notion of poor quality includes only simplices that have a circumradiusto-shortest edge ratio larger than some appropriate bound 9 . The only new edges created by the Delaunay insertion of a vertex are edges connected to (see Figure 2). Because is the circumcenter of some Delaunay simplex , and there were no vertices inside the circumsphere of before was inserted, no new edge can be shorter than the circumradius of . Because has a circumradiusto-shortest edge ratio larger than 9 , every new edge has length at least 9 times that of the shortest edge of . Rupperts Delaunay renement algorithm [13] employs a bound of 9 @ 7 , and Chews second Delaunay renement algorithm [3] employs a bound of 9 @ 6 . Chews rst Delaunay renement algorithm [2] splits any triangle whose circumradius is greater than the length of the shortest edge in the entire mesh, thus achieving a bound of 9 @ 6 , but forcing all triangles to have uniform size. In the same manner, the three-dimensional algorithms of Dey et al. and Chew achieve a bound of 9 @ 7 . With these bounds, every new edge created is at least as long as some other edge already in the mesh. Hence, no vertex is ever inserted closer to another vertex than the length of the shortest edge in the initial triangulation. Delaunay renement must eventually terminate, because the augmented triangulation will run out of places to put vertices.

fg g

(a)

Figure 4: The orthogonal projections of points and sets of points onto facets and segments. ment is initially represented by one subsegment, until it is subdivided.) The bold edges of the tetrahedralization illustrated in Figure 3(b) are subsegments; other edges are not. Similarly, each facet is subdivided into triangular faces called subfacets. All of the triangular faces visible in Figure 3(b) are subfacets, but most of the faces in the interior of the tetrahedralization are not. Any vertex inserted into a segment or facet during Delaunay renement remains there permanently. However, keep in mind that the edges that partition a facet into subfacets are not permanent, are not treated like subsegments, and are subject to ipping according to the Delaunay criterion. Many approaches to tetrahedral mesh generation permanently triangulate the input facets as a separate step prior to tetrahedralizing the interior of a region. The problem with this approach is that these independent facet triangulations may not be collectively ideal for forming a good tetrahedralization. The algorithm discussed herein uses an alternative approach, wherein facet triangulations are rened in conjunction with the tetrahedralization. The tetrahedralization process is not beholden to poor decisions made earlier. 

(b) Figure 3: (a) Any facet of a PLC may contain holes, slits, and vertices; these may support intersections with other facets and segments or allow a user of the nite element method to apply boundary conditions. (b) When a PLC is tetrahedralized, each facet of the PLC is partitioned into triangular subfacets, which respect the holes, slits, and vertices. This idea generalizes without change to higher dimensions. Unfortunately, my description of Delaunay renement thus far has a gaping hole: mesh boundaries have not been accounted for. The aw in the procedure presented above is that the circumcenter of a poor simplex might not lie in the mesh at all. Delaunay renement algorithms, including those of Chew, Ruppert, and Dey et al., are distinguished primarily by how they handle boundaries.

' S & & $  & 43 & $

3 C

3 C3 $ CC  &  B

The input upon which my three-dimensional Delaunay renement algorithm operates is called a piecewise linear complex (PLC). See Miller, Talmor, Teng, Walkington, and Wang [10] for a denition in  . In three dimensions, a PLC is a set of vertices, segments, and facets. A segment is a constraining edge that must be represented by a sequence of contiguous edges in the nal mesh. A facet is a constraining planar surface that must be represented by a set of triangular faces in the nal mesh. A facet can be quite complicated in shape; it is a polygon (not necessarily convex), possibly augmented by holes, slits, and vertices in its interior. Figure 3(a) demonstrates some of the possibilities. A piecewise linear complex  is required to have the following properties. First,  contains both endpoints of each segment of  . Similarly, facets of  are segment-bounded: for any facet in  , every edge and vertex of the facet must be a segment or vertex of  . Second,  is closed under intersection. For example, if two facets of  intersect at a line segment, that line segment must be a segment of  . Third, if a segment of  intersects a facet of  at a point in the segments interior, then the segment must be entirely contained in the facet. The process of tetrahedral mesh generation necessarily divides each segment into smaller edges called subsegments. (Each seg-

The following sections use the notion of the orthogonal projection of a geometric entity onto a line or plane. Given a facet or subfacet and a point ! , the orthogonal projection proj " #! of ! onto is the point that is coplanar with and lies in the line that passes through ! orthogonally to , as illustrated in Figure 4. The projection exists whether or not it falls in . Similarly, the orthogonal projection proj $ #! of ! onto a segment or subsegment % is the point that is collinear with % and lies in the plane through ! orthogonal to % . Sets of points may be projected as well. If and & are facets, then proj " '& is the set ( proj " #! 0) !12&43 .

`3 A

c3 A

3 A 3A P X X a 0 $ BUb d $ ae $   & 3 S 

In this section, I describe a three-dimensional Delaunay renement algorithm that produces well-graded tetrahedral meshes satisfying any circumradius-to-shortest edge ratio bound greater than two. The algorithm takes a facet-bounded PLC as its input. Tetrahedralized and untetrahedralized regions of space must be separated by facets so that, in the nal mesh, any triangular face not shared by two tetrahedra is a subfacet. The rst step is to form a Delaunay tetrahedralization of the input vertices. Some input segments and facets might be missing (or partly missing) from this mesh. The tetrahedralization is rened by inserting additional vertices into the mesh, using the Bowyer/Watson algorithm to maintain the Delaunay property, until all segments and facets are respected and all constraints on tetrahedron quality are met. Vertex insertion is governed by three rules.

Figure 6: Missing segments are recovered through the same recursive splitting procedure used for encroached subsegments that arent missing. In this sequence of illustrations, the thin line represents a segment missing from the triangulation. Segment recovery is illustrated here in two dimensions for clarity, but operates no differently in three. (a) (b)

croach upon any subsegment, it is not inserted; instead, all the subsegments it would encroach upon are split. A tetrahedron is said to be skinny if its circumradius-to-shortest edge ratio is larger than some bound 9 . Each skinny tetrahedron is normally split by inserting a vertex at its circumcenter, thus eliminating the tetrahedron; see Figure 5(c). However, if the new vertex would encroach upon any subsegment or subfacet, then it is not inserted; instead, all the subsegments it would encroach upon are split. If the skinny tetrahedron is not eliminated as a result, then all the subfacets its circumcenter would encroach upon are split.

(c) Figure 5: Three operations for three-dimensional Delaunay renement. (a) Splitting an encroached subsegment. The original subsegment is encroached because there is a vertex in its diametral sphere. In this example, the two subsegments created by bisecting the original subsegment are not encroached. (b) Splitting an encroached subfacet. The triangular faces shown are subfacets of a larger facet, with tetrahedra (not shown) atop them. In this example, all equatorial spheres (included the two illustrated) are empty after the split. (c) Splitting a skinny tetrahedron. The diametral sphere of a subsegment is the (unique) smallest sphere that encloses the subsegment. Following Ruppert [13], a subsegment is said to be encroached if a vertex other than its endpoints lies inside or on its diametral sphere. A subsegment may be encroached whether or not it actually appears as an edge of the tetrahedralization. It is a property of Delaunay tetrahedralizations that if a subsegment is missing from the tetrahedralization, then it is encroached. Any encroached subsegment that arises is immediately split into two subsegments by inserting a vertex at its midpoint, as illustrated in Figure 5(a). These subsegments may or may not be encroached themselves; splitting continues until no subsegment is encroached.

The equatorial sphere of a triangular subfacet is the (unique) smallest sphere that passes through the three vertices of the subfacet. A subfacet is encroached if a non-coplanar vertex lies inside or on its equatorial sphere. It is a property of Delaunay tetrahedralizations that if a subfacet does not appear in the tetrahedralization, and it is not covered by other faces that share the same equatorial sphere, then it is encroached. (The question of what subfacets should not be missing from the mesh will be considered shortly.) Each encroached subfacet is normally split by inserting a vertex at its circumcenter; see Figure 5(b). However, if the new vertex would en-

Encroached subsegments are given priority over encroached subfacets, which have priority over skinny tetrahedra. These encroachment rules are intended to recover missing segments and facets, and to ensure that all vertex insertions are valid. Although the proof is omitted here, one can show that if there are no encroached subsegments, then each subfacet circumcenter lies in the containing facet; and if there are no encroached subfacets, then each tetrahedron circumcenter lies in the mesh. Because encroached subsegments have priority, the algorithm begins by recovering all the segments that are missing from the initial tetrahedralization. Each missing segment is bisected by inserting a vertex into the mesh at the midpoint of the segment (more accurately, at the midpoint of the place where the segment should be). After the mesh is adjusted to maintain the Delaunay property, the two resulting subsegments may appear in the mesh. If not, the procedure is repeated recursively for each missing subsegment until the original segment is represented by a contiguous linear sequence of edges of the mesh, as illustrated (in two dimensions) in Figure 6. We are assured of eventual success because the Delaunay triangulation always connects a vertex to its nearest neighbor; once the spacing of vertices along a segment is sufciently small, its entire length will be represented. In the engineering literature, this process is sometimes called stitching. When no encroached subsegment remains, missing facets are recovered in an analogous manner. The main complication is that if a facet is missing from the mesh, it is difcult to say what its subfacets are. With segments there is no such problem; if a segment is missing from the mesh, and a vertex is inserted at its midpoint, one knows unambiguously where the two resulting subsegments are. But how may we identify subfacets that do not yet exist? The solution is straightforward. For each facet, it is necessary to maintain a two-dimensional Delaunay triangulation of its vertices, independently from the tetrahedralization in which we hope its subfacets will eventually appear. By comparing the triangles of a facets triangulation against the faces of the tetrahedralization, one can identify subfacets that need help in forcing their way into

Facet

Triangulation

equatorial sphere of f f c F

S s

diametral sphere of s p projF (p)

PLC

Mesh

inside

outside

Figure 8: If a vertex encroaches upon a Delaunay subfacet of a facet , but its projection into the plane containing lies outside , then encroaches upon some subsegment of as well.

Figure 7: The top illustrations depict a rectangular facet and its triangulation. The bottom illustrations depict the facets position as an interior boundary of a PLC, and its progress as it is inserted into the tetrahedralization. Most of the vertices and tetrahedra of the mesh are omitted for clarity. The facet triangulation and the tetrahedralization are maintained separately. Shaded triangular subfacets in the facet triangulation (top center) are missing from the tetrahedralization (bottom center). The bold dashed line (bottom center) represents a tetrahedralization edge that passes through the facet. A missing subfacet is forced into the mesh by inserting a vertex at its circumcenter (top right and bottom right). The vertex is independently inserted into both the triangulation and the tetrahedralization. the mesh. For each triangular subfacet in a facet triangulation, look for a matching face in the tetrahedralization; if the latter is missing, insert a vertex at the circumcenter of the subfacet (subject to rejection if subsegments are encroached), as illustrated in Figure 7. The new vertex is independently inserted into both the facet triangulation and the tetrahedralization. Similarly, the midpoint of an encroached subsegment is independently inserted into the tetrahedralization and into each facet triangulation that contains the subsegment. In essence, the same procedure is used to recover missing segments. However, the process of forming a one-dimensional triangulation is so simple that it passes unnoticed. Which vertices of the tetrahedralization need to be considered in a facet triangulation? If a facet appears in a Delaunay tetrahedralization as a union of faces, then the triangulation of the facet is determined solely by the vertices of the tetrahedralization that lie in the plane of the facet. If a vertex lies near a facet, but is not coplanar with the facet, it may cause a subfacet to be missing (as in Figure 7, bottom center), but it cannot otherwise affect the shape of the triangulation. If a facet appears as a union of faces in a Delaunay tetrahedralization, then those faces form a two-dimensional Delaunay triangulation of the facet. Furthermore, because each facet is segment-bounded, and segments are recovered (in the tetrahedralization) before facets, each facet triangulation can safely ignore vertices that lie outside the facet (coplanar though they may be). The requirements set forth in Section 4 ensure that all of the vertices and segments of a facet must be explicitly identied in the input PLC. The only additional vertices to be considered are those that were inserted in segments to help recover those segments and other facets. The algorithm maintains a list of the vertices on each segment, ready to be called upon when a facet triangulation is initially formed.

Unfortunately, if a facets Delaunay triangulation is not unique because of cocircularity degeneracies, then the facet might be represented in the tetrahedralization by faces that do not match the independent facet triangulation. An implementation must detect these cases and correct the triangulation so that it matches the tetrahedralization. When no encroached subsegment or subfacet remains, every input segment and facet is represented by a union of edges or faces of the mesh. The rst time the mesh reaches this state, all external tetrahedra (lying in the convex hull of the input vertices, but outside the region enclosed by the facet-bounded PLC) are removed prior to splitting any skinny tetrahedra. This measure prevents problems that might arise if superuous skinny tetrahedra are split, such as overrenement and failure to terminate because of spurious small angles formed between the PLC and its convex hull. One further amendment to the algorithm is necessary to obtain the best possible bound on the circumradius-to-shortest edge ratios of the tetrahedra. When several encroached subfacets exist, they should not be split in arbitrary order. If a vertex ! encroaches upon a subfacet of a facet , but the projected point proj " #! does not lie in , then splitting is not the best choice. One can show (with the following lemma) that there is some other subfacet of that is encroached upon by ! and contains proj " #! . (The lemma assumes that there are no encroached subsegments in the mesh, as they have priority.) A better bound is achieved if the algorithm splits rst and delays the splitting of indenitely.

3A

3A

Lemma 1 (Projection Lemma) Let be a subfacet of the Delaunay triangulated facet . Suppose that is encroached upon by some vertex ! , but ! does not encroach upon any subsegment of . Then proj " #! lies in the facet , and ! encroaches upon a subfacet of that contains proj " #! .

3A Proof: First, I prove that proj " 3#! A lies in . Suppose for the sake of contradiction that proj " 3#! A lies outside the facet . Let be the centroid of ; clearly lies inside . Because all facets are segment-bounded, the line segment connecting to proj " # 3 ! A must intersect some subsegment in the boundary of . Let be the

3A

plane that contains and is orthogonal to , as illustrated in Figure 8. Because is a Delaunay subfacet of , its circumcircle (in the plane containing ) encloses no vertex of . However, its equatorial sphere may enclose verticesincluding ! and might not appear in the tetrahedralization. It is apparent that ! and proj " #! lie on the same side of , because the projection is dened orthogonally to . Say that a point is inside if it is on the same side of as , and outside if it is on the same side as ! and proj " #! . The circumcircle of cannot enclose the endpoints of , because is Delaunay in . Furthermore, the circumcenter of lies in ; otherwise, a vertex of would encroach upon some boundary segment of . It follows that

3A 3A

equatorial sphere of f c F

E f e

equatorial sphere of g p g

Figure 9: If a vertex encroaches upon a subfacet of a Delaunay triangulated facet , but does not encroach upon any subsegment of , then encroaches upon the subfacet(s) of that contains proj " .

Figure 10: The radius of each disk illustrated is the local feature size of the point at its center. . Thanks to the Projection Lemma, a dihedral angle of at least this rule ensures that no vertex in the interior of one segment or facet can cause the splitting of a subsegment or subfacet in an incident segment or facet, except that a vertex in the interior of a facet may encroach upon a subsegment of that facet. In theory (and more so in practice), the projection condition can be relaxed somewhat, although the analysis becomes too complicated to present here. Rupperts algorithm and the three-dimensional algorithm described here produce nicely graded meshes, in the sense that a small feature in one part of a mesh does not unreasonably reduce the edge lengths at other, distant parts of the mesh. To formalize this idea, Ruppert introduces a function called the local feature size, dened on all points. Given a PLC  , the local feature size lfs #! at any point ! is the radius of the smallest ball centered at ! that intersects two nonincident features of  (where each of the two features might be a vertex, segment, or facet). Figure 10 illustrates the notion in two dimensions (with no facets) by giving examples of such balls for a variety of points. The function lfs  is continuous and has the property that its directional derivatives (where they exist) are bounded in the range  66! , as the following lemma shows. This property sets a lower bound (within a constant factor to be derived shortly) on the rate at which edge lengths grade from small to large as one moves away from a small feature.

the portion of s equatorial sphere outside lies entirely inside or on the diametral sphere of (as the gure demonstrates). Because ! is inside or on the equatorial sphere of , ! also lies inside or on the diametral sphere of , contradicting the assumption that ! encroaches upon no subsegment of . It follows that proj " #! must be contained in some subfacet of . (The containment is not necessarily strict; proj " #! may fall on an edge interior to , and be contained in two subfacets.) To complete the proof of the lemma, I shall show that ! encroaches upon . If @ the result follows immediately, so assume that @ . Again, let be the centroid of . The line segment connecting to proj " #! must intersect some edge of the subfacet , as illustrated in Figure 9. Let be the plane that contains and is orthogonal to . Say that a point is on the -side if it is on the same side of as . Because the triangulation of is Delaunay, the portion of s equatorial sphere on the -side lies entirely inside or on the equatorial sphere of . The point ! lies on the -side or in (because proj " #! lies in ), and ! lies inside or on the equatorial sphere of , so it must also lie inside or on the equatorial sphere of , and hence encroaches upon .

C 

3A

3A

3A

3A

3A

3A

If ! encroaches upon some subfacet of but no subfacet of contains proj " #! , then ! also encroaches upon some boundary subsegment of . Because encroached subsegments have priority, subsegments encroached upon by ! are split until none remains. The Projection Lemma guarantees that any subfacets of that were encroached upon by ! are eliminated in the process. On the other hand, if proj " #! lies in a subfacet of , and no subsegment is encroached, then splitting is a good choice. As a result, several new subfacets will appear, at least one of which contains proj " #! ; if this subfacet is encroached, then it is split as well, and so forth until the subfacet containing proj " #! is not encroached. The Projection Lemma guarantees that any other subfacets of that were encroached upon by ! are eliminated in the process.

3A

3A

Lemma 2 (Ruppert [13]) For any PLC  , and any two points and , lfs # lfs $" &%(' " ' .

3A

3A

&1&a &a Q  3 $ 43 & $ $ ( & &)( '(23 $

The proof of termination for this algorithm is related to that of Rupperts. Rupperts analysis requires that any two incident input seg  . Similarly, the threements be separated by an angle of at least dimensional analysis requires that the input PLC satisfy the following projection condition: if two constraints (each being a segment or facet) in  intersect, then the orthogonal projection of either one onto the other cannot intersect the interior of the other. (Here, interior is dened to exclude all boundaries, including isolated slits and input vertices in the interior of the facet.) For example, if two convex facets intersect along a segment, they must be separated by

With each vertex , associate an insertion radius 32 equal to the length of the shortest edge connected to immediately after is introduced into the tetrahedralization. Consider what this means in each possible circumstance. If is an input vertex, then 2 is the Euclidean distance between and the input vertex nearest . If is a vertex inserted to split a subsegment or subfacet, then 32 is the distance between and the nearest encroaching mesh vertex. If there is no encroaching vertex in the mesh (a vertex was considered for insertion but rejected as encroaching), then 2 is the radius of the diametral/equatorial sphere of the encroached subsegment/subfacet. If is a vertex inserted at the circumcenter of a skinny tetrahedron, then 2 is the circumradius of the tetrahedron. If a vertex is considered for insertion but rejected because of an

Proof: The ball having radius lfs $" centered at " intersects two nonincident features of  . The ball having radius lfs $" )%0' " ' centered at contains the prior ball, and thus also intersects the same two features. Hence, the smallest ball centered at that contains two nonincident features of  has radius no larger than lfs $" 1%(' " ' .

f 3f A f

3 A

"

3 A

3 A

3 A

ff

! f !

encroachment, its insertion radius is dened the same way, even though it is never actually inserted. Each vertex , whether inserted or rejected, has a parent vertex ! , unless is an input vertex. Intuitively, for any non-input vertex , ! is the vertex that is responsible for the insertion of . If is a vertex inserted to split a subsegment or subfacet, then ! is the encroaching vertex, whether that vertex is inserted or rejected. If there are several encroaching vertices, choose the one nearest . If is a vertex inserted (or rejected) at the circumcenis the most recently inserted ter of a skinny tetrahedron, then ! endpoint of the shortest edge of that tetrahedron.

3f A f 3f f f f f A 3f A f f f
@B9

rv > B rp p rp v rv

p rp rv

3f A

rp = 2rv

(a)
p rp f v rv

(b)

Lemma 3 Let be a vertex of the mesh, and let ! lfs , or 2 parent, if one exists. Then either 2

if is the midpoint of an encroached subsegment or the circumcenter of an encroached subfacet.

if

@ ! 3 f A be its ! 3 f A ! ! , where

is the circumcenter of a skinny tetrahedron;

Proof: If is an input vertex, there is another input vertex a distance of 2 from , so lfs # 2 , and the theorem holds. If is inserted in an encroached subsegment or subfacet, and its parent ! is an input vertex or lies on another subsegment or subfacet, then there are two cases to consider. The case in which lies in a boundary segment of a facet that contains ! will be considered shortly. Otherwise, and ! must lie on nonincident features (because  satises the projection condition), so lfs # 2. If is inserted at the circumcenter of a skinny tetrahedron, then its parent ! @ ! is the most recently inserted endpoint of the shortest edge of the tetrahedron; see Figure 11(a). Hence, the length of the shortest edge of the tetrahedron is at least . Because the tetrahedron is skinny, its circumradius-to-shortest edge ratio exceeds 9 , so its circumradius is 2 9 . If is inserted at the midpoint of an encroached subsegment , and its parent ! is a tetrahedron/subfacet circumcenter that was considered for insertion but rejected because it encroaches upon , then ! lies inside or on the diametral sphere of . Because the tetrahedralization/facet triangulation is Delaunay, the circumsphere/circumcircle centered at ! encloses no vertices, and in particular does not enclose the endpoints of . Hence, # 7 2 ; see Figure 11(b) for an example where the relation is equality. If is inserted at the circumcenter of an encroached subfacet , and its parent ! is a tetrahedron circumcenter that was considered for insertion but rejected because it encroaches upon , then ! lies inside or on the equatorial sphere of . Because the tetrahedralization is Delaunay, the circumsphere centered at ! encloses no vertices, including the vertices of . Recall that the algorithm splits only if contains proj #! . Given these facts, one can show that # 7 2 ; see Figure 11(c) for an example where the relation is equality.

f ! f

rp = 2rv

3f A ! f

(c) Figure 11: The relationship between the insertion radii of a child and its parent. (a) When a skinny tetrahedron is split, the childs insertion times larger than that of its parent. (b) When a radius is at least subsegment is encroached upon by a circumcenter, the childs insertion radius may be a factor of smaller than its parents. (c) When a subfacet is encroached upon by the circumcenter of a skinny tetrahedron, the childs insertion radius may be a factor of smaller than its parents.

3f A

3f A !

 

 

!  !

    " #!$       #!$ " #!$

Tetrahedron Circumcenters

! C B!

Subfacet Circumcenters

Subsegment Midpoints

! C 2!

3A

Figure 12: Dataow diagram illustrating the worst-case relation between a vertexs insertion radius and the insertion radii of the children it begets. If no cycle has a product smaller than one, Delaunay renement will terminate. A guarantee of termination alone is not satisfying when a graded mesh is required. What follows is a proof that each edge of the output mesh has length proportional to the local feature sizes of its endpoints. Hence, edge lengths are determined by local considerations; features lying outside the ball that denes the local feature size of a point can only weakly inuence the sizes of edges that contain that point. Lemma 3 was concerned with the relationship between the insertion radii of a child and its parent; the next lemma is concerned 2 and lfs . For any vertex , with the relationship between lfs 2 dene 2 @ lfs . Think of 2 as the one-dimensional density

Lemma 3 limits how quickly the insertion radius can decrease as one traverses a sequence of descendants of a vertex. If vertices with ever-smaller insertion radii cannot be generated, then edges shorter than existing ones cannot be introduced, and Delaunay renement is guaranteed to terminate. Figure 12 expresses this notion as a dataow graph. Mesh vertices are divided into four classes: input vertices (which are omitted from the gure because they cannot contribute to cycles), vertices inserted into segments, vertices inserted into facet interiors, and vertices inserted at circumcenters of tetrahedra. Labeled arrows indicate how a vertex can cause the insertion of a child whose insertion radius is some factor times that of its parent. If the graph has no cycle whose product is less than one, termination is guaranteed. This goal is achieved by choosing 9 to be at least 7 .

')(% &

' (% & )

' 0% & 1

of vertices near when is inserted, weighted by the local feature size. One would like this density to be as small as possible. 2 # 6 for any input vertex, but 2 tends to be larger for a vertex inserted late.

2
$

@ 3'5 9 % C 7 7 A 9

D f

Lemma 4 Let be a vertex with parent ! (following Lemma 3). Then 2 2

Proof: By Lemma 2, lfs By denition, the # lfs #! insertion radius 2 is ' ! ' if ! is a mesh vertex, whereas if ! is rejected, then 2 ' ! ' . Hence, we have

! !

Theorem 5 Suppose the quality bound 9 is strictly larger than 7 , and the input PLC satises the projection condition. Then there exist xed constants 6 , " 6 , and $ 6 such that, for any vertex inserted (or rejected) at the circumcenter of a skinny ; for any vertex inserted (or rejected) at tetrahedron, 2 # " ; and for the circumcenter of an encroached subfacet, 2 # any vertex inserted at the midpoint of an encroached subsegment, 2 # $ . Hence, the insertion radius of every vertex has a lower bound proportional to its local feature size.

!2 !2D # The result follows by dividing both sides by ! 2 .

@

3f A ! f ! f lfs 3 f A #

# 6 % A % 'f ! '. )

3 f A .0 Suppose that
.

Theorem 6 (Ruppert [13]) For any vertex of the output mesh, 2 the distance to its nearest neighbor is at least lfs . Proof: Inequality 3 indicates that 2 # so Theorem 5 shows that lfs

lfs #!

3 A% 2(! % 2 ! 2 %

!2

added after , then ' ' @ 2 lfs 2 , and the theorem holds. If was added after , apply the theorem to , yielding

f0(

By Lemma 2, lfs

Proof: Consider any non-input vertex ! is an input vertex, then @ lfs

# 6 by Lemma 3. Otherwise, assume for the sake of induction that the lemma is true for ! . Hence, (  "  $ 3 . # First, suppose is inserted or considered for insertion at the 9 . circumcenter of a skinny tetrahedron. By Lemma 3, 2 Therefore, by Lemma 4, 2 # 6 % It follows that one can prove that 2 # if is chosen sufciently large that (  "  $ 3 # 6 % (1)

 f 2
f

'10

f % &

with parent !

3 f A . If !

) . D $ # D 3 E  D ToD provide an example, suppose 9 @ 7 . Then 24 @ D G , D  D 2 " @ 65 , and 2 $ @ 77 6 . Hence, the spacing of vertices is at worst about 7 5 times smaller than the local feature size. Note that as 9 approaches 7 , the values of 2 , 2 " , and 2 $ approach innity.
It follows that

'f '

0(

32( A&%(' f0( ' lfs 3 f A , so  f ( ' lfs 3 f A 2 $ ' f0( ' D '0
lfs

'f '

0(

' (% &

$#

% & 2

$ is larger than and " , $ for any vertex . If was

$#

% ) &

f f

!1

lfs

32( A D
$

%2 &

' 2

2
9

2 2 2

'  &  !" $# % ! D

2

Second, suppose is inserted or considered for insertion at the circumcenter of a subfacet . If its parent ! is an input vertex or lies on a segment or facet not incident to the facet containing , then lfs # 2 , so 2 # 6 and the theorem holds. If ! is the circumcenter of a skinny tetrahedron (rejected for insertion because by Lemma 3, so by Lemma 4, it encroaches upon ), 2 2 # 6 % 7 . It follows that one can prove that 2 # " if " is chosen so that

3f A !

2 2

C 2

'10

Third, suppose is inserted at the midpoint of a subsegment . If its parent ! is an input vertex or lies on a segment or facet not incident to the segment containing , then lfs # 2 , and the theorem holds. If ! is the circumcenter of a skinny tetrahedron or encroached subfacet (rejected for insertion because it encroaches upon ), 2 by Lemma 3, so by Lemma 4, 2 # 6 % $ 7 (  " 3 . It follows that one can prove that 2 # if $ is chosen so that

6 %

C 7 2

# 2
"

(2)

3f A !
2

C 2 

'0 2 2

As Figure 13 shows, the algorithm performs much better in practice. In this example, as soon as all encroached subsegments and subfacets have been eliminated (upper left), the largest circumradius-to-shortest edge ratio is already less than 7 6 . The shortest @ , so the spectre of edge lengths edge length is 6 , and lfs 75 times smaller than the local feature size has not materialized. the As the quality bound 9 decreases, the number of tetrahedra  . in nal mesh increases gracefully until drops below With 9 6 9A@ 6  , the algorithm fails to terminate on this example. Figure 14 offers a demonstration of the grading of a tetrahedralization generated by Delaunay renement. A cube has been truncated at one corner, cutting off a portion whose width is onemillionth that of the cube. Although this mesh satises a bound on circumradius-to-shortest edge ratio of 9 @ 6 7 , reasonably good grading is apparent. For this bound there is no theoretical guarantee, but the worst edge is 5 times shorter than the local feature size at one of  its endpoints. If a bound of 9 @ 7 is applied, the worst edge is (rather than 7 5 ) times smaller than the local feature size at one of its endpoints. After the initial Delaunay tetrahedralization has been formed, my implementation runs in constant time per vertex insertionin practice. The theoretical worst-case running time is surely much time per vertex insertion, even after worse, and may be amortization, in pathological cases. I do not know how to devise such a case, however. In practice, the rst few vertex insertions usually banish any large vertex-free regions, after which further vertex insertions tend to be well-localized and thus execute in constant time.

' C 3

D5

D 63

D3

7 328 A

If the quality bound 9 is strictly larger than 7 , Inequalities 1, 2, and 3 are simultaneously satised by choosing

6 %

C7  ( 2

2

" 3

# 2

9  3 @B bca A & @ (X`a 0 $ BbcBd $ ae $

(3)

2 @ 9 % 9 6  % 7 C 7 

2
"

3 6 % C 7 A 9 % C
9 7

7 

The only known theoretical guarantee about Delaunay renements ability to remove sliver tetrahedra is provided by Chew [4], whose bound on tetrahedron quality is too weak to provide any practical assurance. Nonetheless, it is natural to wonder whether Delaunay renement might be effective in practice. If one inserts a vertex at the circumcenter of any tetrahedron with a small dihedral angle, will the algorithm fail to terminate?

143 vertices, 346 tetrahedra.

 63 , & '0) 1 @ 6 D (H  , 9A@B7 D  & ' @ 6 F HED  7  , '0) 1 @ 6 ,

1009 tetrahedra.

9 @ 6 D 7 , F & '0D )  1 @  6 D 7   , & ' ,  6 D F G 5 5 ,6 334 vertices, '0) 1 @@AE

D G(G , & D'0) 1 @B7 5 D 3  ,D FBH63 9A@B7 D  7 , & D'0) 1 @B7 6 D 5  ,ED 3 & ' @ 65 5 G , '0) 1 @ , & ' @ 65(G G , '0) 1 @ 6 G , 307 vertices, 891 tetrahedra. 1761 vertices, 7383 tetrahedra.
9 @ 6

Figure 15: Meshes created by Delaunay renement with bounds on the smallest dihedral angle . and 6 . The mesh on the right was generated with bounds on both tetrahedron volume and dihedral angle, so that enough tetrahedra were generated to ensure that the mesh on the left wasnt  merely a uke. (The best attainable lower bound drops by 7 7 as a result.) Experiments with very large meshes suggest that a mini  can be obtained reliably. See Chew [4] for some mum angle of 6 intimations as to why slivers might be eliminated so readily. One of his useful observations is that the region in which a tetrahedrons apex must lie so that the tetrahedron can simultaneously have poor dihedral angles and a good circumradius-to-shortest edge ratio is relatively small. Unfortunately, my success in removing slivers is probably due in part to the severe restrictions on input angles I have imposed. Practitioners report that they have the most difculty removing slivers at the boundary of a mesh, especially near small angles. Mesh improvement techniques such as optimization-based smoothing and topological transformations can likely remove some of the imperfections that cannot be removed directly by Delaunay renement.

'0) 1

76

5  

5596 tetrahedra.

 F , & '0) 1 @ 6 D C  , 9A@ 6 D ( D  6C6  , &2  '    @  6 DFBF C H '0) 1 @ 5 , 1397 vertices,

13969 tetrahedra.

9 @ 6 D  52F 6 , D & '0 @ D  5  , ) 1   , &2 ' @AE  6 D 6 G  75 , 3144 '0) 1 @ vertices,

Figure 13: Several meshes of a PLC generated with

different bounds ( ) on circumradius-to-shortest edge ratio. Below each mesh is listed the smallest dihedral angle , the largest dihedral angle , and the shortest edge length . The algorithm does not terminate on this PLC for the bound   .

'

'0 '0)) 11

& $ B 0 3 & $

Figure 14: At left, a mesh of a truncated cube. At right, a crosssection through a diagonal of the top face. As Figure 15 demonstrates, Delaunay renement can succeed for useful dihedral angle bounds. Each of the meshes illustrated was generated by enforcing a lower bound on dihedral angles, rather than a circumradius-to-shortest edge ratio bound. However, the implementation prioritizes poor tetrahedra according to their ratios, and thus slivers are split last. I suspect that the program generates meshes with fewer tetrahedra this way, and that the likelihood of termination is greater. Intuitively, one expects that a vertex inserted at the circumcenter of the tetrahedron with the largest ratio is more likely to eliminate more bad tetrahedra. Both meshes illustrated have dihedral angles bounded between

& '0) 1

Delaunay renement is an effective technique for three-dimensional mesh generation. Its theoretical guarantees on tetrahedron quality and grading make it attractive. However, these guarantees are not entirely satisfying. There is no guarantee that slivers can be eliminated. The bounds on edge length are reassuring, but might not be strong enough for all practical purposes, as the number of tetrahedra is inversely proportional to the cube of the edge length. Although Ruppert uses Theorem 6 to prove the size-optimality of the triangular meshes his algorithm produces, an analogous result is not possible in three dimensions. Fortunately, Delaunay renement algorithms outperform their worst-case bounds. There is a great deal of slack in the analysis; the relationship between the insertion radii of a child and its parent is usually looser than in the worst case. This slack accumulates as one traces a sequence of descendants from an input vertex, and makes it possible to achieve better bounds on tetrahedron quality and edge length than the theory promises. (Miller et al. [10] cannot exploit this slack, and thus require about 8,000 times as many vertices to guarantee a bound of 9A@B7 on the PLC in Figure 13). The main outstanding problem in three-dimensional Delaunay renement is the question of how best to handle small input angles. The difculty is that if vertices lying in segments or facets can encroach upon incident subsegments and subfacets, then additional arcs appear in the dataow diagram of Figure 12, creating cycles of ever-diminishing insertion radii. Fortunately, it is possible to modify the algorithm so that it always terminates, even when the

D3

Figure 16:

A counterexample demonstrating that the threedimensional Delaunay renement algorithm is not size-optimal.

projection condition is not satised. The insertion radius of each vertex can be explicitly computed. A vertex should be rejected if its insertion radius is too small compared to its parents. As a result, some skinny tetrahedra may remain in the nal mesh, but only in the vicinity of small input angles. A few other tricks are needed to ensure that all segments and facets are recovered: segments are recovered by using a technique based on concentric spherical shells suggested by Ruppert [13], then facets are recovered by forming a constrained Delaunay tetrahedralization [15]. The reason why a size-optimality result cannot be established is worth exploring. Gary Miller and Dafna Talmor (private communication) have pointed out the counterexample depicted in Figure 16. Inside this PLC, two segments pass very close to each other without intersecting. The PLC might reasonably be tetrahedralized with a few dozen tetrahedra having bounded circumradius-toshortest edge ratios, if these tetrahedra include a sliver tetrahedron whose four vertices are the endpoints of the two internal segments. However, the best the present Delaunay renement algorithm can promise is to ll the region with tetrahedra whose edge lengths are proportional to the distance between the two segments. Because this distance may be arbitrarily small, the algorithm is not size-optimal. However, if guaranteed bounds could be established for the dihedral angles, and not merely the circumradius-to-shortest edge ratios, then size-optimality might be proven using ideas like those with which Mitchell and Vavasis [11, 12] demonstrate the size-optimality of their octree-based algorithms. Many theoretical improvements can be made to the core algorithm described here. The bound on circumradius-to-shortest edge ratio can be improved to 6 5 by replacing equatorial spheres with narrower protective shells called equatorial lenses. (This technique also requires the use of constrained Delaunay tetrahedralizations.) The bound can be further improved to 6 6 with a technique called range-restricted segment splitting, albeit by sacricing good grading in theory (not in practice). With or without these improvements, a bound as small as 6 can be applied to any tetrahedron that does not intersect the interior of a segment or facet, and thus only boundary tetrahedra are subject to the weaker bound of 7 , 6 5 , or 6 6 . The projection condition can be somewhat weakened; for examC  separation between incident segments sufces. All these ple, a improvements are described in detail elsewhere [14]. The improvements in circumradius-to-shortest edge ratio, in particular, are signicant and further propel the approach described here beyond the pioneering tetrahedral mesh generation algorithms of Dey et al., Chew, and Miller et al.

DH

D3

DH

D3

[2] L. Paul Chew. Guaranteed-Quality Triangular Meshes. Technical Report TR-89-983, Department of Computer Science, Cornell University, 1989. [3] . Guaranteed-Quality Mesh Generation for Curved Surfaces. Proceedings of the Ninth Annual Symposium on Computational Geometry, pages 274280. ACM, May 1993. [4] . Guaranteed-Quality Delaunay Meshing in 3D. Proceedings of the Thirteenth Annual Symposium on Computational Geometry, pages 391393. ACM, June 1997. [5] Tamal Krishna Dey, Chanderjit L. Bajaj, and Kokichi Sugihara. On Good Triangulations in Three Dimensions. International Journal of Computational Geometry & Applications 2(1):7595, 1992. [6] Lori A. Freitag and Carl Ollivier-Gooch. A Comparison of Tetrahedral Mesh Improvement Techniques. Fifth International Meshing Roundtable (Pittsburgh, Pennsylvania), pages 87100. Sandia National Laboratories, October 1996. [7] William H. Frey. Selective Renement: A New Strategy for Automatic Node Placement in Graded Triangular Meshes. International Journal for Numerical Methods in Engineering 24(11):21832200, November 1987. [8] Charles L. Lawson. Software for Surface Interpolation. Mathematical Software III (John R. Rice, editor), pages 161 194. Academic Press, New York, 1977. [9] Gary L. Miller, Dafna Talmor, Shang-Hua Teng, and Noel Walkington. A Delaunay Based Numerical Method for Three Dimensions: Generation, Formulation, and Partition. Proceedings of the Twenty-Seventh Annual ACM Symposium on the Theory of Computing, pages 683692, May 1995. [10] Gary L. Miller, Dafna Talmor, Shang-Hua Teng, Noel Walkington, and Han Wang. Control Volume Meshes using Sphere Packing: Generation, Renement and Coarsening. Fifth International Meshing Roundtable, pages 4761, October 1996. [11] Scott A. Mitchell and Stephen A. Vavasis. Quality Mesh Generation in Three Dimensions. Proceedings of the Eighth Annual Symposium on Computational Geometry, pages 212 221, 1992. [12] . Quality Mesh Generation in Higher Dimensions. Submitted to SIAM Journal on Computing, February 1997. [13] Jim Ruppert. A Delaunay Renement Algorithm for Quality 2-Dimensional Mesh Generation. Journal of Algorithms 18(3):548585, May 1995. [14] Jonathan Richard Shewchuk. Delaunay Renement Mesh Generation. Ph.D. thesis, School of Computer Science, Carnegie Mellon University, Pittsburgh, Pennsylvania, May 1997. Available as Technical Report CMU-CS-97-137. . A Condition Guaranteeing the Existence of Higher[15] Dimensional Constrained Delaunay Triangulations. This proceedings, June 1998. [16] David F. Watson. Computing the -dimensional Delaunay Tessellation with Application to Voronoi Polytopes. Computer Journal 24(2):167172, 1981. [17] Nigel P. Weatherill. Delaunay Triangulation in Computational Fluid Dynamics. Computers and Mathematics with Applications 24(5/6):129150, September 1992.

bcB $ 2C

[1] Adrian Bowyer. Computing Dirichlet Tessellations. Computer Journal 24(2):162166, 1981.