Also note that this is an intersection query. It's still GJK, but it's somewhat simplified since all we care about is whether the shapes intersect, not precisely where. If you want where, you want a variation of this algorithm which I'm not covering, as I haven't got the code for it (never needed it, never written it).

And I'm only going to deal with 3D, because that's the most useful case. It's fairly trivial to adapt the algorithm to 2D, and I'll make a note of how to do that. Don't ask me about 4-or-more-D, I'm taking a geometric approach to this, and I can't visualize that many dimensions.

We've got a few bits of information to store here. I'm going to roughly follow Casey's naming conventions. The newest point under consideration will always be called a. b, c, and d refer to previously considered points. v is our search vector (he calls this D, but that can get confused with the point, so I'm renaming it). n is the number of points left in our simplex from the previous step. That gives us the following:

Don't mind the random function definitions, we'll come back to what those are later. Actually, the body of intersect is up next...

As Casey points out, a lot of important information is encoded in the history of how we got our points to begin with. One of the most important bits, which he misses in his video, is contained here, in the outer loop. So let's take a look.

Do you see it? It's the dot product. If update is called, it will be called with an a such that its dot with d is positive (well, OK, points aren't vectors, so that's actually a-0). That little fact is going to pop up in every subcase of update, and it's going to simplify things tremendously.

The update function's job is to look at the existing simplex (encoded by b, c, d, and n) and decide how to update that given the addition of a. If it ever forms a simplex enclosing the origin, it returns true, and the outer loop exits. It has to handle several cases, each dependent on the number (n) of points left from the last iteration.

In high-level terms, update's main job is to find out what portion of the figure formed by adding a to the existing simplex is closest to the origin, to discard all other parts of the figure, and to update the search direction (v) such that we find a better a in the next iteration.

This case is pretty simple. We've got an old point in b, and adding a forms a line segment. Logically, the origin could be closest to one of the endpoints (a, b) or the line itself. However, we're going to be smart about this. Casey already covers why the feature of our line closest to the origin will never be b, but what about a?

Take a look at the (poorly drawn) diagram. We've got our two points and the origin. We know the origin won't be "behind" b (the shaded region) because we searched away from that region to find it (Casey also covers this case). But can't the origin be "behind" a? Well, no. Take a look at the angle x in that diagram. That's the angle between v (marked 0-b in the picture) and a-0 (remember, these are vectors, draw them tail to tail to find the angle between them), and we already know something about that angle - it's acute (or right, but let's ignore that for the moment). If it weren't, our dot product check in the outer loop would have terminated the search and we wouldn't be trying to update a 1-point simplex. Since x is acute, its supplementary angle is obtuse, which means that both of the other angles in the triangle are acute.

If the origin were "behind" a, then the inner angle in the diagram by a would have to be obtuse, but we already know that can't happen. That means that the origin is closest to the line segment between a and b. That means we can shade away another region of the diagram:

Now remember, we're in 3D here. This diagram is an edge-on view of the space we're searching. Those shaded regions aren't representing areas behind lines, they're representing volumes behind planes. These planes run through the endpoints of our line segment and are perpendicular to the line itself (even if my drawing is a little skewed). This will be important later on.

So adding a point to a point simplex always leaves us with a line, and in this case update does this:

Where cross_aba is:

The cross_aba function is just a useful shortcut to compute a vector perpendicular to a in the direction of b (those are its parameters a and b, not our global GJK state).

So now we have a line segment left over from our last iteration, and adding a new a is going to make a triangle. To the diagram!

OK, so before we fill this in, let's start with what we already know. Coming into this iteration, we have two points that form a line segment. We already know that that line segment defines a slab of space in which the origin can be found, so let's put our two planes on the diagram and mark that space as excluded. We also know that the origin can't be "behind" the line (relative to a), because a was found in the direction of the origin (think about it - that's where we looked for it at the end of the previous case).

Note that a itself is not constrained to this same region, but that doesn't matter for reasons that will soon become apparent. More diagrams!

The main feature here is the addition of another exclusion zone. This one's pretty straightforward. We searched along v (perpendicular to bc) and found a, and we know (by the test in the outer loop) that a is on the far side of the origin, so the origin can't be out past a along our search vector v.

Now here's the neat thing: no matter where we found a, its Voronoi region ends up entirely in the excluded space where we know the origin can't be, so we never need to actually consider it. The only regions we care about are beyond the edges ab, ac, and either above or below the triangle (that's in front of or behind the screen in the diagram above - this is 3D, after all). So the actual code is quite simple:

We compute our vector to the origin (ao) and our edge vectors (ab, ac). From there we compute the triangle normal (abc). From that, we compute vectors perpendicular to the edges (abp, acp), testing to see if the origin is out past each edge as we go. If it is, then we end with a new line simplex, updating b, c, and v appropriately. If it isn't beyond the edges, then it's either above or below the triangle, and we update appropriately (being mindful of the winding, as it will save us some work later on), and our simplex is now a triangle.

There's no easy way to diagram this out in 2D, so I'm not going to. In fact, when I was working this out the first time, I ended up building the amusing little construct you see to the right, along with some folded paper to represent the various planes the space around the figure is divided into. (It's almost, but thankfully not quite, too silly to work.)

Anyway, let's think about what we've got in front of us:

We've got our incoming simplex, which is a triangle, with its vertices stored in b, c, and d. Each of its edges has been tested during its construction, and we know that the origin is not in any edge's Voronoi region. That leaves us with an infinitely long triangular prism (running perpendicular to the triangle's face normal) in which the origin may be found. But is it even infinite? Well, no. We were careful to build our triangle such that it was wound with respect to which side the origin falls on - the origin isn't "below" the triangle, so half of our infinitely tall prism is gone. Further, we have our new point a, which is "above" not only the triangle but also the origin (by the dot product in the outer loop - this is the same as how we got our fourth plane in the previous case). This leaves us a proper triangular prism.

And, as was the case when we were building a triangle, a's Voronoi region is also guaranteed to fall outside that prism, so we needn't consider it.

That leaves us with three new edges, three new faces, and the actual volume within the tetrahedron formed by our triangle bcd and the new point a.

That's a big block of code, but look closely: everything below check_face: is a slightly simplified repeat of the n == 2 case, so ignore that and focus on the top portion.

All we do here is test the planes of the three new faces. If the origin is out in front of a face, we reduce our simplex to that face's edge along the incoming triangle (because that's the one that's already tested) and jump to a simplified version of our n == 2 code (we could just put that code in a function and call it, but it would recompute some values we already have on hand and do an unnecessary above/below test when we already know we're "above"). The check_face: code handles both the cases where the origin is nearest to one of the triangle edges and where it's directly above the face appropriately.

If the origin isn't above any of the face planes, then it's inside the tetrahedron, and we report having found an intersection.

This code is an intersection query. It doesn't tell you where the convex objects intersect. Further, it's intended to be used as a rendering optimization (is this object in view? can I skip rendering that shadow map? etc.), so it isn't even a hundred percent precise (that's the code you've been seeing in particular that's not precise, not the logic of the algorithm itself). It's designed this way for performance reasons, so let's discuss exactly what these limitations are and how to deal with them (because if you run it as is on any nontrivial inputs it will eventually get stuck in an infinite loop).

First off, since we don't care about where the intersection is, we can avoid normalizing a lot of our intermediate values since we're only interested in doing a series of plane side tests. That said, depending on the scale of your world, cross_aba can produce some monstrously large vectors. Your support mapping functions either need to deal with this, or you can normalize v at the end of each update and not worry about it beyond that.

Second, these are floating point numbers we're working with. You can get into situations where rounding errors leave you in an infinite loop as your support mapping bounces back and forth between a small number of nearby points. You could write a bunch of code to detect this, maybe sprinkle in some epsilons (which are not just 0.00001), but in practice it's such a rare case (on the order of one in ten thousand queries) that it's easier to just cap the total number of iterations through the outer loop and conservatively declare that there's an intersection if you blow past that limit (which it almost certainly is if your support mapping is reasonably well behaved).

Finally, take a look back over the three cases up above. The n == 1 case always exits with n equal to 2. And no case after that ever sets n to anything less than 2. So that can be collapsed into the setup code above the loop, leaving only two cases (n == 2, n == 3) for the actual update function to deal with.

Putting that all together and cleaning up a bit, we get:

The algorithm is fairly trivial to reduce to 2D. For one thing, the n == 3 case is entirely gone, since a triangle is enough to enclose the origin in 2D. The cross and cross_aba calls need to be replaced with the appropriate "find me a vector perpendicular to x" formula, and the end of the n == 2 case (where we check if we're above/below the triangle) simply becomes return true.