实时全特效GI（Global Illumination，全局光照）一直是图形学中可望而不可及的圣杯，前人为此作出了无数努力，尽管诸多想法距离圣杯都还有一定的距离，但依然产生了不少有趣的技术。本文讨论的基于球谐函数（Spherical Harmonics）的预计算辐射传输（Precomputed Radiance Transfer）算法就是其中的一项明珠。
从渲染方程开始
考虑渲染方程：
\[L(x \to \Theta) = L_e(x \to \Theta) + \int_{\mathcal S^2}f_s(\Phi \to x \to \Theta)L(x \leftarrow \Phi)\cos\langle N_x, \Phi\rangle d\omega_\Phi\]$L_e(x \to \Theta)$没啥好讨论的，是直接可以获得的场景数据，这里简单起见就假设它不存在好了；我们再把非直接光照也扔掉，于是原式就变成了：
\[L(x \to \Theta) = \int_{\mathcal S^2}f_s(\Phi \to x \to \Theta)L_e(x \leftarrow \Phi)\cos\langle N_x, \Phi\rangle d\omega_\Phi\]我们进一步做一个假设：整个场景中只有一个环境光源。所谓环境光源，可以理解为一个包裹整个场景的、和物体距离非常之大（比如无限远）的球形光源。它的一个特征是从辐射亮度只和方向$\Phi$有关，和被照明点的位置$x$无关，因为不管$x$是多少，和光源的距离比起来也可以忽略。此时，渲染方程被简化为了：
\[L(x \to \Theta) = \int_{\mathcal S^2}f_s(\Phi \to x \to \Theta)V(x \to \Phi)L_E(\Phi)\cos\langle N_x, \Phi\rangle d\omega_\Phi\]其中$V(x \to \Phi)$表示从$x$点朝$\Phi$方向发射的射线是否没有被任何东西挡住，$L_E(\Phi)$是环境光源沿$-\Phi$方向的辐射亮度。
投影与重建
这里不太严格地介绍一下本文用到的投影与重建的思想。给定一组基函数$\beta = {\beta_0, \beta_1, \ldots}$和被投影的函数$f$，
\[c_i = \int f(x)\beta_i(x)dx\]可以用来衡量$f$与$\beta_i$的相似程度。将一系列的$\beta_i$用$c_i$加权求和，可以得到$f$的一个近似表示：
\[\hat f = \sum_{i=0}^n c_i\beta_i\]至于近似水平如何，就要看$f$性质好不好、$\beta_i$表现力如何以及$n$够不够大了。很多人大一学过的傅里叶级数就是非常好的示例。特别地，如果$\beta$满足如下性质：
\[\int \beta_i(x)\beta_j(x)dx = \begin{cases}\begin{aligned} &0, &i \ne j \\ &1, &i = j \end{aligned}\end{cases}\]就称$\beta$是一组单位正交基。单位正交基有一个超～棒的性质，令：
\[\begin{aligned} c_i &= \int f(x)\beta_i(x)dx \\ d_i &= \int g(x)\beta_i(x)dx \end{aligned}\]则：
\[\begin{aligned} &\int f(x)g(x)dx \\ \approx& \int \left(\sum_{i=0}^n c_i\beta_i(x)\right)\left(\sum_{i=0}^n d_i\beta_i(x)\right)dx \\ =& \sum_{i=0}^n c_id_i\int \beta_i(x)\beta_i(x)dx + \sum_{\begin{subarray}{c} i, j \in \{0, \ldots, n\} \\ i \ne j \end{subarray}} c_id_j\int \beta_i(x)\beta_j(x)dx \\ =& \sum_{i=0}^n c_id_i \end{aligned}\]PRT框架
看看上面的简化版渲染方程，再看看“单位正交基的美妙性质”所计算的积分，是不是觉得它们之间有一分相似？我们把$f_s$、$V$和$\cos$三项打包，把$L_E$单独打个包，分别看作$f$和$g$，就得到了：
\[L(x \to \Theta) = \int_{\mathcal S^2}f(\Phi \to x \to \Theta)g(\Phi)d\omega_\Phi\]现令：
\[\begin{aligned} t_k &= \int_{\mathcal S^2}f(\Phi \to x \to \Theta)\beta_k(\Phi)d\omega_\Phi \\ l_k &= \int_{\mathcal S^2}g(\Phi)\beta_k(\Phi)d\omega_\Phi \end{aligned}\]于是$L(x \to \Theta)$的近似版就出现了：
\[\hat L(x \to \Theta) = \sum_{k=0}^n t_kl_k\]这就是整个PRT的框架了——剩下的工作“只不过”是找到合适的$\beta$，以及计算$f$、$g$和$\beta$乘积的积分。PRT的思想就是尽可能地把更多的计算量挪到一劳永逸的预计算阶段，减小实时渲染阶段的负担。
球谐函数的定义
这里要介绍的球谐函数（Spherical Harmonics，SH）就是我们要使用的$\beta$。说实话，SH的定义实在是太蛋疼了，我觉得自己没法说清楚它的来龙去脉，这里直接摆上来吧。首先引入伴随Legendre多项式作为辅助函数：
\[\begin{aligned} &\text{for }l \in \mathbb N, m \in \mathbb N, m \le l, x \in [-1, 1] \\ &~~~~~~~~~~~~~~P_l^m(x) = \begin{cases}\begin{aligned} &x(2m + 1)P_m^m(x), &l = m + 1 \\ &(-1)^m(2m-1)!!(1-x^2)^{m/2}, &l=m \\ &\frac{x(2l-1)P_{l-1}^m(x) - (l+m-1)P_{l-2}^m(x)}{l-m}, &\text{otherwise} \end{aligned}\end{cases} \end{aligned}\]然后球谐函数长这样子：
\[\begin{aligned} &\text{for }l \in \mathbb N, m \in \mathbb Z, |m| \le l \\ &~~~~~~~~~~~~~~y_l^m(\theta, \phi) = \begin{cases}\begin{aligned} &\sqrt 2 K_l^m\cos(m\phi)P_l^m(\cos\theta), &m > 0 \\ &\sqrt 2 K_l^m\sin(-m\phi)P_l^{-m}(\cos\theta), &m < 0 \\ &K_l^0P_l^0(\cos\theta) \end{aligned}\end{cases} \end{aligned}\]其中，
\[K_l^m = \sqrt{\frac{(2l+1)(l-|m|)!}{4\pi(l+|m|)!}}\]是归一化系数。根据$l$的不同，可以把SH分成很多层，层数越高所代表的分量的频率也就越高。比如$l = 0$对应1个SH，就只是个常量函数；$l = 1$对应3个SH，有了一些粗粒度的变化；$l = 2$对应5个SH，代表了具有更高频率的分量，等等：
一般我们用“L阶SH”来代指前L层的SH。
以$(\theta, \phi)$为参数的球谐函数有时不那么好用，可以把它们换成$(x, y, z)$坐标。Wiki上有一些三维欧氏空间上的低阶球谐函数的表达式，可以直接拿过来用。最后，$l$和$m$两个参数终归不太好使，我们可以利用公式$k = l(l+1) + m$将它们编码成一个线性下标，以后记作$y^{(k)}(\theta, \phi)$。
求解SH系数
蒙特卡洛暴力即可，设在球面立体角上进行采样所使用的概率密度函数是$p_{\mathcal S^2}$，总采样数是$N$，估计量走起：
\[\begin{aligned} \hat t_k &= \frac 1 N \sum_{i=1}^N \frac{f_s(\Phi_i \to x \to \Theta)V(x \to \Phi_i)\cos\langle N_x, \Phi_i\rangle y^{(k)}(\theta_{\Phi_i}, \phi_{\Phi_i})}{p_{\mathcal S^2}(\Phi_i)} \\ \hat l_k &= \frac 1 N \sum_{i=1}^N \frac{L_E(\Phi_i)y^{(k)}(\theta_{\Phi_i}, \phi_{\Phi_i})}{p_{\mathcal S^2}(\Phi_i)} \end{aligned}\]当然，$p_{\mathcal S^2}$的选取还是可以讲究一番的，总之是按重要性采样的原则来，这里不多赘述。
旋转SH系数
SH函数还有个美妙的性质，就是我们可以通过旋转其系数来实现对原函数的旋转。这有什么用呢？可以用来高效地旋转我们的环境光。设想要是没有这个性质，那么环境光一转，就得用蒙特卡洛方法重新算$\hat l_k$；而有了这个性质，我们直接把一个线性变换作用到原来的$\hat l_k$上即可，岂不美哉。
这个系数的旋转方法用起来简单，推导起来就没那么简单了。许多paper长篇累牍地讲述如何推导SH系数的旋转矩阵，以及如何减少这一过程需要的计算量。我在这里找到一个非常⑨的做法，下面讨论的算法即来源于此。
首先我们需要一些SH的性质，它们的证明不在本文的讨论范围之内：
- SH系数旋转可以通过一个线性变换完成。
- 每一阶的SH系数旋转可以分别完成。
现在以$l = 2$的SH系数旋转为例进行说明。假设投影函数$P$将笛卡尔坐标表示的三维方向向量投影为2阶SH函数的5个系数，我们希望进行的旋转是个线性变换，用3阶方阵$M$表示，$R$是用来旋转SH系数的5阶方阵，则对任意方向向量$N$，有：
\[RP(N) = P(MN)\]现在随便选5个方向向量$N_0, \ldots, N_4$，于是：
\[R [P(N_0), \ldots, P(N_4)] = [P(MN_0), \ldots, P(MN_4)]\]记$A = [P(N_0), \ldots, P(N_4)]$，若$A$可逆，则$R$可以被表示为：
\[R = [P(MN_0), \ldots, P(MN_4)]A^{-1}\]于是我们得到了计算旋转矩阵$R$的方法——只需要适当地选取$N_0, \ldots, N_4$以保证$A$可逆即可。不过，在实际使用的时候我们往往并不关心$R$本身是多少，而是希望计算出某个球谐系数向量$x$对应的$Rx$。这个就非常好说了：
\[Rx = [P(MN_0), \ldots, P(MN_4)]A^{-1}x\]至此问题已经解决了。给定旋转矩阵$M$和SH系数$x$，旋转后的SH系数可按如下步骤计算：
- 预先选好$N_0, \ldots, N_4$，计算出$A^{-1}$。这一步可以直接硬编码结果。
- 用$M$旋转$N_0, \ldots, N_4$，得到$N_0’, \ldots, N_4’$。
- 用$P$计算出$N_0’, \ldots, N_4’$对应的SH系数，以它们为列构造一个5x5矩阵$S$。
- 计算$S(A^{-1}x)$并作为结果输出。
用这样的思路，可以推得任意阶SH系数的旋转方法。下面简单地罗列了一组前3阶的$N$和$A$：
l = 0
此时的SH是常量函数，旋转不会影响其系数值。
l = 1
结合Wiki上的SH函数公式，$P$可以表示为：
\[P\left([x, y, z]^T\right) = \frac 1 2\sqrt{\frac{3}{\pi}}\left[\begin{matrix} y/r \\ z/r \\ x/r \end{matrix}\right] ~~ \text{where } r = \sqrt{x^2 + y^2 + z^2}\]现取$N_0 = [1, 0, 0]^T, N_1 = [0, 1, 0]^T, N_2 = [0, 0, 1]^T$，则：
\[A^{-1} = \left[P(N_0), P(N_1), P(N_2)\right]^{-1} = \left(\frac 1 2\sqrt{\frac{3}{\pi}}\left[\begin{matrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{matrix}\right]\right)^{-1} = 2\sqrt{\frac{\pi}{3}}\left[\begin{matrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{matrix}\right]\]l = 2
首先根据SH给出$P$。设$r = \sqrt{x^2 + y^2 + z^2}$，则：
\[P\left([x, y, z]^T\right) = \left[ \sqrt{\frac{15}{4\pi}}\frac{xy}{r^2}, \sqrt{\frac{15}{4\pi}}\frac{yz}{r^2}, \sqrt{\frac{5}{16\pi}}\frac{-x^2-y^2+2z^2}{r^2}, \sqrt{\frac{15}{4\pi}}\frac{zx}{r^2}, \sqrt{\frac{15}{16\pi}}\frac{x^2-y^2}{r^2} \right]^T\]设$k = 1/\sqrt{2}$，取：
\[N_0 = [1, 0, 0]^T~~N_1 = [0, 0, 1]^T~~N_2 = [k, k, 0]^T~~N_3 = [k, 0, k]^T~~N_4 = [0, k, k]^T\]通过数值计算可以求得$A^{-1}$为：
\[A^{-1} = \left[\begin{matrix} 0 & k_0 & 0 & -k_0 & k_1 \\ k_0 & 0 & k_2 & -k_0 & k_0 \\ k_1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & k_1 & 0 \\ 0 & k_1 & 0 & 0 & 0 \\ \end{matrix}\right]~~\text{where } \begin{cases} k_0 = 0.91529123286551084 \\ k_1 = 1.83058246573102168 \\ k_2 = 1.5853309190550713 \end{cases}\]实现
我简单地实现了1~5阶SH的投影、系数旋转和重建过程，代码放在这里。随便画个模型看看：
环境光一开始在人像的右侧，上图中从左至右是将环境光绕垂直方向旋转了0度、90度、180度、270度后的结果，用了3阶SH，效果意外地很不错。
拓展
以下是一些更加高级的内容，我将它们实现在了代码中，但在这里不过多介绍。
在之前讨论PRT的时候，我们一开始就舍弃了间接光照，这多多少少会让物体表面看起来暗淡了一些，在材质的反照率较高时尤为明显。如果我们不把间接光照舍弃呢？
考虑将物体上某一点具有的亮度值$L(x \to \Theta)$视为该点对环境光的一种响应，那么对每个立体角微元$d\omega_\Phi$而言，$dL(x \to \Theta)/dL_E(\Phi)$都应该是一个和$L$无关、仅和$\Phi$有关的量$\mathcal I_{x \to \Theta}(\Phi)$。基于此，场景中只存在环境光源时的LTE可以被写作：
\[L(x \to \Theta) = \int_{\mathcal S^2}\mathcal I_{x \to \Theta}(\Phi)L_E(\Phi)d\omega_\Phi\]$\mathcal I_{x \to \Theta}(\Phi)$衡量了$L(x \to \Theta)$对来自$\Phi$方向的环境光的敏感程度。我们可以把$\mathcal I_{x \to \Theta}$和$L_E$分别投影到SH上，从而得到包含了间接光照的结果。将$\mathcal I_{x \to \Theta}$投影到SH空间的过程和path tracing非常相似，我就直接放结果了：
上图中最右侧的图像是用path tracing渲染的参考结果，上面的一排分别是用1~5阶SH投影和重建直接光照得到的结果，下面的一排则是考虑了间接照明的结果，可以“右键->在新窗口中查看图片”来查看大图。在阶数达到4、5阶时，下面的图像和参考图像间的差别已经很难用肉眼分辨了。