SDF Raymarching Soft Shadows in GLSL

This article details the mathematical formulation and GLSL implementation used to calculate soft shadow penumbras in raymarched Signed Distance Fields (SDFs). It covers the geometric derivation introduced by John Hart and popularized by Inigo Quilez, discusses the relationship between surface distance and ray length, and provides both the classic penumbra formula and the improved version that eliminates step-size shadow artifacts.

The Mathematical Foundation

In standard raymarching, a ray is cast from a surface point toward a light source to determine visibility. While a hard shadow only checks whether a ray intersects an occluder, a soft shadow estimates the light source's partial occlusion (the penumbra) by evaluating how closely the ray misses surrounding geometry.

At any distance \(t\) along the light ray, the signed distance field provides the distance \(d(t)\) to the nearest obstacle. This value represents the radius of an empty sphere centered at the ray's current sample point. The maximum angle subtended by this empty sphere from the light ray's origin can be approximated using the ratio of the distance to the surface \(d(t)\) over the distance traveled along the ray \(t\):

\[\tan(\theta) \approx \frac{d(t)}{t}\]

To compute the overall shadow intensity \(S\), we track the minimum ratio observed along the entire ray path from \(t_{\min}\) to \(t_{\max}\):

\[S = \min_{t \in [t_{\min}, t_{\max}]} \left( \frac{k \cdot d(t)}{t} \right)\]

Where:

• \(d(t)\) is the evaluated SDF distance at distance \(t\) along the ray.
• \(t\) is the current distance traveled from the surface origin toward the light source.
• \(k\) is a user-defined scaling factor that controls penumbra sharpness. Larger values of \(k\) produce sharper shadows, while smaller values produce wider, softer penumbras.
• \(S\) is clamped to the range \([0.0, 1.0]\), where \(0.0\) represents full shadow (umbra) and \(1.0\) represents fully lit.

Classic GLSL Implementation

The classic approach evaluates the ratio directly within the sphere-tracing loop:

float calcSoftshadow(in vec3 ro, in vec3 rd, in float mint, in float maxt, in float k) {
    float res = 1.0;
    float t = mint;
    
    for(int i = 0; i < 64 && t < maxt; i++) {
        float h = map(ro + rd * t);
        
        // Ray has hit an occluder
        if(h < 0.001) {
            return 0.0;
        }
        
        // Compute penumbra ratio and update minimum
        res = min(res, k * h / t);
        t += h;
    }
    
    return clamp(res, 0.0, 1.0);
}

The Improved Formulation (Artifact-Free)

The classic formula can introduce visual banding and premature darkening when the sample step skips past the true minimum distance point. Inigo Quilez addressed this issue by treating the current and previous sample spheres geometrically to analytically estimate the closest approach of the ray to the occluder between consecutive steps.

Given the current sphere radius \(d\) and previous sphere radius \(d_{\text{prev}}\), the distance \(y\) from the current sample point back along the ray to the actual point of closest approach is:

\[y = \frac{d^2}{2 \cdot d_{\text{prev}}}\]

The corrected minimum distance \(d_{\text{perp}}\) to the surface along the ray segment between these two samples is derived via the Pythagorean theorem:

\[d_{\text{perp}} = \sqrt{d^2 - y^2}\]

The improved penumbra factor evaluation substitutes \(d_{\text{perp}}\) and the corrected distance along the ray \((t - y)\):

\[S = \min \left( S, \frac{k \cdot \sqrt{d^2 - y^2}}{t - y} \right)\]

Improved GLSL Implementation

float calcSoftshadow(in vec3 ro, in vec3 rd, in float mint, in float maxt, in float k) {
    float res = 1.0;
    float t = mint;
    float ph = 1e10; // Previous distance step
    
    for(int i = 0; i < 64 && t < maxt; i++) {
        float h = map(ro + rd * t);
        
        if(h < 0.001) {
            return 0.0;
        }
        
        // Geometric correction for closest approach
        float y = (i == 0) ? 0.0 : (h * h) / (2.0 * ph);
        float d = sqrt(max(0.0, h * h - y * y));
        
        res = min(res, k * d / max(0.0, t - y));
        ph = h;
        t += h;
    }
    
    return clamp(res, 0.0, 1.0);
}

This improved version preserves continuous penumbra gradients and maintains realistic softness regardless of step size variations.