BVH: implement an optimized self-overlap traversal.

Previously cloth self-collision and other code that wants self-overlap
data was using an ordinary BVH overlap utility function. This works,
but is suboptimal because it fails to take advantage of the fact that
it is comparing two identical trees.

The standard traversal for self-collision essentially devolves into
enumerating all elements of the tree, and then matching them to the
tree starting at the root. However, the tree itself naturally partitions
the space, so overlapping elements are likely to be mostly placed nearby
in the tree structure as well.

Instead, self-overlap can be computed by recursively computing ordinary
overlap between siblings within each subtree. This only considers each
pair of leaves once, and provides significantly better locality.
It also avoids outputting every overlap pair twice in different order.

Differential Revision: https://developer.blender.org/D16915

source/blender/blenkernel/intern/collision.c

@@ -1521,8 +1521,9 @@ static bool cloth_bvh_obj_overlap_cb(void *userdata,
 
 static bool cloth_bvh_self_overlap_cb(void *userdata, int index_a, int index_b, int UNUSED(thread))
 {
-  /* No need for equal combinations (eg. (0,1) & (1,0)). */
-  if (index_a < index_b) {
+  /* This shouldn't happen, but just in case. Note that equal combinations
+   * (eg. (0,1) & (1,0)) would be filtered out by BLI_bvhtree_overlap_self. */
+  if (index_a != index_b) {
     ClothModifierData *clmd = (ClothModifierData *)userdata;
     struct Cloth *clothObject = clmd->clothObject;
     const MVertTri *tri_a, *tri_b;
@@ -1599,8 +1600,8 @@ int cloth_bvh_collision(
   if (clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_SELF) {
     bvhtree_update_from_cloth(clmd, false, true);
 
-    overlap_self = BLI_bvhtree_overlap(
-        cloth->bvhselftree, cloth->bvhselftree, &coll_count_self, cloth_bvh_self_overlap_cb, clmd);
+    overlap_self = BLI_bvhtree_overlap_self(
+        cloth->bvhselftree, &coll_count_self, cloth_bvh_self_overlap_cb, clmd);
   }
 
   do {

source/blender/blenlib/BLI_kdopbvh.h

@@ -73,9 +73,11 @@ typedef struct BVHTreeRayHit {
 } BVHTreeRayHit;
 
 enum {
-  /* Use a priority queue to process nodes in the optimal order (for slow callbacks) */
   BVH_OVERLAP_USE_THREADING = (1 << 0),
   BVH_OVERLAP_RETURN_PAIRS = (1 << 1),
+  /* Use a specialized self-overlap traversal to only test and output every
+   * pair once, rather than twice in different order as usual. */
+  BVH_OVERLAP_SELF = (1 << 2),
 };
 enum {
   /* Use a priority queue to process nodes in the optimal order (for slow callbacks) */
@@ -163,6 +165,9 @@ bool BLI_bvhtree_update_node(
     BVHTree *tree, int index, const float co[3], const float co_moving[3], int numpoints);
 /**
  * Call #BLI_bvhtree_update_node() first for every node/point/triangle.
+ *
+ * Note that this does not rebalance the tree, so if the shape of the mesh changes
+ * too much, operations on the tree may become suboptimal.
  */
 void BLI_bvhtree_update_tree(BVHTree *tree);
 
@@ -191,6 +196,12 @@ BVHTreeOverlap *BLI_bvhtree_overlap(const BVHTree *tree1,
                                     unsigned int *r_overlap_num,
                                     BVHTree_OverlapCallback callback,
                                     void *userdata);
+/** Compute overlaps of the tree with itself. This is faster than BLI_bvhtree_overlap
+ *  because it only tests and returns each symmetrical pair once. */
+BVHTreeOverlap *BLI_bvhtree_overlap_self(const BVHTree *tree,
+                                         unsigned int *r_overlap_num,
+                                         BVHTree_OverlapCallback callback,
+                                         void *userdata);
 
 int *BLI_bvhtree_intersect_plane(BVHTree *tree, float plane[4], uint *r_intersect_num);
 

source/blender/blenlib/intern/BLI_kdopbvh.c

@@ -95,6 +95,7 @@ BLI_STATIC_ASSERT((sizeof(void *) == 8 && sizeof(BVHTree) <= 48) ||
 typedef struct BVHOverlapData_Shared {
   const BVHTree *tree1, *tree2;
   axis_t start_axis, stop_axis;
+  bool use_self;
 
   /* use for callbacks */
   BVHTree_OverlapCallback callback;
@@ -1231,6 +1232,62 @@ static bool tree_overlap_traverse_num(BVHOverlapData_Thread *data_thread,
   return false;
 }
 
+/** Calls the appropriate recursive overlap traversal function. */
+static void tree_overlap_invoke_traverse(BVHOverlapData_Thread *data,
+                                         const BVHNode *node1,
+                                         const BVHNode *node2)
+{
+  if (data->max_interactions) {
+    tree_overlap_traverse_num(data, node1, node2);
+  }
+  else if (data->shared->callback) {
+    tree_overlap_traverse_cb(data, node1, node2);
+  }
+  else {
+    tree_overlap_traverse(data, node1, node2);
+  }
+}
+
+/** Self-overlap traversal with callback. */
+static void tree_overlap_traverse_self_cb(BVHOverlapData_Thread *data_thread, const BVHNode *node)
+{
+  for (int i = 0; i < node->node_num; i++) {
+    /* Recursively compute self-overlap within each child. */
+    tree_overlap_traverse_self_cb(data_thread, node->children[i]);
+
+    /* Compute overlap of pairs of children, testing each one only once (assume symmetry). */
+    for (int j = i + 1; j < node->node_num; j++) {
+      tree_overlap_traverse_cb(data_thread, node->children[i], node->children[j]);
+    }
+  }
+}
+
+/** Self-overlap traversal without callback. */
+static void tree_overlap_traverse_self(BVHOverlapData_Thread *data_thread, const BVHNode *node)
+{
+  for (int i = 0; i < node->node_num; i++) {
+    /* Recursively compute self-overlap within each child. */
+    tree_overlap_traverse_self(data_thread, node->children[i]);
+
+    /* Compute overlap of pairs of children, testing each one only once (assume symmetry). */
+    for (int j = i + 1; j < node->node_num; j++) {
+      tree_overlap_traverse(data_thread, node->children[i], node->children[j]);
+    }
+  }
+}
+
+/** Calls the appropriate recursive self-overlap traversal. */
+static void tree_overlap_invoke_traverse_self(BVHOverlapData_Thread *data_thread,
+                                              const BVHNode *node)
+{
+  if (data_thread->shared->callback) {
+    tree_overlap_traverse_self_cb(data_thread, node);
+  }
+  else {
+    tree_overlap_traverse_self(data_thread, node);
+  }
+}
+
 int BLI_bvhtree_overlap_thread_num(const BVHTree *tree)
 {
   return (int)MIN2(tree->tree_type, tree->nodes[tree->leaf_num]->node_num);
@@ -1243,20 +1300,20 @@ static void bvhtree_overlap_task_cb(void *__restrict userdata,
   BVHOverlapData_Thread *data = &((BVHOverlapData_Thread *)userdata)[j];
   BVHOverlapData_Shared *data_shared = data->shared;
 
-  if (data->max_interactions) {
-    tree_overlap_traverse_num(data,
-                              data_shared->tree1->nodes[data_shared->tree1->leaf_num]->children[j],
-                              data_shared->tree2->nodes[data_shared->tree2->leaf_num]);
-  }
-  else if (data_shared->callback) {
-    tree_overlap_traverse_cb(data,
-                             data_shared->tree1->nodes[data_shared->tree1->leaf_num]->children[j],
-                             data_shared->tree2->nodes[data_shared->tree2->leaf_num]);
+  const BVHNode *root1 = data_shared->tree1->nodes[data_shared->tree1->leaf_num];
+
+  if (data_shared->use_self) {
+    /* This code matches one outer loop iteration within traverse_self. */
+    tree_overlap_invoke_traverse_self(data, root1->children[j]);
+
+    for (int k = j + 1; k < root1->node_num; k++) {
+      tree_overlap_invoke_traverse(data, root1->children[j], root1->children[k]);
+    }
   }
   else {
-    tree_overlap_traverse(data,
-                          data_shared->tree1->nodes[data_shared->tree1->leaf_num]->children[j],
-                          data_shared->tree2->nodes[data_shared->tree2->leaf_num]);
+    const BVHNode *root2 = data_shared->tree2->nodes[data_shared->tree2->leaf_num];
+
+    tree_overlap_invoke_traverse(data, root1->children[j], root2);
   }
 }
 
@@ -1273,9 +1330,12 @@ BVHTreeOverlap *BLI_bvhtree_overlap_ex(
   bool overlap_pairs = (flag & BVH_OVERLAP_RETURN_PAIRS) != 0;
   bool use_threading = (flag & BVH_OVERLAP_USE_THREADING) != 0 &&
                        (tree1->leaf_num > KDOPBVH_THREAD_LEAF_THRESHOLD);
+  bool use_self = (flag & BVH_OVERLAP_SELF) != 0;
 
   /* 'RETURN_PAIRS' was not implemented without 'max_interactions'. */
   BLI_assert(overlap_pairs || max_interactions);
+  /* Self-overlap does not support max interactions (it's not symmetrical). */
+  BLI_assert(!use_self || tree1 == tree2 && !max_interactions);
 
   const int root_node_len = BLI_bvhtree_overlap_thread_num(tree1);
   const int thread_num = use_threading ? root_node_len : 1;
@@ -1293,6 +1353,10 @@ BVHTreeOverlap *BLI_bvhtree_overlap_ex(
     return NULL;
   }
 
+  if (UNLIKELY(use_self && tree1 != tree2)) {
+    use_self = false;
+  }
+
   const BVHNode *root1 = tree1->nodes[tree1->leaf_num];
   const BVHNode *root2 = tree2->nodes[tree2->leaf_num];
 
@@ -1308,6 +1372,7 @@ BVHTreeOverlap *BLI_bvhtree_overlap_ex(
   data_shared.tree2 = tree2;
   data_shared.start_axis = start_axis;
   data_shared.stop_axis = stop_axis;
+  data_shared.use_self = use_self;
 
   /* can be NULL */
   data_shared.callback = callback;
@@ -1317,7 +1382,7 @@ BVHTreeOverlap *BLI_bvhtree_overlap_ex(
     /* init BVHOverlapData_Thread */
     data[j].shared = &data_shared;
     data[j].overlap = overlap_pairs ? BLI_stack_new(sizeof(BVHTreeOverlap), __func__) : NULL;
-    data[j].max_interactions = max_interactions;
+    data[j].max_interactions = use_self ? 0 : max_interactions;
 
     /* for callback */
     data[j].thread = j;
@@ -1329,16 +1394,11 @@ BVHTreeOverlap *BLI_bvhtree_overlap_ex(
     settings.min_iter_per_thread = 1;
     BLI_task_parallel_range(0, root_node_len, data, bvhtree_overlap_task_cb, &settings);
   }
+  else if (use_self) {
+    tree_overlap_invoke_traverse_self(data, root1);
+  }
   else {
-    if (max_interactions) {
-      tree_overlap_traverse_num(data, root1, root2);
-    }
-    else if (callback) {
-      tree_overlap_traverse_cb(data, root1, root2);
-    }
-    else {
-      tree_overlap_traverse(data, root1, root2);
-    }
+    tree_overlap_invoke_traverse(data, root1, root2);
   }
 
   if (overlap_pairs) {
@@ -1377,6 +1437,23 @@ BVHTreeOverlap *BLI_bvhtree_overlap(
                                 BVH_OVERLAP_USE_THREADING | BVH_OVERLAP_RETURN_PAIRS);
 }
 
+BVHTreeOverlap *BLI_bvhtree_overlap_self(
+    const BVHTree *tree,
+    uint *r_overlap_num,
+    /* optional callback to test the overlap before adding (must be thread-safe!) */
+    BVHTree_OverlapCallback callback,
+    void *userdata)
+{
+  return BLI_bvhtree_overlap_ex(tree,
+                                tree,
+                                r_overlap_num,
+                                callback,
+                                userdata,
+                                0,
+                                BVH_OVERLAP_USE_THREADING | BVH_OVERLAP_RETURN_PAIRS |
+                                    BVH_OVERLAP_SELF);
+}
+
 /** \} */
 
 /* -------------------------------------------------------------------- */

source/blender/draw/intern/mesh_extractors/extract_mesh_vbo_mesh_analysis.cc

@@ -347,7 +347,7 @@ static void statvis_calc_intersect(const MeshRenderData *mr, float *r_intersect)
     data.mlooptri = mr->mlooptri;
     data.epsilon = BLI_bvhtree_get_epsilon(tree);
 
-    BVHTreeOverlap *overlap = BLI_bvhtree_overlap(tree, tree, &overlap_len, bvh_overlap_cb, &data);
+    BVHTreeOverlap *overlap = BLI_bvhtree_overlap_self(tree, &overlap_len, bvh_overlap_cb, &data);
     if (overlap) {
       for (int i = 0; i < overlap_len; i++) {
         const MPoly *f_hit_pair[2] = {

source/blender/editors/uvedit/uvedit_select.c

@@ -4427,7 +4427,7 @@ static int uv_select_overlap(bContext *C, const bool extend)
   BLI_bvhtree_balance(uv_tree);
 
   uint tree_overlap_len;
-  BVHTreeOverlap *overlap = BLI_bvhtree_overlap(uv_tree, uv_tree, &tree_overlap_len, NULL, NULL);
+  BVHTreeOverlap *overlap = BLI_bvhtree_overlap_self(uv_tree, &tree_overlap_len, NULL, NULL);
 
   if (overlap != NULL) {
     GSet *overlap_set = BLI_gset_new_ex(overlap_hash, overlap_cmp, __func__, tree_overlap_len);