基於預計算的全局光照（Global Illumination Based On Precomputation）

• 2021 年 9 月 29 日
• 筆記
• 圖形學

基於影像的光照(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}\right)=\int_{\Omega^{+}} L_{i}\left(x, \omega_{i}\right) f_{r}\left(x, \omega_{i}, \omega_{r}\right)\left(n \cdot \omega_{i}\right) \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}\right) \approx \frac{\int_{\Omega_{f_{r}}} L_{i}\left(x, \omega_{i}\right) \mathrm{d} \omega_{i}}{\int_{\Omega_{f_{r}}} \mathrm{~d} \omega_{i}} \cdot \int_{\Omega^{+}} f_{r}\left(x, \omega_{i}, \omega_{r}\right) \cos \theta_{i} \mathrm{~d} \omega_{i}\)

對渲染方程拆分成兩個部分（環境光積分、BRDF積分）後就可以通過預計算的方式（後面兩節會介紹如何預計算）分別減少這些積分的運行時開銷，總結這種方法的好處是：

• Split Sum 方法和原始蒙特卡洛方法的影像效果幾乎一模一樣
• 由於不用對環境貼圖進行多重取樣，性能開銷大大減低了

過濾環境貼圖

環境光積分：\(\frac{\int_{\Omega_{f_{r}}} L_{i}\left(x, \omega_{i}\right) \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\theta_{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}}\right)} = \frac{1}{\pi \alpha^{2}\cos^4 \theta_h } \cdot \exp {\left(\frac{\cos^2 \theta_h -1}{a^{2}\cos^2 \theta_h }\right)}\)

可以預想到，該BRDF 的積分結果依賴於三個參數：

• \(F_0\)（Fresnel項係數）
• \(\alpha\) （粗糙度 roughness）
• \(\theta_v\) （反射方向與法線的夾角，實際上決定了 \(\theta_{vh}\)、\(\theta_{h}\)）

我們可以把 \(F_0\) 項拆出來，讓BRDF 積分拆成兩個積分，但是這些積分都減少了一個依賴的參數：

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

這樣，對於相同的材質（相同的BRDF），我們就可以針對剩餘兩個依賴的參數 \(\alpha\) 、\(\theta_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_1\cdot 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}(\theta, \varphi)=\left\{\begin{array}{rl}\sqrt{2} \operatorname{Re}\left(Y_{l}^{m}\right) & m>0 \\ \sqrt{2} \operatorname{Im}\left(Y_{l}^{m}\right) & m<0 \\ Y_{l}^{0} & m=0\end{array}=\left\{\begin{array}{cc}\sqrt{2} K_{l}^{m} \cos m \varphi P_{l}^{m}(\cos \theta) & m>0 \\ \sqrt{2} K_{l}^{m} \sin |m| \varphi P_{l}^{|m|}(\cos \theta) & m<0 \\ K_{l}^{0} P_{l}^{0}(\cos \theta) & m=0\end{array}\right.\right.\)

其中，

\(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}\right)=\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}\right)=\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}\right)=\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}\right)=\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}\right)=\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}\right)=\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\neq j)\)

• 通過 投影（Projection）可以很方便得到 SH 係數（SH coefficients）

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

• 通過係數向量組與基函數組的點積（前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} = M\cdot \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)d\omega\)

這樣我們就可以預計算好 \(M\) ，等到知道 \(\vec{b}\) 後，乘起來就能得到 \(\vec{c}\)

推導：

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

\(\begin{align} c_i &= \int_{\Omega} a(\omega)b(\omega)B_i(\omega) d\omega \\ &= \int_{\Omega} M_i(\omega)b(\omega)d\omega \\ &= \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 的物體。

球諧光照（Spherical Harmonic Lighting）效果：

圖分別為：diffuse物體的unshadowed效果、diffuse物體的shadowed&interreflection效果、glossy物體的unshadowed效果、glossy物體的shadowed&interreflection效果

Diffuse 物體的球諧光照

物體的渲染方程：

\(L(r)=\int_{\Omega^+} L(\mathbf{\omega})\cdot \rho \cdot V(\mathbf{\omega}) \max (0, n\cdot \mathbf{\omega}) \mathrm{d} \mathbf{\omega}\)

其中，\(\rho\) 為 BRDF 項；\(V\) 為 Visibility，表示不被遮擋的程度，往往表現為自遮擋產生的陰影現象。

• 由於物體是 diffuse 的，因此它的 BRDF 將是一個常數 \(\rho\) （無論從哪個方向觀察都得到相同的BRDF值）
• 對於 lighting 部分 \(L(\mathbf{\omega})\)，原本需要對 Environment Map 進行查詢，而現在可以換成使用 SH 函數去表示：\(L(\mathbf{\omega}) \approx \sum l_{i} B_{i}(\mathbf{\omega})\)
• 對於 transfer function 部分 \(T(\omega)=V(\mathbf{\omega}) \max (0, n\cdot \mathbf{\omega})\) ， 也可以換成使用 SH 函數去表示：\(T(\omega) \approx \sum T_j B_j(\omega)\)

代入渲染方程整理後得：

\(\begin{align} L(r) &\approx \frac{\rho}{\pi} \int_{\Omega^+} \sum_{i=1}^{n^2}l_{i}B_i(\omega) \sum_{j=1}^{n^2}T_j\mathrm{B}_{j}(\mathbf{\omega}) \mathrm{d} \mathbf{\omega} \\ &= \frac{\rho}{\pi}\sum_{i=1}^{n^2} \sum_{j=1}^{n^2} l_{i}T_j\int_{\Omega^+} B_i(\omega) \mathrm{B}_{j}(\mathbf{\omega})d\omega \\ &= \frac{\rho}{\pi}\sum_{i=1}^{n^2} l_{i}T_i\int_{\Omega^+} B_i(\omega) \mathrm{B}_{i}(\mathbf{\omega})d\omega \\ &= \frac{\rho}{\pi}\sum_{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(\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, n\cdot \mathbf{\omega}) \mathrm{B}_{i}(\mathbf{\omega}) \mathrm{d} \mathbf{\omega}\)

此時，我們可以理解成係數 \(l_i\) 代表了環境光照（lighting）的資訊，而係數 \(T_i\) 代表了光路傳輸（light transfer）的資訊

運行時渲染過程：

1. 在 vertex shding 階段計算頂點的 SH 顏色 \(L(r) \approx \frac{\rho}{\pi} \sum_{i=1}^{n^2} l_{i} T_{i} = \frac{\rho}{\pi}( \vec{l}\cdot \vec{T})\)

2. 在 pixel/fragment shading 階段得到插值後的 SH 顏色即為該像素的顏色

Glossy 物體的球諧光照

物體的渲染方程：

\(L(r)=\int_{\Omega^+} L(\mathbf{\omega})\cdot \rho(r,\mathbf{\omega}) \cdot V(\mathbf{\omega}) \cdot \max (0, n\cdot \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(\mathbf{\omega}) \approx \sum l_{i} B_{i}(\mathbf{\omega})\)
• 對於傳輸函數（transfer function） \(T(r,\omega)=\rho(r,\mathbf{\omega}) V(\mathbf{\omega}) \max (0, n\cdot \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 情況下的係數向量。

代入渲染方程整理後得：

\(\begin{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}) \cdot 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}\)

運行時渲染過程：

1. 在 vertex shding 階段計算頂點的 SH 顏色 \(L(r) \approx \sum_{i=1}^{n^2}\sum_{j=1}^{n^2} l_{i} T_{ij} B_j(r) = T\vec{l}\cdot \vec{B(r)}\)
2. 在 pixel/fragment shading 階段得到插值後的 SH 顏色即為該像素的顏色

球諧光照預計算 Transfer

PRT 的 SH Lighting 預計算是針對單個模型物體的預計算，即輸入一個 模型物體 + Environment Map，然後輸出模型中每個頂點的 Lighting SH 係數和 Transfer SH 係數或矩陣。這樣也意味著這種 GI 方式的 Transfer 只考慮了模型自身的幾何遮蔽關係，並不會考慮其它物體。

TODO

球諧函數的旋轉（SH Rotation）

PRT 的一個問題是如果 lighting 部分是預計算的，那就只適用於靜態環境光下的靜態物體渲染；環境光或者物體只要有變化，PRT 就不得不進行重新預計算；但得益於 SH 的旋轉不變性，我們至少可以讓 SH Lighting 適用於動態旋轉的情形而不必重新預計算。

環境光旋轉往往用來實現場景的晝夜輪轉效果，物體旋轉就比較常見了。順便註：環境光旋轉和物體旋轉在 PRT 渲染中是等價的（只是說看相對於哪個東西來看待旋轉而已）

現在給定旋轉 \(R\) ，然後有SH投影函數 \(P\)（輸入一個球面向量，輸出第 \(l\) 層band的 SH 係數向量組）。假設我們有 \(2l+1\) 個任意的球面向量，我們需要想辦法求出能等價於旋轉 \(n_i\) 的矩陣 \(M\) 來旋轉 SH投影函數 \(P\) ：

\(M P\left(n_{i}\right)=P\left(R\left(n_{i}\right)\right), i \in[-l, l]\)

整理得：

\(M\left[P\left(n_{-l}\right), \ldots, P\left(n_{l}\right)\right]=\left[P\left(R\left(n_{-l}\right)\right), \ldots, P\left(R\left(n_{l}\right)\right)\right]\)

記 \(A = [P_{(n−l)}, …, P (n_{l})]\) ，如果矩陣 \(A\) 是可逆的，則：

\(M=\left[P\left(R\left(n_{-l}\right)\right), \ldots, P\left(R\left(n_{l}\right)\right)\right] 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 問題。

這裡列舉一些遊戲的做法：

• 《原神》將室內外的 light probe 分開烘焙，使用一個 volume 標記室內區域，室內外區域分別使用對應的 light probe 數據，同時在門口區域，進行過渡處理
• 對馬島之魂：為每個 light probe 設置室內外標記，記錄權重值 。在四面體網格中計算重心坐標時，會傳入一個表示室內外權重的權重值，根據權重值和四面體上四個點的室內外權重值，重新計算出新的重心坐標。
• 使命召喚無盡戰爭：對每個四面體，記錄每個頂點相對面上的遮擋資訊，計算時考慮是否被遮擋。
• 底特律變人：使用 lightmap grid 八叉樹模式保存的light probe，如果遇到在牆兩側的light probe，則進行再次細分，並將取樣光照的位置進行虛擬偏移：

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[\begin{array}{l}a \\ b \\ c\end{array}\right]=\left[\begin{array}{lll}\vec{P}_{0}-\vec{P}_{3} & \vec{P}_{1}-\vec{P}_{3} & \overrightarrow{P_{2}}-\vec{P}_{3}\end{array}\right]^{-1}\left[\vec{P}-\overrightarrow{P_{3}}\right]\)

\(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 效果：


總結

GI 總體上分為兩派：一派是通過一套統一標準的渲染流程把任何物體的 GI 計算出來，往往計算量極大但效果更加物理更加真實，典型的例子便是離線光線追蹤；另一派是把 GI 的效果看成一個個部分來組成，這樣我們可以選擇其中一些GI效果組合使用，適應不同的性能（往往計算量相對低些，尤其是實時渲染）的同時也能帶來能接受的 GI 效果（雖然往往不是嚴謹的物理正確），典型的例子就是這篇部落格所介紹的方法和另一些完全動態的 GI 方法（如 SSAO）。

個人歸納：

• 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的解碼