games101 作业笔记（作业 6）
这次作业总体难度不大，主要麻烦在找各种函数和变量是什么，互相之间有什么关系，理清这个就很麻烦，都理清了以后，就是按照课上的公式来就行了
首先是 Rander 函数，和上次区别不大，坐标的转换已经帮我们写好了，这次缺少的反而是上次给的。第一步还是定义光线方向，并归一化，没什么区别
Vector3f dir = normalize(Vector3f(x, y, -1)); // Don't forget to normalize this direction!
第二步调用 castRay 函数计算出当前像素的颜色，这里原来的 castRay 是在 Render.cpp 里的，新的框架将其移到了 Scene 类中，因此需要通过 scene.castRay 来调用，而且这次函数的参数也有所不同，一个是因为函数移到了 Scene 类中，因此就不需要再传递 scene 过去了，另一个则是光线单独用一个 Ray 类来管理了，传递的参数从光源起点和方向两个参数，变成传递一个 Ray 的对象了，因此我们需要通过原来这两个参数构造个 ray 传递过去
framebuffer[m++] = scene.castRay(Ray(eye_pos, dir), 0);
然后就是第二个需要修改的函数 Triangle::getIntersection，这个函数前半部分就是上次的 rayTriangleIntersect，虽然这个函数这次框架里也还在，但两个函数好像用途有些许不同，新的这个是需要返回一个 Intersection 对象，我们需要补充的就是构造这个对象部分的代码
不知道为什么原来的代码似乎没处理 t < 0 的情况，可能前面的多个判断已经排除了这种情况吧，不过我还是加上了这句
Ray类中重载了()运算符，可以直接返回光线方程的值，因此我们可以直接调用
// Ray.hpp
Vector3f operator()(double t) const{return origin+direction*t;}
// 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; // 光源与交点的距离，后面用来找最先相交的点用的，用 t 就行了
inter.obj = this; // 相交的三角形对象，就是当前对象 this
inter.m = m; // 交点的材质就是三角形面的材质
然后就是这次的新内容 IntersectP 函数判断光线和包围盒是否相交，基本按照公式就行了，tEnter 初始化为 double 的最小值，需要注意 std::numeric_limits
包围盒是用对角线上两个点表示的，分别是 pMin 和 pMAx，pMin 的 xyz 三个分量就可以代表左底前三个面了，同理 pMax 也是一样。这里唯一需要注意的就是，如果光线的某个分量是往负方向走的，我们需要交换一下 tMin 和 tMax，因为负方向的话就是先和右上后面相交了
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
double tEnter = std::numeric_limits<double>::lowest();
double tExit = std::numeric_limits<double>::max();
for(int i=0; i<3; i++) {
double tMin = (pMin[i] - ray.origin[i]) * invDir[i];
double tMax = (pMax[i] - ray.origin[i]) * invDir[i];
if(!dirIsNeg[i]) std::swap(tMin, tMax);
tEnter = std::max(tEnter, tMin);
tExit = std::min(tExit, tMax);
}
return tEnter < tExit && tExit >= 0;
}
最后就是 BVH 判断相交了，这里我们可以发现，叶子节点只存了一个对象，也就是最后只需要和一个三角形判断是否相交，这个只需要调用我们刚刚写的 Triangle::getIntersection 即可，那么整体的思路就是按照树的先序遍历的方法遍历，直到遇到叶子节点或者一个包围盒不和光线相交的中间节点就返回
然后这里判断包围盒是否相交的这个 IntersectP 函数还需要传递一个判断某个分量是否是正的的数组，需要构造一下。
然后就是左子树和右子树都遍历一遍后，选择交点距离最近的，也就是最先相交的返回，这里记得一定要用临时变量暂存左右子树的结果，不要直接调用的同时判断返回值，然后再调用一次把结果返回，每个节点重复调用一次，最坏情况下计算量可是 2 的指数级的增加。
Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
// TODO Traverse the BVH to find intersection
if(node->bounds.IntersectP(ray, ray.direction_inv, {ray.direction.x > 0, ray.direction.y > 0, ray.direction.z > 0})) {
if(node->left == nullptr && node->right == nullptr)
return node->object->getIntersection(ray);
else {
Intersection isectLeft = getIntersection(node->left, ray);
Intersection isectRight = getIntersection(node->right, ray);
return isectLeft.distance < isectRight.distance ? isectLeft : isectRight;
}
}
return Intersection();
}
我这边渲染出来只用了 4s，如果超过几十秒肯定是有问题的更不用说那些跑了几十分钟的
另外顺带提一下，ppm文件 linux 自带的图片管理器应该就能打开，windows 的话用 PS 之类的也能打开，或者网上随便找个在线 ppm 转 png 的网站转换一下也行
提高
这是我感觉目前提高里最难的一次，作业文档里给的参考资料虽然思路和公式还好理解，但要写成代码其实不大好写，网上参考了很多代码，但都各不相同。
整体的思路其实还好理解，普通的BVH是直接从中间划分的，虽然很方便，但是当物体分布不太均匀，有可能就会出现包围盒重叠部分过多的情况。因此 SAH 就是考虑如何划分才能重叠。
如果理解不了老师给的资料的话，可以参考下这篇文章，我觉得姑且也还算清楚 PBRT-E4.3-层次包围体(BVH)（一）
然后下面是我的代码，首先我定义了桶相关的数据结构和参数，桶数量上资料里说是小于32，网上很多代码似乎都选择了12，我这也就照搬了，有兴趣的可以改成不同大小的看看效果
constexpr int nBuckets = 12;
constexpr double cTrav = .125;
constexpr double cIsect = 1;
struct BucketInfo {
Bounds3 bounds;
int count = 0;
};
然后修改 BVHAccel 构造函数中调用 recursiveBuild 构造 bvh 的部分
root = splitMethod == SplitMethod::NAIVE ? recursiveBuild(primitives) :
recursiveBuildWithSAH(primitives);
当然实际上SAH的这个函数前半部分和原本的是一样的，物体个数是 1 和 2 的时候是一样的，直接从 else 开始修改也是可以的。
网上的各类博文里，有的人认为物体数量较小（比如小于桶的数量）时，用原本的二分可能效率更高，我这边没有尝试过，我觉得也可以试试
然后 else 里面先计算出整个的包围盒，并初始化一些变量
Bounds3 centroidBounds;
for (int i = 0; i < objects.size(); ++i)
centroidBounds = Union(centroidBounds, objects[i]->getBounds().Centroid());
int bestDim = 0;
int bestSplit = 0;
double minCost = std::numeric_limits<double>::max();
然后可以参考老师给的资料中的伪代码，遍历三个轴
先初始化一个桶，然后计算每个物体应该属于哪个桶，Offset 函数可以计算出点在包围盒内的空间，而且范围是在[0, 1]，这个函数也是原本代码里就提供给我们的，然后乘上桶的数量，就可以映射到 [0, nBuckets] 的范围。这里我发现老师的 Vector3f 重载[]的函数是有问题的，所以这里要么修改原本的函数，要么就手动判断，反正其实效果是一样的。
后面就是遍历 nBuckets - 1 种划分方式，找到其中的代价最小的即可
for(int dim = 0; dim < 3 ; ++dim) {
std::array<BucketInfo, nBuckets> buckets;
for(auto &object:objects) {
auto offset = centroidBounds.Offset(object->getBounds().Centroid());
int b = nBuckets * (dim == 0 ? offset.x :
dim == 1 ? offset.y : offset.z);
if(b == nBuckets) b = nBuckets - 1;// 在边界上的要单独判断
buckets[b].count++; // 桶里物体数量+1
buckets[b].bounds = Union(buckets[b].bounds, object->getBounds());
}
for(int i = 0; i < nBuckets - 1; ++i) {
Bounds3 boundsLeft, BoundsRight;
int countLeft = 0, countRight = 0;
for(int j = 0; j <= i; ++j) {
countLeft += buckets[j].count;
boundsLeft = Union(boundsLeft, buckets[j].bounds);
}
for(int j = i + 1; j < nBuckets; ++j) {
countRight += buckets[j].count;
BoundsRight = Union(BoundsRight, buckets[j].bounds);
}
double cost = cTrav + (boundsLeft.SurfaceArea() * countLeft +
BoundsRight.SurfaceArea() * countRight) * cIsect / centroidBounds.SurfaceArea();
if(cost < minCost) {
minCost = cost;
bestDim = dim;
bestSplit = i;
}
}
}
最后就是按照最小代价的方式划分即可，这里用了 STL 自带的划分函数，类似于老师上课讲的快速选择算法
auto middling = std::partition(objects.begin(), objects.end(), [=](auto object) {
auto offset = centroidBounds.Offset(object->getBounds().Centroid());
int b = nBuckets * (bestDim == 0 ? offset.x : bestDim == 1 ? offset.y : offset.z);
if(b == nBuckets) b = nBuckets - 1;
return b <= bestSplit;
});
auto beginning = objects.begin();
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);
完整代码如下
BVHBuildNode* BVHAccel::recursiveBuildWithSAH(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 {
Bounds3 centroidBounds;
for (int i = 0; i < objects.size(); ++i)
centroidBounds =
Union(centroidBounds, objects[i]->getBounds().Centroid());
int bestDim = 0;
int bestSplit = 0;
double minCost = std::numeric_limits<double>::max();
for(int dim = 0; dim < 3 ; ++dim) {
std::array<BucketInfo, nBuckets> buckets;
for(auto &object:objects) {
auto offset = centroidBounds.Offset(object->getBounds().Centroid());
int b = nBuckets * (dim == 0 ? offset.x :
dim == 1 ? offset.y : offset.z);
if(b == nBuckets) b = nBuckets - 1;
buckets[b].count++;
buckets[b].bounds = Union(buckets[b].bounds, object->getBounds());
}
for(int i = 0; i < nBuckets - 1; ++i) {
Bounds3 boundsLeft, BoundsRight;
int countLeft = 0, countRight = 0;
for(int j = 0; j <= i; ++j) {
countLeft += buckets[j].count;
boundsLeft = Union(boundsLeft, buckets[j].bounds);
}
for(int j = i + 1; j < nBuckets; ++j) {
countRight += buckets[j].count;
BoundsRight = Union(BoundsRight, buckets[j].bounds);
}
double cost = cTrav + (boundsLeft.SurfaceArea() * countLeft +
BoundsRight.SurfaceArea() * countRight) * cIsect / centroidBounds.SurfaceArea();
if(cost < minCost) {
minCost = cost;
bestDim = dim;
bestSplit = i;
}
}
}
auto middling = std::partition(objects.begin(), objects.end(), [=](auto object) {
auto offset = centroidBounds.Offset(object->getBounds().Centroid());
int b = nBuckets * (bestDim == 0 ? offset.x :
bestDim == 1 ? offset.y : offset.z);
if(b == nBuckets) b = nBuckets - 1;
return b <= bestSplit;
});
auto beginning = objects.begin();
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;
}
可能因为我电脑处理器还算不错，原本就是渲染时间也就 4s，所以换成 SAH 后看不出提升还是 4s，反而建树的时间偶尔会达到 1s，不过大多数情况似乎还是 0s。加了毫秒显示以后似乎平均都是 4800 毫秒左右，提升似乎不到 1%，如果面数和物体数量再多一些，分散的不均匀一些效果估计会明显一点，这里只有一个物体，所有三角形都是集中在一起的，其实效果不会差太多的。
OωO
