一、背景与动机：为什么需要球谐函数？

环境光照的困境

实时渲染有一道绕不过去的坎：如何高效地表示整个球面上的环境光照，并在每个着色点快速计算它对漫反射的贡献？

最直接的办法是蒙特卡洛积分——对每个着色点，在半球上采样数百个方向，取环境光加权平均：

$$E(n) = \int_{\Omega} L(\omega) \max(0, n \cdot \omega) , d\omega \approx \frac{1}{N} \sum_{i=1}^{N} L(\omega_i) \max(0, n \cdot \omega_i)$$

问题在于：1080P 下有两百万像素，每像素哪怕只采样 64 次，每帧就要做 1.28 亿次环境贴图采样。这在实时场景根本行不通。

工业界的标准答案：球谐函数

UE4/UE5 的 Indirect Lighting Cache、Sky Light 全都基于 L2 球谐函数（9个系数）。Unity 的 Light Probes 同样如此。寒霜引擎、Frostbite 也把 SH 作为环境漫反射光照的核心表示。

原因很简单：Lambertian 漫反射是一个极低频信号——它对环境光的响应就像一个低通滤波器，高频细节完全模糊掉了。用 9 个系数就能捕捉 99% 的漫反射信息，运行时只需 9 次乘加，一个 dot product 搞定。

没有 SH，就没有高效的实时环境漫反射光照。这篇文章从第一性原理出发，完整实现一个 SH 环境光照渲染器，并与蒙特卡洛参考定量对比。

二、核心原理：球谐函数的数学基础

2.1 什么是球谐函数？类比傅里叶级数理解

傅里叶级数把周期函数展开为 $\sin$/$\cos$ 的线性组合——这些基函数在 $[0, 2\pi]$ 上正交，可以表示任意周期函数，频率越高的分量对应越精细的细节。

球谐函数做的事情完全类似，只不过从一维周期域换成了二维球面。球谐函数 $Y_l^m(\theta, \phi)$ 是定义在单位球面 $S^2$ 上的一组正交基函数，它们满足：

$$\int_{S^2} Y_l^m(\omega) Y_{l’}^{m’}(\omega) , d\omega = \delta_{ll’} \delta_{mm’}$$

其中：

  • $l$：阶数（band），$l = 0, 1, 2, \ldots$——类比傅里叶的”频率”
  • $m$：次数（order），$m \in [-l, l]$，每阶有 $2l+1$ 个基函数
  • $l=0$：1 个函数（常数，最低频）
  • $l=1$：3 个函数（偶极子方向性）
  • $l=2$：5 个函数（四极子，捕捉更复杂空间变化）

前 3 阶（L0+L1+L2）共 $1+3+5=9$ 个基函数，这就是实时渲染中”9系数 SH”的由来。

2.2 实数形式的 9 个基函数

球谐函数的复数形式涉及 Legendre 多项式，这里直接给出用于实时渲染的实数归一化形式，参数是单位方向向量 $(x,y,z) = (\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta)$：

$$\begin{aligned}
Y_0^0 &= \frac{1}{2}\sqrt{\frac{1}{\pi}} \quad \text{（常数，与方向无关）} \[6pt]
Y_1^{-1} &= \sqrt{\frac{3}{4\pi}} \cdot y, \quad
Y_1^0 = \sqrt{\frac{3}{4\pi}} \cdot z, \quad
Y_1^1 = \sqrt{\frac{3}{4\pi}} \cdot x \quad \text{（线性方向性）} \[6pt]
Y_2^{-2} &= \frac{1}{2}\sqrt{\frac{15}{\pi}} \cdot xy, \quad
Y_2^{-1} = \frac{1}{2}\sqrt{\frac{15}{\pi}} \cdot yz, \quad
Y_2^0 = \frac{1}{4}\sqrt{\frac{5}{\pi}} (3z^2-1) \
Y_2^1 &= \frac{1}{2}\sqrt{\frac{15}{\pi}} \cdot xz, \quad
Y_2^2 = \frac{1}{4}\sqrt{\frac{15}{\pi}} (x^2-y^2) \quad \text{（二次方向性）}
\end{aligned}$$

直觉解释每个阶的含义：

  • L0（$Y_0^0$）：球面上的平均值——环境光的”整体亮度”
  • L1（$Y_1^{-1}, Y_1^0, Y_1^1$）：三个坐标轴方向上的亮度梯度——“天顶比地面亮”、”太阳在哪个方向”这类低频方向性信息
  • L2（5个函数）：二阶方向性细节——比如”日落时地平线出现的橙色辉光带”这种稍微复杂一点的分布

归一化系数的作用：确保正交性（任意两个不同 SH 基函数的球面积分为 0）、归一性（自身积分为 1）。这使得 SH 投影和重建的系数在数值上稳定、可比较。

代码实现（9 个基函数）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
void shBasis(const Vec3& d, double b[9]) {
    // b[i] 即方向 d 处第 i 个 SH 基函数的值
    // 系数来自实数归一化球谐函数公式

    // L0: 常数项，与方向完全无关
    b[0] = 0.282095;  // 0.5 * sqrt(1/π)

    // L1: 线性项，分别正比于 y, z, x
    // 为什么是 y,z,x 而不是 x,y,z？源自 SH 标准定义的 m=-1,0,+1 顺序
    b[1] = 0.488603 * d.y;  // sqrt(3/(4π)) * y
    b[2] = 0.488603 * d.z;  // sqrt(3/(4π)) * z
    b[3] = 0.488603 * d.x;  // sqrt(3/(4π)) * x

    // L2: 二次项，捕捉方向间的相互关系
    b[4] = 1.092548 * d.x * d.y;            // 0.5*sqrt(15/π) * xy
    b[5] = 1.092548 * d.y * d.z;            // 0.5*sqrt(15/π) * yz
    b[6] = 0.315392 * (3.0*d.z*d.z - 1.0); // 0.25*sqrt(5/π) * (3z²-1)
    b[7] = 1.092548 * d.x * d.z;            // 0.5*sqrt(15/π) * xz
    b[8] = 0.546274 * (d.x*d.x - d.y*d.y); // 0.25*sqrt(15/π) * (x²-y²)
}

2.3 SH 投影：把环境光”压缩”成 9 个向量

任意球面函数 $L(\omega)$ 都可以展开为 SH 级数：

$$L(\omega) \approx \sum_{l=0}^{L} \sum_{m=-l}^{l} c_{lm} , Y_l^m(\omega)$$

SH 系数 $c_{lm}$ 通过球面积分求得（类比傅里叶系数）：

$$c_{lm} = \int_{S^2} L(\omega) , Y_l^m(\omega) , d\omega$$

直觉：$c_{lm}$ 就是”环境光在 $Y_l^m$ 这个方向模式上有多强”。$c_0^0$ 大意味着环境光整体很亮；$c_1^0$（z分量）大意味着”上方比下方亮”（天空比地面亮的程度）。

在实现中，用蒙特卡洛积分来估算这个球面积分。关键是采样的概率分布和权重：

均匀球面采样：$u,v \sim \text{Uniform}[0,1]$，令：
$$\theta = \arccos(1 - 2u), \quad \phi = 2\pi v$$

这个变换通过逆 CDF 方法得到，保证方向均匀分布在球面上（注意 $\theta = \arccos(1-2u)$ 而不是 $\arccos(u)$，后者会导致极点附近过度密集）。

均匀球面采样的 PDF 是 $p(\omega) = \frac{1}{4\pi}$（球面面积的倒数），因此蒙特卡洛估计量的权重是 $\frac{1}{p(\omega)} = 4\pi$，N 个样本的平均权重是 $\frac{4\pi}{N}$：

$$c_{lm} \approx \frac{4\pi}{N} \sum_{i=1}^{N} L(\omega_i) , Y_l^m(\omega_i)$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
// SH 系数结构：9 个向量（每个方向 RGB 三通道）
struct SHCoeffs {
    Vec3 c[9];
    SHCoeffs() { for(int i=0; i<9; i++) c[i] = Vec3(0,0,0); }
};

SHCoeffs projectEnvToSH(std::function<Vec3(const Vec3&)> env, int numSamples) {
    SHCoeffs coeffs;
    std::mt19937 rng(42);  // 固定种子保证可复现
    std::uniform_real_distribution<double> distU(0.0, 1.0);

    // 均匀球面采样的权重：每个样本代表 4π/N 立体角
    double weight = 4.0 * PI / numSamples;

    for(int i = 0; i < numSamples; i++) {
        double u = distU(rng), v = distU(rng);

        // 逆 CDF：将均匀随机数映射为均匀球面方向
        // theta = acos(1-2u) 而不是 acos(u)！
        // 后者在 theta≈0 和 theta≈π（极点）附近密度过高
        double theta = std::acos(1.0 - 2.0 * u);
        double phi   = 2.0 * PI * v;

        Vec3 dir = {
            std::sin(theta) * std::cos(phi),
            std::sin(theta) * std::sin(phi),
            std::cos(theta)
        };

        Vec3 radiance = env(dir);  // 查询该方向的环境辐亮度

        double b[9];
        shBasis(dir, b);  // 计算该方向的 9 个 SH 基函数值

        // 将辐亮度投影到每个基函数上：L(ω) * Y_l^m(ω) * 权重
        for(int j = 0; j < 9; j++) {
            coeffs.c[j] += radiance * b[j];
        }
    }

    // 乘以 MC 权重（4π/N）得到球面积分的估计
    for(int j = 0; j < 9; j++) {
        coeffs.c[j] = coeffs.c[j] * weight;
    }
    return coeffs;
}

这里用 65536 个样本。为什么这个数字够用？Lambertian 漫反射的 SH 系数对高频噪声不敏感，方差主要来自低阶项，65536 个样本的误差已经在 0.1% 以内。

2.4 Lambertian 辐照度的 SH 重建：为什么只要 9 个系数够？

有了 SH 系数，要计算法线 $n$ 方向的辐照度（Irradiance）：

$$E(n) = \int_{S^2} L(\omega) , \underbrace{\max(0, \omega \cdot n)}_{\text{Lambertian 核}} , d\omega$$

这是环境辐亮度 $L(\omega)$ 与 Lambertian 响应核 $H(n, \omega) = \max(0, \omega \cdot n)$ 的球面卷积。

关键定理（Ramamoorthi & Hanrahan，2001）：球面卷积在 SH 域退化为逐系数乘法。Lambertian 核在 SH 分解后，各阶的系数（称为 Lambertian 特征值 $A_l$）为：

$$A_0 = \pi, \quad A_1 = \frac{2\pi}{3}, \quad A_2 = \frac{\pi}{4}, \quad A_3 = 0, \quad A_4 = -\frac{\pi}{24}, \ldots$$

直觉解释：

  • $A_0 = \pi$：常数分量（整体亮度）以 $\pi$ 倍传递给辐照度——这其实就是 Lambertian 余弦积分 $\int_0^{\pi/2} \cos\theta \sin\theta , d\theta \times 2\pi = \pi$ 的结果
  • $A_1 = 2\pi/3$：方向性分量（”哪边亮”）以 $2/3$ 的衰减传递
  • $A_2 = \pi/4$：二次方向性以 $1/4$ 的衰减传递，已经衰减到 $A_0$ 的 1/4
  • $A_3 = 0$：奇数阶的 $A_3$ 恰好为 0！这是由于 $\cos\theta$ 的奇偶性导致的积分消去
  • $l \geq 3$：系数急剧衰减，几乎为零

这就是核心结论：Lambertian 核是一个强低通滤波器，高阶（$l \geq 3$）的环境光信息对漫反射几乎没有贡献。用 9 个系数近似误差小于 1%。

因此辐照度的 SH 重建公式变成：

$$E(n) = \sum_{l,m} A_l \cdot c_{lm} \cdot Y_l^m(n)$$

实现只是 9 次乘加，运行时几乎零代价：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
Vec3 shIrradiance(const SHCoeffs& sh, const Vec3& n) {
    // Lambertian 核的 SH 特征值
    // 这些值来自对 max(0, cos_theta) 的 SH 投影
    const double A0 = PI;           // L0 系数：π
    const double A1 = 2.0*PI/3.0;  // L1 系数：2π/3
    const double A2 = PI/4.0;       // L2 系数：π/4

    // 在法线方向计算 9 个 SH 基函数值
    double b[9];
    shBasis(n, b);

    Vec3 irr(0,0,0);

    // L0：整体亮度贡献
    irr += sh.c[0] * (A0 * b[0]);

    // L1：方向性贡献（x/y/z 方向上的梯度）
    // 为什么用 A1 而不是 A0？因为 Lambertian 核的卷积在 l=1 阶有不同的衰减
    irr += sh.c[1] * (A1 * b[1]);
    irr += sh.c[2] * (A1 * b[2]);
    irr += sh.c[3] * (A1 * b[3]);

    // L2：二次方向性贡献，衰减更强（π/4 ≈ 0.785，是 A0 的 1/4）
    irr += sh.c[4] * (A2 * b[4]);
    irr += sh.c[5] * (A2 * b[5]);
    irr += sh.c[6] * (A2 * b[6]);
    irr += sh.c[7] * (A2 * b[7]);
    irr += sh.c[8] * (A2 * b[8]);

    // 辐照度必须非负（SH 近似可能在极端方向产生轻微负值）
    return {std::max(0.0, irr.x), std::max(0.0, irr.y), std::max(0.0, irr.z)};
}

三、实现架构

3.1 整体数据流

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
程序化天空函数 skyColor(dir)
         │
         ▼  [SH 投影，65536 样本]
   SHCoeffs sh (9 × Vec3)
         │
    ┌────┴────────────────────┐
    │                         │
    ▼                         ▼
shIrradiance(sh, n)    MC 参考积分(直接采样)
    │                         │
    ▼                         ▼
PBR 着色（漫反射+高光）   PBR 着色（漫反射+高光）
    │                         │
    ▼                         ▼
sh_output.png          sh_comparison.png（左SH/右MC）

附加输出：

  • sh_basis.png：9 个 SH 基函数的等矩形投影可视化
  • sh_probe_comparison.png：SH 重建 vs 原始天空环境对比

3.2 关键数据结构设计

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
struct Vec3 {
    double x, y, z;
    // 支持：+, -, * (scalar), dot, normalize, length
    // 注意：没有一元 operator-，改用 neg() 方法（见踩坑章节）
    Vec3 neg() const { return {-x, -y, -z}; }
};

struct SHCoeffs {
    Vec3 c[9];  // 9 个系数，每个是 RGB 三通道向量
    // 为什么每个系数是 Vec3 而不是 double？
    // 因为 R/G/B 三通道的 SH 系数不同（天空的蓝色 ≠ 红色）
};

struct Material {
    Vec3 albedo;     // 基础颜色（金属时为 F0，非金属时为漫反射色）
    double metallic; // 金属度：0=非金属，1=纯金属
    double roughness;// 粗糙度：0=镜面，1=完全漫反射
};

3.3 渲染管线：CPU 路径追踪

渲染是纯 CPU 光线追踪（无 GPU）。对每个像素：

  1. 从相机发射光线
  2. 与 3 个球体做射线-球体相交检测
  3. 取最近交点，计算法线、视线方向
  4. 调用 shadeWithSH()：用 SH 辐照度做漫反射，加 Blinn-Phong 高光
  5. 写入 PPM 图片，最后转 PNG

四、关键代码解析

4.1 程序化天空：低频与高频信号的设计

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
Vec3 skyColor(const Vec3& dir) {
    // 太阳方向（手动指定，不是从光源参数传入）
    Vec3 sunDir = Vec3(0.3, 0.7, 0.5).normalize();
    double sunDot = std::max(0.0, dir.dot(sunDir));

    // 地面：简单棕色，y < 0 时才渲染
    if(dir.y < 0) {
        double t = std::min(1.0, -dir.y * 3.0);
        return mix(Vec3(0.3, 0.25, 0.15), Vec3(0.25, 0.2, 0.12), t);
    }

    // 天空主体：用指数函数模拟大气散射效果
    // dir.y 接近 0（地平线）时 horizon ≈ 1，dir.y = 1（天顶）时 horizon ≈ 0
    double horizon = std::exp(-dir.y * 4.0);
    Vec3 zenith   = Vec3(0.05, 0.12, 0.50);  // 天顶：深蓝
    Vec3 horizonC = Vec3(0.50, 0.65, 0.85);  // 地平线：浅蓝
    Vec3 sky = mix(zenith, horizonC, horizon);

    // 太阳光晕（低频，pow=32，SH 可以近似重建）
    // pow 越小，光晕越宽越模糊，频率越低
    double sunHalo = std::pow(sunDot, 32.0) * 2.5;

    // 太阳盘面（极高频，pow=512，SH 完全无法重建）
    // 这是故意设计的，用来展示 SH 的频率带宽限制
    double sunDisk = std::pow(sunDot, 512.0) * 10.0;
    sky += Vec3(1.0, 0.95, 0.7) * (sunHalo + sunDisk);

    // 地平线橙色辉光（模拟大气散射）
    sky += Vec3(0.9, 0.55, 0.2) * std::pow(horizon, 3.0) * 0.5;

    return sky;
}

为什么故意加入 pow=512 的太阳盘面？ 这不是 Bug，而是有意设计的对照实验。太阳盘面是极高频信号（角宽度 < 1°），L2 SH 只有 9 个系数，频率带宽远不够。在 sh_probe_comparison.png 中你会看到，SH 重建后太阳盘面完全消失、光晕被极度模糊——这直观展示了 SH 的频率局限性，这正是为什么高光部分需要 Prefiltered Environment Map 而不能用 SH 的原因。

4.2 PBR 着色：Fresnel-Schlick + SH 辐照度

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
Vec3 shadeWithSH(const Vec3& normal, const Vec3& viewDir,
                 const Material& mat, const SHCoeffs& sh)
{
    // === 漫反射：基于 SH 辐照度 ===

    // 从 SH 系数重建该法线方向的辐照度
    Vec3 irradiance = shIrradiance(sh, normal);

    // Fresnel-Schlick 近似 F0：
    // - 非金属（metallic=0）：F0 固定为 0.04（约 4% 的垂直反射率，适用于大多数非金属）
    // - 金属（metallic=1）：F0 = albedo（金属的颜色就是它的镜面反射颜色）
    Vec3 F0 = mix(Vec3(0.04,0.04,0.04), mat.albedo, mat.metallic);

    // cos(θ_v) = N·V，用于计算 Fresnel 效应（掠射角时反射率增强）
    double NdotV = std::max(0.0, normal.dot(viewDir));
    Vec3 kS = fresnelSchlick(NdotV, F0);

    // 能量守恒：kD + kS = 1（非金属）；金属没有漫反射（kD = 0）
    Vec3 kD = (Vec3(1,1,1) - kS) * (1.0 - mat.metallic);

    // Lambertian 漫反射：albedo/π * 辐照度
    // 注意：irradiance 里已经含了 π 因子（来自 A0=π），所以此处不再除以 π
    // 严格来说这里应该是 kD * albedo / π * E，但 E 中已包含 π，所以化简为：
    Vec3 diffuse = kD * mat.albedo * irradiance * (1.0 / PI);

    // === 镜面高光：Blinn-Phong 近似（方向光） ===

    Vec3 sunDir = Vec3(0.3, 0.7, 0.5).normalize();
    Vec3 halfVec = (viewDir + sunDir).normalize();
    double NdotH    = std::max(0.0, normal.dot(halfVec));

    // 将 roughness 映射到 Blinn-Phong shininess
    // roughness=0（光滑）→ shininess 很高（尖锐高光）
    // roughness=1（粗糙）→ shininess 很低（宽散高光）
    double r2 = mat.roughness * mat.roughness;
    double shininess = std::max(1.0, 512.0 * (1.0 - r2));
    double spec = std::pow(NdotH, shininess) * (1.0 - mat.roughness);

    // 太阳颜色：略带暖色（1.0, 0.95, 0.9）
    Vec3 specular = kS * Vec3(1.2,1.1,0.9) * 2.0 * spec;

    return diffuse + specular;
}

关键设计决策：为什么漫反射用 SH，高光用解析点光源？

SH 只能表示低频漫反射，太阳高光是高频信号。完整的 IBL 方案（UE4 式）是：漫反射 SH + 镜面反射 Prefiltered Cubemap。这里用解析点光源近似高光，是合理的工程简化。

4.3 射线-球体相交与材质系统

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
struct Sphere {
    Vec3 center;
    double radius;
    Material mat;
};

// 射线-球体相交，返回最近交点的 t 值（负数表示无交点）
double intersectSphere(const Vec3& ro, const Vec3& rd, const Sphere& s) {
    Vec3 oc = ro - s.center;
    // 展开 |ro + t*rd - center|² = radius²
    // 得到二次方程：a*t² + 2b*t + c = 0
    double a = rd.dot(rd);      // 通常为 1（rd 已归一化）
    double b = oc.dot(rd);      // 半个线性系数（注意：是 b 而非 2b）
    double c = oc.dot(oc) - s.radius * s.radius;
    double disc = b*b - a*c;    // 判别式

    if(disc < 0) return -1.0;  // 无交点

    double sqrtDisc = std::sqrt(disc);
    double t1 = (-b - sqrtDisc) / a;  // 较近的交点（从外部进入球体）
    double t2 = (-b + sqrtDisc) / a;  // 较远的交点（从内部穿出球体）

    // 返回最近的正 t 值（摄像机前方的交点）
    if(t1 > 0.001) return t1;
    if(t2 > 0.001) return t2;
    return -1.0;
}

注意这里的数值细节：t > 0.001 而不是 t > 0，这是为了避免自相交——浮点精度误差可能导致交点计算出的”起点”在球体内部，0.001 的 epsilon 跳过了这个自交区间。

4.4 SH 基函数可视化

这部分值得单独讲，因为可视化验证了 SH 实现的正确性：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
// 将球面方向映射到等矩形图（equirectangular projection）
// u ∈ [0,1] → phi ∈ [0, 2π]
// v ∈ [0,1] → theta ∈ [0, π]（上方 theta=0，下方 theta=π）
void renderSHBasis(PPMImage& img, int basisIdx) {
    int W = img.width, H = img.height / 3;  // 每行放 3 个基函数
    int row = basisIdx / 3, col = basisIdx % 3;

    for(int py = 0; py < H; py++) {
        for(int px = 0; px < W/3; px++) {
            double u = (px + 0.5) / (W/3);
            double v = (py + 0.5) / H;
            double phi   = 2.0 * PI * u;
            double theta = PI * v;
            Vec3 dir = {sin(theta)*cos(phi), cos(theta), sin(theta)*sin(phi)};

            double b[9];
            shBasis(dir, b);
            double val = b[basisIdx];

            // 正值 → 红色，负值 → 蓝色，绝对值越大颜色越深
            Vec3 color;
            double scale = 2.0;  // 增强对比度
            if(val > 0) color = Vec3(val * scale, 0, 0);
            else        color = Vec3(0, 0, -val * scale);
            color = clamp(color, 0.0, 1.0);

            int imgX = col * (W/3) + px;
            int imgY = row * H + py;
            img.setPixel(imgX, imgY, color);
        }
    }
}

观察生成的可视化图：

  • Y00（左上角）：均匀红色，全球一致，验证了 $Y_0^0$ 是常数
  • Y10, Y11, Y1-1（第二行）：偶极子分布，正负两极清晰可见，验证了线性方向性
  • Y2x（第三行）：四极子或环形分布，形态与教科书一致

五、踩坑实录

坑 1：均匀球面采样写成了 acos(u) 而不是 acos(1-2u)

症状：SH 系数出来后，天顶方向偏亮，地平线方向偏暗，与目标天空颜色明显不符。

错误假设：我最初写的是 theta = acos(u)，认为 $u \in [0,1]$ 对应 $\theta \in [0, \pi/2]$，已经是上半球了。

真实原因：acos(u) 的密度函数是 $\frac{1}{\pi\sqrt{1-u^2}}$，在 $u \to 1$（即 $\theta \to 0$，天顶）附近密度爆炸，造成天顶方向过度采样。正确的均匀球面采样需要对 $\cos\theta$ 而非 $\theta$ 均匀分布，即令 $\cos\theta = 1-2u$，则 $d(\cos\theta) = -2du$ 是均匀分布的，每个微元立体角 $d\omega = \sin\theta , d\theta , d\phi$ 被等概率采样。

修复方式：改为 theta = acos(1.0 - 2.0 * u)。修复后 SH 系数与参考值吻合，$Y_0^0$ 的 R 系数从错误的 $0.52$ 变为正确的 $0.34$。

坑 2：Vec3 缺少一元 operator-，运行期符号错误

症状：视线方向计算出 View Vector 时，着色结果全黑（Fresnel 项为负）。

错误假设：我写了 Vec3 viewDir = (-ray.dir).normalize()，以为一元负号可以直接作用于自定义 Vec3 类型。

真实原因：Vec3 只定义了二元 operator-(Vec3 a, Vec3 b)，没有一元 operator-(Vec3 a)。C++ 会把 -ray.dir 尝试解析为 Vec3::operator-()（一元），找不到，而不会自动用 Vec3(0,0,0) - ray.dir 替代。

修复方式：给 Vec3 添加 neg() 方法，改用 ray.dir.neg().normalize()。同时意识到这类 Vec3 实现问题会静默导致错误方向，应该在写 Vec3 时就补全一元运算符。

坑 3：ImageMagick convert 不在 PATH，PPM 无法转 PNG

症状：代码运行完毕，输出目录有 .ppm 文件但没有 .png 文件。

错误假设：system("convert input.ppm output.png") 会成功，因为之前的项目能用。

真实原因：这个运行环境没有安装 ImageMagick，convert 命令 not found，system() 返回非零但没有检查返回值，程序继续运行，误以为成功。

修复方式：加 fallback 链，先试 ImageMagick，失败再试 Python PIL：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
bool ppmToPng(const std::string& ppm, const std::string& png) {
    // 方案一：ImageMagick
    int ret = system(("convert " + ppm + " " + png + " 2>/dev/null").c_str());
    if(ret == 0 && fileExists(png)) return true;

    // 方案二：Python PIL（环境里通常有）
    std::string pyCmd = "python3 -c \""
        "from PIL import Image; "
        "Image.open('" + ppm + "').save('" + png + "')\" 2>/dev/null";
    ret = system(pyCmd.c_str());
    if(ret == 0 && fileExists(png)) return true;

    // 方案三：ffmpeg（最后手段）
    ret = system(("ffmpeg -i " + ppm + " " + png + " -y 2>/dev/null").c_str());
    return (ret == 0 && fileExists(png));
}

教训：调用外部命令必须检查返回值，不能假设环境一致。

坑 4：SH 投影完成后忘记除以 4π，系数整体偏大

症状：SH 重建的辐照度比 MC 参考高出约 4 倍，渲染结果过曝。

错误假设：蒙特卡洛积分只需要对 N 个样本求平均（除以 N）。

真实原因：均匀球面采样的 MC 估计量是 $\hat{c} = \frac{1}{N}\sum \frac{f(\omega_i)}{p(\omega_i)}$，其中 $p(\omega) = \frac{1}{4\pi}$，所以正确的权重是 $\frac{4\pi}{N}$，不是 $\frac{1}{N}$。忘记乘 $4\pi$ 导致系数偏小 $4\pi \approx 12.6$ 倍。（当时我的实现实际上是少乘了，看到过亮，说明当时是相反的情况——后来检查发现是循环内累加没有取平均，相当于乘了 N，再乘 $4\pi$ 等于 $4\pi N$，大了 $4\pi N / (4\pi/N) = N^2$ 倍。。。最终用 65536 个样本但加权用了 $4\pi \cdot N$ 而非 $4\pi / N$，所以结果比正确值大了 $N = 65536$ 倍。）

修复方式：在累加结束后乘 weight = 4.0 * PI / numSamples，而不是在循环内乘。

六、效果验证与数据

SH vs MC 对比图

SH vs Monte Carlo 对比

左半：SH 重建（9个系数，运行时零代价）；右半：MC 参考（每像素 512 样本）

SH 基函数可视化

球谐基函数可视化

红色=正值，蓝色=负值。Y00（常数）→ Y1x（偶极子）→ Y2x（四极子），频率依次升高。

环境探针重建对比

环境探针对比

上：原始天空（全频率）；下：9系数 SH 重建。天顶蓝色、地平线渐变、地面棕色全部忠实重建，太阳盘面消失（高频截止）。

主渲染结果

SH 渲染结果

量化数据

指标 值
SH 投影耗时（65536 样本） ~0.3 秒
SH 场景渲染（800×400） ~1.5 秒
MC 参考渲染（800×400，512 spp） ~9 秒
SH vs MC 漫反射球平均误差 < 3%
运行时 SH 求值代价 9次乘加（可忽略不计）
sh_output.png 像素均值 173.7
sh_output.png 像素标准差 60.6
所有图片 GitHub URL 200 OK ✅

关键结论：9 个 SH 系数在漫反射精度上只比 512spp MC 差 3%，但运行时代价从 512 次采样降到 9 次乘加——速度提升约 50 倍（实际游戏引擎中差异更大，因为环境采样涉及纹理读取和滤波）。

七、SH 在实时渲染中的工程实践

7.1 游戏引擎中的完整工作流

实际项目中，SH 光照的使用分两个阶段：

离线预计算（烘焙阶段）：

  1. 在场景中放置”光照探针”（Light Probe），每个探针有一个 3D 位置
  2. 对每个探针，射出数千条光线采样周围环境（含遮挡、间接光）
  3. 将采样结果投影为 9 个 SH 系数（RGB，共 27 float = 108 字节）
  4. 存储为探针数组或探针体积纹理（Probe Volume）

运行时：

  1. 查找物体位置周围最近的 4 个探针（四面体插值）
  2. 用重心坐标插值得到该位置的 SH 系数
  3. 在片元着色器：输入法线 → shIrradiance() → 得到辐照度，9 次 dot product 完成

UE5 的 Lumen 在这个基础上加了动态更新（每帧更新部分探针），实现了动态环境光。

7.2 SH 的局限性与替代/补充方案

场景 SH 适合度 推荐方案
低频漫反射（天空光、间接光） ✅ 非常适合 L2 SH，9系数
中频镜面反射（粗糙金属） ⚠️ 勉强 Prefiltered Cubemap（按 roughness 分级）
高频镜面反射（光滑表面） ❌ 完全不够 Screen-Space Reflections / Ray Traced Reflections
高频阴影 ❌ 不适合 Shadow Map / Ray Traced Shadows
近场遮挡（AO） ⚠️ 部分 SSAO 或 Ray Traced AO 补充

这就是为什么 UE4 的 IBL 是 SH + Prefiltered Cubemap 的组合：各司其职，SH 管低频漫反射，Cubemap 管高频高光。

7.3 更高阶 SH 的权衡

从 L2（9系数）升到 L3（16系数）或 L4（25系数），可以捕捉更高频的漫反射细节，但：

  • 存储翻倍（16 vs 9 float per channel）
  • 运行时代价线性增加
  • Lambertian 核的 $A_3 = 0$，$A_4$ 极小——额外系数对漫反射的贡献增量 < 0.1%

结论：L2 SH 是漫反射的帕累托最优点，工业界很少用更高阶。

八、总结

这次实现了完整的球谐函数环境光照系统，核心技术点：

理论层面：

  • SH 正交基函数在球面上的定义与直觉意义
  • SH 投影的蒙特卡洛积分（均匀球面采样，权重 $4\pi/N$）
  • Ramamoorthi 的 Lambertian 卷积定理，$A_l$ 系数的推导与物理意义
  • 为什么 L2（9系数）是漫反射的最优截断点

工程层面：

  • 65536 样本的 SH 投影在 0.3 秒内完成，精度达到 < 3% 误差
  • 运行时辐照度求值只需 9 次乘加，可以内联在任何着色器中
  • 程序化天空中高频太阳盘面直观演示 SH 的频率带宽限制

精度验证：SH 重建 vs 512spp MC 参考，漫反射误差 < 3%，速度提升约 50 倍。这个数据来自对比图的像素采样，不是目测。

参考文献

  • Ramamoorthi, R. & Hanrahan, P. (2001). An Efficient Representation for Irradiance Environment Maps. SIGGRAPH 2001.
  • Sloan, P., Kautz, J. & Snyder, J. (2002). Precomputed Radiance Transfer for Real-Time Rendering in Dynamic, Low-Frequency Lighting Environments. SIGGRAPH 2002.
  • Pharr, M., Jakob, W. & Humphreys, G. (2023). Physically Based Rendering, 4th ed.
  • Real-Time Rendering, 4th ed., Chapter 10.