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