delaunay.tex

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

  C = A + \lambda (B-A)\;,
\end{equation}
where
\begin{equation}
  \lambda = \frac{(F_y - E_y)\,(E_x - A_x) - (F_x - E_x)\,(E_y-A_y)}{
    \det \left|\begin{array}{cc}
        B_x - A_x & B_y - A_y \\
        F_x - E_x & F_y - E_y
    \end{array}\right|}\;.
\end{equation}
The value of $\lambda$ should range between $0$ and $1$.
\par
If, at some stage of the incremental construction, a boundary edge
$AB$ fails the Delaunay \emph{InCircle} test for the circle $CABD$,
then I simply split it into two edges by adding the point of
intersection into the triangulation.  The rest of the process is very
much like the process of edge validation in the original incremental
algorithm.

\subsection{Triangulation of Height Fields}

Often, a velocity field (or other object that we want to triangulate)
is defined on a regular Cartesian grid. One way to perform a
triangulation in this case is to select a smaller subset of the
initial grid points, using them as the input to a triangulation
program. We need to select the points in a way that preserves the main
features of the original image, while removing some unnecessary
redundancy in the regular grid description.

\plot{sphere}{width=6in,height=2in}{Illustration of Garland's
  algorithm for triangulation of height fields. The left plot shows
  the input image of a sphere, containing 100 by 100 pixels. The
  middle plot shows 500 pixels, selected by the algorithm and
  triangulated. The right plot is the result of local plane
  interpolation of the triangulated surface.}
\par
\cite{height} surveyed different approaches
to this problem and proposed a fast version of the incremental
\emph{greedy insertion} algorithm. Their algorithm adds points
incrementally, selecting at each step the point with the maximum
interpolation error with respect to the current triangulation. Though
a straightforward implementation of this idea would lead to an
unacceptably slow algorithm, Garland and Heckbert have discovered
several sources for speeding it up. First, we can take advantage of
the fact that only a small area of the current triangulation gets
updated at each step. Therefore, it is sufficient to recompute the
interpolation error only inside this area. Second, the maximum
extraction procedure can be implemented very efficiently with a
priority queue data structure.

\inputdir{.}

\sideplot{opengl}{width=2in,height=2in}{An image from the previous
  example, as rendered by the OpenGL library. The shades on polygonal
  (triangulated) sides are exaggerated by a simulation of the direct
  light source.}
\par
Figure \ref{fig:sphere} illustrates this algorithm with a simple
example. The original image (the left plot) contained 10,000 points,
laid out on a regular rectangular grid. The algorithm selects a
smaller number of points and immediately triangulates them (the middle
plot).  The image can be reconstructed by local plane interpolation
(the right plot.) The reconstruction quality can be further improved
by increasing the number of triangles. Figure \ref{fig:opengl} shows
the same image as rendered by the OpenGL graphics library
\cite[]{opengl}.

\inputdir{tri}

\plot{marmousi}{width=6in}{Applying the height triangulation algorithm
  to the Marmousi model. The left plot shows a smoothed and windowed
  version of the Marmousi model. The middle plot is a result of
  10,000-point triangulation, followed by linear interpolation. The
  right plot shows the difference between the two images.}
\par
Figure \ref{fig:marmousi} shows an application of the height
triangulation algorithm to the famous Marmousi model. The left plot
shows a smoothed and windowed version of the Marmousi, plotted on a
501 by 376 computational grid. In the middle plot, 10,000 points from
the original 188,376 were selected for triangulation and interpolated
back to the rectangular grid. The right plot demonstates the small
difference between the two images.

\subsection{Mesh Refinement}

One the main properties of Delaunay triangulation is that, for a
given set of nodes, it provides the maximum smallest angle among
all possible triangulations. It is this property that supports the wide
usage of Delaunay algorithms in the mesh generation problems.
However, it doesn't guarantee that the smallest angle will always be small.
In fact, for some point distributions, it is impossible to avoid
skinny small-angle triangles. The remedy is to add additional
nodes to the triangulation so that the quality of the triangles is
globally improved. This problem has become known as
\emph{mesh refinement} \cite[]{ruppert}.

\plot{hole}{width=6in,height=1.5in}{An illustration of Rivara's
  refinement algorithm. The left plot shows an input to the algorithm:
  a valid Delaunay triangulation with ``skinny'' triangles. Three
  other plots show successive applications of the refinement
  algorithm.}
\par
One of the recently proposed mesh refinement algorithms is Rivara's
\emph{backward longest-side refinement} technique \cite[]{rivara}.  The
main idea of the algorithm is to trace the LSPP (longest-side
propagation path) for each refined triangle. The LSPP is an ordered
list of triangles, connected by a common edge, such that the longest
triangle edge is strictly increasing. After tracing the LSPP, we
bisect the longest edge and insert its midpoint into the
triangulation. Rivara's algorithm is remarkably efficient and easy to
implement. In comparison with alternative methods, it has the
additional advantage of being applicable in three dimensions.
\par
Figure \ref{fig:cerveny} demonstrates an application of different
triangulation techniques to a simple layered model, borrowed from the
Seismic Unix demos (where it is attributed it to V.\v{C}erven\'{y}.)
Another model from \cite{hale} is used in Figure \ref{fig:susalt}.

\plot{cerveny}{width=6in,height=4in}{A comparison of different
  triangulation techniques on a simple layered model. The top left
  plot shows the original model; the top right plot, the result of
  noncomforming triangulation; the two bottom plots, conforming
  triangulation and an additional mesh refinement.}

\plot{susalt}{width=6in,height=4in}{A comparison of different
  triangulation techniques on a simple salt-type model. The top left
  plot shows the original model; the top right plot, the result of
  noncomforming triangulation; the two bottom plots, conforming
  triangulation and an additional two-step mesh refinement.}

\subsection{Implementation Details}
Edge operations form the basis of the incremental algorithm.
Therefore, it is convenient to describe triangulation with
edge-oriented data structures \cite[]{stolfi}. I have followed the idea
of \cite{hansen} of associating with each edge two pointers to the
end points and two pointers to the adjacent triangles.  The triangle
structure is defined by three pointers to the edges of a triangle.
Additionally, as required by the point location algorithm, each
triangle has a pointer to its ``children.'' This pointer is NULL when
the triangle belongs to the current Delaunay triangulation.
\par
Many researchers have pointed out that the geometric
primitives used in triangulation must be robust with respect to
round-off errors of finite-precision calculation. To assure the
robustness of the code, I used the adaptive-precision predicates of
\cite{shewchuk}, available as a separate package from the
\texttt{netlib} public-domain archive.

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: