In this article I derive a ray-sphere intersection algorithm from it's math. While certain knowledge in mathematics and programming will be assumed, I am trying to show all the steps of the process.

Ray-sphere intersection is used to determine wether and where a ray hits a sphere. This is useful in selecting 3D objects with the cursor (well, spheres in this case), calculating bullet hits in games, etc., but you knew that. Why else would you look up such a specific topic?

So let's start with what we have: a sphere (your object, a hitbox, ...) and a ray (a bullet, your cursor, ...). What we we'd like to know is, if they intersect and what the intersection point(s) are, if any.

Definitions

A sphere is described as: \(0 = |P - C| - r\) where
\(P\) is an arbitrary point on the spheres hull, a 3-component vector
\(C\) is the spheres center location, a 3-component vector
\(r\) is the spheres radius, a scalar

So for any point on the sphere you subtract the center position. This always leaves a vector with the length of the radius. Subtracting the radius, you'll get zero. Makes sense, right?

A ray is described as: \(R = O + \hat D*t\) where

\(R\) is an arbitrary point on the ray, a 3-component vector
\(O\) is the rays origin, a 3-component vector
\(\hat{D}\) is direction of the ray, a 3-component unit vector
\(t\) is the length of the ray, a skalar

Now we could have included \(t\) in \(\hat{D}\) and sacrificed \(\hat{D}\)'s status as unit vector, but we need that distinction later. The equation should make sense as well, basically your line equation but in 3 dimensions.

Doing the Math

Intersection means, we want to find the point(s) where those to objects meet. Meaning where they have a point in common. Before, we determined that \(P\) is an arbitrary point on the spheres and \(R\) an arbitrary point on the ray. So in saying \(P=R\), we intersected those two equations. Now we need to see if it has a solution.

Intersecting gives us the following equation. We shall solve it for \(t\), the only unknown variable and the distance of intersection(s).
\(|O + \hat{D}t - C| - r = 0\)

So let's begin. Add \(r\) to both sides and remove the absolute by squaring both sides.
\((O + \hat{D}t - C)^2 = r^2\)

Expand the left term. Note that this is the dot product!
\((O + \hat{D}t - C) ⋅ (O + \hat{D}t - C) = r^2\)

Multiply the terms on the left side (each with the others)
\(O^2 + O\hat{D}t - OC + O\hat{D}t + \hat{D}t^2 - C\hat{D}t - OC -C\hat{D}t + C^2 = r^2\)

Add the terms that can be added
\(O^2 + 2O\hat{D}t - 2OC + \hat{D}t^2 - 2C\hat{D}t + C^2 = r^2\)

Rearrange to \(at^2 + bt + c = 0\)
\(\hat{D}t^2 + 2O\hat{D}t - 2C\hat{D}t + O^2 - 2OC + C^2 - r^2 = 0\)

Still re-arranging. Making a single \(b\) term by pulling out \(t\)
\(\hat{D}t^2 + 2\hat{D}t ⋅ (O - C) + O^2 - 2OC + C^2 - r^2 = 0\)

Making the \(c\) term simpler. Binomial theorem!
\(\hat{D}t^2 + 2\hat{D}t ⋅ (O - C) + (O - C)^2 - r^2 = 0\)
"a" term: \(\hat{D}t^2\)
"b" term: \(2\hat{D}t ⋅ (O - C)\)
"c" term: \((O - C)^2 - r^2\)

Now we could try to solve for \(t\) via the normal quadratic formula, but there is a convenient optimization: a unit vector squared equals 1. This is because:

Our equation is thus only

and we can apply the reduced quadratic formula

Doing so, we also see that the factor \(2\) in \(p\) cancels out with the divisor in the \(\frac{p}{2}\) term and we are left with a not too ugly equation for the intersection distances.

From the equation you see there are multiple cases, determined by the term under the root (the discriminant \(\Delta\)).

• \(\Delta < 0 \rightarrow\) No real solution exists, the root would be negative. Objects do not intersect.
• \(\Delta = 0 \rightarrow\) The ray grazed the sphere in just a single point, meaning \(t_1 = t_2\)
• \(\Delta > 0 \rightarrow\) This is a direct hit, so there is a distinguished entry and exit point.

With the distance known, you can then calculate the intersection points by inserting \(t\) back into the ray equation: \(R = O + \hat{D}*t\)

Writing the Code

Before you think about implementing, think about what you want from the function, what its inputs will be. Then draft it's signature. I am using C++ with the GLM vector math library in this example, but feel free to try in a language of your choice. Anyways, this is what I came up with:

int RaySphereIntersection(
    const vec3& rayPos, const vec3& rayDir, 
    const vec3& spherePos, float sphereRadius, 
    float& dist1, float& dist2
);

The first two arguments define our ray. Passed by const reference to avoid unnecessary copies. Remember that rayDir is supposed to be a unit vector, so make sure it's normalized. The next two arguments are our sphere.

The output shall be the two possible intersection distances. I did not include the calculation of the hit locations in the function as there are use cases where it's not needed and would just be performance degrading. Better calculate it separately and on demand. As return value I chose the number of actual intersections. This allows for concise intersection checks in the code.

(!) Note that the distances can be negative if the the ray is pointing away from the object.

Now, without further ado, the code. Where possible, results were cached to avoid unnecessary calculations. The discriminant is checked for being negative, if it is, 0 (as in zero solutions) is returned. For the other cases, the calculations are done and the amount of solutions is returned. Note that in the last line I compare against a very small and positive value, but not zero. Due to the how floating point numbers behave the result might not get to be exactly zero and as a consequence I gave a certain margin to make sure "1" is returned when the difference is reasonably small.

int RaySphereIntersection(
    const vec3& rayPos, const vec3& rayDir, 
    const vec3& spherePos, float sphereRadius,
    float& dist1, float& dist2)
{
    vec3 o_minus_c = rayPos - spherePos;

    float p = dot(rayDir, o_minus_c);
    float q = dot(o_minus_c, o_minus_c) - (sphereRadius * sphereRadius);

    float discriminant = (p * p) - q;
    if (discriminant < 0.0f)
    {
        return 0;
    }

    float dRoot = sqrt(discriminant);
    dist1 = -p - dRoot;
    dist2 = -p + dRoot;

    return (discriminant > 1e-7) ? 2 : 1;
}

Now, it's not optimal yet. Neither in performance, nor in being generalized. It stays however close to it's roots and serves as base for further exploration.

Usage

Finally, an example on how to use the code. Also, I'd like to emphasise again that the direction vector must be normalized, otherwise the code won't work correctly.

float d1, d2;
if(RaySphereIntersection(rayOrigin, normalize(rayDirection), 
   spherePos, radius, d1, d2) > 0)
{
    // they intersect. d1 & d2 are valid now.
}

Assignments

• Can you optimise RaySphereIntersection() so it returns earlier in case of no intersections?
• Does your compiler optimize the calculation of dist1 and dist2 by swapping the operands?
• Write a function to calculate the intersection points.

Next Post Previous Post