games101 作业6 详解SAH
代码分析
作业6 的代码框架在作业五的基础上进行了进一步地丰富与改进
一开始还是创建场景 添加物体以及灯光 这次的物体只有一个用objloader读取的兔子的3d模型 但是在将该模型转换为meshtriangle数据结构的时候，也为该模型中的各个三角形片元创建BVH（Bounding volume hierarchies）的加速结构 这里虽然也为场景中的各个物体也构建了BVH 但是目前场景中只有一个物体
MeshTriangle bunny(Utils::PathFromAsset("model/bunnyAssignment6/bunny.obj"));
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();
MeshTriangle(const std::string& filename){
.......
std::vector<Object*> ptrs;
for (auto& tri : triangles)
ptrs.push_back(&tri);
bvh = new BVHAccel(ptrs);
}
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;
root = recursiveBuild(primitives);
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);
}
而BVH的构建过程就是一颗树的递归构建的过程 划分BVH节点的方法 就是在最长的维度上 对所有的三角形进行排序 找出位于中间的这个三角形的位置 左边的就是左节点 右边就是右节点，划分结束的终止条件就是只有一个物体
对应课件的ppt：
pbrt图示：
BVHBuildNode* BVHAccel::recursiveBuild(std::vector<Object*> objects)
{
BVHBuildNode* node = new BVHBuildNode();
// Compute bounds of all primitives in BVH node
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;
}
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 {
//用每个三角形对应的box的中点 拼成一个box 然后找出这个box最长的维度 在这个维度上进行划分
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;
}
之后就是光线追踪castRay 求交 然后进行着色计算 和作业5类似
但是这次求交要利用我们构建的BVH加速结构进行加速 这也是我们要完成的部分
分析与解决
与BVH求交
这里主要是对构建的BVH加速结构进行一个遍历 和节点的boundingbox求交 如果没有交点就返回 如果有 则查看该节点是不是叶子节点 不是叶子节点 则继续遍历该节点的子节点 如果是叶子节点 则对其中的物体进行相交测试 找出最近的交点返回：
Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
// TODO Traverse the BVH to find intersection
Intersection hit1;
Intersection hit2;
std::array<int, 3> dirIsNeg{ 0,0,0 };
for (int i = 0; i < 3; i++)
{
if (ray.direction[i] < 0) {
dirIsNeg[i] = 1;
}
}
if (!node->bounds.IntersectP(ray, ray.direction_inv, std::move(dirIsNeg)))
{
return Intersection{};
}
if (node->left == nullptr && node->right == nullptr) {
return node->object->getIntersection(ray);
}
hit1 = getIntersection(node->left, ray);
hit2 = getIntersection(node->right, ray);
return hit1.distance > hit2.distance ? hit2 : hit1;
}
boundingbox求交
boundingbox是个长方体 并且与各个面与xyz轴平行 我们只需要求出和每组平面 的tenter 与 texit 并从enter中找出最大值 exit中找出最小值 这样可以防止错误的将交点找到长方体的延展面上 同时也要考虑dir光线为负的情况 这样的话会和坐标值较小的平面先相交：
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
float txmin = (pMin - ray.origin).x * invDir.x;
float txmax = (pMax - ray.origin).x * invDir.x;
float tymin = (pMin - ray.origin).y * invDir.y;
float tymax = (pMax - ray.origin).y * invDir.y;
float tzmin = (pMin - ray.origin).z * invDir.z;
float tzmax = (pMax - ray.origin).z * invDir.z;
//如果方向某个维度是负的 会和pMax对应维度上的面上先相交
if (dirIsNeg[0]) {
std::swap(txmin, txmax);
}
if (dirIsNeg[1]) {
std::swap(tymin, tymax);
}
if (dirIsNeg[2]) {
std::swap(tzmin, tzmax);
}
float tenter = std::max(std::max(txmin, tymin), tzmin);
float texit = std::min(std::min(txmax, tymax), tzmax);
return tenter < texit && texit >= 0;
}
三角形求交
和作业5中的实现类似
这里需要返回交点的属性 比如材质 法线等等：
inline Intersection Triangle::getIntersection(Ray ray)
{
Intersection inter;
//如果法线与光线异侧 说明光线是从背侧射入
if (dotProduct(ray.direction, normal) > 0)
return inter;
double u, v, t_tmp = 0;
// D E2 = S1
Vector3f pvec = crossProduct(ray.direction, e2);
// S1 . E1
double det = dotProduct(e1, pvec);
if (fabs(det) < EPSILON)
return inter;
double det_inv = 1. / det;
//S
Vector3f tvec = ray.origin - v0;
u = dotProduct(tvec, pvec) * det_inv;
if (u < 0 || u > 1)
return inter;
//S2
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.happened = true;
inter.coords = ray(t_tmp);
inter.normal = normal;
inter.distance = t_tmp;
inter.obj = this;
inter.m = this->m;
return inter;
}
SAH
SAH 全称是Surface Area Heuristic 是在构建BVH过程中 进行划分的一种方法 之前我们也提到了一种简单的划分方法 之前的方法在下图所示的情况下:
会选择图中蓝线这样的划分 最终boundingbox会有部分重叠 从而降低光线求交的效率
SAH提出了一个函数来量化地估计划分区域的代价 以将一个boundingbox划分为AB两个为例：
其中t-trav是遍历树节点的代价 这里固定为一个常数0.125 pbrt的解释是我们将光线求交的代价设置为1 遍历效率大概是求交的1/8 重要的是它们之间的相对关系 而不是数值
t-isect代表着求交的代价 具体设置为每个区域内图元的数目
pa pb代表着光线穿过子节点的概率 使用A的boundingbox面积 与 总boundingbox面积的比值来确定
具体流程概括如下：
1.首先我们和之前的划分算法一样 找出最长的维度
2.这里我们定义一定数目的bucket 在整体的boundingbox上进行划分 计算出每个物体在哪个bucket的索引
3.每个bucket都代表着一种划分方法 我们遍历计算各方法的代价 然后找出最小代价对应的bucket索引：
4.借助最小索引将物体分为两堆
代码：参考https://github.com/mmp/pbrt-v3/blob/master/src/accelerators/bvh.cpp
BVHBuildNode* BVHAccel::SAHrecursiveBuild(std::vector<Object*> objects)
{
BVHBuildNode* node = new BVHBuildNode();
// Compute bounds of all primitives in BVH node
Bounds3 bounds;
for (int i = 0; i < objects.size(); ++i)
bounds = Union(bounds, objects[i]->getBounds());
//假如物体数目过少 就直接使用原来的划分方法 没必要
if (objects.size() <= 2) {
node = recursiveBuild(objects);
return node;
}
else {
Bounds3 centroidBounds;
for (int i = 0; i < objects.size(); ++i)
centroidBounds =
Union(centroidBounds, objects[i]->getBounds().Centroid());
int dim = centroidBounds.maxExtent();
//计算索引
constexpr int nBuckets = 8;
BucketInfo buckets[nBuckets];
for (int i = 0; i < objects.size(); ++i) {
int b = nBuckets *
centroidBounds.Offset(
objects[i]->getBounds().Centroid())[dim];
if (b == nBuckets) b = nBuckets - 1;
buckets[b].count++;
buckets[b].bounds =
Union(buckets[b].bounds, objects[i]->getBounds());
}
//计算每个划分的代价
float cost[nBuckets - 1];
for (int i = 0; i < nBuckets - 1; ++i) {
Bounds3 b0, b1;
int count0 = 0, count1 = 0;
for (int j = 0; j <= i; ++j) {
b0 = Union(b0, buckets[j].bounds);
count0 += buckets[j].count;
}
for (int j = i + 1; j < nBuckets; ++j) {
b1 = Union(b1, buckets[j].bounds);
count1 += buckets[j].count;
}
cost[i] = 0.125f +
(count0 * b0.SurfaceArea() +
count1 * b1.SurfaceArea()) /
bounds.SurfaceArea();
}
//找出最小代价
float minCost = cost[0];
int minCostSplitBucket = 0;
for (int i = 1; i < nBuckets - 1; ++i) {
if (cost[i] < minCost) {
minCost = cost[i];
minCostSplitBucket = i;
}
}
float leafCost = objects.size();
int mid = 0;
//如果包含的物体过多 或者缩小划分代价 则进行划分
if (objects.size() > maxPrimsInNode || minCost < leafCost) {
for (int i = 0; i < objects.size(); ++i) {
int b = nBuckets *
centroidBounds.Offset(
objects[i]->getBounds().Centroid())[dim];
if (b == nBuckets) b = nBuckets - 1;
// 满足条件的元素放到数组前半部分
if (b <= minCostSplitBucket) {
std::swap(objects[i], objects[mid]);
mid++;
}
}
}
auto beginning = objects.begin();
auto middling = objects.begin() + mid;
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 = SAHrecursiveBuild(leftshapes);
node->right = SAHrecursiveBuild(rightshapes);
node->bounds = Union(node->left->bounds, node->right->bounds);
}
return node;
}
结果展示
BVH:
BVH-SAH:
