\input TR10GBWideFormat \input QuadMacros \begfile{GBb08.tex of November 1, 1983 4:41 pm --- Stolfi} \chapter{\Ch8}{Closest-Point Problems} \shorttitle{Closest-point problems} \def\omitfigures{T} % Things to do: % fix numbers of theorems etc. % change symbol for Delaunay triangulation to T instead of D\star % Replace ``outer boundary'', ``exterior boundary'' by ``perimeter'' % "Subdivision" or "simple subdivision"? % additional macros \def\E{\Escr} % \E(p,q), \E K = Bisector of p and q (resp K) \def\C{\Cscr} % \C(p,q) = Locus of points closer to p than to q \def\V{\Vscr} % \VP = Voronoi of P, \VP K = region of it def. by K \def\D{\Dscr} % \DP = Delaunay of P, \DP K = region of it def. by K \def\T{\Tscr} % \TP = A Delaunay triangulation \def\L{\Lscr} % The lower portion of \hull{\Lambda(P)} \def\H{\Hscr} % A convex hull \def\F{\Fscr} % A facet of a convex hull. Maybe should use F instead. \def\plc#1{{\langle #1 \rangle}} % \plc{X,Y,Z\hom W} = Plane Xx+Yy+Zz+W=0. \def\hom{\relvv } % Separates weight in homogeneous coordinates. \def\posor{\mathop{\char P\char O\char S}} \def\negor{\mathop{\char N\char E\char G}} \def\plane{\mathop{\char p\char l\char a\char n\char e}} \def\circle{\mathop{\char c\char i\char r\char c\char l\char e}} \def\outcircle{\mathop{\char O\char u\char t\char C\char i\char r\char c\char l\char e}} \def\base{{$e$}} \def\newbase{{$e\pr$}} \def\root{{$\mathop{\char r \char o\char o\char t}$}} % reference macros \def\Sham{Sh} \def\Kirki{Ki1} \def\Kirkii{Ki2} \def\ShamHo{SH} \def\GuiSto{GuiSto} \def \GrSi{GS} \def\Baum{Ba} \def\Braid{Br} \def\IyKa{IK} \def\MuPr{MP} \def\MaSu{MS} \def\Kle{Kel} \def\Knuth{Kn} \def\Even{Ev} \def\Lee{Le} \def\Har{Har} \def\PH{PH} \def\Sei{Sei} % Variables of the program \def\lcand{\hbox{\it lcand}} \def\rcand{\hbox{\it rcand}} \def\basel{\hbox{\it basel}} \def\ldi{\hbox{\it ldi}} \def\rdi{\hbox{\it rdi}} % Begin % Begin {\it Where is the nearest post-office?\/} We may consider this question as classical instance of the so-called {\it post-office problem}: given the coordinates of a finite set $P$ of points (the {\it sites}) and an arbitrary point $q$ (the {\it query}), find the site in $P$ that is closest to $q$. \figspace3cm To make this problem precise, we must define the nature of sites and queries and what is meant by ``closest''. In other words, we must specify the space $\Sscr$ to which the points belong, and the metric $\dist$ used to measure distances between them. Unless stated otherwise, we will assume $\Sscr$ is the plane $\reals^2$, and $\dist$ is the familiar Cartesian metric. Generalizations to spaces of higher dimension and to other metrics will be discussed briefly in section \Sec6. The straightforward algorithm for the post-office problem computes $\dist(p,q)$ for each site $p$, and selects the site for which the distance is minimum. If there are $n$ sites this algorithm takes $\Theta(n)$ time; since all sites must be examined, the algorithm is optimal. \section{\Ch8.\Sec1. The Voronoi and Delaunay diagrams} % Begin There are however several variants of the basic problem for which this na\um{\i}ve algorithm is no longer the best. An example is the {\it batched} version, in which many distinct query points are to be tested against a fixed set $P$ of sites. Another example is given by certain {\it on-line} applications where the set $P$ is given in advance of the queries, and the the goal is to minimize the response time to each query, rather than the total time. In such situations, a better strategy is to pre-process $P$ into a suitable data structure that allows each query to be anwered more efficiently. One possible way to achieve this is to precompute, for each site $p$ of $P$, the set of all queries $q$ for which $p$ is the closest site. This set is called the {\it Voronoi region} associated with $p$, and the collection of those regions is the {\it Voronoi diagram} of $P$. The post-office problem then is reduced to that of locating the Voronoi region that contains a given query point $q$. The figure below shows ten sites ($\bullet$) and their Voronoi regions (bounded by solid lines).\numfig5cm{[\Fig1]} From this example we can see that each Voronoi region is a convex polygon containing the corresponding site. Indeed, the Voronoi diagram of $n$ sites is a simple subdivision of the plane with $n$ convex faces and $O(n)$ edges and vertices (these and other properties of Voronoi diagrams will be proved in the following sections). As we saw in chapter \Ch5, we can build from such a subdivision a data structure that allows each query to be answered in $O(\log n)$ time. The construction of the data structure takes $O(n)$ time and $O(n)$ space. We will show in what follows that construction of the Voronoi diagram takes $O(n\log n)$ time and $O(n)$ space. Therefore, by combining those results we get an algorithm for the on-line post-office problem that requires only $O(n\log n)$ preprocessing time, $O(n)$ space, and $O(\log n)$ time per query. This algorithm also solves the batched vesion, with $n$ sites and $m$ queries, in $O((m+n)\log n)$ total time. \footprob{\Prob{\Sbsec{\Sec1.0}.1}}{Give an $O(m\log m + n\log n)$ algorithm for the batched post-office problem that does not require the construction of Kirkpatrick's structure (or its equivalent).\lp {\it Ans.}: Compute the Voronoi of $P$, sort the query points, and use line sweep.} If we connect two sites by a straight line segment whenever their corresponding Voronoi regions have a common edge, we obtain the so-called {\it Delaunay diagram} of the set $P$. Therefore, every Delaunay edge corresponds to a Voronoi edge, and vice versa. The Delaunay diagram for the ten sites of the previous example is given below (solid lines), superimposed on their Voronoi diagram (broken lines).\numfig5cm{[\Fig2]} This example illustrates an interesting property of the Delaunay diagram: its edges do not cross. Moreover, the ordering of the Delaunay neighbors around a given site $p$ agrees with the ordering of their Voronoi regions around the region of $p$. These facts are all the more surprising if we consider that a Delaunay edge need not intersect the corresponding Voronoi edge, as the figure shows. Together, they imply that the Delaunay diagram is a subdivision of the plane that is topologically the dual of the Voronoi diagram. These assertions will be made precise and proved in the following sections. \note{Ramblings on the importance of Delaunay} The Voronoi diagram is an important tool, but certainly not a panacea. In spaces of dimension $d\geq 3$, the Voronoi diagram may have up to $\Theta(n^{\lceil d/2\rceil})$ elements (Klee et al., [\Kle, \Sei]). Even in the plane, it is a somewhat bulky structure whose construction, storage, and manipulation may be hard to justify in some applications. Other techniques based on hashing or on general data structures like quad trees may be more efficient in such cases, especially if the input points are randomly chosen from a smooth probability distribution (see section \Sec?). \note{Anything else?} The post-office problem has been studied by many famous mathematicians in the last two centuries \note{give names and dates.} ? Voronoi introduced the diagram that carries his name in ????; ? Thiessen and ? Dirichlet discovered the same idea independently at about the same time (????), and for this reason some researchers (particularly in Europe) call it the {\it Dirichlet} or {\it Thiessen} {\it tesselation}. It has inspired many important concepts of computational geometry, and has a surprisingly large (and still growing) set of applications. \note{Fill with references} \subsection{\Ch8.\Sbsec{\Sec1.1}. The Voronoi diagram} The Voronoi and Delaunay diagrams play a central role in this chapter, so we will spend the rest of this section in the study of their properties. First, let us establish a notation for the elements (faces, edges, and vertices) of Voronoi diagrams. \definition{\Def1}{Let $P\subset \Sscr$ be a collection of sites, and $K$ a subset of $P$. The {\it Voronoi region of $K$}, denoted by $\VPK$, consists of all points $x$ of $\Sscr$ that are equidistant from all the sites in $K$, and closer to them than to any other site of $P$. The function $\VP$ is called the {\it Voronoi diagram} of $P$.} We will write $\V K$ instead of $\VPK$ when $P$ is known from the context. For any site $p$, the region $\V\set{p}$ is obviously the Voronoi face associated with $p$, i.e. the set of all points $x$ for which $p$ is the (only) closest site. The interpretation of $\V K$ when $K$ has more than one site is also quite natural. The expression $\V\set{p1, p2}$ (with $p1\neq p2$) denotes the edge between the faces $\V\set{p1}$ and $\V\set{p2}$, or $\empty$ if the two faces are not adjacent. The region $\V\set{p1, p2, \ldots, pm}$, with $m\geq 3$, consists of the single vertex (if any) that is common to $\V\set{p1},\V\set{p2}, \ldots, \V\set{pm}$ and to no other face of the Voronoi.\figspace5cm These statements will be proved in the following sections. \footprob{\Prob{\Sbsec{\Sec1.1}.1}}{What is $\VP\empty$?\lp {\it Ans.}: $\VP\empty$ is $\Sscr$ if $P=\empty$, and empty otherwise.} \footprob{\Prob{\Sbsec{\Sec1.1}.2}}{(a) Describe $\VP$ when $\Sscr=\reals$ (i.e., the real line). (b) Give an algorithm for the post-office problem that requires $O(n\log n)$ preprocessing time, $O(n)$ space, and $O(\log n)$ query time. (c) Show that these are lower bounds for the post-office problem in any Cartesian metric space (i.e., $\reals^k$) with any metric of the sort $\dist(x, y)=f(x-y)$.\lp {\it Ans.}: (a) List of intervals, defined by midpoints of consecutive sites. (b) Sort the sites and do binary search. (c) Consider $n$ points in a line, and show that to know the Voronoi one needs to know the order of the points on the line.} \subsection{\Ch8.\Sbsec{\Sec1.2}. Nature of Voronoi regions} We begin by introducing a couple of auxiliary concepts. For any two points $p$ and $q$, let us denote by $\E(p, q)$ the set of all points of the plane that are equidistant from $p$ and $q$. Let us also denote by $\C(p,q)$ the set of those points which are closer to $p$ than to $q$. As we know, $\E(p, q)$ is the {\it bisector} of $p$ and $q$, that is, the straight line perpendicular to the segment $pq$ and passing through its midpoint, while $\C(p,q)$ is the open half-plane delimited by $\E(p,q)$ and containing $p$. See exercise \Prob{\Sbsec{\Sec1.2}.2} for a proof of this fact.\figspace3cm \footprob{\Prob{\Sbsec{\Sec1.2}.1}}{What are $\E(p,p)$ and $\C(p, p)$?\lp {\it Ans.:} $\E(p,p)=\Sscr$ and $\C(p,p)=\empty$.} \footprob{\Prob{\Sbsec{\Sec1.2}.2}}{Prove that $\E(p,q)$ is the line perpendicular to the segment $pq$ and passing through its middle point, and that $\C(p,q)$ is the half-plane delimited by $\E(p,q)$ and contining $p$.\lp {\it Ans.}: Let $x$ be a point of the plane, $l$ be the line through $p$ and $q$, and $o$ be the midpoint of $pq$. Let $y$ be the orthogonal projection of $x$ onto $l$. (see figure). Since $xy$ is perpendicular to $l$, we have\lp \boxdisp{\eqalign{\bar{xp}\leq\bar{xq} \rmiff\bar{xp}^2\leq\bar{xq}^2\cr \null\rmiff\bar{xy}^2 + \bar{yp}^2\leq\bar{xy}^2+\bar{yq}^2\cr \null\rmiff\bar{yp}^2\leq\bar{yq}^2\cr \null\rmiff\bar{yp}\leq\bar{yq}.\cr}}\lp By analogous reasoning we conclude that $\bar{xp}\geq\bar{xq}\rmiff\bar{yp}\geq\bar{yq}$, and therefore\lp \boxdisp{\eqalign{x\in E(p, q)\rmiff\bar{xp}=\bar{xq}\rmiff\bar{yp}=\bar{yq}\cr x\in C(p,q)\rmiff\bar{xp}<\bar{xq}\rmiff\bar{yp}<\bar{yq}.\cr}}\lp Let $\pi$ be the line perpendicular to $l$ and passing through $o$; let $L$ and $R$ be the half-planes bounded by $\pi$ and containing $p$ and $q$, respectively. Obviously $\bar{yp}=\bar{yq}$ happens if and only if $y=o\in\pi$, and $\bar{yp}<\bar{yq}$ if and only if $y\in L$. Since $xy$ is parallel to $\pi$, we have $y\in\pi\rmiff x\in \pi$, and $x\in L\rmiff y\in L$. We conclude that $\E(p,q)=\pi$ and $\C(p, q)=L$.} If $P$ consists of only two sites $p$ and $q$, its Voronoi diagram is quite simple: the faces are the two half-planes $\V\set{p}=\C(p,q)$ and $\V\set{q}=\C(q,p)$, and the only edge is $\V\set{p, q}=\E(p, q)$, as follows immediately from the definitions (see figure above). If $P$ consists of three sites, we will have either of the two cases below, depending on whether they are collinear or not.\figspace3cm Note that in the second case the single Voronoi vertex (which is the circumcenter of the three points) may lie inside or outside the triangle determined by them, depending on whether the latter is acute or obtuse. The figure below shows the the topologically distinct Voronoi diagrams of four sites:\figspace5cm Let us consider in more detail the last of these examples. The face surrounding site $p$ consists of the points that are closer to $p$ than to $q$, $r$, or $s$, that is, $$\V\set{p}=\C(p,q)\inte\C(p,r)\inte\C(p,s).$$ The edge between $p$ and $q$ consists of all points that are equidistant from those two sites, but closer to them than to $r$ or $s$; that is, $$\V\set{p,q}=\E(p,q)\;\inte\;\C(p,r)\inte\C(q,r)\inte\C(p,s)\inte\C(q,s).$$ Finally, the Voronoi vertex $v$ is equidistant from the three sites $p,q$, and $r$, and closer to them than to $s$. This means $s$ is the only point in $$\V\set{p,q,r}=\E(p,q)\inte\E(p,r)\inte\E(q,r)\;\;\inte\;\;\C(p,s)\inte\C(q,s)\inte\C(r,s).$$ Note that the three bisectors $\E(p,q)$, $\E(p,r)$, and $\E(q,r)$ intersect at a common point, the circumcenter of the triangle $pqr$. More generally, in the Voronoi diagram of an arbitrary collection $P$ of sites, any region can be expressed as the intersection of bisectors and half-planes of the form $\E(p,q)$ and $\C(p,q)$, respectively: \lemma{\Th1}{For any nonempty subset $K$ of $P$, $$\VPK=\inter{p,q\in K} \E(p,q)\;\inte\;\inter{{\scriptstyle r\in K}\atop{\scriptstyle s\in P-K}} \C(r, s).\eqno(\Eq2)$$} Lemma \Th1 trivially follows from the definitions and has some important consequences: \theorem{\Th3}{For any site $p\in P$, the region $\V\set{p}$ is an open convex polygon containing $p$.} \proof{According to lemma \Th1, we have $$\V\set{p}=\inter{s\in P-\set{p}} \C(p, s).\eqno(\Eq2)$$ So, the Voronoi region of $\set{p}$ is the (non-empty) intersection of a finite number of half-planes $\C(p, s)$, and therefore by theorem \Th{\Ch3.?} it is a convex polygon\foot\dag{Note that the polygon may be unbounded; in particular, if $p$ is the only site in $P$, $\V\set{p}$ is the whole plane.}. Since each $\C(p,s)$ is an open set containing $p$, the same is true of their intersection $\V\set{p}$.} \theorem{\Th4}{For any two distinct sites $p, q\in P$, the region $\V\set{p,q}$ is an open interval (possibly empty or unbounded) of the bisector $\E(p,q)$.} \footprob{\Prob{\Sbsec{\Sec1.1}.3}}{Prove theorem \Th4.\lp {\it Ans.}: From lemma \Th1,\lp \boxdisp{\V\set{p, q}=\E(p, q)\;\inte\;\inter{{\scriptstyle r\in \set{p, q}}\atop{\scriptstyle r\in P-\set{p, q}}} \C(r, s).\boxeqno(\Eq3)}\lp The second term of equation (\Eq3) is the intersection of open half-planes, and therefore is an open convex set by theorem \Th{\Ch3.?}. The region $\V\set{p,q}$ is the intersection of this set and the line $\E(p,q)$, and therefore is a convex subset of $\E(p,q)$, i.e. an open interval.} \theorem{\Th5}{If $K$ contains three or more sites of $P$, the region $\V K$ is either empty or contains a single point, which is the center of a circle passing through all sites in $K$.} \footprob{\Prob{\Sbsec{\Sec1.1}.4}}{Prove theorem \Th5.\lp {\it Ans.}: According to lemma \Th1, the region $\V K$ is contained in the set\lp \boxdisp{E=\inter{p,q\in K} \E(p,q).}\lp If $K$ contains three distinct points $p, q, r$, then $E$ is contained in the intersection of the three bisectors $\E(p,q)$, $\E(p, r)$ and $\E(q, r)$. These three lines are all distinct, and therefore their intersection is at most a single point. The same then is true of $E$, and consequently of $\V K$. Finally, any point in $\V K$ is by definition equidistant from all sites in $K$, and therefore is the center of a circle passing through all those sites.} Every point of the plane belongs to exactly one region $\V K$, so the Voronoi diagram is a partition of the plane into a finite number of convex polygons ({\it faces}), open line segments ({\it edges}), and isolated points ({\it vertices}). Combining this fact with the three theorems above, we can conclude the following: \theorem{\Th6}{A point $x$ of the plane belongs to a face, an edge, or a vertex of the Voronoi diagram $\VP$ depending on whether there are one, two, or more than two sites in $P$ at minimum distance from $x$.} \subsection{\Ch8.\Sbsec{\Sec1.3}. Unbounded regions} Since the Voronoi diagram covers the whole plane with a finite number of regions, one or more of them must be unbounded. These regions are easily characterized: \theorem{\Th7}{A Voronoi region $\VP\set{p}$ of a site $p\in P$ is unbounded if and only if $p$ lies on the boundary of\/\ $\hull(P)$.} \proof{Let $p$ be on the boundary of $\hull(P)$; by theorem \Th{\Ch3.?HC}, there is a supporting line $l$ of $P$ passing through $p$. Let $r$ be the ray perpendicular to $l$, with origin $p$, and on the side of $l$ opposite to $P$.\figspace4cm For any point $x\in r$, and any point $q\in P-\set p$, we obviously have $\dist(x,p)\leq\dist(x,q)$, so $x\in Vp$. This implies $r\subset \V\set{p}$, so $V\set{p}$ is unbounded. \pp Conversely, let's assume $\V\set{p}$ is unbounded. By theorem \Th{\Ch3.?R} there must be some ray $r$ emanating from $p$ that is entirely in $\V\set{p}$. Let $q$ be any point of $P-\set{p}$, and let $q1$ be the projection of $q$ onto the line of $r$. \figspace3cm \pp Suppose $q1$ is on the ray $r$ and distinct from $p$; then, by taking a point $x$ on $r$ sufficiently far away from $p$, we can make $\dist(x,p)>\dist(x,q)$, contradicting the fact that $r\subset \V\set{p}$. We conclude that $q1$ is on the ray opposite to $r$, that is, $q$ is in the half-plane bounded by the perpendicular $l$ to $r$ at $p$ and opposite to $r$. Since this holds for all $q\in P-\set{p}$, $l$ is a supporting line of $P$ and, by lemma \Th{\Ch3.?H}, $p$ is on the boundary of $\hull P$.} Note that $p$ doesn't have to be a vertex of the hull for $\V\set{p}$ to be infinite; it may lie on an edge instead. A similar argument shows that a Voronoi edge $\V\set{p,q}$ is unbounded if and only if $p$ and $q$ are consecutive sites on the boundary of $\hull(P)$, that is, if and only if the segment $pq$ is contained in an edge of $\hull(P)$ and is free of other sites. The discussion that follows will be simplified if we add to the plane a dummy point at infinity $\omega$, which transforms it into a compact manifold (the {\it extended plane}) topologically equivalent to the surface of a sphere. By convention, we will consider $\omega$ to be an additional Voronoi vertex, incident to all unbounded edges and faces of the Voronoi. \note{State also that from now on we will consider only set $P$ and Voronoi diagrams with {\em at least two sites}.} With this convention every edge has two (not necessarily distinct) endpoints, every face of the Voronoi is bounded by a closed chain of edges and vertices, and so forth. Topological considerations (see problem \Prob{\Sbsec{\Sec6.1}.6}) readily show that the endpoints of every edge (if any) are Voronoi vertices, and therefore the boundary of a face is the union of edges and vertices. In brief, we have: \theorem{\Th{7B}}{The vertices, edges and faces of the Voronoi diagram constitute a simple subdivision of the extended plane} This way of extending the plane to a compact manifold is not entirely without problems, since not all properties of Voronoi vertices are valid for $\omega$. In particular, $\omega$ does not satisfy theorem \Th6 and many of its corollaries. \subsection{\Ch8.\Sbsec{\Sec1.4}. The circle property} To complete the characterization of Voronoi regions, we will give a condition for $\V K$ to be non-empty. Let us say that a circle is {\it $P$-free} if it contains no points of $P$ in its interior. The following theorem is an immediate consequence of the definitions: \theorem{\Th8}{A point $x$ of the plane is in a Voronoi region $\VPK$ if and only if $x$ is the center of a $P$-free circle whose boundary contains $K$ but no other sites of $P$.} So, $\VP K$ is nonempty if and only if all sites in $K$, and only those, lie on the boundary of some $P$-free circle. The two special cases below are particularly important: \theorem{\Th{8A}}{Two Voronoi faces $\V\set{p}$ and $\V\set{q}$ have a common edge if and only if there is a $P$-free circle passing through $p$ and $q$ but through no other site of $P$.\figspace3cm} \theorem{\Th{8B}}{Three or more regions $\V\set{p1}, \V\set{p2}, \ldots, \V\set{pm}$ have a finite vertex $v$ in common if and only if there is a $P$-free circle with center $v$ passing through $p1, p2, p3, \ldots, pm$.\figspace3cm} \subsection{\Ch8.\Sbsec{\Sec1.5}. The Delaunay diagram} \note{This subsection needs to be shortened. For each lemma and theorem in here, consider (a) replacing it by an unjustified statement, (b) replacing it by a justified statement, (c) stating it but leaving the proof as an exercise, or (d) keeping it as it is.} As we remarked before, the Delaunay diagram is derived from the Voronoi in a very straightforward way. \definition{\Def{1A}}{The {\it Delaunay diagram} of $P$ is a graph with the given sites as vertices and straight line segments as edges, such that two sites are connected by an edge if and only if the corresponding Voronoi faces share a common Voronoi edge.} A characterization of Delaunay edges readily follows from theorem \Th{8A}: \theorem{\Th{8C}}{Two sites $p,q$ are connected by a Delaunay edge if and only if there is a $P$-free circle passing through $p$ and $q$ and through no other site of $P$.\figspace3cm} We can say that such a circle is a ``certificate of Delaunayhood'' for that edge. Lemma \Th{8C} has the following corollary: \theorem{\Th9}{No two Delaunay edges intersect\foot1{\note{This seems to overlap in part the proof of the angle/quadrialteral property of $\incircle$.}}.} \proof{Suppose two distinct Delaunay edges $pq$ and $rs$ intersect at some point $o$. By theorem \Th{8A}, there is a $P$-free circle $C$ having $pq$ for chord and passing through no other site. There is also a similar circle $D$ for the edge $rs$. \pp Let us consider first the case where the four sites are collinear. Since the two open segments $pq$ and $rs$ intersect, at least one of them must contain an endpoint of the other, say $r\in pq$. But then any circle $C$ with chord $pq$ will have $r$ in its interior, a contradiction. \pp Now let us consider the case where the four points are not collinear; since the open segments intersect, it follows all four sites are distinct. Therefore $p$ and $q$ must lie outside the circle $D$. On the other hand, the intersection point $o$ lies in $rs$, and therefore inside $D$. \figspace3cm \pp So, $D$ must intersect the segments $po$ and $qo$ at points $p\pr$ and $q\pr$, which are both interior to $C$ and lie on opposite sides of the line $rs$. Therefore, $D$ intersects $C$ at four (or more) distinct points. Since three points determine a unique circle, $C$ and $D$ must coincide, contradicting the fact that $C$ passes through $p$ and $q$ but $D$ does not.} This shows the Delaunay diagram is a graph properly embedded in the plane. Another property of the Delaunay diagram is the following: \lemma{\Th{9A}}{Let $X$ be a subset of $P$, and let $u\in X$ and $v\in P-X$ be such that the distance $\dist(u, v)$ is minimum. Then $uv$ is an edge of the Delaunay diagram of $P$.} \proof{Consider the circle $C$ with diameter $uv$, and suppose there was some site $w$ other than $u$ and $v$ inside $C$ or on its boundary:\figspace3cm \pp It is clear from the figure that $\dist(u,w)$ and $\dist(w,v)$ would be both strictly less than $\dist(u,v)$. The site $w$ would have to be in $X$ or in $P-X$; but either case contradicts the minimality of $\dist(u,v)$.} Theorem \Th{9A} has the following consequence: \theorem{\Th{9B}}{The Delaunay diagram is a connected graph.} Removal of the Delaunay edges and vertices breaks the plane into one or more open and connected subsets, the {\it faces} of the Delaunay. Since Delaunay edges are finite in number and length, exactly one of them is unbounded, the {\it exterior face}; all other faces (the {\it interior} ones) are bounded. As in the case of the Voronoi, it is convenient to consider the Delaunay diagram as embedded in the extended plane, and include the point at infinity $\omega$ into the exterior face. From theorem \Th{9B} we can then conclude that every face of the Delaunay is a connected subset of the extended plane, bounded by a single ring of edges and vertices. Therefore, every face is homeomorphic to a disk. In brief, we have the following: \theorem{\Th{9C}}{The Delaunay diagram is a simple subdivision of the extended plane.} Definition \Def{1A} immediately establishes a one-to-one correspondence between the edges of the two diagrams: the segment $pq$ is a Delaunay edge if and only if $\V\set{p,q}$ is non-empty, i.e., is a Voronoi edge. Since $\V\set{p,q}$ is contained in the bisector of $pq$, the two edges are perpendicular to each other\foot\dag{As is revealed by the examples given, the two edges need not cross each other. {\bf Gabriel graph?}.}. Thanks to this fact, we can extend the one-to-one correspondence unambiguously also to {\it directed} edges, in the following way: if $e$ is a directed edge of either diagram, then we define $e\rot$ to be the corresponding edge of the other, oriented $90\deg$ counterclockwise from $e$.\figspace3cm This definition is justified, since it establishes the formal duality (in the sense of chapter \Ch5) between the Voronoi and Delaunay diagrams. This fact follows mainly from lemma \Th{10} below: \lemma{\Th{10}}{The neighbors of a site $p$ in the Delaunay diagram are exactly the sites whose Voronoi regions are adjacent to $\V\set{p}$, in the same counterclockwise order.} \proof{Let $\V\set{p1}, \V\set{p2}, \ldots, \V\set{pm}$ be the Voronoi regions adjacent to $\V\set{p}$, in counterclockwise order. By definition, the Delaunay neighbors of $p$ are precisely $p1, p2, \ldots, pm$. All we have to show is that those sites occur in that order around $p$.\figspace4cm \pp Let $ei$, for $1\leq i\leq m$, be the Voronoi edge $\V\set{p, pi}$, directed counterclockwise as seen from $p$. Let also $e\pri$ denote the Delaunay edge $ppi$, directed away from $p$; that is, $e\pri=ei \rot^{-1}$. Now consider the direction vector $di$ of each Voronoi edge $ei$; since $\V\set{p}$ is a convex polygon, the points $d1, d2, \ldots, dm$ will lie on unit circle in precisely that order. But the direction of $e\pri$ is exactly $90\deg$ to the left of $di$. It follows the Delaunay edges $e\pr1, e\pr2, \ldots, e\prm$ (and therefore their destination sites $p1, p2, \ldots, pm$) lie in that order around $p$.} Stated in the notation of edge algebras, lemma \Th{10} says that $e\rot^{-1}\onext\rot=e\lnext$ for any directed Voronoi edge $e$. The next result can now be proved quite easily: \theorem{\Th{12}}{The Delaunay and Voronoi diagrams, as subdivisions of the extended plane, are dual to each other.} Among the many corollaries of this result, the following is noteworthy: \lemma{\Th{11}}{The sites $p1, p2, \ldots, pm$ are the vertices of an internal (resp. external) Delaunay face, in counterclockwise order, if and only if $\V\set{p1}, \V\set{p2}, \ldots, \V\set{pm}$ are the Voronoi faces incident to a finite (resp. infinite) Voronoi vertex, in counterclockwise order.} This gives us a valuable characterization of Delaunay faces: \theorem{\Th{13B}}{Three or more sites $p1, p2, p3, \ldots, pm$ lie on the same interior face of the Delaunay if and only if there is a $P$-free circle passing through all of them.\figspace3cm} From theorem \Th{13B}, and the fact that Delaunay edges do not cross, we conclude that every interior face of the Delaunay is an open convex polygon whose vertices lie on a common circle. As for the exterior face, we have the following: \theorem{\Th{14}}{The exterior face of the Delaunay is the exterior of $\hull(p)$} \proof{By theorem \Th{11}, a site $p$ is incident to the exterior face $F$ of the Delaunay if and only if $\V\set{p}$ is incident to the vertex at infinity $\omega$, i.e., if $\V\set{p}$ is unbounded. By theorem \Th7, this happens if and only if $p$ lies on the boundary of $\hull(P)$. From this and the fact that Delaunay edges do not cross the theorem follows.} We conclude immediately that \theorem{\Th{15}}{Any two consecutive sites on the boundary of $\hull(P)$ are adjacent in the Delaunay diagram of $P$.} The duality between the two diagrams suggests a uniform notation for Delaunay elements: \definition{\Def2}{For any subset $K$ of $P$, we define the {\it Delaunay region} $\DP K$ as the dual of the Voronoi region $\VP K$ if this is nonemty, and $\empty$ if $\VP K$ is empty.} As before, we write just $\D K$ when the set $P$ is implied by the context. For any site $p$, the region $\D \set{p}$ is just $\set{p}$, a Delaunay vertex. The region $\D\set{p,q}$ is the open segment $pq$ if the two sites $p, q$ are Delaunay neighbors, and is empty if they are not. For any set $K$ of three or more sites, $\D K$ is the interior face of $\DP$ whose vertices are exactly the sites in $K$, or empty if there is no such face. The next lemma is simply theorem \Th8 recast in this notation, and summarizes theorems \Th{8C} and \Th{13B}: \theorem{\Th{13}}{$\DP K$ is nonempty if and only if all sites in $K$, and only those, lie on the boundary of some $P$-free circle.} \subsection{\Ch8.\Sbsec{\Sec1.6}. Delaunay diagrams as convex hulls} The quadratic mapping we introduced in section \Sec{\Ch{I}.8} gives an invaluable insight into the nature of Delaunay diagrams. Consider the set $\lambda(P)$ consisting of the images on the paraboloid $\Lambda$ of all the sites in $P$. As we have seen, an oriented circle $C$ on the plane is mapped to the intersection of $\Lambda$ with some oriented plane $\pi$, and the left side of $C$ maps to the part of $\Lambda$ that lies to the right of $\piC$. If $C$ is finite and counterclockwise oriented, the left side is the interior, and the plane $\pi$ is oriented ``upwards'' (that is, its normal vector points towards increasing $z$ coordinates. Conversely, any upwards-oriented plane $\pi$ corresponds to a finite counterclockwise oriented circle $C$ on the $xy$ plane with the above properties. \figspace5cm A site $p$ of $P$ will be inside a circle $C$ if and only if $\lambda(p)$ lies below the corresponding plane $\pi$. We conclude that {\it a circle $C$ is $P$-free if and only if its corresponding plane $\pi$ supports $\lambda(P)$ from below.} In particular, let $C$ be the circumcircle of an interior face $F$ of $\DP$; since $C$ is $P$-free and passes through the vertices $p1, p2, \ldots, pm$ of $F$, its corresponding plane $\pi$ supports $\lambda(P)$ from below and passes through the points $\lambda(p1), \lambda(p2), \ldots, \lambda(pm)$. The polygon $\langle \lambda(p1), \lambda(p2), \ldots, \lambda(pm)\rangle$ on $\pi$ is strictly convex, since the same is true of its projection $F$ on the $xy$-plane. We conclude that {\it $\langle \lambda(p1), \lambda(p2), \ldots, \lambda(pm)\rangle$ is a face on the underside of the three-dimensional convex hull of $\lambda(P)$.} \figspace5cm The converse is also true, namely the vertices of any face on the underside of $\H=\hull{\lambda(P)}$ are the images of the sites on an interior face of $\DP$. Moreover, recall that two sites $p,q$ are adjacent on the Delaunay if and only if $pq$ is the chord of a $P$-free circle passing though no other sites. This happens if and only if there is a plane supporting $\lambda(P)$ from below, and passing through $\lambda(p)$, $\lambda(q)$, and no other point of $\lambda(P)$; in other words, if and only if $\lambda(p)\lambda(q)$ is an edge of the underside of $\H$. To complete the analysis, we have only to observe that every point of $\lambda(P)$ is a vertex on the underside of $\H$ (as follows from the strict concavity of the paraboloid $\Lambda$). In conclusion, {\it the Delaunay diagram of a set $P$ is isomorphic to the ``underside'' of the three-dimensional convex hull of $\lambda(P)$.} More precisely, \theorem{\Th{18}}{Let $\H$ be the three-dimensional convex hull of $\lambda(P)$. Any non-empty Delaunay region $\DP K$ is the projection on the $x,y$ plane of a facet $F$ on the underside of $\H$, whose vertices are the points in $\lambda(K)$. \figspace5cm} \note{Dual: Voronoi diagram is isomorphic to intersection of the half-spaces dual to the points in $\lambda(P)$. The face normals of the convex hull are dual to the Voronoi vertices.} \note{Mumble about importance of theorem \Th{18} for understanding Delaunays} \note{A thought: if this correspondence were stated before (say, just after theorem \Th{8C}), maybe it would be unnecessary to prove that the Delaunay is a subdivision, etc.} \section{\Ch8.\Sec2. Delaunay triangulations} % Begin The Delaunay and Voronoi diagrams contain the same geometrical and structural information, and we can efficiently obtain one from the other. In principle, any algorithm using or constructing a Voronoi diagram can be viewed as working on the Delaunay, and vice-versa. However, the Delaunay diagram is in a sense the more ``natural'' of the two; for example, its vertices are the original data points, and therefore we need not worry about introducing new vertices and points at infinity. For this reason, algorithms generally become simpler (and possibly more efficient) when described in the Delaunay form. The algorithms are further simplified when all interior faces of the Delaunay diagram are triangles. If the sites are independently sampled from a continuous two-dimensional distribution, the probability of four or more sites lying on the same circle (and, therfore, of a Delaunay face having four or more vertices) is zero. In addition, the rounding errors inherent in floating-point arithmetic make it pointless to check whether four points are exactly cocircular\foot\dag{\note{Actually, rounding errors are much nastier than that. Virtually every theorem of geometry and algebra must be revised and adequately qualified before it can be applied to floating-point computations. In the floating-point world, to get anything like a correct algorithm one must employ some special care and be content with ``backwards'' correctness (namely, there is a small perturbation of the input sites for which the output is the correct Delaunay diagram. However, our algorithm does not guarantee that.}}. In practical appliactions and in informal discourse we are therefore justified in assuming all interior faces of the Delaunay diagram to be triangular, and treating any exceptions as degenerate cases. Moreover, rather than handling those exceptions explicitly in the algorithms, it is generally simpler to work with a ``triangularized'' version of the Delaunay diagram, rather than with the diagram itself. This means breaking every Delaunay face with $k>3$ vertices into $k-2$ triangles by means of $k-3$ non-intersecting diagonal edges. The result will be a {\it triangular subdivision}, a straight-line subdivision of the plane whose interior faces are all triangles. \figspace 4cm A triangular subdivision vith vertex set $P$ and convex perimeter is called a {\it triangulation of $P$}. Since the perimeter of a Delaunay diagram is always convex, the same is true of all its triangularized versions. This justifies the following definition: \definition{\Def3}{A {\it Delaunay triangulation} of a point set $P$ is any triangular subdivision of the plane obtained by adding straight non-intersecting diagonal edges to the faces of the Delaunay diagram of $P$.} Unlike the Delaunay diagram, the Delaunay triangulation is not unique, a property that must be borne in mind when designing algorithms for its construction. However, all Delaunay triangulations of a given set $P$ of $n$ sites have the same number of elements: if $h$ is the number of sites that lie on the boundary of $\hull{P}$, then every triangulation of $P$ (Delaunay or not) has exactly $3(n-1)-h$ edges and $2(n-1)-h$ interior faces. Observe that those are obvious upper bounds for the number of edges and interior faces of the Delaunay diagram of $P$, and therefore also for the number of edges and finite vertices of the Voronoi. \digress{A general triangular subdivision need not have convex perimeter, as the examples below illustrate. Recall however that a subdivision must be connected and have at least one edge (and therefore, if the edges are straight, at least two vertices.) \figspace 3cm We warn the reader that triangular subdivisions with non-convex perimeter frequently arise as intermedite stages in the bulding of Dealunay triangulations. \dp A set $P$ whose sites lie all on a straight line has a unique (Delaunay) triangulation, consisting of a straight path with $n-1$ edges and no interior faces. \figspace 2cm In order to make the formulas above valid also for this special case, we must define $h$ by counting each site ``as many times as it is encountered when we go around the boundary of $\hull{P}$''. In other words, when all sites lie on the same line we must count twice each point, except for the first and last on the line. This will give $h=2(n-1)$, satisfying the formulas.} \digress{If we consider only its combinatorial or topological structure, a Delaunay triangulation $\T$ always has a well defined dual. O focurse, the dual subdivision can be drawn on top of $\T$ in the ``canonical'' way, with each dual vertex inside the corresponding primal face, and each dual edge crossing only its primal edge. In this subdivision all vertices have degree 3 (except for the dual of the exterior face of $\T$). \dp It is tempting to define the ``geometric dual'' of $\T$ as a specific subdivision that is {\it geometrically} related to it in the same way that the Voronoi diagram is related to the Delaunay. That is, the geometric dual of $\T$ would be a straight-line subdivision $\Wscr$ of the plane whose veritices are the circumcenters of the faces of $\T$, and such that two vertices are connected by an edge if and only if the corresponding faces of $\T$ share a common edge. \dp Unfortunately, any two triangles of $\T$ that are part of the same face of the Delaunay diagram have the same circumcircle, and therefore their dual vertices in $\Wscr$ (which should be topologically distinct) would turn out t be coincident. Furthermore, an edge of $\T$ common to those two triangles would correspond to an edge of $\Wscr$ of length zero. Clearly, $\Wscr$ is not a subdivision according to our definitions. It may be possible to extend or modify our definitions to accomodate dummy edges and coincident vertices, but it seems likely that $\Wscr$ will always be a concept harder to understand and more awkward to manipulate than $\T$. This is yet another argument for working with the Delaunay diagram instead of the Voronoi.} \footprob{\Prob{\Sbsec{\Sec2.0}.1}}{(a) Can we view any Delaunay triangulation as a correct Delaunay diagram computed for a perturbed version of the input sites? More precisely, can we find, for every Delaunay triangulation $\T$ of a set $P$ of $n$ sites, another set $P\pr$ of $n$ sites whose Delaunay diagram $\D{P\pr}$ is isomorphic to $\T$? If so, try to give useful upper bounds for the ``perturbation'' $\dist(P, P\pr)$. If not, characterize the triangulations that can be Delaunay diagrams.\lp {\it Ans.}: I dont know, but this may be easier to think after $\lambda$-mapping. Seems this question is equivalent to asking if every triangulation of the faces of a convex polyhedron circumscribed in a sphere can be realized by displacing its vertices on the surface of that sphere.} \subsection{\Ch8.\Sbsec{\Sec2.0}. Characterization} The following lemma characterizes Delaunay triangulations: \theorem{\Th{15}}{A triangulation $\T$ of a set $P$ is a Delaunay triangulation of $P$ if and only if every interior face of $\T$ is inscribed in a $P$-free circle.} \proof{The ``only if'' part is immediate from lemma \Th{13B}, so we only have to prove the converse. So, assume every interior face of $\T$ is inscribed in a $P$-free circle. Let $C$ be a $P$-free circle circumscribed on any interior face (triangle) $F$ of $\T$, and let $K$ be the set of all sites in $P$ that lie on the boundary of $C$. Since $K$ has three or more points, it determines an interior Delaunay face $\DP K$; since the vertices of $F$ are in $K$, the triangle $F$ is contained in $\DP K$, and the edges of $F$ are either edges or diagonals of $\DP K$. The exterior face of $\T$ is the complement of $\hull{P}$, and therefore is exactly the exterior face of $\DP$. We conclude $\T$ is obtained from $\DP$ by the addition of diagonal edges so as to subdivide its faces into triangles, that is, $\T$ is a Delaunay triangulation of $P$.} The quadratic mapping $\lambda$ helps in the visualization of this concept, too. It is easy to see that a we can obtain a Delaunay triangulation by projecting on the $xy$-plane the underside of $\hull{\lambda(P)}$, after decomposing every face of the latter into triangles by means of non-intersecting diagonals.\figspace4cm \subsection{\Ch8.\Sbsec{\Sec2.1}. The angle property} Let $P$ consist of four sites $a,b,c,d$, defining in this order the vertices of a convex quadrilateral $Q$. This set admits of at most two triangulations, consisting of the sides $ab$, $bc$, $cd$, and $da$ of $Q$, and exactly one of the diagonals, $ac$ or $bd$:\figspace3cm Let us consider the conditions for the first triangulation above to be a Delaunay one. By theorem \Th{15}, the circles $abc$ and $cda$ must be $P$-free, that is, $\incircle(a,b,c,d)$ and $\incircle(c,d,a,b)$ must both be false. From the properties of $\incircle$ we know these two are equivalent. Furthermore, by lemma \Th{20C}, $\incircle (a,b,c,d)$ is false if and only if $\hat a+\hat c \geq 180\deg \geq \hat b+\hat d$. Similarly, for the second triangulation above to be Delaunay, $\incircle(b,c,d,a)$ and $\incircle(d,a,b,c)$ must be both false, which is equivalent by lemma \Th{20C} to $\hat b+\hat d \geq 180\deg\geq \hat a+\hat c$. We therefore conclude that {\it a Delaunay triangulation of the vertices of $Q$ must use the diagonal of $Q$ connecting the two opposite internal angles with largest sum}. If $Q$ is simple but not convex, the only possible triangulation of its vertices uses all four sides and {\it both} diagonals of $Q$. So, the above rule is trivially true even in this case.\figspace3cm This rule can be easily extended to Delaunay triangulations of more than four vertices: \definition{\Def{11}}{Let $prq$ and $qrs$ be two adjacent interior faces of an arbitrary triangulation. The vertices $prqs$ (in that order) constitute the {\it edge quadrilateral} of the edge $pq$. We say that $e$ satisfies the {\it angle property} if it connects the opposite interior angles of $Q$ having the largest sum, that is, if $\hat p+\hat q \geq \hat s+\hat r$.\figspace3cm The {\it strict} angle property is similarly defined, with $>$ instead of $\geq$.} From the definition it folows immediately that an edge quadrilateral $Q=prqs$ is a simple polygon, and the diagonal $pq$ is interior to $Q$. Note however that $Q$ need not be convex, so the other diagonal $rs$ may lie outside $Q$ or on its boundary. The following conditions are equivalent to the angle property: \indent 1. $\hat p+\hat q \geq 180\deg$; \indent 2. $\hat s+\hat r \leq 180\deg$; \indent 3. $\lnot \incircle(p,r,q,s)$; \indent 4. $\lnot \incircle(q,s,p,r)$; \indent 5. $\lnot \outcircle(r,q,s,p)$; \indent 6. $\lnot \outcircle(s,p,r,q)$. The following result is an immediate consequence of theorem \Th{15}: \lemma{\Th{21A}}{Let $e$ be any interior edge of an arbitrary triangulation $\T$ of $P$. If $e$ does not have the angle property, then no Delaunay triangulation of $P$ uses $e$. Furthermore, if $e$ does not have the strict angle property, then $e$ is not in the Delaunay diagram of $P$.} So, the failure of an edge of $\T$ to satisfy the angle condition is witness to the ``non-Delaunayhood'' of that edge. The converse is not true, however: an edge of $\T$ may pass the angle test, and yet there may be no Delaunay triangulation of $P$ containing $e$. This happens because the angle test is a very local one: it takes into account only the four sites in the edge quadrialteral of $e$, whereas to establish the ``Delaunayhood'' of $e$ we must consider {\em all} sites of $P$. The figure below shows a counterexample:\figspace3cm In view of this, the following result may seem surprising: \lemma{\Th{19}}{A triangulation $\T$ of $P$ is a Delaunay triangulation if and only if every interior edge of $\T$ has the angle property.} \proof{The ``only if'' part is immediate from theorem \Th{21A}, so let us prove the ``if'' part. Let $\T$ be a triangulation of $P$ whose edges all pass the angle test. If $\T$ is not a Delaunay triangulation, it must have at least one interior face $F$ whose circumcircle $C$ strictly contains some site $r$ of $P$. Let the $psq$ (in counterclockwise order) be the vertices of $F$, so that the $pq$ is the edge nearest to $r$.\figspace3cm \pp Now assume $F$ and $r$ have been chosen, among all possible such pairs, so as to maximize the angle $qrp$. We claim that $prq$ is an internal face of the triangulation $\T$, and indeed it is the right face of the edge $pq$. To see why, suppose the right face is $F\pr=puq$ with $u\neq q$. The site $u$ cannot be inside the circle $prq$, for then $\angle puq>\angle prq$. It follows that either $pu$ intersects $rq$, or $uq$ intersects $pr$. \pp In the first case we have $\angle pru > \angle prq$. Furthermore, we have $\outcircle(p,r,q,u)$, and threfore $\incircle(u,q,p,r)$, that is, $r$ is inside the circumcircle of $F\pr$. But then $F\pr$ and $r$ would have been a better choice than $F$ and $r$, a contradiction. A similar reasoning holds for the second case, so $u=r$ and $prq$ is the right face of $pq$. \pp So, we conclude that if $\T$ is not a Delaunay triangulation, there is some edge quadrilateral $prqs$ such that $r$ is inside the circumcircle of $qsp$, that is, $\incircle(q,s,p,r)=\incircle(p,r,q,s)$ is true. But this contadicts the assumption that all edges (in particular $pq$) satisfy the angle test.} In conclusion, the angle property is only a necessary conditition for the ``Delaunayhood'' of a single edge, but is both necessary and sufficient in determining the ``Delaunayhood'' of the whole triangulation. \vfill\eject \footprob{\Prob{\Sbsec{\Sec2.1}.1}}{Assume you are given a pointer to a quad-edge data structure in which the {\fx Data} field of each primal directed edge points to a record containing the (finite) coordinates of a vertex, ostensibly the origin of that edge. (a) Write an algorithm to test whether the data structure is a consistent representation of some straight-line subdivision of the plane. You can assume the structure is consistent with some edge algebra, but nothing can be said {\it a priori} about its genus and the consistency of the {\fx Data} links. (b) Write another algorithm, to be run after (a), that tests whether the subdivision is a triangulation (with convex perimeter). Try to get $O(n)$ time and $O(n)$ working storage for a subdivision with $n$ edges. Is it possible to do it in $O(n)$ time and $O(1)$ storage?} \footprob{\Prob{\Sbsec{\Sec2.1}.2}}{Write an algorithm that tests whether a given triangulation with convex outer boundary is a Delaunay triangulation. Assume the triangulation is given by quad-edge structure that passes the tests of problem \Prob{\Sbsec{\Sec2.1}.1}. {\it Hint\/}: use lemma \Th{19}.} \footprob{\Prob{\Sbsec{\Sec2.1}.3}}{Let us define a {\it $C$-subdivision} as a convex subdivision of the plane with convex outer boundary and such that all vertices of each interior face $F$ lie on a common circle $CF$. Clearly, every Delaunay diagram is a $C$-subdivison. (a) Exhibit a $C$-subdivision that is not a Delaunay diagram. (b) Write an algorithm that tests whether a quad-edge data structure represents a $C$-subdivision. (c) Extend the angle property to $C$-subdivisions, as follows: we say that an edge $e=pq$ of a $C$-subdivision, with left and right faces $F$ and $G$, satisfies the strict angle property if the internal angles at $p$ and $q$ of the polygon $F\uni G$ add up to more than $180\deg$. Prove that this is the same as saying that no vertices of $F$ lie on or inside $CG$, or equivalently that no vertices of $G$ lie on or inside $CF$. (d) Prove that a $C$-subdivision is a Delaunay diagram if and only if every edge in it satisfies the strict angle property. (e) Write an algorithm that tests whether a $C$-subdivision a Delaunay diagram.} \footprob{\Prob{\Sbsec{\Sec2.1}.4}}{Prove the following lemma: a triangulation of $P$ is Delaunay if and only if every edge is the chord of a $P$-fre circle.} \footprob{\Prob{\Sbsec{\Sec2.1}.5}}{Prove the following lemma: a triangulation of $P$ is Delaunay if and only if every edge $e$ with edge quadrilateral $prqs$ is the chord of an $\set{r,s}$-free circle.} \footprob{\Prob{\Sbsec{\Sec2.1}.6}}{Prove that (a) the only $C$-subdivision of a set $P$ with minimal set of edges is $\DP$, and (b) every $C$-subdivision of $P$ with maximal set of edges is a Delaunay triangulation of $P$.} \subsection{\Ch8.\Sbsec{\Sec2.2}. The edge swapping rule.} If $e$ does {\it not} have the angle property, then properties 1--6 are all false. In particular, $\hat p$ and $\hat q$ are both less that $180\deg$, which means the diagonal $rs$ of the quadrilateral $prqs$ lies entirely inside it. Also, $\hat r+\hat s>180$, which implies the edge $rs$ now satisfies the (strict) angle property. This fact will be referred to as the {\it edge swapping rule}: in an arbitrary triangulation $\T$, replacement of an interior edge $e=pq$ that does not have the angle property by the opposite diagonal $rs$ of its edge quadrilateral $prqs$ will produce another triangulation $\T\pr$, and this new edge will have the angle property. \note{Make the convention that exterior edges always satisfy the angle property, so ``interior edge'' can be replaced by ``edge'' almost everywhere.} It is instructive to consider the meaning of the angle test and the swapping rule in terms of the quadratic mapping $\lambda$. Any triangulation of $P$ is the projection on the $xy$ plane of a $z$-monotone polyhedral surface $\L$ with triangular faces whose projection covers $\hull(P)$, and whose vertices are the points in $\lambda(P)$. An edge $pq$ with quadrilateral $prqs$ has the strict angle property if and only if $s$ is outside the (clockwise oriented) circle $prq$, that is, $\lambda(s)$ is above the plane $\lambda(p)\lambda(r)\lambda(q)$. This is also equivalent to saying that the line $\lambda(r)\lambda(s)$ passes above the line $\lambda(p)\lambda(q)$; that is, the two triangles $\lambda(p)\lambda(r)\lambda(q)$ and $\lambda(q)\lambda(s)\lambda(p)$ are ``folded upwards'' at the diagonal $\lambda(p)\lambda(q)$.\figspace3cm The weak angle property is of course similar, but the two triangles are allowed to be coplanar. The swapping rule then can be interpreted as saying that if a pair of adjacent triangles $\lambda(p)\lambda(r)\lambda(q)$ and $\lambda(q)\lambda(s)\lambda(p)$ of $T\ast$ are ``folded downwards'' around their common edge, we can replace them by the triangles $\lambda(s)\lambda(p)\lambda(r)$ and $\lambda(r)\lambda(q)\lambda(s)$. The resulting polyhedral surface will still be $z$-monotone (thanks to the special properties of the paraboloid $\Lambda$), and the two new faces will obviously be folded upwards. These facts correspond to our previous conclusion that the swap of $pq$ does not introduce crossing edges, and that the new edge $rs$ satisfies the strict angle property. The two triangles that were swapped constitute a ``notch'' in the infinite polyhedron consisting of all points above the surface $\L$. The swapping rule then allows us to fill in that notch with the tetrahedron $\lambda(p)\lambda(q)\lambda(r)\lambda(s)$. Indeed, the angle property being false is equivalent to $$\incircle(p,r,q,s)\;\eqv\;\NEG(\lambda(p),\lambda(r),\lambda(q),\lambda(s))\;\eqv\;\Delta(\lambda(p),\lambda(r),\lambda(q),\lambda(s))< 0,\eqno(\Eq{67X})$$ As we recall, the signed volume of the tetrahedron $\lambda(p)\lambda(r)\lambda(q)\lambda(s)$ is ${1\over 6}\Delta(\lambda(p),\lambda(r),\lambda(q),\lambda(s))$, so condition (\Eq{67X}) can be interpreted as meaning that the tetrahedron $\lambda(p)\lambda(r)\lambda(q)\lambda(s)$ is a ``notch'' in $\L$. \footprob{\Prob{\Sbsec{\Sec2.2}.1}}{Give the code (in terms of {\fx Splice} and other quad-edge functions) for the edge swap operation.} \subsection{\Ch8.\Sbsec{\Sec2.3}. Angle signature.} \note{Define angle signature of a triangulation, and prove that the edge swaps decrease it. Use then lemma \Th{19} to prove that a triangulation has minimal angle signature iff it is Delaunay} \note{Somewhere: caveat on sets $P$ lying on a single line (the exterior face is the complement of a degenerate convex polygon (a line segment), and there are no interior faces).} \section{\Ch8.\Sec{2A}. Applications of Voronoi and Delaunay diagrams} % Begin Let's postpone for the moment the question of how to construct the Voronoi diagram, and let's discuss how it can be used to solve the post-office problem and its relatives. \subsection{\Ch8.\Sbsec{\Sec{2A}.1}. The post-office problem} {\bf JUNK:} To find the site $p$ that is closest to a given $x\in \Sscr$, all we have to do is to locate the region $\V\set{p}$ that contains $x$. The Lee-Preparata algorithm (??? 5.???) can be used for this purpose, since $\VP$ is a convex (and therefore monotone) subdivision of the plane. Its query time is not optimal ($O(\log n)^2$), but is significantly beter than that of exhaustive search. Kirkpatrick's algorithm (5.???) has a better asymptotic performance, just $O(\log n)$ time per query. As it is formulated, Kirkpatrick's algorithm requires the Voronoi regions to be subdivided even further, so as to obtain a triangulation. Two possible ways to implement this subdivision are shown below\foot1{Note that these subdivisions are {\it not} related to the Delaunay triangulation $\TP$ in any simple way: the edges here connect Voronoi vertices, not data points.}. \figspace 4cm Given $\VP$, we can construct either of these triangulations (and the data structure required by algorithm 5.???) in linear time, so the preprocessing cost is still dominated by the computation of the Voronoi diagram. {\bf ??? Is there a location algorithm that does not require the full Voronoi, only the Delaunay triang.?} \subsection{\Ch8.\Sbsec{\Sec{2A}.2}. Other closest-point problems} {\bf JUNK:} A common variant of the post-office problem is the {\it nearest neighbor problem}: for each point $p\in P$, find the point $q\in P$ that is closest to it. The straightforward solution, computing all $n-1$ distances from each $p$ and selecting the minimum, gives a $\Theta(n^2)$ algorithm. If we were to construct the Voronoi diagram for $P-\set{p}$ and locate $p$ in it, the total work would be even worse, $\Theta(n^2\log n)$ in the worst case. However, it is the case that $\VP$ already contains enough information to allow us to solve the nearest-neighbor problem in linear time, as shown by the following result: \lemma{8.11}{The point $q\in P$ that is closest to a given point $p\in P$ is adjacent to $p$ in the Delaunay diagram for $P$, and the midpoint of $pq$ is on the Voronoi edge $\V\set{pq}$.} \proof{Let $C$ be the circle with diameter $pq$, and $C\pr$ be the one with center $p$ and radius $pq$ \figspace 3cm Since $q$ is the closest neighbor of $p$, there are no points of $P$ interior to $C\pr$; therefore, $p$ and $q$ are the only points of $P$ on $C$. It follows that the center $o$ of $C$, which is the midpoint of $pq$, is a point of the edge $e{pq}$ and is not a Voronoi vertex, so $p$ and $q$ are adjacent in $\DP$.} Therefore, to find the closest neighbor of $p$ we have to look only at the vertices $q$ adjacent to it in $\DP$, and select the closest one. This may take $\Theta(n)$ time for a single vertex $p$, since its degree can be as high as $n-1$. The total number of edges, however, is at most $3n-6$, so the computation of all nearest neighbors can be done in only $O(n)$ time. Another way to look at this result is to observe that, although the {\it maximum} degree of a point in $\DP$ is $n-1$, the {\it average} degree is $O(1)$ (more precisely, $6\left(1-{2\over n}\right)$). The computation of the Voronoi, requiring $O(n\log n)$ time, will be the limiting step in this algorithm. A slight variation of this algorithm also solves the {\it closest pair problem}, namely finding the pair $p,q\in P$ whose distance is minimal. \note{Closest $X$--$Y$ pair, MST} \subsection{\Ch8.\Sbsec{\Sec{2A}.3}. Circle problems} {\bf JUNK:} Another application of the Voronoi diagram is in finding the largest $P$-free circle $C$. If no constraints are imposed on $C$, the problem has no solution: we can always place a circle of arbitrary radius far enough from $P$ so that it includes no points of $P$. If the center of $C$ is required to lie in the covex hull of $P$, the problem becomes meaningful; using theorem 8.???, we can easily prove that the center of $C$ is either a Voronoi vertex interior to the hull of $P$, or the intersection of a Voronoi edge with an edge of the hull. It turns out that all those points, and therefore the largest $C$, can be found in $O(n)$ time {\bf I guess...}. \subsection{\Ch8.\Sbsec{\Sec{2A}.2}. Good triangulations} Many applications, which at first seem unrelated to the post-office problem, require the construction of triangulations having as vertices a specified set of points. Two examples are interpolation of two-argument functions (surface fitting) on irregularly distributed data points, and finite-element alnalysis of two-dimensional structures. In those problems, the numerical accuracy and stability of the result are strongly dependent on the triangulation used. As a general rule, faces that are more-or-less equilateral are preferrable to long and thin ones. If we accept the angle signature of the triangulation as a formal measure of its quality, we conclude that a triangulation of $P$ is optimal if and only if it is a Delaunay one. \digress{It was once thought that a Delaunay triangulation had minimum total edge length among all triangulations with same vertices. Several proofs of this ``fact'' were published before a counter-example was found.} \section{\Ch8.\Sec3. Construction of Voronoi and Delaunay diagrams} % Begin Since the Delaunay diagram is topologically equivalent to the underside of the three-dimensional convex hull of $\lambda(P)$, any algorithm that computes the latter can be used to construct the former. The computation of the images $\lambda(P)$ takes only $O(n)$ time (more precisely, $2n$ multiplications and $n$ additions), so we can construct the Delaunay asymptotically as fast as we can compute three-dimensional hulls. In particular, the algorithm of Preparata and Hong described in chapter \Ch4 allows us to compute the Delaunay diagram in $O(n\log n)$ time. From the Delaunay diagram (in most reasonable representations) we can obtain the convex hull of $P$ in $O(n)$ further time, by walking around the exterior face. This gives immediately a $\Omega(n\log n)$ lower bound for contruction of $\DP$, so the algorithm derived from the Preparata-Hong one is asymptotically optimal in the worst case. Actually, this reduction does not require the full power of three-dimensional convex hull algorithms, for we need only the faces and edges that are visible from below, and we know {\it a priori} that all points in $\lambda(P)$ will be vertices of that part of the hull. With these assumptions the algorithms can be considerably simplified. \note{plug here the discussion about convex hull algorithms being usable for Delaunay construction (from section \Sec1)} \note{Jarvis's algorithm (or rather its $k$-dimensional generalization, the ``gift-wrapping algorithm) corresponds to an algorithm published by McLain. Check also the Voronoi algorithms of Watson, Brostow, Bowyer, Brassel and Reif, Greene and Sibson, etc., and look at them as convex-hull algorithms.} \note{The merge step in the $O(n\log n)$ Voronoi algorithm of Kirkpatrick, based on random splitting, under this correspondece should give an $O(n)$ hull-merging algorithm for convex polyhedra (with the assumption that all vertices lie on the paraboloid, so none has to be discarded). By dualization, we should get an $O(n)$ intersection algorithm for convex polyhedra (with similar assumptions). Investigate!} {\bf JUNK:} To compute the edges and vertices of a Voronoi region $Vp$, we can start with $Vp\gets \Sscr$ and then perform $Vp\gets Vp\inte C(p,q)$ for all $q$ in $P-\set{p}$. After each step, $Vp$ will be essentially a convex polygon, except that some of its vertices may be at infinity. We can represent such a region by the list of its vertices, in order of polar angle around the point $p$, and use the algorithms of chapter 4 to compute each intersection. With this data structure, each intersection can be found in $O(\log n)$ time, so in $O(n\log n)$ we can compute a single region, and $O(n^2\log n)$ the whole diagram. In the worst case, the running time may be actually $\Theta(n^2\log n)$. We clearly need better algorithms. \subsection{\Ch8.\Sbsec{\Sec3.1}. The edge swapping algorithm} Theorem \Th{19} and the edge swapping rule suggest an algorithm for transforming an arbitrary triangulation $\T$ of $P$ into a Delaunay one: if some edge of $\T$ does not satisfy the angle property, then swap it, and continue repeating this process until the angle test is satisfied by all edges. More precisely, \algorithm{\Alg0}{Transforms a triangulation into a Delaunay one by repeated edge swaps.} \algbody{ \comm{The input is a triangulation $\T$ of the $n$ sites. On output, $\T$ will be a Delaunay triangulation. Uses a set $S$ to store the ``dubious'' edges, the ones which as far as we know may not satisfy the angle property.} \step{\Step1}{Let $S$ be the set of all edges of $\T$.} \step{\Step2}{While $S\neq\empty$, do} \begblock \step{\Step3}{Remove some edge $e=pq$ from $S$.\lp If $e$ is an interior edge, and does not satisfy the angle property, then} \begblock \step{\Step4}{Let $prqs$ be the edge quadrilateral of $e$.\lp Swap $e$, and set $S\gets S\uni \set{pr, rq, qs, st}$} \endblock \unstep{Repeat from step 2.} \endblock } The correctness of algorithm \Alg0 is easily established. The key step is showing that every time we get to step \Step2, the current subdivision $\T$ is indeed a triangulation of $P$, and any interior edge of the triangulation that is not in $S$ satisfies the angle property. The first time we get to step \Step2 this is trivially true. As for the other times, note first that the edge $e$ can be deleted from $S$ without affecting the invariant, whether step \Step4 is executed or or not. Observe also that \Step4 replaces the triangles $prq$ and $qsp$ by $prs$ and $sqr$, but does not affect any other triangles. So, the only edges that can have their angle property destroyed by step \Step4 are the ones on those triangles, i.e. $pr$, $rq$, $qs$, and $sp$. The swapping rule tells us that the new edge $rs$ satisfies the angle test, and the resulting graph is still a proper triangulation of $P$. When the algorithm terminates we have $S=\empty$, and then from this invariant we conclude all edges satisfy the angle property. At any time, the number of edges in the triangulation is at most $3n-6$, so step \Step4 must be executed at least once every $3n-6$ iterations. The angle signature of the triangulation is strictly reduced by each edge swap. Since there is only a finite number of triangulations of $P$, we conclude algorithm \Alg0 eventally terminates. \note{Discuss performance of this algorithm. Should we prove that any pair of vertices is conneced/disconnected at most once, so the total number of edge swaps is $O(n^2)$, and so is the running time of the algorithm? If not, leave as exercise?} \note{Explore the concept of {\it edge signature:} for each edge $e=pq$, the sum of the angles $\hat p$ and $\hat q$ of its edge quadrilateral. Does this help in showing $O(n^2)$ running time?} \note{Mumble some: it doesn't matter whether we distinguish or not the edge $pq$ from $qp$. Does it matter much if $S$ is a bag rather than a set?} \subsection{\Ch8.\Sbsec{\Sec3.2}. Adding sites in sorted order} \note{Maybe this should go in the subdivision chapter} To apply algorithm \Alg0, we need first construct an arbitrary triangulation of the set $P$. One possible method is to start from a triangulation of minimal size (say, two sites connected by an edge), and connect to it the remaining sites, one at a time, in such a way that after each insertion the current graph is a triangulation. This process is easier to implement if the sites are inserted in order of increasing $x$-coordinate (with ties broken by $y$-coordinate). Let us denote by $P\pr$ the set of sites inserted so far, and by $\T\pr$ the current triangulation. The next site $s$ will always lie outside the convex hull of $P\pr$, that is, in the exterior face of $\T\pr$. We can build a triangulation $\T\prr$ for the set $P\prr=P\pr\uni \set{s}$ by simply connecting $s$ to all sites on the exterior face of $\T\pr$ that are visible from $s$.\figspace4cm These sites are precisely those on the outer boundary of $\T\pr$ that lie between the common tangents of $s$ and $\T\pr$, on the side nearest to $s$. The sorting of the input points is helpful here, too: it is easy to see that the site $t$ that was last inserted in $\T\pr$ is visible from $s$. Therefore, we have only to start at $t$ and walk around the outer boundary of $\T\pr$ in both directions, connecting each site we encounter to $s$, until the connecting edge becomes tangent to $\T\pr$. The insertion of $s$ then takes $O(k)$ time, where $k$ is the number of edges connecting $s$ to $\T\pr$. This can be as high as $i-1$ for the $i$-th vertex. It would seem that the total time for inserting $n$ vertices would be quadratic in $n$, but this is not so. The final triangulation $\T$ has only $O(n)$ edges, and no edges are deleted in the process; we conclude that the $n$ insertions take only $O(n)$ time. Therefore, the time spent in building $\T$ is dominated by initial sorting of the sites, which may require $\Theta(n\log n)$ time in the worst case. Since from any triangulation we can get the convex hull of $P$ in linear time, this time is optimal. Anyway, the cost of building $\T$ is overwhelmed by that of transforming it into a Delaunay trinagulation by algorithm \Alg0. We have shown an $O(n^2)$ upper bound for this process, and that bound is attained for the class of inputs illustrated below:\note{References!} \figspace4cm Later on we will see that a Delaunay triangulation for $P$ can be constructed in $O(n\log n)$ time, so the only practical advantage of the method outlined above is its relative simplicity. In fact, the above method can be further simplified (without exceeding the $O(n^2)$ total time bound) if we intersperse the edge-swapping process with the construction of $\T$. The idea is to apply the edge-swapping algorithm \Alg0 to each intermediate triangulation, so that at every stage we have a {\it Delaunay} triangulation for the sites added so far. More precisely, at each step we have a {\em Delaunay} triangulation $\T\pr$ for the set $P\pr$ of the sites added so far, and we want to construct from it a Delaunay triangulation $\T\prr$ for the set $P\prr=P\pr\uni\set{s}$, where $s$ is the next site in the sorted list. We begin by connecting $s$ to $\T\pr$ as explained above, resulting in a triangulation $\T\star$ of $P\prr$ that is not necessarily a Delaunay one. This triangulation consists of all sites, edges and interior faces of $\T\pr$, plus zero or more new faces $sp0p1, sp1p2, \ldots, sp{m-1}pm$, where $p0, p1, \ldots, pm$ are the neighbors of $s$, in counterclockwise order around $s$. All edges in the Delaunay triangulation $\T\pr$ were known to satisfy the angle condition. What can we say about $\T\star$? are all its edges suspect? Clearly not: the only edges for which the condition may have become false are those adjacent to the new triangles, namely $p0p1, p1p2, \ldots, p{m-1}pm$. The new edges $sp0, sp1, \ldots, spm$ will always satisfy the angle condition. To see why, observe that for each $pi$ there is a supporting line $li$ of $\T\pr$ passing through $pi$ and strictly to the left of $s$. Therefore, there is a $P\prr$-free circle $Ci$, tangent to $li$, that passes through $s$ and $pi$ and through no other site of $P\prr$. We conclude $spi$ is an edge of the Delaunay diagram of $P\prr$, and (by lemma \Th{21}) it satisfies the angle condition in {\it any} triangulation of $P\prr$. So, in step \Step1 of algorithm \Alg0 we can initialize $S$ to $\set{p0p1, p1p2, \ldots, p{m-1}pm}$ only, and proceed to the main loop (steps \Step2--\Step4). At this point, the following statements about the current triangulation are obviously true: \begitems \itemno{1.}Every edge $e=pq$ in $S$ has $pqs$ as its left face. \itemno{2.}Every edge that was not in $\T\pr$ is incident to $s$, and is an edge of the Delaunay diagram of $P\prr$. \enditems Now consider what happens on the execution of steps \Step3--\Step4. If the edge $e=pq$ lies on the exterior face, or satisfies the angle condition, it is simply ignored, and we return to step 2. Otherwise $e$ is deleted from the triangulation, and replaced by the new edge $sr$.\figspace3cm At this point algorithm \Alg0 would add to the bag $S$ the four edges $pr$, $rq$, $qs$ and $sp$, on the grounds that the swapping of $e$ may have falsified their angle properties. However, by statement 2 above and lemma \Th{21} we know that $qs$ and $sp$ still are ``good'', so we only have to stack $pr$ and $rq$. Now observe that statement 1 above is still true, since the stacked edges $pr$ and $rq$ have left faces $prs$ and $rqs$. We will show now that statement 2 too is still true. The only non-obvious part is to show that the new edge $rs$ is a Delaunay edge for $P\prr$. The right face $prq$ of the deleted edge $e$ was an internal face of the triangulation. The vertices $p$, $q$, and $r$ are all distinct from $s$. By statement 2 the edges $pr$, $rq$ and $qp$ were all original edges of $\T\pr$, and therefore $prq$ was one of its interior faces. So, its circumcircle $C$ was $P\pr$-free. On the other hand, $C$ contains $s$, for $e$ failed the angle test. Then consider the circle $Cr$ passing through $s$ and tangent internally to $C$ at $r$: this circle is $P\prr$-free, and indeed passes through no other site of $P\prr$ other than $r$ and $s$. This proves $rs$ is an edge of the Delaunay diagram of $P\prr$.\figspace3cm Therefore, when we return to step 2 statements 1 and 2 above are still true. An inductive argument shows that the above reasoning applies to all iterations of steps \Step2--\Step4. \note{Plan 1: discuss now that these observations also hold if $s$ is not in order, except that the construction of $\T\star$ from $\T\pr$ is more complex. Then give algorithm \Alg1 for a general $s$, leaving insertion out.} \note{Plan 2: present algorithm \Alg1 for the particular case of $s$ in sorted $x$-order, and discuss general case later. \note{Adopted this one.}} These simplifications are implemented in algorithm \Alg1 below: \algorithm{\Alg1}{Incremental construction of a Delaunay triangulation.} \algbody{ \comm{The input is a Delaunay triangulation $\T$ for the sites in $P\pr$, and a site $s$ not in $P\pr$ such that $s$ follows all vertices of $P\pr$ in $x$-order (disambiguated by $y$). On output, $\T$ will have been transformed into a Delaunay triangulation of $P\prr=P\pr\uni\set{s}$. The set $S$ keeps track of the ``dubious'' edges, as in algorithm \Alg0.} \step{\Step0}{Connect $s$ to $\T$, transforming it into a triangulation of $P\prr=P\pr\uni\set{s}$.} \step{\Step1}{Let $sp0p1, sp1p2, \ldots, sp{m-1}pm$ be the interior faces of $T$ incident to $s$, in counterclockwise order. Set $S\gets\set{p0p1, p1p2, \ldots, p{m-1}pm}$.} \step{\Step2}{While $S\neq\empty$, do} \begblock \step{\Step3}{Remove some edge $e=pq$ from $S$. If $e$ is an interior edge, and does not satisfy the angle property, then} \begblock \step{\Step4}{Let $prqs$ be the edge quadrilateral of $e$. Swap $e$, and set $S\gets S\uni \set{pr, rq}$} \endblock \endblock } We can siplify this algorithm even further. If the set $S$ is implemented by a stack, and in steps \Step1 and \Step4 the edges are pushed into $S$ in the given order, then the contents of $S$ at any time, from the bottom up, is a path $x0,x1, x2, \ldots, xk$ starting from the site $x0=p0$. Furthermore, by condition 1 the left face of each edge $xi x{i+1}$ is the triangle $xi x{i+1} s$. These triangles constitute a ``fan-shaped'' region whose ``hub'' is the site $s$. The effect of steps \Step3 and \Step4 on this region is particularly simple: at each iteration, the last triangle is either removed from the region (i.e., the top edge $e$ of $S$ is ignored), or it is replaced by two triangles ($e$ is swapped and two edges are pushed into $S$).\figspace4cm Note that every edge $xi x{i+1}$ on the stack, except the first one, is folowed in counterclockwise order around $xi$ by the edges $xi s$ and $xi x{i-1}$. The stack $S$ is then superfluous: given its top edge $e$, we can get the one just below it by computing $e\onext\onext\sym$. \footprob{\Prob{\Sbsec{\Sec3.1}.1}}{Code algorithm \Alg1 without using the stack $S$.} The total number of iterations of steps \Step2--\Step4 is easily determined. If in step \Step0 we introduce $m$ new edges connecting the site $s$ to the rest of $\T$, then the initial size of $S$ is $m-2$. Every time $s$ gains a new incident edge in step \Step4, two other edges are pushed onto $S$. Every iteration of steps \Step2--\Step4 deletes exactly one edge from $S$. Therefore, if the final degree of $s$ in $\T$ is $d$, then the total number of iterations is $2d-m-2$. If $P$ contains $n$ sites, the number of edges in the triangulation $\T$ is $O(n)$, and so are the degree of $s$ and the running time of algorithm \Alg1. \subsection{\Ch8.\Sbsec{\Sec3.3}. Adding sites in random order} Algorithm \Alg1 is essentially a procedure for updating a Delaunay triangulation to account for the insertion of a new site $s$. However, we based its development on the assumption that the sites were going to be inserted in order of increasing $x$ coordinate. What modifications in algorithm \Alg1 are required to remove this restriction? The somewhat surprising answer is that only step \Step0 has to be modified. If the sites are not inserted in increasing $x$-order, the new site $s$ may lie anywhere on the plane. We then have to distinguish three main cases: \itemno{1.}$s$ lies in the exterior face of $\T$. This is essentially the situation we had in algorithm \Alg1: we connect $s$ to one or more consecutive sites $p0, p1, \ldots, pm$ on the outer boundary of $\T$. The edges $sp0$ and $spm$ are the tangents from $s$ to the outer boundary of $\T$. \itemno{2.}$s$ lies in an interior face of $\T$. We connect $s$ to the three vertices $p0, p1, p2, p3=p0$ of that face.\figspace3cm \itemno{3.}$s$ lies on some edge $e$ of $\T$. We first delete $e$, and connect $s$ to its two endpoints. The edge $e$ may belong to zero, one, or two internal faces of $\T$. The third site in each of those faces must also be connected to $s$.\figspace3cm \digress{There is of course a fourth case, namely $s$ coincides with some other site $u$ in $P\pr$. The Voronoi diagram of $P$ was defined under the assumption that $P$ was a {\it set}, that is, the sites were all distinct. This assumption may not be valid in typical applications of the incremental method, so we have to provide for this case, too. Obviously, if $s$ is already in $P\pr$, its insertion does not change $P\pr$ nor its Delaunay triangulations, so in this case we can skip the rest of the algorithm\foot\dag{For some applications, we may have to keep track of the multiplicity of each site in $P$, that is, how many times it was inserted. Furthermore, the sites usually contain other information besides their position, so \note{etc}}} In all three cases above, the insertion of $s$ creates a star-shaped or ``fan-shaped'' region consisting of zero or more new internal faces $sp0p1, sp1p2, \ldots, sp{m-1}pm$ (where $pm$ may or may not be $p0$). Every internal face of of the new triangulation $\T\star$ that is not one of these was present also in the initial Delaunay triangulation $\T\pr$. The same arguments used above show that $spi$ is an edge of the Delaunay diagram of $P\prr$, and therefore will satisfy the angle property in any triangulation of $P\prr$. Therefore, the only edges of $\T\pr$ that may have had their angle property falsified are $p0p1, p1p2, \ldots, p{m-1}pm$. The rest of the correctness proof of algorithm \Alg1 then carries through unchanged. The analysis of the running time requires some revision. It is still true that steps \Step1--\Step4 take only $O(d)$ time, where $d=O(n)$ is the degree of $s$ in the final triangulation $\T\prr$. However, step \Step0 now has the additional burden of locating $s$ in $\T\pr$. Moreover, if $s$ is found to lie in the exterior face, the determination of the sites visible from it is not as simple as before. We know it is possible to locate $s$ in $O(\log n)$ time, given a suitable search data structure like that of Kirkpatrick. Unfortunately, Kirkpatrick's search structure takes $\Theta(n)$ time to build from scratch, and we know of no simple algorithm for updating it, after each call to algorithm \Alg1, in less than $\Theta(n)$ worst-case time. The same applies to all other efficient point location techniques, so there is no much justification for using them in algorithm \Alg1. The exhaustive enumeration of all interior faces of $\T\pr$ takes $O(n)$ time too, so it is a reasonable solution. In view of the above, when $s$ falls in the exterior face of $\T\pr$ we find all sites visible from it by a a simple sequential traversal of the outer boundary of $\T\pr$. Step \Step0 (and the entire algorithm) will then take $O(n)$ time in the worst case. In a sense, this time is optimal: as we remarked before, the $d$ new edges added by algorithm \Alg1 must be in any Delaunay triangulation of $P\prr$, and $d$ may be $\Theta(n)$ in the worst case, so any algorithm for adding a site to a Delaunay triangulation will take that much time in the worst case. In fact, if we had a search data structure for $\T\pr$ that could be updated in $O(\log n + d)$ time to a search data structure for $\T\prr$, and allowed us to do point location in $O(\log n)$ time, then algorithm \Alg1 run in $O(\log n + d)$ time, and would be truly optimal. \note{The following go in the point location chapter:} \note{Discuss alternate algorithm where we locate the site closest to $s$, delete all edges that are on the way, and connect it to $s$, then proceed as before. After all, locating the closest site is supposedly what Delaunay diagrams are for!} What about site deletion? Does there exist an $O(n)$ algorithm that, given a Delaunay triangulation $\T\prr$ for $P\prr$ and a site $s$ in $P\prr$, computes a Delaunay triangulation $\T\pr$ for $P\pr=P\prr-\set{s}$? Surprisingly, the question is still open. Removal of a site $s$ and its incident edges from $\T\prr$ leaves a star-shaped ``hole'', and to build $\T\pr$ it is sufficient to triangulate that hole by adding diagonal edges that satisfy the angle property.\figspace3cm \footprob{\Prob{\Sbsec{\Sec3.3}.1}}{Let $\T\star$ be the graph resulting from the removal of an arbitrary site $s$ from a Delaunay triangulation $\T\prr$ of $P\prr$. Prove that there is always a Delaunay triangulation $\T\pr$ of $P\pr=P\prr-\set{s}$ that contains $\T\star$. That is, the only edges we have to delete from $\T\prr$ to build $\T\pr$ are those incident on $s$.} Unfortunately, we know of no algorithm for doing this that takes less than $\Theta(d\log d)$ time, where $d$ is the number of sites on the boundary of the hole (i.e., the degree of $s$). Since $d$ is $\Theta(n)$ in the worst case, this is asymptotically as bad as computing the Delaunay triangulation of $P\pr$ from scratch. In fact, we know of no algorithm faster than $O(n\log n)$ for computing the Delaunay diagram of $n$ sites that are the vertices of a convex polygon. The counterclockwise ordering of the sites around that polygon can be obtained from the Delaunay diagram in linear time; so, if the sites are given in arbitrary order, the $\Theta(n\log n)$ lower bound for sorting $n$ real numbers readily gives a similar bound for this problem. This argument however does not apply if the sites are given in order. In this case there is only an exponential number of triangulations, so straightforward leaf-counting methods in the decision-tree model also fail. We regard this as a major open problem in this area. If we build a Delaunay triangulation of $n$ sites by inserting them one by one with algorithm \Alg1, the total time will be $O(n^2)$, even if we use a naive point location method. Therefore, the preliminarty sorting of the sites by $x$-coordinate does not improve the worst-case time; its main value is the simplification it introduces in step \Step0. \note{??? Conjecture: if points are sorted by $x$ and inserted in middle order (middle point first, then one-quater and three quarters, then one eight, ..., the WORST case is $O(n\log n)$. \note{Probably false.}} \footprob{\Prob{\Sbsec{\Sec3.3}.2}}{(a) Prove or diasprove this conjecture: For any set $P$ of $n$ sites there is an ordering $s1, s2, \ldots, sn$ of them such that if we construct a Delaunay triangulation of $P$ by inserting the sites in that order using algorithm \Alg1, the total time is $O(n\log n)$ time (exclusive of point location). That is, for that particualr ordering the average number of new edges created at each insertion is $O(\log n)$. (b) If the conjecture is true, how fast can we find such an ordering? (c) Let $\bar T(P)$ be the total running time to build a Delaunay triangulation of $P$ by algorithm \Alg1, averaged over all orderings of the sites. What is the maximum of $\bar T(P)$ over all $P$ of size $n$? This might give a useful probabilistic algorithm for Delaunay triangulation.} \subsection{\Ch8.\Sbsec{\Sec3.4}. Expected performance of incremental algorithms} \note{Point location in a triangulation by walking along edges toward query point. $O(n)$ worst case, mumblably $O(n^{1/2})$ average for uniform distribution. Also $O(n^{1/d})$ for $d$ dimensions.} \note{Can we get an $O(n)$ average time if we do point location by bucket sort? Or at least an $O(n\log n)$ average time if we keep an unbalanced kirkpatrician search structure?} \subsection{\Ch8.\Sbsec{\Sec3.5}. Divide-and-conquer algorithms} Another flavor of algorithm for the construction of Delaunay diagrams can be obtained from the standard ``divide and conquer'' paradigm: divide the set $P$ in two halves $X$ and $Y$, costruct Delaunay diagrams $\DX$ and $\DY$ for each half, and ``merge'' them into the complete diagram $\DP$. This may involve deleting some $X$--$X$ or $Y$--$Y$ edges, and will certainly require adding some $X$--$Y$ or {\it cross} edges. However, we do not have to introduce any new $X$--$X$ or $Y$--$Y$ edges: \lemma{\Th{102}}{Let $X$ be a subset of $P$. Then every edge of $\DP$ both of whose endpoints are in $X$ is an edge of $\DX$} \proof{The edge $e$ is a chord of a $P$-free circle, which is also an $X$-free circle.} This result is quite valuable in the design of divide-and-conquer algorithms, as we will see below. If we manage to do the merge step in $O(n)$ time, then the total time for constructing $\DP$ will be $O(n\log n)$. \subsection{\Ch8.\Sbsec{\Sec3.6}. Splitting the sites by a separating line} The algorithm will be greatly simplified if $X$ and $Y$ are {\it linearly separable}, that is, if they lie strictly on opposite sides of some line $l$. We may look at this process (via the quadratic mapping $\lambda$) as a divide-and-conquer method for computing the three-dimensional convex hull: we divide the set of points $\lambda(P)$ in two subsets $\lambda(X)$ and $\lambda(Y)$, compute the hulls of each subset, and merge them into a single hull for $\lambda(P)$. The splitting is done by a plane parallel to the $z$-axis and containing the splitting line $l$. The merge step consists in wrapping a polyhedral ``sleeve'' around the two hulls. \figspace6cm This is quite similar to the classic Preparata-Hong three-dimensional convex hull algorithm [\PH]. It is in fact a special case of the latter, but it is made considerbly simpler because of two reasons. First, we know that {\it a priori} that all points of $\lambda(P)$ lie on the convex hull, so in the merge step we never have to worry about eliminatingvertices and edges of the two partial hulls that fall inside the connecting sleeve. Second, we are interested only on the ``underside'' of $\hull{\lambda(P)}$, and therefore we need to compute only the underside of that sleeve. In turn, this implies that the initial setup and the termination conditions for the merge step will have to be modified. Considering the importance of Delaunay triangulations, we believe it is worth describing these modifications in full. This Delaunay construction algorithm turns out to be the dual of a Voronoi construction algorithm first described by Shamos and Hoey [\ShamHo]. It is also equivalent to (but somewhat simpler than) a Delaunay construction procedure of D. T. Lee [\Lee]. These two algorithms were developed independently from that of Preparata-Hong, and it took quite some time for their equivalence to be recognized. Unlike those two, our algorithm need not compute straight line intersections, unless the coordinates of Voronoi vertices are required output. The only geometric primitives used by our algorithm are the $\incircle$ and $\posor$ predicates. \note{Have we said something about the quad-edge data structure?} In the development of the merge pass, we will be guided by \Th{102}, which says that we only have to add edges between $X$ and $Y$ (and possibly delete old ones). Actually, since we are interested in constructing Delaunay triangulations, rather than Delaunay diagrams, we must extend theorem \Th{102} accordingly: \theorem{\Th{103}}{Let $X, Y$ linearly separable sets of points, and $\TX, \TY$ be Delaunay triangulations of $X$ and $Y$. Then there is a Delaunay triangulation $\TP$ of $P=X\uni Y$ such that \begitems \itemno{1.} Every edge of $\TP$ whose endpoints are both in $X$ (resp. in $Y$) is an edge of $TX$ (resp. $TY$); and \itemno{2.} Every edge of $TX$ (resp. $TY$) that is not in $\TP$ is not a chord of any $P$-free circle. \enditems} \proof{We build a triangulation $\TP$ satisfying this theorem by taking every edge of $\TX$ or $\TY$ whose two endpoints lie on the same face of $\DP$, and then adding enough $X$--$Y$ edges to complete the triangulation. We will show that this excludes only ``definitely bad'' edges of $\TX$ and $\TY$. \pp Let then $F$ be an interior face of $\DP$. Since $F$ is convex, and $X$ and $Y$ are linearly separable, the vertices of $F$ belonging to $X$ must occur in consecutive positions around $F$, and the same must be true of those belonging to $Y$. That is, $F=\seq{x1, x2, \ldots, xm, y1, y2, \ldots, yn}$, with $xi\in X$ and $yi\in Y$ for all $i$ (note that $m$ or $n$ can be zero). \figspace3cm \pp Observe that the circumcircle of $F$ is $P$-free, and therefore also $X$-free; and observe also that the sites of $X$ that lie on the boundary of that circle are precisely $\set{x1, x2, \ldots, xm}$. Therefore, if $m=2$ the edge $x1x2$ is an edge of $\DX$, and therefore is in $\TX$. If $m\geq 3$, the polygon $FX=\seq{x1, x2, \ldots, xm}$ is a face of $\DX$; it follows that the edges of $TX$ whose endpoints are in $\set{x1, x2, \ldots, xm}$ include all sides of $FX$, and also a set of diagonals of $FX$ that decompose it into triangles. \pp In the same way, if $n=2$ the edge $y1y2$ is in $\TY$; if $n\geq 3$, the edges of $\TY$ with endpoints in $\set{y1, y2, \ldots, yn}$ are the sides of $FY=\seq{y1, y2, \ldots, yn}$ plus a set of diagonals that triangulate $FY$. \pp Therefore, if we take all edges of $\TX$ and $\TY$ whose endpoints are both vertices of $F$, we will obtain almost an almost complete triangulation of those vertices. To complete the triangulation $\TP$, we have to add only the edges $xmy1$ and $ynx1$, and possibly one of the diagonals $x1y1$ or $xmyn$, depending on the values of $m$ and $n$: \capfig4cm{$m=0$\hfil $m=1$\hfil $m\geq 2$ and $n\geq2$\hfilneg} The cases $n=0$ and $n=1$ are symmetrical to the above. These extra edges are $X$--$Y$ edges, and, since all faces of $\DP$ are triangulated, the result is a Delaunay triangulation satisfying part 1 of the theorem. \pp An edge $e$ of $\TX$ is not included in $\TP$ only if its endpoints do not belong to the same of $\DP$. By theorem \Th{13B}, that means there is no $P$-free circle passing through that edge. The same applies to $\TY$, and this completes our proof.} \footprob{\Prob{\Sbsec{\Sec3.6}.1}}{Generalize theorem \Th{103}: If $\set{X1, X2, \ldots, Xm}$ is a partition of $P$, and $\Ti$ is a Delaunay triangulation of $Xi$ for each $i$, there is a Delaunay triangulation $\TP$ of $P$ such that every edge of $\TP$ with both endpoints in some $Xi$ is an edge of $\Ti$, and any edges of $\T^i$ not in $\TP$ is not a chord of any $P$-free circle.} \footprob{\Prob{\Sbsec{\Sec3.6}.2}}{Show that theorem \Th{103} remains valid if the sets $X$ and $Y$ are circularly separable. {\it Hint\/}: use the quadratic mapping $\lambda$?} Assume then we have partitioned $P$ into two subsets $L$ and $R$, respectively to the left and to the right of some oriented line $l$. Assume also we have recursively computed Delaunay triangulations $\TL$ and $\TR$ for them. We must now merge the two into a Delaunay triangulation of $P$. According to the above theorem, we have to delete only the ``definitely bad'' $L$--$L$ and $R$--$R$ edges (if any), and add only $L$--$R$ edges. As we will shortly see, it is possible to both delete the bad edges and add the new ones in a single pass, by a very simple and efficient algorithm whose running time is $O(\max\set{\abs{L}, \abs{R}})$. All $L$--$R$ edges we have to add, which we call the {\it cross} edges, intersect the line $l$, and therefore can be lineraly ordered along it. We can talk about successive cross edges, the first cross edge, and so forth. The algorithm we are about to present will add the cross edges one by one in that order. Before adding each new cross edge, the algorithm will delete enough ``definitely bad'' $L$--$L$ and $R$--$R$ edges in such a way that the graph remains at all times a triangular subdivision of the plane. However, the subdivision will have convex perimeter only at the end of the merge pass. What is the structure of the cross edges? Notice that any two consecutive edges crossing $l$ must belong to the same interior triangle of $\TP$. Therefore, those two edges must have a common endpoint on one side of $l$, and they must be connected by an existing $L$--$L$ or $R$--$R$ edge on the other side of $l$. \figspace 4cm These facts have the following important consequence. Suppose we have just determined a cross edge $e$ (which we assume to be oriented from right to left). Let $u\in R$ be its origin and $v\in L$ be its destination. The next cross edge $e\pr$ after $e$ will either go from $u$ to a neighbor $w$ of $v$ in $\TL$, or go from some neighbor $w$ of $u$ in $\TR$ to $v$. In either case, the new vertex $w$ must lie on the right side of the line supporting the edge $e$. \figspace 4cm We can intuitively describe the behavior of the algorithm 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. By induction, the most recent cross edge $e$ is a chord of some $P$-free circle $C$, which established the Delaunayhood of that edge. Consider continuously transforming the circle $C$ into other circles having $e$ as a chord and lying further and further into the left half-plane of $e$. As we know, there is only a single degree of freedom, since the center of the circle is constrained to lie on the bisector of $e$. The transformed circle will remain $P$-free for a while, but unless $e$ is the last cross edge, at some point its circumference will touch a new site $w$. The next cross edge $e\pr$ will then connect either $u$ to $w$ or $w$ to $v$, depending on whether $w$ is in $L$ or in $R$. \figspace 5cm The computation of $w$ proceeds as follows. As we remarked, $w$ lies in the left half-plane of $e$, and is a neighbor of either $u$ or $v$. The algorithm enumerates the neighbors of $u$ in $\TR$, in clockwise order, and selects among them the site $wR$ which will be hit first by the expanding circle. It then enumerates the neighbors of $v$ in $\TL$, in counterclockwise order, and determines which of them ($wL$) will be hit first by the circle. A final test decides which of $wL$ and $wR$ will be encountered first. In principle, we could continue the recursive subdivision of $P$ until we get to sets consisting of a single site. However, the Delaunay diagram of a single site is not a subdivision of the plane as defined in [\GuiSto], and therefore it cannot be represented as a quad-edge data structure. This restriction could be circumvented by appropraiately extending the definitions\foot\ddag{\note{This is not mathematically clean but is quite feasible. We did that in the Cedar implementation.}}. However, if we did so we would have many more special cases to deal with in the merge algorithm. Therefore, we will apply the algorithm only to sets of with two or more sites; in particular, we wil only use the divide-and-conquer method when each of $L$ and $R$ has at least two elements. So, we have to construct the triangulation explicitly when $P$ has only two or three sites. This is a relatively simple task (see problem \Prob{\Sbsec{\Sec3.6}.????}). The first cross edge is essentially one of the outer common tangents to the perimeters of $\TL$ and $\TR$. For clarity, we will describe the computuation of the first cross-edge separately, as algorithm \Alg{7B} below. To help in this computation, the main algorithm always returns a ``root'' edge of the of the constructed triangulation, an edge on the boundary of $\TP$ whose right face is the exterior one of $\TP$. \figspace 4cm In algorithm \Alg7 we use the expression ``$p$ is valid'' to mean that $p$ is a site strictly on the right side of the edge $e$. This predicate is equivalent to $\POS(e\dest, e\org, p)$. \figspace 3cm \algorithm{\Alg7}{Divide-and-conquer Delaunay triangulation with line split.} \algbody{ \step{1}{If $\abs{P} \leq 3$, construct the Delaunay triangulation $\TP$ of $P$ directly and return an arbitrary ``root'' edge of it. Else,} \step{2}{Split $P$ by a line $l$ into two subsets $L$ and $R$ (with at least two vertices each).} \step{3}{Recursively compute the Delaunay triangulations $\TL$ and $\TR$ for $L$ and $R$.} \step{4}{Determine the endpoints of the first cross edge (algorihm \Alg{7B}). Connect those sites by an edge $e$, oriented from $R$ to $L$.} \step{5}{\sna{Merge loop} Repeat the following steps:} \begblock \comm{Enumerate and delete the edges out of $e\dest$, in counterclockwise order, locating the first $L$ site ($a\dest$) that would be hit by the rising bubble:} \step{6}{Set $a\gets e\sym\onext$, and $a\pr\gets a\onext$. If $a\dest$ and $a\pr\dest$ are both valid, then, while $\incircle(e\dest, e\org, a\dest, a\pr\dest)$, do} \begblock \step{7}{Delete the edge $a$. Set $a\gets a\pr$, and $a\pr\gets a\onext$.} \endblock \comm{Now enumerate and delete the edges out of $e\org$, in clockwise order, locating the first $R$ hit ($b\dest$):} \step{8}{Set $b\gets e\oprev$, $b\pr\gets b\oprev$. If $b\dest$ and $b\pr\dest$ are both valid, then, while $\incircle(b\dest, e\dest, e\org, b\pr\dest)$, do} \begblock \step{9}{Delete the edge $b$. Set $b\gets b\pr$, and $b\pr\gets b\oprev$.} \endblock \comm{Check for termination:} \step{10}{If neither $a\dest$ nor $b\dest$ are valid, terminate the algorithm (returning $e$ as the ``root'' edge).} \comm{Now compare $a\dest$ and $b\dest$, and connect the next cross edge to the winner:} \step{11}{If $b\dest$ is not valid, or $a\dest$ is valid and $\incircle(e\dest, e\org, b\dest, a\dest)$, then} \begblock \step{12}{Connect $e\org$ to $a\dest$, and set $e$ to that new edge.} \endblock \unstep{Else} \begblock \step{13}{Connect $b\dest$ to $e\dest$, and set $e$ to that new edge.} \endblock \endblock } When Algorithm \Alg{7B} is called, we know the root edges $rL$ of $\TL$ and $rR$ of $\TR$, returned by the recursive calls of algorithm \Alg7 in step \Step3. The endpoints of the first cross edge are then computed by walking around $\TL$ and $\TR$, in a manner very similar to O'Rourke's algorithm for computing the intersection of convex polygons. The predicate ``$p$ is ahead of $e$'', used in this algorithm, means that the point $p$ lies ahead of the point $e\org$ on the oriented line supporting the directed edge $e$. The predicate ``$p$ is behind $e$ is similarly defined. \figspace 4cm {\bf STOPPED HERE November 1, 1983 2:49 pm} \algorithm{\Alg{7B}}{Common outer tangent of two linearly separable triangulations.} \algbody{ \comm{On input, $\TL$ and $\TR$ are two triangulations, respectively on the left and on the right side of a separating line $l$. We are given an edge $rL$ from the perimeter of $\TL$, with the exterior face at its right, and a similarly placed one $rR$ from $\TR$.} \comm{The output consists of two edges $f$ and $g$ with the same properties, and such that no sites lie on the right side of the line from $f\org$ to $g\org$, or on the segment connecting those two points.} \algfig 2.5cm \comm{\hbox{\titlefont ???}} \step{1}{Set $f\gets rL$, $g\gets rR$.} \step{2}{Repeat the following steps:} \begblock \step{3}{If $g\org$ is to the left of $f$, or ahead of $f$, then set $f\gets f\rprev$; else,} \step{4}{If $f\org$ is to the left of $g$, or behind $g$, then set $g\gets g\lnext$; else} \step{5}{Terminate the loop and go to step 4.} \endblock } We will now elaborate of the above algorithm. Recall first that the number of edges in a triangulation is at most linear on the number of vertices. 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$. We now formally state the lemmas that prove the correctness of the above algorithm. \lemma{9.2}{When the {\lcand} iteration stops, the circumcircle of the triangle defined by {\lcand} and {\basel} is free of other $L$ points. Furthermore all deleted edges are not in the final Delaunay triangulation.} \proof{Let $A$ and $X$ denote respectively the left and right endpoints of {\basel}, and, as in lemma 8.3, let $Ni$, $1ik$, denote the Delaunay neighbors of $A$ in $L$ that are above {\basel}. We prove our lemma inductively on the sequence of cross edges discovered. 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 $ANiN{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. The inclusion of $X$ in the circle $ANmN{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. 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 NkAN{k+1}$ is convex then $X$ is obviously outside the circle $ANkN{k+1}$, so the $\incircle$ test fails. If the angle $\angle NkAN{k+1}$ is concave, then two things happen. First of all the circle $ANkN{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 $ANkN{k+1}$. So in this case the $\incircle$ test fails also. Finally note that the main loop starts with our ascending circle being the halfplane below the lower outer common tangent of $L$ and $R$, and ends with that circle being the halfplane above the upper such tangent.} A symmetric statement holds for the {\rcand} iteration. Furthermore, the final choice between {\lcand} and {\rcand} is again done by the $\incircle$ test, namely {\rcand} will be chosen if \prg{InCircle[lcand\pdest, lcand\porg, rcand\porg, rcand\pdest]} is true. It now follows that: \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.} \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.} 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. \digress{It is instructive to see what those concepts correspond to in the dual domain. Suppose we have computed the Voronoi diagrams for $L$ and $R$. Each diagram covers the whole plane, even though $L$ and $R$ are separated by the line $l$. By the dual of lemma \Th{102}, all edges of $\V P$ on the boundary between two $L$ regions or two $R$ regions are also edges of $\V L$ or $\V R$. Therefore, in the merge step we have only to find the Voronoi edges separating an $L$ region from an $R$ region (and possibly delete some parts of $\V L$ and $\V R$). \dp It turns out that the Voronoi edges to be added constitute a simple path $\pi$ that is monotone in the $y$-direction. In fact, those edges are precisely the duals of the Delaunay cross-edges, and the conclusions we have reached about their structure simply say those cross edges form the dual of a simple path. It can be shown that the Voronoi diagram of $P$ coonsists of those parts of $\V L$ that lie to the left of $\pi$, plus those parts of $\V R$ to the right of $\pi$, plus $\pi$ itself. See [\ShamHo] for more details. \figspace 4cm \dp We would like to remark again that the Voronoi version of the algorithm, which may be easier to understand at an informal level, is noticeably harder to formalize and to prove correct. For example, notice that the dual of a Delaunay {\it triangulation} (which is what the algorithm constructs) may have edges of zero length, and therefore two or more vertices that are topologically distinct and yet have the same position in the plane. \note{Say this earlier in the chapter?} } {\bf STOPPED HERE October 27, 1983 6:40 pm} \note{Separate the common tangent computation?} \note{Decide what root edge(s) is best to return, and fix the tangent computation accordingly. Perhaps return first and last, as in original prog?} \note{Recheck degenerate cases.} \note{Perhaps leave details for program bleow, or delete the latter entirely.} 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{. . .}''. \vfill\eject {\prog \pcomm{create a first cross edge ({\basel}) by computing the lower common tangent of $L$ and $R$} DO WHILE CCW[ldi\porg, ldi\pdest, rdi\porg] DO ldi _ ldi\plnext OD; IF CCW[rdi\pdest, rdi\porg, ldi\porg] THEN rdi _ rdi.Rprev ELSE EXIT; OD; \null basel \pgets Connect[rdi\psym, ldi, left]; \null \pcomm{initialize {\lcand} to the counterclockwise first edge incident to {\basel$\dest$} and above {\basel}} \pcomm{symmetrically for {\rcand}} lcand \pgets basel\psym\ponext; rcand \pgets basel\poprev; \null \pcomm{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 \pcomm{this is the merge iteration} DO \pcomm{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 \pcomm{symmetrically, advance {\rcand} clockwise} . . . \null \pcomm{if both {\lcand} and {\rcand} are invalid then {\basel} is the upper common tangent} IF NOT ValidL[lcand] AND NOT ValidR[rcand] THEN EXIT; \null \pcomm{connect to either {\lcand} or {\rcand}} \pcomm{if both are valid, then choose the appropriate one using the $\incircle$ test} \pcomm{if {\lcand} is the winner then make the next cross edge close the triangle with {\lcand} and {\basel},} \pcomm{else make the next cross edge close the triangle with {\rcand} and {\basel}} \pcomm{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 \pcomm{connect to the right} basel \pgets Connect[rcand, basel\psym, left]; rcand \pgets basel\psym\plnext; END ELSE . . . \pcomm{connect to the left} OD \endprog} \vfill\eject 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$. 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]. \subsection{MORE JUNK:} To simplify the algorithm, we will always split the set $P$ by a line parallel to the $y$-axis, which corresponds to splitting $\lambda(P)$ by a plane parallel to the $y$ and $z$ axes\foot\dag{In practice, it may be better to have a more flexible choice for the splitting line. This may allow us to improve the average-case performance of the algorithm, as we will discuss later on.\note{do it!}}. It is advantageous to start out sorting the points of $P$ according to increasing $x$-coordinate, so that all future splits can be done in constant time. Points with the same $x$-coordinate require special care: such ties must be resolved by sorting in increasing $y$-coordinate. This simplifies considerably the handling of degenerate cases. \footprob{\Prob{\Sbsec{\Sec3.6}.???}}{Let $L$ and $R$ be as defined above. Show that there is a line $l$ separating $L$ from $R$ and {\em containing no point of $P$}.} \footprob{\Prob{\Sbsec{\Sec3.6}.??}}{Modify algorithm \Alg??? so that it constructs the Delaunay diagram of $P$, rather than a Delaunay triangulation.} \footprob{\Prob{\Sbsec{\Sec3.6}.????}}{Write an algorithm that constructs A Delaunay triangulation for a set $P$ with only two or three sites. Be sure to handle the case of three collinear sites. Return the answer as an edge of the triangulation whose right face is the exterior one.} \subsection{\Ch8.\Sbsec{\Sec3.7}. Splitting the sites at random} \subsection{MORE JUNK:} Unfortunately, theorem \Th{103} does not hold in general if $X$ and $Y$ are not linearly separable. 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. \footprob{\Prob{\Sbsec{\Sec3.7}.1}}{Exhibit two Delaunay triangulations $\TX$ and $\TY$ with disjoint (but not linearly separable) vertex sets $X$ and $Y$ such that there is no Delaunay triangulation of $P=X\uni Y$ satisfying condition 1 of theorem \Th{103}.} \section{\Ch8.\Sec4. Lower bounds} % Begin \note{Sorting lower bound.} \note{Do we gain by having points sorted in one or more directions.} \section{\Ch8.\Sec5. Variants of the Vornoi and Delaunay} \subsection{\Ch8.\Sbsec{\Sec5.1}. Line and polygon Voronoi.} \note{Toussaint and others. see also Hill's manuscript.} \subsection{\Ch8.\Sbsec{\Sec5.2}. $k$-point Voronoi.} \note{Definition.} \note{Circumcircle property.} \note{Correspondence to stellations of the Delaunay polytope.} \note{smallest circle enclosing $k$ points} \subsection{\Ch8.\Sbsec{\Sec5.3}. The relative neighborhood graph.} \note{???} \subsection{\Ch8.\Sbsec{\Sec5.4}. The Gabriel graph.} \note{???} \subsection{\Ch8.\Sbsec{\Sec5.5}. Furthest-point Voronoi.} \note{Top of Delaunay polytope.} \note{min spanning circle} \note{$m$ spanning circles} \subsection{\Ch8.\Sbsec{\Sec5.6}. Weighted Voronoi.} \section{\Ch8.\Sec6. Voronoi diagram in other spaces} % Begin The concept of Voronoi diagram is easily extended to other spaces. The definition of Voronoi regions can be applied to any space $\Sscr$: if $K$ is a subset of the sites in $P$, then $\VPK$ is the set of all $x\in \Sscr$ that are equidistant from all sites in $K$, and strictly closer to them than to any site in $P-K$. Lemma \Th1 also remains valid: $$\VPK = \inter{p,q\in K} \E(p,q)\;\inte\;\inter{{\scriptstyle r\in K}\atop{\scriptstyle s\in P-K}} \C(r, s),\eqno(\Eq{201})$$ where $\E(p,q)$ and $\C(p,q)$ are as defined in section \Sbsec{\Sec1.2}. This formula suggests that the topological nature of the Voronoi region $\VPK$ is determined in large part by the nature of the sets $\E(p,q)$ and $\C(p,q)$. The first term of equation (\Eq{201}) is the set $$\E K = \inter{p,q\in K} \E(p,q)$$ of all points of $\Sscr$ that are equidistant from all sites in $K$. Roughly speaking, if $\Sscr$ is a $d$-dimensional space, then $\E (p,q)$ will be a $(d-1)$-dimensional surface, and in general $\E K$ will be a $(d-\leftv K\rightv+1)$-dimensional subset of $\Sscr$. Each $\C(p,q)$ will be a half-space bounded by the surface $\E (p,q)$. Therefore, we can expect $\VPK$ to be a $(d-\leftv K\rightv+1)$-dimensional ``polytope'' whose facets are pieces of bisectors $\E(p,q)$. This analysis is of course a gross simplification, that ignores multitude of degenerate cases possible in spaces of higer dimensions and with non-Euclidean metrics. These cases are discussed in somewhat more detail in the following sections. \footprob{\Prob{\Sbsec{\Sec6.0}.1}}{(a) Prove that in $\Sscr=\reals^k$, $VP(p)$ is simply connected (or empty) under any metric. What about $Vp(K)$ for $\leftv K\rightv > 1$?} \footprob{\Prob{\Sbsec{\Sec6.0}.2}}{(a) Characterize the metrics of $\Sscr=\reals^k$ for which $VP(K)$ is always convex. (b) Same for star-shapedness.} \subsection{\Ch8.\Sbsec{\Sec6.1}. Euclidean spaces of higer dimensions.} In an Euclidean space of dimension $d\geq 2$, the set $\E(p,q)$ is the hyperplane of $\Sscr$ perpendicular to the segment $pq$ and passing through its middle point. The set $C(p,q)$ is the half-space of $\Sscr$ delimited by $E(p,q)$ and contining $p$. The proof of these assertions is a trivial generalization of \Prob{\Sbsec{\Sec1.2}.2}. \theorem{\Th34}{Let $P$ is a finite set of points of an Euclidean space $\Sscr$. For any nonempty subset $K$ of $P$, $VP(K)$ is either empty or an open convex polytope in the subspace $E(K)$.} \proof{We can rewrite equation (\Eq2) as $$\eqalign{VP(K)=E(K)\;\inte\;\inter{{p\in K}\atop{q\in P-K}} C(p,q)\cr\vsk4 \null=\inter{{p\in K}\atop{q\in P-K}}(E(K)\inte C(p,q)).\cr}$$ Each term $E(K)\inte C(p,q)$ is either empty or is an open half-space of $E(K)$; therefore, the intersection of all those terms $V(K)$ is either empty or an open convex polytope of $E(K)$.} From what we have seen so far, we can conclude that a Voronoi diagram is a partition of $\Sscr$ into open convex polytopes of varying dimensions. In order to show that $VP$ is a convex subdivision of $\Sscr$, as \note{should have been} defined in chapter 5, we have only to show that the proper incidence relationships hold for the elements of $VP$; for example, that the endpoints of an edge (if any) are vertices, that the boundary of a face consists of vertices and edges only, and so forth. These properties follow readily from the next lemma: \lemma{\Th{35}}{The closure of a voronoi region $\V P K$ is $$\union{P\subset X\subset K} \VPX.$$} \theorem{\Th{36}}{The Voronoi diagram of a finite set $P$ in Euclidean space $\reals^k$ is a simple subdivision of $\reals^k$.} Further insight into the meaning of theorem \Th{??} can be gained by considering the three-dimensional case. The region $V(p)$ corresponding to a single site $p$ is an open (an possibly unbounded) convex polyhedron containing $p$, called a {\it cell} of the diagram. If $K$ consists of two distinct sites $p$ and $q$, then $E(K)$ is the bissecting plane of those two points, and the Voronoi region $V(K)$ is either empty or an open convex polygon (a {\it face}) contained in that plane. If $K$ consists of three non-collinear sites $p,q$ and $r$, the bisector $E(K)$ is the line perpendicular to the plane $pqr$ and passing through their circumcenter. The region $V(K)$ will be a convex open subset of that line, i.e. an {\it edge} of the diagram. If $K$ consists of four non-coplanar sites, then $E(K)$ is the center of the circumscribed sphere, and $V(K)$ is either empty or that single point (a {\it vertex} of $VP$). In general, if $K$ is any subset of $P$ containing four non-coplanar sites, $V(K)$ is either empty or a vertex. If all sites in $K$ are coplanar, but at least three of them are not collinear, then $V(K)$ is either empty or an edge. If all sites in $K$ are collinear, but $K$ contains at least two distinct points, then $V(K)$ is either empty or a face of the diagram. Finally, if $K$ contains a single point, then $V(K)$ is a cell. These results can be expressed in a more general form: \theorem{\Th{37}}{If $K\subset P\subset \reals^k$ is nonempty and has affine rank $r$, then the Voronoi region $VP(K)$ is either empty or has dimension $(k-r)$.} \footprob{\Prob{\Sbsec{\Sec6.1}.7}}{Prove theorem \Th{37}.\lp {\it Ans.}: Recall that the affine rank of a set $K$ is the dimension of the smallest affine space containing $K$. This is equivalent to saying that $K$ is contained in some $r$-dimensional affine subspace $X$ of $\Sscr$, and there are $r+1$ sites $p0, p1, \ldots, pr$ in $K$ such that the vectors $p1-p0, p2-p0, \ldots, pr-p0$ are linearly independent.\lp Each bisector $E(p0, pi)$ is an hyperplane perpendicular to the vector $pi-p0$; therefore, the intersection $E\pr=\interi E(p0, pi)$ of those $r$ bisectors is an affine subspace perpendicular to all those vectors, and therefore to the linear subspace generated by them.\lp \note{must prove that dim $E\pr=k-r$ exactly.} The latter has dimension $r$, the dimension of $E\pr$ is at most $k-r$. Since $E(K)\subset E\pr$, the same is true of $E(K)$.\lp \note{must show that every other bisector is either parallel to $E\pr$ or contains it, so the intersection is either empty or $E\pr$.}} \footprob{\Prob{\Sbsec{\Sec6.1}.8}}{Prove that the Delaunay diagram (where $D(K)=\hull(K)$ if $V(K)\neq\empty$) in a Cartesian space is always the dual of the Voronoi diagram, even if the metric is not Euclidean, or give a counterexample.} \footprob{\Prob{\Sbsec{\Sec6.1}.9}}{Suppose the sites (but not necessarily the query points) lie on some $m$-dimensional affine subspace of $\Sscr=\reals^k$. Show that the post-office problem in $\Sscr$ with Euclidean metric can be reduced to the post office problem in $\reals^m$, with only $O(nk)$ extra preprocessing time and $O(k)$ extra query time.} \footprob{\Prob{\Sbsec{\Sec6.1}.10}}{(a) Suppose it is known the sites are the vertices of a convex $k$-dimensional polytope. Can this knowledge be exploited to give faster post-office algorithms (both with and without preprocessing) than the general case, in the Euclidean metric? (b) What about $l1$ metric? (c) Can we do better if the sites are known to lie on the surface of a $k$-dimensional sphere?} \footprob{\Prob{\Sbsec{\Sec6.1}.11}}{Show that in Eucludean spaces of arbitrary dimension a region $V(K)$ is unbounded if and only if it is nonempty and $\hull(K)$ is contained in the boundary of $\hull(P)$. In particular, show that if $P$ is in general position (any subset of $k+1$ sites has affine rank $k$) then a region $V(K)$ is unbounded if and only if $\hull(K)$ is a facet of $\hull(P)$.} To complete the characterization of Voronoi regions, we will give a condition for $V(K)$ and $D(K)$ to be nonempty. Let's define a {\it $P$-free ball} as a $k$-dimensional ball of $\Sscr$ that contains no sites of $P$ in its interior. Then the following holds: \theorem{\Th{37}}{A point $x$ belongs to $VP(K)$ if and only if there is an open $P$-free ball $S$ centered at $x$ such that the its boundary contains $K$ but no other sites of $P$.} \proof{Immediate from the definitions} So, $V(K)$ is nonempty if and only if all sites in $K$, and only those, lie on the boundary of some $P$-free ball. \definition{\Def2}{Let $P$ be a finite set of points (sites) in $\Sscr$, and $K$ be any subset of $P$. If the Voronoi region $VP(K)$ is nonempty, the corresponding {\it Delaunay region} is the open convex hull of $K$; otherwise $DP(K)$ is $\empty$. The collection $\rset{D{K}\relv K\subset P}$ of all Delaunay regions is called the {\it Delaunay diagram} of $P$.} The open convex hull of a set $K$ is just the convex hull of $K$ minus its boundary in the affine space spanned by $K$; that is, the open hull of a singleton set $\set{p}$ is $\set{p}$ itself, the open hull of two points $\set{p,q}$ is the open line segment $pq$, and so forth. The definition above provides a natural correspondence betwen Voronoi and Delaunay regions, which, as we will see, establishes the duality between $VP$ and $DP$. In fact, this definition (and the duality it defines) remains valid in any Euclidean space of arbitrary dimension. We can already see (by theorem \Th?) that a non-empty Delaunay region $D(K)$ has dimension $r$ if and only if the corresponding Voronoi region $V(K)$ has dimension $k-r$. \note{Remember to prove that $DP$ is a simple subdivision of $\Sscr$, dual of $VP$. Is it so for any metric?} By theorem 8.?, $\hull(K)$ is a Delaunay region if and only if there is a $P$-free $k$-dimensional ball passing through the sites of $K$ and through no other sites. In the particular case of a planar Delaunay, we have \note{mention that in a ``general'' set of sites every $k+1$ points determine a single sphere, and $k+2$ or more points never lie on the same sphere.} \theorem{8.?}{The Voronoi diagram $VP$ is a simple subdivision of $\Sscr$.} \proof{Since the regions $VP(K)$ are pairwise disjoint, we can conclude from the previous lemma that the closure $\bar{V(K)}$ of a region $V(K)$ is the union of Voronoi regions. Since the regions are convex, by theorem 5.? {\bf that doesn't exist yet} we conclude that $VP$ is a simple subdivision of $\Sscr$.} \note{By a {\bi subdivision} of a topological space $X$ we mean a partition of $X$ into $k+1$ collections $\Fscr0, \Fscr1, \ldots, \Fscrk$ of subsets, with the properties that (1) each element of $\Fscri$ is homeomorphic to the $i$-dimensional ball $\spherei$, and (2) the closure of each element of $\Fscri$ is the union of elements of the collections $\Fscrj$ with $jJ˜âJ˜WJ˜ÚJ˜ÄJ˜²J˜®J˜ˆJ˜´J˜uJ˜HJ˜´J˜=J˜ÆJ˜íJ˜SJ˜ðJ˜J˜xJ˜XJ˜¾J˜ºJ˜mJ˜†—˜EJ˜¬J˜óJ˜˜J˜¹J˜ÜJ˜ƒJ˜¸J˜NJ˜¯——˜.˜J˜ßJ˜âJ˜¯J˜„J˜ÚJ˜þJ˜ªJ˜¨J˜®J˜ØJ˜ìJ˜”—˜4J˜J˜ J˜‡J˜üJ˜ÆJ˜k—˜CJ˜ J˜¿J˜¬J˜½J˜³J˜©J˜J˜"—˜6J˜ J˜“—˜:J˜áJ˜Ú——˜D˜J˜¶J˜ÕJ˜ÔJ˜yJ˜¯J˜œJ˜ J˜ÆJ˜“—˜?J˜¹J˜ZJ˜òJ˜…J˜ÕJ˜öJ˜ÂJ˜•—˜@J˜7J˜ÍJ˜•J˜ñJ˜¢J˜•J˜îJ˜ÆJ˜†J˜ÀJ˜J˜ÉJ˜£J˜»J˜‹J˜‘J˜µJ˜áJ˜”J˜@J˜IJ˜¿J˜…J˜ŸJ˜WJ˜ÞJ˜›—˜@J˜ J˜•J˜«J˜‰J˜ÿJ˜ðJ˜ñJ˜•J˜ÔJ˜ÉJ˜7J˜ôJ˜áJ˜ùJ˜ÅJ˜ÁJ˜äJ˜ÕJ˜¯—˜RJ˜ÆJ˜³—˜AJ˜òJ˜€J˜XJ˜Ø—˜LJ˜˜J˜µJ˜ÁJ˜ºJ˜ÄJ˜÷J˜µJ˜”J˜´J˜óJ˜J˜J˜ŒJ˜³J˜€J˜ÜJ˜ÍJ˜áJ˜„J˜J˜³J˜µJ˜ÉJ˜NJ˜¡J˜§J˜ÝJ˜,J˜UJ˜µJ˜³J˜TJ˜ÛJ˜ŽJ˜ÝJ˜§J˜ÝJ˜×J˜žJ˜ØJ˜J˜¥J˜¬J˜ J˜ãJ˜,J˜0J˜”J˜!J˜OJ˜­J˜ J˜åJ˜ J˜¸ J˜ÎJ˜J˜åJ˜¨J˜˜J˜—˜AJ˜J˜dJ˜•J˜†——˜#˜J˜J˜EJ˜——˜:˜=J˜9—˜6J˜J˜J˜?J˜,—˜DJ˜ —˜6J˜ —˜;J˜!J˜J˜—J˜5—˜6˜J˜¬J˜ÒJ˜òJ˜²J˜™—˜IJ˜ÓJ˜ÂJ˜øJ˜¦J˜fJ˜J˜çJ˜³J˜°J˜7J˜£J˜üJ˜ðJ˜ãJ˜—J˜ŽJ˜…J˜¸J˜'J˜sJ˜îJ˜ØJ˜pJ˜ÚJ˜—J˜MJ˜¹J˜ê—˜IJ˜/J˜—˜>J˜GJ˜·—˜9J˜¿J˜ÿJ˜¬J˜gJ˜VJ˜–——˜6˜J˜J˜——˜ ˜J˜¯J˜‚J˜QJ˜–J˜KJ˜‡J˜«J˜´J˜…J˜œJ˜²J˜¤J˜˜J˜gJ˜ž——J˜J˜—…—¢ u