games101 作业6

题目

在之前的编程练习中，我们实现了基础的光线追踪算法，具体而言是光线传输、光线与三角形求交。我们采用了这样的方法寻找光线与场景的交点：遍历场景中的所有物体，判断光线是否与它相交。在场景中的物体数量不大时，该做法可以取得良好的结果，但当物体数量增多、模型变得更加复杂，该做法将会变得非常低效。因此，我们需要加速结构来加速求交过程。在本次练习中，我们重点关注物体划分算法 Bounding Volume Hierarchy (BVH)。本练习要求你实现 Ray-Bounding
Volume 求交与 BVH 查找。首先，你需要从上一次编程练习中引用以下函数：• Render() in Renderer.cpp: 将你的光线生成过程粘贴到此处，并且按照新框架更新相应调用的格式。
• Triangle::getIntersection in Triangle.hpp: 将你的光线-三角形相交函数粘贴到此处，并且按照新框架更新相应相交信息的格式。

在本次编程练习中，你需要实现以下函数：

//Bounds3.hpp
IntersectP(const Ray& ray, const Vector3f& invDir,const std::array<int, 3>& dirIsNeg) ;

说明：这个函数的作用是判断包围盒 BoundingBox 与光线是否相交，你需要按照课程介绍的算法实现求交过程。

//BVH.cpp
getIntersection(BVHBuildNode* node, const Ray ray);

说明： 建立 BVH 之后，我们可以用它加速求交过程。该过程递归进行，你将在其中调
用你实现的 Bounds3::IntersectP

代码框架解读和知识点讲解

为了能够更容易理解作业的答案，我们先对代码框架整体进行一个讲解以及相关知识点的剖析。

作业框架

整个流程如下：

1. 创建并初始化Scene。
2. 根据obj文件构建MeshTriangle；（MeshTriangle是Object的子类，保存三角形mesh,并在内部根据obj文件中的数据构建BVH，同时支持光线求交，返回包围盒等）
3. Scene中添加物体MeshTriangle和光源Light
4. Scene中构建BVH。(Scene的BVH中叶子节点中保存的是物体(MeshTriangle), 而MeshTriangle中的BVH叶子节点中保存的是三角形,物体和光线求交时会调用物体的getIntersection，而getIntersection内部又会调用物体的BVH，最终递归调用到三角形和光线求交。因为MeshTriangle和Triangle都集成自Object，均有getIntersection，而Triangle的getIntersection 方法会最后返回三角形和光线求交的结果)。

    Scene scene(1280, 960);const auto path = std::string(MODELS_PATH) + "/bunny/bunny1.obj";MeshTriangle bunny(path);scene.Add(&bunny);scene.Add(std::make_unique<Light>(Vector3f(-20, 70, 20), 1));scene.Add(std::make_unique<Light>(Vector3f(20, 70, 20), 1));scene.buildBVH();

6. 开始渲染，遍历framebuffer中的每个像素点，为每个像素点生成 射线Ray

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 i = 0; i < scene.width; ++i) {// generate primary ray directionfloat x = (2 * (i + 0.5) / (float)scene.width - 1) *imageAspectRatio * scale;float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;Vector3f dir = Vector3f(x, y, -1); // Don't forget to normalize this direction!dir = normalize(dir);Ray ray(eye_pos,dir,0);framebuffer[m++] = scene.castRay(ray, 0);}UpdateProgress(j / (float)scene.height);}} 

7. 射线和Scene求交，即和Scene中的每个物体求交，进一步和物体中的每个三角形求交，得到Intersection，Intersection保存了是否相交(happened),交点的坐标(coords)，射线的时间t(distance),以及相交的三角形(obj)、材质(m)和法向量(normal).
   注意：这里的bvh->Intersect 将是本次作业的核心。

Intersection Scene::intersect(const Ray &ray) const{return this->bvh->Intersect(ray);
}
Vector3f Scene::castRay(const Ray &ray, int depth) const{if (depth > this->maxDepth) {return Vector3f(0.0,0.0,0.0);}Intersection intersection = Scene::intersect(ray);...
}

8. 根据求交结果，如果相交，则根据相交点的材质，纹理坐标、法线等信息计算像素点的颜色值。材质的不同，计算的方式不同，本次作业默认使用材质类型为：DIFFUSE_AND_GLOSSY, 所以采用Phong 着色方式计算。其他两种REFLECTION_AND_REFRACTION, REFLECTION在后续作业中使用。

Vector3f Scene::castRay(const Ray &ray, int depth) const
{if (depth > this->maxDepth) {return Vector3f(0.0,0.0,0.0);}Intersection intersection = Scene::intersect(ray);Material *m = intersection.m;Object *hitObject = intersection.obj;Vector3f hitColor = this->backgroundColor;
//    float tnear = kInfinity;Vector2f uv;uint32_t index = 0;if(intersection.happened) {Vector3f hitPoint = intersection.coords;Vector3f N = intersection.normal; // normalVector2f st; // st coordinateshitObject->getSurfaceProperties(hitPoint, ray.direction, index, uv, N, st);switch (m->getType()) {case REFLECTION_AND_REFRACTION:{// 递归处理反射 折射...}case REFLECTION:{// 递归处理反射...}default:{// Phong 着色...}

知识点剖析

1. 如何表示一条射线？

$$r(t)=o+t\ast d(0<=t<\infty )$$

其中$r(t)$表示射线上任意一点，$o$表示射线的起点，$d$ 表示射线的方向，$t$是标量，可以认为是射线从起点出发的时间长度，或者看作某一点到起点的远近程度。
代码表示如下：

struct Ray{//Destination = origin + t*directionVector3f origin;Vector3f direction, direction_inv; // 保存方向的倒数是为了加速运算double t;//transportation time,double t_min, t_max;
}

2. 如何表示一个三维包围盒？

使用两个点表示，包围盒中所有物体在三个方向的最小值和最大值
$$P_{min}=(x_{min},y_{min},z_{min}),P_{max}=(x_{max},y_{max},z_{max})$$
代码表示如下：

class Bounds3
{
public:Vector3f pMin, pMax; Vector3f Diagonal();// 对角线向量int maxExtent();// 最长轴索引double SurfaceArea();// 表面积Vector3f Centroid();//中心点Bounds3 Intersect(const Bounds3& b);// 交集inline Bounds3 Union(const Bounds3& b1, const Bounds3& b2)// 两个包围盒求并集inline Bounds3 Union(const Bounds3& b, const Vector3f& p)// 点和包围盒求并集Vector3f Offset(const Vector3f& p); // p点相对包围盒的pMin的归一化值bool Overlaps(const Bounds3& b1, const Bounds3& b2);// 是否有交集bool Inside(const Vector3f& p, const Bounds3& b);// 判断一个点是否在包围盒内部(包含边界)bool IntersectP(const Ray& ray, const Vector3f& invDir,const std::array<int, 3>& dirisNeg); // 判断光线是否和包围盒相交  dirIsNeg表示ray的每个方向是否为正方向
}

3. 包围盒如何和一条射线求交？

(1). 计算每个轴的交点参数t
根据射线方程$r(t)=o+t\ast d(0<=t<\infty )$，对于x轴，和最小面$x=pMin.x$的交点
$$P_{minX}=o+t_{minX}\ast d$$
由于$P_{minX}.x=o.x+t_{minX}\ast d.x=pMin.x$,所以$$t_{minX}=(pMin.x-o.x)/d.x$$
同理可得
$$\begin{gathered} t_{minX}=(pMin.x-o.x)/d.x \\ {t}_{maxX}=(pMax.x-o.x)/d.x \\ {t}_{minY}=(pMin.y-o.y)/d.y \\ {t}_{maxY}=(pMax.y-o.y)/d.y \\ {t}_{minZ}=(pMin.z-o.z)/d.z \\ {t}_{maxZ}=(pMax.z-o.z)/d.z \end{gathered}$$
(2) 处理射线方向
如果射线方向为负，则$t_{minX}$ 和$t_{maxX}$ 的意义反了，需要交换，保证$t_{minX}$ 始终是靠近的面，$t_{maxX}$是远离的面。y、z 轴同理。
(3) 计算整体进入和离开参数$t_{enter},t_{exit}$。
光线进入即光线进入到包围盒轴平面的交集区，也就是最晚到达靠近射线起点的包围盒轴平面，即$t_{enter}=max(t_{minX}，t_{minY}，t_{minZ})$,光线离开即光线离开到包围盒轴平面的交集区，也就是最早离开远离射线起点的包围盒轴平面。即$t_{exit}=max(t_{maxX}，t_{maxY}，t_{maxZ})$。进入需要考虑三个方向全部进入，而离开只要有一个方向离开即可。如PPT中所示的2D示例，3D同理。

(4) 判断是否有交
$t_{enter}<t_{exit}$：三轴区间有重叠，说明射线穿过包围盒
$t_{enter}>=ray.t_{min}且t_{exit}<=ray.t_{max}$：交点在射线的有效区间内

注意

• 射线平行于某轴时，需判断起点是否在包围盒该轴的范围内。
• 若不在，则必不相交。
• 若在，则该轴对 t 区间无影响，$t_{min}=-\infty ,t_{max}=+\infty$。

4. 如何递归构造BVH

可参考博客
什么是BVH?
包围体层次结构(BVH)是一个包含包围盒和几何对象的树结构。每个节点都有一个包围盒。所有子节点的包围盒都包含在父节点的包围盒中。只有子节点持有具体的几何对象。如下图所示。

注意：作业中为了简化处理，叶子节点只能包含一个物体。
构造BVH的方法有很多种！常见的方式有：中点划分，等量划分、SAH(也是作业中的提升部分)等。
作业中提供了最简单的实现方式。
(1) 中点划分
下面采用源码+注释的方式讲解原理

BVHBuildNode* BVHAccel::recursiveBuild(std::vector<Object*> objects)
{BVHBuildNode* node = new BVHBuildNode();// 计算所有objects的包围盒Bounds3 bounds;for (int i = 0; i < objects.size(); ++i)bounds = Union(bounds, objects[i]->getBounds());// 如果包围盒中的物体只有一个，则更新叶子节点的包围盒，以及物体，并退出递归。if (objects.size() == 1) {// Create leaf _BVHBuildNode_node->bounds = objects[0]->getBounds();node->object = objects[0];node->left = nullptr;node->right = nullptr;return node;}// 如果包围盒中的物体只有两个，则分别递归更新左子树和右子树，同时更新包围盒，由于不是叶子节点，不需要更新object。else if (objects.size() == 2) {node->left = recursiveBuild(std::vector{objects[0]});node->right = recursiveBuild(std::vector{objects[1]});node->bounds = Union(node->left->bounds, node->right->bounds);return node;}// 否则，从分割轴对所有物体进行排序，然后二分，一半物体作为node的左子树，一半作为右子树。// 分割轴指：所有物体的中心点组成的包围盒的最长轴。else {// 计算所有物体的中心点组成的包围盒Bounds3 centroidBounds;for (int i = 0; i < objects.size(); ++i)centroidBounds =Union(centroidBounds, objects[i]->getBounds().Centroid());// 获得包围盒的最长轴作为分割轴。int dim = centroidBounds.maxExtent();// 根据分割轴，对所有物体进行某一个方向上的位置排序，并存放在objects容器中。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;}// 对排序后的物体进行分割，找到中间物体middling auto beginning = objects.begin();auto middling = objects.begin() + (objects.size() / 2);auto ending = objects.end();// 左子树存储beginning 到middling (不包含middling )auto leftshapes = std::vector<Object*>(beginning, middling);// 左子树存储middling 到ending auto rightshapes = std::vector<Object*>(middling, ending);assert(objects.size() == (leftshapes.size() + rightshapes.size()));// 递归构建左子树和右子树node->left = recursiveBuild(leftshapes);node->right = recursiveBuild(rightshapes);// 更新当前节点的包围盒：即左子树和右子树包围盒的并集。node->bounds = Union(node->left->bounds, node->right->bounds);}return node;
}

(2)SAH 表面积启发式算法（Surface Area Heuristic）
参考文章
构建BVH 其实就三步：
1.计算所有物体的质心边界。（中心点组成的包围盒）
2.选择一个轴（x,y 或者z轴）。
3.根据选择的轴，和一定的策略，对所有物体分为两部分，分别作为左节点和右节点。
然后递归处理。

SAH和上述的中点划分方式不同在于 选择轴的方式不同，划分物体的策略不同。
SAH核心点：
1.代价公式：
$$cost=Ct+(leftArea/totalArea)\ast leftCount\ast Ci+(rightArea/totalArea)\ast rightCount\ast Ci$$
其中
Ct: 经过该节点的遍历代价（通常设为1）
Ci: 叶子节点交点测试代价（通常设为1）
leftArea/rightArea: 左/右子包围盒表面积
leftCount/rightCount: 左/右子物体数
totalArea： 所有物体包围盒总面积
2.原理
1.首先按照每个轴，对所有物体进行分桶（也就是将所有物体分配到$bucketnum$个桶中，每个桶都记录了存放了多少个物体以及这些物体的包围盒）
分桶：就是根据坐标的某个分量，比如x，计算该物体应该在哪个桶中，计算如下：
$$index=(cente{r}_x-mi{n}_x)/(ma{x}_x-mi{n}_x)\ast bucketnum$$
其中 $cente{r}_x$为该物体的包围盒中心点的x坐标。
2.计算每个分割点的代价，找到最小代价的分割（分割点就是桶的索引）
代码说明可能比较清楚。
$cost(j)=cost[0,j]+cost[j+1,bucketnum],j=[1,bucketnum]$

// 计算每个分割点的SAH代价
for (int i = 1; i < nBuckets; i++)
{BucktInfo leftBucket, rightBucket;for (int j = 0; j < i; j++) {leftBucket.count += buckets[j].count;leftBucket.bounds = Union(leftBucket.bounds, buckets[j].bounds);}for (int j = i; j < nBuckets; j++) {rightBucket.count += buckets[j].count;rightBucket.bounds = Union(rightBucket.bounds, buckets[j].bounds);}if (leftBucket.count == 0 || rightBucket.count == 0)continue;float leftArea = leftBucket.bounds.SurfaceArea();float rightArea = rightBucket.bounds.SurfaceArea();float totalArea = bounds.SurfaceArea();// SAH costfloat cost = 1 + (leftArea / totalArea) * leftBucket.count + (rightArea / totalArea) * rightBucket.count;if (cost < minCost) {minCost = cost;bestSplitIndex = i;bestSplitAxis = axis;}
}

3. 遍历每个轴，即可得到代价的最小的轴以及分割点。
4. 根据代价的最小的轴以及分割点，对所有物体进行划分。（可以具体看代码）

题解

本次作业主要是练习光线和场景中的物体包围盒求交，利用BVH 加速求交的过程。

1.IntersectP 光线和包围盒求交

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// dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x>0),int(y>0),int(z>0)], use this to simplify your logic// TODO test if ray bound intersects// Calculate intersection t values for x, y, z planes  float tMinX = (pMin.x - ray.origin.x) * invDir.x;float tMaxX = (pMax.x - ray.origin.x) * invDir.x;if (!dirIsNeg[0]) std::swap(tMinX, tMaxX);float tMinY = (pMin.y - ray.origin.y) * invDir.y;float tMaxY = (pMax.y - ray.origin.y) * invDir.y;if (!dirIsNeg[1]) std::swap(tMinY, tMaxY);float tMinZ = (pMin.z - ray.origin.z) * invDir.z;float tMaxZ = (pMax.z - ray.origin.z) * invDir.z;if (!dirIsNeg[2]) std::swap(tMinZ, tMaxZ);// Find the largest tMin and smallest tMax  float tEnter = std::max(tMinX, std::max(tMinY, tMinZ));float tExit = std::min(tMaxX, std::min(tMaxY, tMaxZ));// Check if the ray intersects the bounding box  bool res = (tEnter <= tExit) && (tExit >= ray.t_min) && (tEnter <= ray.t_max);return res;
}

2 BVHAccel::getIntersection 光线和BVH的某个节点求交，并返回交点信息

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{// TODO Traverse the BVH to find intersectionIntersection isect;// Check if the ray intersects the bounding box of the node  if (!node->bounds.IntersectP(ray, ray.direction_inv, { ray.direction.x > 0, ray.direction.y > 0, ray.direction.z > 0 }))return isect;// If the node is a leaf, check intersection with the object  if (node->left == nullptr && node->right == nullptr) {return node->object->getIntersection(ray);}// Otherwise, traverse the left and right children  Intersection leftIsect = getIntersection(node->left, ray);Intersection rightIsect = getIntersection(node->right, ray);// Return the closer intersection  if (leftIsect.happened && (!rightIsect.happened || leftIsect.distance < rightIsect.distance))return leftIsect;return rightIsect;}

3 SAH 构建BVH

BVHBuildNode* BVHAccel::recursiveBuildSAH(std::vector<Object*> objects)
{BVHBuildNode* node = new BVHBuildNode();// Compute bounds of all primitives in BVH nodeBounds3 bounds;for (int i = 0; i < objects.size(); ++i)bounds = Union(bounds, objects[i]->getBounds());if (objects.size() == 1) {// Create leaf _BVHBuildNode_node->bounds = objects[0]->getBounds();node->object = objects[0];node->left = nullptr;node->right = nullptr;return node;}else if (objects.size() == 2) {node->left = recursiveBuild(std::vector{ objects[0] });node->right = recursiveBuild(std::vector{ objects[1] });node->bounds = Union(node->left->bounds, node->right->bounds);return node;}// 计算所有物体的质心边界Bounds3 centroidBounds;for (int i = 0; i < objects.size(); ++i)centroidBounds =  Union(centroidBounds, objects[i]->getBounds().Centroid());// 选择分割轴int dim = centroidBounds.maxExtent();// Initialize buckets for SAH //•	SAH 代价公式：cost = Ct + (leftArea/totalArea) * leftCount * Ci + (rightArea/totalArea) * rightCount * Ci/*Ct: 经过该节点的遍历代价（通常设为1）Ci: 叶子节点交点测试代价（通常设为1）leftArea/rightArea: 左/右子包围盒表面积leftCount/rightCount: 左/右子物体数*/const int nBuckets = 12; // Number of buckets for SAHstruct BucktInfo{int count = 0;Bounds3 bounds;};float minCost = std::numeric_limits<float>::infinity();int bestSplitIndex = -1;int bestSplitAxis = dim;if (objects.size() >4) {for (int axis = 0; axis < 3; axis++){std::vector<BucktInfo> buckets(nBuckets);float mincoord = centroidBounds.pMin[axis];float maxcoord = centroidBounds.pMax[axis];float extent = maxcoord - mincoord + 1e-5f;// 分桶for (int i = 0; i < objects.size(); ++i) {float centroid = objects[i]->getBounds().Centroid()[axis];int bucketIndex = std::min(nBuckets - 1, static_cast<int>((centroid - mincoord) / extent * nBuckets));buckets[bucketIndex].count++;buckets[bucketIndex].bounds = Union(buckets[bucketIndex].bounds, objects[i]->getBounds());}// 计算每个分割点的SAH代价for (int i = 1; i < nBuckets; i++){BucktInfo leftBucket, rightBucket;for (int j = 0; j < i; j++) {leftBucket.count += buckets[j].count;leftBucket.bounds = Union(leftBucket.bounds, buckets[j].bounds);}for (int j = i; j < nBuckets; j++) {rightBucket.count += buckets[j].count;rightBucket.bounds = Union(rightBucket.bounds, buckets[j].bounds);}if (leftBucket.count == 0 || rightBucket.count == 0)continue;float leftArea = leftBucket.bounds.SurfaceArea();float rightArea = rightBucket.bounds.SurfaceArea();float totalArea = bounds.SurfaceArea();// SAH costfloat cost = 1 + (leftArea / totalArea) * leftBucket.count + (rightArea / totalArea) * rightBucket.count;if (cost < minCost) {minCost = cost;bestSplitIndex = i;bestSplitAxis = axis;}}}}std::vector<Object*> leftshapes, rightshapes;if (bestSplitIndex != -1){// If a good split was found, use itint dim = bestSplitAxis;float mincoord = centroidBounds.pMin[dim];float maxcoord = centroidBounds.pMax[dim];float extent = maxcoord - mincoord + 1e-5f;for (auto& obj : objects) {float centroid = obj->getBounds().Centroid()[dim];int bucketIndex = std::min(nBuckets - 1, static_cast<int>((centroid - mincoord) / extent * nBuckets));if (bucketIndex < bestSplitIndex)leftshapes.push_back(obj);elserightshapes.push_back(obj);}// 防止极端情况，leftshapes 或者 rightshapes 为空if (leftshapes.empty() || rightshapes.empty()) {leftshapes = std::vector<Object*>(objects.begin(), objects.begin() + objects.size() / 2);rightshapes = std::vector<Object*>(objects.begin() + objects.size() / 2, objects.end());}}if(bestSplitIndex == -1 || leftshapes.empty() || rightshapes.empty()) {std::sort(objects.begin(), objects.end(), [bestSplitAxis](auto f1, auto f2) {return f1->getBounds().Centroid()[bestSplitAxis] <f2->getBounds().Centroid()[bestSplitAxis];});// If no good split was found, use the naive methodleftshapes = std::vector<Object*>(objects.begin(), objects.begin() + objects.size() / 2);rightshapes = std::vector<Object*>(objects.begin() + objects.size() / 2, objects.end());}assert(objects.size() == (leftshapes.size() + rightshapes.size()));node->left = recursiveBuild(leftshapes);node->right = recursiveBuild(rightshapes);node->bounds = Union(node->left->bounds, node->right->bounds);return node;
}

作业答案

本次作业的答案放在的git仓库中：作业地址

参考资料

【Unity Graphics】BVH（Bounding volume hierarchy）朴素的构建算法（2）
BVH with SAH (Bounding Volume Hierarchy with Surface Area Heuristic)
如何将Surface Area Heuristic(表面面积启发式)运用到BVH的构建中