A Brute-Force Constrained Delaunay Triangulation?

Question

I need to create a triangle mesh from a set of points. The set has very few points so it doesn't need to be fast or optimised (I will deal with 100 points maximum). The mesh needs to be a constrained "delaunay triangulation". In the image below I showed (on the left) the set of points I start from (blue and red dots). I also know the connections between these points (the outline in black). The mesh needs to look like the example on the right (including the edges in grey that form outside and inner triangles).

I can't use libraries.

I looked at many different algorithms. They are many and it's easy to be confused. I would like to know if there is a naive and thus hopefully simpler algorithm I can use in order to produce the mesh on the right? Brute force approach is fine (ps: I can do a Delaunay triangulation).

Answer 1

Thanks for all the answers.

I went through the process of developing a solution to this problem so thought I would share my own experience, hoping people facing the problem will find the insight useful.

So from my own experience implementing an algorithm I came to the conclusion:

1. There is not really quick way to this problem. It is not reasonable to think it can be achieved in just 50 lines of code. In fact the routine that I wrote (C++) is about 400 to 500 lines (hard to tell with comments). So reasonably compact yet challenging and it took me a 2 to 3 days to get the logic right (it can be tricky).

2. I found the algorithm propose by Sloan in "A FAST ALGORITHM FOR GENERATING CONSTRAINED DELAUNAY TRIANGULATIONS" to be perfectly well suited for the problem at hand. The reality when it comes to Delaunay triangulation which was a new subject for me, is that there seems to be a lot of different algorithms approach and this research is pretty old. So for a new comer it's really hard to know where to start.

   2.1. It's hard to know which algorithm is recent, simple in its comprehension and fast and simple to implement.

   2.2 Generally once you have understood the principle it's mostly a matter of coding the logic in the most efficient way (and it seems that this what most of the algorithms/papers are fighting above).

   2.3 I found the paper from Sloan to be understandable, very well explained. If you follow the logic and the instructions, then anyone can really implement a constrained Delaunay triangulation.

So in conclusion:

1. I recommend the Sloan paper because it includes an explanation on how to create a Delaunay triangulation followed by a constrained triangulation if necessary.

2. To answer my own question there is not really brute force to this problem. Implementing this technique just requires to implement the full logic and most implementation must require more or less the same amount of work

3. With the nuance, that I wasn't looking after much optimisations because my point sets are really small. So I am sure many algorithms are better than the one described by Sloan; they probably propose optimized data structures and algorithms optimized to minimize steps such point insertion in triangulation etc.

So anyway Sloan worked. A small image to illustrate the answer and make it more attractive;-)

EDIT

This is production code so helas I can't share that... I could lead me to be fired. The process is very simple though. You look for the intersection between a segment (your constraint) and all edges in the model. Then for each intersected edge, you swap the diagonal between the 2 triangles that this edges belongs to. If the new diagonal intersects the segment still, then add the new diagonal back onto the stack of intersected edges for this segment. If the new diagonal doesn't intersect the segment then add it to the stack of newly created edges. Keep processing the stack of intersected edges until it's empty.

Then once this is finished you need to process the list of new added edges. For each one of them, check that the Delaunay triangulation criterion is respected. If not swap the diagonal of the triangle this edge belongs to. Simple ...

This is just the paper ...

Point Set

26.9375 10.6875
32.75 9.96875
31.375 4.875
27.6562 2.0625
23.9375 -0.75
18.1562 -0.75
10.875 -0.75
6.60938 3.73438
2.34375 8.21875
2.34375 16.3125
2.34375 24.6875
6.65627 29.3125
10.9688 33.9375
17.8438 33.9375
24.5 33.9375
28.7188 29.4062
32.9375 24.875
32.9375 16.6562
32.9375 16.1562
32.9062 15.1562
8.15625 15.1562
8.46875 9.6875
11.25 6.78125
14.0312 3.875
18.1875 3.875
21.2812 3.875
23.4687 5.5
25.6562 7.125
8.46875 19.7812
27 19.7812
26.625 23.9688
24.875 26.0625
22.1875 29.3125
17.9062 29.3125
14.0312 29.3125
11.3906 26.7188
8.75 24.125

These are x/y/z coordinates (z=0)

Segments:

0 1
1 3
3 5
5 7
7 9
9 11
11 13
13 15
15 17
17 19
19 20
20 22
22 24
24 26
26 0
28 29
29 31
31 33
33 35
35 28

Indices start at 0 (0 -> first vertex in vertex list)

Answer 2

I tried it with alpha shapes with good results for a few shapes https://concavehull.codeplex.com/ but it's nowhere near the original constrained delaunay triangulation.

Here is my alpha-shape algorithm:https://alphashape.codeplex.com.

Answer 3

A simple approach seems to be to implement a ear clipping algorithm. Without optimisation as in hash grids or quad trees. For ear clipping you just check every three consecutive vertices a,b, and c. If b is convex and no other vertex of the polygon lies inside the triangle abc then you can clip this triangle reducing the boundary of the polygon by one vertex, b.

Additionally you have to store the neighbourhood relations. Thus, reference from each triangle its, at most three, neighbours.

When the triangulation is finished you convert it to the constrained Delaunay triangulation (CDT). This can be done by edge flipping. Therefore you have to check for every triangle the circumcircle. If no vertex of a neighbouring triangle lies inside the triangle is conform to the CDT otherwise flip the edge of the triangle where the violation occurs.

Edit due to @Betterdev in the comments blow: Possible holes in the input polygon can be added to the initial boundary by adding a bridge. As a preprocessing one can connect a vertex of a hole to a vertex of the boundary by a "double" edge. This is always possible and makes each hole part of the main polygon boundary; and works well with ear clipping. Storing the neighbour through these bridges is vital to the flipping however.