Delaunay triangulation
www.comp.lancs.ac.uk
en.wikipedia.org
jonathanpuckey.com
yagosweb.blogspot.com
Delaunay triangulation
• Triangulation aiming to preserve the triangles
to be as equilateral as possible (in such a
representation, each triangle represents the
local value on the surface in the best way)
• It is unique
– Independent on the starting point or the
orientation of the input dataset
– If 4 and more points are not lying on a circle
Delaunay triangulation
• Input: P = {p1, p2, …, pn}
• Output: Triangulation T for P
• Definition of triangulation T for P represents the
space division into the set of m triangles T = {t1,
t2, …, tm} which fulfill:
– Two arbitrary triangles can share maximally one edge
– The union of all triangles from T forms the convex hull
of P
– None of the triangles contains another point from P
Active Edge List (AEL)
• Data structure often used for construction of DT
• Contains the topology of the DT triangles
• Lets consider two adjacent triangles ti, tj from DT,
sharing one edge marked as eij in ti and as eji in tj
• Each edge eij (Active Edge) in ti triangle oriented
counter-clockwise keeps:
– Pointer to the following edge ei+1 in ti
– Pointer to edge eji from the adjacent triangle tj
Active Edge List (AEL)
• Except for edges lying on the convex hull H, each
edge e from DT is represented twice (as eij and
eji), with different orientations
• These doubled edges are called twin edges
• Each triangle is then described by a triplet of
edges (eij, ei+1j,ei+2j) with counter-clockwise
orientation and forming a Circular List
• The list of all such edges forms the Active Edge
List
Active Edge List (AEL)
VD construction – algorithms
• Direct construction:
– Local switching
– Incremental approach
– Divide and conquer
• Indirect construction:
– Via Voronoi diagram
Local switching
• Modifying of a general triangulation to DT
• Based on switching the “illegal” edges in
adjacent triangles forming a convex quad
• Complexity O(n2)
Local switching
Edge legalization
• Edge flip = swapping the quad diagonals
• The resulted triangles are both legal = locally
optimal according to the selected criterion
Edge legalization
• Typical criteria:
– Minimization of the maximal angle
– Vertices lying inside a circumscribed circle of the
triangle
– Minimal/maximal triangle height v
– Minimal/maximal area of triangle S
– …
Incremental approach
• Can be used in 2D and 3D
• Incremental addition of points into already
created DT
• For already existing Delaunay edge e = p1p2 we
search for such a point p, which has the minimal
Delaunay distance dD(p1p2 , p) from p1p2
• Each Delaunay edge is oriented, the point p is
searched only on the left side from this edge
• We use the test for orientation of the triangle
vertices if it is counter-clockwise (determinant
test)
Incremental approach
• We add edges of triangle (p1 , p2 , p) to DT
• If such a point p does not exist (the examined
edge lies on the convex hull), we change the
edge orientation and repeat the search
• Complexity O(n2)
Delaunay distance
• Let k(S, r) be a circle and l a line intersecting
with k in points a, b and p point lying on k
• Delaunay distance of point p from edge a,b is
marked as dD(h, p)
Incremental approach
• When constructing we can use the modified
AEL structure:
– It contains edges e for whose we are searching for
points p, it doesn’t store the topology model
Incremental approach
Incremental approach
Incremental approach
Pseudocode
Pseudocode
Incremental insertion method
• Uses so-called simplex (bounding triangle)
• Frequent method for DT construction
• Complexity O(n2)
• Principle:
– In each step we add one point to DT and perform
the legalization of DT
Incremental insertion method
• Input: set P = {p0, p1, …, pn} of points in a plane
• Select p0 as a point with the highest y-axis
value (or also the x-axis)
• We add two other points p-1 (sufficiently low
and far away to the right) and p-2 (sufficiently
high and far away to the left) so that P lies
inside the triangle p0 p-1 p-2
Incremental insertion method
• We create the DT sets {p-2 ,p-1, p0, p1, …, pn}
and at the end we remove all edges containing
points p-2 and p-1
• DT for the set {p-2 ,p-1, p0 } is the triangle {p-2 ,
p-1, p0 }
Incremental insertion method
• We don’t want to determine the exact
position of p-2, p-1 , so for determining the
position of pj wrt. the oriented line we use the
following equivalence:
Step 7 – finding the triangle containing
p
• The most computationally demanding step (it
is not efficient to search for p in all triangles)
• The most common methods:
– Walking method (heuristic method, O(n2))
– DAG tree (ternary tree construction, O(n log n))
Walking method
• By traversing the adjacent triangles we are
gradually approaching the searched triangle ti
• We are testing the mutual position of p and
edge eij in AEL.
• Point p lies on the left side from all edges of
the searched triangle
Divide and conquer
• Input set of points is divided into smaller
parts, each of them is triangulated separately
• Resulting triangulations are merged and
legalized
Assignment
• Implement the Delaunay triangulation using
the incremental approach
Useful details for implementation
• We have to be able to determine the
circumscribed circle = circle containing three
vertices
• We can do this in the following way:
– Create a class RealPoint(float x, float y)
• Its distance method calculates the distance between
points p1 and p2:
– sqrt((p1.x - p2.x)2 + (p1.y - p2.y)2)
Useful details for implementation
• Class Circle is determined by its center
(RealPoint c) and radius (float r)
• Testing if a point p lies inside a circle:
– Method inside
• if (c.distanceSq(p) < r2) return true;
where distanceSq = (p1.x - p2.x)2 + (p1.y - p2.y)2
Useful details for implementation
• Calculating the circle with three points lying on it
(RealPoint p1, p2, p3):
– Method circumCircle(p1, p2, p3)
cp = crossproduct (p1, p2, p3);
if (cp <> 0) {
p1Sq = p1.x2 + p1.y2;
p2Sq = p2.x2 + p2.y2;
p3Sq = p3.x2 + p3.y2;
num = p1Sq *(p2.y - p3.y) + p2Sq *(p3.y - p1.y) +
p3Sq *(p1.y - p2.y);
cx = num / (2.0 * cp);
num = p1Sq *(p3.x - p2.x) + p2Sq*(p1.x - p3.x) +
p3Sq*(p2.x - p1.x);
cy = num / (2.0f * cp); c.set(cx, cy);
c.set(cx, cy);
r = c.distance(p1);
Useful details for implementation
• crossproduct (p1, p2, p3)>
u1 = p2.x() - p1.x();
v1 = p2.y() - p1.y();
u2 = p3.x() - p1.x();
v2 = p3.y() - p1.y();
return u1 * v2 - v1 * u2;