基础题部分

根据教程PDF，首先需要引用如下函数（在作业5的基础上稍作修改）：

renderer() in Renderer.cpp

解说见注释：

// The main render function. This where we iterate over all pixels in the image,
// generate primary rays and cast these rays into the scene. The content of the
// framebuffer is saved to a file.
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(-1, 5, 10);int m = 0;for (uint32_t j = 0; j < scene.height; ++j) {for (uint32_t i = 0; i < scene.width; ++i) {// generate primary ray direction// TODO: Find the x and y positions of the current pixel to get the// direction//  vector that passes through it.// Also, don't forget to multiply both of them with the variable// *scale*, and x (horizontal) variable with the *imageAspectRatio*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)); // Don't forget to normalize this direction!// 在本次代码框架中，castRay()函数放在了Scene类中// 并且castRay()的传入参数修改为(const Ray &ray, int depth)// 所以要先声明Ray ray(const Vector3f& ori, const Vector3f& dir)Ray ray(eye_pos, dir);// 调用scene.castRay()// 并将返回结果(Vector3f类型的像素RGB颜色)存入framebuffer中framebuffer[m++] = scene.castRay(ray, 0);}UpdateProgress(j / (float)scene.height);}UpdateProgress(1.f);// save framebuffer to fileFILE* 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 * clamp(0, 1, framebuffer[i].x));color[1] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].y));color[2] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].z));fwrite(color, 1, 3, fp);}fclose(fp);
}

Triangle::getIntersection in Triangle.hpp

解说见注释：

// TODO find ray triangle intersection
inline Intersection Triangle::getIntersection(Ray ray)
{Intersection inter;if (dotProduct(ray.direction, normal) > 0)return inter;// u,v为三角形重心坐标，即PPT中的b1,b2// t_tmp是PPT中的tdouble 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;// 如果有intersection，则赋值给Intersection inter并将inter返回// Intersection的定义：见Intersection.hppinter.happened = true; // 光线和该三角形有交叉点inter.coords = ray(t_tmp); // 利用公式o+td，求得intersection（交点）坐标inter.normal = normal; // 传入该三角形的法线inter.distance = t_tmp; // 与光源点的距离，用t表示即可inter.obj = this; // 指针指向自身inter.m = m; // 传入该三角形材质(material)return inter;
}

然后实现BVH的功能：

IntersectP() in the Bounds3.hpp

这里使用PPT中的公式：

如果光线方向是相反的（即xyz轴相应的方向向量为负数），则要交换tmin和tmax，这里闫老师的PPT没讲，需要特别注意！

代码和详解如下：

// TODO test if ray bound intersects
inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,const std::array<int, 3>& dirIsNeg) const
{// invDir: ray direction(x,y,z), invDir=(1.0/x,1.0/y,1.0/z), use this because Multiply is faster that Division// invDir就是1/Dir，因为乘法的运算比除法快，故使用invDir来加快运算速度// dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x<0),int(y<0),int(z<0)], use this to simplify your logic// 上面这一句注释在源代码中写错了，应该是[int(x<0),int(y<0),int(z<0)]float tmin_x = (pMin.x - ray.origin.x) * invDir.x;float tmin_y = (pMin.y - ray.origin.y) * invDir.y;float tmin_z = (pMin.z - ray.origin.z) * invDir.z;float tmax_x = (pMax.x - ray.origin.x) * invDir.x;float tmax_y = (pMax.y - ray.origin.y) * invDir.y;float tmax_z = (pMax.z - ray.origin.z) * invDir.z;// swap t_min and t_max if direction is negative// 如果[int(x<0),int(y<0),int(z<0)]，需要交换tmin和tmax// 这里PPT上面没有讲，需要特别注意if (dirIsNeg[0]) std::swap(tmin_x, tmax_x);if (dirIsNeg[1]) std::swap(tmin_y, tmax_y);if (dirIsNeg[2]) std::swap(tmin_z, tmax_z);float t_enter = std::max(tmin_x, std::max(tmin_y, tmin_z));float t_exit = std::min(tmax_x, std::min(tmax_y, tmax_z));if (t_enter < t_exit && t_exit >= 0) return true;else return false;
}

getIntersection()in BVH.cpp

这里参考PPT中的伪代码即可轻松完成：

代码如下：

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{// TODO Traverse the BVH to find intersectionstd::array<int, 3> dirIsNeg = {int(ray.direction.x<0), int(ray.direction.y<0), int(ray.direction.z<0)};bool flag_insert = node->bounds.IntersectP(ray, ray.direction_inv, dirIsNeg);Intersection intersection; // default: happened = false// if (ray misses node.bbox) return;if (!flag_insert) return intersection;// if (node is a leaf node)if (node->left == nullptr && node->right == nullptr){// test intersection with all objs;// return closest intersection;return node->object->getIntersection(ray);}// hit1 = Intersect(ray, node.child1);Intersection hit1 = getIntersection(node->left, ray);// hit2 = Intersect(ray, node.child2);Intersection hit2 = getIntersection(node->right, ray);// return the closer of hit1, hit2return (hit1.distance < hit2.distance ? hit1 : hit2);
}

渲染结果如下：

BVH用时：4秒钟

附加题部分：SAH加速的实现

这里主要参考了如下资料：

闫老师给的CMU PPT链接（完整PPT可供下载）

中文讲解：BVH with SAH (Bounding Volume Hierarchy with Surface Area Heuristic) - lookof - 博客园

这里主要讲几个重点，其中，成本cost的计算方式为：

伪代码：

中文讲解为：

将以上步骤简化：

cost(A,B) = p(A)*(n in A) + p(B)*(m in B)

其中：

p(A) = s(A) / s(C)

p(A) = s(B) / s(C)

因为分母同为s(C)，所以我们对比的时候无需对比分母，故cost(A,B)可以简化为：

cost(A,B) = s(A)*(n in A) + s(B)*(m in B)

分别在xyz轴中寻找上式的最小值即可。

代码实现：

首先，往BVH类的声明（BVH.hpp）中加入以下函数：

BVHBuildNode* BVHAccel::recursiveSAHBuild(std::vector<Object*> objects);

然后，在BVH.cpp中：

1）声明一个结构体，用于记录cost（划分成本的计算结果）和index（在何处划分）

// my struct
struct Cost
{int index;double cost_value;
};

2）声明一个函数compute_cost，用于计算xyz轴中按照任意一个轴进行划分时的最小成本，返回一个Cost类型的变量（划分方案&最小成本）

// my func
// 参数解释：
// objects为未划分的primitives（三角形）列表
// b为分成几个buckets
// p为一个bucket中有几个primitives（即三角形）
// dim为维度，0为x轴，1为y轴，2为z轴
Cost compute_cost(std::vector<Object*> objects, int b, int p, int dim)
{// 先将xyz轴按照升序排列switch (dim){case 0: // for x axisstd::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {return f1->getBounds().Centroid().x <f2->getBounds().Centroid().x;});break;case 1: // for y axisstd::sort(objects.begin(), objects.end(), [](auto f1, auto f2){return f1->getBounds().Centroid().y <f2->getBounds().Centroid().y;});break;case 2: // for z axisstd::sort(objects.begin(), objects.end(), [](auto f1, auto f2){return f1->getBounds().Centroid().z <f2->getBounds().Centroid().z;});break;}// 记录cost的列表std::vector<Cost> cost_list;// 循环计算不同划分方案的成本Bounds3 bucketA;// 将buckets从小到大逐个加入划分方案for (int i = 0; i < b-1; ++i) {for (int j = 0; j < p; ++j) {// 将bucket中的元素从小到大加入包围盒A中bucketA = Union(bucketA, objects[p*i+j]->getBounds());}// 将剩余的元素从小到大加入包围盒B中Bounds3 bucketB;for (int k = p*(i+1); k < objects.size(); ++k) {bucketB = Union(bucketB, objects[k]->getBounds());}// 计算成本double sA = bucketA.SurfaceArea(); // 获取包围盒A的表面积double sB = bucketB.SurfaceArea();double nA = (i+1)*p; // 包围盒A中的primitives（三角形）的个数double nB = objects.size() - nA;double cost_val = sA*nA + sB*nB; // compute costCost cost;cost.index = p*(i+1)-1;cost.cost_value = cost_val;cost_list.emplace_back(cost);}// 按照升序排列（第一个为cost的最小值）std::sort(cost_list.begin(), cost_list.end(), [](auto f1, auto f2){return f1.cost_value < f2.cost_value;});// 返回cost最小值（和对应的划分方案index）return cost_list[0];
}

3）声明一个函数get_partition_axis，对比xyz轴的各自最小划分成本，方便后续选择划分成本最小的轴进行划分：

// my func
Cost get_partition_axis(Cost cost_x, Cost cost_y, Cost cost_z, int& dim)
{if (cost_x.cost_value < cost_y.cost_value && cost_x.cost_value < cost_z.cost_value) {dim = 0;return cost_x;} else if (cost_y.cost_value < cost_z.cost_value) {dim = 1;return cost_y;} else {dim = 2;return cost_z;}
}

4）最后是BVHAccel::recursiveSAHBuild函数的实现：

BVHBuildNode* BVHAccel::recursiveSAHBuild(std::vector<Object*> objects)
{BVHBuildNode* node = new BVHBuildNode();// num of primitives in a bucket// cannot be smaller than 2!int p = 5;if (objects.size() <= p) {node = recursiveBuild(objects);return node;}int b = objects.size() / p; // num of bucketsif (objects.size() % p != 0) b++;Cost cost_x;Cost cost_y;Cost cost_z;Cost final_cost;int dim = 0; // 默认按照x轴划分，防bugcost_x = compute_cost(objects, b, p, 0);cost_y = compute_cost(objects, b, p, 1);cost_z = compute_cost(objects, b, p, 2);final_cost = get_partition_axis(cost_x, cost_y, cost_z, dim);switch (dim) {case 0:std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {return f1->getBounds().Centroid().x <f2->getBounds().Centroid().x;});break;case 1:std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {return f1->getBounds().Centroid().y <f2->getBounds().Centroid().y;});break;case 2:std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {return f1->getBounds().Centroid().z <f2->getBounds().Centroid().z;});break;}auto beginning = objects.begin();auto middling = objects.begin() + final_cost.index; // SAH changedauto ending = objects.end();auto leftshapes = std::vector<Object*>(beginning, middling);auto rightshapes = std::vector<Object*>(middling, ending);assert(objects.size() == (leftshapes.size() + rightshapes.size()));node->left = recursiveSAHBuild(leftshapes);node->right = recursiveSAHBuild(rightshapes);node->bounds = Union(node->left->bounds, node->right->bounds);return node;
}

需要运行SAH的时候，记得需要修改如下两处：

第一处，修改默认splitMethod的值：

// BVHAccel Declarations
inline int leafNodes, totalLeafNodes, totalPrimitives, interiorNodes;
class BVHAccel {public:// BVHAccel Public Typesenum class SplitMethod { NAIVE, SAH };// 从这里开始// BVHAccel Public Methods// BVHAccel(std::vector<Object*> p, int maxPrimsInNode = 1, SplitMethod splitMethod = SplitMethod::NAIVE);// 把上面注释掉的那一句，改成下面这一句BVHAccel(std::vector<Object*> p, int maxPrimsInNode = 1, SplitMethod splitMethod = SplitMethod::SAH);// 到这里为止Bounds3 WorldBound() const;~BVHAccel();Intersection Intersect(const Ray &ray) const;Intersection getIntersection(BVHBuildNode* node, const Ray& ray)const;bool IntersectP(const Ray &ray) const;BVHBuildNode* root;// BVHAccel Private MethodsBVHBuildNode* recursiveBuild(std::vector<Object*>objects);BVHBuildNode* recursiveSAHBuild(std::vector<Object*>objects); // my func// BVHAccel Private Dataconst int maxPrimsInNode;const SplitMethod splitMethod;std::vector<Object*> primitives;
};

第二处，在BVHAccel::BVHAccel() （BVH.cpp）中：

BVHAccel::BVHAccel(std::vector<Object*> p, int maxPrimsInNode, SplitMethod splitMethod): maxPrimsInNode(std::min(255, maxPrimsInNode)), splitMethod(splitMethod),primitives(std::move(p))
{time_t start, stop;time(&start);if (primitives.empty())return;// 从这里开始if (splitMethod == SplitMethod::NAIVE) {root = recursiveBuild(primitives);} else if (splitMethod == SplitMethod::SAH) {root = recursiveSAHBuild(primitives); // for SAH}// 到这里结束time(&stop);double diff = difftime(stop, start);int hrs = (int)diff / 3600;int mins = ((int)diff / 60) - (hrs * 60);int secs = (int)diff - (hrs * 3600) - (mins * 60);printf("\rBVH Generation complete: \nTime Taken: %i hrs, %i mins, %i secs\n\n",hrs, mins, secs);
}

然后大功告成，可以编译运行啦！

SAH运行之后的结果：耗时3秒（比BVH的4秒少1秒）

【Games101 作业6 + 附加题】渲染兔子 BVH SAH 代码相关推荐

【Games101 作业5】光线追踪渲染小球
本次作业中,有2个函数需要修改: 其中第二个比较简单,我们先讲第二个. rayTriangleIntersect():在Triangle.hpp中按照PPT公式来实现即可. 需要注意的点是:公式中的 ...
RT-thread 柿饼UI demo(文本浏览+电子相册) ---- 暨柿饼入门课第一周作业附加题
一.题目要求完成下图的应用制作,范进中举复制群文件内的范进中举.txt内的内容. 二.实现过程 2.1.整体思路整体的框架是在一个page中放入三个button控件(负责控制三个显示界面的切换)和 ...
SLAM十四讲第一讲作业附加题：ORB-SLAM2
SLAM十四讲第一讲大作业附加题:ORB-SLAM2 (供参考)@[素履以往] 借鉴 https://blog.csdn.net/weixin_37661692/article/details/8 ...
GAMES101作业7-路径追踪实现过程代码框架超全解读
目录 Path Tracing算法过程讨论蒙特卡洛积分直接光照 direct illumination 间接光照 indirect illumination 编辑合成全局光照解决一些存在的问 ...
2021年人工神经网络第四次作业-第五题：危险品识别
简介: 通过对于物品X射线数据集合的整理,挑选出15类体积比较大的物品,训练LeNet网络进行识别. 关键词: X射线,危险品识别,LeNet,Paddle #mermaid-svg-wZUMACG ...
网页附加题写出下图的html,附加题（写HTML文件）：根据给定的博客名单，自动生成HTML网页...
收集了学生CSDN博客地址很久了,但一直没来得及整理成贺利坚老师的完美班级网页名册.今天突然想,一共有6个班学生,如果手动写的话,太费事了.我们程序员,就是让费事不费脑的工作自动化,即使是第一次花很多 ...
GAMES101作业5-从头到尾理解代码Whitted光线追踪
目录 Whitted Ray-Tracing Whitted光线追踪 What Why How 1 发射主射线primary ray 实现步骤 (1)定义相机 (2)计算primary主射线的方向 R ...
python123第一周作业答案程序题_[python爬虫]第一周作业_顾静
第一次作业第一题 a = 10,b = 3 计算下面c的值及输出数据类型 1.c = a/b - a 2.c = a/b * a 3.c = 0.1 * a//b - a 4.c = a//b + ...
GAMES101作业7提高-实现微表面模型你需要了解的知识
目录微表面材质模型微平面理论 Microfacet Theory BSDF(浅浅的提一下) 微表面BRDF的实现 Cook-Torrance BRDF 漫反射的BRDF 镜面反射的BRDF 1 法 ...

【Games101 作业6 + 附加题】渲染兔子 BVH SAH 代码