二叉树BVH
BVH是空间切分技术之一，除了BVH之外还有kdtree、octree。
下面先以静态场景为例，讲解BVH的生成算法。Note：演示代码是用pbrt源码改的测试版，这是为了理解代码和方便调试。
静态场景生成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
写作不易，您的支持是我写作的动力！
