Spherical Harmonics for Irradiance Probes

Diffuse illumination is a convolution. The light arriving at a surface point is the incoming radiance integrated against a cosine lobe, and that integration is a smoothing operation with a very specific frequency response: the cosine lobe has almost no energy above the second spherical harmonic band. Everything above that band contributes nothing to the result. Storing more than a few coefficients per probe is therefore wasted work, and this is the observation that made spherical harmonic irradiance practical for real-time rendering.

The basis and its ordering

Spherical harmonics are the spherical analogue of a Fourier series. A function on the sphere is expanded as a sum over bands l = 0, 1, 2, …, where band l contains 2l+1 basis functions. The first band is a constant, the second is a set of three linear functions of direction (effectively x, y, z), and the third is quadratic.

For rendering, the coefficient count matters. Bands 0 through 2 contain 1 + 3 + 5 = 9 coefficients. Bands 0 through 3 would need 16, and 0 through 4 need 25.

The basis functions are usually indexed by a single integer and computed from a direction. The first nine, in the convention used by most graphics code, are proportional to 1; to y, z, x; to xy, yz, 3z²−1, xz, x²−y². The constants in front of each vary with the normalization convention, which is the single most common source of brightness errors in a probe-based pipeline. A common real-valued normalization gives

$$Y_{00} = \frac{1}{2}\sqrt{\frac{1}{\pi}}, \quad Y_{1,-1} \propto y, \quad Y_{10} \propto z, \quad Y_{11} \propto x$$

with the projection coefficient of a function g being the integral of g times the basis function over the sphere. If the projection and the reconstruction use different constants, the result is scaled by the product of the two, which is usually a factor like π or 4π, yielding a plausible-looking image that is uniformly too bright or too dark.

import math

# Real-valued SH basis, bands 0..2, for a unit direction (x, y, z).
def sh_basis(x, y, z):
    return [
        0.282095,                                  # band 0
        0.488603 * y, 0.488603 * z, 0.488603 * x,  # band 1
        1.092548 * x * y, 1.092548 * y * z,
        0.315392 * (3.0 * z * z - 1.0),
        1.092548 * x * z,
        0.546274 * (x * x - y * y),                # band 2
    ]

# Project a set of radiance samples into coefficients.
def project(samples, directions):
    coeffs = [0.0] * 9
    for L, (x, y, z) in zip(samples, directions):
        b = sh_basis(x, y, z)
        for i in range(9):
            coeffs[i] += L * b[i]
    n = len(samples)
    return [c * (4.0 * math.pi / n) for c in coeffs]   # sphere solid angle / count

The weight 4π/n is the solid angle per sample when the samples are distributed uniformly over the sphere. For samples concentrated in a cube map’s texels, the correct weight is each texel’s solid angle, which varies across the faces; using a uniform weight biases the result toward the poles of the cube’s parameterization.

The cosine convolution in coefficient space

Convolution on the sphere becomes multiplication in the coefficient domain, one factor per band. Convolving the environment with the clamped cosine kernel therefore scales coefficient l by a factor that depends only on l:

$$\hat{A}_l = 2\pi \int_0^{\pi/2} P_l(\cos\theta), \cos\theta, \sin\theta, d\theta$$

Evaluating this for the first few bands gives approximately 3.14159, 2.09440, 0.78540, 0, −0.13090, and so on. The exact values are π, 2π/3, π/4, 0, −π/24. The integration stops at the horizon because the clamped cosine is zero below it, and that half-range integral is exactly what produces the sequence. The coefficient at l = 3 vanishes: a clamped cosine has no third-order component, which is why an irradiance map needs only nine coefficients and why fitting more of them can only add noise.

The negative value at l = 4 is worth noticing. A clamped cosine is not positive-definite on the sphere in the coefficient sense, and truncating the series at l = 2 of a function that has l = 4 content produces ringing. This is the Gibbs phenomenon in spherical coordinates: overshoot near sharp features, which in lighting terms means a dark halo around a bright window or a bright fringe along a wall edge.

That l = 3 coefficient being exactly zero has a further consequence: any information the environment carries at l = 3 is invisible to a diffuse surface, so it can be removed without changing the diffuse result and without any visible artifact.

Ringing and what to do about it

Ringing appears whenever the reconstructed irradiance goes negative. Because irradiance is a nonnegative quantity, a negative value is a signal that the representation has overshot, and clamping to zero is the minimal fix. Clamping changes the energy in the reconstructed function, and it produces flat spots where the true function is a shallow gradient.

The standard remedy is windowing: attenuate the higher-band coefficients before reconstruction so that the truncated series decays without overshoot. A cosine window applied to bands 1 and 2 reduces the sharpness of the reconstruction slightly while eliminating the negative lobe. The window’s strength is a trade-off between accuracy at smooth features and suppression of overshoot at sharp ones. For most indoor environments a mild window is invisible; for a scene dominated by a single bright source it visibly dims the environment, because the window is attenuating the very directional information that conveys the source.

A second approach is to convolve the environment with a wider kernel before truncating, which is physically motivated if the receiver has some roughness or if the surface is not purely Lambertian. But doing so changes the answer for a truly diffuse surface, which is a correctness problem rather than a filtering one.

/* Evaluate irradiance in direction n from nine coefficients, with a
   band-1/2 window w1, w2 applied before reconstruction. */
static float irradiance(const float c[9], const float n[3],
                        float w1, float w2) {
    float b[9];
    b[0] = 0.282095f;
    b[1] = 0.488603f * n[1];  b[2] = 0.488603f * n[2];  b[3] = 0.488603f * n[0];
    b[4] = 1.092548f * n[0] * n[1];
    b[5] = 1.092548f * n[1] * n[2];
    b[6] = 0.315392f * (3.0f * n[2] * n[2] - 1.0f);
    b[7] = 1.092548f * n[0] * n[2];
    b[8] = 0.546274f * (n[0] * n[0] - n[1] * n[1]);

    const float a0 = 3.141593f, a1 = 2.094395f, a2 = 0.785398f;
    float e = c[0] * b[0] * a0;
    for (int i = 1; i < 4; i++) e += c[i] * b[i] * a1 * w1;
    for (int i = 4; i < 9; i++) e += c[i] * b[i] * a2 * w2;
    return e < 0.0f ? 0.0f : e;
}

Probes in space, and what a probe cannot know

A single probe stores the irradiance arriving at one point. Applying it across a surface assumes the incoming radiance is the same at every point, which is true only for distant illumination. For a probe near a wall, the wall blocks half the environment, and the coefficients must have been captured with that occlusion baked in.

The usual compromise places probes on a coarse grid and interpolates coefficients between them, with a visibility test or a parallax correction applied at lookup. Interpolating in coefficient space is equivalent to interpolating the reconstructed functions, because the reconstruction is linear in the coefficients, which is convenient. The error it introduces is a blur across occlusion boundaries: a probe on one side of a corner leaks light from the other side, and the leak appears as a soft bright region in the corner.

The band-limited nature sets a hard floor on spatial detail. Two probes separated by a distance at which the lighting changes faster than band 2 can represent will produce a visible interpolation artifact regardless of how many probes are added. The fix is more angular detail in the probe, not more probes; the sampling and pdf machinery for doing that accurately is the same one used for direct lighting.

What this representation is not for

Spherical harmonics are poor at representing specular reflections. A sharp reflection is close to a delta function in angle, requiring a very high band count; truncating such an environment to nine coefficients turns a mirror into a smeared gradient. The same is true for shadows with hard edges and for small bright sources. Specular image-based lighting solves that problem with a separate representation, a prefiltered map and a BRDF lookup table, because the specular integral needs a different decomposition than the diffuse one.

They are also poor at representing a radiance field that must be re-rendered from many viewpoints on a moving object, because the coefficients are tied to the environment around the probe, and a rotating object sees a time-varying convolution that would have to be recomputed. The rotation of a band-limited representation can be done exactly for the low bands, which is what makes prefiltered environment maps rotatable with a small matrix multiply, but the rotation must respect the same normalization convention as everything else in the pipeline.

The narrow, correct role is diffuse illumination from a moderately smooth environment, where nine coefficients and a per-material albedo produce the right answer cheaply. Using them outside that role is not a failure of the basis; it is a failure to check whether the function being represented lives in the space.