ccdVec3PointTriDist2 computes different distance with / without witness points.

[hongkai-dai]

In ccdVec3PointTriDist2

libccd/src/ccd/vec3.h

Lines 193 to 202 in 7931e76

	/**
	* Returns distance^2 of point P from triangle formed by triplet a, b, c.
	* If witness vector is provided it is filled with coordinates of point
	* from which was computed distance to point P.
	*/
	CCD_EXPORT ccd_real_t ccdVec3PointTriDist2(const ccd_vec3_t *P,
	const ccd_vec3_t *a,
	const ccd_vec3_t *b,
	const ccd_vec3_t *c,
	ccd_vec3_t *witness);

, with the following input

P: 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
A: 0.06999979436939574029 -0.00000001708219171670 1.02781464659664711903
B: 0.06999979436939574029 -0.00000001708219171670 -0.81718535340335285433
C: -0.40000020563060423306 -0.00000002075754440556 1.02781464312938330963

the returned distance is different when we pass the witness points, versus setting witness=null.

When we set witness to null, the returned squared distance is 0. When we ask it to compute the witness point, the returned squared distance is 0.00000000000000031080. As we can see, all three points A, B and C have negative y value, hence the point P (origin) cannot have 0 distance to the triangle. ccdVec3PointTriDist2 computes the wrong value when the witness point is not passed in.

Inside the implementation of ccdVec3PointTriDist2, it wants to compute the distance from a point P to a triangle ABC. The way it does that is first to write the projection of point P on the plane coinciding with ABC as Q = A + s * AB + t * AC. In this test example, the point Q is inside the triangle ABC, so we only need to compute |PQ|² as the distance from P to ABC. The implementation takes different cod path to compute |PQ|².

1. When we ask the code to compute the witness point Q, the code first compute Q as Q = A + s * AB + t * AC, and then compute the squared norm of the vector PQ, as implemented in

   libccd/src/vec3.c

   Lines 182 to 189 in 7931e76

   	if (witness){
   	ccdVec3Scale(&d1, s);
   	ccdVec3Scale(&d2, t);
   	ccdVec3Copy(witness, x0);
   	ccdVec3Add(witness, &d1);
   	ccdVec3Add(witness, &d2);
   	dist = ccdVec3Dist2(witness, P);

2. When we don't ask the code to compute the witness point Q, the code expands the expression |PQ|² as |PQ|² = |PA + s * AB + t * AC|² = |PA|² + s²|AB|² + t² |AC|² + 2 st AB.dot(AC) + 2s * AB.dot(PA) + 2t * AC.dot(PA), it evaluates the summation of the quantities |PA|², s²|AB|², etc, as in the code

   libccd/src/vec3.c

   Lines 190 to 197 in 7931e76

   	}else{
   	dist = s * s * v;
   	dist += t * t * w;
   	dist += CCD_REAL(2.) * s * t * r;
   	dist += CCD_REAL(2.) * s * p;
   	dist += CCD_REAL(2.) * t * q;
   	dist += u;
   	}

   . This summation can cause numerical errors up to epsilon, when some of the quantities (like |PA|²) could be in the order of 1. The summation result end up being inaccurate. Hence although the squared distance should be non-zero, the summation result is 0.

The code path 1 (when computing the witness point Q) is numerically more robust. I think this is a bug as the distance result should be the same, with or without computing the witness point.

[sherm1]

It's probably worth mentioning that this bug is causing a failure in FCL as reported by @wxmerkt in flexible-collision-library/fcl#390 (running in double precision). In double precision distance tolerances are around 10⁻¹⁵ so a distance of 10⁻⁸ is still significant (by 7 orders of magnitude or so!).

[Mikez2015]

Hello, Has anyone fixed this issue?

[SeanCurtis-TRI]

I'm pretty sure that has gone untouched for the last year. The best strategy may simply be to create a PR. @hongkai-dai has characterized the problem well and a code change and unit test should be pretty straightforward.