# Mesh 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 (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 invertable 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 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 algorithm is very hard (no polynomial time algorithm is know yet). The relevant search term is ["Graph Isomorphism"]( https://en.wikipedia.org/wiki/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. Instead 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 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/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 obtain a matching vertex for each vertex, 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.