PRT(Precomputed Radiance Transfer【2002】)原理实现

声明

本文源自对Games202课程，作业2的总结。

参考

手把手教你写GAMES202作业：GAMES202-作业2： Precomputed Radiance Transfer（球谐函数）
GAMES 202 作业2
Games202课程
个人Blog 课程总结：Games202(P6、P7)环境光照与PRT全局光照

实现目标

实现材质为Diffuse的，基于球谐函数的PRT预处理。
分别实现无阴影、有阴影、有多次光线弹射的预处理。

注意：

材质类型为Diffuse，不考虑Glossy的物体

文章目录

声明
参考
实现目标
一、架构基础
- --需要提供光线追踪算法--
- --需要提供球谐函数基函数--
- - 1. 使用硬编码
  - 2. 公式求解基函数
二、预计算核心算法实现
- 1. 预计算环境光
- 2. 预计算传输项球谐函数
- - 无阴影函数
  - 有阴影函数
  - 有阴影且多次反射函数
- 3. 保存计算得到的L和T的球谐系数
三、使用预计算数据进行渲染
四、基于漫反射材质的PRT总结
- 什么时候会用到快速的场景切换呢？
- 与光照烘培的区别
五、基于Glossy物体的PRT
六、游戏界环境是否使用PRT来做全局光照
- 那我们就没必要学习吗？？

一、架构基础

–需要提供光线追踪算法–

需要实现函数：
光线与物体是否相交：scene->rayIntersect(Ray3f ray)
光线与物体相交信息：scene->rayIntersect(Ray3f ray,Intersection its)

–需要提供球谐函数基函数–

1. 使用硬编码

保存前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)) * yreturn -0.488603 * d.y();}double HardcodedSH10(const Eigen::Vector3d& d) {// sqrt(3/(4pi)) * zreturn 0.488603 * d.z();}double HardcodedSH1p1(const Eigen::Vector3d& d) {// -sqrt(3/(4pi)) * xreturn -0.488603 * d.x();}double HardcodedSH2n2(const Eigen::Vector3d& d) {// 0.5 * sqrt(15/pi) * x * yreturn 1.092548 * d.x() * d.y();}double HardcodedSH2n1(const Eigen::Vector3d& d) {// -0.5 * sqrt(15/pi) * y * zreturn -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 * zreturn -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 * zreturn 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));}

2. 公式求解基函数

当阶数较大时，无法使用硬编码编写，则需要使用公式直接计算。

首先根据球谐函数定义，有：
得到 $A_l^m$ (归一化系数)

double kml = sqrt((2.0 * l + 1) * Factorial(l - abs(m)) /(4.0 * M_PI * Factorial(l + abs(m))));

求解 $P_l^m(cos\theta)$ ，根据球谐函数如下性质，求得 $cos\theta$ 在任意基函数的值

注：
-. 当 $P_m^m(x)$ 中 $m = 0$ 时，值为1，不是0！（这里公式错了）
-. 因为基函数可以不需要求解 $1)^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 > 0if (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 neededreturn 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 wellreturn 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 Pmlreturn pmm1;
}

当 $m < 0$ 时,根据对称性质
在这里插入图片描述
即：偶数取反，奇数不变。（因为基函数可以不考虑符号正负，故可以认为如下等式成立） $P_l^{-m} = P_l^m$

求解 $e^{im\phi} = cos(m\phi) + i·sin(m\phi)$ 。【为什么乘以 $\sqrt 2$ ，网上没有找到证明，但是确实乘以 $\sqrt 2$ 才能得到正确结果！】
- 当m>0时， $e^{im\phi} = \sqrt{2} * cos(m\phi)$
- 当m<0时， $e^{im\phi} = \sqrt{2} * sin(-m\phi)$
- 当m=0时， $e^{im\phi} = 1$

综上得出最终代码：

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));
}

二、预计算核心算法实现

1. 预计算环境光

$SH系数_{(l,m)} = \sum_i^{sampleNum} Le(i) * SH基函数_{(l,m)}(Dir(i)) * dA$
其中：

$\sum_i^{sampleNum} = \sum_{天空盒k}^6 \sum_y^{height} \sum_x^{weight}$
$i$ 为采样像素点
$D i r (i)$ 为该像素的方向向量
$d A$ 为像素立体角

其中计算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;

2. 预计算传输项球谐函数

$\int_{\Omega}\quad V(i)\quad BRDF(i,o) \quad max(n\cdot w_i,0)\quad di$
注意：若材质为Diffuse，则 $BRDF(i,0)=1/\pi$ ;

将函数展开到球谐函数上，函数有：

无阴影函数
有阴影函数
有阴影且多次反射函数

这些函数都是离散的，所以需要采样来得到。

无阴影函数

输入：顶点法线 $n$ 、采样方向 $w_i$
输出： $\frac{1}{\pi} max(n\cdot w_i，0)$

有阴影函数

输入：顶点位置 $v$ ，顶点法线 $n$ 、采样方向 $w_i$
处理： $Visibility =在v点向w_i方向发射射线，检测是否碰到物体$
输出： $\frac{1}{\pi} max(n\cdot w_i，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;}

3. 保存计算得到的L和T的球谐系数

保存为txt文档

三、使用预计算数据进行渲染

将预计算好的光照球谐系数作为Uniform参数传入Shader
将预计算好的顶点传输项系数作为顶点数据传入Shader

使用顶点着色方式着色

//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可被替代为光照烘培，即直接计算每个点在当前环境下的颜色值，并记录在一张纹理上。

因此，基于漫反射的PRT相比于光照烘培，只是多了环境光可旋转，可改变这样一个条件，而付出的代价则是9倍（3阶球谐记录）的存储空间。

那么PRT有什么优势呢？？？？？

五、基于Glossy物体的PRT

每个顶点多记录一维数据（视口俯仰角度），根据摄像机方向和物体法线方向的夹角，得到该夹角下的球谐系数，再计算得到该夹角下的颜色值。

如何得到Glossy物体系数矩阵？？？
之后再补充！！

这相比于烘培就有很大的优势了，可以实时得到光泽反射的效果（区别于镜面反射，即可以反射一定的全局光照，但只是一个很模糊的结果），这是光照烘培做不到的。

但开销同样巨大，一般情况下，光泽反射需要的阶数更高，通常需要16 * 16（4阶）或25 * 25（5阶）倍的显存开销。

六、游戏界环境是否使用PRT来做全局光照

现代实时渲染 一般不使用 PRT来做全局，
一是因为对显存开销巨大
二是因为现代GPU性能相比于预处理年代【2002】已经好了很多倍，不需要如此大量的预处理，很多计算可以放如GPU中实时计算。
三是因为PRT应用了球谐函数做拟合，本身就是一种很大的近似，而现代工业中，如果需要做全局光照的，一般是3A游戏，而PRT的效果显然达不到。