Some thoughts on constrained Delaunay triangulations

[DanielVandH]

Need to implement constrained Delaunay triangulation at some point. The method here https://www.sciencedirect.com/science/article/pii/S0925772115000322?via%3Dihub would be good. Once this is implemented, de Berg's method will probably be removed.

Some thoughts:

• The interface for this will be important. It is important to be able to support:
  • General vertices.
  • Specified edges. These are just edges that will appear in the triangulation, and could be e.g. a matrix where each column defines an edge. For example, E = [1 2 3; 4 5 6] defines the edges [1, 4], [2, 5], and [3, 6] that have to appear in the triangulation. Could also allow E = [[1, 4], [2, 5], [3, 6]] for the same thing.
  • Specified holes. Holes will be a list of vertices v such that v[begin] == v[end], listed in clockwise order. Holes will be special as all edges that are inside the hole are removed. How this is done is to be decided, but the adjacent map should make this relatively simple. Note that the clockwise order is important to be consistent with all triangles being counterclockwise.
  • A specification for the boundary. Need to somehow handle a boundary loop as being different to a hole, and to allow the boundary to simply be the convex hull. Will need to think about this a bit more. Maybe we simply have a rule that boundaries are given counterclockwise and holes clockwise, so that if no counterclockwise boundary is given then we default to the convex hull. This could be extended to allow for the triangulation of multiply-connected domains, where there are multiple exterior boundaries.

With this interface, we could define three methods:

triangulate(pts) # Unconstrained Delaunay
triangulate(pts, edges) # Constrained Delaunay, no holes (this is a planar straight line graph)
triangulate(pts, edges, loops) # Constrained Delaunay, with holes (this is a piecewise linear complex)

(Maybe a PSLG or PLC struct needs to be defined?)

The returned value for all of these could still be a Triangulation. Triangulation would be updated to be something like

struct Triangulation
    points
    adjacent
    adjacent2vertex
    graph
    edges # better named needed probably - this will just be the list of constrained edges, not all edges in the triangulation - edges(graph) serves that purpose
    loops
    boundary_nodes
end

(The boundary_nodes field is just something extra to store boundaries for both the exterior and for the interiors.)

• With these new features, assigning a value to edges that have no associated triangle is more complicated. e.g. if (i, j) is a boundary edge, currently adj(i, j) = 0 so that (i, j) is a boundary edge and (i, j, 0) is a ghost triangle. When we have interior holes, this is more delicate. Perhaps if we have a set of loops loops = [v1, v2, v3] (say), then adj(i, j) = -1 if the edge is a part of v1, -2 if it is a part of v2, and -3 if it is a part of v3, and so on. (Note that this drops support for non-one-based-indexing.) This is probably the simplest choice, though we would need to adjust DefaultAdjacentValue to something other than -4 (de Berg being removed will free the indices -1, -2, -3.

• Point location will also be more complicated. This is where I struggled the most with the development of the unconstrained code (so many corner cases...), so this will most likely require the most thought again.

• My Gmsh code could be extended as a good way to test all this, providing the interfaces I give above for triangulate except for generate_mesh (except entirely in terms of loops, since the main point there is dealing with boundaries and holes).

Will probably have to sit on these ideas for a while before getting the time to complete it, and to just think about the interface in general. The point location work needs a lot of thought as mentioned.

[DanielVandH]

My Gmsh code could be extended as a good way to test all this, providing the interfaces I give above for triangulate except for generate_mesh (except entirely in terms of loops, since the main point there is dealing with boundaries and holes).

This has now been done (not on this master yet, need to incorporate the interface). Just provide a vector of vector of vectors to provide a list of boundaries, with boundaries split by segments, and the first boundary is the outer-most boundary. Very simple to do. Will be a good test case for point location, etc.

Does give me the idea to remove loops from the above, and just let boundary_nodes be such that boundary_nodes[m][n] is the list of boundary nodes for the nth segment of the mth boundary. The outer-most boundary is boundary_nodes[m]. If a vector of floats was provided, boundary_nodes is just a vector (rather than a vector of vector of vectors).

Example mesh:

[DanielVandH]

Some progress for this is in https://github.com/DanielVandH/DelaunayTriangulationDev.jl, where I've also rewritten how the predicates are implemented (with each predicate returning a Certificate type defined via EnumX.jl). No constrained triangulation yet, but everything has been rewritten to remove some assumptions that will no longer be true with multiple boundaries (e.g. BoundaryIndex is the only boundary index). Will be ported over eventually, probably once the interface is fully finalised + constrained triangulation is added, aiming for a 1.0 release.

[DanielVandH]

See #32 for some progress

[DanielVandH]

Still thinking about this. https://github.com/DanielVandH/DelaunayTriangulation.jl/tree/convex-triangulation has some work on triangulation convex polygons, i.e. vertex deletion, which involves code that will later go into constrained triangulation. Just have to keep working on it so that I understand how it works in the presence of cocircular/collinear points.