Solving Quadratic Ray Intersections Without Losing the Far Root

A sphere is the cheapest curved primitive a ray tracer can support: three floats for the centre, one for the radius, and a quadratic to solve. The quadratic is where the trouble starts. The formula every algebra course teaches has a catastrophic cancellation for exactly the configurations that a renderer meets constantly — a ray shot at a far-away object, where the two roots are close in magnitude and their difference is what carries the useful information.

The quadratic that a sphere produces

Substituting $o + t d$ into $|p - c|^2 = r^2$ and writing $f = o - c$ gives

$$(d \cdot d),t^2 + 2 (f \cdot d),t + (f \cdot f - r^2) = 0$$

so with normalised $d$ the coefficients are $a = 1$, $b = 2(f \cdot d)$, $c = f \cdot f - r^2$. The discriminant is $b^2 - 4ac$; its sign partitions the ray’s line into the three cases that matter. Positive means two distinct roots and the line passes through the sphere. Zero means a tangent, where the hit has zero extent and shading is undefined because the geometric normal depends on a direction that is degenerate. Negative means the line misses, and the renderer moves on. Only the sign of the discriminant is needed for that decision, and the sign survives rounding far better than the value does — a discriminant that computes as $-10^{-9}$ when the true value is zero is a tangent treated as a miss, which is visually identical.

The remaining work is extracting the roots. With the doubled form, the roots are

$$t = \frac{-b/2 \pm \sqrt{(b/2)^2 - ac}}{a}, \qquad \text{here} \qquad t = -(f \cdot d) \pm \sqrt{(f \cdot d)^2 - (f \cdot f - r^2)}.$$

Call $h = f \cdot d$ and $k = f \cdot f - r^2$. Then $t = -h \pm \sqrt{h^2 - k}$, and everything problematic is visible in one expression. For a sphere far away, $f$ is large, so $h$ is large, so $h^2$ is enormous while $k = |f|^2 - r^2$ is only large. The discriminant $h^2 - k$ is a difference of nearly equal magnitudes whenever the ray nearly grazes — the grazing case being the silhouette, which every visible sphere has. The discriminant loses roughly half its significant digits relative to $h^2$, because the subtraction of two numbers of magnitude $h^2$ retains only the digits where they differ.

Worse, the subsequent $\pm$ step forms one root by adding $\sqrt{h^2-k}$, which is close to $|h|$, to a term of magnitude $|h|$. If $h > 0$, the root $-h + \sqrt{h^2-k}$ is a difference of nearly equal magnitudes, and if the sphere is small and far the true root can be orders of magnitude smaller than either operand. The absolute error in the square root, which is at the level of a few ulps of $h^2$, is then comparable to the root itself. The computed near root can be off by a fractional amount, or even come out negative when the true value is positive.

The stable form and why it is only one addition longer

The fix is the textbook one, and it always surprises people that it is not what is taught: for the root formed by a difference, do not compute it directly. Instead compute the root formed by the sum, which has no cancellation, and recover the other root from the product of the roots, $c/a$ (or $t_1 t_2 = (|f|^2 - r^2) = k$ in the normalised case).

#include <math.h>

typedef struct { float x, y, z; } vec3;
static inline float dot3(vec3 a, vec3 b) { return a.x*b.x + a.y*b.y + a.z*b.z; }

/* Nearest positive hit of normalized ray (o, d) with sphere (c, r). */
int ray_sphere(vec3 o, vec3 d, vec3 c, float r, float t_min, float t_max, float *t)
{
    vec3 f = (vec3){o.x - c.x, o.y - c.y, o.z - c.z};
    float h = dot3(f, d);                       /* half the linear coefficient */
    float k = dot3(f, f) - r * r;               /* constant term, a == 1 */
    float disc = h * h - k;
    if (disc < 0.0f) return 0;                  /* line misses entirely */
    float s = sqrtf(disc);

    /* Form the root that adds magnitudes; no cancellation there. */
    float t0, t1;
    if (h > 0.0f) { t1 = -h - s; t0 = k / t1; } /* t1 is definitely the far root */
    else          { t0 = -h + s; t1 = k / t0; } /* t0 is definitely the near root */
    if (t0 > t_min && t0 < t_max) { *t = t0; return 1; }
    if (t1 > t_min && t1 < t_max) { *t = t1; return 1; }
    return 0;                                   /* both behind, or past t_max */
}

The branch on the sign of h selects the root whose computation is a sum rather than a difference. In both branches the recovered root comes from a division by a quantity of comparable magnitude, so the relative error stays near a few ulps instead of degrading with the ratio of the roots. The extra work is one comparison and one divide; the divide is the expensive part, but it is also the operation that keeps the result meaningful, and it costs less than the wrong answer does.

Two details in that code are not obvious from the algebra. k / t1 can overflow or produce a non-finite value when t1 is very small in magnitude, which happens when the origin is almost exactly on the sphere’s surface; the guard against t1 == 0 is a guard against division by zero, and it belongs before the division, not after checking the result for infinity. And the ordering of the two checks matters only for rays that hit both roots inside the interval, which happens when the origin is inside the sphere: then t0 is negative and t1 positive, so the near-root test fails on t_min and the far root is returned. A renderer that needs the exit point of a solid sphere uses the same code with the roles reversed.

What the discriminant sign hides

A more insidious failure is the one the discriminant cannot report. When $h^2 - k$ is computed as a small positive number from operands of much larger magnitude, its value is dominated by rounding noise, so a computed $+10^{-9}$ can correspond to a true value of $-10^{-9}$. The result is a spurious pair of roots, both of which land within a hair of the silhouette point, so the sphere gains a ring of pixels whose hit points are essentially random. Nothing in the arithmetic signals this; the test simply reports a hit at a slightly wrong place.

The standard response is to compare the discriminant against a tolerance scaled by the magnitudes involved rather than against zero, and to treat the ambiguous band as either a tangent or a miss consistently. Consistency is what matters: if the primary-ray test says hit and the shadow-ray test says miss near the silhouette, the pixel has a lit rim. A tolerance of a few ulps of $h^2$, applied identically in every code path that evaluates the same geometry, keeps primary and secondary rays agreeing.

The roots themselves also need an interval test, and this is where spheres differ from triangles. A ray that starts inside a sphere has a negative near root, and a renderer must decide whether the surface counts as hit from the inside. For opaque surfaces the near root is discarded on $t_{\min}$ and the far one returned as a back-face hit; for transparent surfaces both roots are needed and the ray continues inside the medium. That decision is the subject of ray parameter intervals and the epsilon offset problem, because the interval is not just a filter but the mechanism that prevents a ray from immediately re-hitting the surface it left.

General quadrics and the same cancellation

A sphere is one member of a family. Ellipsoids, cylinders, cones, and paraboloids are all quadrics, and all of them produce the same three coefficients from a general implicit form. The difference is that the coefficients are derived from a symmetric matrix rather than from a dot product, so they take more arithmetic, but the conditioning analysis is identical: whenever the primitive’s silhouette is in view, the discriminant is a difference of large nearly equal terms, and the same root-recovery trick applies, with the product of roots replaced by $c/a$ in the general case rather than $k$.

The place where general quadrics diverge is bounding. A sphere has a trivial bound — the same three floats plus a radius — and a box that encloses it, so it can be dropped into the slab-tested hierarchy with a bounding box and no special cases. An arbitrary quadric can be unbounded within any box large enough to contain its visible part (a plane is an unbounded quadric; so is a cone), which forces the renderer either to clip the primitive to a region explicitly or to treat it as a planar approximation inside a bound. Analytic quadrics are therefore usually reserved for primitives that are boundable: spheres, axis-aligned cylinders, and ellipsoids with a stored inverse transform.

There is a further numeric hazard specific to transformed quadrics. An ellipsoid is often stored as a unit sphere plus a 4×4 matrix that scales and rotates it, and the natural implementation transforms the ray into the sphere’s local space, solves, and then maps the parameters back. That mapping assumes the transform’s inverse is well conditioned. A matrix that squashes one axis by a large factor has an inverse that stretches it by the same factor, and a ray direction transformed by that inverse can lose most of its precision in the squashed component. Solving in world space with the full quadratic form avoids the inverse but costs more per test; solving in local space costs less but pushes the error into the transform. The choice has to be made per primitive type, and the two paths must not be mixed for the same surface, because their t values differ in the last bits and a shadow ray tested through the other path will not agree with the primary one at the silhouette.

Degenerate inputs that the algebra does not exclude

$a$ can be zero. That happens when the ray direction is a zero vector, which a renderer should never produce but a transform can: a direction scaled by a matrix with a zero row comes out zero, and then the quadratic is linear, with a single root $t = -c/b$. The formula above divides by $a$ implicitly by assuming it is 1, so a degenerate direction silently becomes a “sphere at infinity” and either misses or returns garbage. A guard that rejects near-zero direction length at the point where the ray is created is cheaper than a guard in every intersection routine.

$c = 0$ is not degenerate: it means the origin lies exactly on the sphere, and one root is zero while the other is $-2h$. With a strictly positive $t_{\min}$ the zero root is filtered, which is the desired behaviour for a ray leaving a surface. $k/t_1$ in that case yields exactly zero, so the recovery is well behaved.

The radius can be zero, giving a point-sphere whose discriminant is $(f \cdot d)^2 - f \cdot f$, which for a normalised direction is negative unless the origin lies exactly on the point. That is correct behaviour under the algebra and the reason point lights are not modelled as zero-radius spheres: an intersection test against a point has no surface to shade.

Finally, the ordering of the two roots is part of the contract, not an implementation detail. Callers of an intersection routine assume that a reported hit is the nearest one along the ray within the interval, and code that walks a list of primitives relies on that to keep the running t_max shrinking. The stable form preserves the ordering, since it returns t0 when it is in range and t1 otherwise, and the sign-based branch determines which is which before any comparison. A version that computed both roots with the naive formula and then sorted them numerically would still be ordering-correct and numerically wrong; the two properties are independent, and only one of them is repairable by sorting.