基於Spherical Harmonics的實時全局照明
本文主要內容：頻率與基函數、球面諧波函數、預計算輻射傳輸、球面均勻采樣、構建正交基、SH實現、球面諧波函數的旋轉及實現。
1. 頻率與基函數
1.1 傅里葉級數展開
任何一個函數可以寫成常數和一系列基函數(不同頻率sin和cos項)的線性組合，基函數數量越多越接近於原函數的形狀：
1.2 頻率
在空間上圖像信號數值的變化是否劇烈，任何一張圖(也就是二維函數)的頻率,也就是頻域上對應的內容可以用一張頻譜表示出來：
頻譜最中心處是低頻內容，我們可以做一個filtering(濾波)，從而去除一系列頻率上的內容，我們對這張圖用一個低通濾波器，從而把高頻的內容去除掉：
卷積：卷積是一個模糊操作，在圖上取任意一點，取點的周圍一定區域內的像素值進行加權平均並將結果寫回這個點，這就是卷積。
在空間域上做卷積等於在函數上做卷積；等於在頻域上做原圖頻譜和卷積核頻譜 的乘積操作：
對於任意的product integral(兩個函數先乘積在積分)，我們將其認為是做了一個卷積操作。
空間域上的兩個信號\(f(x)\)和\(g(x)\)進行一個卷積，等於在頻域上讓兩個信號相乘；如果兩個信號有一個信號是低頻的，那么頻域上相乘后得到的結果也是低頻的，最終相乘在積分的結果也是低頻。的可以總結為：積分之后的頻率取決於積分前最低的頻
1.3 基函數
一個函數可以描繪成其他函數的線性組合，如\(f(x)\)可以描繪成一系列的\(B_i\)函數乘以各自對應的系數最終再相加在一起，這一系列的函數\(B_i\)就是基函數。
2. 球諧函數
SH(Spherical Harmonics)是一系列基函數，系列中的每個函數都是二維函數，並且每個二維函數都是定義在球面上的，階數越高，頻率越高。
下圖是對SH的可視化，與一維的傅里葉一樣，SH也存在不同頻率的函數，但不同頻率的函數個數也不同，頻率越高所含有的基函數越多：
\(l\)表示的是階數，通常第\(l\)階有\(2l+1\)個基函數，前\(n\)階有\(n^2\)個基函數。
越偏白的藍色地方值越大，越黑的地方值越小；而黃色中則表示偏白的地方表示其絕對值大，偏黑的地方表示絕對值小；藍色表示正，黃色表示負。
2.1 投影
由於一個函數\(f(\omega)\)可以由一系列基函數和系數的線性組合表示，那么怎么確定基函數前面的系數，這就需要通過投影操作：
我們知道函數\(f(\omega)\),通過對應的基函數\(B_i(\omega)\)進行投影操作，從而求出各基函數對應的系數\(C_i\)，與以下操作是同一個道理：在空間中想描述一個向量，可以xyz三個坐標來表達，把xyz軸當做三個基函數，把向量投影到xyz軸上，得到三個系數就是三個坐標。
2.2 重建
知道基函數對應的系數，就能用系數和基函數恢復原來的函數。
由於基函數的階可以是無限個的，越高的階可恢復的細節就越好，但一方面是因為更多的系數會帶來更大的存儲壓力、計算壓力，而描述變化比較平滑的環境漫反射部分，用3階SH就足夠。
對於重建diffuse來說這里的\(f(\omega)\)對應環境光照。
2.3 球諧函數的正交性
旋轉一個基函數之后，得到的函數就不再是一個基函數(因為基函數有嚴格的朝向等限制)，但是旋轉球諧函數等價於同階基函數的線性組合。這意味着可以簡單地進行投影、重建、旋轉。
2.4 內容補充
diffuse BRDF類似於一個低通濾波器，使用一些低頻信息就可以恢復出原始內容。因為積分之后的頻率取決於積分前最低的頻率，當diffuse BRDF使用低頻信息即可恢復內容時，也就意味着無論光照項是多么復雜，其本應該用多高頻的基函數去表示，但我們希望得到的是其與BRDF之積的積分，所以可以使用比較低頻的基函數去描述燈光。
前3階球諧基函數：
3. 預計算輻射傳輸
PRT(Precomputed Radiance Transfer)思想：將光照切割成Lighting, Light Transport兩部分：
首先，因為在Diffuse的情況下，BRDF幾乎是一個常數，可以把它提到外面：
而Lighting項可以寫成基函數的形式：
同樣的，Light transport項也可以寫成基函數求和的形式：
\(l_i\)、\(l_q\)是常數，可以提出來：
因為球諧的正交性，只有\(i = q\)的時候，\(B_{i}\left(\mathbf{i}\right) B_{q}\left(\mathbf{i}\right)\)相乘才有值，所以這個函數的復雜度仍是O(n)。
下面將基函數形式的Lighting項帶入渲染方程：
對於積分中的部分來說，相當於light transport乘與一個基函數，這就成了lighting transport投影到一個基函數的系數。因為對light transport做了預計算，visibility成為了常量，只能對靜止物體進行計算。
一個SH基函數旋轉后都可以被同階的SH基函數線性組合得到，因此可以解決光源旋轉問題。
現在將lighting這個球面函數，通過SH的基函數用一堆系數來表示，這些系數排成一行也就是組成了向量，因此光照變成了一個向量：
投影到SH空間：
從SH空間重建：
4. 球面均勻采樣
4.1 Hammersley采樣
Hammersley這是一種均勻分布的2D隨機采樣，他將十進制轉換成二進制，再將二進制轉換到[0,1]之間的小數，這一過程被稱作 Radical Inverse。具體表示方法見下圖：
Hammersley采樣點集合：
\(N\)是我們的采樣點總數，\(\phi(i)\)就是Radical Inverse之后得到的小數，代碼實現如下：
float RadicalInverse( uint bits )
{
//reverse bit
//高低16位換位置
bits = (bits << 16u) | (bits >> 16u);
//A是5的按位取反
bits = ((bits & 0x55555555) << 1u) | ((bits & 0xAAAAAAAA) >> 1u);
//C是3的按位取反
bits = ((bits & 0x33333333) << 2u) | ((bits & 0xCCCCCCCC) >> 2u);
bits = ((bits & 0x0F0F0F0F) << 4u) | ((bits & 0xF0F0F0F0) >> 4u);
bits = ((bits & 0x00FF00FF) << 8u) | ((bits & 0xFF00FF00) >> 8u);
return float(bits) * 2.3283064365386963e-10;
}
vec2 Hammersley(uint i,uint N)
{
return vec2(float(i) / float(N), RadicalInverse(i))
}
4.2 2D坐標到球坐標的映射函數
對於球面來說，如果用方位角\((\theta, \phi)\)直接拿Hammersley均勻采樣，會出現兩極密度高，中間密度低的情況：
通過求球面概率密度分布函數的反函數，可以對球面均勻采樣的映射函數：
現在，在2D平面坐標\((\xi_{x},\xi_{y})\)做Hammersley均勻采樣，然后再映射回球坐標系即可得到均勻的分布：
PS：我們可以很容易的得到 \(cos\theta=1-2 \xi_{x}\)
4.3 球坐標系與直角坐標系的轉換
直角坐標系\((x,y,z)\)轉球坐標系\((r, \theta, \varphi)\)：
球坐標系\((r, \theta, \varphi)\)轉直角坐標系\((x,y,z)\)：
現在我們可以結合2D坐標到球坐標的映射函數寫出輸入\((\xi_{x},\xi_{y})\)，輸出\((x,y,z)\)的方法：
vec3 HemisphereSample_uniform(float u, float v)
{
float phi = v * 2.0 * PI;
float cosTheta = 1.0 - u; // float cosTheta = u;//一回事兒
float sinTheta = sqrt(1.0 - cosTheta * cosTheta);
return vec3(cos(phi) * sinTheta, sin(phi) * sinTheta, cosTheta);
}
注意這里的\(1-2 \xi_{x}\)被改為了\(1- \xi_{x}\)，所以是半球空間。球面空間只需要將float cosTheta = 1.0 - u;改為float cosTheta = 1.0 - 2u;
5. 根據單位法向量構建正交基
根據反射定律和折射定律，我們知道如何求解反射方向和折射方向。用蒙特卡洛方法渲染時，需要生成與法向量或反射方向呈一定概率分布的方向采樣：
通過上文中HemisphereSample_uniform函數可以生成一個直角坐標系下的向量，而且是在正上半球（結合Hammersley可以生成正上半球的均勻采樣向量）。我們可以把它當做是生成在切線空間的，因此我們還需要計算一個正交基，將它轉換回世界空間：
將基下的向量與正交基投影，可以得到原向量，下面提供正交基的一種計算方法：
// [ Duff et al. 2017, "Building an Orthonormal Basis, Revisited" ]
float3x3 GetTangentBasis( float3 TangentZ )
{
const float Sign = TangentZ.z >= 0 ? 1 : -1;
const float a = -rcp( Sign + TangentZ.z );
const float b = TangentZ.x * TangentZ.y * a;
float3 TangentX = { 1 + Sign * a * Pow2( TangentZ.x ), Sign * b, -Sign * TangentZ.x };
float3 TangentY = { b, Sign + a * Pow2( TangentZ.y ), -TangentZ.y };
return float3x3( TangentX, TangentY, TangentZ );
}
6. 實現
6.1 編碼球面紋理信息
以光照探針所在點為相機位置，記錄周圍環境到Cubemap：
之后的操作有兩種。閆老師提到，濾波后一次采樣等於沒濾波多次采樣。所以我們可以不濾波，直接用隨機向量多次采樣來計算漫反射；也可以先濾波，然后進行一次采樣。下面采用先濾波，在進行一次采樣的方法。
6.2 預計算SH系數\(l_{i}\)
投影\(L(\mathbf{i})\)到SH空間：
6.1.2 濾波環境貼圖
float3 PrefilterEnvMap(float3 NormalDir)
{
float3 FilteredColor = 0;
float Weight = 0;
// 根據法線得到正交基
float3x3 TangentToWorld = GetTangentBasis(NormalDir);
const uint NumSamples = 4096;
for( uint i = 0; i < NumSamples; i++ )
{
// 計算二維隨機向量
float2 R = Hammersley(i, NumSamples);
// 通過球面映射轉換到直角坐標系得到隨機向量
float3 E = HemisphereSample_uniform(R.x, R.y);
// 把隨機向量轉換為世界空間
float3 H = mul(E, TangentToWorld);
// 得到cos項
float NoL = saturate(dot(NormalDir, H));
if( NoL > 0 )
{
FilteredColor += AmbientCubemap.Sample(H, 0).rgb * NoL;
Weight += NoL;
}
}
return FilteredColor / max( Weight, 0.001 );
}
可以用一個低分辨率的Cubemap來存儲，因為漫反射頻率非常低。
6.1.2 求出Spherical Harmonics
前3階基函數，可以用計算的結果替代常量計算部分：
float[9] SHcosineLobe(Vector3 normal)
{
float[] sh = new float[9];
float x = normal.x;
float y = normal.y;
float z = normal.z;
sh[0] = 1.0f / 2.0f * Sqrt(1.0f / PI);
sh[1] = Sqrt(3.0f / (4.0f * PI)) * z;
sh[2] = Sqrt(3.0f / (4.0f * PI)) * y;
sh[3] = Sqrt(3.0f / (4.0f * PI)) * x;
sh[4] = 1.0f / 2.0f * Sqrt(15.0f / PI) * x * z;
sh[5] = 1.0f / 2.0f * Sqrt(15.0f / PI) * z * y;
sh[6] = 1.0f / 4.0f * Sqrt(5.0f / PI) * (-x * x - z * z + 2 * y * y);
sh[7] = 1.0f / 2.0f * Sqrt(15.0f / PI) * y * x;
sh[8] = 1.0f / 4.0f * Sqrt(15.0f / PI) * (x * x - z * z);
return sh;
}
計算一次采樣的SH結果：
float3[9] GetSHSampleColor(float3 Dir)
{
float3[9] Result;
float3 color = AmbientCubemap.Sample(Dir, 0);
float[9] SHValue = SHcosineLobe(Dir);
for(uint i = 0; i < 9; ++i)
{
Result[i] = SHValue[i] * Color;
}
return Result;
}
下面需要生成很多隨機采樣方向，計算球諧向量，然后取平均。這里可以用CS的Barrier操作，只分一個線程組，等待組內所有采樣都記錄下之后，進行平均，類似Mipmap一樣一級一級兩兩平均：
// 采樣次數
const THREAD_COUNT = 256;
// 暫存顏色
float3[THREAD_COUNT] CollectedData;
float3[9] GetResult(uint id)
{
// 權重
Weight = 4.0f * PI;
// 計算一個采樣方向，類似上文中的半球空間
float3 SampleWay = SphereSample_uniform(Hammersley(id, THREAD_COUNT)).xyz;
// 采樣出球諧顏色
float3[9] Col = GetSHSampleColor(SampleWay);
[unroll]
for(uint i = 0; i < 9; ++i)
{
// 把第i項球諧顏色存儲
CollectedData[id] = Col[i] * Weight;
uint threadID = THREAD_COUNT;
while(threadID > 0)
{
threadID /= 2;
bool bInside = id < threadID;
float3 data = 0.0f;
// 等當前線程組所有線程都把第i項存到CollectedData中
GroupMemoryBarrierWithGroupSync();
if(isInside)
{
// 取均值
data = CollectedData[id * 2] + CollectedData[id * 2 + 1];
data *= 0.5f;
}
//等待所有可取均值的線程完成
GroupMemoryBarrierWithGroupSync();
if(isInside)
{
// 存儲均值
CollectedData[id] = data;
}
}
if(id == 0)
{
// 最終結果
Col[i] = CollectedData[0];
}
}
}
代碼中的Weight其實是蒙特卡洛方法中的\(\frac{1}{p\left(X_{i}\right)}\)
因為是球面均勻采樣所以PDF是\(\frac{1}{4\pi}\)。
現在我們得到了最終的9個SH Color。
6.3 運行時從SH空間重建光照
float3 GetDiffuseSH(float3[9] SHResult, float3 Normal)
{
float[9] Bi = SHcosineLobe(Normal);
float3 Col = 0.0f;
for(uint i = 0; i < 9; ++i)
{
Col += SHResult[i] * Bi[i];
}
return Col;
}
得到效果
7. 球諧的旋轉
PRT 的一個問題是如果 lighting 部分是預計算的，那就只適用於靜態環境光下的靜態物體渲染；環境光或者物體只要有變化，PRT 就不得不進行重新預計算；但得益於 SH 的旋轉不變性，我們至少可以讓 SH Lighting 適用於動態旋轉的情形而不必重新預計算。
環境光旋轉和物體旋轉在 PRT 渲染中是等價的，只是說看相對於哪個東西來看待旋轉而已。
7.1 SH系數旋轉的性質
- SH系數的旋轉是對SH系數的線性變換
- 對每一階SH系數的旋轉可以分別進行
第一條性質意味着給定SH系數\({t}=\left(t_{0}, \ldots, t_{k-1}\right)\)，我們希望將某個旋轉\(R\)作用到它上面，則存在\(k\)階方陣\(M_R\)使得\(M_R{t}\)就正好是我們所需要的新SH系數。
第二條性質意味着我們可以分別處理每一階的SH系數。比如我們使用了前3階SH來進行投影，那么就有\(l=0\) 對應的1個系數、 \(l=1\) 對應的3個系數以及\(l=2\) 對應的5個系數，加起來共9個SH系數。\(l=0\)的SH是個常量函數，不需要處理； \(l=1\)時的3個系數可以用一個3階方陣進行變換；\(l=2\) 時的5個系數則可以用一個5階方陣進行變換。
也就是說運行時根據光源旋轉解出每階(\({i}\))對應的\(M_{R}\)，就能在\(l_{i}\)不變的情況下得到\(L_i\)
7.2 方陣\(M_R\)的計算方法
以\(l=1\)為例（更高階可類推）：
設\(R\)是我們希望進行的旋轉操作的旋轉矩陣， \(P\)是求出任意三維方向向量所對應的5個SH值的投影函數，則對任意方向向量 \(x\)， \(M_R\)應滿足：
先用\(R\) 來變換\(x\)，再將變換結果投影到SH上，等價於先把 \(x\) 投影到SH上，之后再用 \(R\)去變換投影結果。
假設隨便帶入3個方向向量\({n}_{0},{n}_{1},{n}_{2}\)，這個等式都應該成立：
因為\({P}({n}_{i})\)和\({P}({R}{n}_{i})\)都是列向量，我們可以把這一坨式子寫成矩陣形式：
記\({A}=\left[{P}\left({n}_{0}\right), {P}\left({n}_{1}\right), {P}\left({n}_{2}\right)\right]\)如果\({A}\)是可逆的，那么我們可以很容易地把\(M_R\)找出來：
也就是說，只要\({n}_{0},{n}_{1},{n}_{2}\)選取的值使\({A}\)可逆，就可以解出\(M_R\)。
7.3 實現方法
7.3.1 預計算
預計算\({A}^{-1}\)，選取\(2 l+1\)個方向向量 \({n}_{0}, \ldots,{n}_{2l}\)使得 \({A}=\left[{P}\left({n}_{0}\right), \ldots, {P}\left({n}_{2 l}\right)\right]\) 可逆。求出\({A}^{-1}\)並存儲。
同時存儲向量 \({n}_{0},...,{n}_{2l}\)
7.3.2 運行時
- 計算主光源的旋轉矩陣\(R\)，將\(R\)應用於 \({n}_{0}, \ldots,{n}_{2l}\)，得到\({n}_{0}^{\prime}, \ldots, {n}_{2 l}^{\prime}\)
- 用\(P\)計算出\({n}_{0}^{\prime}, \ldots, {n}_{2 l}^{\prime}\)對應的SH投影值，得到\(2 l+1\)個\(2 l+1\)維列向量，將這些列向量構造成一個\(2 l+1\)階方陣\({S}\)
- 設待變換的SH系數用\(2 l+1\)維列向量 \(x\) 表示，計算\({S}{A}^{-1}x\) 即可。可以寫成\({S}({A}^{-1}x)\)來節省一些計算量。
這里就相當於，對於向量\({n}_{0},...,{n}_{2l}\)，已知\({A}^{-1}\)，這兩個在運行時已經變成了常量，只需要帶入變量\(R\)得到\(M_R\)即可。然后根據：
\({M}_{R} {P}({x})={P}({R} {x})\)，將\(M_R\)與\({P}({x})\)相乘即可得到旋轉后的\(L_i\)。
7.3.4 偽代碼
float3 GetDiffuseSH(float3[9] SHResult, float3 Normal, float3x3 RotationMatrix)
{
// 1. 運行時步驟1
float3 L1_n0 = mul(RotationMatrix, L1N0);
float3 L1_n1 = mul(RotationMatrix, L1N1);
float3 L1_n2 = mul(RotationMatrix, L1N2);
float3 L2_n0 = mul(RotationMatrix, L2N0);
float3 L2_n1 = mul(RotationMatrix, L2N1);
float3 L2_n2 = mul(RotationMatrix, L2N2);
float3 L2_n1 = mul(RotationMatrix, L2N3);
float3 L2_n2 = mul(RotationMatrix, L2N4);
// 2.運行時步驟2
float3 S_L1_n0 = SHcosineLobe_L1(L1_n0);
float3 S_L1_n1 = SHcosineLobe_L1(L1_n1);
float3 S_L1_n2 = SHcosineLobe_L1(L1_n2);
float[5] S_L2_n0 = SHcosineLobe_L2(L2_n0);
float[5] S_L2_n1 = SHcosineLobe_L2(L2_n1);
float[5] S_L2_n2 = SHcosineLobe_L2(L2_n2);
float[5] S_L2_n3 = SHcosineLobe_L2(L2_n3);
float[5] S_L2_n4 = SHcosineLobe_L2(L2_n4);
float[9] S1 = float3x3(
S_L1_n0.x, S_L1_n1.x, S_L1_n2.x,
S_L1_n0.y, S_L1_n1.y, S_L1_n2.y,
S_L1_n0.z, S_L1_n1.z, S_L1_n2.z,
);
float[25] S2 = float[25]{
S_L2_n0[0],S_L2_n1[0],S_L2_n2[0],S_L2_n3[0],S_L2_n4[0],
S_L2_n0[1],S_L2_n1[1],S_L2_n2[1],S_L2_n3[1],S_L2_n4[1],
S_L2_n0[2],S_L2_n1[2],S_L2_n2[2],S_L2_n3[2],S_L2_n4[2],
S_L2_n0[3],S_L2_n1[3],S_L2_n2[3],S_L2_n3[3],S_L2_n4[3],
S_L2_n0[4],S_L2_n1[4],S_L2_n2[4],S_L2_n3[4],S_L2_n4[5],
};
// 3.運行時步驟3
float[9] MR_L1 = mul(S1, A_Inv_L1);
float[25] MR_L2 = mul(S2, A_Inv_L2);
float[9] Bi = SHcosineLobe(Normal);
float L0 = Bi[0];
float3 L1 = float3(Bi[1], Bi[2], Bi[3]);
L1 = mul(MR_L1, L1);
float[5] L2 = float[5]{Bi[4], Bi[5], Bi[6], Bi[7], Bi[8]};
L2 = mul(MR_L2, L2);
return SHResult[0] * L0 +
SHResult[1] * L1[0] + SHResult[2] * L1[1] + SHResult[3] * L1[2] +
SHResult[4] * L2[0] + SHResult[5] * L2[1] + SHResult[6] * L2[2] + SHResult[7] * L2[4] + SHResult[8] * L2[5];
}
8. 游戲引擎中的實現
Unity中的Light Probe，Tetrahedral Tessellations方案。UE4的Indirect Light Cache和VLM均是基於球諧的，只不過在差值方式，存儲等方面有區別。這里先不寫了，文章有點長
引用
主要來源 GAMES202
截圖來自於GAMES 202 PPT和Unity引擎內
http://corysimon.github.io/articles/uniformdistn-on-sphere/
https://www.cnblogs.com/KillerAery/p/15335369.html
https://airguanz.github.io/2018/11/20/SH-PRT.html
https://zhuanlan.zhihu.com/p/51267461
https://zhuanlan.zhihu.com/p/103715075
https://zhuanlan.zhihu.com/p/49746076
https://zhuanlan.zhihu.com/p/362460950
https://zhuanlan.zhihu.com/p/373697912
https://zhuanlan.zhihu.com/p/149217557
https://zhuanlan.zhihu.com/p/351071035