基于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. 运行时步骤1float3 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