Geometric Data Structures For Computer Graphics: Gabriel Zachmann & Elmar Langetepe

Geometric Data Structures for Computer Graphics ∗
Gabriel Zachmann & Elmar Langetepe

Informatik II/I
University of Bonn, Germany
email: {zach,langetep}@cs.uni-bonn.de
http://web.cs.uni-bonn.de/~zach
http://web.cs.uni-bonn.de/I/staff/langetepe.html
13th August 2003
Course Description
In recent years, methods from computational geometry have been widely adopted by the computer
graphics community yielding elegant and efficient algorithms. This course aims at familiarizing practi-
tioners in the computer graphics field with a wide range of data structures from computational geometry.
It will enable attendees to recognize geometrical problems and select the most suitable data structure
when developing computer graphics algorithms.
The course will focus on algorithms and data structures that have proven to be versatile, efficient,
fundamental and easy to implement. Thus practitioners and researchers will benefit immediately from
this course for their everyday work.
∗ You can find the slides and the multimedia material for this tu-
torial at http://cg.cs.uni-bonn.de/course/
2 Zachmann/Langetepe: Geometric Data Structures for CG
Prerequisites
Participants of this course should be familiar with the basic principles of computer graphics and the type
of problems in the area. Familiarity with computer programming is not required.
The intended audience are practitioners working in 3D computer graphics (VR, CAD/CAM, entertain-
ment, animation, etc.) and students from both computer graphics and computational geometry, possibly
working on a master or PhD thesis.
Syllabus
Our goal is to familiarize practitioners and researchers in computer graphics with some very versatile and
ubiquitous geometric data structures, enable them to readily recognize geometric problems during their
work, modify the algorithms to their needs, and hopefully make them curious about further powerful
treasures to be discovered in the area of computational geometry.
In order to achieve these goals in an engaging yet sound manner, the general concept throughout the
course is to present each geometric data structure in the following way: first, the data strucure will be
defined and described in detail; then, some of its fundamental properties will be highlighted; after that,
one or more computational geometry algorithms based on the data structure will be presented; and finally,
a number of recent, representative and practically relevant algorithms from computer graphics will be
described in detail, showing the utilization of the data structure in a creative and enlightening way.
We have arranged the topics in roughly increasing degree of difficulty. The hierarchical data structures
are ordered by increasing flexibility, while the non-hierarchical topics build on each other. Finally, the last
topic presents a generic technique for dynamizing geometric data structures.
Detailed informations about the topics discussed and the corresponding presenters are given in the
schedule below.
Component 1. Introduction
Introduction – Zachmann/Langetepe.
Component 2. Quadtree/Octree
Construction, complexity, balancing, navigation – Langetepe.
Terrain visualization, iso-surface generation – Zachmann.
Component 3. Bounding volume hierarchies
Definition, construction, hierarchical collision detection – Zachmann.
Component 4. Voronoi diagrams / Delaunay triangulations
Definition, construction, generalizations, applications – Langetepe.
NURBS tesselation, texture synthesis – Zachmann.
Component 5. Distance fields
Definition, modelling, morphing – Zachmann.
Component 6. Dynamization of geometric data structures
Definition, amortized Insert and Delete – Langetepe.
Siggraph 2003 Tutorial 16

Zachmann/Langetepe: Geometric Data Structures for CG 3
Contents
1 Introduction 5
2 Quadtrees and K-d-Trees 5

2.1 Quadtrees and Octrees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.2 Complexity and Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.3 Neighbor Finding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 K-d-Trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.3 Height Field Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Isosurface Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.5 Ray Shooting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5.1 3D Octree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.5.2 5D Octree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3 BSP Trees 15
3.1 Rendering Without a Z-Buffer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Representing Objects with BSPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3 Boolean Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.4 Construction Heuristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.4.1 Convex objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4.2 Cost driven heuristic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4.3 Non-uniform queries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4.4 Deferred, self-organizing BSPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4 Bounding Volume Hierarchies 22

4.1 Construction of BV Hierarchies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.1.1 Bottom-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.1.2 Insertion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.1.3 Top-down . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.1.4 Construction criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.1.5 The criterion for collision detection . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2 Collision Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5 Voronoi Diagrams 28
5.1 Definitions and Elementary Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.1.1 Voronoi Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.1.2 Delaunay Triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.2 Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.3 Generalization of the Voronoi Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.3.1 Voronoi Diagram and Delaunay Triangulation in 3D . . . . . . . . . . . . . . . . . 32
5.3.2 Constrained Voronoi diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.3.3 Other Types of Generalizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.4 Applications of the Voronoi Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.4.1 Nearest Neighbor or Post Office Problem . . . . . . . . . . . . . . . . . . . . . . . 34
5.4.2 Motion planning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5.4.3 Other Applications of the Voronoi Diagram in 2D . . . . . . . . . . . . . . . . . . . 36
5.5 Texture Synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.6 Shape Matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

6 Distance Fields 39
6.1 Computation and representation of DFs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
6.1.1 Propagation method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
6.1.2 Projection of distance functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.2 Applications of DFs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
6.2.1 Morphing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
6.2.2 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
7 Dynamization of Geometric Data Structures 44

7.1 Model of the Dynamization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
7.2 Amortized Insert and Delete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
7.2.1 Amortized Insert: Binary Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
7.2.2 Amortized Delete: Occasional Reconstruction . . . . . . . . . . . . . . . . . . . . . 47
7.2.3 Amortized Insert and Amortized Delete . . . . . . . . . . . . . . . . . . . . . . . . 48
7.3 Worst-Case sensitive Insert and Delete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
7.4 A Simple Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

1 Introduction A natural generalization of the one-dimenional

search tree to k dimensions is shown in Section 2.2.
In recent years, methods from computational ge- The k-d-tree is efficient for axis-parallel rectangu-
ometry have been widely adopted by the computer lar range queries.
graphics community. Many solutions draw their The quadtree description was adapted from de Berg
elegance and efficiency from the mutually enrich- et al.23 and the k-d-tree introduction was taken
ing combination of such geometrical data structures from Klein.55
with computer graphics algorithms.
This course imparts a working knowledge of a 2.1 Quadtrees and Octrees
number of essential geometric data structures and
explains their elegant use in several representative, 2.1.1 Definition
diverse, and current areas of research in computer A quadtree is a rooted tree so that every internal
graphics such as terrain visualization, texture syn- node has four children. Every node in the tree cor-
thesis, modelling, and tesselation. Our selection responds to a square. If a node v has children, their
of data structures and algorithms consists of well- corresponding squares are the four quadrants, as
known concepts, which are both powerful and easy shown in Figure 1.
to implement, ranging from quadtrees to BSP trees
and bounding volume trees. We do not try to pro-
vide an exhaustive survey of the topics touched
upon here — this would be far beyond the scope
NE NW SW SE
of this course. Neither do we try to present leading

edge algorithms because we feel that for practical
purposes a good trade-off between simplicity and
efficiency is important.
Figure 1: An example of a quadtree.
Participants will learn to recognize geometric prob-
lems and develop a thorough understanding of suit- Quadtrees can store many kinds of data. We will
able algorithms. As a prerequisite, they should al- describe the variant that stores a set of points and
ready be familiar with the basic principles of com- suggest a recursive definition. A simple recursive
puter graphics and the type of problems in the area. splitting of squares is continued until there is only
The general concept throughout the course is to one point in a square. Let P be a set of points.
present each geometric data structure in the follow- The definition of a quadtree for a set of points
ing way: first, the data strucure will be defined and in a square Q = [x1 Q : x2Q ] × [y1 Q : y2 Q ] is as
described in detail; then, some of its fundamen- follows:
tal properties will be highlighted; after that, one
or more computational geometry algorithms based • If |P| ≤ 1 then the quadtree is a single leaf
on the data structure will be presented; and finally, where Q and P are stored.
a number of recent, representative and practically
relevant algorithms from computer graphics will be • Otherwise let Q NE , Q NW , QSW and QSE denote
described in detail, showing the utilization of the the four quadrants. Let xmid := (x1 Q + x2 Q )/2
data structure in a creative and enlightening way. and ymid := (y1Q + y2Q )/2, and define
PNE := {p ∈ P : p x > xmid ∧ py > ymid },

2 Quadtrees and K-d-Trees PNW := {p ∈ P : p x ≤ xmid ∧ py > ymid },
PSW := {p ∈ P : p x ≤ xmid ∧ py ≤ ymid },
Within this section we will present some funda- PSE := {p ∈ P : p x > xmid ∧ py ≤ ymid }.
mental geometric data structures. The quadtree consists of a root node v, Q is
In Section 2.1, we introduce the quadtree struc- stored at v. In the following, let Q(v) denote
ture. Its definition and complexity, the recursive
the square stored at v. Furthermore v has four
construction scheme and a standard application are
children: The X-child is the root of the quadtree
presented. Quadtrees and octrees have applications
of the set PX , where X is an element of the set
in mesh generation as shown in Section 2.3, 2.4, 2.5.
{NE, NW, SW, SE}.

2.1.2 Complexity and Construction

The recursive definition implies a recursive con- q
struction algorithm. Only the starting square has
to be chosen adequately. If the split operation can-
not be performed well, the quadtree is unbalanced.
Despite this effect, the depth of the tree is related
to the distance between the points.
Theorem 1 The depth of a quadtree for a set P of

Figure 2: The square q has many west neighbors.
points in the plane is at most log( s/c) + 32 , where c
is the smallest distance between any two points in
P and s is the side length of the initial square. For convenience, we extend the neighbor search.
The cost of the recursive construction and the The given node can also be internal, that is, v and
complexity of the quadtree depends on the depth v0 should be adjacent corresponding to the given
of the tree. direction and should also have the same depth. If
there is no such node, we want to find the deepest
Theorem 2 A quadtree of depth d, which stores a node whose square is adjacent.
set of n points, has O((d + 1)n) nodes and can be The algorithm works as follows. Suppose we
constructed in O((d + 1)n) time. want to find the north neighbor of v. If v happens
to be the SE- or SW-child of its parent, then its
Due to the degree 4 of internal nodes, the total north neighbor is easy to find, it is the NE- or NW-
number of leaves is one plus three times the num- child of its parent, respectively. If v itself is the NE-
ber of internal nodes. Hence it suffices to bound the or NW-child of its parent, then we proceed as fol-
number of internal nodes. lows. Recursively find the north neighbor of µ of
Any internal node v has one or more points in- the parent of v. If µ is an internal node, then the
side Q(v). The squares of the node of a single depth north neighbor of v is a child of µ; if µ is a leaf, the
cover the initial square. Therefore, at every depth north neighbor we seek for is µ itself.
we have at most n internal nodes that gives the This simple procedure runs in time O(d + 1).
node bound.
The most time-consuming task in one step of the Theorem 3 Let T be quadtree of depth d. The
recursive approach is the distribution of the points. neighbor of a given node v in T a given direction,
The amount of time spent is only linear in the as defined above, can be found in O(d + 1) time.
number of points and the O((d + 1)n) time bound
holds. Furthermore, there is also a simple procedure
The 3D equivalent of quadtrees are octrees. The that constructs a balanced quadtree out of a given
quadtree construction can be easily extended to oc- quadtree T. This can be done in time O(d + 1)m
trees in 3D. The internal nodes of octrees have eight and O(m) space if T has m nodes. For details
sons and the sons correspond to boxes instead of see Berg et al.23
squares. Similar results hold for octrees as well.
2.1.3 Neighbor Finding 2.2 K-d-Trees

A simple application of the quadtree of a point set The k-d-tree is a natural generalization of the one-
is neighbor finding, that is, given a node v and a dimenional search tree.
direction, north, east, south or west, find a node v0 Let D be a set of n points in Rk . For conve-
so that Q(v) is adjacent to Q(v0 ). Normally, v is nience let k = 2 and let us assume that all X- and
a leaf and v0 should be a leaf as well. The task is Y-coordinates are different. First, we search for a
equivalent to finding an adjacent square of a given
square in the quadtree subdivision.
Obviously, one square may have many such neigh-
bors, as shown in Figure 2.

Y
i can be computed efficiently. We simply have to
g
compute all nodes v with:
d
k
5
b
R(v) ∩ Q 6= ∅.
c
f
Additionally we have to test whether the points
inside the subtree of v are inside Q.
j The efficiency of the k-d-tree with respect to
e
a
h
range queries depends on the depth of the tree. A
X
balanced k-d-tree can be easily constructed. We
5 10
sort the X- and Y-coordinates. With this order
we recursively split the set into subsets of equal
x=5
size in time O(log n). The construction runs in
x<5 5<x time O(n log n). Altogether the following theorem
holds:
y=4 y=5
y<4 4<y y<5 5<y
Theorem 4 A balanced k-d-tree for n points in
x=2 x<3 x=3 3<x x=8 x=7 the plane can be constructed in O(n log n) and
x<2 2<x x<7
x<8 8<x 7<x
needs O(n) space. A range query with an √ axis-
a b d j g
parallel rectangle can be answered in time O( n +
y=3 y=2 y=6
y<3 3<y y<2 2<y y<6 6<y a), where a denotes the size of the answer.
e c h f k i
Figure 3: A k-d-tree for k = 2 and a rectangular 2.3 Height Field Visualization

range query. The nodes correspond to split lines.
A special area in 3D visualization is the rendering
of large terrains, or, more generally, of height fields.
split-value s of the X-coordinates. Then we split D A height field is usually given as a uniformly-
by the split-line X = s into subsets: gridded square array h : [0, N − 1]2 → R, N ∈ I,
of height values, where N is typically in the order
D<s = {(x, y) ∈ D; x < s} = D ∩ {X < s} of 16,384 or more (see Figure 4). In practice, such
D>s = {(x, y) ∈ D; x > s} = D ∩ {X > s}. a raw height field is often stored in some image file
format, such as GIF. A regular grid is, for instance,
For both sets we proceed with the Y-coordinate one of the standard forms in which the US Geo-
and split-lines Y = t1 and Y = t2 . We repeat logical Survey publishes their data, known as the
the process recursively with the constructed sub- Digital Elevation Model (DEM).32
sets. Thus, we obtain a binary tree, namely the 2-d- Alternatively, height fields can be stored as tri-
tree of the point set D, as shown in Figure 3. Each angular irregular networks (TINs) (see Figure 5).
internal node of the tree corresponds to a split-line. They can adapt much better to the detail and fea-
For every node v of the 2-d-tree we define the rect- tures (or lack thereof) in the height field, so they
angle R(v), which is the intersection of halfplanes can approximate any surface at any desired level of
corresponding to the path from the root to v. For accuracy with fewer polygons than any other rep-
the root r, R(r) is the plane itself; for the sons resentation.64 However, due to their much more
of r, say le f t and right, we produce to halfplanes complex structure, TINs do not lend themselves as
R(le f t) and R(right) and so on. The set of rectan- well as more regular representations to interactive
gles {R(l) : l is a leaf} gives a partition of the plane visualization.
into rectangles. Every R(l) has exactly one point of The problem in terrain visualization is a low view-
D inside. point directed at the horizon: then, there are a few
This structure supports range queries of axis- parts of the terrain that are very close, while the
parallel rectangles. For example, if Q is an axis- majority of the visible terrain is at a larger distance.
parallel rectangle, the set of sites v ∈ D with v ∈ Q Close parts of the terrain should be rendered with
high detail, while distant parts should be rendered

Figure 4: A height field approximated by a grid.15 Figure 5: The same height field approximated by
a TIN.
T-vertices!
4 8
Figure 6: In order to use quadtrees for defining a Figure 7: A quadtree defines a recursive subdivi-
height field mesh, it should be balanced. sion scheme yielding a 4-8 mesh. The dots denote
the newly added vertices. Some vertices have de-
gree 4, some 8 (hence the name).
with very little detail in order to maintain a high depicted in Figure 6. If we render that as-is, then
frame rate. there will be a gap (a.k.a. crack ) between the top
In order to solve this problem, a data structure is left square and the fine detail squares inside the top
needed that allows to quickly determine the desired right square. The vertices causing this problem are
level of detail in each part of the terrain. Quadtrees called T-vertices. Triangulating them would help in
are such a data structure, in particular, since they theory, but in practice this leads to long and thin
seem to be a good compromise between the sim- triangles which have problems on their own. The
plicity of non-hierarchical grids and the good adap- solution is, of course, to triangulate each node.
tivity of TINs. The general idea is to construct Thus, a quadtree offers a recursive subdivision
a quadtree over the grid, and then traverse this scheme to define a triangulated regular grid (see
quadtree top-down in order to render it. At each Figure 7): start with a square subdivided into two
node, we decide whether the detail offered by ren- right-angle triangles; with each recursion step, sub-
dering it is enough, or if we have to go down fur- divide the longest side of all triangles (the hy-
ther. pothenuse) yielding two new right-angle triangles
One problem with quadtrees (and quadrangle- each (hence this scheme is sometimes referred to
based data structures in general) is that nodes are as “longest edge bisection”).65 This yields a mesh
not quite independent of each other. Assume we where all vertices have degree 4 or 8 (except the
have constructed a quadtree over some terrain as

red
red
blue
level 0 1 2 3 4
blue
Figure 8: The 4-8 subdivision can be generated Figure 9: The red quadtree can be stored in the
by two interleaved quadtrees. The solid lines con- unused “ghost” nodes of the blue quadtree.
nect siblings that share a common father.
border vertices), which is why such a mesh is often if error(i) < τ then
called a 4-8 mesh. return
This subdivision scheme induces a directed acyclic end if
graph (DAG) on the set of vertices: vertex j is a if Bi outside viewing frustum then
child of i if it is created by a split of a right angle return
end if
at vertex i. This will be denoted by an edge (i, j).
submesh( j, cl )
Note that almost all vertices are created twice (see V += pi
Figure 7), so all nodes in the graph have 4 children submesh( j, cr )
and 2 parents (except the border vertices).
During rendering, we will choose cells of the sub- where error(i) is some error measure for vertex i,
division at different levels. Let M0 be the fully and Bi is the sphere around vertex i that completely
subdivided mesh (which corresponds to the origi- encloses all descendant triangles.
nal grid) and M be the current, incompletely subdi- Note that this algorithm can produce the same
vided mesh. M corresponds to a subset of the DAG vertex multiple times consecutively; this is easy
of M0 . The condition of being crack-free can be re- to check, of course. In order to produce one strip,
formulated in terms of the DAGs associated with the algorithm has to copy older vertices to the cur-
M0 and M: rent front of the list at places where it makes a
“turn”; again, this is easy to detect, and the inter-
M is crack-free ⇔ ested reader is referred to.65
M does not have any T-vertices ⇔ One can speed up the culling a bit by noticing
that if Bi is completely inside the frustum, then we
∀(i, j) ∈ M : (i0 , j) ∈ M, where parent(j) = {i, i0 }
do not need to test the child vertices any more.
(1)
We still need to think about the way we store our
In other words: you cannot subdivide one trian- terrain subdivision mesh. Eventually, we will want
gle alone, you also have to subdivide the one on to store it as a single linear array for two reasons:
the other side. During rendering, this means that if 1. The tree is complete, so it really would not make
you render a vertex, then you also have to render sense to store it using pointers.
all its ancestors (remember: a vertex has 2 parents).
Rendering such a mesh generates (conceptually) 2. We want to map the file that holds the tree into
a single, long list of vertices that are then fed into memory as-is (for instance, with Unix’ mmap
the graphics pipeline as a single triangle strip. The function), so pointers would not work at all.
pseudo-code for the algorithm looks like this (sim-
plified): We should keep in mind, however, that with cur-
rent architectures, every memory access that can
submesh(i,j)

physical
space ª ⊕
ª ⊕ ⊕ ⊕
node cell
⊕ ª
?
ª ⊕
ª ⊕ ª ⊕
computational
space
⊕ ª
Figure 10: A scalar field is often given in the form Figure 11: Cells straddling the isosurface are tri-
of a curvilinear grid. By doing all calculations in angulated according to a lookup table. In some
computational space, we can usually save a lot of cases, several triangulations are possible, which
computational effort. must be resolved by heuristics.
not be satisfied by the cache is extremely expensive During rendering we need to calculate the indices
(this is even more so with disk accesses, of course). of the child vertices, given the three vertices of a
The simplest way to organize the terrain vertices triangle. It turns out that by cleverly choosing the
is a matrix layout. The disadvantage is that there indices of the top-level vertices this can be done as
is no cache locality at all across the major index. efficiently as with a matrix layout.
In order to improve this, people often introduce The interested reader can find more about this
some kind of blocking, where each block is stored topic in Lindstrom et al.,64 Lindstrom and Pas-
in matrix and all blocks are arranged in matrix or- cucci,65 Balmelli et al.,8 Balmelli et al.,7 and many
der, too. Unfortunately, Lindstrom and Pascucci65 others.
report that this is, at least for terrain visualization,
worse than the simple matrix layout by a factor 10! 2.4 Isosurface Generation
Enter quadtrees. They offer the advantage that
vertices on the same level are stored fairly close One technique (among many others) of visualizing
in memory. The 4-8 subdivision scheme can be a 3-dimenional volume is to extract isosurfaces and
viewed as two quadtrees which are interleaved (see render those as a regular polygonal surface. It can
Figure 8): we start with the first level of the “red” be used to extract the surfaces of bones or organs in
quadtree that contains just the one vertex in the medical scans, such as MRI or CT.
middle of the grid, which is the one that is gen- Assume for the moment that we are given a
erated by the 4-8 subdivision with the first step. scalar field f : R3 → R. Then the task of finding
Next comes the first level of the “blue” quadtree an isosurface would “just” be to find all solutions
that contains 4 vertices, which are the vertices gen- (i.e., all roots) of the equation f (~x ) = t.
erated by the second step of the 4-8 subdivision Since we live in a discrete world (at least in com-
scheme. This process repeats logically. Note that puter graphics), the scalar field is usually given in
the blue quadtree is exactly like the red one, except the form of a curvilinear grid : the vertices of the
it is rotated by 45°. When you overlay the red and cells are called nodes, and we have one scalar and a
the blue quadtree you get exactly the 4-8 mesh. 3D point stored at each node (see Figure 10). Such
Notice that the blue quadtree contains nodes that a curvilinear grid is usually stored as a 3D array,
are outside the terrain grid; we will call these nodes which can be conceived as a regular 3D grid (here,
“ghost nodes”. The nice thing about them is that the cells are often called voxels ).
we can store the red quadtree in place of these ghost The task of finding an isosurface for a given value
nodes (see Figure 9). This reduces the number t in a curvilinear grid amounts to finding all cells of
of unused elements in the final linear array down which at least one node (i.e., corner) has a value
to 33%. less than t and one node has a value greater than

y
x=0 m0 m0
y=0 1
Figure 12: Octrees offer a simple way to compute Figure 13: Volume data layout should match the
isosurfaces efficiently. order of traversal of the octree.
t. Such cells are then triangulated according to a edge will be visited exactly four times during the
lookup table (see Figure 11). So, a simple algorithm complete procedure. Therefore, when we visit an
works as follows:66 compute the sign for all nodes edge for the first time, we compute the vertex of the
(⊕ , > t , ª , < t); then consider each cell in isosurface on that edge, and store the edge together
turn, use the eight signs as an index into the lookup with the vertex in a hash table. So whenever we
table, and triangulate it (if at all). need a vertex on an edge, we first try to look up
Notice that in this algorithm we have only used that edge in the hash table. Our observation also
the 3D array — we have not made use at all of the allows us to keep the size of the hash table fairly
information exactly where in space the nodes are low: when an edge has been visited for the fourth
(except when actually producing the triangles). We time, then we know that it cannot be visited any
have, in fact, made a transition from computational more; therefore, we remove it from the hash table.
space (i.e., the curvilinear grid) to computational
space (i.e., the 3D array). So in the following, we 2.5 Ray Shooting
can, without loss of generality, restrict ourselves to
consider only regular grids, that is, 3D arrays. Ray shooting is an elementary task that frequently
The question is, how can we improve the exhaus- arises in ray tracing, volume visualization, and in
tive algorithm. One problem is that we must not games for collision detection or terrain following.
miss any little part of the isosurface. So we need a The task is, basically, to find the earliest hit of a
data structure that allows us to discard large parts given ray when following that ray through a scene
of the volume where the isosurface is guaranteed to composed of polygons or other objects.
not be. This calls for octrees. A simple idea to avoid checking the ray against
The idea is to construct a complete octree over the all objects is to partition the universe into a regular
cells of the grid96 (for the sake of simplicity, we will grid (see Figure 14). With each cell we store a list
assume that the grid’s size is a power of two). The of objects that occupy that cell (at least partially).
leaves point to the lower left node of their associ- Then, we just walk along the ray from cell to cell,
ated cell (see Figure 12). Each leaf ν stores the min- and check the ray against all those objects that are
imum νmin and the maximum νmax of the 8 nodes stored with that cell.
of the cell. Similarly, each inner node of the octree In this scheme (and others), we need a technique
stores the min/max of its 8 children. called mailboxes that prevents us from checking the
Observe that an isosurface intersects the volume ray twice against the same object:38 every ray gets
associated with a node ν (inner or leaf node) if and a unique ID (we just increment a global variable
only if νmin ≤ t ≤ νmax . This already suggests holding that ID whenever we start with a new ray);
how the algorithm works: start with the root and during traversal, we store the ray’s ID with the ob-
visit recursively all the children where the condi- ject whenever we have performed an intersection
tion holds. At the leaves, construct the triangles as test with it. But before doing an intersection test
usual. with an object, we look into its mailbox whether or
This can be accelerated further by noticing that not the current ray’s ID is already there; if so, then
if the isosurface crosses an edge of a cell, then that

Figure 14: Ray shooting can be implemented ef- Figure 15: The same scenario utilizing an octree.
ficiently with an octree.
we know that we have already performed the inter- those nodes and leaves that are stabbed by the
section test in an earlier cell. ray.
In the following, we will present two methods
Here, we will describe a top-down method.79 The
that both utilize octrees to further reduce the num-
idea is to work only with the ray parameter in order
ber of objects considered.
to decide which children of a node must be visited.
Let the ray be given by
2.5.1 3D Octree
A canonical way to improve any grid-based method ~x = ~p + td~
is to construct an octree (see Figure 15). Here, the and a voxel v by
octree leaves store lists of objects (or, rather, point-
ers to objects). Since we are dealing now with poly- [xl , xh ] × [yl , yh ] × [zl , zh ]
gons and other graphical objects, the leaf rule for
In the following, we will describe the algorithm as-
the octree construction process must be changed
suming that all di > 0; later, we will show that the
slightly:
algorithm works also for all other cases.
1. maximum depth reached; or, First of all, observe that if we already have the
2. only one polygon/object occupies the cell. line parameters of the intersection of the ray with
the borders of a cell, then it is trivial to compute the
We can try to better approximate the geometry of
line intervals half-way in between (see Figure 16):
the scene by changing the rule to stop only when
there are no objects in the cell (or the maximum 1 l
depth is reached). tm
α = (t + tαh ) , α ∈ {x, y, z} (2)
2 α
How do we traverse an octree along a given ray?
Like in the case of a grid, we have to make “hori- So, for 8 children of a cell, we need to compute only
zontal” steps, which actually advance along the ray. three new line parameters. Clearly, the line inter-
With octrees, though, we also need to make “verti- sects a cell if and only if max{til } < min{thj }.
cal” steps, which traverse the octree up or down. The algorithm can be outlined as follows:
All algorithms for ray shooting with octrees can traverse( v, tl , th )
be classified into two classes: compute tm
determine order in which sub-cells are hit by the ray
• Bottom-up: this method starts at that leaf in the
for all sub-cells vi that are hit do
octree that contains the origin of the ray; from traverse( vi , tl |tm , tm |th )
there it tries to find that neighbor cell that is end for
stabbed next by the ray, etc.
where tl |tm means that we construct the lower bound-
• Top-down: this method starts at the root of the ary for the respective cell by passing the appropri-
octree, and tries to recurse down into exactly ate components from tl and tm .

tyh
thx
tm l
tm
y
y > tx
tm
x
tlx
tly
tm l
y < tx
Figure 16: Line parameters are trivial to compute Figure 17: The sub-cell that must be traversed first
for children of a node. can be found by simple comparisons. Here, only
the case tlx > tly is depicted.
In order to determine the order in which sub- y 3 7

cells should be traversed, we first need to determine
which sub-cell is being hit first by the ray. In 2D, 2 6
this is accomplished by two comparisons (see Fig-
ure 17). Then, the comparison of tm m
x with ty tells
5
us which cell is next. 0 4

In 3D, this takes a little bit more work, but is es-
z
sentially the same. First, we determine on which x
side the ray has been entering the current cell by
Table 1. Figure 18: Sub-cells are numbered according to
Next, we determine the first sub-cell to be vis- this scheme.
ited by Table 2 (see Figure 18 for the numbering
scheme). The first column is the entering side de- makes use of octrees to adaptively decompose the
termined in the first step. The third column yields problem.
the index of the first sub-cell to be visited: start The underlying technique is a discretization of
with an index of zero; if one or both of the con- rays, which are 5-dimenional objects. Consider a
ditions of the second column hold, then the corre- cube enclosing the unit sphere of all directions. We
sponding bit in the index as indicated by the third can identify any ray’s direction with a point on
column should be set. that cube, hence it is called direction cube (see Fig-
Finally, we can traverse all sub-cells according to ure 19). The nice thing about it is that we can now
Table 3, where “ex.” means the exit side of the ray perform any hierarchical partitioning scheme that
for the current sub-cell. works in the plane, such as an octree: we just apply
If the ray direction contains a negative compo- the scheme individually on each side.
nent(s), then we just have to mirror all tables along Using the direction cube, we can establish a one-
the respective axis (axes) conceptually. This can be to-one mapping between direction vectors and points
implemented efficiently by an XOR operation. on all 6 sides of the cube, i.e.,
2.5.2 5D Octree S2 ↔ [−1, +1]2 × {+x, −x, +y, −y, +z, −z}
In the previous, simple algorithm, we still walk We will denote the coordinates on the cube’s side
along a ray every time we shoot it into the scene. by u and v.
However, rays are essentially static objects, just like
the geometry of the scene! This is the basic obser-
vation behind the following algorithm.3,6 Again, it

current exit side

sub-cell YZ XZ XY
Side condition index bits
0 4 2 1
tm
z < tlx 0
XY 1 5 3 ex
tm
y < tlx 1
2 6 ex 3
max{til } Side tm
x < ty
l 0 3 7 ex ex
XZ 4 ex 6 5
tz < tly
m 2
tlx YZ 5 ex 7 ex
tly XZ tm l
y < tx 1 6 ex ex 7
YZ
tlz XY tz < tlx
m 2 7 ex ex ex
Table 1: Determines the enter- Table 2: Determines the first Table 3: Determines the traver-
ing side. sub-cell. sal order of the sub-cells.
d~
u u
v + =
v
Figure 19: With the direction cube, we can dis- Figure 20: A uv interval on the direction cube
cretize directions, and organize them with any hi- plus a xyz interval in 3-space yield a beam.
erarchical partitioning scheme.
Within a given universe B = [0, 1]3 (we assume the remaining two uv-intervals define a cone in 3-
it is a box), we can represent all possibly occurring space. Together (more precisely, their Minkowski
rays by points in sum) they define a beam in 3-space that starts at
the cell’s box and extends in the general direction
R = B × [−1, +1]2 of the cone (see Figure 20).
(3)
× {+x, −x, +y, −y, +z, −z} Since we have now defined what a 5D cell of
the octree represents, it is almost trivial to define
which can be implemented conveniently by 6 copies how objects are assigned to sub-cells: we just com-
of 5-dimenional boxes. pare the bounding volume of each object against the
Returning to our goal, we now build six 5-dimen- sub-cells 3D beam. Note that an object can be as-
sional octrees as follows. Associate (conceptually) signed to several sub-cells (just like in regular 3D
all objects with the root. Partition a node in the octrees). The test whether or not an object inter-
octree, if sects a beam could be simplified further by enclos-
ing a beam with a cone, and then checking the ob-
1. there are too many objects associated with it; jects bounding sphere against that cone. This just
and increases the number of false positives a little bit.
2. the node’s cell is too large. Having computed the six 5D octrees for a given
If a node is partitioned, we must also partition its scene, ray tracing through that octree is almost
set of objects and assign each subset to one of the trivial: map the ray onto a 5D point via the direc-
children. tion cube; start with the root of that octree which
Observe that each node in the 5D octree defines a is associated to the side of the direction cube onto
beam in 3-space: the xyz-interval of the first three which the ray was mapped; find the leaf in that oc-
coordinates of the cell define a box in 3-space, and

tree that contains the 5D point (i.e., the ray); check orientation and position of a splitting plane can be
the ray against all objects associated with that leaf. chosen arbitrarily. To get a feeling for a BSP tree,
By locating a leaf in one of the six 5D octrees, Figure 23 shows an example for a set of objects.
we have discarded all objects that do not lie in the The definition of a BSP (short for BSP tree) is
general direction of the ray. But we can optimize fairly straight-forward. Here, we will present a re-
the algorithm even further. cursive definition. Let h denote a plane in Rd , h+
First of all, we sort all objects associated with a and h− denote the positive and negative half-space,
leaf along the dominant axis of the beam by their resp.
minimum (see Figure 21). If the minimum coordi-
nate of an object along the dominant axis is greater Definition 1 (BSP tree)
than the current intersection point, then we can Let S be a set of objects (points, polygons, groups
stop — all other possible intersection points are far- of polygons, or other spatial objects), and let S(ν)
ther away. denote the set of objects associated with a node ν.
Second, we can utilize ray coherence as follows. Then the BSP T(S) is defined by
We maintain a cache for each level in the ray tree
1. If |S| ≤ 1, then T is a leaf ν which stores S(ν) :=
that stores the leaves of the 5D octrees that were
S.
visited last time. When following a new ray, we
first look into the octree leaf in the cache whether 2. If |S| > 1, then the root of T is a node ν; ν stores
it is contained therein, before we start searching for a plane hν and a set S(ν) := {x ∈ S|x ⊆ hν }
it from the root. (this is the set of objects that lie completely
Another trick (that works with other ray acceler- inside hν ; in 3D, these can only be polygons,
ation schemes as well) is to exploit the fact that we edges, or points). ν also has two children T −
do not need to know the first occluder between a and T + ; T − is the BSP for the set of objects
point on a surface and a light source. Any occluder S− := {x ∩ h− +
ν |x ∈ S}, and T is the BSP for
suffices to assert that the point is in shadow. So + +
the set of objects S := {x ∩ hν |x ∈ S}.
we also keep a cache with each light source which
stores that object (or a small set) which has been an This can readily be turned into a general algorithm
occluder last time. for constructing BSPs. Note that a splitting step
Finally, we would like to mention a memory op- (i.e., the construction of an inner node) requires us
timization technique for 5D octrees, because they to split each object into two disjoint fragments if it
can occupy a lot of memory. It is based on the ob- straddles the splitting plane of that node. In some
servation that within a beam defined by a leaf of the applications though (such as ray shooting), this is
octree the objects at the back (almost) never inter- not really necessary; instead, we can just put those
sect with a ray emanating from that cell (see Fig- objects into both subsets.
ure 22). So we store objects with a cell only if they Note that with each node of the BSP a convex
are within a certain distance. Should a ray not hit cell is associated (which is possibly unbounded): the
any object, then we start a new intersection query “cell” associated with the root is the whole space,
with another ray that has the same direction and a which is convex; splitting a convex region into two
starting point just behind that maximum distance. parts yields two convex regions. In Figure 23, the
Obviously, we have to make a trade-off between convex region of one of the leaves has been high-
space and speed here, but when chosen properly, lighted as an example.
the cut-off distance should not reduce performance With BSPs, we have much more freedom to place
too much while still saving a significant amount of the splitting planes than with k-d trees. However,
memory. this also makes that decision much harder (as al-
most always in life). If our input is a set of poly-
gons, then a very common approach is to choose
3 BSP Trees one of the polygons from the input set and use
this as the splitting plane. This is called an auto-
BSP trees (short for binary space partitioning trees ) partition (see Figure 24).
can be viewed as a generalization of k-d trees. like While an auto-partition can have Ω(n2 ) frag-
k-d trees, BSP trees are binary trees, but now the ments, it is possible to show the following in 2D.23,75

1 2 3 4
Figure 21: By sorting objects with in each 5D leaf, Figure 22: By truncating the beam (or rather, the
we can often stop checking ray intersection quite list of objects) we can save a lot of memory usage
early. of a 5D octree, while reducing performance only
insignificantly.
h1 near
polygons
h3 h3 h2
h4 far
h2 h4 polygons
h1
Figure 23: An example BSP tree for a set of ob-

jects.
Figure 25: BSP trees are an efficient data structure

encoding visibility order of a set of polygons.
However, all of these examples producing quadratic

Figure 24: Left: an auto-partition. Right: an ex- BSPs violate the principle of locality : polygons are
ample configuration of which any auto-partition small compared to the extent of the whole set. In
must have quadratic size. practice, no BSPs have been observed that exhibit
the worst-case quadratic behavior.69
Lemma 1
Given a set S of n line segments in the plane, the 3.1 Rendering Without a Z-Buffer
expected number of fragments in an auto-partition
BSP trees were introduced to computer graphics
T(S) is in O(n log n); it can be constructed in time
by Fuchs et al.36 At the time, hidden-surface re-
O(n2 log n).
moval was still a major obstacle towards interactive
In higher dimensions, it is not possible to show computer graphics, because a z-buffer was just too
a similar result. In fact, one can construct sets of costly in terms of memory.
polygons such that any BSP (not just auto-partitions) In this section, we will describe how to solve this
must have Ω(n2 ) many fragments (see Figure 24 problem, not so much because the application itself
for a “bad” example for auto-partitions). is relevant today, but because it nicely exhibits one
of the fundamental “features” of BSP trees: they

3
enable efficient enumeration of all polygons in vis- 1
ibility order from any point in any direction.1 out out
6 in 1 2 3
A simple algorithm to render a set of polygons 6 out
in 5 4
with correct hidden-surface removal, and without out
in
5
out
in out out in out
a z-buffer, is the painter’s algorithm : render the 7 7
out 2 in out
scene from back to front as seen from the current 4
viewpoint. Front polygons will just overwrite the Figure 26: Each leaf cell of BSP representation
contents of the frame buffer, thus effectively hid- of an object is completely inside or completely
ing the polygons in the back. There are polygon outside.
configurations where this kind of sorting is not al-
ways possible, but we will deal with that later.
How can we efficiently obtain such a visibility or- A B
der of all polygons? Using BSP trees, this is al-

most trivial: starting from the root, first traverse
the branch that does not contain the viewpoint,
then render the polygon stored with the node, then
traverse the other branch containing the viewpoint
(see Figure 25).
∪ ∩ \ ª
For sake of completeness, we would like to men-
tion a few strategies to optimize this algorithm. Figure 27: Using BSPs, we can efficiently com-
First of all, we should make use of the viewing pute these boolean operations on solids.
direction by skipping BSP branches that lie com-
pletely behind the viewpoint.
Furthermore, we can perform back-face culling 3.2 Representing Objects with BSPs
as usual (which does not cause any extra costs). We BSPs offer a nice way to represent volumetric polyg-
can also perform view-frustum culling by testing onal objects, which are objects consisting of poly-
all vertices of the frustum against the plane of a gons that are closed, i.e., they have an “inside” and
BSP node. an “outside”. Such a BSP representation of an ob-
Another problem with the simple algorithm is ject is just like an ordinary BSP for the set of poly-
that a pixel is potentially written to many times gons (we can, for instance, build an auto-partition),
(this is exactly the pixel complexity ), although only except that here we stop the construction process
the last write “survives”. To remedy this, we must (see Definition 1) only when the set is empty. These
traverse the BSP from front to back. But in order to leaves represent homogeneous convex cells of the
actually save work, we also need to maintain a 2D space partition, i.e., they are completely “in” our
BSP for the screen that allows us to quickly discard “out”.
those parts of a polygon that fall onto a screen area Figure 26 shows an example for such a BSP rep-
that is already occupied. In that 2D screen BSP, we resentation. In this section, we will follow the con-
mark all cells either “free” or “occupied”. Initially, vention that normals point to the “outside”, and
it consists only of a “free” root node. When a new that the right child of a BSP node lies in the positive
polygon is to be rendered, it is first run through half-space and the left child in the negative half-
the screen BSP, splitting it into smaller and smaller space. So, in a real implementation that adheres to
convex parts until it reaches the leaves. If a part these conventions, we can still stop the construction
reaches a leaf that is already occupied, nothing hap- when only one polygon is left, because we know
pens; if it reaches a free leaf, then it is inserted be- that the left child of such a pseudo-leaf will be “in”
neath that leaf, and this part is drawn on the screen. and the right one will be “out”.
Given such a representation, it is very easy and
efficient, for instance, to determine whether or not
a given a point is inside an object. In the next sec-
tion, we will describe an algorithm for solving a
1 Actually, the first version of Doom used exactly this algorithm slightly more difficult problem.
to achieve its fantastic frame rate (at the time) on PCs even with-
out any graphics accelerator.

H First, we need to define some nomenclature:

T
T − , T + = left and right child of T, resp.

R(T) = region of the cell of node T (which is convex)
T ⊕ , T ª =portion of T on the posi-
tive/negative side of H, resp.
Finally, we would like to define a node T by the

Figure 28: The fundamental step of the construc- tuple (HT , p T , T − , T + ), where H is the splitting
tion is this simple operation, which merges a BSP plane, p is the polygon associated with T (with
and a plane. p ⊂ H).
R(T)
H H The pseudo-code below is organized into 8 cases,
H
HT P HT HT
of which the code details 4 (see Figure 29):
pT pT
split-tree( T, H, P ) → (T ª , T ⊕ )
{P = H ∩ R(T)}
"anti-parallel on" "pos./pos." "mixed" case T is a leaf :
"leaf"
return (T ª , T ⊕ ) ← (T, T)
Figure 29: The main building block of the algo- case “anti-parallel” and “on” :
rithm consists of these four cases (plus analogous return (T ª , T ⊕ ) ← (T + , T − )
ones). case “pos./pos.” :
(T +ª , T +⊕ ) ← split-tree(T + , H)
+ → → T ª ← (HT , p T , T − , T +ª )
T ⊕ ← T +⊕
Figure 30: Computation of boolean operations is case “mixed” :
based on a general merge operation. (T +ª , T +⊕ ) ← split-tree(T + , H, P ∩ R(T + ))
(T −ª , T −⊕ ) ← split-tree(T − , H, P ∩ R(T − ))
T ª ← (HT , p T ∩ H − , T −ª , T +ª )
3.3 Boolean Operations
T ⊕ ← (HT , p T ∩ H + , T −⊕ , T +⊕ )
In solid modeling, a very frequent task is to com-
return (T ª , T ⊕ )
pute the intersection or union of a pair of objects. end case
More generally, given two objects A and B, we want
to compute C := A op B, where op ∈ {∪, ∩, \, ª} This might look a little bit confusing at first sight,
(see Figure 27). This can be computed efficiently but it is really pretty simple. A few notes might be
using the BSP representation of objects.68,69 Fur- in order.
thermore, the algorithm is almost the same for all The polygon P is only needed in order to find
of these operations: only the elementary step that the case applying at each recursion. Computing
processes two leaves of the BSPs is different. P ∩ R(T + ) might seem very expensive. However,
We will present the algorithm for boolean oper- it can be computed quite efficiently by computing
ations bottom-up in three steps. The first step is P ∩ HT+ , which basically amounts to finding the two
a sub-procedure for computing the following sim- edges that intersect with HT . Please see Chin18 for
ple operation: given a BSP T and a plane H, con- more details on how to detect the correct case.
struct a new BSP T̂ whose root is H, such that It seems surprising at first sight that function
T̂ − , T ∩ H − , T̂ + , T ∩ H + (see Figure 28). split-tree does almost no work — it just tra-
This basically splits a BSP tree by a plane and then verses the BSP tree, classifies the case found at each
puts that plane at the root of the two halves. Since recursion, and computes p ∩ H + and p ∩ H − .
we will not need the new tree T̂ explicitly, we will The previous algorithm is already the main build-
describe only the splitting procedure (which is the ing block of the overall boolean operation algo-
bulk of the work anyway). rithm. The next step towards that end is an al-
gorithm that performs a so-called merge operation
on two BSP trees T1 and T2 . Let Ci denote the set
of elementary cells of a BSP, i.e., all regions R(L j )

split merge merge
combine
Figure 31: A graphical depiction of the merge step in the algorithm for boolean operations on objects
represented by BSP trees.
of tree Ti where L j are all the leaves. Then the Operation T1 Result
merge of T1 , T2 yields a new BSP tree T3 such that
C3 = {c1 ∩ c2 |c1 ∈ C1 , c2 ∈ C2 , c1 ∩ c2 6= ∅} (see in T1
∪
Figure 30). out T2
The merge operation consists of two cases. The in T2
∩
first, almost trivial, case occurs when one of the two out T1
operands is a leaf: then at least one of the two re-
in T2c
gions is homogenous, i.e., completely inside or out- \
out T1
side. In the other case, both trees are inhomoge-
nous over the same region of space: then, we just in T2c
ª
split one of the two with the splitting plane from out T2
the root of the other, and we obtain two pairs of
BPSs, that are smaller, and still cover the same re- Furthermore, we would like to point out that the
gions in space; those two pairs can be merged recur- merge function is symmetric: it does not matter
sively (see Figure 31). The following pseudo-code whether we partition T2 with H1 or, the other way
describes this recursive procedure more formally: round, T1 with H2 — the result will be the same.
merge( T1 , T2 ) → T3
if T1 or T2 is a leaf then 3.4 Construction Heuristics
perform the cell-op as required by the boolean op-
eration to be constructed (see below) One can prove that in general, auto-partitions have
else the same complexity as other partitionings.23,75 In
(T2ª , T2⊕ ) ← split-tree(T2 , H1 , . . .) addition, it has been proven that “fat” objects (i.e.,
T3− ← merge(T1− , T2ª ) objects with bounded aspect ration) and “unclut-
tered” scenes admit a BSP with linear size.21,22
T3+ ← merge(T1+ , T2⊕ ) However, for practical applications, the “hidden”
T3 ← (H1 , T3− , T3+ ) constants do matter. So, the construction strategy
end if should produce “good” BSPs. Depending on the ap-
The function cell-op is the only place where the plication itself, however, the definition of exactly
semantic of the general merge operation is special- what a “good” BSP tree is can be quite different.
ized. When we have reached that point, then we Basically, there are two kinds of applications:
know that one of the two cells is homogeneous, so
• Classification: These are applications where the
we can just replace it by the other node’s sub-tree
BSP is used to determine the inside/outside sta-
suitably modified according to the boolean opera-
tus of a point, or whether or not a line intersects
tion. The following table lists the details of this
with an object.
function (assuming that T1 is the leaf):

In this case, we try to optimize the balancedness T has been visited (dito for P(T + )). This probabil-
of the BSP. ity depends, of course, on the application the BSP
is going to be used for. For instance, in the case of
• Visibility: These are applications where the BSP inside/outside queries
is used to sort the polygons in “visibility order”
such as rendering without z-buffer. Vol(T − )
P(T − ) =
Consequently, we try to minimize the number Vol(T)
of splits, i.e., the size of the BSP.
Obviously, trying to optimize Equation 4 globally
is prohibitively expensive. Therefore, a local heuris-
3.4.1 Convex objects
tic must be employed to guide the splitting process.
As an example, consider an convex object. In that A simple heuristic is the following:69 estimate the
caese, an auto-partition has size O(n), takes time cost of the subtrees by the number of polygons, and
O(n2 ) to construct, and is a linear tree. This does add a penalty for polygons that are split by the cur-
not seem to be a very smart BSP (although it is per- rent node, so that
fectly suitable for visibility ordering).
If we allow arbitrary splitting planes, then we C(T) = 1 + |S− |α + |S+ |α + βs (5)
can balance the BSP much better. The construction
time will be where S is the set of polygons associated with a
node, s the set of polygons split by a node, and
T(n) = n + 2T( n2 + αn) ∈ O(n1+δ ) , 0 < α < n α, β are two parameters that can be used to make
2
the BSP more balanced or smaller. (69 reports that
where α is the fraction of polygons that are split α = 0.8, . . . , 0.95 and β = 41 , . . . , 43 are usually
at each node. The following table shows the actual good values to start from.) Again, this is yields an
complexities for some values of α: optimization process, but now it only a local pro-
cess.
α 0.05 0.2 0.4 If the BSP is to be an auto-partition, then a very
quick approximization of the local optimum yields
n1.15 n2 n7
very good results: just sort the polygons according
to their size and evaluate Equation 5 for the first k
As mentioned above, the question now is how to
polygons. Then, choose the one that produces the
choose the splitting planes. We propose the fol-
least costly BSP (subtree), the rationale being that
lowing simple heuristic:2 Compute a representative
the probability that a polygon gets split is propor-
vertex for each polygon (barycenter, bbox center,
tional to its size, so we try to get rid of them as
etc.). Determine a plane, such that on both sides
early as possible.
are about the same number of points, and such that
An even simpler heuristic was proposed by.36
all points are far away from the plane (obviously,
Randomly choose k polygons from S and select the
this is an optimization).
one that produces the least number of splits. They
report that k = 5 yielded near-optimal BSPs (for
3.4.2 Cost driven heuristic visibility ordering).
In order to derive construction criteria, obviously,
one needs to define the quality of a BSP. An abstract 3.4.3 Non-uniform queries
measure is the cost of a BSP T
In the previous section, we assumed that all queries
C(T) = 1 + P(T − )C(T − ) + P(T + )C(T + ) (4) are uniformly distributed over a certain domain.
This is a valid assumption if nothing is known
where P(T − ) is the probability that the left subtree about the distribution. On the other hand, if we
T − will be visited under the condition that the tree know more about it, then we should make use
of this and construct the BSP such that frequent
queries can be answered most quickly.3
2 For the case of convex objects, a heuristic has already been pro-
posed by.87 However, we believe that their heuristic has some
flaws. 3 The same principle is underlying the Huffman encoding scheme.

Indeed, quite often we know more about the queries. 3.4.4 Deferred, self-organizing BSPs
For instance, in ray tracing, the starting points are
Now, what should we do if we just don’t know the
usually not uniformly distributed in space; for in-
query distribution, or if it is too expensive to mea-
stance, they usually do not emanate from the inte-
sure it by experiments? The answer is, of course,
rior of objects. Also, the prominent polygons of an
to defer the complete construction of the BSP, in
object are hit more frequently than those that are
other words, we only build as much of the BSP that
within cavities or completely inside.
is absolutely necessary. In addition, we keep a his-
According to,1 we can estimate the costs of a
tory (in some form) of all the queries we have seen
query by
so far, so whenever we have to build a new part of
C(query) = # nodes visited the BSP we base the construction on that history
≤ depth(BSP) · # stabbed leaf cells (as a best guess of the probability distribution of all
queries to come).2
So, according to this, we should minimize the num- As an example, let us consider the problem of de-
ber of stabbed leaf cells before a polygon hit occurs. tecting intersections between a 3-dim. line segment
The factors influencing the probability that a ray and a set of polygons.4
hits a polygon are Since a BSP is now no longer completely con-
structed before we use it, the nodes must store ad-
1. If the angle between a ray and a polygon is large, ditional information:
then the probability is large.
1. the polygon P defining the splitting plane; or, if
2. If the polygon is large (relative to the total size
it is a preliminary leaf,
of the object/universe), then the probability is
large. 2. a list L ⊆ S of polygons associated with it, for
which the subtree has not been built yet.
3. . . .
Let ω(l) denote the density of all rays l over some The algorithm for answering queries with a ray R
domain D; this could be measured, or it could be now also triggers the construction of the BSP:
derived from the geometry. Let S be the set of poly- testray(R,ν)
gons over which the BSP is to be built. Assign to if ν is a leaf then
each polygon p a score for all P ∈ Lν do
Z if R intersects P then
return hit
score(p) = w(S, p, l)ω(l)dl end if
D
end for
where the weight w is defined as else
ν1 ← child of ν that contains startpoint of R
Area(p) ν2 ← other child
w(S, p, l) = sin2 (n p , rl )
Area(S) testray(R,ν1 )
if no hit in ν1 then
and n p is the normal of p and rl is the direction of testray(R,ν2 )
l. end if
So the algorithm for constructing a BSP adapted end if
to a given distribution is the following random-
Since the line segment is finite (in particular, if it is
ized greedy strategy: sort all polygons according to
short), this algorithm can usually stop much earlier,
score(p); choose randomly one out of the “top k“
but the details have been omitted here.
and split the set S. Thus, polygons with a higher
The initial BPS is just one node (the root) with
probability of getting hit by the ray tend to end up
L = S, i.e., all polygons of the object or scene.
higher in the tree.
Now we need to fill in the open issues, namely:
The BSP thus constructed has been named cus-
tomized BSP by.1 The authors report that it usu- 4 Sometimes, this is also called “collision detection”, in particular
ally has about 2× as many polygons as its “obliv- in the game programming industry, because we are only inter-
ious” version, but the queries have to visit only ested in a yes/no answer, while in ray tracing we want to deter-
2 − −10× less polygons. mine the earliest intersection.

1. when do we split a preliminary leaf? (and which

one?)
2. how do we split?
AABB sphere DOP OBB spherical shell
For the “when” question: we keep an access counter
per node, which gets incremented every time that
node is traversed during a query. Whenever one
of the counters is over a certain threshold, then we
split that node (leaf). The threshold can be an abso- convex hull prism cylinder intersection
of other BVs
lute value or relative to the sum of all counters.
Figure 32: Some of the most commonly used BVs,
For the “how” question: we keep a counter per
and some less often used ones.
polygon P ∈ Lν , which is incremented every time
an intersection with P is found. We sort Lν accord-
ing to this counter; this can be done incrementally met.5 Here, objects can be anything from points
whenever one of the counters changes. If a split is to complete graphical objects. With BV hierarchies,
to be performed for ν, then we use the first polygon almost all queries, which can be implemented with
from Lν as the splitting plane. space partitioning schemes, can also be answered,
It turns out that many polygons are never hit by too. Example queries and operations are ray shoot-
a line segment. Therefore, with this algorithm, a ing, frustum culling, occlusion culling, point loca-
BSP subtree will never be “wasted” on these poly- tion, nearest neighbor, collision detection.
gons, and they will be stored at the end of the lists
Definition 2 (BV hierarchy)
at the leaves.
Let O = {o1 , . . . , on } be a set of elementary objects.
There are other ways to organize the polygon
A bounding volume hierarchy for O, BVH(O), is
lists at the leaves: move the polygon currently be-
defined by
ing hit to the front of the list; or, swap it with the
one in front of it. However, according to,2 the strat- 1. If |O| = e, then BVH(O) := a leaf node that
egy that sorts the list seemed to work best. stores O and a BV of O;
According to,2 the performance gain in their ex-
periments was a factor of about 2–20. 2. If |O| > e, then BVH(O) := a node ν with
n(ν) children ν1 , . . . , νn , where each child νi is a
BV hierarchy
S BVH(Oi ) over a subset Oi ⊂ O,
4 Bounding Volume such that Oi = O. In addition, ν stores a BV
of O.
Hierarchies
The definition mentions two parameters. The
Like the previous hierarchical data structures, bound- threshold e is often set to 1, but depending on the
ing volume hierarchies (BVHs) are mostly used to application, the optimal e can be much larger. Just
prevent performing an operation exhaustively on like sorting, when the set of objects is small, it is
all objects. Like with previously discussed hier- often cheaper to perform the operation on all of
archical data structures, one can improve a huge them, because recursive algorithms always incur
range of applications and queries using BVHs, such some overhead.
as ray shooting, point location queries, nearest- Another parameter in the definition is the arity.
neighbor search, view frustum and occlusion culling, Mostly, BV hierarchies are constructed as binary
geographical data bases, and collision detection (the trees, but again, the optimum can be larger. And
latter will be discussed in more detail below). what is more, as the definition suggests, the out-
Often times, bounding volume (BV) hierarchies degree of nodes in a BV hierarchy does not neces-
are described as the opposite of spatial partitioning sarily have to be constant, although this often sim-
schemes, such as quadtrees or BSP trees: instead of plifies implementations considerably.
partitioning space, the idea is to partition the set Effectively, these two parameters, e and n(ν),
of objects recursively until some leaf criterion is control the balance between linear search/operation
5 However, we will argue at the end that BV hierarchies are just at

the other end of a whole spectrum of hierarchical data structures.

(which is exhaustive), and a maximally recursive d is any metric, very often just the Euclidean dis-
algorithm. tance.) Let
There are more design choices possible according
to theSdefinition. For inner nodes, it only requires diam(G) = max d(g, f )
g, f ∈G
that Oi = O; this means, that the same object
o ∈ O could be associated with several children. be the diameter of G.
Depending on the application, the type of BVs, and Then we can define tightness
the construction process, this may not always be
avoidable. But if possible, you should always split h(B, G)
τ := .
the set of objects into disjoint subsets. diam(G)
Finally, there is, at least, one more design choice:
the type of BV used at each node. Again, this See Figure 33 for an illustration.
does not necessarily mean that each node uses the Since the Hausdorff distance is very sensitive to
same type of BV. Figure 32 shows a number of the outliers, one could also think of other definitions
most commonly used BVs. The difference between such as the following one:
OBBs6 and AABBs is that OBBs can be oriented ar-
bitrarily (hence “oriented bounding boxes”). Definition 4 (Tightness by volume)
DOPs53,57,99 are a generalization of AABBs: basi- Let C(ν) b the set of children of a node ν of the
cally, they are the intersection of k slabs. Prisms BV hierarchy. Let Vol(ν) be the volume of the BV
and cylinders have been proposed by Barequet et stored with ν.
al.9 and Weghorst et al.,94 but they seem to be too Then, we can define the tightness as
expensive computationally. A spherical shell is the
intersection of a shell and a cone (the cone’s apex Vol(ν)
τ := .
coincides with the sphere’s center), and a shell is ∑ν0 ∈C(ν) Vol(ν0 )
the space between two concentric spheres. Finally,
one can always take the intersection of two or more Alternatively, we can define it as
different types of BVs.52
There are three characteristic properties of BVs: Vol(ν)
τ := ,
∑ν0 ∈L(ν) Vol(ν0 )
• tightness,
• memory usage, where L(ν) is the set of leaves beneath ν.
• number of operations needed to test the query
object against a BV. Getting back to the tightness definition based on
the Hausdorff distance, we observe a fundamental
Often, one has to make a trade-off between these difference between AABBs and OBBs:40
properties: generally, the type of BV that offers
• The tightness of AABBs depends on the orien-
better tightness also requires more operations per
tation of the enclosed geometry. What is worse
query and more memory.
is that the tightness of the children of an AABB
Regarding the tightness, one can establish a the-
enclosing a surface of small curvature is almost
oretical advantage of OBBs. But first, we need to
the same as that of the father.
define tightness.40
The worst case is depicted in Figure 34. The
Definition 3 (Tightness by Hausdorff distance) tightness of the father is τ = h/d, while the
Let B be a BV, G some geometry bounded by B, i.e., h0 h/2
tightness of a child is τ 0 = d/2 = d/2 = τ.
g ⊂ B. Let
• The tightness of OBBs does not depend on the
h(B, G) = max min d(b, g) orientation of the enclosed geometry. Instead,
b∈B g∈G
it depends on its curvature, and it decreases ap-
be the directed Hausdorff distance, i.e., the maxi- proximately linearly with the depth in the hier-
mum distance of B to the nearest point in G. (Here, archy.
Figure 35 depicts the situation for a sphere. The
Hausdorff distance from an OBB to an enclosed

h φ Œ/2
h(B, G) d
diam(G) h
d
G
h0
Figure 33: One way to de- Figure 34: The tightness of an Figure 35: The tightness of an
fine tightness is via the directed AABB remains more or less OBB decreases for deeper lev-
Hausdorff distance. constant throughout the levels els in a OBB hierarchy for small
of a AABB hierarchy for surfaces curvature surfaces.
of small curvature.
spherical arc is h = r(1 − cos φ), while the di- simplification is to just approximate each object by
ameter of the arc is d = 2r sin φ. Thus, the its center (baryenter or bounding box center), and
tightness for an OBB bounding a spherical arc then deal only with sets of points during the con-
1−cos φ struction. Of course, when the BVs are finally com-
of degree φ is τ = 2 sin φ , which approaches 0
linearly as φ → 0. puted for the nodes, then the true extents of the
objects must be considered.
This makes OBBs seem much more attractive In the following we will describe algorithms for
than AABBs. The price of the much improved each construction strategy.
tightness is, of course, the higher computational ef-
fort needed for most queries per node when travers- 4.1.1 Bottom-up
ing an OBB tree with a query.
In this class, we will actually describe two algo-
rithms.
4.1 Construction of BV Hierarchies
Let B be the set of BVs on the top-most level of
Essentially, there are 3 strategies to build BV trees: the BV hierarchy that has been constructed so far.80
For each bi ∈ B find the nearest neighbor bi0 ∈ B;
• bottom-up,
let di be the distance between bi and bi0 . Sort B with
• top-down,
respect to di . Then, combine the first k nodes in B
• insertion
under a common father; do the same with the next
From a theoretical point of view, one could pursue k elements from B, etc. This yields a new set B0 ,
a simple top-down strategy, which just splits the and the process is repeated.
set of objects into two equally sized parts, where Note that this strategy does not necessarily pro-
the objects are assigned randomly to either subset. duce BVs with a small “dead space”: in Figure 36,
Asymptotically, this usually yields the same query the strategy would choose to combine the left pair
time as any other strategy. However, in practice, (distance = 0), while choosing the right pair would
the query times offered by such a BV hierarchy are result in much less dead space.
by a large factor worse. The second strategy is less greedy in that it com-
During construction of a BV hierarchy, it is con- putes a tiling for each level. We will describe it
venient to forget about the graphical objects or first in 2D.62 Again, let B be the set of BVs on the
primitives, and instead deal with their BVs and top-most level so far constructed, with |B| = n.
consider those as the atoms. Sometimes, another

Figure 36: A simple greedy strategy can produce Figure 37: A less greedy strategy combines BVs
much “dead space”. by computing a “tiling”.
The algorithm first computes the center ci for each to be performed on the BV tree. See below for a few
bi ∈ B. Then, it sorts B along the x-axis √ with re- criteria.
spect to c x . Now, the set B is split into n/k vertical
i Usually, algorithms in this class have complexity
“slices” (again with respect to cix ). Now, each slice O(n log n).
is sorted according to ciy and subsequently split into
√
n/k “tiles”, so that we end up with k tiles (see 4.1.3 Top-down
Figure 37). Finally, all nodes in a tile are combined This scheme is the most popular one. It seems to
under one common father, its BV is combined, and produce very good hierarchies while still being very
the process repeats with a new set B0 . efficient, and usually it can be implemented easily.
In Rd it works quite
√ similarly: we just split each The general idea is to start with the complete set
d
slice repeatedly by n/k along all coordinate axes. of elementary BVs, split that into k parts, and cre-
ate a BV tree for each part recursively. The splitting
4.1.2 Insertion is guided by some heuristic or criterion that (hope-
fully) produces good hierarchies.
This construction scheme starts with an empty tree.
Let B be the set of elementary BVs. The following
pseudo-code describes the general procedure: 4.1.4 Construction criteria
1: while |B| > 0 do In the literature, there is a vast number of crite-
2: choose next b ∈ B ria for guiding the splitting, insertion, or merging,
3: ν := root during BV tree construction. (Often, the authors
4: while ν 6= leaf do endow the thus constructed BV hierarchy with a
5: choose child ν0 , new name, even though the BVs utilized are well
so that insertion of b into ν0 causes minimal known.) Obviously, the criterion depends on the
increase in the costs of the total tree
application for which the BV tree is to be used. In
6: ν := ν0
7: end while
the following, we will present a few of these crite-
8: end while ria.
For ray tracing, if we can estimate the probabil-
All insertion algorithms only vary step 2 and/or 5. ity that a ray will hit a child box when it has hit the
Step 2 is important because a “bad” choice in the father box, then we know how likely it is, that we
beginning can probably never be made right after- need to visit the child node when we have visited
wards. Step 5 depends on the type of query that is the father node. Let us assume that all rays em-
anate from the same origin (see Figure 38). Then,

θν for α = x, y, z do
θ ν0 sort B along axis α with respect to the BV centers
find
ν0 ½
Area(b1 , . . . , b j )
kα = min j+
j=0...n Area(B)
ν
¾
Area(b j+1 , . . . , bn )
(n − j)
Area(B)
Figure 38: The probability of a ray hitting a child
end for
box can be extimated by the surface area. choose the best kα
where Area(b1 , . . . , b j ) denotes the surface area of
we can observe that the probability that a ray s hits the BV enclosing b1 , . . . , b j .
a child box ν0 under the condition that it has hit the If the query is a point location query (e.g., is a
father box ν is given point inside or outside the object), then the
volume instead of the surface area should be used.
θ ν0 Area(ν0 ) This is because the probability that a point is con-
P(s hits ν0 |s hits ν) = ≈ (6)
θν Area(ν) tained in a child BV, under the condition that it is
contained in the father BV, is proportional to the
where Area denotes the surface area of the BV, and
ratio of the two volumes.
θ denotes the solid angle subtended by the BV. This
For range queries, and for collision detection, the
is because for a convex object, the solid angle sub-
volume seems to be a good probability estimation,
tended by it, when seen from large distances, is ap-
too.
proximately proportional to its surface area.39 So, a
A quite different splitting algorithm does not
simple strategy is to just minimize the surface area
(explicitely) try to estimate any probabilities. It
of the BVs of the children that are produced by a
just approximates each elementary BV/object by
split.6
its center point. It then proceeds as follows. For
A more elaborate criterion tries to establish a cost
a given set B of such points, compute its princi-
function for a split and minimize that. For ray trac-
pal components (the Eigenvectors of the covariance
ing, this cost function can be approximated by
matrix); choose the largest of them (i.e., the one ex-
Area(ν1 ) Area(ν2 ) hibiting the largest variance); place a plane orthog-
C(ν1 , ν2 ) = C(ν1 ) + C(ν2 ) onal to that principal axis and through the barycen-
Area(ν) Area(ν)
(7) ter of all points in B; finally, split B into two sub-
where ν1 , ν2 are the children of ν. The optimal split sets according to the side on which the point lies.
B = B1 ∪ B2 minimizes this cost function: (This description is a slightly modified version of
Gottschalk et al.40 ) Alternatively, one can place the
C(B1 , B2 ) = min C(B0 , B \ B0 ) splitting plane through the median of all points, in-
B0 ∈P (B) stead of the barycenter. This would lead to balanced
trees, but not necessarily better ones.
where B1 , B2 are the subsets of elementary BVs (or
objects) assigned to the children. Here, we have as-
sumed a binary tree, but this can be extended to 4.1.5 The criterion for collision detection
other arities analogously. In the following, we will describe a general cri-
Of course, such a minimization is too expensive terion that can guide the splitting process of top-
in practice, in particular, because of the recursive down BV hierarchy construction algorithms, such
definition of the cost function. So, Fussell and Sub- that the hierarchy produced is good in the sense of
ramanian,37 Müller et al.,67 and Beckmann et al.10 fast collision detection100 (see Section 4.2).
have proposed the following approximation algo- Let C(A, B) be the expected costs of a node pair
rithm: (A, B) under the condition that we have already de-
termined during collision detection that we need to
6 For the insertion scheme, the strategy is to choose that child traverse the hierarchies further down. Assuming
node whose area is increased least.39

B the following that AABBs are used as BVs. How-

ever, similar arguments should hold for all other
kinds of convex BVs.
A The event of box A intersecting box B is equiva-
B1 lent to the condition that B’s “anchor point” is con-
tained in the Minkowski sum A ⊕ B. This situa-
A1 tion is depicted in Figure 39.7 Because B1 is a child
d
of B, we know that the anchor point of B1 must
lie somewhere in the Minkowski sum A ⊕ B ⊕ d,
where d = anchor(B1 ) − anchor(B). Since A1 is
inside A and B1 inside B, we know that A1 ⊕ B1 ⊂
locus of A ⊕ B ⊕ d. So, for arbitrary convex BVs the prob-
possible ability of overlap is
anchor points
Vol(A1 ⊕ B1 )
P(A1 , B1 ) =
Vol(A ⊕ B ⊕ d)
d (11)
Vol(A1 ⊕ B1 )
=
Figure 39: By estimating the volume of the Vol(A ⊕ B)
Minkowski sum of two BVs, we can derive an es-
timate for the cost of the split of a set of polygons In the case of AABBs, it is safe to assume that
associated with a node. the aspect ratio of all BVs is bounded by α. Conse-
quently, we can bound the volume of the Minkowski
sum by
binary trees and unit costs for an overlap test, this
can be expressed by 2p
Vol(A) + Vol(B) + Vol(A) Vol(B) ≤
C(A, B) = 4 + ∑ P(Ai , Bj ) · C(Ai , Bj ) (8) α
i,j=1,2 Vol(A ⊕ B) ≤
p
where Ai , Bj are the children of A and B, resp., and Vol(A) + Vol(B) + 2α Vol(A) Vol(B) (12)
P(Ai , Bj ) is the probability that this pair must be
So we can estimate the volume of the Minkowski
visited (under the condition that the pair (A, B) has
sum of two boxes by
been visited).
An optimal construction algorithm would need Vol(A ⊕ B) ≈ 2(Vol(A) + Vol(B))
to expand (8) down to the leaves:
yielding
C(A, B) =P(A1 , B1 ) + P(A1 , B1 )P(A11 , B11 )
Vol(A1 ) + Vol(B1 )
+ P(A1 , B1 )P(A12 , B11 ) + . . . + P(A1 , B1 ) ≈ (13)
Vol(A) + Vol(B)
P(A1 , B2 ) + P(A1 , B2 )P(A11 , B21 )
+... Since Vol(A) + Vol(B) has already been com-
(9) mitted by an earlier step in the recursive construc-
tion, Equation 10 can be minimized only by min-
and then find the minimum. Since we are inter- imizing Vol(A1 ) + Vol(B1 ). This is our criterion
ested in finding a local criterion, we approximate for constructing restricted boxtrees.
the cost function by discarding the terms corre- Construction algorithm. According to the cri-
sponding to lower levels in the hierarchy, which terion derived above, each recursion step will try
gives to split the set of polygons so that the cost func-
¡ tion (10) is minimized. This is done by trying to
C(A, B) ≈ 4 1 + P(A1 , B1 ) + . . .
¢ (10) find a good splitting for each of the three coordinate
+ P(A2 , B2 )
7 In the figure, we have chosen the lower left corner of B as
Now we will derive an estimate of the probability its anchor point, but this is arbitrary, of course, because the
P(A1 , B1 ). For sake of simplicity, we will assume in Minkowski sum is invariant under translation.

axes, and then selecting the best one. Along each sphere,46,74 OBBs,40 DOPs,57,99 AABBs,58,88,98 and
axis, we consider three cases: both subsets form convex hulls,31 to name but a few.
lower boxes with respect to its parent, both are up- Given two hierarchical BV volume data struc-
per boxes, or one upper and one lower box. tures for two objects A and B, almost all hierar-
In each case, we first try to find a good “seed” chical collision detection algorithms implement the
polygon for each of the two subsets, which is as following general algorithm scheme:
close as possible to the outer border that is perpen- traverse(A,B)
dicular to the splitting axis. Then, in a second pass, if A and B do not overlap then
we consider each polygon in turn, and assign it to return
that subset whose volume is increased least. In or- end if
der to prevent “creeping greediness”,8 we run alter- if A and B are leaves then
natingly through the array of polygons. After good return intersection of primitives
splitting candidates have been obtained for all three enclosed by A and B
axes, we just pick the one with least total volume of else
the subsets. for all children A[i] and B[j] do
traverse(A[i],B[j])
The algorithm and criterion we propose here could
end for
also be applied to construct hierarchies utilizing end if
other kinds of BVs, such as OBBs, DOPs, and even
convex hulls. We suspect that the volume of AABBs This algorithm quickly “zooms in” on pairs of close
would work fairly well as an estimate of the volume polygons. The characteristics of different hierarchi-
of the respective BVs. cal collision detection algorithms lie in the type of
This algorithm has proven to be geometrically BV used, the overlap test for a pair of nodes, and
robust, since there is no error propagation. There- the algorithm for construction of the BV trees.
fore, a simple epsilon guard for all comparisons suf- The algorithm outlined above is essentially a si-
fices. multaneous traversal of two hierarchies, which in-
It can be shown that this algorithm is, under cer- duces a so-called recursion tree (see Figure 40).
tain assumptions,9 in O(n), where n is the number Each node in this tree denotes a BV overlap test.
of polygons.100 Leaves in the recursion tree denote an intersection
test of the enclosed primitives (polygons); whether
or not a BV test is done at the leaves depends on
4.2 Collision Detection how expensive it is, compared to the intersection
Fast and exact collision detection of polygonal ob- test of primitives.
jects undergoing rigid motions is at the core of During collision detection, the simultaneous traver-
many simulation algorithms in computer graphics. sal will stop at some nodes in the recursion tree. Let
In particular, all kinds of highly interactive appli- us call the set of nodes, of which some children are
cations such as virtual prototyping need exact col- not visited (because their BVs do not overlap), the
lision detection at interactive speed for very com- “bottom slice” through the recursion tree (see the
plex, arbitrary “polygon soups”. It is a fundamen- dashed lines in Figure 40).
tal problem of dynamic simulation of rigid bodies, One idea is to save this set for a given pair of ob-
simulation of natural interaction with objects, and jects.63 When this pair is to be checked next time,
haptic rendering. we can start from this set, going either up or down.
Bounding volume trees seem to be a very effi- Hopefully, if the objects have moved only a little
cient data structure to tackle the problem of colli- relative to each other, the number of nodes that
sion detection for rigid bodies. All kinds of differ- need to be added or removed from the bottom slice
ent types of BVs have been explored in the past: is small. This scheme is called incremental hierar-
chical collision detection.
8 This happens, for instance, when the sequence of polygons hap-
pens to be ordered such that each polygon increases one subset’s 5 Voronoi Diagrams
volume just a little bit, so that the other subset never gets a
polygon assigned to.
9 These assumption have been valid in all practical cases we have For a given set of sites inside an area the Voronoi
encountered so far. diagram is a partition of the area into regions of

A A1 1
B C B2 B3 C2 C3 2 3
D E F G D4 D5 E4 E5 F6 F7 G6 G7 4 5 6 7
D6 D7 E6 E7 F4 F5 G4 G5
Figure 40: The recursion tree is induced by the simultaneous traversal of two BV trees.
B a simple incremental construction approach which

can be easily generalized to 3D (see Section 5.3.1).
A We present generalizations of the Voronoi dia-
B1
gram and the Delaunay Triangulation in Section 5.3.
In Section 5.3.1 transformations to three dimen-
sions are given and in Section 5.3.2 the concept of
constrained Voronoi diagrams is introduced. A col-
A1
lection of other interesting generalizations is pre-
sented in Section 5.3.3.
A2
B2 In Section 5.4 applications of the Voronoi dia-
gram and the Delaunay triangulation in 2D and 3D
are shown. First, in Section 5.4.1 we discuss the
Figure 41: Hierarchical collision detection can famous post office problem. Then a simple appli-
discard many pairs of polygons with one BV cation of the Voronoi diagram of line segments for
check. Here, all pairs of polygons in A1 and B2 motion planning is presented in Section 5.4.2. Fi-
can be discarded. nally, in Section 5.4.3 a collection of 2D-applications
is shown.
Note, that we can only sketch many of the sub-
the same neighborship. The Voronoi diagram and jects here. For further details and further literature
its dual have been used for solving numerous prob- see one of the surveys mentioned above. The fig-
lems in many fields of science. ures are taken from Aurenhammer and Klein.5
We will concentrate on its application to geo-
metric problems in 2D and 3D. For an overview
of the Voronoi diagram and its dual in computa- 5.1 Definitions and Elementary
tional geometry one may consult the surveys by Properties
Aurenhammer,4 Bernal,11 Fortune34 and Auren-
5.1.1 Voronoi Diagram
hammer and Klein.5 Additionally, chapters 5 and 6
of Preparata and Shamos77 and chapter 13 of Edels- Let S a set of n ≥ 3 point sites p, q, r, . . . in the
brunner29 could be consulted. plane. In the following we assume that the points
We start in Section 5.1 with the simple case of are in general position, i.e., no four of them lie on
the Voronoi diagram and the Delaunay triangula- the same circle and no three of them on the same
tion of n points in the plane, under the Euclidean line.
distance. Additionally we mention some of the ele- For points p = (p1 , p2 ) and x = (x1 , x2 ) let
mentary structural properties that follow from the d(p, x) denote their Euclidean distance. By pq we
definitions. denote the line segment from p to q. The closure of
In Section 5.2 different algorithmic schemes for a set A will be denoted by A.
computing the structures are mentioned. We present

Definition 5 For p, q ∈ S let vertices; they belong to the common boundary of

three or more Voronoi regions.
B(p, q) = {x | d(p, x) = d(q, x)} There is an intuitive way of looking at the Voronoi
diagram. For any point x in the plane we can ex-
be the bisector of p and q. B(p, q) is the perpendic- pand the circle C(r) with center x and radius r by
ular line through the center of the line segment pq. increasing r continuously. We detect three cases de-
It separates the halfplane pending on which event occurs first:
D(p, q) = {x | d(p, x) < d(q, x)} • If C(r) hits one of the n sites, say p, then x ∈
VR(p, S).
containing p from the halfplane D(q, p) contain-
ing q. We call • If C(r) hits two sites p and q simultaneously x
\ belongs to the Voronoi edge of p and q.
VR(p, S) = D(p, q)
q∈S,q6= p
• If C(r) hits three sites p, q and r simultaneously
x is the Voronoi vertex of p, q and r.
the Voronoi region of p with respect to S. Finally,
We will enumerate some of the significant prop-
the Voronoi diagram of S is defined by
erties of Voronoi diagrams.
[
V(S) = VR(p, S) ∩ VR(q, S). 1. Each Voronoi region VR(p, S) is the intersec-
p,q∈S,p6=q tion of at most n − 1 open halfplanes containing
the site p. Every VR(p, S) is open and convex.
An illustration is given in Figure 42. It shows
Different Voronoi regions are disjoint.
how the plane is decomposed by V(S) into Voronoi
regions. Note that it is convenient to imagine a 2. A point p of S lies on the convex hull of S iff its
simple closed curve Γ around the “interesting” part Voronoi region VR(p, S) is unbounded.
of the Voronoi diagram.
3. The Voronoi diagram V(S) has O(n) many edges
and vertices. The average number of edges in
the boundary of a Voronoi region is less than 6.
The Voronoi diagram is a simple linear structure

Γ
and provides for a partition of the plane into cells
of the same neighborship. We omit the proofs and
refer to the surveys mentioned in the beginning.
Note, that the Voronoi edges and vertices build
a graph. Therefore the diagram normally is repre-
sented by a graph of linear size. For example the
diagram can be represented by a doubly connected
edge list DCEL, see de Berg et al.,23 or with the help
of an adjacency matrix.
5.1.2 Delaunay Triangulation

We consider the dual graph of the Voronoi diagram,
Figure 42: A Voronoi diagram of points in the Eu- the so called Delaunay triangulation. In general, a
clidean plane. triangulation of S is a planar graph with vertex set
S and straight line edges, which is maximal in the
The common boundary of two Voronoi regions sense that no further straight line edge can be added
belongs to V(S) and is called a Voronoi edge, if it without crossing other edges. The triangulation of
contains more than one point. If the Voronoi edge a point set S has not more than O(|S|) triangles.
e borders the regions of p and q then e ⊂ B(p, q)
holds. Endpoints of Voronoi edges are called Voronoi

Definition 6 The Delaunay triangulation DT(S) is • Sweep

the dual Graph of the Voronoi diagram. The edges
are convenient for the construction of the Voronoi
of DT(S) are called Delaunay edges.
diagram or the Delaunay triangulation, respectively.
Obviously, the Delaunay triangulation DT(S) is They can also be generalized to other metrics and
a triangulation of S, an example is shown in Fig- sites other than points, for example line segments
ure 43. or polygonal chains. The result of the algorithms is
We present two equivalent definitions of the De- stored in a graph of linear size, see above.
launay triangulation. They are applied for the com- All these approaches run in deterministic time
putation of the diagram and give also rise to gener- O(n log n). We explain a simple Incremental con-
alization, for example if the dual of a Voronoi dia- struction technique which runs in O(n log n) ex-
gram is no longer well-defined. pected time and computes the Delaunay triangula-
tion. The presentation is adapted from Klein and
1. Two points p, q of S give rise to a Delaunay edge Aurenhammer.5 The technique can easily be gen-
iff a circle C exists that passes through p and q eralized to the three dimensional case as we will see
and does not contain any other site of S in its in Section 5.3.1.
interior or boundary. Simple incremental construction: The inser-
tion process is described as follows: We construct
2. Three points of S give rise to a Delaunay trian- DT = DT({p , . . . , p , p }) by inserting the site
i 1 i−1 i
gle iff their circumcircle does not contain a point p into DT . We follow Guibas and Stolfi41 and
i i−1
of S in its interior. construct DTi by exchanging edges, using Law-
son’s60 original edge flipping procedure, until all
edges invalidated by pi have been removed.
w
V(S)
It is helpful to extend the notion of triangle to
the unbounded face of the Delaunay triangulation.
DT(S) s
If pq is an edge of the convex hull of S we call the
supporting outer halfplane H not containing S an
r
p v infinite triangle with edge pq. Its circumcircle is H
itself, the limit of all circles through p and q whose
center tend to infinity within H. As a consequence,
each edge of a Delaunay triangulation is now adja-
q
cent to two triangles.
Figure 43: Voronoi diagram and Delaunay trian- Those triangles of DTi−1 whose circumcircles con-
gulation. tain the new site, pi , are said to be in conflict with
pi . According to the (equivalent) definition of the
DTi , they will no longer be Delaunay triangles.
5.2 Computation Let qr be an edge of DTi−1 , and let T(q, r, t) be
the triangle adjacent to qr that lies on the other
The construction of the Voronoi diagram has time side of qr than pi ; see Figure 44. If its circumcir-
complexity Θ(n log n). The lower bound Ω(n log n) cle C(q, r, t) contains pi then each circle through
can be achieved by the following reductions. q, r contains at least one of pi , t. Consequently, qr
cannot belong to DTi , due to the (equivalent) def-
• A reduction to the convex hull problem is given
inition. Instead, pi t will be a new Delaunay edge,
by Shamos.85
because there exists a circle contained in C(q, r, t)
• A reduction to the e-closeness problem is given that contains only pi and t in its interior or bound-
by Djidjev and Lingas26 and by Zhu and Mirza- ary. This process of replacing edge qr by pi t is called
ian.101 an edge flip.
The necessary edge flips can be carried out effi-
The well-known computation paradigms ciently if we know the triangle T(q, s, r) of DTi−1
that contains pi , see fig. Figure 45. The line seg-
• Incremental construction,
ments connecting pi to q, r, and s will be new De-
• Divide-and-Conquer and
launay edges, by the same argument from above.

t
5.3 Generalization of the Voronoi
Diagram
q
5.3.1 Voronoi Diagram and Delaunay
r Triangulation in 3D
C(pi,t) We will see that incremental construction is also

pi
appropriate for the 3D case. The following descrip-
tion was adapted from Aurenhammer and Klein.5
C(q,r,t) Let S be a set of n point sites in 3D. The bisec-
tor of two sites p, q ∈ S is the perpendicular plane
Figure 44: If triangle T(q, r, t) is in conflict with pi
through the midpoint of the line segment pq. The
then former Delaunay edge qr must be replaced
region VR(p, S) of a site p ∈ S is the intersec-
by pi t.
tion of halfspaces bounded by bisectors, and thus
is a 3-dimenional convex polyhedron. The bound-
Next, we check if e. g. edge qr must be flipped. If ary of VR(p, S) consists of facets (maximal subsets
so, the edges qt and tr are tested, and so on. We within the same bisector), of edges (maximal line
continue until no further edge currently forming segments in the boundary of facets), and of vertices
a triangle with, but not containing pi , needs to be (endpoints of edges). The regions, facets, edges, and
flipped, and obtain DTi . vertices of V(S) define a cell complex in 3D.
Two task have to be considered: This cell complex is face-to-face: if two regions
have a non-empty intersection f , then f is a face
1. Find the triangle of DTi−1 that is in conflict
(facet, edge, or vertex) of both regions. As an ap-
with pi .
propriate data structure for storing a 3-dimenional
2. Perform all flips starting from this triangle. cell complex we mention the facet-edge structure in
It can be shown that the second task is bounded Dobkin and Laszlo.27
by the degree of pi in the new triangulation. If the Complexity: The number of facets of VR(p, S)
triangle of DTi−1 containing pi is known, the struc- is at most n − 1, at most one for each site q ∈
tural work needed for computing DTi from DTi−1 S \ {p}. Hence, by the Eulerian polyhedron for-
is proportional to the degree d of pi in DTi . mula, the number of edges and vertices of VR(p, S)
So we yield an obvious O(n2 ) time algorithm is O(n), too. This shows that the total number of
for constructing the Delaunay triangulation of n components of the diagram V(S) in 3D is O(n2 ).
points: we can determine the triangle of DTi−1 In fact, there are configurations S that force each
containing pi within linear time, by inspecting all pair of regions of V(S) to share a facet, thus¡ achiev-
¢
candidates. Moreover, the degree of pi is trivially ing their maximum possible number of n2 ; see,
bounded by n. e.g., Dewdney and Vranch.25 This fact sometimes
The last argument is too crude. There can be sin- makes Voronoi diagrams in 3D less useful com-
gle vertices in DTi that do have a high degree, but pared to 2-space. On the other hand, Dwyer28
their average degree is bounded by 6. showed that the expected size of V(S) in d-space is
With a special implementation using a directed only O(n), provided S is drawn uniformly at ran-
acyclic graph (DAG), also called Delaunay tree due dom in the unit ball. This result indicates that high-
to Boissonnat and Teillaud,12 we can detect the tri- dimenional Voronoi diagrams will be small in many
angles of DTi−1 which are in conflict with pi in practical situations.
O(log i) expected time. In analogy to the 2-dimenional case, the Delau-
Altogether we get the following result: nay triangulation DT(S) in 3D is defined as the ge-
ometric dual of V(S). It contains a tetrahedron for
Theorem 7 The Delaunay triangulation of a set each vertex, a triangle for each edge, and an edge
of n points in the plane can be easily incremen- for each facet, of V(S). Equivalently, DT(S) may
tally constructed incrementally in expected time be defined using the empty sphere property, by in-
O(n log n), using expected linear space. The av- cluding a tetrahedron spanned by S as Delaunay iff
erage is taken over the different orders of inserting its circumsphere is empty of sites in S. The circum-
the n sites. centers of these empty spheres are just the vertices

q q
q
t t
t
pi pi pi
F T
s s
s r SF r
SF r
(i) (ii) (iii)
Figure 45: Updating DTi−1 after inserting the new site pi . In (ii) the new
Delaunay edges connecting pi to q, r, s have been added, and
edge qr has already been flipped. Two more flips are neces-
sary before the final state shown in (iii) is reached.
of V(S). DT(S) is a partition of the convex hull

of S into tetrahedra, provided S is in general posi-
tion. Note that the edges of DT(S) may form the
complete graph on S.
Simple incremental construction: Among the
various proposed methods for constructing V(S)
in 3D, incremental insertion of sites (compare Sec-
p
q p
q
tion 5.2) is most intuitive and easy to implement.
Basically, two different techniques for integrating a
new site p into V(S) have been applied. The more
obvious method first determines all facets of the re-
gion of p in the new diagram, V(S ∪ {p}), and then
deletes the parts of V(S) interior to this region; see
Figure 46: The Voronoi diagram of a point set in e.g. Watson,93 Field,33 and Tanemura et al.86 In-
L1 and L2 . Note, that there are structural differ- agaki et al.47 describe a robust implementation of
ences. this method.
In the dual environment, this amounts to detect-
ing and removing all tetrahedra of DT(S) whose
circumspheres contain p, and then filling the ’hole’
with empty-sphere tetrahedra with p as apex, to
obtain DT(S ∪ {p}). An example of an edge flip in
3D is shown in Figure 48. Joe,49 Rajan,78 and Edels-
brunner and Shah30 follow a different and numeri-
cally more stable approach. Like in the planar case,
after having added a site to the current Delaunay
triangulation, certain flips changing the local tetra-
hedral structure are performed in order to achieve
local “Delaunayhood”. The existence of such a se-
quence of flips is less trivial, however. Joe48 demon-
strated that no flipping sequence might exist that
turns an arbitrary tetrahedral triangulation for S
into DT(S).
Figure 47: An Euclidean Voronoi diagram of line A complete algorithm with run time O(n2 ) can
segments. be found in Shah.30

Figure 48: Two–into–three tetrahedra flip for five sites.
5.3.2 Constrained Voronoi diagrams nay triangulation DT(S, L) have been proposed in
Lee and Lin,61 Chew,16 Wang and Schubert92 and
For many applications it is not useful to consider
Wang.91 In Seidel81 and Kao and Mount51 one
neighborship relations without any constraints. For
will find good implementation schemes. An appli-
example, a small distance between two cities may
cation of DT(S, L) to quality mesh generation can
become unimportant if there is a (natural or artifi-
be found in Chew.17
cial) border line between them and the border line
itself should not serve as a new site. To overcome
such problems the concept of constrained Voronoi 5.3.3 Other Types of Generalizations
diagrams was introduced.61 We simply list some of the generalization schemes
Let S be a set of n points in the plane, and let and show examples of some intuitive ones.
L be a set of non-crossing line segments spanned
by S, thus we have |L| ≤ 3n − 6. The segments • Different metrics
in L may be considered as obstacles. The bounded – L1 , a comparison of L1 and L2 is shown in
distance between two points x and y in the plane is Figure 46
defined as follows: – L∞
½
d(x, y) if xy ∩ L = ∅ – Convex distance functions
b(x, y) =
∞ otherwise • Different space
where d stands for the Euclidean distance. In the – On trees and graphs
resulting constrained Voronoi diagram VD(S, L), – Higher dimensions
regions of sites that are close but not visible from
• Weights
each other are separated by segments in L, an ex-
• More general sites
ample is shown in Figure 49.
The exact dual of VD(S, L) may be no longer a – Line segments (see Figure 47).
full triangulation of S (even if S is included) but – Polygonal chains
we can modify VD(S, L) in order to dualize it into • Farthest point Voronoi diagram
a "near to Delaunay" triangulation DT(S, L) that • K-th order Voronoi digram
also includes L. • Colored objects
For every line segment l we proceed as follows:
All sites of the clipped regions to the left of a seg-
5.4 Applications of the Voronoi
ment l build a special Voronoi diagram to the right
and vice versa. The neighborship within these spe- Diagram
cial diagrams lead to additional edges inserted into 5.4.1 Nearest Neighbor or Post Office Problem
the dual of VD(S, L). Note that the endpoints of
the line segments are always neighbors in the new We consider the well-known post office problem.
diagrams and therefore the line segments itself are For a set S of sites in the plane and an arbitrary
inserted (see Figure 49). query point x we want to compute the point of S
Algorithms for computing the constrained Voro- closest to x efficiently.
noi diagram VD(S, L) and the constrained Delau-

Figure 49: The constrained Voronoi diagram VD(S, L) and its dual.
In the field of computational geometry there is a Theorem 8 Given a set S of n point sites in the
general technique for solving such query problems. plane, one can, within O(n2 ) time and storage,
One tries to decompose the query set into classes construct a data structure that supports nearest
so that every class has the same answer. Now for a neighbor queries: for an arbitrary query point x,
single answer we only have to determine its class. its nearest neighbor in S can be found in time
This technique is called locus approach. O(log n).
The Voronoi diagram represents the locus ap-
proach for the post office problem. The classes cor- The simple technique can be easily extended to
respond to the regions of the sites. For a query 3D. There are also more efficient approaches, i.e.,
point x we want to determine its class/region and Edelsbrunner29 constructs a O(log n) search struc-
return its owner. ture for the Voronoi diagram in linear time and
To solve this task a simple technique can be ap- with linear space.
plied. We draw a horizontal line through every ver-
tex of the diagram and sort the lines in O(n log n) 5.4.2 Motion planning
time, see Figure 50. The lines decompose the di- We present a simple application of Voronoi dia-
agram into slabs. For every slab we sort the set grams for computing collision-free paths in polyg-
of crossing edges of the Voronoi diagram in linear onal environments. Let us assume that a two di-
time. Altogether we need O(n2 ) time for the sim- mensional environment with polygonal obstacles
ple construction. is given. For simplicity we assume that a circle-
shaped robot is given and a collision-free path from
s to t should be computed (see Figure 51).
We can solve this problem by computing the
Voronoi diagram VD(S) of the corresponding line
segments. This gives rise to a reasonable solu-
tion. Let us assume that the robot moves between
x
to obstacle line segments l1 and l2 , for safetyness

it should maximize its distance to both segments.
That is, for every position x we want to maximize
|xli | = min |xy|, i = 1, 2 .

y∈li
Figure 50: After constructing the slabs, a query
point x can be located quickly. This goal is reached if the robot moves along the
biscetor of l1 and l2 . Altogether the robot should
For a query point x we locate its slab in O(log n) try to move along the Voronoi diagram of the line
time and afterwards its region in O(log n) time by segments of the polygonal scene. For every seg-
binary search. ment of VD(S) we can easily compute whether the
robot fits between the corresponding obstacles.

t' t
yt
s' s ys
l
Figure 51: Computing a collision-free path from s to t for a circle-shaped robot in

the presence of polygonal obstacles.
Now we can solve our problem as follows: • Nearest Neighbor Search

• Compute the Voronoi diagram VD(S) of the set – O(n) for all nearest neighbors of the sites
of line segments S. – O(k log2 n) expected time for k-th nearest
• Consider the Voronoi regions of s and t in VD(S). neighbors of query point x
• Consider the closest point s0 ∈ VD(S) to s, re- • Minimum Spanning Tree and
spectively t0 ∈ VD(S) to s. TSP-Heuristic, O(n log n)
• Delete the segments of VD(S) where the robot • Largest empty circle, O(n)
does not fit between the obstacles. • Smallest enclosing circle (square with fixed ori-
• Consider the graph G of the remaining diagram entation), O(n)
together with s0 and t0 . • Smallest color spanning circle (square with fixed
• Compute a path from s0 to t0 in G by standard orientation), O(nk), where k is the number of
graph algorithms. colors
Altogether we have the following result: • Localization problems, see Hamacher42

• Clustering of objects, see Dehne and Noltemeier24
Theorem 9 There is a collision-free path of a circle-
shaped robot with radius r in a polygonal environ- All these results stem more or less from the lin-
ment of a set line segments S if and only if the ra- ear complexity of the diagram. As we have already
dius r is smaller than |ss0 | und |tt0 | and a collision- mentioned the complexity of the diagrams in three
free movement from s0 to along the segments of dimension is also linear in many practical situa-
VD(S) exists. tions. Thus many of the presented problems can be
solved in three dimensions with almost the same
5.4.3 Other Applications of the Voronoi time bound. We will present some special applica-
Diagram in 2D tions for 3D.
There are many different geometrical applications

of the Voronoi diagram and its dual. Here we sim- 5.5 Texture Synthesis
ply list some of them, together with some perfor- Textures are visual detail of the rendered geometry,
mance results, provided that the diagram is given: which have become very important in the past few
years, because the cost of rendering with texture is
• Closest Pair of sites, O(n) the same as the cost without texture. Virtually all

N(p)
|| · ||
|| · ||
Td Id
p
|| · ||
T I N(p) N(pi )
p pi
random border T0 I0
Figure 52: The texture synthesis algorithm pro- Figure 53: Using an image pyramid, the texture
ceeds in scan line order through the texture and synthesis process becomes fairly robust against
considers only the neighborhood around the cur- different scales of detail in the sample images.
rent pixel as shown.
Figure 54: Some results of the texture synthesis algorithm.95 In each pair, the image on the left is the
original one, the one on the right is the (partly) synthesized one.
real-world objects have texture, so it is extremely Of course, images not satisfying these criteria can
important to render them in synthetic worlds, too. be used as textures as well (such as façades), but if
Texture synthesis generally tries to synthesize you want to synthesize such images, then a statis-
new textures, either from given images, from a tical or stochastic approach is probably not feasible.
mathematical description, or from a physical model. In the following, we will describe a stochastic
Mathematical descriptions can be as simple as a algorithm that is very simple, very efficient, and
number of sine waves to generate water ripples, works remarkably well.95 Given a sample image, it
while physical models try to describe the physi- does not, like most other methods, try to compute
cal or biological effects and phenomena that lead to explicitly the stochastic model. Instead, it uses the
some texture (such as patina or fur). In all of these sample image itself, which implicitly contains that
“model-based” methods, the knowledge about the model already.
texture is in the model and the algorithm. The We will use the following terminology:
other class of methods starts with one or more
images; then they try to find some statistical or I = Original (sample) image
stochastic description (explicitely or implicitely) of T = New texture image
these, and finally it generates a new texture from pi = Pixel from I
the statistic.
Basically, textures are images with the following p = Pixel from T to be generated next
properties: N(p) = Neighborhood of p (see Figure 52)
1. Stationary: if a window with the proper size is Initially, T is cleared to black. The algorithm
moved about the image, the portion inside the starts by adding a suitably sized border at the left
window always appears the same. and the top, filled with random pixels (this will be
thrown away again at the end). Then, it performs
2. Local: each pixel’s color in the image depends the following simple loop in scan line order (see
only on a relatively small neighborhood. Figure 52):

1: for all p ∈ T do by content” approach, where a query is specified by

2: find the pi ∈ I that minimizes |N(p) − N(pi )|2 providing a (possibly crude) shape, for which the
3: p := pi database is to return best matches.10 The funda-
4: end for mental step here is the matching of shapes, i.e., the
Well, the search in line 2 is exactly a nearest-neighbor calculation of a similarity measure.
search! This can be performed efficiently with the Almost all approaches perform the following steps:
algorithm presented in Section 5.4.1: if N(p) con-
1. Define a transformation function that takes a
tains k pixels, then the points are just 3k-dimen-
shape and computes a so-called feature vector
ional vectors of RGB values, and the distance is just
in some high dimensional space, which (hope-
the Euclidean distance.
fully) captures the shape in its essence. Nat-
Obviously, all pixels of the new texture are deter-
urally, those transformation functions are pre-
ministically defined, once the random border has
ferred that are invariant under rotation and/or
been filled. The shape of the neighborhood N(p)
translation and tessellation.
can be chosen arbitrarily, it must just be chosen
such that all but the current pixel are already com- 2. Define a similarity measure d on the feature
puted. Likewise, other “scans” of the texture are vectors, such that if d( f 1 , f 2 ) is large, then the
possible and sensible (for instance a spiral scan or- associated shapes s1 , s2 do not look similar. Ob-
der), they must just match the shape of N(p). viously, this is (partly) a human factors issue. In
The quality of the texture depends on the size of almost all algorithms, d is just the Euclidean dis-
the neighborhood N(p). However, the optimal size tance.
itself depends on the “granularity” in the sample
image. In order to make the algorithm indepen- 3. Compute a feature vector for each shape in the
dent, we can synthesize an image pyramid (see Fig- database and store them in a data structure that
0 1
ure 53). First, we generate a pyramid I , I , . . . , I d allows for fast nearest-neighbor search.
0
for the sample image I . Then, we synthesize the 4. Given a query, i.e., a shape, compute its feature
texture pyramid T 0 , T 1 , . . . , T d level by level with vector, and retrieve the nearest neighbor from
the above algorithm, starting at the coarsest level. the database. Usually, the system also retrieves
The only difference is that we extend the neigh- all k nearest neighbors. Often times, you are
borhood N(p) of a pixel p over k levels as depicted not interested in the exact k nearest neighbors
by Figure 53. Consequently, we have to build a but only in approximate nearest neighbors (be-
nearest-neighbor search structure for each level, cause the feature vector is an approximation of
because as we proceed downwards in the texture the shape anyway).
pyramid, the size of the neighborhood grows.
Of course, now we have replaced the parameter of The main difference among most shape matching
the best size of the neighborhood by the parameter algorithms is, therefore, the transformation from
of the best size per level and the best number of shape to feature vector.
levels to consider for the neighborhood. However, So, fast shape retrieval essentially requires a fast
as Wei and Levoy95 report, a neighborhood of 9 × (approximate) nearest neighbor search. We could
9 (at the finest level) across 2 levels seems to be stop our discussion of shape matching here, but for
sufficient in almost all cases. sake of completeness, we will describe a very simple
Figure 54 shows two examples of the results that algorithm (from the plethora of others) to compute
can be achieved with this method. a feature vector.71
The general idea is to define some shape func-
5.6 Shape Matching tion f (P1 , . . . , Pn ) → R, which computes some ge-
ometrical property of a number of points, and then
As the availability of 3D models on the net and in evaluate this function for a large number of ran-
databases increases, searching for such models be- dom points that lie on the surface of the shape. The
comes an interesting problem. Such a functionality resulting distribution of f is called a shape distri-
is needed, for instance, in medical image databases, bution.
or CAD databases. One question is how to specify a
query. Usually, most researchers pursue the “query 10This idea seems to originate from image database retrieval,
where it was called QBIC = “query by image content”.

10000 In other words, DS tells for any point p the distance

cube
cylinder to the closest point on the surface S.
8000 pillow
sphere
torus
We can augment the DF further to obtain a vec-
tor distance field VS by storing a vector to the clos-
6000 50
est
¯ point ¯ p ∈ S for each point x. So, DS (x) =
density
¯VS (x)¯. This is a vector field. Another vector field

4000
that is often used in the context of DFs is the gradi-
2000 ent of the distance field, or just gradient field.
Figure 56 shows a distance field for a simple poly-
0
0 1 2 3 4 5 6 7
gon in the plane. The distance of a point from the
distance surface is color-coded (the distance for points inside
Figure 55: The shape distribution of a number of the polygon is not shown). Figure 57 shows a vec-
different simple objects. tor distance field (for the same surface).
Apparently, distance fields (DFs) have been “in-
vented” in many different areas, such as computa-
For the shape function, there are a lot of possi- tional physics,82,84 robotics,54 GIS (see Figures 58
bilities (your imagination is the limit). Examples and 59 for an example), and image processing.73
are: Not surprisingly, DFs come under many other names,
• f (P1 , P2 ) = |P1 − P2 |; such as distance maps and potential fields. The
• f (P1 ) = |P1 − P0 |, where P0 is a fixed point, process of, or the algorithm for computing a DF is
such as the bounding box center; sometimes called distance transform.11
• f (P1 , P2 , P3 ) = ∠(P1 P2 , P1 P3 ); Distance fields have close relationships to isosur-
• f (P1 , P2 , P3 , P4 ) = volume of the tetrahedron faces and implicit functions. When regarded as an
between the four points. implicit function, the isosurface for the iso-value
0 of a DF is exactly the original surface (see Fig-
Figure 55 shows the shape distributions of a few ure 56). However, the converse is not true in gen-
simple objects with the distance between two points eral, i.e., the DF of a surface defined by an implicit
as shape function. function is not necessarily identical to the original
implicit function.
There is another data structure related to distance
6 Distance Fields fields, namely Voronoi diagrams (see Section 5).
Given a vector distance field for a set of points,
Distance fields can be viewed as a way to represent edges, and polygons (not necessarily connected),
surfaces (and solids). They are a very powerful data then all points in space, whose vectors point to the
structure, because they contain a huge amount of same feature (point, edge, or oplygon), are in the
information (compared to just the surface itself). same Voronoi cell (see Figures 56 and 57). (We
This is why they usually take fairly long to com- could also regard the vector of the DF as a kind of
pute. Consequently, this computation can usually ID of the respective Voronoi cell.)
be performed only as a preprocessing step. Thus,
they are difficult to adapt for deformable geometry.
6.1 Computation and representation
Definition 5 (Distance field) of DFs
Let S be a surface in R3 . Then, the distance field of
For special surfaces S, we may be able to compute
a surface S is a scalar function DS : R3 → R such DS analytically. In general, however, we have to
that for all p ∈ R3 , discretize DS spatially, i.e., store them in a 3D voxel
© ª grid, octree, or other space partitioning data struc-
DS (p) = sgn(p) · min d(p, q)|q ∈ S , (14)
ture. Voxel grids and octrees are the most com-
where monly used data structures for storing DFs, and
(
−1, if p inside
sgn(p) = 11Sometimes, this term bears the connotation of producing inac-
+1, if p outside
curate distance fields.

Figure 56: Example of a distance field in the plane Figure 57: The vector distance field for the same
for a simple polygon (thick black lines). The dis- polygon as shown on the left. Only a few vectors
tance field inside the polygon is not shown. The are shown, although (in theory) every point of the
dashed lines show the Voronoi diagram for the field has a vector.
same polygon. The thin lines show isosurfaces
for various isovalues.
Figure 58: A network of roads described in the Figure 59: The distance map of these roads. The
plane by a set of edges. distance of each point in the plane from a road
is color coded. It could be used, for example,
to determine the areas where new houses can be
built.
we will describe algorithms for them in the follow- terpolation. Other interpolation methods could, of
ing. More sophisticated representations try to store course, be used as well.
more information at each voxel in order to be able The simplest discrete representation of a DF is
to extract distances quickly from the field.45 the 3D voxel grid. However, for most applications
Since each cell in a discrete distance field stores this is too costly not only in memory utilization but
one signed distance to the surface for only one also in computational efforts. So, usually DFs are
“representative”, the distance of other points must represented using octrees (see Figure 60), because
be interpolated from these values. One simple they are simple to construct and offer fairly good
method is to store exact distances for the nodes (i.e., adaptivity and allow for easy algorithms.35 This
the corners of the voxels), and generate all other representation is called an adaptively sampled dis-
distances in the interior of voxels by trilinear in- tance field (ADF). Actually, the kind of hierarchical

√ √ √ √ √ √
part of matrix 10 3 8 3 10 3 6 5 6 3
√ √ √ √ √ √ √ √
for backward scan 3 6 5 6 3 6 3 2 3 6
√ √ √ √ √ √ √ √
6 5 2 5 6 5 2 1 2 5
√ √ √ √ √ √ √ √
3 6 5 6 3 6 3 2 3 6
√ √ √ √ √ √ √ √
layer 10 10 8 10 10 3 6 5 6 3
-2 -1
-2 0
√ √ √ √
8 5 2 5 4
-1 √ √ √ √
8 2 1 2 5
0
2 1 0 1 2
√ √ √ √
1 5 2 1 2 8
√ √ √ √
5 5 2 5 8
2
1 2
√ √ √ √ √ √
3 6 5 6 3 10 3 8 3 10
Figure 60: Representing distance fields by octrees √ √
6 3
√ √
2 3
√
6 3
√ √
6 5
√
6 3
√ √ √ √ √ √ √ √
has both memory and computational advantages. part of matrix
5 2
√ √
6 3
1
√ √
2 3
2
√
5
6 3
6 5 2
√ √
6 5
√
5 6
6 3
√ √ √ √ √ √ √ √
for forward scan 3 6 5 6 3 10 10 8 10 10
space partitioning is a design parameter of the al- Figure 61: By convoluting a “distance matrix”
gorithm; for instance, BSPs (see Section 3) proba- with an initially binary distance field (0 or ∞ ),
bly offer much more efficient storage at the price of one gets an approximate distance field. Here a
more complex interpolation. 5 × 5 × 5 example of such a matrix is shown.
Constructing ADFs works much like construct-
ing regular octrees for surfaces, except that here, More formally speaking, since we want to com-
we subdivide a cell if the distance field is not well pute the DF only for a discrete set of points in space,
approximated by the interpolation function defined we can reformulate Equation 14 as follows
so far by the cell corners. ©
The following two algorithms for computing DFs D̃S (x, y, z) = min D(x + i, y + j, z + k)
i,j,k∈Z3
produce only flat representations (i.e., grids), but ª (15)
they can be turned into an ADF by coalescing cells +d M (i, j, k)
in a bottom-up manner if the resulting cell still de- where d M is the distance of a node (i, j, k) from
scribes the DF well enough.12 the center node (0, 0, 0). Actually, this is already
a slight approximation of the exact DF.
6.1.1 Propagation method We can now further approximate this by not con-
sidering the infinite “neighborhood” (i, j, k) ∈ Z3 ,
This method bears some resemblance to region grow-
but instead only a local one I ⊂ Z. Since a prac-
ing and flood filling algorithms and is also called
tical computation of D̃S would compute each node
chamfer method. It can produce only approxima-
in some kind of scan order, we choose I such that
tions of the exact DF.
it contains only nodes that have already been com-
The idea is to start with a binary DF, where all
puted by the respective scan. Thus, d M can be pre-
voxels intersecting the surface are assigned a dis-
computed and stored conveniently in a 3-dimen-
tance of 0, and all others are assigned a distance of
sional matrix (see Figure 61).
∞ . Then, we somehow “propagate” these known
This process looks very similar to conventional
values to voxels in the neighborhood, taking advan-
convolution, except that here we perform a mini-
tage of the fact that we already know the distances
mum operation over a number of sums, while the
between neighboring voxels.13 Note that these dis-
conventional convolution performs a sum over a
tances only depend on the local neighborhood con-
number of products.
stellation, no matter where the “propagation front”
Obviously, a single scan will not assign meaning-
currently is.
ful distance values to all nodes, no matter what or-
der the scan preforms. Therefore, we need at least
two scans: for instance, first, a 3D scanline order
from top to bottom (“forward”), and second, back-
wards from bottom to top. In each pass, only one
half of the 3D matrix d M needs to be considered
12Usually, it is faster to compute an ADF top-down, though. (see Figure 61). It could still happen, of course, that
13This idea has been proposed as early as 1984,13 and maybe even there are still nodes with distance ∞ in the respec-
earlier. tive part of the chamfer matrix.

Figure 62: The distance func- Figure 63: More complex sites Figure 64: The distance func-
tion of a point site in the plane have a bit more complex dis- tion of sites in 3D is different
is a cone.44 tance functions. for the different slices of the
volume.
The time complexity of this method is in O(m), The key technique applied in this approach is to
m = number of voxels. However, the result is only embed the problem into more dimensions than are
an approximation of the distance field; in fact, the inherent in the input/output itself. This is a very
accuracy of the distance values can be pretty low.14 general technique that often helps to get a better
In order to further improve accuracy, one can and simpler view on a complex problem.
make a few more passes, not only from top to bot- The method described in the following was pre-
tom and back, but also from left to right, etc. An- sented by Hoff et al.,44 and, from a different point
other option is to increase the size of d M . of view, by Sethian.83
Analogously, an approximate vector distance field Let us consider the distance field of a single point
can be computed.50 In that case, the matrix d M con- in 2D. If we consider this a surface in 3D, then we
tains vectors instead of scalars. In order to improve get a cone with the apex in that point. In other
the accuracy, we can proceed in two steps: first, for words, the distance field of a point in a plane is just
each voxel in the neighborhood of the surface (usu- the orthogonal projection of suitable cone onto the
ally just the 33 -neighborhood) we compute vec- plane z = 0.
tors to the exact closest point on the surface; then, Now, if there are n point sites in the plane, we
this shell is propagated by several passes with the get n cones, so each point in the plane will be “hit”
vector-valued matrix d M . The time for computing by n distance values from the projection (see Fig-
the exact distance shell can be optimized by utiliz- ure 62). Obviously, only the smallest one wins
ing an octree or BV hierarchy for closest point com- and gets “stored” with that point — and this is ex-
putation, and initializing each search with the dis- actly what the z-buffer of graphics hardware was
tance (and vector) to the closest point of a nieghbor designed for.
voxel. These cones are called the distance function for
a point. The distance function for other features
6.1.2 Projection of distance functions (line segments, polygons, curves) is a little bit more
complicated (see Figure 63).
The previous method can be applied to all kinds of So, computing a discrete DF for a set of sites in
surfaces S. However, it can produce only very gross the plane can be done by rendering all the asso-
approximations. In the following, we describe a ciated distance functions (represented as polygonal
method that can produce much better approxima- meshes) and reading out the z-buffer (and possibly
tions. However, it can be applied only to polygonal the frame buffer if we also want the site IDs, i.e.,
surfaces S. discretized Voronoi regions). If we want to com-
pute a three-dimenional DF, then we proceed slice-
by-slice. One noteworthy catch, though, is that the
14Furthermore, the distance values are not even an upper bound
distance functions of the sites changes from slice to
on the real distance. slice. For instance, the distance function for a point

that lies not in the current slice is a hyperboloid,

with the point site being coincident with the vertex
of the hyperboloid (see Figure 64).
6.2 Applications of DFs

Due to their high information density, distance fields
have a huge number of applications. Aomng them
are motion planning in robotics,54,59 collision de-
tection, shape matching,70 morphing,19 volumet-
ric modeling,14,35,97 navigation in virtual environ-
ments,90 reconstruction,56 offset surface construc-
tion,76 and dynamic LODs, to name but a few. In
the following, we will highlight two easy ones of
them.
6.2.1 Morphing
Figure 65: By specifying correspondences, the
One interesting application of DFs is morphing, i.e.,
user can identify feature points that should transi-
the problem of finding a “smooth” transition from
tion into each other during the morph.20
a shape S ⊂ Rd to a shape T ⊂ Rd , i.e., we are
looking for a shape transformation M(t, Rd ), t ∈
[0, 1], such that M(0, S) = S, M(1, S) = T. to each other, there will be in-between models that
Sometimes, there is a bit of confusion about terms: almost vanish.
sometimes, morphing refers to any smooth transi- Therefore, in addition to the DF interpolation, we
tion, sometimes it refers only to those transitions usually want to split the morphing into two steps
that are bijective, continuous, and have a continu- (for each t): first, S is warped by some function
ous inverse. In the latter case, the term morphing W(t, x) into S0 , then DS0 and DT are interpolated,
is equivalent to the notion of homeomorphism in i.e., M(t, S) = tDW(t,S) + (1 − t)DT .
the area of topology, and the term metamorpho- The idea of the warp is that it encodes the ex-
sis is used to refer to the broader definition. The pertise of the user who usually wants a number of
difference between a homeomorphism and a meta- prominent feature points from S to transition onto
morphosis is that the former does not allow to cut other feature points on T. For instance, she wants
or glue the shape during the transition (this would the corresponding tips of nose, fingers, and feet of
change the topology). two different characters transition into each other.
In the following, we will describe a simple method Therefore, the user has to specify a small number
for metamorphosis of two shapes given as volumet- of correspondences (p0,i , p1,i ), where p0,i ∈ S, and
ric models; it was presented by Cohen-Or et al.20 p1,i ∈ T (see Figure 65). Then, we determine a
The nice thing about a volumetric representation is warp function W(t, x) such that W(1, p0,i ) = p1,i .
that it can naturally handle genus changes (i.e., the Since the warp function should distort S as little
number of “holes” in S and T is different). as possible (just “align” them), we can assemble it
In the simplest form, we can just interpolate lin- from a rotation, then a translation, and finally an
early between the two DFs of S and T yielding a elastic part that is not a linear transformation:
new DF ¡ ¢¡ ¢
W(t, x) = (1 − t)I + tE R(a, tθ)x + tc (17)
M(t, S) = J(t, S, T) = tDS + (1 − t)DT . (16)
where E is the elastic transformation, R(a, θ) is a
In order to obtain a polygonal representation, we rotation about axis a and angle θ, and c is the trans-
can compute the isosurface for the value 0. lation.
This works quite well, if S and T are “aligned” The rotation and translation of the warp function
with each other (in an intuitive sense). However, can be determined by least-squares fitting, while
when morphing two sticks that are perpendicular the elastic transformation can be formulated as a

2 2
1 1
0 0
-1 -1
-2 2 -2
1
0
-1
-2
Figure 66: Example of differ- Figure 67: Example of a 3D dif- Figure 68: Close-up view show-
ence of two DFs. On the left, ference operation.14 ing the adaptivity.
the two DFs for tool and ob-
ject are shown superimposed;
on the right, the result is shown.
scattered data interpolation problem. One method Figure 66 shows an example of the difference oper-
that has proven quite robust is the approach using ation.
radial basis functions. The interested reader will If ADFs are used for representation of the DFs,
find the details in Cohen-Or et al.20 then a recursive procedure has to be applied to
the object’s ADF in order to compute the resulting
6.2.2 Modeling ADF. It is similar to the generation of ADFs (see
above), except that here the subdivision is governed
A frequently used modeling paradigm in CAD is by the compliance of new object ADF with the tool
volumetric modeling, which means that objects are, ADF: a cell is further subdivided if it contains cells
conceptually, defined by specifying for each point from the tool ADF that are smaller, and if the tri-
whether or not it is a member of the object. the ad- linear interpolation inside the cell does not approx-
vantage is that it is fairly simple to define Boolean imate well enough the ADF of the tool.
operations on such a representation, such as union,
subtraction, and intersection. Examples of vol-
umetric modeling is constructive solid geometry 7 Dynamization of Geometric
(CSG) and voxel-based modeling. The advantage
of using distance fields for volumetric modeling is
Data Structures
that they are much easier to implement than CSG,
We present a generic approach for the dynamiza-
but provide much better quality than simple voxel-
tion of an arbitrary static geometric data structure.
based approaches.
Often a simple static data structure is sufficient if
Computing the Boolean operation on two DFs
the set of represented geometric objects will have
basically amounts to evaluating the appropriate op-
few changes over time. Once created, the static
erator on the two fields. The following tables gives
structure mostly has to cope with data queries due
these operators for a field DO (for the object) and a
to its geometric intention. If the set of objects
field DT (for the tool):
varies very much over time, there is need for more
complex dynamic structures which allow efficient
Boolean op. operator on DF insertion and deletion of objects.
¡ ¢
Union DO∪T = min DO , DT For example a one dimensional sorted array of
¡ ¢ a fixed set M is sufficient for x is Element of M
Intersection DO∩T = max DO , DT
¡ ¢ queries. But if the set M has many changes over
Difference DO−T = max DO , −DT
time, a dynamic AVL-tree would be more likely.
The AVL-tree implementation is more complex

since rotations of the tree has to be considered for The dynamization module should export a dy-
insertion and deletion of objects. Additionally, the namic abstract (geometric) data type TDyn with
AVL-tree dynamization was invented for the spe- the following operations:
cial case of the one-dimenional search. We want
to show that it is possible to dynamize a simple Build(W, D): Build the structure W of type
static data structure indirectly but also efficiently TDyn with data objects in
in a general setting. Once this generic approach the set D.
is implemented it can be used for many static data Query(W, q): Gives the answer (objects of D) to
structures. a query to W with query object q.
The generic approaches presented here are not
Extract(W, D): Collects all data objects D of W
optimal against a dynamization adapted directly to
in a single set and returns a
a single data structure, but they are easy to imple-
pointer this set.
ment and efficient for many applications.
In Section 7.1 we formalize the given problem Erase(W): Delete the complete data
and define some requirements. In Section 7.2 we structure W from the storage.
present methods allowing insertion and deletion in Insert(W, d): Insert object d into W.
amortized efficient time. For many applications this
Delete(W, d): Delete d out of W.
is already efficient enough. Within this section
the dynamization technique is explained in detail Note, that the new operations Delete and Insert
and the amortized cost of the new operations are are necessary since we have a dynamic data type
shown. Similar ideas for the worst-case sensitive now.
approach are sketched in Section 7.3. The effort Additionally, we introduce some cost functions
of the dynamization itself is amortized over time. for the operations of the abstract dynamic and the
For details see Klein55 or the work of Overmars72 abstract static data type. For example let BV (n) de-
and van Kreveld.89 We present a simple example in notes the time function for the operation Build(V, D)
Section 7.4. of TStat. The notations are fully presented in Fig-
ure 69. The cost functions depend on the imple-
7.1 Model of the Dynamization mentation of TStat. Note, that the cost function of
TDyn will depend on the cost functions of TStat
Let us assume that TStat is a static abstract (geo-
together with the efficiency of the our general dy-
metric) data type. Since we have a geometric data
namization.
structure, we assume that the essential motivation
of TStat is a query operation on the set of stored
Modul
objects D. i.e., for a query object q the answer is Dynamize
always a subset of D which might be empty.
We want to define the generic dynamization by a
module that imports the following operations from K
Export Import
TStat: ®
Comp.: ADT TDyn
build(V, D): Build the structure V of type
ADT TStat Comp.:
TStat with all data objects in the
BW (n) Build(W, D)
set D.
QW (n) Query(W, q) build(V, D) BV (n)
query(V, q): Gives the answer (objects of D) to EW (n) Extract(W, D) query(V, q) QV (n)
a query to V with query object q. IW (n) Insert(W, d) extract(V, D) EV (n)
extract(V, D): Collects all data objects D of V in DW (n) Delete(W, d)
a single set an returns a pointer to Space.: SW (n) Space: SV (n)
this set.
erase(V): Delete the complete data structure Figure 69: Dynamization in the generic sense.
V from the storage.
In order to guarantee some bounds for the corre-
sponding cost functions of TDyn the cost functions

of TStat must not increase arbitrarily. On the other Build(W, D):

hand for the proof of some time bounds we need Compute binary representation of n = |D|.
some kind of monotonic behavior in the functions,
they should not oscillate. Altogether we define the Decompose D into sets Di with
following requirements which are fulfilled in many |Di | = 2i w.r.t. the representation of n.
cases: Compute build(Vi , Di ) for every Di .
1. QV (n) and EV (n) increase
√ monotonically in n; In principle, the binary structure Wn can be con-
examples: 1, log n, n, n, n log n, n2 , 2n . structed as quick as the corresponding structure V.
BV (n)
2. n and SVn(n) increase monotonically in n; ex- Lemma 10
amples n, n log n, n2 , 2n .
BW (n) ∈ O(BV (n)).
3. For all f ∈ {QV , BV , EV , SV } there is a constant
C
√ ≥ 1, so that f (2n) ≤ C f (n); examples: 1, Computing the binary representation of n and the
n, n, n2 , and also log n with n > 1, as well as decomposition into Di can be done in linear time O(n).
the products of this functions, but not 2n . The operation build(Vi , Di ) needs BV (2i ) time.
We have i ≤ l = blog nc and therefore we con-
4. EV (n) ≤ 2 · BV (n). clude:
Moreover, we assume that the that the query op- blog nc blog nc
BV (2i )
eration can be decomposed, i.e., for a decomposition ∑ BV (2i ) = ∑ 2i
V = V1 ∪ V2 ∪ · · · ∪ Vj of the data set V the results i=0 i=0 2i
of the single operations query(Vi , d) lead to the so- blog nc
BV (n)
lution of query(V, d). This is true for many kinds ≤ ∑ 2i
n
of range queries. i=0
BV (n)
≤ 2log n
7.2 Amortized Insert and Delete n
∈ O(BV (n)).
7.2.1 Amortized Insert: Binary Structure
BV (n)
We used the fact that n increases monotoni-
The very first idea is to implement Insert(W, d)
cally.
and Delete(W, d) directly by
Altogether we have
Insert(W, d) : Extract(W, D); Build(W, D ∪ {d}) BW (n) ∈ O(n + BV (n)) = O(BV (n))
Delete(W, d) : Extract(W, D); Build(W, D \ {d}). since BV (n) is at least linear.
Similar results hold for some other operations.
This throw-away implementation is not very ef-
We can prove
ficient. Therefore we distribute the n data objects of
the static structure V among several structures Vi . EW (n) ≤ log n EV (n)
If a new element has to be inserted we hope that
only a single structure Vi may be concerned. Let the results of extraxt(Vi ) for at most log n struc-
tures Vi .
Additionally,
n = al 2l + al−1 2l−1 + . . . + a1 2 + a0 mit ai ∈ {0, 1}.
QW (n) ≤ log n Qv (n)
Then al al−1 . . . a1 a0 is the binary representation
of n. For every ai = 1 we build a structure Vi which holds if we assume that the query can be decompose
has 2i elements. The collection of these structures as well, see the requirements.
is a representation of Wn which is called binary It is also easy to see that
structure (see Figure 70). To build up the binary
blog nc
structure Wn we proceed as follows:
SW (n) ≤ ∑ SV (2i ) ∈ O(SV (n)).
i=0

V4 V4
4
2
V3
23
V2
22
V1
21
V0
20
W23 W24
Figure 70: The binary structure Wn contains the structure Vi if ai = 1 holds for the
binary representation of n. For examples see n = 23 (left) and n = 24
(right).
Therefore it remains to analyse IW (n). Note, that this is not an expected value and that
As we have seen from Figure 70 sometimes the W is a function of s, i.e., the length of the operation
whole structure of Wn is destroyed when Wn+1 was sequence. The current data set may have a number
build up. In this example we had to perform the of elements n ≤ s.
following tasks: For the Insert(W, d) operation one can prove
µ ¶
extract(V0 , D0 ); extract(V1 , D1 ); extract(V2 , D2 ); log s
I W (s) ∈ O BV (s) .
D := D0 ∪ D1 ∪ D2 ∪ {d}; s
build(V3 , D);
Note, that except insertions there are only queries
In general we have to build up Vj and extract and in s. Queries do not change the size of the data set
erase Vj−1 , Vj−2 , . . ., V0 only if ai = 1 holds for and so we can also replace s by the number of in-
i = 0, 1 . . . , j − 1 and a j = 0 holds (in the binary sertions here.
representation of the current n). Altogether, the results are presented in the fol-
In this special case we have lowing theorem.
Ã !
j−1
Theorem 11 A static abstract data type as pre-
IW (n) ≤ ∑ EV (2i ) + Cj + BV (2 j )
sented in Figure 69 can be dynamized by means
i=0
of the binary structure in a dynamic abstract data
≤ EV (2 j ) + Cj + BV (2 j ) type TDyn so that the operation Insert(W, d) is
³ ´
∈ O BV (2 ) .j performed in amortized time
µ ¶
log s
For a long sequence of insertions many of them I W (s) ∈ O BV (s) .
are performed without extreme reconstructions. Thus s
the effort for all Insert(W, d) is amortized over Let n be the size of the current data set. We have
time.
Generally, let Work be an arbitrary operation SW (n) ∈ log n Sv (n)
with cost function W. For a sequence of s differ- BW (n) ∈ BV (n))
ent operations let Work be applied k times. If
QW (n) ∈ log n Qv (n).
total cost of k Work operationen
≤ W(s)
k 7.2.2 Amortized Delete: Occasional
Reconstruction
holds for a monotonically increasing cost function
W, we say that the operation Work is performed in Assume that we did not have implemented the Insert
amortized time W(s). operation, yet.

If we have to delete an object we can not choose the binary structure for the insertion. The opera-
its location beforehand. Therefore the deletion of tion weak.delete is only available for the structures
an objects is much more difficult than the insertion. Vi and we have to extend it to W in order to apply
Deletion may cause fundamental reconstruction. the result of Section 7.2.2. If Weak.Delete(W, d) is
For many data structures it is easier to simply applied, d should be marked as deleted in W. But
mark an object as deleted. Physically the object re- we do not know in which of the structures Vi the
mains in the structure but does no longer belong to element d lies. Therefore in addition to the binary
the data set D. These objects have bad influence on structure we construct a balanced searchtree T that
the running time of all operations although they stores this information. For every d ∈ W there is a
are no longer necessary. Therefore from time to pointer to the structure Vi with d ∈ Vi , an example
time we have to reconstruct the data structure for is shown in Figure 71.
the actual data set. The additional cost of the search tree T is covered
First of all, for TStat we introduce an addi- as follows. Query operations are not involved. For
tional operation weak.delete(V, d) with cost func- Weak.Delete(W, d) there is an additional O(log n)
tion WDV (n). We simply want to construct a for searching the corresponding Vi and for marking
strong delete function with an acceptable amortized d as deleted in Vi .
time bound for TDyn. If an object d has to be inserted we have to up-
Therefore we use weak.delete(V, d) until D has date T. The object d gets an entry for its struc-
only the half size of V. Then we erase V and build ture Vj in T, this is done in time O(log n) and it
a new structure V out of D. The cost of the oc- will not affect the time bound for the insertion.
casional reconstruction is amortized over the pre- But furthermore if V0 , . . . , Vj−1 has to be erased
ceeding delete-operations. This gives the following the corresponding objects should point to Vj after-
result. wards. This can be efficiently realized by collecting
the pointers of T to Vi in a list for every Vi . We
Theorem 12 A static abstract data type as pre- collect the pointers and change them to "‘Vj "’. This
sented in Figure 69 with an additional operation operation is already covered by time O(BV (2 j )) for
weak.delete(V, d) and with additional cost func- constructing Vj .
tion WDV (n) can be dynamized by means of oc- Altogether we conclude:
casional reconstruction in a dynamic abstract data
type TDyn so that Theorem 13 A static abstract data type as pre-
sented in Figure 69 with an additional operation
BW (r) = BV (r) weak.delete(V, d) and with additional cost function
EW (r) ∈ O (EV (r)) WDV (n) can be dynamized by means of binary
QW (r) ∈ O (QV (r)) structure, searchtree T and occasional reconstruc-
tion in a dynamic abstract data type TDyn so that
SW (r) ∈ O (SV (r))
µ ¶ the amortized time for insertion reads
B (s)
DW (s) ∈ O WDV (s) + V , µ
B (s)
¶
s I W (s) ∈ O log s V ,
s
holds. The size of the current actual data set is de-
noted with r and s denotes the length of the opera- and the amortized time for deletion reads
tion sequence. µ ¶
B (s)
DW (s) ∈ O log s + WDV (s) + V .
We omit the proof here. s
For the rest of the operations we have
7.2.3 Amortized Insert and Amortized Delete
BW (r) = BV (r)
In the preceeding sections we have discussed Insert
and Delete separately. Now we want to show how EW (r) ∈ O (log r EV (r))
to combine the two approaches. QW (r) ∈ O (log r QV (r))
A static abstract data type with a weak delete im- SW (r) ∈ O (SV (r)) .
plementation is given. As in Section 7.2.1 we use

1 2 3 4 5 6 7 8 V3
9 10 11 12 V2
T
13 14 V1
15 V0
Figure 71: A structure W15 with a searchtree T storing pointers to Vi for every d ∈ Vi .
The size of the current actual data set is denoted 7.4 A Simple Example
with r and s denotes the length of the operation
For convenience, we take a simple example from
sequence.
Section 2 and apply Theorem 14, thus implement-
ing worst-case sensitive insertion and deletion.
7.3 Worst-Case sensitive Insert and In Section 2.2 an easy implementation of the
Delete static k-d-tree was presented with
√ Sk-d (n) = O(n)
and query time Qk-d (n) = O( n + a), where a
In the last section we have seen that it is easy to
represents the size of the answer (see Theorem 4).
amortize the cost of Insert and Delete analytically
Obviously, weak.delete(k-d, x) can be implemented
over time. The main idea for the construction of
in O(log n) time, thus we obtain WDk-d (n) =
the dynamic data structure was given by the binary
O(log n). Additionally we have Bk-d (n) = O(n log n).
structure of Wn which has fundamental changes
Let Dyn(k-d) denote the dynamic variant based
from time to time but the corresponding costs were
upon the statically implemented k-d-tree.
amortized. Now we are looking for the worst-
Application of Theorem 14 results in:
case cost of Insert and Delete. The idea is to dis-
tribute the construction of Vj itself over time, i.e., Build(Dyn(k-d), D) ∈ O(n
the structure Vj should be finished if Vj−1 , Vj−2 , . . . √log n)
Query(Dyn(k-d), q) ∈ O(³ n log´n + a)
V0 has to be erased.
We only refer to the result of this approach. The Insert(Dyn(k-d), d) ∈ O log2 n
ideas are very similar to the ideas of the preceed- Delete(Dyn(k-d), d) ∈ O (log n)
ing sections. Technically, a modified binary repre- Space O(n).
sentation is used in order to distribute the effort
of the reconstruction over time. For the interested
reader we refer to Klein55 or van Kreveld89 and
References
Overmars.72
[1] S. Ar, B. Chazelle, and A. Tal, Self-
Customized BSP Trees for Collision Detection,
Theorem 14 A static abstract data type as pre- Computational Geometry: Theory and Applica-
sented in Figure 69 with an additional operation tions, 15 (2000), pp. 91–102.
weak.delete(V, d) and with additional cost function
[2] S. Ar, G. Montag, and A. Tal, Deferred,
WDV (n) can be dynamized in a dynamic abstract
Self-Organizing BSP Trees, in Eurographics, Sept.
data type TDyn so that 2002, pp. 269–278. http://www.ee.technion.
ac.il/~ayellet/papers.html.
Build(W, D) ∈ O(BV (n))
Query(W, q) ∈ O(log [3] J. Arvo and D. B. Kirk, Fast Ray Trac-
³ n · QV (n))
´ ing by Ray Classification, in Computer Graphics
log n
Insert(W, d) ∈ O n BV (n) (SIGGRAPH ’87 Proceedings), M. C. Stone, ed.,
³ ´
BV (n) vol. 21, July 1987, pp. 55–64.
Delete(W, d) ∈ O log n + WDV (n) + n
Space O(SV (n)). [4] F. Aurenhammer, Voronoi diagrams: A sur-
vey of a fundamental geometric data structure,
ACM Comput. Surv., 23 (1991), pp. 345–405.
Here n denotes the number of relevant, stored
data objects.

[5] F. Aurenhammer and R. Klein, [15] P. J. C. Brown, Selective Mesh Refine-

Voronoi Diagrams, in Handbook of Com- ment for Rendering, PhD dissertation, Em-
putational Geometry, J.-R. Sack and J. Ur- manuel College, University of Cambridge, Feb.
rutia, eds., Elsevier Science Publishers B.V. 1998. http://www.cl.cam.ac.uk/Research/
North-Holland, Amsterdam, 2000, pp. 201– Rainbow/publications/pjcb/thesis/.
290. http://wwwpi6.fernuni-hagen.de/ [16] L. P. Chew, Constrained Delaunay triangula-
Publikationen/tr198.pdf. tions, Algorithmica, 4 (1989), pp. 97–108.
[6] J. Avro and D. Kirk, A survey of ray trac- [17] , Guaranteed-quality mesh generation for
ing acceleration techniques, in An Introduction to curved surfaces, in Proc. 9th Annu. ACM Sympos.
Ray Tracing, A. Glassner, ed., Academic Press, San Comput. Geom., 1993, pp. 274–280.
Diego, CA, 1989, pp. 201–262. ISBN 0-12-286160- [18] N. Chin, Partitioning a 3D Convex Polygon
4. with an Arbitrary Plane, in Graphics Gems III,
[7] L. Balmelli, J. Kovacevic, and D. Kirk, ed., Academic Press, 1992, chapter V.2,
M. Vetterli, Quadtrees for Embedded Sur- pp. 219–222.
face Visualization: Constraints and Efficient Data [19] D. Cohen-Or, A. Solomovici, and
Structures, in Proc. of IEEE International Confer- D. Levin, Three-Dimensional Distance Field
ence on Image Processing (ICIP), vol. 2, Oct. 1999, Metamorphosis, ACM Transactions on Graphics,
pp. 487–491. 17 (1998), pp. 116–141. http://visinfo.zib.
[8] L. Balmelli, T. Liebling, and M. Vet- de/EVlib/Show?EVL-1998-152.
terli, Computational Analysis of 4-8 Meshes [20] , Three-Dimensional Distance Field Meta-
with Application to Surface Simplification using morphosis, ACM Transactions on Graphics, 17
global Error, in Proc. of the 13th Canadian Confer- (1998), pp. 116–141. ISSN 0730-0301.
ence on Computational Geometry (CCCG), Aug. [21] M. de Berg, Linear Size Binary Space Parti-
2001. tions for Fat Objects, in Proc. 3rd Annu. European
[9] G. Barequet, B. Chazelle, L. J. Sympos. Algorithms, vol. 979 of Lecture Notes
Guibas, J. S. B. Mitchell, and A. Tal, Comput. Sci., 1995, pp. 252–263.
BOXTREE: A Hierarchical Representation for [22] , Linear size binary space partitions for un-
Surfaces in 3D, Computer Graphics Forum, 15 cluttered scenes, Algorithmica, 28 (2000), pp. 353–
(1996), pp. C387–C396, C484. ISSN 0167-7055. 366.
[10] N. Beckmann, H.-P. Kriegel, [23] M. de Berg, M. van Kreveld,
R. Schneider, and B. Seeger, The M. Overmars, and O. Schwarzkopf,
R∗ -tree: An efficient and robust access method Computational Geometry: Algorithms and Ap-
for points and rectangles, in Proc. ACM SIGMOD plications, Springer-Verlag, Berlin, Germany,
Conf. on Management of Data, 1990, pp. 322–331. 2nd ed., 2000.
[11] J. Bernal, Bibliographic notes on Voronoi di- [24] F. Dehne and H. Noltemeier, A com-
agrams, tech. rep., National Institute of Standards putational geometry approach to clustering prob-
and Technology, Gaithersburg, MD 20899, 1992. lems, in Proc. 1st Annu. ACM Sympos. Comput.
[12] J.-D. Boissonnat and M. Teillaud, Geom., 1985, pp. 245–250.
On the randomized construction of the Delaunay [25] A. K. Dewdney and J. K. Vranch, A
tree, Theoret. Comput. Sci., 112 (1993), pp. 339– convex partition of R3 with applications to Crum’s
354. http://www.inria.fr/cgi-bin/wais_ra_ problem and Knuth’s post-office problem, Utilitas
sophia?question=1140. Math., 12 (1977), pp. 193–199.
[13] G. Borgefors, Distance transformations [26] H. Djidjev and A. Lingas, On computing
in arbitrary dimensions, in Computer. Vision, the Voronoi diagram for restricted planar figures,
Graphics, Image Processing, vol. 27, 1984, in Proc. 2nd Workshop Algorithms Data Struct.,
pp. 321–345. vol. 519 of Lecture Notes Comput. Sci., 1991,
pp. 54–64.
[14] P.-T. Bremer, S. D. Porumbescu,
F. Kuester, B. Hamann, K. I. Joy, [27] D. P. Dobkin and M. J. Laszlo, Prim-
and K.-L. Ma, Virtual Clay Modeling us- itives for the manipulation of three-dimensional
ing Adaptive Distance Fields, in Proceedings of subdivisions, Algorithmica, 4 (1989), pp. 3–32.
the 2002 International Conference on Imaging [28] R. A. Dwyer, Higher-dimensional Voronoi di-
Science, Systems, and Technology (CISST 2002), agrams in linear expected time, Discrete Comput.
H. R. A. et al., ed., vol. 1, Athens, Georgia, 2002. Geom., 6 (1991), pp. 343–367.

[29] H. Edelsbrunner, Algorithms in Combina- in SIGGRAPH 96 Conference Proceedings,

torial Geometry, vol. 10 of EATCS Monographs H. Rushmeier, ed., ACM SIGGRAPH, Aug. 1996,
on Theoretical Computer Science, Springer- pp. 171–180. held in New Orleans, Louisiana,
Verlag, Heidelberg, West Germany, 1987. 04-09 August 1996.
[30] H. Edelsbrunner and N. R. Shah, In- [41] L. J. Guibas and J. Stolfi, Primitives for
cremental topological flipping works for regular the manipulation of general subdivisions and the
triangulations, Algorithmica, 15 (1996), pp. 223– computation of Voronoi diagrams, ACM Trans.
241. Graph., 4 (1985), pp. 74–123.
[31] S. A. Ehmann and M. C. Lin, Accurate [42] H. W. Hamacher, Mathematische Lö-
and Fast Proximity Queries Between Polyhedra sungsverfahren für planare Standortprobleme,
Using Convex Surface Decomposition, in Com- Verlag Vieweg, Wiesbaden, 1995.
puter Graphics Forum, vol. 20, 2001, pp. 500–510. [43] I. Hoff, Kenneth E., T. Culver,
ISSN 1067-7055. J. Keyser, M. Lin, and D. Manocha,
[32] A. A. Elassal and V. M. Caruso, USGS Fast Computation of Generalized Voronoi Dia-
digital cartographic data standards - Digital Ele- grams Using Graphics Hardware, pp. 375–376.
vation Models, Technical Report Geological Sur- [44] K. E. Hoff III, T. Culver, J. Keyser,
vey Circular 895-B, US Geological Survey, 1984. M. Lin, and D. Manocha, Fast Com-
[33] D. A. Field, Implementing Watson’s algo- putation of Generalized Voronoi Diagrams Us-
rithm in three dimensions, in Proc. 2nd Annu. ing Graphics Hardware, LA, California, Aug.8–
ACM Sympos. Comput. Geom., 1986, pp. 246– 13 1999, pp. 277–286. http://www.cs.unc.edu/
259. ~geom/voronoi/.
[34] S. Fortune, Voronoi diagrams and Delaunay [45] J. Huang, Y. Li, R. Crawfis, S. C. Lu,
triangulations, in Computing in Euclidean Geom- and S. Y. Liou, A complete distance field rep-
etry, D.-Z. Du and F. K. Hwang, eds., vol. 1 of Lec- resentation, in Proceedings of the conference on
ture Notes Series on Computing, World Scientific, Visualization 2001, 2001, pp. 247–254. ISBN 0-
Singapore, 1st ed., 1992, pp. 193–233. 7803-7200-X.
[35] S. F. Frisken, R. N. Perry, A. P. [46] P. M. Hubbard, Real-Time Collision Detec-
Rockwood, and T. R. Jones, Adap- tion and Time-Critical Computing, in SIVE 95,
tively Sampled Distance Fields: A General Rep- The First Worjshop on Simulation and Interaction
resentation of Shape for Computer Graphics, in in Virtual Environments, no. 1, Iowa City, Iowa,
Siggraph 2000, Computer Graphics Proceedings, July 1995, University of Iowa, pp. 92–96.
K. Akeley, ed., Annual Conference Series, 2000, [47] H. Inagaki, K. Sugihara, and
pp. 249–254. http://visinfo.zib.de/EVlib/ N. Sugie, Numerically robust incremental
Show?EVL-2000-61. algorithm for constructing three-dimensional
[36] H. Fuchs, Z. M. Kedem, and B. F. Voronoi diagrams, in Proc. 4th Canad. Conf.
Naylor, On Visible Surface Generation by Comput. Geom., 1992, pp. 334–339.
a Priori Tree Structures, in Computer Graphics [48] B. Joe, 3-Dimensional Triangulations from Lo-
(SIGGRAPH ’80 Proceedings), vol. 14, July 1980, cal Transformations, SIAM J. Sci. Statist. Com-
pp. 124–133. put., 10 (1989), pp. 718–741.
[37] D. Fussell and K. R. Subramanian, [49] , Construction of Three-Dimensional De-
Fast Ray Tracing Using K-D Trees, Technical Re- launay Triangulations Using Local Transforma-
port TR-88-07, U. of Texas, Austin, Dept. Of Com- tions, Comput. Aided Geom. Design, 8 (1991),
puter Science, Mar. 1988. pp. 123–142.
[38] A. S. Glassner, ed., An Introduction to Ray [50] M. W. Jones and R. A. Satherley,
Tracing, Academic Press, 1989. Shape representation using space filled sub-voxel
[39] J. Goldsmith and J. Salmon, Automatic distance fields, pp. 316–325.
Creation of Object Hierarchies for Ray Trac- [51] T. C. Kao and D. M. Mount, Incre-
ing, IEEE Computer Graphics and Applications, 7 mental construction and dynamic maintenance of
(1987), pp. 14–20. constrained Delaunay triangulations, in Proc. 4th
[40] S. Gottschalk, M. Lin, and Canad. Conf. Comput. Geom., 1992, pp. 170–175.
D. Manocha, OBB-Tree: A Hierarchical
Structure for Rapid Interference Detection,

[52] N. Katayama and S. Satoh, The SR-tree: Detail Rendering of Height Fields, in SIGGRAPH
An Index Structure for High-Dimensional Near- 96 Conference Proceedings, H. Rushmeier, ed.,
est Neighbor Queries, in Proc. ACM SIGMOD ACM SIGGRAPH, Aug. 1996, pp. 109–118. held
Conf. on Management of Data, 1997, pp. 369–380. in New Orleans, Louisiana, 04-09 August 1996.
[53] T. L. Kay and J. T. Kajiya, Ray Trac- [65] P. Lindstrom and V. Pascucci, Visual-
ing Complex Scenes, in Computer Graphics (SIG- ization of Large Terrains Made Easy, in Proc. IEEE
GRAPH ’86 Proceedings), D. C. Evans and R. J. Visualization, San Diego, 2001.
Athay, eds., vol. 20, Aug. 1986, pp. 269–278. [66] W. E. Lorensen and H. E. Cline,
[54] R. Kimmel, N. Kiryati, and A. M. Marching Cubes: A High Resolution 3D Surface
Bruckstein, Multi-Valued Distance Maps for Construction Algorithm, in Computer Graphics
Motion Planning on Surfaces with Moving Obsta- (SIGGRAPH ’87 Proceedings), M. C. Stone, ed.,
cles, IEEE Transactions on Robotics and Automa- vol. 21, July 1987, pp. 163–169.
tion, 14 (1998), pp. 427–436. [67] G. Mller, S. Schfer, and W. D. Fell-
[55] R. Klein, Algorithmische Geometrie, Addison- ner, Automatic Creation of Object Hierarchies
Wesley, Bonn, 1997. http://www.oldenbourg. for Radiosity Clustering, Computer Graphics Fo-
de/cgi-bin/rotitel?T=24382. rum, 19 (2000). ISSN 0167-7055.
[56] R. Klein, A. Schilling, and [68] B. Naylor, J. Amanatides, and
W. Straer, Reconstruction and simplification W. Thibault, Merging BSP Trees Yields Poly-
of surfaces from contours, in Pacific Graphics, hedral Set Operations, in Computer Graphics
1999, pp. 198–207. http://cg.cs.uni-bonn. (SIGGRAPH ’90 Proceedings), F. Baskett, ed.,
de/publications/publication.asp?id=31. vol. 24, Aug. 1990, pp. 115–124.
[57] J. T. Klosowski, M. Held, J. S. B. [69] B. F. Naylor, A Tutorial on Binary Space
Mitchell, H. Sowrizal, and K. Zikan, Partitioning Trees, ACM SIGGRAPH ’96 Course
Efficient Collision Detection Using Bounding Vol- Notes 29, (1996).
ume Hierarchies of k-DOPs, IEEE Transactions on
[70] M. Novotni and R. Klein, A Geometric
Visualization and Computer Graphics, 4 (1998),
Approach to 3D Object Comparison, in Interna-
pp. 21–36.
tional Conference on Shape Modeling and Appli-
[58] T. Larsson and T. Akenine-Mller, cations, May 2001, pp. 167–175.
Collision Detection for Continuously Deforming
[71] R. Osada, T. Funkhouser,
Bodies, in Eurographics, 2001, pp. 325–333. short
B. Chazelle, and D. Dobkin, Matching
presentation.
3D models with shape distributions, in Pro-
[59] J.-C. Latombe, Robot Motion Planning, ceedings of the International Conference on
Kluwer Academic Publishers, Boston, 1991. Shape Modeling and Applications (SMI-01),
[60] C. L. Lawson, Software for C1 surface inter- B. Werner, ed., Los Alamitos, CA, May 7–11 2001,
polation, in Math. Software III, J. R. Rice, ed., Aca- pp. 154–166.
demic Press, New York, NY, 1977, pp. 161–194. [72] M. H. Overmars, The Design of Dynamic
[61] D. T. Lee and A. K. Lin, Generalized De- Data Structures, vol. 156 of Lecture Notes Com-
launay Triangulation for Planar Graphs, Discrete put. Sci., Springer-Verlag, Heidelberg, West Ger-
Comput. Geom., 1 (1986), pp. 201–217. many, 1983.
[62] S. Leutenegger, J. Edgington, and [73] D. W. Paglieroni, Distance Transforms:
M. Lopez, STR : A Simple and Efficient Algo- Properties and Machine Vision Applications,
rithm for R-Tree Packing, in Proceedings of the CVGIP: Graphical Models and Image Processing,
13th International Conference on Data Engineer- 54 (1992), pp. 56–74.
ing (ICDE’97), Washington - Brussels - Tokyo, [74] I. J. Palmer and R. L. Grimsdale, Colli-
Apr. 1997, pp. 497–507. ISBN 0-8186-7807-0. sion Detection for Animation using Sphere-Trees,
[63] T.-Y. Li and J.-S. Chen, Incremental 3D Computer Graphics Forum, 14 (1995), pp. 105–
Collision Detection with Hierarchical Data Struc- 116. ISSN 0167-7055.
tures, in Proc. VRST ’98, Taipei, Taiwan, Nov. [75] M. S. Paterson and F. F. Yao, Effi-
1998, ACM, pp. 139–144. cient binary space partitions for hidden-surface
[64] P. Lindstrom, D. Koller, W. Rib- removal and solid modeling, Discrete Comput.
arsky, L. F. Hughes, N. Faust, and Geom., 5 (1990), pp. 485–503.
G. Turner, Real-Time, Continuous Level of

[76] B. A. Payne and A. W. Toga, Dis- [88] G. van den Bergen, Efficient Collision
tance field manipulation of surface models, IEEE Detection of Complex Deformable Models using
Computer Graphics and Applications, 12 (1992), AABB Trees, Journal of Graphics Tools, 2 (1997),
pp. 65–71. pp. 1–14.
[77] F. P. Preparata and M. I. Shamos, [89] M. J. van Kreveld, New Results on Data
Computational Geometry: An Introduction, Structures in Computational Geometry, ph.D.
Springer-Verlag, New York, NY, 1985. dissertation, Dept. Comput. Sci., Utrecht Univ.,
Utrecht, Netherlands, 1992.
[78] V. T. Rajan, Optimality of the Delaunay tri-
angulation in Rd , in Proc. 7th Annu. ACM Sym- [90] M. Wan, F. Dachille, and A. Kauf-
pos. Comput. Geom., 1991, pp. 357–363. man, Distance-field based skeletons for virtual
navigation, in Proceedings of the conference on
[79] J. Revelles, C. Urena, and M. Las- Visualization 2001, 2001, pp. 239–246. ISBN 0-
tra, An Efficient Parametric Algorithm for Oc- 7803-7200-X.
tree Traversal, in WSCG 2000 Conference Pro-
ceedings, University of West Bohemia, Plzen,
[91] C. A. Wang, Efficiently updating the con-
strained Delaunay triangulations, BIT, 33 (1993),
Czech Republic, 2000. http://visinfo.zib.de/
pp. 238–252.
EVlib/Show?EVL-2000-307.
[92] C. A. Wang and L. Schubert, An optimal
[80] N. Roussopoulos and D. Leifker, algorithm for constructing the Delaunay triangu-
Direct spatial search on pictorial databases using
lation of a set of line segments, in Proc. 3rd Annu.
packed R-trees, in Proceedings of ACM-SIGMOD
ACM Sympos. Comput. Geom., 1987, pp. 223–
1985 International Conference on Management of
232.
Data, May 28–31, 1985, LaMansion Hotel, Austin,
Texas, S. Navathe, ed., New York, NY 10036, USA, [93] D. F. Watson, Computing the n-Dimensional
1985, pp. 17–31. ISBN 0-89791-160-1. http:// Delaunay Tesselation with Applications to
www.acm.org/pubs/articles/proceedings/ Voronoi Polytopes, Comput. J., 24 (1981),
mod/318898/p17-roussopoulos/ pp. 167–172.
p17-roussopoulos.pdf;http://www.acm. [94] H. Weghorst, G. Hooper, and D. P.
org/pubs/citations/proceedings/mod/ Greenberg, Improved Computational Meth-
318898/p17-roussopoulos/. ods for Ray Tracing, ACM Transactions on Graph-
ics, 3 (1984), pp. 52–69.
[81] R. Seidel, Constrained Delaunay triangula-
tions and Voronoi diagrams with obstacles, Tech- [95] L.-Y. Wei and M. Levoy, Fast Texture
nical Report 260, IIG-TU Graz, Austria, 1988. Synthesis Using Tree-Structured Vector Quanti-
zation, in Siggraph 2000, Computer Graphics Pro-
[82] J. A. Sethian, An Analysis of Flame Propaga- ceedings, K. Akeley, ed., Annual Conference Se-
tion, PhD thesis, Dept. of Mathematics, University ries, 2000, pp. 479–488. http://visinfo.zib.
of California, Berkeley, 1982. http://visinfo. de/EVlib/Show?EVL-2000-87.
zib.de/EVlib/Show?EVL-1982-2.
[96] J. Wilhelms and A. V. Gelder, Octrees
[83] , Level Set Methods, Cambridge University for Faster Isosurface Generation Extended Ab-
Press, 1996. http://visinfo.zib.de/EVlib/ stract, in Computer Graphics (San Diego Work-
Show?EVL-1996-149. shop on Volume Visualization), vol. 24, Nov. 1990,
[84] , Level set methods and fast marching pp. 57–62.
methods, Cambridge University Press, 1999. [97] C. Youngblut, R. E. Johnson,
[85] M. I. Shamos, Computational Geometry, S. H. Nash, R. A. Wienclaw, and
ph.D. thesis, Dept. Comput. Sci., Yale Univ., New C. A. Will, Different Applications of
Haven, CT, 1978. Two-Dimensional Potential Fields for Vol-
ume Modeling, Tech. Rep. UCAM-CL-TR-
[86] M. Tanemura, T. Ogawa, and 541, University of Cambridge, Computer
W. Ogita, A New Algorithm for Three- Laboratory, 15 JJ Thomson Avenue, Cam-
Dimensional Voronoi Tesselation, J. Comput. bridge CB3 0FD, United Kingdom, Aug. 2002.
Phys., 51 (1983), pp. 191–207. http://www.cl.cam.ac.uk/users/nad/pubs/.
[87] E. Torres, Optimization of the Binary Space [98] G. Zachmann, Real-time and Exact Collision
Partition Algorithm (BSP) for the Visualization of Detection for Interactive Virtual Prototyping, in
Dynamic Scenes, in Eurographics ’90, C. E. Van- Proc. of the 1997 ASME Design Engineering Tech-
doni and D. A. Duce, eds., Sept. 1990, pp. 507–518. nical Conferences, Sacramento, California, Sept.
1997. Paper no. CIE-4306.

[99] , Rapid Collision Detection by Dynamically China, Nov.11–13 2002, pp. 121–128. http://
Aligned DOP-Trees, in Proc. of IEEE Virtual Real- www.gabrielzachmann.org/.
ity Annual International Symposium; VRAIS ’98, [101] B. Zhu and A. Mirzaian, Sorting does not
Atlanta, Georgia, Mar. 1998, pp. 90–97. always help in computational geometry, in Proc.
3rd Canad. Conf. Comput. Geom., Aug. 1991,
[100] , Minimal Hierarchical Collision Detection, pp. 239–242.
in Proc. ACM Symposium on Virtual Reality
Software and Technology (VRST), Hong Kong,

Geometric Data Structures For Computer Graphics: Gabriel Zachmann & Elmar Langetepe

Uploaded by

Copyright:

Available Formats

Geometric Data Structures For Computer Graphics: Gabriel Zachmann & Elmar Langetepe

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Geometric Data Structures For Computer Graphics: Gabriel Zachmann & Elmar Langetepe

Uploaded by

Copyright:

Available Formats

Geometric Data Structures for Computer Graphics ∗

Gabriel Zachmann & Elmar Langetepe

13th August 2003

Siggraph 2003 Tutorial 16

2 Quadtrees and K-d-Trees 5

4 Bounding Volume Hierarchies 22

Siggraph 2003 Tutorial 16

7 Dynamization of Geometric Data Structures 44

Siggraph 2003 Tutorial 16

1 Introduction A natural generalization of the one-dimenional

of this course. Neither do we try to present leading

PNE := {p ∈ P : p x > xmid ∧ py > ymid },

Siggraph 2003 Tutorial 16

2.1.2 Complexity and Construction

Theorem 1 The depth of a quadtree for a set P of

2.1.3 Neighbor Finding 2.2 K-d-Trees

Siggraph 2003 Tutorial 16

Figure 3: A k-d-tree for k = 2 and a rectangular 2.3 Height Field Visualization

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

In order to determine the order in which sub- y 3 7

us which cell is next. 0 4

Siggraph 2003 Tutorial 16

current exit side

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

Figure 23: An example BSP tree for a set of ob-

Figure 25: BSP trees are an efficient data structure

However, all of these examples producing quadratic

Siggraph 2003 Tutorial 16

der of all polygons? Using BSP trees, this is al-

Siggraph 2003 Tutorial 16

H First, we need to define some nomenclature:

T − , T + = left and right child of T, resp.

Finally, we would like to define a node T by the

Siggraph 2003 Tutorial 16

split merge merge

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

1. when do we split a preliminary leaf? (and which

5 However, we will argue at the end that BV hierarchies are just at

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

B the following that AABBs are used as BVs. How-

Siggraph 2003 Tutorial 16

Siggraph 2003 Tutorial 16

B a simple incremental construction approach which

Siggraph 2003 Tutorial 16

Definition 5 For p, q ∈ S let vertices; they belong to the common boundary of

The Voronoi diagram is a simple linear structure

5.1.2 Delaunay Triangulation