I'm trying to implement a collision detection system to help solve constrained non-linear optimization problems in robotics as my master thesis.

Basically, what I need out of the system is the penetration depth and contact normal between each pair of potentially colliding bodies. I was thinking of using the GJK algorithm as a basis to implement the narrow phase of the system on GPUs using CUDA/OpenCL but wasn't sure if this line of thinking would get me very far due to the inherent lack of ‘numerical robustness’ of GJK using 32 bit floating point representations (the de facto standard on GPUs). There is also the question of code divergence (each parallel thread running GJK may not run the same instructions).

As a starting point, we would only be interested in determining collisions between capsules, spheres and cylinders. Keeping this in mind as well, is there a better way to do this? Are there documented ways to overcome issues with GJK and numerical stability on 32 bit floating point representations? Are there simpler or more efficient ways of implementing narrow phase collision detection systems on GPUs?

Thanks very much for your time and help! I'm totally new to this field and any help would be greatly appreciated, especially from those with prior experience implementing similar systems ?

Don't have much to offer in the way of a GJK implementation, but if the numerical instability with single precision floating point is due to the a limitation of the range of values, then GPUs do have support for double precision float( unless we are talking about GPUs 8+ years ago ).

GJK will not give you penetration depth. GJK computes the distance between two disjoint convex shapes. If you want penetration depth you need SAT, EPA or whatever you like.

GJK is indeed a challenge to get robust in 32bit precision. But if you stick to convex polyhedra you will be fine. You can also use spheres and capsules, but you use the GJK with the sphere center/inner capsule segment and add the radius later. Here are two presentations I gave some years ago that might be useful. I would express the cylinder as a polyhedron. If you don't want convex hull, I recommend to replace cylinders with boxes.

http://media.steampowered.com/apps/valve/2015/DirkGregorius_Contacts.pdf

http://media.steampowered.com/apps/valve/2013/DGregorius_GDC2013.zip

For GJK I recommend Erin Catto's presentation here. It also comes with sample source. It is written in a way that it straight forward to port to 3D.

https://box2d.org/files/ErinCatto_GJK_GDC2010.pdf

https://box2d.org/files/ErinCatto_GJK_Source.zip

This is a very ambitious project for a master thesis. That will keep you busy for a while! Just as a warning ?

@cgrant Thanks a lot for your reply! Indeed, you are right that modern GPUs do support 64 bit floating point, however I believe the throughput for 64 bit precision is terrible (something like 1/32 ) of the throughout of 32 bit versions for the graphics card I have (GTX 1080) and that just won't cut it :(

@Dirk Gregorius Thanks a lot for your inputs Dirk and for sharing the amazing presentations! Apologies for my mistake, I meant the EPA algorithm used in conjunction with GJK is also prone to numerical instability using 32 bit values. Am I right in my understanding?

Thanks for the heads up about the effort involved! At this stage, I'm just doing some reconnaissance and the implementation bits are yet to be finalized so there is still a lot of room for changes in design and implementation strategies.

If the implementation is very challenging and the applications we're looking into are such that we're willing to sacrifice a fair bit of accuracy for performance, would you also say that it can be an option just to use simple intersection tests (Sphere-Sphere, Capsule-Sphere etc) especially if all bodies involved in the scene are represented using these primitives? The reason I ask is that for our purposes (motion planning for robotics), we're primarily concerned about the robot not running into obstacles and we set up the problem as a constrained optimization problem in which one of the constraints to be satisfied is that the penetration depth between any two rigid bodies is <0. Even if these obstacles are not a true representation of real world obstacles (due to a bounding volume like a sphere, capsule etc) being used to model the object, could the simpler calculation of penetration depth between these bounding volumes be sufficient for our purposes? This would probably be over-conservative and sub-optimal but lead to very quick computation times especially when used on a GPU.

Thank you very much for your time and help! Your thoughts and inputs are much appreciated ?

I think for such a project I would use spheres, capsules and boxes. For boxes you can easily implement a specialized SAT. I would start with a stupic BF implementation first though and once you have some tests you can start optimizing. A good example for an optimized box vs box SAT can be found in ODE (dBoxBox iirc). For the other ones you can implement what I show in my presentation. Since you are in robotics you probably also don't need capsule stacking, so you can skip the case where you try to create two contact points for capsule vs capsule. My advice for you is to always start with simplest stuff so you can iterate on working solutions!

I am not aware of any publicly available implementation of EPA that deals with all the numerical problems. You would probably need a third brute-force algorithm in the case EPA fails. GJK/EPA is IMO a bad choice, but you can definitely ship with it and make it work. It is your choice.

I am not aware of any publicly available implementation of EPA that deals with all the numerical problems.

See Prof. Stephen Cameron's code from the 1990s:

http://www.cs.ox.ac.uk/people/stephen.cameron/distances/

That solved the numerical problems. It's not free code; I paid royalties to use it when I sold a product back in the 1990s. You can look at it on the Oxford site, but the comments are stripped. Prof. Cameron died last year, so licensing the code will be tough. I found that numerical bug and worked with him on fixing it, and have a copy with comments, but can't distribute it.

The problem is terminating the search for the closest set of points. The metric for that has a small difference between two large values numerical problems and can get into an infinite loop.

@Dirk Gregorius Thank you very much for your help and suggestions! That clears a lot of things for me and as suggested, I will work on incrementally improving the solution starting from extremely simple stuff. Lovely presentations! Thanks a lot for sharing them, your advice is very helpful for students like me!

@Nagle Thank you very much for sharing this link and the useful information! I'm very sorry to hear about his death, may his soul RIP. I will have a look at the code in detail and see if we can work on creating a working solution with GJK and EPA, lest its not overkill for our particular use-case. It will probably take me sometime to understand the uncommented code, I'll get back to you in case I don't understand something, thanks a ton for the description of what the problem is! Is the code which is online free of numerical errors? Is there any way to use the code for research (non-profit) without paying a fee?

John, this is GJK code, right? GJK is well understood these days in my opinion. I was talking about Expanding Polytope Algorithm (EPA).