BVH层级包围体原理及实现代码

在本章中，我们将了解 Kay 和 Kajiya 使用和描述的BVH方法。 将边界体块分组为我们自己分组的更大体积等的想法非常简单，也很容易理解。

推荐：用 NSDT设计器 快速搭建可编程3D场景。

1、层级包围体块

在图 2 中，我们显示了一组与其各自包围体相关联的对象（为简单起见，这里是一个框）。 通过将这些边界框合并在一起（通常通过邻近度），我们获得了代表对象组的更大边界体块。

图 2：包围体层次结构示例。 边界框 C 中包含的对象不需要进行测试。

很容易看出，如果光线不与这些较大的组中的任何一个相交，那么我们就可以避免测试这些组所包围的体积，可能会同时拒绝许多对象。 显然，这节省了大量渲染时间。

原理非常简单，但为了提高效率，对象和体积必须按邻近程度进行分组。 如果它们不按邻近程度分组，则会出现问题，如图 3 所示。红色茶壶的两个边界框即使彼此相距很远，也会被分组在一起。 在图 3 中，射线与两个包围体相交，并且必须测试四个茶壶是否与射线相交。 在图 1 中，仅测试其中两个茶壶的交叉点（D 和 E）。

图 3：空间被细分为更小的区域。

为了按邻近度对对象进行分组，Kay 和 Kajiya 建议将它们插入空间分区数据结构中。 这种结构将空间划分为子区域，并且对象通常根据其位置插入到这些子区域中。 该过程如图 3 所示。

如果我们创建一个包围场景中所有对象的框，那么我们可以将该框细分为八个大小相等的子框（在图 4 和 5 中，此过程以 2D 形式进行说明）。 然后可以将场景的每个对象插入（按照它们添加到场景中的顺序）它们重叠的子框中。 然后我们可以假设一个子体素中包含的所有对象彼此非常接近（至少比其他子框中的对象更接近）。

图 4：二维的八叉树是四叉树。 该节点分为四个单元，对象被插入到它们重叠的单元中。 细胞细分的过程可以根据需要重复多次。 例如，直到每个单元格中有一个对象或当我们达到用户定义的最大深度时。

如图 4 所示，一个对象可能与多个子框或单元格重叠。 在这种情况下，Kay 和 Kajiya 任意将物体插入其中一个重叠的单元格中。 我们选择将它们插入到对象边界框质心所在的单元格中（由每个茶壶中心的点表示）。

论文提出了两种空间数据结构：中值割和八叉树。 第一个与您可能听说过的 k-d 树空间分区数据结构有关。 我们不会在本课中介绍这些结构。 后者，即八叉树，将用于我们实现 Kay 和 Kajiya 的论文。 如前所述，它将一个立方体划分为八个单元，这些单元本身又被细分等。当所有对象都插入到八叉树中或达到用户定义的最大深度时，递归过程就会停止。

2、建立层次结构

八叉树数据结构在另一课中详细介绍。 不过，简而言之，对象的创建和插入过程如下。 首先，我们计算场景中所有对象的边界体积，然后计算整个场景边界体积，这是这些体积组合的结果。

从场景的整体体积中，我们可以定义八叉树的大小（位于整体体积中心的立方体，其尺寸是沿着 x、y 和 z 轴的任何体积范围的最大值）。 这个立方体代表八叉树的顶部节点，即它的根。

当我们在这个八叉树中插入对象时，我们以自上而下的方式遍历这棵树。 如果我们要插入对象的节点是叶子，并且叶子还没有包含任何对象，那么我们将对象插入到该节点中。

如果节点已经包含一个对象（除非我们已达到用户定义的最大深度），则当前节点将被分割为八个单元格，并且该节点包含的对象以及新对象将被插入到单元格中（注意对象 可能重叠多个单元格）它们重叠。

图 5：插入四叉树中的对象。

重复此过程，直到所有对象都插入到八叉树中。

第二步，这次以自下而上的方式遍历八叉树。 我们从叶子开始，计算它们所包含的对象的整体边界框（叶子可以有一个或多个对象）。 然后，通过组合其子边界体积来计算叶子上方节点的整体体积框。 重复这个过程直到我们到达根（此时应该与我们之前计算的整个场景的范围相同）。

在这两个过程结束时，我们获得了场景的表示形式，作为包围体的层次结构，其中包围体根据邻近度进行分组（节点叶子中包含的对象被分组在一起，等等）。 请注意，赭石实际上并不用作加速结构（不要将此技术与用作加速结构的八叉树混淆 - 我们将在将来编写有关此技术的课程）。

八叉树实际上对于光线相交测试没有真正的好处。 它仅用于计算和构建此卷层次结构。

图 6：与四叉树相交。 如果某个单元格相交，则测试该单元格的子级是否存在相交。 这个过程一直持续到叶节点相交为止。 然后测试该叶节点中包含的几何形状。
BVH 的交集非常简单。 我们从根节点开始，测试射线是否与该节点的包围体相交。

通常，如果节点的包围体与射线相交，我们会在树中向下一层并测试与该节点子节点的包围体的相交。 如果节点是叶子，我们测试射线是否与其包含的任何对象相交。 该解决方案实施起来非常简单，但是可以进一步优化。

当我们测试射线与节点子节点的包围体的相交时，假设多个单元的包围体已相交，我们应该首先研究相交距离最小的节点的子节点，因为可见对象最有可能包含在这些单元中。 这个想法如图 6 所示（我们用光线与单元边界相交处的点来表示与单元边界体积的交点。请记住，实际上，我们与单元边界体积有一个交点，由图中的球体表示，而不是单元本身）。

正如你所看到的，光线与根节点（黑色）相交，这导致测试根的子节点（红色）。 单元按照创建的顺序进行测试：0、1、2、3。射线与单元 0 和 1 的包围体相交，因此接下来测试这些节点的子节点。 然而问题是，单元 0 的子单元将在单元 1 的子单元之前进行测试，即使单元 1 与边界体积的相交距离 (t1) 小于与单元 0 的体积的相交距离 (t0）。 在测试单元 0 中的子项之前先测试单元 1 中的子项会更有效，因为这样可以更快地找到与茶壶 A 的交集。

Kay 和 Kajiya 建议使用一个列表，其中插入被射线穿过的节点的子节点，并根据它们的相交距离（从最小到最大）进行排序。 然后，我们从此列表中删除第一个节点并测试其子节点。 如果这些子级中的任何一个与射线相交，它们就会将自己插入到列表中。 它们将被插入到列表中现有元素的前面，因为它们的相交距离必然小于已插入的任何节点的相交距离。

通过遵循此过程，我们确保始终首先测试层次结构中距离射线原点最近的节点。 如果这些节点不包含可见对象，我们将继续测试距离射线原点较远的节点等，直到列表为空。

while (priority_queue is not empty) { 
    Node node = get node from the priority_queue with highest priority (and remove it) 
    if node is a leaf 
        for each object contained in node 
            if object is intersected by ray 
                 keep this object as potential visible object 
        if an object was intersected 
            return 
    else 
        // keep traversing the hierarchy by looking down the node's children
        for each children in node 
            if node's children is intersected by ray at distance t 
                insert node into priority_queue (use intersection distance t to define priority of the node in the list) 
}

图 7：第一个相交包围体不一定包含最近的相交对象。

我们需要注意最后一个细节。 如果查看图 7，你会发现我们无法使用到边界体积的相交距离来确定这些体积中的哪些体积包含可见对象。 即使到 B 的包围体的距离tvb小于到 C 的包围体的距离tvc，光线将相交的对象是 C 而不是 B。我们只能确定，当到相交对象的距离小于到列表中下一个节点的包围体的相交距离时，我们可以停止层次结构的遍历。

如图7所示，优先级列表中的节点应按以下顺序VB、VC、VA（因为tvb < tvc < tva）。 当我们与 VB 包含的对象（B 的边界体积）相交时，我们找到对象 B的相交距离tb。但是tb不小于tvc，因此我们还必须测试与对象 C 的交集，我们发现tc低于
tb，因此 C 成为潜在的相交对象而不是 B。但是因为tc低于tva，我们不需要测试与 A 的交集，也不需要测试与 VA 之后列表中的任何节点的交集。 因此，我们可以停止层次遍历并返回 B 作为可见对象。

如果我们在光线追踪器中使用所有这些技术，我们现在会得到以下结果：

Render time                                 : 1.61 (sec)
Total number of triangles                   : 16384
Total number of primary rays                : 307200
Total number of ray-triangles tests         : 41341952
Total number of ray-triangles intersections : 59017
Total number of ray-volume tests            : 1531064

渲染速度现在快了 1.8 倍（与前一章的渲染相比）。 射线体积测试的数量减少了 6.4（证明包围体层次结构非常有效），射线三角形测试的数量减少了 1.95，射线三角形相交的数量减少了 1.89。

3、BVH实现代码

完整的源代码可在本课程的最后一章中找到。

BVH 类变得更加复杂。 我们向该类添加了一个 OctreeNode 和一个 Octree 类。 有兴趣了解八度音程的读者可以参考以下部分，他们可以在其中找到有关此主题的课程。 八叉树节点保存八个指向其子节点的其他八叉树节点的指针。 类构造函数将场景范围作为输入参数（第 38 行）。

该范围用于计算和设置根节点的尺寸和位置（其质心。第 40-49 行）。 新对象（更准确地说是它们的边界体积，但请记住 Extents 类保存指向封闭对象的指针）从根节点插入（第 52-53 行）。

如果当前节点是一个空叶子（尚未在该叶子中插入任何对象），则将对象插入该叶子中并返回。

如果该节点是叶子并且包含至少一个或多个对象，但我们已经达到八叉树最大深度，我们仍然将该对象插入到该节点中（第 73 行）。

然而，如果我们还没有达到最大深度，我们将节点标记为“内部”，并重新插入节点已经保存的对象以及八叉树中的新对象（第77-82行）。 如果该节点不是叶子，我们检查该对象应该插入该节点的哪些子节点（第 86-93 行）。 如果子节点尚不存在，我们将创建它（第 97-98 行），最后将对象插入到该子节点中（第 99 行）。

构建卷层次结构的函数很简单。 我们从根开始，沿着树向下走，直到到达叶节点。 该节点的包围体是通过组合它所指向的所有对象的包围体来计算的。 然后，当我们再次向上移动时，当前节点的包围体是通过组合其子节点的包围体来计算的（第 130 行）。

class BVH : public AccelerationStructure
{
    static const uint8_t kNumPlaneSetNormals = 7;
    static const Vec3f planeSetNormals[kNumPlaneSetNormals];
    struct Extents
    {
        Extents()
        {
            for (uint8_t i = 0; i < kNumPlaneSetNormals; ++i)
                d[i][0] = kInfinity, d[i][1] = -kInfinity;
        }
        void extendBy(const Extents &extents)
        {
            for (uint8_t i = 0; i < kNumPlaneSetNormals; ++i) {
                if (extents.d[i][0] < d[i][0]) d[i][0] = extents.d[i][0];
                if (extents.d[i][1] > d[i][1]) d[i][1] = extents.d[i][1];
            }
        }
        bool intersect(
            const float *precomputedNumerator, const float *precomputeDenominator, 
            float &tNear, float &tFar, uint8_t &planeIndex);
        float d[kNumPlaneSetNormals][2]; // d values for each plane-set normals
        const Object *object; // pointer contained by the volume (used by octree)
    };
    Extents *extents;
    struct OctreeNode
    {
        OctreeNode *child[8];
        std::vector<const Extents *> data;
        Extents extents;
        bool isLeaf;
        uint8_t depth; // just for debugging
        OctreeNode() : isLeaf(true) { memset(child, 0x0, sizeof(OctreeNode *) * 8); }
        ~OctreeNode() { for (uint8_t i = 0; i < 8; ++i) if (child[i] != NULL) delete child[i]; }
    };
    struct Octree
    {
        Octree(const Extents &extents) : root(NULL)
        {
            float xdiff = extents.d[0][1] - extents.d[0][0];
            float ydiff = extents.d[1][1] - extents.d[1][0];
            float zdiff = extents.d[2][1] - extents.d[2][0];
            float dim = std::max(xdiff, std::max(ydiff, zdiff));
            Vec3f centroid(
                (extents.d[0][0] + extents.d[0][1]),
                (extents.d[1][0] + extents.d[1][1]),
                (extents.d[2][0] + extents.d[2][1]));
            bounds[0] = (Vec3f(centroid) - Vec3f(dim)) * 0.5f;
            bounds[1] = (Vec3f(centroid) + Vec3f(dim)) * 0.5f;
            root = new OctreeNode;
        }
        void insert(const Extents *extents)
        { insert(root, extents, bounds[0], bounds[1], 0); }
        void build()
        { build(root, bounds[0], bounds[1]); }
        ~Octree() { delete root; }
        struct QueueElement
        {
            const OctreeNode *node; // octree node held by this node in the tree
            float t; // used as key
            QueueElement(const OctreeNode *n, float thit) : node(n), t(thit) {}
            // comparator is > instead of < so priority_queue behaves like a min-heap
            friend bool operator < (const QueueElement &a, const QueueElement &b) { return a.t > b.t; }
        };
        Vec3f bounds[2];
        OctreeNode *root;
    private:
        void insert(
            OctreeNode *node, const Extents *extents, 
            Vec3f boundMin, Vec3f boundMax, int depth)
        {
            if (node->isLeaf) {
                if (node->data.size() == 0 || depth == 16) {
                    node->data.push_back(extents);
                }
                else {
                    node->isLeaf = false;
                    while (node->data.size()) {
                        insert(node, node->data.back(), boundMin, boundMax, depth);
                        node->data.pop_back();
                    }
                    insert(node, extents, boundMin, boundMax, depth);
                }
            } else {
                // insert bounding volume in the right octree cell
                Vec3f extentsCentroid = (
                    Vec3f(extents->d[0][0], extents->d[1][0], extents->d[2][0]) +
                    Vec3f(extents->d[0][1], extents->d[1][1], extents->d[2][1])) * 0.5;
                Vec3f nodeCentroid = (boundMax + boundMin) * 0.5f;
                uint8_t childIndex = 0;
                if (extentsCentroid[0] > nodeCentroid[0]) childIndex += 4;
                if (extentsCentroid[1] > nodeCentroid[1]) childIndex += 2;
                if (extentsCentroid[2] > nodeCentroid[2]) childIndex += 1;
                Vec3f childBoundMin, childBoundMax;
                Vec3f boundCentroid = (boundMin + boundMax) * 0.5;
                computeChildBound(childIndex, boundCentroid, boundMin, boundMax, childBoundMin, childBoundMax);
                if (node->child[childIndex] == NULL)
                      node->child[childIndex] = new OctreeNode, node->child[childIndex]->depth = depth;
                insert(node->child[childIndex], extents, childBoundMin, childBoundMax, depth + 1);
            }
        }
        void computeChildBound(
            const uint8_t &i, const Vec3f &boundCentroid,
            const Vec3f &boundMin, const Vec3f &boundMax,
            Vec3f &pMin, Vec3f &pMax) const
        {
            pMin[0] = (i & 4) ? boundCentroid[0] : boundMin[0];
            pMax[0] = (i & 4) ? boundMax[0] : boundCentroid[0];
            pMin[1] = (i & 2) ? boundCentroid[1] : boundMin[1];
            pMax[1] = (i & 2) ? boundMax[1] : boundCentroid[1];
            pMin[2] = (i & 1) ? boundCentroid[2] : boundMin[2];
            pMax[2] = (i & 1) ? boundMax[2] : boundCentroid[2];    
        }
        // bottom-up construction
        void build(OctreeNode *node, const Vec3f &boundMin, const Vec3f &boundMax)
        {
            if (node->isLeaf) {
                // compute leaf node bounding volume
                for (uint32_t i = 0; i < node->data.size(); ++i) {
                    node->extents.extendBy(*node->data[i]);
                }
            }
            else {
                for (uint8_t i = 0; i < 8; ++i)
                    if (node->child[i]) {
                        Vec3f childBoundMin, childBoundMax;
                        Vec3f boundCentroid = (boundMin + boundMax) * 0.5;
                        computeChildBound(i, boundCentroid, boundMin, boundMax, childBoundMin, childBoundMax);
                        build(node->child[i], childBoundMin, childBoundMax);
                        node->extents.extendBy(node->child[i]->extents);
                    }
                }
            }
        }
    };
    Octree *octree;
public:
    BVH(const RenderContext *rcx);
    const Object* intersect(const Ray<float> &ray, IsectData &isectData) const;
    ~BVH();
};

在 BVH 类的 intersect 方法中（第 37 行），我们首先测试光线是否与整个场景的包围体（赭石根节点的范围）相交。

如果确实如此，我们用该节点以及从射线原点到其包围体的相交距离初始化一个priority_queue列表。 我们重载了运算符 <（上面第 63 行），使其表现得像最小堆（默认情况下，它的表现像最大堆，将具有最高键值的元素设置为最高优先级。我们想要相反）。

其余代码与上面给出的伪代码类似。 我们获取列表顶部的第一个节点（我们也将其从列表中删除。第 72-73 行）。 如果该节点是叶子，则我们测试光线是否与该节点包含的任何对象相交，并跟踪相交最小距离（第 79 行）。 如果它是内部节点，我们测试射线是否与该节点的子节点的包围体相交，当相交时，我们将子节点添加到列表中（使用相交距离作为设置元素在列表中的优先级的键。第93行）。

重复此过程直到列表为空，但如果列表第一个节点的相交距离大于到相交对象的距离（第 71 行），则可以更快停止。

BVH::BVH(const RenderContext *rcx) : AccelerationStructure(rcx), extents(NULL), octree(NULL) 
{ 
    Extents sceneExtents; 
    extents = new Extents[rcx->objects.size()]; 
    for (uint32_t i = 0; i < rcx->objects.size(); ++i) { 
        for (uint8_t j = 0; j < kNumPlaneSetNormals; ++j) { 
            rcx->objects[i]->computeBounds(planeSetNormals[j], extents[i].d[j][0], extents[i].d[j][1]); 
        } 
        extents[i].object = rcx->objects[i]; 
        sceneExtents.extendBy(extents[i]); 
    } 
    // create hierarchy                                                                                                                                                                                     
    octree = new Octree(sceneExtents); 
    for (uint32_t i = 0; i < rcx->objects.size(); ++i) { 
        octree->insert(extents + i); 
    } 
    octree->build(); 
} 
 
inline bool BVH::Extents::intersect( 
    const float *precomputedNumerator, const float *precomputeDenominator, 
    float &tNear, float &tFar, uint8_t &planeIndex) 
{ 
    __sync_fetch_and_add(&numRayVolumeTests, 1); 
    for (uint8_t i = 0; i < kNumPlaneSetNormals; ++i) { 
        float tn = (d[i][0] - precomputedNumerator[i]) / precomputeDenominator[i]; 
        float tf = (d[i][1] - precomputedNumerator[i]) / precomputeDenominator[i]; 
        if (precomputeDenominator[i] < 0) std::swap(tn, tf); 
        if (tn > tNear) tNear = tn, planeIndex = i; 
        if (tf < tFar) tFar = tf; 
        if (tNear > tFar) return false;  //test for an early stop 
    } 

    return true; 
} 
 
const Object* BVH::intersect(const Ray<float> &ray, IsectData &isectData) const 
{ 
    const Object *hitObject = NULL; 
    float precomputedNumerator[BVH::kNumPlaneSetNormals], precomputeDenominator[BVH::kNumPlaneSetNormals]; 
    for (uint8_t i = 0; i < kNumPlaneSetNormals; ++i) { 
        precomputedNumerator[i] = dot(planeSetNormals[i], ray.orig); 
        precomputeDenominator[i] = dot(planeSetNormals[i], ray.dir);; 
    } 
#if 0 
    float tClosest = ray.tmax; 
    for (uint32_t i = 0; i < rc->objects.size(); ++i) { 
        __sync_fetch_and_add(&numRayVolumeTests, 1); 
        float tNear = -kInfinity, tFar = kInfinity; 
        uint8_t planeIndex; 
        if (extents[i].intersect(precomputedNumerator, precomputeDenominator, tNear, tFar, planeIndex)) { 
            IsectData isectDataCurrent; 
            if (rc->objects[i]->intersect(ray, isectDataCurrent)) { 
                if (isectDataCurrent.t < tClosest && isectDataCurrent.t > ray.tmin) { 
                    isectData = isectDataCurrent; 
                    hitObject = rc->objects[i]; 
                    tClosest = isectDataCurrent.t; 
                } 
            } 
        } 
    } 
#else 
    uint8_t planeIndex = 0; 
    float tNear = 0, tFar = ray.tmax; 
    if (!octree->root->extents.intersect(precomputedNumerator, precomputeDenominator, tNear, tFar, planeIndex) 
        || tFar < 0 || tNear > ray.tmax) 
        return NULL; 
    float tMin = tFar; 
    std::priority_queue<BVH::Octree::QueueElement> queue; 
    queue.push(BVH::Octree::QueueElement(octree->root, 0)); 
    while(!queue.empty() && queue.top().t < tMin) { 
        const OctreeNode *node = queue.top().node; 
        queue.pop(); 
        if (node->isLeaf) { 
            for (uint32_t i = 0; i < node->data.size(); ++i) { 
                IsectData isectDataCurrent; 
                if (node->data[i]->object->intersect(ray, isectDataCurrent)) { 
                    if (isectDataCurrent.t < tMin) { 
                        tMin = isectDataCurrent.t; 
                        hitObject = node->data[i]->object; 
                        isectData = isectDataCurrent; 
                    } 
                } 
            } 
        } 
        else { 
            for (uint8_t i = 0; i < 8; ++i) { 
                if (node->child[i] != NULL) { 
                    float tNearChild = 0, tFarChild = tFar; 
                    if (node->child[i]->extents.intersect(precomputedNumerator, precomputeDenominator, 
                        tNearChild, tFarChild, planeIndex)) { 
                        float t = (tNearChild < 0 && tFarChild >= 0) ? tFarChild : tNearChild; 
                        queue.push(BVH::Octree::QueueElement(node->child[i], t)); 
                    } 
                } 
            } 
        } 
    } 
#endif 

    return hitObject; 
}

通常，我们应该在 BVH 中插入的对象是三角形而不是网格。 为了尽快完成基本部分，我们将在本课程的后续版本中解释如何完成此操作。 但在下一章中，我们将展示如何将三角形插入网格加速结构中。 作为练习，你可以尝试修改本课程中找到的 BVH 代码，以支持三角形的插入，以网格实现为例。

---

原文链接：BVH原理及实现代码 — BimAnt