Collision detection bug for certain face * face collisions.

[thomcc]

I believe I've found a bug in your physics or collision code. It's possible you're already aware of this, since after checking out older versions of the code, it seems to be present in all of them. The result of this bug is objects clipping through the geometry entirely.

An easily reproducible test case can be found by modifying your testphys/testphys.cpp code to only have a single cube that starts at {0.0f, 0.0f, 5.5f} (Here's a gist that has those changes available). When you run the program, it the cube will pass entirely through the world geometry as it falls.

I did a bit of sleuthing in a debugger to try and figure out the issue on my own. I did not have a ton of luck, but some. What I've found is that this seems to come from src.v in calcpoints ending up as a zero vector. Normalizing a zero vector ends poorly, and the rest of the code can't cope (and the simple solution of checking before normalizing and substituting a unit vector does not solve the problem). The zero comes from PlaneProjectOf returning zero on line 128 of gjk.h, but while I believe I understand that code, I don't understand it well enough to know how that should be handled (my guess is that that code is fine, and the caller should be checking for the zero vector).

A few things I've found do solve (well, avoid) the problem that might help track it down:

• Changing the starting position so that the x and y positions are not 0, 0. For example, going from {0.0f, 0.0f, 5.5f} to {0.0f, 0.1f, 5.5f} avoids it.
• Certain starting z positions don't trigger it. E.g. 1.9 won't, even though 1.85 and 1.95 will. Nothing below 0.7 seems to trigger it (until the cube begins to start inside the world geometry, at which point it always happens).
• Changing the orientation so that the collision does not occur as face * face avoids it.

Probably worth mentioning that I'm on a mac and using clang, and am using a modified version of glwin.h (you can see it here) that uses SDL to allow running from a mac. I don't believe this to be related in any way, since it happens at all optimization levels and otherwise works perfectly, but I obviously could be wrong.

P.S. This is a great project. Very straightforward code, not bloated or overengineered. Extremely pleasant to read through and learn from. Thanks for sharing it.

[melax]

Hi,
You are doing nothing wrong, the code is not properly handling the cases you presented. I easily reproduced the issue you described.

GJK works best (and is easiest) when there is a gap between things. When things are penetrating, its more difficult to write numerically robust code to find a the ideal plane between them.

Also... Because I had recently been doing various interaction demos with hand sized objects, it looks like some of the tune-able parameters (such as target gaps between things) are set for a smaller scale.

So after something has fallen for 5 meters, it will move a greater distance (about 12cm per frame) than the constant physics_driftmax (currently 0.03 or 3cm per frame). As a result, with steps this big, when the physics code first starts generating contacts (between the box and the ground) it might be already after they are sticking inside of each other. This is why increasing the drop height may cause tunneling, and another increase will cause it to work again.

The ContactPatch function is the one that takes this distance parameter so that it only looks within a given range. But, as mentioned, this current value isn't big enough. Rather than just increasing this parameter (currently drift_max) to another arbitrary number, we can calculate what it should be. At about line 420 in physics.h, Instead of:

  for(auto &c : ContactPatch(SupportFunc(rb,shape), SupportFunc(*cell), physics_driftmax) )

This could take into account the speed of the object to see if there is a potential that it may come into contact during the physics update.

    float distance_range = std::max(physics_driftmax, magnitude(rb->linear_momentum) *physics_deltaT / rb->mass);  // dont need to create potential contacts if beyond this range
    for (auto  cell : cells)
      for(auto &c : ContactPatch(SupportFunc(rb,shape), SupportFunc(*cell), distance_range) )
        contacts_out.push_back(PhysContact(rb, NULL, c));  

I tested this change, and a box with your specs for size, starting position, and height. It doesn't fall through the floor anymore.

Not a concern for the small box currently, but I just noticed another potential problem is how the code approximates the contact manifold in ContactPatch() To keep the physics engine as state-less as possible (immediate mode), contact information isn't cached/saved from previous frames. So to generate all the potential contacts it pulls the shapes apart and rolls one of them by a few degrees about each axis invoking gjk for each. This too is sensitive to size. I just realized that forgot to divide by 2 when generating the quat, so its rolling 4 degrees (instead of 2 as mentioned in the comment). It currently pulls things apart only by 0.1m (yet another hard coded fixed number found in code). This means that an 4 degree roll of anything much larger than a 1m box (1.4 face diagonal) could mean this function is temporarily putting the shapes in a penetrating position thus could result is less robust code being used for this spread of contact points.

I appreciate you trying out the various test cases. It does point to where the problems are.
When the falling box x,y is at exactly 0,0, there are a lot of points on the box and the world geometry that are all coincidentally coplanar - something we see in programmer art, but not as often with artist content. These are the sorts of corner cases that need to be dealt with.

Admittedly, in real applications, penetration happens and its nice to be able to deal with that better. Better expanding-polytope-algorithm is on my todo list. One workaround is to use the history (velocity), and backtrack to a point where things are properly separated and contacts can be generated robustly.

I'll have a break in a couple of weeks where i'll hopefully get the chance to go over some of this code to make it more robust and readable.

If you have some time, you might want to look for articles by Gino van den Bergen. He fully explains GJK, its challenges, and how to properly make it robust.

Thanks again,

[thomcc]

Very interesting, thanks for taking the time to respond!

I actually had continued trying to figure out the issue after submitting this, err, issue. I ended up adding a block that checked for a zero next.v after the nextSimplex call in Separated, and used calcpoints to generate a Contact that we'd return, after clearing the separation and coming up with a normal value -- the idea was that since clearly since next.v is the origin, the simplex contains the origin, and so we can fake a hitInfo and quit.

It sort of worked, but not in a convincing way at all, the bounce caused by the collision response would be a significantly higher (~2x) when coming from these Contacts than otherwise. Amusingly, clearing the separation was done so that the check in ContactPatch wouldn't fail, but changing the input to ContactPatch seems much more sensible in retrospect!

I'll definitely take a look into the articles you mentioned. Most of what I've used to understand the collision code had been from the description of the algorithm here, which was very helpful for understanding, although was clearly describing a different variant of the algorithm (one which doesn't track the closest point to the origin, give separation, etc). But more concrete information on the subject is definitely something I'm interested it.

Anyway, thanks again!

[melax]

Dealing with penetrating cases better...

Background

Even with good intersection avoidance, there will still be those problematic times when things end up getting stuck inside each other. When things are touching or penetrating, gjk's regular main loop terminates finding that the origin lies within the last simplex - a tetrahedron. From here, the expanding polytope algorithm (proposed initially by Gino I believe), is a way to find the best candidate plane. This EPA algorithm is harder to make robust. Even Dirk has suggested that a more brute force approach might be the way to go in this case. see http://www.gamedev.net/topic/649946-epa-expanding-polytope-algorithm/

Admittedly, the code I had for this was not thorough and would only work in the simplest of cases. When expanding the polytope, i would take the face closest to the origin and look for a minkowski support vert above that, but then only turned that triangle into 3 triangles ignoring if the new vert happened to also be above other faces on the current polytope.

Gino's has a paper discussing EPA "Proximity Queries and Penetration Depth
Computation on 3D Game Objects" http://movement.stanford.edu/courses/cs468-01-fall/Papers/van-den-bergen.pdf The paper does explain that you cant simply turn one triangle into 3 (figure 6), but rather he proposes tessellation (figure 7).

imho, what you really want to do here, is an actual convex hull expansion since each time you expand a face, it might affect more than just the current face, and perhaps even faces beyond the adjacent neighboring faces. Also, its better to minimize extra tessellation. Ideally one polygon per plane of each hull if using an n-gon rep. Whatever the contact situation, the minkowski sum itself should be a fairly large (non-flat) regular volume. Consequently some of the robustness issues that make convex hull calculation difficult wont be an issue. Now, its not necessary to compute the entire minkowski sum's convex hull, simple iterate/search based on the closest triangle/face/plane to the origin and simply terminate when you first hit a faces that is determined to be on the surface of the minkowski sum.

Results

Today i adapted the convex hull generation code for a first-pass-implementation to test this idea. Specifically replaced the greedy metric of most distant vertex to minkowski support above closest face, and also changed the termination condition to return first closest face (plane) it finds that doesn't have a supportmap point above it. Simple testing is showing this new version is handling deeper convex-convex penetration cases that would fail before. Still, without more testing, I'm not claiming that its bullet proof. Here's an example screenshot showing the difference:

With this new EPA, the proposed (minimal penetration) plane now has less distance by which opposing geometries cross. Furthermore, it actually corresponds to something that is on the minkowski sum. Full minkowski sum generated here only for visualization purposes.