GAMES101-现代计算机图形学学习笔记（作业07）
Assignment 07
原课程视频链接以及官网
b站视频链接: link.
课程官网链接: link.
作业
作业描述
在之前的练习中，我们实现了 Whitted-Style Ray Tracing 算法，并且用 BVH 等加速结构对于求交过程进行了加速。在本次实验中，我们将在上一次实验的基础上实现完整的 Path Tracing 算法。至此，我们已经来到了光线追踪版块的最后一节内容，实现光线追踪。
本次代码的流程为：
你需要从上一次编程练习中直接拷贝以下函数到对应位置：
• Triangle::getIntersection in Triangle.hpp: 将你的光线-三角形相交函数粘贴到此处，请直接将上次实验中实现的内容粘贴在此。
• IntersectP(const Ray& ray, const Vector3f& invDir, const std::array<int, 3>& dirIsNeg) in the Bounds3.hpp: 这个函数的作用是判断包围盒 BoundingBox 与光线是否相交，请直接将上次实验中实现的内容粘贴在此处，并且注意检查 t_enter = t_exit 的时候的判断是否正确。
• getIntersection(BVHBuildNode* node, const Ray ray)in BVH.cpp: BVH 查找过程，请直接将上次实验中实现的内容粘贴在此处.
在本次实验中，你只需要修改这一个函数:
• castRay(const Ray ray, int depth) in Scene.cpp: 在其中实现 Path Tracing 算法
思路
上节课的代码：
Triangle::getIntersection
inline Intersection Triangle::getIntersection(Ray ray)
{
Intersection inter;
if (dotProduct(ray.direction, normal) > 0)
return inter;
double u, v, t_tmp = 0;
Vector3f pvec = crossProduct(ray.direction, e2);
double det = dotProduct(e1, pvec);
if (fabs(det) < EPSILON)
return inter;
double det_inv = 1. / det;
Vector3f tvec = ray.origin - v0;
u = dotProduct(tvec, pvec) * det_inv;
if (u < 0 || u > 1)
return inter;
Vector3f qvec = crossProduct(tvec, e1);
v = dotProduct(ray.direction, qvec) * det_inv;
if (v < 0 || u + v > 1)
return inter;
t_tmp = dotProduct(e2, qvec) * det_inv;
if (t_tmp < 0)
{
return inter;
}
inter.distance = t_tmp;
inter.coords = ray(t_tmp);
inter.happened = true;
inter.m = m;
inter.normal = normal;
inter.obj = this;
return inter;
}
Bounds3::IntersectP
inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir, const std::array<int, 3>& dirIsNeg) const
{
// 光线进入点
float tEnter = -std::numeric_limits<float>::infinity();
// 光线离开点
float tExit = std::numeric_limits<float>::infinity();
for (int i = 0; i < 3; i++)
{
float min = (pMin[i] - ray.origin[i]) * invDir[i];
float max = (pMax[i] - ray.origin[i]) * invDir[i];
// 坐标为负的话，需要进行交换
if (!dirIsNeg[i])
{
std::swap(min, max);
}
tEnter = std::max(min, tEnter);
tExit = std::min(max, tExit);
}
return tEnter <= tExit && tExit >= 0;
}
BVHAccel::getIntersection
Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
Intersection inter;
// 光线方向
float x = ray.direction.x;
float y = ray.direction.y;
float z = ray.direction.z;
// 判断坐标是否为负
std::array<int, 3> dirsIsNeg{ int(x > 0),int(y > 0),int(z > 0) };
// 判断结点的包围盒与光线是否相交
if (node->bounds.IntersectP(ray, ray.direction_inv, dirsIsNeg) == false) return inter;
if (node->left == nullptr && node->right == nullptr)
{
inter = node->object->getIntersection(ray);
return inter;
}
// 递归判断子节点是否存在与光线相交的情况
auto hit1 = getIntersection(node->left, ray);
auto hit2 = getIntersection(node->right, ray);
if (hit1.distance < hit2.distance)
return hit1;
return hit2;
}
本次作业的代码：
①实现着色过程：
castRay(const Ray ray, int depth) in Scene.cpp: 是用来实现 Path Tracing 算法的，它会用到以下几个函数/变量：
intersect(const Ray ray) in Scene.cpp: 求一条光线与场景的交点
sampleLight(Intersection pos, float pdf) in Scene.cpp: 在场景的所有光源上按面积 uniform 地 sample 一个点，并计算该 sample 的概率密度
sample(const Vector3f wi, const Vector3f N) in Material.cpp: 按照该材质的性质，给定入射方向与法向量，用某种分布采样一个出射方向
pdf(const Vector3f wi, const Vector3f wo, const Vector3f N) in Material.cpp: 给定一对入射、出射方向与法向量，计算 sample 方法得到该出射方向的概率密度
eval(const Vector3f wi, const Vector3f wo, const Vector3f N) in Material.cpp: 给定一对入射、出射方向与法向量，计算这种情况下的 f_r 值
RussianRoulette in Scene.cpp: P_RR, Russian Roulette 的概率
需要注意的几点：
- 我们将光照对物体表面某一点的贡献分为光源和其他反射物，当光源直接能够打到物体上时，不需要再考虑其他反射物的贡献，因为这根光线直接击中了光源
- 需要判断每一根光线是否击中光源
- 没有击中光源才需要通过 Russian Roulette 来计算其他反射物的光照贡献（间接光照）
Vector3f Scene::castRay(const Ray &ray, int depth) const
{
Intersection inter = intersect(ray);
if (inter.happened)
{
// 如果射线第一次打到光源，直接返回
if (inter.m->hasEmission())
{
if (depth == 0)
{
return inter.m->getEmission();
}
else return Vector3f(0, 0, 0);
}
Vector3f L_dir(0, 0, 0);
Vector3f L_indir(0, 0, 0);
// 随机 sample 灯光，用该 sample 的结果判断射线是否击中光源
Intersection lightInter;
float pdf_light = 0.0f;
sampleLight(lightInter, pdf_light);
// 物体表面法线
auto& N = inter.normal;
// 灯光表面法线
auto& NN = lightInter.normal;
auto& objPos = inter.coords;
auto& lightPos = lightInter.coords;
auto diff = lightPos - objPos;
auto lightDir = diff.normalized();
float lightDistance = diff.x * diff.x + diff.y * diff.y + diff.z * diff.z;
Ray light(objPos, lightDir);
Intersection light2obj = intersect(light);
// 如果反射击中光源
if (light2obj.happened && (light2obj.coords - lightPos).norm() < 1e-2)
{
Vector3f f_r = inter.m->eval(ray.direction, lightDir, N);
L_dir = lightInter.emit * f_r * dotProduct(lightDir, N) * dotProduct(-lightDir, NN) / lightDistance / pdf_light;
}
if (get_random_float() < RussianRoulette)
{
Vector3f nextDir = inter.m->sample(ray.direction, N).normalized();
Ray nextRay(objPos, nextDir);
Intersection nextInter = intersect(nextRay);
if (nextInter.happened && !nextInter.m->hasEmission())
{
float pdf = inter.m->pdf(ray.direction, nextDir, N);
Vector3f f_r = inter.m->eval(ray.direction, nextDir, N);
L_indir = castRay(nextRay, depth + 1) * f_r * dotProduct(nextDir, N) / pdf / RussianRoulette;
}
}
return L_dir + L_indir;
}
return Vector3f(0, 0, 0);
}
② 多线程
我们可以将成像平面进行分块，然后分块地计算 Path Tracing
需要修改 Renderer::Render(const Scene& scene)
void Renderer::Render(const Scene& scene)
{
std::vector<Vector3f> framebuffer(scene.width * scene.height);
float scale = tan(deg2rad(scene.fov * 0.5));
float imageAspectRatio = scene.width / (float)scene.height;
Vector3f eye_pos(278, 273, -800);
int m = 0;
// 射线数量
int spp = 16;
std::cout << "SPP: " << spp << "\n";
int process = 0;
// 创造匿名函数，为不同线程划分不同块
auto castRayMultiThreading = [&](uint32_t rowStart, uint32_t rowEnd, uint32_t colStart, uint32_t colEnd)
{
for (uint32_t j = rowStart; j < rowEnd; ++j) {
int m = j * scene.width + colStart;
for (uint32_t i = colStart; i < colEnd; ++i) {
// generate primary ray direction
float x = (2 * (i + 0.5) / (float)scene.width - 1) *
imageAspectRatio * scale;
float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;
Vector3f dir = normalize(Vector3f(-x, y, 1));
for (int k = 0; k < spp; k++) {
framebuffer[m] += scene.castRay(Ray(eye_pos, dir), 0) / spp;
}
m++;
process++;
}
// 互斥锁，用于打印处理进程
std::lock_guard<std::mutex> g1(mutex_ins);
UpdateProgress(1.0*process / scene.width / scene.height);
}
};
int id = 0;
constexpr int bx = 5;
constexpr int by = 5;
std::thread th[bx * by];
int strideX = (scene.width + 1) / bx;
int strideY = (scene.height + 1) / by;
// 分块计算光线追踪
for (int i = 0; i < scene.height; i += strideX)
{
for (int j = 0; j < scene.width; j += strideY)
{
th[id] = std::thread(castRayMultiThreading, i, std::min(i + strideX, scene.height), j, std::min(j + strideY, scene.width));
id++;
}
}
for (int i = 0; i < bx*by; i++) th[i].join();
UpdateProgress(1.f);
//for (uint32_t j = 0; j < scene.height; ++j) {
// for (uint32_t i = 0; i < scene.width; ++i) {
// // generate primary ray direction
// float x = (2 * (i + 0.5) / (float)scene.width - 1) *
// imageAspectRatio * scale;
// float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;
// Vector3f dir = normalize(Vector3f(-x, y, 1));
// for (int k = 0; k < spp; k++){
// framebuffer[m] += scene.castRay(Ray(eye_pos, dir), 0) / spp;
// }
// m++;
// }
// UpdateProgress(j / (float)scene.height);
//}
//UpdateProgress(1.f);
// save framebuffer to file
FILE* fp = fopen("binary.ppm", "wb");
(void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);
for (auto i = 0; i < scene.height * scene.width; ++i) {
static unsigned char color[3];
color[0] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].x), 0.6f));
color[1] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].y), 0.6f));
color[2] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].z), 0.6f));
fwrite(color, 1, 3, fp);
}
fclose(fp);
}
并且在该 cpp 中引入头文件和定义 mutex 变量：
#include <thread>
#include <mutex>
std::mutex mutex_ins;
这里参考了 链接: https://zhuanlan.zhihu.com/p/350405670 的实现
③ 实现 Microfacet 材质
Microfacet 材质 就是考虑了微表面的材质，它在计算 BRDF 的时候，在r反射项引入了三个项来更加真实的计算 specular reflectance：
在 blinn-phong 模型中，我们这样计算光照：
而在 Microfacet 材质 中，计算方式变成了这样：
其中
r
s
r_{s}
rs 就是我们引入微表面模型来计算 specular reflectance 的项：
r
s
=
D
∗
G
∗
F
4
∗
(
n
⃗
⋅
l
⃗
)
∗
(
n
⃗
⋅
v
⃗
)
r_{s}=\frac{D * G * F}{4 *(\vec{n} \cdot \vec{l}) *(\vec{n} \cdot \vec{v})}
rs=4∗(n
⋅l
)∗(n
⋅v
)D∗G∗F
考虑能量守恒定律，所以 d + s = 1 d + s = 1 d+s=1
且菲涅尔项就表示了光在反射部分的能量占比，所以通过计算 s = f r e s n e l ( . . . ) s = fresnel(...) s=fresnel(...) , 就可以求出 d = 1 − s d = 1-s d=1−s
至于 G 和 D ，我们可以通过不同的方式进行计算，这里有一些参考：
http://simonstechblog.blogspot.com/2011/12/microfacet-brdf.html
https://learnopengl-cn.github.io/07%20PBR/02%20Lighting/.
我参考的是:
D
(
n
,
h
,
α
)
=
α
2
π
(
(
n
⋅
h
)
2
(
α
2
−
1
)
+
1
)
2
D(n, h, \alpha)=\frac{\alpha^{2}}{\pi\left((n \cdot h)^{2}\left(\alpha^{2}-1\right)+1\right)^{2}}
D(n,h,α)=π((n⋅h)2(α2−1)+1)2α2
G s u b ( n , v , k ) = n ⋅ v ( n ⋅ v ) ( 1 − k ) + k G ( n , v , l , k ) = G s u b ( n , v , k ) G s u b ( n , l , k ) G_{sub}(n, v, k)=\frac{n \cdot v}{(n \cdot v)(1-k)+k} \\ G(n, v, l, k)=G_{s u b}(n, v, k) G_{s u b}(n, l, k) Gsub(n,v,k)=(n⋅v)(1−k)+kn⋅vG(n,v,l,k)=Gsub(n,v,k)Gsub(n,l,k)
其中 α \alpha α 表示粗糙度， n n n 表示法线， h h h 表示半程向量， k k k 是基于粗糙度的项。
实现该材质，主要是通过修改 Material::eval 函数来实现的：
Vector3f Material::eval(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
switch(m_type){
case DIFFUSE:
{
// calculate the contribution of diffuse model
float cosalpha = dotProduct(N, wo);
if (cosalpha > 0.0f) {
Vector3f diffuse = Kd / M_PI;
return diffuse;
}
else
return Vector3f(0.0f);
break;
}
case Microfacet:
{
// Disney PBR 方案
float cosalpha = dotProduct(N, wo);
if (cosalpha > 0.0f) {
float roughness = 0.35;
Vector3f V = -wi;
Vector3f L = wo;
Vector3f H = normalize(V + L);
// 计算 distribution of normals: D
float D = DistributionGGX(N, H, roughness);
// 计算 shadowing masking term: G
float G = GeometrySmith(N, V, L, roughness);
// 计算 fresnel 系数: F
float F;
float etat = 1.85;
fresnel(wi, N, etat, F);
Vector3f nominator = D * G * F;
float denominator = 4 * std::max(dotProduct(N, V), 0.0f) * std::max(dotProduct(N, L), 0.0f);
Vector3f specular = nominator / std::max(denominator, 0.001f);
// 能量守恒
float ks_ = F;
float kd_ = 1.0f - ks_;
Vector3f diffuse = 1.0f / M_PI;
// 因为在 specular 项里已经考虑了折射部分的比例：F，所以折射部分不需要再乘以 ks_ （ks_ * Ks * specular）
return Ks * specular + kd_ * Kd * diffuse;
}
else
return Vector3f(0.0f);
break;
}
}
}
还需要的一些改动:
在 Material.hpp 中添加函数声明
float DistributionGGX(Vector3f N, Vector3f H, float roughness)
{
float a = roughness * roughness;
float a2 = a * a;
float NdotH = std::max(dotProduct(N, H), 0.0f);
float NdotH2 = NdotH * NdotH;
float nom = a2;
float denom = (NdotH2 * (a2 - 1.0) + 1.0);
denom = M_PI * denom * denom;
return nom / std::max(denom, 0.0000001f); // prevent divide by zero for roughness=0.0 and NdotH=1.0
}
float GeometrySchlickGGX(float NdotV, float roughness)
{
float r = (roughness + 1.0);
float k = (r*r) / 8.0;
float nom = NdotV;
float denom = NdotV * (1.0 - k) + k;
return nom / denom;
}
float GeometrySmith(Vector3f N, Vector3f V, Vector3f L, float roughness)
{
float NdotV = std::max(dotProduct(N, V), 0.0f);
float NdotL = std::max(dotProduct(N, L), 0.0f);
float ggx2 = GeometrySchlickGGX(NdotV, roughness);
float ggx1 = GeometrySchlickGGX(NdotL, roughness);
return ggx1 * ggx2;
}
在 Material.hpp 里添加枚举类型：
enum MaterialType { DIFFUSE, Microfacet};
在 Material.hpp 支持 Microfacet 材质的采样方式：
Vector3f Material::sample(const Vector3f &wi, const Vector3f &N){
switch(m_type){
case DIFFUSE:
{
// uniform sample on the hemisphere
float x_1 = get_random_float(), x_2 = get_random_float();
float z = std::fabs(1.0f - 2.0f * x_1);
float r = std::sqrt(1.0f - z * z), phi = 2 * M_PI * x_2;
Vector3f localRay(r*std::cos(phi), r*std::sin(phi), z);
return toWorld(localRay, N);
break;
}
case Microfacet:
{
// uniform sample on the hemisphere
float x_1 = get_random_float(), x_2 = get_random_float();
float z = std::fabs(1.0f - 2.0f * x_1);
float r = std::sqrt(1.0f - z * z), phi = 2 * M_PI * x_2;
Vector3f localRay(r*std::cos(phi), r*std::sin(phi), z);
return toWorld(localRay, N);
break;
}
}
}
float Material::pdf(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
switch(m_type){
case DIFFUSE:
{
// uniform sample probability 1 / (2 * PI)
if (dotProduct(wo, N) > 0.0f)
return 0.5f / M_PI;
else
return 0.0f;
break;
}
case Microfacet:
{
// uniform sample probability 1 / (2 * PI)
if (dotProduct(wo, N) > 0.0f)
return 0.5f / M_PI;
else
return 0.0f;
break;
}
}
}
sphere.hpp 中修改圆的相交判定，因为圆的相交判定精度不够，并且需要对其中的 evalDiffuseColor 函数进行返回，尽管它没有被引用到，如果不返回会报错：
Intersection getIntersection(Ray ray){
Intersection result;
result.happened = false;
Vector3f L = ray.origin - center;
float a = dotProduct(ray.direction, ray.direction);
float b = 2 * dotProduct(ray.direction, L);
float c = dotProduct(L, L) - radius2;
float t0, t1;
if (!solveQuadratic(a, b, c, t0, t1)) return result;
if (t0 < 0) t0 = t1;
if (t0 < 0) return result;
// 相交判定修改
if (t0 > 0.5) {
result.happened = true;
result.coords = Vector3f(ray.origin + ray.direction * t0);
result.normal = normalize(Vector3f(result.coords - center));
result.m = this->m;
result.obj = this;
result.distance = t0;
}
return result;
}
Vector3f evalDiffuseColor(const Vector2f &st)const {
//return m->getColor();
return {};
}
main.cpp 中添加球体:
Material* m = new Material(Microfacet, Vector3f(0.0f));
m->Ks = Vector3f(0.45, 0.45, 0.45);
m->Kd = Vector3f(0.3, 0.3, 0.25);
Sphere sphere1(Vector3f(150, 100, 300), 100, m);
MeshTriangle floor("./models/cornellbox/floor.obj", white);
//MeshTriangle shortbox("./models/cornellbox/shortbox.obj", m);
//MeshTriangle tallbox("./models/cornellbox/tallbox.obj", white);
MeshTriangle left("./models/cornellbox/left.obj", red);
MeshTriangle right("./models/cornellbox/right.obj", green);
MeshTriangle light_("./models/cornellbox/light.obj", light);
scene.Add(&floor);
scene.Add(&sphere1);
// scene.Add(&shortbox);
// scene.Add(&tallbox);
scene.Add(&left);
scene.Add(&right);
scene.Add(&light_);
修改 C++ 标准为 C++ 17：
消除 _CRT_SECURE_NO_WARNINGS 错误：
结果
SPP 16:
SPP 256：
Microfacet 材质 （SPP 16）：
参考：
http://simonstechblog.blogspot.com/2011/12/microfacet-brdf.html
https://learnopengl-cn.github.io/07%20PBR/02%20Lighting/
https://zhuanlan.zhihu.com/p/350405670