作者：i_dovelemon

日期：2020-01-04

主題：Rendering Equation,Irradiance Environment Map,Spherical Harmonic

引言

        在實時圖形渲染中，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 的對比：

兩種方法對比

        我們假設原始 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

[2] Cubemap Texel Solid Angle

[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