a little more about ... Delaunay Triangulations: Last time we defined Delaunay triangulations and presented a number of their properties. Today we present an efficient algorithm for constructing Delaunay triangulations. Recall that we are given a set of n sites P in the plane. The Delaunay triangulation is a straight line planar graph whose internal faces are triangles such that (assuming general position) three sites p_i, p_j, and p_k are in a common triangular face if and only if their circumcircle contains no other sites.

The algorithm that we present is a randomized incremental algorithm. As with any such algorithm, the idea is to insert sites in random order, one at a time. With each insertion we update the triangulation appropriately. The issues involved with the analysis will be showing that the number of structural changes in the diagram is not very large (it will be O(1) in the expected case for each insertion). As with other incremental algorithm, we need some way of keeping track of where newly inserted sites are to be placed in the diagram. As in the trapezoidal decomposition algorithm, we could create a history graph data structure in which to store the sites. And doing so would provide us with a data structure for doing point locations in the final triangulation. Instead, we will opt for a simpler bucketing method, which will be described later. In this case, we will need to argue that the number of times that a site is rebucketed on average is not too large (it will be O(logn) in the expected case). The result will be an O(n log n) algorithm.

Incremental update: When we add the next site, p_i , the problem is to convert the current Delaunay triangulation into a new Delaunay triangulation containing this site. This will be done by creating a non Delaunay triangulation containing the new site, and then incrementally ``fixing'' this triangulation to restore the Delaunay properties. The fundamental changes will be: (1) adding a site to the middle of a triangle, and creating three new edges, and (2) performing an edge flip. Both of these operations can be performed in O(1) time, assuming that the triangulation is maintained, say, as a DCEL.

Figure 63: Basic triangulation changes.

The algorithm that we will describe has been known for many years, but was first analyzed by Guibas, Knuth, and Sharir. We start with an initial triangulation. Guibas, Knuth and Sharir suggest starting with the triangle of three sites ``at infinity''. (Not just any enclosing triangle will work, since we need to be sure that the newly added vertices will not affect the structure of the circumcircles, and hence the interior structure of the triangulation). This guarantees that all sites to be added, will lie within some triangle of the existing triangulation. The sites are added in random order. When a new site p is added, we find the triangle of the current triangulation that contains this site (we will see how later), insert the site in this triangle, and join this site to the three surrounding vertices. This creates three new triangles, pab, pbc, and pca, each of which may or may not satisfy the empty circle condition. How do we test this? For each of the triangles that have been added, we check the vertex of the triangle that lies on the opposite side of the edge that does not include p. (If there is no such vertex, because this edge is on the convex hull, then we are done.) If this vertex fails the incircle test, then we swap the edge (creating two new triangles that are adjacent to p). This replaces one triangle that was incident to p with two new triangles. We repeat the same test with these triangles. An example is shown in the figure below.

The following is a description of the algorithm (Guibas, Knuth, and Sharir give a nonrecursive version). The current triangulation is kept in a global data structure. The edges in the following algorithm are actually pointers to the DCEL.

Insert(p):
  Find the triangle  containing p;
  Insert edges pa, pb, and pc into triangulation:
    SwapTest(ab);                             // Fix the surrounding edges
    SwapTest(bc);
    SwapTest(ca);

SwapTest(ab):
  if (ab is an edge on the exterior face) return
  Let d be the vertex to the right of edge ab;
    if (inCirc(p, a, b, d) )                // d violates the incircle test
       Flip edge ab for pd;
       SwaptTest(ad); // Fix the new suspect edges
       SwaptTest(db);

Figure 64: Point insertion.

There is only one major issue in establishing the correctness of the algorithm. When we performed empty circle tests, we only tested (1) triangles containing the site p, and (2) only sites that lay on the opposite side of an edge of such a triangle. To establish (1), observe that it suffices to consider only triangles containing p because since p is the only newly added site, it is the only site that can cause a violation of the empty circle property. To establish (2) we argue that if for every site d, which is opposite from p along some edge ab, lies outside the circumcircle of pab, then all these circumcircles are empty. A complete proof takes some effort, but here is a simple justification. What could go wrong? It might be that d lies outside the circumcircle, but there is some other site (e.g. a vertex e of a triangle adjacent to d that lies inside the circumcircle). This is illustrated in the following figure. We claim that this cannot happen. It can be shown that if e lies within the circumcircle of pab, then a must lie within the circumcircle of bde. (The argument is a exercise in geometry.) However, this violates the assumption that the initial triangulation (before the insertion of p) was a Delaunay triangulation.

Figure 65: Proof of sufficiency of testing neighboring sites.

As you can see, the algorithm is very simple. The only things that need to be implemented are the DCEL (or other data structure) to store the triangulation, the incircle test, and locating the triangle that contains p. The first two tasks are straightforward. The point location involves a little thought.

Point Location: The point location can be accomplished by one of two means. One approach is to maintain a history graph point location data structure, just as we did in the trapezoid decomposition case. A simpler approach is based on the idea of maintaining the uninserted sites in a set of buckets. Think of each triangle of the current triangulation as a bucket that holds the sites that lie within this triangle and have yet to be inserted. Whenever an edge is flipped, or when a triangle is split into three triangles through point insertion, some old triangles are destroyed and are replaced by a constant number of new triangles. When this happens, we lump together all the sites in the buckets corresponding to the deleted triangles, create new buckets for the newly created triangles, and reassign each site into its new bucket. Since there are a constant number of triangles created, this process requires O(1) time per site that is rebucketed.

Analysis: To analyze the expected running time of algorithm we need to bound two things:

(1) how many changes are made in the triangulation on average with the addition of each new site, and (2) how much effort is spent in rebucketing sites. As usual, our analysis will be in the worst case (for any point set) but averaged over all possible insertion orders. We argue first that the expected number of edge changes with each insertion is O(1) by a simple application of backwards analysis. First observe that (assuming general position) the structure of the Delaunay triangulation is independent of the insertion order of the sites so far. Thus, any of the existing sites is equally likely to have been the last site to be added to the structure. Suppose that some site p was the last to have been added. How much work was needed to insert p? Observe that the initial insertion of p involved the creation of three new edges, all incident to p. Also, whenever an edge swap is performed, a new edge is added to p. These are the only changes that the insertion algorithm can make. Therefore the total number of changes made in the triangulation for the insertion of p is proportional to the degree of p after the insertion is complete.

Thus, by a backwards analysis, the expected time to insert the last point is equal to the average degree of a vertex in the triangulation. (The only exception are the three initial vertices at infinity, which must be the first sites to be inserted.) However, from Euler's formula, we know that the average degree of a vertex in any planar graph is at most 6. (To see this, recall that a planar graph can have at most 3n edges, and the sum of vertex degrees is equal to twice the number of edges, which is at most 6n.) Thus, (irrespective of which insertion this was) the expected number of edge changes for each insertion is O(1).

Next we argue that the expected number of times that a site is rebucketed (as to which triangle p it lies in) is O(logn). Again this is a standard application of backwards analysis. Consider the i th stage of the algorithm (after i sites have been inserted into the triangulation). Consider any one of the remaining n - i sites. We claim that the probability that this site changes triangles is at most 3/i (under the assumption that any of the i points could have been the last to be added).

To see this, let q be an uninserted site and let D be the triangle containing q after the ith insertion. As observed above, after we insert the ith site all of the newly created triangles are incident to this site. Since T is incident to exactly three sites, if any of these three was added last, then T would have come into existence after this insertion, implying that q required (at least one) rebucketing. On the other hand, if none of these three was the last to have been added, then the last insertion could not have caused q to be rebucketed. Thus, (ignoring the three initial sites at infinity) the probability that q required rebucketing after the last insertion is exactly 3/i. Thus, the total number of points that required rebucketings as part of the last insertion is (n-i)3/i. To get the total expected number of rebucketings, we sum over all stages, giving

 sum(i=1..n) (n-i)3/i <= sum(i=1..n) n(3/i)  = 3n sum(i=1..m) 1/i
=3n ln n.

Thus, the total expected time spent in rebucketing is O(n log n), as desired. There is one place in the proof that we were sloppy. (Can you spot it?) We showed that the number of points that required rebucketing is O(n log n), but notice that when a point is inserted, many rebucketing operations may be needed (one for the initial insertion and one for each each additional edge flip). We will not give a careful analysis of the total number of individual rebucketing operations per point, but it is not hard to show that the expected total number of individual rebucketing operations will not be larger by more than a constant factor. The reason is that (as argued above) each new insertion only results in a constant number of edge flips, and hence the number of individual rebucketings per insertion is also a constant. But a careful proof should consider this.