games101 作业6
这是games101 现代计算机图形学 作业06 by coolwx
感谢https://blog.csdn.net/miyu1994/article/details/107006010/ 这位大佬的文章，因为我一开始真的写不出来SAH（启发式搜索），看了这位大神的算法，我最终才理解了他的算法，其实还是比较简单的，关于SAH，我应该会写一篇笔记，来记录我的感想，再次感谢以上这位大神的算法，因为我一开始真的想不出这个启发式算法，所以代码我会使用他的代码，然后添上我的注释，便于理解（哎，还是太菜，没写出来）。
作业要求：
在之前的编程练习中，我们实现了基础的光线追踪算法，具体而言是光线传
输、光线与三角形求交。我们采用了这样的方法寻找光线与场景的交点：遍历场景
中的所有物体，判断光线是否与它相交。在场景中的物体数量不大时，该做法可以
取得良好的结果，但当物体数量增多、模型变得更加复杂，该做法将会变得非常低
效。因此，我们需要加速结构来加速求交过程。在本次练习中，我们重点关注物体
划分算法 Bounding Volume Hierarchy (BVH)。本练习要求你实现 Ray-Bounding
Volume 求交与 BVH 查找。
首先，你需要从上一次编程练习中引用以下函数：
• Render() in Renderer.cpp: 将你的光线生成过程粘贴到此处，并且按照新框
架更新相应调用的格式。
• Triangle::getIntersection in Triangle.hpp: 将你的光线-三角形相交函数
粘贴到此处，并且按照新框架更新相应相交信息的格式。
在本次编程练习中，你需要实现以下函数：
• IntersectP(const Ray& ray, const Vector3f& invDir,
const std::array<int, 3>& dirIsNeg) in the Bounds3.hpp: 这个函数的
作用是判断包围盒 BoundingBox 与光线是否相交，你需要按照课程介绍的算
法实现求交过程。
• getIntersection(BVHBuildNode* node, const Ray ray)in BVH.cpp: 建
立 BVH 之后，我们可以用它加速求交过程。该过程递归进行，你将在其中调
用你实现的 Bounds3::IntersectP.
其他，这次作业里面有两个坑。
一个是castRay方法现在在Scene.cpp里面，然后castRay函数发生了一点变化，ray作为参数会传进去，其实ray也就定义了direction和origin两个参数。
细读这次的所有代码，也算是对整个光线追踪有了更深刻的理解吧，可惜我本人还是太菜了。
第二个坑是我一开始没画出来我佛了，后来也是参考了这位大神的代码，发现我对于Bounds包围盒中的求解中，没有考虑方向中有负数的情况，直接佛了。。
首先进入Triangle.h
修改getInsection函数，需要注意的是，这里需要理解以下参数
coords：光线经过的那一点的坐标，这里使用ray（double）（注意，我们在Ray.h中重载了（），用于计算t时刻的光线坐标）
happend：光线和三角形是否相交？
m：材质
normal：交点处三角形的法线
distance：光线经过的时间，也就是我们算出来的t_tmp。
object：事实上你进入Objects.h,就会发现基类Object（接口，里面一堆纯虚函数），这里作为一个指针，指向本对象，也就是this（框架的编程还是很漂亮的），指代我们的三角形。
（注意我们作业5中我自己写的求交，感觉low爆了，还是框架的写的好看）
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);//s1
double det = dotProduct(e1, pvec);
if (fabs(det) < EPSILON)
return inter;//no solution
double det_inv = 1. / det;
Vector3f tvec = ray.origin - v0;//s
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 intersection
if(t_tmp<0)
return inter;
inter.distance = t_tmp;
inter.happened = true;
inter.m = m;
inter.obj = this;
inter.normal = normal;
inter.coords = ray(t_tmp);
return inter;
}
然后进入Renderer.cpp，修改Render函数，注意这里，我们上面两行代码计算x和y就是作业5的答案，然后下面对比作业5的代码出现了变化，因为框架修改了
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 float 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* Vector3f dir = Vector3f(x, y, -1); dir = normalize(dir); Ray ray(eye_pos,dir); framebuffer[m++] = scene.castRay(ray,0); // Don't forget to normalize this direction! } 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 * 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); }
进入Bounds3.cpp ，修改IntersectP函数，此函数用于求包围盒是否和光线相交，这里特别要注意负数的情况，一开始我想后面这个数组有什么用，然后一直没用到，最后卡住了，佛了，其实这个就是判断大小用的，算最大最小的时候注意负数的情况
inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
const std::array<int, 3>& dirIsNeg) const
{
float t_Min_x = (pMin.x - ray.origin.x)*invDir[0];
float t_Min_y = (pMin.y - ray.origin.y)*invDir[1];
float t_Min_z = (pMin.z - ray.origin.z)*invDir[2];
float t_Max_x = (pMax.x - ray.origin.x)*invDir[0];
float t_Max_y = (pMax.y - ray.origin.y)*invDir[1];
float t_Max_z = (pMax.z - ray.origin.z)*invDir[2];
if(!dirIsNeg[0])
{
float t = t_Min_x;
t_Min_x = t_Max_x;
t_Max_x = t;
}
if(!dirIsNeg[1])
{
float t = t_Min_y;
t_Min_y = t_Max_y;
t_Max_y = t;
}
if(!dirIsNeg[2])
{
float t = t_Min_z;
t_Min_z = t_Max_z;
t_Max_z = t;
}
float t_enter = std::max(t_Min_x,std::max(t_Min_y,t_Min_z));
float t_exit = std::min(t_Max_x,std::min(t_Max_y,t_Max_z));
if(t_enter<t_exit&&t_exit>=0)
return true;
else
return false;
// 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
}
最后进入BVH.cpp,修改我们的getInsection,注意，我们这里是为了求交，根据闫神所讲的，肯定要到叶子节点才算求交，而中间如果和包围盒都没交点，那么和三角形一定无交点，简单二叉树遍历就行。
Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const { Intersection intersect; Vector3f invdir(1./ray.direction.x,1./ray.direction.y,1./ray.direction.z); std::array<int, 3> dirIsNeg; dirIsNeg[0] = ray.direction.x>0; dirIsNeg[1] = ray.direction.y>0; dirIsNeg[2] = ray.direction.z>0; if(!node->bounds.IntersectP(ray,invdir,dirIsNeg)) { return intersect; } if(node->left == nullptr && node->right==nullptr) { return node->object->getIntersection(ray); } Intersection h1 = getIntersection(node->left,ray); Intersection h2 = getIntersection(node->right,ray); return h1.distance<h2.distance?h1:h2; return intersect; // TODO Traverse the BVH to find intersection }
其实到这里，我们的任务已经完成了，效果如下：
但是其实还没完，注意我们的作业提高要求：
SAH 查找：自学 SAH(Surface Area Heuristic) , 正确实现 SAH 加速，并且提交结果图片，并在 README.md 中说明 SVH 的实现方法，并对比 BVH、SVH 的时间开销。(可参考 http://15462.courses.cs.cmu.edu/fall2015/lecture/acceleration/slide_024，也可以查找其他资料)
点进去ppt，发现是这个东西
告诉我们是一个启发性算法，这里的思想是如果有两个包围盒，B包围A，那么光线过B也过A的概率就是SA/SB，
注意，我们在这一次作业的时候，观察前面的rebuild方法，是递归地将objects数组分成左右两个部分，那么SA可以说代表左边的objects形成的总包围盒的面积，而SB可以说是代表右边的划分总的包围盒的面积，SN代表我们的最大的包围盒的面积，我们的想法是启发式搜索（哇本人看到启发式搜索就头疼），这里Ctrav表示当前划分一次的花费，这里取1/8（这个值应该是可以版的，具体看https://blog.csdn.net/qq_39300235/article/details/106999514，这篇文章说的很好）
我们的目标是分割这两个数组（搜索），所花费的cost越少越好。（注意，我们一开始划分数组的时候，方法是就是取最大的跨度维度，然后从中间划分）
但是启发式搜索目的是不要让我们从中间划分，而是划分的时候让cost（也就是上面的式子C越小越好），找出这种划分方法就行。
https://blog.csdn.net/qq_39300235/article/details/106999514这篇文章其实很清楚，但是里面一堆桶的说法很麻烦，记住，我们的目的就是找出某个划分方法，然后让cost最小，这样形成的BVH结构不会有很大的深度。
我截取这篇文章的核心思想就是在每一个维度（x，y，z）分别计算一种划分（左，右），目的是求左右之间的中间值
这道题坑的地方是其实他 在bound3s中写了一个offset（），这个东西是就是用来求上面那篇文章的桶的标号的。
还有，这里面有一个bid = bid -1 ，其实是我们分了12个桶（数组最大数11），但是根据offset是可以算出12的（12棵树中间的间隔其实有13个，那么就把12的地方换成11就行，也就是桶的大小最后一个桶会大一点）
感谢https://blog.csdn.net/miyu1994/article/details/107006010/大神的代码，我这里直接用了他的代码，不过标记了一点注释，这样比较容易理解。
1 BVHBuildNode* BVHAccel::recursiveBuild(std::vector<Object*> objects) 2 { 3 BVHBuildNode* node = new BVHBuildNode(); 4 5 // Compute bounds of all primitives in BVH node 6 Bounds3 bounds; 7 for (int i = 0; i < objects.size(); ++i) 8 bounds = Union(bounds, objects[i]->getBounds()); 9 if (objects.size() == 1) { 10 // Create leaf _BVHBuildNode_ 11 node->bounds = objects[0]->getBounds(); 12 node->object = objects[0]; 13 node->left = nullptr; 14 node->right = nullptr; 15 return node; 16 } 17 else if (objects.size() == 2) { 18 node->left = recursiveBuild(std::vector{objects[0]}); 19 node->right = recursiveBuild(std::vector{objects[1]}); 20 21 node->bounds = Union(node->left->bounds, node->right->bounds); 22 return node; 23 } 24 else 25 { 26 Bounds3 centroidBounds; 27 //算出最大的包围盒（通用的，因为两个方法都要用到） 28 for (int i = 0; i < objects.size(); ++i) 29 centroidBounds = 30 Union(centroidBounds, objects[i]->getBounds().Centroid()); 31 32 std::vector<Object*> leftshapes; 33 std::vector<Object*> rightshapes; 34 35 switch (splitMethod)//这里注意了在BVH.h里面有个枚举类，构造函数中的初始将决定是普通方法，还是SAH 36 { 37 case SplitMethod::NAIVE: 38 { 39 int dim = centroidBounds.maxExtent();//算出最大的跨度对应的值，x为0，y为1，z为2 40 int index = objects.size() / 2; 41 switch (dim) 42 //排序，针对最大的跨度排序 43 { 44 case 0: 45 std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) { 46 return f1->getBounds().Centroid().x < 47 f2->getBounds().Centroid().x; 48 }); 49 break; 50 case 1: 51 std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) { 52 return f1->getBounds().Centroid().y < 53 f2->getBounds().Centroid().y; 54 }); 55 break; 56 case 2: 57 std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) { 58 return f1->getBounds().Centroid().z < 59 f2->getBounds().Centroid().z; 60 }); 61 break; 62 } 63 64 auto beginning = objects.begin(); 65 auto middling = objects.begin() + index; 66 auto ending = objects.end(); 67 //递归算法，枢轴是中间元素。 68 leftshapes = std::vector<Object *>(beginning, middling); 69 rightshapes = std::vector<Object *>(middling, ending); 70 } 71 break; 72 case SplitMethod::SAH: 73 { 74 float nArea = centroidBounds.SurfaceArea();//算出最大的 75 76 int minCostCoor = 0; 77 int mincostIndex = 0; 78 float minCost = std::numeric_limits<float>::infinity(); 79 std::map<int, std::map<int, int>> indexMap; 80 //indexmap用于记录x，y，z（前一个int代表x，y，z，后一个map代表那个轴对应的map） 81 //遍历x，y，z的划分 82 for(int i = 0; i < 3; i++) 83 { 84 int bucketCount = 12;//桶的个数，这里定了12个桶，就是在某一个轴上面划分了12个区域 85 std::vector<Bounds3> boundsBuckets; 86 std::vector<int> countBucket; 87 for(int j = 0; j < bucketCount; j++) 88 { 89 boundsBuckets.push_back(Bounds3()); 90 countBucket.push_back(0); 91 } 92 93 std::map<int, int> objMap; 94 95 for(int j = 0; j < objects.size(); j++) 96 { 97 int bid = bucketCount * (centroidBounds.Offset(objects[j]->getBounds().Centroid()))[i];//算出对应x，y。z上的id值，这里【i】代表x，y，z 98 if(bid > bucketCount - 1)//实质是可以划分13个区域的，将最后一个区域合并。 99 { 100 bid = bucketCount - 1; 101 } 102 Bounds3 b = boundsBuckets[bid]; 103 b = Union(b, objects[j]->getBounds().Centroid()); 104 boundsBuckets[bid] = b; 105 countBucket[bid] = countBucket[bid] + 1; 106 objMap.insert(std::make_pair(j, bid)); 107 } 108 109 indexMap.insert(std::make_pair(i, objMap)); 110 //对于每一个划分，计算他所对应的花费，方法是对于桶中的每一个面积，计算他的花费，最后进行计算 111 //找出这个划分。 112 for(int j = 1; j < boundsBuckets.size(); j++) 113 { 114 Bounds3 A; 115 Bounds3 B; 116 int countA = 0; 117 int countB = 0; 118 for(int k = 0; k < j; k++) 119 { 120 A = Union(A, boundsBuckets[k]); 121 countA += countBucket[k]; 122 } 123 124 for(int k = j; k < boundsBuckets.size(); k++) 125 { 126 B = Union(B, boundsBuckets[k]); 127 countB += countBucket[k]; 128 } 129 130 float cost = 1 + (countA * A.SurfaceArea() + countB * B.SurfaceArea()) / nArea;//计算花费 131 //找出这个花费。 132 if(cost < minCost) 133 { 134 minCost = cost; 135 mincostIndex = j; 136 minCostCoor = i; 137 } 138 } 139 } 140 //加入左右数组，这里很重要，具体还是看那篇博客 141 for(int i = 0; i < objects.size(); i++) 142 { 143 if(indexMap[minCostCoor][i] < mincostIndex) 144 { 145 leftshapes.push_back(objects[i]); 146 } 147 else 148 { 149 rightshapes.push_back(objects[i]); 150 } 151 } 152 } 153 break; 154 dafault: 155 break; 156 } 157 158 assert(objects.size() == (leftshapes.size() + rightshapes.size())); 159 //递归计算，同普通方法 160 node->left = recursiveBuild(leftshapes); 161 node->right = recursiveBuild(rightshapes); 162 163 node->bounds = Union(node->left->bounds, node->right->bounds); 164 } 165 166 return node; 167 }
最后终于完成了。
这一次还是很有挑战性的，里面涉及启发式搜索，BVH，bounds求交
