(∩^o^)⊃━☆ﾟ.*･｡

0x00 To begin with

这篇文章将会包含以下内容：

  • 部分课程内容回顾
  • GAMES202 作业 2

For reference👇：

master 分支上是 202 全部作业汇总，working 分支用于提交代码。需要原始作业可下载 master 分支的包

0x01 前置准备 Warm Up

渲染方程 Render Equation

James T. Kajiya（没错就是那个提出 Kajiya-Kay 头发着色模型的 Kajiya）在 1986 年提出了渲染方程（The Render Equation），原论文里长这样：

$$ I(x,x')=g(x,x')\left[ \epsilon(x,x')+\int_{S}\rho(x,x',x'')I(x',x'')dx'' \right] $$

当然现在看得更舒服的形式长这样：

$$ L_{o}(p,\omega_{o})=L_{e}(p,\omega_{o})+\int_{\Omega^+}L_{i}(p,\omega_{i})f_{r}(p,\omega_{i},\omega_{o})\cos\theta_{i}\space d\omega_{i} $$

其中：

  • $L_{o}(p,\omega_{o})$：要求出射方向的辐射率（Radience）
  • $L_{e}(p,\omega_{o})$：自发光项
  • $\Omega^+$：积分域，仅对表面正半球积分
  • $L_{i}(p,\omega_{i})$：Lighting，某一入射方向的辐射率
  • $f_{r}(p,\omega_{i},\omega_{o})$：BRDF，描述某一入射方向的光线在出射方向的贡献
  • $\cos \theta_{i}$：表面法向与某一入射方向形成的夹角余弦

所以渲染方程其实就描述了在当前 Shading Point 上，给定入射光与出射方向（观察方向），在出射方向上能够接收到的 Radiance，任何在表面上进行的光线传播的 Shading 都需要解这个方程。

Split Sum

对于环境光，通常使用 CubeMap 或者等距柱形投影将图像投射到无限远的球面上，要获得环境光的 Shading（目前暂不考虑 Shadowing），就需要求解渲染方程，核心是对于 Lighting 部分的求解。

先抛开渲染方程，想象一下当表面是理想镜面的情况：对于某个 Shading Point，其反射行为是镜面反射，那么我们可以很容易地求出射方向对应的入射方向，根据这个方向查询环境光的颜色。

更 General 一点，当表面比较粗糙的情况，能够反射到出射方向的入射角在镜反射方向上会在某个角度 $\phi$ 范围内抖动分布，最终反射到出射方向的 Lighting 受到这个范围内所有采样点颜色的影响。这个过程可以类比于对原始的环境贴图进行滤波，且粗糙度越大，滤波核越大，这样就可以仍使用镜反射方向查询滤波后的图像。这些滤波后的图像当然可以预计算，形成一个模糊层级，然后渲染时根据不同粗糙度使用不同图像。

但是目前 Lighting 项在积分内部，也就是说我们要对整个半球面上的入射方向都进行一次查询，还是比较 Costy 的，有没有办法把 Lighting 从积分里独立出来？下面有一种近似方法：

$$ \int_{\Omega}f(x)g(x)dx\approx \frac{\int_{\Omega_{G}}f(x)dx}{\int_{\Omega_{G}}dx} \cdot \int_{\Omega}g(x)dx $$

这个近似在 $g(x)$ 贡献比较小或者比较平滑（低频）时是相对准确的。漫反射 BRDF 就很好地满足这个性质，所以可以对渲染方程应用这个近似：

$$ L_{o}(p,\omega_{o})\approx\frac{\int_{\Omega_{f_{r}}}L_{i}(p,\omega_{i})d\omega_{i}}{\int_{\Omega_{f_{r}}}d\omega_{i}}\cdot \int_{\Omega^+}f_{r}(p,\omega_{i},\omega_{o})\cos \theta_{i} \space d\omega_{i} $$

乘号左边的部分代表 Lighting 在 BRDF 范围内的平均 Radience，可以认为是对出射方向有贡献的入射方向范围内在环境贴图上采样颜色的平均值。通过对环境贴图的预滤波，可以通过一次查询得到这个值。

球谐函数 Spherical Harmonics

球谐函数

维基百科上球谐函数是下面这个式子在 $m,l$ 不同取值时的解：

$$ Y^m_{l}(\theta, \phi)=(-1)^m\sqrt{ \frac{2l+1}{4\pi}\frac{(l-|m|)!}{(l+|m|)!} }P^m_{l}(\cos \theta)e^{im\phi} $$

实数形式下可视化出来是长这样的：

看上去是不是很像物理/量子化学里面的电子云啥的（这俩还真有关系），但在图形学里面，我们不需要理解其物理意义，可以简单将其理解为一组定义在球面上的两两正交的基函数，有了它我们可以使用傅里叶变换的思想去拟合任何的球面函数。

基函数（Basis Function）：一组两两之间相互独立的函数，即在彼此之上的投影为 0，可以类比线性代数的基向量。一组基函数张成的空间是其经过线性组合能够表示的所有函数的集合

$$ \mathcal{F} = \left\{ f(x) = \sum_{i=1}^n c_{i}B_{i}(x) \,\middle|\, c_{i} \in \mathbb{R} \right\} $$

可以通过将原函数投影到某个基函数上获取该基函数前的系数 $c_{i}$：

$$ c_{i}=\int_{\Omega}f(x)B_{i}(x)\,dx $$

实际使用时需要把积分离散化为一系列求和的形式，给定一个函数 $f(\omega_{i})$，可以通过投影求得其在 $Y^m_{l}(\omega_{i})$ 上的系数：

$$ c^m_{l}\approx\sum_{i=1}^Nf(\omega_{i})\cdot Y^m_{l}(\omega_{i})\cdot \omega_{i} $$

然后就可以用类似泰勒展开的样子用基函数表示原函数：

$$ \begin{aligned} f(x)=&\sum^{\infty}_{l=0}\sum^l_{m=-l}c^m_{l}\cdot Y^m_{l}(\omega_{i}) \\ =& \, c^0_{0}Y^0_{0}(\omega_{i})+c^1_{-1}Y^1_{-1}(\omega_{i})+c^1_{0}Y^1_{0}(\omega_{i})+\,\dots \end{aligned} $$

对于低频的函数，只取前几项进行近似，就能达到 Adequate 的效果。如果 BRDF 是粗糙的，就相当于一个低通滤波器，其和 Lighting 的乘积的积分仍然是低频的。

PRT（Precomputerd Radience Transfer）

球谐函数用一组基函数拟合渲染方程中的 Lighting 项，解决了 Shading 的问题，如果需要 Shadowing，就需要加入可见性项 $V(p,\omega_{i})$，渲染方程变为：

$$ L_{o}(p,\omega_{o})=\int_{\Omega} \underbrace {L_{i}(p,\omega_{i})}_{Lighting} \, \underbrace {f_{r}(p,\omega_{i},\omega_{o})}_{BRDF} \, \underbrace {V(p,\omega_{i})}_{Visibility} \, \max(0,\cos \theta_{i})\,d\omega_{i} $$

暴力算法就是朝各个方向分别对 Lighting，BRDF 和 Visibility 计算，然后再逐项相乘，开销太大。Sloan 在 2002 年提出了PRT（Precomputed Radience Transfer），将积分内部分为 Lighting 和 Light Transport（光线传输）两个部分：

$$ L_{o}(p,\omega_{o})=\int_{\Omega} \underbrace {L_{i}(p,\omega_{i})}_{Lighting} \, \underbrace {f_{r}(p,\omega_{i},\omega_{o}) \, V(p,\omega_{i}) \, \max(0,\cos \theta_{i})}_{Light \,Transport} \, d\omega_{i} $$

假设 BRDF 项是 Diffuse 的，BRDF 项是一个定值 $\frac{\rho}{\pi}$，将 Lighting 表示为系数和基函数的形式，然后交换积分与求和的位置，就可以将积分内部再次表示为在基函数上的投影：

$$ \begin{aligned} L(p,\omega_{o})=& \frac{\rho}{\pi} \int_{\Omega}\sum l_{i}B_{i}(\omega_{i})V(p,\omega_{i})\max(0,\cos \theta_{i}) \, d\omega_{i} \\ \approx& \frac{\rho}{\pi} \sum l_{i} \underbrace{ \int_{\Omega}B_{i}(\omega_{i})V(p,\omega_{i})\max(0,\cos \theta_{i}) }_{Precomputed} \, d\omega_{i} \\ \approx& \frac{\rho}{\pi} \sum l_{i}T_{i} \end{aligned} $$

注意这里是对 $\omega_{i}$ 进行积分，对每一个基函数上的投影进行求和

其中 $l_{i}$ 是 Lighting 项在基函数上投影的一组系数，$T_{i}$ 是 Light Transport 项在基函数上投影的系数。这样就把多个函数乘积的积分转变为一个简单的点乘。

然而这些预计算需要一些前提：

  • Visibility 是固定的，即要求场景中的物体不能动
  • Lighting 是“固定”的，Lighting 的成分不能动态变化，但 Lighting 可以旋转（相当于对基函数张成的线性空间进行旋转）

很合理，就像烘焙的光照信息是静态的一样。

0x02 预计算球谐系数 Precompute SH Coef

光照项系数 Lighting

要算 SH 系数，就要把环境光照投影到球谐函数上，球谐函数可以调用 sh::EvalSH 获取，Lighting 项可以通过预存的 images 数组获取，剩下一个立体角 $d\omega_{i}$：

$$ d\omega_{i}=\sin \theta \, d\theta \, d\phi $$

其中 $\theta$ 是极角，$\phi$ 是方位角。对于每一个环境贴图上的像素，需要计算其投影到单位球面的面积，框架给了 CalcArea 这个函数计算这个面积。

先把 texel 从 $[0,res-1]$ 映射到以 CubeMap 当前面中心为原点的 $\left[ \frac{1}{res}-1,1-\frac{1}{res} \right]$ 坐标范围内。然后通过四个角点的差值计算面积，每个点所对应的面积可以使用一个公式求解：

$$ A(x,y)=\arctan\left( \frac{xy}{\sqrt{ x^2+y^2+1 }} \right) $$

所以在 PrecomputeCubemapSH 中，只需要为每个 SH 基函数上累加每个 texel 的贡献，就可以得到 Lighting 项的球谐系数。

无阴影传输项系数 UnShadowed

传输项球谐系数的计算在 preprocess 函数中，assignmet sheet 里给的计算传输项伪代码比较长，但框架里主要关注到 shFunc 这个 lambda 表达式，该表达式返回需要投影的函数在给定方向上的函数值，后续会传入 sh::ProjectFunction 进行计算。

对于无阴影的情况，要投影的函数就是简单的 cos 项：

$$ M_{DU}=\max(N_{x}\cdot \omega_{i},0) $$

这里 H 后来我放 if 外面了，免得下面再写一遍

有阴影传输项系数 Shadowed

$$ M_{DU}=V(\omega_{i})\max(N_{x}\cdot \omega_{i},0) $$

要投影的函数再乘上一个 Visibility 项即可：

如果想看这步的结果，就把 prt.xml 里改成 shadowed，然后运行下 nori.exe。

1<string name="type" value="shadowed" />

记得改完代码 build 一下，不然 .exe 还是上一次的结果

反射传输项系数 Inter-reflection

对于具有相互反射的传输项，其方程变为：

$$ L_{DI}=L_{DS}+\frac{\rho}{\pi}\int_{S}\hat{L}(x',\omega_{i})(1-V(\omega_{i}))\max(N_{x}\cdot \omega_{i},0)\,d\omega_{i} $$

我们已经计算了每个顶点的直接可见光源的 Lightning（$L_{DS}$），还需要计算间接光的 Lighting，其中：

  • $\hat{L}(x',\omega_{i})$：从 $\omega_{i}$ 到点 $x'$ 的 Radiance
  • $1-V(\omega_{i})$：被遮挡的部分

所以这个式子对每个顶点计算来自直接可见光源的 Lighting，然后遍历所有可能的 $\omega_{i}$ 计算间接光的 Lighting，对于每个交点返回重心坐标插值后的球谐系数，然后以当前交点为原点，迭代计算传输项系数。

因为觉得 lambda 有点好玩，就拿 lambda 写了一个递归，其实也可以另外封装一个函数，图片太长不好截，直接放代码：

 1 if (m_Type == Type::Interreflection)
 2        {
 3            // for each vertex cast a ray to the scene
 4            // if occluded return interpolated sh coef
 5            // then recursively call this function
 6            for(int i = 0;i < mesh->getVertexCount(); i++)
 7            {
 8                const Point3f &v = mesh->getVertexPositions().col(i);
 9                const Normal3f &n = mesh->getVertexNormals().col(i);
10
11                // iterate lambda
12                // need to pass pos, normal, scene, original coeff and current bounce
13                std::function<std::unique_ptr<std::vector<double>>(Eigen::MatrixXf*, const Point3f&, const Normal3f&, const Scene*, int)> iterateSHFunc;
14                iterateSHFunc = [&](Eigen::MatrixXf* TransportSHCoeffs, const Point3f &pos, const Normal3f &normal, const Scene* s, int bounce) -> std::unique_ptr<std::vector<double>>
15                {
16                    // Step 1: allocate and assign 0 to coeffs
17                    std::unique_ptr<std::vector<double>> coeffs(new std::vector<double>());
18                    coeffs->assign(SHCoeffLength, 0.0);
19                
20                    // Step 2: judge by bounce to terminate
21                    if (bounce > m_Bounce)
22                        return coeffs;
23                
24                    // Step 3: generate sample_side^2 uniformly and stratified samples over the sphere
25                    const int sample_side = static_cast<int>(floor(sqrt(m_SampleCount)));
26                    std::random_device rd;
27                    std::mt19937 gen(rd());
28                    std::uniform_real_distribution<> rng(0.0, 1.0);
29                    for (int t = 0; t < sample_side; t++)
30                    {
31                        for (int p = 0; p < sample_side; p++)
32                        {
33                            double alpha = (t + rng(gen)) / sample_side;
34                            double beta = (p + rng(gen)) / sample_side;
35                            double phi = 2.0 * M_PI * beta;
36                            double theta = acos(2.0 * alpha - 1.0);
37                
38                            // Step 4: for each wi cast a ray to the scene
39                            Eigen::Array3d d = sh::ToVector(phi, theta);
40                            const auto wi = Vector3f(d.x(), d.y(), d.z());
41                            float H = wi.dot(normal); 
42                
43                            Ray3f ray = Ray3f(pos, wi);
44                            Intersection its;
45                            if(H > 0.0 && s->rayIntersect(ray, its))
46                            {
47                                // Step 5: if occluded, get the intersected point and its normal
48                                Point3f intersectPos = its.p;
49                                Vector3f baryCoord = its.bary;
50                                Point3f idx = its.tri_index;
51                                MatrixXf normals = mesh->getVertexNormals();
52                                Normal3f intersectNormal = Normal3f(normals.col(idx.x()).normalized() * baryCoord.x() +
53                                                                        normals.col(idx.y()).normalized() * baryCoord.y() +
54                                                                        normals.col(idx.z()).normalized() * baryCoord.z()).normalized();
55                                auto nextBounceCoeffs = iterateSHFunc(TransportSHCoeffs, intersectPos, intersectNormal, scene, bounce + 1);
56                
57                                // Step 6: accumulate the coeffs
58                                for (int k = 0; k < SHCoeffLength; k++)
59                                {
60                                    auto interpolateSH = (TransportSHCoeffs->col(idx.x()).coeffRef(k) * baryCoord.x() +
61                                                                TransportSHCoeffs->col(idx.y()).coeffRef(k) * baryCoord.y() +
62                                                                TransportSHCoeffs->col(idx.z()).coeffRef(k) * baryCoord.z());
63                
64                                    (*coeffs)[k] += (interpolateSH + (*nextBounceCoeffs)[k]) * H;
65                                }
66                            }
67                        }
68                    }
69                    // Step 7: return the coeffs
70                    double weight = (sample_side * sample_side);
71                    for (unsigned int c = 0; c < coeffs->size(); c++) {
72                        (*coeffs)[c] /= weight;
73                    }
74                
75                    return coeffs;
76                };
77                
78                // Step 8: call the iterateSHFunc
79                auto interreflectCoeffs = iterateSHFunc(&m_TransportSHCoeffs, v, n, scene, 1);
80                for (int j = 0; j < SHCoeffLength; j++)
81                {
82                    m_TransportSHCoeffs.col(i).coeffRef(j) += (*interreflectCoeffs)[j];
83                }
84            }
85        }

这里使用 std::function 包装了一下这个 lambda，不然其内部不能调用自己；生成采样点的过程可以照着 sh::ProjectFunction 里面抄

在调用时可以修改 m_Bounce 的值控制光线反弹次数（lambda 迭代次数），我这里仍然是默认值 1，之前改了一次 2 给我跑了一个多小时，(๑╹◡╹๑)。

0x03 实时球谐光照 Realtime SH Lighting

运行 nori.exe 后，计算得到的球谐系数将保存在 prt/scenes/cubemap/ 对应环境贴图文件夹下，需要手工先把 light.txt 和 transport.txt 放到 homework2/assets/cubemap/ 对应的环境贴图文件夹下。

engine.js 中的第 88-114 行已经帮我们完成了解析这些文件的工作，最终 Lighting 和 Transport 项的球谐系数被保存在 precomputeL 和 precomputeLT 两个数组中，这两个数组包含在 envmap 中所有环境纹理的球谐系数，其形状为：

$$ \begin{aligned} precomputeL =& EnvMapNum*3*(SHOrder+1)^{2} \\ precomputeLT =& EnvMapNum*VertexCount*(SHOrder+1)^{2} \end{aligned} $$

可以通过 guiParams.envmapId 获取当前 GUI 中激活的环境纹理索引，以从以上数组中获取对应的球谐系数。

新建 PRT 材质 New PRT Material

在 src/material/ 下新建PRTMaterial.js，编写一个最简单的材质，传的统一变量仅为预计算光照项的各分量（RGB），attribute（aPrecomputeLT）按照 assignment sheet 里的要求传入：

 1class PRTMaterial extends Material {
 2
 3    constructor(vertexShader, fragmentShader) {
 4        super({
 5            // PRT
 6            'uPrecomputeL[0]': { type: 'precomputeL', value: null },
 7            'uPrecomputeL[1]': { type: 'precomputeL', value: null },
 8            'uPrecomputeL[2]': { type: 'precomputeL', value: null },
 9        }, ['aPrecomputeLT'], vertexShader, fragmentShader, null);
10    }
11}
12
13async function buildPRTMaterial(vertexPath, fragmentPath) {
14
15
16    let vertexShader = await getShaderString(vertexPath);
17    let fragmentShader = await getShaderString(fragmentPath);
18
19    return new PRTMaterial(vertexShader, fragmentShader);
20
21}

在 Shader 中，uPrecomputeL 是一个类型为 mat3 的数组，索引对应 RGB 分量的球谐系数，需要逐一绑定

utils/tools.js 中提供了用于将 precomputeL 分解为三个 mat3 类型的数组的函数 getMat3ValueFromRGB，每帧的 render 中需要对其进行绑定，WebGLRenderer.js 中：

1let precomputeL_RGBMat3 = getMat3ValueFromRGB(precomputeL[guiParams.envmapId]);
2for(let j = 0; j < 3 ; j++){
3	if(k == `uPrecomputeL[${j}]`){
4		gl.uniformMatrix3fv(this.meshes[i].shader.program.uniforms[k],false,precomputeL_RGBMat3[j]);
5	}
6}

之后可以先在 index.html 中引入：

1<script src="src/materials/PRTMaterial.js" defer></script>

然后在 loadOBJ.js 中增加 PRT 材质的 case：

1case 'PRTMaterial':
2	material = buildPRTMaterial("./src/shaders/PRTShader/PRTVertex.glsl", "./src/shaders/PRTShader/PRTFragment.glsl");
3	break;

编写着色器 Write Shader

在 shaders/PRTShader/ 下新建 PRT 的顶点和片元着色器文件，作业中需要我们在顶点着色器中完成对 Lighting 和 Transport 的球谐系数的点乘，结果 vColor 经插值后传至片元着色器，直接输出。

所以片元着色器很简单：

 1// prtFragment.glsl
 2#ifdef GL_ES
 3precision mediump float;
 4#endif
 5
 6varying highp vec3 vColor;
 7
 8void main(void) {
 9   gl_FragColor = vec4(vColor, 1.0);
10}

顶点着色器中需要完成对 vColor 的计算，看一下传过来的 uPrecomputeL[3] 和 aPrecomputeLT 长什么样子：

emm，都是 mat3 没办法简单点乘，那就写一个函数，把它们都分成三个 vec3，然后点乘求和即可。

 1// prtVertex.glsl
 2attribute vec3 aVertexPosition;
 3attribute vec3 aNormalPosition;
 4attribute mat3 aPrecomputeLT;
 5
 6uniform mat4 uModelMatrix;
 7uniform mat4 uViewMatrix;
 8uniform mat4 uProjectionMatrix;
 9uniform mat3 uPrecomputeL[3];
10
11varying highp vec3 vNormal;
12varying highp vec3 vColor;
13
14float LoT (mat3 L, mat3 LT) {
15  vec3 L0 = L[0];
16  vec3 L1 = L[1];
17  vec3 L2 = L[2];
18  vec3 LT0 = LT[0];
19  vec3 LT1 = LT[1];
20  vec3 LT2 = LT[2];
21  return dot(L0, LT0) + dot(L1, LT1) + dot(L2, LT2);
22}
23
24void main(void) {
25
26  vNormal = (uModelMatrix * vec4(aNormalPosition, 0.0)).xyz;
27
28  gl_Position = uProjectionMatrix * uViewMatrix * uModelMatrix *
29                vec4(aVertexPosition, 1.0);
30
31  // rgb color from precomputed light transport
32  float r = LoT(uPrecomputeL[0],aPrecomputeLT);
33  float g = LoT(uPrecomputeL[1],aPrecomputeLT);
34  float b = LoT(uPrecomputeL[2],aPrecomputeLT);
35    vColor = vec3(r, g, b);
36}

结果 Outcome

最后在 engine.js 中调用 loadOBJ 加载 202 娘，顺便把还有一张没在 GUI 上的环境纹理放上去：

 1var envmap = [
 2	'assets/cubemap/GraceCathedral',
 3	'assets/cubemap/Indoor',
 4	'assets/cubemap/Skybox',
 5	// Begin TOP changes Add envmap CornellBox
 6	'assets/cubemap/CornellBox'
 7	// End TOP changes
 8];
 9// ...
10
11	loadOBJ(renderer, 'assets/mary/', 'mary', 'PRTMaterial', boxTransform);
12// ...
13
14function createGUI() {
15	const gui = new dat.gui.GUI();
16	const panelModel = gui.addFolder('Switch Environemtn Map');
17	// Begin TOP changes Add envmap CornellBox case in UI
18	panelModel.add(guiParams, 'envmapId', { 'GraceGathedral': 0, 'Indoor': 1, 'Skybox': 2, 'CornellBox':3}).name('Envmap Name');
19	// End TOP changes
20	panelModel.open();
21}

Inter-reflection m_Bounce = 1

0x04 球谐旋转 SH Rotation

球谐函数具有两个很好的性质：

  • 旋转不变性：旋转球谐展开后的函数，相当于旋转这些球谐函数
  • 独立的线性旋转：对于每层球谐函数上的系数，可以分别进行旋转，并且该旋转是线性的

把球谐函数和基向量类比就比较清楚了：旋转一个三维空间中的向量，相当于旋转坐标基，在这组新的坐标基下，其系数和旋转前其在原坐标系的系数相等。那么要求旋转后在原坐标系下的系数，只需要将新的坐标基使用原来的坐标基表示，然后重新整理系数即可。

对球谐每阶 $l$，有 $2l+1$ 维的球谐系数组成的向量，假设存在 $2l + 1$ 个任意的三维法向 $n$，旋转函数 $R$，能够返回该阶球谐投影系数的函数 $P$，需要求得旋转整个球面函数的矩阵 $M$，有：

$$ MP(n_{i})=P(R(n_{i}))\,,\,i \in[-l,l] $$

LHS：对整个球面函数进行旋转；RHS：先旋转方向 $n_i$，再代入球面函数

记矩阵 $A=[P(n_{-l}),\dots,P(n_{l})]$，若其可逆，则等式两边同时右乘一个 $A^{-1}$，等式变为：

$$ M=[P(R(n_{-l})),\dots,P(R(n_{l}))]A^{-1} $$

Lighting 项有 RGB 三个分量，所以这里的 $P(R(n_{i}))$ 实际上是一个三维向量

总结一下思路：

  1. 对于环境纹理的旋转，可以获取其旋转矩阵 $R$；
  2. 在每阶 $l$ 选取 $2l+1$ 个单位方向向量，获取其在该阶上的球谐系数构成矩阵 $A$；
  3. 对该矩阵求逆得到 $A^{-1}$；
  4. 使用旋转矩阵依次旋转每阶的单位方向向量，获取这些旋转后的向量在该阶上的球谐系数 $S$；
  5. 解得 $M=SA^{-1}$

主要的实现在 utils/tools.js 中，直接上代码：

  1function getRotationPrecomputeL(precompute_L, rotationMatrix){
  2
  3	// Step 1: Calc each band's rotation matrix 
  4	// Caution: for each dir the Rotation matrix is a inverse of the rotation matrix 
  5    let r = mat4Matrix2mathMatrix(rotationMatrix);
  6
  7    let shRotateMatrix3x3 = computeSquareMatrix_3by3(r);
  8    let shRotateMatrix5x5 = computeSquareMatrix_5by5(r);
  9
 10	// Step 2: Initialize rotated L    //      R   G   B
 11    let result = [];                   // SH0
 12    for(let i = 0; i < 9; i++){        // SH1
 13        result[i] = [];                // SH2
 14    }                                  // ...
 15
 16	// Step 3: Calc SH coeffs for RGB component of precompute_L
 17    for(let i = 0; i < 3; i++){
 18        let L_SH_R_3 = math.multiply([precompute_L[1][i], precompute_L[2][i], precompute_L[3][i]], shRotateMatrix3x3);
 19        let L_SH_R_5 = math.multiply([precompute_L[4][i], precompute_L[5][i], precompute_L[6][i], precompute_L[7][i], precompute_L[8][i]], shRotateMatrix5x5);
 20		
 21		let L_SH_R_3_arr = L_SH_R_3.toArray();
 22		let L_SH_R_5_arr = L_SH_R_5.toArray();
 23
 24        result[0][i] = precompute_L[0][i]; // Y1^{0} remains the same
 25		result[1][i] = L_SH_R_3_arr[0];
 26        result[2][i] = L_SH_R_3_arr[1];
 27        result[3][i] = L_SH_R_3_arr[2];
 28        result[4][i] = L_SH_R_5_arr[0];
 29        result[5][i] = L_SH_R_5_arr[1];
 30        result[6][i] = L_SH_R_5_arr[2];
 31        result[7][i] = L_SH_R_5_arr[3];
 32        result[8][i] = L_SH_R_5_arr[4];
 33    }
 34
 35	return result;
 36}
 37
 38function computeSquareMatrix_3by3(rotationMatrix){ // 计算方阵SA(-1) 3*3 
 39	
 40	// 1、pick ni - {ni}
 41	let n1 = [1, 0, 0, 0]; let n2 = [0, 0, 1, 0]; let n3 = [0, 1, 0, 0];
 42
 43	// 2、{P(ni)} - A  A_inverse
 44	let pSHn1 = SHEval(n1[0], n1[1], n1[2], 3);
 45	let pSHn2 = SHEval(n2[0], n2[1], n2[2], 3);
 46	let pSHn3 = SHEval(n3[0], n3[1], n3[2], 3);
 47	
 48	let A = math.matrix([
 49		[pSHn1[1], pSHn2[1], pSHn3[1]], // Y1^{-1} 
 50		[pSHn1[2], pSHn2[2], pSHn3[2]], // Y1^{0}
 51		[pSHn1[3], pSHn2[3], pSHn3[3]]  // Y1^{1}
 52	]);
 53	let A_inverse = math.inv(A);
 54
 55	// 3、用 R 旋转 ni - {R(ni)}
 56	let Rn1 = math.multiply(rotationMatrix, n1);
 57	let Rn2 = math.multiply(rotationMatrix, n2);
 58	let Rn3 = math.multiply(rotationMatrix, n3);
 59
 60	// get array from math matrix
 61	let Rn1_arr = Rn1.toArray();
 62	let Rn2_arr = Rn2.toArray();
 63	let Rn3_arr = Rn3.toArray();
 64
 65	// 4、R(ni) SH投影 - S
 66	let pSHRn1 = SHEval(Rn1_arr[0], Rn1_arr[1], Rn1_arr[2], 3);
 67	let pSHRn2 = SHEval(Rn2_arr[0], Rn2_arr[1], Rn2_arr[2], 3);
 68	let pSHRn3 = SHEval(Rn3_arr[0], Rn3_arr[1], Rn3_arr[2], 3);
 69
 70	let S = math.matrix([
 71		[pSHRn1[1], pSHRn2[1], pSHRn3[1]], // Y1^{-1}
 72		[pSHRn1[2], pSHRn2[2], pSHRn3[2]], // Y1^{0}
 73		[pSHRn1[3], pSHRn2[3], pSHRn3[3]]  // Y1^{1}
 74	]);
 75
 76	// 5、S*A_inverse
 77	let SA_inverse = math.multiply(S, A_inverse);
 78	return SA_inverse;
 79}
 80
 81function computeSquareMatrix_5by5(rotationMatrix){ // 计算方阵SA(-1) 5*5
 82	
 83	// 1、pick ni - {ni}
 84	let k = 1 / math.sqrt(2);
 85	let n1 = [1, 0, 0, 0]; let n2 = [0, 0, 1, 0]; let n3 = [k, k, 0, 0]; 
 86	let n4 = [k, 0, k, 0]; let n5 = [0, k, k, 0];
 87
 88	// 2、{P(ni)} - A  A_inverse
 89	let pSHn1 = SHEval(n1[0], n1[1], n1[2], 3);
 90	let pSHn2 = SHEval(n2[0], n2[1], n2[2], 3);
 91	let pSHn3 = SHEval(n3[0], n3[1], n3[2], 3);
 92	let pSHn4 = SHEval(n4[0], n4[1], n4[2], 3);
 93	let pSHn5 = SHEval(n5[0], n5[1], n5[2], 3);
 94
 95	let A = math.matrix([
 96		[pSHn1[4], pSHn2[4], pSHn3[4], pSHn4[4], pSHn5[4]], // Y2^{-2} 
 97		[pSHn1[5], pSHn2[5], pSHn3[5], pSHn4[5], pSHn5[5]], // Y2^{-1}
 98		[pSHn1[6], pSHn2[6], pSHn3[6], pSHn4[6], pSHn5[6]], // Y2^{0}
 99		[pSHn1[7], pSHn2[7], pSHn3[7], pSHn4[7], pSHn5[7]], // Y2^{1}
100		[pSHn1[8], pSHn2[8], pSHn3[8], pSHn4[8], pSHn5[8]]  // Y2^{2}
101	]);
102	let A_inverse = math.inv(A);
103
104	// 3、用 R 旋转 ni - {R(ni)}
105	let Rn1 = math.multiply(rotationMatrix, n1);
106	let Rn2 = math.multiply(rotationMatrix, n2);
107	let Rn3 = math.multiply(rotationMatrix, n3);
108	let Rn4 = math.multiply(rotationMatrix, n4);
109	let Rn5 = math.multiply(rotationMatrix, n5);
110
111	// get array from math matrix
112	let Rn1_arr = Rn1.toArray();
113	let Rn2_arr = Rn2.toArray();
114	let Rn3_arr = Rn3.toArray();
115	let Rn4_arr = Rn4.toArray();
116	let Rn5_arr = Rn5.toArray();
117
118	// 4、R(ni) SH投影 - S
119	let pSHRn1 = SHEval(Rn1_arr[0], Rn1_arr[1], Rn1_arr[2], 3);
120	let pSHRn2 = SHEval(Rn2_arr[0], Rn2_arr[1], Rn2_arr[2], 3);
121	let pSHRn3 = SHEval(Rn3_arr[0], Rn3_arr[1], Rn3_arr[2], 3);
122	let pSHRn4 = SHEval(Rn4_arr[0], Rn4_arr[1], Rn4_arr[2], 3);
123	let pSHRn5 = SHEval(Rn5_arr[0], Rn5_arr[1], Rn5_arr[2], 3);
124
125	let S = math.matrix([
126		[pSHRn1[4], pSHRn2[4], pSHRn3[4], pSHRn4[4], pSHRn5[4]], // Y2^{-2}
127		[pSHRn1[5], pSHRn2[5], pSHRn3[5], pSHRn4[5], pSHRn5[5]], // Y2^{-1}
128		[pSHRn1[6], pSHRn2[6], pSHRn3[6], pSHRn4[6], pSHRn5[6]], // Y2^{0}
129		[pSHRn1[7], pSHRn2[7], pSHRn3[7], pSHRn4[7], pSHRn5[7]], // Y2^{1}
130		[pSHRn1[8], pSHRn2[8], pSHRn3[8], pSHRn4[8], pSHRn5[8]]  // Y2^{2}
131	]);
132
133	// 5、S*A_inverse
134	let SA_inverse = math.multiply(S, A_inverse);
135	return SA_inverse;
136}

有一些要注意的点，旋转天空盒的矩阵是如下代码生成并传入 Shader 的：

 1// WebGLRenderer.js
 2let cameraModelMatrix = mat4.create();
 3mat4.fromRotation(cameraModelMatrix, timer, [0, 1, 0]);
 4
 5if (k == 'uMoveWithCamera') { // The rotation of the skybox
 6	gl.uniformMatrix4fv(
 7		this.meshes[i].shader.program.uniforms[k],
 8		false,
 9		cameraModelMatrix);
10}

对整个环境纹理进行旋转，对某一个方向 $n_{i}$ 来说，旋转后的 Lighting 相当于对这个方向应用旋转矩阵的逆，即在原始天空盒的 $R^{-1}n_{i}$ 方向的 Lighting 的值，所以传入计算球谐函数的旋转矩阵是天空盒旋转矩阵的逆（由于旋转矩阵是单位阵，所以等价于转置）。

tools.js 有一个函数 mat4Matrix2mathMatrix，按照名字应该是把 WebGL 的矩阵格式转换为 math.js 的矩阵格式：

 1function mat4Matrix2mathMatrix(rotationMatrix){
 2
 3	let mathMatrix = [];
 4	for(let i = 0; i < 4; i++){
 5		let r = [];
 6		for(let j = 0; j < 4; j++){
 7			r.push(rotationMatrix[i*4+j]);
 8		}
 9		mathMatrix.push(r);
10	}
11	return math.matrix(mathMatrix)
12
13}

看起来好像是把东西都掏出来再塞回去，不过这俩矩阵存储方式是不一样的，WebGL 是按行优先存储，而 math.js 是按列优先，如果做下打印：

1console.log("WebGL mat4:", cameraModelMatrix);
2console.log("math matrix:", mat4Matrix2mathMatrix(cameraModelMatrix));

会发现这俩看起来似乎一样，但是在内存中存储的方式有区别。WebGL 在内存中以一维数组形式排布，math.js 的矩阵在内存中按行以二维数组存储：

$$ WebGLMatrix=\begin{bmatrix} m_{00} & m_{10} & m_{20} &m_{30} \\ m_{01} & m_{11} & m_{21} &m_{31} \\ m_{02} & m_{12} & m_{22} &m_{32} \\ m_{03} & m_{13} & m_{23} &m_{33} \end{bmatrix},\quad mathMatrix=\begin{bmatrix} \left[ m_{00},m_{01},m_{02},m_{03} \right] \\ \left[ m_{10},m_{11},m_{12},m_{13} \right] \\ \left[ m_{20},m_{21},m_{22},m_{23} \right] \\ \left[ m_{30},m_{31},m_{32},m_{33} \right] \end{bmatrix} $$

所以 mat4Matrix2mathMatrix 这个函数实际上就是把 WebGL 矩阵转到 math.js 的矩阵，还做了一次转置，相当于直接拿到了旋转矩阵的逆，可以在 getRotationPrecomputeL 中直接调用这个函数获得天空盒旋转矩阵的逆，然后再传入 computeSquareMatrix_xbyx 计算球谐系数的旋转矩阵。另外，由于该矩阵是 math.js 封装的 matrix，调用 math.multiply 方法的返回值还是 matrix（即便是矩阵乘向量），所以不能通过下标访问，可以使用 ._data[i] 或者先用 .toArray 方法拿到数组，再按下标访问。

最后将 WebGLRenderer.js 第 53 和第 63 行取消注释，不出意外就能看到模型的着色正确跟随天空盒旋转！