improved decimation (elalish#1147)

* improved decimation

* only collapse intersected Boolean edges

* ensure retains verts are not removed

* collapse all short edges

* fix compile

* gotten part

* cleanup

* pca's patch

include/manifold/manifold.h

@@ -311,6 +311,7 @@ class Manifold {
   Manifold Warp(std::function<void(vec3&)>) const;
   Manifold WarpBatch(std::function<void(VecView<vec3>)>) const;
   Manifold SetTolerance(double) const;
+  Manifold Simplify(double tolerance = 0) const;
   ///@}
 
   /** @name Boolean

src/boolean_result.cpp

@@ -856,7 +856,7 @@ Manifold::Impl Boolean3::Result(OpType op) const {
 
   UpdateReference(outR, inP_, inQ_, invertQ);
 
-  outR.SimplifyTopology();
+  outR.SimplifyTopology(nPv + nQv);
   outR.RemoveUnreferencedVerts();
 
   if (ManifoldParams().intermediateChecks)

src/edge_op.cpp

@@ -48,26 +48,31 @@ struct ShortEdge {
   VecView<const Halfedge> halfedge;
   VecView<const vec3> vertPos;
   const double tolerance;
+  const int firstNewVert;
 
   bool operator()(int edge) const {
-    if (halfedge[edge].pairedHalfedge < 0) return false;
+    const Halfedge& half = halfedge[edge];
+    if (half.pairedHalfedge < 0 ||
+        (half.startVert < firstNewVert && half.endVert < firstNewVert))
+      return false;
     // Flag short edges
-    const vec3 delta =
-        vertPos[halfedge[edge].endVert] - vertPos[halfedge[edge].startVert];
+    const vec3 delta = vertPos[half.endVert] - vertPos[half.startVert];
     return la::dot(delta, delta) < tolerance * tolerance;
   }
 };
 
 struct FlagEdge {
   VecView<const Halfedge> halfedge;
   VecView<const TriRef> triRef;
+  const int firstNewVert;
 
   bool operator()(int edge) const {
-    if (halfedge[edge].pairedHalfedge < 0) return false;
+    const Halfedge& half = halfedge[edge];
+    if (half.pairedHalfedge < 0 || half.startVert < firstNewVert) return false;
     // Flag redundant edges - those where the startVert is surrounded by only
     // two original triangles.
     const TriRef ref0 = triRef[edge / 3];
-    int current = NextHalfedge(halfedge[edge].pairedHalfedge);
+    int current = NextHalfedge(half.pairedHalfedge);
     TriRef ref1 = triRef[current / 3];
     bool ref1Updated = !ref0.SameFace(ref1);
     while (current != edge) {
@@ -92,9 +97,15 @@ struct SwappableEdge {
   VecView<const vec3> vertPos;
   VecView<const vec3> triNormal;
   const double tolerance;
+  const int firstNewVert;
 
   bool operator()(int edge) const {
-    if (halfedge[edge].pairedHalfedge < 0) return false;
+    const Halfedge& half = halfedge[edge];
+    if (half.pairedHalfedge < 0) return false;
+    if (half.startVert < firstNewVert && half.endVert < firstNewVert &&
+        halfedge[NextHalfedge(edge)].endVert < firstNewVert &&
+        halfedge[NextHalfedge(half.pairedHalfedge)].endVert < firstNewVert)
+      return false;
 
     int tri = edge / 3;
     ivec3 triEdge = TriOf(edge);
@@ -106,7 +117,7 @@ struct SwappableEdge {
       return false;
 
     // Switch to neighbor's projection.
-    edge = halfedge[edge].pairedHalfedge;
+    edge = half.pairedHalfedge;
     tri = edge / 3;
     triEdge = TriOf(edge);
     projection = GetAxisAlignedProjection(triNormal[tri]);
@@ -295,7 +306,7 @@ void Manifold::Impl::CleanupTopology() {
  * Rather than actually removing the edges, this step merely marks them for
  * removal, by setting vertPos to NaN and halfedge to {-1, -1, -1, -1}.
  */
-void Manifold::Impl::SimplifyTopology() {
+void Manifold::Impl::SimplifyTopology(int firstNewVert) {
   if (!halfedge_.size()) return;
 
   CleanupTopology();
@@ -315,43 +326,43 @@ void Manifold::Impl::SimplifyTopology() {
   {
     ZoneScopedN("CollapseShortEdge");
     numFlagged = 0;
-    ShortEdge se{halfedge_, vertPos_, epsilon_};
+    ShortEdge se{halfedge_, vertPos_, epsilon_, firstNewVert};
     s.run(nbEdges, se, [&](size_t i) {
-      CollapseEdge(i, scratchBuffer);
+      const bool didCollapse = CollapseEdge(i, scratchBuffer);
+      if (didCollapse) numFlagged++;
       scratchBuffer.resize(0);
-      numFlagged++;
     });
-  }
 
 #ifdef MANIFOLD_DEBUG
-  if (ManifoldParams().verbose && numFlagged > 0) {
-    std::cout << "found " << numFlagged << " short edges to collapse"
-              << std::endl;
-  }
+    if (ManifoldParams().verbose && numFlagged > 0) {
+      std::cout << "collapsed " << numFlagged << " short edges" << std::endl;
+    }
 #endif
+  }
 
-  {
+  while (1) {
     ZoneScopedN("CollapseFlaggedEdge");
     numFlagged = 0;
-    FlagEdge se{halfedge_, meshRelation_.triRef};
+    FlagEdge se{halfedge_, meshRelation_.triRef, firstNewVert};
     s.run(nbEdges, se, [&](size_t i) {
-      CollapseEdge(i, scratchBuffer);
+      const bool didCollapse = CollapseEdge(i, scratchBuffer);
+      if (didCollapse) numFlagged++;
       scratchBuffer.resize(0);
-      numFlagged++;
     });
-  }
+    if (numFlagged == 0) break;
 
 #ifdef MANIFOLD_DEBUG
-  if (ManifoldParams().verbose && numFlagged > 0) {
-    std::cout << "found " << numFlagged << " colinear edges to collapse"
-              << std::endl;
-  }
+    if (ManifoldParams().verbose && numFlagged > 0) {
+      std::cout << "collapsed " << numFlagged << " colinear edges" << std::endl;
+    }
 #endif
+  }
 
   {
     ZoneScopedN("RecursiveEdgeSwap");
     numFlagged = 0;
-    SwappableEdge se{halfedge_, vertPos_, faceNormal_, tolerance_};
+    SwappableEdge se{halfedge_, vertPos_, faceNormal_, tolerance_,
+                     firstNewVert};
     std::vector<int> edgeSwapStack;
     std::vector<int> visited(halfedge_.size(), -1);
     int tag = 0;
@@ -365,13 +376,13 @@ void Manifold::Impl::SimplifyTopology() {
         RecursiveEdgeSwap(last, tag, visited, edgeSwapStack, scratchBuffer);
       }
     });
-  }
 
 #ifdef MANIFOLD_DEBUG
-  if (ManifoldParams().verbose && numFlagged > 0) {
-    std::cout << "found " << numFlagged << " edges to swap" << std::endl;
-  }
+    if (ManifoldParams().verbose && numFlagged > 0) {
+      std::cout << "swapped " << numFlagged << " edges" << std::endl;
+    }
 #endif
+  }
 }
 
 // Deduplicate the given 4-manifold edge by duplicating endVert, thus making the
@@ -549,16 +560,17 @@ void Manifold::Impl::RemoveIfFolded(int edge) {
   }
 }
 
-// Collapses the given edge by removing startVert. May split the mesh
-// topologically if the collapse would have resulted in a 4-manifold edge. Do
-// not collapse an edge if startVert is pinched - the vert will be marked NaN,
-// but other edges may still be pointing to it.
-void Manifold::Impl::CollapseEdge(const int edge, std::vector<int>& edges) {
+// Collapses the given edge by removing startVert - returns false if the edge
+// cannot be collapsed. May split the mesh topologically if the collapse would
+// have resulted in a 4-manifold edge. Do not collapse an edge if startVert is
+// pinched - the vert would be marked NaN, but other edges could still be
+// pointing to it.
+bool Manifold::Impl::CollapseEdge(const int edge, std::vector<int>& edges) {
   Vec<TriRef>& triRef = meshRelation_.triRef;
   Vec<ivec3>& triProp = meshRelation_.triProperties;
 
   const Halfedge toRemove = halfedge_[edge];
-  if (toRemove.pairedHalfedge < 0) return;
+  if (toRemove.pairedHalfedge < 0) return false;
 
   const int endVert = toRemove.endVert;
   const ivec3 tri0edge = TriOf(edge);
@@ -567,7 +579,7 @@ void Manifold::Impl::CollapseEdge(const int edge, std::vector<int>& edges) {
   const vec3 pNew = vertPos_[endVert];
   const vec3 pOld = vertPos_[toRemove.startVert];
   const vec3 delta = pNew - pOld;
-  const bool shortEdge = la::dot(delta, delta) < tolerance_ * tolerance_;
+  const bool shortEdge = la::dot(delta, delta) < epsilon_ * epsilon_;
 
   // Orbit endVert
   int current = halfedge_[tri0edge[1]].pairedHalfedge;
@@ -594,20 +606,14 @@ void Manifold::Impl::CollapseEdge(const int edge, std::vector<int>& edges) {
       if (!ref.SameFace(refCheck)) {
         refCheck = triRef[edge / 3];
         if (!ref.SameFace(refCheck)) {
-          return;
-        } else {
-          // Don't collapse if the edges separating the faces are not colinear
-          // (can happen when the two faces are coplanar).
-          if (CCW(projection * pOld, projection * pLast, projection * pNew,
-                  epsilon_) != 0)
-            return;
+          return false;
         }
       }
 
       // Don't collapse edge if it would cause a triangle to invert.
       if (CCW(projection * pNext, projection * pLast, projection * pNew,
               epsilon_) < 0)
-        return;
+        return false;
 
       pLast = pNext;
       current = halfedge_[current].pairedHalfedge;
@@ -654,6 +660,7 @@ void Manifold::Impl::CollapseEdge(const int edge, std::vector<int>& edges) {
   UpdateVert(endVert, start, tri0edge[2]);
   CollapseTri(tri0edge);
   RemoveIfFolded(start);
+  return true;
 }
 
 void Manifold::Impl::RecursiveEdgeSwap(const int edge, int& tag,

src/impl.h

@@ -329,9 +329,9 @@ struct Manifold::Impl {
 
   // edge_op.cpp
   void CleanupTopology();
-  void SimplifyTopology();
+  void SimplifyTopology(int firstNewVert = 0);
   void DedupeEdge(int edge);
-  void CollapseEdge(int edge, std::vector<int>& edges);
+  bool CollapseEdge(int edge, std::vector<int>& edges);
   void RecursiveEdgeSwap(int edge, int& tag, std::vector<int>& visited,
                          std::vector<int>& edgeSwapStack,
                          std::vector<int>& edges);

src/manifold.cpp

@@ -415,6 +415,25 @@ Manifold Manifold::SetTolerance(double tolerance) const {
   return Manifold(impl);
 }
 
+/**
+ * Return a copy of the manifold simplified to the given tolerance, but with its
+ * actual tolerance value unchanged. If no tolerance is given, the current
+ * tolerance is used for simplification.
+ */
+Manifold Manifold::Simplify(double tolerance) const {
+  auto impl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+  const double oldTolerance = impl->tolerance_;
+  if (tolerance == 0) tolerance = oldTolerance;
+  if (tolerance >= oldTolerance) {
+    impl->tolerance_ = tolerance;
+    impl->CreateFaces();
+    impl->SimplifyTopology();
+    impl->Finish();
+  }
+  impl->tolerance_ = oldTolerance;
+  return Manifold(impl);
+}
+
 /**
  * The genus is a topological property of the manifold, representing the number
  * of "handles". A sphere is 0, torus 1, etc. It is only meaningful for a single

test/boolean_complex_test.cpp

@@ -912,6 +912,8 @@ TEST(BooleanComplex, CraycloudBool) {
   Manifold m2 = ReadMesh("Cray_right.glb");
   Manifold res = m1 - m2;
   EXPECT_EQ(res.Status(), Manifold::Error::NoError);
+  EXPECT_FALSE(res.IsEmpty());
+  res = res.Simplify();
   EXPECT_TRUE(res.IsEmpty());
 }
 

test/boolean_test.cpp

@@ -138,6 +138,14 @@ TEST(Boolean, Cubes) {
 #endif
 }
 
+TEST(Boolean, Simplify) {
+  Manifold cube = Manifold::Cube().Refine(10);
+  Manifold result = cube + cube.Translate({1, 0, 0});
+  EXPECT_EQ(result.NumTri(), 1928);
+  result = result.Simplify();
+  EXPECT_EQ(result.NumTri(), 20);
+}
+
 TEST(Boolean, NoRetainedVerts) {
   Manifold cube = Manifold::Cube(vec3(1), true);
   Manifold oct = Manifold::Sphere(1, 4);

test/properties_test.cpp

@@ -77,14 +77,14 @@ TEST(Properties, Tolerance) {
 }
 
 TEST(Properties, ToleranceSphere) {
-  const int n = 100;
+  const int n = 1000;
   Manifold sphere = Manifold::Sphere(1, 4 * n);
   EXPECT_EQ(sphere.NumTri(), 8 * n * n);
 
   Manifold sphere2 = sphere.SetTolerance(0.01);
-  EXPECT_LT(sphere2.NumTri(), 15000);
-  EXPECT_NEAR(sphere.Volume(), sphere2.Volume(), 0.02);
-  EXPECT_NEAR(sphere.SurfaceArea(), sphere2.SurfaceArea(), 0.01);
+  EXPECT_LT(sphere2.NumTri(), 2500);
+  EXPECT_NEAR(sphere.Volume(), sphere2.Volume(), 0.05);
+  EXPECT_NEAR(sphere.SurfaceArea(), sphere2.SurfaceArea(), 0.05);
 #ifdef MANIFOLD_EXPORT
   if (options.exportModels) ExportMesh("sphere.glb", sphere2.GetMeshGL(), {});
 #endif