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::min() 其实是正的最小值,应该用 lowest,等价于 -max()。tExit初始化为最大值,当然两个用正无穷和负无穷其实也行。

包围盒是用对角线上两个点表示的,分别是 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,如果超过几十秒肯定是有问题的更不用说那些跑了几十分钟的

image-20220709165157978

另外顺带提一下,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%,如果面数和物体数量再多一些,分散的不均匀一些效果估计会明显一点,这里只有一个物体,所有三角形都是集中在一起的,其实效果不会差太多的。