作者：i_dovelemon
日期：2020-01-04
主題：Rendering Equation,Irradiance Environment Map,Spherical Harmonic
ChangeLog-2020/01/11: 添加 Light Probe Blend 相關描述
引言
在實時圖形渲染中，Global Illumination 是聖杯級的效果。為了實現這個效果，前輩們開發了很多的技術。但是這些技術大都只能用於靜態物體上，對於動態的物體卻不能很好的支持。所以，為了讓動態的物體也有一點 GI 的效果，開發出了一系列的技術。今天，我們就來介紹其中一種技術：Diffuse Irradiance Environment Map。在游戲開發領域，一般稱之為 Light Probe(注：當然 Light Probe 能夠實現更多的效果，Diffuse 的 GI 是其中一種)。
Diffuse Irradiance Environment Map 是基於 Environment Map 來實現的。所以，它不會考慮陰影和模型本身的光照影響。同時，我們也只探討光照中的 Diffuse 部分，即 Lambert BRDF 部分。
文章中會存在大量的渲染相關的術語，諸如 irradiance，radiance，solid angle 等等。我們假設你已經了解了這些基礎性的概念知識，如果不是，PBRT [文獻1] 是一個很好的參考資料。
本文將主要從兩個方面來講述：一個是傳統的計算 Diffuse Irradiance Environment Map 的方法，我們稱之為 Brute force；另外一種是基於 Spherical Harmonic 的方法。
背景知識
我們回顧下渲染方程，可以知道一個點在半球范圍里面受到的 irradiance 為：
$E(\vec{n})=\int_{\Omega(\vec{n})}^{ }L(\vec{w})(\vec{n}\cdot\vec{w})d\vec{w} \ \ \ \ \ \ \ \ (1)$
也就是說，對於一個固定的 Environment Map (即 $L(\vec{w})$ 相同)來說，irradiance 只和 normal 有關。所以，我們可以通過預計算，將 Environment Map 對應的 Irradiance Environment Map 保存為一個和 normal 相映射的形式，然后通過頂點的 normal 來獲取對應的 irradiance 信息。獲取到 irradiance 信息之后，帶入下面的公式，就能夠得到最終需要顯示的顏色值：
$B(\vec{p},\vec{n})=f(\vec{p})E(\vec{n}) \ \ \ \ \ \ \ \ (2)$
其中 $f(\vec{p})$ 表示的是 Diffuse 的 BRDF。
BruteForce 方法
公式(1)中計算 irradiance 的方法，是一個在半球范圍里面積分的形式，這種方式不存在解析解，沒有辦法直接去計算得到。但是，由於光照環境是通過 Environment Map 來表達的，我們可以將公式(1)轉化為離散的形態，如下公式所示：
$E(\vec{n})=\sum_{i=0}^{N-1}L(\vec{w})(\vec{n}\cdot\vec{w})d(\vec{w})\ \ \ \ \ \ \ \ (3)$
其中，N 表示的是整張 Environment Map 上的所有像素的數量；$L(\vec{w})$ 表示的是在 $\vec{w}$ 方向上的 radiance；$d(\vec{w})$ 表示的是在 $\vec{w}$ 方向上像素的 solid angle。
這樣，我們就有了一個方法來實際計算一個 normal 方向上的 irradiance 的值了。
我們知道了如何計算一個 normal 對應的 irradiance 的值，那么這個值該怎么保存了？很明顯的，我們可以利用另外一張 Cubemap 來保存各個 normal 計算出來的對應的 irradiance 的值，而這個新的 Cubemap 就是 Diffuse Environment Irradiance Map。以下是整個過程的偽代碼：
for pixel_iem in IrradianceEnvironmentMap n = GetNormal(pixel_iem) irradiance = 0 for pixel in EnvironmentMap L = GetRadiance(pixel) w = GetRadianceDir(pixel) dw = GetTexelSolidAngle(pixel) irradiance += L * max(0, dot(n,w)) * dw pixel_iem = irradiance
Cubemap Texel Solid Angle
上面的代碼中，唯一可能比較難計算的是：GetTexelSolidAngle。[文獻2] 中詳細的解釋了如何定義這個函數，以及該函數背后的數學原理，這里給出實際的代碼，不再贅述：
private static float AreaElement(float x, float y) { return Mathf.Atan2(x * y, Mathf.Sqrt(x * x + y * y + 1.0f)); } private static float TexelCoordSolidAngle(float x, float y, int size) { // Scale up to [-1,1] range (inclusive), offset by 0.5 to point to texel center float u = 2.0f * (x + 0.5f) / size - 1.0f; float v = 2.0f * (y + 0.5f) / size - 1.0f; float invRes = 1.0f / size; // Project area for this texel float x0 = u - invRes; float y0 = v - invRes; float x1 = u + invRes; float y1 = v + invRes; return AreaElement(x0, y0) - AreaElement(x0, y1) - AreaElement(x1, y0) + AreaElement(x1, y1); }
我們知道，solid angle 在整個球上的積分值為 $4\pi$。前面我們將公式(1)轉化成了離散的形式，這樣就導致所有像素的 solid angle 的總和與 $4\pi$ 存在一定的誤差，所以需要進行修正。修正之后的偽代碼如下所示：
for pixel_iem in IrradianceEnvironmentMap n = GetNormal(pixel_iem) irradiance = 0 totalSolidAngle = 0 for pixel in EnvironmentMap L = GetRadiance(pixel) w = GetRadianceDir(pixel) dw = GetTexelSolidAngle(pixel) irradiance += L * max(0, dot(n,w)) * dw totalSolidAngle += dw pixel_iem = irradiance * 4 * PI / totalSolidAngle
好了，至此我們就得到了一張 Diffuse Irradiance Environment Map。在渲染的時候，我們只要通過像素的 normal 來采樣 Irradiance Environment Map 就可以得到對應的 irradiance。然后帶入公式(2)中，得到最終需要顯示的顏色值。
To PI or not to PI?
這里有一個容易引起困惑的地方。我們知道，Lambert 光照模型的 BRDF 如下所示：
$f = \frac{c_{diff}}{\pi}\ \ \ \ \ \ \ \ (4)$
而有游戲開發經驗的同學就知道，在游戲里面我們定義 Diffuse 的光照模型如下所示：
$B=c_{diff}*c_{light}*max(0,dot(\vec{n},\vec{w}))\ \ \ \ \ \ \ \ (5)$
這里卻沒有 $\pi$ 相關的值。這是因為在傳統的游戲里面，我們定義的 $c_{light}$ 並不是以光學輻射度的單位來定義的，而是以一種對美術更加友好的定義方式：當一個純白的 Lambert 表面被一束平行於表面 normal 的光所照射時所呈現的顏色為 $c_{light}$。也就是說，傳統游戲開發中定義的 $c_{light}$，實際上是真實光學輻射度單位輸入除以 $\pi$ 之后的結果，所以公式(5)中就不存在 $\pi$。
說這么多的意思是，我們定義 Environment Map 是以真實的輻射度單位來保存的，也就是說在計算最終顏色的時候，我們需要自己除以 $\pi$ 來保證結果的正確性，即將公式(4)帶入公式(2)中計算最終的顏色值，即：
$B=\frac{c_{diff}}{\pi}E(\vec{n})\ \ \ \ \ \ \ \ (6)$
這里為了簡化 shader 中的計算，我們將 $\pi$ 的計算放在了 Diffuse Irradiance Environment Map 里，即：
for pixel_iem in IrradianceEnvironmentMap n = GetNormal(pixel_iem) irradiance = 0 totalSolidAngle = 0 for pixel in EnvironmentMap L = GetRadiance(pixel) w = GetRadianceDir(pixel) dw = GetTexelSolidAngle(pixel) irradiance += L * max(0, dot(n,w)) * dw totalSolidAngle += dw pixel_iem = irradiance * 4 * PI / totalSolidAngle pixel_iem = pixel_iem / PI
關於 $\pi$ 的詳細討論可以參考[文獻3]。
結果
以下是一些通過 BruteForce 方法計算出來的 Diffuse Irradiance Environment Map 和原始 Environment Map 的對比結果圖，Diffuse Irradiance Environment Map 大小是 32x32：
Spherical Harmonic 方法
Spherical Harmonic 是信號處理里面的一種變換方法。和 Fourier 變換相似，都是將信號轉化到頻域中去，以此來更加精簡的表達原始復雜的信號。不同的是，Spherical Harmonic 更加適合用來處理球面相關的信號。而渲染相關的問題，都是在一個球面范圍里面進行，所以選擇使用 SH 的方法。關於 SH 的描述，[文獻4] 講解的非常詳細，這里就不再贅述。神奇的地方在於，BruteForce 的方法得到的是最終的 Diffuse Irradiance Environment Map，而 SH 的方法得到的是 SH 系數（一般是9個系數）。然后在實際渲染的時候，我們根據這9個系數，重建原始的信號，得到對應的 irradiance。
Prefilter
根據[文獻5]中的描述，我們知道如果使用 SH coefficient 的表示方法來編碼 Environment Map 的話，將使用如下的公式：
$L(\theta,\phi)=\sum_{l,m}^{ }L_{lm}Y_{lm}(\theta,\phi)\ \ \ \ \ \ \ \ (7)$
而同樣的，使用 SH 編碼 Irradiance Environment Map 的話，將使用如下的公式：
$E(\theta,\phi)=\sum_{l,m}^{ }E_{lm}Y_{lm}(\theta,\phi)\ \ \ \ \ \ \ \ (8)$
同時定義：
$A=(\vec{n}\cdot\vec{w})$
$A(\theta)=\sum_l^{ }A_lY_{l0}(\theta)\ \ \ \ \ \ \ \ (9)$
根據上面的定義，我們得到：
$E_{lm}=\sqrt{\frac{4\pi}{2l+1}}A_lL_{lm}\ \ \ \ \ \ \ \ (10)$
引入新的變量：
$\hat{A}_l=\sqrt{\frac{4\pi}{2l+1}}A_l\ \ \ \ \ \ \ \ (11)$
將公式(9)(10)(11)帶入公式(8)，得到：
$E(\theta,\phi)=\sum_{l,m}^{ }\hat{A}_lL_{lm}Y_{lm}(\theta,\phi)\ \ \ \ \ \ \ \ (12)$
公式(12)中，$\hat{A}_l$是可以預先計算出來的，$Y_{lm}(\theta,\phi)$ 通過帶入 normal，也能夠計算出來，只有 $L_{lm}$ 是未知的。所以，我們 Prefilter 操作的目的就是計算出 $L_{lm}$ 的值。
[文獻5]中講述了我們只需要3階的 SH 系數，就能夠很好的表達信號，所以我們只需要計算出來 $l <= 2$ 的 SH 系數即可。
根據文獻[4]中的描述，計算 SH 系數的方式就是將信號投影到對於的基向量上去即可，即：
$L_{lm}=\int_{\Omega}^{ }L(\vec{w})Y_{lm}(\vec{w})d(\vec{w})\ \ \ \ \ \ \ \ (13)$
同樣的，我們將這個積分形式的方程，轉化為離散的形式[文獻6]，如下所示：
$L_{lm}=\sum_{i=0}^{N-1}L(\vec{w})Y_{lm}(\vec{w})d(\vec{w})\ \ \ \ \ \ \ \ (14)$
這樣，我們就能夠通過計算，得到3階球諧的9個系數。但是我們知道光的單位是有RGB三個部分組成，每一個部分可以單獨的進行 SH 系數的求解，所以最終的結果是9個RGB系數。以下是求解這些系數的偽代碼：
foreach sh_coefficient sh_coefficient = 0 totalSolidAngle = 0 for pixel in EnvironmentMap L = GetRadiance(pixel) sh = GetSHBais(pixel) dw = GetTexelSolidAngle(pixel) sh_coefficient += L * sh * dw totalSolidAngle += dw sh_coefficient = sh_coefficient * 4 * PI / totalSolidAngle
計算過程十分簡單，唯一需要注意的點是GetSHBais 函數的實現。這個函數的定義可以通過預先計算得到，如下所示[文獻5]：
至此，Prefilter 的工作就完成了。
Rendering
通過公式(12)，我們在知道了 SH 系數的情況下，就可以重建原始的 irradiance 信號。由於公式(12)中只有 $L_{lm}$ 是未知的，其他兩個部分都是可以通過預計算得到，所以合並預計算的部分，我們得到根據 SH 系數重建信號的公式[文獻5]：
這樣，在知道了 irradiance 的情況下，帶入到公式(2)中，得到：
$B=\frac{c_{diff}}{\pi}E\ \ \ \ \ \ \ \ (15)$
結果
以下是通過 SH 的方法得到的 Diffuse Irradiance Environment Map 和原始 Environment Map 的對比：
Light Probe Blend
上面的描述中，我們都是從一個特定的點去繪制環境貼圖，然后以此環境貼圖為周圍光照的 radiance 描述來構建 light probe。對於 reflection environment map 來說，我們可以通過簡單的數學計算，就能夠得到一個 localized 的效果（見文獻[8]）。但是對於 light probe 來說，由於需要在半球范圍里面進行積分運算，無法使用文獻[8]中描述的方式實現 localized 的效果。但是，我們可以通過對兩張在不同的點生成的 light probe 進行插值，以此來構建他們之間某個點的 light probe，從而實現一定程度上的 localized 的效果。
對於 SH 表達的 light probe 來說，在兩個 light probe 之間進行插值就是簡單對 SH coefficient 進行插值即可（見文獻[9]）:
當然在實際的游戲開發中，你不可能只對兩個 light probe 進行插值。為了較好的表達整個場景的效果，我們需要在場景中擺放很多個 light probe。然后根據被渲染物體所在的位置選擇一個或多個 light probe，根據一定的權重進行插值。這部分的知識也十分復雜，可以參考文獻[9]中的具體描述。文獻[9]描述了常見的系統設計方案，以及它所存在的問題，同時給出了 Unity 對此進行的改進和它所使用的方案。感興趣的讀者可以去了解。
兩種方法對比
我們假設原始 Environment Map 的尺寸是 NxNx6，而 BruteForce 方法計算得到的 Diffuse Irradiance Environment Map 的尺寸為 MxMx6，那么對於 BruteForce 的方法來說，就是一個 O(NxNxMxM) 的操作。而對於 SH 方法來說，它的計算時間為 O(9xNxN)。兩個方法在Prefilter上面，SH 的速度大大提高。同時，對於 BruteForce 方法來說，得到的結果是一張 Cubemap，在渲染的時候需要進行采樣，而 SH 的方法則是通過一些簡單的計算得到最終的結果。
以下是兩種方式得到的 Diffuse Irradiance Environment Map 的對比：
可以看到，通過 SH 方式得到的結果和 BruteForce 的方法得到的結果誤差非常小。
結論
如果在實際使用過程中，你需要使用 Diffuse Irradiance Environment Map，也是建議通過先求 SH 系數，然后重建 Diffuse Irradiance Environment Map，這樣的方法比 BruteForce 來計算得到 Diffuse Irradiance Environment Map 的速度要快的多。
當然除了這里提到的方法，還有很多其他的方法來計算 Diffuse Irradiance Environment Map，比如[文獻7]中，使用 Rieman 積分的方式，加速 BruteForce 方法來得到結果。
這里只是介紹了基礎的知識，在實際項目開發過程中還需要處理諸多的問題，比如：Light Probe Auto Layout，Light Probe Blend 等等復雜的問題，后面有機會會專門講解這方面的知識。
本文的配套代碼可以在這里獲取得到：https://github.com/idovelemon/UnityProj/tree/62eff639347645f380d651dd80b8720010f6097b/IrradianceEnvironmentMap。值得注意的是，學術界對球面坐標系的定義是 Z 軸朝上，而在 Unity 里面是 Y 軸朝上，實際代碼實現的時候需要轉化下方向。
參考文獻
[1] Physically Based Rendering : From Theory to Implementation
[3] PI or not to PI in game lighting equation
[4] Spherical Harmonic Lighting:The Gritty Details
[5] An Efficient Representation for Irradiance Environment Maps
[6] Real-Time Computation of Dynamics Irradiance Environment Maps
[7] GraphicsLab Project 之 IBL - Diffuse 光照
[9] Light Probe Interpolation using Tetrahedral Tessellations