I recently came across Casey's wonderful video explaining almost all of the algorithm.
https://www.youtube.com/watch?v=Qupqu1xe7Io
I started my own implementation but have gotten lost.
It's mentioned in the (10 year old) video that discussion could take place on "the forums", which I assume to be referring to the MollyRocket forums which I don't think are up anymore.
I was hoping some of you who have implemented this in the past could give me some tips.
1) How do you check if a tetrahedron contains the origins. I did some Googling, but what I found was riddled with math terms I couldn't follow.
2) Do you perform this check first and only if it fails check which of the remaining faces/edges/single point point closest to the origin or does doing those checks somehow lead into determining if the tetrahedron contains the origin?
My thoughts are you can see if all the negative normals of the 3 unchecked faces point in the direction of the origin and be done with it, but I don't know if that's missing some cases or is much less efficient than what you really want to be doing.
Any help would be greatly appreciated.
GJK Implementation
Casey Muratori
801 posts
/ 1 project
Casey Muratori is a programmer at Molly Rocket on the game 1935 and is the host of the educational programming series Handmade Hero.
GJK Implementation
1) How do you check if a tetrahedron contains the origins.You are correct in your thoughts on this one - you know the four planes of the tetrahedron, so you need only evaluate the plane equation for each and "and" the results together.
2) Do you perform this check first and only if it fails check which of the remaining faces/edges/single point point closest to the origin or does doing those checks somehow lead into determining if the tetrahedron contains the origin?Well, as far as I know, the optimal way to structure the checks is still not quite known. Unlike the other cases where we have worked them out definitively, nobody (that I know of) has proven the optimal way to structure the checks (or even which checks to do, I suppose!)
That said, the usual way I do it is to check the three plane questions (you don't need to check the fourth, since you "came from" that direction already). This gives you a "bitfield", if you will, composed of 3 bits, each being set to positive if the origin lies outside the corresponding plane. You then branch to one of 8 code paths corresponding to the 8 combinations of bits, since each one lets you know something about the simplex.
One of those cases (no bits set) means that you are inside the simplex, so you have found an intersection.
Hope that helps,
- Casey
GJK Implementation
Another method for checking if a tetrahedron (with vertices A, B, C, D) contains a point (P, here P = 0) is checking if there exists a convex combination of the vertices yielding the points:
(A,B,C,D are 3d vectors, x,y,z,u are scalar variables).
As you can see, a convex combination is a linear combination with two additional properties. Ax + By + Cz + Du = P and x + y + z + u = 1 together give you a 4x4 linear equation system which you would need to solve. It should have a unique solution (if A, B, C, D span an actual tetrahedron and do not lie all on the same plane), then you can check if x,y,z,u >= 0.
This is probably something you would not want to do in an efficient implementation. But it gives theoretical insight on the number of if-statements you could need.
If you think further about this, you can solve the entire problem of deciding which part of the tetrahedron to choose from the solution of this equation system, more explicitly from the signs of x,y,z,u. We already saw one case: if x,y,z,u >= 0, then take the whole tetrahedron (ABCD). If only three have sign >= 0 (say x, z, u), then take the side consisting of these 3 vertices (ACD). If only two have sign >= 0 (say x, u), then take the corresponding edge (AD). If only one has positive (or non-negative) sign, say x, then the corresponding vertex is the right one (A). The case that all of them are negative or zero is not posible, since x + y + z + u = 1. Also we can apply the trick used by Casey in the video to further cut down the number of cases: Each of the simplices returned contains the vertex A, since from the previous step of the algorithm we can already deduce that x can not be negative or zero.
This cuts the number of if statements needed from 6, not to 5, not to 4, but to 3 if statements (if y >= 0, if z >= 0, if u >= 0, each yielding a different case).
But you can immediately see that I cheated. While tests for y,z,u >= 0 are very cheap the computation of y,z,u are not. They can be computed explicitly using the Cramer rule (which I think is what the papers mentioned in the video have done) or using Gauss elimination, even without doing a division. The Cramer rule is very computation intensive (at least for 4x4) and Gauss elimination without division could cost you precision of the result, as far as I see it (I am no numerical expert). The number of additions / multiplications you need for Gauss should be not too big (I estimate ~60 multiplications for the calculation of the signs of y,z,u), this should not be to bad compared to the calculations used for the plane equations.
With this one could also try to reduce the number of if statements in the triangular case. Casey lists 5 different results of the function, but two of them seem very similar (the cases "origin is above the triangle" and "origin is below the triangle"). Maybe we only need one of them ("origin is above or below the triangle"), so this reduces the number of cases to 4 and two if statements could theoretically suffice.
1 | Ax + By + Cz + Du = P where x + y + z + u = 1 and x,y,z,u >= 0 |
(A,B,C,D are 3d vectors, x,y,z,u are scalar variables).
As you can see, a convex combination is a linear combination with two additional properties. Ax + By + Cz + Du = P and x + y + z + u = 1 together give you a 4x4 linear equation system which you would need to solve. It should have a unique solution (if A, B, C, D span an actual tetrahedron and do not lie all on the same plane), then you can check if x,y,z,u >= 0.
This is probably something you would not want to do in an efficient implementation. But it gives theoretical insight on the number of if-statements you could need.
If you think further about this, you can solve the entire problem of deciding which part of the tetrahedron to choose from the solution of this equation system, more explicitly from the signs of x,y,z,u. We already saw one case: if x,y,z,u >= 0, then take the whole tetrahedron (ABCD). If only three have sign >= 0 (say x, z, u), then take the side consisting of these 3 vertices (ACD). If only two have sign >= 0 (say x, u), then take the corresponding edge (AD). If only one has positive (or non-negative) sign, say x, then the corresponding vertex is the right one (A). The case that all of them are negative or zero is not posible, since x + y + z + u = 1. Also we can apply the trick used by Casey in the video to further cut down the number of cases: Each of the simplices returned contains the vertex A, since from the previous step of the algorithm we can already deduce that x can not be negative or zero.
This cuts the number of if statements needed from 6, not to 5, not to 4, but to 3 if statements (if y >= 0, if z >= 0, if u >= 0, each yielding a different case).
But you can immediately see that I cheated. While tests for y,z,u >= 0 are very cheap the computation of y,z,u are not. They can be computed explicitly using the Cramer rule (which I think is what the papers mentioned in the video have done) or using Gauss elimination, even without doing a division. The Cramer rule is very computation intensive (at least for 4x4) and Gauss elimination without division could cost you precision of the result, as far as I see it (I am no numerical expert). The number of additions / multiplications you need for Gauss should be not too big (I estimate ~60 multiplications for the calculation of the signs of y,z,u), this should not be to bad compared to the calculations used for the plane equations.
With this one could also try to reduce the number of if statements in the triangular case. Casey lists 5 different results of the function, but two of them seem very similar (the cases "origin is above the triangle" and "origin is below the triangle"). Maybe we only need one of them ("origin is above or below the triangle"), so this reduces the number of cases to 4 and two if statements could theoretically suffice.
GJK Implementation
Or you can do it on one go if you have a decent 4x4 matrix inversion code.
Take a matrix ABCD which is constructed so that each collumn is the coordinates A, B, C and D with the 4th component 1. Invert it and then multiply it with the P which works because the 4x4 equation can be represented by ABCD*xyzu = P so after inverting and left multiplying you get ABCD' * ABCD * xyzu = I * xyzu = ABCD' * P
Then you can just check all the results are positive.
Take a matrix ABCD which is constructed so that each collumn is the coordinates A, B, C and D with the 4th component 1. Invert it and then multiply it with the P which works because the 4x4 equation can be represented by ABCD*xyzu = P so after inverting and left multiplying you get ABCD' * ABCD * xyzu = I * xyzu = ABCD' * P
Then you can just check all the results are positive.
GJK Implementation
Thanks for the pointers everyone.
For now I'm sticking to the geometric solution since that makes the most sense to me.
The bitfield and branch part makes perfect sense and the cases where a single bit is set really just means we're reiterating the triangle case, right?
The cases where 2 triangles normals point towards the origin is tricky though. I've given it a go, but my code is making either 5 or 7 total checks and I'm not seeing the way to make it 6 through every path. That means either I'm performing the checks in an unoptimal order or I've made mistakes in my checks.
Could anyone give this a once-over and tell me if there are any blatant mistakes I've made?
This is my thought process.
The simplex has 3 faces (ABC, ACD, ADB), 3 edges (AB, AC, AD) and one point (A) to check for being closest to the origin.
We just tested the BCD triangle and moved closer to the origin, so checking those is useless.
To help myself visualize the simplex, my mental image is the newest point A to be at the top above triangle BCD.
So it looks like a less MSPainty version of this:
If 2 triangle normals point towards the origin, then they have to have a shared edge in the middle.
For my code I refer to the points of this arrangement as T for top, L for "left" (clockwise) of the middle shared edge, M for the shared vertex, and R for "right" (counter-clockwise).
For example, if the normals of ABC and ACD are ones that point toward to origin, then
T = A;
L = B;
M = C;
R = D;
Which would look like so:
All that exposition out of the way, this is some pseudo code of how I'm running my checks
So after the 3 if's to determine the plane test, there are paths in this code that can add an additional 2 to 4 tests.
More importantly, are these checks right or am I misunderstanding the geometry?
Again, really appreciate the hints, folks.
For now I'm sticking to the geometric solution since that makes the most sense to me.
The bitfield and branch part makes perfect sense and the cases where a single bit is set really just means we're reiterating the triangle case, right?
The cases where 2 triangles normals point towards the origin is tricky though. I've given it a go, but my code is making either 5 or 7 total checks and I'm not seeing the way to make it 6 through every path. That means either I'm performing the checks in an unoptimal order or I've made mistakes in my checks.
Could anyone give this a once-over and tell me if there are any blatant mistakes I've made?
This is my thought process.
The simplex has 3 faces (ABC, ACD, ADB), 3 edges (AB, AC, AD) and one point (A) to check for being closest to the origin.
We just tested the BCD triangle and moved closer to the origin, so checking those is useless.
To help myself visualize the simplex, my mental image is the newest point A to be at the top above triangle BCD.
So it looks like a less MSPainty version of this:
If 2 triangle normals point towards the origin, then they have to have a shared edge in the middle.
For my code I refer to the points of this arrangement as T for top, L for "left" (clockwise) of the middle shared edge, M for the shared vertex, and R for "right" (counter-clockwise).
For example, if the normals of ABC and ACD are ones that point toward to origin, then
T = A;
L = B;
M = C;
R = D;
Which would look like so:
All that exposition out of the way, this is some pseudo code of how I'm running my checks
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 | TL = L - T;
TR = R - T;
TLM = cross(TL, TM);
TMR = cross(TM, TR);
// check the middle edge first to have an even triangle check on either branch path
if(dot(cross(TLM, TM), TO) < 0)
{
if(dot(cross(TL, TLM), TO) >= 0)
{
// we're on the "left" of TL according to the plane TLM
if(dot(TL, TO) >= 0)
{
// it's TL, TLxAOxTL
}
else
{
if(dot(TM, TO) >= 0)
{
// it's TM, TMxTOxTM
}
else
{
// it's T, TO
}
}
}
else
{
// we're in between TL and TM on plane TLM
// inside triangle
// can't be negative TLM or else plane test would be different
// it's TLM, +TLM
}
}
else
{
// we're to the "right" of TM according to the plane TLM
if(dot(cross(TMR, TR), TO) >= 0)
{
// we're to the "right" of TR according to the plane TMR
if(dot(TR, TO) >= 0)
{
// it's TR, TRxTOxTR
}
else
{
// it's T, TO
}
}
else
{
// we're in between TM and TR
if(dot(cross(TM, TMR), TO) >= 0)
{
// we're to the "left" of M according to the plane TMR
if(dot(TM, TO) >= 0)
{
// it's TM, TMxTOxTM
}
else
{
// it's T, TO
}
}
else
{
// we're in between TM and TR according to plane TMR
// can only be positive face normal
// it's TMR, +TMR
}
}
}
|
So after the 3 if's to determine the plane test, there are paths in this code that can add an additional 2 to 4 tests.
More importantly, are these checks right or am I misunderstanding the geometry?
Again, really appreciate the hints, folks.