作业要求

在之前的编程练习中，我们实现了基础的光线追踪算法，具体而言是光线传输、光线与三角形求交。我们采用了这样的方法寻找光线与场景的交点：遍历场景中的所有物体，判断光线是否与它相交。在场景中的物体数量不大时，该做法可以取得良好的结果，但当物体数量增多、模型变得更加复杂，该做法将会变得非常低效。因此，我们需要加速结构来加速求交过程。在本次练习中，我们重点关注物体划分算法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& dirIsNeg) in the Bounds3.hpp: 这个函数的作用是判断包围盒BoundingBox与光线是否相交，你需要按照课程介绍的算法实现求交过程。 • getIntersection(BVHBuildNode* node, const Ray ray)in BVH.cpp: 建立BVH之后，我们可以用它加速求交过程。该过程递归进行，你将在其中调用你实现的Bounds3::IntersectP High-Level Request：自学 SAH(Surface Area Heuristic) , 正确实现 SAH 加速，并且提交结果图片，并在 README.md 中说明 SVH 的实现方法，并对比 BVH、SVH 的时间开销。(可参考 http://15462.courses.cs .cmu.edu/fall2015/lecture/acceleration/slide_024，也可以查找其他资料。这次的作业就是上次作业的一个升级，从性能开销方便做加速优化。

知识点回顾： BVH加速算法

层次包围体技术 (BVH) 指的是将所有包围体分层逐次地再次包围，获得一个更大的包围体，直到包围住所有物体。实际上，它是一个树形结构，因此可以仿照树的结构，将两个或三个小的包围体包围成一个更大的包围体，以此类推。如图，从大节点到小节点，我们逐步套索并分类全部的物体，我们先判定与规范大包围盒是否相交，再逐个随着树深度深入寻找相交情况。另外，对于静止物体和运动物体，我们推荐使用不同的两套处理逻辑进行相应的实现。（对于大量动态物体，更推荐使用静态网格体对各部分进行筛选与判断。）

BVH树的构建

在传统的BVH构建方法当中，我们首先创建根节点，再递归的创建好相应的子节点。

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());
    // 当size为1时，无需再细分，即为叶子节点
    if (objects.size() == 1) {
        // 具体代码，篇幅考虑省略
    }
    // 当size为2时，直接二分即可，无需采用其他策略
    else if (objects.size() == 2) {
        // 具体代码，篇幅考虑省略
    }
    else {
        Bounds3 centroidBounds;
        for (int i = 0; i < objects.size(); ++i)
            centroidBounds = Union(centroidBounds, objects[i]->getBounds().Centroid());

        int dim = centroidBounds.maxExtent();
        // 根据dimension的最长轴，选择对应轴进行包围盒划分
        switch (dim) {
        case 0:  // x轴最长
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().x <
                       f2->getBounds().Centroid().x;});
            break;
        case 1:  // y轴最长
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().y <
                       f2->getBounds().Centroid().y;});
            break;
        case 2:  // z轴最长
            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;
}

从逻辑上看，这个逻辑就是：

初始化节点：创建一个新的BVH节点。
计算包围盒：计算输入物体集合的包围盒。
叶子节点处理：如果物体集合只有一个物体，将该物体作为叶子节点，停止递归。
两物体特殊情况处理：如果物体集合有两个物体，直接将它们分为左右子节点。
中间节点处理：
- 计算物体的质心包围盒。
- 找出质心包围盒最长的轴作为分割轴。
- 根据分割轴对物体进行排序。
- 将物体分为两组，递归构建左右子树。
- 更新当前节点的包围盒为左右子树包围盒的并集。
返回节点：返回构建好的BVH节点。

知识点回顾：相交判定

一般而言，光线的定义方式为点向式，那么确定线上一点只需要知道距离起点的距离t就可以唯一确定光线路径上的某一点。以与aabb包围盒观察的时候的x，y转写著来的四个t值（两组进入+离开）。如图所示，与其相交的最核心因素就是存在t同时满足x y两个方向均位于与aabb相交的范围内。观察上述三幅图可以得出，只要发生区间交叠，光线与平面就能相交，那么区间交叠出现的条件便是：光线进入平面处的最大t值小于光线离开平面处的最小t值那么问题就变成了如何求光线进入平面处的最大t值以及光线离开平面处的最小t值其中：光线的参数方程为R(t) = O + t * Dir 一般平面方程为aX+bY+cZ+d=0，因为AABB的六个面分别平行于XY、XZ、YZ平面，所以平面的方程为X=d，Y=d，Z=d 光线与垂直于x轴的两个面相交时，t = (d - O.x) / Dir.x 光线与垂直于y轴的两个面相交时，t = (d - O.y) / Dir.y 光线与垂直于z轴的两个面相交时，t = (d - O.z) / Dir.z （除以距离角度量保证进行了相应方向的投影。）如果为三维实现，我们判定交界的问题就被转化为筛选最大下界以及最小上界的大小关系；如果其中的下界小于上界并且t值均为正值（这是一个易于忽略的问题，因为光线需保证其射线性且我们统一好了光线t参数的值，为了保证鲁棒性在此声明，避免发生意料外的错误）,我们就认为光线与包围盒成功相交。

实现

移植Render（）

我们发现render方法中仍然是逐像素的计算部分没有给出。我们可以将其获得视角视线的部分改为ndcX,ndcY的签名，并由此向下继续进行；在框架5中，我们的光线投射效果如下：

Vector3f castRay(
        const Vector3f &orig, const Vector3f &dir, const Scene& scene,
        int depth)

在这里的框架中，我们需要调用封装在场景类中的光线投射方法：

Vector3f Scene::castRay(const Ray &ray, int depth)

阅读方法中，我们深度设置为零，就可以将其渲染的视线平面上；关于前一个封装接口Ray，我们在计算时候已经设定好了视线方向，这显然是ray中参与的其中一个参数；（这个光线投射效果是递归的，因此从0开始记参数很合理)

Ray(const Vector3f& ori, const Vector3f& dir, const double _t = 0.0)

我们根据上述构造函数初始化对应类，根据框架中计算的方向和我们自己设定的ori，初始化构造一个ray传入参数，就完成了render（）方法的移植。

移植Triangle::getIntersection（）

首先Intersection创建时默认为False，那么我们在判断相交的时候，只需要逐个判定不相交的情况，如果通过了所有检测，那么我们将其置为相交逻辑这个方法其实已经规范的基本实现完成了。一共完成了以下思路：

判断光线是否从三角形正面射入，如果不是则直接返回；
验证光线是否与三角形平面平行，若平行也返回
判断光线交点是否位于三角形内部，若不在内部同样返回；
确保交点在射线的正方向上，若不在正方向也返回。
如果所有条件都满足，最后会将交点的信息赋值一个对象并返回，表示发生了相交，返回的点即发生相交的位置。我们需要补充的部分就是对应的交点属性部分，参考类进行视线。

inter.happened = true;
inter.coords = ray(t_tmp);
inter.normal = normal;
inter.distance = t_tmp;
inter.obj = this;
inter.m = m;

补全IntersectP（）

在这个方法中，我们考虑的时包围盒与光线的相交关系。首先，我们记录光源距离包围盒的最近、最远端位置。其中光线是利用源头、方向向量以及步长t及进行定义，在这里我们主要是目的是得到对应的t。

 float x_t_min = (pMin.x - ray.origin.x) * invDir.x;
 float x_t_max = (pMax.x - ray.origin.x) * invDir.x;
 float y_t_min = (pMin.y - ray.origin.y) * invDir.y;
 float y_t_max = (pMax.y - ray.origin.y) * invDir.y;
 float z_t_min = (pMin.z - ray.origin.z) * invDir.z;
 float z_t_max = (pMax.z - ray.origin.z) * invDir.z;

其中invDir部分是对应分量的倒数。在加速计算的逻辑下，用倒数乘法加速计算是一种常见的思想。

由于t的定义会受到源头到包围盒的方向影响（可能在world space的作用下，pMin的值小于管线源头，也就是光源可能更靠近坐标原点的情况下某个t值的大小关系产生的问题，所以在这里给出三个判定位置，如果正向，我们把实际上的minmax进行一个调换。

    if (!dirIsNeg[0]) std::swap(x_t_min, x_t_max);
    if (!dirIsNeg[1]) std::swap(y_t_min, y_t_max);
    if (!dirIsNeg[2]) std::swap(z_t_min, z_t_max);

之后，我们要计算光线进入的最晚时间以及光线离开的最早时间并进行比较，完成这一功能的书写。（详情见上文：相交判断的思路）

实现getIntersection（）

我们已经做好了对包围盒的逐层判断，那么我们需要建立对应的树，从根节点向下逐步划分进行相交判定再实现相应的相交逻辑。主要的逻辑就是，如果发生了相交那么向下遍历左右节点，若达到叶子节点，则执行相应部分的渲染逻辑。需要补全的部分如下：

首先对方向进行初步判定，用于传入对应的Intersect中。
创建相交结构，并向下判断，如果并未存在相交如果为空树那么直接返回，由于默认构造为false直接判断为未相交。
在叶子节点处，我们已经确定对应的物体很大可能出现相交的逻辑，那么我们调取里面对应的物体与光线相交（此处与上篇作业部分逻辑接洽）
在中间节点的时候我们采取先序遍历，当中间节点的左右子节点都与光线相交时，它会比较左右子节点返回的交点距离，选择更近的那个交点作为自己的返回值。这样可以确保在递归过程中，每一步都能传递回当前最优的交点信息，提高遍历效率。中间节点的这一逻辑是BVH树高效遍历的关键环节，无法被其他更简单的方式完全替代。

std::array<int, 3> disIsNeg = { ray.direction.x > 0, ray.direction.y > 0, ray.direction.z > 0 };
    Intersection inter;

    if (!node || !node->bounds.IntersectP(ray, ray.direction_inv, disIsNeg)) {
        return inter;
    }

    //若达到叶子节点时
    if (!node->left && !node->right) {
        return node->object->getIntersection(ray);
    }

    //若达到中间节点
    Intersection leftInter = getIntersection(node->left, ray);
    Intersection rightInter = getIntersection(node->right, ray);
    return leftInter.distance < rightInter.distance ? leftInter : rightInter;

效果（FirstLevel）

蓝色背景渲染的兔子（Low Polygen）

提高任务：实现SAH加速结构

实现 SAH 加速，通过修改 recursiveBuild （）进行实现；

实现逻辑

这里使用的是基于分桶（bucketing）的方法。代码段会对每个轴 (x, y, z) 进行操作：

对于每个轴，首先初始化一个桶（bucket）：创建一个大小为B的桶数组。B通常比较小，例如小于32。
然后计算每一个物体p的质心（centroid），看看这个物体落在哪个桶里。
- 将物体p的包围盒与桶b的包围盒做并集操作，也就是扩展桶b的包围盒，使其能够包含物体p。
- 增加桶b中的物体计数。
对于每个可能的划分平面（总共有B-1个），使用表面积启发式（SAH）公式评估其成本。
执行成本最低的划分（如果找不到有效的划分，就将当前节点设为叶子节点）。原本需要对所有物体的每一种可能划分进行评估，现在只需要对B-1个划分进行评估。因此，分桶方法可以在构建BVH时，有效地降低计算复杂度，提高算法的效率。 SAH定义了一个花费模型，用于估算射线-对象相交的计算成本。这个模型假设，当射线穿过一个节点（或空间区域）时，将产生一定的花费（表示为 $C$ ）。这个花费与光线穿过子节点的概率 $p_x$ 成正比。另外，也需要考虑额外的花费（表示为 $C_{trav}$ ），比如计算射线与包围盒的相交等。 划分策略：SAH的目标是找到一种空间划分方式，使得整体花费最小。假设有A个物体被划分到x子节点，B个物体被划分到y子节点，且假设穿过子节点的概率p与该节点的包围盒大小成正比。那么，空间划分的总花费C可以近似为：

其中， $S_A$ 、 $S_B$ 分别表示x、y子节点的表面积， $S_N$ 表示整个节点的表面积， $N_A$ 、 $N_B$ 分别表示x、y子节点中的物体数量， $C_{\text {isect }}$ 表示射线与物体相交的计算成本。 空间划分：为了实现上述的优化目标，利用桶分类即可加速分区。

具体实现

通过将物体分组到多个桶中，探索分割点的候选位置，从而找到最优的分割点以最小化光线与BVH树遍历中的计算成本。
具体来说，桶分类的逻辑包括：计算总表面积用于归一化左右子节点概率，根据物体质心在最大扩展轴上的位置对物体进行排序，将整个物体列表均匀分成10个桶并依次尝试每两个桶之间的分割点计算SAH成本，最后记录所有候选分割点中SAH成本最小的那个点作为最优分割点。
这种做法在减少计算成本的同时，尽可能地找到能降低SAH成本的分割点，但可能存在一些优化空间，比如动态调整桶的数量、预计算累积表面积或使用更稀疏的分割点采样等。

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 {
        Bounds3 centroidBounds;
        for (int i = 0; i < objects.size(); ++i)
            centroidBounds =
            Union(centroidBounds, objects[i]->getBounds().Centroid());
        int maxExtentDimension = centroidBounds.maxExtent();
        switch (SplitMethod::SAH) {
        case SplitMethod::NAIVE: {

            switch (maxExtentDimension) {
            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);
            break;
        }
        case SplitMethod::SAH: {
            float totalSurfaceArea = centroidBounds.SurfaceArea();
            int bucketSize = 10;
            int optimalSplitIndex = 0;
            float minCost = std::numeric_limits<float>::infinity(); //最小花费

            // Sort objects based on their centroid along the max extent dimension
            switch (maxExtentDimension) {
            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;
            }

            // Find the split that minimizes the SAH cost
            //在计算开销的逻辑中，我们要寻找一种合理的桶建立方式。
            for (int i = 1; i < bucketSize; i++) {
                auto beginning = objects.begin();
                auto middling = objects.begin() + (objects.size() * i / bucketSize);
                auto ending = objects.end();
                auto leftshapes = std::vector<Object*>(beginning, middling);
                auto rightshapes = std::vector<Object*>(middling, ending);
                //求左右包围盒:
                Bounds3 leftBounds, rightBounds;
                for (int k = 0; k < leftshapes.size(); ++k)
                    leftBounds = Union(leftBounds, leftshapes[k]->getBounds().Centroid());
                for (int k = 0; k < rightshapes.size(); ++k)
                    rightBounds = Union(rightBounds, rightshapes[k]->getBounds().Centroid());
                float SA = leftBounds.SurfaceArea(); //SA
                float SB = rightBounds.SurfaceArea(); //SB
                float cost = 0.125 + (leftshapes.size() * SA + rightshapes.size() * SB) / totalSurfaceArea; //计算花费
                if (cost < minCost) {//如果花费更小，记录当前坐标值
                    minCost = cost;
                    optimalSplitIndex = i;
                }
            }
            //找到optimalSplitIndex后的操作等同于BVH
            auto beginning = objects.begin();
            auto middling = objects.begin() + (objects.size() * optimalSplitIndex / bucketSize);//划分点选为当前最优桶的位置
            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);//返回pMin pMax构成大包围盒

        }
        }
    }
    return node;
}

为了进行两种结构的加速比较，我们进行一个switch开关通过控制enum来切换加速结构的选择，经测速度由5.0s——4.0s左右，当然在我的好机器上效果又有出入。这种toy级别的加速处理需要在后续深入了解后解决。

You'll Never
Walk Alone

GAMES101(6)-加速结构

作业要求

知识点回顾： BVH加速算法

BVH树的构建

知识点回顾：相交判定

实现

移植Render（）

移植Triangle::getIntersection（）

补全IntersectP（）

实现getIntersection（）

效果（FirstLevel）

提高任务：实现SAH加速结构

相关知识

实现逻辑

具体实现

You'll NeverWalk Alone

作业要求

知识点回顾： BVH加速算法

BVH树的构建

知识点回顾： 相交判定

实现

移植Render（）

移植Triangle::getIntersection（）

补全IntersectP（）

实现getIntersection（）

效果（FirstLevel）

提高任务：实现SAH加速结构

相关知识

实现逻辑

具体实现

You'll Never
Walk Alone

知识点回顾：相交判定