How to choose the points to triangulate a domain

Question

Choosing Mathematica to expose the problem, defined the points:

points = {{1, 0}, {2, 2.5}, {3, 0}, {2.5, 3.5}, {4, 5}, {2.5, 4.5}, 
          {2, 7}, {1.5, 4.5}, {0, 5}, {1.5, 3.5}, {1, 0}};

it's possible to define the following boundary:

boundary = BoundaryMeshRegion[points, Line[Range[Length[points]]]];

and finally it's possible to obtain the following triangulation:

TriangulateMesh[boundary, Method -> "ConstrainedQuality", MaxCellMeasure -> 0.1]

Mesh Plot

where, if I am not mistaken, this triangulation:

• must have at least the imposed ones as nodes;

• must satisfy the Delaunay condition;

• must satisfy the condition on the area (at most equal to 0.1).

If, on the one hand, I think I understand the constraints to which the triangles of the mesh must undergo, what isn't absolutely clear to me is how to generate the additional points of the first attempt, which if I haven't misunderstood are called Steiner points.

If we were dealing only with star domains, an idea would be to distribute these points along the sides of smaller polygons obtained through homothety with respect to the centroid of the known points of the boundary, while I would not know how to obtain an equally general rule by also incorporating non-star domains.

So, the question arises: how would you solve this problem? How is it customary to generate such internal points and then proceed with their connection obtaining a triangulation?

Thank you!

Answer 1

One approach is outlined in Pfeifle and Seidel, Triangular BSplines for Blending and Filling of Polygonal Holes, Section 3.

My summary may not be watertight or 100% follow the paper, but roughly speaking:

1. select a maximum edge length $\epsilon$.
2. Refine domain edges so that no edge is longer than $\epsilon$
3. Triangulate the domain
4. Wherever a triangle's edges are too long, add it's centroid as an interior point. When all new ints have been added, flip edges to achieve a Delaunay triangulation. Repeat step 4 until no edge is longer than $\epsilon$.

You can also do a relaxation at each step that uniformizes the distribution of points. E.g. each point is moved towards the average position of its neighbours.

There are variations and improvements possible, and it's not hard to come up with pathological cases, but this is a good starting point. It may be desirable to let $\epsilon$ vary to allow the density of the interior and boundary points to vary depending on the complexity of the nearby boundary.