ccdMPRPenetration() bug, wrong dir penetration depth

[Mikez2015]

Hello.

ccdMPRPenetration() return wrong dir penetration depth.

Please, see this video, bright green stick - it's dir penetration:
https://youtu.be/vN_u3Ig1C2c

I think this is due to this problem: ccdVec3PointTriDist2 computes different distance with / without witness points. #55

Really need help.. I hope the project is not abandoned.

[Mikez2015]

Bug in ccdVec3PointTriDist2

[Mikez2015]

Bug there:

if ((ccdIsZero(s) || s > CCD_ZERO)
            && (ccdEq(s, CCD_ONE) || s < CCD_ONE)
            && (ccdIsZero(t) || t > CCD_ZERO)
            && (ccdEq(t, CCD_ONE) || t < CCD_ONE)
            && (ccdEq(t + s, CCD_ONE) || t + s < CCD_ONE))

If we replace this condition with:
if (1)

issue is solved.

[SeanCurtis-TRI]

Just as a reality check, these are the lines of code to which you are referring, yes?

https://github.com/danfis/libccd/blob/master/src/vec3.c#L176-L180

The solution will have to be more nuanced than just setting that entire branch of logic to true. :) That block code only works if the query point projects onto the triangle and shouldn't be exercised unless that's true. I can imagine in your case it seems to work because the point in questions is positioned right at the threshold of successfully projecting.

The true bug is going to be more subtle. It's going to be related to the tolerances built into the successful projection test and what happens in the alternate branch when that threshold is crossed. Ideally, the two branches should, at the limit, produce the same answer as we approach the threshold; but they are apparently not doing so.

Looking at your video, I hypothesize that what we're seeing is that when the "has projected" test fails and we evaluate the distance from the point P to the nearest point on each triangle edge Qi, the closest pair (P, Qi)'s distance is essentially epsilon in length which means the direction of the vector between them is largely defined by rounding error.

Ultimately, if you could provide some code that reproduces this (as done in the issue you linked to), it would aid considerably in resolving this. (Although, it clearly increases the burden on you, it definitely increases the likelihood of the issue being addressed).

[Mikez2015]

Just as a reality check, these are the lines of code to which you are referring, yes?
Yes these lines

The solution will have to be more nuanced than just setting that entire branch of logic to true. :) That block code only works if the query point projects onto the triangle and shouldn't be exercised unless that's true. I can imagine in your case it seems to work because the point in questions is positioned right at the threshold of successfully projecting.

The true bug is going to be more subtle. It's going to be related to the tolerances built into the successful projection test and what happens in the alternate branch when that threshold is crossed. Ideally, the two branches should, at the limit, produce the same answer as we approach the threshold; but they are apparently not doing so.

I understand that this is not the best solution, but I needed at least some solution that would work. I was also surprised that it works, but it is. I checked on different models - boxes, cylinders, pyramids, moved them from different sides and from different angles, this is my strange solution that really works well. I saw some problems with the direction of the penetration vector only when using spheres.
Of course, I would like to know a truly mathematically correct solution, and not such as mine. But I'm not good at geometry.

Looking at your video, I hypothesize that what we're seeing is that when the "has projected" test fails and we evaluate the distance from the point P to the nearest point on each triangle edge Qi, the closest pair (P, Qi)'s distance is essentially epsilon in length which means the direction of the vector between them is largely defined by rounding error.

I also think that this is due to calculation errors, but it is definitely not related to float precision, since I tried the option with double too.
Unfortunately, to fix the ccdVec3PointTriDist2 function myself, I do not have enough knowledge of geometry.

Ultimately, if you could provide some code that reproduces this (as done in the issue you linked to), it would aid considerably in resolving this. (Although, it clearly increases the burden on you, it definitely increases the likelihood of the issue being addressed).

Which part of the code do you mean? You won’t read thousands of lines of code. I am completely sure that the problem is in the ccdVec3PointTriDist2 function, which is called from the findPenetr function, which is called from the ccdMPRPenetration function. Most likely the problem is exactly what you wrote "I hypothesize that what we're seeing is that when the" has projected "test fails and we evaluate the distance from the point P to the nearest point on each triangle edge Qi, the closest pair (P, Qi) 's distance is essentially epsilon in length which means the direction of the vector between them is largely defined by rounding error. " This is the only logical reason for the error.
I would be very grateful if you could help fix this, like many other libccd users.

[SeanCurtis-TRI]

I checked on different models

You'd have to check various configurations as well. You could change numerous models and not run into the situation where this change would explode if you just don't configure them correctly.

it is definitely not related to float precision

I didn't say it was related to "float precision". I said it was related to tolerances. You'll note that ccdIsZero actually tests the absolute value of the quantity under test against an epsilon. That epsilon value tends to scale with the precision available. But it's still a tolerance and that can lead to funny artifacts if your two branches don't have proper continuity.

One of the factors that is important, from a geometry sense, is the scale of the problem. The default epsilon value works best when all of the parameters are approximately unit magnitude. Ultimately, epsilon should really scale with respect to the size of the problem (otherwise, for large triangles, the absolute value of epsilon falls below rounding error, and for small triangles, the absolute value becomes a larger portion of the valid numerical space).

I agree with you, the complete diagnosis of the code is complex. And the solution more so. I hadn't intended to suggest you should implement the solution. I was merely documenting why your apparent solution is not a general solution that could/should be included in a PR. Resolving what appears to be a very real bug is going to be more complex than just setting the test to true.

Which part of the code do you mean?

I mean a piece of reproducible code that would allow someone who has no access to your code base the ability to see the bug. That would the single greatest contribution contained in this issue. It makes the bug definitive and empowers someone to address it.

In an ideal world, you should be able to reproduce this in just a few lines of code.

• Configure a triangle
• Select a point
• Call ccdVec3PointTriDist2
• perturb the point by a small amount (epsilon)
• Call ccdVec3PointTriDist2 again and show that the result jumps discontinuously.

The challenge is finding the triangle and point that is causing you problems because it is, ultimately, being generated inside libccd. Still, build enough things in debug build and you should be able to capture those parameters.

[Mikez2015]

The challenge is finding the triangle and point that is causing you problems because it is, ultimately, being generated inside libccd. Still, build enough things in debug build and you should be able to capture those parameters.

Thx a lot, I will try it.
I will continue to fight this error and write the process here.

[Mikez2015]

ccdVec3PointTriDist2

The function incorrectly considers the distance from the zero point to the triangle with the following coordinates:

A(1.650970, -0.102323, -2.994949);
B(-1.613953, -0.102323, 2.534958);
C(-1.613953, -0.102323, -2.994949);

The function returns the distance to the triangle equal to 0.143770, but the correct distance is 0.102323.

A couple of pictures for clarity:

I myself will continue to search for errors in this function, but I will be grateful for any help.

[Mikez2015]

One of the factors that is important, from a geometry sense, is the scale of the problem. The default epsilon value works best when all of the parameters are approximately unit magnitude. Ultimately, epsilon should really scale with respect to the size of the problem (otherwise, for large triangles, the absolute value of epsilon falls below rounding error, and for small triangles, the absolute value becomes a larger portion of the valid numerical space).

I tried to change the value of CCD_EPS in different directions and in different orders, but this did not produce any result.

[Mikez2015]

I began to research this comparison on the correctness of work:
https://github.com/danfis/libccd/blob/master/src/vec3.c#L176-L180

I found out that if the function does not work properly, the condition is not met
(ccdIsZero(t) || t > CCD_ZERO)
, because the variable "t" becomes negative.

Is it possible that a variable "t" is being calculated incorrectly?
https://github.com/danfis/libccd/blob/master/src/vec3.c#L173

[SeanCurtis-TRI]

This is a great effort; you've got the right idea.

However, I've looked at your example and I think you have an error. Here's the scenario you've created:

• Two triangles, lying on the y = -0.102323 plane.
• The query point is the origin (0, 0, 0).
• The distance from origin to triangle is 0.102323 if and only if the origin's projection onto the y = -0.10232 plane lies inside the triangle.

         +z
         ^
         ┆
    B    ┆
    │ ╲  ┆O
┄┄┄┄┄┄┄┄╲┼┄┄┄┄┄┄┄┄┄┄┄> +x
    │    ┆╲
    │    ┆  ╲
    C────┆────A
         ┆

Figure 1: "Bad" triangle drawn on the y = -0.102323 plane. The picture is merely illustrative showing the relative position of the triangle vertices.

For the "erroneous" triangle, you assume the correct distance is 0.102323 which implies the origin projects to a point inside the second triangle. However, that appears not to be the case. Consider the vector from A to B (AB). To be inside the triangle, the vector from A to O (AO) must lie on AB's "left" side. By my calculations, it doesn't. Here's some python code to illustrate that.

import numpy as np
def det(a, b):
  return a[0] * b[1] - a[1] * b[0]
# We're going to turn this into a 2D problem; because the y-value doesn't matter
# for testing the *projection* of the origin on the y = -0.102323 plane.
A = np.array((1.650970, -2.994949))
B = np.array((-1.613953, 2.534958))
p_AB = B - A
p_AO = -A
# for origin to project onto the triangle, this value must be <= 0.
det(p_AO, p_AB)
# evaluates to 0.6485673141370007

In fact, if I take this a step further, the projection of the origin lies 0.100994 units outside the triangle. And that means the origin is 0.14377027 units away from the triangle.

Unfortunately, this means the "error" condition isn't really an error condition. We'll have to keep digging.

[Mikez2015]

For the "erroneous" triangle, you assume the correct distance is 0.102323 which implies the origin projects to a point inside the second triangle. However, that appears not to be the case.

Stop, sorry, I'm a little confused. Let us take one moment to uttle. I correctly understood that you are claiming that the point
P (0,0)
does not belong to the triangle

A(1.650970,-2.994949);
B(-1.613953,  2.534958);
C(-1.613953,  -2.994949);

?

[SeanCurtis-TRI]

I was worried that might be a bit opaque. Sorry about that. I'm happy to walk through it.

I correctly understood that you are claiming that the point P (0,0) does not belong to the triangle

Correct. Everything else in the message was about showing that to be the case.

[Mikez2015]

I was worried that might be a bit opaque. Sorry about that. I'm happy to walk through it.

I correctly understood that you are claiming that the point P (0,0) does not belong to the triangle

Correct. Everything else in the message was about showing that to be the case.

The function ccdVec3PointTriDist2 thinks it belongs. I have to double-check this with another simple algorithm. Thanks for the help, I will write about the result later.

[Mikez2015]

Correct. Everything else in the message was about showing that to be the case.

I double-checked, you are absolutely right.
Then it turns out that the issue is that the function ccdVec3PointTriDist2 incorrectly determines whether the point belongs to the triangle.

So the error is at the beginning of the function, in this range:

https://github.com/danfis/libccd/blob/master/src/vec3.c#L142-L180