GAMES 202 - 作业 2
作业 2: Precomputed Radiance Transfer
原课程视频链接以及官网
b站视频链接: https://www.bilibili.com/video/av887241709/.
课程官网链接: https://sites.cs.ucsb.edu/~lingqi/teaching/games202.html.
作业 2: Precomputed Radiance Transfer
总览
本次作业的重点是实现 针对环境光照和漫反射传输项的预计算球谐系数，这一部分通过 c++ 实现；第二部分是实时球谐光照计算，根据前面 求解得到的球谐系数计算环境光照，这一部分在 webgl 上实现
预计算球谐系数
预计算球谐系数需要分别将环境光照和漫反射传输项投影至球面上用球谐系数表示
环境光照
将环境光照投影至球谐系数表示主要是通过 prt.cpp 的 ProjEnv::PrecomputeCubemapSH() 函数来实现的
首先球谐系数的求解如下：
c i = ∫ Ω f ( ω ) B i ( ω ) d ω c_{i}=\int_{\Omega} f(\omega) B_{i}(\omega) \mathrm{d} \omega ci=∫Ωf(ω)Bi(ω)dω
它可以通过黎曼积分来近似求解：
c i = ∑ Ω f ( ω ) B i ( ω ) d ω c_{i}=\sum_{\Omega} f(\omega) B_{i}(\omega) \mathrm{d} \omega ci=Ω∑f(ω)Bi(ω)dω
B i ( ω ) B_{i}(\omega) Bi(ω) 表示 basic function， f ( ω ) f(\omega) f(ω) 为原函数，这里就是环境贴图的像素值
首先在环境贴图里，每个单位像素都对应了一个方向向量，整张贴图的方向实际上就是一个球面函数，将其投影在球谐系数的过程就可以描述为：将方向向量投影至球谐 basic function，然后乘以环境贴图的像素值，并累加，就得到了球谐系数，与上述方程一致
所以代码可以写成：
// 在环境贴图上每个像素
for (int i = 0; i < 6; i++)
{
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)
{
// 拿到该像素的方向
Eigen::Vector3f dir = cubemapDirs[i * width * height + y * width + x];
int index = (y * width + x) * channel;
// 该像素的光照
Eigen::Array3f Le(images[i][index + 0], images[i][index + 1],
images[i][index + 2]);
// 像素所代表的矩形区域投影到单位球面的面积
float wOmege = CalcArea(x, y, width, height);
// 投影
for (int l = 0; l <= SHOrder; l++)
{
for (int m = -l; m <= l; m++)
{
int k = sh::GetIndex(l, m);
double basisFunc = sh::EvalSH(l, m, dir.cast<double>().normalized());
// 对于 cubemap 中的每一处光线，都去累加它们把灯光 Le，以 wOmege面积投影在 basisFunc 上的系数
// 得到的结果一系列近似了环境光球面的 SH 系数
SHCoeffiecents[k] += Le * wOmege * basisFunc;
}
}
}
}
}
return SHCoeffiecents;
至于 ProjEnv::PrecomputeCubemapSH() 中其他的代码，实际上是在确定每个像素的方向而已：
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);
}
}
}
漫反射传输项
漫反射传输项分为 Diffuse Unshadowed、Diffuse Shadowed 和 Diffuse Inter-reflection(bonus)，我们需要将其投影至 SH space
这一部分的代码在 PRTIntegrator::preprocess() 中，因为课程已经帮我们写好了投影过程：
// sh::ProjectFunction 包含了投影过程
auto shCoeff = sh::ProjectFunction(SHOrder, shFunc, m_SampleCount);
所以我们只需要判断相交/阴影等情况即可
Diffuse 投影至 SH space 的总体代码：
for (int i = 0; i < mesh->getVertexCount(); i++)
{
// 获取当前 mesh 的坐标和法线
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());
double H = wi.dot(n);
if (m_Type == Type::Unshadowed)
{
if (H > 0)
{
return H;
}
return 0;
}
else
{
Ray3f ray(v, wi.normalized());
// 没有hit到其他地方的射线就对最终光照信息有贡献
if (H > 0 && !scene->rayIntersect(ray))
{
return H;
}
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];
}
}
Diffuse Unshadowed
当射线在上半球时，就投影至 SH space，描述为代码就有：
if (m_Type == Type::Unshadowed)
{
if (H > 0)
{
return H;
}
return 0;
}
Diffuse Shadowed
如果射线在上半球，且当前pos 发出的射线没有跟其他东西相交，描述为代码有：
else
{
Ray3f ray(v, wi.normalized());
// 没有hit到其他地方的射线就对最终光照信息有贡献
if (H > 0 && !scene->rayIntersect(ray))
{
return H;
}
return 0;
}
Diffuse Inter-reflection(bonus)
这一部分在 http://www.cse.chalmers.se/~uffe/xjobb/Readings/GlobalIllumination/Spherical%20Harmonic%20Lighting%20-%20the%20gritty%20details.pdf. 中的 34 页说的较为清楚了，这里就先不实现了：
实时球谐光照计算
实时球谐光照计算的关键就在于如何重构球谐光照，根据球谐函数的意义我们有：两个函数的乘积在球面空间上的积分值与它们的球谐系数向量组的点积相同
∫ S f 1 ( s ) f 2 ( s ) d s = ∑ i = 0 N c 1 i ⋅ c 2 i \int_{S} f_1(s) f_2(s) d s = \sum_{i=0}^{N} c1_{i} \cdot c2_{i} ∫Sf1(s)f2(s)ds=i=0∑Nc1i⋅c2i
假设 f 1 ( s ) f_1(s) f1(s) 投影至 SH 的系数向量为 c1， f 2 ( s ) f_2(s) f2(s) 投影至 SH 的系数向量为 c2，那么 f 1 ( s ) f_1(s) f1(s) 乘以 f 2 ( s ) f_2(s) f2(s) 的积分就等于它们在 SH 上的系数点乘；在渲染方程中， f 1 ( s ) f_1(s) f1(s) 就是环境光照项， f 2 ( s ) f_2(s) f2(s) 就是漫反射传输项，所以重构球谐光照的伪代码就可以描述如下：
for i = 1: 系数个数
light += light_SHCoeffiecents[i] * diffuse_SHCoeffiecents[i];
end
由于这里的 light 的 SH 系数是用三通道 RGB 编码，所以在重构时，只需要用每个通道上组成向量与 diffuse 的 SH 系数相乘即可：
// light SH coeff: uPrecomputeLR R channel,uPrecomputeLG G channel, uPrecomputeLB B channel
// diffuse SH coeff: vPrecomputeLT
vec3 SHrecon(mat3 uPrecomputeLR, mat3 uPrecomputeLG, mat3 uPrecomputeLB, mat3 vPrecomputeLT){
vec3 result = vec3(uPrecomputeLR[0][0], uPrecomputeLG[0][0], uPrecomputeLB[0][0]) * vPrecomputeLT[0][0] +
vec3(uPrecomputeLR[0][1], uPrecomputeLG[0][1], uPrecomputeLB[0][1]) * vPrecomputeLT[0][1] +
vec3(uPrecomputeLR[0][2], uPrecomputeLG[0][2], uPrecomputeLB[0][2]) * vPrecomputeLT[0][2] +
vec3(uPrecomputeLR[1][0], uPrecomputeLG[1][0], uPrecomputeLB[1][0]) * vPrecomputeLT[1][0] +
vec3(uPrecomputeLR[1][1], uPrecomputeLG[1][1], uPrecomputeLB[1][1]) * vPrecomputeLT[1][1] +
vec3(uPrecomputeLR[1][2], uPrecomputeLG[1][2], uPrecomputeLB[1][2]) * vPrecomputeLT[1][2] +
vec3(uPrecomputeLR[2][0], uPrecomputeLG[2][0], uPrecomputeLB[2][0]) * vPrecomputeLT[2][0] +
vec3(uPrecomputeLR[2][1], uPrecomputeLG[2][1], uPrecomputeLB[2][1]) * vPrecomputeLT[2][1] +
vec3(uPrecomputeLR[2][2], uPrecomputeLG[2][2], uPrecomputeLB[2][2]) * vPrecomputeLT[2][2];
return result;
}
理解完这一点之后，就可以根据课程的流程继续做下去了：
① 创建并调用一种使用预计算数据的材质（materials 文件夹）
class PRTMaterial extends Material {
constructor(vertexShader, fragmentShader) {
super({
'uPrecomputeLR': { type: 'updatedInRealTime', value: null },
'uPrecomputeLG': { type: 'updatedInRealTime', value: null },
'uPrecomputeLB': { type: 'updatedInRealTime', 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);
}
uPrecomputeLR、uPrecomputeLG、uPrecomputeLB 在更新时传入，放在 WebGLRenderer.js 即可：
if (k == 'uPrecomputeLR') {
gl.uniformMatrix3fv(
this.meshes[i].shader.program.uniforms[k],
false,
precomputeL_RGBMat3[0]);
}
if (k == 'uPrecomputeLG') {
gl.uniformMatrix3fv(
this.meshes[i].shader.program.uniforms[k],
false,
precomputeL_RGBMat3[1]);
}
if (k == 'uPrecomputeLB') {
gl.uniformMatrix3fv(
this.meshes[i].shader.program.uniforms[k],
false,
precomputeL_RGBMat3[2]);
}
aPrecomputeLT 是通过 MeshRender.js 里的 draw() 传入:
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);
}
其中步长为 36 ,因为每个顶点的漫反射传输项的 SH 系数是 9 个 float, 4x9=36; offset 为 ii * 12,因为在 shader 里, aVertexPosition 包含 3个float,相比 aPrecomputeLT 先传入
创建好了后就可以在 loadOBJ.js 里加载对应的材质了:
switch (objMaterial) {
case 'PhongMaterial':
material = buildPhongMaterial(colorMap, mat.specular.toArray(), light, Translation, Scale, "./src/shaders/phongShader/phongVertex.glsl", "./src/shaders/phongShader/phongFragment.glsl");
shadowMaterial = buildShadowMaterial(light, Translation, Scale, "./src/shaders/shadowShader/shadowVertex.glsl", "./src/shaders/shadowShader/shadowFragment.glsl");
break;
// TODO: Add your PRTmaterial here
case 'PRTMaterial':
material = buildPRTMaterial("./src/shaders/prtShader/prtVertex.glsl", "./src/shaders/prtShader/prtFragment.glsl");
break;
case 'SkyBoxMaterial':
material = buildSkyBoxMaterial("./src/shaders/skyBoxShader/SkyBoxVertex.glsl", "./src/shaders/skyBoxShader/SkyBoxFragment.glsl");
break;
}
PRT 材质的 shader:
顶点着色器
attribute vec3 aVertexPosition;
attribute mat3 aPrecomputeLT;
uniform mat4 uModelMatrix;
uniform mat4 uViewMatrix;
uniform mat4 uProjectionMatrix;
varying highp mat3 vPrecomputeLT;
void main(void) {
gl_Position = uProjectionMatrix * uViewMatrix * uModelMatrix * vec4(aVertexPosition, 1.0);
vPrecomputeLT = aPrecomputeLT;
}
片段着色器
#ifdef GL_ES
precision mediump float;
#endif
uniform mat3 uPrecomputeLR;
uniform mat3 uPrecomputeLG;
uniform mat3 uPrecomputeLB;
varying highp mat3 vPrecomputeLT;
vec3 SHrecon(mat3 uPrecomputeLR, mat3 uPrecomputeLG, mat3 uPrecomputeLB, mat3 vPrecomputeLT){
vec3 result = vec3(uPrecomputeLR[0][0], uPrecomputeLG[0][0], uPrecomputeLB[0][0]) * vPrecomputeLT[0][0] +
vec3(uPrecomputeLR[0][1], uPrecomputeLG[0][1], uPrecomputeLB[0][1]) * vPrecomputeLT[0][1] +
vec3(uPrecomputeLR[0][2], uPrecomputeLG[0][2], uPrecomputeLB[0][2]) * vPrecomputeLT[0][2] +
vec3(uPrecomputeLR[1][0], uPrecomputeLG[1][0], uPrecomputeLB[1][0]) * vPrecomputeLT[1][0] +
vec3(uPrecomputeLR[1][1], uPrecomputeLG[1][1], uPrecomputeLB[1][1]) * vPrecomputeLT[1][1] +
vec3(uPrecomputeLR[1][2], uPrecomputeLG[1][2], uPrecomputeLB[1][2]) * vPrecomputeLT[1][2] +
vec3(uPrecomputeLR[2][0], uPrecomputeLG[2][0], uPrecomputeLB[2][0]) * vPrecomputeLT[2][0] +
vec3(uPrecomputeLR[2][1], uPrecomputeLG[2][1], uPrecomputeLB[2][1]) * vPrecomputeLT[2][1] +
vec3(uPrecomputeLR[2][2], uPrecomputeLG[2][2], uPrecomputeLB[2][2]) * vPrecomputeLT[2][2];
return result;
}
void main(void){
vec3 result = SHrecon(uPrecomputeLR, uPrecomputeLG, uPrecomputeLB, vPrecomputeLT);
gl_FragColor = vec4(result,1);
}
环境光球谐旋转
球谐函数有两个性质利于我们推导球谐旋转:
针对性质一
对于某层 band
l
l
l, 假设我们有
2
l
+
1
2l+1
2l+1 3D normalized 法向量
n
n
n , 旋转函数
R
R
R, 能够返回该 band 上球谐投影系数的函数
P
P
P，根据上述两条性质存在一个该 band上 SH coefficient vector 的旋转矩阵
M
M
M 根据上述性质我们得出:
M
P
(
n
i
)
=
P
(
R
(
n
i
)
)
,
i
∈
[
−
l
,
l
]
MP\left(n_{i}\right)=P\left(R\left(n_{i}\right)\right), i \in[-l, l]
MP(ni)=P(R(ni)),i∈[−l,l]
整理后有:
M
[
P
(
n
−
l
)
,
…
,
P
(
n
l
)
]
=
[
P
(
R
(
n
−
l
)
)
,
…
,
P
(
R
(
n
l
)
)
]
M\left[P\left(n_{-l}\right), \ldots, P\left(n_{l}\right)\right]=\left[P\left(R\left(n_{-l}\right)\right), \ldots, P\left(R\left(n_{l}\right)\right)\right]
M[P(n−l),…,P(nl)]=[P(R(n−l)),…,P(R(nl))]
假设物体的坐标为
x
x
x,那么对物体进行旋转就可以描述如下:
M
x
=
[
P
(
R
(
n
−
l
)
)
,
…
,
P
(
R
(
n
l
)
)
]
A
−
1
x
Mx=\left[{P}\left({R}\left({n}_{-l}\right)\right), \ldots, {P}\left({R}\left({n}_{l}\right)\right)\right] {A}^{-1}x
Mx=[P(R(n−l)),…,P(R(nl))]A−1x
所以针对 band
l
l
l, 计算旋转矩阵
M
M
M 的步骤就可以整理成:
过程用代码描述如下(这里是 band 为 1 的例子):
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 pSH1 = SHEval(n1[0], n1[1], n1[2], 3);
let pSH2 = SHEval(n2[0], n2[1], n2[2], 3);
let pSH3 = SHEval(n3[0], n3[1], n3[2], 3);
let A_mat = math.matrix([[pSH1[1], pSH1[2], pSH1[3]],
[pSH2[1], pSH2[2], pSH2[3]],
[pSH3[1], pSH3[2], pSH3[3]]]);
let A_inverse_mat = math.inv(A_mat);
// 3、用 R 旋转 ni - {R(ni)}
let n1_rot = vec4.create();
let n2_rot = vec4.create();
let n3_rot = vec4.create();
vec4.transformMat4(n1_rot, n1, rotationMatrix);
vec4.transformMat4(n2_rot, n2, rotationMatrix);
vec4.transformMat4(n3_rot, n3, rotationMatrix);
// 4、R(ni) SH投影 - S
let pRotSH1 = SHEval(n1_rot[0], n1_rot[1], n1_rot[2], 3);
let pRotSH2 = SHEval(n2_rot[0], n2_rot[1], n2_rot[2], 3);
let pRotSH3 = SHEval(n3_rot[0], n3_rot[1], n3_rot[2], 3);
let S_mat = math.matrix([[pRotSH1[1], pRotSH1[2], pRotSH1[3]],
[pRotSH2[1], pRotSH2[2], pRotSH2[3]],
[pRotSH3[1], pRotSH3[2], pRotSH3[3]]]);
// 5、S*A_inverse
let result = math.multiply(S_mat, A_inverse_mat);
return result;
}
针对性质二
需要注意的是由于针对球谐的旋转是可以线性变化的,且每一个 band 的旋转是独立的,所以对于不同的 band 只需要分别进行旋转,再组合起来即可:
// band 1
let rotMatBand1 = computeSquareMatrix_3by3(rotationMatrix);
// band 2
let rotMatBand2 = computeSquareMatrix_5by5(rotationMatrix);
let result = [];
for(let i = 0; i < 3; i++)
{
// 针对不同的 band 去旋转
let rotSHBand1 = math.multiply(rotMatBand1, [precompute_L[i][1], precompute_L[i][2], precompute_L[i][3]]);
let rotSHBand2 = math.multiply(rotMatBand2, [precompute_L[i][4], precompute_L[i][5], precompute_L[i][6],
precompute_L[i][7], precompute_L[i][8]]);
// 重新组合在一起
result[i] = mat3.fromValues(precompute_L[i][0], rotSHBand1._data[0], rotSHBand1._data[1],
rotSHBand1._data[2], rotSHBand2._data[0], rotSHBand2._data[1],
rotSHBand2._data[2], rotSHBand2._data[3], rotSHBand2._data[4]);
}
结果
旋转结果1:
旋转结果2:
项目地址：
链接: https://github.com/Crocs512/GAMES202-HW.