GAMES202作业2
作业要求
物体在不同光照下的表现不同，PRT(Precomputed Radiance Transfer) 是一个计算物体在不同光照下表现的方法。光线在一个环境中，会经历反射，折射，散射，甚至还会物体的内部进行散射。为了模拟具有真实感的渲染结果，传统的Path Tracing 方法需要考虑来自各个方向的光线、所有可能的传播形式并且收敛速度极慢。PRT 通过一种预计算方法，该方法在离线渲染的 Path Tracing 工具链中预计算 lighting 以及 light transport 并将它们用球谐函数拟合后储存，这样就将时间开销转移到了离线中。最后通过使用这些预计算好的数据，我们可以轻松达到实时渲染严苛的时间要求，同时渲染结果可以呈现出全局光照的效果。预计算部分是 PRT 算法的核心，也是其局限性的根源。因为在预计算 lighttransfer 时包含了 visibility 以及 cos 项，这代表着实时渲染使用的这些几何信息已经完全固定了下来（P.S. 我们可以基于球谐函数的旋转性质让光源旋转起来，本轮作业的提高部分中将简单的涉及到这一点）。所以 PRT 方法存在的限制包括：
- 不能计算随机动态场景的全局光照
- 场景中物体不可变动
本次作业的工作主要分为两个部分：cpp 端的离线预计算部分以及在 WebGL框架上使用预计算数据部分。
根据上面要求，我们把作业内容分成两部分
基于球面谐波函数计算Light项和Transform项
在WebGL使用预计算的数据实时计算光照信息
基于球面谐波函数预计算Light项和Transport项
为什么需要预计算
- 首先我们需要了解预计算辐射度传输（PRT）
预计算辐射度传输（Precomputed Radiance Transfer，PRT） 是一类通过预计算辐射度传输的基于物理的物体渲染方法。所谓辐射度传输，可以理解成物体自身的阴影/AO和表面互反射等有关于光路传输的信息。因此 PRT 往往适用于动态光照下的静态物体/静态材质
-
PRT 可以实现静态环境光照：不仅预计算 transfer，还顺便预计算环境光部分（可以称为 radiance）；这样就能以更低性能开销实现运行时基于物理的渲染，只是光照动态性会有所限制（但不一定完全静态
-
PRT 可以实现动态环境光照：仅预计算 transfer 部分 这里需要补充一下的是，在PRT之前，实时渲染的一大难题就是解决阴影和互反射问题，caustics现象就更不用说了，难上加难。主要难点在于这些物理现象要求对每个点的半球面求解入射光线的积分（渲染方程），这实在是一个很难实时完成的任务。PRT算法的核心有两个
-
将环境光照用基底编码，如球谐光照方法；
-
将入射光映射为包含阴影、反射等信息的transferred radiance，预先计算、存储起来，在实时渲染时将其直接应用到实时光照上；
如果不做PRT的话，假如一张光照贴图像素是6 * 64 * 64，对某个shading point，要得到结果颜色，至少得算6 * 64 * 64次，计算量太大。
所以预计算可以在保证效果的情况下大幅度减少渲染所消耗的时间。
如何进行预计算
预计算可以分为对Lighting部分的预计算和对Litght transport部分的预计算。
而要把这两个部分独立开来，第一时间想到的办法应该是用两张贴图去存取数据。
但仔细想想就能发现问题，就光照而言，可以用一张通用的环境光贴图，根据入射角度来查询，但是环境中每一个点的遮挡关系都是不一样的，也就是说如果用贴图的方式，每个点都需要一张遮挡关系的贴图。
也就是上面这种遮挡贴图每个点都需要一张，这样会占用大量的内存。
为了解决这个问题，球面谐波函数就出现了。
- 球面谐波函数
对于任意一个函数来说，都可以用一系列的基函数的线性组合来表示，最典型的例子就是信号处理中所提到的傅里叶变换，傅里叶变换将f(x)表示成另一系列基函数（各种频率的正弦谐波）的线性组合：
而 球谐（Spherical Harmonics，SH） 便是定义在球面上的一系列2D基函数，它与2D傅里叶序列有点相似，但非常适合球面函数 f(ω)（即参数为单位球面向量）。
一般任何式子都可以拆解为多个系数乘对应频率的球面谐波函数的基函数的和
其中的系数ci可以理解为原函数在对应基函数上的投影或者滤波
上图是球谐函数的前n阶的图像化展示，每一阶上SH的频率是相同的，并且阶级越高频率越高，能表示的信息越详细。 如果要完全复原一个任意函数，我们需要无穷阶的 SH；而如果我们只要复原一个任意函数的近似（换句话说只重建出该函数的低频信息），那么我们完全可以只需要前几阶的 SH（包含 l=0,l=1,...,l=n 每一阶的所有基函数）。
球面谐波函数也可以理解为一种滤波，把原函数的不同频率信息分开记录到不同频率的球面谐波函数中
SH 基函数的形式比较复杂，这里不作过多描述，接下来介绍用SH作为球面函数的基函数的几点好处：
- 基函数之间具有正交性（orthonormal），不同的基函数点积为0
- 通过 投影（Projection） 可以很方便得到 SH 系数（SH coefficients)
- 通过系数向量组与基函数组的点积（前n阶的基函数/系数共有 n*n 个）可以很方便重建球面函数
- 支持插值，对 SH 系数的插值相当于对重建的函数值的插值。
- 旋转不变性（rotational invariance），对函数 f的旋转 RSH等价于对 f(ω)的自变量的旋转 R3D，可以理解为光照旋转即为旋转基函数的球面图，而旋转过后的球面函数（基函数）可以被同阶的基函数的线性组合得到
用不同阶的基函数得到的结果图
- 结合上面球面谐波函数的性质，我们可以得出一个结论，任意一个二维函数的积分，都可以用球面谐波函数来描述，越高阶的球面谐波函数与原函数越拟合。
在Diffuse的物体上，Light是二维函数，View项也是二维函数，BRDF用常数替代。因为我们不需要考虑观察方向，diffuse物体中无论从哪个方向观察某一个点，对应的BRDF值是不变的。
我们只需要把Lighting部分和View和Cos组成的Transport部分分别用球面谐波函数表示即可。
编译问题
如果你用的是vs code，需要下载cmake 和mingw。当你的电脑里有vs2019或者更高的版本的时候你需要注意，直接在命令行里面用cmake ..后很可能没办法直接make。因为可能你用cmake ..调用的是vs自带的编译器MSVC，这个你可以在cmake过程中的信息看到。
如果用的是MSVC编译的是没办法用mingw的make的
你需要用vs打开nori.sln然后去生成项目
如果生成项目遇到报错
<lambda_9ed74708f63acbd4deb1a7dc36ea3ac3>::operator()
在课程论坛的讨论区可以找到答案
MSVC 对于代码中的中文字符支持有问题（应该是会吞换行符），需要启用 utf-8 编译选项：
在 prt/CMakeLists.txt 112 行添加：
target_compile_options(nori PUBLIC /utf-8) # MSVC unicode support
预计算
首先我们对渲染方程进行拆解
先把L(i)转换成用球面谐波方程表示
这个li是球面谐波函数的系数，也就是公式里的ci换了个字母，ci的计算方式如下，可以理解为原函数在各阶球面谐波函数上的投影的积分。
我们把原式化简
我们把Bi也就是基函数留在积分里，li提出积分外，在积分内的Bi**（ViCos）可以看做是球面谐波函数的系数，因为他满足计算系数ci的公式
我们把它当做系数T算，最终化简得到L(o)的结果可以变成两个球面谐波函数的各阶系数乘积的和，也就是把每一阶的liTi都加起来就可以得到光照颜色了，只要我们预计算出这两个值便可以极快的得到光照结果了。
对Lo预计算
下面的函数用于将环境光贴图的所有像素中包含的光照，都投影到基函数上去然后算出每阶的球面谐波函数的系数，然后将这些系数记录在SHCoeffiecents中，目前用的是3阶球面谐波函数，3阶球面谐波函数对应的l为2，共有九个系数
//prt.cpp
template <size_t SHOrder>
std::vector<Eigen::Array3f> PrecomputeCubemapSH(const std::vector<std::unique_ptr<float[]>> &images,
const int &width, const int &height,
const int &channel)
{
std::vector<Eigen::Vector3f> cubemapDirs;
cubemapDirs.reserve(6 * width * height);
for (int i = 0; i < 6; i++)
{
Eigen::Vector3f faceDirX = cubemapFaceDirections[i][0];
Eigen::Vector3f faceDirY = cubemapFaceDirections[i][1];
Eigen::Vector3f faceDirZ = cubemapFaceDirections[i][2];
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)
{
float u = 2 * ((x + 0.5) / width) - 1;
float v = 2 * ((y + 0.5) / height) - 1;
Eigen::Vector3f dir = (faceDirX * u + faceDirY * v + faceDirZ).normalized();
cubemapDirs.push_back(dir);
}
}
}
constexpr int SHNum = (SHOrder + 1) * (SHOrder + 1);
std::vector<Eigen::Array3f> SHCoeffiecents(SHNum);
for (int i = 0; i < SHNum; i++)
SHCoeffiecents[i] = Eigen::Array3f(0);
float sumWeight = 0;
//i对应6个cubemap
for (int i = 0; i < 6; i++)
{
//y和x分别对应cubmap中的像素
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)
{
// TODO: here you need to compute light sh of each face of cubemap of each pixel
// TODO: 此处你需要计算每个像素下cubemap某个面的球谐系数
//将cubemap上的每个点转化成其对应的方向向量dir
Eigen::Vector3f dir = cubemapDirs[i * width * height + y * width + x];
int index = (y * width + x) * channel;
//这里的Le是从环境光贴图获取的光照量，可以理解为渲染方程中的Li
Eigen::Array3f Le(images[i][index + 0], images[i][index + 1],
images[i][index + 2]);
//delta_s是用来计算cubemap上像素对应立体角权重的，因为不同角度的像素转化成立体角的大小不一样
auto delta_s = CalcArea(x, y, width, height);
//ShOrder是球面谐波函数的阶数，此时我们用的球面谐波函数阶数为3
for (int j = 0; j <= SHOrder; j++)
{
//第n阶函数对应应该有2n+1个系数，所以从k=-当前阶数开始
for (int k = -j; k <= j; k++)
{
//EvalSH的作用是求出第j阶第k个基函数在某个方向的值，下面求的方向是dir的方向
auto basic_sh_proj = sh::EvalSH(j, k, Eigen::Vector3d(dir.x(), dir.y(), dir.z()).normalized());
//我们需要把六个方向cubemap上所有点的光照都投影到第j阶第k个基函数里面，同样把投影值累加到第j阶第k个系数里面
//3阶球面谐波函数基函数一共有9个，所以每个系数我们都得算一遍
SHCoeffiecents[sh::GetIndex(j, k)] += Le * basic_sh_proj * delta_s;
}
}
}
}
}
return SHCoeffiecents;
}
}
关于CalcArea，这个函作用是归一化，使得每次投影区域的大小都是一样的，如果没有CalcArea函数，越是角落的纹素投影到球体上的单位角越小。
对Transport预计算
在计算Transport之前，我们需要先弄明白其中的3中不同情况，我们可以根据是否存在自阴影以及是否存在二次反射分出三种情况
无自遮挡 Diffuse Unshadowed
存在自遮挡Diffuse Shadowed
存在二次反射 Diffuse Inter-reflection
- Diffuse Unshadowed
如果不考虑自遮挡，则不需要考虑visibility项，所以Transport项其实可以理解为只计算max(Nx*wi,0)，其实也就是算Cos项
- Diffuse Shadowed
如果考虑自遮挡，则需要同时考虑visibility和Cos项，也就通过比较深度来决定是返回0还是Cos值
- Diffuse Inter-reflection
考虑反射的话函数式考虑自遮挡的，其主要核心就是在原本考虑自遮挡的情况，加上间接光的能量。二次光源能量计算与第一次相同，下面为公式的伪代码演算
//prt.cpp
Lo1=L(i1)*V(i1)*Cos1
Lo2=L(i2)*V(i2)*Cos2
LoSum=Lo1+Lo2*V(i1)*Cos1
因为最后要计算的是Transport的系数，所以我们这么理解：最终Transport系数=原Transport系数1+(2次光源对应点的Transport系数2Cos项2)(原Transport系数1*Cos项1)
Transport=Transport1+Transport2*Cos1
Transport1=V(i1)*Cos1
Transport2=V(i2)*Cos2
以上三种情况的代码实现
for (int i = 0; i < mesh->getVertexCount(); i++)
{
const Point3f &v = mesh->getVertexPositions().col(i);
const Normal3f &n = mesh->getVertexNormals().col(i);
auto shFunc = [&](double phi, double theta) -> double {
Eigen::Array3d d = sh::ToVector(phi, theta);
const auto wi = Vector3f(d.x(), d.y(), d.z());
//Cos代表的是wi方向与法线的点积
double Cos = wi.normalized().dot(n.normalized());
if (m_Type == Type::Unshadowed)
{
// TODO: here you need to calculate unshadowed transport term of a given direction
// TODO: 此处你需要计算给定方向下的unshadowed传输项球谐函数值
//当Transport的类型是无自遮挡的时候，Transport函数只包含Cos
return std::max(Cos,0.0);
}
else
{
// TODO: here you need to calculate shadowed transport term of a given direction
// TODO: 此处你需要计算给定方向下的shadowed传输项球谐函数值
//当Transport的类型是自遮挡的时候，Transport包含了Cos和View
//通过RayIntersect发射一根射线，v是发射点的位置，wi是发射的方向，如果不存在交点则代表该方向没有自遮挡，直接返回Cos即可
if (!scene->rayIntersect(Ray3f(v, wi.normalized())))
return std::max(Cos, 0.0);
//如果存在交点，不考虑二次反射的情况下，该点被遮挡，所以返回的View为0，自然View*Cos=0
else return 0;
}
};
auto shCoeff = sh::ProjectFunction(SHOrder, shFunc, m_SampleCount);
for (int j = 0; j < shCoeff->size(); j++)
{
m_TransportSHCoeffs.col(i).coeffRef(j) = (*shCoeff)[j]/M_PI;
}
}
//此处是多次射线的判断，若m_Type == Type::Interreflection则上面的自遮挡的条件会成立，下面的二次反射得到的系数会累加到shCoeff内
if (m_Type == Type::Interreflection)
{
// TODO: leave for bonus
for (int i = 0; i < mesh->getVertexCount(); i++)
{
//仿照上面的过程，对每个顶点做一次预计算
const Point3f& v = mesh->getVertexPositions().col(i);
const Normal3f& n = mesh->getVertexNormals().col(i);
//该函数返回的是多次弹射带来的额外系数，这里用了2次bounds，函数里是是直接迭代了两次，并没有使用递归的方法
auto shFunc = [&](double phi, double theta) -> double {
Eigen::Array3d d = sh::ToVector(phi, theta);
const auto wi = Vector3f(d.x(), d.y(), d.z());
double Cos1 = wi.dot(n);
Ray3f ray(v, wi.normalized());
// 与场景求交，若存在交点，则说明该方向有自遮挡
Intersection its;
if (Cos1> 0.0 && scene->rayIntersect(ray, its))
{
//这里对交点进行处理
//首先找这个交点位于的三角形片段，通过getVertexNormals()这个函数获取这个片段的三个顶点的各个数据
auto normal21 = mesh->getVertexNormals().col(its.tri_index.x());
auto normal22 = mesh->getVertexNormals().col(its.tri_index.y());
auto normal23 = mesh->getVertexNormals().col(its.tri_index.z());
//bary1是交点的重心坐标
const Vector3f& bary2 = its.bary;
//插值得到交点的法线
auto normal2 = bary2.x() * normal21 + bary2.y() * normal22 + bary2.z() * normal23;
//得到交点坐标
auto vertex2 = its.p;
//这里的shFunc1不考虑自遮挡情况，因为上面假设光线只弹射两次。
auto shFunc1 = [&](double phi, double theta) -> double {
Eigen::Array3d d2 = sh::ToVector(phi, theta);
const auto wi2 = Vector3f(d2.x(), d2.y(), d2.z());
double Cos2 = wi2.dot(normal2);
Ray3f ray1(vertex2, wi2.normalized());
// 与场景求交，跟上面同理
Intersection its;
if (Cos1 > 0.0 && scene->rayIntersect(ray1, its))
{
return 0;
}
return Cos2;
};
// shCoeff1是额外的Transport系数
auto shCoeff1 = sh::ProjectFunction(SHOrder, shFunc1, m_SampleCount);
//定义总Transport的SH系数
double Transport_sh_Csum=0.0;
for (int j = 0; j < shCoeff1->size(); j++)
{
Transport_sh_Csum+=Cos1* (*shCoeff1)[j];
}
return Cos1 * Transport_sh_Csum;
}
return 0;
};
// for each 采样光照方向1
auto shCoeff = sh::ProjectFunction(SHOrder, shFunc, m_SampleCount);
for (int j = 0; j < shCoeff->size(); j++)
{
m_TransportSHCoeffs.col(i).coeffRef(j) += (*shCoeff)[j];
}
}
}
- 在实现完预计算代码后，我们得把数据导出
把prt.xml填写在命令参数里
在WebGL使用预计算的数据实时计算光照信息
根据作业的说明，我们需要自定义三个函数来实现PRT渲染
-
prtVertex.glsl
-
prtFragment.glsl
-
PRTMaterial.js
同时我们还有修改一些函数让所有函数能串联在一起
PRTMaterial.js
首先我们在materials文件夹下仿照Material写一个用于PRT的PRTMaterial
这里面的uPrecomputeL[0]~uPrecomputeL[2]分别对应的是L的RGB的系数
//PRTMaterial.js
class PRTMaterial extends Material {
constructor(vertexShader, fragmentShader) {
super({
'uPrecomputeL[0]': { type: 'precomputeL', value: null},
'uPrecomputeL[1]': { type: 'precomputeL', value: null},
'uPrecomputeL[2]': { type: 'precomputeL', value: null},
},
['aPrecomputeLT'],
vertexShader, fragmentShader, null);
}
}
async function buildPRTMaterial(vertexPath, fragmentPath) {
let vertexShader = await getShaderString(vertexPath);
let fragmentShader = await getShaderString(fragmentPath);
return new PRTMaterial(vertexShader, fragmentShader);
}
在index.html补上PRTMaterial.js
//index.html
<script src="src/materials/Material.js" defer></script>
<script src="src/materials/ShadowMaterial.js" defer></script>
<script src="src/materials/PhongMaterial.js" defer></script>
<script src="src/materials/SkyBoxMaterial.js" defer></script>
<!-- Edit Start -->
<script src="src/materials/PRTMaterial.js" defer></script>
<!-- Edit End -->
在loadOBJ.js支持新PRT材质的加载
//engine.js
// TODO: load model - Add your Material here
// loadOBJ(renderer, 'assets/bunny/', 'bunny', 'addYourPRTMaterial', boxTransform);
// loadOBJ(renderer, 'assets/bunny/', 'bunny', 'addYourPRTMaterial', box2Transform);
// Edit Start
let maryTransform = setTransform(0, -35, 0, 20, 20, 20);
loadOBJ(renderer, 'assets/mary/', 'mary', 'PRTMaterial', maryTransform);
// Edit End
在WebGLRender.js中循环给precomputel实时的值
//WebGLRenderer.js
if (k == 'uMoveWithCamera') { // The rotation of the skybox
gl.uniformMatrix4fv(
this.meshes[i].shader.program.uniforms[k],
false,
cameraModelMatrix);
}
// Bonus - Fast Spherical Harmonic Rotation
//let precomputeL_RGBMat3 = getRotationPrecomputeL(precomputeL[guiParams.envmapId], cameraModelMatrix);
// Edit Start
let Mat3Value = getMat3ValueFromRGB(precomputeL[guiParams.envmapId])
for(let j = 0; j < 3; j++){
if (k == 'uPrecomputeL['+j+']') {
gl.uniformMatrix3fv(
this.meshes[i].shader.program.uniforms[k],
false,
Mat3Value[j]);
}
}
// Edit End
prtVertex.glsl
从上面的公式我们其实可以知道，其实只需要把Light的系数和Transport的系数相乘并将各阶结果累加起来就可以得到最后的颜色结果。
但是要搞清楚代码是如何实现的我们就得先了解Light和Transport的系数是如何保存和传输的。
我们先打开一份light.txt
一共有9行数据，每行数据有3个浮点数，再根据prt的代码我们可以得知，3阶Light的球面谐波函数的基函数一共有九个，每行的3个数分别代表的是RGB三个通道的系数。
我们再打开transport.txt
从上面我们不难看出，我们需要对每个点都预计算9个球面谐波函数的T系数。
我们可以从MeshRender.js208行里的 draw()函数看出每次传入prtVertex.glsl里的aPrecomputeLT的数据量是36,36刚好是9个flaot的数据量，所以每个顶点中的aPrecomputeLT包含了Transport的9个系数
for (var ii = 0; ii < 3; ++ii) {
gl.enableVertexAttribArray(this.shader.program.attribs['aPrecomputeLT'] + ii);
// void gl.vertexAttribPointer(index, size, type, normalized, stride, offset);
gl.vertexAttribPointer(this.shader.program.attribs['aPrecomputeLT'] + ii, 3, gl.FLOAT, false, 36, ii * 12);
}
从engine.js的第106行可以看出Light.txt的分割方式是换行分割，也就是Light.txt内的信息被分为了9组存在precomputeL内的
precomputeL[i] = val.split(/[(\r\n)\r\n]+/);
所以我们再WebGLRenderer64行加上这个逻辑。Mat3Value是包含了RGB三个3×3矩阵的数组
let Mat3Value = getMat3ValueFromRGB(precomputeL[guiParams.envmapId])
for(let j = 0; j < 3; j++){
if (k == 'uPrecomputeL['+j+']') {
gl.uniformMatrix3fv(
this.meshes[i].shader.program.uniforms[k],
false,
Mat3Value[j]);
}
}
这样我们传入prtVertex内uPrecomputeL数组的数据分别是Light的RGB各自对应的9个系数
//prtVertex.glsl
attribute vec3 aVertexPosition;
attribute vec3 aNormalPosition;
attribute mat3 aPrecomputeLT;
uniform mat4 uModelMatrix;
uniform mat4 uViewMatrix;
uniform mat4 uProjectionMatrix;
uniform mat3 uPrecomputeL[3];
varying highp vec3 vNormal;
varying highp mat3 vPrecomputeLT;
varying highp vec3 vColor;
//L与LT的点积，因为BRDF传的是一个值，而L传的是有RGB三个维度的颜色，所以我们需要把L分为LR、LG、LB分别计算。因为球面谐波函数一共有9个系数，所以我们用mat3类型的三维矩阵记录。
vec3 L_dot_LT(mat3 LR,mat3 LG,mat3 LB,mat3 LT)
{
vec3 result=vec3(LR[0][0],LG[0][0],LB[0][0])*LT[0][0]+
vec3(LR[0][1],LG[0][1],LB[0][1])*LT[0][1]+
vec3(LR[0][2],LG[0][2],LB[0][2])*LT[0][2]+
vec3(LR[1][0],LG[1][0],LB[1][0])*LT[1][0]+
vec3(LR[1][1],LG[1][1],LB[1][1])*LT[1][1]+
vec3(LR[1][2],LG[1][2],LB[1][2])*LT[1][2]+
vec3(LR[2][0],LG[2][0],LB[2][0])*LT[2][0]+
vec3(LR[2][1],LG[2][1],LB[2][1])*LT[2][1]+
vec3(LR[2][2],LG[2][2],LB[2][2])*LT[2][2];
return result;
}
void main(void) {
// 无实际作用，避免aNormalPosition被优化后产生警告
vNormal = (uModelMatrix * vec4(aNormalPosition, 0.0)).xyz;
mat3 LR=uPrecomputeL[0];
mat3 LG=uPrecomputeL[1];
mat3 LB=uPrecomputeL[2];
vColor=L_dot_LT(LR,LG,LB,aPrecomputeLT);
gl_Position = uProjectionMatrix * uViewMatrix * uModelMatrix * vec4(aVertexPosition, 1.0);
}
LR、LG、LB、LT分别有9个系数
由此通过L_dot_LT函数分别将LR、LG、LB系数与LT一一相乘并且累加便可以得到最终的颜色结果。
- 这里说明一下，需要更改背景为CornellBox需要修改一下engine.js的代码
//engine.js
function createGUI() {
const gui = new dat.gui.GUI();
const panelModel = gui.addFolder('Switch Environemtn Map');
// Edit Start
panelModel.add(guiParams, 'envmapId', { 'GraceGathedral': 0, 'Indoor': 1, 'Skybox': 2, 'CornellBox': 3}).name('Envmap Name');
// Edit End
panelModel.open();
}
//engine.js
var envmap = [
'assets/cubemap/GraceCathedral',
'assets/cubemap/Indoor',
'assets/cubemap/Skybox',
// Edit Start
'assets/cubemap/CornellBox',
// Edit End
];
环境光球面谐波函数的旋转
在实现球面谐波函数的旋转之前，我们需要先了解SH的性质
- SH系数的旋转是堆SH系数的线性变换
- 对每一阶SH系数的旋转可以分别进行
第一条性质意味着当我们给定SH的系数，例如SH=(SH0,......,SHk-1)，则当其旋转的时候，存在一个k阶矩阵M使得M*SH得到的结果正好是新的系数。
第二条性质代表我们可以分别处理每一阶的SH系数。我们使用3阶SH来投影，l=0对应1个系数，l=1对应3个系数，l=2对应5个系数。l=0的SH是个常量函数，我们不需要处理，当l=1的时候，我们需要一个3阶矩阵旋转3个系数，l=5时我们需要一个5阶矩阵来旋转这5个系数。
我们可以根据光源的旋转解出每阶对应的M，即可在光照大小不变的情况下求出Li。
根据上面的推倒，我们可以得出公式
M*P(x)=P(Rx)
其中R是旋转矩阵，P是求任意三维方向向量对应的多个SH值的投影函数。
我们用一阶SH来演示，先随机定义三个方向向量n0,n1,n2
M*P(n0)=P(Rn0)
M*P(n1)=P(Rn1)
M*P(n2)=P(Rn2)
当我们把P(ni)和P(Rni)都写成列向量的形式的时候，我们可以把上面式子组合成一条式子
M*[P(n0),P(n1),P(n2)]=[P(Rn0),P(Rn1),P(Rn2)]
记A=[P(n0),P(n1),P(n2)]，当A可逆，代表M可解
M=[P(Rn0),P(Rn1),P(Rn2)]*A-1
接下来我们需要修改tools.js里的代码，需要注意的是，在外界getRotationPrecomputeL(precompute_L, rotationMatrix)这条函数的时候，传入的的rotationMatrix并不是天空盒的旋转矩阵，而是旋转矩阵的逆矩阵，所以我们需要进行一个invert操作把其逆转换回来。
还有一点需要注意的是，mat4Matrix2mathMatrix(rotationMatrix)这条函数作用是吧glMAtrix和math两个数学库里的矩阵相互转化，但是这两个数学库前者矩阵是列优先，后者是行优先，mat4Matrix2mathMatrix的转换并没有考虑这一点，所以我们需要把最后的renturn值加上转置
function getRotationPrecomputeL(precompute_L, rotationMatrix){
let rotationMatrix_inverse = mat4.create()
//传入的逆矩阵，所以要逆转换回来
mat4.invert(rotationMatrix_inverse, rotationMatrix)
//把原矩阵转换成shRotateMatrix3x3能用的格式
let r = mat4Matrix2mathMatrix(rotationMatrix_inverse)
let shRotateMatrix3x3 = computeSquareMatrix_3by3(r);
let shRotateMatrix5x5 = computeSquareMatrix_5by5(r);
let result = [];
for(let i = 0; i < 9; i++){
result[i] = [];
}
for(let i = 0; i < 3; i++){
//这里要注意先后顺序,因为系数是行向量
let L_SH_R_3 = math.multiply([precompute_L[1][i], precompute_L[2][i], precompute_L[3][i]], shRotateMatrix3x3);
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);
result[0][i] = precompute_L[0][i];
result[1][i] = L_SH_R_3._data[0];
result[2][i] = L_SH_R_3._data[1];
result[3][i] = L_SH_R_3._data[2];
result[4][i] = L_SH_R_5._data[0];
result[5][i] = L_SH_R_5._data[1];
result[6][i] = L_SH_R_5._data[2];
result[7][i] = L_SH_R_5._data[3];
result[8][i] = L_SH_R_5._data[4];
}
return result;
}
function computeSquareMatrix_3by3(rotationMatrix){ // 计算方阵SA(-1) 3*3
// 1、pick ni - {ni}
let n1 = [1, 0, 0, 0]; let n2 = [0, 0, 1, 0]; let n3 = [0, 1, 0, 0];
// 2、{P(ni)} - A A_inverse
let Pn1=SHEval(n1[0],n1[1],n1[2],3);
let Pn2=SHEval(n2[0],n2[1],n2[2],3);
let Pn3=SHEval(n3[0],n3[1],n3[2],3);
//记A为[Pn1,Pn2,Pn3],Pn1、Pn2、Pn3都是列向量
let A=math.matrix(
[
[Pn1[1],Pn2[1],Pn3[1]],
[Pn1[2],Pn2[2],Pn3[2]],
[Pn1[3],Pn2[3],Pn3[3]],
]);
let A_inverse=math.inv(A);
// 3、用 R 旋转 ni - {R(ni)}
let n1_r=math.multiply(rotationMatrix,n1);
let n2_r=math.multiply(rotationMatrix,n2);
let n3_r=math.multiply(rotationMatrix,n3);
// 4、R(ni) SH投影 - S
let Pn1_r=SHEval(n1_r[0],n1_r[1],n1_r[2],3);
let Pn2_r=SHEval(n2_r[0],n2_r[1],n2_r[2],3);
let Pn3_r=SHEval(n3_r[0],n3_r[1],n3_r[2],3);
//记S为[Pn1_r,Pn2_r,Pn3_r],Pn1_r、Pn2_r、Pn3_r都是列向量
let S=math.matrix(
[
[Pn1_r[1],Pn2_r[1],Pn3_r[1]],
[Pn1_r[2],Pn2_r[2],Pn3_r[2]],
[Pn1_r[3],Pn2_r[3],Pn3_r[3]],
]);
// 5、S*A_inverse
return math.multiply(S,A_inverse);
}
function computeSquareMatrix_5by5(rotationMatrix){ // 计算方阵SA(-1) 5*5
// 1、pick ni - {ni}
let k = 1 / math.sqrt(2);
let n1 = [1, 0, 0, 0]; let n2 = [0, 0, 1, 0]; let n3 = [k, k, 0, 0];
let n4 = [k, 0, k, 0]; let n5 = [0, k, k, 0];
// 2、{P(ni)} - A A_inverse
let Pn1 = SHEval(n1[0], n1[1], n1[2], 3);
let Pn2 = SHEval(n2[0], n2[1], n2[2], 3);
let Pn3 = SHEval(n3[0], n3[1], n3[2], 3);
let Pn4 = SHEval(n4[0], n4[1], n4[2], 3);
let Pn5 = SHEval(n5[0], n5[1], n5[2], 3);
let A = math.matrix(
[
[Pn1[4], Pn2[4], Pn3[4], Pn4[4], Pn5[4]],
[Pn1[5], Pn2[5], Pn3[5], Pn4[5], Pn5[5]],
[Pn1[6], Pn2[6], Pn3[6], Pn4[6], Pn5[6]],
[Pn1[7], Pn2[7], Pn3[7], Pn4[7], Pn5[7]],
[Pn1[8], Pn2[8], Pn3[8], Pn4[8], Pn5[8]],
]);
let A_inverse = math.inv(A);
// 3、用 R 旋转 ni - {R(ni)}
let n1_r = math.multiply(rotationMatrix, n1);
let n2_r = math.multiply(rotationMatrix, n2);
let n3_r = math.multiply(rotationMatrix, n3);
let n4_r = math.multiply(rotationMatrix, n4);
let n5_r = math.multiply(rotationMatrix, n5);
// 4、R(ni) SH投影 - S
let Pn1_r = SHEval(n1_r[0], n1_r[1], n1_r[2], 3)
let Pn2_r = SHEval(n2_r[0], n2_r[1], n2_r[2], 3)
let Pn3_r = SHEval(n3_r[0], n3_r[1], n3_r[2], 3)
let Pn4_r = SHEval(n4_r[0], n4_r[1], n4_r[2], 3)
let Pn5_r = SHEval(n5_r[0], n5_r[1], n5_r[2], 3)
let S = math.matrix(
[
[Pn1_r[4], Pn2_r[4], Pn3_r[4], Pn4_r[4], Pn5_r[4]],
[Pn1_r[5], Pn2_r[5], Pn3_r[5], Pn4_r[5], Pn5_r[5]],
[Pn1_r[6], Pn2_r[6], Pn3_r[6], Pn4_r[6], Pn5_r[6]],
[Pn1_r[7], Pn2_r[7], Pn3_r[7], Pn4_r[7], Pn5_r[7]],
[Pn1_r[8], Pn2_r[8], Pn3_r[8], Pn4_r[8], Pn5_r[8]],
]);
// 5、S*A_inverse
return math.multiply(S,A_inverse);
}
function mat4Matrix2mathMatrix(rotationMatrix){
let mathMatrix = [];
for(let i = 0; i < 4; i++){
let r = [];
for(let j = 0; j < 4; j++){
r.push(rotationMatrix[i*4+j]);
}
mathMatrix.push(r);
}
return math.transpose(mathMatrix);
}
参考文章：
GAMES202-作业2： Precomputed Radiance Transfer（球谐函数）
图形学渲染基础(7) 实时环境光照(Real-time Environment Mapping)
GAMES202笔记：基于SH的GI方案
【推荐】国内首个AI IDE，深度理解中文开发场景，立即下载体验Trae
【推荐】编程新体验，更懂你的AI，立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包，你的智能百科全书，全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell：AI 加持，快人一步
· 全程不用写代码，我用AI程序员写了一个飞机大战
· DeepSeek 开源周回顾「GitHub 热点速览」
· 记一次.NET内存居高不下排查解决与启示
· MongoDB 8.0这个新功能碉堡了，比商业数据库还牛
· 白话解读 Dapr 1.15：你的「微服务管家」又秀新绝活了