CIS677: Notes for lecture 6: Delaunay triangulations, Convexity

CIS677
Notes for lecture 6:
Delaunay triangulations, Convexity

Outline

Voronoi diagrams and Delaunay triangulations
- edge flipping and angle optimality
- uses: interpolation, geometric graphs
- divide and conquer algorithm
- sweep algorithm
Convexity
- convex, affine, linear, conical combinations
- Caratheodory's theorem
- polarity and Duality
- the Farkas Lemma
- convex hulls and intersections of halfspaces

Edge Flipping and Angle Optimality

Suppose a triangulation is not locally Delaunay, so some triangles abc and abd have site d inside the circumcircle of abc.

The quadrilateral acbd is convex, and has two triangulations:

T = {abc, abd}, with edge ab, the original one, and
T' = {cda, cdb}, with edge cd.

Facts:

T' is Delaunay (since the other one isn't)
Angle increase: The minimum angle of T' (of the six angles) is larger than the minimum angle of T.

Edge flip: The transformation of the triangulation
deleting {abc, abd} and inserting {cda, cdb},
that is, swapping cd for ab.

Edge flip:	The transformation of the triangulation deleting {abc, abd} and inserting {cda, cdb}, that is, swapping cd for ab.

Note: if {a,b,c,d} are co-circular, both triangulations are Delaunay and have the same smallest angle.

The edge flipping algorithm

...for computing a Delaunay triangulation:

	Find any triangulation of S;
	repeat
		find non-locally-Delaunay pair of triangles
		flip edge
	until the triangulation is locally Delaunay

Minimum angle of a triangulation: the smallest angle of any triangle in the triangulation

Theorem: The edge flipping algorithm terminates, and the Delaunay triangulation maximizes the minimum angle

Proof: Given a triangulation T, let

Angle(T) = (a₁, a₂,...) be the sequence of measures of angles of triangles of T, sorted into nondecreasing order.

Angle(T) =	(a₁, a₂,...) be the sequence of measures of angles of triangles of T, sorted into nondecreasing order.

Suppose T is not Delaunay, and so not locally Delaunay, and so there is a triangulation T' resulting from an edge flip.

Then Angle(T') is lexicographically larger than Angle(T).

Hence if T is not Delaunay, Angle(T') is larger than Angle(T),
so the Delaunay triangulation is lexicographically maximum.

Since the number of triangulations is finite,
and each triangulation generated is distinct,
the algorithm terminates.
QED

Other properties and uses of the Delaunay triangulation

Interpolation

The Delaunay triangulation is frequently used for "piecewise-linear" interpolation of scattered data: given a set of points in the plane, and function values at those points, the interpolated value of the function at a point p is the value found using linear interpolation within the triangle containing p.
That is, the graph of the interpolant function is:
for each triangle abc of the Delaunay triangulation, lift it to the 3-dimensional triangle
[a f(a)][b f(b)][c f(c)].

What is the best triangulation for this interpolation method?

It depends on the function being interpolated, which is unknown in general (except at the sites). However:

The Delaunay triangulation minimizes the maximum error when the function is
f(x,y) = x² + y² or (x-a)²+(y-b)²
for any fixed a,b.

The Delaunay triangulation minimizes the "roughness" of the interpolant, the integral of the (squared) length of the gradient.

For some functions, the Delaunay triangulation is not the best; better triangulations are often found by local improvement using edge flipping.

Geometric graphs

For a given set S of sites, the geometric graph G(S) on S has:

the sites as vertices, and
has all edges {a,b} between sites a and b,
with weight equal to the distance between the a and b.

For a given set S, the following undirected graphs have vertex set S, and weighted edges {a,b}, for a,b in S, with weight d(a,b). Let:

G(S) have all edges.

RNG(S) , the Relative Neighborhood Graph, have edges between sites a and b when no third site is closer to a and to b than d(a,b). That is, the shaded lune in the figure is empty.

GG(S) , the Gabriel graph, have edges between sites a and b when no site is closer than a and b to the midpoint of segment ab. The shaded circle is empty.

MST(S) be the minimum spanning tree of G(S).

G(S)	have all edges.
RNG(S)	, the Relative Neighborhood Graph, have edges between sites a and b when no third site is closer to a and to b than d(a,b). That is, the shaded lune in the figure is empty.
GG(S)	, the Gabriel graph, have edges between sites a and b when no site is closer than a and b to the midpoint of segment ab. The shaded circle is empty.
MST(S)	be the minimum spanning tree of G(S).

Fact: the edge sets of these graphs have
MST(S) subset RNG(S) subset GG(S) subset DG(S) subset G(S).

Proof sketch: Recall that an edge e of a graph is in no MST of the graph if e is the heaviest edge in some cycle of the graph. If an edge {a,b} is not in RNG(S) due to some site c, then {a,b} is heaviest in the cycle a-b-c-a. For the rest, proof by picture. QED

Spanners

Given a graph G and parameter t, a t-spanner G' of G has the same vertices, and perhaps fewer edges, but the distance between a and b in G' with a factor of t of the distance between a and b in G.

Fact: DG(S) is a 2.43-spanner of G(S).

Proof omitted.

Computing Voronoi diagrams and Delaunay triangulations

From one, it's easy to get the other in O(n) time.
The Delaunay triangulation is generally easier to compute.

Some algorithms

edge flipping
- not so slow, really
- useful for triangulations besides Delaunay
divide-and-conquer
- O(n log n) worst case
- split the points in half by a line
- compute Delaunay triangulation of each piece
- merge
sweep
- O(n log n) worst case
- use sweep line paradigm in some way
randomized incremental
- O(n log n) expected
- yes, for this too.

Divide and conquer

The hard step is the merge:

Given sets L and R, S= L union R, separated by a vertical line, and their triangulations DG(L) and DG(R), compute DG(S).

There are RR edges, LL edges, and cross edges.

All new Delaunay edges in the merge are cross edges: if an RR edge is Delaunay now, it was before.

Some triangles of DG(L) are deleted: those with circumcircles containing points of R.
...and symmetrically for DG(R).

The cross edges are ordered along the split line, and consecutive edges share a vertex.

The cross edges are built from bottom to top, starting from a cross edge that is a convex hull edge.

Conceptually, maintain an empty disk on the current cross edge {a,b}, pushing it up, but keeping a and b on the bounding circle.
The first site this "rising bubble" hits gives the next cross edge.
That site is on an edge incident to a or b. To find it,
walk around the edges incident to a, deleting those bounding triangles conflicting with b. Find the first edge {a,a'} not deleted. Similarly walk around the edges incident to b, finding edge {b,b'}. Either b or b' is hit first by the "rising bubble", yielding the next cross edge {a,b'} or {a',b}.

Since O(n) work is done in the merge, O(n log n) is needed overall.
Dwyer has shown that O(n log log n) is needed for points uniform in the square, using bucketing.
Midway through a merge; edge ac is deleted

Sweep

Fortune has shown that there is a sweep line algorithm for Delaunay triangulations/Voronoi diagrams.

There are at least four viewpoints of the sweep algorithm.

When viewed as a Delaunay triangulation algorithm:
Sweep a vertical line from left to right.
For a point p on the sweep line,
the tangent empty circle at p:
is empty: has no sites inside;
is tangent to the sweep and left of it;
contains a site (on the circle, not inside)

S(p): the sites on the tangent empty circle of p.
The points on the sweep form equivalence classes based on S(p).
S(p)=a, S(p')={a,b}

Intervals of points have the same S(p) with |S(p)|=1, separated by points with |S(p)|=2, corresponding to Delaunay edges.

Sweep Events:

site events:	The sweep meets a site
circle events:	the sweep contains a point with S(p)=3; yields a Delaunay circle; corresponds to the shrinking down of an interval

Suppose the circle event has S(p)={a,b,c},
the previous interval had label {a}, bounded by points with labels {a,b} and {a,c}
Then output Delaunay edge bc;
New label {b,c} has possible events with its neighbors on the sweep: schedule these.

*A circle event, and just before*

At a site event, suppose site a is met by the sweep,
in an interval labeled by {b}, and bounded by {b,c} and {b,d}
Then a new interval, labeled by {a} and bounded by {a,b} on both sides, will be seen.
Output Delaunay edge ab.
Check for circle events for {a,b} and {b,c}, and for {a,b} and {b,d}.
Some scheduled circle events will not actually occur: for example,
if a site event splits an interval, that interval cannot shrink to a point.
if a circle event shrinks an interval to a point, that interval cannot expand to shrink a neighboring interval to a point.

A site event

Such "false alarm" events can be detected and removed when real events occur, and there are a bounded number of them.

Variations of Voronoi diagrams

Farthest-point Voronoi diagrams
- sites are points
- VR(p)==points with p farthest site
- circle emptiness of closest point --> outside of circle emptiness
- As will be shown:
  - ordinary Voronoi diagram == lower hull of points in 3d
  - farthest-point VD == upper hull of those points
The Voronoi diagram of line segments:
- the sites are line segments
- bisectors are parabolas or lines
- used for analysis of shapes
- part of this diagram is also known as the medial axis.
order k Voronoi diagrams
- sites are points
- VR(B), for B subset S == points with k nearest sites = B.

Convexity Overview

We next turn to the problem of computing convex hulls of point sets in higher dimensions; to do this, many basic facts about convex sets must be known. While actually proving these facts rigorously is a course in itself, proving some of them will motivate introducing some machinery that is useful in its own right.

Some of the facts needed are these:

Either a point p is in conv(S), or there is a subset T of S with p and S on opposite sides of the hyperplane through T.
The convex hull of S is the convex hull of the vertices of the convex hull of S.
A face of the convex hull of S is the convex hull of S', for some subset S' of S.
A set is the convex hull of a set of points iff it is the intersection of a finite set of halfspaces iff it is the intersection of all halfspaces bounded by hyperplanes through d points in S, with all points of S on one side.

As an example of the last fact: recall that the Voronoi region of a site is a polygon, if bounded: the intersection of all the halfspaces containing the site and bounded by the perpendicular bisectors of the site and other sites.

But it is also the convex hull of its vertices. This fact holds for any bounded convex polyhedral sets. They are both: the convex hulls of their vertices, and the intersection of the halfspaces bounded by their facets (edges in 2d).

The convex hull problem is to convert between these presentations: given a set S of points, find the halfspaces whose intersection is conv(S).

We'll next review basic ideas of convexity, then state and/or prove some basic facts about convex sets.

Combinations and closures

Lecture 1 gave a definition for convex combination; here we generalize it slightly, and define also three other combinations: linear, affine and positive.

Let S be a set of points, and z_p a coefficient for each point p in S. Then a point of the form

sum_{p in S} z_p p

is a combination of the points of S, of the kind shown in the table, depending on the properties of the coefficients. The table also gives the name of the corresponding sets closed under such combination, and the closure operation applied to the set.

*Combinations of point sets*
combination	conditions on z_p	closed set	closure
linear	none	linear subspace	span(S)
positive, conical	z_p >= 0	cone	cone(S)
affine	sum_{p in S} z_p = 1	affine subspace, flat	aff(S)
convex	z_p >= 0, sum_{p in S} z_p = 1	convex set	conv(S)

Note that a convex combination is affine and positive: conv(S) is the intersection of cone(S) with aff(S).

Dependence

A set S is affinely dependent if some point p in S is in aff(S-{p}), and is affinely independent if it is not affinely dependent.

For example, three collinear points are affinely dependent.

There are similar ideas of independence for the other forms:

linear independence for linear combinations,
"in convex position" for convex combinations.

dim(P) : ==dim(aff(P)) == one less than the size of the largest affinely independent set in P.

Also, dim(P) == the linear dimension of P-p, for any point p in P. That is, translate aff(P) to contain the origin; the resulting linear subspace has dimension dim(P).

Examples

(need linear independence for span(S), affine independence for the rest.)

Two distinct points S;

span(S) is a plane through the origin,

cone(S) is a planar "wedge" with apex at the origin,

aff(S) is the line through the points, and

conv(S) is the line segment between the points.

Three distinct points:

span(S) is a 3-dimensional hyperplane through the origin.

cone(S) is a cone with apex at the origin, with three sides.

aff(S) is the plane through the points, and

conv(S) is the triangle with the three points as vertices.

A few more definitions, for sets P and Q:

relint(P) : (mentioned in Lecture 5) the relative interior of P, relint(P), is the interior of P, as a subset of aff(P),

P+Q : for sets P and Q, the Minkowski sum P+Q is the set {p+q | p in P, q in Q}.

simplex : conv(S) for an affinely independent set S.

relint(P)	: (mentioned in Lecture 5) the relative interior of P, relint(P), is the interior of P, as a subset of aff(P),
P+Q	: for sets P and Q, the Minkowski sum P+Q is the set {p+q \| p in P, q in Q}.
simplex	: conv(S) for an affinely independent set S.

Intersections

Convex sets, flats, cones, and linear subspaces can all be described as intersections of halfspaces, or equivalently, as solution sets of systems of linear inequalities.

P(H):	for a set H of closed halfspaces, P(H) is the intersection of the halfspaces in H.

Since each halfspace is convex, and the intersection of convex sets is convex, P(H) is convex: a polyhedral set or convex polyhedron.

However, not all P(H) can be described as conv(S) for some finite point set S. For example, unbounded Voronoi regions cannot be.

This makes the connections between the H-polyhedra and V-polyhedra messy; the correct relation is that every P(H) can be written as conv(S)+cone(Y), for some sets S and Y.

Rather than deal with such awkwardness, it's better to move to cones.

span(S)	is a plane through the origin,
cone(S)	is a planar "wedge" with apex at the origin,
aff(S)	is the line through the points, and
conv(S)	is the line segment between the points.

span(S)	is a 3-dimensional hyperplane through the origin.
cone(S)	is a cone with apex at the origin, with three sides.
aff(S)	is the plane through the points, and
conv(S)	is the triangle with the three points as vertices.

CIS677 Notes for lecture 6: Delaunay triangulations, Convexity