Simplification & subdivision

include/geometrycentral/surface/mutation_manager.h

@@ -128,20 +128,23 @@ class MutationManager {
   // to `tSplit`.
   // In general, both the `tSplit` and the `newVertexPosition` are used; `tSplit` is necessary to allow callbacks to
   // interpolate data; if called with either the other will be inferred.
-  void splitEdge(Edge e, double tSplit);
-  void splitEdge(Edge e, Vector3 newVertexPosition);
-  void splitEdge(Edge e, double tSplit, Vector3 newVertexPosition);
+  // Returns the new halfedge which points away from the new vertex (so he.vertex() is the new vertex), and is the same
+  // direction as e.halfedge() on the original edge. The halfedge direction of the other part of the new split edge is
+  // also preserved.
+  Halfedge splitEdge(Edge e, double tSplit);
+  Halfedge splitEdge(Edge e, Vector3 newVertexPosition);
+  Halfedge splitEdge(Edge e, double tSplit, Vector3 newVertexPosition);
 
   // Collapse an edge.
-  // Returns true if the edge could actually be collapsed.
-  bool collapseEdge(Edge e, double tCollapse);
-  bool collapseEdge(Edge e, Vector3 newVertexPosition);
-  bool collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition);
-
-  // Split a face (i.e. insert a vertex into the face)
-  void splitFace(Face f, const std::vector<double>& bSplit);
-  void splitFace(Face f, Vector3 newVertexPosition);
-  void splitFace(Face f, const std::vector<double>& bSplit, Vector3 newVertexPosition);
+  // Returns the new vertex if the edge could be collapsed, and Vertex() otherwise
+  Vertex collapseEdge(Edge e, double tCollapse);
+  Vertex collapseEdge(Edge e, Vector3 newVertexPosition);
+  Vertex collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition);
+
+  // Split a face (i.e. insert a vertex into the face), and return the new vertex
+  Vertex splitFace(Face f, const std::vector<double>& bSplit);
+  Vertex splitFace(Face f, Vector3 newVertexPosition);
+  Vertex splitFace(Face f, const std::vector<double>& bSplit, Vector3 newVertexPosition);
 
   // ======================================================
   // ======== High-level mutations

include/geometrycentral/surface/quadric_error_simplification.h

@@ -0,0 +1,37 @@
+#pragma once
+
+#include "geometrycentral/surface/manifold_surface_mesh.h"
+#include "geometrycentral/surface/mutation_manager.h"
+#include "geometrycentral/surface/vertex_position_geometry.h"
+
+#include "geometrycentral/utilities/elementary_geometry.h"
+
+#include <queue>
+
+namespace geometrycentral {
+namespace surface {
+class Quadric {
+public:
+  Quadric();
+  Quadric(const Eigen::Matrix3d& A_, const Eigen::Vector3d& b_, double c_);
+  Quadric(const Quadric& Q1, const Quadric& Q2);
+  double cost(const Eigen::Vector3d& v);
+  Eigen::Vector3d optimalPoint();
+
+  Quadric operator+=(const Quadric& Q);
+
+protected:
+  Eigen::Matrix3d A;
+  Eigen::Vector3d b;
+  double c;
+
+  friend Quadric operator+(const Quadric& Q1, const Quadric& Q2);
+};
+
+Quadric operator+(const Quadric& Q1, const Quadric& Q2);
+
+void quadricErrorSimplify(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo, double tol = 0.05);
+void quadricErrorSimplify(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo, double tol, MutationManager& mm);
+
+} // namespace surface
+} // namespace geometrycentral

include/geometrycentral/surface/subdivide.h

@@ -0,0 +1,19 @@
+#pragma once
+
+#include "geometrycentral/surface/manifold_surface_mesh.h"
+#include "geometrycentral/surface/mutation_manager.h"
+#include "geometrycentral/surface/vertex_position_geometry.h"
+
+namespace geometrycentral {
+namespace surface {
+
+// Subdivide to form quad mesh
+// Note: these do not take MutationManager arguments since MutationManager only supports triangular meshes at the moment
+void linearSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo);
+void catmullClarkSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo);
+
+void loopSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo);
+void loopSubdivide(ManifoldSurfaceMesh& mesh, VertexPositionGeometry& geo, MutationManager& mm);
+
+} // namespace surface
+} // namespace geometrycentral

src/CMakeLists.txt

@@ -46,6 +46,8 @@ SET(SRCS
   surface/tufted_laplacian.cpp
   surface/flip_geodesics.cpp
   surface/transfer_functions.cpp
+  surface/quadric_error_simplification.cpp
+  surface/subdivide.cpp
   #surface/detect_symmetry.cpp
   #surface/mesh_ray_tracer.cpp
 
@@ -105,13 +107,15 @@ SET(HEADERS
   ${INCLUDE_ROOT}/surface/mesh_graph_algorithms.h
   ${INCLUDE_ROOT}/surface/mesh_ray_tracer.h
   ${INCLUDE_ROOT}/surface/parameterize.h
+  ${INCLUDE_ROOT}/surface/quadric_error_simplification.h
   ${INCLUDE_ROOT}/surface/rich_surface_mesh_data.h
   ${INCLUDE_ROOT}/surface/rich_surface_mesh_data.ipp
   ${INCLUDE_ROOT}/surface/polygon_soup_mesh.h
   ${INCLUDE_ROOT}/surface/signpost_intrinsic_triangulation.h
   ${INCLUDE_ROOT}/surface/signpost_intrinsic_triangulation.ipp
   ${INCLUDE_ROOT}/surface/simple_idt.h
   ${INCLUDE_ROOT}/surface/simple_polygon_mesh.h
+  ${INCLUDE_ROOT}/surface/subdivide.h
   ${INCLUDE_ROOT}/surface/surface_centers.h
   ${INCLUDE_ROOT}/surface/surface_mesh.h
   ${INCLUDE_ROOT}/surface/surface_mesh.ipp

src/surface/manifold_surface_mesh.cpp

@@ -1263,6 +1263,64 @@ Vertex ManifoldSurfaceMesh::collapseEdgeTriangular(Edge e) {
   modificationTick++;
 }
 
+Face ManifoldSurfaceMesh::removeEdge(Edge e) {
+  if (e.isBoundary()) {
+    throw std::runtime_error("not implemented");
+  }
+
+  // Halfedges/edges/faces that will be removed
+  // (except first face)
+  std::vector<Halfedge> toRemove{e.halfedge(), e.halfedge().twin()};
+  std::vector<Halfedge> ringHalfedges;
+  for (Halfedge heStart : toRemove) {
+    Halfedge he = heStart.next();
+    while (he != heStart) {
+      // The one-ring must not contain any other copies of e, or we cannot remove the edge
+      if (he.edge() == e) {
+        return Face();
+      }
+      ringHalfedges.push_back(he);
+      he = he.next();
+    }
+  }
+
+  // If both faces are the same, we cannot remove the edge. This should have been caught above
+  if (toRemove[0].face() == toRemove[1].face()) {
+    return Face();
+  }
+  Face keepFace = toRemove[0].face();
+
+
+  // Record these before we break pointers
+  Vertex src = e.halfedge().vertex();
+  Vertex dst = e.halfedge().twin().vertex();
+  Halfedge altSrcHedge = e.halfedge().twin().next();
+  Halfedge altDstHedge = e.halfedge().next();
+
+  // Hook up next and face refs for the halfedges along the ring
+  size_t N = ringHalfedges.size();
+  for (size_t i = 0; i < N; i++) {
+    heNextArr[ringHalfedges[i].getIndex()] = ringHalfedges[(i + 1) % N].getIndex();
+    heFaceArr[ringHalfedges[i].getIndex()] = keepFace.getIndex();
+  }
+
+  // only update vHalfedgeArr if needed to avoid disturbing boundary halfedges
+  if (src.halfedge().edge() == e) {
+    vHalfedgeArr[src.getIndex()] = altSrcHedge.getIndex();
+  }
+  if (dst.halfedge().edge() == e) {
+    vHalfedgeArr[dst.getIndex()] = altDstHedge.getIndex();
+  }
+
+  fHalfedgeArr[keepFace.getIndex()] = ringHalfedges[0].getIndex();
+
+  // Actually delete all of the elements
+  deleteElement(toRemove[1].face());
+  deleteEdgeBundle(e);
+
+  modificationTick++;
+  return keepFace;
+}
 
 
 bool ManifoldSurfaceMesh::removeFaceAlongBoundary(Face f) {

src/surface/mutation_manager.cpp

@@ -59,16 +59,16 @@ bool MutationManager::flipEdge(Edge e) {
   return true;
 }
 
-void MutationManager::splitEdge(Edge e, double tSplit) {
+Halfedge MutationManager::splitEdge(Edge e, double tSplit) {
   Vector3 newPos{0., 0., 0.};
   if (geometry) {
     VertexData<Vector3>& pos = geometry->vertexPositions;
     newPos = (1. - tSplit) * pos[e.halfedge().tailVertex()] + tSplit * pos[e.halfedge().tipVertex()];
   }
-  splitEdge(e, tSplit, newPos);
+  return splitEdge(e, tSplit, newPos);
 }
 
-void MutationManager::splitEdge(Edge e, Vector3 newVertexPosition) {
+Halfedge MutationManager::splitEdge(Edge e, Vector3 newVertexPosition) {
 
   double tSplit = -1;
   GC_SAFETY_ASSERT(geometry, "must have geometry to split by position");
@@ -80,10 +80,10 @@ void MutationManager::splitEdge(Edge e, Vector3 newVertexPosition) {
     tSplit = pointLineSegmentNeaestLocation(newVertexPosition, posTail, posTip);
   }
 
-  splitEdge(e, tSplit, newVertexPosition);
+  return splitEdge(e, tSplit, newVertexPosition);
 }
 
-void MutationManager::splitEdge(Edge e, double tSplit, Vector3 newVertexPosition) {
+Halfedge MutationManager::splitEdge(Edge e, double tSplit, Vector3 newVertexPosition) {
 
   // Invoke before callbacks
   for (EdgeSplitPolicy* policy : edgeSplitPolicies) {
@@ -102,11 +102,13 @@ void MutationManager::splitEdge(Edge e, double tSplit, Vector3 newVertexPosition
   for (EdgeSplitPolicy* policy : edgeSplitPolicies) {
     policy->afterEdgeSplit(newHeFront, newHeBack, tSplit);
   }
+
+  return newHeFront;
 }
 
 // Collapse an edge.
-// Returns true if the edge could actually be collapsed.
-bool MutationManager::collapseEdge(Edge e, double tCollapse) {
+// Returns the new vertex if the edge could be collapsed, and Vertex() otherwise
+Vertex MutationManager::collapseEdge(Edge e, double tCollapse) {
   Vector3 newPos{0., 0., 0.};
   if (geometry) {
     // Find the nearest tCoord
@@ -116,7 +118,7 @@ bool MutationManager::collapseEdge(Edge e, double tCollapse) {
   return collapseEdge(e, tCollapse, newPos);
 }
 
-bool MutationManager::collapseEdge(Edge e, Vector3 newVertexPosition) {
+Vertex MutationManager::collapseEdge(Edge e, Vector3 newVertexPosition) {
 
   double tCollapse = -1;
   GC_SAFETY_ASSERT(geometry, "must have geometry to split by position");
@@ -130,8 +132,7 @@ bool MutationManager::collapseEdge(Edge e, Vector3 newVertexPosition) {
 
   return collapseEdge(e, tCollapse, newVertexPosition);
 }
-
-bool MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition) {
+Vertex MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertexPosition) {
 
   // Invoke before callbacks
   // TODO need to handle possiblity that collapse fails -- check before calling
@@ -140,7 +141,7 @@ bool MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertexPo
   }
 
   Vertex newV = mesh.collapseEdgeTriangular(e);
-  if (newV == Vertex()) return false;
+  if (newV == Vertex()) return Vertex();
 
   if (geometry) {
     VertexData<Vector3>& pos = geometry->vertexPositions;
@@ -152,11 +153,11 @@ bool MutationManager::collapseEdge(Edge e, double tCollapse, Vector3 newVertexPo
     policy->afterEdgeCollapse(newV, tCollapse);
   }
 
-  return true;
+  return newV;
 }
 
 // Split a face (i.e. insert a vertex into the face)
-void MutationManager::splitFace(Face f, const std::vector<double>& bSplit) {
+Vertex MutationManager::splitFace(Face f, const std::vector<double>& bSplit) {
   Vector3 newPos = Vector3::zero();
   if (geometry) {
     size_t iV = 0;
@@ -167,16 +168,17 @@ void MutationManager::splitFace(Face f, const std::vector<double>& bSplit) {
     }
   }
 
-  splitFace(f, bSplit, newPos);
+  return splitFace(f, bSplit, newPos);
 }
 
-void MutationManager::splitFace(Face f, Vector3 newVertexPosition) {
+Vertex MutationManager::splitFace(Face f, Vector3 newVertexPosition) {
   // TODO
   throw std::runtime_error("Face split based on vertex position not implemented yet");
+  return Vertex();
 }
 
 
-void MutationManager::splitFace(Face f, const std::vector<double>& bSplit, Vector3 newVertexPosition) {
+Vertex MutationManager::splitFace(Face f, const std::vector<double>& bSplit, Vector3 newVertexPosition) {
   // Invoke before callbacks
   for (FaceSplitPolicy* policy : faceSplitPolicies) {
     policy->beforeFaceSplit(f, bSplit);
@@ -192,6 +194,8 @@ void MutationManager::splitFace(Face f, const std::vector<double>& bSplit, Vecto
   for (FaceSplitPolicy* policy : faceSplitPolicies) {
     policy->afterFaceSplit(newV, bSplit);
   }
+
+  return newV;
 }
 
 
@@ -290,7 +294,7 @@ MutationPolicyHandle MutationManager::registerPolicy(MutationPolicy* policyObjec
 }
 
 void MutationManager::removePolicy(const MutationPolicyHandle& toRemove) {
-  // Remove from all lists 
+  // Remove from all lists
   removeFromVector(vertexRepositionPolicies, toRemove.policy);
   removeFromVector(edgeFlipPolicies, toRemove.policy);
   removeFromVector(edgeSplitPolicies, toRemove.policy);