games101——作业6


总览

在之前的编程练习中,我们实现了基础的光线追踪算法,具体而言是光线传输、光线与三角形求交。我们采用了这样的方法寻找光线与场景的交点:遍历场景中的所有物体,判断光线是否与它相交。在场景中的物体数量不大时,该做法可以取得良好的结果,但当物体数量增多、模型变得更加复杂,该做法将会变得非常低效。因此,我们需要加速结构来加速求交过程。在本次练习中,我们重点关注物体划分算法 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<int, 3>& dirIsNeg) in the Bounds3.hpp: 这个函数的作用是判断包围盒 BoundingBox 与光线是否相交,你需要按照课程介绍的算法实现求交过程。
  • getIntersection(BVHBuildNode* node, const Ray ray) in BVH.cpp: 建立 BVH 之后,我们可以用它加速求交过程。该过程递归进行,你将在其中调用你实现的Bounds3::IntersectP.

开始实现

编译运行

基础代码只依赖于 CMake,下载基础代码后,执行下列命令,就可以编译这个项目:

mkdir build
cd ./build
cmake ..
make

在此之后,你就可以通过 ./Raytracing 来执行程序。

代码框架

我们修改了代码框架中的如下内容:

  • Material.hpp: 我们从将材质参数拆分到了一个单独的类中,现在每个物体实
    例都可以拥有自己的材质。
  • Intersection.hpp: 这个数据结构包含了相交相关的信息。
  • Ray.hpp: 光线类,包含一条光的源头、方向、传递时间 t 和范围 range.
  • Bounds3.hpp: 包围盒类,每个包围盒可由 pMinpMax 两点描述(请思考为什么)。Bounds3::Union 函数的作用是将两个包围盒并成更大的包围盒。与材质一样,场景中的每个物体实例都有自己的包围盒.- BVH.hpp: BVH 加速类。场景 scene 拥有一个 BVHAccel 实例。从根节点开始,我们可以递归地从物体列表构造场景的 BVH.

代码框架解析

作业6的代码框架与作业5大体一致,主要有几处修改的地方。

一是将光线的源头、方向、传递时间 t 和范围 range 封装成了 Ray 光线类。因此,与光线有关的代码,需要有所修改。

struct Ray{
    //Destination = origin + t*direction
    Vector3f origin;
    Vector3f direction, direction_inv;
    double t;//transportation time,
    double t_min, t_max;

    Ray(const Vector3f& ori, const Vector3f& dir, const double _t = 0.0): origin(ori), direction(dir),t(_t) {
        direction_inv = Vector3f(1./direction.x, 1./direction.y, 1./direction.z);
        t_min = 0.0;
        t_max = std::numeric_limits<double>::max();

    }

    Vector3f operator()(double t) const{return origin+direction*t;}

    friend std::ostream &operator<<(std::ostream& os, const Ray& r){
        os<<"[origin:="<<r.origin<<", direction="<<r.direction<<", time="<< r.t<<"]\n";
        return os;
    }
};

二是交点也被封装成了 Intersection 结构体,其主要包括,是否相交happened、交点坐标coords、交点法向量normal、交点与相机距离 distance,相交物体 obj,材质 m

struct Intersection
{
    Intersection(){
        happened=false;
        coords=Vector3f();
        normal=Vector3f();
        distance= std::numeric_limits<double>::max();
        obj =nullptr;
        m=nullptr;
    }
    bool happened;
    Vector3f coords;
    Vector3f normal;
    double distance;
    Object* obj;
    Material* m;
};
#en

三是因为引入了BVH,所以在判断光线与物体是否有交点前,先对光线与包围框是否有交点进行判断。这里需要注意的是,空间中只有一个物体,而物体是由多个三角形网格组成。所以这里其实要构建两个 BVH,一个是针对空间中物体的 BVH,另一个是针对单一物体的组成网格的 BVH。运行代码,也可以看出,执行了两次 BVH 的构建。
那么在之后判断光线与物体中网格相交实际要执行三个步骤,第一步是判断光线是否与该物体包围框相交,如果相交,第二步判断是否与该物体中的三角形相交,而第二步由可以拆分成两个步骤,一是判断是否与三角形所在包围框相交,如果相交,二就是判断是否与三角形相交。

下面主要说明一下 BVH 的构建过程:

  • 首先对于所有物体寻找一个包围框
  • 递归地将物体的集合划分进两个子集中,然后重新计算子集的包围框
  • 在合适的时候停止,比如一个包围框的物体个数为5个(代码里为1),在叶子节点存储物体信息
    那么如何划分节点,一般地做法是选择一个维度去划分,一般是选最长的轴进行划分。在物体的中间位置划分节点 —— 每个子集有相同数量的物体,其对应的代码部分如下
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);

作业代码

Render

需要注意的地方在于,eyepos 的位置在 ( − 1 , 5 , 10 ) (-1, 5, 10) (1,5,10)。所以其实屏幕的中心坐标就不在 ( 0 , 0 , − 1 ) (0,0,-1) (0,0,1) 了,而是位于 ( − 1 , 5 , 9 ) (-1,5,9) (1,5,9)。所以其实只是相机源点发生了改变,而射线方向其实与作业5求法相同。

Vector3f dir = normalize(Vector3f(x, y, -1));
Ray ray(eye_pos, dir);
framebuffer[m++] = scene.castRay(ray, 0);

Triangle::getIntersection

因为将光线封装为 Ray 类,所以相应光线-三角形相交函数也需要有所修改:

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.origin + t_tmp * ray.direction;
    inter.normal = this->normal;
    inter.distance = dotProduct(t_tmp * ray.direction, t_tmp * ray.direction);
    inter.obj = this;
    inter.m = this->m;

    return inter;
}

IntersectP

这个函数的作用是判断包围盒 BoundingBox 与光线是否相交。其原理为:求光线与轴对齐包围盒交点 —— 光线进入所有平面内才算进入包围盒内,光线只要出了一个平面就算出了包围盒,所以 t e n t e r = m a x ( t m i n ) t_{enter}=max(t_{min}) tenter=max(tmin), t e x i t = m i n ( t m a x ) t_{exit}=min(t_{max}) texit=min(tmax) ,只要 t e n t e r < t e x i t & & t e x i t > = 0 t_{enter}<t_{exit} \&\& t_{exit}>=0 tenter<texit&&texit>=0 要求 t e x i t > = 0 t_{exit}>=0 texit>=0 是因为 t e x i t < 0 t_{exit}<0 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
    float t_enter;
    float t_exit;
    Vector3f t_enter_v3f = (pMin - ray.origin) * invDir;
    Vector3f t_exit_v3f = (pMax - ray.origin) * invDir;

    // for(int i=0; i<3; ++i) {
    //     if(!dirIsNeg[i])
    //         std::swap(t_enter_v3f[i], t_exit_v3f[i]);
    // }
    if(!dirIsNeg[0])
        std::swap(t_enter_v3f.x, t_exit_v3f.x);
    if(!dirIsNeg[1])
        std::swap(t_enter_v3f.y, t_exit_v3f.y);
    if(!dirIsNeg[2])
        std::swap(t_enter_v3f.z, t_exit_v3f.z);

    t_enter = std::max(t_enter_v3f.x, std::max(t_enter_v3f.y, t_enter_v3f.z));
    t_exit = std::min(t_exit_v3f.x, std::min(t_exit_v3f.y, t_exit_v3f.z));

    if (t_enter < t_exit && t_exit >= 0)
        return true;
    else
        return false;       
}

这里需要注意的是,如果一个轴光线的方向为负数,那么其离包围盒远端距离越小,所以这与光线方向为正,正好相反,所以需要 swap 一下

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 t_enter;
    float t_exit;
    Vector3f t_enter_v3f = (pMin - ray.origin) * invDir;
    Vector3f t_exit_v3f = (pMax - ray.origin) * invDir;

    // for(int i=0; i<3; ++i) {
    //     if(!dirIsNeg[i])
    //         std::swap(t_enter_v3f[i], t_exit_v3f[i]);
    // }
    if(!dirIsNeg[0])
        std::swap(t_enter_v3f.x, t_exit_v3f.x);
    if(!dirIsNeg[1])
        std::swap(t_enter_v3f.y, t_exit_v3f.y);
    if(!dirIsNeg[2])
        std::swap(t_enter_v3f.z, t_exit_v3f.z);

    t_enter = std::max(t_enter_v3f.x, std::max(t_enter_v3f.y, t_enter_v3f.z));
    t_exit = std::min(t_exit_v3f.x, std::min(t_exit_v3f.y, t_exit_v3f.z));

    if (t_enter < t_exit && t_exit >= 0)
        return true;
    else
        return false;       
}

getIntersection

其算法伪代码如下:

其思想就是如果光线与该节点包围盒没有交点,那么必定与包围盒内的物体没有交点。

如果有交点,如果当前节点是叶节点,那么就与包围盒内的物体判断是否有交点。如果不是叶节点,即中间节点,那么就递归其左右子树,取一个最近的交点返回。

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    Intersection isect;
    // TODO Traverse the BVH to find intersection
    if(!node->bounds.IntersectP(ray, ray.direction_inv, std::array<int, 3> {ray.direction.x>0, ray.direction.y>0, ray.direction.z>0}))
        return isect;
    if(node->object != nullptr) // leaf node
        return node->object->getIntersection(ray);

    // internal node
    Intersection isect_left, isect_right;
    isect_left = getIntersection(node->left, ray);
    isect_right = getIntersection(node->right, ray);

    return isect_left.distance <= isect_right.distance ? isect_left : isect_right;
}

结果如下图所示


进阶代码

进阶代码要求实现 SAH(Surface Area Heuristic),其具体原理可以看这篇文章

简单来说如果图元分布不均匀,BVH 得到的划分结果会很差,比如下面的左图包围框有很大的重叠面积,那么射线与这两个包围框相交的概率肯定会比右图的大,计算就会增加。

SAH 就是考虑求交代价与包围盒面积,给出一个总的代价公式,我们寻找划分位置使得这个总代价最小。

这里的 S ( A ) S(A) S(A) 表示左子树包围框面积, S ( B ) S(B) S(B) 表示右子树包围框面积, S ( C ) S(C) S(C) 表示当前节点包围框面积, a a a表示左子树物体个数(认为求交代价为1), b b b表示右子树物体个数, 0.125 0.125 0.125 为遍历划分代价。

其代码如下

BVHBuildNode* BVHAccel::recursiveBuild_SAH(std::vector<Object*> objects)
{
    BVHBuildNode* node = new BVHBuildNode();
    // std::cout<<objects.size()<<std::endl;

    // Compute bounds of all primitives in SAH node
    Bounds3 bounds;
    for (int i = 0; i < objects.size(); ++i)
        bounds = Union(bounds, objects[i]->getBounds());
    if (objects.size() == 1) {
        // Create leaf _SAHBuildNode_
        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_SAH(std::vector{objects[0]});
        node->right = recursiveBuild_SAH(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;
        }

        int B_size = 10;
        int minCostIdx = 1;
        float minCost = std::numeric_limits<float>::infinity();
        float SC = bounds.SurfaceArea();
        // float SC = centroidBounds.SurfaceArea();
        
        auto beginning = objects.begin();
        auto ending = objects.end();
        
        for(int i=1; i<B_size; ++i) {// search minCost and minCostIdx
            auto middling = objects.begin() + std::max(1, (int)(objects.size() * i / B_size));

            auto leftshapes = std::vector<Object*>(beginning, middling);
            auto rightshapes = std::vector<Object*>(middling, ending);

            assert(objects.size() == (leftshapes.size() + rightshapes.size()));

            Bounds3 leftbounds, rightbounds;

            for(int j=0; j<leftshapes.size(); ++j) {
                leftbounds = Union(leftbounds, leftshapes[j]->getBounds());
                // leftbounds = Union(leftbounds, leftshapes[j]->getBounds().Centroid());
            }

            for(int j=0; j<rightshapes.size(); ++j) {
                rightbounds = Union(rightbounds, rightshapes[j]->getBounds());
                // rightbounds = Union(rightbounds, rightshapes[j]->getBounds().Centroid());
            }

            float SA = leftbounds.SurfaceArea();
            float SB = rightbounds.SurfaceArea();
            float c_A_B = leftshapes.size() * SA / SC + rightshapes.size() * SB / SC + 0.125f;

            if(c_A_B < minCost) {
                minCost = c_A_B;
                minCostIdx = i;
            }
        }

        auto middling_minCost = objects.begin() + std::max(1, ((int)objects.size() * minCostIdx / B_size));
        auto leftshapes_minCost = std::vector<Object*>(beginning, middling_minCost);
        auto rightshapes_minCost = std::vector<Object*>(middling_minCost, ending);

        assert(objects.size() == (leftshapes_minCost.size() + rightshapes_minCost.size()));
        
        // std::cout<<leftshapes_minCost.size()<<" "<<rightshapes_minCost.size()<<std::endl;
        node->left = recursiveBuild_SAH(leftshapes_minCost);
        node->right = recursiveBuild_SAH(rightshapes_minCost);

        node->bounds = Union(node->left->bounds, node->right->bounds);
    }

    return node;
}

这里考虑九个划分位置,即左边 1/10个物体,右边 9/10,左边 2/10,右边 8/10 …,然后面积我这里用的是整个包围框的面积,也有用质心包围框的面积,我感觉都可以,下面比较下两种方法的时间,可见 SAH 比 BVH 快了 1秒。

  • 8
    点赞
  • 12
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值