Volumetric Shadow Raymarching with CSM in GLSL

This article details the algorithm for volumetric shadow raymarching through cascaded shadow map (CSM) frustums in GLSL. It covers the end-to-end rendering pipeline, including view-ray formulation, dynamic cascade selection, phase function evaluation, and numerical integration via transmittance accumulation.

The Algorithm: Cascaded Volumetric Raymarching

Volumetric shadowing—often referred to as god rays or crepuscular rays—is computed by numerically integrating light in-scattering along the camera’s view ray through a participating medium. When utilizing Cascaded Shadow Maps (CSM), the algorithm dynamically selects the correct depth slice along the ray to determine whether a given point in space is occluded.

The overall algorithm proceeds through four primary stages:

  1. Ray Setup and Bounds Determination: Identify the camera ray origin, direction, and ray march limits based on the scene's depth buffer.
  2. Ray Stepping and Jittering: Step along the ray using uniform, exponential, or equiangular distribution, typically offset by blue noise to prevent undersampling artifacts.
  3. Cascade Indexing and Visibility Evaluation: For each sample position, select the appropriate cascade layer, project the point into that cascade's light space, and query the shadow map.
  4. Scattering and Transmittance Accumulation: Apply the medium's phase function and compute light attenuation using the Beer-Lambert law.

Step 1: Ray Setup and Jittering

The view ray is defined from the camera position \(\mathbf{x}_0\) along the normalized ray direction \(\mathbf{v}\) up to a maximum distance \(t_{\text{max}}\), defined either by the scene geometry depth at that pixel or a fixed view distance:

\[t_{\text{max}} = \min(t_{\text{depth}}, t_{\text{far}})\]

To break up regular stepping artifacts (banding), apply an interleaved gradient noise or blue noise offset \(J \in [0, 1)\) to the sample positions:

float stepSize = t_max / float(NUM_STEPS);
float rayStartOffset = texture(u_BlueNoise, gl_FragCoord.xy / u_NoiseResolution).r;

Step 2: Dynamic Cascade Selection

As the ray marches away from the camera, the distance along the view vector increases. Because CSM partitions the view frustum into multiple depth intervals (cascades), each sample point's view-space depth \(z_{\text{view}}\) determines which shadow cascade must be sampled.

In GLSL, cascade selection is typically executed by comparing the sample point's linear view depth against an array of cascade split distances:

int getCascadeIndex(float linearDepth) {
    int cascadeIndex = NUM_CASES - 1;
    for (int i = 0; i < NUM_CASES - 1; ++i) {
        if (linearDepth < u_CascadeSplits[i]) {
            cascadeIndex = i;
            break;
        }
    }
    return cascadeIndex;
}

Step 3: Coordinate Transformation and Shadow Evaluation

Once the cascade index \(c\) is determined, the world-space sample point \(\mathbf{x} = \mathbf{x}_0 + t \cdot \mathbf{v}\) is transformed into that cascade's light-space projection coordinates:

vec4 lightSpacePos = u_LightViewProjMatrices[cascadeIndex] * vec4(samplePos, 1.0);
vec3 projCoords = lightSpacePos.xyz / lightSpacePos.w;
projCoords = projCoords * 0.5 + 0.5; // Transform to [0, 1] UV coordinates

The visibility factor \(V(\mathbf{x})\) is retrieved by comparing the projected depth against the depth recorded in the 2D texture array:

float shadowDepth = texture(u_CascadedShadowMaps, vec3(projCoords.xy, float(cascadeIndex))).r;
float visibility = (projCoords.z - u_ShadowBias <= shadowDepth) ? 1.0 : 0.0;

Step 4: Transmittance and Scattering Accumulation

At each step, light scattered toward the eye is weighted by the Henyey-Greenstein phase function \(P(\theta)\), which models directional scattering based on the angle between the view direction and the light direction \(\mathbf{L}\):

\[P(\theta) = \frac{1 - g^2}{4\pi (1 + g^2 - 2g \cos\theta)^{3/2}}\]

The total accumulated in-scattering \(L_{\text{in}}\) and transmittance \(T\) are updated iteratively via the Beer-Lambert law:

vec3 lightDir = normalize(u_SunDirection);
float cosTheta = dot(-rayDir, lightDir);
float g = 0.6; // Forward scattering coefficient
float phase = (1.0 - g * g) / (4.0 * 3.14159265 * pow(1.0 + g * g - 2.0 * g * cosTheta, 1.5));

vec3 accumulatedScattering = vec3(0.0);
float opticalDepth = 0.0;

for (int i = 0; i < NUM_STEPS; ++i) {
    float t = (float(i) + rayStartOffset) * stepSize;
    if (t >= t_max) break;

    vec3 samplePos = cameraPos + rayDir * t;
    float linearDepth = t * dot(rayDir, cameraForward);

    int cascadeIndex = getCascadeIndex(linearDepth);
    vec4 lightSpacePos = u_LightViewProjMatrices[cascadeIndex] * vec4(samplePos, 1.0);
    vec3 projCoords = (lightSpacePos.xyz / lightSpacePos.w) * 0.5 + 0.5;

    float shadowDepth = texture(u_CascadedShadowMaps, vec3(projCoords.xy, float(cascadeIndex))).r;
    float visibility = (projCoords.z - u_ShadowBias <= shadowDepth) ? 1.0 : 0.0;

    if (visibility > 0.0) {
        float stepOpticalDepth = u_AbsorptionExtinction * stepSize;
        opticalDepth += stepOpticalDepth;
        float stepTransmittance = exp(-opticalDepth);

        accumulatedScattering += u_LightColor * phase * stepTransmittance * stepSize;
    }
}

The resulting accumulatedScattering value is added to the scene's color buffer during composite or tone-mapping passes.