Spherical Harmonic Irradiance Evaluation in GLSL

Spherical harmonic irradiance evaluation is an efficient technique used in real-time graphics to compute diffuse indirect lighting from distant environment radiance. By projecting environment lighting into low-order spherical harmonic (SH) basis functions—typically order 3, totaling 9 coefficients per color channel—the expensive integral of cosine-weighted incident radiance over a hemisphere is reduced to a simple polynomial evaluation. In GLSL, this evaluation takes a normalized surface normal vector and the 9 precomputed SH coefficients to reconstruct diffuse irradiance in a handful of multiply-add operations per pixel.

Theoretical Foundations

Diffuse Lambertian reflection acts as a low-pass filter on incident radiance. According to the Funk-Hecke theorem, convolving incident radiance with a clamped cosine lobe (\(\max(\mathbf{n} \cdot \mathbf{\omega}, 0)\)) dampens higher-frequency signals. Ravi Ramamoorthi and Pat Hanrahan demonstrated that using just three bands of spherical harmonics (degrees \(l = 0, 1, 2\), yielding 9 scalar basis functions) captures more than 99% of diffuse irradiance energy.

The reconstructed irradiance \(E(\mathbf{n})\) along a surface normal \(\mathbf{n} = (x, y, z)\) is expressed as:

\[E(\mathbf{n}) = \sum_{l=0}^{2} \sum_{m=-l}^{l} \hat{A}_l L_{lm} Y_{lm}(\mathbf{n})\]

Where:

  • \(L_{lm}\) are the spherical harmonic coefficients of the incident lighting environment (stored as RGB triplets).
  • \(Y_{lm}(\mathbf{n})\) are the real spherical harmonic basis functions evaluated at \(\mathbf{n}\).
  • \(\hat{A}_l\) are the cosine lobe convolution factors:
    • \(\hat{A}_0 = \pi\)
    • \(\hat{A}_1 = \frac{2\pi}{3}\)
    • \(\hat{A}_2 = \frac{\pi}{4}\)

From Spherical Harmonics to Polynomials

Evaluating raw trigonometric basis functions on the GPU is computationally inefficient. Because the normal \(\mathbf{n} = (x, y, z)\) is a unit Cartesian vector, the products \(\hat{A}_l Y_{lm}(\mathbf{n})\) can be expanded into Cartesian polynomials:

  • Band 0 (\(l=0\)): Constant term.
  • Band 1 (\(l=1\)): Linear terms proportional to \(y\), \(z\), and \(x\).
  • Band 2 (\(l=2\)): Quadratic terms proportional to \(xy\), \(yz\), \((3z^2 - 1)\), \(xz\), and \((x^2 - y^2)\).

By folding the convolution factors \(\hat{A}_l\) and the normalization constants of \(Y_{lm}\) directly into the coefficients, the evaluation collapses into a quadratic polynomial:

\[E(\mathbf{n}) = c_1 L_{00} - c_3 L_{20} + 2 c_2 (L_{1-1}y + L_{10}z + L_{11}x) + 2 c_4 (L_{2-1}yz + L_{21}xz + L_{2-2}xy) + c_3 L_{20} (3z^2) + c_5 L_{22} (x^2 - y^2)\]

The scalar constants derived from the basis definitions are:

  • \(c_1 = 0.282095 \times \pi \approx 0.886227\)
  • \(c_2 = 0.488603 \times \frac{2\pi}{3} \approx 1.023326\)
  • \(c_3 = 0.315392 \times \frac{\pi}{4} \approx 0.247708\)
  • \(c_4 = 1.092548 \times \frac{\pi}{4} \approx 0.858086\)
  • \(c_5 = 0.546274 \times \frac{\pi}{4} \approx 0.429043\)

GLSL Implementation

Coefficients are typically passed to shaders via a uniform buffer as an array of nine vec3 values representing the RGB channels for each SH index (\(L_0\) to \(L_8\)).

// Incident lighting SH coefficients: L00, L1-1, L10, L11, L2-2, L2-1, L20, L21, L22
uniform vec3 u_SH[9];

vec3 evaluateSHIrradiance(vec3 normal) {
    // Ensure normal is normalized
    vec3 n = normalize(normal);

    // Pre-derived constant factors (A_l * Y_lm constants)
    const float c1 = 0.429043;
    const float c2 = 0.511664;
    const float c3 = 0.743125;
    const float c4 = 0.886227;
    const float c5 = 0.247708;

    // Linear components
    vec3 linear = 2.0 * c2 * (u_SH[1] * n.y + u_SH[2] * n.z + u_SH[3] * n.x);

    // Quadratic polynomial components
    vec3 quadratic = c1 * u_SH[4] * (n.x * n.y) +
                    c1 * u_SH[5] * (n.y * n.z) +
                    c3 * u_SH[6] * (n.z * n.z) - c5 * u_SH[6] +
                    c1 * u_SH[7] * (n.x * n.z) +
                    c1 * u_SH[8] * (n.x * n.x - n.y * n.y);

    // Constant scalar component (Band 0 combined with Band 2 offset)
    vec3 constant = c4 * u_SH[0];

    // Total diffuse irradiance
    return max(constant + linear + quadratic, vec3(0.0));
}

Matrix Formulation Alternative

An alternative implementation rewrites the evaluation as a quadratic surface form:

\[E(\mathbf{n}) = \mathbf{v}^T M \mathbf{v}\]

Where \(\mathbf{v} = (x, y, z, 1)^T\) and \(M\) is a symmetric \(4 \times 4\) matrix constructed for each color channel. In GLSL, this is expressed using dot products:

uniform mat4 u_SH_R;
uniform mat4 u_SH_G;
uniform mat4 u_SH_B;

vec3 evaluateSHMatrix(vec3 normal) {
    vec4 v = vec4(normalize(normal), 1.0);
    return max(vec3(
        dot(v, u_SH_R * v),
        dot(v, u_SH_G * v),
        dot(v, u_SH_B * v)
    ), 0.0);
}

While the matrix form simplifies the shader code to three matrix-vector multiplications, pre-baking the constants into the 9 polynomial terms directly requires fewer overall arithmetic operations, making the polynomial approach the preferred choice for execution speed. In both variations, clamping with max(..., 0.0) prevents negative irradiance artifacts resulting from SH ringing.