It’s been quite a while since I posted on the blog (a) at all, and (b) about a topic remotely related to the name of my domain and subsequently my book. Thanks to a gentle nudge from my friend and colleague Pål-Kristian Engstad I’ll try to rectify the situation with this post about how one could go about optimizing a sphere-triangle intersection test, suitable for a SIMD-implementation (such as on the PS3 SPUs).
This article isn’t so much about producing useful code, but more about the steps that you should be able to identify and perform when optimizing vector-heavy code. I also briefly discuss how to deal with separating axis tests when one of the objects is a sphere (which doesn’t have any obvious features).
Testing sphere vs. triangle
There are several different methods we can use to determine the intersection between a sphere and a triangle.
One is to “inflate” or “grow” the triangle by the sphere (forming the Minkowski sum of them) while simultaneously shrinking the sphere to a point. The test then becomes that of testing if the sphere center (the point) lies inside the spherically inflated triangle. In practical terms, this resolves into testing containment of the sphere center in one of 3 spheres (at the triangle vertices), 3 cylinders (aligned to the triangle edges), or a triangular prism (centered around the triangle itself).
Another option is to use the method of testing for a separating axis. Since this method is usually applied to (convex) polygons and polytopes only, and not for spheres, I thought it is a more interesting solution to elaborate on here, since I get to touch on two things at once.
Before continuing, let’s briefly review the separating axis test, to see where it comes from and how it works.
Justifying separating-axis testing
If two convex objects, C1 and C2, are not intersecting there must exist a pair of points P1 and P2, one from each object, such that no other pair of points are closer than these points. Given such a pair, we can insert a separating plane between P1 and P2, with normal P2 - P1. We call P2 - P1 a separating axis.
Since there is an infinite number of points on either object, which is not practical for testing, we classify points into a finite number sets we call features. We do this in such a way that all points of a feature share the same test axis for separation. This way, it is sufficient to form candidate axes for all pairs of features from each object, and see if we can find object separation on any of those axes.
If, after examining all candidate axes, we cannot find separation, we must conclude that C1 and C2 are intersecting.
For convex polyhedra (aka polytopes) features are the vertices, edges, and faces of the polytopes.
Separating axis testing for a sphere
From the above we conclude that a triangle has three kinds of features that are involved in testing: vertices (3), edges (3), and a face (1). But what about a sphere? Well, if you followed the justification of the test, it should be clear a sphere has only one feature, but it varies from test to test: the point on the sphere surface closest to each given triangle feature!
If we define the problem as testing the intersection between a triangle T, defined by the vertices A, B, and C, and a sphere S, given by its center point P and a radius r, the axes we need to test are therefore the following:
- normal to the triangle plane
- perpendicular to AB through P
- perpendicular to BC through P
- perpendicular to CA through P
- through A and P
- through B and P
- through C and P
We can now proceed with the actual testing, below.
Testing if sphere lies outside the triangle plane
The sphere S cannot intersect the triangle T if the sphere center P lies further away from the plane of T than r units. This can be expressed as follows:
Translate problem so sphere is centered at origin A = A - P B = B - P C = C - P Compute a vector normal to triangle plane (V), normalize it (N) V = cross(B - A, C - A) N = V / sqrt(dot(V, V)) Compute distance d of sphere center to triangle plane d = abs(dot(A, N)) separated = d > r
Since both d and r are positive, we can square both sides of the test; this allows us to remove the abs() too. N is only used once, so we can substitute the expression in the calculation of d, and move out the constant sqrt() expression from within the dot product of A and N:
A = A - P B = B - P C = C - P V = cross(B - A, C - A) d = dot(A, V) / sqrt(dot(V, V)) separated = d^2 > r^2
d^2 is equivalent to dot(A, V)^2 / dot(V, V) and since dot(V, V) is always positive we can multiply by dot(V, V) on both sides of the inequality, giving the final form:
A = A - P B = B - P C = C - P V = cross(B - A, C - A) d = dot(A, V) e = dot(V, V) separated = d * d > r * r * e
Testing if sphere lies outside a triangle vertex
For the axis through the sphere center P and the vertex A we have separation if:
Translate problem so sphere is centered at origin A = A - P B = B - P C = C - P Compute distance between sphere center and vertex A d = sqrt(dot(A, A)) The plane through A with normal A (”A - P”) separates sphere iff: (1) A lies outside the sphere, and separated1 = d > r (2) if B and C lie on the opposite side of the plane w.r.t. the sphere center separared2 = dot(B - A, A) > 0 separared3 = dot(C - A, A) > 0 separated = separated1 & separated2 & separated3
(Here we assume that A is the feature that separates T and S. As such, if B or C are closer, in the projection, then A cannot possibly be the separating feature. B or C might be, but we will find that out when testing for B and C as the potential separating feature, so we do not need to test that case here. The idea is to not do any early outs in this code.)
As before we square terms to get rid of the sqrt() and also use the linearity (or distributive properties) of the dot product to remove some vector subtractions, to give the optimized pseudocode:
A = A - P B = B - P C = C - P aa = dot(A, A) ab = dot(A, B) ac = dot(A, C) separated = (aa > r * r) & (ab > aa) & (ac > aa)
If we do the same thing for the remaining two axes and combine the tests we arrive at the final form for testing all vertices:
A = A - P B = B - P C = C - P aa = dot(A, A) ab = dot(A, B) ac = dot(A, C) bb = dot(B, B) bc = dot(B, C) cc = dot(C, C) rr = r * r separateda = (aa > rr) & (ab > aa) & (ac > aa) separatedb = (bb > rr) & (ab > bb) & (bc > bb) separatedc = (cc > rr) & (ac > cc) & (bc > cc) separated = separateda | separatedb | separatedc
Testing if sphere lies outside a triangle edge
To determine if there is a separating plane through the triangle edge AB with a normal parallel to the axis through the sphere center P and perpendicular to AB, we first need to compute the point Q on the line AB to which P projects. This Q is really just the projection of P onto the plane and must, by construction, fall onto the line AB (which is why we project P onto AB). We compute Q as follows:
AB = B - A t = dot(P - A, AB) / dot(AB, AB) Q = A + t * AB
(For a full justification of why Q is computed like this, see Section 5.1.2 of RTCD.)
To test separation we now need to see if Q is further than r from P and if C lies on the opposite side of the plane through AB with normal PQ.
Compute Q AB = B - A t = dot(P - A, AB) / dot(AB, AB) Q = A + t * AB Test separation QP = P - Q QC = C - Q separated1 = sqrt(dot(QP, QP)) > r separated2 = dot(QP, QC) < 0 separated = separated1 & separated2
We can get rid of the sqrt() by squaring both sides of the first inequality (as both sides are positive) and we can prepare to get rid of the division by multiplying both sides of the inequalities by the positive number e^2:
AB = B - A d = dot(P - A, AB) e = dot(AB, AB) Q = A + (d / e) * AB QP = P - Q QC = C - Q separated1 = dot(QP * e, QP * e) > r * r * e * e separated2 = dot(QP * e, QC * e) < 0 separated = separated1 & separated2
We finalize getting rid of the division through the following transformation:
AB = B - A d = dot(P - A, AB) e = dot(AB, AB) Q = A * e + d * AB QP = P * e - Q QC = C * e - Q separated = [dot(QP, QP) > r * r * e * e] & [dot(QP, QC) < 0]
Doing all three edge tests together results in the following pseudocode:
AB = B - A d1 = dot(P - A, AB) e1 = dot(AB, AB) Q1 = A * e1 + d1 * AB QP1 = P * e1 - Q1 QC = C * e1 - Q1 separated1 = [dot(QP1, QP1) > r * r * e1 * e1] & [dot(QP1, QC) < 0] BC = C - B d2 = dot(P - B, BC) e2 = dot(BC, BC) Q2 = B * e2 + d2 * BC QP2 = P * e2 - Q2 QA = A * e2 - Q2 separated2 = [dot(QP2, QP2) > r * r * e2 * e2] & [dot(QP2, QA) < 0] CA = A - C d3 = dot(P - C, CA) e3 = dot(CA, CA) Q3 = C * e3 + d3 * CA QP3 = P * e3 - Q3 QB = B * e3 - Q3 separated3 = [dot(QP3, QP3) > r * r * e3 * e3] & [dot(QP3, QB) < 0] separated = separated1 | separated2 | separated3
Yet again we translate there problem so the sphere is at the origin (P = 0), which, after some shuffling to remove unary minuses and making all comparisons greater-than, turns the pseudocode into the final version:
A = A - P B = B - P C = C - P AB = B - A BC = C - B CA = A - C d1 = dot(A, AB) e1 = dot(AB, AB) d2 = dot(B, BC) e2 = dot(BC, BC) d3 = dot(C, CA) e3 = dot(CA, CA) Q1 = A * e1 - d1 * AB QC = C * e1 - Q1 Q2 = B * e2 - d2 * BC QA = A * e2 - Q2 Q3 = C * e3 - d3 * CA QB = B * e3 - Q3 separated1 = [dot(Q1, Q1) > r * r * e1 * e1] & [dot(Q1, QC) > 0] separated2 = [dot(Q2, Q2) > r * r * e2 * e2] & [dot(Q2, QA) > 0] separated3 = [dot(Q3, Q3) > r * r * e3 * e3] & [dot(Q3, QB) > 0] separated = separated1 | separated2 | separated3
All together now
Finally, let’s combine all three (or 7) tests together in one final block of pseudocode:
A = A - P B = B - P C = C - P rr = r * r V = cross(B - A, C - A) d = dot(A, V) e = dot(V, V) sep1 = d * d > rr * e aa = dot(A, A) ab = dot(A, B) ac = dot(A, C) bb = dot(B, B) bc = dot(B, C) cc = dot(C, C) sep2 = (aa > rr) & (ab > aa) & (ac > aa) sep3 = (bb > rr) & (ab > bb) & (bc > bb) sep4 = (cc > rr) & (ac > cc) & (bc > cc) AB = B - A BC = C - B CA = A - C d1 = ab - aa d2 = bc - bb d3 = ac - cc e1 = dot(AB, AB) e2 = dot(BC, BC) e3 = dot(CA, CA) Q1 = A * e1 - d1 * AB Q2 = B * e2 - d2 * BC Q3 = C * e3 - d3 * CA QC = C * e1 - Q1 QA = A * e2 - Q2 QB = B * e3 - Q3 sep5 = [dot(Q1, Q1) > rr * e1 * e1] & [dot(Q1, QC) > 0] sep6 = [dot(Q2, Q2) > rr * e2 * e2] & [dot(Q2, QA) > 0] sep7 = [dot(Q3, Q3) > rr * e3 * e3] & [dot(Q3, QB) > 0] separated = sep1 | sep2 | sep3 | sep4 | sep5 | sep6 | sep7
It should be noted that throughout I’ve been deliberately sloppy with the separation criteria (i.e. greater-than vs. greater-equal, etc) to allow for the final version of the code to perform the same comparison for several tests in parallel. (Using SIMD instructions, the sixteen scalar greater-than comparisons can be packed into four SIMD comparisons.) Also, if implemented as derived, the resulting code is not the greatest example of floating-point robustness since we have cross-multiplied and squared terms to get rid of divisions and square roots in our quest for faster code. Finally, like Knuth, I have only “proved” the pseudocode to be correct, I haven’t actually tested it, so beware bugs. Caveat emptor, in other words!
Aside: there is no “separating-axis theorem”
As a side note, I feel compelled to point out that there is no such thing as a “separating-axis theorem” despite people talking about it in academic papers and even writing Wikipedia articles about it! Why? Because it is not a theorem and no one has ever proven it in print, nor would anyone attempt to, because it is a trivial consequence of the separating hyperplane theorem, which is a real proven and fundamental theorem in convex analysis! The correct vernacular is “separating-axis test.” Get it right.