实时渲染-05：球谐与辐射传输

环境阴影

环境预积分解决了无阴影环境着色。加入可见性 $V(p,\boldsymbol\omega_i)$ 后，问题变得更复杂：天空各个方向都可能提供光，而每个点对这些方向的遮挡不同。

将环境离散为许多方向光，可以为每盏光建立 Shadow Map，再累加结果；方向数增加，阴影生成和查询成本也随之增加。改为每像素随机采样少量方向，则需要回答这些方向是否被挡住，采样不足还会产生噪声。环境图中有一个很亮的小窗口时，均匀采样甚至可能一直错过主要贡献。

实际实现常先为太阳或最亮的一两个方向计算阴影，再用较便宜的方法处理其余环境。环境遮蔽 Ambient Occlusion，AO 适合近似均匀环境下的局部遮挡，不能准确替代任意方向光照的阴影。Imperfect Shadow Maps 以较粗的阴影表示减小多光源成本，Lightcuts 按贡献聚合大量光源；实时光线追踪直接查询几何可见性。这些选择各有预算与误差来源。

这里沿另一条思路处理：如果场景固定，可否提前记录每个点对环境各方向的响应，之后只更换光照？在回答之前，先理解方向函数如何压缩。

频率与滤波

函数变化缓慢，主要由低频成分描述；细密纹理和尖锐变化需要更高频率。傅里叶展开 Fourier Expansion 将一个周期信号分解为不同频率的正弦、余弦之和，例如

$$f(x)=a_0+\sum_{k=1}^{\infty} \bigl[a_k\cos(kx)+b_k\sin(kx)\bigr].$$

每个系数决定对应频率的贡献。有限项重建只能保留其中一部分，截去高频通常使细节变少。

右图展示移到中心后的频谱幅度，中心附近表示缓慢变化的亮暗，离中心越远表示更快的变化。十字方向的强分量与图像中的方向结构、边界等有关。亮度显示通常还经过对数缩放，因此不能直接把频谱图像的显示亮度当成原图像素能量。

低通滤波 Low-Pass Filtering 保留低频、减弱高频。下图只保留中心一小片频谱，左侧头发、衣服和手部的细节随之变模糊，而大范围明暗仍然存在。

空间域的卷积 Convolution 写作

$$h(y)=\int f(x)g(y-x)\,\mathrm dx.$$

它将核 $g$ 移到位置 $y$，对附近的 $f$ 加权汇总。卷积定理 Convolution Theorem 给出 $\widehat h=\widehat f\widehat g$：若核的高频响应很小，乘法就削弱了输入的高频。因此“局部平均”和“频域低通”是同一个操作的两种观察方式。

光照积分中也有两个函数相乘再积分，但需要留意变量。$\int f(x)g(x)\,\mathrm dx$ 只是一个标量内积，不能直接谈它的输出频率。只有保留法线、观察方向等输出参数，再分析核随参数移动或旋转的关系，才可以说明积分形成了什么滤波。

基函数

基函数 Basis Functions 是一组用于组合其他函数的函数。取固定的 $B_i$，用不同系数表示不同信号：

$$f(x)\approx\sum_{i=0}^{N-1}c_iB_i(x).$$

正弦、余弦是一种选择，多项式 $1,x,x^2,\ldots$ 是另一种选择。基函数固定后，信号的存储便可由大量采样值改为少量系数。不过，只有 $f$ 位于有限基张成的空间中时，有限展开才完全准确；一般情况下是投影近似。

将函数当成向量，可以定义内积 $\langle f,g\rangle=\int fg$。若基函数正交归一 Orthonormal，满足 $\langle B_i,B_j\rangle=\delta_{ij}$，则

$$c_j=\langle f,B_j\rangle.$$

理由是把展开代入：$\langle\sum_i c_iB_i,B_j\rangle=\sum_i c_i\delta_{ij}=c_j$。投影所做的工作就是提取信号沿每个基方向的分量。非正交基也能使用，但需要通过 $G_{ij}=\langle B_i,B_j\rangle$ 构成的 Gram 矩阵解系数，不能直接套用这一积分。

球谐函数

环境光按方向变化，每个单位方向都对应单位球面上的一点。因此需要球面上的基函数。球谐函数 Spherical Harmonics，SH 按角向变化的快慢分带，可以看成球面方向信号的频率展开。

记实 SH 为 $Y_l^m(\boldsymbol\omega)$，$l=0,1,2,\ldots$ 是频带，$m=-l,\ldots,l$ 区分同一带中的函数。第 $l$ 带有 $2l+1$ 项，前 $n$ 带共有

$$\sum_{l=0}^{n-1}(2l+1)=n^2$$

项。包含 $l=0,1,2$ 的截断有 $1+3+5=9$ 项，常称二阶 SH；“前五个频带”则对应 $25$ 项。写代码时应明确阶数是否包含最高带，避免差一带。

$l=0$ 在所有方向保持常数；$l=1$ 已能表达相反方向的变化；更高带的正负区域分布更细。两种颜色表示函数的正负，离球心的距离表示幅度大小。它们不是材质颜色或真实几何，也不能只凭“花瓣数”确定任意函数的频谱。

球面方向用极角 $\theta$、方位角 $\phi$ 表示，积分测度为 $\mathrm d\omega=\sin\theta\,\mathrm d\theta\,\mathrm d\phi$。正交归一关系定义在完整球面 $S^2$：

$$\int_{S^2}Y_l^mY_{l'}^{m'}\,\mathrm d\omega =\delta_{ll'}\delta_{mm'}.$$

于是环境光系数为

$$L_l^m=\int_{S^2}L_i(\boldsymbol\omega)Y_l^m(\boldsymbol\omega)\,\mathrm d\omega, \qquad L_i\approx\sum_{l,m}L_l^mY_l^m.$$

RGB 光照对三个通道分别投影。用 Cube Map 离散积分时，需要乘每个纹素代表的立体角；六个面上的等面积纹素并不对应相等立体角。

漫反射卷积

Lambertian 表面的环境反射为 $L_o=\rho E(\mathbf n)/\pi$，其中

$$E(\mathbf n)=\int_{S^2}L_i(\boldsymbol\omega) \max(0,\mathbf n\cdot\boldsymbol\omega)\,\mathrm d\omega.$$

背面方向通过 $\max(0,\cdot)$ 归零，所以这和接收半球上的积分相同。对某一法线，响应是以法线为轴、从正前方到切向逐渐下降的截断余弦核 Clamped-Cosine Kernel。转动法线，相当于在球面上转动这一平滑核；环境中的小细节被宽范围平均，漫反射响应因此比原环境平滑。

对这种只依赖夹角的核，球面卷积在 SH 域中变成逐频带缩放：

$$E_l^m=A_lL_l^m, \qquad A_l=2\pi\int_0^1\mu P_l(\mu)\,\mathrm d\mu.$$

$P_l$ 是 Legendre 多项式，例如 $P_0=1$、$P_1=\mu$、$P_2=(3\mu^2-1)/2$；$\mu$ 是入射余弦。代入可得

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

同一频带的所有 $m$ 乘相同 $A_l$，因为核绕自身轴旋转不变。高带响应总体较弱，前三带适合近似漫反射；$A_4$ 仍非零，说明这并不是严格的九项精确表示。

九项近似

只保留 $l=0$，环境仅剩整体平均，无法区分不同法线。

增加 $l=1$，四项近似已经表现出上下、左右方向的宽范围亮暗和颜色变化，但局部形状仍有差别。

保留到 $l=2$，九项近似进一步恢复中等尺度结构。

这组三幅图对应的示例 RMS 误差约为 $25\%$、$8\%$、$1\%$。课件另引用 Basri 与 Jacobs（2001）关于不同照明下平均误差小于 $3\%$ 的结果，与这组示例的 RMS 数字应分开理解；平均误差不等于任意像素都满足同一误差上界。强方向性光照和截断重建中的负值，都需要结合具体应用评估。

九项实 SH 是常数、一次项和二次项的组合，因此 $E(\mathbf n)$ 可以整理为

$$E(\mathbf n)=\mathbf n_h^{\mathsf T}M\mathbf n_h, \qquad \mathbf n_h=(n_x,n_y,n_z,1)^{\mathsf T}.$$

$M$ 是由环境 SH 系数和上述卷积系数组成的对称 $4\times4$ 矩阵，每个颜色通道各一份。这里末尾的 $1$ 用来容纳常数项和一次项，与位置投影无关。

// M 已由九项环境 SH 系数构造；输入是单位法线。
float evaluateIrradiance(vec3 normal, mat4 M) {
    vec4 direction = vec4(normal, 1.0);
    return dot(direction, M * direction);
}

这仍是无阴影漫反射环境着色。若有一个遮挡物挡住半边天空，每个点的响应核便不再只是转动的余弦，不能仅靠这份通用矩阵处理。

传输预计算

预计算辐射传输 Precomputed Radiance Transfer，PRT 将固定场景如何响应环境光提前保存。将渲染方程写成光照与光传输 Light Transport 的乘积：

$$L_o(p,\boldsymbol\omega_o)=\int_{S^2} L_i(\boldsymbol\omega_i)T(p,\boldsymbol\omega_i,\boldsymbol\omega_o) \,\mathrm d\omega_i.$$

对于只有一次表面反射的 Lambertian 接收点，可取

$$T(p,\boldsymbol\omega_i)=\frac{\rho(p)}\pi V(p,\boldsymbol\omega_i)\max(0,\mathbf n\cdot\boldsymbol\omega_i).$$

$T$ 包含该点的材质、朝向和遮挡，背面方向为零。环境可变化，固定场景的 $T$ 保持不变。朴素地使用 $6\times64\times64$ Cube Map，每点要处理 $24\,576$ 个方向贡献；将光照代入 SH 展开后，得到

$$L_o(p)\approx\sum_iL_i^{\rm SH} \underbrace{\int_{S^2}B_i(\boldsymbol\omega)T(p,\boldsymbol\omega)\,\mathrm d\omega}_{t_i(p)}.$$

每个 $t_i$ 是场景对第 $i$ 个基函数照明的响应，可以提前计算。渲染时只需 $\sum_iL_i^{\rm SH}t_i$。这里将 $\rho/\pi$ 包含在传输中；若实现把反照率放在点积外，预计算定义也要相应改变。

还可以同时将 $T=\sum_jt_jB_j$ 展开，此时

$$\int_{S^2}\sum_iL_i^{\rm SH}B_i\sum_jt_jB_j\,\mathrm d\omega =\sum_{i,j}L_i^{\rm SH}t_j\delta_{ij} =\sum_iL_i^{\rm SH}t_i.$$

看似需要 $N^2$ 项的双重求和，通过正交性只剩 $N$ 项。物理接收半球已包含在置零的 $T$ 中；不能把 SH 的正交关系直接改写为半球上的关系。

基照明响应

传输系数的直观意义，是用一组已知照明分别照场景，记录每个点的结果。完整投影可以逐点追踪许多方向的可见性，也可以离线计算更复杂的反射路径。

左边三幅球面分别描述三个带正负值的基照明；右边是模型在相应照明下的响应。红蓝用来显示系数的正负，真实材质并没有变成红蓝。基函数允许负值，因为它们是线性展开的中间量；最终物理环境是许多基函数加权组合的结果。

// 离线：保存固定场景对每个基函数的响应。
for each shading point p:
    transport[p] = zeros(N)
    for each sampled direction omega with density pdf:
        response = albedo[p] / PI
                 * visibility(p, omega)
                 * max(dot(normal[p], omega), 0)
        for i in 0 .. N - 1:
            transport[p][i] += basis(i, omega) * response / pdf
    transport[p] /= sampleCount

// 每帧：环境已投影为 lightSH，逐颜色通道计算。
for each shading point p:
    outgoing[p] = dot(lightSH, transport[p])

可逐顶点存储传输，着色后插值颜色，也可插值传输系数再在片元处求值。漫反射点积对传输是线性的，因此固定光照下两者在线性插值意义上相容；法线、材质或高频遮挡变化仍会受网格与插值分辨率限制。

固定与旋转

PRT 将成本移到离线阶段，代价是固定条件。几何位置、法线、遮挡物或材质发生变化，原传输便可能失效；只换环境光照则只需替换光照系数。每点 $25$ 项的漫反射传输仍有存储成本，大场景需要压缩或分区。

SH 的一个有用性质是旋转封闭：转动环境后，第 $l$ 带的 $2l+1$ 个系数通过一个小矩阵相互混合，不会跨到其他频带：

$$\mathbf L_l'=D_l(R)\mathbf L_l.$$

$R$ 是旋转，$D_l$ 是该频带对应的旋转矩阵。旋转不会增加频率，因此可以直接处理系数，无需每次重采样整张环境图。对于孤立的刚性物体，也可以将世界光照转到物体局部坐标中查询；物体与其他物体的相对位置变化会改变遮挡，不能仅靠转动系数解决。

多次反弹

固定场景中，光照强度与最终响应保持线性。传输可以包含环境光经过多次反射、再到达当前点的贡献，运行时仍以当前环境系数和已存响应组合。

左图不计可见性，腋下和褶皱偏亮；中图加入阴影，遮挡区域变暗；右图再加入 Interreflection，其他表面反弹来的光使部分暗处重新获得照明。环境光本身是直接照明，只有传输包含表面间反弹，结果才包含相应间接光。

路径记号可以明确预计算范围：$L$ 表示光源，$E$ 表示观察端，$D$ 表示漫反射事件，$G$ 表示 Glossy，$S$ 表示理想镜面；$|$ 表示选择，$*$ 表示重复零次或多次。$LGE$ 是一次 Glossy 反射，$L(D|G)^*E$ 允许漫反射或 Glossy 的多次反弹，$LS^*(D|G)^*E$ 还允许前面的镜面路径，并可能形成 Caustics，即光经反射或折射集中后的亮区。

纳入更多路径通常增加离线计算难度和数据表示要求。运行时无需逐条重走这些已编码路径，但并不意味着预计算、内存或更新成本消失。半透明、毛发、动态场景和动态材质都需要扩展传输表示，基础固定场景 PRT 不直接覆盖这些变化。

光泽传输

漫反射出射辐亮度与观察方向无关，传输可存为一个向量。Glossy 反射有宽度有限的方向瓣，随着观察方向变化会看到不同响应；传输同时依赖入射方向和出射方向，是两个球面方向构成的四维函数。

先沿入射方向展开，再沿出射方向展开：

$$T(\boldsymbol\omega_i,\boldsymbol\omega_o) \approx\sum_i\sum_jt_{ij}B_i(\boldsymbol\omega_i)B_j(\boldsymbol\omega_o).$$

将光照代入，对入射方向积分，得到出射系数

$$o_j=\sum_iL_i^{\rm SH}t_{ij}, \qquad L_o(\boldsymbol\omega_o)\approx\sum_jo_jB_j(\boldsymbol\omega_o).$$

图中一排色块是一个方向函数的系数向量，方阵记录“每个入射基如何贡献到每个出射基”。矩阵乘法先产生出射方向函数，之后还要在实际观察方向上重建；不能把输出系数向量本身当成像素 RGB。

$N$ 个基函数需要每点约 $N^2$ 个传输数：$16$ 项对应 $16\times16=256$ 项，运行时也从点积转为矩阵乘法。再考虑表面位置的两个自由度，局部传输可视为六维量；各向异性、空间变化 BRDF 可以纳入，但存储和预计算相应增加。

接近镜面的材质反射范围窄，需要更高角向频率。提高 SH 阶数可以改善表示，代价是矩阵迅速增长；也可以换基函数或直接追踪近镜面路径。低阶 SH 的不足是截断精度有限，不是 SH 数学上无法表示高频。

小波表示

小波 Wavelets 使用不同位置与尺度上的局部基函数。与覆盖整个球面的 SH 相比，局部变化可以集中在少量系数中。课程以 Haar Wavelet 为例，在 Cube Map 的各个面上进行二维变换。

一维的一对值 $a,b$ 可改写为平均分量 $s=(a+b)/\sqrt2$ 与差异分量 $d=(a-b)/\sqrt2$。二维对子块沿两个轴处理，得到一个低频子块和水平、垂直、对角细节子块，再递归分解低频子块。低频逐级集中到左上角，其他子块记录不同尺度的细节；只有保存全部系数时，才能通过逆变换完整恢复输入。

白色为正、黑色为负，灰色为零。左上常数基记录整体平均；其余基只在特定区域取非零值，水平、垂直或交叉的正负分区记录不同方向的差异。变小、移位的分区使局部细节不必影响整张图。

若大部分系数很小，可以只保留绝对值较大的项。变换本身是线性的，但“根据输入选出最大系数”属于非线性近似 Nonlinear Approximation。课程示例中保留约 $0.1\%$–$1\%$ 的项，不能将这一比例当成任意环境都达到相同精度的保证。计算内积时，要保持光照与传输基索引及归一化一致，不能把各自保留的系数按数组顺序直接点乘。

这张图截取一个完整 Cube Map 面的系数布局。递归细分的左上区域保存较粗的结构，较大子块记录细尺度差异；大部分位置接近灰色的零值，只有少量位置有明显响应。它已经是系数图，不能按原教堂空间位置直接解释色块。

右图茶壶上的反射结构和阴影细节比左图更清楚，说明局部稀疏表示可以在有限预算下保留更高频率。图中“All frequency”描述的是方法可保留多尺度细节，并不保证有限系数重建绝对无误差。

小波的旋转不像 SH 那样方便：转动环境后，局部基所在的位置和尺度关系改变，往往需要先恢复 Cube Map，再旋转、变换和筛选系数。标准 JPEG 主要使用离散余弦变换 DCT，不能因二者都有变换压缩步骤就把它称为 Haar 小波压缩。

其他表示也有各自取舍：Spherical Gaussian 以球面上的平滑瓣拟合方向分布，Zonal Harmonics 使用绕轴对称的球谐成分，分片常数表示则直接对方向区域分区。进一步把光照、BRDF、可见性分别展开，需要处理多个基函数乘积的积分，预计算量和动态自由度都会改变。选择表示时，应同时考虑频率、旋转、存储与哪些场景量需要实时变化。