PRT预计算辐射传输方法
PRT（Precomputed Radiance Transfer）技术是一种用于实时渲染全局光照的方法。它通过预计算光照传输来节省时间，并能够实时重现面积光源下3D模型的全局光照效果。
由于PRT方法的局限，它不能计算随机动态场景的全局光照，场景中物体也不可变动。
Basic Idea
- 光的传输与背景光照内容本身无关，因此两部分拆开计算。
- 环境光照使用球谐函数拟合，来近似表示
Diffuse Case
对于完全粗糙表面，brdf系数是一个常量，因此：
Precompute
1.球谐函数表达环境光照，离线计算和保存
公式：
\[l_i = \int_{\Omega}{L(\omega)B_i(\omega)d\omega}
\]
code:
float CalcPreArea(const float& x, const float& y)
{
return std::atan2(x * y, std::sqrt(x * x + y * y + 1.0));
}
float CalcArea(const float& u_, const float& v_, const int& width,
const int& height)
{
// transform from [0..res - 1] to [- (1 - 1 / res) .. (1 - 1 / res)]
// ( 0.5 is for texel center addressing)
float u = (2.0 * (u_ + 0.5) / width) - 1.0;
float v = (2.0 * (v_ + 0.5) / height) - 1.0;
// shift from a demi texel, mean 1.0 / size with u and v in [-1..1]
float invResolutionW = 1.0 / width;
float invResolutionH = 1.0 / height;
// u and v are the -1..1 texture coordinate on the current face.
// get projected area for this texel
float x0 = u - invResolutionW;
float y0 = v - invResolutionH;
float x1 = u + invResolutionW;
float y1 = v + invResolutionH;
float angle = CalcPreArea(x0, y0) - CalcPreArea(x0, y1) -
CalcPreArea(x1, y0) + CalcPreArea(x1, y1);
return angle;
}
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);
constexpr int Order = SHOrder;
cout << "SHNum:" << SHNum << endl;
cout << "Order:" << Order << endl;
std::vector<Eigen::Array3f> SHCoeffiecents(SHNum);
for (int i = 0; i < SHNum; i++)
SHCoeffiecents[i] = Eigen::Array3f(0);
float sumWeight = 0;
//float u_ = 0, v_ = 0;
int bindex = 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]);
//Le[0] = 0.5f; Le[1] = 0.5f; Le[2] = 0.5f;//for test
// 从这里开始
//u_ = x / (float)width;
//v_ = y / (float)height;
float area = CalcArea(x, y, width, height);
// 对于每个基函数
// L^2 + M+L+1 = index
bindex = 0;
for (int L = 0; L <= Order; L++) {
//cout << "l:" << L << endl;
for (int M = -L; M <= L; M++) {
//cout << "m:" << M << endl;
//cout << "bindex:" << bindex << endl;
//if (bindex < SHNum)
{
// 对于每个颜色
Eigen::Vector3d dird = Eigen::Vector3d(dir.x(), dir.y(), dir.z());
dird.normalize();
float sh = sh::EvalSH(L, M, dird) * area;
SHCoeffiecents[bindex][0] += Le[0] * sh;
SHCoeffiecents[bindex][1] += Le[1] * sh;
SHCoeffiecents[bindex][2] += Le[2] * sh;
}
bindex++;
}
}
}
}
}
return SHCoeffiecents;
}
结果：
light.txt
三列分别对应r,g,b的9个基函数系数数据。
2.离线计算光传播并保存
T有9个分量，对应于分别使用2阶球谐函数的9个基函数作为光照，在式子中的积分结果。
code:
std::unique_ptr<std::vector<double>> ProjectFunction(
int order, const SphericalFunction& func, int sample_count) {
CHECK(order >= 0, "Order must be at least zero.");
CHECK(sample_count > 0, "Sample count must be at least one.");
// This is the approach demonstrated in [1] and is useful for arbitrary
// functions on the sphere that are represented analytically.
const int sample_side = static_cast<int>(floor(sqrt(sample_count)));
std::unique_ptr<std::vector<double>> coeffs(new std::vector<double>());
coeffs->assign(GetCoefficientCount(order), 0.0);
// generate sample_side^2 uniformly and stratified samples over the sphere
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;
// See http://www.bogotobogo.com/Algorithms/uniform_distribution_sphere.php
double phi = 2.0 * M_PI * beta;
double theta = acos(2.0 * alpha - 1.0);
// evaluate the analytic function for the current spherical coords
double func_value = func(phi, theta);
// evaluate the SH basis functions up to band O, scale them by the
// function's value and accumulate them over all generated samples
for (int l = 0; l <= order; l++) {
for (int m = -l; m <= l; m++) {
double sh = EvalSH(l, m, phi, theta);
(*coeffs)[GetIndex(l, m)] += func_value * sh;
}
}
}
}
// scale by the probability of a particular sample, which is
// 4pi/sample_side^2. 4pi for the surface area of a unit sphere, and
// 1/sample_side^2 for the number of samples drawn uniformly.
double weight = 4.0 * M_PI / (sample_side * sample_side);
for (unsigned int i = 0; i < coeffs->size(); i++) {
(*coeffs)[i] *= weight;
}
return coeffs;
}
for (int i = 0; i < mesh->getVertexCount(); i++)
{
const Point3f& v = mesh->getVertexPositions().col(i);
const Normal3f& n = mesh->getVertexNormals().col(i);
double dot = 0.0f;
double* dd = ˙
bool intersect = false;
auto shFunc = [&](double phi, double theta) -> double {
Eigen::Array3d d = sh::ToVector(phi, theta);
const auto wi = Vector3f(d.x(), d.y(), d.z());
if (m_Type == Type::Unshadowed)
{
cout << "Unshadowed" << endl;
// TODO: here you need to calculate unshadowed transport term of a given direction
// TODO: 此处你需要计算给定方向下的unshadowed传输项球谐函数值
//*dd = 0.1f;//有奇怪的编译器优化
//cout << "dd:" << *dd << endl;
*dd = wi.x() * n.x() + wi.y() * n.y() + wi.z() * n.z();// wi.dot(n);
//cout << "dd:" << *dd << endl;
//cout << "wi:" << wi.x()<<" "<<wi.y()<<" " << wi.z() << endl;
//cout << "n:" << n.x() << " " << n.y() << " " << n.z() << endl;
if (*dd > 0) {
//cout << "dd:" << *dd << endl;
return *dd;
}
return 0;
}
else
{
// TODO: here you need to calculate shadowed transport term of a given direction
// TODO: 此处你需要计算给定方向下的shadowed传输项球谐函数值
cout << "Shadowed" << endl;
intersect = (*scene).rayIntersect(Ray3f(v, wi));
//cout << "intersect:"<< intersect << endl;
if (!intersect) {
*dd = wi.x() * n.x() + wi.y() * n.y() + wi.z() * n.z();// wi.dot(n);
//cout << "dd:" << *dd << endl;
//cout << "wi:" << wi.x()<<" "<<wi.y()<<" " << wi.z() << endl;
//cout << "n:" << n.x() << " " << n.y() << " " << n.z() << endl;
if (*dd > 0) {
//cout << "dd:" << *dd << endl;
return *dd;
}
}
return 0;
}
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];
}
}
result: transport.txt
行数为模型顶点数，行数据为每个顶点处的对应每个基函数光照的积分结果。
Runtime Render
code:
在顶点着色器中计算：
attribute mat3 aPrecomputeLT;
attribute vec3 aVertexPosition;
attribute vec3 aNormalPosition;
attribute vec2 aTextureCoord;
uniform mat4 uModelMatrix;
uniform mat4 uViewMatrix;
uniform mat4 uProjectionMatrix;
uniform float uPrecomputeL[27];
varying highp vec2 vTextureCoord;
varying highp vec3 vFragPos;
varying highp vec3 vNormal;
varying highp vec3 vColor;
void main(void) {
vFragPos = (uModelMatrix * vec4(aVertexPosition, 1.0)).xyz;
vNormal = (uModelMatrix * vec4(aNormalPosition, 0.0)).xyz;
gl_Position = uProjectionMatrix * uViewMatrix * uModelMatrix *
vec4(aVertexPosition, 1.0);
vTextureCoord = aTextureCoord;
float r = 0.0;
float g = 0.0;
float b = 0.0;
int index = 0;
for(int i=0;i<3;i++){
for(int j=0;j<3;j++){
r += aPrecomputeLT[i][j]*uPrecomputeL[index];
g += aPrecomputeLT[i][j]*uPrecomputeL[index+1];
b += aPrecomputeLT[i][j]*uPrecomputeL[index+2];
index +=3;
}
}
//r=0.0;g=0.0;b=0.0;
//r = uPrecomputeL[26];
//r = aPrecomputeLT[0][0];
vColor = vec3(r,g,b);
}
其中uPrecomputeL是长度为3*9的float数组，对应于light.txt的数据。
aPrecompute是3x3矩阵，存储了transport.txt对应顶点的9个数据。
result:
优点：运行时计算非常简单和快速
缺点：
- 物体位置旋转不能发生变化，否则光传播发生改变
- 在shadow和inter-reflection case中考虑了物体的自遮挡和内部反射，但不会考虑物体之间的遮挡与反射，可以认为是一种局部(local)光照方法。
