Ray-sphere intersection

From DmWiki

To find the intersection of a sphere and a ray, you will need to begin with the following information.

  • Sphere with radius r positioned at point c.
  • Ray (line segment) beginning at a point p and travelling in direction d.

Algebraic Method

The most straightforward way to intersect a ray with a sphere is to do it algebraically. You can express a sphere in implicit form like this:

(\mathbf{x} - \mathbf{c}) \cdot (\mathbf{x} - \mathbf{c}) = r^2

This equation is satisfied when the point x lies exactly on the sphere. In order to intersect with a ray, use the parametric equation for the ray:

\mathbf{x}(t) = \mathbf{o} + t \mathbf{d}

where, for rays, we require t \geq 0. Then we simply set the x in the first equation equal to the x in the second equation.

(\mathbf{o} + t\mathbf{d} - \mathbf{c}) \cdot (\mathbf{o} + t\mathbf{d} - \mathbf{c}) = r^2

Solving for t, this expands into a quadratic equation: At^2 + Bt + C = 0\, where

A = \mathbf{d} \cdot \mathbf{d},
B = 2(\mathbf{o} - \mathbf{c}) \cdot \mathbf{d},
C = (\mathbf{o} - \mathbf{c}) \cdot (\mathbf{o} - \mathbf{c}) - r^2

This can be solved by the usual means, i.e. the quadratic formula.

This method can be highly optimized; a very fast ray-sphere intersection routine (due to davepermen (http://www.devmaster.net/forums/index.php?showuser=134)) is below. Note that it assumes d is a unit vector.

float intersectRaySphere(const Ray &ray, const Sphere &sphere) {
	Vec dst = ray.o - sphere.o;
	float B = dot(dst, ray.d);
	float C = dot(dst, dst) - sphere.r2;
	float D = B*B - C;
	return D > 0 ? -B - sqrt(D) : std::numeric_limits<float>::infinity();

Geometric Method

The primary advantage of the geometric ray-sphere intersection routine is that it can detect lack of intersection earlier than the algebraic routine. It can also deal with rays that have d not unit length without the need of an additional square root calculation, though if the same ray is used for more than two or thee intersection computations it is still faster to normalize d in advance. Depending on the relative speed of comparisons and mathematical operations on the underlying hardware, the geometric method is not necessarily faster than the algebraic method, though in general it does perform better.

First, we find a vector from the ray origin to the center of the sphere,

\mathbf{o_c} = \mathbf{c} - \mathbf{p}

and compute the square of its length

l^2_{oc} = \mathbf{o_c} \cdot \mathbf{o_c}.

If l^2_{oc} is greater than r2, the the ray originates outside the sphere. We now find the parametric t value that corresponds to the ray's closest approach to the center of the sphere,

t_{ca} = \frac{\mathbf{o_c} \cdot \mathbf{d}}{\mathbf{d} \cdot \mathbf{d}}.

(note: if d is a unit vector, the division by d · d can be omitted.)

If tca is less than zero, the sphere center lies behind the ray, which means we have no intersection (or the camera is inside the sphere). Otherwise, we can use the pythagorean theorem to find the length of the half cord defined by the ray and the sphere:

l^2_{hc} = \frac{r^2 - l^2_{oc}}{\mathbf{d}\cdot\mathbf{d}} + {t_{ca}}^2

(note: again, we can leave out the division by d · d if we know d is a unit vector)

If l^2_{hc} is less than zero, the ray misses the sphere. Otherwise the t of intersection is

t = t_{ca} - \sqrt{l^2_{hc}}.

If l^2_{oc} was less than r2 the ray originated inside the sphere. Depending on the application, this case can be treated in several ways; if the point where the ray exited the far side of the sphere is desired, we compute tca and l^2_{hc} as above, though we skip the check on tca < 0, and compute t as

t = t_{ca} + \sqrt{l^2_{hc}}.

A full implementation of the geometric ray-sphere intersection is given by

float intersectRaySphere(const Ray &ray, const Sphere &sphere) {
	Vec oc = sphere.c - ray.p;
	float l2oc = dot(oc,oc);
	if (l20c < sphere.r2) { // starts inside of the sphere
		float tca = dot(oc, ray.d) / dot(ray.d, ray.d);			// omit division if ray.d is normalized
		float l2hc = (sphere.r2 - l20c) / dot(ray.d, ray.d) + tca*tca;  // division
		return tca + sqrt(l2hc);
	} else {
		float tca = dot(oc, ray.d);
		if (tca < 0) // points away from the sphere
			return std::numeric_limits<float>::infinity();
		float l2hc = (sphere.r2 - l20c)/dot(ray.d, ray.d) + (tca*tca);	// division
		return l2hc > 0 ? tca - sqrt(l2hc) : std::numeric_limits<float>::infinity();

DevMaster navigation