Hamburger Bounds for Moment Shadow Maps in GLSL

This article explains how to evaluate four statistical depth moments using Hamburger bounds to compute shadow factors in Moment Shadow Maps (MSM). Moment Shadow Mapping utilizes a four-component vector of trigonometric or power depth moments to reconstruct a sharp upper bound on occlusion while eliminating common artifacts like light leaking. By leveraging Hamburger's theorem and solving a Hankel linear system via Cholesky decomposition, developers can evaluate shadow coverage in GLSL with closed-form, per-fragment matrix operations.

Theoretical Foundation

In Moment Shadow Mapping, the shadow map records the first four power moments of depth \(z \in [0, 1]\):

\[b = (b_1, b_2, b_3, b_4)^T = \left( \mathbb{E}[z], \mathbb{E}[z^2], \mathbb{E}[z^3], \mathbb{E}[z^4] \right)^T\]

Given a surface fragment at depth \(z_f\), the goal is to evaluate the maximum probability that an occluder exists in front of the fragment (\(z < z_f\)). Hamburger's theorem provides a sharp upper bound for this cumulative distribution function using a \(3 \times 3\) positive-semidefinite Hankel matrix of moments:

\[M = \begin{bmatrix} 1 & b_1 & b_2 \\ b_1 & b_2 & b_3 \\ b_2 & b_3 & b_4 \end{bmatrix}\]

Applying Hamburger’s bound requires calculating the extremal distributions supported by the moment vector. In practice, this is solved by decomposing \(M\) and finding the roots of a quadratic polynomial that yields the bound on occlusion.

Mathematical Formulation

1. Moment Regularization

Because filtering and precision limits can cause \(M\) to become negative or singular, the moment vector \(b\) must first be regularized using a small bias parameter \(\alpha\) (typically \(10^{-4}\) to \(10^{-3}\)):

\[b' = (1 - \alpha) b + \alpha \begin{bmatrix} 0.5 \\ 0.5 \\ 0.5 \\ 0.5 \end{bmatrix}\]

2. Cholesky Factorization

To invert the system efficiently, compute the lower-triangular Cholesky factorization \(M = L L^T\):

\(L_{11} = 1\)
\(L_{21} = b'_1\)
\(L_{31} = b'_2\)
\(L_{22} = \sqrt{b'_2 - L_{21}^2}\)
\(L_{32} = \frac{b'_3 - L_{21} L_{31}}{L_{22}}\)
\(L_{33} = \sqrt{b'_4 - L_{31}^2 - L_{32}^2}\)

3. Forward Substitution

Solve the linear equation \(L y = v\), where \(v = (1, z_f, z_f^2)^T\):

\(y_1 = 1\)
\(y_2 = \frac{z_f - L_{21}}{L_{22}}\)
\(y_3 = \frac{z_f^2 - L_{31} - L_{32} y_2}{L_{33}}\)

4. Root Finding and Bound Evaluation

The vector \(y\) defines the coefficients of a polynomial whose roots \(z_1\) and \(z_2\) define the points of support for the extremal distribution. By computing the squared norm and determining where \(z_f\) falls relative to the roots, the bound reduces to an evaluation of occlusion weight:

\[G(z_f) = \frac{y_1^2 + y_2^2}{y_1^2 + y_2^2 + y_3^2}\]

When evaluating the condition with the roots:

If \(z_f \le z_1\), the surface is fully lit (\(G = 0\)).
If \(z_1 < z_f \le z_2\), the surface is partially shadowed according to the lower extremal point.
If \(z_f > z_2\), the surface is evaluated against the upper bound.

GLSL Implementation

The following complete GLSL function implements the regularization, Cholesky factorization, forward substitution, and Hamburger bound evaluation for four moments:

// Evaluates shadow intensity in [0, 1], where 0 is fully lit and 1 is fully shadowed
float evaluateHamburgerBound(vec4 moments, float fragmentDepth, float bias) {
    // Step 1: Regularization to avoid non-positive definite Hankel matrices
    vec4 b = mix(moments, vec4(0.5, 0.5, 0.5, 0.5), bias);

    // Fragment depth vector: v = (1.0, z, z^2)
    float z = fragmentDepth;

    // Step 2: Compute Cholesky factorization L of Hankel matrix
    // L11 is 1.0, L21 = b.x, L31 = b.y
    float L21 = b.x;
    float L31 = b.y;

    float L22 = sqrt(max(0.0, b.y - L21 * L21));
    float invL22 = 1.0 / max(L22, 1e-6);
    float L32 = (b.z - L21 * L31) * invL22;
    
    float L33 = sqrt(max(0.0, b.w - L31 * L31 - L32 * L32));
    float invL33 = 1.0 / max(L33, 1e-6);

    // Step 3: Forward substitution: Solve L * y = (1.0, z, z^2)
    float y1 = 1.0;
    float y2 = (z - L21) * invL22;
    float y3 = (z * z - L31 - L32 * y2) * invL33;

    // Step 4: Evaluate the polynomial roots and Hamburger bound
    // Vector c = (1, -z1 - z2, z1 * z2) derived from the null-space computation
    float ySquareSum = y1 * y1 + y2 * y2;
    float totalNorm = ySquareSum + y3 * y3;

    // Determine quadratic polynomial coefficients: c0 + c1*z + c2*z^2
    float c1 = -y2 * invL22 + (L32 * invL22) * (y3 * invL33);
    float c2 = -y3 * invL33;
    float c0 = 1.0 - L21 * c1 - L31 * c2;

    // Solve quadratic equation: c2 * z^2 + c1 * z + c0 = 0
    float desc = c1 * c1 - 4.0 * c2 * c0;
    float sqrtDesc = sqrt(max(0.0, desc));
    float z1 = (-c1 - sqrtDesc) / (2.0 * c2);
    float z2 = (-c1 + sqrtDesc) / (2.0 * c2);

    // Order the roots
    if (z1 > z2) {
        float temp = z1;
        z1 = z2;
        z2 = temp;
    }

    // Step 5: Classify occlusion based on root intervals
    if (z <= z1) {
        return 0.0; // Fully in light
    } else if (z <= z2) {
        // Intermediate bound
        return (ySquareSum) / totalNorm;
    } else {
        // Fully behind both occluder approximations
        return 1.0;
    }
}

Practical Considerations

Numerical Stability: When \(L_{22}\) or \(L_{33}\) approach zero, the moment distribution represents a single plane (Dirac delta distribution). Enforcing minimum thresholds via max(L, 1e-6) prevents division-by-zero artifacts.
Texture Formats: Four moments require a 4-channel texture. Using RGBA16F can work, but RGBA32F is recommended to avoid catastrophic cancellation during the \(L_{33}\) root calculations.
Filtering: Because moments are linear combinations of depth expectations, textures can be blurred using standard hardware linear filtering, Gaussian blurs, or compute-based convolutions prior to evaluating the Hamburger bound.