- Irradiance Environment Map基本原理
Irradiance Environment Map（也叫Irradiance Map或Diffuse Environment Map），屬於Image Based Lighting技術中的一種。
Irradiance Map的詳細定義可參考GPU Gems 2 Chapter 10.“Real-Time Computation of Dynamic Irradiance Environment Maps”。簡單說來就是一種用於近似Environment Diffuse Lighting的方法。想象一個場景中有k個方向光，方向分別為d1…dk，光照強度為i1…ik，對於一個法線和Diffuse Color分別為n和c的Lambert表面，其光照強度為：
對於Environment Lighting，我們可以用一個Cube Map來表示，Cube Map里的每一個texel就是一個方向光，光強度為texel的值，方向為texel的location。這樣就能通過一個Cube Map來表示任意的Environment Lighting。一般把這個Cube Map叫做Light Probe。
對於Lambert表面，其光照強度只和法線n和光照方向l相關，所以給定一個Light Probe，可以計算出所有可能的法線方向的光照，然后存儲到一個Cube Map里，渲染時，只需要使用法線n去這個Cube Map里索引就能得到Environment Lighting，這個存儲着光照的Cube Map就叫Irradiance Map。計算的偽代碼如下：
diffuseConvolution(outputEnvironmentMap, inputEnvironmentMap) { for_all {T0: outputEnvironmentMap} sum = 0 N = envMap_direction(T0)
for_all {T1: inputEnvironmentMap} L = envMap_direction(T1) I = inputEnvironmentMap[T1] sum += max(0, dot(L, N)) * I
T0 = sum return outputEnvironmentMap }
對於每一個法線n都需要去遍歷所有的光線方向，算法復雜度為O(NM)，N為Light Probe的大小，M為Irradiance Map的大小。
- Spherical Harmonics
由於Diffuse光照本身是變化很緩慢的低頻數據，所以可以使用SH來加速計算。把算法分為兩步：
1. 把Light Probe投影到SH上，求解出SH系數存儲下來。
2. 將Light Probe的SH和Diffuse Transfer的SH做卷積即可求出Irradiance Map。
具體做法如下：
//將LightProbe投影到SH上 for each texel of the lightProbe { lightSample = texelRadiance; weight += texelSolidAngle; //計算光照方向 l = texelDirection; //根據光照方向求出SH基函數 SHBasis = calculateSHBasis(l); //累加SH系數 lightSH += lightSample*SHBasis*texelSolidAngle; } lightSH = lightSH*4*PI/weight; for each texel of the irradianceMap { //法線方向 n = texelDirection; // 求出cosine lobe的SH diffuseSH = calculateDiffuseSH(n); // 用cosine lobe的SH和light probe的SH做卷積 irradiance = dotSH(diffuseSH, lightSH); // lambert brdf irradiance *= 1/PI; texelValue = irradiance; }
使用SH來計算的話，Light Probe和Irradiance Map只需要分別遍歷一遍，所以算法復雜度為O(KN+KM)，N為Light Probe的大小，M為Irradiance Map的大小，其中K為SH系數的個數，對於Diffuse光照，使用3階的SH函數就能獲得不錯的近似結果，3階的SH有9個系數，所以K遠小於N和M。
因為Diffuse光照本身是低頻的，所以輸出的Irradiance Map可以使用較小的分辨率，那么整個算法的開銷主要是在第一步——把Light Probe投影到 SH上。
- 使用GPU計算Light Probe SH
GPU Gems 2 的Chapter 10介紹了使用pixel shader來計算Light Probe SH的方法，使用SM5.0的Compute Shader來計算可以獲得更大的加速比。
觀察求解Light Probe SH的過程——遍歷所有的texel，對於每個texel，求解出SH，然后累加，最后累加的結果就是SH系數。如果使用並行的算法，偽代碼如下：
g_mutex; for each texel of the lightProbe { lightSample = texelRadiance; //計算光照方向 l = texelDirection; //根據光照方向求出SH基函數 SHBasis = calculateSHBasis(l); //累加SH系數 g_mutex.lock(); lightSH += lightSample*SHBasis*texelSolidAngle; weight += texelSolidAngle; g_mutex.unlock(); } lightSH = lightSH*4*PI/weight;
但是對於GPU的Compute Shader，只能同步一個Group里的thread，對於不同Group的thread無法同步，所以無法使用一個加鎖的全局變量不斷累加的方法。
既然無法一次求出所有的累加結果，那么就先求出每個Group的累加結果，然后根據GroupID寫入到輸出的Buffer，然后把這個Buffer作為輸入，重復之前的操作，直到輸出Buffer的Size為1時就求出了結果。比如對於一個512x512的Cube Map，固定Thread Group大小為8x8，那么我們分配64x64個Group，其輸出Buffer大小為64x64（每一個Group輸出一個結果），運行Compute Shader計算結果輸出到Buffer，這時數據就縮小到了64x64，然后重復之前的操作，下一輪的數據大小就變成了8x8（64/8）。一直重復這個操作，直到輸出的Buffer大小為1時就求解出了結果。
這個算法的思路和HDR渲染中求解場景的平均亮度是一樣的，在求解平均亮度時，每次把Texture的Size縮小到1/4做Down Sample，直到Texture大小為1x1時就求出了平均亮度。
求解每個Group的SH累加結果的算法參考Nvidia的 “Optimizing Parallel Reduction in CUDA”， Parallel Reduction的思路就是一個遞歸的tree-based approach，如下圖
對於Shader代碼，使用循環來模擬這個過程，具體做法是設置一個步長step，把相隔step個步長的數據相加，然后step乘以2，重復這個過程，直到step大於N，N為輸入數據的大小，循環累加的代碼如下：
for (uint s = 1; s < groupthreads; s *= 2) // stride: 1, 2, 4, 8, 16,32, 64, 128 { int index = 2 * s * GI; if (index < (groupthreads)) sharedMem[index] += sharedMem[index + s]; GroupMemoryBarrierWithGroupSync(); }
算法的流程如下圖：
完整的Shader代碼如下：
#define THREAD_SIZE_X 8 #define THREAD_SIZE_Y 8 #define GROUP_THREADS THREAD_SIZE_X*THREAD_SIZE_Y Texture2D<float> g_InputBuffer : register(t0); RWTexture2D<float> g_OutputBuffer : register(u0); groupshared float g_ShareMem[GROUP_THREADS]; [numthreads(THREAD_SIZE_X, THREAD_SIZE_Y, 1)] void ReductionCS(uint3 Gid : SV_GroupID, uint3 DTid : SV_DispatchThreadID, uint3 GTid : SV_GroupThreadID, uint GI : SV_GroupIndex) { //加載數據到share memory中 uint Idx = DTid.y*g_InputBuffer.Length.x + DTid.x; g_ShareMem[GI.xy] = g_InputBuffer[Idx]; GroupMemoryBarrierWithGroupSync(); // 循環累加所有數據 [unroll] for (uint s = 1; s < GROUP_THREADS; s *= 2) // stride: 1, 2, 4, 8, 16, 32, 64, 128 { int index = 2 * s * GI; if (index < GROUP_THREADS) g_ShareMem[index] += g_ShareMem[index + s]; GroupMemoryBarrierWithGroupSync(); } if (GI == 0) { //寫入結果 g_OutputBuffer[Gid.xy] = g_ShareMem[0]; } }
NV的paper中還提到了可以進一步優化，GPU在運行Thread Group時，會把Thread划分為Warp，Nvidia的GPU中一個Warp包含32個Thread，這些線程是由SIMD32處理器同步運行的，所以如果線程數目小於32時，可以去掉GroupMemoryBarrierWithGroupSync() 的調用來提升性能。（對於AMD的GPU，把線程划分為Wavefronts，和NV的Warp對應，AMD的每個Wavefronts中包含64個同步執行的Thread）
- 運行結果
運行的參考對象是DirectX ToolKit中的SHProjectCubeMap，這個函數使用CPU來計算Light Probe SH。測試使用的Light Probe是一個512x512的R16G16B16A16_FLOAT格式的HDR Cube Map，測試結果發現在Release模式下使用Compute Shader可以比DXTK的CPU版本快10倍以上，加速比很高。
- 其他
1. Light Probe一般都是HDR格式的，所以生成的Irradiance Map也是HDR格式的，那么使用Irradiance Map計算光照時需要使用HDR渲染。
2. Diffuse Lighting一般用3階的SH足矣，對於HDR的Light Probe，5階的SH會更准確（4階的Diffuse Transer SH系數為0，所以無需使用）。當然也意味着更多的計算量。
3. 使用GPU計算的時候，計算的結果參考了DXTK的函數，但是發現和DXTK的結果有些許偏差。仔細調試后發現是浮點累加誤差導致的，DXTK的結果並不准確，在其SHProjectCubeMap源碼實現中，是使用float類型不斷累加，比如對於求解solid angle weight，最后的結果會累加到幾百萬，但是對於單個的weight數據，其值可能只有1左右，而float有效數字只有6位，所以在累加的過程中會有誤差，數據到幾百萬時，誤差可能會達到幾十。理論上來講，solid angle weight是使用texture uv求解的，那么對於Cube Map的每一個face，其累加結果應該都是相同的，DXTK因為使用float累加，浮點誤差導致其每個face的solid angle weight累加結果並不相同，所以DXTK的結果並不准確，只能作為參考值。
實際上相比於DXTK，GPU算法的准確度更高，因為調試的時候發現GPU求解出來的結果，每一個face的solid angle累加結果都是相同的，不會像DXTK那樣有誤差。因為在GPU算法中，累加是分層求解的（參考前面的算法圖解），每一層中的節點數據范圍都差不多，這樣浮點的誤差就不會像直接累加那么大，所以GPU求解的速度和精度都好於DXTK的SHProjectCubeMap。
參考資料：
http://http.developer.nvidia.com/GPUGems2/gpugems2_chapter10.html
http://developer.download.nvidia.com/assets/cuda/files/reduction.pdf
http://diaryofagraphicsprogrammer.blogspot.com/2014/03/compute-shader-optimizations-for-amd.html
https://seblagarde.wordpress.com/2012/06/10/amd-cubemapgen-for-physically-based-rendering/