【GAMES101】作业6 加速结构

一.作业描述
二.作业解析
- **1. 递归构造BVH树分析**
- **2. Render函数**
- **3. Triangle::getItersection函数**
- **4. 射线与AABB包围盒相交判断**
- **5. BVH::getItersection函数**
- **6. SAH优化思路**

一.作业描述

二.作业解析

1. 递归构造BVH树分析

看视频学习后过了太久，BVH树怎么构造已经忘记了，这里重新根据作业提供的代码整理了一下思路，重新学习了。后面要递归BVH树去判断相交，所以了解一下BVH树的结构也是很必要的。

BVH树节点node属性：

Bound：包围盒；
left，right：左右指针；
object：Object类；

递归实现：
1.递归出口：
当场景中物体总数=1时，包围盒为当前物体的包围盒，左右节点为NULL，返回当前节点；

2.当场景中物体总数=2时，取两个物体的边界并集作为当前节点的包围盒，左右节点分别再进行一次递归；

3.一般递归思路：
取当前所有物体的包围盒并集作为当前节点的总包围盒大小，求出总包围盒大小三维（x，y，z）中长度最大的一个维度dim，并将所有物体的质心坐标在维度dim上的值进行排序，将排序后的前1/2物体和后1/2的物体分别作为左右子集继续进行递归。

归纳一下，BVH树叶节点中存储比较重要的信息就是包围盒bounds。
在遍历BVH树判断相交时，如果发现当前节点的左右指针均为空，就说明当前节点的包围盒大小已经不可再分了。

递归构造的代码在这里mark一下。

BVHBuildNode* BVHAccel::recursiveBuild(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;}else {Bounds3 centroidBounds;for (int i = 0; i < objects.size(); ++i)centroidBounds =Union(centroidBounds, objects[i]->getBounds().Centroid());int dim = centroidBounds.maxExtent();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() + (objects.size() / 2);auto 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 = recursiveBuild(leftshapes);node->right = recursiveBuild(rightshapes);node->bounds = Union(node->left->bounds, node->right->bounds);}return node;
}

2. Render函数

Render仅有些数据结构上的简单变动。
作业代码实现：

//Render函数：
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 directionfloat x = (2 * (i + 0.5) / (float)scene.width - 1) *imageAspectRatio * scale;float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;// 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*// Don't forget to normalize this direction!Vector3f dir = normalize(Vector3f(x, y, -1)); Ray ray=Ray(eye_pos,dir,0);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);
}

3. Triangle::getItersection函数

求三角形与射线是否相交的getItersection函数的基础判断部分已经实现了，相比作业5返回值出现了一点变化，之前是返回true，这次时返回结构体Intersection类型的inter，return前给各个属性赋上计算出的值即可。

Intersection属性:

happened：是否有交点；
coords：交点坐标；
normal：交点处平面法向量；
obj：相交物体；
m:相交点面材质；

如果判断相交，则修改inter中的属性即可。
如下公式所示，根据射线ray的定义，t代表的是射线初始点到交点的距离。
计算出t的值以后，可根据origin+t*dir求出相交点坐标。

依次为inter中的属性赋值，具体代码如下。

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;// TODO find ray triangle intersectioninter.happened=true;inter.coords=ray(t_tmp);inter.normal=normal;inter.obj=this;inter.distance=t_tmp;inter.m=this->m;return inter;
}

4. 射线与AABB包围盒相交判断

二维判断方法如课程中所讲，最终结果算出近端tnear取从x，y两个维度txmin和tymin中的最大值，tfar取x，y两个维度txmax和tymax中的最小值。
三维判断与二维判断思路一致，只是多了一个维度。

概括一下判断相交思路：
如果texit<0，那么包围盒在射线”后面“，没有相交；
如果tenter<texit且texit>=0，射线与包围盒相交；

实现代码：

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 intersectsauto Vmin=(pMin-ray.origin)*invDir;auto Vmax=(pMax-ray.origin)*invDir;if(dirIsNeg[0]==1)std::swap(Vmin.x,Vmax.x);if(dirIsNeg[1]==1)std::swap(Vmin.y,Vmax.y);if(dirIsNeg[2]==1)std::swap(Vmin.z,Vmax.z);float tmin=fmax(Vmin.x,fmax(Vmin.y,Vmin.z));float tmax=fmin(Vmax.x,fmin(Vmax.y,Vmax.z));if(tmin<tmax&&tmax>=0)return true;else return false;
}

5. BVH::getItersection函数

加速后的求交的判断步骤为：
深度优先遍历BVH树，判断射线是否与当前节点的轴对称包围盒AABB相交，如果不相交，则直接结束；如果相交，判断当前节点的包围盒是否可以继续划分（即当前节点是否有左右子节点），若可以继续划分，返回左右求交获得的intersection中最近的intersection。

代码如下：

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{// TODO Traverse the BVH to find intersectionIntersection inter,interL,interR;std::array<int, 3> dirIsNeg;if(ray.direction.x<0)dirIsNeg[0]=1;if(ray.direction.y<0)dirIsNeg[1]=1;if(ray.direction.z<0)dirIsNeg[2]=1;bool flag=node->bounds.IntersectP(ray,ray.direction_inv,dirIsNeg);if(!flag)return inter;if(node->left==nullptr&&node->right==nullptr){return inter=node->object->getIntersection(ray);}interL=getIntersection(node->left,ray);interR=getIntersection(node->right,ray);return interL.distance<interR.distance?interL:interR;}

6. SAH优化思路

BVH存在着包围盒重叠的问题，然而划分时要尽可能减少划分后两部分包围盒重叠的体积，因为重叠的体积越大，射线穿过重叠区域的可能性越大，遍历两个子树的可能性就越高，计算消耗越多。

此外，如果空间中物体分布得很不均匀，比如右边存在大量物体，而左边只有少量物体，此时使用BVH就会得到很差的划分结果。

SAH（Surface Area Heuristic）基于表面积对划分进行评估，这种方法通过对求交代价和遍历代价进行评估，给出了每一种划分的代价（Cost），而我们的目的便是去寻找代价最小的划分。
内容参考：PBRT-E4.3-层次包围体(BVH)
↓↓↓↓↓↓↓↓↓↓

c(A,B)为当前划分的代价；
p(A),p(B)为射线落在划分的两个子包围盒的概率；
t(i),t(j)为对每个物体求交的代价，t_trav为遍历代价。

SAH考虑了物体在空间中的分布，也考虑到了划分后两部分包围盒重叠的体积问题，优化了性能。

这里就不放具体的代码了，我尝试了：
B站 SAH和知乎 SAH
这两种方法大差不差，处理细节上有些不同，但是我的虚拟机处处透露着离谱的气息，单BVH的render都用了15s，别人描述的都是6s，SAH优化以后时间还多了1s（实属反向优化了），后面越来越慢…就当是先对SAH做了个简单的了解吧。

实现效果：