International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Vol XXXV, Part B2. Istanbul 2004 
points can be directly on the circle). The two-dimensional DT is 
illustrated in Figure ! by the dashed lines. The Delaunay 
triangulation is popular for modelling surfaces because among 
all the possible triangulations of a set of points, it creates one 
where the minimum angle in each triangle is maximized 
(triangles are as equilateral as possible), thus being useful for 
interpolation. The generalization to three dimensions of the 
Delaunay triangulation is the Delaunay tetrahedralization: each 
triangle becomes a tetrahedron that satisfies the empty 
circumsphere rule. The DT is unique for a set of points, except 
when there are degenerate cases in the set (if five or more points 
are cospherical in 3D). In these cases, an arbitrary choice must 
be made among all the possible solutions. The number of 
tetrahedra in a DT constructed with # points depends on the 
configuration of these points, and can be up to O(7). 
  
Figure !. Two-dimensional VD (bold lines) and DT (dashed 
lines). 
Most of the properties of the 2D VD/DT generalize to 3D, 
except that the minimum angle in each Delaunay tetrahedra is 
not maximized. There can indeed be almost ‘flat Delaunay 
tetrahedra. These tetrahedra, called slivers, have their four 
vertices almost lying on a plane and thus have a volume of 
nearly zero. For many applications where the Delaunay 
tetrahedralization is used, e.g. to perform simulation in 
engineering or when the tetrahedra are used to perform 
interpolation directly, these tetrahedra are bad and must be 
removed. Here, one might wonder why use them if their 
properties are not good? First, it should be said that in most 
cases the Delaunay tetrahedralization has a tendency to favour 
equilateral tetrahedra over slivers. Second, the Voronoi diagram 
is not affected by them; the Voronoi cells in 3D will still be 
Tround' (ie. relatively spherical) even if the DT has many 
slivers. Third, many GIS operations (e.g. spatial analysis 
functions) use the properties of the VD, and if only one 
tetrahedron is not Delaunay, then the VD is corrupted. 
Both the VD and the DT represent the same thing, just from a 
different viewpoint. The duality between the two structures in 
three dimensions is simple: each polyhedron becomes a point 
and each line becomes a face, and vice-versa. For example, a 
aa 
fs 
  
Figure 2. A Voronoi cell in 3D. The edges are the Delaunay 
edges joining the generator to its natural neighbours. 
704 
Delaunay tetrahedron becomes a Voronoi vertex (its position is 
the centre of the circumsphere around the tetrahedron): a 
Delaunay edge becomes a (convex) Voronoi face; and a 
Delaunay triangular face becomes an edge spanned by the two 
Voronoi vertices that are dual to the two tetrahedra sharing the 
face. For example, in Figure 2, the number of edges joining the 
generator is equal to the number of faces of the Voronoi cell. 
3. 3D VD/DT ALGORITHMS AND DATA 
STRUCTURES 
As mentioned in the previous section, both the VD and the DT 
are geometrically equivalent. By having one structure, its dual 
can always be constructed. Because it is easier, from an 
algorithmic and data structure point-of-view, to manage 
tetrahedra over arbitrary polyhedra (they have a constant 
number of vertices and neighbours), we construct, store and 
modify a VD by working only on its dual. The VD is extracted 
from a DT in O(n) time, n being the number of data points in 
the set. 
We first describe in this section basic operations needed to 
construct and modify a Delaunay tetrahedralization and then 
discuss some possible data structures that can be used to 
efficiently store the DT and/or the VD. 
3.1 Flipping in 3D 
A flip is a local topological operation that modifies the 
configuration of adjacent tetrahedra in a tetrahedralization. If 
we consider five points (a, b, c, d, e! in R^, there exist three 
ways to tetrahedralize them: either with two, three or four 
tetrahedra, depending on their configuration in space. Figure 3 
shows one such configuration: the point e is inside a tetrahedron 
abcd. Figure 4 shows the other configuration where the 
polyhedron abcde is tetrahedralized with either 2 or 3 
tetrahedra. Based on this, we can define different kinds of flips. 
A flipl4 is the operation that will insert e inside a tetrahedron 
abcd (splitting it into 4 tetrahedra), and a flip41 is the inverse 
operation that will delete € and merge together the 4 tetrahedra. 
À flip23 transforms a tetrahedralization of 2 tetrahedra by one 
with 3, and a flip32 is the inverse operation. 
  
  
a 
3 PEN 
| flipl4 / BN 
vd À / | | D 
"à X re / À 
/ T i 
; pdt v S A 
" / > i i M x 
b dés pde 
e pu et De dc uu 
Ha x yet 
d d 
3.2 Point Location 
The point location problem involves finding which tetrahedron 
in a DT contains a query point x. This is needed for different 
operations, for example to insert a new point in a DT or to 
interpolate, as it is explained in Section 5. The method we 
describe here, called "walking, was discussed in the earliest 
papers about the construction of triangulations in two 
dimensions (Gold et al., 1977; Green and Sibson, 1978). 
Interna 
Its gen: 
method 
triangle 
move ic 
the san 
continu 
such ne 
simple 
if a poi 
extra st 
Mücke 
3.3 Ce 
Many € 
One so 
constru 
dimensi 
one dir 
Implem 
dimensi 
1996). ; 
et al. 
paradig; 
ofa D 
deletion 
or simp 
Algoritl 
'Increme 
construi 
added c 
each ins 
steps ar 
identifie 
previous: 
contain 
increme 
develop 
tetrahed 
the hole 
every v 
tetrahed 
algorith 
tetrahed 
algorith 
encount 
Another 
numeric 
is kept 
tetrahed 
of flips. 
contains 
added t