数理解析研究所講究録
1068 巻 1998 年 110-120

110

Optimality and Integer Programming Formulations of
Triangulations in General Dimension
Akira
Hiroshi$\mathrm{T}\mathrm{a}\mathrm{j}\mathrm{i}\mathrm{m}\mathrm{a}^{*}$$\mathrm{I}\mathrm{m}\mathrm{a}\mathrm{i}^{\uparrow}$($\uparrow \mathrm{t}\mathrm{a}\mathrm{j}\mathrm{i}\mathrm{m}\mathrm{a}@\mathrm{t}\mathrm{r}1.\mathrm{i}\mathrm{b}\mathrm{m}.\mathrm{c}\mathrm{o}$(imai@is.

$\mathrm{s}.\mathrm{u}$-tokyo.

. jp)$\mathrm{a}\mathrm{c}$. jp)

Abstract
The properties of triangulations in two and three dimensions are computationally in. Three IP formulations of triangulations are
vestigated by using integer programming
introduced, two based on the stable set problem, and the other based on the set partitioning
problem. Some properties that are interesting from a theoretical or practical point of vieware considered as objective functions for . Finally, some computational results are given.
This approach allows three-dimensional triangulations to be treated in a flexible andefficient way, and has led to the discovery of some interesting properties of three-dimensional
triangulations.
$(\mathrm{I}\mathrm{P})$$\mathrm{I}\mathrm{P}$1

Introduction

Triangulation is a major topic in computational geometry, and there have been many studies ofvarious approaches to triangulating a configuration of points, such as Delaunay triangulation, theminimum weight triangulation, and the lexicographic triangulation. Most of these studies havefocused on triangulations in two dimensions, and less attention has been paid to triangulations
in higher dimensions.
As a result of the recent advances in the performance of computers, however, the number ofapplications using triangulation in three dimensions is growing. For the finite element method(FEM), a three-dimensional polyhedron has to be divided into meshes. There are many meshingtechniques: one approach is to distribute points inside the polyhedron and generate meshes by
triangulating them. This approach is gaining popularity because of its simplicity and independence of dimensions [14]. Another major application is volume rendering, which is a method forvisualizing the results of FEM, or semi-transparent objects such as clouds and flames. Unlike
ordinary computer graphics methods that visualize the surface of materials, volume rendering
meshes the three-dimensional space [15].
To contribute to these applications, the properties of triangulations in three dimensions
should be further investigated. In three dimensions, some properties that are valid in two
dimensions do not hold. Bistellar flips (Figure 1) show that the number of triangles is notconstant in three dimensions. Sch\"onhardt’s polytope (Figure 2), which can be obtained bytwisting a prism in such a way that the side faces become concave, shows that not all polyhedra
can be triangulated in three dimensions. Some results on triangulatability are given in [13].
In two dimensions, Delaunay triangulation is optimal in many respects; for example, itmaximizes the minimum angle, and minimizes the largest circumsphere and the largest minimum
enclosing sphere. Many applications, such as meshing algorithms for FEM, depend on this*Tokyo Research Laboratory, IBM Research.
\dagger Department of Information Science, University of Tokyo.

111

.
$1;:$

’
$;$

,

;

$|$

’
$|4$

$\int)$$1;\mathfrak{l},$‘$!1;;;;\mathrm{t};$Figure 1: Bistellar flips in 2 and 3 dimensions

Figure 2: Sch\"onhardt’s polytope

characteristic. But in three dimensions, it only minimizes the largest minimum enclosing sphere[12], and does not necessarily lead to optimality in other respects. In other words, Delaunay
triangulation is not so good in three or more dimensions. Thus, it is necessary to obtain a good
triangulation for each individual application, namely, to optimize the triangulation.
In this paper we investigate the optimality of triangulations in two and three dimensions
formulations. We first give three formulations in general
by using integer programming
dimension. Two of them are based on the stable set problem, and the other is based on theset partitioning problem. We then compare these formulations in size, and in other respects
related to the efficiency of the procedure for solving the problem. We consider several objectivefunctions in both two and three dimensions, which are interesting from a theoretical or practicalpoint of view. We also give some computational results. Through the experiments, we foundsome interesting properties of triangulations in three dimensions.
As a preliminary step, we define some special terms. A -simplex is a -dimensional polytopethat is the convex hull of $d+1$ affinely independent points. For example, a line segment, atriangle, and a tetrahedron correspond to a 1-simplex, a 2-simplex, and a 3-simplex, respectively.
An -face of a -simplex is an -simplex $(0\leq\dot{i}\leq d)$ that is the convex hull of a subset of thevertices of the -simplex. In particular, a $(d-1)$ -face is called a facet. Two -simplices intersectwhen the intersection is non-empty and is not a face of at least one of the two simplices. In
this paper, especially for IP formulations, we consider the division of the convex hull of a point
configuration $A$ of points in -dimensional space using -simplices, but we use the term
$triangulati_{\mathit{0}}n$ for convenience. We assume that $A$ is a configuration in general position (no $d+1$
points lie on the same $(d-1)$ -dimensional hyperplane).
This paper is organized as follows. In Section 2, we give two kinds of IP formulationsof triangulations in general dimension, then compare and investigate them. In Section 3, weconsider objective functions. In Section 4, we give some computational results. Section 5summarizes this paper.
$(\mathrm{I}\mathrm{P})$

$d$

$d$

$\dot{i}$$d$$\dot{i}$

$d$

$d$

$n$

2

$d$

$d$

Formulations

In this section, we give IP formulations by regarding triangulation as the stable set problem(Formulations and 1), or the set partitioning problem (Formulation 2). Especially for theformer, we consider two cases: that in which each -simplex corresponds to an element (Formulation ), and that in which each -simplex $(0\leq i\leq d)$ corresponds to an element (Formulation
$0$

$d$

$0$

1).

$\dot{i}$112

Formulation : the stable set problem of d-simplices

2.1

$0$

corresponds to the set of -simplices, andIn the intersection graph $G(V, E)$ of -simplices,an arc is defined between two nodes when the corresponding two -simplices intersect. A triangulation has no pair of simplices that intersect and thus corresponds to a stable set; further, it
is maximal, because another -simplex surely intersects with some -simplices that constitutethe triangulation. On the basis of this observation, Imai [10] gave an algorithm for enumerating
triangulations in general dimension.
As mentioned in Section 1, not all polytopes can be triangulated in dimensions higher than
two. Thus, there exist cases in which regions like Sch\"onhardt’s polytope remain untriangulated
inside the convex hull, and the maximal stable set is not a triangulation. Therefore, we have toensure that the stable set becomes a triangulation by imposing the condition that the sum ofof the convex hull.
of -simplex is equal to the volume
the
$V$

$d$

$d$

$d$

$d$

$d$

$\dot{\mathrm{v}}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}v_{i}$$d$minimize$\mathrm{s}.\mathrm{t}$$V$$\dot{i}$(1)$\sum_{i}c_{l}Xi$.
$x_{i}$$x_{i}+x_{j}$$\in$$\{0,1\}$$\leq$

1 (for$=$$i,$$j\mathrm{s}.\mathrm{t}$. -simplex and -simplex$d$$\dot{i}$

$d$

$j$

intersect)$V$$\sum_{i}v_{i}x_{i}$Here, checking the intersections among -simplices is essentially redundant. For example, consider triangulations of the convex hull of four points on a plane. We have four triangles, andfour pairs of triangles that intersect. On the other hand, the four intersections are derived froman intersection of two diagonals. In other words, if we explicitly handle lower-dimensional facesof -simplices, the formulation can be more efficient.
$d$

$d$

2.2

Formulation 1: the stable set problem of -simplices $(i\leq d)$$\dot{i}$In the two-dimensional case, Kyoda et al. [11] formulated the minimum weight triangulation(MWT) problem as the stable set problem on the intersection graph of edges (1-simplices). MWT
is a famous problem for which it is not known whether a solution can be obtained in polynomial. In two dimensions, the maximal stable sets of edges are always maximum, and thenumber of edges $M$ is constant. Thus MWT is obtained by using the following formulation:$\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}[9]$minimize$\mathrm{s}.\mathrm{t}$(2)$\sum_{\dot{*}}\mathrm{q}xi$.
$x_{i}$$x_{i}+x_{j}$$\in$$\{0,1\}$$\leq$

1 (for$=$$M$$\dot{i},j\mathrm{s}.\mathrm{t}$. edge and edge$\dot{i}$$j$intersect)$\sum_{i}x_{i}$2.2.1

Extending to higher dimensions

The above formulation (2) is valid only in two dimensions, and the cost and variables are assigned
only to edges, so we need to extend it to bigher dimensions and assign variables to d-simplices.
We first consider the simplest example in three dimensions, a bistellar flip defined on the convex
hull of five points. In the lattice-like structure in Figure 3, the -th layer corresponds to $(\dot{i}-1)-$
simplices, and straight arcs among layers correspond to face relations of simplices in different$\dot{i}$113

dimensions. For example, edges $(0,1)$ , $(0,4)$ , and $(1,4)$ are facets of, and necessary for triangle(0,1,4). Thus, the lattice is interpreted as a poset when we assign a variable to each simplex
in such a way that the variable is equal to 1 if the simplex appears in the triangulation, andotherwise. We will call the relations defined by the poset lattice constraints.
The lattice itself is independent of the configuration, whereas information on intersections
depends on the configuration. In Figure 3, edge $(0,4)$ and triangle (1,2,3) intersect. On theother hand, point (4) and tetrahedron (0,1,2,3) intersect in Figure 4.
only one of them can appear in the trianguWhen edge $(0,4)$ and triangle (1,2,3)
lation. This can be represented as a constraint that the sum of the corresponding variables isless than or equal to 1. The following condition is necessary to make the formulation efficient:$0$$\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{t}\rangle$

Proposition 1 All the vertices of simplices are assumed to be in general position. Then, two-simplices intersect if and only if they have a -face and a -face, respectively, that intersectand satisfy $k+j=d$ .
$k$

$d$

$j$

Proof: Owing to the space limitation, this proof is omitted and will be given in a full paper.
We obtain a set of -simplices that do not intersect by imposing the intersection constraints
only on the pairs of simplices whose sum of dimensions is equal to , and the lattice constraints.
The volume constraint is again necessary, as in formulation (1), for the set to be a triangulation.
$d$

$d$

minimize$\mathrm{s}.\mathrm{t}$(3)$\sum_{i}c_{?}x_{i}^{d}$.
$x_{i}^{k}$$\in$$\{0,1\}(0\leq k\leq d)$$x_{i}^{k}+x_{j}^{d-k}$$\leq$

1 (for$\dot{i},j\mathrm{s}.\mathrm{t}$$x_{i}^{kk1}-x_{j^{+}}$$\geq$$0$(for$i,j\mathrm{s}.\mathrm{t}$

$=$

$V$

. -simplex and $(d-k)$ -simplex intersect, $0\leq k\leq d$ ). -simplex isa facet of $(k+1)$ -simplex $0\leq k\leq d-1$ )$k$$\dot{i}$$k$$\dot{i}$

$j$

$j,$

$\sum_{i}v_{i}x_{i}^{d}$2.3

Formulation 2: the set partitioning problem

We use the term chambers for the minimal cells obtained by dividing the convex hull with all
the possible $(d-1)$ -simplices. Obtaining a triangulation can be treated as an instance of the set
partitioning problem by regarding each chamber as an element, and each -simplex as a subset
[1].
$d$

minimize$\mathrm{s}.\mathrm{t}$(4)$cx^{d}$.
$x_{i}^{d}$$Ax^{d}$

2.3.1

$\in$

{0,1}

$=$

1 $A_{ij}=\{$

1$0$-simplex
otherwise$d$$j$contains chamber$\dot{i}$Cocircuit form constraints

A large amount of geometric computation is necessary for handling chambers, and therefore
be a
we consider another type of constraint called cocircuit form constraints instead. Let
containing . Let
are defined by a hyperplane$(d-1)$ -simplex. Two half-spaces$e_{k}$$\mathcal{H}_{e_{k}}^{+},\mathcal{H}_{e_{k}}^{-}$$\mathcal{H}_{e_{k}}$$e_{k}$

114

Figure 3: Face relations among simplices

Figure 4: Intersection of a -simplex and a 3-simplex
$0$

115

defined by an element
be a -simplex
constraint holds for all $(d-1)$ -simplices [7].
$d$$e_{k}\cup\{a\}$$t_{i}$$a$of the point set $A$ and$\pm 1$$t_{i}=e_{k^{\cup\{a}} \}\sum_{\in aA\cap \mathcal{H}_{e}^{+}k},xi-\sum_{=tiek\cup\{a\},a\in A\cap \mathcal{H}_{e_{k}}^{-}}X_{i}=\{$$0$$e_{k}$. The following

is on the boundaryof the convex hullotherwise$e_{k}$For $(d-1)$ -simplices on the boundary, the above equation is the same as the set partitioning
constraint for the adjacent chamber. For interior $(d-1)$ -simplices, it corresponds to the differenceof the constraints for the chambers on both sides.
De Loera et al. [7] showed that the cocircuit form constraints are sufficient to define theaffine hull of the characteristic vectors of -simplices that form a triangulation. From now on,
we will call the formulation using the cocircuit form constraints Formulation 2.
$d$

2.4

Comparison of formulations

In this section, we compare Formulations $0,1$ , and 2. Table 1 shows the size of each formulation.
Although Formulation 1 uses fewer constraints than Formulation , Formulation 2 is still themost compact.
$0$

2.4.1

Polytopes

We denote the polytope that corresponds to the linear programming relaxation of formulation ,to the subspace that corresponds
respectively. The projection of, and
1, and 2 by
to the variables for -simplices is denoted by . Then, the following relations hold among thesepolytopes:$0$$P_{0},$$P_{1}’$$P_{2}$$P_{1}’$$d$$P_{1}$

Proposition 2 $P_{0}=P_{1}$ .
Proof: Owing to the space limitation, this proof is omitted and will be given in a full paper.

Lemma 1 For the convex hull of $d+3$ points in dimensions, there exist two triangulations
inside the convex hull andthat share no -simplex. Further, there exists a -simplex that
does not belong to either of the two triangulations.
$d$

$d$

$d$

$l_{i}es$

Proposition 3 $P_{2}\subset P_{0}(n\geq d+3)$ .
Proof outline: By focusing on a chamber, we can observe that both the volume and intersection.
constraints are satisfied if the set partitioning condition is satisfied. Thus
From Lemma 1, we can configure the overlap of two triangulations, both of which are assigned
a weight of . We then perturb it slightly by using another simplex so that it still satisfies the.
constraints of Formulation . It then violates some chamber constraints, and$P_{2}\subseteq P_{0}$$\frac{1}{2}$$0$$P_{2}\subset P_{0}$$\square$

116

2.4.2

Spanning triangulations

All the formulations above allow triangulations that do not use points inside the convex hull.
However, triangulations are often assumed to use all the points: we call them spanning triangulations. In order to avoid non-spanning triangulations, we eliminate non-empty d-simplices
in advance. For random point sets, the expected number of empty -simplices is $O(n^{d})$ [3],which is significantly smaller than the upper bound[3].
the boundary of the convex hull. The lower bound is$d$$(_{d-1})$

2.5

Cutting planes

As cutting planes for the problem of formulation (2), Kyoda et al. [11] applied clique cuts andodd-cycle cuts, both of which are well known for the stable set problem, and convex polygon
cuts which use geometric information on triangulations. Formulations and 1 are based on
the stable set problem, and these cutting planes can cut off fractional solutions of their linear
relaxations. In this section, we investigate the effectiveness of cutting planes withrespect to Formulation 2.
The convex polygon cut is based on the property that the number of edges is constant fortriangulations in two dimensions. Thus, this cut is valid only in two dimensions. Although it
was given as a constraint on the number of edges, it can also be described as a constraint on
the number of triangles, as follows:$0$$\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{g}_{\mathrm{T}}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{g}$denote a configuration of points, $C$ denote its convex hull, and $M$ denotethe cardinality of the spanning triangulations of V. Further let denote the $\dot{i}ncident$ vector of a
, and letdenote thespanning triangulation of a point configuration
vertices.
as
their
the
elements
have
only
to
that
dimensions of , which correspond triangles
ofThe following inequality holds:

Definition 1 Let$\mathcal{V}$

$x$

$V$

$A(A\supseteq \mathcal{V}, A\backslash \mathcal{V}\cap C=\emptyset)$$\mathcal{V}$$x$

Convex polygon cut$\sum_{i\in V}Xi\leq M$

Here we focus on an example in which Formulation 2 has a fractional solution [7]. A set ofvertices 1, . . . , 5 of a regular pentagon and its center constitutes a point set in two dimensions,with coordinates $x\{123\}=$and has 20 triangles and 16 triangulations. The vector$x_{\{2\mathrm{s}4} \}=x_{\{345\}}=x_{\{145\}}=x_{\{125\}}=x_{\{013\}}=x_{\{024\}}=x\{035\}=x_{\{014}\}=X\{025\}=\frac{1}{2}$, andall other coordinates $=0$ , satisfies the constraints of Formulation 2. The subgraph of theintersection graph of triangles induced by the nodes corresponding to the triangles with value
and without vertex , constitutes an odd-cycle of size 5. The corresponding constraint$0$$x\in \mathcal{R}^{20}$$\frac{1}{2}$$0$$x\{123\}+x\{234\}+X\{345\}+x\{145\}+X\{125\}\leq 2$

cuts off . Thus, odd-cycle cuts can be effective cutting planes. On the other hand, the convex
polygon cut $\sum_{ijk}x_{i}jk\leq M(=5)$ is satisfied by and redundant. For convex polygon cuts, westate the following:
$x$

$x$

Conjecture 1 Any convex polygon cut

3

$iS$

valid for

$P_{2}$.

Objective Functions

In this section, we introduce some objective functions that are interesting from a theoreticalor practical point of view. Objective functions can be categorized into two types. One issimple summation, which only requires us to solve IP on the basis of one of the formulations we

117

introduced in Section 2. The other is bottleneck optimization, such as minimizing the maximum.
If we try to solve such problems solely by using , the branch-and-bound procedure just goesprogressively deeper and the lower bound never improves, and thus we cannot obtain a solution.
Binary search is one way to avoid this numerical instab.i.lity.
$\mathrm{I}\mathrm{P}$Minimize the maximum aspect ratio. In applications such as FEM, we should avoid flatness of triangles for the sake of computational stability and accuracy. The most straight, which is the ratio of the radii of theforward index of the flatness is the aspect ratio$(\mathrm{A}\mathrm{R})$circumscribing and inscribing spheres.

Maximize the minimum angle. Delaunay triangulation in two dimensions maximizes theminimum angle. In three dimensions, the angle itself has two varieties; the dihedral angle
and the solid angle. There are several meshing algorithms that avoid some of the badtoo large or too small–angles (see [4], for example), but none for triangulating a point
configuration. We can obtain a solution by specifying the $\min-\max/\max-\min$ angle as theobjective function.
Minimize the number of triangles. The number of triangles varies in three dimensions, andis an interesting objective function related to the following open problem:
Open Problem 1 Is there a polynomial-time algorithm for triangulating an arbitrary
convex polyhedron with the minimum number of tetrahedra ?
$f\mathit{5}f$

Minimize the weight. MWT in two dimensions is an interesting problem, as we mentioned
in Section 2.2. In three dimensions, the weight is extended from the edge length to thesurface area of the triangle.

4

Computational Experiments

This section gives the results of computational experiments based on Formulation 2. All themodel 990, using the IBM Optimization Subroutine
experiments were done on an IBM
Library to solve mathematical programming problems. We also used the library CGAL [6] forhandling geometric objects, and for obtaining Delaunay triangulation, we used the program
for three dimensions.
Triangle [16] for two dimensions, and the program$\mathrm{R}\mathrm{S}/6000$$\mathrm{D}\mathrm{e}\mathrm{W}\mathrm{a}\mathrm{l}\mathrm{l}[8]$4.1Problem size and computational time

We measured the size and the CPU time for obtaining the minimum weight (spanning) triangulation of uniformly distributed point sets in two and three dimensions (Tables 2 and 3). The
size of solvable problems is quite small in three dimensions.
We obtained integer solutions by means of the linear programming relaxations in all thecases, whereas in [11], based on the stable set problem, we obtained fractional solutions in manycases and applied cutting planes and the branch-and-bound method.
We experimented with several kinds of objective functions, including the ones in Section
3, and encountered fractional solutions after solving the relaxed linear programming problems
only when the objective function was the cardinality, that is, when the problem was unweighted.
Formulation 2 is fairly good practically, but still, the cutting planes we mentioned in Section
2.5 would be necessary to solve unweighted cases efficiently.

118

Table 2: Size and CPU time of the minimum weight triangulation in$2\mathrm{D}$Size and CPU time of the minimum weight triangulation in$\mathrm{T}\mathrm{a}\mathrm{b}\mathrm{l}\underline{\mathrm{e}3.}$$3\mathrm{D}$: CPU time for solving the IP problem
Others: CPU time for the rest$\mathrm{I}\mathrm{P}$

4.2

Optimality

Figures 5, 6, and 7 show the cardinality, weight and maximum aspect ratio, respectively, of
triangulations in three dimensions. We used 10 point sets of cardinality 10, 20, and 30 randomly, and AR in these figures represent Delaugenerated in a unit cube with different seeds.
nay triangulation, the minimum weight triangulation, the minimum cardinality triangulation,and the triangulation with the minimum maximum aspect ratio, respectively.
From Figure 5 we can observe that Delaunay triangulation tends to contain more tetrahedra
than the other triangulations. Figure 6 shows that Delaunay triangulation is not at all close
to the minimum weight triangulation in three dimensions. Figure 7 shows that Delaunay triangulation avoids flat tetrahedra. There is still a gap between Delaunay triangulation and the
optimal solution, which is obtained by the triangulation with the minimum maximum aspectratio.
$\mathrm{D},$5$\mathrm{W},$$\mathrm{C}$Conclusions and Remarks

In this paper, we have introduced IP formulations of triangulations of a point configurationin general dimension, and shown through computational experiments that the properties oftriangulations differ significantly between two and three dimensions.
Another very important property of triangulations in three dimensions is that not all polyhedra can be triangulated [13]. Through the use of IP and maximization of the volume that can betriangulated, it may be possible to gain some insight into triangulatability. Among the formulations introduced in this paper, the one using cocircuit form constraints was the most efficientfor triangulations of a point configuration. But it is still an open question which formulation is
the most efficient for triangulations of non-convex polytopes.
Our experiments were based on linear programming, and it would be interesting to compare
our approach with graph-based algorithms for the stable set problem ([2] for example).

References[1] T. V. Alekseyevskaya. Combinatorial bases in systems of simplices and chambers. Discrete
Mathematics, 157:15-37, 1996.

119

Figure 5: Cardinality of triangulations

Figure 7: Minimum maximum aspect ratio of
triangulations

Figure 6: Weight of triangulations

120[2] E. Balas and J. Xue. Minimum weighted coloring of triangulated graphs, with application
to maximum weighted vertex packing and clique finding in arbitrary graphs. SIAM Journalof Computing, $20(2):209^{-}221$ , 1991.
[3] I. B\’ar\’any and Z. F\"uredi. Empty simplices in euclidian spaces. Canad. Math. Bull., 30:436445, 1987.

[4] M. Bern, P. Chew, D. Eppstein, and J. Ruppert. Dihedral bounds for mesh generation inhigh dimensions. In 6th A CM-SIAM Symp. Discrete Algorithms, pages 189-196, 1995.
[5] M. Bern and D. Eppstein. Mesh generation and optimal triangulation. In Computing in
Euclidean Geometry, volume 4 of Lecture Notes Series on Computing, pages 47-123. World
Scientific, 2nd edition, 1995.
[6] CGAL.

$\mathrm{h}\mathrm{t}\mathrm{t}\mathrm{p}://\mathrm{w}\mathrm{w}\mathrm{w}.\mathrm{c}\mathrm{S}.\mathrm{r}\mathrm{u}\mathrm{u}.\mathrm{n}\mathrm{l}/\mathrm{C}\mathrm{G}\mathrm{A}\mathrm{L}/$.

[7] J. A. De Loera, S. Hosten, F. Santos, and B. Sturmfels. The polytope of all triangulationsof a point configuration. Documenta Mathematica, 1:103-119, 1996.

[8]$\mathrm{D}\mathrm{e}\mathrm{W}\mathrm{a}\mathrm{l}\mathrm{l}$. http:.

$//\mathrm{m}\mathrm{i}\mathrm{l}\mathrm{e}\mathrm{s}.\mathrm{C}\mathrm{n}\mathrm{u}\mathrm{c}\mathrm{e}.\mathrm{c}\mathrm{n}\mathrm{r}.\mathrm{i}\mathrm{t}/\mathrm{c}\mathrm{g}/\mathrm{S}\mathrm{W}\mathrm{O}\mathrm{n}\mathrm{T}\mathrm{h}\mathrm{e}\mathrm{w}_{\mathrm{e}\mathrm{b}}.\mathrm{h}\mathrm{t}\mathrm{m}1$[9] M. R. Garey and D. S. Johnson. Computers and Intractability: A.Guide to the Theory of$NP$ -Completeness. W. H. Freeman and Company, 1979.
[10] H. Imai and K. Imai. Triangulation and convex polytopes. In RIMS Kokyuroku, pages 149166. Research Institute for Mathematical Sciences, Kyoto University, 1996. (in Japanese).

[11] Y. Kyoda, K. Imai, F. Takeuchi, and A. Tajima. A branch-and-cut approach for minimum
weight triangulation. In Proceedings of the 8th Annual International Symposium on Algo), volume 1350 of Lecture Notes in Computer Science,rithms and Computation (ISAA
pages 384-393. Springer Verlag, 1997.
$C’ \mathit{9}7$[12] V. T. Rajan. Optimality of the Delaunay triangulation in
Computational Geometry, pages 357-363, 1991.

$R^{d}$. In Proc. 7th A CM Symp.

[13] J. Ruppert and R. Seidel. On the difficulty of triangulating three-dimensional nonconvex
polyhedra. Discrete.& Computational Geometry, 7:227-253, 1992.
[14] K. Shimada. Physically-based automatic mesh generation. Simulation, $12(1):11-20$ , 1993.
(in Japanese).
[15] P. Shirley and A. Tuchman. A polygonal approximation to direct scalar volume rendering.
Computer Graphics, $24(5):63^{-70}$ , 1990.

[16] Triangle. http:$//\mathrm{w}\mathrm{w}\mathrm{w}.\mathrm{C}\mathrm{S}.\mathrm{C}\mathrm{m}\mathrm{u}.\mathrm{e}\mathrm{d}\mathrm{u}/\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{k}\mathrm{e}/\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}.\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{h}.\mathrm{h}\mathrm{t}\mathrm{m}\mathrm{l}$.