Hello. I recently implemented the GJK algorithm for collision detection in my 3D engine, basing my implementation on a great video located here. The algorithm is robust, and seems to work for all convex objects and all cases so far. Now I need to determine the contact point, normal, and penetration depth of the collision. EPA seems like a good algorithm for this, because it uses the same principles of the Minkowski Sum and support mappings that GJK does.

Here are the only two resources I have been able to find that give much detail about EPA...

http://graphics.stanford.edu/courses/cs468-01-fall/Papers/van-den-bergen.pdf

http://uu.diva-portal.org/smash/get/diva2:343820/FULLTEXT01

The idea behind EPA is very simple. It starts with the tetrahedron simplex that GJK terminates with. It then projects the origin onto each triangle of the simplex. Lets call this point v. It then finds the triangle whose v is closest to the origin, and calls the support function using v as the direction. The triangle is subdivided, and new triangles are creating using the result of the support function as a new vertex. This is repeated over and over, so the algorithm essentially expands the simplex inside the Minkowski Sum. Once the distance between v and the support point is less than a certain threshold, we have found the penetration depth, and can easily determine the contact point and normal.

The only issue I'm having with EPA is the subdivision of the triangles. The only two resources I've found give two completely different methods for subdivision....

The first says to subdivide the edges of the triangle and use the support point to create 6 new triangles. The problem with this is, what the hell do you do with the adjacent triangles? They will have to be subdivided too, or else the mesh will have holes in it. This is a very ugly problem to solve, and he does not give much detail about the rest of the subdivision process.

The second one gives a more detailed explanation, but is completely different. Instead of subdividing the single triangle, we find every triangle whose normal is facing toward the support point. We then delete all of those triangles, and then create new triangles that all converge at the support point to fill the hole. This is what I'm trying to implement right now. It just doesn't seem like it would be very efficient. EPA is an iterative algorithm, and this must be done many times per frame.

I've looked and looked, and those were the only two decent resources I could find, and I cannot come up with a better method myself. Bullet uses EPA and I've looked at their source code, but it's utterly unreadable. So my question is, does anybody who has experience with EPA the best way of doing this? Or can anybody point me to a place with a different method? I'm not sure who came up with EPA, and I'm unable to find an official paper from the original author documenting the correct way to do it. Any help would be greatly appreciated.

Edit: I'm using method two now, and now I understand the methods they use to make it more efficient. I will post an explanation of the algorithm with sample code when I'm done

I think EPA was introduced by Gino v.d. Bergen. Personally I find it very difficult to make it numerically stable. The other issue you need to deal with is when GJK terminates with a simplex that is not a tetrahedron or even only gives you a fuzzy answer whether it is just very close or penetrating. I recommend testing objects with very different dimensions to see if you run into problems. In practice people use a brute-force fallback where they sample a fixed number possible directions on a unit sphere if things go bad. And they will go bad...

So why not use SAT instead?

My GJK implementation always terminates with a tetrahedron simplex if there is a collision, so that's not a problem. I almost have EPA working now. I added debugging code to draw the polytope and let me go through the algorithm step by step while rotating the view to inspect the polytope and make sure everything is working correctly. It expands the polytope correctly for several iterations, and then I start getting incorrect triangle normals. I'm almost positive I handle the triangle windings correctly, but I'll keep checking to find out what the problem is. As for the numerical stability... This is my first time writing a physics engine, so I don't expect everything to work 100% all the time

And I've never actually looked into using SAT in 3D. How is the performance? Does it work curved shapes? And does it compute the contact info?

Edit: I got it working now. It finds the correct penetration depth and normal. It seems to work very well. I've tested it with boxes and spheres of different sizes. I'll add other types of shapes soon. The only thing I'm confused about is how to find the contact points relative to each object