【cg】【球谐光照】PRT与PRTProbe
前言
书接前文【cg】【球谐光照】球谐函数。
上回书基本上讨论完了我们在实现本文所需要的理论知识，本文是对前面几篇理论知识的综合应用。
其中核心应用的是球谐函数的投影与重建，以及球谐函数的正交完备性和旋转不变性。
本文需要用到的在上回书中提到的公式列举如下。
\[ \left\{ \begin{array}{l} P_0(x) = 1 \\ P_1(x) = x \\ P_2(x) = \frac{1}{2}(3x^2 -1) \\ P_3(x) = \frac{1}{2}(5x^3 -3x) \\ \end{array} \right. \tag{4.6} \]
\[ \int_S Y_l^m(\theta, \varphi) Y_k^n(\theta, \varphi) d\omega = \delta_{lk} \delta_{mn} \tag{5.0} \]
\[ f(\theta, \varphi) = \sum_{k=0}^{k<n} c_k \cdot b_k \tag{5.3} \]
\[ c_k = \int_S f(\theta, \varphi) \cdot b_k(\theta, \varphi) d\omega \tag{5.5} \]
\[ C_l^{m'} = \sum\limits_{m = -l}^{m = l} c_l^m \cdot D_{m', m}^l(\alpha, \beta, \gamma) \tag{5.13} \]
PRT
预计算光照传输(Precomputed Radiance Transfer)[7]是指将渲染方程分为光照部分和传输部分，利用球谐函数的性质，分别预计算其球谐系数，以达到在运行时省略掉积分运算，快速计算出物体表面辐射照度的全局光照方法。
在PBR[8]里我们了解过对于法线方向为\(\vec{n}\)的某个物体表面\(p\)的渲染方程如下。
\[ L(p, \omega_0) = \int_{\Omega} L(p, \omega_i) \cdot f(p, \omega_0, \omega_i) \cdot (\vec{n} \cdot \omega_i) \cdot \mathrm{d} \omega_i \tag{6.1} \]
或
\[ L(p, \omega_0) = \int_{S} L(p, \omega_i) \cdot f(p, \omega_i, \omega_0) \cdot \max(0, (\vec{n} \cdot \omega_i)) \cdot \mathrm{d} \omega_i \tag{6.2} \]
其中\(L(p, \omega_i)\)为从\(\omega_i\)处输入的光源，即为上面所说的光照部分，其余的便为传输部分，记为\(T\)，则有
\[ T(p, \omega_i, \omega_0) = f(p, \omega_i, \omega_0) \cdot \max(0, \vec{n} \cdot \omega_i) \tag{6.3} \]
此时渲染方程可表示为
\[ L(p, \omega_0) = \int_{S} L(p, \omega_i) \cdot T(p, \omega_i, \omega_0) \cdot \mathrm{d} \omega_i \tag{6.4} \]
根据\((5.3)(5.5)\)，我们分别将两部分使用球谐函数展开，并求其系数。
光源部分
设其球谐系数为\(l_k\)，则有
\[ L(p, \omega_i) = \sum_{k=0}^{k<n} l_k \cdot b_k(\omega_i) \tag{6.5} \]
其中球谐系数由\((5.5)\)，需要在球面积分得到。
\[ l_k = \int_{S} L(p, \omega_i) \cdot b_k(\omega_i) \mathrm{d} \omega_i \tag{6.6} \]
传输部分
设传输部分的球谐系数为\(t_k\)，同理有
\[ T(p, \omega_i, \omega_0) = \sum_{k=0}^{k<n} t_k \cdot b_k(\omega_i) \tag{6.7} \]
同理有其系数也需要做球面积分。
\[ t_k = \int_{S} T(p, \omega_i, \omega_0) \cdot b_k(\omega_i) \mathrm{d} \omega_i \tag{6.8} \]
将\((6.5)(6.7)\)代入\((6.4)\)，可得
\[ \begin{array}{l} L(p, \omega_0) &=& \int_{S} \sum\limits_{p=0}^{p<n} l_p \cdot b_p(\omega_i) \cdot \sum\limits_{q=0}^{q<n} t_q \cdot b_q(\omega_i) \cdot \mathrm{d} \omega_i \\ &=& \sum\limits_{p=0}^{p<n} \sum\limits_{q=0}^{q<n} l_p \cdot t_q \cdot \underbrace{\int_{S} b_p(\omega_i) b_q(\omega_i) \mathrm{d} \omega_i}_{= \delta_{pq}} \\ &=& \sum\limits_{k=0}^{k<n} \sum\limits_{k=0}^{k<n} l_p \cdot t_q \cdot \delta_{pq} \\ &=& \sum\limits_{k=0}^{k<n} l_k \cdot t_k \end{array} \tag{6.9} \]
可以看到，上述推导利用了球谐函数的正交完备性，即公式\((5.0)\)，将渲染方程中的积分项消去了，使得结果只剩下光照部分的球谐系数与传输部分的乘积的累加值，如果我们能够预计算出这系数球谐系数，那么在运行时只需要做一个简单的点积即可得到最终的辐射照度值。
需要注意的是，漫反射(diffuse)的brdf是\(\frac{1}{\pi}\)，与观察方向\(\omega_0\)无关，所以\((6.8)\)也与观察方向无关，即无论在哪个观察方向，我们求得的\(t_k\)值是一样的。
而对于光泽反射(glossy)和镜面反射(specular)而言，其brdf与观察方向是有关的，不同的观察方向，求得的\(t_k\)值是不同的，即球谐系数\(t_k\)是关于观察方向\(\omega_0\)的函数。
\[ t_k(\omega_0) = \int_{S} T(p, \omega_i, \omega_0) \cdot b_k(\omega_i) \mathrm{d} \omega_i \tag{6.10} \]
在实际应用中，我们无法求出所有观察方向对应的\(t_k\)，所以会按一定的规则将观察方向离散化而有限个数，只对这有限个方向计算球谐系数，然后再通过在这有限个结果中插值得到任意观察方向的球谐系数或最终辐射照度的近似值。
PRT Probe
在实际应用中，我们无法对场景中所有的物体表面都将它们的光照的系数和传输的系数求出来，这种计算量是相当大的。
我们通常会在场景中选择一些代表，我们只计算这些代表处的球谐系数和光照结果，这些代表就是我们经常谈到的光照探针(light probe)。
如图所示，当我们需要点a处的光照结果时，我们会计算以probe所在位置为位置，以点a的法线方向为法线方向的假想表面的光照结果。计算点b的时候同理。
所以对于一个probe来说，我们需要计算朝向四面八方的无数个法线位置对应的光照结果，这样才有能力代表它周围朝向不同法线的物体表面。在实际应用中我们通常离散化出有限数量个法线方向进行计算，并通过法线方向进行插值得到任意法线方向的近似结果。
除此之外，由于一个probe周围不同的物体表面都有不同的材质属性，如果我们在计算光照的过程中需要使用到这些属性，比如计算glossy时需要粗糙度，比如计算specular时需要粗糙度、金属度等，probe无法决定使用哪位表面的材质属性进行计算，所以probe一般很难计算glossy和specular，一般我们只计算diffuse部分。
对于漫反射的光源部分的球谐系数(6.6)，因为它是与法线方向无关的，所以对于任意法线方向来说，光源部分的球谐系数都是一样的，总共只需要计算一次即可。
对于漫反射的传输部分的球谐系数，我们将其brdf值\(f(p, \omega_i, \omega_0) = \frac{1}{\pi}\)代入到\((6.3)(6.8)\)，得
\[ t_k = \frac{1}{\pi} \int_{S} max(0, \vec{n} \cdot \omega_i) \cdot b_k(\omega_i) \mathrm{d} \omega_i \tag{6.11} \]
我们当然可以离散化有限个法线方向分别使用上述公式求出每个法线方向对应的系数，但是利用球谐函数的旋转不变性，我们有一个更好的做法。
首先，我们将这个公式写回带\((l、m)\)的形式，即
\[ t_l^m = \frac{1}{\pi} \int_{S} max(0, \vec{n} \cdot \omega_i) \cdot Y_l^m(\omega_i) \mathrm{d} \omega_i \tag{6.12} \]
我们先只求法线方向为正z轴的球谐系数，即\(\vec{n_0} = (0, 0, 1)\)的球谐系数。显然此时\(\vec{n} \cdot \omega_i = \cos\theta_i\)，那么上述公式变为
\[ t_l^m = \frac{1}{\pi} \int_{S} \cos\theta_i \cdot Y_l^m(\omega_i) \mathrm{d} \omega_i \tag{6.13} \]
由于目标函数与方位角\(\varphi\)无关，所以上述公式与m无关，则上述公式可退化为
\[ \begin{array}{l} t_l^m = t_l^0 &=& \frac{1}{\pi} \int\limits_{S} \cos\theta_i \cdot Y_l^0(\omega_i) \mathrm{d} \omega_i \\ &=& \frac{1}{\pi} \int\limits_{\varphi = 0}^{\varphi = 2\pi} [ \int\limits_{\theta=0}^{\theta=\frac{\pi}{2}} \cos\theta \cdot Y_l^0(\theta) \cdot \sin\theta \mathrm{d} \theta ] \mathrm{d} \varphi \\ &=& \frac{1}{\pi} \sqrt{\frac{2l+1}{4\pi}} \int\limits_{\varphi = 0}^{\varphi = 2\pi} [ \int\limits_{\theta=0}^{\theta=\frac{\pi}{2}} \cos\theta \sin\theta \cdot P_l^0(\cos\theta) \cdot \mathrm{d} \theta ] \mathrm{d} \varphi \\ \end{array} \tag{6.14} \]
对于其他的任意法线方向\(\vec{n}\)，我们都可以通过旋转操作\(R(\alpha, \beta, \gamma)\)将\(\vec{n_0} = (0, 0, 1)\)旋转到这个法线方向上，那么此时，由上面讨论过的球谐函数的旋转不变性，可通过代入\((5.13)\)得新的球谐系数为
\[ \begin{array}{l} T_l^{m'} &=& \sum\limits_{m = -l}^{m = l} t_l^m \cdot D_{m', m}^l(\alpha, \beta, \gamma) \\ &=& \sum\limits_{m = -l}^{m = l} t_l^m \cdot D_{0, m}^l(\theta, 0, 0) \\ &=& \sum\limits_{m = -l}^{m = l} t_l^m \cdot [\sqrt{\frac{4\pi}{2l+1}} Y_l^m(\vec{n}) \cdot \delta_{m'm}] \\ &=& t_l^{m'} \cdot \sqrt{\frac{4\pi}{2l+1}} Y_l^{m'}(\vec{n}) \\ \end{array} \tag{6.15} \]
将\((6.14)\)代入\((6.15)\)，得
\[ T_l^{m'} = \frac{1}{\pi} Y_l^{m'}(\vec{n}) \int\limits_{\varphi = 0}^{\varphi = 2\pi} [ \int\limits_{\theta=0}^{\theta=\frac{\pi}{2}} \cos\theta \sin\theta \cdot P_l^0(\cos\theta) \cdot \mathrm{d} \theta ] \cdot \mathrm{d} \varphi \tag{6.16} \]
这就是传输部分旋转后的球谐系数的表达式，其中\(P_l^0(\cos\theta)\)为前面提到的勒让德多项式，如果我们代入勒让德多项式的表达式\((4.6)\)，那么\((6.16)\)的定积分的具体值是可以求出来的。
比如
\[ \begin{array}{l} T_0^{m'} &=& \frac{1}{\pi} Y_0^{m'}(\vec{n}) \int\limits_{\varphi = 0}^{\varphi = 2\pi} [ \int\limits_{\theta=0}^{\theta=\frac{\pi}{2}} \cos\theta \sin\theta \cdot \mathrm{d} \theta ] \cdot \mathrm{d} \varphi \\ &=& \frac{1}{\pi} Y_0^{m'}(\vec{n}) \cdot \frac{1}{2} \int\limits_{\varphi = 0}^{\varphi = 2\pi} \mathrm{d} \varphi \\ &=& \frac{1}{\pi} Y_0^{m'}(\vec{n}) \cdot \frac{1}{2} \cdot 2\pi \\ &=& \frac{1}{\pi} \cdot \frac{2\pi}{2} \cdot Y_0^{m'}(\vec{n}) \end{array} \]
再比如
\[ \begin{array}{l} T_1^{m'} &=& \frac{1}{\pi} Y_1^{m'}(\vec{n}) \int\limits_{\varphi = 0}^{\varphi = 2\pi} [ \int\limits_{\theta=0}^{\theta=\frac{\pi}{2}} \cos^2\theta \sin\theta \cdot \mathrm{d} \theta ] \cdot \mathrm{d} \varphi \\ &=& \frac{1}{\pi} Y_1^{m'}(\vec{n}) \int\limits_{\varphi = 0}^{\varphi = 2\pi} [-\frac{1}{3}\cos^3\theta]_0^{\frac{\pi}{2}} \mathrm{d} \varphi \\ &=& \frac{1}{\pi} Y_1^{m'}(\vec{n}) \cdot \frac{1}{3} \int\limits_{\varphi = 0}^{\varphi = 2\pi} \mathrm{d} \varphi \\ &=& \frac{1}{\pi} Y_1^{m'}(\vec{n}) \cdot \frac{1}{3} \cdot 2\pi \\ &=& \frac{1}{\pi} \cdot \frac{2\pi}{3} \cdot Y_1^{m'}(\vec{n}) \end{array} \]
至此，对于Probe的任意法线方向的漫反射，其光照部分的球谐系数只需要算一次积分\((6.6)\)，其传输部分的球谐系数甚至都不用算积分，直接使用公式\((6.16)\)求解定积分的值即可。
ue4本身封装了传输部分的球谐系数的函数，如下所示，可以看到，与我们计算出来的结果是一致的。
1 |
|
实现相关
光源部分球谐系数
由\((6.6)\)，我们使用有限数量次采样，使用蒙特卡洛积分将该定积分转换为有限数量次的均值求出近似值。
\[ l_k \approx \frac{1}{N} \sum \frac{L(p, \omega_i) \cdot b_k(\omega_i)}{pdf(\omega_i)} \]
当然也可以使用黎曼和，以步进极角和方向角的方式求这个定积分的近似值。
\[ l_k \approx \frac{1}{N_{\theta} N_{\varphi}} \sum^{N_\varphi} \sum^{N_\theta} L(p, \omega_i) \cdot b_k(\omega_i) \cdot \sin\theta_i \]
所以关于每一个积分项的代码如下 ，其中Radiance便是入射的光源\(L(\omega_i)\)。我们通常使用computer shader去计算这个累加结果。
1 | FThreeBandSHVector SHBasisCoefficients = SHBasisFunction3(TexelNormal); |
传输部分球谐系数
根据\((6.16)\)，我们不需要为计算传输部分的球谐系数做任何积分计算，只需要调用\((6.16)\)在ue中的实现即可。
1 | FThreeBandSHVector DiffuseTransferCofficients = CalcDiffuseTransferSH3(TexelNormal, 1) / PI; |
Irradiance
根据\((6.9)\)，只需要将上面两个球谐系数做一次点积，就可以得到最终的辐射照度结果。
1 | float3 Irradiance = max(half3(0, 0, 0), DotSH(LightCoefficients, DiffuseTransferCofficients)); |
显然，对于一个probe来说，我们可以求得任意一个法线方向对应的irradiance结果，但是如果我们要把所有的结果都存下来，那至少每个probe都需要一个cubemap，这显然是不现实的。我们这里分别尝试了两种方法来解决这个问题。
1.再次利用球谐函数将irradiance结果展开成球谐系数并只保存球谐系数，在采样时利用系数重建irradiance值。
该方法需要使用\((5.5)\)积分求系数。
该方法存储的是球谐系数，所以在场景绘制采样时需要采样更多次，如三阶球谐至少需要采样3次才可将9个系数采样出来。
2.只求有限个法线方向的irradiance结果，其它法线方向通过已知的法线插值出最终结果。
这里选取直角坐标系坐标轴上的6个方向(+x、-x、+y、-y、+z、-z)进行存储，其他法线方向的结果由这些方向插值而来。
我们可以对法线方向进行正八面体编码，以利用2D贴图的插值算法自动插值出结果，且只需要1次采样。
当然我们这里只个有6个顶点的值，所以最终一个probe存储的结果如下图所示。
我们将最终所有probe的结果按照ID保存在一张贴图上，最终生成的整个场景的irradiance结果如下所示。
