Mesh Comparison

For mesh related unit tests, a comparison is made between the output of the code and a previously computed snapshot mesh which is known to be correct. This page explains how the mesh comparison is done.

Equality vs. Isomorphism

Say that you are given two algorithms that output a 3d mesh, how would you test that the two algorithms produce the same result? Naively, you would go through every vertex/edge/face and check that it is the same. Sadly, this doesn't really do what you want. Consider these two meshes consisting of a single quadrilateral (the numbers are the indices of the vertices/edges):

graph TD
 0((0)) ---|0| 1((1)); 
 2((2)) ---|1| 3((3));
 0 ---|2| 2;
 1 ---|3| 3;
 0'((1)) ---|0| 1'((0)); 
 2'((2)) ---|1| 3'((3));
 0' ---|2| 2';
 1' ---|3| 3';

Clearly, we should consider them as being the same. However, if we just naively check everything based on indices, we would come to the conclusion that they are different. Indeed, edge 2 of mesh 1 goes between vertices 0 and 2, while edge 2 of mesh 2 goes between vertices 1 and 2. In this case, we can see quite quickly that re-indexing the vertices of mesh 2 as follows would result in an identical mesh to mesh 1:

0 -> 1
1 -> 0
2 -> 2
3 -> 3

This is exactly the problem we have to solve: "Can you find a re-indexing of the vertices, edges, faces of mesh 1 such that it becomes identical to mesh 2?"

Here is a more mathematical way to formulate the problem. Let \(M_1 = (V_1, E_1, F_1)\) and \(M_2 = (V_2, E_2, F_2)\) be the two 3D meshes. We need to find invertible maps \(f_V \colon V_1 \to V_2\), \(f_E \colon E_1 \to E_2\) and \(f_F \colon F_1 \to F_2\) such that:

• If \(e = (v, v') \in E_1\) is an edge in the first mesh, then \(f_E(e) =(f_V(v), f_V(v')) \in E_2\) is an edge in the second mesh.

• If \(f = (e^1, e^2, \dots, e^n) \in F_1\) is a face in the first mesh, then \(f_F(f)=(f_E(e^1),f_E(e^2),\dots , f_E(e^n)) \in F_2\) is a face in the second mesh.

In this case we would say that the two meshes are isomorphic, which literally means "same shape".

In the example above we would have \(f_V(0) = 1, f_V(1) = 0, f_V(2) = 2\) and \(f_V(3) = 3\), while \(f_E(e_i) = e_i\). Then, indeed:

• \(f_E((0,1)) = f_E(e_0) = e_0 = (1,0) = (f_V(0), f_V(1))\)
• \(f_E((2,3)) = f_E(e_1) = e_1 = (2,3) = (f_V(2), f_V(3))\)
• \(f_E((0,2)) = f_E(e_2) = e_2 = (1,2) = (f_V(0), f_V(2))\)
• \(f_E((1,3)) = f_E(e_3) = e_3 = (0,3) = (f_V(1), f_V(3))\)

The algorithm

The general version of this problem is very hard (no polynomial time algorithm is known yet). The relevant search term is "Graph Isomorphism". Luckily when working with 3D meshes, we usually have more data available: mesh attributes. For our use-case, if the resulting positions of the 3D mesh got changed, then the mesh would be considered different. With this we can come up with an algorithm that can test for the existence of a mesh isomorphism relatively quickly.

Canonical/sorted indices

In practice, it turns out to be quite annoying to construct the isomorphism from \(M_1\) to \(M_2\) immediately. Insteadh we construct two isomorphisms to a "canonical" mesh \(M\). Then, if we know that \(M_1\) is isomorphic to \(M\), and \(M_2\) is isomorphic to \(M\) we can conclude that \(M_1\) and \(M_2\) are isomorphic as well by transitivity. The "canonical/sorted" mesh represents a mesh with the same vertices, edges, and faces as \(M_1\) or \(M_2\), except that the vertices, edges, and faces have been sorted based on the mesh attributes. If the sorted version of \(M_1\) does not equal the sorted version of \(M_2\), then we know that they are not isomorphic!

To keep track of this canonical mesh, we store 6 maps per domain:

• to_sorted1, to_sorted2 are maps going from the indices of \(M_1\) and \(M_2\) respectively to the index in the "sorted" mesh.
• from_sorted1, from_sorted2 are the inverse maps of to_sorted1 and to_sorted2. Given an index in the sorted mesh, they tell us what the index is in the original mesh.

• set_ids, set_sizes. When sorting, it's possible that two indices have the exact same attribute value, e.g., two vertices with the same position. In this case we can not yet distinguish one from the other. We treat them as being in the same "set". What the maps mean, is best seen with an example. Say we have an attribute with the following values:

  Index	Attribute Value
  0	3
  1	2
  2	1
  3	2
  4	2
  5	1
  6	3
  7	1
  8	1
  9	4

  After (stable) sorting, the table would look like this:

  Index	Attribute Value	Set ID	Set size
  2	1	0	4
  5	1	0	4
  7	1	0	4
  8	1	0	4
  1	2	4	3
  3	2	4	3
  4	2	4	3
  0	3	7	2
  6	3	7	2
  9	4	9	1

  The "ID" of a "set" is the index of its first element. The size is simply the number of elements in that set. Initially all the elements are in the same set. We don't need two copies of these maps, since the two meshes should result in the same sets, if they are isomorphic.

Here's a schematic diagram of the from/to_sorted maps:

flowchart LR
subgraph mesh1[Mesh 1]
v0((0)) ---|0| v1((1))
v2((2)) ---|1| v3((3))
v0 ---|2| v2
v1 ---|3| v3
end
subgraph mesh2[Mesh 2]
v0'((3)) ---|2| v1'((0))
v2'((1)) ---|3| v3'((2))
v0' ---|1| v2'
v1' ---|0| v3'
end
subgraph sorted_mesh[Sorted Mesh]
v0''((2)) ---|3| v1''((3))
v2''((0)) ---|0| v3''((1))
v0'' ---|1| v2''
v1'' ---|2| v3''
end
mesh1 -->|to_sorted1| sorted_mesh
mesh2 -->|to_sorted2| sorted_mesh
sorted_mesh -->|from_sorted1| mesh1
sorted_mesh -->|from_sorted2| mesh2

Concretely the maps would have these values for the vertex domain

index	from_sorted1	to_sorted1	from_sorted2	to_sorted2
0	2	2	1	3
1	3	3	2	0
2	0	0	3	1
3	1	1	0	2

and these values for the edge domain

index	from_sorted1	to_sorted1	from_sorted2	to_sorted2
0	1	3	3	2
1	2	0	1	1
2	3	1	0	3
3	0	2	2	0

Step 1 (simple verifications)

First thing we check is if the meshes have the same number of vertices, edges and face, and if they have the same type of attributes (== same name, domain, and type).

Step 2 (sorting)

The goal of the sorting is to make the set sizes as small as possible. We first sort the vertices, since edges and faces depend on the vertex ordering. We go through all the attributes on the domain, and sort the vertices based on the values of the attribute. Concretely:

• Look up an attribute A.
• Per set, we sort using the attribute values of A (this updates the from_sorted maps).
• We verify that the values of A agree on both meshes, i.e., that A1[from_sorted1[i]] == A2[from_sorted2[i]]
• We update the set ID's and set sizes.

After sorting the vertices we calculate the to_sorted maps from the from_sorted maps. For the other domains we do the same thing, except for some special case attributes:

• For edges, we sort the .edge_verts attribute ourselves, since the attribute values depend on the indices of the vertices. To do this, we create an ordered edge from each original edge as follows: if the edge has vertices v1 and v2, the ordered edge will be given vertices set_ids[to_sorted[v1]] and set_ids[to_sorted[v2]]. Recall that vertices in the same set were vertices which we couldn't distinguish yet. So, instead of treating edges as going from one specific vertex to another specific vertex, we treat them as going from one set of vertices to another set of vertices. We can then sort the edges based on the ordered edges.
• For the corners, we handle the .corner_vert and .corner_edge attributes separately. Instead of sorting based on the corner's vertex or edge directly, we sort based on the corner's vertex set or edge set using the set_ids[to_sorted[...]] composition like for sorting edges.

Step 3 (ensuring all sets have 1 element)

Once all the sets have a single element we are done, since we then can make the maps \(M_1 \to M_2\) as from_sorted2[to_sorted1[...]]. So, this final step's goal is to ensure that all the sets have a single element. Once we know that all the vertex sets have a single element, we can ensure that the other domain sets will also have single elements by sorting them again (since duplicate edges and faces are not allowed). To ensure that all the vertex sets have a single element, we make use of the vert_to_edge_map:

• Take a vertex in mesh 1, we need to find exactly 1 matching vertex in Mesh 2.
• The only possibilities are the vertices in the same set.
• For each vertex in the set, compare its edges to the edges of the vertex from Mesh 1. If they agree, add it as a candidate.
• If there are no candidates at the end, it means that the two meshes have different mesh topologies. If there is exactly one candidate, then that is our matching vertex. If there is more than one candidate, then we have a lot of overlapping vertices and edges in the original mesh. We make the simplifying assumption that this never happens (for the meshes we're interested in).

As such we match all the vertices one-by-one, and reduce all the vertex sets to size 1. After that, we recalculate the inverse maps and sort the edges, corners and faces one last time to ensure that the mesh topologies are the same.