optimizing collider

[pca006132]

I'm trying to optimize the collider by adding a finer collision detection at the end to check, e.g., edge-face collision instead of just relying on the bounding box. With a very simple implementation, there seems to be decent speed improvement (30% overall? From 26618ms to 17685ms) because this eliminates most of the false collision pairs, but there is a failing case with non-manifold result:

My patch (mainly for debugging now):

diff --git a/src/collider/include/manifold/collider.h b/src/collider/include/manifold/collider.h
index 37a0cb91..7d2a47ba 100644
--- a/src/collider/include/manifold/collider.h
+++ b/src/collider/include/manifold/collider.h
@@ -154,18 +154,25 @@ struct CreateRadixTree {
   }
 };
 
-template <typename T, const bool selfCollision, typename Recorder>
+struct NaiveRefinement {
+  bool operator()(size_t a, size_t b) { return true; }
+};
+
+template <typename T, const bool selfCollision, typename Recorder,
+          typename Refinement>
 struct FindCollisions {
   VecView<const T> queries;
   VecView<const Box> nodeBBox_;
   VecView<const std::pair<int, int>> internalChildren_;
   Recorder recorder;
+  Refinement refinement;
 
   int RecordCollision(int node, const int queryIdx) {
     bool overlaps = nodeBBox_[node].DoesOverlap(queries[queryIdx]);
     if (overlaps && IsLeaf(node)) {
       int leafIdx = Node2Leaf(node);
-      if (!selfCollision || leafIdx != queryIdx) {
+      if ((!selfCollision || leafIdx != queryIdx) &&
+          refinement(queryIdx, leafIdx)) {
         recorder.record(queryIdx, leafIdx);
       }
     }
@@ -330,45 +337,53 @@ class Collider {
   }
 
   template <const bool selfCollision = false, const bool inverted = false,
-            typename T>
-  SparseIndices Collisions(const VecView<const T>& queriesIn) const {
+            typename T, typename Refinement=collider_internal::NaiveRefinement>
+  SparseIndices Collisions(const VecView<const T>& queriesIn,
+                           Refinement refinement = collider_internal::NaiveRefinement()) const {
     ZoneScoped;
     using collider_internal::FindCollisions;
+    using collider_internal::NaiveRefinement;
     // note that the length is 1 larger than the number of queries so the last
     // element can store the sum when using exclusive scan
-    if (queriesIn.size() < collider_internal::kSequentialThreshold) {
+    // if (queriesIn.size() < collider_internal::kSequentialThreshold) {
       SparseIndices queryTri;
       for_each_n(
           ExecutionPolicy::Seq, countAt(0), queriesIn.size(),
           FindCollisions<T, selfCollision,
-                         collider_internal::SeqCollisionRecorder<inverted>>{
-              queriesIn, nodeBBox_, internalChildren_, {queryTri}});
-      return queryTri;
-    } else {
-      // compute the number of collisions to determine the size for allocation
-      // and offset, this avoids the need for atomic
-      Vec<int> counts(queriesIn.size() + 1, 0);
-      Vec<char> empty(queriesIn.size(), 0);
-      for_each_n(
-          ExecutionPolicy::Par, countAt(0), queriesIn.size(),
-          FindCollisions<T, selfCollision, collider_internal::CountCollisions>{
-              queriesIn, nodeBBox_, internalChildren_, {counts, empty}});
-      // compute start index for each query and total count
-      manifold::exclusive_scan(counts.begin(), counts.end(), counts.begin(), 0,
-                               std::plus<int>());
-      if (counts.back() == 0) return SparseIndices(0);
-      SparseIndices queryTri(counts.back());
-      // actually recording collisions
-      for_each_n(
-          ExecutionPolicy::Par, countAt(0), queriesIn.size(),
-          FindCollisions<T, selfCollision,
-                         collider_internal::ParCollisionRecorder<inverted>>{
-              queriesIn,
-              nodeBBox_,
-              internalChildren_,
-              {queryTri, counts, empty}});
+                         collider_internal::SeqCollisionRecorder<inverted>,
+                         Refinement>{
+              queriesIn, nodeBBox_, internalChildren_, {queryTri}, refinement});
       return queryTri;
-    }
+    // } else {
+    //   // compute the number of collisions to determine the size for allocation
+    //   // and offset, this avoids the need for atomic
+    //   Vec<int> counts(queriesIn.size() + 1, 0);
+    //   Vec<char> empty(queriesIn.size(), 0);
+    //   for_each_n(
+    //       ExecutionPolicy::Par, countAt(0), queriesIn.size(),
+    //       FindCollisions<T, selfCollision, collider_internal::CountCollisions,
+    //                      Refinement>{queriesIn,
+    //                                  nodeBBox_,
+    //                                  internalChildren_,
+    //                                  {counts, empty},
+    //                                  refinement});
+    //   // compute start index for each query and total count
+    //   manifold::exclusive_scan(counts.begin(), counts.end(), counts.begin(), 0,
+    //                            std::plus<int>());
+    //   if (counts.back() == 0) return SparseIndices(0);
+    //   SparseIndices queryTri(counts.back());
+    //   // actually recording collisions
+    //   for_each_n(
+    //       ExecutionPolicy::Par, countAt(0), queriesIn.size(),
+    //       FindCollisions<T, selfCollision,
+    //                      collider_internal::ParCollisionRecorder<inverted>,
+    //                      Refinement>{queriesIn,
+    //                                  nodeBBox_,
+    //                                  internalChildren_,
+    //                                  {queryTri, counts, empty},
+    //                                  refinement});
+    //   return queryTri;
+    // }
   }
 
   static uint32_t MortonCode(vec3 position, Box bBox) {
diff --git a/src/manifold/src/impl.cpp b/src/manifold/src/impl.cpp
index 590e5186..d9596047 100644
--- a/src/manifold/src/impl.cpp
+++ b/src/manifold/src/impl.cpp
@@ -21,6 +21,7 @@
 #include "manifold/hashtable.h"
 #include "manifold/parallel.h"
 #include "manifold/svd.h"
+#include "manifold/tri_dist.h"
 #include "mesh_fixes.h"
 
 namespace {
@@ -45,13 +46,13 @@ struct Transform4x3 {
   vec3 operator()(vec3 position) { return transform * vec4(position, 1.0); }
 };
 
+template <bool calculateTriNormal>
 struct AssignNormals {
   VecView<vec3> faceNormal;
   VecView<vec3> vertNormal;
   VecView<const vec3> vertPos;
   VecView<const Halfedge> halfedges;
   const double precision;
-  const bool calculateTriNormal;
 
   void operator()(const int face) {
     vec3& triNormal = faceNormal[face];
@@ -657,9 +658,14 @@ void Manifold::Impl::CalculateNormals() {
     faceNormal_.resize(NumTri());
     calculateTriNormal = true;
   }
-  for_each_n(policy, countAt(0), NumTri(),
-             AssignNormals({faceNormal_, vertNormal_, vertPos_, halfedge_,
-                            precision_, calculateTriNormal}));
+  if (calculateTriNormal)
+    for_each_n(policy, countAt(0), NumTri(),
+               AssignNormals<true>({faceNormal_, vertNormal_, vertPos_,
+                                    halfedge_, precision_}));
+  else
+    for_each_n(policy, countAt(0), NumTri(),
+               AssignNormals<false>({faceNormal_, vertNormal_, vertPos_,
+                                     halfedge_, precision_}));
   for_each(policy, vertNormal_.begin(), vertNormal_.end(),
            [](vec3& v) { v = SafeNormalize(v); });
 }
@@ -703,11 +709,35 @@ SparseIndices Manifold::Impl::EdgeCollisions(const Impl& Q,
         QedgeBB[e] = Box(vertPos[edges[e].first], vertPos[edges[e].second]);
       });
 
+  auto newPrecision = std::max(Q.precision_, precision_);
   SparseIndices q1p2(0);
+  std::cout << std::setprecision(15);
+  std::cout << "precison: " << newPrecision << std::endl;
+  auto edgeCollisionRefinement = [&](size_t edge, size_t face) {
+    std::array<vec3, 3> tri = {vertPos_[halfedge_[3 * face].startVert],
+                               vertPos_[halfedge_[3 * face + 1].startVert],
+                               vertPos_[halfedge_[3 * face + 2].startVert]};
+    double distance = DistanceTriangleTriangleSquared(
+        {vertPos[edges[edge].first], vertPos[edges[edge].second],
+         vertPos[edges[edge].second]},
+        tri);
+
+    if (std::sqrt(distance) <= 1e3 * newPrecision &&
+        std::sqrt(distance) > 2 * precision_)
+      std::cout << std::sqrt(distance) << " "
+                << glm::length(vertPos[edges[edge].first] -
+                               vertPos[edges[edge].second])
+                << " " << glm::length(tri[1] - tri[0]) << " "
+                << glm::length(tri[2] - tri[1]) << " "
+                << glm::length(tri[0] - tri[2]) << std::endl;
+    return std::sqrt(distance) <= 2 * newPrecision ;
+  };
   if (inverted)
-    q1p2 = collider_.Collisions<false, true>(QedgeBB.cview());
+    q1p2 = collider_.Collisions<false, true>(QedgeBB.cview(),
+                                             edgeCollisionRefinement);
   else
-    q1p2 = collider_.Collisions<false, false>(QedgeBB.cview());
+    q1p2 = collider_.Collisions<false, false>(QedgeBB.cview(),
+                                              edgeCollisionRefinement);
 
   if (inverted)
     for_each(policy, countAt(0_uz), countAt(q1p2.size()),

Error case:

[ RUN      ] Hull.Fail
precison: 1.87068878173828e-10
1.21131336592732e-07 24.0340458178366 2.91621236208739 8.72515856104347 5.82655334387241
1.60093420138062e-07 98.4000055854313 2 2.37783235835963 3.10710262535084
1.27515246202283e-07 35.0823071951458 2.91621304233482 52.5861043176949 54.2472584862783
5.7979293203925e-08 35.0823071951458 57.2387040806607 61.8186622666715 11.6061857296013
5.7979293203925e-08 35.0823071951458 5.82655211654379 57.2387040806607 54.2472584862783
6.01224874383682e-08 64.038470103086 8.7251509509721 20.0860930597135 11.6061996531557
6.01224874383682e-08 64.038470103086 22.8393674546187 40.5268325638132 20.0860930597135
8.30666705269674e-08 68.4000182061008 5.82655211654379 57.2387040806607 54.2472584862783
4.59502335852199e-08 68.3311380429698 11.6061972164029 22.8393674546187 11.6062039995614
4.59502335852199e-08 68.3311380429698 5.82655984703431 11.6061972164029 5.82655039677679
1.58178382649028e-07 37.6813484766753 64.7383344623228 40.5268527590104 46.9816724057741
precison: 1.87068878173828e-10
1.07584822699139e-07 14.4638554785791 6.13133058510376 26.8829098138667 32.4126130511599
2.10663738739981e-08 8.72516286239364 12.2132860069237 12.2133057474357 24.0340458178366
1.60093420239025e-07 3.10710262535084 27.8207861013618 90.7172594241457 98.4000055854313
1.89429547827839e-08 20.0860930597135 35.0823071951458 12.2132921676034 45.0030041820258
1.89429547827839e-08 20.0860930597135 32.4126225820928 45.0030041820258 64.038470103086
1.45857789497073e-07 61.8186622666715 65.9350614266107 68.3311380429698 15.220430808056
5.7979293203925e-08 57.2387040806607 35.0823071951458 12.2132921676034 45.0030041820258
6.79503788179453e-08 57.2387040806607 32.4126225820928 45.0030041820258 64.038470103086
1.45857789497073e-07 57.2387040806607 65.9350614266107 68.3311380429698 15.220430808056
1.45857789497073e-07 5.82655211654379 65.9350614266107 68.3311380429698 15.220430808056
1.27515246202283e-07 54.2472584862783 35.0823071951458 12.2132921676034 45.0030041820258
1.45857789497073e-07 30.0000121370557 65.9350614266107 68.3311380429698 15.220430808056
1.45857789497073e-07 17.2923873012639 65.9350614266107 68.3311380429698 15.220430808056
1.45857789497073e-07 63.9554176289687 65.9350614266107 68.3311380429698 15.220430808056
1.45857789497073e-07 2.91621643495683 65.9350614266107 68.3311380429698 15.220430808056
1.45857789497073e-07 2.91621506355526 65.9350614266107 68.3311380429698 15.220430808056
1.45857789497073e-07 5.82654760265774 65.9350614266107 68.3311380429698 15.220430808056
1.45857789497073e-07 8.72515244254667 65.9350614266107 68.3311380429698 15.220430808056
unknown file: Failure
C++ exception with description "Error in file: /home/pca006132/code/manifold/src/manifold/src/boolean_result.cpp (278): 'edgePos.size() % 2 == 0' is false: Non-manifold edge! Not an even number of points." thrown in the test body.

And I verified that there are some edge-face pairs printed that has to go into the collision result, but I am not sure why their distance is larger than 2*precisioon.

@elalish any idea?

[pca006132]

Btw, the changes to CalculateNormals is another optimization that improves the performance of itself by 10%, which may be helpful when importing large meshes.

[elalish]

Hmm, any idea how much bigger? One minor thing: I generally convert sqrt(d) > p into d > p * p, since it's both faster and has less rounding error (probably). Usually the issue is with inverting nearly singular matrices, since that explodes the rounding error, but that can come in more subtle ways too. If the only error is in Hull, I'd guess something is creating more rounding error than it should.

[elalish]

This looks both a bit sketchy and inefficient:

double distance = DistanceTriangleTriangleSquared(
        {vertPos[edges[edge].first], vertPos[edges[edge].second],
         vertPos[edges[edge].second]},
        tri);

These are exactly the kind of geometric checks that so frequently go wrong. I would guess we'll have to carefully rewrite this function to guarantee bounded floating-point error.

[pca006132]

Sure, I will open a PR for this. I just use this as this is apparently the simplest implementation. I was also thinking about floating-point error, but it seems unlikely to me that the error would be this large.

[pca006132]

I was thinking more about efficiency, because this function is indeed taking a large amount of time in the profiling result (and function calls in Boolean3 are now taking a lot less time!), but I feel that line-triangle distance is not something I can write in several minutes, so I just keep using this triangle-triangle distance :)

[elalish]

I bet if we look inside that tri-tri distance it's really just three tri-edge distances, or something close. Could probably use some cleanup in any case.

[pca006132]

Sure. Anyway, based on the preliminary results, I think it will be worth it to put some effort into this optimization.