光线求交加速算法：边界体积层次结构(Bounding Volume Hierarchies)1-BVH引入

光线求交加速算法：边界体积层次结构(Bounding Volume Hierarchies)1

 

BVH引入

光线和物体求交的加速算法中，最常见的是物体(图元)细分和空间细分。边界体积层次结构（BVH）是一种基于图元(primitives)细分的光线相交加速方法，其中，图元被划分为不相交集的层次结构。 （相反，空间细分通常将空间划分为不相交集的层次结构。）下图显示了一个简单场景的边界框层次结构：

其中，图元存储在叶中，每个节点存储其所有子节点图元的边界框。 因此，当光线穿过树时，只要它不与节点的边界相交，就可以跳过该节点下方的子树。

对于一个二叉树形式的BVH，其在叶子结点中存储一个单一的图元，该树节点总数为2n-1，其中n是图元数。将有n个叶节点和n-1个内部节点。如果叶子存储多个图元，则需要更少的节点。BVH比kd树更快速地构建，kd树通常比BVH提供更快的射线相交测试，但是建造时间要长得多。另一方面，与kd树相比，BVH通常在数值上更健壮，并且由于舍入误差而导致的交集遗漏更少。

 

构建BVH

BVH的构建分为三个阶段。

1. 计算有关每个图元的边界信息并将其存储在数组中，该数组将在树构建期间使用。
2. 使用相应的分割算法来构建树，本文主要介绍四种分区算法(SAH, HLBVH, Middle, EqualCounts) ，构建的结果是一棵二叉树，其中每个内部节点持有指向其子节点的指针，每个叶子节点持有对一个或多个基元的引用。
3. 最后该树被转换为更紧凑(因此效率更高)的无指针表示，以在渲染期间使用。

我们根据这三个阶段写出BVH构建代码：

// BVH构造函数，构造BVH2叉树
BVHAccel::BVHAccel(std::vector<std::shared_ptr<Primitive>> p,
                   int maxPrimsInNode, SplitMethod splitMethod)
    : maxPrimsInNode(std::min(255, maxPrimsInNode)),
      splitMethod(splitMethod),
      primitives(std::move(p)) {

    // 利用图元构造BVH2叉树，从以下三个阶段进行。(+)表示代码未展开 (=)表示展开后
    // 1.+初始化图元信息数组 
    // 2.+利用图元信息并使用相应的图元分割算法构建BVH
    // 3.+树的节点用深度优先排列表示
}

对于要存储在BVH中的每个图元，我们将其边界框的质心，其完整边界框及其索引存储在BVHPrimitiveInfo结构体的图元数组中。该结构体如下：

struct BVHPrimitiveInfo {
    BVHPrimitiveInfo(size_t primitiveNumber, const Bounds3f &bounds)
        : primitiveNumber(primitiveNumber), bounds(bounds),
          centroid(.5f * bounds.pMin + .5f * bounds.pMax) { }
    size_t primitiveNumber;//图元索引
    Bounds3f bounds;//图元边界框
    Point3f centroid;//图元边界框的质心
};

我们初始化图元信息数组的代码则为：

//1.=初始化图元信息数组
std::vector<BVHPrimitiveInfo> primitiveInfo(primitives.size());
for (size_t i = 0; i < primitives.size(); ++i)
    primitiveInfo[i] = { i, primitives[i]->WorldBound() };

接下来使用相应的图元分区算法以及图元信息构建BVH树，我们贴出代码：

// 2.=使用图元信息与相应的图元分割算法构建BVH树
    MemoryArena arena(1024 * 1024);
    int totalNodes = 0;
    std::vector<std::shared_ptr<Primitive>> orderedPrims;
    orderedPrims.reserve(primitives.size());
    BVHBuildNode *root;//构建完成树的根节点
    if (splitMethod == SplitMethod::HLBVH)
        root = HLBVHBuild(arena, primitiveInfo, &totalNodes, orderedPrims);
    else
        root = recursiveBuild(arena, primitiveInfo, 0, primitives.size(),
                              &totalNodes, orderedPrims);
    primitives.swap(orderedPrims);//orderedPrims是构造后排序过的图元，我们用其覆盖原图元

其中，我们把图元分区算法用枚举类型表示：

//SAH:表面积启发式法。HLBVH：线性边界体积层次结构法。Middle：中点分割法。EqualCounts：等数量分割法
enum class SplitMethod { SAH, HLBVH, Middle, EqualCounts };

HLBVH方法我们利用函数HLBVHBuild()构建，其他三个方法在函数recursiveBuild()中区分构建。MemoryArena为PBRT中内存管理类，用以优化内存分配。

每个BVHBuildNode代表BVH的一个节点。 所有节点都存储一个Bounds3f，它表示该节点下所有子级的边界框。 每个内部节点在其子级中存储指向其两个子级的指针。 内部节点还记录了坐标轴，图元沿该坐标轴进行了划分，以分配给它们的两个子节点。叶子节点需要记录其中存储了哪些图元，索引从firstPrimOffset到firstPrimOffset+nPrimitives之间的图元都为该叶子结点内的图元。内部结点则不保存图元，仅保存边界框信息。BVHBuildNode结构如下：

struct BVHBuildNode {
    // 构造叶子结点
    void InitLeaf(int first, int n, const Bounds3f &b) {
        firstPrimOffset = first;
        nPrimitives = n;
        bounds = b;
        children[0] = children[1] = nullptr;
    }
    //构造内部结点
    void InitInterior(int axis, BVHBuildNode *c0, BVHBuildNode *c1) {
        children[0] = c0;
        children[1] = c1;
        bounds = Union(c0->bounds, c1->bounds);//Union表示联合c1和c2的边界框
        splitAxis = axis;
        nPrimitives = 0;
    }
    Bounds3f bounds;
    BVHBuildNode *children[2];
    int splitAxis, firstPrimOffset, nPrimitives;
};

除了用于分配节点内存的MemoryArena和保存图元信息的数组BVHPrimitiveInfo之外，recursiveBuild()还将范围[start，end）作为参数。它负责返回由从primaryInfo [start]到包括primitiveInfo [end-1]的范围表示的图元子集的BVH。如果此范围仅涵盖单个图元，则递归已触底并创建叶节点。否则，将使用一种分区算法在该范围内对数组的元素进行分区，并相应地对该范围内的数组元素进行重新排序，以使[start，mid）和[mid，end）的范围代表已分区的子集。如果分区成功，则将这两个图元集依次传递给递归调用，这些递归调用本身将返回指向当前节点的两个子节点的指针。totalNodes跟踪已创建的BVH节点的总数。orderedPrims数组用于存储已排序的图元引用。代码如下：

BVHBuildNode *BVHAccel::recursiveBuild(MemoryArena &arena,
        std::vector<BVHPrimitiveInfo> &primitiveInfo, int start,
        int end, int *totalNodes,
        std::vector<std::shared_ptr<Primitive>> &orderedPrims) {
    BVHBuildNode *node = arena.Alloc<BVHBuildNode>();
    (*totalNodes)++;
    // +计算所有BVH图元的边界框
    int nPrimitives = end - start;
    if (nPrimitives == 1) {
        // +创建叶子结点
    } else {
        // +计算图元质心的边界，选择分割的坐标轴
        // +将图元集分割成两个子集
    }
    return node;
}

首先展开计算边界框的代码：

// =计算所有SVH图元的边界框
Bounds3f bounds;
for (int i = start; i < end; ++i)
    bounds = Union(bounds, primitiveInfo[i].bounds);

在叶节点，将与叶重叠的图元附加到orderedPrims数组，并初始化叶节点对象：

// =创建叶子结点
int firstPrimOffset = orderedPrims.size();
for (int i = start; i < end; ++i) {
    int primNum = primitiveInfo[i].primitiveNumber;
    orderedPrims.push_back(primitives[primNum]);//由于前序遍历，所以保存是有序的
}
node->InitLeaf(firstPrimOffset, nPrimitives, bounds);
return node;

对于内部结点我们按什么方式为图元分区呢？我们通常选择一个坐标轴，图元覆盖该坐标轴具有最大的范围，如下图所示：

此处进行分区的总体目标是选择一个基本图元分区，该图元与两个所得基本图元集的边界框没有太多的重叠-如果存在实质性的重叠，遍历两个子树可能比更有效地分割图元的集合需要更多的计算。在讨论表面积启发式方法(SAH)时，将很快更加严格地找到有效的原始分区。

我们这里先写出寻找最大范围的坐标轴：

// =计算图元质心的边界，选择分割的坐标轴
Bounds3f centroidBounds;
for (int i = start; i < end; ++i)
    centroidBounds = Union(centroidBounds, primitiveInfo[i].centroid);
int dim = centroidBounds.MaximumExtent();

如果所有质心点都在同一位置（即质心边界的体积为零），则递归停止，并使用图元创建叶节点； 在这种(异常)情况下，此处的分割方法均无效。 否则，将使用所选方法对图元进行分区，并将其传递给recursiveBuild()，对其进行两个递归调用。代码如下：

// =将图元集分割成两个子集
int mid = (start + end) / 2;
if (centroidBounds.pMax[dim] == centroidBounds.pMin[dim]) {
    // +创建叶子结点
} else {
    // +根据不同的分割方法，分割图元集
    node->InitInterior(dim,
                       recursiveBuild(arena, primitiveInfo, start, mid,
                                      totalNodes, orderedPrims),
                       recursiveBuild(arena, primitiveInfo, mid, end,
                                      totalNodes, orderedPrims));
}

这里创建叶子结点代码与上部分创建叶子结点完全相同，在递归调用之前我们需要完成图元集分区，一个简单的splitMethod是Middle，它首先计算沿分割轴的图元质心的中点。根据它们的质心是在中点之上还是之下，将图元分为两组。使用C ++标准库函数std :: partition()可以轻松完成此分区，该函数接受数组中的一系列元素和比较函数，并对数组中的元素进行排序。std :: partition()返回一个指针，该指针指向该谓词具有假值的第一个元素，该元素将转换为对primitiveInfo数组的偏移量，以便我们可以将其传递给递归调用:

switch (splitMethod) {
    case SplitMethod::Middle: {
        // 中点法
        Float pmid = (centroidBounds.pMin[dim] + centroidBounds.pMax[dim]) / 2;
        BVHPrimitiveInfo *midPtr = std::partition(
            &primitiveInfo[start], &primitiveInfo[end - 1] + 1,
            [dim, pmid](const BVHPrimitiveInfo &pi) {
                return pi.centroid[dim] < pmid;
            });
        mid = midPtr - &primitiveInfo[0];
        // 对于许多带有大重叠边界框的素数，可能无法分区； 在这种情况下将不执行break跳出，而是进行下一种分区法(等数量分区法)
        if (mid != start && mid != end) break;
        }
        //...其他分区方法
}

如果所有图元都具有大的重叠边界框，则此拆分方法可能无法将图元分为两组。在这种情况下，执行将落入SplitMethod :: EqualCounts方法来重试。EqualCounts方法将图元划分为两个相等大小的子集，以使它们中的n个的前半部分是沿选定轴的质心坐标值最小的n / 2个，而后半部分是质心坐标值最大的半个子集。尽管这种方法有时效果很好，下图(b)的情况则这种方法的效果较差：

如图，a)对于某些图元分布，例如此处所示的分布，基于质心沿所选轴的中点(粗蓝线)进行分割的效果很好。(用虚线显示两个结果原始组的边界框。)(b)对于这种分布，中点是次优的选择；两个结果边界框基本重叠。(c)如果将(b)中的同一组图元沿着此处所示的线进行拆分，则生成的边界框会更小并且完全不重叠，从而在渲染时会带来更好的性能。

我们列出EqualCounts方法的代码：

switch (splitMethod) {
    //中点法
    case SplitMethod::Middle:{
        //...
    }
    //等数量法
    case SplitMethod::EqualCounts: {
        // 利用一半数量分区
        mid = (start + end) / 2;
        std::nth_element(&primitiveInfo[start], &primitiveInfo[mid],
                         &primitiveInfo[end - 1] + 1,
                         [dim](const BVHPrimitiveInfo &a,
                               const BVHPrimitiveInfo &b) {
                             return a.centroid[dim] < b.centroid[dim];
                         });
        break;
    }
}

通过标准库调用std :: nth_element()也可以轻松实现此方案。 它需要一个开始，中间和结束指针以及一个比较函数。 它对数组进行排序，以使中间指针处的元素是数组完全排序后将存在的元素，并且使得中间一个元素之前的所有元素小于中间元素，而后面一个元素的所有元素都小于中间元素 它比它更大。 这种排序可以在O(n)时间内完成，元素数量为n，这比对数组进行完全排序的O(nlogn)更有效。

EqualCounts和Middle是比较简单的方法，但同时会有比较大的不稳定性，比如一些特殊情况会使分区之后的遍历性能大大降低，如上述提到的图片示例。后面用用两篇介绍两种更好的方法，这里先不提及它们。

 

紧凑的BVH遍历

构建BVH树后，可以将其转换为紧凑的表示形式-这样做可以提高缓存，内存以及整体系统性能。 最终的BVH存储在内存中的线性数组中。 原始树的节点以先根深度优先布局，这意味着每个内部节点的第一个子节点在内存中紧随其后。 在这种情况下，只需显式存储每个内部节点的第二个子节点的偏移量。 相关的树形拓扑和内存中节点顺序排列的关系如下图：

利用LinearBVHNode结构存储遍历BVH所需的信息。除了每个节点的边界框外，对于叶节点，它还存储节点中图元的偏移量和图元计数。对于内部节点，它存储到第二个子节点的偏移量，以及在构建层次结构时图元沿哪个坐标轴进行分区。此信息在下面的遍历中使用，以尝试沿射线从前到后的顺序访问节点。LinearBVHNode结构如下：

struct LinearBVHNode {
    Bounds3f bounds;
    union {
        int primitivesOffset;    // 叶结点
        int secondChildOffset;   // 内部结点
    };
    uint16_t nPrimitives;  // 当是内部结点时，值为0
    uint8_t axis;          // 内部结点的分割轴(x,y,z其中之一)
    uint8_t pad[1];        // 确保整个结构大小为32整数倍
};

对该结构进行填充以确保其大小为32个字节。 这样做可以确保：如果分配的节点使第一个节点与高速缓存行对齐，则随后的节点都不会跨越高速缓存行。

我们通过flattenBVHTree()方法将生成的树转换为LinearBVHNode表示形式，该方法执行深度优先遍历并将节点以线性顺序存储在内存中：

// 3.=树的节点用深度优先排列表示
nodes = AllocAligned<LinearBVHNode>(totalNodes);
int offset = 0;
flattenBVHTree(root, &offset);

BVHAccel类再添加私有成员变量：指向LinearBVHNodes数组的指针：

//BVHAccel类 私有成员变量private：
LinearBVHNode *nodes = nullptr;

将树展平为线性表示很简单； * offset参数将当前偏移量跟踪到BVHAccel的私有成员nodes数组中：

int BVHAccel::flattenBVHTree(BVHBuildNode *node, int *offset) {
    LinearBVHNode *linearNode = &nodes[*offset];
    linearNode->bounds = node->bounds;
    int myOffset = (*offset)++;//每次调用flattenBVHTree执行偏移
    if (node->nPrimitives > 0) {//叶子结点
        linearNode->primitivesOffset = node->firstPrimOffset;
        linearNode->nPrimitives = node->nPrimitives;
    } else {//当内部结点时，创建内部flattened BVH结点
        linearNode->axis = node->splitAxis;
        linearNode->nPrimitives = 0;
        flattenBVHTree(node->children[0], offset);
        linearNode->secondChildOffset = flattenBVHTree(node->children[1], offset);
    }
    return myOffset;
}

在内部节点上，进行递归调用以展平两个子树。 第一个节点直接进行递归调用，第二个节点的偏移量需要通过递归返回。

BVH的遍历会在光线求交时进行，求交代码如下：

bool BVHAccel::Intersect(const Ray &ray,
        SurfaceInteraction *isect) const {
    bool hit = false;
    Vector3f invDir(1 / ray.d.x, 1 / ray.d.y, 1 / ray.d.z);
    int dirIsNeg[3] = { invDir.x < 0, invDir.y < 0, invDir.z < 0 };
    // +跟随ray穿过BVH节点以找到相交的图元
    return hit;
}

每当Intersect()中的while循环开始迭代时，currentNodeIndex会将偏移量保存到要访问的节点的节点数组中。它以值0开头，代表树的根。仍然需要访问的节点存储在nodesToVisit[]数组中，该数组充当栈。toVisitOffset保存到栈中下一个空闲元素的偏移量(toVisitOffset在访问内部结点时总++，相当于入栈，当访问叶子结点时则--，相当于出栈):

// =跟随ray穿过BVH节点以找到相交的图元 
int toVisitOffset = 0, currentNodeIndex = 0;
int nodesToVisit[64];
while (true) {
    const LinearBVHNode *node = &nodes[currentNodeIndex];
    // +对BVH节点检查射线
}

在每个节点处，我们检查射线是否与节点的边界框相交。如果是这样，我们将访问该节点，如果它是叶节点，则测试其与图元的相交，如果是内部节点，则对其子级进行相同处理。 如果未找到相交，则从nodesToVisit []中检索下一个要访问的节点的偏移量(如果堆栈为空，则遍历完成):

// =对BVH节点检查射线
if (node->bounds.IntersectP(ray, invDir, dirIsNeg)) {
    if (node->nPrimitives > 0) {
        // +光线和BVH叶结点中的图元相交 
    } else {
        // +将远BVH节点(第二个子结点)放在nodesToVisit堆栈上，前进到下一个近节点
    }
} else {
    if (toVisitOffset == 0) break;
    currentNodeIndex = nodesToVisit[--toVisitOffset];
}

如果当前节点是叶子，则必须测试射线与它内部的图元的相交。然后从nodesToVisit堆栈中找到下一个要访问的节点。即使在当前节点中找到一个相交，也要访问其余的节点，以防其中一个产生更近的相交。但是，如果找到相交点，则射线的tMax值将更新为相交距离。这使得可以有效地丢弃比交叉点更远的所有剩余节点：

// =光线和BVH叶结点中的图元相交 
for (int i = 0; i < node->nPrimitives; ++i)
    if (primitives[node->primitivesOffset + i]->Intersect(ray, isect))
        hit = true;
if (toVisitOffset == 0) break;
currentNodeIndex = nodesToVisit[--toVisitOffset];

对于射线撞击的内部节点，需要访问其两个子节点。最好是在访问第二个孩子之前先访问射线穿过的第一个孩子(近的结点)，以防光线在第一个孩子中相交，从而可以更新射线的tMax值，从而减少射线的测试次数。

一种有效方法，是使用光线方向矢量的符号作为坐标轴，在该坐标轴上为当前图元划分图元节点：如果符号为负，则应该在第一个孩子之前访问第二个孩子，因为进入第二个孩子的子树的图元在分区点的近侧。(相反为正号方向)。这样做很简单：首先将要访问的节点的偏移量复制到currentNodeIndex，将另一个节点的偏移量添加到nodesToVisit堆栈中。由于内存中节点的深度优先布局，第一个孩子紧随当前节点之后。

// =将远BVH节点(第二个子结点)放在nodesToVisit堆栈上，前进到下一个近节点
if (dirIsNeg[node->axis]) {
   nodesToVisit[toVisitOffset++] = currentNodeIndex + 1;
   currentNodeIndex = node->secondChildOffset;
} else {
   nodesToVisit[toVisitOffset++] = node->secondChildOffset;
   currentNodeIndex = currentNodeIndex + 1;
}

遍历过程图解如下图(红点为图元所在位置)：

值得注意的是：即使光线已经找到和图元的交点，依旧会继续遍历其他边界框和光线有交点的结点，保证交点为最近值。

最后我们给出完整的BVHaccel类：

class BVHAccel : public Aggregate {
  public:
    enum class SplitMethod { SAH, HLBVH, Middle, EqualCounts };

    BVHAccel(std::vector<std::shared_ptr<Primitive>> p,
             int maxPrimsInNode = 1,
             SplitMethod splitMethod = SplitMethod::SAH);
    Bounds3f WorldBound() const;
    ~BVHAccel();
    bool Intersect(const Ray &ray, SurfaceInteraction *isect) const;
    bool IntersectP(const Ray &ray) const;

  private:
    BVHBuildNode *recursiveBuild(
        MemoryArena &arena, std::vector<BVHPrimitiveInfo> &primitiveInfo,
        int start, int end, int *totalNodes,
        std::vector<std::shared_ptr<Primitive>> &orderedPrims);
    BVHBuildNode *HLBVHBuild(
        MemoryArena &arena, const std::vector<BVHPrimitiveInfo> &primitiveInfo,
        int *totalNodes,
        std::vector<std::shared_ptr<Primitive>> &orderedPrims) const;
    BVHBuildNode *emitLBVH(
        BVHBuildNode *&buildNodes,
        const std::vector<BVHPrimitiveInfo> &primitiveInfo,
        MortonPrimitive *mortonPrims, int nPrimitives, int *totalNodes,
        std::vector<std::shared_ptr<Primitive>> &orderedPrims,
        std::atomic<int> *orderedPrimsOffset, int bitIndex) const;
    BVHBuildNode *buildUpperSAH(MemoryArena &arena,
                                std::vector<BVHBuildNode *> &treeletRoots,
                                int start, int end, int *totalNodes) const;
    int flattenBVHTree(BVHBuildNode *node, int *offset);

    const int maxPrimsInNode;
    const SplitMethod splitMethod;
    std::vector<std::shared_ptr<Primitive>> primitives;
    LinearBVHNode *nodes = nullptr;
};