games101 作业6
这是games101 现代计算机图形学 作业06 by coolwx
感谢https://blog.csdn.net/miyu1994/article/details/107006010/ 这位大佬的文章,因为我一开始真的写不出来SAH(启发式搜索),看了这位大神的算法,我最终才理解了他的算法,其实还是比较简单的,关于SAH,我应该会写一篇笔记,来记录我的感想,再次感谢以上这位大神的算法,因为我一开始真的想不出这个启发式算法,所以代码我会使用他的代码,然后添上我的注释,便于理解(哎,还是太菜,没写出来)。
作业要求:
其他,这次作业里面有两个坑。
一个是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 }
其实到这里,我们的任务已经完成了,效果如下:
但是其实还没完,注意我们的作业提高要求:
![](https://img2020.cnblogs.com/blog/1729737/202102/1729737-20210205020417833-1388539937.jpg)
告诉我们是一个启发性算法,这里的思想是如果有两个包围盒,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求交