Variance Shadow Mapping with Chebyshev Bounds in GLSL

Variance Shadow Mapping (VSM) resolves standard shadow map aliasing and filtering limitations by treating depth distribution as a statistical probability problem rather than a binary depth test. By storing the first and second moments of depth in a texture, an application can filter shadow maps using standard hardware features like bilinear filtering, mipmapping, and Gaussian blurs. Inside a GLSL fragment shader, Chebyshev's inequality uses these stored moments to calculate an upper bound on the probability that a given surface fragment is occluded, producing smooth, continuous shadow transitions.

Storing Depth Moments

Traditional shadow mapping stores only the closest light-space depth (\(z\)), making linear filtering mathematically invalid because the depth comparison is non-linear. VSM overcomes this by writing two values into a floating-point texture (such as RG32F or RG16F):

  1. The light-space depth: \(z\)
  2. The squared light-space depth: \(z^2\)

When this texture is filtered across a footprint using bilinear filtering or blurring passes, the sampled values represent the expected depth \(E[x]\) and the expected squared depth \(E[x^2]\) of the occluders covering that region.

The Mathematics of Cantelli's (Chebyshev's) Inequality

In the receiver's fragment shader, the current fragment has a light-space depth \(t\). If \(t \le E[x]\), the fragment is closer to the light than the average occluder, meaning it is completely unshadowed.

When \(t > E[x]\), the fragment might be in shadow. VSM approximates the non-shadowed percentage using a one-sided variant of Chebyshev's inequality, also known as Cantelli's inequality. It states that for any distribution with mean \(\mu\) and variance \(\sigma^2\):

\[P(x \ge t) \le p_{\text{max}}(t) = \frac{\sigma^2}{\sigma^2 + (t - \mu)^2}\]

In this equation:

  • The mean \(\mu = E[x]\) is retrieved directly from the texture's red component.
  • The variance \(\sigma^2 = E[x^2] - (E[x])^2\) is computed using both the red and green components.
  • The distance \((t - \mu)\) represents the difference between the receiver depth and the mean occluder depth.
  • \(p_{\text{max}}(t)\) serves as an upper bound approximation of visibility (the lit factor).

GLSL Implementation

The following GLSL function demonstrates the computation of the visibility factor inside a fragment shader:

float ChebyshevUpperBound(vec2 moments, float receiverDepth) {
    // Surface is fully lit if the receiver depth is before the occluder depth
    if (receiverDepth <= moments.x) {
        return 1.0;
    }

    // Compute variance: Var(X) = E(X^2) - E(X)^2
    float variance = moments.y - (moments.x * moments.x);

    // Prevent numerical instability and reduce shadow acne
    float minVariance = 0.00002;
    variance = max(variance, minVariance);

    // Compute distance to mean
    float d = receiverDepth - moments.x;

    // Apply Cantelli's inequality
    float pMax = variance / (variance + d * d);

    return pMax;
}

// Fragment shader execution
uniform sampler2D shadowMap;
in vec4 vLightSpacePosition;

void main() {
    vec3 projCoords = vLightSpacePosition.xyz / vLightSpacePosition.w;
    projCoords = projCoords * 0.5 + 0.5; // Map to [0, 1] range

    vec2 moments = texture(shadowMap, projCoords.xy).rg;
    float visibility = ChebyshevUpperBound(moments, projCoords.z);

    // Apply visibility to direct lighting
    vec3 color = ambientColor + directColor * visibility;
    gl_FragColor = vec4(color, 1.0);
}

Mitigation of Light Bleeding

Because Chebyshev's inequality gives an upper bound rather than the exact probability, complex geometry with multiple overlapping occluders can cause the computed variance to become excessively large. When \(\sigma^2\) is large, \(p_{\text{max}}\) remains high even when \(t \gg \mu\), creating an artifact known as "light bleeding," where shadows erroneously fade out.

To reduce light bleeding in GLSL, developers apply a non-linear remap function to clamp low probabilities to zero:

float ReduceLightBleeding(float pMax, float minThreshold) {
    // Rescale pMax from [minThreshold, 1.0] to [0.0, 1.0]
    return clamp((pMax - minThreshold) / (1.0 - minThreshold), 0.0, 1.0);
}

By tweaking the threshold value (typically between \(0.1\) and \(0.3\)), light bleeding is suppressed in dense geometry while retaining soft shadow borders elsewhere.