collapse: add fixup callback restore halfedge data

include/geometrycentral/surface/manifold_surface_mesh.h

@@ -2,9 +2,13 @@
 
 #include "geometrycentral/surface/surface_mesh.h"
 
+#include <functional>
+
 namespace geometrycentral {
 namespace surface {
 
+using EdgeCollapseFixupCallback = std::function<void(Halfedge, Halfedge)>;
+
 class ManifoldSurfaceMesh : public SurfaceMesh {
 
 public:
@@ -68,7 +72,42 @@ class ManifoldSurfaceMesh : public SurfaceMesh {
 
   // Collapse an edge. Returns the vertex adjacent to that edge which still exists. Returns Vertex() if not
   // collapsible. Assumes triangular simplicial complex as input (at least in neighborhood of collapse).
-  Vertex collapseEdgeTriangular(Edge e);
+  //
+  // However, note that when we edge collapse, we also remove the degenerate face(s) incident to the edge,
+  // which results in removing further edge(s):
+  //  |        x           |
+  // edge degenerate-face edge
+  //
+  //  |        x           +
+  // edge degenerate-face deleted-edge
+  //
+  // This destroys halfedge data that we still need in the neighborhood post collapse:
+  //
+  //  he1 | he2        x            he3 | he4
+  //     edge      degenerate-face     edge
+  //
+  //  he1 | he2        x            -- | --
+  //     edge      degenerate-face    deleted-edge
+  //
+  //  he1 |  he2
+  //     edge (after)
+  //
+  // he2 and he3 become completely unnecessary post collapse, but we've lost he4, which disrupts invariants the user may have wanted to maintain, such as parameterization boundaries (e.g., if he1 and he4 were on opposite sides of a UV boundary).
+  //
+  // Although we'd like to do so, we cannot directly perform
+  //
+  //  he1 | --         x             -- | he4
+  //     edge      degenerate-face     edge
+  //
+  // because this violates the implicit sibling property of manifold meshes (a part of which seems to assume that abs(he - he.twin) = 1).
+  //
+  // Also, the data may be stored in parallel structures outside this mesh object.
+  // Therefore, we should allow the user to provide a callback, |fixup|, that performs, essentially, for all halfedge/corner containers C,
+  //
+  // C[he2] = C[he4] (and analogs)
+  //
+  // after collapse.
+  Vertex collapseEdgeTriangular(Edge e, EdgeCollapseFixupCallback fixup = {});
 
   // Removes a vertex, leaving a high-degree face. If the input is a boundary vertex, preserves an edge along the
   // boundary. Return Face() if impossible (generally because doing so would make a manifold mesh nonmanifold).

src/surface/manifold_surface_mesh.cpp

@@ -1037,7 +1037,7 @@ Vertex ManifoldSurfaceMesh::insertVertex(Face fIn) {
 }
 
 
-Vertex ManifoldSurfaceMesh::collapseEdgeTriangular(Edge e) {
+Vertex ManifoldSurfaceMesh::collapseEdgeTriangular(Edge e, EdgeCollapseFixupCallback fixup) {
   /*  must maintain these
       std::vector<size_t> heNextArr;    // he.next(), forms a circular singly-linked list in each face
       std::vector<size_t> heVertexArr;  // he.vertex()
@@ -1124,6 +1124,12 @@ Vertex ManifoldSurfaceMesh::collapseEdgeTriangular(Edge e) {
     deleteElement(vA);
     deleteElement(fA);
 
+    if (fixup) {
+      // We destroyed heA2.edge(), but that also destroys heC2.
+      // heC2 (aka heC1.next) shows up in the next configuration as heA1.
+      fixup(heA1, heC2);
+    }
+
     return vB;
   }
 
@@ -1208,6 +1214,15 @@ Vertex ManifoldSurfaceMesh::collapseEdgeTriangular(Edge e) {
       deleteElement(fA);
       deleteElement(fB);
 
+      if (fixup) {
+        // We destroyed heB1.edge() and heA2.edge() along with their halfedges,
+        // but that also destroys heD1 and heC2 halfedges.
+        // They become heB2 and heA1 in the new configuration.
+        // We need to restore their data.
+        fixup(heB2, heD1);
+        fixup(heA1, heC2);
+      }
+
       return vB;
     }
 
@@ -1255,6 +1270,13 @@ Vertex ManifoldSurfaceMesh::collapseEdgeTriangular(Edge e) {
       deleteElement(fA);
       deleteElement(fB);
 
+      if (fixup) {
+        // We destroyed e, heB1, and heA2 along with their halfedges,
+        // and A2 and B1 are gone. However, that leaves us to reassign
+        // C1, because destroying B1 also destroyed C1, and made it become B2.
+        // We need to fix up data in B2 to match C1.
+        fixup(heB2, heC1);
+      }
       return vB;
     }