基础题部分

根据教程PDF，首先需要引用如下函数（在作业5的基础上稍作修改）：

renderer() in Renderer.cpp

解说见注释：

// The main render function. This where we iterate over all pixels in the image,
// generate primary rays and cast these rays into the scene. The content of the
// framebuffer is saved to a file.
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
            // 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*
            float x = (2 * (i + 0.5) / (float)scene.width - 1) * imageAspectRatio * scale;
            float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;
            
            // 从这里开始
            Vector3f dir = normalize(Vector3f(x, y, -1)); // Don't forget to normalize this direction!

            // 在本次代码框架中，castRay()函数放在了Scene类中
            // 并且castRay()的传入参数修改为(const Ray &ray, int depth)
            // 所以要先声明Ray ray(const Vector3f& ori, const Vector3f& dir)
            Ray ray(eye_pos, dir);
            
            // 调用scene.castRay()
            // 并将返回结果(Vector3f类型的像素RGB颜色)存入framebuffer中
            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);    
}

Triangle::getIntersection in Triangle.hpp

解说见注释：

// TODO find ray triangle intersection
inline Intersection Triangle::getIntersection(Ray ray)
{
    Intersection inter;

    if (dotProduct(ray.direction, normal) > 0)
        return inter;
    
    // u,v为三角形重心坐标，即PPT中的b1,b2
    // t_tmp是PPT中的t
    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;

    // 如果有intersection，则赋值给Intersection inter并将inter返回
    // Intersection的定义：见Intersection.hpp
    inter.happened = true; // 光线和该三角形有交叉点
    inter.coords = ray(t_tmp); // 利用公式o+td，求得intersection（交点）坐标
    inter.normal = normal; // 传入该三角形的法线
    inter.distance = t_tmp; // 与光源点的距离，用t表示即可
    inter.obj = this; // 指针指向自身
    inter.m = m; // 传入该三角形材质(material)

    return inter;
}

然后实现BVH的功能：

IntersectP() in the Bounds3.hpp

这里使用PPT中的公式：

如果光线方向是相反的（即xyz轴相应的方向向量为负数），则要交换tmin和tmax，这里闫老师的PPT没讲，需要特别注意！

代码和详解如下：

// TODO test if ray bound intersects
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
    // invDir就是1/Dir，因为乘法的运算比除法快，故使用invDir来加快运算速度
    // dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x<0),int(y<0),int(z<0)], use this to simplify your logic
    // 上面这一句注释在源代码中写错了，应该是[int(x<0),int(y<0),int(z<0)]

    float tmin_x = (pMin.x - ray.origin.x) * invDir.x;
    float tmin_y = (pMin.y - ray.origin.y) * invDir.y;
    float tmin_z = (pMin.z - ray.origin.z) * invDir.z;

    float tmax_x = (pMax.x - ray.origin.x) * invDir.x;
    float tmax_y = (pMax.y - ray.origin.y) * invDir.y;
    float tmax_z = (pMax.z - ray.origin.z) * invDir.z;

    // swap t_min and t_max if direction is negative
    // 如果[int(x<0),int(y<0),int(z<0)]，需要交换tmin和tmax
    // 这里PPT上面没有讲，需要特别注意
    if (dirIsNeg[0]) std::swap(tmin_x, tmax_x);
    if (dirIsNeg[1]) std::swap(tmin_y, tmax_y);
    if (dirIsNeg[2]) std::swap(tmin_z, tmax_z);

    float t_enter = std::max(tmin_x, std::max(tmin_y, tmin_z));
    float t_exit = std::min(tmax_x, std::min(tmax_y, tmax_z));

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

getIntersection()in BVH.cpp

这里参考PPT中的伪代码即可轻松完成：

代码如下： 

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

    Intersection intersection; // default: happened = false

    // if (ray misses node.bbox) return;
    if (!flag_insert) return intersection;

    // if (node is a leaf node)
    if (node->left == nullptr && node->right == nullptr)
    {
        // test intersection with all objs;
        // return closest intersection;
        return node->object->getIntersection(ray);
    }

    // hit1 = Intersect(ray, node.child1);
    Intersection hit1 = getIntersection(node->left, ray);
    // hit2 = Intersect(ray, node.child2);
    Intersection hit2 = getIntersection(node->right, ray);
    // return the closer of hit1, hit2
    return (hit1.distance < hit2.distance ? hit1 : hit2);
}

渲染结果如下：

 BVH用时：4秒钟

附加题部分：SAH加速的实现

这里主要参考了如下资料：

闫老师给的CMU PPT链接（完整PPT可供下载）

中文讲解：BVH with SAH (Bounding Volume Hierarchy with Surface Area Heuristic) - lookof - 博客园

这里主要讲几个重点，其中，成本cost的计算方式为：

伪代码：

 中文讲解为：

将以上步骤简化：

cost(A,B) = p(A)*(n in A) + p(B)*(m in B)

其中：

p(A) = s(A) / s(C)

p(A) = s(B) / s(C)

因为分母同为s(C)，所以我们对比的时候无需对比分母，故cost(A,B)可以简化为：

cost(A,B) = s(A)*(n in A) + s(B)*(m in B)

分别在xyz轴中寻找上式的最小值即可。

代码实现：

首先，往BVH类的声明（BVH.hpp）中加入以下函数：

BVHBuildNode* BVHAccel::recursiveSAHBuild(std::vector<Object*> objects);

然后，在BVH.cpp中：

1）声明一个结构体，用于记录cost（划分成本的计算结果）和index（在何处划分）

// my struct
struct Cost
{
    int index;
    double cost_value;
};

2）声明一个函数compute_cost，用于计算xyz轴中按照任意一个轴进行划分时的最小成本，返回一个Cost类型的变量（划分方案&最小成本）

// my func
// 参数解释：
// objects为未划分的primitives（三角形）列表
// b为分成几个buckets
// p为一个bucket中有几个primitives（即三角形）
// dim为维度，0为x轴，1为y轴，2为z轴
Cost compute_cost(std::vector<Object*> objects, int b, int p, int dim)
{
    // 先将xyz轴按照升序排列
    switch (dim)
    {
    case 0: // for x axis
        std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
            return f1->getBounds().Centroid().x <
                   f2->getBounds().Centroid().x;
        });
        break;
    
    case 1: // for y axis
        std::sort(objects.begin(), objects.end(), [](auto f1, auto f2){
            return f1->getBounds().Centroid().y <
                   f2->getBounds().Centroid().y;
        });
        break;
    
    case 2: // for z axis
        std::sort(objects.begin(), objects.end(), [](auto f1, auto f2){
            return f1->getBounds().Centroid().z <
                   f2->getBounds().Centroid().z;
        });
        break;
    }

    // 记录cost的列表
    std::vector<Cost> cost_list;

    // 循环计算不同划分方案的成本
    Bounds3 bucketA;
    // 将buckets从小到大逐个加入划分方案
    for (int i = 0; i < b-1; ++i) {
        for (int j = 0; j < p; ++j) {
            // 将bucket中的元素从小到大加入包围盒A中
            bucketA = Union(bucketA, objects[p*i+j]->getBounds());
        }
        // 将剩余的元素从小到大加入包围盒B中
        Bounds3 bucketB;
        for (int k = p*(i+1); k < objects.size(); ++k) {
            bucketB = Union(bucketB, objects[k]->getBounds());
        }
        // 计算成本
        double sA = bucketA.SurfaceArea(); // 获取包围盒A的表面积
        double sB = bucketB.SurfaceArea();
        double nA = (i+1)*p; // 包围盒A中的primitives（三角形）的个数
        double nB = objects.size() - nA;
        double cost_val = sA*nA + sB*nB; // compute cost
        Cost cost;
        cost.index = p*(i+1)-1;
        cost.cost_value = cost_val;
        cost_list.emplace_back(cost);
    }

    // 按照升序排列（第一个为cost的最小值）
    std::sort(cost_list.begin(), cost_list.end(), [](auto f1, auto f2){
        return f1.cost_value < f2.cost_value;
    });

    // 返回cost最小值（和对应的划分方案index）
    return cost_list[0];
}

3）声明一个函数get_partition_axis，对比xyz轴的各自最小划分成本，方便后续选择划分成本最小的轴进行划分：

// my func
Cost get_partition_axis(Cost cost_x, Cost cost_y, Cost cost_z, int& dim)
{
    if (cost_x.cost_value < cost_y.cost_value && cost_x.cost_value < cost_z.cost_value) {
        dim = 0;
        return cost_x;
    } else if (cost_y.cost_value < cost_z.cost_value) {
        dim = 1;
        return cost_y;
    } else {
        dim = 2;
        return cost_z;
    }
}

4）最后是BVHAccel::recursiveSAHBuild函数的实现：

BVHBuildNode* BVHAccel::recursiveSAHBuild(std::vector<Object*> objects)
{
    BVHBuildNode* node = new BVHBuildNode();

    // num of primitives in a bucket
    // cannot be smaller than 2!
    int p = 5;

    if (objects.size() <= p) {
        node = recursiveBuild(objects);
        return node;
    }

    int b = objects.size() / p; // num of buckets
    if (objects.size() % p != 0) b++;

    Cost cost_x;
    Cost cost_y;
    Cost cost_z;
    Cost final_cost;
    int dim = 0; // 默认按照x轴划分，防bug
    cost_x = compute_cost(objects, b, p, 0);
    cost_y = compute_cost(objects, b, p, 1);
    cost_z = compute_cost(objects, b, p, 2);
    final_cost = get_partition_axis(cost_x, cost_y, cost_z, dim);

    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() + final_cost.index; // SAH changed
    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 = recursiveSAHBuild(leftshapes);
    node->right = recursiveSAHBuild(rightshapes);
    node->bounds = Union(node->left->bounds, node->right->bounds);

    return node;
}

需要运行SAH的时候，记得需要修改如下两处：

第一处，修改默认splitMethod的值：

// BVHAccel Declarations
inline int leafNodes, totalLeafNodes, totalPrimitives, interiorNodes;
class BVHAccel {

public:
    // BVHAccel Public Types
    enum class SplitMethod { NAIVE, SAH };

    // 从这里开始
    // BVHAccel Public Methods
    // BVHAccel(std::vector<Object*> p, int maxPrimsInNode = 1, SplitMethod splitMethod = SplitMethod::NAIVE);
    // 把上面注释掉的那一句，改成下面这一句
    BVHAccel(std::vector<Object*> p, int maxPrimsInNode = 1, SplitMethod splitMethod = SplitMethod::SAH);
    // 到这里为止
    Bounds3 WorldBound() const;
    ~BVHAccel();

    Intersection Intersect(const Ray &ray) const;
    Intersection getIntersection(BVHBuildNode* node, const Ray& ray)const;
    bool IntersectP(const Ray &ray) const;
    BVHBuildNode* root;

    // BVHAccel Private Methods
    BVHBuildNode* recursiveBuild(std::vector<Object*>objects);
    BVHBuildNode* recursiveSAHBuild(std::vector<Object*>objects); // my func

    // BVHAccel Private Data
    const int maxPrimsInNode;
    const SplitMethod splitMethod;
    std::vector<Object*> primitives;
};

第二处，在BVHAccel::BVHAccel() （BVH.cpp）中：

BVHAccel::BVHAccel(std::vector<Object*> p, int maxPrimsInNode, SplitMethod splitMethod)
    : maxPrimsInNode(std::min(255, maxPrimsInNode)), splitMethod(splitMethod),
      primitives(std::move(p))
{
    time_t start, stop;
    time(&start);
    if (primitives.empty())
        return;
    
    // 从这里开始
    if (splitMethod == SplitMethod::NAIVE) {
        root = recursiveBuild(primitives);
    } else if (splitMethod == SplitMethod::SAH) {
        root = recursiveSAHBuild(primitives); // for SAH
    }
    // 到这里结束
    
    time(&stop);
    double diff = difftime(stop, start);
    int hrs = (int)diff / 3600;
    int mins = ((int)diff / 60) - (hrs * 60);
    int secs = (int)diff - (hrs * 3600) - (mins * 60);

    printf(
        "\rBVH Generation complete: \nTime Taken: %i hrs, %i mins, %i secs\n\n",
        hrs, mins, secs);
}

然后大功告成，可以编译运行啦！

SAH运行之后的结果：耗时3秒（比BVH的4秒少1秒）