Exact Boolean: speed up when there are many separate components.

Use bounding box tests quickly tell that two components cannot
have a containment relation between each other. This change
cut about 0.6s off a test with 25 big icospheres.

3 files changed, 127 insertions, 63 deletions

diff --git a/source/blender/blenlib/BLI_mesh_intersect.hh b/source/blender/blenlib/BLI_mesh_intersect.hh
index 7000349c5af..6b8e79f376c 100644
--- a/source/blender/blenlib/BLI_mesh_intersect.hh
+++ b/source/blender/blenlib/BLI_mesh_intersect.hh
@@ -29,6 +29,7 @@
 
 #  include "BLI_array.hh"
 #  include "BLI_double3.hh"
+#  include "BLI_float3.hh"
 #  include "BLI_index_range.hh"
 #  include "BLI_map.hh"
 #  include "BLI_math_mpq.hh"
@@ -340,6 +341,63 @@ class IMesh {
 std::ostream &operator<<(std::ostream &os, const IMesh &mesh);
 
 /**
+ * A Bounding Box using floats, and a routine to detect possible
+ * intersection.
+ */
+struct BoundingBox {
+  float3 min{FLT_MAX, FLT_MAX, FLT_MAX};
+  float3 max{-FLT_MAX, -FLT_MAX, -FLT_MAX};
+
+  BoundingBox() = default;
+  BoundingBox(const float3 &min, const float3 &max) : min(min), max(max)
+  {
+  }
+
+  void combine(const float3 &p)
+  {
+    min.x = min_ff(min.x, p.x);
+    min.y = min_ff(min.y, p.y);
+    min.z = min_ff(min.z, p.z);
+    max.x = max_ff(max.x, p.x);
+    max.y = max_ff(max.y, p.y);
+    max.z = max_ff(max.z, p.z);
+  }
+
+  void combine(const double3 &p)
+  {
+    min.x = min_ff(min.x, static_cast<float>(p.x));
+    min.y = min_ff(min.y, static_cast<float>(p.y));
+    min.z = min_ff(min.z, static_cast<float>(p.z));
+    max.x = max_ff(max.x, static_cast<float>(p.x));
+    max.y = max_ff(max.y, static_cast<float>(p.y));
+    max.z = max_ff(max.z, static_cast<float>(p.z));
+  }
+
+  void combine(const BoundingBox &bb)
+  {
+    min.x = min_ff(min.x, bb.min.x);
+    min.y = min_ff(min.y, bb.min.y);
+    min.z = min_ff(min.z, bb.min.z);
+    max.x = max_ff(max.x, bb.max.x);
+    max.y = max_ff(max.y, bb.max.y);
+    max.z = max_ff(max.z, bb.max.z);
+  }
+
+  void expand(float pad)
+  {
+    min.x -= pad;
+    min.y -= pad;
+    min.z -= pad;
+    max.x += pad;
+    max.y += pad;
+    max.z += pad;
+  }
+};
+
+/** Assume bounding boxes have been expanded by a sufficient epsilon. */
+bool bbs_might_intersect(const BoundingBox &bb_a, const BoundingBox &bb_b);
+
+/**
  * The output will have duplicate vertices merged and degenerate triangles ignored.
  * If the input has overlapping co-planar triangles, then there will be
  * as many duplicates as there are overlaps in each overlapping triangular region.
diff --git a/source/blender/blenlib/intern/mesh_boolean.cc b/source/blender/blenlib/intern/mesh_boolean.cc
index 00d53a010b0..2b286de5120 100644
--- a/source/blender/blenlib/intern/mesh_boolean.cc
+++ b/source/blender/blenlib/intern/mesh_boolean.cc
@@ -1863,6 +1863,7 @@ static Vector<ComponentContainer> find_component_containers(int comp,
                                                             const IMesh &tm,
                                                             const PatchesInfo &pinfo,
                                                             const TriMeshTopology &tmtopo,
+                                                            Array<BoundingBox> comp_bb,
                                                             IMeshArena *arena)
 {
   constexpr int dbg_level = 0;
@@ -1888,6 +1889,12 @@ static Vector<ComponentContainer> find_component_containers(int comp,
     if (dbg_level > 0) {
       std::cout << "comp_other = " << comp_other << "\n";
     }
+    if (!bbs_might_intersect(comp_bb[comp], comp_bb[comp_other])) {
+      if (dbg_level > 0) {
+        std::cout << "bounding boxes don't overlap\n";
+      }
+      continue;
+    }
     int nearest_tri = NO_INDEX;
     int nearest_tri_close_vert = -1;
     int nearest_tri_close_edge = -1;
@@ -1962,6 +1969,51 @@ static Vector<ComponentContainer> find_component_containers(int comp,
 }
 
 /**
+ * Populate the per-component bounding boxes, expanding them
+ * by an appropriate epsilon so that we conservatively will say
+ * that components could intersect if the BBs overlap.
+ */
+static void populate_comp_bbs(const Vector<Vector<int>> &components,
+                              const PatchesInfo &pinfo,
+                              const IMesh &im,
+                              Array<BoundingBox> &comp_bb)
+{
+  const int comp_grainsize = 16;
+  /* To get a good expansion epsilon, we need to find the maximum
+   * absolute value of any coordinate. Do it first per component,
+   * then get the overall max. */
+  Array<double> max_abs(components.size(), 0.0);
+  parallel_for(components.index_range(), comp_grainsize, [&](IndexRange comp_range) {
+    for (int c : comp_range) {
+      BoundingBox &bb = comp_bb[c];
+      double &maxa = max_abs[c];
+      for (int p : components[c]) {
+        const Patch &patch = pinfo.patch(p);
+        for (int t : patch.tris()) {
+          const Face &tri = *im.face(t);
+          for (const Vert *v : tri) {
+            bb.combine(v->co);
+            for (int i = 0; i < 3; ++i) {
+              maxa = max_dd(maxa, fabs(v->co[i]));
+            }
+          }
+        }
+      }
+    }
+  });
+  double all_max_abs = 0.0;
+  for (int c : components.index_range()) {
+    all_max_abs = max_dd(all_max_abs, max_abs[c]);
+  }
+  constexpr float pad_factor = 10.0f;
+  float pad = all_max_abs == 0.0 ? FLT_EPSILON : 2 * FLT_EPSILON * all_max_abs;
+  pad *= pad_factor;
+  for (int c : components.index_range()) {
+    comp_bb[c].expand(pad);
+  }
+}
+
+/**
  * The cells and patches are supposed to form a bipartite graph.
  * The graph may be disconnected (if parts of meshes are nested or side-by-side
  * without intersection with other each other).
@@ -2003,19 +2055,23 @@ static void finish_patch_cell_graph(const IMesh &tm,
   }
   int tot_components = components.size();
   Array<Vector<ComponentContainer>> comp_cont(tot_components);
-  for (int comp : components.index_range()) {
-    comp_cont[comp] = find_component_containers(
-        comp, components, ambient_cell, tm, pinfo, tmtopo, arena);
-  }
-  if (dbg_level > 0) {
-    std::cout << "component containers:\n";
-    for (int comp : comp_cont.index_range()) {
-      std::cout << comp << ": ";
-      for (const ComponentContainer &cc : comp_cont[comp]) {
-        std::cout << "[containing_comp=" << cc.containing_component
-                  << ", nearest_cell=" << cc.nearest_cell << ", d2=" << cc.dist_to_cell << "] ";
+  if (tot_components > 1) {
+    Array<BoundingBox> comp_bb(tot_components);
+    populate_comp_bbs(components, pinfo, tm, comp_bb);
+    for (int comp : components.index_range()) {
+      comp_cont[comp] = find_component_containers(
+          comp, components, ambient_cell, tm, pinfo, tmtopo, comp_bb, arena);
+    }
+    if (dbg_level > 0) {
+      std::cout << "component containers:\n";
+      for (int comp : comp_cont.index_range()) {
+        std::cout << comp << ": ";
+        for (const ComponentContainer &cc : comp_cont[comp]) {
+          std::cout << "[containing_comp=" << cc.containing_component
+                    << ", nearest_cell=" << cc.nearest_cell << ", d2=" << cc.dist_to_cell << "] ";
+        }
+        std::cout << "\n";
       }
-      std::cout << "\n";
     }
   }
   if (dbg_level > 1) {
diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc
index 97f856476c5..5b7b6d9a929 100644
--- a/source/blender/blenlib/intern/mesh_intersect.cc
+++ b/source/blender/blenlib/intern/mesh_intersect.cc
@@ -744,62 +744,12 @@ std::ostream &operator<<(std::ostream &os, const IMesh &mesh)
   return os;
 }
 
-struct BoundingBox {
-  float3 min{FLT_MAX, FLT_MAX, FLT_MAX};
-  float3 max{-FLT_MAX, -FLT_MAX, -FLT_MAX};
-
-  BoundingBox() = default;
-  BoundingBox(const float3 &min, const float3 &max) : min(min), max(max)
-  {
-  }
-
-  void combine(const float3 &p)
-  {
-    min.x = min_ff(min.x, p.x);
-    min.y = min_ff(min.y, p.y);
-    min.z = min_ff(min.z, p.z);
-    max.x = max_ff(max.x, p.x);
-    max.y = max_ff(max.y, p.y);
-    max.z = max_ff(max.z, p.z);
-  }
-
-  void combine(const double3 &p)
-  {
-    min.x = min_ff(min.x, static_cast<float>(p.x));
-    min.y = min_ff(min.y, static_cast<float>(p.y));
-    min.z = min_ff(min.z, static_cast<float>(p.z));
-    max.x = max_ff(max.x, static_cast<float>(p.x));
-    max.y = max_ff(max.y, static_cast<float>(p.y));
-    max.z = max_ff(max.z, static_cast<float>(p.z));
-  }
-
-  void combine(const BoundingBox &bb)
-  {
-    min.x = min_ff(min.x, bb.min.x);
-    min.y = min_ff(min.y, bb.min.y);
-    min.z = min_ff(min.z, bb.min.z);
-    max.x = max_ff(max.x, bb.max.x);
-    max.y = max_ff(max.y, bb.max.y);
-    max.z = max_ff(max.z, bb.max.z);
-  }
-
-  void expand(float pad)
-  {
-    min.x -= pad;
-    min.y -= pad;
-    min.z -= pad;
-    max.x += pad;
-    max.y += pad;
-    max.z += pad;
-  }
-};
-
 /**
  * Assume bounding boxes have been expanded by a sufficient epsilon on all sides
  * so that the comparisons against the bb bounds are sufficient to guarantee that
  * if an overlap or even touching could happen, this will return true.
  */
-static bool bbs_might_intersect(const BoundingBox &bb_a, const BoundingBox &bb_b)
+bool bbs_might_intersect(const BoundingBox &bb_a, const BoundingBox &bb_b)
 {
   return isect_aabb_aabb_v3(bb_a.min, bb_a.max, bb_b.min, bb_b.max);
 }