nbul 2004 
position is 
edron); a 
€: and a 
y the two 
haring the 
oining the 
oi cell. 
id the DT 
>, Its dual 
from an 
manage 
constant 
store and 
extracted 
points in 
eeded to 
and then 
used to 
ifies the 
zation. If 
ust three 
or four 
Figure 3 
rahedron 
here the 
2 or 3 
of flips. 
rahedron 
? Inverse 
trahedra. 
a by one 
ahedron 
different 
T or to 
hod we 
earliest 
in two 
  
  
  
Figure 4. Flips 23 and 32. 
Its generalization to three dimensions is straightforward as the 
method uses only the adjacency relationships between the 
triangles. The idea is: starting from a given tetrahedron f, we 
move to one of its neighbours f, if the query point x and /; are on 
the same side of the triangular face shared by / and f,. We 
continue walking from tetrahedron to tetrahedron until /, has no 
such neighbour, which means that /, contains x. The method is 
simple to implement as only one function — one that determines 
if a point is left or right of a plane in 3D — is needed and no 
extra storage is required. It is also very efficient in practice, as 
Mücke et al. (1999) show. 
3.3 Construction Algorithms 
Many different algorithms can be used to construct a 3D VD. 
One solution, as described in Brown (1979), involves firstly 
constructing the convex hull of the set of points in (d + 1) 
dimensions — 4D in our case — and then projecting the result 
one dimension lower to get the Delaunay tetrahedralization. 
Implementations of convex hull algorithms in higher 
dimensions are readily available, e.g. Qhull (Barber et al, 
1996). Another solution is using the DeWall algorithm (Cignoni 
et al, 1998), which is based on the divide-and-conquer 
paradigm. These algorithms might be useful for the construction 
of a DT, but local modifications (insertion of a new point, 
deletion or movement of one) are either slow and complicated, 
or simply impossible. 
Algorithms that allow local modifications are called 
‘incremental insertion algorithms’ and they proceed as follow to 
construct a DT. Starting from a valid DT, each point of a set is 
added one at a time and the tetrahedralization is updated after 
each insertion. To insert a single point x in a DT, the following 
steps are required. First, the tetrahedron that contains x must be 
identified with the point location algorithm described in the 
previous section. Then, all the tetrahedra whose circumspheres 
contain x must be deleted and replaced by new ones. The first 
increment insertion algorithm, valid in any dimensions, was 
developed by Watson (1981). His idea is simple: all the 
tetrahedra that 'conflict' with x are deleted from the DT and then 
the hole thus created is filled by creating edges linking x to 
every vertex of the hole (they prove that the new resulting 
tetrahedra are guaranteed to be Delaunay). Although the 
algorithm is simple to implement, the fact that the 
tetrahedralization is temporarily destroyed can corrupt the 
algorithm. Field (1986) explains the problems that are 
encountered when implementing the method. 
Another incremental insertion algorithm, due to Joe (1991), is 
numerically more stable because a complete tetrahedralization 
is kept during the whole updating process. The conflicting 
tetrahedra are deleted and replaced by new ones by a sequence 
of flips. The first step is the insertion of x in the tetrahedron that 
contains it by using a flip!4. The four new tetrahedra are then 
added to a stack that will control the sequence of flips to 
705 
perform to restore the 'Delaunayness' in the tetrahedralization. 
Each tetrahedron on the stack must be tested against its 
neighbours, if it is not Delaunay then a flip — a flip23 or a 
flip32, depending on the configuration of adjacent tetrahedra — 
will destroy some tetrahedra and replace them by other ones 
(the new ones are then pushed on the stack). The algorithm 
terminates when the stack is empty. The time complexity of this 
algorithm, and of Watson's algorithm, is O(7) for a set of n 
points, not just for the insertion of a single point. This is worst- 
case optimal since a DT of n points can theoretically have up to 
O(n’) tetrahedra. 
3.4 Deletion Algorithms 
The deletion of a single vertex in a Delaunay tetrahedralization 
is often simply referred to as the ‘inverse of the incremental 
insertion algorithm’. Like the insertion operation, it is a local 
operation that involves modifying only some adjacent tetrahedra 
of a DT. Figure 5 illustrates a two-dimensional example where 
the vertex x is deleted from a Delaunay triangulation. Although 
the problem appears to be simple, it is a much more difficult 
task to implement than the insertion of a point, and very few 
algorithms can be found in the literature. 
   
7 S. 7 ^ Ha TL 
Figure 5. Left: x is the vertex to be deleted in a 2D Delaunay 
triangulation. Right: re-triangulation of the polygon. 
The most elegant algorithm, which is valid in any dimensions, is 
by Devillers (2002). In two dimensions, the method involves 
deleting all the triangles incident to the vertex x and re- 
triangulating the hole by using a priority queue of the potential 
triangles that could be used to fill the hole. Devillers' algorithm 
states that the potential triangle having the smallest power — the 
power is a geometric function defined in Aurenhammer (1987) 
— with respect to x is guaranteed to be Delaunay. Because the 
operation is local, the number of edges & incident to a vertex is 
usually used to analyse deletion algorithms. Devillers' method 
has a time complexity of O(k log A) in two dimensions. 
Although possible, the implementation of the algorithm in 3D 
requires many modifications to handle the degenerate cases, 
and, to our knowledge, has not been implemented yet. Because 
the number of tetrahedra present in a Delaunay 
tetrahedralization of m points varies depending on the 
configuration of the points, the time complexity of the method 
in 3D is O(r log &), where ¢ is the number of tetrahedra needed 
to fill the hole. A simpler solution involves keeping a list of all 
potential tetrahedra and testing (Delaunay empty circumsphere 
test) them against each vertex forming the hole. The resulting 
algorithm is slower — a time complexity of O(r &) — and an 
implementation can be found in CGAL (Devillers and Teillaud, 
2003). 
However, these methods temporarily destroy the 
tetrahedralization and some problems can arise. For this reason, 
we have developed a method that uses the flips described in 
Section 3.1 and an algorithm similar to the one implemented in 
CGAL. The methods works fine for most cases and we are 
currently working on making it fully robust for all the 
degenerate cases.