四 实时环境光照 Real-Time Environment Mapping¶
本课程为闫令琪老师的 GAMES202-高质量实时渲染
此为个人的课堂笔记，如有疏漏，请多指正
1 环境光照着色 Shading from Environment Lighting¶
基于图像的渲染 Image-Based Lighting (IBL)：通过 Spherical map 或者 Cube map 这些环境光照图直接计算环境光照并给定点着色的技术
计算实时环境光照的关键在于：已经有了环境光照，即已经有了 Based Image，如何以最快的速度确定某个着色点应该如何着色
由于不论何种光照，计算着色都需要解渲染方程 Rendering equation，因此如何以近似的方法，在尽量保证质量的情况下，降低解渲染方程的复杂度，便是实时光照需要重点解决的问题
在暂时不考虑阴影问题（即 可见性\(V\)）时，要解决的渲染方程如下：
1）渲染方程的通用解法¶
蒙特卡洛积分 Monte Carlo integration，详见 GAMES101-光线追踪（蒙特卡洛积分与路径追踪）
问题：需要大量的样本才能将结果收敛到期望值，开销过大，速度很慢
因此，希望找到一个避免大量采样的方法来减少开销
2）环境光照渲染方程的近似解法 Spilt Sum 分离求和法¶
Spilt Sum 指的是将一个大的和式近似地拆解为多个和式计算，从而降低计算复杂度
（1）对 光照部分\(L\) 进行预滤波 PreFilter¶
已知现象：
- 如果 BRDF项 的光泽度较高 glossy，那么积分时的积分率（支持集）support 就会特别小
- 如果 BRDF项 的漫反射程度较高 diffuse，那么积分函数就会特别平滑 smooth
根据以上现象，可以利用“实时阴影-2 阴影映射技术背后的数学”中提及的用于估算定积分的经典的近似不等式，将 光照项\(L\) 从积分中提到外面来，从而将渲染方程近似为以下方程：
其效果相当于将 光照项\(L\) 的特定区域进行滤波 Filter，具体到环境光照所使用的 Spherical map 或者 Cube map，就是将对应环境光照图进行了模糊处理。并且这个 Filter 操作可以提前进行，即 PreFilter。甚至可以进行多次滤波核大小不同的 PreFilter 得到多张环境光照图，在最终使用时只需对两张环境光照图进行三线性插值即可，类似于 MipMap 的原理
Spilt Sum 的第一步最终实现的效果就是将多次的采样变为一次对 理想镜面反射方向\(r\) 的直接采样：
（2）对 BRDF 项的所有可能组合情况进行预计算 Precompute¶
以下内容以微表面模型 Microfacet Model 作为例子进行讨论。微表面模型的具体内容可参考 GAMES101-材质与外观（微表面模型）
a）初始情况¶
微表面模型的 BRDF 如下：
其中 \(F\)、\(G\) 和 \(D\) 分别是菲涅尔方程、几何遮蔽函数和法线分布函数。其中 \(h\) 是介于 \(\omega_i\) 和 \(\omega_o\) 之间的半程向量 (half vector)
其中 菲涅尔项\(F\) 如下：
在该 BRDF 模型上直接对所有可能组合情况进行预计算，需要考虑 RGB、roughness、入射角度\(\omega_i\)（这里表现为半程向量 \(h\)） 共五维的情况来打表，计算和存储量为海量，无法实现
b）对 菲涅尔项\(F\) 和 微表面分布函数\(D\) 的近似¶
对于 菲涅尔项\(F\) 可以进行石里克近似 Schlick's approximation：
其中，\(R_0=\left(\frac{\eta_1-\eta_2}{\eta_1+\eta_2}\right)^2\)是基底颜色（基础反射率），其中 \(\eta_1\) 和 \(\eta_2\) 是着色点所处表面两侧介质的折射率
对于 微表面分布函数\(D\) 的分布可以用 Beckmann 模型描述：
其中 \(\alpha\) 定义了 roughness，\(h\) 是介于 \(\omega_i\) 和 \(\omega_o\) 之间的半程向量（可以被认为与 入射角\(\omega_i\) 相关）
通过以上近似，就可以将最开始的五维降低为三维：基底颜色\(R_0\)、粗糙度roughness、入射角度\(\omega_i\)
c）将 菲涅尔项\(F\) 的 基础反射率\(R_0\) 提取到定积分外部¶
通过如下变换将 \(R_0\) 提取到定积分外部：
于是定积分的部分不再依赖于 基础反射率\(R_0\)，在预计算时只需要考虑 粗糙程度roughness 和 入射角度\(\theta\) 即可
通过以上过程将一个五维的预计算降低至二维，使得其消耗很小，从而能在实时光照中实现对 BRDF项 的所有可能组合情况进行预计算
综上，通过 Spilt Sum 分离求和法成功实现对环境光照渲染方程的近似求解，最终实现的效果对比如下：
2 环境光照阴影 Shadow from Environment Lighting¶
如果在着色时需要考虑环境光的可见性，则实现实时环境光照是一个相当困难的问题：
- 因为环境光照来自四面八方，所以如果把环境光照下的着色看作多光绘制 (many-light rendering)，则每一个光源都要生成一张阴影图，那么阴影图的数量将极其庞大
- 而如果把该问题看作抽样问题，则环境光照的 可见性项\(V\) 可能是任意复杂度，不能根据 Split Sum 方法估计积分结果。因为 光照项\(L\) 的支撑集是整个半球（即 support 很大，因为是四面八方的环境光）、BRDF项 也可能并不平滑（在 光照项\(L\) 的 support 很大的情况下想要用 Split Sum 方法则必须保证另一项 BRDF 的积函数十分平滑，然而这并不一定，比如材质很 glossy 的时候），以上约束导致并不能利用 Split Sum 方法将 可见性项\(V\) 拆解到定积分的外面
因此工业界的一般解决方案是选取环境中最亮的那个光源（例如太阳）或前几个光源生成阴影图，然后由此生成阴影
3 预计算辐射亮度传输 Precomputed Radiance Transfer (PRT)¶
之前提到环境光照阴影是很难在实时实现渲染的，但是 预计算辐射亮度传输PRT 技术是这个问题的一个解决方案
在了解 PRT 前需要了解的前置基础知识：
- 频域 Frequency Domain 和滤波 Filtering
- 球谐函数 Spherical Harmonics 及其基函数 Basis functions
1）频域 Frequency Domain 和滤波 Filtering¶
详细内容参考 GAMES101-光栅化（深度测试与抗锯齿）
在此基础上对于渲染方程的定积分的一个更通用的理解：
两个函数 \(f(x)\) 和 \(g(x)\) 相乘的结果再进行积分，这个操作可以认为是一个卷积/滤波。而 \(f(x)\) 和 \(g(x)\) 相乘后再积分可以认为是时域上两个信号相乘后再卷积，其结果就是频域上的两个信号相乘（参考卷积定理 Convolution Theorem）
而如果这两个频域上的信号有一个是低频 Low frequency 的，则它们相乘的结果也是低频的，也就对应着卷积后的结果是平滑 smooth 的。积分后结果的频率取决于两个信号中更低的那个的频率（两个频谱图相乘，如果一个地方没信号，相乘的结果当然也就没信号了）
2）球谐函数 Spherical Harmonics 及其基函数 Basis functions¶
（1）基函数 Basis functions¶
一个函数 \(f(x)\) 可以描述为其他一系列函数的线性组合，这些组合的函数就称为函数 \(f(x)\) 的基函数 Basis functions
例如，傅里叶级数展开的一系列函数就是一套基函数
（2）球谐函数 Spherical Harmonics¶
球谐函数 Spherical Harmonics：指的是一系列定义在球面上的二维基函数
不同阶的 SH基函数 如下：
每一阶 SH 有 2L+1 个基函数，前 n 阶 SH 共有 n^2^ 个基函数。L 越大则代表频率越大
对于一个二维球面函数，如果想要展开成一系列 SH基函数 的线性组合，计算每一个基函数前面的系数 \(c_i\)（这个过程可以被称为投影，即将一个二维球面函数投影到每一个 SH基函数 上）的公式如下：
注：投影操作一般是通过点乘实现的，而这个积分的本质其实就是点乘
这样就可以得到每一个二维球面函数关于 SH基函数 的线性组合表达式了。不过一般情况下，为了避免过于复杂的计算量，并且保证实现效果的前提下，只会使用到前4阶 (L = 0~3) 的 SH基函数
（3）球谐函数的一些重要性质¶
- orthonormal（正交性）：SH是一组正交基，SH的一个基函数投影到另一个任意基函数上的值为0：\(\int_{\Omega}B_i(i)\cdot B_j(i)\mathrm{d}i=1(i=j);\int_{\Omega}B_i(i)\cdot B_j(i)\mathrm{d}i=0(i\ne j)\)
- simple projection/reconstruction（投影运算简单）：\(c_i=\int_{\Omega}f(\omega)B_i(\omega)\mathrm{d}\omega\)
- simple rotation（旋转运算简单）：任何一个被一组SH基函数线性组合表示的二维球面函数发生旋转时，可以等价为对每一个基函数进行旋转，而每一个基函数旋转后的结果，都可以被其同阶的基函数的线性组合表示（因此可以很方便地进行预计算——提前打表）
- simple convolution（卷积运算简单）
（4）球谐函数的使用实例：环境光下的漫反射着色计算¶
UC San Diego 的 Ravi Ramamoorthi 教授的博士论文：Precomputation-Based Rendering (Ravi Ramamoorthi)
主要思想：当计算漫反射时，渲染方程的 BRDF项 仅需要前三阶共9个 SH基函数 即可较为准确地描述出来，即只需要前三阶的低频信息即可。而考虑到渲染方程的积分可以视为频域上的两个信号 \(L\)项 和 BRDF项 相乘，若其中的 BRDF项 为低频的，则结果也是低频的，此处的 BRDF项 就类似于一个低通滤波器，因此对于 光照\(L\)项 也就根本不需要高频信息来描述，仅需要前三阶的低频信息即可
在 BRDF项 仅使用前3阶 SH基函数 描述的同时，对 光照\(L\)项 也仅使用前3阶 SH基函数 描述的实现效果：
综上所述：
- 当想要尽量完整地描述光照信息时，使用足量的 SH基函数 进行描述
- 当只需要低频的信息时（比如漫反射下的 BRDF项 和与其关联的 \(L\)项），使用低频的 SH基函数 描述即可实现
有关球谐函数的两个个比较易懂的解释参考：
球谐函数介绍（Spherical Harmonics） - 知乎 (zhihu.com)
球面高斯介绍（Spherical Gaussian） - 知乎 (zhihu.com)
3）预计算辐射亮度传输 Precomputed Radiance Transfer (PRT)¶
对于环境光照计算的完整渲染方程如下：
如果想要计算这个积分，则对于每一个着色点 shading point，需要计算6（Map的六个图） * 64 * 64（分辨率）次，计算量极为庞大
预计算辐射亮度传输 PRT 的核心思想是假设场景中的只有光照会发生变化，将绘制方程的被积函数拆分成光照 (lighting) 和光线传输 (light transport) 这两部分，并分别预计算两者的纹理图像，再把图像从空间域（时域）转换到频域，最终把着色时绘制方程中的定积分计算转换成向量的点积 (Diffuse Case)，或者转换成向量与矩阵的乘法 (Glossy Case)
（1）预计算光照项 (lighting) 和光线传输项 (light transport)¶
对于光照项 Lighting，可以用 SH基函数 的线性组合来近似表示
对于光线传输项 Lighting transport，在场景只有光照情况会发生改变时，可以被视为一个二维球面函数，也可以用 SH基函数 的线性组合来近似表示
即，将光照项 Lighting 和光线传输项 Lighting transport 的信号从时域（空间域）转换到频域：
（2）PRT 的 Diffuse Case¶
A）Diffuse Case 下对渲染方程的第一种拆解方式¶
在漫反射模型中，BRDF项 可以认为是一个 常数\(\rho\)，与光线 入射方向\(\omega_i\) 和 视角方向\(\omega_o\) 均没有关系，因此在渲染方程中可以拆解到积分外面：
对于光照项 Lighting，可以用 SH基函数 的线性组合来近似表示（p 是对基函数的计数，\(\omega_i\) 是光线入射方向）：
在 PRT 中，可以交换积分和求和得到基本近似的结果，因此渲染方程转换为如下表示：
此时，积分项 \(\int_{\Omega^+}B_p(\omega_i)cos\theta_i V(\omega_i) \mathrm{d}\omega_i\) 的计算结果很凑巧就是光线传输项 Lighting transport \(f(\omega_i) = cos\theta_i V(\omega_i)\) 用 SH基函数 线性表示时的系数 \(t_p\)
（注：这里这个“凑巧”的推导过程如下：\(f(\omega_i)\) 是一个二维球面函数，其可以用 SH基函数 近似线性表示为 \(f(\omega_i)\approx\sum\limits_{p}t_p\cdot B_p(\omega_i)\)，而想要求每个基函数对应的系数 \(t_p\)，就需要求 \(f(\omega_i)\) 在这个基函数下面的投影，也就是求下面这个定积分：\(t_p=\int_{\Omega^+} f(\omega_i)·B_p(\omega_i)\mathrm{d}\omega_i\)，将 \(f(\omega_i)\) 展开就是 \(t_p=\int_{\Omega^+} B_p(\omega_i)cos\theta_i V(\omega_i) \mathrm{d}\omega_i\)，就是上面这个积分项）
最终得到渲染方程的近似表示如下：
这样，渲染方程就近似简化为了对一个点积进行求和，其中 \(l_i\) 和 \(T_i\) 可以被预计算打表，只要对于每一个着色点提前对球面函数求在每一个基函数上的投影以得到每一个基函数的系数即可
以上分析的完整过程如下：
这个简化过程的缺点在于：除了光照之外（仅限于环境光照发生的旋转），场景必须是不可变的。因为一旦场景发生了更改，\(t_p=\int_{\Omega^+} B_p(\omega_i)cos\theta_i V(\omega_i) \mathrm{d}\omega_i\) 中的这个定积分结果就会发生变化，\(t_p\) 的值就会发生变化，提前打表的结果就无效了。那么为什么类似情况的光照（仅限于环境光照发生的旋转）的 \(L(\omega_i)\) 不受影响呢？因为根据球谐函数旋转运算简单 simple rotation 的性质，当光照信息这个球面函数发生旋转时，相当于每一个基函数都发生对应旋转，而基函数旋转后可以直接表示为另一套基函数的组合，因此只需要额外打一张旋转角度-基函数系数的表然后在光照信息发生旋转时代入计算就可以解决该问题
注：对于可预计算的 \(t_p\)（对应下图的 \(T_i\)），可以把其中的 \(B_p(\omega_i)\)（对应下图的\(B_i(i)\)）视为一套以基函数作为光照的另类的光照项 Lighting，这样 \(t_p\) 就是对这套光照项的预计算，最后的渲染方程就是将这套预计算的每个值分别乘一个系数 \(l_p\) 然后求和的结果：
B）Diffuse Case 下对渲染方程的第二种拆解方式¶
对于 Diffuse Case，渲染方程的光照传输项中的 BRDF 部分和视角方向 \(\omega_o\) 无关，此时光照项 Lighting 和光线传输项 Lighting transport 可以都用 SH基函数 的线性组合来近似表示：
然后将其代入渲染方程得到下式：
虽然这个双重求和公式看起来是 \(O(n^2)\) 的复杂度，但是根据球谐函数的正交性特性， \(B_p(\omega_i)·B_q(\omega_i)\) 的值只有在 \(p\) 和 \(q\) 是相同情况下（即同一个基函数）才不是0。这意味着这个双重求和只需要求 \(p=q\) 的情况，又降低为 \(O(n)\) 复杂度了
以上分析的完整过程如下：
（3）PRT 的 Glossy Case¶
Glossy 的物体相对于 Diffuse 的物体，其渲染方程最大的一点区别在于：BRDF项 不再是一个常数，而是和光线入射方向 \(\omega_i\) 和视角方向 \(\omega_o\) 直接相关。光线入射方向 \(\omega_i\) 在之前的讨论中一直是看作一个可变量处理，因此不影响复杂度；但是视角方向 \(\omega_o\) 在 Glossy Case 中是一个新增的可变量，这就导致以 PRT 的方法近似得到的渲染方程（以第一种拆解方式为例）新增了一个维度 \(\omega_o\)：
对此的一个直观的理解就是：之前 Diffuse Case 无论摄像机视角怎么移动，着色结果都不会发生变化，因此 \(T_i\) 可以直接打表进行一个一维的预计算；而在 Glossy Case 中一旦摄像机视角 \(\omega_o\) 发生移动，\(t_p\) 作为一个关于 \(\omega_o\) 的二维函数（\(\theta\) 和 \(\phi\)）就会发生变化，如果仍然想要为 \(t_p\) 进行打表，需要为场景中的每一个可能的摄像机视角（\(\theta\) 和 \(\phi\)）进行打表，这显然是一个极其大的计算和存储量
进一步，考虑到 \(t_p(\omega_o)\) 本身也是一个二维球面函数，也可以利用 SH基函数 的线性组合来表示（注意由于 \(t_p\) 本身就是一个根据 \(p\) 变化的值，已经有了一个维度，这里再进行一次近似变换需要多出一个维度，所以前面的系数是一个二维矩阵 \(t_{pk}\)，而不是之前的一维向量）：
代入到公式中可以得到如下近似渲染方程：
其实现的计算效果如下：
这样，渲染方程就近似简化为了对一个向量和矩阵的乘积进行求和，每一个着色点需要提前对球面函数 \(L(\omega_i)\) 求在每一个基函数上的投影以得到每一个基函数的系数 \(l_p\)，然后对于这一维的每一种情况再求一次对球面函数 \(t_p(\omega_o)\) 在每一个基函数上的投影以得到每一个基函数的系数 \(t_{pk}\)
以上分析的完整过程如下：
Glossy Case 的代价相对 Diffuse Case 大幅提升，体现在
- 存储空间上，每个点都需要额外存储一个二维的 transport matrix，而不是之前的一个一维向量 \(t_p\)
- 时间复杂度上，以使用前4阶 SH基函数 的情况为例，Diffuse Case 的每个点只需要计算1次 size 为16的向量点乘，而 Glossy Case 的每个点需要将计算1次 size 为16的向量和 size 为16*16的矩阵的乘积
4）PRT 的优缺点¶
PRT 的优势：
- 对于任何光线路径，从没有反射的 LE (Light-Eye)，到简单的直接光照 L(D|G)E (Light-Diffuse/Glossy-Eye)，再到复杂的数次 bounce 的间接光照 L(D|G)·E 以及涉及到镜面反射的 LS·(D|G)·E (Light-Specular-Diffuse/Glossy-Eye)，都可以将 L 和 E 之间的部分看作是 Light Transport项，然后利用 PRT 的思想进行对应的预计算（虽然情况越复杂预计算所需的空间和时间开销越大，但是不论多长时间其都是独立于实时渲染之外的）
- 预计算带来的实时渲染效率的提高十分显著
PRT 的缺陷：
- PRT 的预计算要求场景不可动（但是视角和光源可动）
- PRT 只适用于比较低频的信息，对于高频的信息（例如材质特别 Glossy，接近于镜面反射的）实现效果不好，因为高频信息需要特别高阶的 SH基函数 才能够近似。因此，虽然根据 PRT 的优势1，它其实任何光照（包括实时全局光照）都可以实现，不仅限于实时环境光照，但是实际应用中考虑到开销问题，PRT 的主要应用场景还是在大部分都是低频信息的实时环境光照
- 大量的预计算数据，预计算需要的时间和空间的开销很大
4）补充内容：小波函数 Wavelet¶
小波函数 Wavelet 和球谐函数类似，也是一系列的基函数。区别于球谐函数是定义在球面上的，小波函数可以理解为定义在图像块上的，并且不同小波基函数的定义域还不相同
小波函数相对于球谐函数的优势在于：它支持表示全频率的信息，对于高频信息的还原效果显著优于球谐函数，因此对于高频的光照和阴影信息都可以较好地还原
但是对应的，小波函数也有其劣势：不支持快速旋转。一旦发生旋转的光照只能被认为是一个新的光照，所有计算过程需要全部重新走一遍（对应的，球谐函数由于有 simple rotation 的性质，对于光照的旋转十分支持，不需要重新计算所有数据）
小波函数的具体内容在此不做额外展开
本篇笔记主要参考了以下两篇博客，感谢 WC Yang 和 zhiwei 两位大佬的分享：
《GAMES202：高质量实时渲染》2 实时环境光照：Split Sum、PRT - 知乎 (zhihu.com)
2023年7月 ziao德