GAMES202笔记：基于SH的GI方案

基于Spherical Harmonics的实时全局照明
本文主要内容：频率与基函数、球面谐波函数、预计算辐射传输、球面均匀采样、构建正交基、SH实现、球面谐波函数的旋转及实现。

1. 频率与基函数

1.1 傅里叶级数展开

任何一个函数可以写成常数和一系列基函数(不同频率sin和cos项)的线性组合，基函数数量越多越接近于原函数的形状：

\[f(x)=\frac{A}{2}+\frac{2 A \cos (t \omega)}{\pi}-\frac{2 A \cos (3 t \omega)}{3 \pi}+\frac{2 A \cos (5 t \omega)}{5 \pi}-\frac{2 A \cos (7 t \omega)}{7 \pi}+\cdots \]

1.2 频率

在空间上图像信号数值的变化是否剧烈，任何一张图(也就是二维函数)的频率,也就是频域上对应的内容可以用一张频谱表示出来：

频谱最中心处是低频内容，我们可以做一个filtering(滤波)，从而去除一系列频率上的内容，我们对这张图用一个低通滤波器，从而把高频的内容去除掉：

卷积：卷积是一个模糊操作，在图上取任意一点，取点的周围一定区域内的像素值进行加权平均并将结果写回这个点，这就是卷积。
在空间域上做卷积等于在函数上做卷积；等于在频域上做原图频谱和卷积核频谱 的乘积操作：

对于任意的product integral(两个函数先乘积在积分)，我们将其认为是做了一个卷积操作。

\[\int_{\Omega} f(x) g(x) \mathrm{d} x \]

空间域上的两个信号\(f(x)\)和\(g(x)\)进行一个卷积，等于在频域上让两个信号相乘；如果两个信号有一个信号是低频的，那么频域上相乘后得到的结果也是低频的，最终相乘在积分的结果也是低频。的可以总结为：积分之后的频率取决于积分前最低的频

1.3 基函数

\[f(x)=\sum_{i} c_{i} \cdot B_{i}(x) \]

一个函数可以描绘成其他函数的线性组合，如\(f(x)\)可以描绘成一系列的\(B_i\)函数乘以各自对应的系数最终再相加在一起，这一系列的函数\(B_i\)就是基函数。

2. 球谐函数

SH(Spherical Harmonics)是一系列基函数，系列中的每个函数都是二维函数，并且每个二维函数都是定义在球面上的，阶数越高，频率越高。

下图是对SH的可视化，与一维的傅里叶一样，SH也存在不同频率的函数，但不同频率的函数个数也不同，频率越高所含有的基函数越多：

\(l\)表示的是阶数，通常第\(l\)阶有\(2l+1\)个基函数，前\(n\)阶有\(n^2\)个基函数。
越偏白的蓝色地方值越大，越黑的地方值越小；而黄色中则表示偏白的地方表示其绝对值大，偏黑的地方表示绝对值小；蓝色表示正，黄色表示负。

2.1 投影

由于一个函数\(f(\omega)\)可以由一系列基函数和系数的线性组合表示，那么怎么确定基函数前面的系数，这就需要通过投影操作：

\[c_{i}=\int_{\Omega} f(\omega) B_{i}(\omega) \mathrm{d} \omega \]

我们知道函数\(f(\omega)\),通过对应的基函数\(B_i(\omega)\)进行投影操作，从而求出各基函数对应的系数\(C_i\)，与以下操作是同一个道理：在空间中想描述一个向量，可以xyz三个坐标来表达，把xyz轴当做三个基函数，把向量投影到xyz轴上，得到三个系数就是三个坐标。

2.2 重建

知道基函数对应的系数，就能用系数和基函数恢复原来的函数。
由于基函数的阶可以是无限个的，越高的阶可恢复的细节就越好，但一方面是因为更多的系数会带来更大的存储压力、计算压力，而描述变化比较平滑的环境漫反射部分，用3阶SH就足够。

\[c_{i}=\int_{\Omega} f(\omega) B_{i}(\omega) \mathrm{d} \omega \]

对于重建diffuse来说这里的\(f(\omega)\)对应环境光照。

2.3 球谐函数的正交性

\[\begin{aligned} &\int_{\Omega} B_{i}(\mathbf{i}) \cdot B_{j}(\mathbf{i}) \mathrm{di}=\mathbf{1} \quad(\mathbf{i}=\boldsymbol{j}) \\ &\int_{\Omega} B_{i}(\mathbf{i}) \cdot B_{j}(\mathbf{i}) \mathrm{di}=\mathbf{0} \quad(\mathbf{i} \neq j) \end{aligned} \]

旋转一个基函数之后，得到的函数就不再是一个基函数(因为基函数有严格的朝向等限制)，但是旋转球谐函数等价于同阶基函数的线性组合。这意味着可以简单地进行投影、重建、旋转。

2.4 内容补充

diffuse BRDF类似于一个低通滤波器，使用一些低频信息就可以恢复出原始内容。因为积分之后的频率取决于积分前最低的频率，当diffuse BRDF使用低频信息即可恢复内容时，也就意味着无论光照项是多么复杂，其本应该用多高频的基函数去表示，但我们希望得到的是其与BRDF之积的积分，所以可以使用比较低频的基函数去描述灯光。
前3阶球谐基函数：

\[\begin{aligned} &Y_{0}^{0}(\theta, \varphi)=\frac{1}{2} \sqrt{\frac{1}{\pi}} \\ &Y_{1}^{-1}(\theta, \varphi)=\frac{1}{2} \sqrt{\frac{3}{2 \pi}} \sin \theta e^{-i \varphi} \\ &Y_{1}^{0}(\theta, \varphi)=\frac{1}{2} \sqrt{\frac{3}{\pi}} \cos \theta \\ &Y_{1}^{1}(\theta, \varphi)=\frac{-1}{2} \sqrt{\frac{3}{2 \pi}} \sin \theta e^{i \varphi} \\ &Y_{2}^{-2}(\theta, \varphi)=\frac{1}{4} \sqrt{\frac{15}{2 \pi}} \sin ^{2} \theta e^{-2 i \varphi} \\ &Y_{2}^{-1}(\theta, \varphi)=\frac{1}{2} \sqrt{\frac{15}{2 \pi}} \sin \theta \cos \theta e^{-i \varphi} \\ &Y_{2}^{0}(\theta, \varphi)=\frac{1}{4} \sqrt{\frac{5}{\pi}}\left(3 \cos ^{2} \theta-1\right) \\ &Y_{2}^{1}(\theta, \varphi)=\frac{-1}{2} \sqrt{\frac{15}{2 \pi}} \sin \theta \cos \theta e^{i \varphi} \\ &Y_{2}^{2}(\theta, \varphi)=\frac{1}{4} \sqrt{\frac{15}{2 \pi}} \sin ^{2} \theta e^{2 i \varphi} \end{aligned} \]

3. 预计算辐射传输

PRT(Precomputed Radiance Transfer)思想：将光照切割成Lighting, Light Transport两部分：

\[L(\mathbf{0})=\int_{\Omega} \underbrace{L(\mathbf{i})}_{\text {lighting }} \underbrace{V(\mathbf{i}) \rho(\mathbf{i}, \mathbf{0}) \max (0, \boldsymbol{n} \cdot \mathbf{i})}_{\text {light transport }} \mathrm{di} \]

首先，因为在Diffuse的情况下，BRDF几乎是一个常数，可以把它提到外面：

\[L(\mathbf{0})=\rho \int_{\Omega} L(\mathbf{i}) V(\mathbf{i}) \max (0, \boldsymbol{n} \cdot \mathbf{i}) \mathrm{d} \mathbf{i} \]

而Lighting项可以写成基函数的形式：

\[L(\mathbf{i}) \approx \sum l_{i} B_{i}(\mathbf{i}) \]

同样的，Light transport项也可以写成基函数求和的形式：

\[T\left(\mathbf{i}\right) \approx \sum_{} l_{q} B_{q}\left(\mathbf{i}\right) \]

\(l_i\)、\(l_q\)是常数，可以提出来：

\[L(\mathbf{0})=\sum_{i} \sum_{q} l_{i} l_{q} \int_{\Omega^{+}} B_{i}\left(\mathbf{i}\right) B_{q}\left(\mathbf{i}\right) \mathrm{d} \mathbf{i} \]

因为球谐的正交性，只有\(i = q\)的时候，\(B_{i}\left(\mathbf{i}\right) B_{q}\left(\mathbf{i}\right)\)相乘才有值，所以这个函数的复杂度仍是O(n)。
下面将基函数形式的Lighting项带入渲染方程：

\[\begin{aligned} &{L(\mathbf{0}) \approx \rho \sum l_{i} \int_{\Omega} \mathrm{B}_{i}(\mathbf{i}) V(\mathbf{i}) \max (0, \boldsymbol{n} \cdot \mathbf{i}) \mathrm{di}} \\ &L(\mathbf{0}) \approx \rho \sum l_{i} T_{i} \end{aligned} \]

对于积分中的部分来说，相当于light transport乘与一个基函数，这就成了lighting transport投影到一个基函数的系数。因为对light transport做了预计算，visibility成为了常量，只能对静止物体进行计算。
一个SH基函数旋转后都可以被同阶的SH基函数线性组合得到，因此可以解决光源旋转问题。
现在将lighting这个球面函数，通过SH的基函数用一堆系数来表示，这些系数排成一行也就是组成了向量，因此光照变成了一个向量：

投影到SH空间：

\[l_{i}=\int_{\Omega} L(\mathbf{i}) \cdot B_{i}(\mathbf{i}) \mathrm{di} \]

从SH空间重建：

\[L(\mathbf{i}) \approx \sum l_{i} B_{i}(\mathbf{i}) \]

4. 球面均匀采样

4.1 Hammersley采样

Hammersley这是一种均匀分布的2D随机采样，他将十进制转换成二进制，再将二进制转换到[0,1]之间的小数，这一过程被称作 Radical Inverse。具体表示方法见下图：

Hammersley采样点集合：

\[p_{i}=\left(x_{i}, y_{i}\right)=(i / N, \phi(i)) \]

\(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均匀采样，会出现两极密度高，中间密度低的情况：

通过求球面概率密度分布函数的反函数，可以对球面均匀采样的映射函数：

\[\left\{\begin{array}{l} \theta=\arccos \left(1-2 \xi_{x}\right) \\ \phi=2 \pi \xi_{y} \end{array}\right.,\xi_{x}, \xi_{y} \sim \text { Uniform }[0,1] \]

现在，在2D平面坐标\((\xi_{x},\xi_{y})\)做Hammersley均匀采样，然后再映射回球坐标系即可得到均匀的分布：

PS：我们可以很容易的得到 \(cos\theta=1-2 \xi_{x}\)

4.3 球坐标系与直角坐标系的转换

直角坐标系\((x,y,z)\)转球坐标系\((r, \theta, \varphi)\)：

\[\begin{gathered} r=\sqrt{x^{2}+y^{2}+z^{2}} \\ \theta=\arccos \frac{z}{r} \\ \varphi=\arctan \frac{y}{x} \end{gathered} \]

球坐标系\((r, \theta, \varphi)\)转直角坐标系\((x,y,z)\)：

\[x=r \sin \theta \cos \varphi\\ y=r \sin \theta \sin \varphi\\ z=r \cos \theta \]

现在我们可以结合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空间：

\[l_{i}=\int_{\Omega} L(\mathbf{i}) \cdot B_{i}(\mathbf{i}) \mathrm{di} \]

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)}\)

\[\int f(x) \mathrm{d} x=\frac{1}{N} \sum_{i=1}^{N} \frac{f\left(X_{i}\right)}{p\left(X_{i}\right)} \]

因为是球面均匀采样所以PDF是\(\frac{1}{4\pi}\)。
现在我们得到了最终的9个SH Color。

6.3 运行时从SH空间重建光照

\[L(\mathbf{i}) \approx \sum l_{i} B_{i}(\mathbf{i}) \]

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\)应满足：

\[{M}_{R} {P}({x})={P}({R} {x}) \]

先用\(R\) 来变换\(x\)，再将变换结果投影到SH上，等价于先把 \(x\) 投影到SH上，之后再用 \(R\)去变换投影结果。
假设随便带入3个方向向量\({n}_{0},{n}_{1},{n}_{2}\)，这个等式都应该成立：

\[\left\{\begin{array}{l} {M}_{R} {P}\left({n}_{0}\right)={P}\left({R} {n}_{0}\right) \\ {M}_{R} {P}\left({n}_{1}\right)={P}\left({R} {n}_{1}\right) \\ {M}_{R} {P}\left({n}_{2}\right)={P}\left({R} {n}_{2}\right) \\ \end{array}\right. \]

因为\({P}({n}_{i})\)和\({P}({R}{n}_{i})\)都是列向量，我们可以把这一坨式子写成矩阵形式：

\[{M}_{R}\left[{P}\left({n}_{0}\right), {P}\left({n}_{1}\right), {P}\left({n}_{2}\right)\right]=\left[{P}\left({R} {n}_{0}\right), {P}\left({R} {n}_{1}\right), {P}\left({R} {n}_{2}\right)\right] \]

记\({A}=\left[{P}\left({n}_{0}\right), {P}\left({n}_{1}\right), {P}\left({n}_{2}\right)\right]\)如果\({A}\)是可逆的，那么我们可以很容易地把\(M_R\)找出来：

\[{M}_{R}=\left[{P}\left({R} {n}_{0}\right), {P}\left({R} {n}_{1}\right), {P}\left({R} {n}_{2}\right)\right] {A}^{-1} \]

也就是说，只要\({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 运行时

1. 计算主光源的旋转矩阵\(R\)，将\(R\)应用于 \({n}_{0}, \ldots,{n}_{2l}\)，得到\({n}_{0}^{\prime}, \ldots, {n}_{2 l}^{\prime}\)
2. 用\(P\)计算出\({n}_{0}^{\prime}, \ldots, {n}_{2 l}^{\prime}\)对应的SH投影值，得到\(2 l+1\)个\(2 l+1\)维列向量，将这些列向量构造成一个\(2 l+1\)阶方阵\({S}\)
3. 设待变换的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