delaunay.tex

国外免费地震资料处理软件包

\append{Incremental DELAUNAY TRIANGULATION and related problems}

Delaunay triangulation \cite[]{delaunay2,sibson,stolfi} is a fundamental
geometric construction, which has numerous applications in different
computational problems. For a given set of nodes (points on the
plane), Delaunay triangulation constructs a triangle tessellation of
the plane with the initial nodes as vertices. Among all possible
triangulations, the Delaunay triangulation possesses optimal
properties, which make it very attractive for practical applications,
such as computational mesh generation. One of the most well-known
properties is maximizing the minimum triangulation angle. In three
dimensions, Delaunay triangulation generalizes naturally to a
tetrahedron tessellation.
\par
Several optimal-time algorithms of Delaunay triangulation (and its
counterpart--Voronoi diagram) have been proposed in the literature.
The divide-and-conquer algorithm \cite[]{shamos,stolfi} and the
sweep-line algorithm \cite[]{fortune} both achieve the optimal $O (N
\log N)$ worst-case time complexity. Alternatively, a family of
incremental algorithms has been used in practice because of their
simplicity and robustness. Though the incremental algorithm can take
$O (N^2)$ time in the worst case, the expectation time can still be $O
(N \log N)$, provided that the nodes are inserted in a random order
\cite[]{knuth}.
\par
The incremental algorithm consists of two main parts:
 \begin{enumerate}
 \item Locate a triangle (or an edge), containing the inserted point.
 \item Insert the point into the current triangulation, making the
   necessary adjustments.
 \end{enumerate}
\par
The Delaunay criterion can be reduced in the second step to a simple
\emph{InCircle} test \cite[]{stolfi}: if a circumcircle of a triangle
contains another triangulation vertex in its circumcenter, then the
edge between those two triangles should be ``flipped'' so that two new
triangles are produced. The testing is done in a recursive fashion
consistent with the incremental nature of the algorithm. When a new
node is inserted inside a triangle, three new triangles are created,
and three edges need to be tested. When the node falls on an edge,
four triangles are created, and four edges are tested. In the case of
test failure, a pair of triangles is replaced by the flip operation
with another pair, producing two more edges to test. Under the
randomization assumption, the expected total time of point insertion
is $O (N)$.  Randomization can be considered as an external part of
the algorithm, provided by preprocessing.
\par
\cite{knuth} reduce the point location step to an efficient $O (N
\log N)$ procedure by maintaining a hierarchical tree structure: all
triangles, occurring in the incremental triangulation process, are
kept in memory, associated with their ``parents.'' One or two point
location tests (\emph{CCW} tests) are sufficient to move to a lower
level of the tree. The search terminates with a current Delaunay
triangle.
\par
To test the algorithmic performance of the incremental construction, I
have profiled the execution time of my incremental triangulation
program with the Unix \texttt{pixie} utility. The profiling result,
shown in Figures~\ref{fig:itime} and~\ref{fig:ctime}, complies
remarkably with the theory: $O (N \log N)$ operations for the point
location step, and $O (N)$ operations for the point insertion step.
The experimental constant for the insertion step time is about $8.6$.
The experimental constant for the point location step is $4$.  The CPU
time, depicted in Figure~\ref{fig:time}, also shows the expected $O (N
\log N)$ behavior.

\inputdir{.}

\sideplot{itime}{width=2.5in}{The number of point insertion operations
(\emph{InCircle} test) plotted against the number of points.}

\sideplot{ctime}{width=2.5in}{Number of point location operations
  (\emph{CCW} test) plotted against the number of points.}

\sideplot{time}{width=2.5in}{CPU time (in seconds per point) plotted
  against the number of points.}
\par
A straightforward implementation of Delaunay triangulation would
provide an optimal triangulation for any given set of nodes. However,
the quality of the result for unfortunate geometrical
distributions of the nodes can be unsatisfactory. In the rest of this
appendix, I describe three problems, aimed at improving the
triangulation quality: conforming triangulation, triangulation of
height fields, and mesh refinement.  Each of these problems can be
solved with a variation of the incremental algorithm.

\subsection{Conforming Triangulation}

\inputdir{tri}

In the practice of mesh generation, the input nodes are often
supplemented by boundary edges: geologic interfaces, seismic rays, and
so on. It is often desirable to preserve the edges so that they appear
as edges of the triangulation \cite[]{SEG-1994-0502}. One possible approach is
\emph{constrained} triangulation, which preserves the edges, but only
approximately satisfies the Delaunay criterion \cite[]{lee,chew}. An
alternative, less investigated, approach is \emph{conforming}
triangulation, which preserves the ``Delaunayhood'' of the
triangulation by adding additional nodes \cite[]{hansen} (Figure
\ref{fig:conform}).  Conforming Delaunay triangulations are difficult
to analyze because of the variable number of additional nodes. This
problem was attacked by \cite{edels}, who suggested an algorithm
with a defined upper bound on added points. Unfortunately,
Edelsbrunner's algorithm is slow in practice because the number of
added points is largely overestimated.  I chose to implement a
modification of the simple incremental algorithm of Hansen and Levin.
Although Hansen's algorithm has only a heuristic justification and
sets no upper bound on the number of inserted nodes, its simplicity is
attractive for practical implementations, where it can be easily
linked with the incremental algorithm of Delaunay triangulation.
\par
The incremental solution to the problem of conforming triangulation
can be described by the following scheme:
 \begin{itemize}
 \item First, the boundary nodes are triangulated.
 \item Boundary edges are inserted incrementally.
 \item If a boundary edge is not present in the triangulations, it is
   split in half, and the middle node is inserted into the triangulation. This
   operation is repeated for the two parts of the original boundary
   edge and continues recursively until all the edge parts 
   conform.
 \item If at some point during the incremental process, a boundary edge
   violates the Delaunay criterion (the \emph{InCircle} test), it is
   split to assure the conformity.
 \end{itemize}

\plot{conform}{width=4in,height=2in}{An illustration of conforming triangulation.
  The left plot shows a triangulation of 500 random points; the
  triangulation in the right plot is conforming to the embedded
  boundary.  Conforming triangulation is a genuine Delaunay
  triangulation, created by adding additional nodes to the original
  distribution.}
\par
To insert an edge $AB$ into the current triangulation, I use the
following recursive algorithm:
 \begin{quote}
 Function \textbf{InsertEdge} ($AB$)
 \begin{enumerate}
 \item Define $C$ to be the midpoint of $AB$.
 \item Using the triangle tree structure, locate triangle $\mathcal{T} = DEF$
   that contains $C$ in the current triangulation.
 \item \textbf{If} $AB$ is an edge of $\mathcal{T}$ \textbf{then return}.
 \item \textbf{If} $A$ (or $B$) is a vertex of $\mathcal{T}$ (for example, $A = D$)
   {\bf then} define $C$ as an intersection of $AB$ and $EF$.
 \item {\bf Else} define $C$ as an intersection of $AB$ and an
   arbitrary edge of $\mathcal{T}$ (if such an intersection exists).
 \item Insert $C$ into the triangulation.
 \item {\bf InsertEdge} ($CA$).
 \item {\bf InsertEdge} ($CB$).
 \end{enumerate}
 \end{quote}
\par
The intersection point  of edges $AB$ and $EF$ is given by the formula
\begin{equation}