Spherical harmonics for view-dependent colour

A real surface doesn't look the same from every side: a glossy car door catches a highlight, a leaf glows when backlit. A radiance field has to store colour as a function of viewing direction, c(d) for unit vectors d, and it has to do so for millions of primitives with a handful of numbers each. 3D Gaussian Splatting (Kerbl et al. 2023) does it the way real-time graphics has for twenty years: expand c(d) in real spherical harmonics (SH) and keep the first 16 coefficients per colour channel.

This lesson builds the exact evaluator the 3DGS reference code uses — constants, signs, offset and clamp — then goes the other way: projecting a function on the sphere onto SH, what truncation does to it, and the classic use that made SH famous in games, irradiance environment maps.

1. A Fourier basis for the sphere

On a circle, any periodic function is a sum of cos(mθ) and sin(mθ); keep the low frequencies and you get a smooth approximation. Spherical harmonics are the same idea on the sphere. They are grouped in bands l = 0, 1, 2, …; band l holds 2l + 1 functions Y_lm, m = −l … l, each a homogeneous polynomial of degree l in the components of the unit vector (x, y, z):

• band 0: a constant (1 function),
• band 1: x, y, z up to scale (3),
• band 2: xy, yz, 3z² − 1, xz, x² − y² (5),
• band 3: seven cubics (7).

Bands 0 … d together hold 1 + 3 + 5 + … + (2d + 1) = (d + 1)² functions: 4, 9, 16 for degrees 1, 2, 3. They are orthonormal over the sphere, ∫ Y_j(ω) Y_k(ω) dω = δ_jk, which is what makes projection a dot product (section 5).

The constants are normalisations: Y_00 = 1/(2√π) = 0.28209… because the sphere has area 4π; band 1 is √(3/4π)·(y, z, x) = 0.48860…·(y, z, x); band 2's ½√(15/π) = 1.09254…, ¼√(5/π) = 0.31539…, ¼√(15/π) = 0.54627…; band 3's ¼√(35/2π), ½√(105/π), ¼√(21/2π), ¼√(7/π), ¼√(105/π). You never derive them at run time — you copy them.

2. The 3DGS convention, exactly

The CUDA rasterizer (computeColorFromSH in cuda_rasterizer/forward.cu) and the Python helper (eval_sh in utils/sh_utils.py) agree on this order and these signs:

result = SH_C0 * sh[0];
result = result - SH_C1 * y * sh[1] + SH_C1 * z * sh[2] - SH_C1 * x * sh[3];
result = result + SH_C2[0] * xy * sh[4] + SH_C2[1] * yz * sh[5]
       + SH_C2[2] * (2.0f * zz - xx - yy) * sh[6]
       + SH_C2[3] * xz * sh[7] + SH_C2[4] * (xx - yy) * sh[8];
result = result + SH_C3[0] * y * (3.0f * xx - yy) * sh[9] + SH_C3[1] * xy * z * sh[10]
       + SH_C3[2] * y * (4.0f * zz - xx - yy) * sh[11]
       + SH_C3[3] * z * (2.0f * zz - 3.0f * xx - 3.0f * yy) * sh[12]
       + SH_C3[4] * x * (4.0f * zz - xx - yy) * sh[13]
       + SH_C3[5] * z * (xx - yy) * sh[14] + SH_C3[6] * x * (xx - 3.0f * yy) * sh[15];

with

SH_C0 = 0.28209479177387814f;
SH_C1 = 0.4886025119029199f;
SH_C2[] = { 1.0925484305920792f, -1.0925484305920792f, 0.31539156525252005f,
           -1.0925484305920792f, 0.5462742152960396f };
SH_C3[] = { -0.5900435899266435f, 2.890611442640554f, -0.4570457994644658f,
             0.3731763325901154f, -0.4570457994644658f, 1.445305721320277f,
            -0.5900435899266435f };

Two details trip everyone up:

• Order within band 1 is y, z, x. It's m = −1, 0, 1, and the m < 0 functions are the "sine-like" ones (y is sin φ on the equator).
• The signs are the Condon–Shortley phase (−1)^m: the textbook real SH (Wikipedia's table) times −1 for odd m. That is the minus on y and x in band 1, and the negative entries SH_C2[1], SH_C2[3], SH_C3[0, 2, 4, 6]. Flipping signs keeps the basis orthonormal, so any sign convention works for training — but a checkpoint trained with one renders wrongly with another, and the result still looks plausible. Test against reference values.

#include <print>

constexpr float C0 = 0.28209479177387814f;
constexpr float C1 = 0.4886025119029199f;

// Bands 0 and 1 in the 3DGS order and signs.
void basis1(float x, float y, float z, float out[4]) {
    out[0] = C0;
    out[1] = -C1 * y;
    out[2] = C1 * z;
    out[3] = -C1 * x;
}

int main() {
    const float axes[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
    const char* names[3] = {"+x", "+y", "+z"};
    for (int a = 0; a < 3; ++a) {
        float y[4];
        basis1(axes[a][0], axes[a][1], axes[a][2], y);
        std::println("{}: {:+.4f} {:+.4f} {:+.4f} {:+.4f}", names[a], y[0], y[1], y[2], y[3]);
    }
}

On +x only Y_0 and Y_3 = −C1 survive; on +z only the m = 0 functions are nonzero — a quick test that catches most order and sign bugs.

3. From coefficients to a pixel colour

A Gaussian stores 16 RGB coefficients. To shade it for a camera:

• Direction. dir = normalize(pos − campos): from the camera towards the Gaussian (the reference uses the Gaussian's mean, not the pixel's ray — one colour per splat per frame). Forget to normalise and band l scales with distance^l: colour changes as you walk.
• The +0.5. Colours are stored shifted: the point cloud is initialised with sh[0] = RGB2SH(rgb) = (rgb − 0.5) / C0, so all-zero coefficients mean mid-grey and the DC term is centred on zero. Evaluation adds the 0.5 back.
• The clamp is max(colour, 0) after the offset; no upper clamp. The kernel remembers which channels were clamped so the backward pass gives them zero gradient.

#include <algorithm>
#include <print>

constexpr float C0 = 0.28209479177387814f;

float rgb_to_sh0(float rgb) { return (rgb - 0.5f) / C0; }
float sh0_to_rgb(float sh0) { return std::max(sh0 * C0 + 0.5f, 0.0f); }

int main() {
    for (float rgb : {0.0f, 0.25f, 0.5f, 1.0f}) {
        float c = rgb_to_sh0(rgb);
        std::println("rgb {:.2f} -> sh0 {:+.4f} -> rgb {:.2f}", rgb, c, sh0_to_rgb(c));
    }
}

Training: one band at a time

3DGS optimises the DC term from the start and raises the active SH degree by one every 1000 iterations (oneupSHdegree) until it reaches 3. Early on the geometry is wrong, and high bands would happily explain geometric errors as "view-dependent colour"; growing the degree slowly makes the model explain what it can with geometry first. The coefficients are stored for degree 3 from the start; evaluation just uses the first (d + 1)². The learning rate for the higher bands is 20× smaller than for the DC term.

Memory: most of a splat is SH

A 3-million-splat scene is ~700 MB of float32, ~575 MB of it SH. Compression work (quantising SH, codebooks, dropping to degree 0 for diffuse splats, distilling SH into a tiny MLP) starts here; so does streaming — send DC first, refine bands later.

4. Projection: from a function to coefficients

Orthonormality turns "best approximation" into "dot product":

c_k = ∫_{S²} f(ω) Y_k(ω) dω

and Σ_k c_k Y_k with k < (d+1)² is the least-squares best band-limited approximation of f. With N directions spread uniformly over the sphere, each stands for 4π/N steradians:

c_k ≈ (4π / N) Σ_i f(ω_i) Y_k(ω_i)

Random uniform directions converge like 1/√N. A Fibonacci sphere (points on a golden-angle spiral, evenly spaced in z) is deterministic and far more accurate for smooth integrands — ideal for tests:

#include <cmath>
#include <numbers>
#include <print>

constexpr float C0 = 0.28209479177387814f;
constexpr float C1 = 0.4886025119029199f;

int main() {
    const int n = 2000;
    const double golden = std::numbers::pi * (3.0 - std::sqrt(5.0));
    double gram[4][4] = {};
    for (int i = 0; i < n; ++i) {
        double z = 1.0 - (2.0 * i + 1.0) / n, r = std::sqrt(1.0 - z * z), phi = golden * i;
        double x = r * std::cos(phi), y = r * std::sin(phi);
        double b[4] = {C0, -C1 * y, C1 * z, -C1 * x};
        for (int j = 0; j < 4; ++j)
            for (int k = 0; k < 4; ++k) gram[j][k] += 4.0 * std::numbers::pi / n * b[j] * b[k];
    }
    for (auto& row : gram) std::println("{:+.4f} {:+.4f} {:+.4f} {:+.4f}", row[0], row[1], row[2], row[3]);
}

The Gram matrix comes out as the identity: the basis is orthonormal, signs and all.

Band-limiting rings

A function that is itself a sum of bands ≤ d comes back exactly. Anything sharper is truncated, and truncated Fourier series ring (the Gibbs phenomenon). The clamped cosine max(cos θ, 0) — the shape of diffuse lighting — has coefficients C0·π, C1·2π/3 and (π/4)·… on the m = 0 functions of bands 0–2; band 3's is exactly zero (all odd bands above 1 vanish). Its degree-3 reconstruction overshoots to 1.0625 at the pole and dips below zero behind it:

That is the practical ceiling of degree-3 SH: soft view dependence (sheen, broad highlights, Fresnel-ish rims) works; a mirror reflection or a tight specular highlight cannot be represented and gets "baked" as blur — or, in 3DGS, faked with extra Gaussians that appear only from some directions.

Rotations stay inside a band

Rotating f mixes the coefficients within each band (by a (2l+1) × (2l+1) orthogonal matrix, the Wigner D-matrix) but never across bands. So the energy per band Σ_m c_lm² is rotation-invariant — a cheap sanity test for a projector, and the reason an object's SH lighting can be rotated without re-projecting. It also means 3DGS splats must rotate their SH when you rotate a scene (a detail many editing tools get wrong: the colours of a rotated scene look "off" until the SH is rotated too).

5. The older use: irradiance environment maps

Before radiance fields, SH were the trick for diffuse lighting from an environment map (Ramamoorthi & Hanrahan 2001, An Efficient Representation for Irradiance Environment Maps). Irradiance at normal n is the radiance L convolved with the clamped cosine. Convolution with a rotationally symmetric kernel is a per-band scale in SH:

E(n) = Σ_{l,m} Â_l · L_lm · Y_lm(n), Â_0 = π, Â_1 = 2π/3, Â_2 = π/4, Â_3 = 0, …

Â_l falls off fast (and is 0 for odd l > 1), so 9 coefficients (bands 0–2) reproduce irradiance with an average error around 1 %. Project the environment once (27 floats for RGB), then shading a pixel is a 9-term polynomial in the normal — still how many engines do ambient diffuse light probes. The formula doesn't care about sign conventions as long as you project and evaluate with the same basis.

#include <algorithm>
#include <cmath>
#include <numbers>
#include <print>

// A "sky" that is 1 above the horizon and 0 below; irradiance on an upward
// and a sideways surface, from its projection onto bands 0-2 (m = 0 terms only:
// the sky is symmetric about z).
int main() {
    const double pi = std::numbers::pi;
    const int n = 4000;
    double c0 = 0, c2 = 0, c6 = 0;
    for (int i = 0; i < n; ++i) {
        double z = 1.0 - (2.0 * i + 1.0) / n;   // uniform in z = uniform on the sphere
        double L = z > 0 ? 1.0 : 0.0;
        double w = 4 * pi / n;
        c0 += w * L * 0.28209479177387814;
        c2 += w * L * 0.4886025119029199 * z;
        c6 += w * L * 0.31539156525252005 * (3 * z * z - 1);
    }
    auto E = [&](double nz) {
        return pi * c0 * 0.28209479177387814 + 2 * pi / 3 * c2 * 0.4886025119029199 * nz +
               pi / 4 * c6 * 0.31539156525252005 * (3 * nz * nz - 1);
    };
    std::println("up {:.3f} (exact {:.3f}), sideways {:.3f} (exact {:.3f})", E(1.0), pi, E(0.0), pi / 2);
}

6. Implementation notes

• Keep one sh_basis(degree, dir, out) with early returns per band; both evaluation and projection call it, so they can never disagree.
• Use (d + 1)² everywhere — 2d + 1 is the size of the last band only.
• float is plenty: the basis values are bounded (|Y| < 1.1 for l ≤ 3) and the sums are short. Normalise with a reciprocal square root in hot loops.
• A zero-length pos − campos (camera exactly at the Gaussian) gives NaN; real renderers never hit it because such splats are culled by the near plane.

In Python (NumPy / PyTorch). sh_utils.eval_sh(deg, sh, dirs) takes sh of shape [..., 3, (deg+1)**2] and unit dirs [..., 3] and is written exactly like the C++ above, vectorised over all Gaussians; the model calls it as colors = torch.clamp_min(eval_sh(...) + 0.5, 0.0) and initialises features_dc = RGB2SH(colors). In NumPy, build the [N, 16] basis matrix once and the colours are one matmul with the [16, 3] coefficients.

In an interview. Expect: "how many floats per splat does SH take and why?" (16 × 3 = 48, degree 3), "why not just store an RGB?" (view dependence — specular, sheen), "what can't degree 3 represent?" (sharp reflections; ringing), "how do you rotate a splat scene?" (rotate positions and covariances and the SH coefficients, band by band), and "how would you compress them?" (quantise / codebook the 45 higher-band floats, prune to lower degree where view dependence is weak).