[Cyan]<PreISLCedar>Stolfi>Papers>SubdCY.tex!1

If we are given only two sites, then the associated Voronoi regions are simply the two (open) half-planes delimited by the bisector of the two sites. More generally, when $n$ sites are given, the region associated with a particular site $p$ will be the intersection of all half-planes containing $p$ and delimited by the bisectors between $p$ and the other sites. It follows that the Voronoi regions are (possibly unbounded) convex polygons, whose edges are portions of inter-site bisectors, and whose vertices are circumcenters of triangles defined by three of the sites. An example Voronoi diagram for a small collection of sites is shown in fig. 7.1.\par

As mentioned above, most of the time we will actually deal with a dual of the Voronoi subdivision, commonly called the Delaunay triangulation. As the name implies, this is a planar subdivision whose vertices are the $n$ sites and whose edges are straight line segments that connect every pair of sites associated with regions that share a common edge. If no four of the sites happen to be cocircular, then this dual is in fact a triangulation and, in any case, additional edges can be introduced to make it one. Although our algorithms can handle the degeneracies implied by four cocircular points, in this section only we will assume that such degeneracies do not arise, so we can speak unambiguously of the Delaunay triangulation of the $n$ sites. If $S$ denotes our collection of sites, then $D(S)$ will denote the Delaunay triangulation of $S$.\par

We say that a circle is {\it point-free} if none of the given sites is contained in its interior. It is known that the circumcircle of each Delaunay triangle is point-free. Actually the converse is also true. If the circumcircle of a triangle defined by three of the sites is point-free, then that triangle is a Delaunay triangle. And if it is possible to find a point-free circle passing through two sites, then these sites will be connected by a Delaunay edge. For a discussion of these facts see [\Lee]. We will speak of these circles as {\it witnesses} to the ‘‘Delaunayhood’’ of the corresponding triangle or edge respectively. For completeness we also state here a couple of other obvious properties of the Delaunay triangulation which will be useful in the sequel.\par

In particular this implies that $D$ should be inside the circle $ABC$ if the points $A$, $B$, and $C$ define a counterclockwise oriented triangle, and {\it outside}, if they define a clockwise oriented one. (In case $A$, $B$, and $C$ are collinear we interpret the line as a circle by adding a point at infinity). If $A$, $B$, $C$, and $D$ are cocircular, then our predicate returns false. Notice that the test is equivalent to asking whether $\angle ABC + \angle CDA >\angle BCD + \angle DAB$. Another equivalent form of it is given below, based on the coordinates of the points.\par

\lemma{8.1}{The test $\InCircle(A,B,C,D)$ is equivalent to
$$\Dscr(A,B,C,D) = \det{\vcenter{
\halign{$\ctr{#}$\quad–$\ctr{#}$\quad–
$\ctr{#}$\quad–$\ctr{#}$\cr
xǍ–yǍ–xǍ↑2+yǍ↑2–1\cr
\noalign{\vskip 2.5pt}
xB̌–yB̌–xB̌↑2+yB̌↑2–1\cr
\noalign{\vskip 2.5pt}
xČ–yČ–xČ↑2+yČ↑2–1\cr
\noalign{\vskip 2.5pt}
xĎ–yĎ–xĎ↑2+yĎ↑2–1\cr}}} < 0.$$}\par

\proof{We consider the following mapping from points in the plane to points in space:
$$\lambda:\ (x,y)\mapsto(x,y,x↑2+y↑2),$$
which lifts each point on the $x,y$-plane onto the paraboloid of revolution $z=x↑2+y↑2$. See fig. 8.2 for an illustration. We first show that $A$, $B$, $C$, and $D$ are cocircular if and only if $\lambda(A)$, $\lambda(B)$, $\lambda(C)$, and $\lambda(D)$ are coplanar, a rather amazing fact.\par

Suppose first that $A$, $B$, $C$, and $D$ are cocircular. If we have the degenerate case where they are collinear, then $\Dscr(A,B,C,D)$ is zero, as we can see by expanding it by the third column. But $\Dscr(A,B,C,D)$ is also the (signed) volume of the tetrahedron defined by $\lambda(A)$, $\lambda(B)$, $\lambda(C)$, and $\lambda(D)$. Since the volume is zero, the points must be coplanar. Otherwise let $(p,q)$ denote the center and $r$ the radius of the circle passing through the points $A$, $B$, $C$, $D$. We must have
$$(xǍ-p)↑2+(yǍ-q)↑2=r↑2,$$
or equivalently
$$-2p\cdot xǍ-2q\cdot yǍ+1\cdot (xǍ↑2+yǍ↑2)+(p↑2+q↑2-r↑2)\cdot 1=0.\eqno(1)$$
This relation also holds for points $B$, $C$, and $D$, and therefore we have a linear dependence among the columns of the determinant $\Dscr(A,B,C,D)$, which implies that its value is zero. So again we can conclude that $\lambda(A)$, $\lambda(B)$, $\lambda(C)$, and $\lambda(D)$ are coplanar.\par

Now conversely suppose that $\lambda(A)$, $\lambda(B)$, $\lambda(C)$, and $\lambda(D)$ are coplanar. If all of $A$, $B$, $C$, and $D$ are collinear, then we are done. So suppose, without loss of generality, that $A$, $B$, and $C$ are not collinear. As above, let $(p,q)$ denote the center and $r$ the radius of the circumcircle of triangle $ABC$. Then $A$, $B$, and $C$ satisfy equation (1) above. Since $A$, $B$, and $C$ are {\it not} collinear, the corresponding three rows of $\Dscr(A,B,C,D)$ are linearly independent. But all four rows are linearly dependent, since the determinant is zero. So the last row can be expressed as a linear combination of the first three, and therefore point $D$ satisfies (1) as well , i.e., it is on the circle $ABC$.\par

The above result shows that planar sections of the paraboloid of revolution $z=x↑2+y↑2$ project onto circles in the $x,y$-plane. The paraboloid is a surface that is convex upwards, and therefore, in a section of it with a plane, the part below the plane projects to the interior of the corresponding circle in the $x,y$-plane, and the part above the plane to the exterior. From this and the standard right-handed orientation convention for the sign of volumes the lemma follows. Notice that this establishes an interesting correspondance between circular queries in the plane and half-space queries in 3-space.}\par

As a side note we remark that Ptolemy’s theorem in Euclidean plane geometry does not lead to a useful implementation of the {\InCircle} test, as we always have
$$AB\times CD+BC\times AD \geq BD\times AC,$$
with equality only when the four points are cocircular. In fact, the quantity one obtains by rendering $AB\times CD+BC\times AD - BD\times AC$ radical-free is essentially the square of the determinant $\Dscr(A,B,C,D)$ above.\par

\lemma{8.2}{If $A$, $B$, $C$, $D$ are any four non-cocircular points in the plane, then transposing any adjacent pair in the predicate $\InCircle(A,B,C,D)$ will change the value of the predicate from true to false, or vice versa. In particular, the boolean sequence $\InCircle(A,B,C,D)$, $\InCircle(B,C,D,A)$, $\InCircle(C,D,A,B)$, $\InCircle(D,A,B,C)$ is either T(rue), F(alse), T, F, or F, T, F, T.}\par

The last two lemmas show that in $\InCircle(A,B,C,D)$ all four points play a symmetric role, even though from the definition $D$ seems to be special. If $A$, $B$, $C$, $D$ define a counterclockwise oriented quadrilateral and $\InCircle(A,B,C,D)$ is true, then $\InCircle(C,B,A,D)$ is false, so reversing the orientation flips the value of the predicate. Whenever we consider applying this test to a (rooted) quadrilateral in the sequel, we will take the quadrilateral to be counterclockwise oriented.\par

For a convex quadrilateral $ABCD$ there are two ways to triangulate it. We can either draw diagonal $AC$ or diagonal $BD$. By lemma 7.2 the sides $AB$, $BC$, $CD$, and $DA$ are always Delaunay edges. If $\InCircle(A,B,C,D)$ is false, then circle $ABC$ is point-free, so the diagonal $AC$ completes the Delaunay triangulation for the four points $A$, $B$, $C$, $D$ (else the diagonal $BD$ does). More generally, for any collection of points, all the convex hull edges are Delaunay edges. For each non-hull Delaunay edge $BD$, if the two triangles $ABD$ and $BDC$ on its two sides form a convex quadrilateral, then $\InCircle(A,B,C,D)$ is true. A converse of this also holds [\Lee]. This is commonly expressed in terms of the ‘‘swapping rule’’, an implementation of which was presented in section 5. This rule states that if there exists a convex quadrilateral $ABCD$ that has been triangulated with diagonal $AC$, but $\InCircle(A,B,C,D)$ is true, then this diagonal should be erased and diagonal $BD$ should be drawn. Success of the swapping rule on an edge ($AC$ in the previous sentence) is a witness for non-Delaunayhood. The above converse then states that a triangulation stationary under the swapping rule is the Delaunay triangulation. In fact, the swapping rule provides a new necessary and sufficient condition for the Delaunayhood of an edge, namely that the rule should fail for any quadrilateral defined by the endpoints of the edge and any two other points, one on each side.\par

A useful intuition about the {\InCircle} test comes from the following observation. Given two points in the plane $X$ and $Y$, the set of circles passing through $X$ and $Y$ forms a one parameter family. For example, we can think of these circles as the family $\left\{Cť\right\}$, where $t$ denotes the signed distance of the center of such a circle along the bisector of $XY$, measured from the midpoint of $XY$. Thus $C{̌-}̄$ denotes the circle corresponding to the half-plane to the left of the line $XY$, $C0̌$ denoted the circle with diameter $XY$, and $Č$ denotes the circle corresponding to right half-plane of $XY$. Note that the portion of these circles to the left of $XY$ strictly decreases (by proper inclusion) as $t$ increases, while the portion to the right of $XY$ strictly increases. See fig. 8.3. The next lemma uses this intuition to derive another property of the Delaunay triangulation that will be useful to us in the next section.\par

\lemma{8.3}{Let $A$ be a site and consider any line $\lscr$ through $A$. For convenience of terminology, assume that $\lscr$ is horizontal. Let $N1̌, N2̌, \ldotss, Nǩ$ denote the Delaunay neighbors of $A$ lying above $\lscr$, listed in counterclockwise order. If $X$ is any point of $\lscr$ to the right of $A$, let $\Gammaǐ$ denote the portion of the circumcircle of $AXNǐ$ that is above $\lscr$. Then the sequence $\left\{\Gammaǐ\right\}$ is unimodal, in the sense that there is some $j$, $1 j k$, such that for $1 i<j$ we have $\Gammaǐ\supset\Gamma{̌i+1}$, while for $j i<k$ we have $\Gammaǐ\subset\Gamma{̌i+1}$.}\par

Consider, as in fig. 8.4, three consecutive Delaunay neighbors $Nǐ$, $N{̌i+1}$, and $N{̌i+2}$ of $A$. The point $N{̌i+2}$ is outside the circle $ANǐN{̌i+1}$ and points $Nǐ$ and $N{̌i+2}$ are on opposite sides of $AN{̌i+1}$. So to get to the circle $AN{̌i+1}N{̌i+2}$ from $ANǐN{̌i+1}$, while always passing through $A$ and $N{̌i+1}$, we must expand on the side of $N{̌i+2}$, and therefore we must contract on the side of $Nǐ$, where $Xǐ$ and $X{̌i+1}$ also lie. This proves our assertion.\par

To prove the lemma now, note that as long as the $Xǐ$ are to the right of $X$, then $X$ is inside the circumcircle of $ANǐN{̌i+1}$, or equivalently, $N{̌i+1}$ is inside the circumcircle of $ANǐX$. After the $Xǐ$ move to the left of $X$, then $N{̌i+1}$ is outside the circumcircle of $ANǐX$. Thus the $\Gammaǐ$ behave as stated.}\par

\proof{Without loss of generality let us assume that in the convex angle $\angle XAY$ the ray $AY$ is counterclockwise of $AX$. Suppose the lemma is false inside the convex angle $\angle XAY$. Let $N$ be the Delaunay neighbor of $A$ outside the circle $AXY$ which is the first to be encounterd as we sweep this angle in counterclockwise fashion. Let also $M$ denote the predecessor of $N$, as in fig. 8.5. Note that $N$ is also outside the circle $AMY$ since $X$, $M$, and $N$ are all on the same side of $AY$. Therefore $Y$ is inside the circle $AMN$, contradicting Delaunayhood. A similar argument proves the lemma outside $\angle XAY$.}\par

In this section we use the tools we have developed so far to describe, analyze, and prove correct a divide and conquer algorithm for computing the Delaunay triangulation of $n$ points in the plane. Topologically the quad-edge data structure gives us the dual for free, so by associating some relevant geometric information with our face nodes, e.g., the coordinates of the corresponding Voronoi vertices, we are simultaneously computing the Voronoi diagram as well. Our algorithm is the dual of that proposed by Shamos and Hoey [\ShamHo] and, like theirs, runs in time $O(n \lg n)$ and uses linear storage. In [\Lee] D. T. Lee proposes a dual algorithm which is similar to but more complex than ours. Unlike both of these, our algorithm need not compute straight line intersections unless the coordinates of Voronoi vertices are needed. The only geometric primitives used by our algorithm are the {\InCircle} test and the predicate $\CCW(A,B,C)$, which is true if the points $A$, $B$, and $C$ form a counterclockwise oriented triangle. The latter is a standard tool in geometric algorithms, and is equivalent to reporting on which side of a line a point lies. It is equivalent to $\InCircle(A,B,C,D)$, for $D$ chosen as the barycenter of triangle $ABC$.\par

As one might expect, in the divide and conquer algorithm we start by partitioning our points into two halves, the left half ($L$) and the right half ($R$), which are separated in the $x$-coordinate. We next recursively compute the Delaunay triangulation of each half. Finally we need to marry the two half triangulations into the Delaunay triangulation of the whole set.\par

We now elaborate on this brief description in stages. First of all it is advantageous to start out by sorting our points in increasing $x$-coordinate. When there are ties we resolve them by sorting in increasing $y$-coordinate and throwing away duplicate points. This makes all future splittings constant time operations. After splitting in the middle and recursively triangulating $L$ and $R$, we must consider the merge step. Note that this may involve deleting some $L-L$ or $R-R$ edges and will certainly require adding some $L-R$ or {\cross} edges. By lemma 7.3, however, no new $L-L$ or $R-R$ edges will be added.\par

What is the structure of the cross edges? All these edges must cross a line parallel to the $y$ axis and placed at the splitting $x$ value. This establishes a linear ordering of the cross edges, so we can talk about successive cross edges, the bottom-most cross edge, etc. The algorithm we are about to present will produce the cross edges incrementally, in ascending $y$-order. See fig. 9.1.\par

Lemma 9.1 has the following important consequence. Let us call the current cross edge the base, and write its directed variant going from right to left as {\basel}. The successor to {\basel} will either be an edge going from the left end-point of {\basel} to one of the $R$-neighbors of the right end-point lying ‘‘above’’ {\basel}, or symmetrically, it will be an edge from the right end-point of {\basel} to one of the $L$-neighbors of the left end-point lying ‘‘above’’ {\basel}. In the program below edges from the left end-point of {\basel} to its candidate $L$-neighbors will be held in the variable {\lcand}, and their symmetric counterparts in {\rcand}. See fig. 9.2.\par

We can intuitively view what happens by imagining that a circular bubble weaves its way in the space between $L$ and $R$ and in so doing gives us the cross edges. Inductively we have a point free circle circumscribing the triangle defined by {\basel} and the previous cross edge. Consider continuously transforming this circle into other circles having {\basel} as a chord, but lying further into the half-plane above {\basel}. As we remarked, there is only a single degree of freedom, as the center of the circle is constrained to lie on the bisector of {\basel}. See fig. 9.3. Our circles will be point free for a while, but unless {\basel} is the upper common tangent of $L$ and $R$, at some point the circumference of our transforming circle will encounter a new point, belonging either to $L$ or $R$, giving rise to a new triangle with a point-free circumcircle. The new $L-R$ edge of this triangle is the next cross edge determined by the body of the main loop below.\par

In more detail, edge {\lcand} is is computed so as to have as destination the first $L$ point to be encountered in this process, and {\rcand} the first $R$ point. A final test chooses the point among these two that would be encountered first. We start the cross edge iteration by computing the lower common tangent of $L$ and $R$, which defines the first cross edge. The edge {\ldi} is the clockwise convex hull edge out of the rightmost $L$ point, and {\rdi} is the counterclockwise convex hull edge out of the leftmost $R$ point. We assume these are returned by the recursive calls of the triangulation procedure to $L$ and $R$. In the program below prose in $\{\}$ is comments; certain portions of code that are just symmetric variants of others shown are elided and are indicated by ‘‘\prg{. . .}’’.\par
\vskip0pt plus1000pt
\hsize\pagewidth pt\eject

$\{${\rm initialize {\lcand} to the counterclockwise first edge incident to {\basel$\dest$} and above {\basel} }$\}$
$\{${\rm symmetrically for {\rcand} }$\}$
lcand \pgets basel\psym\ponext; rcand \pgets basel\poprev;
\null
$\{${\rm let \prg{ValidL} and \prg{ValidR} be procedures which test whether {\lcand} and {\rcand} are above {\basel} }$\}$
ValidL[l] : CCW[basel\porg, basel\pdest, l\pdest];
ValidR[r] : CCW[basel\porg, basel\pdest, r\pdest];
\null
$\{${\rm this is the merge iteration }$\}$
DO
\ \ $\{${\rm advance {\lcand} counterclockwise, deleting the old {\lcand} edge, until the {\InCircle} test fails }$\}$
\ \ IF ValidL[lcand] AND ValidL[lcand\ponext] THEN
\ \ WHILE InCircle[lcand\pdest, lcand\ponext\pdest, lcand\porg, basel\porg] DO
\ \ \ \ t ← lcand\ponext;
\ \ \ \ DeleteEdge[lcand];
\ \ \ \ lcand ← t;
\ \ \ \ OD;
\null

\ \ $\{${\rm connect to either {\lcand} or {\rcand} }$\}$
\ \ $\{${\rm if both are valid, then choose the appropriate one using the {\InCircle} test }$\}$
\ \ $\{${\rm if {\lcand} is the winner then make the next cross edge close the triangle with {\lcand} and {\basel}, }$\}$
\ \ $\{${\rm else make the next cross edge close the triangle with {\rcand} and {\basel} }$\}$
\ \ $\{${\rm advance {\basel} and reinitialize {\lcand} or {\rcand} }$\}$
\ \ IF NOT ValidL[lcand] OR
\ \ \ \ (ValidR[rcand] AND InCircle[lcand\pdest, lcand\porg, rcand\porg, rcand\pdest])
\ \ \ \ THEN BEGIN $\{${\rm connect to the right }$\}$
\ \ \ \ basel \pgets Connect[rcand, basel\psym, left];
\ \ \ \ rcand \pgets basel\psym\plnext;
\ \ \ \ END
\ \ ELSE . . . $\{${\rm connect to the left }$\}$
\ \ OD
\endprog}

\vfill\hsize\textwidth pt\eject

We now elaborate on the above program. Recall first that the number of vertices, edges, and faces in a triangulation are all linearly related. Also, the lower common tangent computation takes linear time as either {\ldi} or {\rdi} advances at each step. What about the cost of the {\lcand} computation? We can account for this loop by charging every iteration to the edge being deleted. Similarly, iterations of the {\rcand} loop can be charged to deleted edges. The rest of the body of the main loop requires constant time and may be charged to the $L-L$ or $R-R$ edge closing the next triangle. This shows that the overall cost of the merge pass is linear in the size of $L$ and $R$.\par

Assuming that the algorithm has worked correctly up to now, then the circumcircle of the triangle defined by {\basel} and the previous cross edge is the circumcircle of a Delaunay triangle and therefore it is point-free. As in our earlier discussion, suppose that we ‘‘push’’ this circle upwards while always making it pass through $A$ and $X$. At some point it will encounter a new $L$ point $N$ for the first time. This point will be above {\basel}. We wish to show that $N$ is the point that the {\lcand} iteration will discover. To see this note first of all that $N$ must be a Delaunay neighbor of $A$, since our circle is a Delaunayhood witness. The {\lcand} iteration proceeds until the first time $X$ falls outside the circle $ANǐN{̌i+1}$, as shown in the code. By the unimodality lemma, this gives us the neighbor $N=Nǐ$ which determines the circle $AXNǐ$ with the smallest extent above {\basel}. This circle now forms the basis for the next iteration through the main loop.\par

The inclusion of $X$ in the circle $ANm̌N{̌m+1}$ for $m = 1,2,\ldotss, m-1$ establishes that $ANm̌$ is not a Delaunay edge in the final triangulation, since each circle with $ANm̌$ as chord contains either $N{̌m+1}$ or $X$ in its interior. Thus the edge deletions performed are valid. It is easy to check that the cross edges drawn, together with the remaining edges of $L$ and $R$ define a triangulation of the $n$ sites.\par

This is essentially the proof, modulo expending some care on certain boundary cases. We must first of all show that the {\lcand} iteration, if started at all, never gets to neighbors on the bottom side of {\basel}. For if the iteration was invoked, then there must be at least two neighbors above {\basel}. We claim that the moment {$\lcand\ \onext$} crosses to the bottom of {\basel}, the {\fx WHILE} {\InCircle} test fails. See fig. 9.4. In this case we had that $X$ was inside the circle $AN{̌k-1}Nǩ$. If the angle $\angle NǩAN{̌k+1}$ is convex then $X$ is obviously outside the circle $ANǩN{̌k+1}$, so the {\InCircle} test fails. If the angle $\angle NǩAN{̌k+1}$ is concave, then two things happen. First of all the circle $ANǩN{̌k+1}$ reverses orientation, and secondly, since $N{̌k+1}$ is outside the circle $AN{̌k-1}Nǩ$ but on the same side of $ANǩ$ as $N{̌k-1}$ and $X$, $X$ now falls inside the circle $ANǩN{̌k+1}$. So in this case the {\InCircle} test fails also.\par

\lemma{9.3}{The triangle defined by {\basel} and the next cross edge as computed by the above algorithm has a circumcircle free of both $L$ and $R$ points, and thus is a Delaunay triangle. Furthermore, after the completion of the cross edge iteration, all remaining $L$ or $R$ triangles are also Delaunay triangles for the full configuration.}\par

\proof{The first statement is obvious from the above discussion. The second statement is most easily proven by showing that each remaining $L--L$ or $R--R$ edge is stationary under the swapping rule. The only possible suspect edges are those bordering triangles containg cross edges. It follows that these edges were {\lcand} or {\rcand} at some point in the algorithm and were not deleted. So the {\InCircle} test of the inner iterations failed there, and this establishes that these edges are safe under the swapping rule.}\par

These lemmas complete the proof that the above algorithm correctly computes the Delaunay triangulation. The algorithm will work with cocircular sites and other degenerate cases. When both {\lcand} and {\rcand} are equally good, it arbitrarily favors {\rcand}. In practice floating-point errors in the computation of the {\InCircle} test usually interfere with any effort to handle degenerate cases in a consistent way.\par

It is worthwhile to comment on a way to view this algorithm as operating in three-dimensional space, on the lifted images of our sites under the map $\lambda$ discussed in the proof of lemma 8.1. These lifted images are points on a convex surface, and therefore they define a convex polyhedron, corresponding to their convex hull. The discussion in the proof of lemma 8.1 has also established that the ‘‘downwards’’ looking faces of this polyhedron are in a one-to-one correspondance with the Delaunay faces of the sites. The upwards facing faces of the polyhedron correspond to triangles of our collection of sites whose circumcircles enclose all the sites. These are of course the faces of the dual of the {\it furthest-point} Voronoi diagram [\Sham] for our collection of sites. Note that when our cross edge iteration ends, the ascending circle is the half-plane above the upper common tangent of $L$ and $R$. The half-plane below this tangent is a circle containing all the sites. If we now let this circle contract until it hits the first site, while always having as chord the last cross edge, we will have produced the next cross edge of the dual of the furthest point Voronoi. In fact, if our cross edge iteration is left to continue, but by using the furthest-point versions of $L$ and $R$ and reversing the sense of the {\InCircle} test, it will compute all the cross edges of this dual and will cycle back to the lower common tangent of $L$ and $R$.\par

From the above discussion it is apparent that our divide-and-conquer algorithm is computing the convex hull of the lifted images of the sites. It is in fact exactly the Preparata-Hong [\PH] algorithm for computing the convex hull of $n$ points in three dimensions. If the {\InCircle} test is replaced by a ‘‘negative volume’’ test, as obtained by substituting in $\Dscr(A,B,C,D)$ the third column by the $z$ coordinates of the points, then the code given above implements the Preparata-Hong algorithm! The reader may have fun verifying that our expanding circles passing through a chord become rotating supporting planes around the lifted image of the chord. Thus we are computing the ‘‘sleeve’’ discussed in [\PH].\par

The algorithm of the previous section assumes that all points are known at the beginning of time. For many applications we are interested in dynamic algorithms, which allow us to update our diagrams as new points are added or deleted. The machinery we have developed so far yields also simple and elegant insertion algorithms. Deletions present problems.\par

It will simplify our discussion to assume that our points are always within some large convex polygon (say a triangle) whose vertices are considered to be among the given points. Thus the entire interior of this polygon will have been triangulated, and when a new point is given, we can be sure that it will fall within one of the existing triangular faces. The algorithm for updating the Delaunay triangulation starts out by locating the new point in one of the triangular faces of the subdivision. This can be done by Kirkpatrick’s method in time $O(\lg n)$ [\Kirki]. Our concern is with how the triangulation is to be updated.\par

WHILE the stack is nonempty DO
\ \ BEGIN
\ \ \ \ $\{${\rm each edge on the stack forms a triangle with the new point }$\}$
\ \ remove an edge from the stack;
\ \ consider the quadrilateral whose diagonal
\ \ \ \ is this edge;
\ \ if this edge passes the swapping rule test,
\ \ \ \ then forget about it,
\ \ \ \ else apply the swapping rule and stack,
\ \ \ \ \ \ in ccw order, the two edges not
\ \ \ \ \ \ incident to the new point;
\ \ OD;
\endprog}

At any one time the set of triangles with the new point $X$ as a vertex form a star-shaped polygon around $X$. This is because the polygon can only get larger with each application of the swap rule on a bounding edge, and if it does grow, then the two new stacked sides are in the angular sector with vertex $X$ and delimited by the deleted edge. Furthermore, because we stack edges in a consistent order, the stacked edges always form a contiguous section of the perimeter of that polygon, with the rest of the perimeter consisiting of the forgotten edges.\par

\proof{Let $X$ be the new point and suppose that the swap rule was just applied to quadrilateral $XLMN$ where it replaced diagonal $LN$ with diagonal $XM$, as in fig. 10.2, Then $X$ must be interior to the circle $LMN$ and by the same argument as in the proof of lemma 10.1 it follows that $XM$ is a Delaunay edge. If the swap rule failed on $XLMN$, then $X$ is outside the circle $LMN$, and so edge $LN$ is Delaunay and can be forgotten.}\par

\proof{One of two things can happen at each application of the swapping rule. If the swap succeeds, then the polygon grows past the old bounding edge, so that edge can never be on its boundary again. If it fails, then that bounding edge is frozen and the algorithm will never again touch any edge in the angular sector with vertex $X$ and delimited by that edge.}\par

\proof{From lemma 10.2 we know that all edges leading to the new point $X$ and introduced during the above algorithm are Delaunay edges. The only other edges bounding new triangles are the edges of the star-shaped polygon discussed in the proof of lemma 10.3. At the termination of the loop the stack is empty and so all the bounding edges have been forgotten at some point. Thus there are no suspect edges left.}\par

The above lemmas imply that this loop correctly computes the final Delaunay triangulation, and that it does so in at most $O(n)$ time. More precisely, if $k$ edges need to be added or deleted to update the Delaunay structure, then our algorithm works in time $O(k)$, since each edge added by the algorithm must be added, and each edge that is deleted must be deleted.\par

We saw that the forgotten edges will always define a path around the new point, which will eventually close. If $e$ is an edge on the path, then the next edge is $e\ \lnext\ \oprev$. This observation gives us a way to implement the above algorithm without any auxiliary storage, such as the stack, by a method that is similar to that used in the previous section. Here {\basel} starts out (arbitrarily) as one of the edges connecting the new point to one of the vertices of its enclosing triangle. By convention we think of the new point as being the left end point of {\basel}. Then {\lcand} will always be invalid (non-existent, in fact) and we will consider for {\rcand} the various edges out of the other endpoint of {\basel}. The successively found {\rcand}$\,$s, as we proceed counterclockwise around the new point, will correspond to the forgotten edges of the previous algorithm. Note, however, that the incremental algorithm ‘‘connects ahead’’ after each deletion, while the algorithm of the previous section would connect all cross edges in strict counterclockwise order.\par

Something quite general is happening here. We have a method with which, given any two Delaunay triangulations $L$ and $R$ (not necessarily linearly separated), and a cross edge between them, we are able to find the next cross edge (if one exists) on a specified side of the original one. Thus we have another way to look at Kirkpatrick’s [\Kirkii] linear time merge of two arbitrary Voronoi subdivisions.\par

The above arguments show that it is possible to insert a new site into the Delaunay structure in total time $O(k)$, if $k$ updates need to be made. Unfortunately we know of no $O(k)$ algorithm for handling the deletion of a site that leaves an untriangulated face of $k$ sides. Our best algorithm has asymptotic complexity $O(k\lg k)$, which in the worst case $k=\Theta(n)$ is as bad as rebuilding the subdivision from scratch. We do not know of a linear algorithm even if we assume that the deletion of the site leaves a convex face. We regard the handling of deletions as the major open problem in this area.\par

In this paper we have presented a new data structure for planar subdivisions which simultaneously represents the subdivision, its dual, and its mirror image. Our {\quadedge} structure is both general (it works for subdivisions on any two-dimensional manifold) and space efficient. We have shown that two topological operations, both simple to implement, suffice to build and dismantle any such structure.\par

We have also shown how, by using the {\quadedge} structure and the {\InCircle} primitive, we can get compact and efficient Voronoi/Delaunay algorithms. The {\InCircle} test is shown to be of value both for implementing and reasoning about such algorithms. The code for these algorithms is sufficiently simple that we have practically given all of it in this paper.\par