Determine which delaunayTriangulation Constraints correspond to input constraints?

Here's a minimal example:

tri = delaunayTriangulation([0 0;1 0;0 1;1 1],[1 4;2 3]);

Warning: Intersecting edge constraints have been split, this may have added new points into the triangulation.

triplot(tri)

adds an extra point at the intersected constraints at (0.5,0.5). Correspondingly, tri.Constraints splits the input 2 constraints into 4 constraints.

Is there a way to retrieve which output constraints are mapped to which input constraints? I'd like to rely on the book-keeping of the algorithm rather than try to recompute this myself geometrically (as the numerics can get tricky in harder examples).

Tools like https://www.cs.cmu.edu/~quake/triangle.html and https://wias-berlin.de/software/tetgen/ allow tagging the input vertices and constraints with markers. Is there a similar feature?

[Matt J]

Is there a way to retrieve which output constraints are mapped to which input constraints?

Don't you mean the other way around? Which input constraints are mapped to which output constraints (inputs map to outputs)?

[Alec Jacobson]

Yeah yeah

Answer 1

[Matt J]

If you convert the final triangulation to a graph, then it seems to me that shortestpath should give you the intermediate nodes on the straight line between a given pair of input edge nodes. The intermediate steps correspond to the output constraints.

[Matt J]

Downloading spatialgraph2D can help with the implementation:

V=[0 0;1 0;0 1;1 1];
C=[1 4;2 3];

tri = delaunayTriangulation(V,C);

Warning: Intersecting edge constraints have been split, this may have added new points into the triangulation.

s=tri.ConnectivityList(:,1:3);
t=tri.ConnectivityList(:,[2,3,1]);

G=simplify(graph(s,t));
G=spatialgraph2D(G,tri.Points(:,1),tri.Points(:,2)).G;


for i=1:height(C)
    
  E=shortestpath(G,C(i,1), C(i,2))      %output edge constraints

end

E = 1×3

1 5 4

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

E = 1×3

2 5 3

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Answer 2

[Steven Lord]

I have not looked at the spatialgraph2D tool Matt J suggested, but I suspect you can do the same thing using functionality included in MATLAB already.

V=[0 0;1 0;0 1;1 1];
C=[1 4;2 3];

tri = delaunayTriangulation(V,C);

Warning: Intersecting edge constraints have been split, this may have added new points into the triangulation.

s=tri.ConnectivityList(:,1:3);
t=tri.ConnectivityList(:,[2,3,1]);

The coordinate data is in the Points property of the delaunayTriangulation. Let's turn it into a table with two variables, one the X coordinates and one the Y coordinates.

coordTable = array2table(tri.Points, ...
    VariableNames = ["x", "y"]);

I specify that there are no edge weights (the empty array as third input, []) and store the coordinate information as "custom attributes" in the Nodes property of the graph by specifying it as the fourth input when I construct G.

G=simplify(graph(s, t, [], coordTable))

G = 
  graph with properties:

    Edges: [8×1 table]
    Nodes: [5×2 table]

If I'd omitted the coordinates table, G.Nodes would be 5-by-0.

GNoCoordinates = simplify(graph(s, t, []))

GNoCoordinates = 
  graph with properties:

    Edges: [8×1 table]
    Nodes: [5×0 table]

Now you can use that coordinate information when plotting the graph to put the nodes in the same place they would be if you plotted the delaunayTriangulation. The extra margin around the graph plot is due to different axes limits; the coordinates are the same.

plot(G, XData = G.Nodes.x, YData = G.Nodes.y)

figure
triplot(tri)

●

Because there are multiple shortest paths between the pair of edges that were constrained, the one that shortestpath finds may not be the ones that include the nodes delaunayTriangulation added. Valid shortest paths from 1 to 4 in G include (1, 2, 4), (1, 3, 4), and (1, 5, 4) and similar for 2 to 3.

for i=1:height(C)
    E=shortestpath(G,C(i,1), C(i,2))      %output edge constraints
end

E = 1×3

1 2 4

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

E = 1×3

2 1 3

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

You could extract the subgraph containing the nodes that form the constrained edges and the nodes introduced by delaunayTriangulation. For the first constraint:

constrainedEdge1 = C(1, :);
addedNodes = (height(V)+1):numnodes(G);

addedNodes are the nodes in G that don't correspond to elements of V, the nodes that were added by delaunayTriangulation when it split the constraints.

S = subgraph(G, [constrainedEdge1, addedNodes]);

In this case only nodes 1 and 4 (the constrained edge) and 5 (the node added by delaunayTriangulation) and edges between two of those three nodes are in S. Nodes 2 and 3 in G are not present in S as you can see in the plot.

figure
plot(S, XData = S.Nodes.x, YData = S.Nodes.y)

However, this subgraph extraction does renumber the nodes (node 3 in S was node 5 in G and similarly for node 2 in S and node 4 in G.) Adding node names to G would help, either by manipulating G directly or by modifying coordTable before using it to create G in the first place. Let's give node 1 name "1", node 2 name "2", etc. then re-extract the subgraph as before.

G.Nodes.Names = string(1:numnodes(G)).';
S = subgraph(G, [constrainedEdge1, addedNodes]);

When plotting I'll use the node names as the labels for the nodes.

figure
plot(S, ...
    XData = S.Nodes.x, ...
    YData = S.Nodes.y, ...
    NodeLabel = S.Nodes.Names);

Index into the names using the EndNodes of the edges to obtain the list of edges in S (using the names from the original data):

E = S.Nodes.Names(S.Edges.EndNodes)

E = 2×2 string array
    "1"    "5"
    "4"    "5"

Or in double:

D = double(E)

D = 2×2

1 5 4 5

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Matt J]

Because there are multiple shortest paths between the pair of edges that were constrained, the one that shortestpath finds may not be the ones that include the nodes delaunayTriangulation added

True, but only because you didn't add weights. Had you weighted the edges by their Euclidean lengths, there would be only one unique shortest path, and it would coincide with the straight line path added in by the delaunayTriangulation constraints.