The triangulation of a set of given two-dimensional points is a well-researched problem and an impressive variety of efficient triangulation methods can be found in the literature [116] [95] [120] [121] [122] [123] [124] [125] [126] [127] [128] [129] [130] [131] [132]. This is only exceeded by the massive scope of applications of the Delaunay triangulation (and the Voronoi tessellation) [115] [95] [133] [134] [135] [136].
A classical approach is the divide-and-conquer method which is described in detail in [95]. A direct practical realization of the algorithm faces two severe problems that had to be overcome.
The first disadvantages of the initial method are the limiting assumptions that no more than three Delaunay nodes lie on a circle (are cocircular) and that no more than two Delaunay nodes are collinear. In these degenerate cases the behavior of the classical algorithm proposed in [95] is unpredictable. It is clear that these degeneracies will occur when a tensor product grid or any other grid that exhibits a partly regular structure is triangulated. So the degeneracies are inevitable and the classical divide-and-conquer algorithm must be modified and extended. Appropriate numerical tolerances must be used to detect the special point constellations.
The second necessary improvement with respect to the original divide-and-conquer algorithm is the use of the quadtree (which was originally needed for the boundary refinement and grid point relaxation steps) for the recursive decomposition of the point set in disjoint subsets. The classical method uses a purely vertical recursive decomposition of the point cloud. This limits the application of the original algorithm to small point sets because due to the repeated vertical partitioning, the point sets approach vertical lines which can in the limit not be distinguished from degenerate cases. This leads to both numerical instabilities and an excessive overhead in ``garbage'' Voronoi edges.