Stochastic Raymarched Soft Shadows in GLSL Compute

This article examines the implementation and evaluation of stochastic ray distribution for compute-based raymarched soft shadows using GLSL. It covers the mathematical formulation of jittered ray directions toward area lights, low-discrepancy sampling techniques, the integration of blue noise to minimize visible clustering, and the raymarching execution model within compute shaders. Finally, it addresses divergence management and spatial filtering strategies required to reconstruct clean penumbras from low-sample-count stochastic distributions.

Area Light Formulation and Ray Jittering

Standard raymarched soft shadows often approximate penumbras using distance-field estimators along a single unjittered ray directed at a point light. Stochastic soft shadows, by contrast, physically model area light sources (such as disks, spheres, or rectangles) by distributing secondary shadow rays across the solid angle subtended by the light.

To evaluate a spherical light source with radius \(R\) at position \(\mathbf{P}_L\) from a surface hit point \(\mathbf{x}\), the target sample point \(\mathbf{S}\) on the light is determined stochastically:

  1. Generate a two-dimensional pseudo-random or quasi-random sample \(\mathbf{\xi} \in [0, 1)^2\).
  2. Map \(\mathbf{\xi}\) onto the surface of a disk oriented toward \(\mathbf{x}\), using concentric disk mapping to avoid areal distortion.
  3. Construct an orthonormal basis (tangent, bitangent, normal) aligned with the vector \(\mathbf{L} = \mathbf{P}_L - \mathbf{x}\).
  4. Offset the center \(\mathbf{P}_L\) along the tangent and bitangent axes by the mapped disk coordinates scaled by \(R\).
  5. Normalize the resulting vector \(\mathbf{D} = \mathrm{normalize}(\mathbf{S} - \mathbf{x})\) to form the primary direction for the shadow raymarch.

Low-Discrepancy and Blue Noise Sampling in GLSL

Evaluating stochastic rays using naive white noise introduces high-frequency, uncorrelated grain that spatial and temporal filters struggle to reconstruct. Modern compute-based pipelines evaluate ray distribution using spatio-temporal blue noise or low-discrepancy sequences like Halton or Sobol offsets.

Blue noise pushes sampling errors into higher frequencies, making the resulting penumbra noise significantly easier to eliminate with a lightweight bilateral blur. In GLSL, this is implemented by sampling a precomputed, tileable blue noise texture indexable by screen-space coordinates and temporal frame index:

// Fetch blue noise offset per pixel
vec2 blueNoise = texelFetch(u_BlueNoiseTex, (ivec2(gl_GlobalInvocationID.xy) + ivec2(u_FrameOffset * 17)) % 128, 0).rg;

// Uniform concentric disk sampling
vec2 diskSample = concentricSampleDisk(blueNoise);

// Build orthonormal basis
vec3 N_L = normalize(u_LightPos - worldPos);
vec3 up = abs(N_L.z) < 0.999 ? vec3(0.0, 0.0, 1.0) : vec3(1.0, 0.0, 0.0);
vec3 T_L = normalize(cross(up, N_L));
vec3 B_L = cross(N_L, T_L);

// Final ray direction toward jittered area light point
vec3 targetLightPos = u_LightPos + (T_L * diskSample.x + B_L * diskSample.y) * u_LightRadius;
vec3 shadowRayDir = normalize(targetLightPos - worldPos);

Compute Raymarching Evaluation

The compute shader executes the raymarch through a signed distance field (SDF) or a voxel grid. Because the ray direction is already stochastically distributed across the light surface, binary visibility (0.0 for occluded, 1.0 for unoccluded) can be evaluated directly rather than tracking step-wise distance heuristics:

layout(local_size_x = 8, local_size_y = 8, local_size_z = 1) in;

float evaluateShadowRay(vec3 ro, vec3 rd, float maxDist) {
    float t = 0.05; // Offset to avoid self-shadowing
    for (int i = 0; i < MAX_STEPS && t < maxDist; ++i) {
        float d = mapSceneSDF(ro + rd * t);
        if (d < 0.001) {
            return 0.0; // Fully occluded
        }
        t += max(d, 0.01); // March forward by safe distance
    }
    return 1.0; // Unoccluded
}

By casting only 1 to 4 samples per pixel (SPP) per frame, evaluation remains performant across complex distance fields.

SIMD Divergence and Compute Optimizations

Compute shaders execute threads in warps (NVIDIA) or wavefronts (AMD). Stochastic raymarching risks performance degradation due to branch divergence: neighboring threads tracing slightly different ray paths will step different numbers of times, causing thread serialization.

To mitigate SIMD divergence:

  • Depth and Normal Rejection: Pixels marked as skybox, background, or facing away from the light source (\(\mathbf{N} \cdot \mathbf{L} \le 0\)) skip raymarching entirely using an early branch before calculating stochastic offsets.
  • Hierarchical Step Snapping: Clamping minimum step sizes prevents paths from getting stuck in geometric crevices that only a fraction of the warp encounters.
  • Subgroup Uniformity: Where available, compute shaders can leverage GL_KHR_shader_subgroup_ballot to dynamically prune threads that have already found occluders, though standard fixed-step loops with early break statements are often sufficient if the divergence latency is hidden by global memory latency.

Denoising Integration

Because low-SPP stochastic evaluations output noisy 1-bit or low-precision visibility masks, an edge-aware filtering pass is mandatory. The compute pipeline outputs the stochastic shadow buffer to a dedicated render target, which is then processed by:

  1. Temporal Accumulation (TAA/TSR): Reprojecting the previous frame's shadow mask via motion vectors to blend samples over time, using exponential moving averages.
  2. Spatial Denoising: An edge-avoiding à-trous wavelet filter or cross-bilateral filter that uses scene normals and linear depth buffers to blur the stochastic noise across the penumbra without bleeding shadows across geometric boundaries.