Simplification of the algorithm using the DtRef method

The Bowyer–Watson algorithm can be simplified by leveraging the topological model of the triangulation and its navigation functions, such as traversing the triangles sharing a common vertex.

To simplify the algorithm, we replace the removal of triangles and the triangulation of the hole with the transformation of non-conforming triangles around the insertion point, which leads to the same result.

The modified algorithm only creates triangles, thereby avoiding memory fragmentation issues and calls to the garbage collector.

This method is referred to as DtRef in the remainder of the document.

The algorithm of the DtRef method is :

1/ Create the "super triangle"
2/ Insert the scatter points
    For each point (p):
      - find the triangle (t) containing point (p)
      - split triangle (t) into three at point (p)
      // Keep Delaunay's conformity
      - As long as there is a non-conforming triangle (ts) around vertex (p):
         - flip the triangle (ts) and its relevant adjacent triangle
      - End of loop (ts)
    End of loop (p)
3/ Remove the triangles at the edge of the fictitious envelope
4/ Remove abnormal triangles on the boundary

Insert process including edge cases

Edge cases (see point status) arising from calculation precision must be taken into account, which somewhat complicates the algorithm.

The triangle search process iterates through all triangles and, for each one, determines the point's status to handle edge cases:

The detailed insertion step of the DtRef method becomes:

2/  Insert the scatter points
    For each point (p):
      For each triangle (t)
        - evaluate the state of point (p) for trianle(t):
          - case near vertex:
             - exclude point (p)
             - next point (p)
          - case near segment and near angular sector:
            - set point (p) aside
            - next point (p)
          - case "inside":
            - split triangle (t) into three at point (p)
            - As long as there is a non-conforming triangle (ts) around vertex (p) :
               - flip the triangle (ts) and its relevant adjacent triangle
            - End of loop (ts)
            - next point (p)
      End of loop (t)
    End of loop (p)

"Near segment and near angular sector" points are set aside for processing after the insertion phase.

Processing them may result in:

Sample processing

Illustration of each sample processing step of the DtRef method. The triangles conforming to the Delaunay triangulation are shaded in blue:

Step 1 Create the super triangle	Step 2 Insert the scatter points, conformity is 100%
Step 3 Remove the triangles at the edge of the fictitious envelope	Step 4 Remove abnormal triangles on the boundary

The final step removes boundary triangles whose longest side exceeds a maximum length or whose ratio (height/base length) is below a threshold:

The DtFlipRef method

The "flip principle" is applied to the DtRef method by eliminating the maintenance of Delaunay conformity during each point insertion, replacing it with iterations performed after all points have been inserted.

The method described below is referred to as DtFlipRef in the remainder of the document. It is a simplified version of the DtRef method with the Flip principle.

Insert process including edge cases

The detailed insertion Flip insert step (1.2) with edge cases is:

1.2 Insert the scatter points
    For each point (p):
      For each triangle (t)
        - evaluate the state of point (p) for trianle(t):
          - case near vertex:
             - exclude point (p)
             - next point (p)
          - case near segment and near angular sector:
            - set point (p) aside
            - next point (p)
          - case "inside":
            - split triangle (t) into three at point (p)
            - next point (p)
      End of loop (t)
    End of loop (p)

"Near segment and near angular sector" points are set aside for processing after the insertion phase.

Flipping triangles

The Flipping step (2) is performed through successive iterations over the set of triangles while processing the points that were set aside:

2/ Flip adjacent triangles that do not meet the Delaunay criterion
   Do:
      2.1/ Process non-compliant triangles:
           For each triangle (t) in the triangulation:
           - If non-conforming flip the triangle (t) and its relevant adjacent triangle
           End of loop (t)
      2.2/ Process the points set aside
   Continue until 100% conformity is reached

The status of the set-aside points is redetermined when the triangles are transformed to reduce them to two cases:

The detailed process of the points set aside is:

2.2/ Process the points set aside
     For each point(p) set aside:
      For each triangle (t):
      - evaluate the status of point (p)
      - "inside" case: split triangle (t)
      - "near vertex" case: exclude point (p)
      - "near segment" case: insert point(p) on segment
      - "near angular sector" case: set point (p) aside
      End of loop (t)
     End of loop (p)

Sample processing

Illustration of each sample processing step of the DtFlipRef method.

Triangles that do not conform to the Delaunay triangulation are shaded in green:

Step 1.1 Create the super triangle	Step 1.2 Insert the scatter points, conformity is 40%
Step 2 Flip adjacent triangles that do not meet the Delaunay criterion	Step 3 Remove the triangles at the edge of the fictitious envelope
Step 4 Remove abnormal triangles on the boundary

DtRef and DtFlipRef metrics for the test dataset

Here is the table of metrics produced by the DtRef and DtFlipRef methods using the test set (981 points).

Metrics N. v and N. s are empty, indicating good dataset quality.

Similarly, the metric Split indicates the absence of collinear points.

Therefore, this dataset does not allow for testing edge cases.

Comparative analysis of the DtRef and DtFlipRef methods

Both methods yield similar metrics:

In fact, these methods area heavily dependent on the distribution of the points. The execution time for this phase can vary significantly depending on the dataset.

These methods are useful for creating and validating the implementation of the fundamental triangulation model (relationships between triangles, point status, split segment, split triangle , shift methods, etc.).

Triangle subdivision

Sorting orders the points and subdivides the triangles by depth, as shown in the illustrations below for the two methods:

DtRefTree division of triangles after one iteration	DtFlipRefTree division of triangles after one iteration
DtRefTree division of triangles after two iterations	DtFlipRefTree division of triangles after two iterations

Triangles that do not meet the Delaunay criterion are light green.

Parallel processing

These two methods pave the way for parallelizing the triangulation calculation.

For instance, one can limit the depth of triangle subdivision:

Merging the triangles requires ensuring the alignment of adjacent triangles along their boundaries.

You can find publications on the "divide and conquer" principle for computing a Delaunay triangulation, such as this one:

DtRefTree and DtFlipRefTree metrics for the test dataset

Updating the metrics table with the two new methods DtRefTree and DtFlipRefTree:

Comparative analysis of the DtRefTree and DtFlipRefTree methods

The point insertion steps of DtRefTree and DtFlipRefTree becomes independent of their distribution.

Sorting the points yields an immediate gain:

Furthermore, the same difference as before is observed between the two methods, just as with the DtRef and DtFlipRef methods:

Convex hull precision

Triangle subdivision with convex hull

The initial triangulation created with the convex hull comprises 14 triangles.

While the previous method starts with only one triangle, this one uses a finer initial triangulation with 14 triangles and conforms more closely to the perimeter of the point set.

The subdivision of the triangles is more precise, as illustrated by the figures below:

Triangulation calculated on the convex hull	Delaunay conformal triangulation
Division of triangles after one iteration	Division of triangles after two iterations

DtCvxTree and DtFlipCvxTree metrics for the test dataset

The convex hull is calculated in a process (Cvx) separate from the triangulation calculation.

The "Σ ms" column displays the total duration of the method.

Updating the metrics table with the two new methods DtCvxTree and DtFlipCvxTree:

Comparative analysis of the DtCvxTree and DtFlipCvxTree methods

Observations on metrics:

Calculation of a convex hull

The calculation is performed at a specific resolution, and the algorithm eliminates points that coincide with the selected vertices.

Here is a simple and efficient algorithm for calculating the convex hull of a set of points that uses the "divide and conquer" principle:

 Initialize the convex hull with the extreme points (x_min, y_max, x_max, y_min).
 Distribute each point located outside the convex hull onto its corresponding projected segment.
 For each non-empty segment (s):
 - Determine the highest point
 - Split the segment (s) at that point
 - Redistribute the points of (s) across the two segments
 End loop (s)

The complexity of the algorithm is O(n log n).

The convex hull of the extreme points and its exterior points

Convex after one iteration

The CvxSurface structure

To achieve this, the convex hull calculation is used to accumulate within a structure (named CvxSurface) all the hulls of the points, proceeding from the outside toward the center; each hull is determined using the points located inside the preceding one. The final contour may consist of three or more points, or it may be flattened into two points or even a single point.

The calculation method for the CvxSurface structure is as follows:

Calculate the convex hull (cvx1) of the set of points
While (cvx1) contains points
 Calculate the convex hull (cvx2) of the points inside (cvx1)
 Continue with (cvx2), which becomes (cvx1)

The CvxRing structure

A second structure (named CvxRing) is derived from the CvxSurface and consists of adjacent concentric rings delimited by two successive convex hulls (an inner one and an outer one), plus a final solid ring with no inner contour.

The structure CvxSurface classifies the 981 points into 50 contours.

The CvxRing structure is composed of 49 rings and a solid last one.

The rings (represented here by alternating shades of green and light blue) cover the space and pass through every point without intersecting.

Triangulation cursor

The area of ​​each ring, bounded by its two boundaries, can be triangulated using triangles adjacent to either the outer (green) boundary or the inner (blue) boundary (Triangulation of rings). To triangulate a ring, a cursor is used on three successive points of its outer (Ti) and inner (Vi) contours. The cursor defines a quadrilateral [Ti, Ti+1, Vi+1, Vi] which, as far as possible, must exhibit an overlap between its outer segment [Ti, Ti+1] and its inner segment [Vi, Vi+1]. The first two vertices define either the triangle adjacent to the outer contour [Ti,Ti+1,Vi] or the triangle adjacent to the inner contour [Vi+1,Vi,Ti] (both being consistently oriented clockwise).

Optimal triangle

As for the third vertices Ti+2 and Vi+2 they serve to select the optimal triangle by measuring the overlap and the elongation length of the two subsequent quadrilaterals. The cursor starts at the vertex (T1) with the minimum x-coordinate on the outer contour and the vertex (V1) closest to the inner contour At each step, the cursor determines which triangle to select outer [Ti,Ti+1,Vi] or inner [Vi+1,Vi,Ti] and advances by one vertex along the contour adjacent to the triangle

Choice of triangle

The choice of triangle is made based on the analysis of the vertices (4 cases):

The inner triangle [Vi+1,Vi,Ti] is invalid (negative area). Select the outer triangle [Ti,Ti+1,Vi] The vertex Vi+1 is included in the outer triangle [Ti,Ti+1,Vi]. Select the inner triangle [Vi+1,Vi,Ti]. Only one of the two quadrilaterals [Ti,Ti+1,Vi+2,Vi+1] and [Ti+1,Ti+2,Vi+1,Vi], or neither of them, may exhibit an overlap between the outer and inner segments. Select the triangle preceding the quadrilateral with the overlap. Finally, choose the triangle that minimizes the quadrilateral elongation length.

Triangulation of the Cvx Ring structure

The CvxRing structure thus allows the triangulation to be partitioned ring by ring:

The triangulation of each ring primarily produces flat triangles, however:

The transformation of the triangles is carried out in five steps:

4/ Transform the triangles on each ring
   For each ring (r) other than the last one:
   - Flip triangles that do not meet the base-max criterion on the inner boundary of the ring (r)
   End of loop (r)
   For each ring (r) other than the last one:
   - Flip triangles that do not meet the Delaunay criterion on the inner boundary of the ring (r)
   End of loop (r)
   For each ring (r):
   - Flip triangles that do not meet the Delaunay criterion on the outer boundary of the ring (r)
   End of loop (r)
   - Flip triangles that do not meet the Delaunay criterion on the last ring

Following the calculation of the CvxSurface structure, the calculation method for the CvxRing structure is as follows:

Compute the (CvxRing) structure from the CvxSurface structure:
1/ Create the rings for each contour
2/ Triangulate each ring
3/ Establish relationships between the triangles at the edges of each ring
4/ Transform the triangles on each ring

Several steps can be parallelized:

Triangulate each ring

Transform triangles on each ring, conformity is 100%

The algorithm with the CvxRing structure

Once the CvxRing structure has been calculated, the triangulation can be formed by accumulating the triangles from the rings.

The new algorithm can be summarized as follows:

0/ Calculate the two structures CvxSurface and CvxRing (which can be processed during a separate upstream step)
1/ Create the triangulation by accumulating the rings (after this step, the triangulation is Delaunay compliant)
2/ Remove abnormal triangles on the boundary

Edge cases

Edge cases are naturally handled by the method that does not require calculating the point status:

DtCvxRing metrics for the test dataset

Updating the metrics table with the new method DtCvxRing with calculations of the CvxSurface and CvxRing structures prior to the triangulation calculation:

Analysis of the DtCvxRing method

Observations on metrics:

In conclusion:

The second dataset

The second dataset is a topographic survey with 569 points.

The third dataset

The third dataset contains 3,839 points (1 of which are excluded due to their proximity to another point).

The fourth dataset

The fourth dataset contains 9,811 points (1 of which are excluded due to their proximity to another point).

Method for evaluating correlation equations for the number of operations

By treating the measurements (Σ op, np) obtained for each method as a series, probabilistic techniques (calculation of standard deviation, variance, and variation) are applied to evaluate the correlation equations.

Subsequently, to evaluate a correlation equation Σ op=f(np), its percentage variation (standard deviation divided by the mean value) is calculated across the five datasets.

A percentage variation of less than 15% is accepted as an indicator of the equation's plausibility.

Evaluation of the complexity of the DtRef method

For the DtRef method, as expected given its O(n²) complexity, we observe that the value 2*(np)² is close to the value Σ op, as shown by the calculation for the last dataset:

2*(np)² = 2*(3,839)² = 29475,842  op = 29291,302.

To confirm this equation for this method and DtFlipRef, we calculate its percentage variation:

The value of the variation is less than 15%, confirming the equation.

To proceed, we represent this equation on a graph with for each method:

It can be seen from the graph that:

It is assumed that:

Complexity evaluation of the DtRefTree, DtFlipRefTree, DtCvxTree and DtFlipCvxTree methods

We generalize the DtRef equation to determine the equations for the other four by calculating:

Its selection of the best correlation for each method leads to the following equations :

Following this calculation, we observe:

Next, the (Σ op / np^√2) ratio of the equation is determined for the three methods, followed by the variation:

Complexity evaluation of the DtCvxRing method

For the DtCvxRing method, and in order to confirm the assumed O(n log n) complexity, we evaluate the following ratio R(p):

A bisection search on the interval [1, 100] for (p) :

Summary of the method complexity assessment

The table below summarizes the correlation equation (Σ op ≈ f(np)), the percentage of variation, and the complexity for each method:

To complement this table, the following graph shows:

A scale factor is applied to each axis to fit the 4:3 ratio.

Quantified evaluation criteria

We use the measurements obtained on the final dataset to compare the last five methods.

Each criterion results in a score from 1 to 10.

Evaluation of methods using quantified criteria

Criteria DtRefTree DtFlip.. DtCvxTree DtFlip.. DtCvxRing Memory 1 5 6 10 6 Calculations 1 6 1 5 10 Shifts count 10 2 10 1 6 Σ 12 13 17 16 22

Design criteria

These criteria do not result in a measured score but identify a characteristic classified as a strength.

Evaluation of methods using design criteria

(1) Parallelizing depth-first sorting methods is possible, provided the algorithm undergoes significant modifications.