Ray Tracing in a Weekend Part 4 - Adding a Sphere


Ray vs. Sphere

Part 3 introduced a ray class. Now it’s time to let the ray hit the first object: a sphere.

The equation of a sphere centered at the origin is [1]

\[ x^2 + y^2 + z^2 = R^2 \]

which is essentially derived from Pythagoras’ Theorem extended to three dimensions. [2] For any \( (x,y,z) \), if \( x^2 + y^2 + z^2 = R^2 \) then \( (x,y,z) \) is on the sphere, otherwise it’s not.

The more general equation of a sphere centered at a point \( (x_0,y_0,z_0) \) is

\[ (x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2 = R^2 \]

This can also be expressed as a dot-product:

$$\begin{pmatrix} x \\ y \\ z \end{pmatrix} - \begin{pmatrix} x_0 \\ y_0 \\ z_0 \end{pmatrix} = \begin{pmatrix} x-x_0 \\ y-y_0 \\ z-z_0 \end{pmatrix} $$ $$\begin{pmatrix} x-x_0 \\ y-y_0 \\ z-z_0 \end{pmatrix} \cdot \begin{pmatrix} x-x_0 \\ y-y_0 \\ z-z_0 \end{pmatrix} = (x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2 = R^2 $$

For center \( C = (x0,y0,z0) \) and a point \( P = (x,y,z) \), we can also write:

$$\begin{pmatrix} P - C \end{pmatrix} \cdot \begin{pmatrix} P - C \end{pmatrix} = R^2 $$

The ray \( p(t) = A + t*B \) intersects the sphere where \( p(t) = P \). If the ray hits the sphere, there is a t for which p(t) satisfies the sphere equation.

$$\begin{pmatrix} p(t) - C \end{pmatrix} \cdot \begin{pmatrix} p(t) - C \end{pmatrix} = R^2 $$ $$\begin{pmatrix} A + t*B - C \end{pmatrix} \cdot \begin{pmatrix} A + t*B - C \end{pmatrix} = R^2 $$

We are looking for t, and A, B, C and R are constants. We can rearrange the equation as follows:

$$\begin{pmatrix} t*B + (A-C) \end{pmatrix} \cdot \begin{pmatrix} t*B + (A-C) \end{pmatrix} = R^2 \hspace{1cm} \vert \space \text{let }oc = (A-C) $$ $$ t*B \cdot t*B + 2 * t*B \cdot oc + oc \cdot oc = R^2 \hspace{1cm} \vert \space \text{let }a = (B \cdot B) $$ $$ a * t^2 + 2 * t * (B \cdot oc) + oc \cdot oc - R^2 = 0 \hspace{1cm} \vert \space \text{let }b = 2 * (B \cdot oc) $$ $$ a * t^2 + b * t + oc \cdot oc - R^2 = 0 \hspace{1cm} \vert \space \text{let }c = (oc \cdot oc) - R^2 $$ $$ a * t^2 + b * t + c = 0 \hspace{1cm} $$

This is a quadratic equation that can be solved for t with the Quadratic Formula:

$$ t = \frac{-b \pm \sqrt{b^2 - 4ac} }{2a} $$

The square root portion of this equation can either be positive (meaning there are two real solutions), or zero (meaning there’s one real solution), or negative (meaning there are no real solutions).

This relates directly to the geometry. If there are two real solutions, the ray intersects the squere at two points. If there’s exactly one real solution, it touches the sphere at one point on the surface. If there are no real solutions, then the ray doesn’t touch the sphere at all.

For the hit_sphere() function we only care about the root and whether it’s greater than 0:

=== sphere-example.cpp

bool hit_sphere(const vec3& center, float radius, const ray& r) {
  vec3 oc = r.origin() - center;
  float a = dot(r.direction(), r.direction());
  float b = 2.0 * dot(oc, r.direction());
  float c = dot(oc, oc) - radius * radius;
  float discriminant = b*b - 4*a*c; // for quadratic equation
  return (discriminant > 0); // all we care about is whether the root is not 0
}

vec3 color(const ray& r) {
  if (hit_sphere(vec3(0,0,-1), 0.5, r)) return vec3(1,0,0); // just return red if sphere was hit

  vec3 unit_direction = unit_vector(r.direction());
  float t = 0.5f * (unit_direction.y() + 1.0f);
  return (1.0f - t) * vec3(1.0f, 1.0f, 1.0f) + t * vec3(0.5f, 0.7f, 1.0f);
}

Javascript and 2D Canvas

Here’s also the same example in Javascript.

var example = {
  hit_sphere: function(center, radius, r) {
    oc = r.origin().sub(center);
    var a = r.direction().dot(r.direction());
    var b = 2.0 * oc.dot(r.direction());
    var c = oc.dot(oc) - radius * radius;
    var discriminant = b*b - 4*a*c; // for quadratic equation
    return (discriminant > 0); // all we care about at this point is whether the root is not 0
  },
  color: function(r) {
    if (this.hit_sphere(new vec3(0,0,-1), 0.5, r)) return new vec3(1,0,0); // just return red if sphere was hit

    var unit_direction = r.direction().unit_vector();
    var t = 0.5 * (unit_direction.y + 1.0);
    return (new vec3(1.0, 1.0, 1.0))
      .mul(1.0 - t)
      .add( 
        (new vec3(0.5, 0.7, 1.0)).mul(t) 
      );
  },
  // ...
};

Please check my git-repository at https://github.com/celeph/ray-tracing-in-a-weekend for the complete code, and Peter Shirley’s excellent books and his original code and book repository.