【GAMES101】作业6 加速结构

一.作业描述

二.作业解析

1. 递归构造BVH树分析

看视频学习后过了太久，BVH树怎么构造已经忘记了，这里重新根据作业提供的代码整理了一下思路，重新学习了。后面要递归BVH树去判断相交，所以了解一下BVH树的结构也是很必要的。

BVH树节点node属性：

Bound：包围盒；
left，right：左右指针；
object：Object类；

递归实现：
1.递归出口：
当场景中物体总数=1时，包围盒为当前物体的包围盒，左右节点为NULL，返回当前节点；

2.当场景中物体总数=2时，取两个物体的边界并集作为当前节点的包围盒，左右节点分别再进行一次递归；

3.一般递归思路：
取当前所有物体的包围盒并集作为当前节点的总包围盒大小，求出总包围盒大小三维（x，y，z）中长度最大的一个维度dim，并将所有物体的质心坐标在维度dim上的值进行排序，将排序后的前1/2物体和后1/2的物体分别作为左右子集继续进行递归。

归纳一下，BVH树叶节点中存储比较重要的信息就是包围盒bounds。
在遍历BVH树判断相交时，如果发现当前节点的左右指针均为空，就说明当前节点的包围盒大小已经不可再分了。

递归构造的代码在这里mark一下。

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 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;
}

2. Render函数

Render仅有些数据结构上的简单变动。
作业代码实现：

//Render函数：
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*
            // Don't forget to normalize this direction!
            Vector3f dir = normalize(Vector3f(x, y, -1)); 
            Ray ray=Ray(eye_pos,dir,0);
            framebuffer[m++] = scene.castRay(ray,0);
            

        }
        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);    
}

3. Triangle::getItersection函数

求三角形与射线是否相交的getItersection函数的基础判断部分已经实现了，相比作业5返回值出现了一点变化，之前是返回true，这次时返回结构体Intersection类型的inter，return前给各个属性赋上计算出的值即可。

Intersection属性:

happened：是否有交点；
coords：交点坐标；
normal：交点处平面法向量；
obj：相交物体；
m:相交点面材质；

如果判断相交，则修改inter中的属性即可。
如下公式所示，根据射线ray的定义，t代表的是射线初始点到交点的距离。
计算出t的值以后，可根据origin+t*dir求出相交点坐标。

依次为inter中的属性赋值，具体代码如下。

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);
    double det = dotProduct(e1, pvec);
    if (fabs(det) < EPSILON)
        return inter;

    double det_inv = 1. / det;
    Vector3f tvec = ray.origin - v0;
    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
    inter.happened=true;
    inter.coords=ray(t_tmp);
    inter.normal=normal;
    inter.obj=this;
    inter.distance=t_tmp;
    inter.m=this->m;

    return inter;
}

4. 射线与AABB包围盒相交判断

二维判断方法如课程中所讲，最终结果算出近端tnear取从x，y两个维度txmin和tymin中的最大值，tfar取x，y两个维度txmax和tymax中的最小值。
三维判断与二维判断思路一致，只是多了一个维度。

概括一下判断相交思路：
如果texit<0，那么包围盒在射线”后面“，没有相交；
如果tenter<texit且texit>=0，射线与包围盒相交；

实现代码：

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
    auto Vmin=(pMin-ray.origin)*invDir;
    auto Vmax=(pMax-ray.origin)*invDir;
    if(dirIsNeg[0]==1)std::swap(Vmin.x,Vmax.x);
    if(dirIsNeg[1]==1)std::swap(Vmin.y,Vmax.y);
    if(dirIsNeg[2]==1)std::swap(Vmin.z,Vmax.z);


    float tmin=fmax(Vmin.x,fmax(Vmin.y,Vmin.z));
    float tmax=fmin(Vmax.x,fmin(Vmax.y,Vmax.z));

    if(tmin<tmax&&tmax>=0)return true;
    else return false;
}

5. BVH::getItersection函数

加速后的求交的判断步骤为：
深度优先遍历BVH树，判断射线是否与当前节点的轴对称包围盒AABB相交，如果不相交，则直接结束；如果相交，判断当前节点的包围盒是否可以继续划分（即当前节点是否有左右子节点），若可以继续划分，返回左右求交获得的intersection中最近的intersection。

代码如下：

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    
    // TODO Traverse the BVH to find intersection
    Intersection inter,interL,interR;


    std::array<int, 3> dirIsNeg;
    if(ray.direction.x<0)dirIsNeg[0]=1;
    if(ray.direction.y<0)dirIsNeg[1]=1;
    if(ray.direction.z<0)dirIsNeg[2]=1;

    bool flag=node->bounds.IntersectP(ray,ray.direction_inv,dirIsNeg);
    if(!flag)return inter;
    if(node->left==nullptr&&node->right==nullptr){
    
        return inter=node->object->getIntersection(ray);
    }

    interL=getIntersection(node->left,ray);
    interR=getIntersection(node->right,ray);

    return interL.distance<interR.distance?interL:interR;

}

6. SAH优化思路

BVH存在着包围盒重叠的问题，然而划分时要尽可能减少划分后两部分包围盒重叠的体积，因为重叠的体积越大，射线穿过重叠区域的可能性越大，遍历两个子树的可能性就越高，计算消耗越多。

此外，如果空间中物体分布得很不均匀，比如右边存在大量物体，而左边只有少量物体，此时使用BVH就会得到很差的划分结果。

SAH（Surface Area Heuristic）基于表面积对划分进行评估，这种方法通过对求交代价和遍历代价进行评估，给出了每一种划分的代价（Cost），而我们的目的便是去寻找代价最小的划分。
内容参考：PBRT-E4.3-层次包围体(BVH)
↓↓↓↓↓↓↓↓↓↓

c(A,B)为当前划分的代价；
p(A),p(B)为射线落在划分的两个子包围盒的概率；
t(i),t(j)为对每个物体求交的代价，t_trav为遍历代价。

SAH考虑了物体在空间中的分布，也考虑到了划分后两部分包围盒重叠的体积问题，优化了性能。

这里就不放具体的代码了，我尝试了：
B站 SAH和知乎 SAH
这两种方法大差不差，处理细节上有些不同，但是我的虚拟机处处透露着离谱的气息，单BVH的render都用了15s，别人描述的都是6s，SAH优化以后时间还多了1s（实属反向优化了），后面越来越慢…就当是先对SAH做了个简单的了解吧。

实现效果：