Better faces (#1040)

* refactor * separate propVert from createFaces * working, but very slow * no longer slow * added test * loosen test * merging * update nix actions * fix --------- Co-authored-by: pca006132 <[email protected]>
elalish · Feb 13, 2025 · f2da5ea · f2da5ea
1 parent 447921c
commit f2da5ea
Show file tree

Hide file tree

Showing 6 changed files with 115 additions and 208 deletions.
diff --git a/.github/workflows/manifold.yml b/.github/workflows/manifold.yml
@@ -290,6 +290,5 @@ jobs:
     steps:
     - uses: actions/checkout@v4
     - uses: DeterminateSystems/nix-installer-action@main
-    - uses: DeterminateSystems/magic-nix-cache-action@main
     - run: nix build -L '.?submodules=1#${{matrix.variant}}'
 
diff --git a/flake.lock b/flake.lock
diff --git a/flake.nix b/flake.nix
@@ -38,34 +38,6 @@
                 version = clipper2-src.rev;
                 src = clipper2-src;
               });
-              # https://github.com/NixOS/nixpkgs/pull/343743#issuecomment-2424163602
-              binaryen =
-                let
-                  testsuite = final.fetchFromGitHub {
-                    owner = "WebAssembly";
-                    repo = "testsuite";
-                    rev = "e05365077e13a1d86ffe77acfb1a835b7aa78422";
-                    hash = "sha256-yvZ5AZTPUA6nsD3xpFC0VLthiu2CxVto66RTXBXXeJM=";
-                  };
-                in
-                prev.binaryen.overrideAttrs (_: rec {
-                  version = "119";
-                  src = pkgs.fetchFromGitHub {
-                    owner = "WebAssembly";
-                    repo = "binaryen";
-                    rev = "version_${version}";
-                    hash = "sha256-JYXtN3CW4qm/nnjGRvv3GxQ0x9O9wHtNYQLqHIYTTOA=";
-                  };
-                  preConfigure = ''
-                    if [ $doCheck -eq 1 ]; then
-                      sed -i '/googletest/d' third_party/CMakeLists.txt
-                      rmdir test/spec/testsuite
-                      ln -s ${testsuite} test/spec/testsuite
-                    else
-                      cmakeFlagsArray=($cmakeFlagsArray -DBUILD_TESTS=0)
-                    fi
-                  '';
-                });
             })
           ];
         };

diff --git a/src/impl.cpp b/src/impl.cpp
@@ -97,98 +97,6 @@ struct UpdateMeshID {
   void operator()(TriRef& ref) { ref.meshID = meshIDold2new[ref.meshID]; }
 };
 
-struct CoplanarEdge {
-  VecView<std::pair<int, int>> face2face;
-  VecView<double> triArea;
-  VecView<const Halfedge> halfedge;
-  VecView<const vec3> vertPos;
-  VecView<const TriRef> triRef;
-  VecView<const ivec3> triProp;
-  const int numProp;
-  const double epsilon;
-  const double tolerance;
-
-  void operator()(const int edgeIdx) {
-    const Halfedge edge = halfedge[edgeIdx];
-    const Halfedge pair = halfedge[edge.pairedHalfedge];
-    const int edgeFace = edgeIdx / 3;
-    const int pairFace = edge.pairedHalfedge / 3;
-
-    if (triRef[edgeFace].meshID != triRef[pairFace].meshID) return;
-
-    const vec3 base = vertPos[edge.startVert];
-    const int baseNum = edgeIdx - 3 * edgeFace;
-    const int jointNum = edge.pairedHalfedge - 3 * pairFace;
-
-    if (numProp > 0) {
-      if (triProp[edgeFace][baseNum] != triProp[pairFace][Next3(jointNum)] ||
-          triProp[edgeFace][Next3(baseNum)] != triProp[pairFace][jointNum])
-        return;
-    }
-
-    if (!edge.IsForward()) return;
-
-    const int edgeNum = baseNum == 0 ? 2 : baseNum - 1;
-    const int pairNum = jointNum == 0 ? 2 : jointNum - 1;
-    const vec3 jointVec = vertPos[pair.startVert] - base;
-    const vec3 edgeVec =
-        vertPos[halfedge[3 * edgeFace + edgeNum].startVert] - base;
-    const vec3 pairVec =
-        vertPos[halfedge[3 * pairFace + pairNum].startVert] - base;
-
-    const double length = std::max(la::length(jointVec), la::length(edgeVec));
-    const double lengthPair =
-        std::max(la::length(jointVec), la::length(pairVec));
-    vec3 normal = la::cross(jointVec, edgeVec);
-    const double area = la::length(normal);
-    const double areaPair = la::length(la::cross(pairVec, jointVec));
-
-    // make sure we only write this once
-    if (edgeIdx % 3 == 0) triArea[edgeFace] = area;
-    // Don't link degenerate triangles
-    if (area < length * epsilon || areaPair < lengthPair * epsilon) return;
-
-    const double volume = std::abs(la::dot(normal, pairVec));
-    // Only operate on coplanar triangles
-    if (volume > std::max(area, areaPair) * tolerance) return;
-
-    face2face[edgeIdx] = std::make_pair(edgeFace, pairFace);
-  }
-};
-
-struct CheckCoplanarity {
-  VecView<int> comp2tri;
-  VecView<const Halfedge> halfedge;
-  VecView<const vec3> vertPos;
-  std::vector<int>* components;
-  const double tolerance;
-
-  void operator()(int tri) {
-    const int component = (*components)[tri];
-    const int referenceTri =
-        reinterpret_cast<std::atomic<int>*>(&comp2tri[component])
-            ->load(std::memory_order_relaxed);
-    if (referenceTri < 0 || referenceTri == tri) return;
-
-    const vec3 origin = vertPos[halfedge[3 * referenceTri].startVert];
-    const vec3 normal = la::normalize(
-        la::cross(vertPos[halfedge[3 * referenceTri + 1].startVert] - origin,
-                  vertPos[halfedge[3 * referenceTri + 2].startVert] - origin));
-
-    for (const int i : {0, 1, 2}) {
-      const vec3 vert = vertPos[halfedge[3 * tri + i].startVert];
-      // If any component vertex is not coplanar with the component's reference
-      // triangle, unmark the entire component so that none of its triangles are
-      // marked coplanar.
-      if (std::abs(la::dot(normal, vert - origin)) > tolerance) {
-        reinterpret_cast<std::atomic<int>*>(&comp2tri[component])
-            ->store(-1, std::memory_order_relaxed);
-        break;
-      }
-    }
-  }
-};
-
 int GetLabels(std::vector<int>& components,
               const Vec<std::pair<int, int>>& edges, int numNodes) {
   UnionFind<> uf(numNodes);
@@ -199,19 +107,6 @@ int GetLabels(std::vector<int>& components,
 
   return uf.connectedComponents(components);
 }
-
-void DedupePropVerts(manifold::Vec<ivec3>& triProp,
-                     const Vec<std::pair<int, int>>& vert2vert,
-                     size_t numPropVert) {
-  ZoneScoped;
-  std::vector<int> vertLabels;
-  const int numLabels = GetLabels(vertLabels, vert2vert, numPropVert);
-
-  std::vector<int> label2vert(numLabels);
-  for (size_t v = 0; v < numPropVert; ++v) label2vert[vertLabels[v]] = v;
-  for (auto& prop : triProp)
-    for (int i : {0, 1, 2}) prop[i] = label2vert[vertLabels[prop[i]]];
-}
 }  // namespace
 
 namespace manifold {
@@ -309,77 +204,102 @@ void Manifold::Impl::InitializeOriginal(bool keepFaceID) {
 
 void Manifold::Impl::CreateFaces() {
   ZoneScoped;
-  Vec<std::pair<int, int>> face2face(halfedge_.size(), {-1, -1});
-  Vec<std::pair<int, int>> vert2vert(halfedge_.size(), {-1, -1});
-  Vec<double> triArea(NumTri());
-
-  const size_t numProp = NumProp();
-  if (numProp > 0) {
-    for_each_n(
-        autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
-        [&vert2vert, numProp, this](const int edgeIdx) {
-          const Halfedge edge = halfedge_[edgeIdx];
-          const Halfedge pair = halfedge_[edge.pairedHalfedge];
-          const int edgeFace = edgeIdx / 3;
-          const int pairFace = edge.pairedHalfedge / 3;
-
-          if (meshRelation_.triRef[edgeFace].meshID !=
-              meshRelation_.triRef[pairFace].meshID)
-            return;
-
-          const int baseNum = edgeIdx - 3 * edgeFace;
-          const int jointNum = edge.pairedHalfedge - 3 * pairFace;
-
-          const int prop0 = meshRelation_.triProperties[edgeFace][baseNum];
-          const int prop1 =
-              meshRelation_
-                  .triProperties[pairFace][jointNum == 2 ? 0 : jointNum + 1];
-          if (prop0 == prop1) return;
-
-          bool propEqual = true;
-          for (size_t p = 0; p < numProp; ++p) {
-            if (meshRelation_.properties[numProp * prop0 + p] !=
-                meshRelation_.properties[numProp * prop1 + p]) {
-              propEqual = false;
-              break;
-            }
-          }
-          if (propEqual) {
-            vert2vert[edgeIdx] = std::make_pair(prop0, prop1);
-          }
-        });
-    DedupePropVerts(meshRelation_.triProperties, vert2vert, NumPropVert());
-  }
+  const int numTri = NumTri();
+  struct TriPriority {
+    double area2;
+    int tri;
+  };
+  Vec<TriPriority> triPriority(numTri);
+  for_each_n(autoPolicy(numTri), countAt(0), numTri,
+             [&triPriority, this](int tri) {
+               meshRelation_.triRef[tri].faceID = -1;
+               const vec3 v = vertPos_[halfedge_[3 * tri].startVert];
+               triPriority[tri] = {
+                   length2(cross(vertPos_[halfedge_[3 * tri].endVert] - v,
+                                 vertPos_[halfedge_[3 * tri + 1].endVert] - v)),
+                   tri};
+             });
 
-  for_each_n(autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
-             CoplanarEdge({face2face, triArea, halfedge_, vertPos_,
-                           meshRelation_.triRef, meshRelation_.triProperties,
-                           meshRelation_.numProp, epsilon_, tolerance_}));
-
-  std::vector<int> components;
-  const int numComponent = GetLabels(components, face2face, NumTri());
-
-  Vec<int> comp2tri(numComponent, -1);
-  for (size_t tri = 0; tri < NumTri(); ++tri) {
-    const int comp = components[tri];
-    const int current = comp2tri[comp];
-    if (current < 0 || triArea[tri] > triArea[current]) {
-      comp2tri[comp] = tri;
-      triArea[comp] = triArea[tri];
+  stable_sort(triPriority.begin(), triPriority.end(),
+              [](auto a, auto b) { return a.area2 > b.area2; });
+
+  for (const auto tp : triPriority) {
+    if (meshRelation_.triRef[tp.tri].faceID >= 0) continue;
+
+    meshRelation_.triRef[tp.tri].faceID = tp.tri;
+    const vec3 base = vertPos_[halfedge_[3 * tp.tri].startVert];
+    const vec3 normal = faceNormal_[tp.tri];
+    std::vector<int> interiorHalfedges = {3 * tp.tri, 3 * tp.tri + 1,
+                                          3 * tp.tri + 2};
+    while (!interiorHalfedges.empty()) {
+      const int h =
+          NextHalfedge(halfedge_[interiorHalfedges.back()].pairedHalfedge);
+      interiorHalfedges.pop_back();
+      if (meshRelation_.triRef[h / 3].faceID >= 0) continue;
+
+      const vec3 v = vertPos_[halfedge_[h].endVert];
+      if (std::abs(dot(v - base, normal)) < tolerance_) {
+        meshRelation_.triRef[h / 3].faceID = tp.tri;
+
+        if (interiorHalfedges.empty() ||
+            h != halfedge_[interiorHalfedges.back()].pairedHalfedge) {
+          interiorHalfedges.push_back(h);
+        } else {
+          interiorHalfedges.pop_back();
+        }
+        const int hNext = NextHalfedge(h);
+        interiorHalfedges.push_back(hNext);
+      }
     }
   }
+}
 
-  for_each_n(autoPolicy(halfedge_.size(), 1e4), countAt(0), NumTri(),
-             CheckCoplanarity(
-                 {comp2tri, halfedge_, vertPos_, &components, tolerance_}));
+void Manifold::Impl::DedupePropVerts() {
+  ZoneScoped;
+  const size_t numProp = NumProp();
+  if (numProp == 0) return;
 
-  Vec<TriRef>& triRef = meshRelation_.triRef;
-  for (size_t tri = 0; tri < NumTri(); ++tri) {
-    const int referenceTri = comp2tri[components[tri]];
-    if (referenceTri >= 0) {
-      triRef[tri].faceID = referenceTri;
-    }
-  }
+  Vec<std::pair<int, int>> vert2vert(halfedge_.size(), {-1, -1});
+  for_each_n(
+      autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
+      [&vert2vert, numProp, this](const int edgeIdx) {
+        const Halfedge edge = halfedge_[edgeIdx];
+        const Halfedge pair = halfedge_[edge.pairedHalfedge];
+        const int edgeFace = edgeIdx / 3;
+        const int pairFace = edge.pairedHalfedge / 3;
+
+        if (meshRelation_.triRef[edgeFace].meshID !=
+            meshRelation_.triRef[pairFace].meshID)
+          return;
+
+        const int baseNum = edgeIdx - 3 * edgeFace;
+        const int jointNum = edge.pairedHalfedge - 3 * pairFace;
+
+        const int prop0 = meshRelation_.triProperties[edgeFace][baseNum];
+        const int prop1 =
+            meshRelation_
+                .triProperties[pairFace][jointNum == 2 ? 0 : jointNum + 1];
+        bool propEqual = true;
+        for (size_t p = 0; p < numProp; ++p) {
+          if (meshRelation_.properties[numProp * prop0 + p] !=
+              meshRelation_.properties[numProp * prop1 + p]) {
+            propEqual = false;
+            break;
+          }
+        }
+        if (propEqual) {
+          vert2vert[edgeIdx] = std::make_pair(prop0, prop1);
+        }
+      });
+
+  std::vector<int> vertLabels;
+  const size_t numPropVert = NumPropVert();
+  const int numLabels = GetLabels(vertLabels, vert2vert, numPropVert);
+
+  std::vector<int> label2vert(numLabels);
+  for (size_t v = 0; v < numPropVert; ++v) label2vert[vertLabels[v]] = v;
+  for (auto& prop : meshRelation_.triProperties)
+    for (int i : {0, 1, 2}) prop[i] = label2vert[vertLabels[prop[i]]];
 }
 
 /**

diff --git a/src/impl.h b/src/impl.h
@@ -227,6 +227,7 @@ struct Manifold::Impl {
       InitializeOriginal();
     }
 
+    DedupePropVerts();
     CreateFaces();
 
     SimplifyTopology();
@@ -267,6 +268,7 @@ struct Manifold::Impl {
   }
 
   void CreateFaces();
+  void DedupePropVerts();
   void RemoveUnreferencedVerts();
   void InitializeOriginal(bool keepFaceID = false);
   void CreateHalfedges(const Vec<ivec3>& triVerts);