\append{Incremental DELAUNAY TRIANGULATION and related problems}Delaunay triangulation \cite[]{delaunay2,sibson,stolfi} is a fundamentalgeometric construction, which has numerous applications in differentcomputational problems. For a given set of nodes (points on theplane), Delaunay triangulation constructs a triangle tessellation ofthe plane with the initial nodes as vertices. Among all possibletriangulations, the Delaunay triangulation possesses optimalproperties, which make it very attractive for practical applications,such as computational mesh generation. One of the most well-knownproperties is maximizing the minimum triangulation angle. In threedimensions, Delaunay triangulation generalizes naturally to atetrahedron tessellation.\parSeveral optimal-time algorithms of Delaunay triangulation (and itscounterpart--Voronoi diagram) have been proposed in the literature.The divide-and-conquer algorithm \cite[]{shamos,stolfi} and thesweep-line algorithm \cite[]{fortune} both achieve the optimal $O (N\log N)$ worst-case time complexity. Alternatively, a family ofincremental algorithms has been used in practice because of theirsimplicity 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}.\parThe 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}\parThe Delaunay criterion can be reduced in the second step to a simple\emph{InCircle} test \cite[]{stolfi}: if a circumcircle of a trianglecontains another triangulation vertex in its circumcenter, then theedge between those two triangles should be ``flipped'' so that two newtriangles are produced. The testing is done in a recursive fashionconsistent with the incremental nature of the algorithm. When a newnode 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 oftest failure, a pair of triangles is replaced by the flip operationwith another pair, producing two more edges to test. Under therandomization assumption, the expected total time of point insertionis $O (N)$.  Randomization can be considered as an external part ofthe 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: alltriangles, occurring in the incremental triangulation process, arekept in memory, associated with their ``parents.'' One or two pointlocation tests (\emph{CCW} tests) are sufficient to move to a lowerlevel of the tree. The search terminates with a current Delaunaytriangle.\parTo test the algorithmic performance of the incremental construction, Ihave profiled the execution time of my incremental triangulationprogram with the Unix \texttt{pixie} utility. The profiling result,shown in Figures~\ref{fig:itime} and~\ref{fig:ctime}, compliesremarkably with the theory: $O (N \log N)$ operations for the pointlocation 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 CPUtime, 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.}\parA straightforward implementation of Delaunay triangulation wouldprovide an optimal triangulation for any given set of nodes. However,the quality of the result for unfortunate geometricaldistributions of the nodes can be unsatisfactory. In the rest of thisappendix, I describe three problems, aimed at improving thetriangulation quality: conforming triangulation, triangulation ofheight fields, and mesh refinement.  Each of these problems can besolved with a variation of the incremental algorithm.\subsection{Conforming Triangulation}\inputdir{tri}In the practice of mesh generation, the input nodes are oftensupplemented by boundary edges: geologic interfaces, seismic rays, andso on. It is often desirable to preserve the edges so that they appearas edges of the triangulation \cite[]{SEG-1994-0502}. One possible approach is\emph{constrained} triangulation, which preserves the edges, but onlyapproximately satisfies the Delaunay criterion \cite[]{lee,chew}. Analternative, less investigated, approach is \emph{conforming}triangulation, which preserves the ``Delaunayhood'' of thetriangulation by adding additional nodes \cite[]{hansen} (Figure\ref{fig:conform}).  Conforming Delaunay triangulations are difficultto analyze because of the variable number of additional nodes. Thisproblem was attacked by \cite{edels}, who suggested an algorithmwith a defined upper bound on added points. Unfortunately,Edelsbrunner's algorithm is slow in practice because the number ofadded points is largely overestimated.  I chose to implement amodification of the simple incremental algorithm of Hansen and Levin.Although Hansen's algorithm has only a heuristic justification andsets no upper bound on the number of inserted nodes, its simplicity isattractive for practical implementations, where it can be easilylinked with the incremental algorithm of Delaunay triangulation.\parThe incremental solution to the problem of conforming triangulationcan 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.}\parTo insert an edge $AB$ into the current triangulation, I use thefollowing 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}\parThe intersection point  of edges $AB$ and $EF$ is given by the formula\begin{equation}