Delaunay Triangulation

Motivation

Let us consider a function $f:\mathbb R^2\to\mathbb R$. Its values are known only in the final set of points $P$. One way to approximate such a function is to divide the convex hull of the set $P$ into triangles with vertices at the points of the set $P$ and replace the function $f$ by a function which is linear on these triangles. This means the following: Each point $x$ of a triangle $p_i, p_j, p_k$ is a convex combination of the vertices

and hence we put

The graph of the approximation of the function $f$ is thus formed by triangles in $\mathbb R ^ 2\times\mathbb R$ with vertices $(p_i,f(p_i)$, $(p_j,f(p_j)$ and $(p_k,f(p_k))$.

Figure 10.1
Linear approximation of the function $f$.

This approximation depends substantially on how we do the triangulation of the convex hull of our set. It seems that a good approximation should be one that has no triangles with small angles. We will deal with such triangulations in the following text.

Angle-optimal triangulation

First, we show that for a given set of points $P$ in a plane, all the triangulations of its convex hull have the same number of triangles. This allows us to compare different triangulations lexicographicly according to the sizes of angles.

Let $\mathcal T$ be a triangulation of the convex hull of a set $P$ with $m=2n-2-k$ triangles. The triangles of this triangulation have $3m$ angles. Let us order them into a sequence

On these sequences we define the lexicographic order

if there is $i\in\{1,2,\dots,3m\}$ such that

A triangulation $\mathcal T$ is called angle-optimal if it is maximal in this ordering.

Legal edges and legal triangulations

For the purposes of our algorithm which should create a suitable triangulation we will need the notion of legal edge and legal triangulation. Before we start with a definition we repeat shortly high school geometry, especially what we know about the angles on the arc of a circle.

Consider points $p$ and $q$ on the circle $\mathcal K$ with a center $s$. The points divide the circle into two arcs. Choose a point $a$ on the first of them and a point $b$ on the second one. Denote the angles in the triangles $paq$ and $pbq$ at the vertices $a$ and $b$ as $\alpha$ and $\beta$, respectively. We call them arc angles of the chord $pq$. Arc angles have these two important properties:

• The size of angles $\alpha$ and $\beta$ does not depend on the position of points $a$ and $b$ on the given arcs.
• Their sum is $\alpha+\beta=180^\circ$.

From here we can easily derive that every triangle $pcq$ where the point $c$ lies in the halfplane $paq$ inside the circle $\mathcal K$ has the angle

at the vertex $c.$ Similarly every triangle $pdq$ with the vertex $d$ lying in the halfplane $paq$ outside the circle $\mathcal K$ has the angle

at the vertex $d$.

Figure 10.2
Arc angles

From the properties of arc angles it can be shown that we can circumscribe a circle to a quadrilateral if and only if the sum of opposite angles is equal to $180^\circ$.

Now consider a triangulation $\mathcal T$ and an edge $pq$ of $\mathcal T$ which has two adjacent triangles $pqr$ a $pqt$. Such an edge is called illegal edge of the triangulation $\mathcal T$ if the point $t$ lies inside the circle circumscribed to the triangle $pqr$. It is equivalent to the fact that the sum of angles at the vertices $r$ and $t$ in the quadrilateral $prqt$ is greater than $180^\circ$ and also to the fact that the point $r$ lies inside the circle circumscribed to the triangle $pqt$.

Figure 10.3
Illegal edge $pq$.

The edges of the triangulation $\mathcal T$ which are not illegal are called legal and the triangulation without illegal edges is called legal triangulation. Immediate consequence of this definition is the following characterization of legal triangulations:

The following assertion will play a crucial role in our algorithm.

Figure 10.4
Flip

Proof

Since the edge $pq$ is illegal in the triangulation $\mathcal T$, the sum of the angles at vertices $r$ and $t$ in the quadrilateral $prqt$ is bigger than $180^\circ$ and so the sum of the angles at $p$ and $q$ in the same quadrilateral is less than $180^\circ$. (The sum of all angles in the quadrilateral is $360^\circ$.) That is why the edge $rt$ in the triangulation $\mathcal T'$ is legal.

Denote the angles before and after the flip according to the following figure:

Figure 10.5
Notation.

Since all the other triangles of both triangulations are the same, for

$$\alpha(\mathcal T)\lt\alpha(\mathcal T'),$$

it suffices to prove that

$$\min \{\alpha_i,\ i=1,2,\dots,6\}\lt\min \{\beta_i,\ i=1,2,\dots,6\}.$$

We show subsequently that for every $i\in\{1,\dots ,6\}$ there exists $j\in\{1,\dots ,6\}$ such that $\beta_i \gt \alpha_j.$ Obviously,

$$\beta _1\gt \alpha_1\quad\text{ a }\quad \beta_4\gt \alpha_5.$$

Further, $\beta _2$ is an arc angle over the chord $pr$ in the circle circumscribed to the triangle $prt$. Since the point $q$ lies outside this circle, we get

$$\beta_2\gt\alpha_5.$$

Similarly $\beta _6$ is an arc angle over the chord $pt$ in the circle circumscribed to the triangle $prt$. Hence

$$\beta_6\gt\alpha_4.$$

Analogously, using the circle circumscribed to the triangle $qrt$ we can prove that

$$\beta_3\gt \alpha_1\quad \text{ and }\quad \beta_5\gt \alpha_2.$$

\(\Box\)

Legal triangulations need not to be angle-optimal as the following figure shows on the example of two triangulations of the pentagon $p_1p_2p_3p_4p_5$ to which we can circumscribe a circle. Moreover, $p_1p_2p_3p_5$ is a square and $p_1p_2p_4$ is a isosceles triangle. Both triangulations $\mathcal T$ and $\mathcal T'$ are legal but because $\alpha(\mathcal T)\lt \alpha(\mathcal T')$, the triangulation $\mathcal T$ is not angle-optimal.

Figure 10.6
The example of legal triangulation which is not angle-optimal.

Delaunay graph and triangulation

Besides angle-optimal triangulation and legal triangulation we introduce another notion of Delaunay triangulation. First we define what Delaunay graph is.

The Delaunay graph of the set $P$ has $P$ as the set of vertices. Two vertices $p_i$ and $p_j$ are connected by an edge if there is a circle on which these two points are lying and all the other points of the set $P$ lie outside the circle. It means, that the center of this circle lies on the edge of the Voronoi diagram which separates the cell determined by $p_i$ from the cell determined by $p_j$. That is why Delaunay graph is dual to the Voronoi diagram.

Figure 10.7
The duality between the Voronoi diagram (dashed) and the Delaunay graph (red).

Figure 10.8
Illustration of the proof of Lemma 10.5

Bounded faces of the Delaunay graph are polygons. Let points $p_1,p_2,\dots,$ $p_k$ form the vertices of one of these polygons. This polygon determines a vertex $v$ in the dual graph, i.e. in the Voronoi diagram. Hence the points $p_1,p_2,\dots,p_k$ have the same distance to the point $v$ and lie on one circle with the center $v$.

Figure 10.9
The quadrilateral $p_1p_2p_3p_4$ in the Delaunay graph (black) and the Voronoi diagram (red).

Delaunay triangulation is defined as an arbitrary triangulation of these polygons in the Delaunay graph. The Delaunay triangulation need not be determined uniquely. Nevertheless, the centers of the circles circumscribed to the triangles in the Delaunay triangulation form the set of the vertices of the Voronoi diagram.

Figure 10.10
Delaunay graph (black) and Delaunay triangulation (black and red).

If no four points of the set $P$ lie on one circle then the Delaunay graph is already the Delaunay triangulation. In this case the Delaunay triangulation is determined uniquely. We say that in this case the points of the set $P$ are in general position.

Every Delaunay triangulation can be characterized by the following property.

Using this lemma and the characterization of legal triangulations by Lemma Lemma 10.2, we get immediately that every Delaunay triangulation is legal. To prove that every legal triangulation is Delaunay is a little bit more complicated.

Figure 10.11
To the proof of Theorem 10.7

From here and from Corollary 10.4 we get:

Basic idea of randomized incremental algorithm

One way how to construct Delaunay triangulation is to use algorithm for the Voronoi diagram, then to find the dual graph to it and to divide its polygons into triangles.

We will present a different approach based on the fact that every legal triangulation is Delaunay and which is based on gradual removing of illegal edges.

A naive version of such an algorithm is the following one: Make a triangulation of the convex hull of the set $P$ and go through its edges. When you find an illegal edge, replace it with a legal edge by a flip and proceed by the same way further. This process is finite because after every flip we get a triangulation which is bigger than the previous one in described lexicographic order. The number of triangulations is finite, so this process must sometimes end.

A sophisticated version of this approach is the following randomized incremental algorithm. Its basic idea can be captured in several points:

1. We denote $p_0$ the maximal point from the set $P$ (in lexicographic order firstly with respect to $y$-coordinate and then to $x$-coordinate). We add points $\pj$ and $\pd$ to the points of the set $P$ so that all other points of $P$ lie inside the triangle $p_0\pj\pd$.
2. We order the remaining points randomly into a sequence $p_1,p_2,\dots,p_n$. Gradually for $r=1,2,\dots,n$ we add the point $p_r$ to a legal triangulation $\mathcal T_{r-1}$ of the triangle $p_0\pj\pd$ where the vertices of single triangles are in the set $$p_{-2},p_{-1},p_0,p_1,\dots,p_{r-1}.$$
3. Using a simultaneously constructed search structure we find in which triangle $p_ip_jp_k$ of the triangulation $\mathcal T_{r-1}$ the point $p_r$ lies and if it lies in the interior or on the edge. Then we connect the point $p_r$ by edges with vertices of this triangle (if it lies in the interior) or with vertices of adjacent triangles (if it lies on an edge). This newly created triangulation of the triangle $p_0\pj\pd$ with vertices in the set $$\{p_{-2},p_{-1},p_0,p_1,\dots,p_{r-1},p_r\}$$ can be "legalized" using flips which exchange illegal edges for legal ones. After removing all illegal edges we get a triangulation $\mathcal T_r$.
4. At the end we remove the points $\pj$ and $\pd$ and the edges going from these points out of the legal triangulation $\mathcal T_n$. Due to a suitable choice of the points $\pj$ and $\pd$ the remaining triangulation is the legal and also Delaunay triangulation of the convex hull of the set $P=\{p_0,p_1,p_2,\dots,p_n\}$.

Figure 10.12
The set $P$ in the triangle $p_0\pj\pd$.

Now we describe aforementioned steps in more details. We start with (2) and (3).1

Legalization

Let us consider a legal triangulation $\mathcal T_{r-1}$ described in the step (2) above. Using a search structure we can find a triangle $p_ip_jp_k$ in whose interior lies the point $p_r$. Then we connect point $p_r$ with the vertices $p_i,p_j,p_k$ by edges.

Figure 10.13
The point $p_r$ in the interior of a triangle $p_ip_jp_k$.

This will create a new triangulation. It is important that the newly added edges based on $p_r$ are legal. This can be proved as follows: Let $\mathcal K$ be a circle circumscribed to the triangle $p_ip_jp_k$. Since the triangulation $\mathcal T_{r-1}$ is legal and hence also Delaunay, there is no point from the set $P_{r-1}=\{\pd,\pj, p_0,\dots,p_{r-1}\}$ in the interior of $\mathcal K$. Let $\mathcal L$ be a circle homothetic to $\mathcal K$ which has the cord $p_rp_i$. See Figure 10.14. In the interior of this circle there are no points from $P_{r-1}$. So the edge $p_rp_i$ is an edge of the Delaunay graph for the set $P_r$ and this has to be legal. Analogously, one can prove that also the edges $p_rp_j$ and $p_rp_k$ are legal.

Figure 10.14
To the proof that $p_rp_i$ is legal.

It can happen that some of the edges $p_ip_j$, $p_ip_k$ a $p_jp_k$ become illegal in the new triangulation. Using a flip we can replace them by edges going from the point $p_r$. For example let $p_ip_j$ be illegal. Then it is the edge not only of the triangle $p_ip_jp_r$ but also of an adjacent triangle $p_ip_jp_l$. We replace it by the edge $p_rp_l$ which is already legal. Now some of the edges $p_ip_l$ and $p_jp_l$ can become illegal. That is why we have to replace it by a flip with an edge based in the point $p_r$.

Figure 10.15
The edge $p_ip_j$ is illegal.

We proceed in this way as long as there is no illegal edge. (After every flip the new triangulation is bigger in the lexicographic order than the previous one.) The procedure is captured in the following animation and pseudocode. The input is the point $p_r$, the edge $p_ip_j$ and a legal triangulation $\mathcal T_{r-1}$.

If we find that the point $p_r$ lies on an edge $p_ip_j$, then we replace the edge $p_ip_j$ by the edges $p_rp_i$ and $p_rp_j$ and we create new edges connecting the point $p_r$ with the vertices $p_k$ and $p_l$ of adjacent triangles to the edge $p_ip_j$. Similarly as above one can show that the new edges $p_rp_i$, $p_rp_j$, $p_rp_k$ a $p_rp_l$ are legal. However the edges $p_ip_k$, $p_jp_k$, $p_ip_l$ a $p_jp_l$ can become illegal. This can be corrected by making gradually flips as in the above case.

Figure 10.16
The point $p_r$ on the edge $p_ip_j$.

Search structure

In the course of the algorithm we have a relevant search structure for every triangulation. This is an oriented graph the leaves of which are the triangles of the actual triangulation and inner nodes are the triangles of the previous triangulations.

At the beginning the triangulation $\mathcal T_0$ has the only triangle $p_0\pj\pd$ and the corresponding search structure $\mathcal D_0$ has only one node given by this triangle. For the triangulation $\mathcal T_{r-1}$ we have a search structure $\mathcal D_{r-1}$. We change the triangulation $\mathcal T_{r-1}$ gradually in the way described in the previous section. Simultaneously we change the search structure. To every transition in the triangulation there is a corresponding transition from one search structure to further search structure. This is illustrated in the Figure 10.17.

Figure 10.17
The change of the search structure corresponding to newly added edges inside the triangle $p_ip_jp_k$.

Figure 10.18
The change of the search structure corresponding to a flip.

Using this procedure we create a search structure $\mathcal D_r$ corresponding to the legal triangulation $\mathcal T_r$. The resulting search structure depends only on the random order of the points $p_1,p_2,\dots,p_r$.

Initial choice of points $p_0$, $p_{-1}$, $p_{-2}$

Let us now explain the steps (1) and (4) of our algorithm. We will refer to the geometric idea for better understanding. Important is that the points $\pj$ a $\pd$ will be determined by their properties, we will not compute their concrete coordinates. In the decisions whether some edge is legal or illegal four points are considered. If none of these points is $\pj$ or $\pd$, we have to find the center of the circle circumcsribed to three of the points and decide if the fourth point lies in the interior of the circle. If some of these points is $\pj$ or $\pd$, we will proceed according to special rules described below which follow from the way how the points $\pj$ a $\pd$ are chosen.

As we have already said for the point $p_0$ we choose that one from $P$ which is maximal in the lexicographic order first with respect to the coordinate $y$ and then with respect to the coordinate $x$.

We choose the point $\pj$ to the right of the set $P$ and under it. (More exactly: its $y$-coordinate is smaller than the minimum of $y$-coordinates of the points from $P$ and its $x$-coordinate is bigger than the maximum of $x$-coordinates of the points from $P$.) The choice of $\pj$ is carried out so that

1. it lies outside all circles circumscribed to all triples of points from $P$ which do not lie in one line,
2. all points $p_i$ with $i\gt 0$ lie under the line $p_0\pj$.

We choose the point $\pd$ to the left of the set $P$ and over it so that

1. it lies outside all circles circumscribed to all triples of points from $P\cup\{\pj\}$ which do not lie in one line,
2. all points $p_i$ with $i\gt 0$ lie under the line $p_0\pd$,
3. all points $p_i$ with $i\ge 0$ lie over the line $\pj\pd$.

So all the points $p_i$ with $i\gt 0$ lie in the interior of the triangle $p_0\pj\pd$ and every quadrilateral $\pd p_i \pj p_j$ is not convex which implies that the point $\pj$ lies outside all circles circumscribed to the triangles $\pd p_i p_j$.

Figure 10.19
The set $P$ in the triangle $p_0\pj\pd$.

Under these assumptions every Delaunay (legal) triangulation of $P\cup\{\pj,\pd\}$ consists of

1. the edges going from the point $\pj$ to the vertices of the right hand side of the convex hull of the set $P$,
2. the edges going from the point $\pd$ to the vertices of the left hand side of the convex hull of the set $P$,
3. a Delaunay triangulation of the set $P$.

Consequently we can obtain a Delaunay triangulation of the set $P$ in the way described in the step (4).

The rules of how to work with our algorithm are based on the following lemma:

Instead of doing a proof we demonstrate the situation by a picture.

Figure 10.20
To Lemma 10.9. The points $\pj$ and $\pd$ are chosen in such a way that in the coloured angles there is no other point of the set $P$.

When finding whether an edge is legal or illegal we use the following rules:

• All edges of the triangle $p_0\pj\pd$ are legal.
• The edge $\pd p_j$ with adjacent vertices $\pj$ and $p_k$ is legal.

Figure 10.21
Illustration of the previous statement.

• The edge $\pj p_j$ with adjacent vertices $\pd$ and $p_k$ is legal.
• Let the edge $p_i p_j$ have adjacent vertices $p_k$ and $p_l$, let the quadrilateral $p_ip_kp_jp_l$ be convex and let there be just one negative number among the numbers $i,j,k,l$. Then the edge $p_ip_j$ is legal if and only if $$\min(k.l)\lt\min(i,j).$$

Figure 10.22
$\min(-2,l)\lt \min(i,j)$, the edge $p_ip_j$ is legal.

Figure 10.23
$\min(k,l)\gt \min(-1,j)$, the edge $\pj p_j$ is illegal.

Pseudocode of the algorithm and its expected time complexity

A brief pseudocode of the algorithm described in details in the previous sections can look like this:

Without proof we give the following statement about expected time complexity of the algorithm. See [Berg et al] for the proof.