Merge pull request elalish#137 from elalish/infiniteLoop

Fixed edge collapse

.vscode/settings.json

@@ -86,11 +86,18 @@
         "set": "cpp",
         "memory_resource": "cpp",
         "random": "cpp",
-        "variant": "cpp"
+        "variant": "cpp",
+        "hash_map": "cpp",
+        "hash_set": "cpp",
+        "codecvt": "cpp",
+        "regex": "cpp",
+        "future": "cpp",
+        "shared_mutex": "cpp"
     },
     "C_Cpp.clang_format_fallbackStyle": "google",
     "editor.formatOnSave": true,
     "cmake.configureArgs": [
-        "-DTHRUST_BACKEND=CPP"
+        "-DMANIFOLD_USE_CUDA=OFF",
+        "-DMANIFOLD_PAR=NONE"
     ],
 }

manifold/src/edge_op.cpp

@@ -33,6 +33,17 @@ __host__ __device__ bool Is01Longest(glm::vec2 v0, glm::vec2 v1, glm::vec2 v2) {
   return l[0] > l[1] && l[0] > l[2];
 }
 
+struct DuplicateEdge {
+  const Halfedge* sortedHalfedge;
+
+  __host__ __device__ bool operator()(int edge) {
+    const Halfedge& halfedge = sortedHalfedge[edge];
+    const Halfedge& nextHalfedge = sortedHalfedge[edge + 1];
+    return halfedge.startVert == nextHalfedge.startVert &&
+           halfedge.endVert == nextHalfedge.endVert;
+  }
+};
+
 struct ShortEdge {
   const Halfedge* halfedge;
   const glm::vec3* vertPos;
@@ -110,6 +121,11 @@ namespace manifold {
  * It also performs edge swaps on the long edges of degenerate triangles, though
  * there are some configurations of degenerates that cannot be removed this way.
  *
+ * Before collapsing edges, the mesh is checked for duplicate edges (more than
+ * one pair of triangles sharing the same edge), which are removed by
+ * duplicating one vert and adding two triangles. These degenerate triangles are
+ * likely to be collapsed again in the subsequent simplification.
+ *
  * Note when an edge collapse would result in something non-manifold, the
  * vertices are duplicated in such a way as to remove handles or separate
  * meshes, thus decreasing the Genus(). It only increases when meshes that have
@@ -119,9 +135,26 @@ namespace manifold {
  * removal, by setting vertPos to NaN and halfedge to {-1, -1, -1, -1}.
  */
 void Manifold::Impl::SimplifyTopology() {
-  VecDH<int> flaggedEdges(halfedge_.size());
   auto policy = autoPolicy(halfedge_.size());
+
+  VecDH<Halfedge> halfedge(halfedge_);
+  VecDH<int> idx(halfedge_.size());
+  sequence(policy, idx.begin(), idx.end());
+  sort_by_key(policy, halfedge.begin(), halfedge.end(), idx.begin());
+
+  VecDH<int> flaggedEdges(halfedge_.size());
+
   int numFlagged =
+      copy_if<decltype(flaggedEdges.begin())>(
+          policy, idx.begin(), idx.end() - 1, countAt(0), flaggedEdges.begin(),
+          DuplicateEdge({halfedge.cptrD()})) -
+      flaggedEdges.begin();
+  flaggedEdges.resize(numFlagged);
+
+  for (const int edge : flaggedEdges) DedupeEdge(edge);
+
+  flaggedEdges.resize(halfedge_.size());
+  numFlagged =
       copy_if<decltype(flaggedEdges.begin())>(
           policy, countAt(0), countAt(halfedge_.size()), flaggedEdges.begin(),
           ShortEdge({halfedge_.cptrD(), vertPos_.cptrD(), precision_})) -
@@ -154,6 +187,47 @@ void Manifold::Impl::SimplifyTopology() {
   }
 }
 
+void Manifold::Impl::DedupeEdge(const int edge) {
+  // Orbit endVert
+  const int startVert = halfedge_[edge].startVert;
+  const int endVert = halfedge_[edge].endVert;
+  int current = halfedge_[NextHalfedge(edge)].pairedHalfedge;
+  while (current != edge) {
+    const int vert = halfedge_[current].startVert;
+    if (vert == startVert) {
+      int newVert = vertPos_.size();
+      vertPos_.push_back(vertPos_[endVert]);
+      current = halfedge_[NextHalfedge(current)].pairedHalfedge;
+      int opposite = halfedge_[NextHalfedge(edge)].pairedHalfedge;
+
+      UpdateVert(newVert, current, opposite);
+
+      int newHalfedge = halfedge_.size();
+      int newFace = newHalfedge / 3;
+      int outsideVert = halfedge_[current].startVert;
+      halfedge_.push_back({endVert, newVert, -1, newFace});
+      halfedge_.push_back({newVert, outsideVert, -1, newFace});
+      halfedge_.push_back({outsideVert, endVert, -1, newFace});
+      PairUp(newHalfedge + 2, halfedge_[current].pairedHalfedge);
+      PairUp(newHalfedge + 1, current);
+
+      newHalfedge += 3;
+      ++newFace;
+      outsideVert = halfedge_[opposite].startVert;
+      halfedge_.push_back({newVert, endVert, -1, newFace});
+      halfedge_.push_back({endVert, outsideVert, -1, newFace});
+      halfedge_.push_back({outsideVert, newVert, -1, newFace});
+      PairUp(newHalfedge + 2, halfedge_[opposite].pairedHalfedge);
+      PairUp(newHalfedge + 1, opposite);
+      PairUp(newHalfedge, newHalfedge - 3);
+
+      break;
+    }
+
+    current = halfedge_[NextHalfedge(current)].pairedHalfedge;
+  }
+}
+
 void Manifold::Impl::PairUp(int edge0, int edge1) {
   halfedge_[edge0].pairedHalfedge = edge1;
   halfedge_[edge1].pairedHalfedge = edge0;

manifold/src/face_op.cpp

@@ -162,7 +162,7 @@ void Manifold::Impl::Face2Tri(const VecDH<int>& faceEdge,
     }
   }
   faceNormal_ = std::move(triNormal);
-  CreateAndFixHalfedges(triVerts);
+  CreateHalfedges(triVerts);
 }
 
 /**

manifold/src/impl.cpp

@@ -132,41 +132,6 @@ struct LinkHalfedges {
   }
 };
 
-struct SwapHalfedges {
-  Halfedge* halfedges;
-  const TmpEdge* edges;
-
-  __host__ void operator()(int k) {
-    const int i = 2 * k;
-    const int j = i - 2;
-    const TmpEdge thisEdge = edges[i];
-    const TmpEdge lastEdge = edges[j];
-    if (thisEdge.first == lastEdge.first &&
-        thisEdge.second == lastEdge.second) {
-      const int swap0idx = thisEdge.halfedgeIdx;
-      Halfedge& swap0 = halfedges[swap0idx];
-      const int swap1idx = swap0.pairedHalfedge;
-      Halfedge& swap1 = halfedges[swap1idx];
-
-      const int next0idx = swap0idx + ((swap0idx + 1) % 3 == 0 ? -2 : 1);
-      const int next1idx = swap1idx + ((swap1idx + 1) % 3 == 0 ? -2 : 1);
-      Halfedge& next0 = halfedges[next0idx];
-      Halfedge& next1 = halfedges[next1idx];
-
-      next0.startVert = swap0.endVert = next1.endVert;
-      swap0.pairedHalfedge = next1.pairedHalfedge;
-      halfedges[swap0.pairedHalfedge].pairedHalfedge = swap0idx;
-
-      next1.startVert = swap1.endVert = next0.endVert;
-      swap1.pairedHalfedge = next0.pairedHalfedge;
-      halfedges[swap1.pairedHalfedge].pairedHalfedge = swap1idx;
-
-      next0.pairedHalfedge = next1idx;
-      next1.pairedHalfedge = next0idx;
-    }
-  }
-};
-
 struct InitializeBaryRef {
   const int meshID;
   const Halfedge* halfedge;
@@ -297,7 +262,7 @@ Manifold::Impl::Impl(const Mesh& mesh,
   CheckDevice();
   CalculateBBox();
   SetPrecision();
-  CreateAndFixHalfedges(mesh.triVerts);
+  CreateHalfedges(mesh.triVerts);
   ALWAYS_ASSERT(IsManifold(), topologyErr, "Input mesh is not manifold!");
   CalculateNormals();
   InitializeNewReference(triProperties, properties, propertyTolerance);
@@ -478,23 +443,6 @@ int Manifold::Impl::InitializeNewReference(
  * Create the halfedge_ data structure from an input triVerts array like Mesh.
  */
 void Manifold::Impl::CreateHalfedges(const VecDH<glm::ivec3>& triVerts) {
-  const int numTri = triVerts.size();
-  halfedge_.resize(3 * numTri);
-  VecDH<TmpEdge> edge(3 * numTri);
-  auto policy = autoPolicy(numTri);
-  for_each_n(policy, zip(countAt(0), triVerts.begin()), numTri,
-             Tri2Halfedges({halfedge_.ptrD(), edge.ptrD()}));
-  sort(policy, edge.begin(), edge.end());
-  for_each_n(policy, countAt(0), halfedge_.size() / 2,
-             LinkHalfedges({halfedge_.ptrD(), edge.cptrD()}));
-}
-
-/**
- * Create the halfedge_ data structure from an input triVerts array like Mesh.
- * Check that the input is an even-manifold, and if it is not 2-manifold,
- * perform edge swaps until it is. This is a host function.
- */
-void Manifold::Impl::CreateAndFixHalfedges(const VecDH<glm::ivec3>& triVerts) {
   const int numTri = triVerts.size();
   // drop the old value first to avoid copy
   halfedge_.resize(0);
@@ -506,14 +454,11 @@ void Manifold::Impl::CreateAndFixHalfedges(const VecDH<glm::ivec3>& triVerts) {
   // Stable sort is required here so that halfedges from the same face are
   // paired together (the triangles were created in face order). In some
   // degenerate situations the triangulator can add the same internal edge in
-  // two different faces, causing this edge to not be 2-manifold. We detect this
-  // and fix it by swapping one of the identical edges, so it is important that
-  // we have the edges paired according to their face.
+  // two different faces, causing this edge to not be 2-manifold. These are
+  // fixed by duplicating verts in SimplifyTopology.
   stable_sort(policy, edge.begin(), edge.end());
   thrust::for_each_n(thrust::host, countAt(0), halfedge_.size() / 2,
                      LinkHalfedges({halfedge_.ptrH(), edge.cptrH()}));
-  thrust::for_each(thrust::host, countAt(1), countAt(halfedge_.size() / 2),
-                   SwapHalfedges({halfedge_.ptrH(), edge.cptrH()}));
 }
 
 /**

manifold/src/impl.h

@@ -68,7 +68,6 @@ struct Manifold::Impl {
 
   void ReinitializeReference(int meshID);
   void CreateHalfedges(const VecDH<glm::ivec3>& triVerts);
-  void CreateAndFixHalfedges(const VecDH<glm::ivec3>& triVerts);
   void CalculateNormals();
   void UpdateMeshIDs(VecDH<int>& meshIDs, VecDH<int>& originalIDs,
                      int startTri = 0, int n = -1, int startID = 0);
@@ -109,6 +108,7 @@ struct Manifold::Impl {
 
   // edge_op.cu
   void SimplifyTopology();
+  void DedupeEdge(int edge);
   void CollapseEdge(int edge);
   void RecursiveEdgeSwap(int edge);
   void RemoveIfFolded(int edge);

manifold/src/properties.cpp

@@ -219,12 +219,14 @@ bool Manifold::Impl::IsManifold() const {
   auto policy = autoPolicy(halfedge_.size());
   bool isManifold = all_of(policy, countAt(0), countAt(halfedge_.size()),
                            CheckManifold({halfedge_.cptrD()}));
+  // std::cout << (isManifold ? "" : "Not ") << "Manifold" << std::endl;
 
   VecDH<Halfedge> halfedge(halfedge_);
   sort(policy, halfedge.begin(), halfedge.end());
-  isManifold &= all_of(policy, countAt(0), countAt(2 * NumEdge() - 1),
-                       NoDuplicates({halfedge.cptrD()}));
-  return isManifold;
+  bool noDupes = all_of(policy, countAt(0), countAt(2 * NumEdge() - 1),
+                        NoDuplicates({halfedge.cptrD()}));
+  // std::cout << (noDupes ? "" : "Not ") << "2-Manifold" << std::endl;
+  return isManifold && noDupes;
 }
 
 /**

test/mesh_test.cpp

@@ -781,24 +781,25 @@ TEST(Boolean, Gyroid) {
   Related(result, input, meshID2idx);
 }
 
-TEST(Boolean, DISABLED_Cylinders) {
-  Manifold rod = Manifold::Cylinder(1.0, 0.4, -1.0, 20);
+TEST(Boolean, Cylinders) {
+  Manifold rod = Manifold::Cylinder(1.0, 0.4, -1.0, 12);
   float arrays1[][12] = {
+      {0, 0, 1, 3,    //
+       -1, 0, 0, 3,   //
+       0, -1, 0, 6},  //
+      {0, 0, 1, 2,    //
+       -1, 0, 0, 3,   //
+       0, -1, 0, 8},  //
+
       {0, 0, 1, 1,    //
        -1, 0, 0, 2,   //
        0, -1, 0, 7},  //
       {1, 0, 0, 3,    //
        0, 1, 0, 2,    //
        0, 0, 1, 6},   //
-      {0, 0, 1, 3,    //
-       -1, 0, 0, 3,   //
-       0, -1, 0, 6},  //
       {0, 0, 1, 3,    //
        -1, 0, 0, 3,   //
        0, -1, 0, 7},  //
-      {0, 0, 1, 2,    //
-       -1, 0, 0, 3,   //
-       0, -1, 0, 8},  //
       {0, 0, 1, 1,    //
        -1, 0, 0, 3,   //
        0, -1, 0, 7},  //
@@ -810,15 +811,16 @@ TEST(Boolean, DISABLED_Cylinders) {
        0, -1, 0, 6},  //
   };
   float arrays2[][12] = {
-      {0, 0, 1, 2,    //
-       -1, 0, 0, 2,   //
-       0, -1, 0, 7},  //
       {1, 0, 0, 3,    //
        0, 0, 1, 2,    //
        0, -1, 0, 6},  //
       {1, 0, 0, 4,    //
        0, 1, 0, 3,    //
        0, 0, 1, 6},   //
+
+      {0, 0, 1, 2,    //
+       -1, 0, 0, 2,   //
+       0, -1, 0, 7},  //
       {1, 0, 0, 3,    //
        0, 1, 0, 3,    //
        0, 0, 1, 7},   //
@@ -860,10 +862,24 @@ TEST(Boolean, DISABLED_Cylinders) {
     }
     m2 += Manifold(rod).Transform(mat);
   }
-
   m1 += m2;
 
   EXPECT_TRUE(m1.IsManifold());
   EXPECT_TRUE(m1.MatchesTriNormals());
   EXPECT_LE(m1.NumDegenerateTris(), 12);
 }
+
+TEST(Boolean, Cubes) {
+  Manifold result = Manifold::Cube({1.2, 1, 1}, true).Translate({0, -0.5, 0.5});
+  result += Manifold::Cube({1, 0.8, 0.5}).Translate({-0.5, 0, 0.5});
+  result += Manifold::Cube({1.2, 0.1, 0.5}).Translate({-0.6, -0.1, 0});
+
+  EXPECT_TRUE(result.IsManifold());
+  EXPECT_TRUE(result.MatchesTriNormals());
+  EXPECT_LE(result.NumDegenerateTris(), 0);
+  auto prop = result.GetProperties();
+  EXPECT_NEAR(prop.volume, 1.6, 0.001);
+  EXPECT_NEAR(prop.surfaceArea, 9.2, 0.01);
+
+  if (options.exportModels) ExportMesh("cubes.glb", result.GetMesh(), {});
+}

test/test_main.cpp

@@ -27,8 +27,6 @@ void print_usage() {
 int main(int argc, char **argv) {
   ::testing::InitGoogleTest(&argc, argv);
 
-  int c;
-
   for (int i = 1; i < argc; i++) {
     if (argv[i][0] != '-') {
       fprintf(stderr, "Unknown option: %s\n", argv[i]);