本文源自对Games202课程，作业2的总结。
实现材质为Diffuse的，基于球谐函数的PRT预处理。
分别实现无阴影、有阴影、有多次光线弹射的预处理。
注意：
需要实现函数：
光线与物体是否相交：scene->rayIntersect(Ray3f ray)
光线与物体相交信息：scene->rayIntersect(Ray3f ray,Intersection its)
保存前4阶，共16个球谐基函数。
double HardcodedSH00(const Eigen::Vector3d& d) {
// 0.5 * sqrt(1/pi)
return 0.282095;
}
double HardcodedSH1n1(const Eigen::Vector3d& d) {
// -sqrt(3/(4pi)) * y
return -0.488603 * d.y();
}
double HardcodedSH10(const Eigen::Vector3d& d) {
// sqrt(3/(4pi)) * z
return 0.488603 * d.z();
}
double HardcodedSH1p1(const Eigen::Vector3d& d) {
// -sqrt(3/(4pi)) * x
return -0.488603 * d.x();
}
double HardcodedSH2n2(const Eigen::Vector3d& d) {
// 0.5 * sqrt(15/pi) * x * y
return 1.092548 * d.x() * d.y();
}
double HardcodedSH2n1(const Eigen::Vector3d& d) {
// -0.5 * sqrt(15/pi) * y * z
return -1.092548 * d.y() * d.z();
}
double HardcodedSH20(const Eigen::Vector3d& d) {
// 0.25 * sqrt(5/pi) * (-x^2-y^2+2z^2)
return 0.315392 * (-d.x() * d.x() - d.y() * d.y() + 2.0 * d.z() * d.z());
}
double HardcodedSH2p1(const Eigen::Vector3d& d) {
// -0.5 * sqrt(15/pi) * x * z
return -1.092548 * d.x() * d.z();
}
double HardcodedSH2p2(const Eigen::Vector3d& d) {
// 0.25 * sqrt(15/pi) * (x^2 - y^2)
return 0.546274 * (d.x() * d.x() - d.y() * d.y());
}
double HardcodedSH3n3(const Eigen::Vector3d& d) {
// -0.25 * sqrt(35/(2pi)) * y * (3x^2 - y^2)
return -0.590044 * d.y() * (3.0 * d.x() * d.x() - d.y() * d.y());
}
double HardcodedSH3n2(const Eigen::Vector3d& d) {
// 0.5 * sqrt(105/pi) * x * y * z
return 2.890611 * d.x() * d.y() * d.z();
}
double HardcodedSH3n1(const Eigen::Vector3d& d) {
// -0.25 * sqrt(21/(2pi)) * y * (4z^2-x^2-y^2)
return -0.457046 * d.y() * (4.0 * d.z() * d.z() - d.x() * d.x()
- d.y() * d.y());
}
double HardcodedSH30(const Eigen::Vector3d& d) {
// 0.25 * sqrt(7/pi) * z * (2z^2 - 3x^2 - 3y^2)
return 0.373176 * d.z() * (2.0 * d.z() * d.z() - 3.0 * d.x() * d.x()
- 3.0 * d.y() * d.y());
}
double HardcodedSH3p1(const Eigen::Vector3d& d) {
// -0.25 * sqrt(21/(2pi)) * x * (4z^2-x^2-y^2)
return -0.457046 * d.x() * (4.0 * d.z() * d.z() - d.x() * d.x()
- d.y() * d.y());
}
double HardcodedSH3p2(const Eigen::Vector3d& d) {
// 0.25 * sqrt(105/pi) * z * (x^2 - y^2)
return 1.445306 * d.z() * (d.x() * d.x() - d.y() * d.y());
}
double HardcodedSH3p3(const Eigen::Vector3d& d) {
// -0.25 * sqrt(35/(2pi)) * x * (x^2-3y^2)
return -0.590044 * d.x() * (d.x() * d.x() - 3.0 * d.y() * d.y());
}
double HardcodedSH4n4(const Eigen::Vector3d& d) {
// 0.75 * sqrt(35/pi) * x * y * (x^2-y^2)
return 2.503343 * d.x() * d.y() * (d.x() * d.x() - d.y() * d.y());
}
double HardcodedSH4n3(const Eigen::Vector3d& d) {
// -0.75 * sqrt(35/(2pi)) * y * z * (3x^2-y^2)
return -1.770131 * d.y() * d.z() * (3.0 * d.x() * d.x() - d.y() * d.y());
}
double HardcodedSH4n2(const Eigen::Vector3d& d) {
// 0.75 * sqrt(5/pi) * x * y * (7z^2-1)
return 0.946175 * d.x() * d.y() * (7.0 * d.z() * d.z() - 1.0);
}
double HardcodedSH4n1(const Eigen::Vector3d& d) {
// -0.75 * sqrt(5/(2pi)) * y * z * (7z^2-3)
return -0.669047 * d.y() * d.z() * (7.0 * d.z() * d.z() - 3.0);
}
double HardcodedSH40(const Eigen::Vector3d& d) {
// 3/16 * sqrt(1/pi) * (35z^4-30z^2+3)
double z2 = d.z() * d.z();
return 0.105786 * (35.0 * z2 * z2 - 30.0 * z2 + 3.0);
}
double HardcodedSH4p1(const Eigen::Vector3d& d) {
// -0.75 * sqrt(5/(2pi)) * x * z * (7z^2-3)
return -0.669047 * d.x() * d.z() * (7.0 * d.z() * d.z() - 3.0);
}
double HardcodedSH4p2(const Eigen::Vector3d& d) {
// 3/8 * sqrt(5/pi) * (x^2 - y^2) * (7z^2 - 1)
return 0.473087 * (d.x() * d.x() - d.y() * d.y())
* (7.0 * d.z() * d.z() - 1.0);
}
double HardcodedSH4p3(const Eigen::Vector3d& d) {
// -0.75 * sqrt(35/(2pi)) * x * z * (x^2 - 3y^2)
return -1.770131 * d.x() * d.z() * (d.x() * d.x() - 3.0 * d.y() * d.y());
}
double HardcodedSH4p4(const Eigen::Vector3d& d) {
// 3/16*sqrt(35/pi) * (x^2 * (x^2 - 3y^2) - y^2 * (3x^2 - y^2))
double x2 = d.x() * d.x();
double y2 = d.y() * d.y();
return 0.625836 * (x2 * (x2 - 3.0 * y2) - y2 * (3.0 * x2 - y2));
}
当阶数较大时，无法使用硬编码编写，则需要使用公式直接计算。
首先根据球谐函数定义，有：
得到
A
l
m
A_l^m
Alm (归一化系数)
double kml = sqrt(
(2.0 * l + 1) * Factorial(l - abs(m)) /
(4.0 * M_PI * Factorial(l + abs(m)))
);
double EvalLegendrePolynomial(int l, int m, double x) {
// Compute Pmm(x) = (-1)^m(2m - 1)!!(1 - x^2)^(m/2), where !! is the double factorial.
double pmm = 1.0;
// P00 is defined as 1.0, do don't evaluate Pmm unless we know m > 0
if (m > 0) {
double sign = (m % 2 == 0 ? 1 : -1);
pmm = sign * DoubleFactorial(2 * m - 1) * pow(1 - x * x, m / 2.0);
}
if (l == m) {
// Pml is the same as Pmm so there's no lifting to higher bands needed
return pmm;
}
// Compute Pmm+1(x) = x(2m + 1)Pmm(x)
double pmm1 = x * (2 * m + 1) * pmm;
if (l == m + 1) {
// Pml is the same as Pmm+1 so we are done as well
return pmm1;
}
// Use the last two computed bands to lift up to the next band until l is
// reached, using the recurrence relationship:
// Pml(x) = (x(2l - 1)Pml-1 - (l + m - 1)Pml-2) / (l - m)
for (int n = m + 2; n <= l; n++) {
double pmn = (x * (2 * n - 1) * pmm1 - (n + m - 1) * pmm) / (n - m);
pmm = pmm1;
pmm1 = pmn;
}
// Pmm1 at the end of the above loop is equal to Pml
return pmm1;
}
当
m
<
0
m<0
m<0时,根据对称性质
即：偶数取反，奇数不变。（因为基函数可以不考虑符号正负，故可以认为如下等式成立）
P
l
−
m
=
P
l
m
P_l^{-m} = P_l^m
Pl−m=Plm
求解 e i m ϕ = c o s ( m ϕ ) + i ⋅ s i n ( m ϕ ) e^{im\phi} = cos(m\phi) + i·sin(m\phi) eimϕ=cos(mϕ)+i⋅sin(mϕ)。【为什么乘以 2 \sqrt 2 2 ，网上没有找到证明，但是确实乘以 2 \sqrt 2 2 才能得到正确结果！】
综上得出最终代码：
if (m > 0) {
return kml * sqrt(2.0) * cos(m * phi) *
EvalLegendrePolynomial(l, m, cos(theta));
}
else if (m < 0) {
return kml * sqrt(2.0) * sin(-m * phi) *
EvalLegendrePolynomial(l, -m, cos(theta));
}
else {
return kml * EvalLegendrePolynomial(l, 0, cos(theta));
}
S
H
系
数
(
l
,
m
)
=
∑
i
s
a
m
p
l
e
N
u
m
L
e
(
i
)
∗
S
H
基函
数
(
l
,
m
)
(
D
i
r
(
i
)
)
∗
d
A
SH系数_{(l,m)} = \sum_i^{sampleNum} Le(i) * SH基函数_{(l,m)}(Dir(i)) * dA
SH系数(l,m)=i∑sampleNumLe(i)∗SH基函数(l,m)(Dir(i))∗dA
其中：
其中 计算Cubemap的一个像素对应的立体角的大小原理可参照
Solid Angle of A Cubemap Texel - 计算Cubemap的一个像素对应的立体角的大小
// 计算球谐系数
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: 此处你需要计算每个像素下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]);
// 立体角
float dA = CalcArea(x, y, width, height);
// 计算球谐系数
for (int l = 0; l <= SHOrder; l++) {
for (int m = -l; m <= l; m++) {
// Eigen库不能在Vector3d，Vector3f之间相互赋值
double basic_fun_value = sh::EvalSHSlow(l, m, Eigen::Vector3d(dir.x(), dir.y(), dir.z()).normalized());
SHCoeffiecents[sh::GetIndex(l, m)] += Le * basic_fun_value * dA;
}
}
}
}
}
return SHCoeffiecents;
∫
Ω
V
(
i
)
B
R
D
F
(
i
,
o
)
m
a
x
(
n
⋅
w
i
,
0
)
d
i
\int_{\Omega}\quad V(i)\quad BRDF(i,o) \quad max(n\cdot w_i,0)\quad di
∫ΩV(i)BRDF(i,o)max(n⋅wi,0)di
注意：若材质为Diffuse，则
B
R
D
F
(
i
,
0
)
=
1
/
π
BRDF(i,0)=1/\pi
BRDF(i,0)=1/π;
将函数展开到球谐函数上，函数有：
这些函数都是离散的，所以需要采样来得到。
输入：顶点法线
n
n
n、采样方向
w
i
w_i
wi
输出：
1
π
m
a
x
(
n
⋅
w
i
，
0
)
\frac{1}{\pi} max(n\cdot w_i，0)
π1max(n⋅wi，0)
输入：顶点位置
v
v
v，顶点法线
n
n
n、采样方向
w
i
w_i
wi
处理：
V
i
s
i
b
i
l
i
t
y
=
在
v
点向
w
i
方向发射射线，检测是否碰到物体
Visibility =在v点向w_i方向发射射线，检测是否碰到物体
Visibility=在v点向wi方向发射射线，检测是否碰到物体
输出：
1
π
m
a
x
(
n
⋅
w
i
，
0
)
∗
V
i
s
i
b
i
l
i
t
y
\frac{1}{\pi} max(n\cdot w_i，0) * Visibility
π1max(n⋅wi，0)∗Visibility
需要先进行有阴影函数处理，得到每个顶点的SH系数。
对于每个顶点，递归调用函数computeInterreflectionSH（&m_TransportSHCoeffs, v, n, scene, 1）。
Func(求解顶点的球谐系数T,顶点位置v,顶点法线n,场景指针,当前弹射次数){
if (当前弹射次数>最大弹射次数) return 球谐系数为0;
对每一个顶点在球面采样
for(theta = 0->pi, phi = 0->2*pi){
1. 在采样方向与物体求交，得到求交信息（三角形顶点编号，交点位置，重心坐标）;
2. 插值计算得到相交点球谐函数系数interpolateSH。
3. 递归调用Func(相交点信息)，得到弹射返回的传输项球谐系数nextBouncesCoeffs。
4. 将interpolateSH（当前采样方向直接传输项） + nextBouncesCoeffs（当前采样方向多次弹射传输项）
}
得到采样方向的总传输项归一化(*coeffs)[i] /= sample_side * sample_side;
并将这数据return;
}
源码：
/// <summary>
/// 递归计算相互反射
/// </summary>
/// <param name="directTSHCoeffs">当前顶点传输项球谐系数</param>
/// <param name="pos">顶点位置</param>
/// <param name="normal">顶点法线</param>
/// <param name="scene">场景</param>
/// <param name="bounces">弹射次数</param>
/// <returns></returns>
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);
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++) {
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系数
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());
// 因为这里假设所有表面都是漫反射，所以每个方向的光照都是相同的。
// 因此，间接光照的值： L_间接 = L * T_间接 （L为环境光 * T为传输项）
// 因此，该方向的直接光照 = L * T_直接 + L_间接 = L * T_直接 + L * T_间接
// (注意：这里 (L_间接) 不用乘以 (1- T_直接) 因为我们做光线追踪时，只有(1- T_直接)范围内才有间接光)
// 综上：当前顶点传输项 T = T + T_间接，因此用加号直接相加即可
(*coeffs)[i] += (interpolateSH + (*nextBouncesCoeffs)[i]) * H;
}
}
}
}
// 这里不应该简单的除以sample_side * sample_side，
// 而是应该在采样中乘以立体角，而立体角=4pi/采样数
// 所以应该乘以4pi/sample_side * sample_side
// 源码
//for (unsigned int i = 0; i < coeffs->size(); i++) {
// (*coeffs)[i] /= sample_side * sample_side;
//}
// 修改代码
float invSample = 1 / sample_side * sample_side;
for (unsigned int i = 0; i < coeffs->size(); i++) {
(*coeffs)[i] = 4.0f * M_PI * invSample;
}
return coeffs;
}
保存为txt文档
使用顶点着色方式着色
//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 vColor;
float L_dot_LT(mat3 PrecomputeL, mat3 PrecomputeLT) {
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) {
for(int i = 0; i < 3; i++)
{
vColor[i] = L_dot_LT(aPrecomputeLT, uPrecomputeL[i]);
}
gl_Position = uProjectionMatrix * uViewMatrix * uModelMatrix * vec4(aVertexPosition, 1.0);
}
将颜色直接输出（经过色调映射）
//prtFragment.glsl
#ifdef GL_ES
precision mediump float;
#endif
varying highp vec3 vColor;
vec3 toneMapping(vec3 color){
vec3 result;
for (int i=0; i<3; ++i) {
result[i] = pow(color[i], 0.45);
}
return result;
}
void main(){
vec3 color = toneMapping(vColor);
gl_FragColor = vec4(color, 1.0);
}
可以看出，基于漫反射材质的PRT，本质上就是计算好每个顶点接收四面八方光照的方程，再与光照相乘得到最终结果。
我们已知，PRT要求预计算的场景不可动，只能切换或旋转球谐光源。所以PRT只能用于静态物体；而在多个球谐光源之间插值或转换，也可以理解为切换球谐光源。
所以可以用在：
如果不改变环境光照，PRT可被替代为光照烘培，即直接计算每个点在当前环境下的颜色值，并记录在一张纹理上。
因此，基于漫反射的PRT相比于光照烘培，只是多了环境光可旋转，可改变这样一个条件，而付出的代价则是9倍（3阶球谐记录）的存储空间。
那么PRT有什么优势呢？？？？？
每个顶点多记录一维数据（视口俯仰角度），根据摄像机方向和物体法线方向的夹角，得到该夹角下的球谐系数，再计算得到该夹角下的颜色值。
如何得到Glossy物体系数矩阵？？？
之后再补充！！
这相比于烘培就有很大的优势了，可以实时得到光泽反射的效果（区别于镜面反射，即可以反射一定的全局光照，但只是一个很模糊的结果），这是光照烘培做不到的。
但开销同样巨大，一般情况下，光泽反射需要的阶数更高，通常需要16 * 16（4阶）或25 * 25（5阶）倍的显存开销。
现代实时渲染 一般不使用 PRT来做全局，
一是因为对显存开销巨大
二是因为现代GPU性能相比于预处理年代【2002】已经好了很多倍，不需要如此大量的预处理，很多计算可以放如GPU中实时计算。
三是因为PRT应用了球谐函数做拟合，本身就是一种很大的近似，而现代工业中，如果需要做全局光照的，一般是3A游戏，而PRT的效果显然达不到。
不，PRT开创了通过预计算得到实时全局光照的思想。
PRT 技术或许很少用在工业中，但这种思想则被用在各种技术中！预计算得到实时全局光照的思想，之后可能会更多的应用到VR，数字孪生，开放世界中。
另外两类全局光照的开创性代表是：
这两类则更为我们所熟知！