ESM Exponential Overflow Mitigation in GLSL

Exponential Shadow Maps (ESM) replace standard shadow map depth comparisons with an exponential test to allow linear filtering, but high exponential factors frequently trigger 32-bit floating-point overflow during depth evaluation. This article examines how ESM mitigates this numerical instability in GLSL by reordering operations to evaluate depth differences inside a single exponential function, leveraging negative exponents to prevent float truncation, and using log-space transformations to avoid precision limits in hardware filtering.

The Overflow Problem in Standard ESM

The standard ESM approximation represents shadow testing using the function:

\[f(d, z) = e^{-c(z - d)} = e^{cd} \cdot e^{-cz}\]

where \(z\) is the receiver depth, \(d\) is the occluder depth stored in the shadow map, and \(c\) is an exponential warp factor (commonly between 40 and 80).

In GLSL, standard single-precision floats (float, FP32) have a maximum representable value of approximately \(3.4028 \times 10^{38}\). The standard library function exp(x) overflows into infinity (+inf) or generates NaN values when \(x > 88.7228\). When \(c\) is set to high values to prevent light bleeding on surfaces far from the occluder, computing \(e^{cd}\) directly in the shadow pass or \(e^{cz}\) in the fragment shader rapidly exceeds this threshold, corrupting the shadow output.

Moving Depth Differences Inside the Exponent

The primary mitigation in GLSL is to avoid evaluating individual exponentials for depth values. Instead of computing \(e^{cd}\) and \(e^{-cz}\) separately, the shader computes the difference between occluder and receiver depths before applying the exponential operation:

\[\text{Shadow} = \text{clamp}(e^{-c(z - d)}, 0.0, 1.0)\]

In GLSL fragment code, this is implemented as:

float occluderDepth = texture(shadowMap, shadowCoord.xy).r;
float receiverDepth = shadowCoord.z;

// Evaluate difference prior to exponentiation
float depthDifference = receiverDepth - occluderDepth;
float shadowFactor = clamp(exp(-c * depthDifference), 0.0, 1.0);

When a surface is in shadow (\(z > d\)), the difference \((z - d)\) is positive, making the exponent \(-c(z - d)\) strictly negative. In IEEE 754 floating-point arithmetic, large negative exponents underflow gracefully to zero rather than overflowing to infinity.

When a surface is fully lit (\(z \le d\)), the exponent is positive or zero. While positive values can cause overflow if \(z \ll d\) (such as depth bias mismatches or receiver geometry far in front of an occluder), applying clamp or bounding the depth difference prevents numerical instability:

float depthDifference = max(receiverDepth - occluderDepth, 0.0);
float shadowFactor = exp(-c * depthDifference);

Mitigating Texture Filter Overflow with Log-ESM

While evaluating differences inside exp() protects the lighting pass, linear hardware filtering across shadow map texels requires storing pre-warped values. If \(e^{cd}\) is stored directly in an FP32 texture, texel values will overflow during texture generation.

To resolve this, implementations use Log-ESM. In Log-ESM, depths are retained in linear space or converted into log-space textures. When generating the shadow map with custom blur passes (such as separable Gaussian blurs), the operations are performed using the log-sum-exp trick to avoid intermediate overflow:

\[\ln\left(\sum w_i e^{c d_i}\right) = c d_{\max} + \ln\left(\sum w_i e^{c (d_i - d_{\max})}\right)\]

In this equation, \(d_{\max}\) is the maximum depth in the blur kernel, and \(w_i\) represents the kernel weights. Because \(d_i - d_{\max} \le 0\), the values passed to exp() are guaranteed not to exceed \(1.0\), completely eliminating floating-point overflow during filtering.

In the final GLSL lighting pass, if the texture stores the filtered occluder depth directly (\(d\)), the final lookup remains:

float shadow = clamp(exp(-c * (shadowCoord.z - texture(shadowMap, shadowCoord.xy).r)), 0.0, 1.0);

This formulation guarantees that the shader executes entirely within the stable range of FP32 hardware arithmetic.