目录
任务
C++预计算部分
1.实现prt.cpp 中的 ProjEnv::PrecomputeCubemapSH() 函数。该函数的功能是将cubemap的环境光转化为3阶球谐函数来表示（也就是9个系数）。
2.在 prt.cpp 中的 PRTIntegrator::preprocess(const Scene *scene) 函数中实现预计算过程。
JS实现部分
1.创建并调用一种使用预计算数据的材质。
2.编写该材质对应的 Shader，推荐在 vertexShader 中通过对应 SH 系数之间点积并累加来计算 vColor，接着把 vColor 传递到 fragmentShader中插值后着色。
实现
C++预计算部分
PrecomputeCubemapSH
对于下面的Original三个形状，可以将其看作为某一种颜色在各个方向上的分布。因此对RGB的三个分量，都可以通过分配球谐函数的系数来拟合。而且当阶数越高，环境光照就越接近真实的光。但环境光不需要特别清晰的效果，因此低阶拟合就能实现很好的效果。这里采用l = 2来拟合环境光。也就是说，只需要九个系数，就能把一张贴图的六个面这么复杂的光照给表示出来。
计算环境光的球谐系数，其实就是计算其在各个球谐函数上的投影。下图为其计算公式。任何一个系数，都是各个方向上的光照投影到其对应的球谐函数上的值。
因为图片是像素的，因此采用离散求和的方式计算系数。
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;
for (int i = 0; i < 6; i++)
{
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某个面的球谐系数
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]);
auto delta_w = CalcArea(x,y,width,height);
for(int l = 0;l<=SHOrder;l++){
for(int m = -l;m<=l;m++){
Eigen::Vector3d d_dir = Eigen::Vector3d( dir.x() ,dir.y() , dir.z() );
d_dir = d_dir.normalized();
auto SH = sh::EvalSH(l,m,d_dir);
SHCoeffiecents[ sh::GetIndex(l,m) ] += Le * SH * delta_w;
}
}
}
}
}
return SHCoeffiecents;
}
}preprocess
对于diffuse材质，其BRDF可以变成如下形式。光源已经提前计算成九个系数了，因此在这里是常数，直接提取出来。想要知道每个顶点的光照，再把剩下的部分也预计算就好了。剩下来的积分部分就是我们需要对每个顶点预计算的部分。该积分在这里采用蒙特卡洛方法对周围随机采样求得。
这里的H就是max(0,n·i)，V通过判断光线是否和场景有交点得出。有交点说明被遮挡，V为0。
virtual void preprocess(const Scene *scene) override
{
// Here only compute one mesh
const auto mesh = scene->getMeshes()[0];
// Projection environment
auto cubePath = getFileResolver()->resolve(m_CubemapPath);
auto lightPath = cubePath / "light.txt";
auto transPath = cubePath / "transport.txt";
std::ofstream lightFout(lightPath.str());
std::ofstream fout(transPath.str());
int width, height, channel;
std::vector<std::unique_ptr<float[]>> images =
ProjEnv::LoadCubemapImages(cubePath.str(), width, height, channel);
auto envCoeffs = ProjEnv::PrecomputeCubemapSH<SHOrder>(images, width, height, channel);
m_LightCoeffs.resize(3, SHCoeffLength);
for (int i = 0; i < envCoeffs.size(); i++)
{
lightFout << (envCoeffs)[i].x() << " " << (envCoeffs)[i].y() << " " << (envCoeffs)[i].z() << std::endl;
m_LightCoeffs.col(i) = (envCoeffs)[i];
}
std::cout << "Computed light sh coeffs from: " << cubePath.str() << " to: " << lightPath.str() << std::endl;
// Projection transport
m_TransportSHCoeffs.resize(SHCoeffLength, mesh->getVertexCount());
fout << mesh->getVertexCount() << std::endl;
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());
const double H = 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传输项球谐函数值
if(H > 0.0)return H;
return 0;
}
else
{
// TODO: here you need to calculate shadowed transport term of a given direction
// TODO: 此处你需要计算给定方向下的shadowed传输项球谐函数值
if(H > 0.0 && !scene->rayIntersect(Ray3f( v , wi.normalized())))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] / M_PI;
}
}
if (m_Type == Type::Interreflection)
{
// TODO: leave for bonus
}
// Save in face format
for (int f = 0; f < mesh->getTriangleCount(); f++)
{
const MatrixXu &F = mesh->getIndices();
uint32_t idx0 = F(0, f), idx1 = F(1, f), idx2 = F(2, f);
for (int j = 0; j < SHCoeffLength; j++)
{
fout << m_TransportSHCoeffs.col(idx0).coeff(j) << " ";
}
fout << std::endl;
for (int j = 0; j < SHCoeffLength; j++)
{
fout << m_TransportSHCoeffs.col(idx1).coeff(j) << " ";
}
fout << std::endl;
for (int j = 0; j < SHCoeffLength; j++)
{
fout << m_TransportSHCoeffs.col(idx2).coeff(j) << " ";
}
fout << std::endl;
}
std::cout << "Computed SH coeffs"
<< " to: " << transPath.str() << std::endl;
}JS实现部分
创建材质
按照作业的要求对属性命名，直接抄Bling-Phong材质的代码稍作修改即可
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引入该文件
<script src="src/materials/PRTMaterial.js" defer></script>编写Shader
编写的重点在VertexShader里，FragmentShader使用插值来获取颜色。
使用上面的公式，因为l和T两项都已经预计算好了，只需要将传入的数据在顶点着色器中进行相乘与求和即可。
这里直接将九个系数分给三个vec3来计算，减少代码量。
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 vec3 vColor;
float L_dot_LT(mat3 PrecomputeLT,mat3 PrecomputeL){
vec3 L_0 = PrecomputeL[0];
vec3 L_1 = PrecomputeL[1];
vec3 L_2 = PrecomputeL[2];
vec3 LT_0 = PrecomputeLT[0];
vec3 LT_1 = PrecomputeLT[1];
vec3 LT_2 = PrecomputeLT[2];
return dot(L_0 , LT_0) + dot(L_1 , LT_1) + dot(L_2 , LT_2);
}
void main(void) {
vNormal = (uModelMatrix * vec4(aNormalPosition, 0.0)).xyz;
vColor[0] = L_dot_LT(aPrecomputeLT,uPrecomputeL[0]);
vColor[1] = L_dot_LT(aPrecomputeLT,uPrecomputeL[1]);
vColor[2] = L_dot_LT(aPrecomputeLT,uPrecomputeL[2]);
gl_Position = uProjectionMatrix * uViewMatrix * uModelMatrix *
vec4(aVertexPosition, 1.0);
}FragmentShader部分
#ifdef GL_ES
precision mediump float;
#endif
varying vec3 vColor;
void main(void) { gl_FragColor = vec4(vColor, 1.0); }结果
预计算的预览结果(skybox)
在html中运行的结果
Bonus部分
球谐函数的旋转
根据作业文档里，球谐函数的两个性质
1.旋转不变形性
假设有一个旋转R，对原函数f(x) 的旋转R(f(x)) 与直接旋转f(x) 的自变量是一样的，即
R(f(x)) = f(
R
(
x
))
2.对每层 band 上的 SH coeff icient，可以分别在上面进行旋转，并且这个旋转是线性变化。每个系数不会对其他系数产生影响。
直接沿用作业文档的说明
对于球谐函数上某层 band l
，该
band
上有一个
2
l
+1
维的、由
SH coeff icient 构成的向量。假设我们有 2
l
+ 1
个任意的
3
D normal vector n
，旋转函数
R
， 能够返回该 band
上球谐投影系数的函数
P
，根据上述两条性质存在一个该
band上 SH coeff icient vector
的旋转矩阵
M
使得：
整理一下得到：
记 A
= [
P
(
n
−
l
)
, ..., P
(
n
l
)]
。如果矩阵
A
是可逆矩阵，此时
SH coeff icient vector
的旋转矩阵
M
易得，即：
要实现求M,直接按照上述的流程做就好了。
function getRotationPrecomputeL(precompute_L, rotationMatrix){
let rotationMatrix_inverse = mat4.create();
mat4.invert(rotationMatrix_inverse,rotationMatrix);
let rotationMatrix_math = mat4Matrix2mathMatrix(rotationMatrix_inverse)
let matrix3x3 = computeSquareMatrix_3by3(rotationMatrix_math);
let matrix5x5 = computeSquareMatrix_5by5(rotationMatrix_math);
let result = [];
for(let i =0;i<9;i++){
result[i] = [];
}
for(let i =0;i<3;i++){
let result_3 = math.multiply( [precompute_L[1][i],precompute_L[2][i],precompute_L[3][i]],matrix3x3 );
let result_5 = math.multiply( [precompute_L[4][i],precompute_L[5][i],precompute_L[6][i],precompute_L[7][i],precompute_L[8][i] ],matrix5x5 );
result[0][i] = precompute_L[0][i];
result[1][i] = result_3._data[0];
result[2][i] = result_3._data[1];
result[3][i] = result_3._data[2];
result[4][i] = result_5._data[0];
result[5][i] = result_5._data[1];
result[6][i] = result_5._data[2];
result[7][i] = result_5._data[3];
result[8][i] = result_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 sh1 = SHEval(n1[0],n1[1],n1[2],3);
let sh2 = SHEval(n2[0],n2[1],n2[2],3);
let sh3 = SHEval(n3[0],n3[1],n3[2],3);
let A = math.matrix([
[sh1[1],sh2[1],sh3[1]],
[sh1[2],sh2[2],sh3[2]],
[sh1[3],sh2[3],sh3[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 sh1_r = SHEval(n1_r[0],n1_r[1],n1_r[2],3);
let sh2_r = SHEval(n2_r[0],n2_r[1],n2_r[2],3);
let sh3_r = SHEval(n3_r[0],n3_r[1],n3_r[2],3);
let S = math.matrix([
[sh1_r[1],sh2_r[1],sh3_r[1]],
[sh1_r[2],sh2_r[2],sh3_r[2]],
[sh1_r[3],sh2_r[3],sh3_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 sh1 = SHEval(n1[0],n1[1],n1[2],3);
let sh2 = SHEval(n2[0],n2[1],n2[2],3);
let sh3 = SHEval(n3[0],n3[1],n3[2],3);
let sh4 = SHEval(n4[0],n4[1],n4[2],3);
let sh5 = SHEval(n5[0],n5[1],n5[2],3);
let A = math.matrix([
[sh1[4],sh2[4],sh3[4],sh4[4],sh5[4]],
[sh1[5],sh2[5],sh3[5],sh4[5],sh5[5]],
[sh1[6],sh2[6],sh3[6],sh4[6],sh5[6]],
[sh1[7],sh2[7],sh3[7],sh4[7],sh5[7]],
[sh1[8],sh2[8],sh3[8],sh4[8],sh5[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 sh1_r = SHEval(n1_r[0],n1_r[1],n1_r[2],3);
let sh2_r = SHEval(n2_r[0],n2_r[1],n2_r[2],3);
let sh3_r = SHEval(n3_r[0],n3_r[1],n3_r[2],3);
let sh4_r = SHEval(n4_r[0],n4_r[1],n4_r[2],3);
let sh5_r = SHEval(n5_r[0],n5_r[1],n5_r[2],3);
let S = math.matrix([
[sh1_r[4],sh2_r[4],sh3_r[4],sh4_r[4],sh5_r[4]],
[sh1_r[5],sh2_r[5],sh3_r[5],sh4_r[5],sh5_r[5]],
[sh1_r[6],sh2_r[6],sh3_r[6],sh4_r[6],sh5_r[6]],
[sh1_r[7],sh2_r[7],sh3_r[7],sh4_r[7],sh5_r[7]],
[sh1_r[8],sh2_r[8],sh3_r[8],sh4_r[8],sh5_r[8]],
])
// 5、S*A_inverse
return math.multiply(S,A_inverse);
}计算间接光照的球谐系数
原理其实和之前在101做的光线追踪一样，只不过将光换成了球谐函数的系数表示而已。但是这里对代码框架不是非常熟悉，参照了大佬的作业https://zhuanlan.zhihu.com/p/596050050。对着代码打完一遍后才明白了一点。以下直接上代码与批注
std::unique_ptr<std::vector<double>> computeInterreflectionSH(Eigen::MatrixXf* directTSHCoeffs,
const Point3f& pos,
const Normal3f& normal,
const Scene* scene,
int bounces)
{
//初始化系数
std::unique_ptr<std::vector<double>> coeffs(new std::vector<double>());
coeffs->assign(SHCoeffLength,0.0);
//如果达到了预定的最大次数，直接返回全0
if( bounces > m_Bounce)return coeffs;
//采样区域大小
const int sample_side = static_cast<int>(floor(sqrt(m_SampleCount)));
//初始化随机数生成器
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> rng(0.0,1.0);
//对指定的点进行周围光线的采样
for(int t = 0;t<sample_side;t++){
for(int p = 0; p < sample_side; p++){
//将随机取到的球面坐标值转化为xyz的坐标值
double alpha = (t + rng(gen)) / sample_side;
double beta = (p + rng(gen)) / sample_side;
double phi = 2.0 * M_PI * beta;
double theta = acos(2.0 * alpha - 1.0);
Eigen::Array3d d = sh::ToVector(phi , theta);
const auto wi = Vector3f(d.x(),d.y(),d.z());
double H = wi.normalized().dot(normal);
Intersection its;
//如果发出的光线打到物体，计算间接光照
if(H > 0.0 && scene->rayIntersect(Ray3f(pos,wi.normalized()),its)){
MatrixXf normals = its.mesh->getVertexNormals();
Point3f idx = its.tri_index;
Point3f hitPos = its.p;
Vector3f bary = its.bary;
//这里是根据三角形的重心坐标对法向量进行插值得到对应点的法向量
Normal3f hitNormal =
Normal3f(normals.col(idx.x()).normalized()*bary.x() +
normals.col(idx.y()).normalized()*bary.y() +
normals.col(idx.z()).normalized()*bary.z()
).normalized();
//递归查询下一个点的系数
auto nextBouncesCoeffs = computeInterreflectionSH(directTSHCoeffs,hitPos,hitNormal,scene,bounces + 1);
//因为只有顶点有SH系数，这里对SH系数属性也进行了一次重心坐标插值
for( int i =0; i < SHCoeffLength; i++){
auto interpolateSH = (directTSHCoeffs->col(idx.x()).coeffRef(i) * bary.x() +
directTSHCoeffs->col(idx.y()).coeffRef(i) * bary.y() +
directTSHCoeffs->col(idx.z()).coeffRef(i) * bary.z()
);
(*coeffs)[i] += (interpolateSH + (*nextBouncesCoeffs)[i]) * H;
}
}
}
}
for (unsigned int i = 0;i < coeffs->size(); i++){
(*coeffs)[i] /= sample_side * sample_side;
}
return coeffs;
}实现这个函数后，在Interreflection材质里将调用的代码补上
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);
auto indirectCoeffs = computeInterreflectionSH(&m_TransportSHCoeffs, v, n, scene,1);
for(int j= 0; j < SHCoeffLength ; j++){
m_TransportSHCoeffs.col(i).coeffRef(j) += (*indirectCoeffs)[j];
}
}
}
左图为shadowed材质的结果，右图为interreflection的2次bounce的结果，可以看到部分区域因为计算了间接光照，明亮了很多。
最后附上在CornellBox里球谐旋转的结果
4817
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
