Constrained Delaunay Triangulation

Let be the set of at least distinct points in -dimensional space () from which at least points are not in the same -dimensional hyperplane and let be the set of -dimensional constraining simplices which are mutually not intersecting and which are not crossing any other point of except their vertices that are the points in .

The -dimensional general triangulation of is a partition of the convex hull of into -dimensional mutually not intersecting simplices of nonzero volume whose vertices are the points in . There may be many triangulations of the same set .

The constrained -dimensional triangulation is such a triangulation which contains all constraining simplices from .

The -dimensional Delaunay triangulation of is such a triangulation , if for each of its simplices there is no point from inside -dimensional sphere (defined as the set of points of the same distance from a given point called center) circumscribed to vertices of that simplex. This is known as empty-circumsphere property or Delaunay property. If no points in are co-spherical without any other point inside that sphere, then the Delaunay triangulation is unique (non-degenerated), otherwise there exist several degenerated Delaunay triangulations of .

The constrained -dimensional Delaunay triangulation of is such a constrained triangulation , if for each of its simplices there is no point from inside -dimensional sphere circumscribed to vertices of that simplex , unless the point is on the opposite side of constraining simplex from than or unless the point is in the plane of bounding and simultaneously inside the -dimensional sphere circumscribed to in that plane.

In spaces of , there may exist no constrained (Delaunay or non-Delaunay) triangulation for given sets and . The example of such a case is the Schönhardt polyhedron.

For a 2D set of points, the examples of general triangulation , constrained triangulation , Delaunay triangulation , and constrained Delaunay triangulation are depicted in Figure 1.

The Delaunay triangulation is appealing due to its interesting properties and especially due to the existence of simple and efficient methods [37,38] for its construction. One of the most favourite ones is the Watson's algorithm. It is an incremental algorithm which starts from an initial Delaunay triangulation (formed by a few or by just a single simplex) completely surrounding the convex hull of . Each point of is then inserted (one at a time) into the current Delaunay triangulation to form a new one as a union of two sets of simplices. In the first set, there are all simplices from whose empty-sphere property is not violated by . The second set is composed by new simplices created in the star-shaped cavity defined by all simplices from whose empty-circumsphere property is violated by . Each newly created simplex is formed by and by the -dimensional simplex bounding the cavity.

An important aspect of assembling the Delaunay triangulation is the handling of degeneracy causing the non-uniqueness of the resulting triangulation. Due to the incremental nature of the Watson's algorithm, the degeneracy occurs even if the set is not degenerated. The remedy to this problem consists in changing the order in which the points are inserted, but this is useless if is degenerated. The other way around is to perturb the location of co-spherical points to make them not co-spherical. In reality, the co-spherical points can be reliably detected only if computationally demanding precise floating point arithmetic is adopted. When standard floating point arithmetic is used, even more caution is required, because the co-spherical predicate is not reliable if the points are not only exactly but also almost co-spherical.

On the other hand, the construction of the constrained Delaunay triangulation is not straightforward. In spaces of dimension 3 or more, not every polyhedron can be discretized into constrained Delaunay triangulation without using additional points. Moreover, it was shown [39] that it is an NP-hard problem to identify such a case. Therefore the most common discretization scenarios follow the way that the Delaunay triangulation is built first and then it is subjected to topological and geometrical modifications to enforce the appearance of in the resulting triangulation. This process is known as boundary recovery. In 2D, the boundary recovery (see Figure 2) can be accomplished purely by topological transformations (diagonal swapping). In higher dimensional spaces, however, the recovery algorithm has usually to resort to insertion of additional points to refine the set . This approach has two negative aspects. Firstly, the description of the real geometry of the boundary of the computational domain may not be available, which can result in the undesirable distortion of the geometry, and secondly, the resulting triangulation is not conforming to the original set . This complicates the application of this approach to the separate discretization of individual adjacent regions composing the whole computational domain (for example in parallel processing), unless the additionally inserted points are again removed.

In the present approach, the constrained Delaunay triangulation is constructed from the beginning simultaneously with inserting points from . Once a constraining simplex from appears in the intermediate triangulation, it is not allowed to be removed (during the cavity formation). This requires only a simple change to the Watson's incremental insertion algorithm. Should the cavity include a simplex bounded by it is accepted to form the cavity only if it is on the same side of as the point being inserted. To make this criterion unambiguous, the inserted point must not be located in plane of and simultaneously inside or on -dimensional sphere circumscribed to in that plane. Thus the goal reduces just to ensure that each simplex in appears in the intermediate constrained Delaunay triangulation immediately after all its vertices have been inserted. This is achieved by the a priori boundary recovery performed before the actual point insertion algorithm starts.