ccdVec3PointTriDist2 avoids a division by zero for ill-formed triangl…

…es (fixes danfis#51)
tgn3000 · Dec 22, 2018 · c02958e · c02958e
1 parent 63d3a91
commit c02958e
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 3 deletions.
diff --git a/src/testsuites/vec3.c b/src/testsuites/vec3.c
@@ -270,4 +270,11 @@ TEST(vec3PointTriDist)
     ccdVec3Set(&P0, -.5, .5, 0.);
     assertTrue(ccdVec3Eq(&w, &P0));
     //fprintf(stderr, "dist: %lf\n", dist);
+
+    ccdVec3Set(&a, -0.36715889, 0.288464308, 0.000100158155);
+    ccdVec3Set(&b, -0.0222680569, 0.0171524286, -2.88337469e-06);
+    ccdVec3Set(&c, 3.0792e-01, -2.4249e-01, 3.8363e-05);
+    ccdVec3Set(&P, 0., 0., 0.);
+    dist = ccdVec3PointTriDist2(&P, &a, &b, &c, &w);
+    assertTrue(dist < 0.0000001);
 }
diff --git a/src/vec3.c b/src/vec3.c
@@ -149,7 +149,7 @@ ccd_real_t ccdVec3PointTriDist2(const ccd_vec3_t *P,
     // computed.
 
     ccd_vec3_t d1, d2, a;
-    ccd_real_t u, v, w, p, q, r;
+    ccd_real_t u, v, w, p, q, r, d;
     ccd_real_t s, t, dist, dist2;
     ccd_vec3_t witness2;
 
@@ -164,8 +164,14 @@ ccd_real_t ccdVec3PointTriDist2(const ccd_vec3_t *P,
     q = ccdVec3Dot(&a, &d2);
     r = ccdVec3Dot(&d1, &d2);
 
-    s = (q * r - w * p) / (w * v - r * r);
-    t = (-s * r - q) / w;
+    d = w * v - r * r;
+    if (ccdIsZero(d)){
+        // To avoid division by zero for zero (or near zero) area triangles
+        s = t = -1.;
+    }else{
+        s = (q * r - w * p) / d;
+        t = (-s * r - q) / w;
+    }
 
     if ((ccdIsZero(s) || s > CCD_ZERO)
             && (ccdEq(s, CCD_ONE) || s < CCD_ONE)