Circumvent bug in mjc_BoxBox.

The box-box collider sometimes outputs bad contact points which are far outside one of the two boxes. This happens at configurations which are very rare numerically, but rare configurations (e.g. perfect box alignment) are often induced by the designer and the physics, making them actually not-so-rare.

Since the box collision functions are very difficult to understand, we circumvent the issue by extending the post processing code to detect and remove such contacts.

PiperOrigin-RevId: 589178197
Change-Id: I82b8a888ef04f9196111c960757e3314d0ab6d8e

src/engine/engine_collision_driver.c

@@ -1444,6 +1444,9 @@ static void mj_makeCapsule(const mjModel* m, mjData* d, int f, const int vid[2],
 
 // test two geoms for collision, apply filters, add to contact list
 void mj_collideGeoms(const mjModel* m, mjData* d, int g1, int g2) {
+  // relative distance (1%) outside of which box-box contacts are removed
+  static mjtNum kBoxRemoveMargin = 1.01;
+
   TM_START;
 
   int num, type1, type2, condim;
@@ -1524,16 +1527,44 @@ void mj_collideGeoms(const mjModel* m, mjData* d, int g1, int g2) {
     mjERROR("too many contacts returned by collision function");
   }
 
-  // remove repeated contacts in box-box
+  // remove bad and repeated contacts in box-box
   if (type1 == mjGEOM_BOX && type2 == mjGEOM_BOX) {
     // use dim field to mark: -1: bad, 0: good
     for (int i=0; i < num; i++) {
       con[i].dim = 0;
     }
 
-    // find bad
+    // get box info
+    const mjtNum* pos1 =  d->geom_xpos + 3 * g1;
+    const mjtNum* mat1 =  d->geom_xmat + 9 * g1;
+    const mjtNum* size1 = m->geom_size + 3 * g1;
+    const mjtNum* pos2 =  d->geom_xpos + 3 * g2;
+    const mjtNum* mat2 =  d->geom_xmat + 9 * g2;
+    const mjtNum* size2 = m->geom_size + 3 * g2;
+
+    // find bad: contacts outside one of the boxes
+    for (int i=0; i < num; i++) {
+      // box sizes with margin
+      mjtNum sz1[3] = {size1[0] + margin, size1[1] + margin, size1[2] + margin};
+      mjtNum sz2[3] = {size2[0] + margin, size2[1] + margin, size2[2] + margin};
+      mju_scl3(sz1, sz1, kBoxRemoveMargin);
+      mju_scl3(sz2, sz2, kBoxRemoveMargin);
+      // mark as bad if outside box
+      if (mju_outsideBox(con[i].pos, pos1, mat1, sz1) ||
+          mju_outsideBox(con[i].pos, pos2, mat2, sz2)) {
+        con[i].dim = -1;
+      }
+    }
+
+    // find duplicates
     for (int i=0; i < num-1; i++) {
+      if (con[i].dim == -1) {
+        continue;  // already marked bad: skip
+      }
       for (int j=i+1; j < num; j++) {
+        if (con[j].dim == -1) {
+          continue;  // already marked bad: skip
+        }
         if (con[i].pos[0] == con[j].pos[0] &&
             con[i].pos[1] == con[j].pos[1] &&
             con[i].pos[2] == con[j].pos[2]) {
@@ -1546,6 +1577,7 @@ void mj_collideGeoms(const mjModel* m, mjData* d, int g1, int g2) {
     // consolidate good
     int i = 0;
     for (int j=0; j < num; j++) {
+      // good: maybe copy
       if (con[j].dim == 0) {
         // different: copy
         if (i < j) {

src/engine/engine_collision_primitive.h

@@ -16,6 +16,7 @@
 #define MUJOCO_SRC_ENGINE_ENGINE_COLLISION_PRIMITIVE_H_
 
 #include <mujoco/mjdata.h>
+#include <mujoco/mjexport.h>
 #include <mujoco/mjmodel.h>
 
 // define and extract geom info
@@ -47,32 +48,32 @@ int mjraw_SphereTriangle(mjContact* con, mjtNum margin,
                          const mjtNum* t1, const mjtNum* t2, const mjtNum* t3, mjtNum rt);
 
 // plane collisions
-int mjc_PlaneSphere     (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
-int mjc_PlaneCapsule    (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
-int mjc_PlaneCylinder   (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
-int mjc_PlaneBox        (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_PlaneSphere     (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_PlaneCapsule    (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_PlaneCylinder   (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_PlaneBox        (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
 
 // sphere and capsule collisions
-int mjc_SphereSphere    (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
-int mjc_SphereCapsule   (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
-int mjc_SphereCylinder  (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
-int mjc_CapsuleCapsule  (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_SphereSphere    (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_SphereCapsule   (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_SphereCylinder  (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_CapsuleCapsule  (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
 
 // box collisions: from engine_collision_box.c
-int mjc_CapsuleBox      (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
-int mjc_SphereBox       (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
-int mjc_BoxBox          (const mjModel* m, const mjData* d,
-                         mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_CapsuleBox      (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_SphereBox       (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
+MJAPI int mjc_BoxBox          (const mjModel* m, const mjData* d,
+                               mjContact* con, int g1, int g2, mjtNum margin);
 
 #ifdef __cplusplus
 }

src/engine/engine_util_misc.c

@@ -837,6 +837,26 @@ mjtNum mju_springDamper(mjtNum pos0, mjtNum vel0, mjtNum k, mjtNum b, mjtNum t)
 
 
 
+// return 1 if point is outside box given by pos, mat, size
+int mju_outsideBox(const mjtNum point[3], const mjtNum pos[3], const mjtNum mat[9],
+                   const mjtNum size[3]) {
+  // vector from pos to point, projected to box frame
+  mjtNum vec[3] = {point[0]-pos[0], point[1]-pos[1], point[2]-pos[2]};
+  mju_rotVecMatT(vec, vec, mat);
+
+  // outside
+  if (vec[0] > size[0] || vec[0] < -size[0] ||
+      vec[1] > size[1] || vec[1] < -size[1] ||
+      vec[2] > size[2] || vec[2] < -size[2]) {
+    return 1;
+  }
+
+  // inside
+  return 0;
+}
+
+
+
 // print matrix to screen
 void mju_printMat(const mjtNum* mat, int nr, int nc) {
   for (int r=0; r < nr; r++) {

src/engine/engine_util_misc.h

@@ -78,6 +78,10 @@ MJAPI void mju_decodePyramid(mjtNum* force, const mjtNum* pyramid,
 // integrate spring-damper analytically, return pos(dt)
 MJAPI mjtNum mju_springDamper(mjtNum pos0, mjtNum vel0, mjtNum Kp, mjtNum Kv, mjtNum dt);
 
+// return 1 if point is outside box given by pos, mat, size
+MJAPI int mju_outsideBox(const mjtNum point[3], const mjtNum pos[3], const mjtNum mat[9],
+                         const mjtNum size[3]);
+
 // print matrix
 MJAPI void mju_printMat(const mjtNum* mat, int nr, int nc);