Bounding Volume Hierachy of pbrt 解析(1)

Tags:

二叉树BVH

BVH是空间切分技术之一,除了BVH之外还有kdtree、octree。

下面先以静态场景为例,讲解BVH的生成算法。

静态场景生成BVH

初始化

首先要确保场景里的3d对象都可以算出它的AABB盒(下面简称BBox)(模型空间),然后经过scene graph以及自身的模型矩阵,把包围盒变换到世界空间,称之为world bound。

代码片段——找出待计算的所有3d对象:

std::vector<ObjectID> objs; // 所有3d对象的id列表
std::function<void(Object objScene)> filterFunc;
filterFunc = [&](Object objScene) {
    if (objScene.hasComponent<SpatialData>() && objScene.hasComponent<MeshRef>()) {
        ObjectID objID = objScene.ID();
        objs.push_back(objID);
    }
    auto sgNode = objScene.component<SceneGraphNode>();
    for (auto childObjID : sgNode->children) {
        Object childObj = m_objMgr->get(childObjID);
        filterFunc(childObj); // 树递归
    }
}; 
filterFunc(objScene);

代码片段——预先算出所有物体的world bound:

std::vector<BVHObjInfo> objInfo(bvhAccel->objs.size());
for (size_t i = 0; i < bvhAccel->objs.size(); ++i) {
    ObjectID objID = bvhAccel->objs[i];
    Object obj = m_objMgr->get(objID);
    auto meshRef = obj.component<MeshRef>();
    auto spatialData = obj.component<SpatialData>();
    Mesh& mesh = meshSet->getMesh(meshRef->meshID);
    auto bound = mesh.Bound();
    auto worldBound = (spatialData->o2w)(bound);
    objInfo[i] = { i, worldBound };
}

代码片段——BVH生成流程:

int totalNodes = 0; // 统计总共动态创建了多少个BVHBuildNode
std::vector<ObjectID> orderedObjs;
orderedObjs.reserve(bvhAccel->objs.size()); // 分配空间 
BVHBuildNode *root;
// 2种构造算法,选1
if (splitMethod == BVHAccel::SplitMethod::HLBVH)
    root = HLBVHBuild(bvhAccel, arena, objInfo, &totalNodes, orderedObjs);
else
    root = recursiveBuild(bvhAccel, arena, objInfo, 0, bvhAccel->objs.size(),
        &totalNodes, orderedObjs);
// 构造完毕
bvhAccel->objs.swap(orderedObjs);

// 计算用深度优先遍历表示的BVH二叉树数组
bvhAccel->nodes = AllocAligned<LinearBVHNode>(totalNodes);
int offset = 0;
flattenBVHTree(bvhAccel, root, &offset);

可以发现,BVH的核心算法代码就在recursiveBuild和HLBVHBuild里面了。

相关数据结构

BVHBuildNode:

// 生成BVH的过程中会用到这个struct,代表BVH的一个节点
struct BVHBuildNode {
    BBox bounds; // 该子树的包围盒 
    BVHBuildNode *children[2]; // 左右孩子节点
    Axis splitAxis; // 是切了哪条坐标轴
    int firstObjOffset, nObjs; // 如果是叶子节点的话,这2个变量记录了被包含的3d物体
};

// 把node初始化为叶子节点
void BVHSystem::InitLeaf(BVHBuildNode * node, int first, int n, const BBox &b) {
    node->firstObjOffset = first;
    node->nObjs = n;
    node->bounds = b;
    node->children[0] = node->children[1] = nullptr;
}

// 把node初始化为内部节点
void BVHSystem::InitInterior(BVHBuildNode * node, Axis axis, BVHBuildNode *c0, BVHBuildNode *c1) {
    node->children[0] = c0;
    node->children[1] = c1;
    node->bounds = Union(c0->bounds, c1->bounds);
    node->splitAxis = axis;
    node->nObjs = 0;
}

BVHObjInfo:

struct BVHObjInfo {
    BVHObjInfo() {}
    BVHObjInfo(size_t objNumber, const BBox &bounds)
        : objNumber(objNumber),
        bounds(bounds),
        centroid(.5f * bounds.pMin + .5f * bounds.pMax) {}
    size_t objNumber; // 记录这个3d对象在bvhAccel->objs数组的位置
    BBox bounds;
    Vector3dF centroid; // 包围盒中心坐标
};

recursiveBuild

BVHBuildNode *BVHSystem::recursiveBuild(
        ComponentHandle<BVHAccel> bvhAccel, 
        MemoryArena &arena,
        std::vector<BVHObjInfo> &objInfo, 
        int start,
        int end, 
        int *totalNodes,
        std::vector<ObjectID> &orderedObjs) {
    BVHBuildNode *node = arena.Alloc<BVHBuildNode>();
    (*totalNodes)++; // 计数

    BBox bounds; // [start, end) 区间内所有物体的包围盒的并集包围盒
    for (int i = start; i < end; ++i)
        bounds = Union(bounds, objInfo[i].bounds);

    int nObjs = end - start;// [start, end)区间内的物体数量
    if (nObjs == 1) { // 如果只有一个物体,那就把node设置为叶节点
        int firstObjOffset = orderedObjs.size();
        for (int i = start; i < end; ++i) {
            int objNum = objInfo[i].objNumber;
            orderedObjs.push_back(bvhAccel->objs[objNum]);
        }
        InitLeaf(node, firstObjOffset, nObjs, bounds);
        return node;
    } else {
        ······
        // B
    }

    return node;
}

B:

// 计算[start, end)里的物体包围盒中心坐标的包围盒
BBox centroidBounds;
for (int i = start; i < end; ++i)
    centroidBounds = Union(centroidBounds, objInfo[i].centroid);
//  把边长最长的坐标轴作为切分轴
Axis dim = centroidBounds.MaximumExtent();

// Partition objs into two sets and build children
if (centroidBounds.pMax[dim] == centroidBounds.pMin[dim]) {
    // 在边长最长的轴上,竟然pMax和pMin还相等,说明包围盒中心都重叠了,
    // 把这些物体都归入同个叶子节点
    int firstObjOffset = orderedObjs.size();
    for (int i = start; i < end; ++i) {
        int objNum = objInfo[i].objNumber;
        orderedObjs.push_back(bvhAccel->objs[objNum]);
    }
    InitLeaf(node, firstObjOffset, nObjs, bounds);
    // (和上面那段代码一模一样)
    return node;
}
else {
    ······
    // C
}

到了C,就真的要做空间切分了,需要根据SplitMethod选择切分算法。

C:

int mid = (start + end) / 2; // 默认mid

// 这个switch的作用就是把[start,end)的物体切分开,并把切分位置记到mid
// 从而可以对切出来的2个子树继续build
switch (bvhAccel->splitMethod) {
    case BVHAccel::SplitMethod::Middle: { 
        // 中心切分法
        // pmid是最长边的边中心点
        float pmid =
            (centroidBounds.pMin[dim] + centroidBounds.pMax[dim]) / 2;
        // 把[start, end)里的物体按pmid切分,小于pmid的在前面
        // 注意,partition并不是排序
        BVHObjInfo *midPtr = std::partition(
            &objInfo[start], &objInfo[end - 1] + 1,
            [dim, pmid](const BVHObjInfo &pi) {
            return pi.centroid[dim] < pmid;
        });
        mid = midPtr - &objInfo[0]; // midPtr指向了pmid右侧第一个对象,然后转成数组下标mid
        if (mid != start && mid != end) break;
        // 如果mid等于start或end,说明这些物体要么都在pmid左侧,要么都在右侧,切分无意义
        // 所以不能break,继续进入EqualCounts切分
        // (但按照前面的逻辑来看,这里不应该一个都切分不出来,可能是一个保险措施)
    }
    case BVHAccel::SplitMethod::EqualCounts: {
        // 按照数量对半切
        // nth_element是偏排序,保证比mid小的都排在mid之前,比mid大的都在mid之后,但不保证有序
        mid = (start + end) / 2;
        std::nth_element(&objInfo[start], &objInfo[mid],
            &objInfo[end - 1] + 1,
            [dim](const BVHObjInfo &a,
                const BVHObjInfo &b) {
            return a.centroid[dim] < b.centroid[dim];
        });
        break;
    }
    case BVHAccel::SplitMethod::SAH:
    default: {
        ····
        // D
        break;
    }
}
// 递归构造子树
InitInterior(node, dim,
    recursiveBuild(bvhAccel, arena, objInfo, start, mid,
        totalNodes, orderedObjs),
    recursiveBuild(bvhAccel, arena, objInfo, mid, end,
        totalNodes, orderedObjs));

D,是基于概率、表面积的启发式切分法SAH(Surface Area Heuristic),要用到一条公式:

\[ cost(A, B) = t_{trav} + P_{A} \sum _{i=1}^{N_{A} } t_{isect}(a_{i}) + P_{B} \sum _{i=1}^{N_{B} } t_{isect}(b_{i}) \]

cost(A, B)是ray和这个内部节点做相交计算的时间开销;

\( P_{A} 、 P_{B}\) 分别是ray经过A、B子树的概率;

\( N_{A} 、N_{B}\) 分别是A、B子树的物体的数量;

\( t_{trav} \) 是遍历一个节点的时间开销,\( t_{isect}(a_{i}) 、 t_{isect}(b_{i}) \) 分别是ray和子树某个对象包围盒做相交计算的时间开销。可以简单假设\(t_{trav}、t_{isect}\)都是常数1。那么上式简化成:

\[ cost(A, B) = 1 + P_{A} \sum _{i=1}^{N_{A} } 1 + P_{B} \sum _{i=1}^{N_{B} } 1 = 1 + P_{A} N_{A} + P_{B} N_{B} \]

\( P_{A} 、 P_{B}\) 可以用表面积比来求出:

\[ P_{A} = \frac {S_{A} }{S_{parent } }\]

\[ P_{B} = \frac {S_{B} }{S_{parent } }\]

然后就是应用问题,因为切分的最优位置并不能直接得到,需要对所有可切分点都做这条公式算出cost,然后取那个最小值。

显然这很暴力,所以pbrt作者设计了bucket机制,其实就是把相邻的物体捆绑打包,当作一个整体,再来遍历各个切分点算cost。

if (nObjs <= 2) { // 如果用了SAH但这个节点只有2个物体,就没必要算SAH了,直接切
    mid = (start + end) / 2;
    std::nth_element(&objInfo[start], &objInfo[mid],
        &objInfo[end - 1] + 1,
        [dim](const BVHObjInfo &a,
            const BVHObjInfo &b) {
        return a.centroid[dim] <
            b.centroid[dim];
    });
}
else {
    constexpr int nBuckets = 12;
    BucketInfo buckets[nBuckets];

    // 初始化bucket
    for (int i = start; i < end; ++i) {
        int b = nBuckets *
            centroidBounds.Offset(
                objInfo[i].centroid)[dim];
        if (b == nBuckets) b = nBuckets - 1;
        Assert(b >= 0 && b < nBuckets);
        buckets[b].count++;
        buckets[b].bounds =
            Union(buckets[b].bounds, objInfo[i].bounds);
    }

    // 计算切分开销
    float cost[nBuckets - 1];
    for (int i = 0; i < nBuckets - 1; ++i) {
        BBox b0, b1;
        int count0 = 0, count1 = 0;
        for (int j = 0; j <= i; ++j) {
            b0 = Union(b0, buckets[j].bounds);
            count0 += buckets[j].count;
        }
        for (int j = i + 1; j < nBuckets; ++j) {
            b1 = Union(b1, buckets[j].bounds);
            count1 += buckets[j].count;
        }
        cost[i] = 1 +
            (count0 * b0.SurfaceArea() +
                count1 * b1.SurfaceArea()) /
            bounds.SurfaceArea();
    }

    // 线性遍历,找出cost最小的
    float minCost = cost[0];
    int minCostSplitBucket = 0;
    for (int i = 1; i < nBuckets - 1; ++i) {
        if (cost[i] < minCost) {
            minCost = cost[i];
            minCostSplitBucket = i;
        }
    }
    float leafCost = nObjs;
    // 如果数量比maxObjsInNode大 或者 切分开销比直接弄成叶子节点开销小,就做切分
    if (nObjs > bvhAccel->maxObjsInNode || minCost < leafCost) {
        BVHObjInfo *pmid = std::partition(
            &objInfo[start], &objInfo[end - 1] + 1,
            [=](const BVHObjInfo &pi) {
            int b = nBuckets *
                centroidBounds.Offset(pi.centroid)[dim];
            if (b == nBuckets) b = nBuckets - 1;
            Assert(b >= 0 && b < nBuckets);
            return b <= minCostSplitBucket;
        });
        mid = pmid - &objInfo[0];
    }
    else {
        // 否则,弄成叶子节点就行了(这段代码出现第三次了)
        int firstObjOffset = orderedObjs.size();
        for (int i = start; i < end; ++i) {
            int objNum = objInfo[i].objNumber;
            orderedObjs.push_back(bvhAccel->objs[objNum]);
        }
        InitLeaf(node, firstObjOffset, nObjs, bounds);
        return node;
    }
}
(未经授权禁止转载)
Written on December 20, 2017

博主将十分感谢对本文章的任意金额的打赏^_^