Cloth: Optimization in self collision

15% to 20% improvement in cloth simulation performance with
self-collision.

The idea is to reduce the number of collisions computed by avoiding
overlapping tris with the same combination (eg. (1,0) and (0,1)).

Reviewed By: zeddb

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

1 files changed, 119 insertions, 66 deletions

diff --git a/source/blender/blenkernel/intern/collision.c b/source/blender/blenkernel/intern/collision.c
index 5db42618a9e..5d4753e7ef1 100644
--- a/source/blender/blenkernel/intern/collision.c
+++ b/source/blender/blenkernel/intern/collision.c
@@ -598,8 +598,13 @@ static float compute_collision_point_edge_tri(const float a1[3],
 }
 
 // w3 is not perfect
-static void collision_compute_barycentric(
-    const float pv[3], float p1[3], float p2[3], float p3[3], float *w1, float *w2, float *w3)
+static void collision_compute_barycentric(const float pv[3],
+                                          const float p1[3],
+                                          const float p2[3],
+                                          const float p3[3],
+                                          float *w1,
+                                          float *w2,
+                                          float *w3)
 {
   /* dot_v3v3 */
 #define INPR(v1, v2) ((v1)[0] * (v2)[0] + (v1)[1] * (v2)[1] + (v1)[2] * (v2)[2])
@@ -646,8 +651,13 @@ static void collision_compute_barycentric(
 #  pragma GCC diagnostic ignored "-Wdouble-promotion"
 #endif
 
-DO_INLINE void collision_interpolateOnTriangle(
-    float to[3], float v1[3], float v2[3], float v3[3], double w1, double w2, double w3)
+DO_INLINE void collision_interpolateOnTriangle(float to[3],
+                                               const float v1[3],
+                                               const float v2[3],
+                                               const float v3[3],
+                                               const double w1,
+                                               const double w2,
+                                               const double w3)
 {
   zero_v3(to);
   VECADDMUL(to, v1, w1);
@@ -867,6 +877,31 @@ static int cloth_collision_response_static(ClothModifierData *clmd,
   return result;
 }
 
+static void cloth_selfcollision_impulse_vert(const float clamp_sq,
+                                             const float impulse[3],
+                                             struct ClothVertex *vert)
+{
+  float impulse_len_sq = len_squared_v3(impulse);
+
+  if ((clamp_sq > 0.0f) && (impulse_len_sq > clamp_sq)) {
+    return;
+  }
+
+  if (ABS(vert->impulse[0]) < ABS(impulse[0])) {
+    vert->impulse[0] = impulse[0];
+  }
+
+  if (ABS(vert->impulse[1]) < ABS(impulse[1])) {
+    vert->impulse[1] = impulse[1];
+  }
+
+  if (ABS(vert->impulse[2]) < ABS(impulse[2])) {
+    vert->impulse[2] = impulse[2];
+  }
+
+  vert->impulse_count++;
+}
+
 static int cloth_selfcollision_response_static(ClothModifierData *clmd,
                                                CollPair *collpair,
                                                uint collision_count,
@@ -881,11 +916,8 @@ static int cloth_selfcollision_response_static(ClothModifierData *clmd,
   cloth1 = clmd->clothObject;
 
   for (int i = 0; i < collision_count; i++, collpair++) {
-    float i1[3], i2[3], i3[3];
-
-    zero_v3(i1);
-    zero_v3(i2);
-    zero_v3(i3);
+    float ia[3][3] = {{0.0f}};
+    float ib[3][3] = {{0.0f}};
 
     /* Only handle static collisions here. */
     if (collpair->flag & (COLLISION_IN_FUTURE | COLLISION_INACTIVE)) {
@@ -959,22 +991,25 @@ static int cloth_selfcollision_response_static(ClothModifierData *clmd,
 
         impulse = magtangent / 1.5;
 
-        VECADDMUL(i1, vrel_t_pre, w1 * impulse);
-        VECADDMUL(i2, vrel_t_pre, w2 * impulse);
-        VECADDMUL(i3, vrel_t_pre, w3 * impulse);
+        VECADDMUL(ia[0], vrel_t_pre, w1 * impulse);
+        VECADDMUL(ia[1], vrel_t_pre, w2 * impulse);
+        VECADDMUL(ia[2], vrel_t_pre, w3 * impulse);
+
+        VECADDMUL(ib[0], vrel_t_pre, -u1 * impulse);
+        VECADDMUL(ib[1], vrel_t_pre, -u2 * impulse);
+        VECADDMUL(ib[2], vrel_t_pre, -u3 * impulse);
       }
 
       /* Apply velocity stopping impulse. */
       impulse = magrelVel / 3.0f;
 
-      VECADDMUL(i1, collpair->normal, w1 * impulse);
-      cloth1->verts[collpair->ap1].impulse_count++;
-
-      VECADDMUL(i2, collpair->normal, w2 * impulse);
-      cloth1->verts[collpair->ap2].impulse_count++;
+      VECADDMUL(ia[0], collpair->normal, w1 * impulse);
+      VECADDMUL(ia[1], collpair->normal, w2 * impulse);
+      VECADDMUL(ia[2], collpair->normal, w3 * impulse);
 
-      VECADDMUL(i3, collpair->normal, w3 * impulse);
-      cloth1->verts[collpair->ap3].impulse_count++;
+      VECADDMUL(ib[0], collpair->normal, -u1 * impulse);
+      VECADDMUL(ib[1], collpair->normal, -u2 * impulse);
+      VECADDMUL(ib[2], collpair->normal, -u3 * impulse);
 
       time_multiplier = 1.0f / (clmd->sim_parms->dt * clmd->sim_parms->timescale);
 
@@ -991,9 +1026,13 @@ static int cloth_selfcollision_response_static(ClothModifierData *clmd,
 
         impulse = repulse / 1.5f;
 
-        VECADDMUL(i1, collpair->normal, w1 * impulse);
-        VECADDMUL(i2, collpair->normal, w2 * impulse);
-        VECADDMUL(i3, collpair->normal, w3 * impulse);
+        VECADDMUL(ia[0], collpair->normal, w1 * impulse);
+        VECADDMUL(ia[1], collpair->normal, w2 * impulse);
+        VECADDMUL(ia[2], collpair->normal, w3 * impulse);
+
+        VECADDMUL(ib[0], collpair->normal, -u1 * impulse);
+        VECADDMUL(ib[1], collpair->normal, -u2 * impulse);
+        VECADDMUL(ib[2], collpair->normal, -u3 * impulse);
       }
 
       result = 1;
@@ -1009,42 +1048,29 @@ static int cloth_selfcollision_response_static(ClothModifierData *clmd,
         float repulse = d * 1.0f / time_multiplier;
         float impulse = repulse / 9.0f;
 
-        VECADDMUL(i1, collpair->normal, w1 * impulse);
-        VECADDMUL(i2, collpair->normal, w2 * impulse);
-        VECADDMUL(i3, collpair->normal, w3 * impulse);
+        VECADDMUL(ia[0], collpair->normal, w1 * impulse);
+        VECADDMUL(ia[1], collpair->normal, w2 * impulse);
+        VECADDMUL(ia[2], collpair->normal, w3 * impulse);
 
-        cloth1->verts[collpair->ap1].impulse_count++;
-        cloth1->verts[collpair->ap2].impulse_count++;
-        cloth1->verts[collpair->ap3].impulse_count++;
+        VECADDMUL(ib[0], collpair->normal, -u1 * impulse);
+        VECADDMUL(ib[1], collpair->normal, -u2 * impulse);
+        VECADDMUL(ib[2], collpair->normal, -u3 * impulse);
 
         result = 1;
       }
     }
 
     if (result) {
-      float clamp = clmd->coll_parms->self_clamp * dt;
-
-      if ((clamp > 0.0f) &&
-          ((len_v3(i1) > clamp) || (len_v3(i2) > clamp) || (len_v3(i3) > clamp))) {
-        return 0;
-      }
+      float clamp_sq = clmd->coll_parms->self_clamp * dt;
+      clamp_sq *= clamp_sq;
 
-      for (int j = 0; j < 3; j++) {
-        if (cloth1->verts[collpair->ap1].impulse_count > 0 &&
-            ABS(cloth1->verts[collpair->ap1].impulse[j]) < ABS(i1[j])) {
-          cloth1->verts[collpair->ap1].impulse[j] = i1[j];
-        }
+      cloth_selfcollision_impulse_vert(clamp_sq, ia[0], &cloth1->verts[collpair->ap1]);
+      cloth_selfcollision_impulse_vert(clamp_sq, ia[1], &cloth1->verts[collpair->ap2]);
+      cloth_selfcollision_impulse_vert(clamp_sq, ia[2], &cloth1->verts[collpair->ap3]);
 
-        if (cloth1->verts[collpair->ap2].impulse_count > 0 &&
-            ABS(cloth1->verts[collpair->ap2].impulse[j]) < ABS(i2[j])) {
-          cloth1->verts[collpair->ap2].impulse[j] = i2[j];
-        }
-
-        if (cloth1->verts[collpair->ap3].impulse_count > 0 &&
-            ABS(cloth1->verts[collpair->ap3].impulse[j]) < ABS(i3[j])) {
-          cloth1->verts[collpair->ap3].impulse[j] = i3[j];
-        }
-      }
+      cloth_selfcollision_impulse_vert(clamp_sq, ib[0], &cloth1->verts[collpair->bp1]);
+      cloth_selfcollision_impulse_vert(clamp_sq, ib[1], &cloth1->verts[collpair->bp2]);
+      cloth_selfcollision_impulse_vert(clamp_sq, ib[2], &cloth1->verts[collpair->bp3]);
     }
   }
 
@@ -1112,6 +1138,28 @@ static void cloth_collision(void *__restrict userdata,
   }
 }
 
+static bool cloth_bvh_selfcollision_is_active(const ClothVertex *verts,
+                                              const MVertTri *tri_a,
+                                              const MVertTri *tri_b)
+{
+  /* Ignore overlap of neighboring triangles. */
+  for (uint i = 0; i < 3; i++) {
+    for (uint j = 0; j < 3; j++) {
+      if (tri_a->tri[i] == tri_b->tri[j]) {
+        return false;
+      }
+    }
+  }
+
+  if (((verts[tri_a->tri[0]].flags & verts[tri_a->tri[1]].flags & verts[tri_a->tri[2]].flags) |
+       (verts[tri_b->tri[0]].flags & verts[tri_b->tri[1]].flags & verts[tri_b->tri[2]].flags)) &
+      CLOTH_VERT_FLAG_NOSELFCOLL) {
+    return false;
+  }
+
+  return true;
+}
+
 static void cloth_selfcollision(void *__restrict userdata,
                                 const int index,
                                 const TaskParallelTLS *__restrict UNUSED(tls))
@@ -1129,21 +1177,7 @@ static void cloth_selfcollision(void *__restrict userdata,
   tri_a = &clmd->clothObject->tri[data->overlap[index].indexA];
   tri_b = &clmd->clothObject->tri[data->overlap[index].indexB];
 
-  for (uint i = 0; i < 3; i++) {
-    for (uint j = 0; j < 3; j++) {
-      if (tri_a->tri[i] == tri_b->tri[j]) {
-        collpair[index].flag = COLLISION_INACTIVE;
-        return;
-      }
-    }
-  }
-
-  if (((verts1[tri_a->tri[0]].flags & verts1[tri_a->tri[1]].flags & verts1[tri_a->tri[2]].flags) |
-       (verts1[tri_b->tri[0]].flags & verts1[tri_b->tri[1]].flags & verts1[tri_b->tri[2]].flags)) &
-      CLOTH_VERT_FLAG_NOSELFCOLL) {
-    collpair[index].flag = COLLISION_INACTIVE;
-    return;
-  }
+  BLI_assert(cloth_bvh_selfcollision_is_active(verts1, tri_a, tri_b));
 
   /* Compute distance and normal. */
   distance = compute_collision_point_tri_tri(verts1[tri_a->tri[0]].tx,
@@ -1542,6 +1576,22 @@ static int cloth_bvh_selfcollisions_resolve(ClothModifierData *clmd,
   return ret;
 }
 
+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) {
+    struct Cloth *clothObject = userdata;
+    const MVertTri *tri_a, *tri_b;
+    tri_a = &clothObject->tri[index_a];
+    tri_b = &clothObject->tri[index_b];
+
+    if (cloth_bvh_selfcollision_is_active(clothObject->verts, tri_a, tri_b)) {
+      return true;
+    }
+  }
+  return false;
+}
+
 int cloth_bvh_collision(
     Depsgraph *depsgraph, Object *ob, ClothModifierData *clmd, float step, float dt)
 {
@@ -1602,8 +1652,11 @@ 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, NULL, NULL);
+    overlap_self = BLI_bvhtree_overlap(cloth->bvhselftree,
+                                       cloth->bvhselftree,
+                                       &coll_count_self,
+                                       cloth_bvh_self_overlap_cb,
+                                       clmd->clothObject);
   }
 
   do {