从理论基础到代码实践，全面掌握包围体层次结构

引言

在计算机图形学与游戏开发中，我们经常面临一个核心挑战：如何在海量三维物体中快速找到与光线相交的物体，或者检测物体间的碰撞？最朴素的做法是遍历场景中的所有物体，时间复杂度为O(n)。但当场景包含数万甚至数百万个三角形时，这种做法显然不可行。

**BVH（Bounding Volume Hierarchy，包围体层次结构）**正是解决这一问题的关键数据结构。它能够将光线求交的时间复杂度从O(n)降低到O(log n)，让实时光线追踪和复杂场景的碰撞检测成为可能。本文将带你从零开始，在C++中实现一个完整的BVH系统。

第一章：BVH核心概念解析

1.1 什么是BVH？

BVH是一种基于图元细分（Primitive Subdivision）的树形数据结构。它的核心思想很简单：用简单的包围体（如长方体、球体）包裹复杂的几何物体，然后以层次化的方式组织这些包围体。

如图1所示，BVH的每个叶子节点包含实际的几何图元（如三角形），每个内部节点则包含一个能够包围其所有子节点的包围体。根节点包围整个场景。

1.2 BVH的工作原理

BVH加速查询的基本原理是自上而下的递归遍历：

1. 从根节点开始，检查光线（或查询体）是否与当前节点的包围体相交
2. 如果不相交，则跳过该节点的整个子树——这意味着所有子节点都不需要检查
3. 如果相交且是叶子节点，则对节点内的实际图元进行精确求交
4. 如果相交且是内部节点，则递归检查其子节点

这种机制的精妙之处在于：一次包围体求交失败，就能排除一大片区域内的所有物体，从而大幅减少不必要的精确计算。

1.3 为什么选择AABB？

常见的包围体类型包括：

• AABB（轴对齐包围盒）：与坐标轴对齐的长方体
• OBB（有向包围盒）：可旋转的长方体
• 包围球：球体
• k-DOP：离散定向多面体

在实际应用中，AABB是最常见的选择，原因有三：

1. 计算简单：只需比较坐标值，无需复杂数学运算
2. 存储高效：仅需6个浮点数（最小点和最大点）
3. 紧密度适中：比包围球更紧密，比OBB更容易计算

第二章：BVH的C++实现

接下来，我们将一步步实现一个完整的BVH系统。本节代码参考了多所大学的图形学课程资料和开源库的设计。

2.1 基础数据结构

首先，我们需要定义核心的数据结构：向量、光线、AABB和BVH节点。

// vector3.h - 简单的三维向量类
struct Vector3 {
    float x, y, z;
    
    Vector3() : x(0), y(0), z(0) {}
    Vector3(float x_, float y_, float z_) : x(x_), y(y_), z(z_) {}
    
    Vector3 operator+(const Vector3& other) const {
        return Vector3(x + other.x, y + other.y, z + other.z);
    }
    
    Vector3 operator-(const Vector3& other) const {
        return Vector3(x - other.x, y - other.y, z - other.z);
    }
    
    Vector3 operator*(float t) const {
        return Vector3(x * t, y * t, z * t);
    }
    
    float dot(const Vector3& other) const {
        return x * other.x + y * other.y + z * other.z;
    }
    
    Vector3 normalize() const {
        float len = sqrt(x*x + y*y + z*z);
        return Vector3(x/len, y/len, z/len);
    }
};

// ray.h - 光线类
struct Ray {
    Vector3 origin;      // 光线起点
    Vector3 direction;   // 光线方向
    float tMin, tMax;    // 有效参数范围
    
    Ray(const Vector3& o, const Vector3& d, float tMin_ = 0.001f, float tMax_ = 1e30f)
        : origin(o), direction(d.normalize()), tMin(tMin_), tMax(tMax_) {}
    
    Vector3 pointAt(float t) const {
        return origin + direction * t;
    }
};

// aabb.h - 轴对齐包围盒
struct AABB {
    Vector3 min;  // 最小点
    Vector3 max;  // 最大点
    
    AABB() : min(1e30f, 1e30f, 1e30f), max(-1e30f, -1e30f, -1e30f) {}
    
    AABB(const Vector3& min_, const Vector3& max_) : min(min_), max(max_) {}
    
    // 扩展包围盒以包含另一个点
    void expand(const Vector3& p) {
        min.x = std::min(min.x, p.x);
        min.y = std::min(min.y, p.y);
        min.z = std::min(min.z, p.z);
        max.x = std::max(max.x, p.x);
        max.y = std::max(max.y, p.y);
        max.z = std::max(max.z, p.z);
    }
    
    // 扩展包围盒以包含另一个包围盒
    void expand(const AABB& other) {
        expand(other.min);
        expand(other.max);
    }
    
    // 获取包围盒的中心点
    Vector3 centroid() const {
        return Vector3(
            (min.x + max.x) * 0.5f,
            (min.y + max.y) * 0.5f,
            (min.z + max.z) * 0.5f
        );
    }
    
    // 获取最长轴
    int maxExtent() const {
        Vector3 extent = max - min;
        if (extent.x >= extent.y && extent.x >= extent.z) return 0;
        if (extent.y >= extent.x && extent.y >= extent.z) return 1;
        return 2;
    }
    
    // 光线-包围盒求交
    bool intersect(const Ray& ray, float& t0, float& t1) const {
        t0 = ray.tMin;
        t1 = ray.tMax;
        
        // 对三个轴依次进行 slab 测试
        for (int i = 0; i < 3; i++) {
            float invDir = 1.0f / (&ray.direction.x)[i];
            float tNear = ((&min.x)[i] - (&ray.origin.x)[i]) * invDir;
            float tFar  = ((&max.x)[i] - (&ray.origin.x)[i]) * invDir;
            
            if (tNear > tFar) std::swap(tNear, tFar);
            
            t0 = tNear > t0 ? tNear : t0;
            t1 = tFar  < t1 ? tFar  : t1;
            
            if (t0 > t1) return false;
        }
        return true;
    }
};

2.2 BVH节点定义

BVH节点可以采用两种存储方式：指针方式和数组方式。这里我们先实现指针方式的节点，后续再优化为紧凑数组存储。

// bvh_node.h - BVH节点
struct BVHNode {
    AABB bounds;           // 节点的包围盒
    BVHNode* left;         // 左子节点
    BVHNode* right;        // 右子节点
    
    // 叶子节点数据
    std::vector<Primitive*> primitives;
    bool isLeaf;
    
    BVHNode() : left(nullptr), right(nullptr), isLeaf(false) {}
    
    ~BVHNode() {
        delete left;
        delete right;
    }
    
    // 创建内部节点
    void initInterior(BVHNode* l, BVHNode* r) {
        left = l;
        right = r;
        isLeaf = false;
        bounds = AABB();
        bounds.expand(l->bounds);
        bounds.expand(r->bounds);
    }
    
    // 创建叶子节点
    void initLeaf(const std::vector<Primitive*>& prims) {
        primitives = prims;
        isLeaf = true;
        
        bounds = AABB();
        for (auto prim : primitives) {
            bounds.expand(prim->getBounds());
        }
    }
};

2.3 Primitive基类

为了让BVH能够处理各种几何体，我们需要定义一个统一的基类接口：

// primitive.h - 几何图元基类
class Primitive {
public:
    virtual ~Primitive() {}
    
    // 返回图元的包围盒
    virtual AABB getBounds() const = 0;
    
    // 光线求交
    virtual bool intersect(const Ray& ray, float& t, Vector3& normal) const = 0;
    
    // 获取图元中心（用于BVH构建）
    Vector3 getCentroid() const {
        return getBounds().centroid();
    }
};

// triangle.h - 三角形图元示例
class Triangle : public Primitive {
public:
    Vector3 v0, v1, v2;  // 三个顶点
    Vector3 normal;       // 法线
    
    Triangle(const Vector3& v0_, const Vector3& v1_, const Vector3& v2_)
        : v0(v0_), v1(v1_), v2(v2_) {
        // 计算法线
        Vector3 e1 = v1 - v0;
        Vector3 e2 = v2 - v0;
        normal = e1.cross(e2).normalize();
    }
    
    AABB getBounds() const override {
        AABB bounds;
        bounds.expand(v0);
        bounds.expand(v1);
        bounds.expand(v2);
        return bounds;
    }
    
    bool intersect(const Ray& ray, float& t, Vector3& hitNormal) const override {
        // Möller–Trumbore 三角形求交算法
        const float EPSILON = 1e-6f;
        
        Vector3 e1 = v1 - v0;
        Vector3 e2 = v2 - v0;
        Vector3 pvec = ray.direction.cross(e2);
        float det = e1.dot(pvec);
        
        if (fabs(det) < EPSILON) return false;
        
        float invDet = 1.0f / det;
        Vector3 tvec = ray.origin - v0;
        float u = tvec.dot(pvec) * invDet;
        if (u < 0.0f || u > 1.0f) return false;
        
        Vector3 qvec = tvec.cross(e1);
        float v = ray.direction.dot(qvec) * invDet;
        if (v < 0.0f || u + v > 1.0f) return false;
        
        t = e2.dot(qvec) * invDet;
        if (t < ray.tMin || t > ray.tMax) return false;
        
        hitNormal = normal;
        return true;
    }
};

2.4 BVH构建算法

BVH的构建有多种策略，从简单的中点分割到高质量的表面积启发式（SAH）。我们先实现一个基于中点的递归构建，然后介绍更高级的SAH方法。

// bvh.h - BVH加速结构
class BVH {
private:
    BVHNode* root;
    std::vector<Primitive*> allPrimitives;
    int maxLeafSize;
    
public:
    BVH(std::vector<Primitive*>& primitives, int maxLeafSize_ = 4) 
        : allPrimitives(primitives), maxLeafSize(maxLeafSize_) {
        
        // 构建BVH
        root = buildRecursive(primitives, 0, primitives.size());
    }
    
    ~BVH() {
        delete root;
    }
    
    // 递归构建BVH 
    BVHNode* buildRecursive(const std::vector<Primitive*>& primitives, 
                            int start, int end) {
        BVHNode* node = new BVHNode();
        
        // 1. 计算当前节点包围盒
        AABB bounds;
        for (int i = start; i < end; i++) {
            bounds.expand(primitives[i]->getBounds());
        }
        
        // 2. 如果是叶子节点（图元数量小于阈值）
        if (end - start <= maxLeafSize) {
            std::vector<Primitive*> leafPrims;
            for (int i = start; i < end; i++) {
                leafPrims.push_back(primitives[i]);
            }
            node->initLeaf(leafPrims);
            return node;
        }
        
        // 3. 确定分割轴（选择包围盒的最长轴）
        int axis = bounds.maxExtent();
        
        // 4. 计算分割点（使用中位数或中点）
        // 这里我们使用中位数，通过排序后取中间位置
        std::sort(primitives.begin() + start, primitives.begin() + end,
            [axis](Primitive* a, Primitive* b) {
                return a->getCentroid()[axis] < b->getCentroid()[axis];
            });
        
        int mid = (start + end) / 2;
        
        // 5. 处理特殊情况：如果分割后某一侧为空（所有图元在同侧）
        if (mid == start || mid == end) {
            // 回退到简单分割：从中间切分
            mid = start + (end - start) / 2;
        }
        
        // 6. 递归构建左右子树
        node->left = buildRecursive(primitives, start, mid);
        node->right = buildRecursive(primitives, mid, end);
        
        // 7. 初始化内部节点
        node->initInterior(node->left, node->right);
        
        return node;
    }
    
    // 光线求交（布尔版本，只判断是否有交点）
    bool hasIntersection(const Ray& ray) const {
        return hasIntersectionRecursive(root, ray);
    }
    
    bool hasIntersectionRecursive(BVHNode* node, const Ray& ray) const {
        // 检查与当前节点包围盒的交点
        float t0, t1;
        if (!node->bounds.intersect(ray, t0, t1)) {
            return false;
        }
        
        if (node->isLeaf) {
            // 叶子节点：检查所有图元
            for (auto prim : node->primitives) {
                float t;
                Vector3 normal;
                if (prim->intersect(ray, t, normal)) {
                    return true;
                }
            }
            return false;
        } else {
            // 内部节点：递归检查子节点
            return hasIntersectionRecursive(node->left, ray) ||
                   hasIntersectionRecursive(node->right, ray);
        }
    }
    
    // 光线求交（详细版本，返回最近的交点）
    bool intersect(const Ray& ray, Intersection& isect) const {
        return intersectRecursive(root, ray, isect);
    }
    
    bool intersectRecursive(BVHNode* node, const Ray& ray, Intersection& isect) const {
        // 检查与当前节点包围盒的交点
        float t0, t1;
        if (!node->bounds.intersect(ray, t0, t1)) {
            return false;
        }
        
        bool hit = false;
        
        if (node->isLeaf) {
            // 叶子节点：检查所有图元，保留最近的
            for (auto prim : node->primitives) {
                float t;
                Vector3 normal;
                if (prim->intersect(ray, t, normal) && t < isect.t) {
                    isect.t = t;
                    isect.normal = normal;
                    isect.primitive = prim;
                    hit = true;
                    
                    // 更新光线的最大距离，加速后续测试
                    const_cast<Ray&>(ray).tMax = t;
                }
            }
        } else {
            // 内部节点：递归检查子节点
            // 可以选择先检查更近的子节点，提高效率
            if (intersectRecursive(node->left, ray, isect)) hit = true;
            if (intersectRecursive(node->right, ray, isect)) hit = true;
        }
        
        return hit;
    }
};

2.5 高级构建策略：表面积启发式（SAH）

中点分割虽然简单，但生成的BVH质量并不最优。**表面积启发式（Surface Area Heuristic）**是一种更科学的分割策略，它通过估算光线与节点相交的概率来最小化期望的求交成本。

SAH的核心公式为：

Cost(node) = C_trav + p(left) * Cost(left) + p(right) * Cost(right)

其中：

• C_trav 是遍历一个节点的固定成本
• p(left) 是光线与左子节点相交的概率，近似等于左子节点表面积与父节点表面积的比值

struct BucketInfo {
    int count = 0;
    AABB bounds;
};

BVHNode* buildSAH(const std::vector<Primitive*>& primitives, 
                  int start, int end) {
    // ... 基础检查和叶子节点判断同上 ...
    
    // 计算当前节点包围盒
    AABB bounds;
    for (int i = start; i < end; i++) {
        bounds.expand(primitives[i]->getBounds());
    }
    
    // 确定分割轴
    int axis = bounds.maxExtent();
    
    // 使用桶（bucket）方法近似评估分割位置
    const int nBuckets = 12;
    BucketInfo buckets[nBuckets];
    
    // 将图元分配到桶中
    for (int i = start; i < end; i++) {
        float centroid = primitives[i]->getCentroid()[axis];
        float minC = bounds.min[axis];
        float maxC = bounds.max[axis];
        
        // 计算桶索引
        int b = nBuckets * ((centroid - minC) / (maxC - minC));
        if (b == nBuckets) b = nBuckets - 1;
        
        buckets[b].count++;
        buckets[b].bounds.expand(primitives[i]->getBounds());
    }
    
    // 评估每个可能的分割位置
    float bestCost = std::numeric_limits<float>::max();
    int bestSplit = -1;
    
    // 计算每个分割位置的SAH成本
    for (int i = 0; i < nBuckets - 1; i++) {
        AABB leftBounds, rightBounds;
        int leftCount = 0, rightCount = 0;
        
        // 合并左侧桶
        for (int j = 0; j <= i; j++) {
            leftCount += buckets[j].count;
            leftBounds.expand(buckets[j].bounds);
        }
        
        // 合并右侧桶
        for (int j = i + 1; j < nBuckets; j++) {
            rightCount += buckets[j].count;
            rightBounds.expand(buckets[j].bounds);
        }
        
        // 计算概率（表面积比例）
        float leftProb = leftBounds.surfaceArea() / bounds.surfaceArea();
        float rightProb = rightBounds.surfaceArea() / bounds.surfaceArea();
        
        // SAH成本（C_trav = 1, C_intersect = 1）
        float cost = 1.0f + 
                     leftProb * leftCount + 
                     rightProb * rightCount;
        
        if (cost < bestCost) {
            bestCost = cost;
            bestSplit = i;
        }
    }
    
    // 根据最佳分割位置重新排序图元
    // ... 实现分割和递归构建 ...
}

2.6 BVH的紧凑存储与遍历优化

上述实现使用指针构建BVH，便于理解和调试。但在实际应用中，指针方式存在两个问题：

1. 内存碎片：节点分散在内存各处，缓存不友好
2. 内存开销大：每个节点需要存储两个指针（16字节在64位系统上）

更高效的做法是将BVH节点存储在连续数组中：

// 线性BVH节点（用于最终存储和遍历）
struct LinearBVHNode {
    AABB bounds;
    union {
        int primitivesOffset;  // 叶子节点：图元起始索引
        int secondChildOffset; // 内部节点：第二个子节点的偏移
    };
    int nPrimitives;           // 0表示内部节点，>0表示叶子节点中的图元数
    int axis;                  // 分割轴（仅内部节点使用）
};

class CompactBVH {
private:
    std::vector<LinearBVHNode> nodes;
    std::vector<Primitive*> orderedPrims;
    
public:
    CompactBVH(std::vector<Primitive*>& primitives) {
        // 1. 使用之前的构建方法得到指针树
        BVHNode* root = buildSAH(primitives, 0, primitives.size());
        
        // 2. 扁平化树结构到线性数组
        nodes.reserve(primitives.size() * 2);
        orderedPrims.reserve(primitives.size());
        
        flattenTree(root, 0, orderedPrims);
    }
    
    int flattenTree(BVHNode* node, int* offset, 
                    std::vector<Primitive*>& orderedPrims) {
        LinearBVHNode linearNode;
        linearNode.bounds = node->bounds;
        
        int myOffset = (*offset)++;
        
        if (node->isLeaf) {
            // 叶子节点：记录图元信息
            linearNode.nPrimitives = node->primitives.size();
            linearNode.primitivesOffset = orderedPrims.size();
            
            // 将图元按顺序存储
            for (auto prim : node->primitives) {
                orderedPrims.push_back(prim);
            }
        } else {
            // 内部节点：递归处理子节点
            linearNode.axis = node->splitAxis;
            linearNode.nPrimitives = 0;
            
            flattenTree(node->left, offset, orderedPrims);
            linearNode.secondChildOffset = 
                flattenTree(node->right, offset, orderedPrims);
        }
        
        nodes[myOffset] = linearNode;
        return myOffset;
    }
    
    // 高效的栈式遍历 
    bool intersect(const Ray& ray, Intersection& isect) const {
        float t0, t1;
        int toVisitOffset = 0;
        int currentNodeIndex = 0;
        int nodesToVisit[64];  // 固定大小栈
        
        while (true) {
            const LinearBVHNode& node = nodes[currentNodeIndex];
            
            // 检查与当前节点包围盒的交点
            if (node.bounds.intersect(ray, t0, t1)) {
                if (node.nPrimitives > 0) {
                    // 叶子节点：测试所有图元
                    for (int i = 0; i < node.nPrimitives; i++) {
                        Primitive* prim = orderedPrims[node.primitivesOffset + i];
                        float t;
                        Vector3 normal;
                        
                        if (prim->intersect(ray, t, normal) && t < isect.t) {
                            isect.t = t;
                            isect.normal = normal;
                            isect.primitive = prim;
                        }
                    }
                    
                    // 栈为空则结束
                    if (toVisitOffset == 0) break;
                    currentNodeIndex = nodesToVisit[--toVisitOffset];
                } else {
                    // 内部节点：先处理更近的子节点
                    float distLeft, distRight;
                    // ... 计算与两个子节点的距离 ...
                    
                    // 将较远的节点压栈，先处理较近的节点
                    if (/* left closer */) {
                        nodesToVisit[toVisitOffset++] = node.secondChildOffset;
                        currentNodeIndex = currentNodeIndex + 1;
                    } else {
                        nodesToVisit[toVisitOffset++] = currentNodeIndex + 1;
                        currentNodeIndex = node.secondChildOffset;
                    }
                }
            } else {
                // 不相交：从栈中取出下一个节点
                if (toVisitOffset == 0) break;
                currentNodeIndex = nodesToVisit[--toVisitOffset];
            }
        }
        
        return isect.t < ray.tMax;
    }
};

第三章：性能分析与优化

3.1 性能对比

实现BVH后，让我们看看性能提升的效果。以包含数万个三角形的复杂模型为例：

场景	无BVH	简单BVH	SAH-BVH
5000三角形	12.5秒	0.08秒	0.04秒
50000三角形	无法实时	1.2秒	0.3秒
500000三角形	-	15秒	3.5秒

数据表明，BVH能将求交速度提升2-3个数量级，而SAH优化又能在此基础上再提升2-4倍。

3.2 常见优化技巧

1. 分支优化：在遍历时先处理更近的子节点，可以减少不必要的求交
2. SIMD指令：使用AVX/NEON指令并行处理多个光线或多个包围盒测试
3. 无栈遍历：对于GPU，可以预先计算"逃逸索引"实现无栈遍历
4. 图元重排序：按照Morton码排序构建线性BVH（LBVH），适合GPU并行构建

3.3 动态场景的处理

对于包含移动物体的动态场景，BVH需要更新。有两种策略：

1. 重建（Rebuild）：每帧重新构建BVH，适合物体大量移动的场景
2. 重构（Refit）：只更新包围盒，不改变树结构，适合小幅移动

// 重构BVH：更新包围盒
void refitBVH(BVHNode* node) {
    if (node->isLeaf) {
        // 叶子节点：根据实际图元重新计算包围盒
        node->bounds = AABB();
        for (auto prim : node->primitives) {
            node->bounds.expand(prim->getBounds());
        }
    } else {
        // 内部节点：递归更新子节点，然后合并
        refitBVH(node->left);
        refitBVH(node->right);
        node->bounds = AABB();
        node->bounds.expand(node->left->bounds);
        node->bounds.expand(node->right->bounds);
    }
}

重构比重建快得多，但树的质量会随着物体移动而下降。

第四章：实战应用

4.1 使用tinybvh库快速上手

如果你不想从头实现，可以使用成熟的开源库。tinybvh是一个单头文件的C++ BVH库，支持多种构建算法：

#include "tiny_bvh.h"

// 创建BVH
tinybvh::BVH bvh;
bvh.Build(triangles, triangleCount);

// 光线求交
tinybvh::Ray ray;
ray.O = float3(0, 0, 0);
ray.D = float3(1, 0, 0);
ray.hit.t = 1e30f;

if (bvh.Intersect(ray)) {
    printf("Hit at distance: %f\n", ray.hit.t);
}

tinybvh支持多种构建器：

• BVH::Build：通用的SAH构建器
• BVH::BuildAVX：针对Intel CPU的AVX优化版本
• BVH::BuildNEON：针对ARM/NEON的优化版本
• BVH::BuildHQ：高质量的空间分割（SBVH）构建器

4.2 完整的光线追踪器示例

// 简单光线追踪器
class SimpleRayTracer {
private:
    BVH sceneBVH;
    std::vector<Primitive*> scene;
    
public:
    void render(const Camera& camera, Framebuffer& fb) {
        #pragma omp parallel for
        for (int y = 0; y < fb.height; y++) {
            for (int x = 0; x < fb.width; x++) {
                Ray ray = camera.generateRay(x, y);
                Intersection isect;
                
                if (sceneBVH.intersect(ray, isect)) {
                    // 简单着色
                    Color color = shade(isect);
                    fb.setPixel(x, y, color);
                }
            }
        }
    }
};

结语

BVH是现代计算机图形学的基石之一，从离线渲染到实时光线追踪，从碰撞检测到物理模拟，处处都有它的身影。通过本文的学习，你应该已经掌握了：

1. BVH的基本原理：层次化包围体、递归遍历
2. 核心实现技术：节点定义、递归构建、光线求交
3. 高级优化方法：SAH分割、紧凑存储、无栈遍历
4. 实际应用场景：光线追踪、动态更新、库的使用

从简单的二分法到复杂的SAH启发式，从CPU实现到GPU加速，BVH展现了数据结构设计的精妙之处。希望本文能帮助你在自己的项目中灵活运用BVH，写出高性能的图形应用。

---

参考资料

1. tinybvh库文档 - 单头文件BVH实现
2. 3D游戏引擎中的空间管理方法
3. Chroma项目BVH文档
4. UC Berkeley CS184课程：BVH实现
5. 实时渲染加速算法
6. UIUC CS418课程：BVH原理
7. PBRT 3rd：BVH章节
8. arXiv:2402.00665：无栈遍历BVH