BIB-VERSION:: CS-TR-v2.0
ID:: CORNELLCS//TR93-1327
ENTRY:: 1993-10-14
ORGANIZATION:: Cornell University, Computer Science Department
LANGUAGE:: English
TITLE:: Mesh Generation With Provable Quality Bounds
AUTHOR:: Mitchell, Scott A.
DATE:: February 1993
PAGES:: 161
COPYRIGHT:: Scott A. Mitchell - ALL RIGHTS RESERVED
ABSTRACT::
We consider the problem of generating a triangulation of provable quality for 
two and three dimensional polyhedral regions. That is, we seek a 
triangulation, allowing additional vertices called Steiner points, such that 
the triangular or tetrahedral elements have a bound on their shape. In three 
dimensions we also seek an upper bound on the number of tetrahedra in the 
triangulation. These triangulation algorithms find application in mesh 
generation for finite element methods. The polyhedral region must be bounded 
and well defined, but may have holes of arbitrary complexity.  

In three dimensions, we assume there are no restrictions on where we may place 
Steiner points. Our triangulation is optimal in the following two senses. 
First, our triangulation achieves the best possible aspect ratio up to a 
constant factor, which is our bound on element shape. Second, for any other 
triangulation of the same region into $m$ triangles with bounded aspect ratio, 
our triangulation has size $n = O(m)$. Such a triangulation is desired as an 
initial mesh for a finite element mesh refinement algorithm. Previous three 
dimensional triangulation schemes either worked only on a restricted class of 
input, or did not guarantee well-shaped tetrahedra, or were not able to bound 
the output size. We build on some of the ideas presented in previous work by 
Bern, Eppstein and Gilbert, who have shown how to triangulate a two 
dimensional polyhedral region with holes, with similar quality and optimality 
bounds.

In two dimensions, we assume the restriction that we may introduce Steiner 
points on the polygon's interior, but not on its boundary. Of all 
triangulations satisfying this restriction, our triangulation has the maximum 
minimum angle, up to a constant factor. This is the first known algorithm for 
this problem with provably optimal element shape. The algorithm solves several 
subproblems of mesh generation.
END:: CORNELLCS//TR93-1327
BODY::
Mesh Generation
With Provable Quality Bounds
Scott A. Mitchell
Ph.D Thesis
93-1327
February 1993
Department of Computer Science
Cornell University
Ithaca, NY 14853-7501
MESH GENERATION WITH PROVABLE QUALITY
BOUNDS
A Dissertation
Presented to the Faculty of the Graduate School
of Cornell University
in Partial Fulfillment of the Requirements for the Degree of
Doctor of Philosophy
by
Scott A. ?4itchell
January 1993
Oc Scott A. Mitchell 1993
ALL RIGHTS RESERVED
Mesh Ceneration with Provable Quality Bounds
Scott A. Mitchell, Ph.D.
Cornell University 1993
We consider the problem of generating a triangulation of provable quality for two
and three dimensional polyhedral regions. That is, we seek a triangulation, allow-
ing additional vertices called Steiner points, such that the triangular or tetrahedral
elements have a bound on their shape. In three dimensions we also seek an up-
per bound on the number of tetrahedra in the triangulation. These triangulation
algorithms find application in mesh generation for finite element methods. The
polyhedral region must be bounded and well defined, but may have holes of arbi-
trary complexity.
In three dimensions, we assume there are no restrictions on where we may place
Steiner points. Our triangulation is optimal in the following two senses. First,
our triangulation achieves the best possible aspect ratio up to a constant factor,
which is our bound on element shape. Second, for any other triangulation of the
same region into m triangles with bounded aspect ratio, our triangulation has size
n = 0(m). Such a triangulation is desired as an initial mesh for a finite element
mesh refinement algorithm. Previous three dimensional triangulation schemes ei-
ther worked only on a restricted class of input, or did not guarantee well-shaped
tetrahedra, or were not able to bound the output size. We build on some of the
ideas presented in previous work by Bern, Eppstein, and Gilbert, who have shown
how to triangulate a two dimensional polyhedral region with holes, with similar
quality and optimality bounds.
In two dimensions, we assume the restriction that we may introduce Steiner
points on the polygon's interior, but not on its boundary. Of all triangulations
satisfying this restriction, our triangulation has the maximum minimum angle,
up to a constant factor. This is the first known algorithm for this problem with
provably optimal element shape. The algorithm solves several subproblems of mesh
generation.
Biographical Sketch
Scott Alan Mitchell was born on April 1, 1966 in Madison, Wisconsin. He enjoyed
a happy childhood full of bicycling the backroads, climbing the Sierra mountains,
canoeing the Wisconsin River, traveling to Australia, and walking with his family.
He married his high school sweetheart, Carissa, an aspiring graphic designer. She
dealt with the rigors of a spouse in graduate school gracefully.
After receiving his bachelor's degree in Applied Math, Engineering, and Physics
from the University of Wisconsin-Madison in May, 1988, he packed up his wife and
his bicycles and moved to Ithaca New York. There he had hopes of bicycling up the
hills of the finger lakes region and of completing his doctorate in Applied Math at
Cornell University. While the former goal was met almost immediately, the latter
was realized in January 1993.
Along Scott's journey through graduate school, he and his wife were joined by a
very special little boy, Evan, with whom Scott is carrying on the Mitchell tradition
of long walks and blank stares.
111
For my wife, Carissa.
iv
Acknowledgements
I thank Steve Vavasis for his wonderful guidance during my stay at Cornell. I truly
appreciate his willingness to make time for me and help with even the most droll
of subproblems. His honesty has made him a pleasure to work with.
I thank Joe Mitchell of Cornell and SUNY Stony Brook for sparking my interest
in computational geometry through his well planned courses. I also must thank
him for introducing me to Steve Vavasis. Joe's weekly work sessions were always
stimulating. I thank those who participated in them over the years, Paul Chew,
Christine Piatko, Robert Freimer, Klara Kedem, Erik Wynters, Estie Arkin, Bob
Connelly and Bob's many visitors. I also thank Robert Freimer for continuing the
sessions in Joe's absence.
I thank my fellow Applied Math students at Cornell. Their support through
the rigorous boot camp of analysis and topology courses was invaluable. I thank
Dolores Pendell for sparing me from many administrative hassles. I appreciate the
picnics in the parks and parties at John Guckenheimer's. I will miss our dilapidated
but homey offices in Sage Hall. I would especially like to thank classmates Rick
Wicklin, Keith Saints, Patrick Worfolk, Jianguo Liu and Greg Henry, as well as
Russ Qualls (civil engineering) for making Cornell bearable through the difficult
times.
v
I thank everyone at the Xerox PARC CSL for a productive yet restful summer
1991 in California. Thanks to Marshall Bern for his continued interest in my work.
Concerning the technical development of this thesis, I thank Steve Vavasis for
his indispensable involvement with the three dimensional algorithm. I thank Paul
Chew for discussing various triangulations with me, and for his insights concerning
the development of the two dimensional algorithm. Thanks to Marshall Bern and
John Gilbert of Xerox PARC and David Eppstein of UC Irvine for help with their
earlier work. Thanks to Joe Mitd?ell of SUNY Stony Brook for helping with the
plane sweep algorithm.
vi
Table of Contents
2
Introduction
1.1 The three dimensional algorithm
1.2 The two dimensional algorithm
1.3 Thesis overview
1.4 Some definitions.
Other work
2.1 Computational geometry literature
2.2 Geometric modeling.
2.3 Engineering literature.
2.3.1 Buratynsici [1990J
3 Element shape
4
Quality mesh generation in three dimensions
4.1 Subdivision of P into cubes, the octree
4.1.1			Data structures and definitions			.
4.1.2			Creating the octree . . .
4.1.3			The vertex phase .
4.1.4			The edge phase . . .			. .
4.1.5			The facet phase
4.1.6			The duplication process .
4.1.7			Summary of the octree algorithm
4.2			Relative sizes of adjacent octree boxes
4.2.1			Vertex cone boxes
4.2.2			Edge cone boxes
4.2.3			Facet Cone Boxes
4.3			Warping
4.4			Two dimensional triangulation			.			. .
4.5			Three			dimensional			triangulation
4.5.1			Empty box tetrahedra			. .			.
4.5.2			Vertex box tetrahedra
4.5.3			Edge box tetrahedra .			. .			.
2
6
10
10
12
12
16
16
18
20
26
27
27
29
30
34
36
38
43
45
46
55
56
56
61
66
67
67
67
vii
4.5.4 Facet box tetrahedra
4.5.5 Two-facet box tetrahedra
Aspect ratio of tetrahedra in AocT
Aspect ratio tetrahedra in A*
Optimality of the cardinality of AocT
Conclusions
5
4.6
4.7
4.8
4.9			. . . . . .			. .
Maxmin angle triangulation with no boundary Steiner points
5.0.1			Two dimensional algorithm overview . 			.			. .			.			.
5.1			Constructing wedges around vertices			.			. .			.			.
5.1.1			Constructing a complete wedge around a vertex .			.			.
5.2			Shrinking the polygon Q to R . 			.			. .			.
5.3			Triangulation of 1?. . 			.			. .			.			.
5.4			Matching R to Q . . 			.			.			.
5.4.1			Triangulating a match trapezoid.			. .			.			.
5.4.2			Triangulating a match triangle
5.5			Conclusions.			. .			.			.
5.5.1			Why the two dimensional analog of Chapter			4?			. .
5.5.2			Running time . . . . . . . . 			. .			. .			. .
5.5.3			Cardinality
5.5.4			Open problems . .
67
69
69
79
82
98
99
100
101
125
133
138
139
140
154
155
155
156
157
157
viii
List of Figures
1.1
1.2
3.1
3.2
3.3
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13
4.14
4.15
4.16
4.17
The optimal min angle triangulation of a polygon without Steiner
points (left) and with interior Steiner points (right)
The optimal min angle may be much better than linear in the ratio
of the length of an edge to its closest interfering point. .
The sides and angles of a general triangle
The radius of the inscribed circle. .
The radius of the circumcircle. . . . . .
8
8
23
24
24
Here a box is duplicated for two components, and one duplicate box
is split because it is crowded
No matter how small we split a box U, if it is outside of ex(b) then
the balance condition will not split b. . . . . . . .
Uniting boxes around a vertex box in two dimensions.
Floaters and tethers for a two dimensional polygon
Determining the components of a box. Sweeping the boundary de-
termines two components, and following tethers determines that all
internal floaters belong to the same component, and there are no
internal components
The distance between two superfacets of an edge inside an extended
vertex box I(ex(B)) but outside an edge box is at least kch(R).
The special case of sharing a vertex with two			boxes.
One case in the proof of Lemma 6 . . 			. .
Another case in the proof of Lemma 6
The six cases for triangulating a box facet			. .
Triangulating a quadrilateral in Case C2			(closeup)			. .
The flat (t3) and curved (u3) tetrahedra			. .
A sphere tangent to three facets of a tetrahedron			.
Constructing a large inscribed sphere of U3
Finding the tetrahedron with an angle at most ?			.
An illustration of the conditions stated in Lemma 20.. . . .
An example path in the proof of Lemma 23 . Here all P vertices
represent edges, and edges represent facets, except for the P vertex
J = z. Note b3 is crowded and b1 is split by the balance condition.
31
33
34
39
42
52
57
60
61
64
65
72
74
76
81
86
92
ix
4.18 An illustration of an instance of Case 2 of the proof of Lemma 23
leading to the case in the proof of Lemma 26 where F' and G' have
a common vertex v.			. . .			. . .			. . .94
5.1			The triangles Vi, spokes Si, and rim edges Ri of a wedge at V. . 			101
5.2			The wedge at V in a given triangulation where the sector of the last
triangle Vk = AYVX contains the interfering vertex W111
5.3			The first subcase of where an interfering vertex W may lie relative
to the spiral triangle V'112
5.4			The second subcase of where an interfering vertex W may lie relative
to the spiral triangle V'			113
5.5			The continuous spiral for an optimal wedge114
5.6			The interior side of E relative to V is the same as the interior side
of E1 with respect to V1124
5.7			The interior side of E relative to V is opposite the interior side of
E1 with respect to V1. . . 			125
5.8			Two smallest spirals for the same vertex V meet at X . . 			126
5.9			Welding the two smallest optimal spirals for V.. . . 			128
5.10			Replacing VX by the two spokes VX' and VX". This is the final
transformation for welding the two smallest optimal spirals at V			129
5.11			Welding the two wedges for one P edge and its two vertices V and			U.129
5.12			The interior angle of Q where we welded two smallest spirals for one
P vertex and two different P edges is at least r/2, regardless of the
interior angle at V. . 131
5.13			The smallest interior angle 0 of Q may occur where we welded the
two smallest spirals for one P edge131
5.14			Shrinking Q to R near a vertex V with interior angle < 3r/4 (left)
and interior angle > 37r/4 (right). . . . 			134
5.15			Bounding ?, the angle between VW and E. . . . . . . . 			135
5.16			Bounding the lengths of shrunken edges F1 and F2 meeting at an
angle < 37r/4. . . . . . . 138
5.17			Triangulating a special case trapezoid with layers of boxes. . 			140
5.18			Drawing squares as a preliminary step to triangulating a general
match trapezoid. . . . . . . . . . . . . . . . . 			143
5.19			Triangulating the squares with the help of a center vertex. . . 			148
5.20			How to triangulate a square drawn around three edges149
5.21			There may be many (1/A) square vertices in a match edge for con-
stant shaped squares, or at most a constant number of square ver-
tices for wide and short (1/A) squares. . . 			150
5.22			The worst case of a triangle containing a vertex of E. . . . 			152
x
Chapter 1
Introduction
Triangulation of polyhedral regions is a fundamental geometric problem for flu-
merical analysis. In particular, if one wishes to solve an elliptic boundary value
problem on a three-dimensional domain, it is necessary to discretize the domain
with a mesh. If the domain is sufficiently complicated, then the method of finite
elements is commonly used. In this method the domain is divided into small con-
vex polyhedral regions with a fixed number of faces called elements. A common
choice for an element in three dimensions is the tetrahedron, and in two dimen-
sions the triangle. Thus, a triangulation of the domain is required. Triangulations
are also desired for other applications as well, such as solid modeling, functional
interpolation, and computer graphics.
For numerical stability in the finite element method, it is necessary that the
elements have some bound on their shape. For most applications it is agreed upon
that the closer an element is to equilateral the better. However, which measure
of an element is best is unclear, even restricted to finite element mesh generation.
For information about element shape in numerical analysis, see Babuska and Aziz
2
(1976j. We have chosen to consider the minimum interior angle for two dimensions.
That is, we seek a triangulation whose minimum angle is as large as possible. For
three dimensions, we consider what we call aspect ratio, the ratio of the radius of
the smallest sphere containing the tetrahedron, to the radius of the larges sphere
contained in the tetrahedron. That is, we seek a triangulation whose aspect ratio
is as small as possible. These two measures are very related, as will be described
in Chapter 3.
One question that often arises when designing geometric software is how the
polyhedral regions are represented. One may define such a region as the regularized
intersection and union of half spaces. This is known as the "computational solid
geometry," or CSG, representation. Also, one may define the region in terms of
the lattice of boundary faces, and which side of the d--H 1 dimensional faces is the
interior of the region. This is known as the "B-rep". The composition of half
spaces is mathematically easier to describe. However, it has several algorithmic
drawbacks, most notably having to traverse a tree of possibly large depth to identify
specific faces of the boundary. We have chosen to use the boundary face lattice
representation for our algorithms.
This thesis presents two triangulation algorithms, one for three dimensions and
one for two dimensions. Both algorithms generate a triangulation for a nonconvex
bounded polyhedral domain with holes.
1.1 The three dimensional algorithm
For the three dimensional algorithm, we allow ourselves the freedom to introduce
additional vertices, called Steiner points, in the region's interior, and also on its
3
boundary. The primary goal of the algorithm is to achieve good aspect ratio.
However, care must be taken not to introduce too many faces, as the number
of tetrahedra in the triangulation is also important. For example, the running
time of a finite element method grows with the number of elements. Some very
complicated and large objects, such as an entire airplane, require hundreds of
thousands of elements and mere disk space can become a constraint.
We are able to achieve both of our goals. That is, the triangulation is optimal
in two respects. First, the best possible aspect ratio is achieved for the tetrahedra,
up to a constant factor. Second, the number of tetrahedra is within a constant
factor of the best possible for any triangulation with bounded aspect ratios. This
second constant factor depends on the smallest interior angle of the input region.
Other authors have considered three-dimensional regions, but, to our knowl-
edge, no previous work has simultaneously addressed the problems of optimality,
aspect ratio and of complicated nonconvex regions. Indeed, as far as we know,
there is no previous algorithm to triangulate a nonconvex three-dimensional region
with guaranteed aspect ratio (regardless of the optimality of the triangulation).
Our work is closely based on earlier work by Bern, Eppstein and Gilbert [1990],
who solved the corresponding problem for two-dimensional polyhedral domains.
The main complication in extending the work to three dimensions concerns vertices
of the domain. In particular, for a two-dimensional domain a vertex is adjacent
to only two edges, and there are essentially three cases for vertices (the two edges
make an acute angle, an obtuse angle, or an angle greater than 18o?). In three
dimensions, however, a large number of faces can meet at a single vertex, forming
a complicated cross-section with some convex and some concave angles. Handling
4
vertices properly necessitates a host of changes in the algorithm. The key idea is
that we isolate the vertex in a cube. Then, we triangulate within the cube by first
triangulating in two dimensions the boundary of the cube, and then taking the
convex hull of these triangles and the vertex to form tetrahedra.
Our running time bounds are probably not the best possible; this is discussed
further at the end of Chapter 4. Our optimality condition is slightly stronger
than Bern et. al.'s condition. In particular, they show that overall the number of
triangles they generate is no more than a constant above optimal. The triangles
used in a particular region of the domain, however, could be arbitrarily smaller
than what is needed for the triangulation to achieve good aspect ratio. In our
triangulation, no tetrahedron is ever more than a constant factor smaller than
the tetrahedron at the same location in the best possible triangulation. This is
achieved by using cube sizes that depend on the geodesic distance between disjoint
faces of the input. We are able to use geodesic distance by duplicating boxes when
the intersection of a box with the input forms several connected components. Our
optimality proof is based on characteristic length functions. Similar functions have
been used before by Miller and Thurston [1990], Miller and Vavasis [1990], and
Miller, Teng and Vavasis [1991]. We have not seen them used before, however,
to analyze the construction of a triangulation. We believe the proof method is
particularly powerful and has great potential. For example, it has been adopted
by Rupert[1992], whose two dimensional algorithm does not even use quadtrees.
Our technique is based on an octree partitioning of the domain. That is, we
use a three dimensional grid of cubes (boxes) to divide up to domain and isolate
disjoint faces of the boundary of the region. The boxes containing the boundary
5
must be deformed away from cubes, and all boxes are triangulated. This technique
has been used before, for example, by Carey, Sharna and Wang [1988].
Our algorithm generates an unstructured mesh. That is, the sizes of the cubes
is allowed to vary. Structured meshes, where the size of the elements are fixed,
have several advantages. See Bern and Eppstein [1992]. However, they are ill
suited to input regions that have both very small and very large faces, or where
the equation being solved by a finite element method has a greatly varying gradient.
In both cases, to get the required resolution one must have an inordinate number
of elements. We have chosen to consider problems where it is advantageous to vary
the element size.
The sizes of the cube, and hence the tetrahedra, is completely determined by
the geometry of the input. For some applications, one may wish there to be smaller
elements in a portion of the region. For example, for using finite element methods
to analyze the temperature of an object, it may be known that a certain point
contains a heat source. Since the gradient of temperature will vary greatly near
that point, to achieve an accurate solution small elements are needed near that
source. All one needs to do to insure smaller elements near that source is to
specify that the boxes may be no larger than a particular size near the source.
That is, the user may define a function that takes as an argument a point of the
input region, and returns the size of the largest desired box containing that point.
Our algorithm could refer to this function to easily comply with the user's wishes.
Section 4.1.7 summarizes the algorithm.
6
1.2 The two dimensional algorithm
The problem of triangulating a two dimensional polygon with holes with minimum
angle within a constant factor of optimal and simultaneous cardinality bounds has
been solved by by Bern, Eppstein and Gilbert [1990]. When no Steiner points
are allowed, the triangulation that maximizes the minimum angle is the Delau-
nay triangulation. The numerous algorithms for the Delaunay triangulation are
summarized in Bern and Eppstein [1991].
The two dimensional problem we consider strikes a middle ground between
arbitrary Steiner points and no Steiner points. That is, we seek to triangulate
a polygon with minimum angle within a constant factor of optimal, where we
are allowed Steiner points on the polygon's interior, but not its boundary. This
algorithm may find application in a number of mesh generation subproblems.
This algorithm allows us to triangulate intersecting regions independently. One
application is that it may be used to allow standard mesh generation algorithms to
accept degenerate input. Degenerate faces such as edges that have the interior of
the input region on both sides are often desired to "probe" a numeric solution at
specific locations. Such an edge is also useful on the boundary between two differ-
ently doped regions in the analysis of semiconductors. We could take a degenerate
edge, and fatten it in one direction into a rectangular hole of finite width. Then
any mesh generation algorithm could be used to triangulate the region outside the
hole. Our no boundary Steiner point algorithm can then be used to triangulate
the hole, respecting the mesh already generated.
Consider another related example. In our three dimensional algorithm, we
triangulate each two dimensional facet of the intersection of a cube and the input
7
region before introducing tetrahedra. It does not matter if in this step we introduce
Steiner vertices in the interior of the facet. However, if we introduced Steiner
points on an edge of the boundary of the facet, that edge is shared by facets of
other cubes. We may have already triangulated that facet, in which case we have
ruined its triangulation by introducing a vertex in the middle of an edge of another
triangle. Thus a triangulation that introduces Steiner points only on the regions
interior is desired.
An important feature of the algorithm that makes it useful is that we give
explicit bounds on the minimum angle of the triangulation in terms of geometric
features of the input, namely edge lengths and certain angles. Consider again the
example in the previous paragraph of triangulating the two dimensional intersection
of the cubes of our three dimensional algorithm. If we could give explicit bounds on
the geometric features of the intersection of our cubes with the input region, then
we would not have to resort to the tedious case by case triangulation presented in
Section 4.4. We do not take this approach because we need a triangulation that
maximizes the minimum height instead.
The algorithm may also find application in mesh refinement. That is, suppose
a mesh exists for a given polygon, but after running a finite element method we
discover that the mesh is too coarse near a heat source. We may then "erase" the
mesh in a neighborhood of the heat source. Then a finer mesh could be generated
strictly inside this neighborhood, and our algorithm could be used to triangulate
the gap between the outer and inner mesh. Care must be taken in selecting the
inner mesh, as in general it is not possible to generate an arbitrarily fine inner mesh
while retaining the outer mesh and still have minimum angle within a constant
8
Figure 1.1: The optimal min angle triangulation of a polygon without Steiner
points (left) and with interior Steiner points (right).
Figure 1.2: The optimal min angle may be much better than linear in the ratio of
the length of an edge to its closest interfering point.
factor of optimal. The advantage for this method is that a mesh for the entire
region does not need to be recalculated.
The algorithms that do not add any Steiner points are inadequate for these
tasks. The four point example of Figure 1.1 shows that for a polygon the optimal
min angle when Steiner points are not allowed may be much worse than the optimal
min angle when Steiner points are allowed on the interior.
A key question is what features of a general polygon P determine its optimal
min angle A when Steiner points are allowed on the interior of P, but not its
boundary. Bern Eppstein Gilbert [1990] considers the problem where Steiner points
9
are allowed everywhere in P. There the distance between an edge and a point on
a closed face disjoint from the edge (an interfering point) determines how small
triangles must be in order to have aspect ratio within a constant factor of optimal.
For determining A, because we are not allowed to add Steiner points on ap, it may
be impossible to make triangles this small. For example, in Figure 1.1, there will
always be a triangle with an edge that is the bottom long edge of P. Indeed, as
we vary the distance between the top and bottom edges of Figure 1.1, the optimal
aspect ratio changes proportionally (up to a constant factor).
Hence the optimal min angle for this problem should depend on the maximum
ratio of the length of an edge E to its distance to an interfering point. However,
the optimal min angle is actually more closely dependent on the angle that the
interfering point makes with the closer vertex V of E. That is, if this angle is large,
it may be possible to triangulate with min angle much more than the interfering
point distance to edge length ratio.
For an intuitive reason for why this is so, consider Figure 1.2. The closest
interfering point to edge E is the vertex of F opposite V. If we vary the length
of edge F, the ratio of the length of E to the closest interfering point changes
proportionally. Suppose we also keep the same combinatorial triangulation as we
vary F, but allow ourselves to change the position of the Steiner points. As there
are ten similar triangles with vertex V, if we halve the length of F, we may also
change the lengths of the edges containing V, so that the resulting triangulation
has min angle much more than half of that of the original triangulation. In fact
for this example we could construct a triangulation with min angle about 2--H1/10
times that of the original triangulation.
10
1.3 Thesis overview
In Chapter 2 we briefly review some relevant literature. In Chapter 3 we give more
details on aspect ratio and minimum angle.
In Chapter 4 we present the three dimensional triangulation algorithm. In
Section 4.1 we discuss the formation of the octree. In Section 4.2 we demonstrate
a relation between the boxes of the octree. In Section 4.3 to Section 4.5 we discuss
the steps necessary to produce a triangulation given an octree. In Section 4.6 to
Section 4.8 we provide the optimality proof. We conclude the three dimensional
algorithm in Section 4.9.
In Chapter 5 we present the two dimensional triangulation algorithm. It has
three main steps. First in Section 5.1 we introduce triangles around the boundary
of P, leaving an untriangulated portion Q. Second in Section 5.2 we describe
the formation of another polygon R inside Q, and in Section 5.3 we triangulate
R using a standard mesh generator. Third, in Section 5.4 we describe how to
triangulate between Q and R. The optimality proof relies on a sequence of lemmas
and theorems found in each of these sections. We conclude the two dimensional
algorithm in Section 5.5.
1.4 Some definitions.
We now describe some assumptions about the input and triangulations that apply
to both algorithms.
Face. A face of the input, P, or another polyhedral region, can have either
zero dimensions, called a vertex, one dimension, called an edge, or two dimensions,
called a facet. In some circumstances we want to regard a three dimensional input
11
as a three-dimensional face. Note that the containment relation induces a partial
order on faces.
Conformal. All of the triangulations we consider will be conformaL If we
consider a triangulation as a lattice of faces, then being conformal means that two
faces intersect in a face of the lattice or not at all. For example, there may be no
vertices in the middle (relative interior) of an edge or triangle. In addition, we only
consider triangulations that are conformal to the input, as defined in Section 4.8.
Nondegenerate. We assume that the input region P is nondegenerate in the
following senses. Every facet has the interior of the region on exactly one side. For
every zero-dimensional face (vertex) v, for every small enough open neighborhood
N of v, the set (N n P) --H fvj has exactly one connected component. For the three
dimensional algorithm we also require that every one-dimensional edge is incident
on exactly two two-dimensional faces. These assumptions may be dropped in future
work.
Chapter 2
Other work
This chapter will review selected previous results concerning triangulations. We
focus our attention on maximizing the minimum angle. For a comprehensive survey
of the computational geometry relevant to finite element mesh generation, see Bern
and Eppstein [1991].
2.1 Computational geometry literature
The first question one might wish to settle is the existence of a triangulation.
Preparata and Shamos [1985] have shown that indeed one can triangulate a planar
straight line graph in optimal O(n log n) time, without adding any Steiner points.
However, there is no bound on the shape of a triangle produced by their algorithm.
One of the first algorithms considered that provides a quality bound to triangle
shape is one that produces the well known Delaunay triangulation (DT). For a
point set, the DT is defined as the subgraph of the complete graph, such that each
edge has a circle through its vertices containing no other input vertex. Delaunay
[1934] shows that it may also be described as the planar dual of the Voronoi
12
13
diagram or Dirichiet tessellation, Aurenhammer [1991] has surveyed the literature
of Delaunay triangulations and Voronoi diagrams. Certain steps must be taken if
the input vertices are degenerate (co-circular) to insure that the DT is indeed a
triangulation.
The notion of the DT may be extended to the constrained Delaunay triangu-
lation (CDT) for planar straight line graphs (PSLG), a category which includes
polygons with holes. For a PSLG, the CDT may be defined as the subgraph of
the complete graph on the vertices of the input. An edge ab is in the CDT if the
vertices are visible to each other, and there is a circle through a and b such that
no other input vertex c visible to a or b is contained in the circle. We say that two
points are visible to each other if the edge between them does not cross an input
edge. Lawson [1977] has shown that the CDT maximizes the minimum angle in
the triangulation. Actually, the CDT lexicographically maximizes the angle vector,
the list of angles of the triangulation from smallest to largest. See Lee [1978] or
Lee and Lin [1986] Numerous optimal O(n log n) algorithms have been proposed
for generating the CDT, which are nicely summarized in Bern and Eppstein [1991].
The notion of a CDT can also be used for mesh generation where Steiner
points are allowed. In the 0(n2) algorithm of Chew[1989b], the edges of a PSLG
are divided into edges of about equal length. The CDT of these smaller edges is
computed. While there exists a poorly shaped triangle, we add the center of the
circle through the vertices of the triangle and recompute the CDT. This results
in triangles with angles between 30 and 120 degrees. The algorithm works with
a slightly restricted class of input, and there are no optimality bounds on the
cardinality of the triangulation.
14
Recently Rupert[1992j has developed a mesh generation algorithm based on the
constrained Delaunay triangulation that overcomes these shortcomings. This result
has simultaneous size and quality bounds, similar to those of the earlier work of
Bern, Eppstein and Gilbert [1990] and the three dimensional algorithm presented
in this thesis and previously in Mitchell and Vavasis [1992J. Unlike those works,
Rupert [1992] allows the simultaneous triangulation of the inside and outside of
a polygon. Rupert[1992] is much easier to implement and in practice produces a
triangulation with fewer triangles than Bern, Eppstein and Gilbert [1990].
Despite the power of the CDT in two dimensions, it is unlikely that the CDT will
be very useful in developing three or higher dimensional algorithms. Rupert [1992]
comments on the infeasibility of extending the results to higher dimension. There
are two main problems. The first is that no (constrained) Delaunay triangulation
may exist for general three dimensional input (Sch5nhardt [1928]), but with careful
placement of Steiner points it is conceivable that this may be overcome. The second
more fundamental problem is that Delaunay tetrahedra may have terrible shape:
Even when the angles between edges are large, the angle between a face and an
edge as defined in Chapter 4 may be arbitrarily small. This can happen when four
vertices are nearly coplanar. That is, the results of Lee [1978] and Lee and Lin
[1986] are peculiar to planar graphs. Rajan[1991] shows that CDT minimizes the
maximum radius of a smallest containing sphere in IRd We note aspect ratio is
dependent on the ratio of the radius of the smallest containing sphere to the radius
of the inscribed sphere (Bern [1992], personal correspondence).
One way generate a triangulation with guaranteed element shape is to lay down
a regular subdivision of the plane, and then perturb this subdivision slightly to
15
conform to the input. The subdivision must be fine enough to isolate geometric
features of the input. See for example Baker, Gross and Rafferty [1988], and the
linear time algorithm described in Chew [1989b]. These triangulations may have
output size much larger than optimal, however.
In higher dimensions, Chazelle and Palios [1989] consider the problem of gen-
erating a triangulation with optimal size, but are not interested in angle bounds.
The first two dimensional algorithm to simultaneously guarantee optimal bounds
on triangle quality and mesh cardinality was Bern, Eppstein and Gilbert [1990].
They overcome the mesh cardinality limitations of using a regular subdivision of
the plane by using a quadtree. See Carey, Sharna and Wang [1988] for a discussion
of quadtrees and octrees. Bern, Eppstein and Gilbert [1990] solves many triangu-
lation problems, including that of optimally triangulating a point set of arbitrary
dimension. The general outline for each method is similar, we describe that of the
algorithm that triangulates a two dimensional polygon with no small angles and
optimal cardinality: Start with a large square containing the entire input. Then,
recursively divide any box containing disjoint faces of the input into four smaller
boxes. Element quality is retained by balancing the quadtree, forcing adjacent
boxes to differ in size by no more than a factor of two. Before forming the final
triangulation, faces of the quadtree that are poorly aligned with the input are per-
turbed. Individual boxes are triangulated independently, The mesh size optimality
proof makes use of the minimum angle of the CDT of the input.
Hierarchical decomposition methods are not limited to quadtrees and octrees.
Using a triangular or tetrahedral decomposition often makes the final triangulation
step easier, although the decomposition process can be more complicated. Moore
16
[1992j discusses the properties simplices must have in order to be useful for a recur-
sive subdivision of space. See also Moore and Warren [1990aj. For adaptive mesh
generation for box-shaped regions see also Moore and Warren [1990b]. Mitchell
[1987] and [1988J consider the problem of adaptive mesh generation of complex
regions in regard to finite element error bounds.
2.2 Geometric modeling
In many mesh generation algorithms geometric modeling plays an important role.
For example, in the three dimensional aigorithm presented in this thesis, we must
repeatedly find the intersection between an octree box and the input polygon.
Often mesh generation algorithms use many of the functions found in geometric
modelers, such as intersection. Also, complicated input regions are usually gener-
ated using a CAD system with a geometric modeler. A geometric modeler is also
useful to view to resulting mesh, especially in three dimensions. Sapidis and Peruc-
chio [1989] review existing mesh generation implementations, and look for robust
extensions to three dimensions. They search in vain for a fully automatic mesh
generation algorithm that may be incorporated into an existing CAD system, but
emphasize the desirability of such an algorithm. Here automatic mesh generation
refers to algorithms that require no user input other than the initial polyhedral
region.
2.3 Engineering literature
The engineering literature is rich with finite element mesh generation algorithms.
See for example, the conference proceedings edited by Hauser and Taylor [1986].
17
There are many techniques that do not directly involve computational geometry,
such as conformal mapping. Often even the engineering literature that does in-
volve computational geometry differs significantly from the results developed by
the computational geometry community. Many of these differences arise merely
from style of presentation. For example, most engineering results rely on empir-
ical evidence rather than mathematical proofs to demonstrate the quality of the
mesh. With sufficient effort, algorithms from the computational geometry commu-
nity could provide empirical evidence, and in some cases proofs could be derived
for algorithms from the engineering community.
Sometimes the differences go beyond this, however, and stem from what is
considered a respectable problem in the two communities. For example, in the
engineering literature there are many cases where reasonable element shape is only
achieved for very well shaped input, no mathematical proof of element shape is
possible. What is considered well shaped input is usually defined by the practical
problems that arise, rather than by lower bounds on the degree of vertices and other
mathematical means. Considerable effort has been made to develop heuristics that
improve element shape in practice, but have no provable guarantee of success.
Another example is that in the engineering literature, many algorithms have
been presented that are useful for particular equations that one wishes to solve,
on particularly well shaped domains. As one might expect, these issues have been
largely ignored by computational geometers, who prefer to focus on the geometry
of the input regions. Most computational geometry literature focuses on arbitrar-
ily complex domains. Also, most computational geometry results have ignored
constant factors in mesh size, and often element shape as well. These constant
18
factors can be important in implementations, and so the computational geometry
literature is limited from an engineering standpoint.
One avenue of cooperation between these fields is the following. A computa-
tional geometry result, such as one presented in this thesis, can be used to generate
a mesh. But, either during or after this process, well developed heuristics from the
engineering literature can be applied to improve the mesh quality. As long as the
heuristics do not reduce the quality of the mesh, we would still have a provably
good algorithm. Hopefully, the restrictions of the computational geometry result
would not hinder the heuristics from developing a good mesh in practice. Such
a result would be especially valuable for very large inputs, where verifying and
improving the quality of heuristic results can be tedious for practitioners.
2.3.1 Buraty?ski [1990]
Buratynski [1990] presents several ideas that were useful in developing the three
dimensional algorithm of this thesis. It is also useful as an example of the compu-
tational geometry limitations of a heuristic result. This is not a criticism of that
work, but an illustration of differences in intent.
In Buratynski [1990], the three dimensional input region is assumed to be rep-
resented by its boundary. The method is octree based. Similar to our three di-
mensional algorithm, in both the octree generation and final triangulation steps
input faces are considered on a hierarchal basis. That is first vertices, then edges
and finally facets. In typical mesh generation algorithms, the octree boxes are
deformed to the input. In this method, however, the input geometry is deformed
to the boxes (called bricks by Buratynski [1990]), and then the inverse transform
is performed to both the input and the box to recover the original input. This
19
approach is easier in that there is less to consider when forming the final trian-
gulation of a box intersected with the input. However, mathematically proving a
lower bound on element shape using these transformations would be difficult.
This method works very well in practice on the class of input it is designed for,
namely injection plastic molds. However, it is unclear if it can handle high (> 12)
degree vertices, so there are correctness concerns from a computational geometry
viewpoint. The analysis of the quality of the mesh does not consider some factors
that we consider important in this thesis, but probably are not important for the
class of input the algorithm was designed for. That is, there are no mathematically
proved lower bounds on element quality in terms of the input geometry. There are
limitations on the size of the boxes that can contain a particular input edge, which
limits how close to optimal the cardinality of the mesh may be. Indeed, a lower
bound on the mesh size is not considered.
The algorithm is automatic in the sense that the user does not specify individual
boxes, although there are user specified starting parameters. These parameters
place some limits on how coarse the mesh may be. Also, there are points in the
algorithm that the program must stop and ask the user for additional parameters
if it encounters unusual features of the input. In light of these user inputs, it is
doubtful that nontrivial mathematical bounds on mesh cardinality and element
shape exist. However, these inputs appear to be valuable practitioners.
Chapter 3
Element shape
In two dimensions we seek a triangulation whose minimum angle is as large as
possible, up to a constant factor. In three dimensions, we seek a triangulation
whose aspect ratio is as small as possible, up to a constant factor. Here the
minimum angle of a triangulation is the minimum of the angles of its triangles,
and the aspect ratio of a triangulation is the maximum of the aspect ratios of its
tetrahedra.
In three dimensions, in order to argue about the aspect ratio of triangulations
of P, we wish to define the angle between any two faces of P or a triangulation of
P, that intersect each other non-trivially. Various concepts such as dihedral angle
and solid angle have been used in the past, but we adopt a different approach.
Interior angle. We define the interior angle or angle between two faces f and
g in the cases where f and g are a facet and an edge, two facets, or two edges
that have a common intersection v. We additionally require that one face is not
a subset of the other. We say that two rays ri, r2 with a common endpoint v can
see each other if there are points vi, V2 on ri, r2 each at distance t > 0 from v such
20
21
that the sector vviv2 lies in P.
If f and g are edges, then we measure angles in the usual way for rays, except
we require that the rays can see each other. We always measure angles with respect
to the interior of P, so that f and g could make an angle close to 36O?. This is
also how we define interior angle in two dimensions.
If f is an edge and g is a facet, then we consider all rays r in the plane of g,
such that r lies between the two edges of g that meet at v, and such that f and r
can see each other. Then the angle between f and g is defined to be the minimum
angle between the ray covering f and all rays r meeting the preceding conditions.
We have minimum rather than infimum because faces are assumed to be closed.
If f and g are both facets whose intersection is a vertex, then we let both the
ray for f and the ray for g vary as in the last paragraph, and take the minimum
angle between the rays to be the angle between f and g that can see each other. If
the facets share a common edge, then the interior angle between the facets is the
angle between the rays in the half planes of f and g perpendicular to the common
edge. Once again, we only consider rays that see each other.
We use the symbol "Q" in Chapter 4 to refer to the smallest such interior angle
between any two faces of P. We note that if the interior angle between two facets
is P, then one of them has a subface edge that makes an angle of no more than p
with the other. Thus there are always either two edges, or an edge and a facet,
that have interior angle o between them.
Aspect ratio. A useful measure of the shape of a tetrahedron is its aspectratio,
which we define to be the ratio of the radius B of the smallest containing sphere,
to the radius r of the largest inscribed sphere. In general the smallest containing
22
sphere is not the sphere through the vertices of the simplex (circumscribing sphere).
A simplex with bounded aspect ratio has bounds on the angles made between any
pair of faces that intersect, and not just on the angles between edges.
It can be proved that the aspect ratio of a tetrahedron is bounded above and
below by constant multiples of the reciprocal of the sharpest interior angle (as
defined above) in the tetrahedron. We do not prove this claim because it is implicit
in theorems and lemmas in upcoming sections in Chapter 4. This is also easy
to show in two dimensions. Thus, up to constant factors, aspect ratio and the
reciprocal of min angle are equivalent measures. In Chapter 4 we will make heavy
use of this fact because the min interior angle is so much easier to analyze.
Sphere ratio. Often another measure goes by the name of aspect ratio. It is
what we will call sphere ratio. The sphere ratio of a d dimensional simplex is the
ratio of the radius of the d dimensional sphere through the vertices of the simplex,
to the radius of the largest sphere contained in the simplex. While related to aspect
ratio, the sphere ratio is not equivalent, even up to constant factors.
Bern [1992j (personal correspondence) has proven that in d dimensions the
sphere ratio of a simplex is between A and A2, where A is the aspect ratio defined
above. Bern and Teng [1992] (personal correspondence) have extended this result
to show the equivalence of aspect ratio with several measures. For example, A is
bounded above and below by constant factors times the ratio of the largest distance
between a pair of supporting hyperplanes to the smallest distance between a pair of
supporting hyperplanes, and also by constant factors times the inverse of smallest
solid angle.
A previously known consequence is that if a triangle has no large angles, its
23
(x
Figure 3.1: The sides and angles of a general triangle.
sphere ratio and aspect ratio are within constant factors of each other.
Sphere ratio is unwieldy, even in two dimensions. We present the following
approximation of sphere ratio for two dimensions.
Theorem 1 Let 5 be the sphere ratio of a triangle, a its smallest angle, and ? its
largest angle. Then
3			1			1
<5<4.
8sinosinP --H --H sinasiuP
Proof. Without loss of generality, assume that the longest edge of the triangle,
which is opposite p, has length 1. Let the shortest edge, which is opposite a, have
length a. Denote the remaining angle by ?, and the length of the opposite edge by
c. See Figure 3.1.
We first bound r, the radius of the inscribed circle. From the hypotenuse leg
theorem of elementary trigonometry, the center of the inscribed circle is at the
intersection of the angle bisectors. llence, from Figure 3.2, we have r = l? tan2
2'
r = la tan a2, and l? + la = 1. Eliminating l? and la from these equations, we have
1 _			a
r -cot--H + cot-.
2			2
a
Since by assumption a <?, we have tan 2 <r <tan a2 Since 0 < a <7r/3, we
24
r
1T `(1
Figure 3.2: The radius of the inscribed circle.
Figure 3.3: The radius of the circumcircle.
have ?21 sin Q <tan ?a2<??23 sin ?. Hence
1o+			2
--H sin ? < r < --H sin a.
4			--H			3
We next examine R, the radius of the circumscribing circle. see Figure 3.3.
Drawing the segments from the center of the circumscribing circle to each of the
vertices of the triangle, we obtain two triangles. The first is a similar triangle with
two sides of length R and the remaining side of length a. We define Pa to be the
angle at a base vertex of this triangle. The second triangle is also a similar triangle,
with two sides of length R, but the remaining side is of length c. We define Pc to
25
be the angle at a base vertex of this second triangle. Then
2 coS Pc
We wish to convert the above equality involving Pc and c into inequalities
involving only P. Write P = --H E, with 0 < ? <2r/3. Note Pa + Pc =
We obtain an upper bound on 1? involving ? as follows. Since c> a, we have
Pc> Pa, and Pc < = _ --H ?2E. lience cosPc > cos(?r2 --H 2E) = sin2C > ?Sifl2E = sin2?
2
Since c is the length of the second longest side of the triangle T, we have c < 1.
Hence
--H sinP
We obtain a lower bound on R involving ? as follows. Since Pa is the base
angle of a similar triangle, Pa < Hence Pc = P --H Pa>??T2 --H ?. Hence cos Pc <
cos 2 --H = sin E = sin P. Since c is the length of the second longest side of the
triangle T, from the triangle inequality we have c > 21. Hence
4sinP
Combining the bounds on 1? with those on r, we obtain the following bounds
on the sphere ratio S = Rr of
3			1			1
8sinasinP --H --H sinasinP
Chapter 4
Quality mesh generation in three
dimensions
We show how to triangulate a three dimensional polyhedral region with holes. We
assume there are no restrictions on where we may place Steiner points. Our tri-
angulation is optimal in the following two senses. First, our triangulation achieves
the best possible aspect ratio up to a constant factor. Second, for any other tri-
angulation of the same region into m triangles with bounded aspect ratio, our
triangulation has size n = O(rn). Such a triangulation is desired as an initial
mesh for a finite element mesh refinement algorithm. Previous three dimensional
triangulation schemes either worked only on a restricted class of input, or did not
guarantee well-shaped tetrahedra, or were not able to bound the output size. We
build on some of the ideas presented in previous work by Bern, Eppstein, and
Gilbert, who have shown how to triangulate a two dimensional polyhedral region
with holes, with similar quality and optimality bounds.
26
27
4.1 Subdivision of P into cubes, the octree
4.1.1 Data structures and definitions
Our triangulation is denoted AocT. We assume we are given a three-dimensional
polyhedral region as follows. There is a list of vertices with coordinates specified.
There is also a list of edges and faces. The three lists are mutually linked. The
list of edges for each vertex is ordered as they occur around the vertex. Similar
orderings are assumed for the other lists. We assume that P is connected; if not,
each component could be triangulated separately.
The main data structure we use for our algorithm is an octree. We commonly
refer to each node of the tree as a box. We associate with each box, b, a polyhe-
dral region of JR3 called the embedding of the box and denoted 1(b). During the
generation of the octree, 1(b) is exactly a three dimensional cube. Later boxes are
warped and triangulated, changing their geometric structure.
An octree node is either a leaf, or has eight children. The embedding of the
eight children of a node are the eight cubes obtained by dividing the embedding of
the box in half in each of the three dimensions. We say a node is split if it is not
a leaf. The process of creating the eight children of a node is called splitting.
For clarity, we use the term box face ("box facet," "box edge," etc.) to denote
a face of a box and P face to denote a face of P. Similarly for a? and ?b.
Duplicate. Some boxes in the octree will be duplicated into the original node
and several new nodes, called duplicates or duplicate boxes. We create duplicates
whenever the intersection of the box with P has more than one connected comp?
nent. Each duplicate represents the same geometric cube in JR3, but is associated
with one connected component of P A 1(b). We use the notation P A b to denote
28
the component of P A 1(b) assigned to a particular box 6.
Whenever we split a box 6, if P A 6 is nonconvex a child box U may have more
than one component (i.e., (P A 6) A 6' might have more than one component even
though PA 6 consists of one component). Whenever a box is split, we immediately
determine whether any of its children contains more than one component of P,
and if so, we make duplicates of the child box.
If a P face is incident upon P A 6 for a box 6, we say that a box contains the
face. We maintain pointers between each box and the faces it contains.
Extended Box. For a given box 6, we define the extended box of 6, ex(b), such
that I(ex(b)) is the cube concentric with 1(b) and with each dimension expanded by
a factor of 5. We use the notation PAex(b) to denote the component of PflI(ex(b))
that contains PA 6. The extended box contains only the P faces of PA ex(b). Note
that ex(b) is not a box of the octree, but may be constructed from boxes of the
octree.
Adjacent. Two boxes are called neighbors if their embeddings intersect non-
trivially, and one is not a duplicate of the other. The size of a box is the length of
an edge. We say that box b1 is adjacent to box b2 if they are neighbors and there
is a point of P common to both, i.e. if PA bi A P A 62 $ $. In certain cases such as
when a box contains faces meeting at a reflex angle, a box may be adjacent to two
duplicates of a second box. We say bi, 62 are balance-adjacent if they are neighbors
and there is a point of P common to one of the boxes and the extended box of the
other, i.e. if(PAex(bi))A(PA62) $ ? or if (PAex(b2))fl(PAbi) $ $. Boxes that
are adjacent are balance-adjacent, but not all balance-adjacent boxes are adjacent,
see for example Figure 4.3.
29
Balance Condition. If we split a box, we immediately split other boxes to
maintain the following invariant called the balance condition: No box is balance
adjacent to a box more than twice its size. Certain boxes containing vertices or
edges of P are called protected boxes, and are exempt from being split by the
balance condition later in the algorithm. Nonetheless, the ratio of the size of a
protected box to an adjacent box is bounded above by a constant that depends
linearly on 1/a, as we prove in the next section.
4.1.2 Creating the octree
We generate the octree by selectively splitting and duplicating nodes. The goal
of splitting and duplicating boxes is to make the boxes small enough so that the
intersection of P with the embedding of any box is as simple as possible. However,
boxes should not be made too small, as this would lead to an excessive number of
tetrahedra in the final triangulation.
The octree generation algorithm is divided into three phases, the vertex phase,
the edge phase, and the facet phase.
It is easy to state the conditions under which we duplicate b, namely, whenever
P n 1(b) has more than one component. Nonetheless, the box duplication process
is actually the most computationally complicated part of the octree generation al-
gorithm, because determining components of P n 1(b) is a nontrivial computational
task. The details of the duplication process are described below.
On the other hand, for splitting, the conditions under which a box is split
are more complicated, but the computational tasks are straightforward. We now
describe the various steps of the splitting algorithm.
30
4.1.3 The vertex phase
How finely we split boxes depends on the following definition.
Vertex cone. A vertex cone of a box b is a set of P faces, F1, ..... . , Fk, that
satisfy the following:
(a) The set is precisely one vertex and all of its superfaces. This vertex is called
the apex of the vertex cone.
(b) The apex is in b.
(c) The faces F1,F2,... , Fk are exactly the P faces incident upon P A ex(b).
Vertex crowded. We say that a box b is vertex crowded if the following are
satisfied.
(a) There is a P vertex v in b.
(b) The superfaces of v are not the only faces of P incident upon P A ex(b).
Equivalently, a box is vertex crowded if P A b contains a vertex that is not the apex
of a vertex cone. We say that a face G of P A ex(b) interferes with v E P A b if
v is not a subface of 0. Thus b is crowded if it contains a P vertex v, and ex(b)
contains a face that interferes with v.
We recursively split and duplicate boxes, maintaining the balance condition at
each step, until every box containing a vertex contains exactly one vertex cone.
Clearly this procedure will terminate once the box sizes become a constant factor
smaller than the minimum path length in P between a vertex and a face that does
not contain that vertex.
31
The conclusion of the vertex phase is a special one-time reorganization of the boxes
so that each vertex of P is far from the boundary of the box that contains it. The
value of this property will become clear in Section 4.3 and later sections.
We first make the following observation concerning the boxes at the conclusion
of the vertex phase. Let b be a box containing a P vertex v. We claim that any
box balance-adjacent to b has size at least that of b. We denote the size of b with
the notation h(b). Let b1 be balance-adjacent to b. The balance condition ensures
that h(b1) is either equal to h(b), 2h(b), or h(b)/2. Hence we need only prove that
h(b1) # h(b)/2.
To prove this, suppose that h(b1) = h(b)/2. Let b'1 be the parent box of b1.
Then either U1 is vertex crowded, or it is split because of the balance condition. We
The vertex centering step
The details of the procedure for determining whether a box b should be split are
straightforward given an enumeration of the P faces bounding P A ex(b). Such an
enumeration is obtained from the procedure for determining whether a box should
be duplicated described below.
Figure 4.1: Here a box is duplicated for two components, and one duplicate box is
split because it is crowded.
M
-I
32
may assume that b'1 is split by the balance condition, since if it is vertex crowded we
simply take k = 1 in the following discussion and the proof holds. There is a chain
of boxes ..... . , 61k such that b'k is vertex crowded, Uj is balance-adjacent to b'j+1
(for i ..... .,k--H1), and such that h(Uj) = 2h(Uj+1) These are just the conditions
necessary for the balance condition to propagate from Uk back to Ui Note that
I(Uk) lies interior to I(ex(b)), no matter how large k is, because ex(b) is 5 times
larger than b in every direction, and the sizes of the Uj's form a decreasing geometric
series. Thus, the vertex causing Uk to be crowded is in PA ex(b), contradicting the
fact that b is not vertex crowded.
Thus, every box b with a vertex is surrounded by boxes that are either twice the
size or the same size as b. The first part of the vertex centering step is as follows.
For every box b with a vertex, we split the boxes balance-adjacent to b that are
twice as large, and then we propagate the balance condition (which causes every
other box in the octree to be split at most one time).
We claim that after this procedure, every box balance-adjacent to a box b
with a vertex is the same size as b. This claim is clearly true before the balance
condition gets propagated. Thus, we have to show that when the balance condition
is propagated, this does not cause any box balance-adjacent to a box b with a vertex
to be split.
Consider a setup of boxes that might cause a box b1 balance-adjacent to a
box b to be split by propagation of the balance condition. Suppose b1 is a box
balance-adjacent to b, where b contains a vertex v, such that b1 is the same size as
b. Suppose elsewhere there is a box t with balance-adjacent bk such that bk is the
same size as t, and bk is balance-adjacent to a box bk?1 four times as large, which
33
b
Figure 4.2: No matter how small we split a box b', if it is outside of ex(b) then the
balance condition will not split 6.
is balance-adjacent to a box 6k--H2 twice as large, and so on, up to bi, which is twice
as large as 62. Under these conditions (and only these conditions) propagating the
balance condition from bk would cause 61 to split.
However, we claim that the conditions described in the previous paragraph
cannot occur. This is because such conditions would contradict the fact that every
facet contained in ex(b) has vertex v as a subface (same argument as above).
Thus, after the first part of the vertex centering step, every box has been split
at most once, and every box with a vertex is balance-adjacent to boxes of the same
size. Consider such a box 6. Without loss of generality, 6 has 26 balance-adjacent
boxes surrounding it of the same size, arranged in a 3 x 3 x 3 group. These boxes
may have been duplicated. If necessary, we adjoin empty boxes in the spaces of
this group where no box is balance-adjacent to 6.
Now, the second part of the vertex centering step is as follows. We reorganize
the octree by uniting 6 with some of the boxes of the 3 x 3 x 3 group, and perhaps
some duplicates of these as well. We form a new box B by uniting 6 with seven
34
Figure 4.3: Uniting boxes around a vertex box in two dimensions.
boxes in a 2 x 2 x 2 group, and perhaps some of their duplicates, to obtain the result
that 1(B) is a cube of side length 2h(b), and B contains the component of P A B
containing the P vertex in b. See Figure 4.3 for an illustration of this procedure in
two dimensions. Note that we can always pick the correct 7 neighbors of b, so that
v is distance at least A(b)/2 from any facet of B. Since we have effectively doubled
the size of b, B may be balance-adjacent to boxes four times smaller than B. We
do not take any steps to repair this violation of the balance condition, as the size
ratio of adjacent boxes is still bounded.
From now on, B is protected, meaning that it is never split again during the
algorithm.
4.1.4 The edge phase
The next phase of the octree generation algorithm focuses on edges of P, and is
analogous to the vertex phase.
We define "edge cones" and "edge crowded" and split and duplicate boxes as for
35
vertices. The main difference is that we take an approach different from centering
to ensure that no box face is close to an edge. The splitting and duplicating
algorithm proceeds as if the protected vertex boxes do not exist. In particular, the
extended box of any box in this phase does not extended into a protected box, and
the vertices of P are ignored.
Edge cone. An edge cone of a box I) is a set of P faces, F1, F2, F3, that satisfy
the following:
(a) The set is precisely one edge and its two superfaces. This edge is called the
apex of the edge cone.
(b) The apex is contained in 6, and
(c) The faces F1, F2, F3 are exactly the P faces incident upon P A ex(b).
Edge crowded. We say that a box 6 is edge crowded if it contains an edge
of P, but it and its superfacets are not the only P faces incident on P A ex(b).
Equivalently, 6 is edge crowded if any edge contained in 6 is not the apex of an
edge cone.
Just as in the vertex case, we recursively duplicate an unprotected box for each
edge cone it contains, and split it if it is edge crowded. We maintain the balance
condition immediately after a box is split. Recall that the protected vertex cone
boxes are not changed in this step.
Identifying edge cones is straightforward if we have a list of facets bounding
P A ex(6), which is obtained from the algorithm for duplication described below.
36
Protecting edge cone boxes
We seek the same goal as the centering vertices step but for edges. The fact that
an edge cone box may be balance-adjacent to a protected vertex box prompts us
to take an approach somewhat different from the one used for vertex cones. Here
if a P edge is near a box face G, we surround the face with protected boxes, as in
the vertex case, but do not require them to be the same size. Unlike the vertex
centering phase, we do not eliminate G at this time, because G may be a face of
a protected vertex box. Instead, after the entire octree is generated but before
we triangulate, we warp the cube (see Section 4.3) into a more general shape that
increases the distance between 0 and the P edge, and eliminates any small angles
as well.
We begin by splitting every unprotected box containing a P edge and its
balance-adjacent boxes and enforcing the balance condition. Every box is split
at most once. This ensures by the "crowded" definition that if a box b (other than
a protected vertex box) is balance-adjacent to a box b' containing an edge E, then
b cannot be balance-adjacent to a box containing any other edge E'
We always protect every box containing a P edge, and every box balance-
adjacent to a box containing a P edge, regardless of whether any face of the box
is close to a P edge. These boxes are never split again.
4.1.5 The facet phase
The third phase of the octree generation algorithm focuses on facets of P. We
define "facet cones" and "facet crowded" as in the previous phases. We split and
duplicate boxes as we did for vertices. As for edges, the main complications arise
37
when centering facets, especially when a facet cone box is balance-adjacent to a
vertex or edge cone boxes. The splitting and duplicating algorithm proceeds as if
the protected vertex and edge boxes do not exist. In particular, the extended box
of any box does not extended into a protected box, and the edges and vertices of
P are ignored.
Facet cone. A facet cone of an unprotected box b is a single P facet, F1, that
satisfies the following:
(a) Facet F1 is contained in b.
(b) Facet F1 is the one and only P face incident upon P A ex(b).
For the sake of consistency, we will call P1 the apex of the facet cone.
Facet crowded. We say that a box b is facet crowded if it contains a facet that
is not the only P facet incident upon P A ex(b). Equivalently, b is facet crowded if
it contains a facet that is not the apex of a facet cone.
Just as in the vertex case, we recursively split a box if it is facet crowded. We
maintain the balance condition immediately after a box is split. Recall that the
protected vertex and edge cone boxes are not changed in this step.
Protecting facet cone boxes
We take the same approach as for edges. We begin as in the edge case by splitting
every unprotected box containing a P facet and its balance-adjacent boxes and
enforcing the balance condition. This ensures that a box is balance-adjacent to at
most one facet cone apex. Thus, in the warping and triangulation step to follow
in Section 4.3, we may consider the boxes containing a facet without reference to
38
any other facets, except for boxes balance-adjacent to a protected edge or vertex
box.
Technically, there is no need to protect boxes because there are no more splitting
phases to consider. However, for the sake of consistent terminology, we protect facet
cone boxes and the boxes balance-adjacent to them.
4.1.6 The duplication process
Recall that we duplicate b whenever it is determined that PnI(b) has more than one
component. In this section we explain how to identify the components of P n 1(b).
This same techniques allow us to determine the components of P n 1(ex(b)), which
is necessary for determining whether boxes are crowded, and when to propagate
the balance condition.
Because we believe that the algorithm proposed here is not the most efficient
algorithm possible for the task, we skip over some of the details involved.
The first part of the duplication algorithm is a preprocessing algorithm, which
is run before the octree generation begins. In the preprocessing algorithm we first
identify the separate components of ap. If P is a convex domain (or some other
fairly simple shape) then a? has a single component and no further preprocessing
is necessary. Otherwise, suppose a? has several components. One component is
the exterior component; we call all the other components floaters. The component
for each facet, edge, and vertex can be determined in 0(n) steps with a standard
graph search, where n is the total number of facets edges and vertices of P.
Once the components are identified, we now perform a second procedure in
which we identify a tether for each floater. First, for each floater C, we identify
the point v (necessarily a vertex) with the largest x coordinate on C. This point
39
Figure 4.4: Floaters and tethers for a two dimensional polygon.
is called the base of the tether. From the base we shoot a ray in the positive x
coordinate direction until we encounter the first point ?? on 8P. This point cannot
be on the same floater. Such a point exists, since the ray must eventually pass
through the exterior component. The directed segment from v to v? is called the
tether for C, and v' is called the head of the tether. See Figure 4.4 for a two
dimensional example of the tethering structure of a polygon. We can carry out the
search for v`in a fairly naive fashion: Simply loop over all facets, and for each facet
F carry out a point-in-polygon test for the point where the ray v passes through
the plane of F. This requires time proportional to the number of edges of F. Thus,
the total time to locate v' for a particular tether (looping over all facets) is 0(n).
Therefore, the time to identify all tethers is 0(n2).
Note the following two properties of tethers: (1) Each tether is contained com-
pletely in P. (2) When regarded as a digraph on the components of 8P, the tethers
form an in-tree rooted at the exterior component.
40
Now we discuss the procedure for determining when a box should be duplicated.
We first make a list of all P facets, edges and vertices contained in the box. A
particular P facet F in the box may intersect the box in more than one component
if the facet is nonconvex. We therefore must first determine all the components of
each facet passing through the box. This may be done by sweeping a line segment
across the facet, requiring O(v log v) steps if v is the number of vertices in the facet.
The total time for identifying all components of all facets for a box is O(n log n).
Then we perform a combinatorial graph search to identify the various comp???
nents of a? n 1(b). Let these components be called sheets. There are two kinds
of sheets: those that intersect the boundary of 1(b), and those that are completely
internal to 1(b). Note that a sheet entirely inside 1(b) must be a floater; such a
floater is called an internal floater for b. The other sheets are called external sheets.
Note that a floater that intersects 1(b) is an internal floater if and only if it lies
entirely in 1(b), otherwise it is one or more external sheets.
Clearly any sheet is incident upon a single connected component of P A 1(b)
since the sheet itself is connected. Thus, the remaining task is to determine which
sheets are connected to each other in P A 1(b).
For external sheets, the following algorithm exactly determines how they are
connected. We consider the intersection of each sheet with the boundary of 1(b)
(the surface of a cube). This intersection is a collection of disjoint simple polygons.
We form the list of polygons for all sheets. We call these polygons the "rims." Note
that these polygons break across edges of the cube. We then sweep over this list
of rims, forming a tree indicating which polygons are contained in others. If these
polygons were in the plane, we could sweep with a vertical line. This would require
41
sorting x coordinates of the polygons, hence requiring O(n log n) steps. On the
surface of the cube we can sort by the sum of x, y, z coordinates, which means that
we sweep over the surface of the cube with a closed polygonal segment (initially
the boundary of a triangle).
In the previous paragraphs we made reference to point-in-polygon tests and
plane sweeps. These are well-known tools in the computational geometry litera-
ture. The reader unfamiliar with these tools should consult Edelsbrunner [1987j
or Preparata and Shamos [1985j.
From the sweep we can build a rooted tree indicating the containment hierarchy
among the polygons. If we assume that the first point swept is outside P, then
the following simple rule determines connectivity in P: Every polygon at an even
level i (counting 0 as the top level) is connected in P to its children at odd level
i + 1. The parities are interchanged if the first point swept is inside P.
Thus, from an O(n log n) procedure we can determine exactly which external
sheets are in the same components of P n 1(b). This leaves the problem of deter-
mining the component of P n 1(b) for an internal floater C. This is determined
from the tether (v, v') of C. If the head v' is inside 1(b), then the floater lies in the
same component as the facet containing v' (because the entire tether is in P). This
facet might be a facet of another internal floater, but because the tethers form an
acyclic graph, eventually we will reach a point ?? on an external sheet. See Figure
4.5.
The other case is that the tether's head v' lies outside 1(b). In this case we find
the point v11 where the tether intersects the surface of 1(b). If we knew where in the
polygon tree v" lay, then we could determine the component containing v11 (which
42
v"
Figure 4.5: Determining the components of a box. Sweeping the boundary deter-
mines two components, and following tethers determines that all internal floaters
belong to the same component, and there are no internal components.
is the same component containing the floater). The polygon regions containing v?
are easily determined if, before starting the polygon containment sweep, we first
insert v" into the list of polygons as a singleton polygon. This should be done for
all tethers. This does not change the running time bound of O(n log n).
Thus, in time O(n log n) we can determine all the components of P n 1(b).
Moreover, it is easy to see that the rims allow us to determine which neighboring
boxes are adjacent to a given duplicate box.
Finally, we need to determine the facets adjacent to PAex(b) for testing crowd-
edness. This is done by running the component-determination algorithm on ex(b),
and then saving the component of P n I(ex(b)) that contains a point in PA b. This
also allows us to determine which neighboring boxes are balance-adjacent to b.
43
4.1.7 Summary of the octree algorithm
The octree algorithm has now been completely specified. The octree structure,
including adjacency, duplication, and P face containment properties remain un-
changed by the algorithms that follow.
1. We initialize by setting the octree equal to a single node representing a large
box with all of the faces of P in its interior.
2. Vertex Phase. We split and duplicate boxes until it is determined that no
box is vertex crowded. The balance condition is enforced each time a box is
split.
3. Every box b' balance-adjacent to a box b with a P vertex is split, in the case
that h(b') = 2h(b). The balance condition is enforced.
4. Every vertex is centered in its containing box by merging eight boxes around
a vertex. Every box containing a vertex is now protected, so is never split
again.
5. Edge Phase. We split and duplicate unprotected boxes until it is determined
that no box is edge crowded. The balance condition is enforced each time a
box is split.
6. Every box containing a P edge and its balance-adjacent boxes are split to
further separate edges from each other. The balance condition is enforced.
7. Every box containing an edge and its balance-adjacent boxes are protected,
and never split again.
44
8. Facet Phase. We split and duplicate boxes until it is determined that no box
is facet crowded. The balance condition is enforced each time a box is split.
9. Every box containing a P facet and its balance-adjacent boxes are split to
further separate facets from each other. The balance condition is enforced.
An important feature of the algorithm is that when it terminates, every leaf
box is either empty, lying entirely inside P, or contains a vertex, edge, or facet
cone. Also recall that the running time per box is proportional to n log n, where n
is the number of faces of P.
We now consider the running time of constructing the octree. We believe that
our current running time analysis is suboptimal, so we omit some of the details.
First, there is a slight addition to the algorithm as described so far that gives
us a better bound. If a box is crowded, with 0(n) additional time per box we
determine the closest face G interfering with some F, where n is the number of
faces of all dimension of P. Here for closest, distance is distance within the
box. We immediately split the box nonuniformly down to a small level around the
closest point of F to G, so that G does not interfere with F in one of the descendent
boxes containing F. We thus get at least one uncrowded box every time we test if
a box is crowded. The running time of this split is dependent on the distance in
the box between F and G. This time may be amortized against the boxes created
by the balance condition when balancing the octree after this split.
This allows us to claim that the total number of boxes constructed by the
octree algorithm is bounded by a constant multiple of the number of leaf boxes
that contain a point of P. From the warping and triangulating rules to follow,
each such leaf box leads to at least one tetrahedron. In particular, the number of
45
such leaf boxes is bounded by ?, the size of the output. Finally, the total amount
of time spent on each box is at most n log n where n is the number of faces of
the polyhedral region P. Thus, a bound on the running time is O(?n log n). Note
that this subsumes the 0(n2) preprocessing to find tethers, since ? = ?(n). See
Section 4.9 for more remarks on this bound.
4.2 Relative sizes of adjacent octree boxes
In the last section, a box may have been protected in a vertex or edge phase of
the octree generation algorithm. A protected box is not subject to the balance
condition in subsequent phases, so boxes adjacent to it may be smaller than half
its size. In this section we prove that the ratio of the sizes of two adjacent boxes
cannot be smaller than a constant multiple of ?. Here c is the smallest interior
angle of P as defined in Section 3. The various kinds of boxes must be considered
as separate cases.
Observe that from the centering algorithm, if the vertex was close to any face
of its containing box, that face was deleted. Thus we have the following lemma.
Lemma 1 The distance between a vertex of P and any face of a box is at least
one eighth the size of the box containing the P vertex.
Immediately after boxes are split in the vertex phase, each P vertex is a vertex
cone apex in the box that contains it. That is, the only P faces in the extended box
of a box containing a P vertex are superfaces of that P vertex. The extra splitting
to further separate vertices decreases the size of boxes (and hence extended boxes).
The merging of boxes to center the P vertices increases the size of boxes, but boxes
are no bigger than they were before the extra splitting step. Hence, even after the
46
separating and merging, the superfaces of a P vertex are the only P faces contained
in the extended box of the box containing that P vertex. So each P vertex still
satisfies the definition of a vertex cone apex.
4.2.1 Vertex cone boxes
Theorem 2 Let B be a box with vertex cone apex v. For any box b adjacent to B,
h(b) > k a h(B), where k is a constant.
Proof. The proof divides into a series of lemmas, depending on the phase in which
the parent of b was split. If the parent of b was split in a separating phase, we
note that any box is split at most once in a centering phase, so that the size of b
is bounded by the size of its first ancestor that was split in a crowded box phase.
Let b' be this ancestor. So, we can ignore the separating phases, and the proof of
the theorem divides into three lemmas depending on when b' was split, and one
corollary.
Lemma 2 If b' was split during the vertex box stage, then h(b) > h(B)/2.
Proof. The balance condition ensures that h(b') is at least one half of the size of
the box containing v before the centering phase begins. 1
Lemma 3 If b1 was split during the edge box stage, then h(b) > k.a.h(B), where
k is a constant and a is as in Section 3.
Proof. The proof breaks up into two cases that consider why b' was split. We
first consider the more difficult case, that when b' was split because of the balance
condition. Let d be the edge crowded box that was split causing the balance
47
condition to split b1. We now divide into two subcases depending on the position
of d relative to B.
Subcase 1. If 1(d) lies entirely outside 1(ex(B)) or intersects ?I(ex(B)) the
balance condition gives the result: Now matter how small a box is that intersects
?I(ex(B)), a box that intersects ?I(B) will be split to a size no smaller than half
the distance from ?I(B) to ?I(ex(B)). This is very similar to Figure 4.2. If 1(d)
actually lies outside of 1(ex(B)), then the same bound holds.
Subcase 2. Otherwise 1(d) lies strictly interior to 1(ex(B)). Observe then
that for two boxes, one either contains the other or their intersection is a subset
of their boundaries. So, there is a box at least half the size of d between 1(d) and
?I(ex(B)). If d is split down to a size one half that of this buffer box, for each of
the children d1 of d, I(ex(di)) is a subset of I(ex(B)). Note that d may have to be
split at most twice to achieve this size, so it is sufficient to bound the size of d1,
by showing that a box d that is actually half the size of the buffer box has a lowe?
bound on its size. That is, we may assume that d is actually half the size of the
buffer box, if it was larger then this implies that b is larger.
The P vertex in B can not interfere with the edge in d, so the only way that
a face can interfere with the edge contained in d is for that face to be a facet or
edge of the vertex cone for B. Let f be the edge in d and g the closest face that
interferes with f. Note that g may not be the face that makes the sharpest angle
with f at v because of alignment odities, but of course g makes an angle p with f
at least this sharpest angle.
The distance from f fl ex(d) to v is at least the distance from a face of 1(B) on
the shortest straightline path to v. So by Lemma 1 and trigonometry, we have that
48
dist(f fl ex(d),g fl ex(d)) > k p h(B). The dependence is actually on sin(P/2).
If p is large, then the balance condition between d and B in the vertex centering
step will be the limiting factor to how large d is. Otherwise, sin(P/2) is bounded
above and below by constant multiples of p.
Since g intersected the extended box of d, and f intersected B, we also have
dist(f n ex(d),g A ex(d)) < 3?h(d). Hence we have the result of the lemma for
d. But since the balance condition caused b' to split after d was split, bi must be
smaller than d, and we have h(b') > 2h(d).
The second case is if b' was split because it itself was crowded. If b1 intersects
?I(ex(B)) as well as aI(B), we see that h(b) > h(B) as before. Otherwise, b is a
special case of a box d in the interior of 1(ex(B)) in the previous case, so that the
proof of the last case proves the desired result. I
Lemma 4 (Corollary) If d is an edge protected box such that 1(d) is contained
in I(ex(B)) and further 1(d) A ?I(ex(B)) = ?, then h(d) > k p h(B), where P is
the sharpest interior angle an edge in B makes with any superface of v.
Proof. This was almost proved in the proof of the previous lemma. For the
completion of the proof of the corollary, we need to consider a box strictly interior
to 1(ex(B)) that contains an edge, but whose parent was split because of the
balance condition and not because it was edge crowded. We note that the last time
its ancestor was split because it was crowded it had the desired size bound from
the proof of the theorem. If the balance condition splits this ancestor, it must first
split a balance-adjacent non-crowded box. Once this is done, the children of this
balance-adjacent box form a buffer of two same sized non-crowded boxes between
the ancestor and crowded boxes. The balance condition can never penetrate the
49
inner box of this buffer by splitting crowded boxes balance-adjacent to the outer
box. Also, splitting other boxes balance-adjacent to the ancestor only once will
not cause the ancestor to split again, so we may assume that there are two layers
of same sized boxes surrounding the first generation children of the ancestor, and
so the ancestor is never split again.
Hence we have the result with a constant k half as large as for the case in
the lemma where the box was split because it was edge crowded. The proof of
the corollary for a protected box that does not actually contain an edge but is
protected because it is balance-adjacent to a box containing an edge follows with
a constant one quarter as large, from the balance condition and the fact that it is
balance-adjacent to a box where the result holds. I
Lemma 5 If b' was split during the facet box stage, then h(b) > k a h(B).
Proof. This is the most difficult of the three lemmas comprising the proof of
Theorem 2. Many cases may be handled in the same way as in the edge lemma.
Difficulties arise when considering a box that contains two facets that meet at an
edge, and this is where the corollary to the previous lemma is crucial. The depen-
dence on a in the lemma arises from the possibility of small angles in the following
discussion: if all of the angles considered below are large, the factor of a may be
eliminated from the result. As such, the proof concentrates on the possibility that
the angles considered below are small. Certain anomalies occur in the following dis-
cussion when the angles under consideration are large. Except where noted, these
anomalies are easily handled, and the resulting tedious discussion omitted. in this
proof there will be a series of constants k1, k2These labels are independent of
the labels of other such constants appearing elsewhere.
50
If b1 was split because of the balance condition, consider the crowded facet box
d that caused us to split b1. In fact, we consider d to be the last descendent of this
box that was split because it was facet crowded.
If 1(d) does not intersect 1(ex(B)) or intersects ?I(ex(B)) the same argument
as in the last lemma gives the desired result.
So it remains to consider the case where d is contained in the strict interior of
1(ex(B)). Let f be the facet contained in d, and g the closest interfering facet in
ex(d). If f and g meet only at v and do not share an edge, then the same argument
as in the last lemma, substituting faces with the same symbols, proves this lemma.
We are left with considering a facet crowded box d containing a facet f, whose
closest interfering facet g shares a common edge e with f. See Figure 4.6. We seek
to obtain a lower bound on the distance between the two facets f and g, where
for distance paths are restricted to outside of protected boxes but inside I(ex(B)).
We bound this distance by a constant times oh(B). Then since the box d can only
be facet crowded if both f and ? intersect I(ex(d)), the box sizes are at least a
constant fraction of this distance.
We first make an argument that shows we need only consider a box adjacent
to a box where the corollary to the last lemma holds. As one travels parallel to e
to boxes d farther from the vertex v, but still inside I(ex(B)), the distance within
I(ex(d)) between f and 9 increases. So d is at least half the size of some box
adjacent to B containing f. We have "half" and not equality because of possible
odities in the alignment of the edge with the coordinate (octree) directions leading
to nonuniformity in which boxes are protected. We now consider the first box d1
in the facet phase, such that d' contains f and is adjacent to both B and a box
51
Be, where Be is an edge cone box protected because it contains e or is adjacent to
a box containing e.
We now describe how badly d' can be crowded, that is we seek to describe
the distance between f and g. Because edges are well surrounded by protected
boxes, the distance from e to the boundary of Be shared with d1 is at least a
constant fraction of h(Be). We recall that since Be is an edge protected box, it
is by definition not part of the extended box of d'. Hence, the closest point of
f fl ex(d') to e is on the shared boundary of Be and d1. So we have that the
distance from f fl I(ex(d')) to e A 1(B) is at least kh(Be).
Let e? be the face making sharpest interior angle p with e at v. The ratio of
the size of Be to B is bounded above by a constant times this angle as in Lemma
4. Combining this with the distance results of the last paragraph, the angle at v
between any point of f A I(ex(d')) and face e is at least a constant factor k1 times
the angle p from e to e?, as in Figure 4.6 left.
Let F be an imaginary facet intersecting e and e?. It may be that F is coincident
with either f or g. In Figure 4.6 left F is coincident with f, in Figure 4.6 right F
is distinct from f and g. If e? is a facet, we take F to intersect it forming an edge
that makes angle p with e. Note that F passes between g and f in I(ex(d1)).
The two right triangles at the top of Figure 4.6 left are similar. First, if we
assume that P is less than 7r/4, then the ratio of their sizes is no more than k1/2.
If P > r/4, then it follows from the definition of crowded that h(Be) is at least a
constant fraction of h(B), and the result of the lemma follows immediately from the
fact that e is well centered in Be and f and g make angle at least c at e. Second,
we let P2, be the angle between Ffle2 and the half plane of g with boundary e. Note
52
c
Ye
e2
Yj
Figure 4.6: The distance between two superfacets of an edge inside an extended
vertex box I(ex(B)) but outside an edge box is at least kah(B).
that the larger triangle's leg between F and g is at least k2 Pg' h(B). Combining
these two results shows that the smaller triangle's leg between F and g is at least
k12k2 Pg'?h(B)
In Figure 4.6 the smaller triangle's leg between f and 9 is completely inside Be
and makes a right angle with bnal(B), so that the the distance from gfll(ex(d1))
to Fn 1(ex(d1)) is at least the length of the smaller triangle's leg. If the orientation
of Be with respect to f and g is such that this leg is not entirely inside Be, then this
distance may be at most a constant k3 smaller than the leg's length, because the
facet that e is well centered in Be implies that the angle g can make with ?I(Be) is
bounded below by a constant. Hence dist(gfll(ex(d1)), Ffl1(ex(d1))) > k4.p;.h(B)
for constant k4 = k1k22k3 Similarly for f (recall that Figure 4.6 left is the case where
F is coincident with f, so Pf' happens to be zero). We will prove below that at
least one of or is at least a constant factor times an interior angle of P.
Assume that Pig is at least ka for some constant k. Substituting this in the above
53
yields dist(g A I(ex(d1)), F A I(ex(d1))) > k5 a
Since F passes between f and g, we have a reverse triangle inequality of dis-
tances between sets, that is dist(ffll(ex(d1)),gAI(ex(d1))) > dist(fAI(ex(d1)), PA
I(ex(d1))) + dist(g A I(ex(d1)), PA I(ex(d1))). By the last paragraph, one of these
two distance summands must be at least k5 a h(B). Hence
dist(f A I(ex(d1)),g A I(ex(d'))) > k5 a
Recall our assumption that d' was split because 9 interfered with f. In partic-
ular, f and g are contained in 1(ex(d')), so h(d') is at least a constant fraction of
the distance between f and g in 1(ex(J1)), that is h(d') > kk5 a h(B). From the
definition of facet crowded, as soon as box d' and its descendants are split down
to a size that is a constant fraction of this distance, the descendants are not facet
crowded. Hence d' `s descendants are at least a constant fraction of this distance,
and by construction so is d, b?, and b. That is, for some constants k and k'
h(b) > k'.h(d') > k.a.h(B).
The proof is complete except for showing that one of or Pg' is within a
constant factor of an interior angle of P. We first show that the angle ?1, is an
angle interior to P, but not necessarily an "interior angle" between P faces as
defined in Section 3. Similarly for P9,. If p; + p?' > 7r/4, then the proof is trivial.
Otherwise, we may draw a cone C with apex at v, axis along e, and P A e2 on its
boundary. See Figure 4.6 right for a spherical cross section of this cone, which is
nearly a "top" view of Figure 4.6 left. Since by assumption the P face making the
smallest interior angle with e is e?, the only P faces visible to e within this cone
are f, g, and e?. The angle at v between P A e? and the intersection of this cone
54
with f is an interior angle. But this angle is just P11 Similarly for g.
We finally show that one of P11 or Pg1 is within a constant factor of an interior
angle between P faces. If e? is a facet that meets f at an edge, then the angle P,1
between e2 n F and f fl C may be arbitrarily small, and not bounded by alpha.
We show that even if e2 is a facet meeting f at an edge and also g at an edge, both
angles P11 and P21 may not be small. See Figure 4.6 right. The proof is as follows.
If one of the angles is greater than half of P, there is nothing to prove. Otherwise
they are both less than half P, but then they are a constant fraction of the angle ?j
between e2 n F and e? fl f and the angle ? between e? n F and e? A g, respectively.
The sum of these angles is either a true interior angle between the two P edges of
f A e? and g A e2 meeting at v, or if it is not interior, there is a P edge inside the
triangular cone formed by f, g, and e? but outside the circular cone C, such that
this P edge makes an interior angle with f or g smaller than this sum. In either
case, this sum is at least a.
For a final note in this proof, it is surprising that the result depends on ?1 and
that there is no explicit dependence on ?, where P1 = Pf' + Pg1 What is hidden
is the close relationship between P' and P, which may be shown by the following
mental exercise. For simplicity assume that e? is a subedge of f, i.e. F =
Fix the angle made by the facets f and g at e at a value less than 7r/4. Then P1
depends linearly on P. However, if P is fixed, then P1 depends linearly on the angle
between f and g at e. That is, provided all angles are small compared to 7r/4, P1
is bounded above and below by the product of P and the angle between f and g
at e. In this sense, P1 acts as a lower bound on P, and the explicit dependence on
P is subsumed by dependence on P1. It is important to realize that P1 is (within a
55
constant factor) of a true interior angle, so that the box size has dependence on a,
and not a2. I
This sequence of lemmas concludes the proof of Theorem 2. I
4.2.2 Edge cone boxes
Theorem 3 Let B be a box with edge cone apex v. For any box b adjacent to B,
where k is a constant.
h(b) >k.a.h(B),
Proof. We again consider when b was last split. We know that B was not protected
in the vertex centering phase, since it is an edge protected box.
if b was last split in the vertex crowded phase, regardless of whether B is
protected or not, the balance condition gives us h(b) > h(B)/2. Similarly if b was
split in the edge crowded phase.
So we are left to consider if b was last split in the facet crowded phase. Let f
and g be the two superfacets of e. Let I be the union of the embeddings of all of
the boxes protected because of e or its vertices. Then since e is well centered in
its protected boxes, the distance between a point of f fl ai and g fl ai is at least
k p h(B), where P is the angle between f and g at e. We may now repeat the
argument of Lemma 3, with facets f and g meeting at e taking the place of edge f
and face g meeting at v. That is, a rewording of the proof that edge boxes adjacent
to a vertex box have the correct size proves that facet boxes adjacent to an edge
box have the right size. I
56
4.2.3 Facet Cone Boxes
Theorem 4 Let B be a box with facet cone apex F. For any box b adjacent to B,
h(b) > h(B)/2.
Proof. We know that B was not protected in the vertex and edge protecting phases.
If b was last split in any phase, the balance condition gives us h(b) > h(B)/2. I
If a box is not adjacent to a protected box, then it is empty, and the balance
condition provides us with h(b) > h(B)/2. We may now summanze.
Theorem 5 Let B be any box of the octree whose embedding contains a point of
P, and b an adjacent box. Then
where k is a constant.
4.3 Warping
h(b) > k a
All of the boxes we will now consider are leaf boxes of the octree. If a box is
adjacent to boxes that are smaller, then we replace large faces of the box with the
smaller faces of the adjacent boxes where they intersect. For example, if a box
is completely surrounded by boxes one-half its size, then each of the six sides of
the box is now four square facets with a common vertex and shared edges, and no
longer strictly speaking a cube. We warp the faces of a box away from faces P
when necessary to achieve good aspect ratio when we triangulate. We warp box
faces that are close to P-faces, where the definition of "close" is below.
The warping consists of moving box vertices. We consider adjacent boxes to
share faces of their common boundary, so that when we warp a box face, there is
57
Figure 4.7: The special case of sharing a vertex with two boxes.
a change in all of the adjacent boxes that contain that face. The warping rules are
selected so that no inconsistencies arise.
Note that a protected vertex box (after the centering step described in Section
4.1.3) is always adjacent to boxes the same size or smaller. This means that,
without loss of generality, we do not have to triangulate and warp facets of a
protected vertex box; instead we can triangulate and warp the facets of the adjacent
boxes.
A special note is needed in the case of protected edge boxes when the superfacets
of the edge make an angle close to 3600. If an edge box B is adjacent to two boxes
b and c, such that b and c are not adjacent to each other, but 1(b) n 1(c) is non-
empty, then there may be two duplicate copies of a face that should be shared by
one face of 1(B).
See Figure 4.7. In this case, we retain the two distinct copies of the same face.
Each face is warped separately for b and c. The portion of the face for b is retained
during the triangulation of P A b but ignored when considering P A c. Similarly for
The purpose of warping is to make sure no P face comes too close to a box
58
face. Note that P vertices are already guaranteed to be well separated from box
faces because of the vertex centering procedure described in Section 4.1.3. Thus,
we only need to warp for P edges and P faces.
The warping is done in two passes. The first warping pass is for box edges close
to P edges. We say that a box edge E is close to a P edge F if the distance from
E to F is less than h/8. llere, h is the size of the smallest box sharing the box
edge E. Note that the boxes containing E are all edge-protected, since they are
all adjacent to a box containing F. This means that the sizes of these boxes are
all within a factor of two from each other.
For every close edge E, we move it away from F as follows. We move every
point on the edge, including the endpoints, by distance h/8 in the direction that is
orthogonal to E and F, oriented away from F. (If E and F happen to be parallel,
we move the points of E in the direction orthogonal to E and coplanar with F.)
Even though the boxes after warping have complicated shapes, we still let the
"size" of a box b be its prewarped edge length, and we continue to use notation
h(b) for this size.
The second warping pass is for box vertices close to P facets. We say that a
box vertex v is close to a P facet F if the distance between them is at most h/16,
where h is the size of the smallest box containing v. For such a vertex, we move it
distance h/16 away from the facet in the direction orthogonal to the facet.
In general, because the octree algorithm separates cones, there is never an
inconsistency in these rules that could cause a vertex or edge to be close to two
P facets. There is one exception to this, namely, the case of a protected edge box
containing two facets that make a small interior or exterior angle. If the two facets
59
make a reflex angle, and are both close to the same vertex, then we are in the
special sharing case described before: there are actually two copies of the vertex,
and each is warped separately for the appropriate facet. See Figure 4.7. If the
facets make an acute angle, and are both close to a vertex of the protected edge
box, then we project the vertex orthogonally onto the plane halfway between the
facets (i.e., the plane that bisects the edge cone).
Whenever we warp a vertex in the second warping pass, we carry its adjacent
edges with it. Thus, box edges continue to be straight segments between box
vertices during the warping. The box facets, however, are no longer well defined
because the bounding edges are no longer coplanar.
We make a new definition as follows. Let an 1 vertex be the intersection of a
P edge with a box facet, or the intersection of a P facet with a box edge, or a box
vertex.
In the last paragraph we referred to the intersection of a P edge with a box
facet; the box facet is no longer well defined after warping. Accordingly, we let the
1 vertex be the intersection of the edge with the prewarped box facet.
Lemma 6 After the warping, no two distinct I vertices x, y are closer than min(t, h/16).
Here h is the size of the smallest box containing x, y. Variable t is defined only
for the case when when x, y are on two facets making a sharp interior angle in a
protected edge box, or when one is a box vertex in between two such facets. In this
case, t = k1ha where a is the angle between the facets, and k1 is a constant.
Proof. This proof has numerous cases, all similar. For the sake of brevity, we
prove the lemma in two cases only.
We consider the case that x and y both arise as the intersection of a P facet F
F
60
E1
E2
Figure 4.8: One case in the proof of Lemma 6.
with box edges E1, E2. (This is one of the more complicated cases of the lemma.)
If E1, E2 are disjoint edges then the lemma is clear. Points x, y before the first
warping pass are separated by at least distance h if they lie on disjoint edges. After
the first warping pass x, y must be separated by at least 3h/4. After the second
warping pass, x, y must be separated by at least 5h/8.
Next, consider the other case that E1, E2 have a common box vertex v. (In the
third case that E1, E2 are the same edge, then x, y are not distinct.) This case has
two subcases, namely, v is not considered close to F during the second warping
pass, or that v is considered close. This case is illustrated in Figure 4.8.
First, observe that in either subcase, the angle /xvy is at least 600; this angle
is originally 900, and the warping steps cannot reduce it more than 300.
Now, if v is not close to F during the second warping pass, this means that the
altitude of triangle xvy (from v to the xy leg) is at least h/16, which means by the
angle bound in the last paragraph that the base of this triangle, that is, the xy
leg, is at least h/16 long.
The same argument holds if v is close to F
Another case of the lemma is the case of a protected edge box, where x is the
61
p
x
:2
Figure 4.9: Another case in the proof of Lemma 6.
intersection of a P facet F with a box edge E, and y is a box vertex. First, we
know that, as a result of the first warping pass, y has distance at least h/8 from the
P edge E in the protected edge box, where h is the size of the protected edge box,
or half that size (if the protected edge box is adjacent to a box half its size). Now,
suppose a point on F is distance less than h/16 from y after the first warping pass.
Then we either move y distance h/16 from F, or distance t, where t is half the
distance from zi to Z2. Here, Zi, Z2 are the endpoints of a path though y orthogonal
to the bisector of the two P facets in the box. See Figure 4.9. But since y after
the first warping pass is distance at least h/8 from E, the same holds for Zi, Z2.
This means that the distance between zi and Z2 is at least 2 sin(a/2) h/8, and
the distance from y after the second warping pass is at least tan(a/2) h/8. This
formula has a lower bound of k1ha. I
4.4 Two dimensional triangulation
After warping, we triangulate facets of boxes. This is a preliminary step to the three
dimensional triangulation algorithm of the next section. Suppose two boxes b, ?? are
adjacent, and suppose their intersection is a facet S. if we assume that h(b) >
62
then 5 is a facet of b'. Under these circumstances, we always triangulate 5 by
considering it a facet of b', not b. (If b, b' are the same size, we break ties arbitrarily).
This means that we can always assume that when triangulating a square 5, there
is no box vertex in the interior of 5. There are four box vertices at the corners of
5, and there may be additional box vertices at points along edges of 5.
We now describe the triangulation of a box facet 5. Note that after warping,
the points of 5 are not necessarily coplanar, although the points are nearly coplanar
because the warping distances are small.
The triangulation of a box facet 5 breaks down into four cases, and cases C
and D have two subcases. The cases depend on which P faces pass through the
facet. Note that we only have to triangulate the portion of 5 interior to P, which
is bounded by a closed path of line segments. Call this path of line segments
?. We can also think of ? as a "perturbed polygon" (since the vertices are nearly
coplanar). The six subcases are illustrated in Figure 4.10; the regions 7r are shaded.
The points on the boundary of 7r are all 1 vertices.
Case A, no P faces pass through the facet 5. In this case, we put in a new point
v at the center of the prewarped box facet, and we connect v to every segment along
the boundary of 5.
Case B, two P facets F, C and their common edge E pass through facet 5. Let
v be the 1 vertex where E passes through the prewarped version of facet 5. Then
we connect v to all segments on the boundary of r, triangulating ?.
Case C, one P facet F passes through (two edges of) the facet 5. Let x, y be
the two points where F passes through two edges of 5, and let v be the midpoint
of x, y. Let c be the center of 5 (before warping). Let h be the edge length of 5
63
(before warping).
Subease Ci, v is within h/4 of c. Then we connect v to every segment along
the boundary of ?. This divides the region into triangles.
Subcase C2, v is further than h/4 from c. Then we connect c to every seg-
ment along the boundary of r. This divides polygonal region ? into triangles
and quadrilaterals. Then we insert a diagonal (arbitrarily) into each quadrilateral,
triangulating ?.
Case D, two P facets F, C pass through facet S, but no P edge. Let x, y and
x', y? be points where F, & respectively cross through the boundary of S. Let v, v'
be the midpoints of these segments. Let c be the center of 5 (before warping). Let
h be the edge length of 5 (before warping).
Subcase Di, v or v1is within h/4 of c. Say v is within h/4 of c. Then we
connect v to every segment along the boundary of ?. This divides ir into triangles.
Subcase D2, v and v' are both further than h/4 from c. Then we connect
c to every segment along the boundary of ?. This divides 7r into triangles and
quadrilaterals. Then we insert a diagonal (arbitrarily) into each quadrilateral,
triangulating ?.
We claim that this triangulation has the following property.
Lemma 7 Let 5 be a facet of a box b, triangulated with the above algorithm. Let
h be the size of 5, and let h' be the size of the smallest box adjacent to 5. Let T be
a triangle in this triangulation. Jf two facets pass through 5, let t be as in the last
lemma. Then the radius of the largest inscribed circle in T is at least k2 min(h', t),
where k2 is a constant.
Proof. Standard geometry tells us the radius of the largest inscribed circle in a
64
B
v
C2			DlV!			D2
CI
v
Figure 4.10: The six cases for triangulating a box facet.
triangle T = Axyz is bounded below by
k3. xyl min(I/zxyl, Izzyxi),
where k3 is a constant. This fact is used below.
This proof breaks down into cases. First, we take the case that T has three
vertices pqv, such that p, q are vertices on the boundary of ?, and v is as in the
above cases. This is the type of triangle we get in all cases except C2 and D2.
First, observe that the distance from p to q is at least min(h', t/2); this follows
from the preceding lemma. Thus, we only need a constant lower bound on Zvpq
and Zvqp. But these angle bounds are immediate, because v is near the center of
the facet, and p and q are along the boundary of the facet.
A slightly different analysis is needed for the triangle vx'y' arising in Case Dl,
where x? and yt are vertices arising from the intersection of the P facet C and S.
We know from Lemma 6 that the distance between two I vertices, one from F n 5
65
e1
Figure 4.11: Triangulating a quadrilateral in Case C2 (closeup)
and one from G n 5, is at least t. In fact, any point of F n 5 is at distance at least
kt from G n 5, from the fact that the common edge between F and G is far from
5 after warping. To be technically correct, since the vertices of 5 are not actually
coplanar, F n 5 refers to any triangulation of 5 intersected with F, and similarly
for G. Hence the distance from v to x'is at least kt, which gives us an edge length,
and the distance from v to 0 n 5 is also at least kt, which gives us angle bounds.
A different kind of triangle arises from splitting quadrilaterals in Cases C2 and
D2. For each of the quadrilaterals in these cases, we can see that the two opposite
edges el and e? have length at least k4 min(h', t/2). See Figure 4.11. For el,
this follows from the same argument as in the last paragraph. For e?, this follows
because facet F in 5 is well separated from the center c. Therefore, the length of
e2 is at least a constant fraction of the length of e?. Similarly, it is clear that edges
e? and e? have length at least k4 min(h',t/2).
Similarly, we can show that the angles made by edges e?, e? with el, e2 are
bounded below by constants. Again, this follows because c is centered, whereas
el, e2 are away from the center.
This is enough to show that the inscribed circles of the triangles obtained as a
66
result of inserting a diagonal are at least a constant multiple of min(h', t/2). I
Theorem 6 Let B be a box. For any triangle T on a facet of B, let r be the radius
of the largest inscribed circle. Then r > k a h(B), where k is a constant.
Proof. This is an immediate consequence of the preceding lemma and Theorem 5.
in particular, h' in the previous lemma is at least kah(B) by Theorem 5. Moreover,
tin the previous lemma was shown to be at least k a h(B) in the proof of Lemma
6. I
General results have been obtained for two-dimensional triangulations with
guaranteed inscribed-circle radius bounds. Bern, Edelsbrunner, Eppstein, Mitchell,
and Tan (1991j have a result concerning optimal two-dimensional triangulation of
polygons, assuming new points cannot be introduced. Bern, Dobkin, and Eppstein
[1991J have similar results in the case that new points can be introduced. We
have not been able to incorporate these algorithms in the present version of our
triangulation algorithm because we need to introduce new points, but the new
points have to be introduced in a carefully controlled fashion.
4.5 Three dimensional triangulation
We now describe how to form tetrahedra from the triangles in the last section. We
triangulate on a box by box basis. The details of how we triangulate depends on a
case analysis of what is contained in a box, and how it intersects the box. However,
the general principle is to find a central vertex, and then form one tetrahedron or
prism for each triangle in the box by taking the convex hull of this vertex and the
triangle. The organization of the argument is very similar to the two-dimensional
triangulation in the last section.
67
4.5.1 Empty box tetrahedra
We first describe how to triangulate in three dimensions a box b containing no P
faces. We take as a central vertex v the centroid of b. (The prewarping centroid
may be used). For every triangle on the surface of the box we form a tetrahedron
by taking the convex hull of v and the triangle.
4.5.2 Vertex box tetrahedra
For a box b containing a vertex of P, the three-dimensional triangulation is partic-
ularly easy. We take as the central vertex v the P vertex itself. We form tetrahedra
by taking the convex hull of v each triangle on the surface of 6. The vertex center-
ing phase ensures that the distance from the central vertex to the triangles is at
least a constant times the box size.
4.5.3 Edge box tetrahedra
We now describe the next case, that of triangulating in three dimensions the pro?
tected edge boxes that contain edges. Let edge E be the edge contained in box 6.
We let v be the midpoint of E n 1(6). We now connect v to every triangle on the
surface of 6, creating tetrahedra.
4.5.4 Facet box tetrahedra
Consider a box 6 containing one P facet F. We introduce a central vertex v at the
centroid of the prewarped box. If v is within distance h(b)/4 to F, then we move
v to the closest point of F to v.
We now have three cases. The first case is if v lies in the interior of P, but not
on F. Here we form tetrahedra by taking the convex hull of every triangle of F
68
and v. This requires the creation of a triangulation on the intersection of F with
b, which is a polygonal region. This polygonal region is nearly convex (it would
be convex if the triangles on each the surface of the box were coplanar). To form
a triangulation of F n 1(b), we pick a centrally located point on the F n 1(b) and
connect it to all the edges. Using similar arguments as in the previous section,
it is not hard to show that the largest inscribed radius of any triangle in this tw?
dimensional triangulation is at least a constant fraction of ??, where ?? is the size
of the smallest box adjacent to b.
The second case is that v lies on F. We take the central vertex v, and proceed as
if it were a vertex of P. That is, for each triangle of 1(b), we generate a tetrahedron
by taking the convex hull of it and v.
The third case is if v lies outside of P. We form tetrahedra by taking the
convex hull of v and every triangle of 1(b). These tetrahedra are clipped by F into
"prisms" with nonparallel sides, and a triangular top and bottom. Zero, one or
two of the vertices of the top may also be vertices of the bottom. If two vertices
are shared, then the prism is a tetrahedron. If one vertex is shared, the prism is
split into two tetrahedra by introducing a facet between the shared vertex, and
one distinct vertex of the top and one distinct vertex of the bottom. If no vertices
are shared between the top and bottom, the prism is split into three tetrahedra
by introducing two facets. One facet is between two of the top vertices and the
opposite bottom vertex. The other facet is between the opposite bottom vertex,
one of the other bottom vertices, and the opposite top vertex.
69
4.5.5 Two-facet box tetrahedra
Perhaps the most complicated case is the case of a box b containing two facets
F1, F2 adjacent to an edge E, but the edge itself is not in the box. Let v be the
centroid of box b (prewarped). If v is within distance h(b)/4 from F1 or F2, then
we move it to the closest point on either of these facets (break ties arbitrarily).
Say F1 is closer.
There are several possible arrangements for v, F1, F2. For example, v could be
on F1, or it could be between F1, F2, or it could be outside P, with F1 facing away
from v and F2 facing towards P
For each facet Fi, (i = 1,2) we triangulate Fj A 1(b) provided that Fi does not
contain v and faces v. This means that either one or both of Fj is triangulated.
We use the same algorithm described in the facet box case for triangulating the
facet.
Now we connect v to all the triangles on the boundary of b that are interior to
P. We also connect v to all the triangles of F1, F2.
If v is contained in F1, or if v is between F1, F2, then we have now divided the
box into tetrahedra. If v is outside P with F1 facing away, then by connecting v
to the triangles of F2 we have divided the interior of P into tetrahedra and prisms.
We now split the prisms into tetrahedra as in the facet box case.
4.6 Aspect ratio of tetrahedra in AocT
Our triangulation is optimal in two respects. First, the maximum aspect ratio
among tetrahedra of our triangulation is optimal up to a constant multiple. Second,
compared to all other triangulations of fixed aspect ratio, our triangulation has the
70
minimum number of tetrahedra up to a constant multiple. In this section and the
next section we establish the first optimality property, focusing on the optimality
of the aspect ratio.
Recall that the tetrahedra we form are either the convex hull of a central vertex
and a triangle that does not contain the central vertex on the surface of a cube or
a facet of P, or else a portion of a prism that is divided into three tetrahedra. We
now want to argue about the aspect ratio of all tetrahedra arising in the algorithm
of the previous section.
There are three types of tetrahedra arising in the triangulation. A Type A
tetrahedron has one vertex v centrally located in the box, and the other vertices
on a box face. This type of tetrahedron arises in the cases of triangulating a box
with a vertex, a box with an edge, or a box with a facet in which the vertex v
in the last section lies near the center of the box, is in P, and is not close to the
second facet (if the box has two facets). Also some of the tetrahedra arising from
the case of two facets in a box and v on one of the facets are of this type.
A Type B tetrahedron arises only in the two-facet case in the last section, when
the vertex v is on one facet, and the triangle it is joined to lies on the other P
facet. These tetrahedra arise in cases analogous to one of the triangles in case Di
of Figure 4.10.
A Type C tetrahedron arises from a vertex v in the last section outside P. This
happens only with boxes with one or two facets and no edges. This causes P A b
to be divided into prisms, and then each prism is divided into tetrahedra. These
tetrahedra arise in cases analogous to triangles in cases C2 and D2 in Figure 4.10.
Intuitively, all three types of tetrahedra have aspect ratio at most 1/0, because
71
they involve a base (the base being a triangle) whose inscribed radius is between
kiah(b) and k2h(b), and the apex is well-centered over the base, and is distance
between k3oh(b) and k4h(b) from the base. Here, ..... . , k4 are constants.
The aspect ratio computations are very technical and involved. We first prove
the aspect ratio bounds for Type A tetrahedra using a construction of an inscribed
sphere. The proofs for Type B and C tetrahedra follow from a different analysis
of the same construction. A Type A tetrahedron has a base size between kioh(b)
and k2h(b), and has a height at least k5h(b).
Consider a particular box b and its central vertex v, and a Type A tetrahedron
t3 in this box. In this section we let subscripts denote dimension, so that t2 is one
of the triangles on the surface of P A b, and t3 is the tetrahedron formed by the
convex hull of t2 and v.
Lemma 8 The smallest containing sphere, S3, of t3 has radius R3 at most kh(b),
where k is a constant.
Proof. The smallest containing sphere of t3 has radius less than the radius of a
sphere that encloses the entire box. If the box was unchanged by the warping step,
we see that a sphere with radius equal to one half a diagonal of the box will do. If
the box was warped, then it may be enclosed in a concentric box with side lengths
5/4 times that of the original box by the definition of "close." Thus we can take
k=5$3/4. I
This concludes the analysis of the smallest sphere containing t3. Next, we
analyze the largest inscribed sphere of t3 by actually constructing an inscribed
sphere. Consider the triangle t2 on ?(P A b). We form a new curved triangle U2
on the surface of a sphere as follows. Let a be the closest point of t2 to v. It does
v
72
a
Figure 4.12: The flat (t3) and curved (u3) tetrahedra.
not matter if a is a vertex of t2 or not, although in the figures of this section a is a
vertex. The intersection of t3 and the sphere centered at v with radius dist(v, a) is
the curved triangle U2. We let U3 be the convex hull of U2 and v. Note that U3 is
contained in t3, and has three flat facets, each one contained in a facet of t3 that
intersects v. This is illustrated in Figure 4.12. We will construct a sphere inscribed
in U3 instead of t3.
Lemma 9 The distance from v to its base t2 is at least h(b)/8.
Proof. This is because, by assumption for the case, t3 is a Type A tetrahedron.
Therefore, v is centered in b. I
Lemma 10 Let U2 be defined as in the last paragraph. Then the radius of the
largest inscribed curved circle in U2 is a constant fraction of the radius of the
largest inscribed circle in t2.
Proof. Let C be the largest inscribed circle in t2. We project this circle onto sphere
73
53. This yields an ellipse, whose minor axis is the product of the original inscribed
radius and a factor depending on the angle between t2 and a tangent plane to
U2. We claim that the angles between v and any edge of t2 are bounded below by
constants. This is because, before warping, the vertices of t2 were coplanar with a
box facet, and v is at least distance h/4 from that facet, where h is the edge length
of the box. After t2 is warped, since the warping distances are small, this can only
change its inclination with respect to the original facet by a small constant angle
that is arrived at via a case analysis (which we omit).
This same angle bound applies to the angle between t2 and the tangent plane.
Lemma 11 (Corollary) If v is within kah(b) of the plane of t2, then the radius
of the largest inscribed curved circle in U2 is a constant times a times the radius
of the largest inscribed circle in t2.
Proof. The key observation is that v is still within a constant of the box size to the
farthest point of the triangle. The angle between v and an edge of t2 is bounded
below by a term linear in the ratio of the distances from v to the plane of t2 and
from v to the farthest point of t2. In the lemma this was a constant, here in the
corollary it is bounded by a. This will be used only for the analysis of Type B
tetrahedra. I
In a series of lemmas we derive a lower bound for the radius of the largest
sphere inscribed in t3 by the radius of the largest curved circle inscribed in U2.
Let 53 denote the largest sphere inscribed in U3, r? its radius, and c? its center.
Similarly, we let ?2 denote the largest circle inscribed in U2, r2 its radius, and C2
its center.
74
v
Figure 4.13: A sphere tangent to three facets of a tetrahedron.
Lemma 12 The sphere 53 is tangent to U3 at points equidistant from v.
Proof. Pick a facet f of U3 adjacent to v. Consider the triangle formed by v,
the point of tangency between 53 and f, and the center of 53. Note that this is a
right triangle. Observe that the length of the hypotenuse and the length of one leg
depends only on the position of 53 and not on on which facet was chosen. Therefore,
the lengths of the other leg (the leg between v and the points of tangency) is the
same regardless of which facets is selected. See Figure 4.13. This theorem applies
to arbitrary tetrahedra with any vertex chosen for v. I
Lemma 13 There is a sphere, 5?3, that is tangent to U3 at exactly the same three
points as ?2 is tangent to U2.
Proof. We call these points of tangency a, b, e. We construct 513 Note that 5?3 is
not contained in U3 or t3, but will serve as a guide for finding an inscribed sphere.
75
We first determine uniquely the center, c', of 313 Consider the line vc?. Now
d is the unique point on vc2 on the opposite side of C2 as v, such that Zvac' is
right. Additionally, since 7c? is perpendicular to the boundary of U2, ad is also
perpendicular to the facet, f, of U3 containing v and a in the plane tangent to the
vertex sphere at a. Hence, da is perpendicular to the facet f. Similarly for b and
e. The hypotenuse-leg theorem gives us congruent triangles, and we see that the
distances from c' to each of a, b, e are equal. So, c, is indeed the center of a sphere
with the desired properties, 3'3. I
Lemma 14 (Corollary) Any sphere tangent to the three faces of U3 containing
v has center on Thd, and points of tangency on 7a, vb and ve.
We now seek to find a sphere inscribed in U3, which is almost as large as the
largest inscribed sphere. Let d be the distance from v to the flat triangle t2. We
see that we can chose an inscribed sphere tangent to the three facets of U3 meeting
at v, with radius r, and distance from center to v equal to t? such that
See Figure 4.14 for a picture of the two dimensional case of this construction.
Let 0 be the angle between Thc' and va. It is obvious from the figure that
= csc(O), so that
d
r1 +csc(O)
We note that the largest inscribed sphere has radius at least that of this in-
scribed sphere, r3 > i:. Note that 0 < 0 < ir/2. Hence l+Jc(O) = H$8j?S?fl8O+1 >
sin(0)/2> kO. Hence
> kdO.
v.
76
Figure 4.14: Constructing a large inscribed sphere of U3.
Also from Lemma 10 and Figure 4.14,
Recall the result of Theorem 6,
0 = r2/d.
> kah(b).
Hence r? > kah(b), and we have proved the following theorem for Type A tetra-
hedra.
Theorem 7 (Aspect ratio of AocT) Let B be a box. For any tetrahedron t
inside B, the aspect ratio oft is at most k/a, where k is a true constant, and a is
the sharpest interior angle of P
This theorem is the main result for this section, which the preceding lemmas prove
for Type A tetrahedra. We now wish to establish this theorem for the remaining
types of tetrahedra, Types B and C.
We consider Type B tetrahedra briefly because the analysis is so similar to
that for Type A tetrahedra. We perform the same construction as for Type A
tetrahedra, defining d and 0 as before. Consider a Type B tetrahedron, the convex
hull of vertex v on P facet F and triangle T on the other P facet G in box b. If
77
the distance from v to G is at least kh(b), then the analysis of Type A tetrahedra
suffices. Otherwise, we have d> kah(b) as in Lemma 6. As before, the construction
gives r3 > kr2. Because d may be closer to the plane of T, the result of Lemma 10
does not hold, and instead we use the result of Lemma 11. This is sufficient because
we have a stronger bound on the shortest altitude of T than when considering Type
A tetrahedra, namely alt(T) > kh(b).
We now consider Type C tetrahedra. These tetrahedra arise from prisms. Since
the central vertex v is well centered in the box, and no P facet is close to it (and no
P edges are in the box), the angle ? between the shortest segment from a triangle T
to v and the normal to the plane of T is no more than a constant slightly larger than
ir/4, in contrast with the Type B tetrahedra. We will make several constructions
as for Type A tetrahedra, only we will chose different vertices for "v,,. However,
our choice of "v" will be related to the box's central vertex, so that Lemma 10 will
hold for Type C tetrahedra. Hence we have r? > kr2 > kah(b), and we now need
only specify the constructions.
Consider the various types of prisms as in Section 4.5.4, having a triangular top
and bottom, and one, two, or three quadrilateral sides depending on how many
vertices of the top are shared by the bottom. Prisms formed by clipping the convex
hull of v and a P triangle have zero shared top and bottom vertices, and will be
considered last. Otherwise, if two vertices of the top of a prism are shared with the
bottom, then the prism is a tetrahedron. We take the top to be a box triangle. We
now perform the same construction as for Type A tetrahedra, with the non-shared
bottom vertex as "v".
If one vertex of the top is shared with the bottom, then the prism is split into
78
two tetrahedra. Define the top as in the last type of prism. The bottom is a
projection of the first onto a P facet, with center of projection v far from the top
and bottom. Recalling that ? is bounded by a constant, the bottom altitude is at
least a constant fraction of the top altitude, kch(b). For the tetrahedron that is
the convex hull of the top and a non shared bottom vertex, we chose "v,, to be the
bottom vertex. For the tetrahedron that is the convex hull of the bottom and a
non shared top vertex, we chose "v" to be the top vertex.
If no top vertex is shared with the bottom, then the prism is split into three
tetrahedra. We let the top be the box or P triangle. As in the previous case,
the bottom triangle has altitude at least a constant fraction of the top altitude,
kah(b). The analysis of the one shared vertex case is sufficient for the tetrahedra
with the top as a facet and the tetrahedra with the bottom as a facet.
For the third tetrahedron, the one that is the convex hull of an edge of the
top and an edge of the bottom, we chose the bottom vertex opposite the two top
vertices as "v,,. The triangle of the tetrahedron opposite t in this subcase arises
from introducing a diagonal on one of the quadrilateral sides of the prism. If both
the top and bottom are subsets of P facets, or in general if the box is a protected
edge box not containing an edge, the P edge is at least kh(b) away from t, and
hence t has altitude at least kah(b). If the top is a subset of a box facet, then t has
altitude at least kh(b). In both the case where the top is a subset of a P face and in
the case where the bottom is a subset of a box face, the angle between the planes
of the top and bottom is between 0 and a constant slightly larger than 3ir/4, by the
fact that neither plane is close to the central vertex, the central vertex is outside
of P, and both planes pass through the box. For the definition of angle between
79
the planes, we chose the normal to the bottom pointing towards the interior of the
prism, and the normal to the top pointing towards the exterior of the prism. In
particular, the angle between the shortest segment from t to "v,, and the normal
to t is bounded by a constant, so that again Lemma 10 holds.
4.7
Aspect ratio tetrahedra in A*
In the last section we established an upper bound on the aspect ratio of our trian-
gulation. In this section we show that we are a constant factor from optimal, by
analyzing the aspect ratio of the optimal triangulation.
We wish to show that i/a acts as a lower bound on the aspect ratio of any
triangulation of P, where a is the smallest interior angle of P. Note that if a >
then in any triangulation of P there is trivially a tetrahedron with aspect ratio at
least k/a for some constant k, since aspect ratios are always greater than 1. Hence,
we need only prove the following theorem
Theorem 8 Let P have smallest interior angle a. If a < 7r/2, any three dimen-
sional triangulation of P has a tetrahedron with aspect ratio at least 1/ sin(a).
Before starting the proof, we state the following simple geometric lemma, whose
proof is omitted.
Lemma 15 Let Q be a plane in lR3, and let x be a point on Q. Let rl,r2 be two
rays starting from point x such that plane Q separates them (or one or both ri, r2
may be contained in Q). Then the angle between ri and Q is bounded above by the
angle between ri and r?.
Proof of Theorem 8. We first show that for any three dimensional triangulation
of P, there is some tetrahedron of the triangulation with angle at most a between
80
two edges or an edge and a facet not containing the edge. Consider the faces F
and G of P that intersect at p, with the angle between F and & equal to a. Note
that we assume that F and G can see each other in the sense of Section 3. Recall
that we may always chose F to be an edge. If F makes angle a with any edge, we
let that edge be &, otherwise C is facet. In either case there is some segment in
the relative interior of C along the ray in C making angle a with F. See Figure
4.15. We let Ci be this segment, and note that p is contained in G. Now, let Cm
be a segment joining a point on P and a point on Cj so that Ci, Cm, and F make
a triangle and so that F and Cm make a right angle. Let p' be the intersection of
F and Cm. Let p" be another point on 0m close to p'. Given a three dimensional
triangulation A*, we may chose 0m close enough to p and p" close enough to p',
such that every tetrahedron of A* containing p" also contains p and p'.
Let A* be one of these tetrahedron containing P? Then A* contains P? P'? and
also has an edge, say F1, that is coincident with a subsegment of F. Without loss
of generality, Cm is moved so that P1is a vertex of A* terminating edge F'. Let
CI be the face of A; opposite p1. Then the angle made by CI and F1 is no larger
than a, because CI is on a plane separating F' and Ci as in the last lemma.
We now show that this small angle leads to a large aspect ratio for A* We
first bound the altitude of A; at p', which bounds the radius of the largest sphere
inscribed in A* If C' is a facet then we chose H = C', otherwise we chose Has
the unique facet of A* containing edge C' but not edge F'. That is, we let H be
the facet of A* that intersects P1 at exactly P Then the angle a' between F' and
the plane of H is at most a. It is clear that since tetrahedra are convex, the radius
of the largest sphere inscribed in A* is at most half of the altitude between P' and
81
F
Gj
Figure 4.15: Finding the tetrahedron with an angle at most a.
H, where as before p' is the vertex of F' that is not p. This altitude is sin(a')/IF'I,
where F'l denotes the length of edge ?? Hence we have an upper bound on the
size of the largest sphere inscribed in A* and we now find a lower bound on the
size of the smallest sphere containing A;. Since a containing sphere must have
diameter at least the longest edge, we have that the smallest containing sphere has
radius at least IF'I/2. Combining these two shows that the aspect ratio of A* is
at least 1/sin(a'), which is at least 1/sin(a) as desired. I
Let T be a three dimensional triangulation. Let A(T) be its aspect ratio, defined
to be be the maximum over all tetrahedron t E T, of the aspect ratio of t. We
may now state the theorem concerning the optimality of the aspect ratio of our
triangulation of P.
Theorem 9 (Aspect ratio optimality)
A(AocT) < kA(A*)
82
where k is a true constant, independent of P and a, and AocT is the triangulation
of P arising from our algorithm, and A* is any other triangulation of P.
4.8 Optimality of the cardinality of AocT
In this section we prove that the number of tetrahedra in our triangulation is no
more than a constant factor larger than any other triangulation with a bounded
aspect ratio. The reasoning in this section is as follows. We first prove that
any triangulation with bounded aspect ratio, which we denote A*, has certain
geometric properties concerning how fast the sizes of the tetrahedra can change. We
then derive lower bounds on how small the tetrahedra of our triangulation AoGT
can be. We then compare the two sets of bounds to show that our triangulation is
within a constant factor of optimal.
Conformal. Any triangulation A of P that we consider in this paper must
be conformat to P. This means that the following condition must hold. Let x be
a point on a P face F. Let F' be the lowest dimensional face of A containing
x. Then F' c F. It is easy to see that the triangulation we construct has this
property; vertices of P are covered by vertices of AocT, and edges of P are covered
by edges of AOGT. A more general, but equivalent definition of conformality is the
following lemma, whose proof is obvious from the definition.
Lemma 16 If a triangulation A is conformal to P, then the following condition
is satisfied. Let F, G be two faces of P, and x, y two points on F, 0. Let F', 0' be
the lowest dimensional faces of A containing x and y. Then F' n 0' c F n 0. In
particular, if F, 0 are disjoint, so are F', 0'.
83
In this section there will be a series of positive constants Ci, .2.... depending
only on the maximum aspect ratio A of the triangulation under consideration. For
example, the maximum number of faces incident to a vertex of A* is bounded by
a constant ci.
Characteristic length function. We define f?: P JR to be the character-
istic length function of a triangulation AT. This function is defined as follows. If
x is a point of P, then we define fT(x) to be the longest edge among all tetrahedra
that contain x.
The next group of lemmas focus on a triangulation A* of bounded aspect ratio
and its characteristic length function f*. We will show below that the longest edge
of any tetrahedron containing x is bounded by a constant times the longest edge
of any other tetrahedron containing x, so that any tetrahedron containing x can
be used to bound f*(x) above and below.
Lemma 17 Suppose x, y E P, and suppose the two tetrahedra defining f*(x) and
are T1 and T2 (in particular, x E T1 and y E T2). If T1 and T2 share a
common edge, then
c2?f*(Y)? f*(x) <?c?.f*(y).
Proof. Let e denote the common edge. Note that because of the aspect ratio
bound of A*, the ratio of the longest to shortest side of any simplex in A* is
bounded above and below. Hence
lel < f*(x) < .4 lel.
This also holds for y. 1
84
Wedge. A wedge is a set of tetrahedra, T1,T2,... , Tk, that share a common
vertex.
Note that the bounded aspect ratio of A* bounds the smallest interior angle
of any simplex in A*, so that the number of simplices in a wedge is bounded by a
constant c1
Lemma 18 Given a bounded aspect ratio triangulation A*, two points x, y, and
any two faces F, G containing x, y respectively, define H = F n G. If H # ?, then
Proof. Note that H is a face of A*. Let v be a vertex of H. We consider the
tetrahedra of a maximal wedge with common vertex v. Since no vertex of P is
degenerate, between any two tetrahedra there is a sequence of tetrahedra of the
wedge such that each shares an edge with its successor. Since wedges consist of
a bounded number of tetrahedra, the function f* on any tetrahedron is within a
constant multiple of f* on any another by Lemma 17. This wedge includes F and
G. Finally, there are tetrahedra T1, T2, such that T1 contains x and determines
f*(x) (i.e., the longest edge of T1 is equal to f*(x)), and similarly T2 determines
f*(y). Then T1 and F have a common vertex (since they both contain x) and
therefore the lengths of their edges are within a constant of each other by the same
argument. The same holds for T2 and G. 1
Geodesic distance. We let distp(F, G) denote the geodesic distance in P
between two closed subsets F, C of P. In other words, this is the length of the
shortest path in P between a point of F and a point of G. This distance is always
at least as great as the Euclidean distance.
85
As an example of this notation, we have the following lemma. We omit the proof
of the lemma, which follows from the meaning of aspect ratio for a tetrahedron.
Lemma 19 Let T be a tetrahedron of a bounded aspect ratio triangulation, let x
be a vertex of T, and let F be a subface of T not containing x. Then distp(x, F) >
Similarly, the following lemma follows from elementary trigonometry applied
to bounded aspect ratio.
Lemma 20 Let T be a tetrahedron with bounded aspect ratio. Let F be an edge or
facet of T, and let x be a point on F whose distance from any subface of F is at
least t. Let y be a point lying on any facet G of T, such that G does not contain
F. Then the distance from y to x is at least c?t.
Proof. An illustration of the claim made by this lemmais in Figure 4.16 in the case
that F is a facet. If 0 is a subface of F, then the lemma is trivial (we could take
= 1). Otherwise, let v be the point on the boundary ofF closest to x. (Note that
if F is an edge, v will be one of its endpoints). Consider the triangle Axvy. Let
LFG denote the interior angle between F and 0. Then Zxvy > min(ZFC, ir/2).
But LFG is bounded below by a constant c? depending on the aspect ratio of
the tetrahedron. Moreover, the vx leg of this triangle has length at least t by
assumption. This means that the distance from x to y is at least t sin c?. I
We now wish to establish the following theorem about the distance between
faces in a triangulation. The theorem has two cases, depending on whether the
faces have a common point or not.
G
86
Figure 4.16: An illustration of the conditions stated in Lemma 20.
Theorem 10 Let F,C be two faces of A*. Let H = FnG. Let x,y be two points
such that x E F, y ? C. Then
1. If H $ ? and distp(x, H) > t, then distp(x, y) > c?t.
2. If H = $, then distp(x,y) > ciof*(x).
Proof. This proof is broken down into eight cases, depending on whether we are in
Case 1 or 2 or the theorem, and depending on whether F is a vertex, edge, triangle
or tetrahedron. The cases are labeled Al to D2. Note that if we are in Case 1 of
the theorem, we can assume that H is a proper subset of both F and G. If H is
not a proper subset, say F = H, then distp(x, H) = 0, and the claim is trivial.
Similarly, if G = H then y E H so distp(x, y) > distp(x, H).
In each of the upcoming cases, let p be the shortest path in P from y to x.
Case Al, F is a vertex and H $ ?. This is a trivial case described in the
previous paragraphs, since H must be the same vertex as
Case A2, F is a vertex and H = $. Let T be the first tetrahedron on p such
that T contains F. Then T must be entered on a face F' of T not containing
vertex F (since G does not contain F). Thus, from Lemma 19, we know that
distp(F, F') > c6f*(x).
87
Case Bi, F is an edge and H $ ?. Then G is an edge or a triangle, and H
is a vertex v common to F and G. Note that t is bounded above by f*(x) since
t cannot be longer than the length of F. Let w be the other endpoint of F. We
take two subcases: either x is within distance c7f*(w)/2 of w or it is further than
this distance. If x is within c7f*(w)/2 of w, then we apply Case A2 to get a lower
bound of c6f*(w) on the distance from y to w. This means that the distance from
y to x must be at least c6f*(w)/2, which is at least ciit. (Note that f*(w) >
since any tetrahedron containing x will also contain w.)
The other subcase is that x is distance at least c6f*(w)/2 from w and at least
t from v. Let T be the first tetrahedron on the path p that contains edge F. Say
that tetrahedron is entered at point y' on face ?? Then F' does not contain F.
By Lemma 20 the distance from y' to x is at least
C7 min(c6f*(w)/2, t).
Thus, in all cases, the distance from F to C is at least q?t.
Case B2, F is an edge and H = ?. This is handled with exactly the same
argument as Case Bi, namely, we observe that either x is within C6 f*(w)/2 of an
endpoint w of F, in which case the distance bound can be reduced to A2, or else
x is further than this distance from either endpoint, in which case we use Lemma
20. Again, in this case, we get a distance bound of c12f*(x) on the distance from
x to y.
Case Cl, F is a triangle and H $ ?. Let us take two subcases:
Snbcase Cl (a), H is an edge of triangle F. This means that G is a second facet
containing edge H. Then we first check whether x is within distance min(ci?, 1)t/3
from a point y? on one of the other two edges E1, E2 of triangle F. If so (say y'
88
is on E1), then y' is at least distance 2t/3 from H, so we use Case Bi (applied to
E1 and G) to get a lower bound of 2c12t/3 on the distance from y' to H, and this
yields a lower bound of ci?t/3 on the distance from x to H. Otherwise, x is further
than distance min(ci?, 1)t/3 from all three edges of F. Then we look at the first
tetrahedron T on path p that has F as a facet. This tetrahedron must be entered
by a point y' not on F. Then Lemma 20 tells us that the distance from y' to x is
at least C7 min(ci?, 1)t/3.
Subcase Cl(b), H is a vertex of triangle F. Let E1, E2 be the two edges of
F containing vertex H, and let E3 be the third edge. If x is within distance
min(ci?, 1)t/3 from a point y' on either E1 or E2, we can prove the theorem as in
Subcase C1(a) by applying the argument of Case Bi to E1 and G, or to E2 and G.
Another possibility is that x is within distance c12f*(x)/2 from a point y' E E3, in
which case we can apply Case B2 to E3 and G, since they are disjoint. This gives a
lower bound of c12f*(x) on the distance from y to G. Otherwise, the distance from
x to any of the edges of F is bounded below by either min(ci?, 1)t/3 or c12f*(x)/2.
Let T be the first tetrahedron on p that has F as a facet. Then T is entered from
a point y? not on F. Finally, Lemma 20 gives us a lower bound of the form
C7 min(c12f*(x)/2, min(ci?, 1)t/3)
on the distance from x to y.
Cases Dl--HD2, F is a tetrahedron. These cases are not needed by the upcoming
theorems, so we omit them and leave them to the reader. I
Theorem 11 There exist constants C13 and ci? dependent on A such that for all
x,y Ei P,
f*(x) ? max(q? f*(y),c14 . distp(x,y)).
89
Proof. If x and y are in a common tetrahedron of A*, then the result follows from
Lemma 18.
Otherwise, let xy denote the shortest path in P from x to y. Label the triangles
of A* crossed by xy by Fj,i = 1,2,... ,m. We assume in this proof that xy
intersects each Fi in a single point; the degenerate case of exact alignment by Fi
and xy can be analyzed by considering a slight perturbation to xy. Note any two
consecutive triangles of the sequence are in a common tetrahedron. Let Xj denote
Fj fl xy. By Lemma 18 we know that f*(x) < C5 f*(x1) and f*(xm) < C5
Define
i* =maxti:F?nf\#$1.
If i* = m, then xi and Xrn are in a common wedge, and f*(x1) < C5 f*(xm).
Combining these inequalities shows that f*(x) < Cj3 f*(y). The other case is that
i* <m, and we note that distp(xi, xi*+i) > cio f*(x1) by Theorem 10. Moreover,
distp(x, y) > distp(xi, Xj*+1). Combining these inequalities gives the result that
< Cj4 distp(x,y). I
We now observe the following characterization of our triangulation AocT
Lemma 21 Let x C P be a point lying in octree box b. Then focT(x) > kiah(b),
where k1 is a constant.
Proof. Let T be the tetrahedron of AocT lying in b containing x. Then each
edge of T on the boundary of P A b has length at least k1ah(b) as determined in
previous sections. A case analysis verifies that every tetrahedron we construct in b
has an edge on the boundary of P A b. This includes those tetrahedra arising from
prisms. I
90
We now have a lemma about geodesic distances in boxes. Note that in a
crowded box b, it is possible to have two points x, y in b such that distp(x, y) is
arbitrarily larger than h(b). For example, P might be shaped like a coil inside b.
There are two cases, however, when we can get bounds on geodesic distances.
Lemma 22 Let b be an uncrowded box in the octree, and let x, y be two points in
PAb. Then
distp(x,y) < 2?h(b).
Proof. Recall that in a box, b, P A b is a single component C, so x, y are both in
this component. Either there is a vertex, edge or facet cone in the box, or b c P.
In the latter case, the distance from x to y in b is the same as the distance from x
to y in P, which is bounded by $3h(b).
The other case is that b contains a cone. Let z be a point on the apex of
the cone. Then x and y are both connected to z by a straight-line path in P, so
distp(x,z) < $3h(b) and distp(y,z) < $3h(b). 1
We wish to establish a series of lemmas that prove that the characteristic length
function at a point of A* is bounded by a constant times the size of the octree box
containing that point in our algorithm. The first lemma involves a construction
that is useful for each of the remaining lemmas. The remaining lemmas prove the
result for a point in a protected vertex box, protected edge box, protected facet
box, and finally for an unprotected (empty) box.
Lemma 23 Let b1 be a box split in the vertex, edge or facet phase, and x a point
of P face F in b1. Then there exists a point y of P face G such that
91
1. G does not contain F, and
2. distp(x,y) < 6?h(b1).
Proof. If b1 is split because of the balance condition, let b1 be the crowded box
that caused it to split. There may be a long chain of boxes that split to lead to a
splitting of bi; say the chain is ..... . , b1. In this chain, b1 is crowded, and it causes
b1?1 to split, which causes b1?2 to split, ... which causes b1 to split. Without loss
of generality, we can assume that all of the boxes ..... . , b1?1 are not crowded. For
example, if b1?, 1' < 1 is crowded, then we can replace 1 with 1'.
The sizes of these boxes must be decreasing geometrically, so that h(bj) =
2?--H1h(b1), because otherwise the balance condition would not be enforced from
bj+i to bj. This means that the Euclidean distance from a point in b1 to a point in
b1 is at most 2?h(b). In other words, 1(ex(bj)) c I(ex(bj)) for all j = 1,2,.. .1.
Recall that the balance condition is only enforced between balance-adjacent boxes,
and neighboring boxes are balance-adjacent only if the smaller box contains a point
of the extended box of the larger box. Hence, the single connected component
P A ex(bi) contains P A ex(bj) for all 5 = 1,2,... 1. This line of reasoning is also
correct in the case that b1 is itself crowded, in which case we take 1 = 1.
Recall that x is a point on a P face F in P A b1. We wish to find a face G close
to F and x that does not contain F. We now have two cases, depending on why
b1 was split.
Case 1. The first case is that b1 is crowded because a face interferes with F,
and this face is visible to x. In this case, we let G be the closest interfering face
to F, and y the closest point of G to x. Then we conclude 0 does not contain F,
and distp(x,y) ? 4$3h(b1).
92
b1
Figure 4.17: An example path in the proof of Lemma 23. Here all P vertices
represent edges, and edges represent facets, except for the P vertex J = z. Note
b3 is crowded and b1 is split by the balance condition.
Case 2. The second case consists of one of three subcases. This first subcase
is that b1 is split because of the balance condition. The second subcase is that b1
is crowded because some face in b1 is interfered with, but this face is not F. The
third subcase is that b1 is crowded because some face interferes with F, but this
face is not visible to x. In the second and third subcases 1 = 1. In any subcase,
since b1 is crowded, we may let J denote a P face in b1 such that some P face in
ex(b1) interferes with it. It may be that J is F. Let H be the face interfering with
J such that the geodesic distance between H and J in P A ex(b1) is m?mmum.
Since P A ex(bi) is a single component containing P A b1, there is a shortest path
entirely in P A ex(b1) from x to a point z of J in b1 visible to H. This path may
be considered as a chain of maximal length edges, from x to xi, xi to X2, ... Xm
to z. A cross-sectional view of an example path is given in Figure 4.17.
If the smallest dimensional face containing xi is not F, we choose G to be this
face, and y to be xl. Our choice of C as the smallest dimensional face containing
y ensures that G does not contain P. Since there is always a straightline path in
ex(b1) from x to xi = y, we conclude that distp(x,y) < 3$3h(b1).
93
Otherwise the smallest dimensional face containing xi is F, that is x1 is in the
relative interior of F. Because the path construction uses maximal length edges,
the only way this can happen is if xl is actually z. Hence there is a straightline
path on F from x to z, and J is contained in F. See Figure 4.18. We choose
G to be H, the face interfering with J, and y to be the closest point of H to z.
Since 0 does not contain J, and J in contained in F, it follows that 0 does not
contain F. Since there exist straightline paths from x E b1 to z E b1 to y E bj,
with h(b1) < h(b1), we conclude that distp(x,y) < 6$3h(b1). I
Lemma 24 For any vertex v of P in a box b, f*(v) < ci5h(b).
Proof. Since the protecting and merging steps for a protected vertex box change
the size of any box by at most a factor of two, we can assume b is a box at the end
of a vertex phase. Let b1 be the parent of b.
We apply Lemma 23 with x = v in box b1. For any vertex, if another face does
not contain it, the two are disjoint. Hence the point y on P face 0 is disjoint from
F = v and from the construction distp(v, y) < 4$3h(b). We apply Theorem 10 to
conclude f*(v) < 4?3h(b). I
clo
Lemma 25 (Corollary) For any pointp ina boxbvprotected for vertex v, f*(p) <
cj6h(bv).
Proof. In a protected box Euclidean distance to the apex is the same as geodesic
distance, so distp(p, v) < $3h(bv). From Theorem 11 and Lemma 24, f*(p) <
max(ci? f*(v), ci? distp(p, v)) < max(q?ci?, $3ci4)h(b). I
Lemma 26 For anypointp ofP edge F in a box b, f*(p) < ci7h(b).
94
by
x
v
y
Figure 4.18: An illustration of an instance of Case 2 of the proof of Lemma 23
leading to the case in the proof of Lemma 26 where F' and G' have a common
vertex v.
Proof. As in the vertex case, Lemma 24, since the splitting steps within a pro
tecting step subdivide any box at most once, we may assume that the parent b1
of b was split in the vertex or edge phase. Assume that the box containing p at
the end of the octree algorithm is not vertex protected, since otherwise Lemma 25
applies.
We apply the Lemma 23 with x = p. Unlike the vertex case, Lemma 24, we
may have that G and F intersect nontrivially. To complete the proof we need to
consider the particular triangulation A*. Let F, and C' be the smallest dimensional
faces in A* such that x E P' c F and y E G' c G. If F' and G' are disjoint, we
can apply Theorem 11 as in Lemma 24 and the proof is complete. Otherwise, by
conformality, F' n G, is a P vertex v. Let bv be the box containing v at the end
of the vertex phase. See Figure 4.18. Since v is centered in b,, and by assumption
p is outside bv, distp(p, v) > h(bv)/8. But h(bv) > ci5f*(v) by Lemma 24. Hence
by Theorem 10 Case 1, distp(p, y) > 8c9ci5f*(v). Recalling from Lemma 23 that
distp(p,y) < 6$3h(b1) completes the proof. 1
95
Lemma 27 (Corollary) For any point p in a box bF protected for edge F, f*(p) <
clsh(bF).
Proof. The proof is identical to the proof of Lemma 25. I
Lemma 28 For any point p of P facet F in a box b, f*(p) < ci9h(b).
Proof. As in the previous lemmas, we may assume that the parent b1 of b was
split in the vertex or edge phase, and not a protecting step. Assume that the box
containing p at the end of the edge phase, is not vertex or edge protected, since
otherwise Lemma 25 or Lemma 27 applies.
We apply Lemma 23 with x = p. As in Lemma 26, F and U may not be
disjoint, so we must consider the particular triangulation A*. Let ?? and ?, be
the smallest dimensional faces in A*, such that x E F' c F and y E 0' c 0.
if F' and 0' are disjoint, we can apply Theorem 11. Otherwise, by conformality,
F' n 0' is a P vertex or a subset of a P edge. Let v be the closest point of F' n 0'
to p. Let bv be the box containing v at the end of the edge phase. Since v ?s
well centered in protected boxes, and p is not in a protected box by assumption,
distp(p, v) > h(bv)/8. Depending on whether bv is vertex or edge protected, we
rely on Lemma 25 or Lemma 27 to show h(bv) > min(ci?,cis)f*(v). Hence, by
Theorem 10 Case 1, distp(p, y) > 8c? min(q?, c18)f*(v). Recalling from Lemma 23
that distp(p, y) < 6?h(b1) completes the proof. I
Lemma 29 (Corollary) For any point p of a box bF protected for facet F f*(p) <
c?oh(b?).
Proof. The proof is identical to the proof of Lemma 25. I
96
Lemma 30 For any point p in an unprotected box b (containing no P faces),
< c21h(b).
Proof. As in the previous lemmas, we may assume that the parent b1 of b was
split in the vertex, edge, or facet phase. Either b1 is split by the balance condition,
or b1 is actually crowded, but no P face is in the octant of its child b. By imitating
the path construction of Lemma 23, with the only difference being that x is not
on any P face, we may identify as G the first face we meet on a path from p to
a face in the crowded box bj. That is, we may show that distp(p, y) < 3?h(b1),
and h(by) < h(b1), for some point y on face G in box b? at the end of the octree
generation algorithm. Depending on whether by is a vertex, edge or facet protected
box, we apply one of the previous three corollaries, Lemma 25, Lemma 27 or Lemma
29 to show max(c16,c18,c20)f*(v) < h(by). We apply Theorem 11 to p, so that
< max(ci?f*(y), ci?distp(p, y)) < max(q?/ min(ci?, ci8, c?o), 3?q4)h(b). 1
We may now summarize the results of the previous lemmas.
Theorem 12 Let x E P lie in a box b. Then f*(x) < c22h(b), where c22 is a
constant depending on the aspect ratio bound of A*.
The following result combines the bounds derived for focT and f*.
Theorem 13 For all polytopes P and constants A, there exists a constant c', de-
pendent only on A and c, such that for any z ? P and any triangulation A* of P
with maximum aspect ratio A,
focT(z) >
Proof. This follows directly from Lemma 21 and Theorem 12 I
97
Theorem 14 For all polytopes P and constants A, there exists a constant ??,
dependent only on A and a, such that
I AoCT I < c111 A* I,
where I A* denotes the number of tetrahedra of A* and similarly for I AocT
Proof. Let AT be a triangulation with bounded aspect ratio. We derive some
inequalities that depend on aspect ratio; these inequalities apply to either A* or
AocT By the definition of fT and bounded aspect ratio A, there are constants
c23 and C24 depending only on A such that
c??f?(x)3 < vol(Tj) < c??f?(x)3
for all x in the interior of Tj and for any tetrahedron Tj.
Thus
c2? < fTj 1 dx < C24
fT(x)3
for every tetrahedron. Summing this chain of inequalities over all tetrahedra in
the triangulation, we get
C23?IATI<?WfT(i?)3dX=Zi!TifT(1r)3 dx < c???I AT
This holds for all triangulations, so in particular for AocT and A*. By Theorem
13
c23 I AoCT I < L ?oc??x?3 dx < f? ?????x??3 dx < c??/(c')? I A*
98
4.9 Conclusions
An important open question is the running time. We demonstrated a running
time bound of O(?n log n) for our algorithm, which is inferior to Bern, Epstein
and Gilbert's running time of O(n log n + ?) for two dimensional regions. However,
we have recently found out by private communication with those authors that
this bound does not actually hold for the algorithm that they propose for two
dimensional nonconvex polygons. Therefore, it is not clear what is the best possible
running time for triangulation of nonconvex polyhedral regions in either two or
three dimensions.
It would be interesting to generalize our algorithm to work for higher dimen-
sional regions. Also, the non-degeneracy conditions on the input region could be
relaxed.
Another open question concerns optimizing the tetrahedra for properties other
than aspect ratio. For example, is there an algorithm for triangulating a three-
dimensional nonconvex polyhedral region to maximize inscribed sphere radius of
the tetrahedra?
Finally, we have plans to implement this algorithm. I implemented the two
dimensional analog of this algorithm in C++ while at Xerox PARC.
Chapter 5
Maxmin angle triangulation with
no boundary Steiner points
We show how to triangulate a two dimensional polygon with holes. We assume the
restriction that we may introduce Steiner points on the polygon's interior, but not
on its boundary. Of all triangulations satisfying this restriction, our triangulation
has the maximum minimum angle, up to a constant factor. This is the first known
algorithm for this problem with provably optimal element shape. The algorithm
solves several subproblems of mesh generation.
The largest possible minimum angle of a triangulation of a polygon P is de-
pendent on the geometry of the polygon. Consider an edge E = VU of P, and W
a point of a face of P disjoint from  (interfering with E). Let r be the ratio of
the distance from W to the closer vertex V of E to the length of E, or L(E). Let
w = LUVW. Then we show that any triangulation of P has min angle at most
h(E,W) = O(w/max(1,--Hlogr)).
99
loo
Hence an upper bound on the max min angle of a triangulation of P is g(P), the
minimum of f(E, W) over all E and W of P.
The algorithm we describe in this chapter constructs a triangulation with min-
imum angle at least a constant factor times g(P). Thus our algorithm is optimal,
up to a constant factor.
5.0.1 Two dimensional algorithm overview
First, we triangulate part of P by introducing sequences of triangles called wedges
at vertices of P. The triangles have min angle within a constant factor of the
optimal min angle of a triangulation of P. We call the remaining untriangulated
polygon Q. We may triangulate Q using more traditional techniques, and still
achieve min angle within a constant factor of optimal.
So that we do not need to introduce vertices on the boundary of Q, we shrink
Q to the polygon R, where 1? closely resembles a small copy of Q lying inside Q.
We triangulate R with the two dimensional analog of the algorithm in Chapter 4.
This triangulates R with min angle within a constant factor of optimal, but also
introduces Steiner points on the boundary of R.
Lastly, we match R to Q, that is we triangulate the region of P between Q and
R using the edges of Q as a guide. No new vertices on the boundary of either Q
or R are introduced in this last step.
The proof of the optimality of our algorithm is a step-by-step argument that we
never introduce a triangle with an angle smaller than a constant times the optimal.
101
R3
v1
v			E=s1			U
Figure 5.1: The triangles Vi, spokes Si, and rim edges Ri of a wedge at V.
5.1
Constructing wedges around vertices
In this section we begin the triangulation of P. We form triangles that share an
edge or vertex with ap. We call the portion of P which remains to be triangulated
Q. The boundary of Q is completely composed of faces from the triangles we
introduce in this section. That is, no triangle of a later section will contain an edge
or vertex of ap. We first provide motivation for how we should triangulate.
The following definitions differ slightly from those given in Chapter 4.
Wedge. A wedge is a sequence of triangles such that all triangles contain a
common vertex, and consecutive triangles contain a common edge. The common
vertex V is called the vertex of the wedge, and the wedge is often called the wedge
at V. The edges containing V are called the spokesof the wedge, and are ordered
as in Figure 5.1. The edges opposite V are called the rim edges and are similarly
ordered.
Interfere. For a closed face F of a polygon, we say that the closed polygon
face G interferes with F if F and G are disjoint. We also say that any point of G
interferes with F.
102
Length. We use the notation L(E) to donote the length of the edge E.
Our goal is to produce a triangulation of P with max min angle at least a
constant factor times the optimal min angle. In order to get a lower bound on the
optimal min angle, we consider a class of triangulations. Any triangulation will
have a wedge at every vertex V of P. Let E be a P edge containing V, and W a
P vertex interfering with E. We consider an arbitrary triangulation from the class
of triangulations that have an edge from V to W. We produce an upper bound
on the min angle of the wedge at V with first spoke E and last spoke VW. We
do this by constructing an "optimal" wedge at V. The optimal wedge is closely
related to the "A-optimal" wedge that we define below.
A-optimal wedge. We define the A-optimal wedge to be the wedge that has
optimal min angle for the problem of constructing a wedge at V subject to only
the constraints that its first spoke is E and its last spoke is VW. Note that such
a wedge may intersect the other faces of a? arbitrarily. Since, by assumption, in
the triangulation under consideration there is a wedge at V that also satisfies these
simple spoke constraints, the min angle of the A-optimal wedge is a lower bound
on the min angle of this triangulation. Our optimal wedge is a wedge consisting of
similar triangles that has min angle within a constant factor of the min angle of
the A-optimal wedge. Note that Theorem 18 below shows that there really is an
A-optimal wedge, and not an infinite sequence of wedges with increasing minimum
angle. To summarize, we have proved the following theorem.
Theorem 15 Let P vertex W interfere with P edge E = VU. Then among all
those wedges at V with first spoke E and last spoke VW, there is a wedge with
optimal min angle called the A-optimal wedge. The min angle of any triangtlation
103
ofF that has edge VW is at least that of the A-optimal wedge.
For simplicity we restrict our attention to the case that L(E) > L( VW). To
get a stronger grasp on what the triangles of the A-optimal wedge are like, we have
the following three theorems.
Theorem 16 Let 0 be the angle at V for a triangle T of an A-optimal wedge in
Theorem 15. Then
0 <
Proof. Suppose that the two spokes of T are of equal length (this is the limiting
case). Let ? be the angle at the other vertices of T. If 0 > ?/3, then ? < 0, and
the min angle of T is given by (7r --H 0)/2. Otherwise 0 < 7r/3 and 0 is the min
angle of T. For 0 > 7r/2 the min angle of T is smaller than the min angle for the
triangle with angle at V half that of T. llence, if T had angle at V greater than
ir/2, we could divide it into two and improve the min angle, a contradiction.
If the spokes of T are not the same length, then creating two triangles by
bisecting the angle at V with an appropriate length spoke improves the min angle
even more. 1
This bound allows us to write an approximation of the min angle of the optimal
wedge in terms of a and b.
Theorern 17 Let A be the min angle of a triangle T of an A-optimal wedge at V
in Theorem 15. Let a be the ratio of the length of the shorter spoke to the longer
spoke of T, and 0 the angle of T at V. Then
k1a0 < A < k2a0,
where k1 and k2 are constants.
104
Proof. Let the longer spoke have length 1 and the shorter spoke length a. Let the
rim edge have length c. Let ? be the angle at the vertex of the longer spoke that
is not V. So, either ? = A or 0 = A.
Suppose A = 0 < ?. Then, from the law of sines, c < a. But from the triangle
inequality a + c> 1 so that a > 1/2. By definition a < 1 so laO <0 <2a9.
Otherwise A = < 0. Then, from the law of sines, sin ? = (a/c) sin 0, and
_ a. Since c + a> 1 we have c> 1/2. But also, from Theorem 16, we have that
c < Vm?2+1 < $2. Hence (a/$2) sin 0 <sin(A) = sin ? < 2a sin 0. Noting that ?
is bounded above by 7r/2 shows that the sines are within a constant factor of their
arguments, and completes the proof. I
it is surprising that we may also derive a constant upper bound on the ratio
between successive spokes.
Theorem 18 Let a be the ratio of the lengths of the shorter spoke to the longer
spoke of a triangle T in an A-optimal wedge in Theorem 15. Then
1/4 < a < 1.
Proof. That a < 1 follows from its definition. To show a > 1/4, consider the min
angle approximation aO. Assume the longer spoke has unit length. Divide T into
two similar triangles by introducing a new spoke along the angle bisector of 0. Let
the length of the new spoke be a' = $a. Then the new triangles have min angle
approximation a'0' = $a0/2. This is an improvement over a0 if a > 1/4. Using
the actual min angle instead of the approximation shows a greater improvement.
Hence if a was less than 1/4, then adding more spokes would improve the min
angle of the A-optimal wedge. But since by definition the A-optimal wedge cannot
be improved, this is a contradiction. I
105
Optimal wedge. We now define an optimal wedge to be the wedge that has
the max min angle approximation aTOT from Theorem 17, subject to 0T <
Contrast this with the definition of an A-optimal wedge, which is the wedge that
has the max min angle. Formally, the optimal wedge for a P edge E = vU and P
vertex V, and interfering point W, is the wedge at V that has maximum minimum
aTOT over all triangles T of the wedge, subject to the constraint that the wedge
has first spoke E, last spoke VW, and 0T < 7r/2 for all triangles T of the wedge.
We consider the optimal wedge because it is easier to analyze than the A-optimal
wedge.
Theorem 19 All triangles in an optimal wedge are similar, with spoke lengths
following a geometric progression. The parameters a, the ratio of stccessive spoke
lengths, and 0, the angle of a triangle at V, determine the shape of an optimal
wedge.
Proof. It is sufficient to show this for any two consecutive triangles T and U in
an optimal wedge at V. Suppose that the number of spokes has been determined,
and also the lengths of all spokes except for the spoke shared by both T and U.
We now show that in order to maximize the minimum angle approximation of T
and U, we should chose this shared spoke so that T and U are similar, and these
three spokes form a geometric progression. Let a? be the ratio of the length of the
longer spoke of T to the length of the shorter spoke, and 0T the angle of T at V.
Similarly define au and Ou.
We first prove that in an optimal wedge aTOT = auOu. Suppose aTOT <auOu,
then we could decrease Ou and increase 0T by rotating the shared spoke about V,
and the resulting max min angle approximation would increase, a contradiction.
106
Next we have three cases corresponding to possible ranges of the length of the
shared spoke in the optimal wedge. Let the shared spoke have length ss Let the
spoke unique to T have length 5T, and the spoke unique to U have length su.
Without loss of generality assume that 5T/5U = r < 1. Let the sum of the angles
ofTandUat Vbew.
Case 1, 5T < 5s <su. Then the problem can be stated as
Maximize aTOT
5. t. aTOT = aUOU
aTaU = r
0T+0U =W
It			It
T < 2' 0U < 2
If we ignore the last constraint, we can use the other constraints to solve for 0T in
the objective function. The problem reduces to
Maximize			--H 0T),
which has solution 0T = w/2. By assumption, the number of spokes have been
determined, so that w < It. Hence this solution satisfies the constraint that we
ignored. Solving for the other variables leads to Ou = w/2 and a? = au =
which shows that triangles T and U are similar. Moreover s? = $rsu, and 8T =
so the spokes lengths are in a geometric series. The solution value is $rw/2.
Case 2, 5? < 5T. Then it follows that su > 5s, and the problem can be stated
as
Maximize aTOT
107
5. t. aTOT = au9u
au/a? = r
0T+0U =W
7r 7r
T < 2' 0U
Using similar techniques, this has optimal solution at 0T = w/(1 +r), Ou = wr/(1 +
= 1 and au = r. The solution value is w/(1 + r). The solution is the same
as that of Case 1 for r = 1. For the case that r < 1, the solution value is inferior
to that of Case 1, and so this will not be the optimal solution.
Case 3 S? > Su. Then if follows that s? > 5T? and the problem can be stated
as
Maximize auOu
5. t. aTOT = aUOU
a?Iau = r
0T+0U =W
o?<?? ??<rn.
2'			--H 2
By symmetry with Case 2, the solution here will always be inferior to the Case
1 solution for r < 1, and identical to it for r = 1. I
This allows us to write the optimal wedge for E = VU and W as the solution
to a simple optimization problem. Let r = L( VW)/L(E) < 1, the ratio of the
length of the last spoke to the length of the first spoke. Let a < 1 be the ratio of
the size of a triangle to the one preceding it. Let w be the angle between VW and
E, that is w = /WE = ZWVU. Let the angle at V of any triangle be 0. Then
108
the number of triangles in the wedge is w/O. Our constraint that the first spoke is
E and the last is VW may be written as aW/0 = r. Thus the optimal wedge is the
solution to:
Maximize a0
st awI8 --H r
w/O > 1, integer
In light of Theorem 17 (and its dependence on Theorem 16), the optimal wedge
has min angle within a constant factor of the solution value of the A-optimal wedge
if 0 < r/2. Let us consider the relaxed problem in the case that r < 1:
Maximize aO
5. t. aw18 = r
Since r < 1, we have that a < 1, and the constraint may be rewritten 0 =
w log a/log r. Hence we may rewrite the problem as simply
Maximize a log
From calculus this has solution a = e?1, and solution value
w
--Helogr
The corresponding value of 0 = --Hw/ log r. Noting that the minimum angle is
always positive, and considering separately the cases when 0 is greater than and
less than r/2, this is sufficient to show that a constant factor times 77eWogr is an
upper bound on the min angle of the A-optimal wedge for E and W.
However, for r <e?1, we wish to find a wedge with min angle at least a constant
factor times this upper bound. In general, in the solution to the relaxed problem
109
w/O = --H log r > 0 is non-integer. However, in the case that r <e?1, we take the
relaxed problem solution, and decrease 0 to the nearest value such that w/O is an
integer. We then increase a so that a"?18 --H r is once again satisfied. Since r <e?1
implies w/O > 1, at worst we halve 0, and a is increased, so that there is a solution
to the unrelaxed problem with value at least one half the relaxed solution. That
is, the unrelaxed problem has solution value at least
w
--H2elogr
Moreover, any feasible solution to the unrelaxed problem with w/O greater than
some value has min angle worse than this feasible solution, hence there is an optimal
wedge and not an infinite sequence of wedges approaching an optimal value.
In the case that r > e?1, we merely note that no triangle can have angle at V
larger than w. Hence, no wedge at V can have min angle more than w. In this
case it is easy to construct a wedge at V, consisting of no more than four similar
triangles, with 0 at least a constant times w and no more that 7r/2, and with a
bounded by a constant. Thus we have proved the following theorem.
Theorem 20 For every edge E and interfering vertex W, there exists an optimal
wedge consisting of a finite number of triangles. The min angle of the optimal wedge
is O(w/ max(1, --H log r)). The minimum angle of the optimal wedge is bounded above
and below by constant factors times the minimum angle of the corresponding A-
optimal wedge.
We narrow the definition of wedge slightly to differentiate between wedges that
we construct and wedges of arbitrary triangulations.
Spiral. A wedge that we construct that consists of similar triangles is called a
spiral.
110
Optimal wedges are spirals. In the following discussion, we will be concerned
with spirals only when E is the longest edge of the spiral containing V, and the
interfering point W defining the spiral is closer to the vertex of the spiral than to the
other vertex of E. The first condition will always be satisfied when dist(V? W) <
L(E). The P edge E will always be the first edge of the spiral.
Theorem 21 Given P edge E = VU and interfering point W, if dist(V, W) <
L(E) and dist(V, W) < dist(U, W), then any triangulation must have min angle at
most a constant times that of the optimal wedge with first spoke E and last spoke
vw.
Proof. Note that W is any point interfering with E, and not necessarily a P
vertex. If VW is an edge of the given triangulation, then Theorem 20 applies. So
there is only something to prove in the case when there is no triangulation edge
from V to W. Note that there will still be a wedge at V, but not one with last
spoke VW. Let the sequence of triangles V1, V2,... Vk be the wedge at V such
that E is the first spoke and W is in the infinite sector defined by the spokes of
Vk. See Figure 5.2. The outline of the proof follows.
We transform the wedge at V in a given triangulation to one that includes
the edge VW. We use different transformations depending on whether k > 1 or
k = 1. For k> 1, there is one step that at most halves the minimum angle, while
all other steps do not decrease the minimum angle. For k = 1, the transformation
does not increase the minimum angle. In either case, the transformed wedge must
have minimum angle at most that of the A-optimal wedge, by definition. Hence
any given triangulation must have min angle at most twice that of the A-optimal
wedge. The proof is completed by noting that the A-optimal wedge has min angle
111
Y W
x
v1
v			E			U
Figure 5.2: The wedge at V in a given triangulation where the sector of the last
triangle Vk = AYVX contains the interfering vertex W.
at most a constant factor larger than the min angle of the optimal wedge. We now
describe the transformations.
If k > 1, then let X be the vertex opposite V of the kth (second to last) spoke.
We have two cases which require different transformations.
Case 1, L( VX) < 2dist(V, W). In this case we will ignore Vk and consider
transformations to the wedge Vl,V2,... Vk--Hi By definition, the wedge Vi, V2,... Vk?1
has min angle at most that of the A-optimal wedge at V for last spoke VX, and
hence at most a constant factor times that of the corresponding optimal wedge.
If we place a vertex X' at the midpoint of VX, the optimal wedge at V with
last spoke VX' has min angle approximation at least half that of the optimal
wedge at V with last spoke VX. This can be seen by noting that we have just
halved r in the optimization problem leading to Theorem 20. Now, moving X' to
W increases r and increases w in the constraint aw10 --H so that the optimal wedge
for last spoke VX' has min angle approximation smaller than that of the optimal
wedge for last spoke VW. Hence, the min angle of the given triangulation is at
most a constant times the min angle of the optimal wedge at V with last spoke
112
x
v			E			U
Figure 5.3: The first subcase of where an interfering vertex W may lie relative to
the spiral triangle V'
vw-
Case 2, L( VX) > 2dist(V? W). We consider a series of transformations of the
wedge V1, V2,.. . Vk that changes it into a wedge that contains W. We replace Vk
by the triangle V' = AXVZ, where the angle of Vk and V' at X are the same,
and the angle of V' at Z is ir/2. We wish to show that the min angle of V' is not
smaller than the min angle of Vk. By assumption the altitude of Vk and V' at V
is at most dist(V? W) < L( VX)/2. Hence the angle ? at X is at most it/6, and
the smallest angle of V' occurs at X. Since the angle at X of the two triangles is
the same, the min angle of V' is at least the min angle of Vk.
We have two subcases. The first subcase is when W is not in the infinite sector
LXVZ. See Figure 5.3. We note that the min angle of the wedge Vi,V2,... Vk?1, V'
is at most a constant times that of the optimal wedge at V with last spoke VZ.
Since the angle at X was unchanged when we replaced Vk by V,, and W is in the
sector at V defined by Vk, but not inside Vk itself, it is always the case that line
ZX separates W from V. Since Z is the closest point of line ZX to V, we have
that dist(V W) > dist(V Z), and /WVU > ZZVU. Referring to the constraints
leading to Theorem 20 shows that the min angle approximation of the optimal
113
x
v			E			U
Figure 5.4: The second subcase of where an interfering vertex W may lie relative
to the spiral triangle V'.
wedge at V with last spoke VW is more than the min angle approximation of the
optimal wedge at V with last spoke VZ. Hence we have proved the theorem for
this subcase.
The second subcase is when W is in the infinite sector ZXVZ. Since also
dist(W? V) < dist(X, V)/2, we note that W is in the shaded region of Figure 5.4.
Recall that the smallest angle of V' occurs at X. We transform ?` by moving
moving Z to W. From the law of sines, the new angle at V is at least the old angle
at X. Also the new angle at X is at least the old angle at X. Hence the min
angle of V' is not decreased by this transformation. We have constructed a wedge
Vj,V2,... Vk?1, V' with min angle at least that of the original wedge at V. This
wedge V1,V2,... V??1, V, contains W as the last vertex so its min angle is at most
that of the A-optimal wedge at V to W, and at most a constant factor times that
of the optimal wedge.
If k = 1, then let Y be the third vertex of V1. We transform V1 into V' by
intersecting it with the half plane containing U defined by the line through VW.
The angle at U is the same for both V1 and V'. By assumption W is closer to V
114
Figure 5.5: The continuous spiral for an optimal wedge.
than to U, so the angle at V of V' is at least the angle at U. The angle at the
new vertex Y, of V' is greater than the angle at Y of the old triangle. Hence V'
has min angle at least that of V1. Next, the optimal wedge for E and Y' has min
angle at least a constant factor times that of V', which is at least that of V1. By
examining the constraints leading to Theorem 20 we see that this optimal wedge
has r smaller than r for the optimal wedge for  and W. Hence the optimal wedge
for E and W has min angle at least a constant factor times that of the optimal
wedge for E and Y', which is at least that of V1, and the proof is complete. I
Continuous spiral. An optimal wedge for P edge E and interfering point W
defines a and 0 as in Theorem 19. We define the continnous spiralfor E and W
to be the tp? ?J polar curve p = L(E)a?/0, where the origin is at V, the ? = 0
axis is aligned with E. We allow ? to vary between 0 and the angle that the other
P edge containing V makes with E. See Figure 5.5. Note that the spoke vertices
115
of the optimal wedge lie on the continuous spiral. All other points of the optimal
wedge lie "inside" the continuous spiral. We say that set Y lies inside a continuous
spiral z if for ally = fp?,O < < 2rJ E Y? there is a z = fpz,?z) E Z such that
= ? and PI' < Pz.
Shrinking factor and optimal spiral. We may rewrite a continuous spiral
as p = L(E)(a'I8)?. We define a110 to be the shrinking factor of the continuous
spiral. The shrinking factor will always be less than one for the continuous spirals
we consider. The shrinking factor and L(E) defines a continuous spiral. The
optimal wedge that defines a continuous spiral has max min angle approximation
for all spirals with spoke vertices on a continuous spiral, angle at V no more than
ir/2, and the fixed last spoke. Hence, we may think of a continuous spiral as
defining an optimal wedge, except that we don't know the angle of the last spoke
of the optimal wedge. So we distinguish as the optimal spiral, the corresponding
optimal wedge whose last spoke vertex has ? that of the other P edge containing V.
The optimal spiral for a continuous spiral has min angle approximation at least a
constant times that of any optimal wedge that defines the continuous spiral. This
fact follows from examining the rounding we did when establishing the optimal
spiral from the relaxed optimization problem solution leading the Theorem 20.
In the following discussion, we will take a given continuous spiral, transform it
by increasing its shrinking factor, and consider the optimal spiral it defines. Hence
we wish to explicitly write how the min angle approximation of the optimal spiral
depends on the shrinking factor of the continuous spiral.
If we ignore where the last spoke of an optimal spiral must lie, then we may
write the corresponding optimal spiral in terms of a shrinking factor 5. That is, we
116
may write it as the solution to
Maximize aO
s. t. a110 --H 5
0<-.
2
If we ignore the last constraint, this has solution a = e?1 and 0 = --H1/logs,
and solution value --H 1/e log 5. The constraint that 0 <7r/2 is only important when
5 is large (close to 1), in which case we can find a spiral with min angle within a
constant factor of the angle between the two P edges meeting at V. Imposing the
constraint of where the last spoke must lie, we may find a feasible solution with
larger a, and 0 no less than a constant fraction of the above solution value. Hence
we have proved the following theorem.
Theorem 22 Let w be the angle between two P edges E and F meeting at a
vertex. Given a continuous spiral for edge E with shrinking factor 5, there is
a corresponding optimal spiral with min angle approximation bounded above and
below by constants times min(w, --H1/elogs).
We now begin to define curves that will contain our first set of triangles in the
triangulation of P. We start with the continuous spiral defined in the following
theorem, and transform it in the next definition.
Theorem 23 Given P edge E = VU , then there is a continuous spiral F called
the empty spiral, such that
1. The only points of a? inside F are closer to U than to V, or are in the P
edges containing V.
117
2. The min angle of the optimal spiral defined by F is at least a constant factor
times the optimal min angle of P.
Proof. Each optimal wedge at V for edge E defined by an interfering point in
turn defines a continuous spiral. Of these continuous spirals, let F be the one with
smallest shrinking rate. Note F will be inside all the other continuous spirals. Each
point interfering with E is a vertex of the optimal spiral it defines, and hence lies
on the corresponding continuous spiral. Hence we may shrink F by an infinitesimal
amount so that 1 is satisfied. Item 2 follows from Theorem 21 and Theorem 22.
If there is no interfering point close enough to V to define an optimal wedge,
then we take S to be the spiral with shrinking factor 1, that is p = L(E). The
corresponding optimal spiral has min angle at least a constant factor times the
angle between the P edges meeting at V. 1
Smallest spiral. For each empty spiral F we define a continuous spiral called
the smallest spiral. Suppose the empty spiral is defined by p = L(E)s?. Then the
smallest spiral is defined by p --H AYLE(54)? The smallest spiral shrinks four times
as fast as the empty spiral. Also, the distance from V to the point at p = 0 in the
smallest spiral is only one quarter the length of E.
Theorem 24 For a smallest spiral 5, the optimal spiral defined by it is called the
smallest optimal spiral and denoted T. Let A be the min angle approximation of
1. A is at least a constant factor times the min angle of the optimal spiral for
the empty spiral at V for E.
118
2. For any rim edge J ofT, let d be the distance from J to a smallest optimal
spiral for another vertex. Then d/L(J) > kA for some constant k.
Proof. The first item follows directly from the definition of the smallest spiral and
Theorem 22.
For the second item, we develop a series of lemmas. We first establish the
length of a rim edge.
Lemma 31 For a point Y = ?py, ?y) on a rim edge J,
for some constant k1.
L(J) < kipy,
Proof. Suppose Y is a spoke vertex. The length of the larger rim edge J containing
Y is L(J) = py sin 0/ sin ?, where 0 and ? are defined as for a triangle of an optimal
wedge. Since a> e?1, ? is never more than a constant factor smaller than 0, and
L(J) < kipy, for some constant k1.
If Y is not actually the vertex of J, then let Z be the vertex of J making larger
angle with E. Since 0 <7r/2, py is at least p? cos ?4T. The bound on L(J) in terms
of pz from the last paragraph completes the proof. 1
We next establish the distance from a point of S to any other smallest spiral.
The proof is not very direct. Let F be the empty spiral defining S. For a point Z
of S with ?z in a certain range we are able to show that it is closer to E than to
F by a constant factor, and its distance to F is reasonably large. If this were true
for all Z of all smallest spirals regardless of ?z, then since all other faces of P are
outside of F, this would show that the two smallest spirals for disjoint P edges are
119
well separated. However, it is not true for points with ?z in the remaining range,
and we must argue directly that any other smallest spiral must be far away from
these points. Let s be the shrinking factor of F.
Lemma 32 For a point Z in S with ?z >
dist(Z, F) > 3pz = 3dist(Z, E).
Proof. For any point W of F, the triangle inequality gives us dist(Z, W) > pw --H
pz. For ?z > ir/2, we have pz < L(A4E54(q2) = L(4E)32r But pw > L(E)s2r. I
This generalizes to
Lemma 33 For a point Z of 5, dist(Z, W) > 3pz for all points W of F with
We are also able to show that points Z with small ?? are closer to E than to
Lemma 34 For a point Z with ?z <
for some constant k, and
dist(Z,F)?			pz
k?ax(1,?Thg(s))?
>2
dist(Z, F)			--Hdist(Z, E).
-$3
Proof. We first find the point W at angle ?w such that pw = 2pz. This is
equivalent to 5?W --H s4?z/2, which has solution ?w = 4?z + --Hlo1gs
It might be the case that no such W exists. That is, perhaps the solution
for pw above is larger than 2r. Then for all W? pw > 2pz. Then the triangle
120
inequality yields dist(W? Z) > pz. But also dist(Z, E) = pz sin ?z <?Pz, and we
have proved the lemma for this Z.
So assume there is such a W. The reasoning of the last paragraph proves the
lemma if we restrict F to the points W' of F such that ?? < ?w.
It remains to consider the distance between Z and any point W' with ?wi >
?w. This distance is at least the distance between Z and the ray at V at angle
?w. Let 6 = --H ?z. If 6 > 7r/2, then this distance is pz and the lemma follows
for all such points W' of F.
Otherwise, 6 < ir/2. Note that dist(Z, W') > pz sin 6, or
dist(W', Z) > pz sin(3?z +			1
--Hlogs
Here sin 6> 26 and we may reduce this to
dist(W',Z) >2 pz
r--Hlogs'
which proves the first assertion of the lemma.
For the second assertion note that dist(Z, E) = pz sin ?z. If ?z = 0, the second
assertion is trivial. Otherwise ?z > 0 and
dist(W',Z) > sin3?z
dist(Z, E) --H sin ?z
Recall that we are in the case 6 < 7r/2, so in particular ? < r/6 and the ratio on
the right hand side of the above is at least 2. I
In the previous two lemmas, we considered the distance from a point Z to F.
However, there may be points interfering with E inside F, but closer to the other
vertex of E than to Z. The following lemmas proves bounds similar to those of
the previous lemmas for this case.
121
Lemma 35 If point W interferes with E and is inside S, then for all Z of S
and
dist(Z, W) > dist(Z, V),
dist(Z, W) > $3dist(Z, E).
Proof. By the definition of smallest spiral, pz = L(4E)s4pz. Hence pz < L(E)/4.
The only possible W interfering with E but inside S are closer to the other vertex
of E than to Z. So W lies on the side of the perpendicular bisector of E opposite V.
Hence dist(Z, W) > L(E)/2 for ?z > 7r/2, and dist(Z, W) > L(E)(2 --H cos
otherwise. Thus dist(Z, W) > L(E)/4 > pz and the first inequality is proved.
To prove the second inequality we consider two ranges of ?z. For ?z > ir/2
we have dist(Z, W) > L(E)/2 > 2dist(Z, E). For ?z < r/2 we have dist(Z, E) =
pz sin ?z < L(4E) sin ?z. For ?z = 0 the second inequality is trivial, otherwise
dist(Z,W)			2--Hcos?w
dist(Z,E) --H			sin?z
This has minimum value ? at ?z = r/3. I
Lemma 36 For point Z of S, and point Zi of Si, if one of the following holds
1.			?z < r/3, or			> r/2, or Z1 is inside 5.
2.			?z1 < 7r/3, or			> r/2, or Z is inside Si,
then
dist(Z, Z1) > 0.07 max(dist(Z, ), dist(Zi, E1)).
Proof. Assume 2 holds. Let G be the closest point of E to Z, and similarly define
G1.
122
The easy case is if 1 also holds. In this case we may relabel so that dist(Z, G) >
dist(Zi, Gi). Note that G1 interferes with E, so from Lemma 32, Lemma 34, or
Lemma 35, we have that dist(Z, Ci) > ?dist(Z, G). But then the triangle in-
equality gives dist(Z,Z1) > dist(Z,01) --H dist(Z1,G1) > --H 1)dist(Z,G).
It remains to consider when 1 does not hold, that is when 7r/3 < ?z < ir/2 and
Si is not inside 5. Suppose the lemma were false and dist(Z, Zi) <0.O7dist(Z, G).
From Lemma 32, Lemma 34 or Lemma 35, we have dist(Z1, G1) < ?dist(Zi, G).
Hence dist(Z,Gi) < dist(Zi,Gi) + dist(Z,Zi) < ?dist(Zi,G) +O.O7dist(Z,C).
But dist(Z1,G) < dist(Zi,Z)+dist(Z,G) <--H (1+0.07)dist(Z,G). So dist(Z,Ci) <
+ 0.07) + O.O7)dist(Z, G) <dist(Z, G).
Hence we have that Gi lies inside the circle at Z with radius dist(Z, E). Since
additionally ?z < 7r/2, we have that the angle that Gi makes with E at V is
no more than ?z + r/2 < 2.S?z. Note G1 interferes with E. Hence either G1 is
not inside 5, and from Lemma 33 we have that dist(V? G1) > 4dist(Z, V), or G1
is inside 5 and from the proof of Lemma 35 we have dist(V G1) > 2dist(Z, V).
Hencedist(Z,G1) >--H dist(V,G1)--Hdist(Z,V) > dist(Z,V) > dist(Z,E). But in the
last paragraph we just showed that dist(Z, Gi) <dist(Z, E), a contradiction, and
hence dist(Z,Z1) > O.O7dist(Z,G).
Since we want a bound in terms of max(dist(Z, E), dist(Zi, E1)), it remains to
consider the case that dist(Z1, G1) > dist(Z, G). Using the same arguments as in
the easy case described above when both 1 and 2 were true, we can show that
dist(Z,Z1) > --H 1)dist(Zi,Gi). I
Only one lemma remains to fill in the case when both suppositions 1 and 2 of
the previous lemma are false.
123
Lemma 37 For point Z of S, and point Z1 of Si, if r/3 < ?z < ir/2, Z1 is not
inside 5, r/3 < ?z1 < r/2, and Z is not inside Si, then
dist(Z,Zi) > 2max(pz,pz1).
Proof. Suppose pz > pz1. We consider the angle ?v1 that V1 makes with E at
V. If ?vj < 4?z, then Lemma 33 applies and dist(V V1) > 3pz. Hence from the
triangle inequality dist(Z, Z1) > dist(V, V1) --H dist(V Z) --H dist(V1, Z1) > 2pz and
the lemma is proved restricted to such V1.
It remains to consider V1 such that ?v1 > 4?z. From the previous para-
graph, it is obvious that we may also restrict our attention to Z1 such that
dist(V Zi) < L(E)/2. This avoids certain complications when considering geodesic
distance below. Define S = ?v1 --H ?z. Note S > r. The main consequence of
this is that the shortest geodesic path from Z to V1 passes through V, and so
dist(Z, V1) = pz + dist(V V1). This and the triangle inequality form the basis of
the proofs in the cases below.
Define ?v to be the angle V makes with E1 at V1. We have two main cases
depending on the orientation of the interior side of E1 relative to that for E.
The first case is when the orientation or handedness of the interior side of E1
with respect to V1 is the same as the handedness of E with respect to V. See Figure
5.6. We have two subcases.
The first subcase is if ?? > ?z1. Then the shortest geodesic path from Zi to Z
passes through V and so has distance at least pv.
The second subcase is if ?v < ?z1. But then since V interferes with V1, we
have pz1 <pv. Hence dist(Z,Zi) > dist(Z,Vj) --H dist(Z1,V1) > pz.
The second case is when the handedness of the interior side of E1 with respect
124
z
v
E
z
v1			E1
Figure 5.6: The interior side of  relative to V is the same as the interior side of
E1 with respect to V1.
to V1 is opposite the handedness of E with respect to V. See Figure 5.7. If ?v > ?,
then the fact that ?z1 < ?/2 insures that ?v --H ?z1 > 7r/2. Hence the closest point
to Z on the ray at V1 at angle (kz1 is V1, so dist(Z,Zi) > pz.
It remains to consider the subcase that ?v <7r. But then ?v < 4?z1 <
so that pz1 < 4pv. Hence by the triangle inequality, dist(Z, Z1) > dist(Z, V1) --H
dist(V1,Z1)>?pz+pv--HpvI4>?pz. I
Thus we have succeeded in showing that for a point Z of a smallest spiral S, its
distance to a point of any other smallest spiral is at least a constant factor times
pzA, where A is the min angle of the corresponding smallest optimal spiral. In
particular, this implies that the curves of two smallest spirals S and 5' for different
P vertices do not intersect. Since the point at angle 0 of 5 is not inside 5', no
point of 5' is inside 5.
Now consider a point Y of a rim edge J of the smallest optimal spiral T. Let
125
z
v
E
ZI
Figure 5.7: The interior side of E relative to V is opposite the interior side of E1
with respect to V1.
Ji and J2 be the vertices of J with ?j1 < ?j?. Note py < pj1. Clearly, the closest
point Z of 5 to Y has pj2 < pz < pj1. Recalling that pj1 and pj2 are within
a constant factor of one another, we see that the distance from Y to a smallest
spiral for another vertex is at least a constant times pyA. Its distance to the
corresponding smallest optimal spiral is at least this distance. Recalling Lemma
31 completes the proof of the second item of the theorem.
This completes the proof of the Theorem 24. 1
5.1.1 Constructing a complete wedge around a vertex
We now describe how we construct a wedge around a P vertex V that includes
both of the P edges containing V.
126
v
E
Figure 5.8: Two smallest spirals for the same vertex V meet at X.
Welding two spirals at a vertex
There a smallest spiral for vertex V and each of the two edges containing it. Our
first task is to "weld" these two smallest spirals.
Let S and T be the two smallest spirals for E and F at V, where L(E) >
Let X be the point common to the boundaries of both S and T as in Figure 5.8.
It may be that both S and T have shrinking factor 1, in which case they are
coincident. In this case we simply ignore T and consider S to be the welded spiral.
Hence we may assume that X is unique. If X lies on E we discard S, and if X
lies on F we discard T, and consider the remaining spiral to be the welded spiral.
Otherwise, we now form the smallest optimal spiral 5' for 5 subject to the last
spoke being VX. That is, in Theorem 22 we take w = ?? instead of w =
Here ?F denotes the angle between F and E at V. Similarly define T'. See Figure
5.9. From Theorem 22 we see that 51 has min angle at least a constant factor
times what is optimal for P, unless ?x is less than half the angle 0 of the smallest
127
optimal spiral when we take w = ?F We now show how to fix 5' in this case.
Should this case arise for T', it is fixed analogously.
First note from Theorem 22 that this case can only happen when ?x <?F/2.
This implies that 5' consists of just one triangle. This also implies that this case
cannot arise simultaneously for S? and T'. llence the min angle of T' is within a
constant factor of what is optimal for P. If the min angle of 5' is within a constant
factor of that for T', there is nothing that needs to be done. Otherwise, let Y
be the second to last spoke vertex of ??, and Z the second to last spoke vertex
of 5'. We move X to the point X' along the line ZX such that VX' bisects
ZYVZ. This transformation keeps T' and 5' inside T and 5. The angle of 5' at
Z is unchanged.
The transformation at most halves the angle of T' at V. Because of the constant
lower bound on a for both 51 and T', the transformation decreases a for T' by at
most a constant factor. Also because of the constant upper and lower bounds on a
for both 5' and T', the min angle of 5' is at least a constant factor times the min
angle of T'. llence the min angles of 5? and T' are both within a constant factor
of the optimal achievable for P, and we have fixed the spirals for this case.
We always make an additional transformation to guarantee a large angle be-
tween the rim edges meeting at X, in order to make Theorem 26 true. We replace
the triangles of 5' and T' that contain X by the three triangles of Figure 5.10.
That is, we replace VX by two spokes, VX' and VX", with a rim edge XIX"
between them. See Figure 5.10. The length of the new spokes is the same as that
of the old spoke VX. For each of the new triangles, the angle at V is one third
of the angle at V between the second to last spokes of 5' and T'. This makes at
128
rF
v			E
Figure 5.9: Welding the two smallest optimal spirals for V.
most a constant factor decrease ;n the min angle of S, and T'. This also insures
that Theorem 28 below is true.
Welding two spirals for an edge
The wedge at V is complete except for the triangles containing P edges E and F.
We consider E, the procedure for F is analogous. We rotate the first spoke N of S,
(the one contained in E), uniformly distributing the angle change over all S', until
the angle N makes with E is the angle of 5'. This at most halves the min angle of
5,. The edge E also contains another P vertex U. We complete the same steps in
the construction of the complete wedge for U as we have for V before continuing.
Call the first spoke in that construction M.
Recall that L(N) = L(M) = L()/4. Without loss of generality, we may
assume that ZNE < ZME. Let X be the other vertex of M, and Y the other
vertex of N as in Figure 5.11. We introduce the edges 7U and YX into our
construction. See Figure 5.11.
129
rF
v			E
Figure 5.10: Replacing VX by the two spokes VX' and VX". This is the final
transformation for welding the two smallest optimal spirals at V.
Figure 5.11: Welding the two wedges for one P edge and its two vertices V and U.
It remains to analyze the min angle of the two new triangles we introduced.
For AYUV, we note that the ratio of the length of E to the length of N is 4, and
its angle at V is the same as that of (the transformed) S'. Hence, its min angle is
within at least a constant factor times optimal.
For AYXU, from the law of sines,
L(YU) _			L(N)
sin ZYVU --H sin ZYUV
Since L(YU) > -3L(E), we have that sin LYUV <?31 sin ZYVU. Since by assump-
tion LXUV> ZYVU, we have /XUY = ZXUV - ZYUV> ZYVU - -31ZYVU =
130
-23LYVU. So the min angle of AYXU is at least a constant times that of AYVU,
and hence also at least a constant times that of the optimal triangulation.
This completes the construction of a complete wedge about any vertex of P.
We form the complete wedge about all vertices of P. This triangulates a portion
of P; the region left untriangulated we call Q, and the triangulated region we call
Qc The boundary of Q is composed of edges of wedges opposite the vertex of the
wedge. Combining the min angle theorems of this section yields the following.
Theorem 25 The min angle of QC is at least a constant factor times the min
angle of any other triangulation of P.
For the following sections, we also have a bound on the smallest interior angle
ofQ.
Theorem 26 The interior angle at any vertex of Q is at least O.39r.
Proof. We must analyze three types of vertices, depending on what types of edges
contain it. The interior angle is at least r for a vertex contained in two rim edges of
a spiral. This follows from the fact that spoke lengths follow a geometric series and
from the upper bound on the angle between successive spokes given by Theorem
16.
The interior angle is at least 7r/2 for a vertex contained in one rim edge of a spiral
and one edge XIX" as in Section 5.1.1. This is the limiting value as L( X'X")
approaches zero and the angle at V between successive spokes approaches zero.
See Figure 5.12 for an example of this construction at a small interior angle of P.
The interior angle is at least 0.39r for a vertex contained in an edge between
two spirals for the same P edge but different P vertices, as in Section 5.1.1. The
131
v			E
Figure 5.12: The interior angle of Q where we welded two smallest spirals for one
P vertex and two different P edges is at least ?r/2, regardless of the interior angle
at V.
v
Figure 5.13: The smallest interior angle a of Q may occur where we welded the
two smallest spirals for one P edge.
limiting case is when the angle for the spiral at one of the P vertices approaches
zero, as illustrated in Figure 5.13. The smallest interior angle a is thus at least
7r/2 --H arcsin(1/3) > O.39r. I
Theorem 27 The interior angle at a vertex of Q is at most 3r/2.
Proof. Consider the interior angle PQ between successive rim edges for the same
smallest spiral. These edges lie in similar triangles of a (transformed) smallest
spiral. So the interior angle of Q at that vertex is equal to 2r minus the sum of
the two angles of the similar triangles that do not occur at V. Hence ?? is equal
132
to ? plus the angle at V of a (transformed) smallest optimal spiral. But the angle
at V of any smallest optimal spiral is at most 7r/2. Hence PQ < 37r/2.
It is clear that the interior angle at rim vertices where we welded two smallest
spirals for the same P vertex is smaller than the case just analyzed. Consider
the interior angle WQ at the first rim vertex, where we welded the two smallest
spirals for a P edge. The limiting case is for the smallest optimal spiral S' left
with the larger angle at V. Consider Figure 5.13, which shows the limiting case
as ZXUV approaches 0. For S' we have a < e?1 but L( VY) = L(E)/4 and
L( VX) > L(E)/2. Hence ZVYX is at least the angle of a smallest optimal spiral
vertex at a triangle of S'. From the analysis of the last paragraph, we see that WQ
is at most 2? minus the angle of S' at V, and hence WQ < 37r/2. I
We now turn to establishing the last theorem of this section which will find
application in the next section.
Delta. For a closed face E of a polygon, we define ?(E) to be the minimum
over all closed faces F C E of the geodesic distance from F to the closest face W
interfering with F.
Theorem 28 For any edge E of ?Q, ?(E)/L(E) is at least a constant times
A(Qc), where A(QC) denotes the minimum angle among all triangles we have in-
troduced between P and Q
Proof. From Theorem 24, the ratio of the length of an edge of ?Q to its distance
to a spiral for another P edge is at most A(Qc). The ratio of L(E) to its distance
to a vertex or edge of the same spiral is at most 1/a, where a is the ratio between
consecutive spokes of the spiral. This is bounded by the constant e. A constant
133
bound also holds for the edge introduced between the first spokes of two spirals for
the same P edge in Section 5.1.1.
We now consider the theorem for an edge E that is the last rim edge of a spiral
at V whose closest interfering point is a vertex of a different smallest spiral at
V. That is, consider when we matched the two smallest spirals for the two edges
containing a P vertex in Section 5.1.1. Let E and F be the rim edges in the last
triangles of the spirals. Let E = XY and F = XZ. Let WE be the angle at V
opposite E, and WF the angle at V opposite F.
Without loss of generality, we seek to bound dist(Z, E)/L(E). Note dist(Z, E) >
L( VX) sin yWF. Since by construction WE = WF, dist(Z, E) > L( VX) sin(??).
Since the first spoke of a spiral is the longest, and WE < r/2, we have that
L( VY) > L(E). Hence dist(Z, E)/L(E) is at least the min angle of the triangle
containing E. I
5.2 Shrinking the polygon Q to R
We have already triangulated Qc, it remains to triangulate Q. in order to trian-
gulate Q, we first shrink it. That is, we form another polygon R, contained in Q,
that is approximately what one would get by moving all the edges of Q towards
its interior.
For every edge E of Q we introduce a shrunken edge F parallel to E in the
interior of Q, at distance ?(E)/4 from E. The placement of the vertices of F
depends on the interior angles at the vertices of E. For every vertex V of Q, we
introduce shrunken vertices according to the following two cases.
Case 1. IfZV < 3r/4, we introduce a single shrunken vertex. This vertex is
134
Figure 5.14: Shrinking Q to R near a vertex V with interior angle < 37r/4 (left)
and interior angle > 3r/4 (right).
placed at the intersection of the shrunken edges for the Q edges containing V.
Case 2. If ZV > 37r/4, then we introduce two shrunken vertices and a shrunken
edge between them. Consider a ray from V at 7r/4 radians clockwise from the bisec-
tor of the interior angle at V. We introduce a shrunken vertex at the intersection
of this ray and the shrunken edge for the clockwise Q edge at V. Similarly, we
introduce a shrunken vertex at the intersection of the ray r/4 radians counter-
clockwise from the interior angle bisector at V and the shrunken edge for the
counter-clockwise Q edge at V. We introduce a shrunken edge between the two
shrunken vertices we just introduced. See Figure 5.14.
Delta for R and Q. We extend the definition of the delta function to R and
Q. That is, for a face F of R we define ?(F) in terms of the geodesic distance in
R between F and its subfaces and interfering faces of R. The following theorems
relating S and the lengths of edges are important for later sections.
Theorem 29 Let F be an edge shrunken for Q edge E. Let V be a vertex of E,
and W the corresponding vertex of F. Then the angle b at V between VW and E
is at least a constant times ?(E)/L(E).
Proof. If we are in Case 2 of the above construction, that is the angle between E
and the other Q edge containing V is at least 37r/4, then the angle b is at least
135
12+1			5(E)14
v			v
%(E)14
E
Figure 5.15: Bounding 4), the angle between VW and E.
r/8.
Otherwise, we are in Case 1 of the above construction. Let E2 be the other
Q edge containing V. Let 0 = LEE2. Define the distances 1 and 12 as in Figure
1 We have that I --H			5 2			________
--H 4tanO and 12 = 4sin Hence 1+12 > k1 6(E)+p6(E2) for some
constant k1. Since tan 4) = 4(S1?1)2), we have the bound that
4) > k2?(E) + ?(E2)
If ?(E2) < ?(E), then 4) > k20/2. Let us ignore the fact that 0 is bounded by a
constant for polygon Q, and prove the theorem for any polygon. In any polygon,
we may get a bound for ?(E) in terms of 0 as follows. Suppose L(E) > L(E2). The
vertex of E2 opposite V is disjoint from E, and its distance to E is L(E2) sin 0.
Similarly, if L(E) < L(E2), then note that the vertex of E opposite V is disjoint
from E, and its distance to E2 is L(E) sin 0. Since this vertex is a subface of
E, this distance bounds ?(E). Hence S(E) < min(L(E), L(E2)) sin 0. That is,
0 > ?(E)/L(E), and 4) > k2442SLEE.
Otherwise ?(E) < ?(E2), so that 4) > k24\OsSEE2. But from the last paragraph, we
may also claim that ?(E2) < min(L(E), L(E2)) sin 0. Substituting this expression
for S(E2) in the equation bounding 4) allows us to conclude the theorem I
136
Theorem 30 For any edge F of R arising from shrinking edge E of Q, ?(F) >
?(E)/2. Also, for any vertex or edge F of R arising from shrinking vertex V of Q,
6(F) > 6(V)/2.
Proof. We first show that the shrunken edges are close to the corresponding Q
edges in terms of 6. This leads to the fact that shrunken edges for two disjoint Q
edges are far from one another in terms of 6. We next note that the shrunken edges
for Q vertices have a large length in terms of 6, and are close to the corresponding
Q vertex. Finally we note that shrunken vertices for Q vertices with small interior
angle are far from other faces. This allows us to conclude that that shrunken edges
have bounded 6 even when considering all faces of I?.
Recall that a shrunken edge F is introduced parallel to a Q edge E at distance
6(E)/4. The angle between a shrunken edge F and the ray from one of its endpoints
through the corresponding Q vertex V is at least 7r/4. Hence any point on a
shrunken edge F for Q edge E is at distance at most ?6(E) from the closest
point of E. For an edge F shrunken for vertex V, we note that since 6(V) is at
least 6 for any edge containing V, the distance from a point of F to V is at most
4?26(V). The triangle inequality allows us to conclude that two shrunken edges E
and F for disjoint Q edges are at least distance (1 --H #2)6(E) from each other.
For any shrunken vertex W arising from a Q vertex V with interior angle greater
than 37r/4, the distance from W to V is at most 4F26(V). Hence the distance from
any shrunken edge for a Q edge E to a shrunken edge for a disjoint Q vertex is at
least (1 --H %)6(E). Also, since the shrunken edge for a vertex V with large interior
angle subtends an angle of it/2 at V, its length is at least the distance from V to
the further shrunken edge for a Q edge containing V. Consider a Q vertex V with
137
edges E1 and E2 meeting at an interior angle larger than 37r/4. Let F1 and F2
denote the shrunken edges arising from E1 and E2. Then the distance from F1 to
F2 is at least max(?(Ei),?(E2)).
Using the above arguments, we can also show that for a shrunken vertex W
arising from a Q vertex V with edges E1 and E2 meeting at an interior angle
greater than 3r/4, S(W) > min(?(E1),6(E2)).
In order to bound 6 for a face, it remains to consider the distance from that
face to a Q vertex V with edges E1 and E2 meeting at interior angles less than
37r/4. Let F1 and F2 denote the shrunken edges for E1 and E2 and W the shrunken
vertex for V. Note that the shrunken edges meet at the same angle as the Q edges.
Since this angle is less than 7r, to determine 6 for a shrunken edge that is disjoint
from F1 and F2 it is not necessary to consider its distance to W. However, it may
be the case that the closest disjoint face to F1 is the other vertex of F2. Hence we
need to consider the length of F1 and F2. See Figure 5.16. Assume that the closest
face disjoint to E1 is X, the vertex of E2 that is not V, and similarly assume that
the closest face disjoint from E2 is Y, the vertex of E1 that is not V. Assume that
E1 and E2 are of equal length, and the angle at V is acute. Introduce F2 first. We
see that the line through F2 cuts off from E1 an edge of length 3/4 times that of
E1. If we now introduce F1, we see that F1 is of length 3/4 that of E1. If a similar
situation occurs at the vertex of E1 that is not V, we see that F1 can be as short
as 1/2 the length of E1. Similarly for F2. If there is some other face of Q that is
the closest face disjoint from E1 rather than X, then this only makes F1 longer.
If the angle at V is not acute, then 1 - --H? times the length of E1 is a weak lower
bound on the length of F1. We omit many tedious cases and symmetry arguments,
138
x
F2
E2
F1
Y
v			E1
Figure 5.16: Bounding the lengths of shrunken edges F1 and F2 meeting at an
angle < 37r/4.
and conclude that a lower bound on ? for a shrunken vertex W for Q vertex V is
--H #2) I
5.3 Triangulation of R
We now triangulate the shrunken polygon R, adding Steiner points to its boundary.
We use the two dimensional analog of Chapter 4. Alternatives are discussed in
Section 5.5.1.
A feature of this triangulation algorithm is that the size of a quadtree box
containing an input edge is at least a constant ci times ? for that edge. Also,
the edge of a triangle contained in an input edge is at least a constant C2 times
the size of the quadtree box containing the triangle. Hence a shrunken edge F is
subdivided by the triangulation of R into edges of length at least c1c2?(E).
Balanced. We also note that the quadtree is balanced: Any quadtree box
139
containing an R edge is at least 1/2 as large as any adjacent box. Any quadtree
box containing an 1? vertex is at least a constant C3 times as large as any adjacent
box. Hence, any two adjacent triangulation edges on the boundary of R differ in
size by at most a factor of 1/(c3c2).
The constants cl, C2 and C3 all depend on the smallest interior angle of R.
However, the smallest interior angle of R is equal to the smallest interior angle of
Q, and Theorem 26 bounds this by the constant O.39ir. Hence ci, C2 and C3 are all
bounded below by true constants, namely k1, k2 and k3.
In the upcoming section we will refer to the size of the quadtree box containing
a given triangulation edge C as
5.4 Matching R to Q
We now show how to triangulate the portion of Q outside of R. We do this by
matching an edge or vertex of Q with its shrunken edges to form a trapezoid or
triangle. Each such figure is triangulated with min angle at least a constant times
?(F)/L(F) for some R edge F. We add Steiner points to a figure's interior and to
its sides, but not to any edge of R or Q. Adjacent figures share sides, so that care
must be taken when adding Steiner points to the sides.
Match faces. We introduce an edge from each Q vertex to the corresponding
vertices of R. We call these edges match edges. This partitions Q \ R into match
trapezoids and match triangles. There is one trapezoid for each edge of Q, and
one triangle for each Q vertex with interior angle greater than 3?/4. For each
trapezoid, one of the parallel edges is a Q edge and the opposite R edge is divided
into several smaller edges by the quadtree triangulation of R. Similarly for the
140
F
M2
M1
E
Figure 5.17: Triangulating a special case trapezoid with layers of boxes.
triangles; one vertex is a Q vertex, and the opposite R edge is subdivided into
several smaller edges.
5.4.1 Triangulating a match trapezoid
Consider a particular trapezoid T, with Q edge E and I? edge F. Let d be the
distance from F to E, and recall that d is bounded above and below by constant
factors times 6(E). Consider a set of lines parallel to E and F, the first one halfway
between E and F, the second halfway between the previous line and F, etc. That is,
consider the set of lines M1, .2.... M? parallel to F, at distance d/2, ....... , d/2?
from F, where n = [log(L(E)/d)J.
We triangulate T by first introducing squares with two sides parallel to F,
where those sides lie on these parallel lines. We then triangulate each square
independently, without adding vertices on the boundary of the square. The min
angle of this triangulation will be at least a constant factor time A, the optimal
min angle of any triangulation of Q.
141
A special case
Before describing the general algorithm and proving this min angle bound, we
consider an easier special case. First recall from Section 5.3 that in general q(X)
is at least a constant k1 times d for any point X E F. So suppose that q(X) = k1d
for all X E F. Then we divide AIn into edges of length 2k1d, and triangulate as
in Figure 5.17. The distance between Mn and F is d/2? > dS(E)/L(E). The min
angle of a triangle is within a constant factor of the ratio of the radii of the smallest
containing sphere to the inscribed sphere. The containing sphere for any triangle
has radius within a constant factor of k1d. The inscribed sphere has radius within a
constant factor of d6(E)/L(E). Hence the min angle of the triangles we introduce
is within a constant factor of ?(E)/(k1L(E)). Note that this holds for the triangles
intersecting a match edge, since the angle between a match edge and F is at least
A, the optimal min angle of any triangulation of R. We have drawn Figure 5.17
as if this angle was r/2. We may repeat this process for the lines M2, .3.... Mn.
Each time we triangulate with min angle at least a constant times 6(E)/(k1L(E)).
We will subdivide Mi into a constant number of edges (linearly dependent on k1),
each of length at least a constant times the length of E. Hence we may triangulate
between M1 and E by introducing an edge from each vertex of Mi to the closer
vertex of F (and one more edge for the remaining trapezoid). The resulting min
angle is within a constant factor of ?(E)/(k1L(E)). We have drawn Figure 5.17 in
the case when M1 is divided into four segments.
Drawing squares
We now show how to modify the special case to handle the general case when F
is not subdivided uniformly by the quadtree. Recall that the quadtree is balanced
142
as described in Section 5.3.
We define the level of an edge to be defined in terms of the size of the quadtree
box containing it, but first we need to make a special definition of size for boxes
containing a vertex of F. The quadtree is not necessarily balanced by a factor of
two at boxes containing a vertex of F. We treat a quadtree box bv containing a
vertex V of F as if it had the same size as the quadtree box bF, where bF is the
unique box that shares a point of F with bv. Since the ratio of the true sizes of
these quadtree boxes is bounded by a constant k3 from Section 5.3, this will affect
the min angle of our triangulation by at most a constant factor.
Level. We now define the level of a triangulation edge so that the smallest level
edges have level 1, the second smallest level 2, the third smallest level 3, etc. That
is, the edges arising from the ith smallest quadtree boxes have level i. Formally,
the level of an edge C in F is log2(sq(G)), where s = 2/minx??q(X). We order
the edges from the leftmost vertex of F to the rightmost vertex of F, and place
them in a list S
The algorithm for drawing the squares is iterative, having steps 1,2,... n, where
as before n = ?log(L(E)/d)J, and is the number of parallel guide lines Mi. At the
i step we draw squares around edges of level at most i. For each square, its edge
parallel to and furthest from F will lie on line Mn?j+1. This edge G is then
considered as an edge of level i + 1 for consideration by the next iteration step.
We say that the square is also of level i + 1. The edges we drew the square around
are consecutive edges in S. We remove that subsequence from S and insert G in
its place.
We now describe exactly which edges we draw squares around in iteration step
143
Figure 5.18: Drawing squares as a preliminary step to triangulating a general
match trapezoid.
i. Identify each maximal length subsequence C of consecutive edges in 5, where
each edge in C has level at most ?. Denote the cardinality of C by m. If m > 1,
then for j=0,1,... [m/2J, we draw a square around the 2jth and 2j + 1th edge of
C. That is, we draw a square around each pair of edges in the subsequence starting
with the first two, and skipping the last if the subsequence is of odd cardinality.
The exception to this rule is if m> 1 and the last edge of C is the last edge of 5.
In that case, the last square we draw is drawn around three edges, namely edges
m --H 2, m --H 1 and m of C. See Figure 5.18.
This concludes the definition of the square drawing algorithm. We now consider
the shape of the squares we have drawn.
Elongation. We define the elongation of a triangulation edge to be 1. We
define the elongation of a square b (and its corresponding edge placed in 5) to be
?GEb21(G)
e(b) --H			21(b)			`
where G ranges over all triangulation edges contained in b, and 1(x) denotes the
level of x. If all the subsequences encountered by the algorithm were of even length,
144
then each square would have elongation exactly 1 -
It is easy to show that the elongation of b may also be expressed as
e(b) = ? e(H)21(H)?!(b),
HEb
where H ranges over the edges that b is immediatly drawn around, both trian-
gulation and square edges. It is this formula that will be used in the lemmas
below.
We seek constant upper and lower bounds on the elongation of a square. We
first derive the upper bounds.
Lemma 38 Any edge not containing the rightmost vertex of F has elongation at
most 1.
Proof. We prove this by induction. Before the first iteration step begins, the
lemma is true by definition of elongation for a triangulation edge. Consider a square
b drawn in the ith iteration step. Now b is drawn around two edges. Moreover,
neither of these edges contains the rightmost vertex of F, and hence by induction
have elongation at most 1. Thus the elongation of b is at most 1/2 + 1/2 = 1. 1
Lemma 39 Any edge containing the rightmost edge of F has elongation at most
2.
Proof. Again we prove this by induction. Before the first iteration step begins, the
edge containing the rightmost vertex of F has elongation 1 by definition. Before
the ith iteration step, there is exactly one edge in S that contains the rightmost
edge of F. If we do not draw a square around this edge in the ith iteration step,
then after this step this edge is a triangulation edge of level at least i and has
145
elongation 1. If we do draw a square around this edge, then it is drawn around at
most three edges, each of which can have level at most i. Thus from the previous
lemma the square has elongation at most 1/2 + 1/2 + 2/2 = 2. 1
We now derive the lower bounds. First we need the following two lemmas about
the algorithm.
Lemma 40 If an edge of level i is not drawn around in iteration step i, then it is
a trian?ulation edge.
Proof. If an edge 0 of level i is not drawn around in iteration step i, then it
is adjacent in 5 to an edge H of level i + 1. Since i + 1 is greater than the
current iteration step, i, H is a triangulation edge. By the balance condition,
all triangulation edges adjacent to an edge of level i + 1 have level at least i. A
triangulation edge adjacent to H is contained in 0, hence 0 contains a triangulation
edge of level i. But by assumption 0 is a level i edge, and so it cannot contain any
other triangulation edges, and is a single triangulation edge. 1
Lemma 41 An edge of level i is drawn around on iteration step i or i + 1.
Proof. Suppose an edge 0 is not drawn around on iteration step i. We have two
cases.
The first case is if 0 does not contain the rightmost vertex of F. Then 0's
successor H in 5 is of level greater than i. By the previous lemma 0 is a single
triangulation edge, and H is a triangulation edge of level exactly i + 1. Hence
0 will be in a maximal subsequence in iteration step i + 1, but will not be the
leftmost edge in the subsequence. So in iteration step i + 1 the algorithm will draw
a square around 0.
146
The second case is if G does contain the rightmost vertex of F. The only way
we would not draw a square around such an edge is if the maximal subsequence
containing it has length 1. But then by the same arguments as the last paragraph,
it is preceded by a triangulation edge H of level i +1. Hence G will be in a maximal
subsequence of length at least 2 in iteration step i + 1. So in iteration step i + 1
the algorithm will draw a square around G. 1
We may now prove a lower bound on the elongation of an edge.
Lemma 42 Any ed?e has elongation at least 1/2.
Proof. We prove this by induction. The lemma is trivially true before the first
iteration step. Consider an edge H drawn around by a square at iteration step i.
From the previous lemma, H may be of level i or i --H 1.
Suppose H is of level i. By induction H has elongation at least 1/2, so
e(H)2I(H)?l(b) > 1/4.
Suppose H is of level i --H1. But by Lemma 40, such an edge H is a triangulation
edge and must have elongation 1. Hence e(H)2l(H)?l(b) = 1/4. Every square is
drawn around at least two edges, so the square will have elongation at least 1/4 +
1/4 = 1/2. I
We now describe how the algorithm is modified when the match edges are not
perpendicular to E and F: If one vertex of the square lies on a match edge, we use
the match edge for a side of a "square.,,
After we describe how to triangulate the squares, we will wish to prove that
the min angle of the triangles is at least a constant factor times A, the optimal
min angle of any triangulation of Q. To that end, we have the following two
147
lemmas concerning the angle that a match edge may make with an R edge F that
it intersects.
Lemma 43 The angle between a match edge and F is at least r/4 + kA, where k
is a constant and A is the optimal min angle of any triangulation of Q.
Proof. Suppose we are in Case 2 of the construction of R. Theorem 27 gives an
upper bound of 3r/2 on the interior angle at a vertex of E. From Case 2 of the
construction of R in Section 5.2, the angle between a match edge and E is at most
?34T --H 4 = T2. Hence the angle that F makes with a match edge is at least the
supplement of this, or
If we are in Case 1 of the construction of R the situation is not as nice. Here the
angle at E may be as great as 3r/4--H k1?(E')/?(E) for some constant k1, where E'
is the other R edge intersecting E at the vertex under consideration. The ratio of
?`s is bounded by a constant times the optimal min angle, or k2A. The angle that
F makes with a match edge is at least the supplement of the above, or r/4+k1k2A.
Lemma 44 The angle between a match edge and F is at most r --H kA.
Proof. If we are in Case 2 of the construction of R, we have that the angle between
E and the match edge is at least 3r/8 --H = 7r/8. Hence the angle between F
and that match edge is at most the supplement of this, or 7r/8.
If we are in Case 1 of the construction of R the situation again is not as nice.
Here the angle at E may be as small as k1?(E)/?(E') for some constant k, where
E' is the other 1? edge intersecting E at the vertex under consideration. The ratio
of ?`s is bounded by the optimal min angle times a constant, that is k2A. The
148
Figure 5.19: Triangulating the squares with the help of a center vertex.
angle that F makes with a match edge is at least the supplement of the above, or
ir--Hk1k2A. 1
Triangulation of squares
We now describe how to triangulate the squares. We triangulate between E and
M1 as in the special case described above.
Center vertex. We then triangulate each square by taking the convex hull of
the center vertexof the square and each edge of the boundary of the square as in
Figure 5.19. For squares that were drawn around two edges in 5 we may chose one
of the vertices of these edges for the center vertex. In particular, if the rightmost
edge is of level i, then we use its leftmost vertex as the center vertex. Otherwise,
the leftmost edge is of level? and we use its rightmost vertex as the center vertex.
For the squares drawn around three edges in 5, it is easy to prove that there are
six cases depending on whether these three edges are triangulation edges or edges
of a drawn square. These squares are triangulated as illustrated in Figure 5.20.
The center vertices in squares without any squares inside them is the midpoint
149
A
y
Th
M
Figure 5.20: How to triangulate a square drawn around three edges.
of a horizontal segment drawn halfway between the bottom and top edges. Note
that the leftmost and rightmost vertical edges may have vertices in their interior,
arising from smaller adjacent squares. In that case we introduce an edge from that
square vertex to the opposite vertex of the triangle that contains it in Figure 5.20.
A special situation can arise for the first square drawn around a triangulation
edge containing a vertex of F. For two trapezoids T1 and T2 sharing a match edge
F, the number of levels n may be different. Moreover, the level of the edges G1
and G2 containing a vertex of F may be different for T1 and T2. Thus the first
square b1 drawn around G1 may be smaller than the first square drawn around
G2. That is, the side of b1 that is contained in F may have vertices in its interior,
arising from the vertices of squares drawn in T2. We triangulate as in Figure 5.21
when b1 is drawn around only triangulation edges. In general, we triangulate by
introducing edges between the match edge vertices and the central vertex. The
extension to squares drawn around three edges is analogous. We show in Lemma
150
mJ
Figure 5.21: There may be many (1/A) square vertices in a match edge for constant
shaped squares, or at most a constant number of square vertices for wide and short
(1/A) squares.
45 that there are strict bounds on the distance between the match vertices and the
central vertex, so that the min angle of triangles in b1 is at least a constant factor
times A.
Min angle of triangles in squares
We now set out to prove that our triangulation of the squares is within a constant
factor of what is optimal.
Lemma 45 Any triangle of a square without a match edge for a side and not
intersecting E has min angle at least a constant times ?(E)/L(E).
Proof. Consider the triangles that contain an edge G parallel to F. Since elonga-
tion is between 1/2 and 1 for G, we have k1L(E)2'(0)--H?--H' < L(G) < k1L(E)2!(G)--Hn
The vertex V opposite G makes altitude at least d21(G)-n--H3, and no more than
d2I(G)--Hn. These bounds follow from the bounds on elongation. The exact altitude
value depends on whether the square has level 1(G), 1(G) + 1, or 1(C) + 2. In all
cases this altitude lies inside the triangle. llence the angle of the triangle at the
vertices of 0 is at least
d21(G)--Hn--H3			d
= arctan
k1L(E)21(G)--Hn			8k1L(E)'
arctan --H____________
151
which is within a constant factor of ?(E)/L(E). Also the angle of the triangle at
V is at least
k1L(E)21(G)--Hn--H1			k1L(E)
arctan			d21(G)--Hn			= arctan 2S(E)
which is bounded below by constant.
Consider a triangle that contains an edge G perpendicular to F. Let b be the
square containing triangle. The length of G is within constant factors of d2l(b)--Hn--H2.
The altitude at the vertex opposite G is within a constant factor of k1L(E)2l(b)--Hn--H1,
and hence the triangle has min angle at least arctan(?(E)/L(E)), which completes
the proof. I
Lemma 46 Any triangle intersecting E has min angle at least a constant times
Proof. As shown before, M1 is subdivided into at most a constant number of
triangulation edges, where the constant factor depends linearly on k1. The distance
between E and M1 is d/2. But the lower bound on the angle that a match edge
makes with E may be written in terms of ?(E)/L(E). Hence the length of M1 is
at least a constant times the length of E. The constant upper bound on the angle
that a match edge may make with E, and the fact that d < L(E)/2 insures that
L(M1) < 2L(E). Thus from the construction each triangulation edge in M1 has
length at least ciL(E) and at most c2L(E). A lower bound on the worst triangle
that may be formed is shown in Figure 5.22. Thus any triangle has min angle at
least
arctan(			--H arctan(1+12ciL$))
Since d/L(E) <1/2, we have that the above is at least O.8d/L(E). I
152
c1L(E)12
L(Efl
Figure 5.22: The worst case of a triangle containing a vertex of E.
We now consider triangles in squares that are on the boundary between two
trapezoids.
Lemma 47 Any triangle in a square with a match edge for a side has min angle
at least a constant times A.
Proof. On the boundary between two trapezoids, the subdivision of the match edge
for one trapezoid may not correspond to the subdivision for the other. However,
from Lemma 40 we see that this can only happen inside the first squares drawn
around the R triangulation edge containing the R vertex of the match edge.
For triangles of a square that is not drawn around the triangulation edge con-
taining the match edge R vertex, there is no vertex in the side of the square.
Hence the arguments of the previous two lemmas are sufficient to show that the
triangulation of the square is at least a constant factor times A.
We now consider the squares drawn around the triangulation edge containing
the match edge R vertex. The algorithm to triangulate R centers vertices of R
in a quadtree box. That is, the two triangulation edges G and G' containing a R
vertex V and contained in R edges F and ?? are contained in the same quadtree
box. Thus the square drawing algorithm considers q(G) and q(G') to be within a
constant factor of one another. We have two cases.
153
The first case is if the lengths of E and E' are within constant factors of one
another, then the levels of & and C' are within constant factors of one another.
Hence from Lemma 41 the square drawn around G has at most a constant number
of vertices in its match edge side arising from squares around F'. Similarly for CI.
The second case is if there are more than a constant number of subdivision
vertices in its match edges, we wish to show two things. First we wish to show
that the box has constant height to width ratio in terms of Euclidean distance.
Second we wish to show that L(F)/L(F') is bounded below by A. Thus we may
triangulate as in Figure 5.21 and achieve min angle at least A. Note that in the
case that the match edge makes small angle with F and all the edges drawn around
are triangulation edges we would instead introduce edges from the match vertices
to the lower right vertex of E in Figure 5.21.
We now provide the details for the second case. Edges F and F' arising from
shrinking edges of the same smallest spiral have length ratio bounded by a constant.
This follows from the constant bound on the ratio between similar triangles of the
same smallest spiral. So it remains to consider edges arising from shrinking two
edges E and E' that were used to weld two spirals. Let us consider the case that
E = XIX" was used to weld two spirals for the same vertex. We have already
shown in Theorem 28 that L(E)/L(E') > kA(QC) for some constant k. However,
in this case E is at most a constant times the lengths of the two spokes containing
its vertices, which are equal in length. Also, these spokes are at most one quarter
of ?(V), where V is the P vertex they contain. Hence 6(E) is at least a constant
factor times L(E). We now consider the case that E' was used to weld two spirals
for the same P edge, and E is the first rim edge of a spiral at V. Again from
154
Theorem 28 we have L(E)/L(E') > kA(QC) for some constant k. However, the
only case that L(E)/L(E') is less than a constant is when the smallest optimal
spiral consists of one triangle. That is, when the interior angle at V is small and
the shrinking rate of the smallest spiral is not too small in relation to it. But as in
the previous welding case, then E is between two spokes of (nearly) equal length,
and hence ?(E) can be no less than a constant times L(E). I
5.4.2 Triangulating a match triangle
Triangulating a match triangle T is no more difficult that triangulating a match
trapezoid. The key observation is a bound on the lengths of the two match edges
of a triangle T. Let V be the Q vertex of T. Let F be the R edge of T opposite
V. Let H1 and H2 be the R edges intersecting F, and N1 and N2 the match edges
from Q to F. Assume L(N1) > L(N2). Note L(N1) is within a constant factor
of ?(H1). Similarly for N2. Since the vertex of H2 opposite V interferes with
H1, we have that ?(H1) < L(H2). Hence L(N1) < L(H2), and we conclude that
L(N1)/L(N2) <--H L(H2)/?(H2).
Thus, we may proceed as in the trapezoid case. We introduce lines parallel to
F, spaced as before. Recall that the angle of T at V is 7r/2. Hence, by the last para-
graph, the angle between these lines and a match edge is bounded by a constant
times L(H2)/?(H2). Hence we may build squares as in the trapezoid case, with min
angle at a match edge at least a constant factor times min (6(H1)/L(H1), 6(H2)/L(H2)),
which is at least a constant factor times A.
The algorithm used to triangulate the squares is the same as in the trapezoid
case. We triangulate between M1 and V by introducing an edge from W to each
of the square vertices of M1. The proof of Lemma 46 shows that these triangles
155
have min angle at least a constant times A.
5.5 Conclusions
5.5.1 Why the two dimensional analog of Chapter 4?
The details of the method used in Section 5.3 to triangulate R is not important
for the correctness of the other parts of the algorithm. Any two dimensional mesh
generation algorithm that generates triangles with guaranteed min angle or aspect
ratio may be used, as long as it satisfies the two mentioned features:
1. The length of a triangulation edge G contained in an R edge E must be at
least a constant times 6 for that edge.
2. The ratio of the lengths of two intersecting triangulation edges must be
bounded by a constant.
Any triangulation algorithm with bounded aspect ratio would have to satisfy the
second condition, although the constant could be quite large.
The method of Bern, Eppstein and Gilbert[1991J satisfies the second condition,
and provides us with a quadtree to guide Section 5.4, if the input has no small
interior angles that need to be cut off. This is the case for R because of the lower
bound on the smallest interior angle given by Theorem 26. The first condition
is not satisfied by their algorithm for general input. Normally, the duplication of
boxes in Chapter 4 is needed to guarantee this. The method is akin to Reimann
sheets in Seidel[1988J. llowever, the fact that we shrunk Q to R by moving edges
a constant fraction of 6 towards the interior of R shows that generating small
triangles in one part of R will not necessarily generate small triangles in another
156
part of R close in Euclidean distance but far in geodesic distance. Hence Bern,
Eppstein and Gilbert [1991] would satisfy the first condition for R.
The algorithm of Rupert[1992J will also satisfy both conditions. This is startling
because it is based on Delaunay triangulations and is not quadtree based. That
algorithm recursively divides edges in half, so that the ratio of the lengths of
two intersecting triangulation edges contained in an edge of R is at most 2. If
Rupert[1992J was used, the algorithm presented in Section 5.4 would have to be
modified to handle the fact that there is no quadtree to guide us.
5.5.2 Running time
We now provide an upper bound on the running time of our algorithm, in the
process describing how we identify the smallest optimal spirals and compute the
delta function.
Let n denote the number of edges of P. We start our algorithm by computing
the mediat axis of P in optimal time O(n log n). See Duda and Hart [1973] and
Kirkpatrick [1979]. This allows us to identify on average a constant number of faces
that have points that could define the smallest optimal spiral for a given vertex
and edge. Thus, after the medial axis is computed, we may identify the smallest
spirals for all edges and vertices of P in total time 0(n).
Producing the triangles of the wedges about a vertex (QC) takes time linear
in the number of triangles produced. The medial axis of Q may be used to effi-
ciently compute the delta function for the faces of Q. This may take time at most
O(m log m), where m is the number of faces of Q. After the medial axis of Q is
computed, R takes time 0(m) to produce. Note that m > n, and m may be as
large as on the order of the cardinality of the final triangulation.
157
Let t be the time to triangulate R, and ? the number of triangulation edges
on the boundary of R. Using geometric series, it is easy to prove that triangulat-
ing between Q and R produces O(?) triangles. Triangulating between Q and R
takes time linear in the number of triangles produced. To achieve this bound in
the triangulation algorithm of Section 5.4, at each iteration level we use the pre-
vious iteration's maximal subsequences to guide our determination of the current
iteration's maximal subsequences. Thus we have the following.
Theorem 31 The running time of the atgorithm is O(? log ? + t), where ? is the
number of edges of the finat triangulation, and t is the time taken to triangutate R.
As the algorithm is described, using the two dimensional analog of Chapter 4
to triangulate R, t is no worse than o(?2 log ?). So, t is the limiting factor in the
running time of this algorithm. See Section 4.9 for more comments on t.
5.5.3 Cardinality
The cardinality of the triangulation is not optimal. The construction of the wedges
around the vertices of P triangulating QC produces an optimal number of triangles,
up to a constant factor. However, triangulating R with any algorithm designed to
maximize the minimum angle, allowing Steiner points arbitrarily in R, is bound to
produce too many triangles for certain input.
5.5.4 Open problems
Concerning triangulations with only interior Steiner points, there are several other
measures for which algorithms may be useful. As mentioned in the introduction,
maximizing the minimum height would be useful for the three dimensional mesh
158
generation algorithm. The algorithm as described bounds sphere ratio in light of
Theorem 1 in Chapter 3, but is not within a constant factor of optimal. Other
possible measures include minimizing the maximum angle.
Is there an algorithm for triangulating a three dimensional polytope with max-
imum minimum angle? The difficulty in extending this work to three dimensions
is that we must define an optimal spiral at a vertex of possibly high degree. The
boundary of a three dimensional spiral may be a complicated surface, and not a
simple curve as presented here.
References
F. Aurenhammer [1991], Voronoi diagrams- a survey of a fundamental geometric
data structure. ACM Computing Surveys 23, 345-405.
I. Babuska and A. K. Aziz [1976], On the angle condition in the finite element
method, SlAM J. Num. Anal. 13:214--H226.
B.S. Baker, E. Cross and C.S. Rafferty [1988], Non-obtuse triangulation of poly-
gons. Disc. and Comp. Ceom. 3 (1988)147-168.
M. Bern, D. Dobkin, and D. Eppstein [1991], Triangulating polygons without large
angles, preprint.
M. Bern, H. Edelsbrunner, D. Eppstein, 5. Mitchell, and T. 5. Tan [1991], Edge-
insertion for optimal triangulations, preprint.
M. Bern, D. Eppstein, J. Gilbert [1990] Provably Good Mesh Generation, Xerox
Palo Alto Research Center. Also, Proc. 1990 Symposium on Foundations of
Computer Science.
M. Bern and D. Eppstein [1991], Mesh Generation and Optimal Triangulation,
preprint.
E. K. Buratynski [1990], A fully automatic three-dimensional mesh generator for
complex geometries, Int. J. for Num. Meth. in Engr. vol 30 931-952.
G. Carey, M. Sharna, and K. Wang [1988], A Class of Data Structures for 2-D and
3-D adaptive mesh refinement, Int. J. Numer. Meth. in Engr. 26:2607--H2622.
B. Chazelle, L. Palios [1989], Triangulating a Nonconvex Polytope, Proc. 5th ACM
Symposium on Computational Geometry. 393-399.
L. P. Chew [1989a], Constrained Delaunay Triangulation, Algorithmica 4:97--H108.
159
160
L. P. Chew [1989bJ, Guaranteed-Quality Triangular Meshes, C. 5. Cornell, TR
89--H983.
B. Delaunay [1934], Sur la sphere vide, Bull. Acad. Sci. USSR (VIII), Classe Sci.
Mat. Nat., 793-800.
R. Duda and P. Hart [1973], Pattern Classification and Scene Analysis, Inter-
science, New York.
II. Edelsbrunner [1987], Algorithms in combinatorial geometry, Springer Verlag,
Berlin.
J. Hauser and C. Taylor [1986], Numerical Grid Generation in Computational Fluid
Dynamics (Proceedings), Pineridge Press, Swansea, U.K.
D. Kirkpatrick [1979], Efficient computation of continuous skeletons, Proc. 20th
lEFE Annual Symp. on Foundations of Comput. Sci., pp. 18-27, Oct. 1979.
C.L. Lawson [1977], Software for C1 surface interpolation, in J. Rice, ed., Mathe-
matical Software III, Academic Press 161-194.
D.T. Lee [1978], Proximity and reachability in the plane. Coordinated Science
Laboratory, Univ. Illinois, Tech.Rep. R-831.
D.T. Lee and A. Lin [1986], Generalized Delaunay triangulation for planar graphs.
Disc. and Comp. Geometry 1 (1986), 201-217.
G. L. Miller, S.-ll. Teng, and 5. A. Vavasis [1991], A Unified Geometric Approach
to Graph Separators, Proc. 32nd Symp. Foundations of Comp. Sci.
G. L. Miller and W. Thurston [1990], Separators in Two and Three Dimensions,
Proc. 22nd Symp. on Theory of Computing, 300--H309.
G. L. Miller and 5. A. Vavasis [1990], Density Graphs and Separators, Department
of Computer Science, Cornell, Tech Report 90-1169, and to appear, Symp?
sium on Discrete Algorithms.
W. F. Mitchell [1987], A Comparison of Adaptive Refinement Techniques for El-
liptic Problems, Department of Computer Science, University of Illinois at
Urbana-Champaign, Report No. UIUCDCS-R-87- 1375.
W. F. Mitchell [1988], Unified Multilevel Adaptive Finite Element Methods for
Elliptic Problems, Department of Computer Science, University of Illinois at
Urbana-Champaign, Report No. UIUCDCS-R-88-1446.
161
5. A. Mitchell and 5. A. Vavasis [1992], Quality Mesh Generation in Three Dimen-
sions, Proc. 8th Symp. on Computational Geometry, 212-221.
D. Moore [1992], Simplicial Mesh Generation with Applications, C.S. Cornell, PhD
thesis, May 1992.
D. Moore and J. Warren [1990a]. Adaptive mesh generation I: packing space. Dep.
of Computer Science, Rice University, Tech. Report TR 90-106.
D. Moore and J. Warren [1990b], Multidimensional Adaptive Mesh Generation,
Department of Computer Science, Rice University, Rice COMP TR 90-106.
F.P. Preparata and M.I. Shamos [1985], Computationat Geometry: An Introduc-
tion, Springer Verlag, New York.
V.T. Rajan [1991], Optimality of the Delaunay Triangulaion in Rd. In Proc. 7th
ACM Symp. Comp. Geometry 357-363.
J. Rupert [1992], A New and Simple Algorithm for Quality 2-Dimensional Mesh
Generation, Computer Science Division (EECS), UC Berkeley, Report No.
UCB/CSD 92/694.
N. Sapidis and R. Perucchio [1989], Advanced techniques for automatic finite ele-
ment meshing from solid models. Computer Aided Design, 21 4 248-253.
E. Sch5nhardt [1928], Uber die Zerlegung von Dreieckspolyedern in Tetraeder.
Math. Ananalen 98 309-312.
R. Seidel [1988], Constrained Delaunay triangulations and Voronoi diagrams with
obstacles, Technical Report 260, Inst. for Information Processing, Graz, Aus-
tria.
