基于预计算的实时环境光照（Real-time Environment Lighting Based On Precomputation）

目录

• 概念：环境光照与全局光照（Environment Lighting & Global Illumination）
• 基于图像的光照(Image Based Lighting，IBL)
  • The Split Sum Approximation
    • 过滤环境贴图
    • 预计算BRDF积分
• 预计算辐射度传输（Precomputed Radiance Transfer，PRT）
  • 球谐（Spherical Harmonics，SH）
  • 球谐光照（Spherical Harmonic Lighting）
    • Diffuse 物体的球谐光照
    • Glossy 物体的球谐光照
  • 球谐光照预计算 Interreflection Transfer
    • Diffuse 物体的 Interreflection Transfer
  • 球谐的旋转（SH Rotation）
• 光照探针（Light Probe）
  • Unity 四面体镶嵌（Tetrahedral Tessellations）
  • UE4 间接光缓存（ILC） & 体积光照贴图（VLM）
• 光照贴图（Light Map）
  • Light Map 结合动态光照
  • Lighting Scenarios 实现昼夜天气系统
• 总结
• 参考

概念：环境光照与全局光照（Environment Lighting & Global Illumination）

---

先理清楚环境光和全局光照的概念，一个完整的光照效果往往包含：主光+间接光+环境光。

• 主光：主要的光源，例如 directional light、spot light、area light
• 间接光：任何光照射到物体后反弹出来的光，体现了物体与物体间的光照关系
• 环境光：本质属于间接光，但可以被看成是一种特殊的光源，例如 skylight、background

虽然环境光很多情况下都被当成一种特殊的光源来处理，但我仍然把它归类为间接光，这是因为：

在现实物理世界中，环境光意味着主光照射到远处（即场景外）的物体（包括天空）后反弹出来的光；在实际 Rendering 上，也可以实时地将场景内物体的间接光照信息包含进环境光。

而 全局光照（Global Illumination,GI），严格意义上指主光+完整的间接光下的光照效果，实际上只要有体现间接光效果（即使是部分效果）的光照都可称为 GI。

因此本文将标题改成 Environment Lighting 以强调本文的算法是基于环境光这个角度出发的，但是所属范畴仍然可算 GI 方法。

此外，GI 总体上分为两派：

• 一派是通过一套统一标准的渲染流程把任何物体的 GI 计算出来，往往计算量极大但效果更加物理更加真实，典型的例子便是离线光线追踪
• 另一派是把 GI 的效果看成一个个部分来组成，这样我们可以选择其中一些 GI 效果组合使用，适应不同的性能（往往计算量相对低些，尤其是实时渲染）的同时也能带来能接受的 GI 效果（虽然往往不是严谨的物理正确），典型的例子就是这篇博客所介绍的方法和另一些完全动态的 GI 方法（如 RSM、SSAO、SSRT）

基于图像的光照(Image Based Lighting，IBL)

---

基于图像的光照（IBL），简单的说就是一类通过 环境贴图（Environment Map） 来保存某个物体的环境光信息，从而实现该物体的基于物理的物体渲染方法。

• IBL 可用于 diffuse/glossy/specular 物体（基本上包含大部分物体了）渲染：这取决于环境贴图，环境贴图分辨率越大，那么所能表示高频信息就越多，从而越适合 specular 物体的渲染；而分辨率越低，所表示的信息更多是低频信息，则 diffuse 物体的渲染也足以满足，而且还能节省一定存储空间。

• IBL 可以实现成动态环境光照：实时渲染出动态的环境贴图。有一定开销，一般只用于少量的specular物体。

• IBL 可以实现成静态环境光照：预渲染环境贴图。需要环境光是静态的。

用 IBL 实现静态环境光照就有点类似 Light Map的做法了，区别在于 IBL 存的是静态环境光信息，Light Map 存的是接受静态环境光信息后的着色结果。

为了表示来自四面八方的环境光信息，IBL 使用 Spherical Map 或者 Cube Map 的方式来存储：

The Split Sum Approximation

---

IBL 中最常见的算法便是基于 The Split Sum Approximation 的算法。

我们知道，一般的渲染方程如下：

(L_{r}left(x, omega_{r} ight)=int_{Omega^{+}} L_{i}left(x, omega_{i} ight) f_{r}left(x, omega_{i}, omega_{r} ight)left(n cdot omega_{i} ight) mathrm{d} omega_{i})

在实际渲染的时候，我们当然可以使用蒙特卡洛方法实现该渲染方程，然而这样的开销是巨大的（每个shading point都要做多重采样，而且结果很容易是noisy的）。

为了进一步加快式子，这里有一个很经典的近似公式：

(int_{Omega^+} f(x) g(x) mathrm{d} x approx frac{int_{Omega_{G}} f(x) mathrm{d} x}{int_{Omega_{G}} mathrm{~d} x} cdot int_{Omega^+} g(x) mathrm{d} x)

积分域 (Omega_G) ：在原本 (Omega^+) 的积分域范围内剔除掉 (f(x) = 0) 的地方而剩余的范围。

想让这个公式近似效果比较精确，那么需要满足以下一种或两种条件：

• 积分域 (Omega_G) 比较小
• (g(x)) 比较光滑，即变化不是很大

而我们的观察是：

• 如果 BRDF 是 glossy/specular 的，那么它的 lobe 往往是花瓣状，即只有很小的积分域才能接受环境光。
• 如果 BRDF 是 diffuse 的，那么它的 lobe 往往是均匀的半球状，即无论哪个方向的环境光打进来， (f_r) 函数的输出几乎没多少变化（甚至是个常数）。

于是基于以上理论，Split Sum 方法对渲染方程改造成这样的近似公式来获得渲染的加速：

(L_{r}left(x, omega_{r} ight) approx frac{int_{Omega_{f_{r}}} L_{i}left(x, omega_{i} ight) mathrm{d} omega_{i}}{int_{Omega_{f_{r}}} mathrm{~d} omega_{i}} cdot int_{Omega^{+}} f_{r}left(x, omega_{i}, omega_{r} ight) cos heta_{i} mathrm{~d} omega_{i})

对渲染方程拆分成两个部分（环境光积分、BRDF积分）后就可以通过预计算的方式（后面两节会介绍如何预计算）分别减少这些积分的运行时开销，总结这种方法的好处是：

• Split Sum 方法和原始蒙特卡洛方法的图像效果几乎一模一样
• 由于不用对环境贴图进行多重采样，性能开销大大减低了

过滤环境贴图

环境光积分：(frac{int_{Omega_{f_{r}}} L_{i}left(x, omega_{i} ight) mathrm{d} omega_{i}}{int_{Omega_{f_{r}}} mathrm{~d} omega_{i}})

因为我们已经拥有一张环境贴图（无论是实时的还是预渲染的）来存储环境光信息了，为了计算环境光部分的积分，需要在 (Omega_{fr}) 范围内做多次光线采样。但是，可以有一个几乎等价但避免运行时多重采样的方式：

预先对纹理进行滤波操作（模糊），我们只需要对滤波后的环境贴图采样一次光线方向就能得到积分。

滤波后的环境贴图实际上称为 Irradiance Environment Map ：原来环境贴图中的单位是 radiance (L(mathrm{x}, omega))，贴图积分后的单位则变成了 Irradiance (E(x))。

当然，BRDF Lobe 的形状越尖锐，即环境光积分范围越小，就需要使用模糊程度更低的环境贴图；反之 BRDF Lobe 的形状越粗壮，即环境光积分范围越大，就需要使用模糊程度更高的环境贴图。

我们可以使用 MIPMAP技术 来生成不同Level的环境贴图，通过三线性插值（Trilinear Interportion）的方式来得出任何模糊程度且任何2D位置的环境光滤波结果。

预计算BRDF积分

回顾 Microfacet BRDF 的组成：

• Fresnel项（举例 Schlick's approximation）

(F=F_0 +(1-F_0)(1-(h cdot omega_{r}))^{5} = F_0 +(1-F_0)(1-(cos heta_{vh}))^{5})

• NDF项（举例 Beckmann NDF）

(D(h)=frac{1}{pi alpha^{2}({n} cdot {h})^{4}} cdot exp {left(frac{({n} cdot {h})^{2}-1}{a^{2}({n} cdot {h})^{2}} ight)} = frac{1}{pi alpha^{2}cos^4 heta_h } cdot exp {left(frac{cos^2 heta_h -1}{a^{2}cos^2 heta_h } ight)})

可以预想到，该BRDF 的积分结果依赖于三个参数：

• (F_0)（Fresnel项系数）
• (alpha) （粗糙度 roughness）
• ( heta_v) （反射方向与法线的夹角，实际上决定了 ( heta_{vh})、( heta_{h})）

我们可以把 (F_0) 项拆出来，让BRDF 积分拆成两个积分，但是这些积分都减少了一个依赖的参数：

(egin{align} int_{Omega^{+}} f_{r}left(p, omega_{i}, omega_{o} ight) cos heta_{i} mathrm{~d} omega_{i}approx & F_{0} int_{Omega^{+}} frac{f_{r}}{F}left(1-left(1-cos heta_{vh} ight)^{5} ight) cos heta_{i}mathrm{~d} omega_{i} \ & + int_{Omega^{+}} frac{f_{r}}{F}left(1-cos heta_{vh} ight)^{5} cos heta_{i} mathrm{~d}omega_{i} end{align})

这样，对于相同的材质（相同的BRDF），我们就可以针对剩余两个依赖的参数 (alpha) 、( heta_v) 建立一张 二维查询表：

预计算辐射度传输（Precomputed Radiance Transfer，PRT）

---

预计算辐射度传输（Precomputed Radiance Transfer，PRT） 是一类通过预计算辐射度传输的基于物理的物体渲染方法。所谓辐射度传输，可以理解成物体自身的阴影/AO（AmbientOcclusion）和表面互反射等有关于光路传输的信息。因此 PRT 往往适用于动态光照下的静态物体/静态材质。

• PRT 可以实现动态环境光照：仅预计算 transfer 部分
• PRT 可以实现静态环境光照：不仅预计算 transfer，还顺便预计算环境光部分（可以称为 radiance）；这样就能以更低性能开销实现运行时基于物理的渲染，只是光照动态性会有所限制（但不一定完全静态）

顺便一提，PRT 方法在2002年 Siggraph 会议被 Peter-Pike Sloan 首先提出，并从此掀起了一波 PRT 方法的研究浪潮。当时 Peter-Pike Sloan 论文所实现的 PRT 算法便是基于球谐（SH）的，实际上 PRT 不仅可以通过 SH 去表示低频环境光场景，还可以有其它方法（如Wavelet）。不过本文仍将主要介绍球谐光照（SH Lighting）方法。

球谐（Spherical Harmonics，SH）

---

对于任意函数 (f(x)) ，不管是连续的还是不连续的，我们可以展开成一系列基函数（每项都带某个系数）的线性组合：

(f(x)=sum_{i} c_{i} cdot B_{i}(x))

例如，多项式展开可以看成是一系列基函数（多项式）的线性组合：(f(x)=c_{0} +c_1cdot x^1 + c_2 cdot x^2 + ...)

再例如，傅立叶变换也可以将 (f(x)) 表示成另一系列基函数（各种频率的正弦谐波）的线性组合：

而 球谐（Spherical Harmonics，SH） 便是定义在球面上的一系列2D基函数，它与2D傅里叶序列有点相似，但非常适合球面函数 (f(omega))（即参数为单位球面向量）。

SH 是分阶数的：在第0阶（l=0）有1个基函数（m=0）；第1阶（l=1）有3个基函数（m=-1,0,1）...第n阶（l=n）有2n+1个基函数（m=-l,...,0,...,l）。实际上，阶数越低的基函数所代表的信息就越是低频。如果要完全复原一个任意函数，我们需要无穷阶的 SH；而如果我们只要复原一个任意函数的近似（换句话说只重建出该函数的低频信息），那么我们完全可以只需要前几阶的 SH（包含 l=0,l=1,...,l=n 每一阶的所有基函数）。

SH 的实数形式的基函数如下：

SH 基函数的形式比较复杂，只建议看一眼就跳过。

(y_{l}^{m}( heta, varphi)=left{egin{array}{rl}sqrt{2} operatorname{Re}left(Y_{l}^{m} ight) & m>0 \ sqrt{2} operatorname{Im}left(Y_{l}^{m} ight) & m<0 \ Y_{l}^{0} & m=0end{array}=left{egin{array}{cc}sqrt{2} K_{l}^{m} cos m varphi P_{l}^{m}(cos heta) & m>0 \ sqrt{2} K_{l}^{m} sin |m| varphi P_{l}^{|m|}(cos heta) & m<0 \ K_{l}^{0} P_{l}^{0}(cos heta) & m=0end{array} ight. ight.)

其中，

(K_{l}^{m}=sqrt{frac{(2 l+1)(l-|m|) !}{4 pi(l+|m|) !}})

(p^m_l) 为勒让德多项式（Legendre polynomials）：

• $P_0^0= 1 $
• $P_{m}^{m} =(1-2 m) P_{m-1}^{m-1} $
• $P_{m+1}^{m} =(2 m+1) z P_{m}^{m} $
• $ P_{l}^{m} =frac{(2 l-1) z P_{l-1}^{m}-(l+m-1) P_{l-2}^{m}}{l-m} $

z为球面向量对应球面坐标的z值，该多项式不依赖x值y值。

以下是前3阶的SH实数形式的基函数：

(l=0)：

• (Y_{00}=s=Y_{0}^{0}=frac{1}{2} sqrt{frac{1}{pi}})

(l=1)：

• (Y_{1,-1}=p_{y}=i sqrt{frac{1}{2}}left(Y_{1}^{-1}+Y_{1}^{1} ight)=sqrt{frac{3}{4 pi}} cdot frac{y}{r})
• (Y_{1,0}=p_{z}=Y_{1}^{0}=sqrt{frac{3}{4 pi}} cdot frac{z}{r})
• (Y_{1,1}=p_{x}=sqrt{frac{1}{2}}left(Y_{1}^{-1}-Y_{1}^{1} ight)=sqrt{frac{3}{4 pi}} cdot frac{x}{r})

(l=2)：

• (Y_{2,-2}=d_{x y}=i sqrt{frac{1}{2}}left(Y_{2}^{-2}-Y_{2}^{2} ight)=frac{1}{2} sqrt{frac{15}{pi}} cdot frac{x y}{r^{2}})
• (Y_{2,-1}=d_{y z}=i sqrt{frac{1}{2}}left(Y_{2}^{-1}+Y_{2}^{1} ight)=frac{1}{2} sqrt{frac{15}{pi}} cdot frac{y z}{r^{2}})
• (Y_{2,0}=d_{z^{2}}=Y_{2}^{0}=frac{1}{4} sqrt{frac{5}{pi}} cdot frac{-x^{2}-y^{2}+2 z^{2}}{r^{2}})
• (Y_{2,1}=d_{s z}=sqrt{frac{1}{2}}left(Y_{2}^{-1}-Y_{2}^{1} ight)=frac{1}{2} sqrt{frac{15}{pi}} cdot frac{z x}{r^{2}})
• (Y_{2,2}=d_{x^{2}-y^{2}}=sqrt{frac{1}{2}}left(Y_{2}^{-2}+Y_{2}^{2} ight)=frac{1}{4} sqrt{frac{15}{pi}} cdot frac{x^{2}-y^{2}}{r^{2}})

SH 与一般的基函数相比，具有以下性质：

• 基函数之间具有正交性（orthonormal）

(int_{Omega} B_{i}(omega) cdot B_{j}(omega) mathrm{d} omega=mathbf{1} quad(i=j))

(int_{Omega} B_{i}(omega) cdot B_{j}(omega) mathrm{d} omega=mathbf{0} quad(i eq j))

• 通过 投影（Projection） 可以很方便得到 SH 系数（SH coefficients）

(f_i = int_{Omega}f(omega)cdot B_i(omega)domega)

• 通过系数向量组与基函数组的点积（前n阶的基函数/系数共有 n*n 个）可以很方便重建球面函数

(f(omega)approx sum^{n^2}_{i=1} f_{i} B_{i}(mathbf{omega}))

• product projection：(c(omega)=a(omega)b(omega)) ，已知 (a) 的形式而 (b) 未知，有

(vec{c} = Mcdot vec b)

其中，(vec{c}) 为 (c(omega)) 的SH系数向量组，(vec{b}) 为 (b(omega)) 的SH系数向量组，(M) 为 (n^2) X (n^2) 的矩阵，其元素为

(M_{ij}=int_{Omega}a(omega)B_i(omega)B_j(omega)domega)

这样我们就可以预计算好 (M) ，等到知道 (vec{b}) 后，乘起来就能得到 (vec{c})

推导：

(令 M_{i}(omega)=a(omega)B_i(omega))

(egin{align} c_i &= int_{Omega} a(omega)b(omega)B_i(omega) domega \ &= int_{Omega} M_i(omega)b(omega)domega \ &= int_{Omega} sum_{j=1}^{n^2} (M_{ij} B_j(omega)) b(omega) d omega \ &= sum_{i=1}^{n^2} M_{ij} int_{Omega} B_j(omega)b(omega) d omega \ &= sum_{j=1}^{n^2} M_{ij}cdot b_j end{align})

• 支持插值，对 SH 系数的插值相当于对重建的函数值的插值。

• 旋转不变性（rotational invariance），对函数 (f) 的旋转 (R_{SH}) 等价于对 (f(omega)) 的自变量的旋转 (R_{3D})：

(R_{SH}(f(omega))=f(R_{3D}(omega)))

球谐光照（Spherical Harmonic Lighting）

---

纹理本质上也是一个函数（信号），输入二维坐标，输出对应纹素的RGBA值。传统的环境光信息往往使用环境贴图（Environment Map）表示，这往往需要一个庞大的二维数组存储各个纹素的值，十分耗费空间，而且采样和纹理I/O也有一定开销。

球谐光照（SH Lighting） 的思路：将渲染方程分为两个球面函数（即 lighting function 和 transfer function），这些球面函数将使用 SH 方法来表示：

• 对于环境光（lighting）部分，只需预先存储若干个 SH 系数而不必存储一整张环境贴图。
• 对于传输函数（transfer function） 部分，只需预先存储若干个 SH 系数（diffuse情况下）或者矩阵（glossy情况下），预计算好传输率可以在实时渲染中以极低代价实现自阴影、互反射的效果。

不过 SH Lighting 一般采用3阶SH，因此 SH 所能表示的 lighting 和 transfer 信息是低频的，只适用于 diffuse 和 glossy 的物体而不适用于 specular 的物体。

SH lighting 效果（分别为 diffuse物体的无阴影情况、diffuse物体的阴影&互反射情况、glossy物体的无阴影情况、glossy物体的阴影&互反射情况）：

Diffuse 物体的球谐光照

Diffuse 物体的渲染方程：

(L(r)=int_{Omega^+} L_e(mathbf{omega})cdot ho cdot V(mathbf{omega}) max (0, ncdot mathbf{omega}) mathrm{d} mathbf{omega})

(L_e) 代表环境光； ( ho) 为 BRDF 项；(V) 为 Visibility，表示不被遮挡的程度，往往表现为自遮挡产生的阴影现象。

• 由于物体是 diffuse 的，因此它的 BRDF 将是一个常数 ( ho) （无论从哪个方向观察都得到相同的BRDF值），也因此 diffuse 物体的 (L) 将是一个常量值，而不受参数 (r) 影响
• 对于 lighting 部分 (L(mathbf{omega}))，原本需要对 Environment Map 进行查询，而现在可以换成使用 SH 函数去表示：(L_e(mathbf{omega}) approx sum l_{i} B_{i}(mathbf{omega}))
• 对于 transfer function 部分 (T(omega)=V(mathbf{omega}) max (0, ncdot mathbf{omega})) ， 也可以换成使用 SH 函数去表示：(T(omega) approx sum T_j B_j(omega))

代入渲染方程整理后得：

(egin{align} L(r) &approx ho int_{Omega^+} sum_{i=1}^{n^2}l_{i}B_i(omega) sum_{j=1}^{n^2}T_jmathrm{B}_{j}(mathbf{omega}) mathrm{d} mathbf{omega} \ &= hosum_{i=1}^{n^2} sum_{j=1}^{n^2} l_{i}T_jint_{Omega^+} B_i(omega) mathrm{B}_{j}(mathbf{omega})domega \ &= hosum_{i=1}^{n^2} l_{i}T_iint_{Omega^+} B_i(omega) mathrm{B}_{i}(mathbf{omega})domega \ &= hosum_{i=1}^{n^2} l_{i}T_i end{align})

这样，我们只需要预计算出 (l_i) 和 (T_i)，就能让运行时的diffuse物体渲染速度大大提升。

预计算过程：

1. 对整个环境光信息（环境贴图），预计算 lighting 的 SH 系数向量组 (vec{l})，其中一个元素为：(l_{i}=int_{Omega^+} L_e(mathbf{omega}) B_{i}(mathbf{omega}) mathrm{d} mathbf{omega})

2. 对每个顶点，预计算 light transfer 的 SH 系数向量组 (vec{T})，其中一个元素为：(T_i=int_{Omega^{+}}V(mathbf{omega}) max (0, ncdot mathbf{omega}) mathrm{B}_{i}(mathbf{omega}) mathrm{d} mathbf{omega})

此时，我们可以理解成系数 (l_i) 代表了环境光照（lighting）的信息，而系数 (T_i) 代表了在某个顶点上光路传输（light transfer）的信息；整个模型预计算完成之后将会得到一个 lighting 系数向量和模型顶点数量对应的 transfer 系数向量。

运行时渲染过程：

1. 在 vertex shding 阶段计算顶点的 SH 颜色 (L(r) approx ho sum_{i=1}^{n^2} l_{i} T_{i} = ho( vec{l}cdot vec{T}))

2. 在 pixel/fragment shading 阶段得到插值后的 SH 颜色即为该像素的颜色

Glossy 物体的球谐光照

Glossy 物体的渲染方程：

(L(r)=int_{Omega^+} L_e(mathbf{omega})cdot ho(r,mathbf{omega}) cdot V(mathbf{omega}) cdot max (0, ncdot mathbf{omega}) mathrm{d} mathbf{omega})

对于 glossy 物体，不同视角观察物体表面同一点会有不同的光照。因此 glossy 的 BRDF 将是四维的函数（参数不仅包含(mathbf{omega})，还包含(r)），这次将 BRDF 算入 transfer，则 transfer function 将额外增加一个二维参数 (r)

• 对于环境光（lighting）即 (L(mathbf{omega}))，原本需要对 Environment Map 进行查询，而现在可以换成使用 SH 函数去表示：(L_e(mathbf{omega}) approx sum l_{i} B_{i}(mathbf{omega}))
• 对于传输函数（transfer function） (T(r,omega)= ho(r,mathbf{omega}) V(mathbf{omega}) max (0, ncdot mathbf{omega})) ，换成使用 SH 函数去表示：(T(r,omega) approx sum_{i=1}^{n^2} T_{i}(r)B_i(omega)approx sum_{i=1}^{n^2} sum_{j=1}^{n^2} T_{ij} B_j(r)B_i(omega))

由于额外多了一个二维的参数，transfer 系数将是一个矩阵而非之前 diffuse 情况下的系数向量

代入渲染方程整理后得：

(egin{align} L(r) & approx int_{Omega^+} sum_{i=1}^{n^2} l_{i} B_{i}sum_{j=1}^{n^2}sum_{k=1}^{n^2} T_{jk} B_j(omega)B_k(r)mathrm{d} mathbf{omega} \ &= sum_{i=1}^{n^2}sum_{j=1}^{n^2}sum_{k=1}^{n^2} l_iT_{jk }B_k(r) int_{Omega^+} B_{i}(omega)B_j(omega)mathrm{d} mathbf{omega} \ &= sum_{i=1}^{n^2}sum_{k=1}^{n^2}l_iT_{ik }B_k(r) end{align})

预计算过程：

1. 对整个环境光信息，预计算 lighting 的 SH 系数向量组 (vec{l})，其中一个元素为：(l_{i}=int_{Omega^+} L(mathbf{omega}) B_{i}(mathbf{omega}) mathrm{d} mathbf{omega})
2. 对每个顶点，预计算 light transfer 矩阵 (T)，其中一个元素为：(T_{ij}=int_{Omega^{+}}mathrm{B}_{j}(mathbf{omega}) T_i(omega) mathrm{d} mathbf{omega}) 、(T_i(r) = int_{Omega^+}T(r,omega)B_i(omega)domega)

由于 Glossy 物体需要每个顶点存储一个矩阵 T，这使得空间开销略大，因此不太实用。

运行时渲染过程：

1. 在 vertex shding 阶段计算顶点的 SH 颜色 (L(r) approx sum_{i=1}^{n^2}sum_{j=1}^{n^2} l_{i} T_{ij} B_j(r) = Tvec{l}cdot vec{B(r)})
2. 在 pixel/fragment shading 阶段得到插值后的 SH 颜色即为该像素的颜色

原论文实际上还有个 BRDF 卷积步骤。但是本身用SH表示环境光信息已经足够低频（模糊）了，对这个卷积步骤，个人先存个疑。

球谐光照预计算 Interreflection Transfer

而前面的渲染方程中：

(L_{x}(r)=L^0_{x}(r)=int_{Omega^+} L_{e}(mathbf{omega})cdot ho(r,mathbf{omega}) cdot V(mathbf{omega}) cdot max (0, ncdot mathbf{omega}) mathrm{d} mathbf{omega})

只考虑几何传播关系的因素只有 (V(omega)cdot max(0,ncdotomega))，即只可实现阴影或者无阴影（令 (V=1)）的效果，并不能实现互反射（interreflection）的效果。

所谓 interreflection，实际上就是某个点朝外围环境看时，(V=0) 意味着有几何部分遮挡住了直接的环境光，但实际上遮挡还意味着这些几何部分上的出射光会反射给该点，这就形成了该点间接的环境光照。

增加 1 次互反射的渲染方程（迭代）：

(L^n_x(r) := L^0_x(r)+int_{Omega^+} L^{n-1}_{x'}(mathbf{omega})cdot ho(r,mathbf{omega}) cdot (1-V(mathbf{omega})) cdot max (0, ncdot mathbf{omega}) mathrm{d} mathbf{omega})

(L^n_x) 意味着在某个点 (x) 的 radiance，而 (x') 则是从 (x) 朝 (omega) 方向进行 raycast 而得到相交的其它点，而 (n) 则代表了迭代了多少次。

Diffuse 物体的 Interreflection Transfer

(L^n_x := L^0_x + hoint_{Omega^+} L^{n-1}_{x'} cdot (1-V(mathbf{omega})) cdot max (0, ncdot mathbf{omega}) mathrm{d} mathbf{omega})

将上述式子左右两边同时分离出 (vec{l}) （即对预计算出来的 lighting SH 系数不做改变 ），则意味着我们实现 Interreflection 的效果只需要在所有顶点 (x) 的 transfer 系数向量做以下迭代：

((T^{n}_x)_i := (T^0_x)_i + ho int_{Omega^{+}}(T^{n-1}_{x'})_i cdot (1-V(mathbf{omega})) cdot max (0, ncdot mathbf{omega}) mathrm{d} mathbf{omega})

球谐的旋转（SH Rotation）

PRT 的一个问题是如果 lighting 部分是预计算的，那就只适用于静态环境光下的静态物体渲染；环境光或者物体只要有变化，PRT 就不得不进行重新预计算；但得益于 SH 的旋转不变性，我们至少可以让 SH Lighting 适用于动态旋转的情形而不必重新预计算。

环境光旋转往往用来实现场景的昼夜轮转效果，物体旋转就比较常见了。顺便注：环境光旋转和物体旋转在 PRT 渲染中是等价的（只是说看相对于哪个东西来看待旋转而已）

现在给定旋转 (R) ，然后有SH投影函数 (P)（输入一个球面向量，输出第 (l) 层band的 SH 系数向量组）。假设我们有 (2l+1) 个任意的球面向量，我们需要想办法求出能等价于旋转 (n_i) 的矩阵 (M) 来旋转 SH投影函数 (P) ：

(M Pleft(n_{i} ight)=Pleft(Rleft(n_{i} ight) ight), i in[-l, l])

整理得：

(Mleft[Pleft(n_{-l} ight), ldots, Pleft(n_{l} ight) ight]=left[Pleft(Rleft(n_{-l} ight) ight), ldots, Pleft(Rleft(n_{l} ight) ight) ight])

记 (A = [P_{(n−l)}, ..., P (n_{l})]) ，如果矩阵 (A) 是可逆的，则：

(M=left[Pleft(Rleft(n_{-l} ight) ight), ldots, Pleft(Rleft(n_{l} ight) ight) ight] A^{-1})

这样，无论 (R) 怎么变化，我们都可以即时根据 (R) 求出 (M) ，然后把最新的 (M) 应用到所有预计算好的 (P(omega)) 上（即预计算好的SH系数上），而不必重新预计算 (P(R(omega))) 。

SH 的快速旋转实现：

对于第 (l) 层 band 需要做如下处理：

1. 选取共 (2l + 1) 个 单位向量 (n_i)，这些向量 (n_i) 投影在该层 band 上的基函数就能得到系数 (P (n_i))。共 (2l + 1) 个 (2l+1) 维向量 (P (n_i)) 构成了矩阵 (A)，并求出 (A^{−1})

如何选取单位向量：要保证投影后构成的 (A) 矩阵可逆。

2. 给定旋转 (R) ，对所有 (n_i) 依次做旋转 (R(n_i)) ，这些旋转后的单位向量同样投影在该层 band 上的基函数就能得到系数 (P(R(n_i)))。 共 (2l + 1) 个 (2l + 1) 维向量 (P(R(n_i))) 构成了矩阵 (S)

3. 求出该层 band 上球谐函数的旋转矩阵 (M = SA^{−1})

4. 用 (M) 乘以该层 band 上的 SH 系数向量组就可以得到旋转后的 SH 系数向量组。

(l=0) 时只有1个系数，故不需要处理；(l=1) 时需要处理3个系数，其中的 (A、S、M) 将会是 3X3矩阵；(l=2) 时需要处理5个系数，其中的 (A、S、M) 将会是 5X5矩阵......

最后，将每一层 band 的结果重新拼接起来即可得到完整的旋转后的 SH 系数结果。

// 伪代码：三阶SH旋转
SHCoeffientsAfterRotation(Array coeffients,Rotation rotation){
    Array res;
    // 处理 l = 0
    res[0] = coffients[0];
    // 处理 l = 1
    // step 1
    n1 = [1,0,0],n2 = [0,1,0],n3 = [0,0,1];
    A = [ // SHProject(n,l,m)：向量n投影到SH中l层第m个基函数，输出对应系数
    [SHProject(n1,1,-1),SHProject(n2,1,-1),SHProject(n3,1,-1)],
    [SHProject(n1,1,0),SHProject(n2,1,0),SHProject(n3,1,0)],
    [SHProject(n1,1,1),SHProject(n2,1,1),SHProject(n3,1,1)]
    ];
    A_inv = InvMatrix(A);
    // step 2
    rn1 = rotation*n1,rn2 = rotation*n2,rn3 = rotation*n3;
    S = [
        [SHProject(rn1,1,-1),SHProject(rn2,1,-1),SHProject(rn3,1,-1)],
        [SHProject(rn1,1,0),SHProject(rn2,1,0),SHProject(rn3,1,0)],
        [SHProject(rn1,1,1),SHProject(rn2,1,1),SHProject(rn3,1,1)]
    ];
    // step 3
    M = S*A_inv;
    // step 4
    coeff_l1 = M*[coeffients[1],coeffients[2],coeffients[3]];
    res[1] = coeff_l1[0];
    res[2] = coeff_l1[1];
    res[3] = coeff_l1[2];
    
    // 处理 l = 2    
    // step 1
	k = 1/sqrt(2);
    n1 = [1,0,0],n2 = [0,0,1],n3 = [k,k,0],n4 = [k,0,k],n5 = [0,k,k];
    A = [
[SHProject(n1,2,-2),SHProject(n2,2,-2),SHProject(n3,2,-2),SHProject(n4,2,-2),SHProject(n5,2,-2)],
[SHProject(n1,2,-1),SHProject(n2,2,-1),SHProject(n3,2,-1),SHProject(n4,2,-1),SHProject(n5,2,-1)],
[SHProject(n1,2,0),SHProject(n2,2,0),SHProject(n3,2,0),SHProject(n4,2,0),SHProject(n5,2,0)],
[SHProject(n1,2,1),SHProject(n2,2,1),SHProject(n3,2,1),SHProject(n4,2,1),SHProject(n5,2,1)],
[SHProject(n1,2,2),SHProject(n2,2,2),SHProject(n3,2,2),SHProject(n4,2,2),SHProject(n5,2,2)]
    ];
    A_inv = InvMatrix(A);
    // step 2
    rn1 = rotation*n1,rn2 = rotation*n2,rn3 = rotation*n3,rn4 = rotation*n4,rn5 = rotation*n5;
    S = [
[SHProject(rn1,2,-2),SHProject(rn2,2,-2),SHProject(rn3,2,-2),SHProject(rn4,2,-2),SHProject(rn5,2,-2)],
[SHProject(rn1,2,-1),SHProject(rn2,2,-1),SHProject(rn3,2,-1),SHProject(rn4,2,-1),SHProject(rn5,2,-1)],
[SHProject(rn1,2,0),SHProject(rn2,2,0),SHProject(rn3,2,0),SHProject(rn4,2,0),SHProject(rn5,2,0)],
[SHProject(rn1,2,1),SHProject(rn2,2,1),SHProject(rn3,2,1),SHProject(rn4,2,1),SHProject(rn5,2,1)],
[SHProject(rn1,2,2),SHProject(rn2,2,2),SHProject(rn3,2,2),SHProject(rn4,2,2),SHProject(rn5,2,2)]
    ];
    // step 3
    M = S*A_inv;
    // step 4
    coeff_l2 = M*[coeffients[4],coeffients[5],coeffients[6],coeffients[7],coeffients[8]];
    res[4] = coeff_l2[0];
    res[5] = coeff_l2[1];
    res[6] = coeff_l2[2];
    res[7] = coeff_l2[3];
    res[8] = coeff_l2[4];
    
    return res;
}

这样我们对 lighting 的 SH 系数做 Rotation，就可以实现支持环境光旋转或物体旋转的实时 SH Lighting。

但是注意，不应当对 transfer 的 SH 系数做 Rotation，因为一个物体的 transfer function 只考虑自身的几何遮蔽（自身与自身的遮蔽关系），而非像物体与环境光那样存在相对关系。

SH Rotation 效果图：

光照探针（Light Probe）

---

光照探针（Light Probe） 是一类场景 GI 方案。它在场景里设置若干个Probe点（可以理解成向四面八方探测光照的点），为每个 Probe 预计算出环境光信息后，在运行时通过物体周围的 Probe 信息来插值来得到物体此时受到的环境光。

实际上，PRT也算是一种基于Probe的GI方案：PRT物体上的每个顶点都是一个Probe，它所探测的环境光照信息将预计算为 SH 系数，然后 shading point 可以根据三角形的顶点对 SH 系数做重心插值从而得到该点的 SH 系数（即代表了环境光）。只是 PRT 的每个三角形顶点都是一个Probe，这些Probe综合起来完成了单个物体的 GI 效果；而一般的 Light Probe 方案是在场景中每隔一段空间放置一个Probe，这些Probe综合起来提供了场景中所有物体的 GI 效果。

不过 Light Probe 经常会遇到错误的 Bleeding 问题：受到几何上不应该存在光照关系的 Probe 影响，常见于墙壁遮挡的内外侧。

例如，下面的一个 Probe 生成在室内：

结果由于 Light Probe 方法没有考虑遮挡，直接且错误地进行了对相邻的 Probe 插值，从而室外本该明亮的一侧变暗，室内本该黑暗的一侧发生漏光：

因此一个好的 Probe 放置方案可以尽量减少这种 Bleeding 问题。

Unity 四面体镶嵌（Tetrahedral Tessellations）

Unity 允许在场景中放置自定义位置的 Probe ，而且这些 Probe 将相连成一个个四面体。

为了增加这种GI方案的真实感和避免过多的Probe带来的存储开销，还应当把 Probe 在光照发生明显变化的地方（如明暗交接处）放密集些，而在光照不怎么变化的地方可以稀疏地放置。

预计算过程：

1. 每个 Probe 从自身往四面八方看到的环境光信息烘焙成3阶SH系数（SH Lighting方法）并记录之

Unity 曾经提供了一种动态实时GI方案，只不过现在被废弃了：它的 Probe 是运行时采样并实时计算出最新的2阶 SH 系数。

2. 将这些 Probe 进行三角化（将空间切分成四面体），切分的算法是使用 Delaunay Triangulation 完美三角剖分。

运行时渲染过程：

1. 在渲染物体的 shading point 时，通过四面体插值的方法取周围4个 Probe 的结果的加权平均作为该点受到的 lighting

四面体插值：假设 (a,b,c,d) 分别为四个 Probe 的权重，则

(left[egin{array}{l}a \ b \ cend{array} ight]=left[egin{array}{lll}vec{P}_{0}-vec{P}_{3} & vec{P}_{1}-vec{P}_{3} & overrightarrow{P_{2}}-vec{P}_{3}end{array} ight]^{-1}left[vec{P}-overrightarrow{P_{3}} ight])

(d=1-a-b-c)

渲染效果：

Unity Tetrahedral Tessellations 的弊端：

• 运行时，查找所处四面体的CPU运算量较高：四面体的分布并不是均匀的，GPU很难做到查找，所以Unity需要通过CPU把每个使用GI的物体所处的四面体找到并将查找结果传递给GPU，这就带来了性能瓶颈。

UE4 间接光缓存（ILC） & 体积光照贴图（VLM）

间接光缓存（Indirect Light Cache）：与 Unity Light Probe 方案基本一样，Probe 也是采用烘焙 SH 系数和插值的方案，只是放置 Probe 的方式有所不同。ILC 基本上是在静态物体表面法线向上均匀放置 Probe ，这是基于假设受环境光影响最大的地方都是在靠近静态物体的地方。

当然实际上 ILC 还包含两种 Probe 放置方式：通过手动放置 Lightmass Importance Volume（重要光照范围） 来限制烘焙光照采样范围，或 Lightmass Character Indirect Detail Volume，来增加一段均匀放置 Probe 的区域。

ILC 通过 CPU 寻找物体周围的 Probe：所有的 Probe 的数据将保存到八叉树中以方便物体找到周围的 Probe。

体积光照贴图（Volumetric Light Map，VLM）：还是采用烘焙 SH 系数和插值的方案，也还是放置 Probe 的方式有所不同。VLM 使用网格（Grid）采样点保存 Probe 数据，在静态物体表面附近，会对网格进一步细分。

VLM 将通过 GPU 来寻找物体周围的 Probe ：所有 Probe 的数据将烘焙至贴图中，不同细分层度的网格将使用 level 不同的贴图，这样可以方便地在 GPU 中进行逐像素的三线性插值。

VLM 与 ILC 比较：

• VLM 比 ILC 的 Probe 要多更多，从而 VLM 的 GI 效果更佳，但要存储的数据更多
• VLM 的存储结构决定了可以通过 GPU 算法来寻找 Probe，这比传统的 CPU 算法性能更加客观（即是 VLM Probe 数量要多得多）
• VLM 更适合将 Probe 应用到体积雾效果中，这是因为 ILC 往往不在几乎没什么物体的空间放置 Probe

UE 4.18 以后使用 VLM 代替了 ILC 方法作为 UE4 默认的 Light Probe 方案。

UE4 ILC & VLM 效果图：

光照贴图（Light Map）

---

有关 Light Map 的初步介绍可见 实时渲染基础（4）纹理-光照贴图

Light Map 结合动态光照

---

在 UE4 和 Unity 中，静态物体的全局光照效果都采用了烘焙 Light Map 的方法。

Light Map 方法实际上就是把环境光被分成了静态部分和动态部分，而 Light Map 记录的正是受静态部分环境光后着色结果，因此参与烘焙的物体（位置、形状、材质）和光照都应当是静态的，这样的局限性很大。

为了让静态物体也能受到动态部分环境光的影响，将 Light Map 中对应位置的着色信息解析出来，并与动态光照部分计算出来的着色进行混合。

这里对Unity、UE4所实现的 Light Map 做一个更深入分析：

• Light Map 一般记录的是受静态部分环境光后的着色结果，着色方法不限于Ray Tracing、辐射度、Shadow、环境光遮蔽（Ambient Occlusion，AO）等算法。

比较有意思的是，对于颜色存储的方案，Unity采用了直观的线性颜色存储方式，而UE4则采用了非线性映射的颜色存储方式：暗色区域的颜色可以有更多细节（更高精度）而明色区域的颜色则更少细节（更低精度）。这是因为在实际画面中，我们更容易注意到暗色区域的细节而不容易注意到明亮区域的细节。

• Light Map 也可以直接混合写入正常纹理，这样物体着色只需要采样一个纹理，但会极度影响纹理的重用（有别的物体使用了相同的纹理，但是光照结果却是截然不同的；有些着色算法依赖原生纹理而非混合纹理的数据）

• Light Map 对于带 Normal Map 的物体，可能无法得到正确的结果（Normal Map没有起作用）。这是因为光照烘焙系统只使用了物体的顶点数据（而没有使用 Normal Map）。而 Unity 和 UE4 的 Light Map 在启动了高质量选项后会多额外一张纹理（称为 Directional Light Map ）来存储入射光的主要方向（世界坐标空间）：在对某个 shading point 着色时，通过解析出来的入射光主要方向与 shading point 的 normal（Normal Map 解压出来的，且需变换成世界坐标空间下的）进行点乘，便可以得到对应 Light Map 颜色的贡献权重。

Lighting Scenarios 实现昼夜天气系统

---

Lighting Scenarios ：一个物体使用多套 Light Map，常用于模拟天气系统以及一天中不同时段的光照。实际上，就是烘焙几个关键时间点的 Light Map， 根据时间进行线性插值。

UE4 的 Lighting Scenarios 效果：

总结

---

个人归纳：

• IBL 是半动态的 GI 方法（环境光可动可静、动态物体、材质仅粗糙度参数可变），常用于 Specular 物体的渲染
• PRT 是半动态的 GI 方法（环境光可动可静、即使静态环境光还可支持旋转、物体不可形变、静态材质）， 常用于 Diffuse/Glossy 物体的渲染
• Light Probe 是半动态的 GI 方法（静态环境光、动态物体）， 实时给场景中任何物体提供合适的静态环境光信息
• Light Map 是全静态的 GI 方法（静态环境光、静态物体），预先给场景中静态物体提供受环境光静态部分影响后的着色结果

参考

---

• [1] GAMES202-高质量实时渲染-闫令琪
• [2] 《Real-Time Rendering 4th Edition》
• [3] Precomputed Radiance Transfer for Real-Time Rendering in Dynamic, Low-Frequency Lighting Environments
• [4] Stupid Spherical Harmonics (SH) Tricks
• [5] Probe-Based Global Illumination | 知乎
• [6] 游戏中的全局光照(四) 漫反射GI
• [7] Light probe interpolation using tetrahedral tessellations | GDC2012
• [8] 虚幻引擎学习之路：渲染模块之全局光照明
• [9] UE4 Lightmap的解码