Bounding Volume Hierarchies

本章节将会是光线追踪系列博客中目前为止最难且最复杂的部分，我们之所以在当前的进度下引入BVH，一方面是为了让代码运行更快，同时也是因为实现BVH要求我们重构hittable类的部分代码，当我们在场景中添加长方体等其他物体时，我们就不再需要再回过头考虑hittable的重构了。

当前我们的渲染器的主要时间瓶颈在光线-物体的相交测试上，当我们向场景中发射数以百万计的光线时，我们可以首先对场景中的物体进行排序，这里的排序主要使用两种方法：细分空间、细分对象。后者通常更容易通过代码实现，在大多数模型上运行速度也一样快。

最终，我们将通过BVH将渲染器的复杂度从线性优化为次线性sublinear，也就是O(log n)。

3.1 The Key Idea

为一组primitive创建bounding volume的关键点在于，找到一个可以完全包围所有物体的volume。假设，我们计算出了一个能够完全包围十个物体的球体，那么我们可以得出这样的结论：任何无法与该球体相交的光线，都不可能与球体所包围的物体相交。我们可以用一段伪代码来说明这个结论：

if (ray hits bounding object)
   return whether ray hits bounded objects
else
   return false

3.2 Hierachies of Bounding Volumes

为了实现次线性的复杂度，我们需要为bounding volumes构建层级。例如说，如果我们将物体分为红色与蓝色的两组，然后使用矩形的bounding volume，则我们得到的结果如下所示：

从图中我们可以看出，红色与蓝色两个组都属于紫色bounding volume的一部分，且两个volume可以发生重叠，并且它们之间没有次序之分，我们只会使用左右来区分它们。用伪代码的形式表示：

if (hits purple)
    hit0 = hits blue enclosed objects
    hit1 = hits red enclosed objects
    if (hit0 or hit1)
    	return true and info of closer hit
return false

3.3 Axis-Aligned Bounding Boxes(AABBs)

在BVH中，光线与物体进行相交测试之前，会首先与bounding volume进行相交测试，在众多的volume模型中，AABB的相交测试的性能优势更为明显，同时也可以更好地实现对物体的划分。

值得一提的是，在光线与bounding volume的相交测试中，我们只需要判断是否相交会发生，并不需要再获取任何额外的信息。

我们使用slab的方法来进行AABBs的相交测试。我们以2D空间中的AABB（也就是一个矩形）为例来说明这个方法。从区间的角度来考虑，一个矩形可以通过两组区间来表示，如下图所示，绿色与蓝色所代表的两个范围共同构成了这样一个矩形：

如果一个光线能够与某个区间相交，则该光线必定与区间的边界相交。如下图所示，光线与x₀和x₁这两个区间的边界相交，会分别返回t₀与t₁：

我们再次回到三维空间中思考问题，那么x₀和x₁所表示的就是三维空间中的平面。我们如何判断光线是否与平面相交呢？实际上还是利用光线与平面的参数方程所构造的等式，我们首先回顾光线的参数方程：

\[\mathbf{P}(t) = \mathbf{Q} + t \mathbf{d}\]

这个等式对于xyz三个分量都是成立的，我们不妨以x坐标值为例，当光线与x₀所代表的平面相交时，我们有：

\[x_0 = Q_x + t_0 d_x\]

其中t₀是等式中的未知数，即：

\[t_0 = \frac{x_0 - Q_x}{d_x}\]

同样的，我们可得：

\[t_1 = \frac{x_1 - Q_x}{d_x}\]

当我们为每个区间计算出t₀与t₁时，我们就可以利用slab方法中的核心技巧了：只有当光线与bounding box相交时，构成bounding box的所有区间所对应的[t₀, t₁]才会发生重叠。如下图中更下方的射线所示：

3.4 Ray Intersection with an AABB

我们将前面所提到的实现思路用下面的伪代码表示：

intervalX <- computeIntersectionX(ray, x0, x1);
intervalY <- computeIntersectionX(ray, y0, y1);
intervalZ <- computeIntersectionX(ray, z0, z1);
return overlaps(intervalX, intervalY, intervalZ);

只是我们还需要考虑到光线的行进方向，还是以X轴为例，如果光线朝着-X的方向前进，则我们会得到t₀,>t₁，所以我们还需要进行如下判断（对所有轴向上的分量均成立）：

\[\displaylines{t_{x0} = \min( \frac{x_0 - Q_x}{d_x}, \frac{x_1 - Q_x}{d_x})\\t_{x1} = \max( \frac{x_0 - Q_x}{d_x}, \frac{x_1 - Q_x}{d_x})}\]

接下来，我们再来看看伪代码中overlaps()函数的实现：

bool overlaps(tInterval1, tInterval2)
    tMin <- max(tInterval1.min, tInterval2.min)
    tMax <- min(tInterval1.max, tInterval2.max)
    return tMin < tMax

现在有了清晰的思路，我们就可以动手为我们的渲染器实现BVH了。

首先需要做的是实现一个AABB的类，按照前面的内容，我们可以将AABB视为由三个轴向上的三个区间围成的一个box。此外，AABB类还需要一个hit()函数，用于判断AABB是否会与给定光线发生相交。下面是AABB类的代码：

#ifndef AABB_H
#define AABB_H

#include <queue>

#include "rayTracing.h"

class aabb
{
public:
    interval x, y, z;

    aabb() = default; // aabb is empty by default

    aabb(const interval& x, const interval& y, const interval& z)
        : x(x), y(y), z(z) {}

    aabb(const point3& a, const point3& b)
    {
        // here we treat a and b as extrema for the bounding box, and we sort them manually
        // thus we just dont require them to be particular minimum-maximum order
        x = a[0] <= b[0] ? interval(a[0], b[0]) : interval(b[0], a[0]);
        y = a[1] <= b[1] ? interval(a[1], b[1]) : interval(b[1], a[1]);
        z = a[2] <= b[2] ? interval(a[2], b[2]) : interval(b[2], a[2]);
    }

    const interval& getIntervalOfAxis(int n) const
    {
        if (n == 1) return y;
        if (n == 2) return z;
        return x;
    }

    bool hit(const ray& rayIncoming, interval t) const
    {
        const point3& origin = rayIncoming.origin();
        const vec3& direction = rayIncoming.direction();

        for (int axis = 0; axis < 3; axis++)
        {
            // 1. calculate t0 and t1 for the intersection of incoming ray and certain interval of AABB
            const interval& axisInterval = getIntervalOfAxis(axis);
            const double divisor = 1.0 / direction[axis];
            double t0 = (axisInterval.min - origin[axis]) * divisor;
            double t1 = (axisInterval.max - origin[axis]) * divisor;
            // 2. make sure t0 < t1
            if (t0 >= t1) std::swap(t0, t1);
            // 3. do the overlap test
            // for AABB, there are three [t0, t1] intervals
            // we name the maximum t0 of the intervals as tMin,
            // and the minimum t1 of the interval as tMax
            if (t0 > t.min) t.min = t0;
            if (t1 < t.max) t.max = t1;
            // 4. if tMin > tMax, then ray will hit AABB
            if (t.max <= t.min) return false;
        }
        return true;
    }
};

#endif

3.5 Constructing Bounding Boxes for Hittables

BVH的核心在于，为场景中所有物体的AABB构建一个层级，每个单一的hittable对象都属于这个层级中的一个子节点，或者说是BVH这棵大树的叶子。

由于我们的BVH采用的是划分物体的方法，所以场景中每一个hittable对象都需要一个AABB，为此，我们需要给hittable类的派生类添加新的成员变量，并且在sphere类的构造函数中计算AABB：

class sphere final : public hittable
{
	...
private:
	...
	aabb bbox;
}

由于光线追踪器已经引入了运动模糊的特性，静态球体与动态球体计算AABB的方式略有不同，前者只需要根据球体的球心位置与半径计算即可，而后者我们则需要先分别计算出球体在time=0和time=1时的AABB，再根据这两个AABB计算出一个可以包围这两个AABB的AABB。听起来我们需要一个新的AABB构造函数了。

由于AABB是由三个轴上的区间构成的，所以新的AABB构造函数同样需要一个新的interval构造函数，该构造函数会根据输入的两个区间输出一个包含两个区间的新区间：

interval(const interval& a, const interval& b)
{
    // create the interval tightly enclosing the two input intervals
    min = a.min <= b.min ? a.min : b.min;
    max = a.max >= b.max ? a.max : b.max;
}

然后是新的AABB构造函数：

aabb(const aabb& a, const aabb& b)
{
x = interval(a.x, b.x);
y = interval(a.y, b.y);
z = interval(a.z, b.z);
}

现在我们继续在sphere类的构造函数中计算球体的AABB：

// stationary sphere
sphere(const point3& center, double radius, const shared_ptr<material>& material)
    : initialCenter(center), radius(fmax(0, radius)), material(material), isMoving(false)
{
    vec3 radiusVec3 = vec3(radius, radius, radius);
    bbox = aabb(initialCenter - radiusVec3, initialCenter + radiusVec3);
}

// moving sphere
sphere(const point3& initialCenter, const point3& finalCenter, double radius, const shared_ptr<material>& material)
    : initialCenter(initialCenter), radius(fmax(0, radius)), material(material), isMoving(true)
{
    vec3 radiusVec3 = vec3(radius, radius, radius);
    aabb initialAABB = aabb(initialCenter - radiusVec3, initialCenter + radiusVec3);
    aabb finalAABB = aabb(finalCenter - radiusVec3, finalCenter + radiusVec3);
    bbox = aabb(initialAABB, finalAABB);

    centerVector = finalCenter - initialCenter;
}

最后，为了能够使用hittable对象的AABB，我们在hittable类中新增一个虚函数boundingBox()：

...

#include "aabb.h"

class material;

...

class hittable
{
public:
    virtual ~hittable() = default;

    virtual bool hit(const ray& r, interval tInterval, hitInfo& info) const = 0;
    
    virtual aabb boundingBox() const = 0;
};

sphere类的中的boundingBox()只需要返回成员变量bbox即可：

class sphere final : public hittable
{
	...

    aabb boundingBox() const override {return bbox;}

private:
    ...
    aabb bbox;

    ...
};

3.6 Creating Bounding Boxed of Lists of Objects

由于hittable类新增了boundingBox()函数，我们同样更新hittableList类。由于BVH是一个层级结构，场景中的hittableList对象将作为这个结构的根节点，同样具有一个AABB，且场景中每增加一个hittable对象，hittableList的AABB都会根据新增的hittable的AABB进行范围上的扩展。

#include "aabb.h"
#include "hittable.h"

#include <vector>

class hittableList final : public hittable
{
public:
	...

    void add(const shared_ptr<hittable>& object)
    {
        objects.push_back(object);
        bbox = aabb(bbox, object->boundingBox());
    }

	..

    aabb boundingBox() const override { return bbox; }

private:
    aabb bbox;
};

3.7 The BVH Node Class

BVH在我们的C++工程里同样派生自hittable类，并且在我们的代码设计中，bvh类可以同时包含根节点与子节点：

#ifndef BVH_H
#define BVH_H

#include "rayTracing.h"

#include "aabb.h"
#include "hittable.h"
#include "hittableList.h"

class bvhNode final : public hittable
{
public:
    bvhNode(hittableList list) : bvhNode(list.objects, 0, list.objects.size())
    {
        // There's a C++ subtlety here. This constructor (without span indices) creates an
        // implicit copy of the hittable list, which we will modify. The lifetime of the copied
        // list only extends until this constructor exits. That's OK, because we only need to
        // persist the resulting bounding volume hierarchy.
    }

    bvhNode(std::vector<shared_ptr<hittable>>& objects, size_t start, size_t end)
    {
        // To be implemented later
    }

    bool hit(const ray& r, interval tInterval, hitInfo& info) const override
    {
        if (!bbox.hit(r, tInterval)) return false;

        bool hitLeft = left->hit(r, tInterval, info);
        bool hitRight = right->hit(r, interval(tInterval.min, hitLeft ? info.t : tInterval.max), info);

        return hitLeft || hitRight;
    }

    aabb boundingBox() const override {return bbox;}

private:
    shared_ptr<hittable> left;
    shared_ptr<hittable> right;
    aabb bbox;
};

#endif

3.8 Splitting BVH Volumes

在BVH加速结构中，最复杂的部分就是构建BVH，现在让我们在bvh类的构造函数中完成这项任务。在BVH中，我们需要找到一个尽可能合理的递归划分策略，所以我们的思路是：

• 随机选择一个轴向
• 使用std::sort对物体进行排序
• 将排序后的物体分为两个subtree，每个subtree各包含一半的物体

由于构建BVH的过程是递归的，且涉及到划分，我们需要考虑到一些特殊的情况。当传入的hittableList只有两个元素时，我们直接在两个subtree中各放一个元素然后结束递归。代码中所使用的算法应该尽可能保证平滑进行，所以我们需要避免对null指针的检查，因此，若hittableList中只有一个元素，我们就为每个subtree复制一个元素。

当然，在构建BVH时，如果我们能够显式地检查到hittableList中只包含三个元素（也就是场景中只有三个物体），我们只进行一次递归调用就可以了。只是对于这种特殊情况的优化我们可以放在后面进行。

随机选择轴向的功能需要我们添加一个随机生成整数值的函数：

inline int randomInt(int min, int max)
{
    // returns a random integer in [min, max]
    return static_cast<int>(randomDouble(min, max + 1));
}

下面是用于构建BVH的构造函数，其中boxCompareX、boxCompareY、boxCompareZ我们还没有定义：

bvhNode(std::vector<shared_ptr<hittable>>& objects, size_t start, size_t end)
{
    int axis = randomInt(0, 2);

    auto comparator = axis == 0 ? boxCompareX
        			: axis == 1 ? boxCompareY
            					: boxCompareZ;

    size_t objectSpan = end - start;

    if (objectSpan == 1)
    {
        left = right = objects[start];
    }
    else if (objectSpan == 2)
    {
        left = objects[start];
        right = objects[start + 1];
    }
    else
    {
        std::sort(objects.begin() + start, objects.begin() + end, comparator);
        auto middle = start + objectSpan / 2;
        left = make_shared<bvhNode>(objects, start, middle);
        right = make_shared<bvhNode>(objects, middle, end);
    }

    bbox = aabb(left->boundingBox(), right->boundingBox());
}

3.9 The Box Comparison Function

现在，我们需要定义出std::sort所使用比较函数：

class bvhNode final : public hittable
{
...

private:
    shared_ptr<hittable> left;
    shared_ptr<hittable> right;
    aabb bbox;

    static bool boxCompare(
        const shared_ptr<hittable> a, const shared_ptr<hittable> b, int axisIndex)
    {
        interval aInterval = a->boundingBox().getIntervalOfAxis(axisIndex);
        interval bInterval = b->boundingBox().getIntervalOfAxis(axisIndex);
        return aInterval.min < bInterval.min;
    }

    static bool boxCompareX(const shared_ptr<hittable> a, const shared_ptr<hittable> b)
    {
        return boxCompare(a, b, 0);
    }

    static bool boxCompareY(const shared_ptr<hittable> a, const shared_ptr<hittable> b)
    {
        return boxCompare(a, b, 1);
    }

    static bool boxCompareZ(const shared_ptr<hittable> a, const shared_ptr<hittable> b)
    {
        return boxCompare(a, b, 2);
    }
};

当我们进行到这一步，我们已经可以在场景中使用BVH来加速光线追踪的计算了：

#include "rayTracing.h"

#include "bvh.h"
#include "camera.h"
#include "hittable.h"
#include "hittableList.h"
#include "material.h"
#include "sphere.h"

int main() {

    // World-------------------------------------------------------------------------------------
    hittableList world;

    ...

    // build BVHs
    world = hittableList(make_shared<bvhNode>(world));


    // Render-------------------------------------------------------------------------------------
    camera cam;

    ...

    cam.render(world);
}

现在渲染超快！

3.10 Another BVH Optimization

我们可以进一步优化BVH。与其随机选择一个轴向进行BVH划分，不如直接选择AABB中最长边界所在的轴向作为划分的依据，从而能够进行更深入的划分。

我们首先为aabb类添加一些内容，包括用于返回AABB中最长边界所在轴向的函数。此外，我们再声明两个大小为空和无限大的AABB：

class aabb
{
public:
    ...

    int longestAxisIndex() const
    {
        // return the index of the longest axis of the AABB
        if (x.size() > y.size())
            return x.size() > z.size() ? 0 : 2;
        else
            return y.size() > z.size() ? 1 : 2;
    }
    
    static const aabb empty, universe;
};

const aabb aabb::empty = aabb(interval::empty, interval::empty, interval::empty);
const aabb aabb::universe = aabb(interval::universe, interval::universe, interval::universe);

当然，我们需要为interval类补充size()这个函数：

double size() const { return max - min; }

现在，我们回到BVH的构造函数。首先，我们先声明一个大小为空的AABB，然后将其边界扩展到可以包含所有对象的范围，最后再使用aabb::lonestAxisIndex()得到AABB中最长边界所在的轴向：

class bvhNode final : public hittable
{
public:
    ...
    
    bvhNode(std::vector<shared_ptr<hittable>>& objects, size_t start, size_t end)
    {
        // build thr bounding box of the span of source objects
        bbox = aabb::empty; // bbox = aabb(); should work as well
        for (size_t objectIndex = start; objectIndex < end; objectIndex++)
        {
            bbox = aabb(bbox, objects[objectIndex]->boundingBox());
        }
        int axis = bbox.longestAxisIndex();

		...

        bbox = aabb(left->boundingBox(), right->boundingBox());
    }

	...
};