球谐函数，就是一组在球面上的基函数
其表达式为
$$
y_l^m=
\left\{
\begin{array}{lr}
\sqrt{2}K^{m}_ {l}cos(m\phi)P^m_ {l}(cos\theta), & m>0 \\
\sqrt{2}K^{m}_ {l}sin(-m\phi)P^{-m}_ {l}(cos\theta), & m < 0 \\
K^0_lP^0_ {l}(cos\theta), & m = 0
\end{array}\right.
$$
其中的符号K只是个系数，表示为
$$ K^m_ {l}=\sqrt{\frac{2l + 1}{4\pi}\frac{(l - |m|)!}{(l+|m|)!}} $$
其中的符号$P^m_l$是伴随勒让德多项式（associated LegendrePolynomials）
伴随勒让德多项式是一个可以递归表达的形式
$$ \begin{array}{lr} (l-m)P^m_l(x)=x(2l-1)p^m_ {l-1}(x)-(l+m-1)P^m_ {l-2}(x)\\ P^{m}_ {l}(x)=(-1)^l(2l-1)!!(1-x^2)^{\frac{l}{2}}\\ P^l_ {l+1}=x(2l+1)P^l_ {l}(x)\ \end{array} $$
对于勒让德伴随多项式，$l = 0,1,2......$，$m=0,1..,l$，l表示了伴随勒让德多项式的次数
而对于球谐函数，$l = 0,1,2......$，$m=-l,..,0,1,..l$，其中l被称为Band(段？带？) ，对于给定的l，在l段就会有2l+1个球谐函数，所以，对于l阶(order)的球谐函数，就有$(2l+1)^2$个球谐函数。
也可以用一个数字i来表示某个球谐函数，那么$i=l(l+1)+m$
好吧，对于这里我讲的很“简”，因为这些东西只有在计算的时候才用得到
这边我直接摘抄了《全局光照技术: 从离线到实时渲染》10.1，作者已经把书在github开源了，可以去康康
基函数是球谐函数光照的基础。
可以把基函数理解为类似坐标轴，比如说对于二维的笛卡尔坐标系，其有X轴和Y轴，我们可以用坐标轴的线性组合来描述一个点，或者一个向量。那么基函数就是函数空间中的基，可以用基函数的线性组合描述一个连续的函数。
给定一个定义在T上的函数f(t)以及一组基函数Bi(t)，该函数可以通过下面的积分被投影（projection）成系数
$$c_i = \int f(t)B_i(t)dt$$
同样的，原始的函数可以通过基函数和其对应系数的乘积之和的形重建
$$
f(t)\approx\tilde{f}(t)=\sum_N^{i=0}c_iB_i(t)
$$
可以理解为系数向量和基函数向量的点乘
这里N表示系数的数量，一个连续函数通常需要无穷多个基函数的线性组合才能完全被重建，当基函数的数量有限，或者N小于基函数的数量时，上述基函数的线性组合只是原始函数的一个近似，记为$\tilde{f}(t)$
再引入一个概念正交多项式
正交多项式（orthogonal polynomials）是一组具有一些有趣特性的多项式，即任意两个基函数之间都是相互正交的，所以当对任意两个基函数的乘积进行积分时，只有当两者完全相同时才会取得非零值
而对于标准正交基函数（orthonormal basis functions），两个多项式乘积的积分必须返回0或者1，如下
$$
\int_1^{-1}F_m(x)F_n(x)dx=\left\{
\begin{array}{lr}
0 ，&m!=n
\\
1 ,&m=n
\end{array}\right.
$$
由多项式组成的基函数还具有另一种特征，不同阶的基函数能够表述原函数不同频率域的特征，阶数越高则保留的频率细节越多，但是需要的系数数量也越多，因此这给出一种可能，即使用少量的低阶多项式基函数来近似原函数的低频部分，这即是预计算辐射传输算法的核心思想。
球谐函数是定义在球面上的标准正交基，而光照很显然也可以表示成关于 $\phi,\theta$ 的函数，这代表着我们可以把光照预计算，投影到球谐函数上，然后在使用的时候重建，这样可以用低廉的性能换取极致的画面。
对一个定义在球面上的函数，计算器在球谐函数上的投影，只需要将原函数f和球谐函数y的乘积进行积分。
$$
c_l^m=\int_Sf(s)y_l^m(s)ds
$$
相反的，如果要重建这个函数，只需要将其系数与对应的球谐函数进行相乘求和，也就是求点乘
$$
\tilde{f}(s)=\sum^{n^2}_{i=0}c_iy_i(s)$$
其中，s代表了一个三维的单位向量$(x,y,z)$，当然也可以将其参数化为$(\theta,\phi)$
球谐函数的另一个很神奇的性质就是
$$\int \tilde{a}(s)\tilde{b}(s)\approx\sum_{i=0}^{n^2} a_ib_i$$
两个定义在球面上的函数的乘积的积分可以近似为其分别在球谐函数上的投影的系数向量的点乘
这个性质可以由其的标准正交基性质推出
$$\int \tilde{a}(s)\tilde{b}(s)\approx\int_S\sum^{n^2}{i=0}a_iy_i(s)\sum^{n^2}{j=0}b_jy_j(s)$$
然后想象一下，首先把系数提出积分 ，又只有当两个y相同的时候其乘积才为1，不相同的时候为0，于是可以消掉一大半，最后合并，只剩下了一个求和，是不是很神奇。
SH - visualizer by iq
根据这个性质，可以对渲染方程进行拆分
把积分里的式子拆成两项,光照项和光照传输项
那么，对于每一项，都用投影的公式计算其投影到球谐函数上的系数，在shader里再进行重建就行了。
这里还需要进一步进行讨论，对于diffuse和glossy的物体我们需要做不同的处理
对于diffuse，其brdf是个常数，于是可以把BRDF项提出渲染方程，这样LightTransport项就成为了只关于$V(w_i)cos(\theta)$的函数，极大的减少了计算量
但是对于Glossy的物体，其brdf会根据不同的出射方向和入射方向而变化，那么LightTransport项就成为了一个四维的函数（把wi和wo参数化为($\theta，\phi$)）
那么需要把LightTransport项二次投影到球谐函数上，成为一个系数矩阵
那么最后还原这一步，也就从点乘变成了矩阵乘法......
作业用了nori作为框架
首先的首先，我遇到了个大坑。编译问题
我一开始用的是Mingw-w64进行编译，但无奈花式编译都不通过，然后我用了室友的电脑进行编译....嗯不是我的问题....
接着我用visual studio2022进行编译...嗯，也是不通过
就在我接近绝望的时候，我看官方网站上说要用visual studio 2019，而且语气特别坚定
Be sure to select the Visual Studio 2019 64 bit compiler.
抱着孤注一掷的心态装了.....然后成了
草
还有一个问题是Visual studio上可能遇到缺少')'等这种错误，解决方法是在项目根目录CMakeLists.txt中112行加上一行
target_compile_options(nori PUBLIC /utf-8) # MSVC unicode support
这部分要在prt/src/prt.cpp中完成
首先，需要把一个CubeMap投影到球谐函数上，作为全局光照
在这一步，其实就是在进行L项的计算
需要修改函数：
std::vector<Eigen::Array3f> PrecomputeCubemapSH(const std::vector<std::unique_ptr<float[]>> &images,const int &width, const int &height,const int &channel)
经过观察法线TODO部分在三个for循环内部，最外层循环是对于cubeMap的六个面进行一个大循环，而内部两个循环是对面上每个像素进行遍历。
其实，这里是想让你使用这个公式
$$c_i = \int f(t)B_i(t)dt$$
只不过把积分换成了累加和而已。重新写一边就是
$$c_i = \sum^{pix}_{cubeMap} L_i(w_i)B_i(w_i)dw_i$$
这里的dw也就是每个像素投影到球面上的面积大小，作为平衡权重了。
球谐函数的计算函数在ext/spherical-harmonics/sh/中已经给出，直接使用即可；$w_i$的计算已经给出，也就是$Eigen::Vector3f dir = cubemapDirs[i * width * height + y * width + x];$
$dw$的计算只需要调用函数CalcArea()就可以了
于是乎不难写出下面的代码
// 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]);
double theta,phi;
sh::ToSphericalCoords(dir,&phi,&theta);
double weight = CalcArea(x, y, width, height);
for (int l = 0; l <= SHNum; l++) {
for (int m = -l; m <= l; m++) {
int h = sh::GetIndex(l, m);
double sh = sh::EvalSH(l, m, phi, theta);
SHCoeffiecents[h] += sh * weight * Le.array();
}
}
作业中只用了三阶的球谐函数，效果还不错
BTW,这个代码是可以抄到的，ext/spherical-harmonics/sh/spherical_harmonics.cc中的函数std::unique_ptr<std::vector<Eigen::Array3f>> ProjectEnvironment( int order, const Image& env)就有这一段.....
还有一点需要注意，由于BRDF是$\frac{\rho}{\pi}$（而$\rho$假定为1），那么光照需要再乘上一个系数$\pi$
为了节省计算时间放到此处的三层循环外做
for(auto &&i : SHCoeffiecents){
i = i * INV_PI
}
放到后面计算lightTransport里做也不是不行..
Light Transport项被分为了三种，
Unshadow, shadowed，interrflection
这三种的区别是，unshadowed不会考虑光线被自己和其他物体遮挡的情况，shadowed考虑到了这点，interrflection则进一步考虑了光线多次弹射的情况。
那么对于Unshadowed，其LT项就是cos项，而shadowed项额外需要一个visibility项。
ps: 这个visibility是顶点到场景的可见性，而不是摄像机到顶点的可见性哦
double H = n.dot(wi);
if (m_Type == Type::Unshadowed)
{
// TODO: here you need to calculate unshadowed transport term of a given direction
// TODO: 此处你需要计算给定方向下的unshadowed传输项球谐函数值
if(H < 0)return 0;
return H;
}
else
{
// TODO: 此处你需要计算给定方向下的shadowed传输项球谐函数值
Ray3f mRay(v,wi.normalized());
if(H > 0 && !scene->rayIntersect(mRay)){
return H;
}
return 0;
}
所以shadowed的计算就是比Unshadowed的计算多了一个于场景相交的判断而已
interreflection的部分不想做，感觉就是一个小号的rt积分器...
渲染的结果会显示在屏幕上，也会保存在prt/scenes中的prt.png，球谐函数的系数被保存在了prt/scenes/cubemap里的对应场景文件夹下
需要注意的是此处有四个场景需要构建，需要在prt/scenes/prt.xml中手动修改
天空场景的渲染结果
在./homework2中需要重建环境光，相较来说这部分就比较复杂了
因为不仅需要会一点webgl，还得看得懂作业框架
js真不习惯
此处先构建好我们的shader
放在shaders/PRTMaterial/下
首先在vertex shader里计算光照，再送入frag插值
in prt.vert
#ifdef GL_ES
precision mediump float;
#endif
//light
uniform vec3 uL0;
uniform vec3 uL1;
uniform vec3 uL2;
uniform vec3 uL3;
uniform vec3 uL4;
uniform vec3 uL5;
uniform vec3 uL6;
uniform vec3 uL7;
uniform vec3 uL8;
uniform mat4 uModelMatrix;
uniform mat4 uViewMatrix;
uniform mat4 uProjectionMatrix;
//light transport
attribute mat3 aPrecomputeLT;
varying highp vec3 vFragPos;
attribute vec3 aNormalPosition;
varying highp vec3 vNormal;
attribute vec3 aVertexPosition;
varying highp vec3 albedo;
vec3 getLight(){
vec3 res = vec3(0);
res = res + uL0 * aPrecomputeLT[0][0];
res = res + uL1 * aPrecomputeLT[0][1];
res = res + uL2 * aPrecomputeLT[0][2];
res = res + uL3 * aPrecomputeLT[1][0];
res = res + uL4 * aPrecomputeLT[1][1];
res = res + uL5 * aPrecomputeLT[1][2];
res = res + uL6 * aPrecomputeLT[2][0];
res = res + uL7 * aPrecomputeLT[2][1];
res = res + uL8 * aPrecomputeLT[2][2];
return res;
}
void main(void){
albedo = getLight();
//无实际意义，放着好看
vFragPos = (uModelMatrix * vec4(aVertexPosition, 1.0)).xyz;
//无实际意义，放着好看
vNormal = (uModelMatrix * vec4(aNormalPosition, 0.0)).xyz;
gl_Position = uProjectionMatrix *uViewMatrix* uModelMatrix *vec4(aVertexPosition, 1.0);
}
我此处直接选择用9个uniform vec3来传入球谐函数的系数
其中uL0代表了band0的球谐函数光照系数，
uL1-3是band1，uL4-8是band2。由于光照是三色，所以是三位向量
aPrecomputeLT是每个顶点的LightTransport，其是三维矩阵，m00代表band1，m01,m02,m10代表band1........LightTransport就是一个浮点数，aPrecomputeLT在框架的renderers/MeshRender.js中送入shader
在函数getLight()里，进行对应系数的累加和，也就是公式
$$\int \tilde{a}(s)\tilde{b}(s)\approx\sum_{i=0}^{n^2} a_ib_i$$
最后把光照数据传入了frag shader中
in prt.frag
#ifdef GL_ES
precision mediump float;
#endif
varying highp vec3 albedo;
vec3 gammaCorrection(vec3 color){
return pow(color,vec3(1. / 2.2) );
}
void main(){
gl_FragColor = vec4(gammaCorrection(albedo),1. );
}
这里做了个gamma校正，原因是看网上有人说nori渲染的问题....
不过校正后观感确实可以力
接下来更具shader来构建一个prtMaterial
in materials/PRTMaterial.js
class PRTMaterial extends Material {
constructor(preComputeL,vertexShader, fragmentShader) {
super({
'uL0' : { type: '', value: null },
'uL1' : { type: '', value: null },
'uL2' : { type: '', value: null },
'uL3' : { type: '', value: null },
'uL4' : { type: '', value: null },
'uL5' : { type: '', value: null },
'uL6' : { type: '', value: null },
'uL7' : { type: '', value: null },
'uL8' : { type: '', value: null },
},
['aPrecomputeLT'], vertexShader, fragmentShader, null);
}
}
async function buildPRTMaterial(precomputeL,vertexPath, fragmentPath) {
let vertexShader = await getShaderString(vertexPath);
let fragmentShader = await getShaderString(fragmentPath);
return new PRTMaterial(precomputeL,vertexShader, fragmentShader);
}
父类构造函数的第一个参数，是shader里的Unform变量，框架中对这个的处理可以去renderers/MeshRender.js里查看，发现其有一个根据类型的type对应处理传Uniform逻辑
这里有两种选择：
其一是什么都不填，或者乱填个没有的type，然后这部分uniform的传递会延缓到后续进行，方便进行不同的环境贴图切换和bonus-球谐函数光照的快速旋转
另一种是把使用的贴图和数据写死在此处，做法如下，写死了使用第0张贴图
'uL0' : { type: '3fv', value: precomputeL[0][0] },
那么下一步就是注册我们的shader和材质了
在loads/loadOBJ.js中找到TODO
然后创建自己的材质
// TODO: Add your PRTmaterial here
case 'PRTMaterial':
console.log('loading material')
material = buildPRTMaterial(precomputeL,"./src/shaders/PRTMaterial/prt.vert","./src/shaders/PRTMaterial/prt.frag")
break;
接下来，需要去renderers/MeshRender.js传递uniform了
in render()
render()函数是每帧调用一次，所以很多计算不用重复做
网上有些作业答案真这么干了，每帧每有一个uniform就计算一次，真的牛皮
在函数的开头先把光照系数复制一份
let lightArray = []
lightArray = precomputeL.slice();
然后找到bonus的那个TODO
在下面这么写
if(k.charAt(0)=='u' && k.charAt(1)=='L'){
gl.uniform3fv(this.meshes[i].shader.program.uniforms[k],lightArray[guiParams.envmapId][pos])
}
这样就可以随着不同的场景图切换环境光贴图辣
环境光贴图的选项并没有关于CornellBox的那张，你得自己手动加上.....
而且SkyBox那张渲染的和实际的不同....
球谐函数有旋转不变性，所以如果我们旋转了一个向量然后再把它投影到SH上，我们会得到和先投影再旋转相同的结果，而且是线性的，其的旋转矩阵是对每个band分别进行计算的,然后分别旋转
所以对于i band其应该有一个(i+1)*(i+1)的方阵作为旋转矩阵
这里有几种方法供选择
第零种是最复杂的做法，可以手算球谐函数的旋转矩阵每个参数
根据公式
$M_{ij} = \int_{s}y_i(Rs)y_j(s)ds $
可以把这里看作是将旋转前的球谐函数投影到了旋转后的球谐函数上
计算出每个band的旋转矩阵
第一种是OldFashion的ZXZXZ方法
按照$R_{SH}(\alpha,\beta,\gamma)=Z_\gamma X_{-90}Z_{\beta}X_{90}Z_{\alpha}$的顺序分别计算其球谐函数的旋转矩阵，最后合一起拼成Rsh
你可能问为啥子不直接XYZ
因为球谐光照的旋转矩阵对于Z轴的旋转具有简单的形式，计算量小
(说白了就是按照下图填表，第0band是常数不用动，1 band就是框起来的3x3，以此类推..
想要更高系数的？自己算去
而其中的两个X，第一个是旋转-90，第二个旋转+90
由于常数，所以可以预存取
这个方法在论文 Spherical Harmonic Lighting 中被提到，其还有一个更加快速的Trick是用泰勒展开近似
你也可以去看 《全局光照技术》，书中有讲解
这种方法需要几百次乘法...开销很大
第二种是快速球谐函数旋转方法，算法特别特别简单
由于球谐函数有旋转不变性，所以如果我们旋转了一个向量然后再把它投影到SH上，我们会得到和先投影再旋转相同的结果
ok，那就有了如下的式子
$$R\cdot P(N)=P(M\cdot N)$$
其中，R是球谐函数的旋转矩阵，M是对于向量的旋转矩阵，而N是任意一个标准化的向量（真的是任意），P是球谐函数，负责投影向量到系数向量的
对于i band上的系数向量都有
$$ R[P(N_1),....P(N_{2i+1})] = [P(M\cdot N_1),....P(M\cdot N_{2i+1})] $$
对于等式左侧的那一坨系数向量矩阵，简称为A
$$ R \cdot A = [P(M\cdot N_1),....P(M\cdot N_{2i+1})] $$
$$ R = [P(M\cdot N_1),....P(M\cdot N_{2i+1})]\cdot A^{-1} $$
那么R就是要求的球谐函数旋转矩阵
P(N)是列向量！！！！！！！！！！！
都是列向量！！！！！
大坑，卡了我两天
这里N的选取需要满足两点
代码我直接放下面了
const k = 1./math.sqrt(2)
let vecs = [[1,0,0,0],
[0,0,1,0],
[k,k,0,0],
[k,0,k,0],
[0,k,k,0]]
let rotVecs = []
function computeInvA3x3(){
shN1 = SHEval3(vecs[0][0],vecs[0][1],vecs[0][2])
shN2 = SHEval3(vecs[1][0],vecs[1][1],vecs[1][2])
shN3 = SHEval3(vecs[2][0],vecs[2][1],vecs[2][2])
matA = math.matrix(
[shN1.slice(1,4),shN2.slice(1,4),shN3.slice(1,4)]
)
matA = math.transpose(matA)
return math.inv(matA)
}
function computeInvA5x5(){
shN1 = SHEval3(vecs[0][0],vecs[0][1],vecs[0][2])
shN2 = SHEval3(vecs[1][0],vecs[1][1],vecs[1][2])
shN3 = SHEval3(vecs[2][0],vecs[2][1],vecs[2][2])
shN4 = SHEval3(vecs[3][0],vecs[3][1],vecs[3][2])
shN5 = SHEval3(vecs[4][0],vecs[4][1],vecs[4][2])
matA = math.matrix(
[shN1.slice(4,9),shN2.slice(4,9),shN3.slice(4,9),shN4.slice(4,9),shN5.slice(4,9)]
)
matA = math.transpose(matA)
return math.inv(matA)
}
function computeRotVec(rot){
rotVecs = new Array(vecs.length)
for (let i = 0; i < rotVecs.length; i++) {
element = vecs[i];
let iRot = vec4.create()
vec4.transformMat4(iRot,element,rot)
rotVecs[i] = iRot
}
}
function compute3x3M(){
let invA = computeInvA3x3()
rot1 = SHEval3(rotVecs[0][0],rotVecs[0][1],rotVecs[0][2])
rot2 = SHEval3(rotVecs[1][0],rotVecs[1][1],rotVecs[1][2])
rot3 = SHEval3(rotVecs[2][0],rotVecs[2][1],rotVecs[2][2])
let Pn = math.matrix(
[rot1.slice(1,4),rot2.slice(1,4),rot3.slice(1,4)]
)
Pn = math.transpose(Pn)
let m = math.multiply(Pn,invA)
return m
}
function compute5x5M(){
let invA = computeInvA5x5()
rot1 = SHEval3(rotVecs[0][0],rotVecs[0][1],rotVecs[0][2])
rot2 = SHEval3(rotVecs[1][0],rotVecs[1][1],rotVecs[1][2])
rot3 = SHEval3(rotVecs[2][0],rotVecs[2][1],rotVecs[2][2])
rot4 = SHEval3(rotVecs[3][0],rotVecs[3][1],rotVecs[3][2])
rot5 = SHEval3(rotVecs[4][0],rotVecs[4][1],rotVecs[4][2])
let Pn = math.matrix(
[rot1.slice(4,9),rot2.slice(4,9),rot3.slice(4,9),rot4.slice(4,9),rot5.slice(4,9)]
)
Pn = math.transpose(Pn)
let m = math.multiply(Pn,invA)
return m
}
function computeBand1Light(light){
let lightBand1 = light[guiParams.envmapId].slice(1,4)
let M = compute3x3M()
return math.multiply(M,math.matrix(lightBand1))
}
function computeBand2Light(light){
let lightBand1 = light[guiParams.envmapId].slice(4,9)
let coefRGB = math.matrix(lightBand1)
let M = compute5x5M()
return math.multiply(M,(coefRGB))
}
其中INV_A矩阵只需要一次计算就可以了
在render函数中的开头
//更新旋转后的向量
computeRotVec(cameraModelMatrix);
//计算band1的系数矩阵
let coefBand1 = computeBand1Light(lightArray)
//计算band2的系数矩阵
let coefBand2 = computeBand2Light(lightArray)
在TODO下面
if(k.charAt(0)=='u' && k.charAt(1)=='L'){
let pos = k.charCodeAt(2) - '0'.charCodeAt()
let coef = []
if(pos == 0){
coef = lightArray[guiParams.envmapId][0]
}else if(pos <=3){
//coef = lightArray[guiParams.envmapId][pos]
coef = coefBand1.subset(math.index(pos-1,[0,1,2]))._data[0]
}else{
//coef = lightArray[guiParams.envmapId][pos]
coef = coefBand2.subset(math.index(pos-4,[0,1,2]))._data[0]
}
gl.uniform3fv(this.meshes[i].shader.program.uniforms[k],coef)
}
完美结束
《全局光照技术：从离线渲染到实时渲染》
Simple and Fast Spherical Harmonic Rotation - Filmicworlds
Spherical Harmonic Lighting: The Gritty Details