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$.\parIf, 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 ofintersection into the triangulation.  The rest of the process is verymuch like the process of edge validation in the original incrementalalgorithm.\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 atriangulation in this case is to select a smaller subset of theinitial grid points, using them as the input to a triangulationprogram. We need to select the points in a way that preserves the mainfeatures of the original image, while removing some unnecessaryredundancy 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 approachesto this problem and proposed a fast version of the incremental\emph{greedy insertion} algorithm. Their algorithm adds pointsincrementally, selecting at each step the point with the maximuminterpolation error with respect to the current triangulation. Thougha straightforward implementation of this idea would lead to anunacceptably slow algorithm, Garland and Heckbert have discoveredseveral sources for speeding it up. First, we can take advantage ofthe fact that only a small area of the current triangulation getsupdated at each step. Therefore, it is sufficient to recompute theinterpolation error only inside this area. Second, the maximumextraction procedure can be implemented very efficiently with apriority 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.}\parFigure \ref{fig:sphere} illustrates this algorithm with a simpleexample. The original image (the left plot) contained 10,000 points,laid out on a regular rectangular grid. The algorithm selects asmaller number of points and immediately triangulates them (the middleplot).  The image can be reconstructed by local plane interpolation(the right plot.) The reconstruction quality can be further improvedby increasing the number of triangles. Figure \ref{fig:opengl} showsthe 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.}\parFigure \ref{fig:marmousi} shows an application of the heighttriangulation algorithm to the famous Marmousi model. The left plotshows a smoothed and windowed version of the Marmousi, plotted on a501 by 376 computational grid. In the middle plot, 10,000 points fromthe original 188,376 were selected for triangulation and interpolatedback to the rectangular grid. The right plot demonstates the smalldifference between the two images.\subsection{Mesh Refinement}One the main properties of Delaunay triangulation is that, for agiven set of nodes, it provides the maximum smallest angle amongall possible triangulations. It is this property that supports the wideusage 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 avoidskinny small-angle triangles. The remedy is to add additionalnodes to the triangulation so that the quality of the triangles isglobally 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.}\parOne of the recently proposed mesh refinement algorithms is Rivara's\emph{backward longest-side refinement} technique \cite[]{rivara}.  Themain idea of the algorithm is to trace the LSPP (longest-sidepropagation path) for each refined triangle. The LSPP is an orderedlist of triangles, connected by a common edge, such that the longesttriangle edge is strictly increasing. After tracing the LSPP, webisect the longest edge and insert its midpoint into thetriangulation. Rivara's algorithm is remarkably efficient and easy toimplement. In comparison with alternative methods, it has theadditional advantage of being applicable in three dimensions.\parFigure \ref{fig:cerveny} demonstrates an application of differenttriangulation techniques to a simple layered model, borrowed from theSeismic 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 withedge-oriented data structures \cite[]{stolfi}. I have followed the ideaof \cite{hansen} of associating with each edge two pointers to theend points and two pointers to the adjacent triangles.  The trianglestructure is defined by three pointers to the edges of a triangle.Additionally, as required by the point location algorithm, eachtriangle has a pointer to its ``children.'' This pointer is NULL whenthe triangle belongs to the current Delaunay triangulation.\parMany researchers have pointed out that the geometricprimitives used in triangulation must be robust with respect toround-off errors of finite-precision calculation. To assure therobustness 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: