\lecturetitle{12}{Triangulation 1}{8 October 1997}{Aaron Brown}

% INCLUDE AN OUTLINE OF THE CLASS TOPICS
\section*{Outline}
\begin{itemize} 
  \item Introduction to Triangulation

  \item Basic Geometric Operations
    \begin{enumerate}
    \item Degeneracies and Precision Issues
    \item Line Side Test
    \item In-Circle Test
    \end{enumerate}

  \item Convex Hulls
    \begin{enumerate}
    \item Dual Problem: Convex Polytopes
    \item Algorithms
	\begin{itemize}
	\item Giftwrapping
	\item Graham Scan
	\item Mergehull
	\item Quickhull
	\end{itemize}
    \end{enumerate}

  \item Voronoi Diagrams

  % ugh. not \latex2e, so use ancient font-changing commands
  \item Delaunay Triangulations ({\it dual of Voronoi diagrams})
    \begin{enumerate}
    \item Properties
    \item Algorithms
        \begin{itemize}
	\item Naive Algorithm
	\item Solution via Reduction to 3D Convex Hull
	\item Direct Algorithms
        \end{itemize}
    \end{enumerate}
\end{itemize}

The following three lectures cover triangulation, which is one of the
many problems in the larger field of ``computational geometry''.
Problems in computational geometry have found many applications in
industry.  We will also cover convex-hulls and Voronoi-diagrams
because they are closely related to triangulation.

\section{Introduction to Triangulation}

Triangulation is a process that takes a region of space and divides it
into subregions. The space may be of any dimension (although we will
focus primarily on the two-dimensional problem here), and the
subregions typically are shaped as triangles (2D), tetrahedrons (3D),
etc. The goal of triangulation is to produce a subdivision that obeys
certain ``nice'' properties, for example one in which the output
triangles are as close to equilateral as possible.

Triangulation is used in many real-world applications. One such
application is {\it finite element simulation}, in which a partial
differential equation (PDE) is simulated over a continuous surface or
space by first using a triangulation to break the space into small,
discrete elements, then simulating the PDE discretely over the
elements.  Finite element simulations are common in fluid flow
applications and similar problems; an example that illustrates the
role of the triangulation can be seen at:

\medskip
\indent \indent {\tt http://www.cs.cmu.edu/\~\/scandal/applets/airflow.html}
\medskip

Other important applications of triangulation include:
\begin{itemize}
\item {\it Surface approximation}: a continuous surface is
approximated by a mesh of triangles. This is used in graphics
applications (e.g., face representation), compression of surface
representations, interpolation of surfaces, and construction
applications ({\em e.g.}, optimizing dirt movement to contour a region
for new construction).

\item {\it Nearest neighbor structure identification}: the
triangulation includes data identifying the nearest neighbors of each
point in the triangulated data set. This is potentially useful for
clustering applications, etc.
\end{itemize}

In general, triangulation operates on several different types of
input. The simplest case (and the one that we focus on in these notes)
is when the input is just a set of points which define the vertices of
the output triangulation. Another possibility is that the input data
is a boundary defining a region of space. In this case, the
triangulation must generate the interior points needed to complete the
triangulation. Finally, the input can be specified as a density
function $\varphi(x,y,z)$, which specifies the density of triangles
across the region to be triangulated. This last format is especially
useful in simulation, as it allows the person developing the model to
focus triangle detail (and thus simulation accuracy) on critical
portions of the model.

\section{Basic Geometrical Operations}
\subsection{Degeneracies and Numerical Precision}
Before getting into the triangulation algorithms themselves, we first
need to consider two important issues that arise when these algorithms
are implemented in the real world: handling of Degeneracies and
numerical precision.

Degeneracies occur (in two dimensions) when three or more points are
co-linear, or when four or more points lie on the same circle.  
Special
care is often needed in triangulation algorithms to make them robust with
respect to degeneracies.

Another related problem that crops up in the real world is the issue
of numerical precision.  Since (in general) computers operate with
only a fixed amount of floating-point precision errors are introduced
by most arithmetic operations.  These errors can lead to inconsistent
answers, especially when points are almost degenerate.  For example if
three points $a$, $b$ and $c$ are almost on a line the answer to
whether $a$ is above $b-c$ could differ depending on how you calculate
it.
There are three potential solutions to the numerical
precision problem (and the related degeneracy problem):
\begin{enumerate}
\item Assume that the problem doesn't exist. Many existing
implementations do this, but it is a bad idea.
\item Design the algorithm to be robust with respect to degeneracies
and to lack of precision (e.g., by using exact or adaptive
arithmetic). This is the best solution.
\item Randomly perturb the input data to avoid degeneracies. This may
work in most cases, but the possibility of degeneracies or precision
problems still remains.
\end{enumerate}

Luckily, most triangulation algorithms only use two basic operations
which are sensitive to degeneracies and numerical precision: the {\it
line side test} and the {\it in-circle test}. These operations
encapsulate all of the needed coordinate comparisons, and thus they
are the only routines that need be implemented in infinite-precision
or adaptive-precision arithmetic. The rest of the triangulation
algorithms are primarily combinatorics and graph manipulations.

\subsection{The Line-side Test}
\label{lineside}
Given a point $C$ and a line $AB$, the {\it line-side test} answers
the question: ``is $C$ above or below $AB$?''

\begin{figure}[hbtp]
\begin{center}
\setlength{\unitlength}{.5in}
\begin{picture}(4,2.5)
\put(1,0){\circle*{.125}}
\put(0.9,0){\makebox(0,0)[r]{\it A}}
\put(3.07,2.07){\circle*{.125}}
\put(3.15,2.15){\makebox(0,0)[bl]{\it B}}
\put(1,0){\vector(1,1){2}}
\put(3,.5){\circle*{.125}}
\put(3.1,0.5){\makebox(0,0)[l]{\it C}}
\end{picture}
\end{center}
\end{figure}

The line-side test is easily implemented via a matrix determinant
operation; if the determinant operation is built with
adaptive-precision arithmetic, the line-side test will be stable with
respect to numerical precision. Let the line be defined by points
$A=(x_a,y_a)$ and $B=(x_b,y_b)$, and let the test point be
$C=(x_c,y_c)$. Then let:
\[
D = \left| 
	\begin{array}{ccc}
	x_a & y_a & 1 \\
	x_b & y_b & 1 \\
	x_c & y_c & 1 \\
	\end{array}
    \right|
\]

\noindent Then the point $C$ is ``above'' the line $A-B$ iff $D >
0$.  Notice as well that the area of triangle $ABC$ is given by $D/2$.

\begin{figure}[htbp]
\begin{center}
\setlength{\unitlength}{.5in}
\begin{picture}(3,1.5)
\put(1,0){\circle*{.125}}
\put(3.1,0){\circle*{.125}}
\put(2,1){\circle*{.125}}
\put(1,0){\vector(1,0){2}}
\put(0.9,0){\makebox(0,0)[r]{(0,0)}}
\put(3.2,0){\makebox(0,0)[l]{(2,0)}}
\put(2,1.2){\makebox(0,0)[b]{(1,1)}}
\end{picture}
\end{center}
\caption{Example of applying the line-side test.}
\label{fig2}
\end{figure}

As an example, consider the line $AB$ and the point $C$ with
$A=(0,0)$, $B=(2,0)$, and $C=(1,1)$, as illustrated in Figure~\ref{fig2}.
In this case,
\[
D = \left|
	\begin{array}{ccc}
	0 & 0 & 1 \\
	2 & 0 & 1 \\
	1 & 1 & 1 \\
    	\end{array}
    \right|
  = 2 > 0.
\]

\noindent Thus, $C$ is above $AB$, and the triangle area is
$\frac{2}{2} = 1$.

\subsection{The In-circle Test}
Given three points $A$, $B$, $C$, and $D$, the {\it in-circle test} answers
the question: ``is $D$ in the circle defined by $A$, $B$, and $C$?''

\begin{figure}[hbtp]
\begin{center}
\setlength{\unitlength}{.25in}
\begin{picture}(4,2.5)
\put(2,1){\circle{2}}
\put(2.5,1.2){\circle*{.2}}
\put(2.3,1.0){\makebox(0,0)[ur]{\it D}}
\put(2,2){\circle*{.2}}
\put(2,2.2){\makebox(0,0)[b]{\it B}}
\put(1.2,1.6){\circle*{.2}}
\put(1.0,1.6){\makebox(0,0)[r]{\it A}}
\put(2.867,.5){\circle*{.2}}
\put(3.067,.5){\makebox(0,0)[l]{\it C}}
\end{picture}
\end{center}
\end{figure}

As in the line-side test, the in-circle test can be implemented with a
matrix determinant. Again, if the determinant operation uses
adaptive-precision arithmetic, the in-circle test should be
stable. Notice, however, that $A$, $B$, and $C$ only define a unique
circle if they are not colinear. 

As before, let $A=(x_a, y_a)$, $B=(x_b, y_b)$, $C=(x_c,y_c)$, and
$D=(x_d,y_d)$. Then the particular matrix determinant used for the
in-circle test is:
\[
D = \left| 
	\begin{array}{cccc}
	x_a & y_a & x_a^2 + y_a^2 & 1 \\
	x_b & y_b & x_b^2 + y_b^2 & 1 \\
	x_c & y_c & x_c^2 + y_c^2 & 1 \\
	x_d & y_d & x_d^2 + y_d^2 & 1 \\
	\end{array}
    \right|
\]

\noindent Then the point $D$ is in the circle defined by $ABC$ iff $D >
0$. 

Notice that this algorithm generalizes to higher dimensions (as does
the line-side test). In this case, in $d$ dimensions, $d+1$ points
determine the $d$-dimensional sphere, and $D$ becomes a $(d+2) \times
(d+2)$ determinant.

\section{Convex Hulls}
Given a set of points $S$, the {\it convex hull problem} is to find
the smallest convex region that contains all points in $S$. The boundary
of the region is called the {\it convex hull} of $S$, and is usually
specified by a set of points and their neighbor relation. 

\begin{figure}[htbp]
\begin{center}
\setlength{\unitlength}{0.00041300in}
\begin{picture}(4939,5007)(2639,-6362)
\put(4966,-1816){\circle*{128}}
\put(2701,-4201){\circle*{128}}
\put(3316,-5701){\circle*{128}}
\put(4816,-6301){\circle*{128}}
\put(7066,-5911){\circle*{128}}
\put(7516,-4126){\circle*{128}}
\put(6751,-2251){\circle*{128}}
\put(3841,-4411){\circle*{128}}
\put(6247,-4804){\circle*{128}}
\put(3766,-2401){\circle*{128}}
\put(4300,-3700){\circle*{128}}
\put(3751,-2386){\line( 2, 1){1200}}
\put(4951,-1786){\line( 4,-1){1800}}
\put(6751,-2236){\line( 2,-5){750}}
\put(7501,-4111){\line(-1,-4){450}}
\put(7051,-5911){\line(-6,-1){2250}}
\put(4801,-6286){\line(-5, 2){1500}}
\put(3301,-5686){\line(-2, 5){600}}
\put(2701,-4186){\line( 3, 5){1072.059}}
\put(5221,-3241){\circle*{128}}
\put(4840,-6025){\circle*{128}}
\put(7081,-3961){\circle*{128}}
\put(6404,-2587){\circle*{128}}
\put(4450,-2530){\circle*{128}}
\put(5146,-4975){\circle*{128}}
\put(3505,-5440){\circle*{128}}
\put(3331,-3385){\circle*{128}}
\put(8566,-5146){\vector(-1, 1){925}}
\put(8600,-5146){\makebox(0,0)[l]{Extreme Point}}
\put(8116,-2656){\vector(-2,-1){960}}
\put(8180,-2656){\makebox(0,0)[l]{Convex Hull}}
\put(6112,-1198){\vector(-2,-1){1066}}
\put(6180,-1198){\makebox(0,0)[l]{Extreme Point}}
\end{picture}
\end{center}
\caption{Example of a 2D Convex Hull}
\label{fig4}
\end{figure}

Convex hulls are well-defined in any dimension; a 2-dimensional
example is given in Figure~\ref{fig4}. The 2-D convex hull can be
thought of the shape that results when a rubber band is stretched to
surround a set of nails, one at each point in $S$; similarly, the 3-D
convex hull is similar to the shape that would be formed if a balloon
were stretched to tightly surround the points $S$ in 3-space.

\subsection{Dual: Convex Polytopes}
\begin{figure}[thbp]
\setlength{\unitlength}{0.00041300in}%
%
\begingroup\makeatletter\ifx\SetFigFont\undefined
% extract first six characters in \fmtname
\def\x#1#2#3#4#5#6#7\relax{\def\x{#1#2#3#4#5#6}}%
\expandafter\x\fmtname xxxxxx\relax \def\y{splain}%
\ifx\x\y   % LaTeX or SliTeX?
\gdef\SetFigFont#1#2#3{%
  \ifnum #1<17\tiny\else \ifnum #1<20\small\else
  \ifnum #1<24\normalsize\else \ifnum #1<29\large\else
  \ifnum #1<34\Large\else \ifnum #1<41\LARGE\else
     \huge\fi\fi\fi\fi\fi\fi
  \csname #3\endcsname}%
\else
\gdef\SetFigFont#1#2#3{\begingroup
  \count@#1\relax \ifnum 25<\count@\count@25\fi
  \def\x{\endgroup\@setsize\SetFigFont{#2pt}}%
  \expandafter\x
    \csname \romannumeral\the\count@ pt\expandafter\endcsname
    \csname @\romannumeral\the\count@ pt\endcsname
  \csname #3\endcsname}%
\fi
\fi\endgroup
\begin{center}
\begin{picture}(9174,4983)(2684,-6190)
\put(10319,-6079){(b)}
\put(2926,-1666){\line( 1, 0){4755}}
\put(4951,-5506){\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(4966,-5506){\line( 5, 4){3315}}
\put(4531,-5071){\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(3886,-1658){\vector( 0,-1){173}}
\put(6773,-1673){\vector( 0,-1){173}}
\multiput(3571,-2603)(6.00000,6.00000){21}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(3691,-2483){\vector( 1, 1){0}}
\multiput(4385,-3567)(6.00000,6.00000){21}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(4505,-3447){\vector( 1, 1){0}}
\multiput(5296,-4689)(6.00000,6.00000){21}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(5416,-4569){\vector( 1, 1){0}}
\multiput(4734,-2492)(6.00000,-6.00000){26}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(4884,-2642){\vector( 1,-1){0}}
\multiput(3706,-3287)(6.00000,-6.00000){26}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(3856,-3437){\vector( 1,-1){0}}
\multiput(5802,-1708)(6.00000,-6.00000){26}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(5952,-1858){\vector( 1,-1){0}}
\multiput(5941,-2125)(-5.91667,-5.91667){25}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(5799,-2267){\vector(-1,-1){0}}
\multiput(6953,-3010)(-5.91667,-5.91667){25}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(6811,-3152){\vector(-1,-1){0}}
\multiput(5036,-1394)(-5.91667,-5.91667){25}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(4894,-1536){\vector(-1,-1){0}}
\put(6339,-3790){\vector( 0, 1){180}}
\put(4021,-3797){\vector( 0, 1){180}}
\put(5956,-3174){\vector(-1, 1){157}}
\put(5190,-4134){\vector(-1, 1){157}}
\put(4662,-4783){\vector(-1, 1){157}}
\put(6671,-2278){\vector(-1, 1){157}}
\multiput(6376,-4367)(-6.04762,6.04762){22}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(6249,-4240){\vector(-1, 1){0}}
\put(5506,-6109){(a)}
\multiput(5243,-5274)(-6.04762,6.04762){22}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(5116,-5147){\vector(-1, 1){0}}
\put(9379,-3604){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}b}}}
\multiput(7973,-3092)(-6.04762,6.04762){22}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(7846,-2965){\vector(-1, 1){0}}
\put(5026,-3783){\vector( 0, 1){180}}
\put(9337,-3058){\line( 5,-6){595}}
\put(10993,-1844){\line( 6,-5){840}}
\put(10862,-3764){\line( 4, 5){984}}
\put(10989,-1840){\line(-4,-3){1652}}
\put(3151,-2116){\line( 5,-6){2765}}
\put(4906,-1276){\line( 6,-5){3168}}
\put(3466,-3796){\line( 1, 0){4875}}
\put(4471,-5026){\line( 4, 5){2472}}
\put(6346,-1291){\line(-4,-3){2924}}
\put(9937,-3784){\line( 1, 0){930}}
\put(2684,-1654){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}a}}}
\put(2929,-2164){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}b}}}
\put(4624,-1279){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}c}}}
\put(3124,-3499){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}d}}}
\put(3244,-3874){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}e}}}
\put(4219,-5104){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}f}}}
\put(4699,-5584){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}g}}}
\put(10064,-2229){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}d}}}
\put(11414,-2044){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}c}}}
\put(11499,-3319){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}f}}}
\put(10319,-4094){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}e}}}
\end{picture}
\end{center}
\caption{(a) Example 2D Convex Polytope Problem (b) Solution}
\label{fig5}
\end{figure}

The convex hull problem has a dual formulation in the problem of {\it
convex polytopes}. The convex polytope problem is to find the region
in space that is the intersection of a set of halfspaces (we assume
all the halfspaces include the origin).  This problem
arises in certain linear programming algorithms. Figure~\ref{fig5}a
illustrates an example of a 2-D convex polytope problem. The solution
to a convex polytope problem in general is a list of the
``non-redundant'' half-spaces (e.g., those lines, planes, or
hyperplanes that define the convex region). In the example of
Figure~\ref{fig5}a, the solution is $(b,d,c,f,e)$, illustrated in
Figure~\ref{fig5}b.

The convex hull problem and the convex polytope problems are
duals, as a problem and solution to one can be mapped into a
problem and solution to the other. Thus an algorithm for finding
convex hulls can be used to find convex polytopes as well.

The transformation that takes an instance of the convex polytope
problem to a convex hull problem is the following, in 2D.
Let the $i$th half-space be defined by the line $\ell_i : a_i x + b_i y
= 1$.  In higher dimensions, the equation is of the form:
\[
\ell_i : \sum_{j=1}^d a_{i_j} x_j  = 1.
\]

Then the dual problem is to find the convex hull of the points $p_i$,
where (in 2D):
\[
P_i = \left(\frac{a_i}{a_i^2 + b_i^2}, \frac{b_i}{a_i^2+b_i^2}\right).
\]
\noindent In higher dimensions,
\[
P_i = \left(\frac{a_{i_1}}{a_{i_1}^2 + a_{i_2}^2 + \cdots + a_{i_d}^2},
	\frac{a_{i_2}}{a_{i_1}^2 + a_{i_2}^2 + \cdots + a_{i_d}^2},
	\ldots,
	\frac{a_{i_d}}{a_{i_1}^2 + a_{i_2}^2 + \cdots + a_{i_d}^2},
        \right).
\]

\noindent See Figure~\ref{fig7} for a graphical illustration.

\begin{figure}[hbtp]
\setlength{\unitlength}{0.00043300in}
\begin{center}
\begin{picture}(7920,3420)(5834,-6301)
\put(13289,-3739){\circle*{128}}
\put(12074,-5314){\circle*{128}}
\put(13649,-5314){\circle*{128}}
\put(6089,-4669){\line( 5, 3){2565}}
\put(7514,-2899){\line( 4,-3){2684}}
\put(6254,-3919){\line( 1,-1){2100}}
\put(7229,-5944){\line( 2, 1){3196}}
\put(12076,-3736){\line( 0,-1){1575}}
\put(12076,-5311){\line( 1, 0){1575}}
\put(13651,-5311){\line(-1, 5){315.654}}
\put(13304,-3739){\line(-1, 0){1228}}
\put(12073,-3741){\circle*{128}}
\put(11064,-4639){\makebox(0,0)[lb]{$\Longrightarrow$}}
\put(11834,-3459){\makebox(0,0)[lb]{$p_4$}}
\put(11909,-5584){\makebox(0,0)[lb]{$p_3$}}
\put(13754,-5584){\makebox(0,0)[lb]{$p_1$}}
\put(13544,-3754){\makebox(0,0)[lb]{$p_2$}}
\put(6889,-6094){\makebox(0,0)[lb]{$\ell_1$}}
\put(8409,-6274){\makebox(0,0)[lb]{$\ell_3$}}
\put(7174,-2989){\makebox(0,0)[lb]{$\ell_2$}}
\put(5734,-4774){\makebox(0,0)[lb]{$\ell_4$}}
\end{picture}
\end{center}
\caption{Mapping Convex Polytopes to Convex Hulls}
\label{fig7}
\end{figure}

\subsection{Algorithms for Finding Convex Hulls}
The following is a summary of the most popular algorithms for finding
convex hulls. Notice that the algorithms are classified into
categories that are similar to those used to classify sorting
algorithms. Those marked with a dagger ($\dagger$) are considered further
below. In the time bounds, $n$ is the number of points in the input
set $S$, and $h$ is the number of points in the hull $H$.
The algorithms and bounds we list here are for 2-dimensions although
many algorithms generalize to higher dimensions.

\begin{itemize}
\item Basic algorithms
    \begin{itemize}
    \item Giftwrapping$^\dagger$: $O(nh)$
    \end{itemize}
\item Incremental algorithms
    \begin{itemize}
    \item Graham Scan$^\dagger$: $T($sort$) + O(n)$
    \item Unordered Insertion: $T($sort$) + O(n)$
    \end{itemize}
\item Divide and Conquer algorithms
    \begin{itemize}
    \item ``Merge-based'' D\&C algorithms
	\begin{itemize}
	\item Mergehull$^\dagger$: $T($sort$) + O(n)$
	\end{itemize}
    \item ``Divide-based'' D\&C algorithms
	\begin{itemize}
	\item Quickhull$^\dagger$: $O(n)$ best, $O(n^2)$ worst, $O(n \log n)$ average
	\item Kirkpatrick-Seidel: $O(n \log h)$
	\end{itemize}
    \item Random Sampling: $O(n \log n)$
    \end{itemize}
\end{itemize}

\noindent The theoretical lower-bound on finding a convex hull is $O(n
\log h)$.

\subsubsection{Giftwrapping}
\index{Convex Hull: Giftwrapping}

The idea behind the giftwrapping algorithm (in 2-d) is to move from point to
point along the convex hull. To start, an extreme point is chosen as
the first point on the hull (this can be done easily by choosing the
point in $S$ with largest $x$-value).

Then, on each step, the point in $S$ with the smallest angle $\alpha$
is chosen to be the next point in the hull, where $\alpha$ is defined
as the angle between a ray originating at the current point and
continuing along the ray that connects the previous and current points,
and a ray connecting the current point and the test point. This is
explained much more clearly in a picture in Figure~\ref{fig8}.

\begin{figure}[htbp]
\setlength{\unitlength}{0.00063300in}%
%
\begingroup\makeatletter\ifx\SetFigFont\undefined
% extract first six characters in \fmtname
\def\x#1#2#3#4#5#6#7\relax{\def\x{#1#2#3#4#5#6}}%
\expandafter\x\fmtname xxxxxx\relax \def\y{splain}%
\ifx\x\y   % LaTeX or SliTeX?
\gdef\SetFigFont#1#2#3{%
  \ifnum #1<17\tiny\else \ifnum #1<20\small\else
  \ifnum #1<24\normalsize\else \ifnum #1<29\large\else
  \ifnum #1<34\Large\else \ifnum #1<41\LARGE\else
     \huge\fi\fi\fi\fi\fi\fi
  \csname #3\endcsname}%
\else
\gdef\SetFigFont#1#2#3{\begingroup
  \count@#1\relax \ifnum 25<\count@\count@25\fi
  \def\x{\endgroup\@setsize\SetFigFont{#2pt}}%
  \expandafter\x
    \csname \romannumeral\the\count@ pt\expandafter\endcsname
    \csname @\romannumeral\the\count@ pt\endcsname
  \csname #3\endcsname}%
\fi
\fi\endgroup
\begin{center}
\begin{picture}(5190,4996)(2639,-5062)
\put(4966,-1816){\circle*{108}}
\put(2701,-4201){\circle*{108}}
%\put(3316,-5701){\circle*{108}}
%\put(4816,-6301){\circle*{108}}
%\put(7066,-5911){\circle*{108}}
\put(7516,-4126){\circle*{108}}
\put(6751,-2251){\circle*{108}}
\put(3841,-4411){\circle*{108}}
\put(6247,-4804){\circle*{108}}
\put(4300,-3700){\circle*{108}}
\put(5221,-3241){\circle*{108}}
%\put(4840,-6025){\circle*{108}}
\put(7081,-3961){\circle*{108}}
\put(6404,-2587){\circle*{108}}
\put(4450,-2530){\circle*{108}}
%\put(5146,-4975){\circle*{108}}
\put(3766,-2401){\circle*{108}}
%\put(3505,-5440){\circle*{108}}
\put(4439,-1894){\makebox(0,0)[lb]{$\alpha_2$}}
\put(3331,-3385){\circle*{108}}
\put(7544,-4129){\line(-2, 5){762}}
\put(6782,-2224){\line(-4, 1){1804}}
\put(4978,-1773){\line(-2,-1){1200}}
\put(7334,-4354){\framebox(375,405){}}
\put(6584,-2449){\framebox(375,405){}}
\put(4784,-1999){\framebox(375,405){}}
\put(3569,-2599){\framebox(375,405){}}
\multiput(6779,-2224)(-26.66666,66.66666){11}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4994,-1774)(-71.14288,11.85712){11}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(7559,-4099){\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(7559,-4084)(0.00000,72.00000){11}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(7829,-4159){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}0}}}
\put(7109,-2314){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}1}}}
\put(4949,-1474){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}2}}}
\put(3299,-2434){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}3}}}
\put(6344,-1939){\makebox(0,0)[lb]{$\alpha_1$}}
\put(7334,-3379){\makebox(0,0)[lb]{$\alpha_0$}}
\end{picture}
\end{center}
\caption{Giftwrapping Algorithm}
\label{fig8}
\end{figure}

Each step of the algorithm takes $O(n)$ time, since $|S| = n$ and
potentially all points in $S$ must be examined to determine the one
with minimum angle. Since the algorithm does $h$ steps (one for each
point on the hull), the running time of the algorithm is $O(nh)$.  If
all the points lie on the hull then the algorithm will take $O(n^2)$
time.  Because of this worst-case behavior the algorithm is rarely
used in practice, although it is parallelizable and does generalize to
3D.

Note that this algorithm also requires additional primitives besides
the line-side and in-circle tests.

\subsubsection{Graham Scan}
\index{Convex Hull: Graham Scan}

The Graham Scan algorithm constructs the convex hull in two pieces
(the top and bottom), then connects them together. For each half, it
begins by sorting the points in decreasing order of their $x$
coordinates. Then, starting at the largest (or smallest) $x$, it scans
through the points, building the hull along the way. As it encounters
a new point, it connects it to the previous point, assuming that it
will be in the hull. If adding that point causes the hull to be
non-convex, the algorithm adds a {\it bridge} to the appropriate point
to ensure that the hull remains convex. This is illustrated in
Figure~\ref{fig9}; the dashed lines are the bridges added.

\begin{figure}[thbp]
\setlength{\unitlength}{0.00083300in}%
%
\begingroup\makeatletter\ifx\SetFigFont\undefined
% extract first six characters in \fmtname
\def\x#1#2#3#4#5#6#7\relax{\def\x{#1#2#3#4#5#6}}%
\expandafter\x\fmtname xxxxxx\relax \def\y{splain}%
\ifx\x\y   % LaTeX or SliTeX?
\gdef\SetFigFont#1#2#3{%
  \ifnum #1<17\tiny\else \ifnum #1<20\small\else
  \ifnum #1<24\normalsize\else \ifnum #1<29\large\else
  \ifnum #1<34\Large\else \ifnum #1<41\LARGE\else
     \huge\fi\fi\fi\fi\fi\fi
  \csname #3\endcsname}%
\else
\gdef\SetFigFont#1#2#3{\begingroup
  \count@#1\relax \ifnum 25<\count@\count@25\fi
  \def\x{\endgroup\@setsize\SetFigFont{#2pt}}%
  \expandafter\x
    \csname \romannumeral\the\count@ pt\expandafter\endcsname
    \csname @\romannumeral\the\count@ pt\endcsname
  \csname #3\endcsname}%
\fi
\fi\endgroup
\begin{center}
\begin{picture}(2753,1414)(4790,-4373)
\put(5461,-3552){\circle*{64}}
\put(5928,-4117){\circle*{64}}
\put(6246,-3870){\circle*{64}}
\put(6423,-3173){\circle*{64}}
\put(6881,-3693){\circle*{64}}
\put(7437,-4284){\circle*{64}}
\put(7146,-3155){\circle*{64}}
\put(7420,-4302){\line(-1, 4){282}}
\put(7138,-3174){\line(-1,-2){269}}
\put(6869,-3712){\line(-5, 6){450}}
\put(6419,-3172){\line(-1,-4){172}}
\put(6247,-3860){\line(-6,-5){330}}
\put(5917,-4135){\line(-4, 5){468}}
\put(5449,-3550){\line(-4,-5){488}}
\multiput(7146,-3155)(-72.30000,0.00000){11}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(4972,-4164){\circle*{64}}
\multiput(6423,-3155)(-68.57143,-27.42857){15}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(4790,-4249){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}7}}}
\multiput(6423,-3155)(-32.06667,-64.13333){16}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(7543,-4346){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}0}}}
\put(7243,-3164){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}1}}}
\put(6820,-3905){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}2}}}
\put(6361,-3067){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}3}}}
\put(6228,-4090){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}4}}}
\put(5858,-4302){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}5}}}
\put(5373,-3464){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}6}}}
\end{picture}
\end{center}
\caption{Graham Scan Algorithm}
\label{fig9}
\end{figure}

The running time for the Graham Scan algorithm is $T(n) = T($sort$) +
O(n)$, since the algorithm must take $n$ steps (one for each point in
the scan), and the amortized cost per step can be shown to be
$O(1)$. Since $T($sort$)$ is $\Omega(n \log n)$ for comparison-based
sorts, the Graham Scan algorithm is essentially an $O(n \log n)$
algorithm (unless the data lends itself to a non-comparison-based
sorting algorithm).

Finally, note that this algorithm is inherently sequential, and thus
does not parallelize well. Even on a sequential machine, it is not the
most efficient in practice. Furthermore, it does not generalize to
higher dimensions. However, it does have an interesting historical
note in that it was developed at Bell Labs in response to a real
application, and thus is one of the few hull algorithms that were not
theoretically motivated.

\subsubsection{Mergehull}
\index{Convex Hull: Mergehull}

The Mergehull algorithm is a divide-and-conquer algorithm that is
similar in structure to merge sort.  Mergehull has a trivial divide
phase, and does most of its work during the merge phase.

Mergehull's divide stage is simple: the algorithm merely splits the
input set $S$ into two subsets $S_1$, $S_2$ by one of the
coordinates. It then recursively solves $S_1$ and $S_2$ (e.g., it
computes the convex hull for both $S_1$ and $S_2$).

After the recursion returns, the algorithm must merge the two convex
hulls for $S_1$ and $S_2$ into a single convex hull for
$S$. Conceptually, this involves locating two {\it bridges} between
the hulls---line segments connecting two points in $S_1$ to two points
in $S_2$. Once the bridges are added to connect the two hulls, the
interior segments of the subhulls are discarded. See
Figure~\ref{fig10} for an illustration. Finding the bridges and
discarding the inside edges requires an Overmars dual binary search,
which takes time $O(\log n)$. Note that it is not sufficient just to
place the bridges between the minima and maxima of the two sub-hulls.

\begin{figure}[htbp]
\setlength{\unitlength}{0.00083300in}%
%
\begingroup\makeatletter\ifx\SetFigFont\undefined
% extract first six characters in \fmtname
\def\x#1#2#3#4#5#6#7\relax{\def\x{#1#2#3#4#5#6}}%
\expandafter\x\fmtname xxxxxx\relax \def\y{splain}%
\ifx\x\y   % LaTeX or SliTeX?
\gdef\SetFigFont#1#2#3{%
  \ifnum #1<17\tiny\else \ifnum #1<20\small\else
  \ifnum #1<24\normalsize\else \ifnum #1<29\large\else
  \ifnum #1<34\Large\else \ifnum #1<41\LARGE\else
     \huge\fi\fi\fi\fi\fi\fi
  \csname #3\endcsname}%
\else
\gdef\SetFigFont#1#2#3{\begingroup
  \count@#1\relax \ifnum 25<\count@\count@25\fi
  \def\x{\endgroup\@setsize\SetFigFont{#2pt}}%
  \expandafter\x
    \csname \romannumeral\the\count@ pt\expandafter\endcsname
    \csname @\romannumeral\the\count@ pt\endcsname
  \csname #3\endcsname}%
\fi
\fi\endgroup
\begin{center}
\begin{picture}(3668,3447)(4306,-4879)
\put(4479,-2906){\circle*{70}}
\put(4349,-3623){\circle*{70}}
\put(5073,-4107){\circle*{70}}
\put(5442,-3920){\circle*{70}}
\put(5554,-3196){\circle*{70}}
\put(5453,-2505){\circle*{70}}
\put(4985,-2117){\circle*{70}}
\put(6707,-2130){\circle*{70}}
\put(4675,-2519){\circle*{70}}
\put(6760,-3515){\circle*{70}}
\multiput(5073,-4108)(71.41935,0.00000){32}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(7808,-3437){\circle*{70}}
\put(7931,-2917){\circle*{70}}
\put(7702,-2378){\circle*{70}}
\put(7296,-4117){\circle*{70}}
\put(6696,-2123){\line( 4,-1){1008}}
\put(7704,-2375){\line( 2,-5){224}}
\put(7928,-2935){\line(-1,-4){128}}
\put(7800,-3447){\line(-3,-4){506.880}}
\put(7296,-4125){\line(-5, 6){524.836}}
\multiput(6757,-3507)(-8.00000,8.00000){2}{\makebox(6.6667,10.0000){\SetFigFont{7}{8.4}{rm}.}}
\put(6749,-3499){\line( 0, 1){1376}}
\put(4984,-2105){\line( 6,-5){468}}
\put(5452,-2495){\line( 1,-6){114}}
\put(5566,-3179){\line(-1,-6){123}}
\put(5443,-3917){\line(-2,-1){370}}
\put(5073,-4102){\line(-3, 2){732}}
\put(4341,-3614){\line( 1, 5){145}}
\put(4486,-2889){\line( 1, 2){189}}
\put(4675,-2511){\line( 3, 4){306.120}}
\put(6105,-1444){\line( 0,-1){3423}}
\put(6200,-4800){\it Split}
\multiput(5002,-2114)(71.33333,0.00000){25}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\end{picture}
\end{center}
\caption{Mergehull Algorithm}
\label{fig10}
\end{figure}

If we assume that the points are presorted (to make the subdivision in
the divide phase take constant time) and that tree data structures are used
throughout, the recurrence for the mergehull algorithm is:
\begin{eqnarray*}
T(n) &=& 2 T\left(\frac{n}{2}\right) + O(\log n) \\
     &=& O(n)
\end{eqnarray*}
When added to the time of the original sort, the total time is $O(n \log n)$.

It is interesting to note that the mergehull algorithm parallelizes,
and that the parallel algorithm does well in practice (it obtains a
nearly linear speedup with the number of processors). The
parallelization is rather simple: the set $S$ of input points is first
divided into $P$ regions (rather than in half, as above). Each
of the $P$ processors is then assigned one of the regions, for which
it sequentially generates a convex hull. Then, the $P$ regions are
merged by finding bridges in parallel.

Finally, the mergehull algorithm can be generalized to three or more
dimensions. As the dimensionality increases, however, the merge
becomes more difficult. For example, in 3D, the bridges between
subhulls are no longer line segments, but instead are surfaces that
connect the hulls. Finding these surfaces is more expensive: in 3D,
the merge takes $O(n)$ time, so the total time for mergehull is $O(n
\log n)$, even when the input points are presorted. 
% In general, in $d$ dimensions, $T(n) = O(n \log^{d-2} n)$.

\subsubsection{Quickhull}
\index{Convex Hull: Quickhull}

The Quickhull algorithm is another variant on the recursive
divide-and-conquer method for finding convex hulls. Unlike the
mergehull algorithm that we have just seen, quickhull concentrates its
work on the divide phase of the algorithm, and uses a trivial merge
step. It is thus very similar in structure to the quicksort sorting
algorithm, and we will see below that it has similar time complexity.

Quickhull works (in 2D) by locating three points $A$, $B$, and $C$ on
the convex hull of its input $S$ (with $C$ above both $A$ and $B$),
then recursively refining the hull between $AC$ and $CB$. Since both
halves of the recursion meet at $C$, which is on the hull, the merge
is essentially a no-op. When the algorithm completes, the upper
half-hull between $A$ and $B$ will be complete. The algorithm is then
repeated for the lower half-hull. Figure~\ref{fig11} details the
process graphically.

\begin{figure}[htbp]
\setlength{\unitlength}{0.00083300in}%
%
\begingroup\makeatletter\ifx\SetFigFont\undefined
% extract first six characters in \fmtname
\def\x#1#2#3#4#5#6#7\relax{\def\x{#1#2#3#4#5#6}}%
\expandafter\x\fmtname xxxxxx\relax \def\y{splain}%
\ifx\x\y   % LaTeX or SliTeX?
\gdef\SetFigFont#1#2#3{%
  \ifnum #1<17\tiny\else \ifnum #1<20\small\else
  \ifnum #1<24\normalsize\else \ifnum #1<29\large\else
  \ifnum #1<34\Large\else \ifnum #1<41\LARGE\else
     \huge\fi\fi\fi\fi\fi\fi
  \csname #3\endcsname}%
\else
\gdef\SetFigFont#1#2#3{\begingroup
  \count@#1\relax \ifnum 25<\count@\count@25\fi
  \def\x{\endgroup\@setsize\SetFigFont{#2pt}}%
  \expandafter\x
    \csname \romannumeral\the\count@ pt\expandafter\endcsname
    \csname @\romannumeral\the\count@ pt\endcsname
  \csname #3\endcsname}%
\fi
\fi\endgroup
\begin{center}
\begin{picture}(3353,1465)(4411,-3630)
\put(4728,-3349){\circle*{76}}
\put(7561,-3349){\circle*{76}}
\put(7446,-2881){\circle*{64}}
\put(6281,-2608){\circle*{64}}
\put(5002,-2784){\circle*{64}}
\put(5276,-3212){\circle*{64}}
\put(5410,-2936){\circle*{64}}
\put(5821,-3255){\circle*{64}}
\put(5814,-2731){\circle*{64}}
\put(6534,-3125){\circle*{64}}
\put(5646,-2440){\circle*{76}}
\put(6898,-3591){\circle*{64}}
\put(5584,-2273){\makebox(0,0)[lb]{$C$}}
\put(5478,-3557){\circle*{64}}
\put(4728,-3361){\line( 1, 1){919.500}}
\put(5646,-2440){\line( 2,-1){1922}}
\multiput(4746,-3361)(71.71795,0.00000){40}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5655,-2440)(0.00000,-70.84615){14}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(5725,-3000){\makebox(0,0)[l]{$d_{\mbox{max}}$}}
\put(4623,-3464){\framebox(229,212){}}
\put(5540,-2555){\framebox(229,212){}}
\put(7464,-3473){\framebox(229,212){}}
\put(4411,-3411){\makebox(0,0)[lb]{$A$}}
\put(7764,-3420){\makebox(0,0)[lb]{$B$}}
\end{picture}
\end{center}
\caption{Quickhull Algorithm}
\label{fig11}
\end{figure}

The tricky part in Quickhull is implementing the divide stage. A
typical algorithm (which works well in the average case but does not
guarantee the best time bound) is the following. We start out with two
points $A$, $B$ on the convex hull. These can be found on the first
iteration by taking the points in $S$ with minimum and maximum $x$
coordinate. Then, apply the following algorithm:
\begin{enumerate}
\item Find the furthest point $C$ from the line joining $A$ and
$B$. To do this, the line-side test (described in
Section~\ref{lineside}) can be used, as it returns the area of the
triangle $ABC$, and the triangle with largest area is the one that has
the furthest point $C$. Note that, when building the upper half-hull,
$C$ should be above $AB$, whereas $C$ should be above the inverse
line $BA$ when building the lower half-hull.
\item Discard all points in triangle $ABC$, as they clearly
will not be on the boundary of the final hull (since $ABC$ is contained within the
hull).
\item Recurse on $A,C$ and $C,B$ to fully refine the hull between $A$
and $C$, and between $C$ and $B$.
\item Connect the two resulting hulls at $C$.
\item Repeat for the lower half-hull, and connect it at $A$ and $B$.
\end{enumerate}

The worst-case running time of the quickhull algorithm is $O(n^2)$,
for if the splits (choices of $C$) are bad, the algorithm may only
eliminate one point on each recursion, resulting in a recursion depth
of $n$ and a per-recursion cost of $n$. The input distribution that
forces this worst-case time is if all the input points are distributed
along a semicircle with exponentially decreasing distance from the
left-hand point (see Figure~\ref{fig12}).

\begin{figure}[hbtp]
\centerline{\psfig{figure=../lecture12/badqhull.eps,width=3in}}
\caption{Worst-case Input Distribution for Quickhull}
\label{fig12}
\end{figure}

Despite this poor worst-case bound, quickhull runs in $O(n \log n)$
time for most ``realistic'' point sets.  In fact, in practice it can
often run as fast as $O(n)$, since the algorithm discards points in
step 2, above, reducing the size of the input set for later
recursions.

There is also a technique developed by Kirkpatrick and Seidel that
locates the ideal ``pivot'' point $C$ using two-variable linear
programming.  This variant guarantees the optimal time complexity of $O(n
\log h)$.  
The technique is conceptually similar to the method of using the true
sample median to force quicksort to run in $\Theta(n \log n)$. The
technique involves finding the bridge segment that connects across
the median $x$ coordinate before recursively solving the subhulls
themselves.  This is the same bridge as found by mergehull, but instead
of requiring the solved subhulls for each half it uses linear-programming
to find the bridge directly in $O(n)$ time.

Finally, note that quickhull parallelizes easily, and generalizes to
three dimensions (in which case the running time is $O(n \log^2 h)$
for an algorithm by Edelsbrunner and Shi).

\section{Delaunay Triangulations and Voronoi Diagrams}

We will not consider the problem of generating 
Delaunay triangualations, but we first consider the dual problem
of generating Voronoi Diagrams.

\subsection{Voronoi Diagrams}
A {\it Voronoi diagram} (or Dirichlet tessellation) is a partition of
space into regions $V(p_i)$ based on a set of input points $p_i \in S$ such
that:
\[
V(p_i) = \{ p: |p_i - p| \leq |p_j - p|, \forall \, j \neq i\}.
\]
\noindent 
where $|p_i - p|$ indicates the distance from $p_i$ to $p$ (here we 
assume this distance is Euclidean, but later we will briefly discuss
other distance metrics).
What this property means is that, for a point $p_i$, its
Voronoi region $V(p_i)$ will contain all points that are closer to
$p_i$ than to any other point in $S$.

Voronoi diagrams are well-defined in any dimension. In two dimensions,
the Voronoi regions are convex polygons; in three dimensions, they are
convex polyhedra, and so forth. Figure~\ref{fig13} gives an example of
a 2D Voronoi diagram.

\begin{figure}[htbp]
\begin{center}
\setlength{\unitlength}{0.00083300in}%
\begin{picture}(2219,1525)(3588,-3230)
\thicklines
\put(5761,-2114){\circle*{76}}
\put(4781,-1973){\circle*{76}}
\put(5046,-3093){\circle*{76}}
\put(3634,-2652){\circle*{76}}
\put(4455,-2652){\line(-1,-2){283}}
\put(3970,-1946){\line( 2,-3){476}}
\put(4446,-2660){\line( 3, 1){741}}
\put(5187,-2413){\line( 1, 4){174}}
\put(5196,-2414){\line( 2,-1){592}}
\end{picture}
\end{center}
\caption{A 2D Voronoi Diagram}
\label{fig13}
\end{figure}

Note that there are actually many different ways to represent the
diagram.  One way is to return for each original point the set of Voronoi
neighbors (points that share a Voronoi edge).  Another is to return
the graph structure formed by the Voronoi edges.

\subsection{Delaunay Triangulations}
\begin{figure}[pt]
\begin{center}
\setlength{\unitlength}{0.00083300in}%
%
\begingroup\makeatletter\ifx\SetFigFont\undefined
% extract first six characters in \fmtname
\def\x#1#2#3#4#5#6#7\relax{\def\x{#1#2#3#4#5#6}}%
\expandafter\x\fmtname xxxxxx\relax \def\y{splain}%
\ifx\x\y   % LaTeX or SliTeX?
\gdef\SetFigFont#1#2#3{%
  \ifnum #1<17\tiny\else \ifnum #1<20\small\else
  \ifnum #1<24\normalsize\else \ifnum #1<29\large\else
  \ifnum #1<34\Large\else \ifnum #1<41\LARGE\else
     \huge\fi\fi\fi\fi\fi\fi
  \csname #3\endcsname}%
\else
\gdef\SetFigFont#1#2#3{\begingroup
  \count@#1\relax \ifnum 25<\count@\count@25\fi
  \def\x{\endgroup\@setsize\SetFigFont{#2pt}}%
  \expandafter\x
    \csname \romannumeral\the\count@ pt\expandafter\endcsname
    \csname @\romannumeral\the\count@ pt\endcsname
  \csname #3\endcsname}%
\fi
\fi\endgroup
\begin{picture}(2747,2477)(3642,-3793)
\put(4280,-3234){\circle*{76}}
\put(5198,-3463){\circle*{76}}
\put(5930,-3005){\circle*{76}}
\put(5390,-2440){\circle*{76}}
\put(6105,-2017){\circle*{76}}
\put(5064,-1664){\circle*{76}}
\put(3811,-2573){\circle*{76}}
\put(4107,-1675){\circle*{76}}
\put(4702,-3014){\circle*{76}}
\put(4693,-3005){\circle*{76}}
\put(4525,-2299){\line( 6,-1){858}}
\put(5383,-2442){\line(-6,-5){687.639}}
\put(4702,-3023){\line(-1, 4){185.177}}
\put(4543,-2273){\line(-1,-4){238}}
\put(4305,-3225){\line( 5, 3){390}}
\put(4711,-3031){\line( 6,-5){522}}
\put(5233,-3466){\line(-4, 1){948}}
\put(4711,-3025){\line( 1, 0){1224}}
\put(5935,-3025){\line(-5,-3){705}}
\put(5399,-2458){\line( 1,-1){556}}
\put(4534,-2290){\circle*{76}}
\put(5417,-2440){\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4605,-3181)(0.00000,-75.00000){9}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(5946,-2987){\line( 1, 6){165}}
\put(6111,-1997){\line(-5,-3){705}}
\put(6123,-2017){\line(-3, 1){1050}}
\put(5073,-1667){\line( 2,-5){306}}
\put(5055,-1673){\line(-1, 0){944}}
\put(4111,-1673){\line(-1,-3){303}}
\put(3808,-2582){\line( 3,-4){483}}
\put(3802,-2581){\line( 5, 2){740}}
\put(4543,-2290){\line( 5, 6){535}}
\put(4543,-2281){\line(-3, 4){462}}
\multiput(4208,-2696)(-48.12500,-57.75000){9}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4128,-2140)(-67.71429,11.28571){8}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4120,-2140)(12.00000,-72.00000){9}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4216,-2716)(72.75000,-24.25000){5}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4507,-2813)(52.85714,42.28571){8}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4877,-2517)(22.14286,66.42857){8}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5032,-2052)(-66.00000,33.00000){8}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4111,-2158)(57.50000,43.12500){9}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4571,-1813)(0.00000,69.28571){8}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5505,-1981)(26.00000,65.00000){9}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4508,-2802)(16.50000,-66.00000){7}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4607,-3198)(71.30000,0.00000){11}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5320,-3198)(0.00000,73.33333){7}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5320,-2758)(-62.85714,37.71429){8}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5328,-2749)(63.12500,25.25000){9}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5833,-2547)(-37.66667,62.77778){10}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5494,-1982)(-66.85714,-11.14286){8}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5849,-2555)(75.42857,-12.57143){8}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5328,-3181)(64.88889,-32.44444){10}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(5912,-3473)(-1.00000,0.00000){2}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\end{picture}
\end{center}
\caption{A Delaunay Triangulation (and associated Voronoi diagram)}
\label{fig14}
\end{figure}

The dual graph of a Voronoi diagram is called a {\it Delaunay
triangulation}. Given a Voronoi diagram, the Delaunay triangulation is
formed by connecting the points at the centers of every pair of
neighboring Voronoi regions. See Figure~\ref{fig14} for an example.
Note that the the Delaunay triangulation will always consist of
triangles, since the Delaunay edges join neighboring Voronoi regions,
and two neighbor regions of adjacent Voronoi edges will themselves be
neighbors. Note that a Delaunay edge is always perpendicular to the
Voronoi edge that the two neighbors share, but that the Delaunay edge
does not necessarily cross the corresponding Voronoi edge, as
illustrated in Figure~\ref{fig15}.

\begin{figure}[pt]
\begin{center}
\setlength{\unitlength}{0.00083300in}%
%
\begingroup\makeatletter\ifx\SetFigFont\undefined
% extract first six characters in \fmtname
\def\x#1#2#3#4#5#6#7\relax{\def\x{#1#2#3#4#5#6}}%
\expandafter\x\fmtname xxxxxx\relax \def\y{splain}%
\ifx\x\y   % LaTeX or SliTeX?
\gdef\SetFigFont#1#2#3{%
  \ifnum #1<17\tiny\else \ifnum #1<20\small\else
  \ifnum #1<24\normalsize\else \ifnum #1<29\large\else
  \ifnum #1<34\Large\else \ifnum #1<41\LARGE\else
     \huge\fi\fi\fi\fi\fi\fi
  \csname #3\endcsname}%
\else
\gdef\SetFigFont#1#2#3{\begingroup
  \count@#1\relax \ifnum 25<\count@\count@25\fi
  \def\x{\endgroup\@setsize\SetFigFont{#2pt}}%
  \expandafter\x
    \csname \romannumeral\the\count@ pt\expandafter\endcsname
    \csname @\romannumeral\the\count@ pt\endcsname
  \csname #3\endcsname}%
\fi
\fi\endgroup
\begin{picture}(1428,1140)(4001,-2792)
\thicklines
\put(4824,-1727){\circle*{70}}
\put(5170,-2193){\circle*{70}}
\put(4834,-2749){\circle*{70}}
\put(4825,-1725){\line( 3,-4){357.480}}
\put(5170,-2211){\line(-3,-5){336.088}}
\put(4827,-2767){\line( 0, 1){1042}}
\multiput(5417,-1664)(-60.46154,-40.30769){14}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4631,-2188)(60.46154,-40.30769){14}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\multiput(4639,-2184)(-69.55556,0.00000){10}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\end{picture}
\end{center}
\caption{Example of a Delaunay edge that does not cross its Voronoi edge}
\label{fig15}
\end{figure}

Finally, note that degeneracies can be encountered in forming a
Delaunay triangulation. In particular, if two Voronoi regions
intersect only at a corner, there will be more than one possible
triangulation (recall that a triangulation contains only
non-overlapping triangles). This is illustrated in Figure~\ref{fig16}.

\begin{figure}[phtb]
\begin{center}
\setlength{\unitlength}{0.00083300in}%
%
\begingroup\makeatletter\ifx\SetFigFont\undefined
% extract first six characters in \fmtname
\def\x#1#2#3#4#5#6#7\relax{\def\x{#1#2#3#4#5#6}}%
\expandafter\x\fmtname xxxxxx\relax \def\y{splain}%
\ifx\x\y   % LaTeX or SliTeX?
\gdef\SetFigFont#1#2#3{%
  \ifnum #1<17\tiny\else \ifnum #1<20\small\else
  \ifnum #1<24\normalsize\else \ifnum #1<29\large\else
  \ifnum #1<34\Large\else \ifnum #1<41\LARGE\else
     \huge\fi\fi\fi\fi\fi\fi
  \csname #3\endcsname}%
\else
\gdef\SetFigFont#1#2#3{\begingroup
  \count@#1\relax \ifnum 25<\count@\count@25\fi
  \def\x{\endgroup\@setsize\SetFigFont{#2pt}}%
  \expandafter\x
    \csname \romannumeral\the\count@ pt\expandafter\endcsname
    \csname @\romannumeral\the\count@ pt\endcsname
  \csname #3\endcsname}%
\fi
\fi\endgroup
\begin{picture}(3118,2724)(3493,-4570)
\put(4801,-2761){\circle*{64}}
\put(4801,-3961){\circle*{64}}
\put(6001,-3961){\circle*{64}}
\multiput(4808,-2758)(50.37500,-50.37500){25}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\thicklines
\put(5399,-2158){\line( 0,-1){2400}}
\put(4199,-3358){\line( 1, 0){2400}}
\multiput(4817,-3967)(51.78261,51.78261){24}{\makebox(6.6667,10.0000){\SetFigFont{10}{12}{rm}.}}
\put(6001,-2761){\circle*{64}}
\thinlines
\put(4808,-3958){\framebox(1200,1200){}}
\put(3000,-3100){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}Delaunay Edge}}}
\put(5932,-1879){\vector(-1,-1){477}}
\put(6466,-2661){\vector(-1,-1){406}}
\put(4367,-3075){\vector( 1, 0){661}}
\put(5964,-1990){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}Voronoi Edge}}}
\put(6537,-2749){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}Delaunay Edge}}}
\put(3250,-2875){\makebox(0,0)[lb]{\smash{\SetFigFont{12}{14.4}{rm}Potential}}}
\end{picture}
\end{center}
\caption{A degenerate case for Delaunay triangulation}
\label{fig16}
\end{figure}

\subsection{Properties of Delaunay Triangulations}
Delaunay triangulations have certain ``nice'' properties, including:
\begin{itemize}
\item The Delaunay triangulation of a non-degenerate set of points is
unique. Thus there is only one triangulation that obeys the following
set of properties.
\item A circle through the three points of a Delaunay triangle
contains {\em no} other points. Also, the intersection of the three
corresponding Voronoi edges is at the center of the circle, as shown
in Figure~\ref{fig17}. This property generalizes to spheres in 3D and
so forth.

\begin{figure}[phbt]
\centerline{\psfig{figure=../lecture12/delcircle.eps,width=3in}}
\caption{Circle test for Delaunay triangles}
\label{fig17}
\end{figure}

\item The minimum angle across all the angles in all the triangles in
a Delaunay triangulation is greater than the minimum angle in any
other triangulation of the same points.
\end{itemize}

\noindent These properties tend to result in triangles that are
``fat'' (closer to equilateral) rather than ``skinny''. This is a
property that is highly desirable in applications, especially finite
element simulations where the geometry of the finite elements affects
the accuracy and quality of the simulation.

\subsection{Algorithms for Delaunay Triangulation}
The Delaunay triangulation for a set of points can be found using one
of several types of algorithms: a naive direct algorithm, by reducing
the problem to a convex hull problem and applying a convex hull
algorithm, and by using a direct algorithm. We will consider each of
these cases in turn.

\subsubsection{Naive Algorithm}
A naive way of finding the Delaunay triangulation is the
following: for every three-point-subset of the input set, check if any
other point is in its circle. Iff not, those three points define a
Delaunay triangle.

This algorithm runs in $O(n^4)$ time, since there are $n^3$ triples of
points to check, and each check takes $O(n)$ time (since the in-circle
test must be applied $n$ times).  Although polynomial-time, this
algorithm is quite impractical since, as we will see, there
are algorithms that run in $O(n \log n)$ time (and the constant factors
in the big-O are quite small).

\subsubsection{Solution via Reduction to 3D Convex Hull}

A more practical method for finding the Delaunay triangulation of a
set of $d$-dimensional points uses a reduction of the Delaunay
triangulation problem to the problem of finding a convex hull in $d+1$
dimensions. For the typical 2D triangulation case, a 3D convex hull is
used.

The reduction is as follows.  The original input points $P$ are
projected upward onto a three-dimensional parabolic surface centered
at the origin. Then, the 3D convex hull of those projected points is
found. When the hull is represented as a neighbor graph and projected
back down to the two-dimensional plane, the result is the Delaunay
triangulation of $P$. More precisely, given $P$, we first form:
\[
P' = \left\{ (x,y,x^2+y^2) : (x,y) \in P \right\}.
\]
\noindent Then let $C=$LowerConvexHull$(P')$, represented as a
neighbor graph. Then, the neighbor graph of $C$ is equivalent to
Delaunay($P$).

The reason that this works is related to the properties of the
Delaunay triangulation, in particular, the property that states that
any circle through three points of the input set containing no other
points circumscribes a Delaunay triangle through those three
points. Recall that any slice of a 3D parabolic surface maps onto a
circle in the $x$-$y$ plane. Also, notice that any face of the 3D
convex hull defines a slice through the parabola. Thus the projection
of this slice onto the plane is a circle, and the circle will pass
through the projection of the three points on the plane. Since the
slice is only for a single hull face, there will be no other projected
points in that circle, so the projection of the hull face defines a
Delaunay triangle.  Since this works for all hull faces, the result
faces of the convex hull in 3-D are the Delaunay triangles in 2-D (when
projected down).

Using this reduction to find the Delaunay triangulation gives an
efficient algorithm that runs in $O(n \log n)$. However, there are
still better algorithms in practice, such as some of the direct
algorithms described in the next section.

\subsubsection{Direct Algorithms}

The following algorithms can be used to find either Voronoi diagrams
or Delaunay triangulations, as the two problems are duals of each
other. Those algorithms marked with a dagger ($\dagger$) are easy to
parallelize.
\begin{itemize}
\item Intersecting Halfspaces$^\dagger$: $O(n^2 \log n)$

This algorithm works as follows: for each input point, draw a line to
every other point. Generate a halfspace separating the input point
from every other point (by bisecting the lines drawn in the first
step). Compute the halfspace intersection of these halfspaces (a
convex hull problem running in $O(n \log n)$). This gives a Voronoi
region. Repeat for all $n$ points.
\item Radial Sweep
\item Incremental Construction: $O(n \log n)$ (expected time)
\item Plane Sweep: $O(n \log n)$
\item Divide and Conquer Algorithms
    \begin{itemize}
    \item Merge-based: $O(n \log n)$
    \item Split-based$^\dagger$: $O(n \log n)$
    \end{itemize}
\end{itemize}

In this class we will only discuss the incremental and the split-based
divide-and-conquer algorithms.

