Howto Raytracer: Ray / Sphere Intersection Theory

Categories:

I will show you how to calculate the intersection point of a ray with a sphere.

A ray r(t)r(t) can be represented by a point on the ray ee and the ray's direction dd. Any point on that ray can then be reached by walking tt times in the direction of dd from ee as r(t)=e+tdr(t)=e + t d. The set RR of all points on the ray is then given by: R={r(t)tR}R = {r(t) \mid t \in \mathbb{R}}

A sphere SS can be represented by its center point cc and its radius rr. The set of all points xx whose distance from the center cc of the sphere equals the radius rr, is by definition the set of points on the sphere: S={xc=rxR3}S = { ||x - c|| = r \mid x \in \mathbb{R}^3}. (x||x|| denotes the length of a vector.)

sphere equation

To find the intersection of the ray and the sphere now, we have to find the points that are in both sets. So we check if a point r(t)r(t) on the ray also fullfills the distance equation of the sphere:

r(t)c=r||r(t) - c|| = r

We now simply have to solve this equation. The way we do it is by rewriting the length of the vector as a dot product. For any vector x=(a,b,c)x = (a,b,c), its (euclidean) length is given by x=a2+b2+c2||x|| = \sqrt{a^2 + b^2 + c^2}. The dot product \cdot of a vector xx with itself is the sum of its squared components xx=a2+b2+c2x \cdot x = a^2 + b^2 + c^2.

So in the end we get the following relation:

r(t)c=(r(t)c)(r(t)c)=r(r(t)c)(r(t)c)=r2 \begin{aligned} ||r(t) - c|| = &\sqrt{(r(t) - c) \cdot (r(t) - c)} = r \\ &(r(t) - c) \cdot (r(t) - c) = r^2 \end{aligned}

The dot product has the same distributive and associative properties as the scalar multiplication, so you can do your maths the normal way:

(r(t)c)(r(t)c)=r2(e+tdc)(e+tdc)=r2ee+tedec+ted+t2ddtdccetcd+cc=r2 \begin{aligned} (r(t) - c) \cdot (r(t) - c) &= r^2 \\ (e+td- c) \cdot (e+td - c) &= r^2 \\ e\cdot e + te\cdot d - e\cdot c &\\ +t e \cdot d + t^2 d\cdot d - t d\cdot c &\\ -c \cdot e - t c \cdot d + c \cdot c &= r^2 \end{aligned}

The dot product is commutative (xy=yx)(x \cdot y = y \cdot x), and after rearranging the terms according to our free paramter tt we end up with:

t2dd+t2(edcd)+ee2(ec)+cc=r2t2dd+t2(edcd)+ee2(ec)+ccr2=0t2dd+t2d(ec)+(ec)(ec)r2=0 \begin{aligned} && t^2 d\cdot d & \\ &&+t 2(e \cdot d - c \cdot d) &\\ && +e\cdot e -2(e \cdot c) + c \cdot c &= r^2 \\ t^2 d\cdot d &+t 2(e \cdot d - c \cdot d) &+ e\cdot e -2(e \cdot c) + c \cdot c - r^2 &= 0 \\ t^2 d\cdot d &+t 2d\cdot(e - c) &+ (e-c) \cdot (e-c) - r^2 &= 0 \end{aligned}

Now we reduced it to a quadratic equation in tt. Quadratic equations of the form at2+bt+c=0at^2 + bt + c = 0 have the two solutions: t1,2=b±b24ac2at_{1,2} = \frac{-b\pm\sqrt{b^2-4ac}}{2a}

Applied to our equation we get:

t1,2=2d(ec)±(2d(ec))24dd((ec)(ec)r2)2ddt_{1,2} = \frac{-2d\cdot(e - c) \pm \sqrt{(2d\cdot(e - c))^2-4 d\cdot d ((e-c) \cdot (e-c) - r^2) }}{2 d\cdot d}

If and how many solutions exist depends on the discriminant, the term under the square root, D=(2d(ec))24dd((ec)(ec)r2)D =(2d\cdot(e - c))^2-4 d\cdot d ((e-c) \cdot (e-c) - r^2):

  • D<0D < 0: No solution. The ray misses the sphere
  • D=0D = 0: One solution. The ray touches the sphere in one point.
  • D>0D > 0: Two solution. The ray hits the sphere in two points, one entry hit and the exit hit on the other side

To get the intersections you calculate the tt values and then evaluate the ray r(t)r(t).

Left: No intersection Middle: One intersection Right: Two intersections Left: No intersection
Middle: One intersection
Right: Two intersections

Finally some C# code that implements the ray / sphere intersection test:

public RayTracer.HitInfo Intersect(Ray ray)
{
    RayTracer.HitInfo info = new RayTracer.HitInfo();

    Vector3 eMinusS = ray.origin - center;
    Vector3 d = ray.direction;
    double discriminant = Math.Pow(2 * Vector3.Dot(d, eMinusS), 2) - 4 * Vector3.Dot(d, d) *
                (Vector3.Dot(eMinusS, eMinusS) - Math.Pow(radius, 2.0f));

    if (discriminant < -Mathf.Epsilon)
    {   // 0 hits
        return info;
    }
    else {      // there will be one or two hits
        float front = -2.0f * Vector3.Dot(d, eMinusS);
        float denominator = 2.0f * Vector3.Dot(d, d);
        if (discriminant <= Mathf.Epsilon)
        {   // 1 hit
            info.time = (float)(front + Math.Sqrt(discriminant)) / denominator;  // does not matter if +- discriminant
        }
        else {  // 2 hits
            float t1 = (float)(front - Math.Sqrt(discriminant)) / denominator;  // smaller t value
            float t2 = (float)(front + Math.Sqrt(discriminant)) / denominator;  // larger t value
            if (t2 < 0) // sphere is "behind" start of ray
            {
                return info;    // no hit
            }
            else {  // one of them is in front
                if (t1 >= 0) info.time = t1; // return first intersection with sphere (usual case, smaller t)
                else info.time = t2;        // return second hit (ray's origin is inside the sphere)
            }
        }
    }

    // if we are here, info.time has been set, otherwise the function would have returned
    info.hitPoint = ray.GetPoint(info.time);
    info.normal = (info.hitPoint - center).normalized;
    return info;
}