Sphere tracing: the step that must not overshoot

A raymarched scene is a loop of four lines, and both of the ways it fails are invisible in the code. It renders holes where the surface should be solid, or it costs ten times what it should. Both come from the same quantity: how far the marcher is allowed to move per step.

Every number in this article is recomputed from its inputs by a script that runs on each build of this site. The code is from the tunnel scene in COLD START.

The loopPermalink to “The loop”#

f32 dist = 0.05f;
int hit = 0;
v3  p = ro;
for (int i = 0; i < 72; i++) {
    p = v3_add(ro, v3_scl(rd, dist));
    f32 d = tun_map(p, t);
    if (d < 0.0025f) { hit = 1; break; }
    dist += d * 0.80f;
    if (dist > 70.0f) break;
}

tun_map returns a signed distance: positive in free space, zero on the surface. The march steps along the ray by that value, evaluates again, and repeats. The name comes from the guarantee: if the returned value is the radius of a sphere centred on the current point that is known to be empty, then moving anywhere inside that sphere cannot cross the surface, so the step needs no check on what it passed through.

Three constants in those nine lines decide everything: the iteration cap of 72, the hit threshold of 0.0025, and the factor 0.80 that this article is mostly about.

The guarantee has a conditionPermalink to “The guarantee has a condition”#

The sphere argument holds only if the function underestimates. Formally, for the march to be safe the map $f$ must satisfy

$$|f(\mathbf{p})| \le \operatorname{dist}(\mathbf{p}, \partial\Omega)$$

everywhere, and the standard way to get that property is to require the field to be 1-Lipschitz — that is, $\|\nabla f\| \le 1$. A field whose gradient magnitude reaches $L$ somewhere can overestimate the distance by a factor of up to $L$ there, and the repair is to scale every step by $1/L$.

This is where hand-written distance functions go wrong. The exact distance to a sphere or a box has unit gradient by construction. Almost anything convenient does not: a field built from a radial expression, a domain warp, a cheap min of two shapes at a joint, or a non-uniform scale all inflate the gradient somewhere, and the march inherits the inflation.

The tunnel is not a distance fieldPermalink to “The tunnel is not a distance field”#

The tunnel wall is the zero set of

$$f(\mathbf{p}) = R(z,\theta) - r, \qquad r = \lVert (x,y) - \mathbf{a}(z) \rVert$$

with a bore radius that flutes six times around the circumference and ripples along the axis, around an axis $\mathbf{a}(z)$ that wanders:

$$\mathbf{a}(z) = \big(1.45\sin 0.062z,\; 1.15\cos 0.051z\big), \qquad R = 3.15 + 0.26\sin(6\theta + 0.34z) + 0.12\sin(1.65z - 1.9t)$$

Radially the field is exact: $\partial f/\partial r = -1$, which is the whole reason this form is used. The other two directions are where the cost sits.

The tangential term is the one that surprises people. Six flutes of amplitude 0.26 mean the radius changes by $6 \times 0.26 = 1.56$ per radian, and a radian at the wall is 3.15 units of arc, so the surface climbs at 0.4952 along the wall — half again as steep as anything the axis does.

Summing the three axial contributions gives 0.3937, and combining the three directions:

$$\lVert \nabla f \rVert \le \sqrt{1 + 0.4952^2 + 0.3937^2} = 1.1833$$

so the field overestimates the true distance by up to 18 per cent, and the largest step that is provably safe is

$$\lambda_{\max} = \frac{1}{1.1833} = 0.8451$$

The shipped code uses 0.80, which keeps 5.3 per cent in hand against a bound that assumed every term peaks at once. That is the entire justification for a constant that otherwise looks like it was tuned by eye — and it is why the comment beside it says the factor could be raised once the axis was made to bend more gently, rather than saying it was fiddled with until the holes went away.

What overshoot actually looks likePermalink to “What overshoot actually looks like”#

An unsafe step produces neither a warning nor a crash. It produces a ray that is inside the wall when it next evaluates the map, reads a negative distance, and — depending on how the hit test is written — either registers a hit at the wrong depth, so the shading normal is computed at a point behind the surface, or fails the test entirely and marches on into what should have been solid. On a moving camera the result is a shimmering rash of dark pixels that appears only at certain angles, which is why it is so often misdiagnosed as a precision problem.

The distinguishing symptom is that it worsens where the field is steepest rather than where the geometry is nearest, and that lowering the relaxation factor removes it while lowering the hit threshold does not.

Grazing rays and the step budgetPermalink to “Grazing rays and the step budget”#

The other failure is quiet. Take a ray approaching a surface at incidence angle $\alpha$, measured from the surface. Each step reduces the perpendicular distance by a factor $(1 - \lambda \sin\alpha)$, so from an initial separation $d_0$ the number of steps needed to reach the hit threshold $\varepsilon$ is

$$n = \frac{\ln(\varepsilon / d_0)}{\ln(1 - \lambda\sin\alpha)}$$

With $d_0 = 3.15$, $\varepsilon = 0.0025$ and $\lambda = 0.80$:

Head-on convergence is geometric and takes five steps. At five degrees it takes 99, and the loop is capped at 72. Setting $n = 72$ and solving for the angle gives the break-even incidence:

$$\alpha_{\min} = \arcsin\left(\frac{1 - (\varepsilon/d_0)^{1/72}}{\lambda}\right) = 6.78°$$

Below 6.78 degrees of incidence this marcher does not converge — it runs out of iterations while still short of the surface and reports a miss. That is what puts a dark band along the far wall of a corridor, and no amount of relaxing the step fixes it, because relaxation makes it worse: a smaller $\lambda$ needs more steps, not fewer.

The available fixes are to raise the cap where it is affordable, to widen $\varepsilon$ with distance so a surface two hundred units away is not being resolved to a quarter of a millimetre, or to accept the band and hide it in fog. This scene does the last of the three, which is a legitimate choice as long as it is a choice.

The normal costs six more evaluationsPermalink to “The normal costs six more evaluations”#

Once a hit is found the surface normal comes from the gradient, by central differences:

static v3 tun_normal(v3 p, f32 t)
{
    const f32 e = 0.010f;
    f32 dx = tun_map(V3(p.x + e, p.y, p.z), t) - tun_map(V3(p.x - e, p.y, p.z), t);
    f32 dy = tun_map(V3(p.x, p.y + e, p.z), t) - tun_map(V3(p.x, p.y - e, p.z), t);
    f32 dz = tun_map(V3(p.x, p.y, p.z + e), t) - tun_map(V3(p.x, p.y, p.z - e), t);
    return v3_norm(V3(dx, dy, dz));
}

Six map evaluations for one normal, on top of the march that found the point. The known alternative is the tetrahedron trick, which samples four vertices of a regular tetrahedron instead of six axis-aligned points and combines them with alternating signs, for a third less work at the same order of accuracy. It is not used here because the map is cheap relative to the march that precedes it, and because central differences make the epsilon’s effect obvious while tuning: at $e = 0.010$ the normal is averaged over a span of 0.020 units, which is what quietly smooths the fluting.

The choice of $e$ is a real trade. Too small, and the difference of two nearly equal floats loses its significant digits — the same cancellation that decides whether two runs of a solver agree. Too large, and the normal describes a surface smoother than the one that was hit, so the specular highlight slides across the geometry as the camera moves.

The cost, and why the render is offlinePermalink to “The cost, and why the render is offline”#

At 1920×1080 with two samples per pixel there are 4,147,200 primary rays in a frame. With the 72-step cap and six evaluations for each normal, one frame costs at most 298.6 million map evaluations for the march and 24.88 million for the normals; with the two-subframe motion blur the demo renders at, that is 647 million evaluations of tun_map for one frame of video.

Each of those evaluates two sines for the axis, a square root for the radius, two more transcendentals for the flute phase and one for the ripple. That is the whole argument for rendering offline rather than in real time: at 60 frames per second the budget is 16.7 milliseconds per frame, and this shot misses it by two orders of magnitude. Rendering slowly instead costs nothing but time, and buys motion blur and anti-aliasing that a real-time renderer has to fake.

It is also why the marcher is the only thing in the frame worth optimising, and why the map avoids an atan2.

Six-fold symmetry without a single transcendentalPermalink to “Six-fold symmetry without a single transcendental”#

The flute term needs $\sin(6\theta + \varphi)$, and $\theta$ only ever appears multiplied by six. Recovering the angle with atan2 and feeding it back into a sine pays for two transcendentals to compute something De Moivre’s formula gives directly: for a unit vector $(c_x, c_y) = (\cos\theta, \sin\theta)$,

$$(c_x + i c_y)^6 = \cos 6\theta + i \sin 6\theta$$

so squaring twice and multiplying once produces $\cos 6\theta$ and $\sin 6\theta$ in multiplications alone:

f32 inv = 1.0f / DM_MAX(r, 1e-5f);
f32 cx = q.x * inv, cy = q.y * inv;
f32 x2 = cx * cx - cy * cy, y2 = 2.0f * cx * cy;   /* z^2 */
f32 x4 = x2 * x2 - y2 * y2, y4 = 2.0f * x2 * y2;   /* z^4 */
f32 c6 = x4 * x2 - y4 * y2;                        /* cos 6t */
f32 s6 = x4 * y2 + y4 * x2;                        /* sin 6t */

The audit that runs on every build of this site checks the identity rather than trusting it: it evaluates both forms across the full turn and requires them to agree, which they can only do if the algebra is right. The unit-vector precondition is the part worth remembering — the identity is a statement about a point on the unit circle, so the division by $r$ is not an optimisation but the thing that makes the next six lines true.

What to check when a march misbehavesPermalink to “What to check when a march misbehaves”#

• Bound the gradient before choosing a step factor. Differentiate the map by hand, take the worst case, and use $1/L$. It is ten minutes of algebra, and it replaces a week of tuning by eye.
• Suspect the field, not the epsilon, when holes appear. Holes that move with the camera and vanish under a smaller step factor are overshoot; holes that stay put are a threshold or a budget.
• Count steps, not milliseconds, when it is slow. Writing the iteration count to the framebuffer as a heat map is the single most useful debug view a raymarcher has: it shows which regions burn the budget at grazing incidence, which no profiler will tell you.
• Keep the map free of state. The march calls it hundreds of millions of times per frame and across every thread at once, so anything cached inside it becomes a race, and anything conditional on frame history breaks the property that a scene is a pure function of time — which is what makes a preview and a final render agree.