PBRT笔记(2)——BVH
BVH
构建BVH树分三步：
- 计算每个图元的边界信息并且存储在数组中
- 使用指定的方法构建树
- 优化树，使得树更加紧凑
//BVH边界信息，存储了图元号，包围盒以及中心点
struct BVHPrimitiveInfo {
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;
};
分割
使用，子图元中质心距离最大的轴向作为分割方向。（另一种方法是尝试所有轴，之后再选择效果最好的那个轴作为分割方向。但是在实践中发现当前方案也有着不错的效果）
int nPrimitives = end - start;
//只有一个图元，则生成叶子节点并且返回
if (nPrimitives == 1) {
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;
}
int mid = (start + end) / 2;
//如果多个片元的组成的聚合体是0空间体，则生成叶子节点，并且返回（这是一种不寻常的现象）
if (centroidBounds.pMax[dim] == centroidBounds.pMin[dim]) {
// Create leaf _BVHBuildNode_
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;
}
int mid = (start + end) / 2;
if (centroidBounds.pMax[dim] == centroidBounds.pMin[dim]) {
<创建叶子节点>
} else {
<使用对应方法切割图元>
//使用递归构造BVH树
node->InitInterior(dim,
recursiveBuild(arena, primitiveInfo, start, mid,totalNodes,orderedPrims),
recursiveBuild(arena, primitiveInfo, mid, end,totalNodes, orderedPrims));
}
//中间对半切的方法
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];
if (mid != start && mid != end) break;
}
这里遇到指针相减的操作
如果两个指针指向同一个数组,它们就可以相减,其结果为两个指针之间的元素数目。同理+1则表示内存偏移一个该类型空间。
[end - 1] + 1是会为了回避容器越界访问的问题，但是取了地址再偏移
end操作返回的迭代器指向容器的“末端元素的下一个”，指向了一个不存在的元素
但是如果有多个图元的边界盒处于与左右节点的边界盒重叠的状态，那分割很可能会失败。SplitMethod::EqualCounts排序速度会更快些
case SplitMethod::EqualCounts: {
// Partition primitives into equally-sized subsets
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;
}
case SplitMethod::SAH:
default: {
//图元数量较少时没必要使用SAH
if (nPrimitives <= 2) {
// Partition primitives into equally-sized subsets
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];
});
} else {
// Allocate _BucketInfo_ for SAH partition buckets
PBRT_CONSTEXPR int nBuckets = 12;
BucketInfo buckets[nBuckets];
//计算出当前图元的质心处于第几个桶，centroidBounds为当前节点中所有图元的边界盒
//之后进行BucketInfo统计，并且调整对应桶的边界盒
for (int i = start; i < end; ++i) {
int b = nBuckets *
centroidBounds.Offset(
primitiveInfo[i].centroid)[dim];
if (b == nBuckets) b = nBuckets - 1;
CHECK_GE(b, 0);
CHECK_LT(b, nBuckets);
buckets[b].count++;
buckets[b].bounds =
Union(buckets[b].bounds, primitiveInfo[i].bounds);
}
Float cost[nBuckets - 1];
for (int i = 0; i < nBuckets - 1; ++i) {
Bounds3f b0, b1;
int count0 = 0, count1 = 0;
//第一个桶到[i]桶，所有的图元数与边界盒
for (int j = 0; j <= i; ++j) {
b0 = Union(b0, buckets[j].bounds);
count0 += buckets[j].count;
}
//[i+1]桶到最后一个桶，所有的图元数与边界盒
for (int j = i + 1; j < nBuckets; ++j) {
b1 = Union(b1, buckets[j].bounds);
count1 += buckets[j].count;
}
//bounds变量为所有图元的边界盒
//通过面积模型计算，相交开销设为1
cost[i] = 1 +
(count0 * b0.SurfaceArea() +
count1 * b1.SurfaceArea()) /
bounds.SurfaceArea();
}
//计算出最小消耗的切割位置与消耗的量
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 = nPrimitives;
//图元数超过最大值（255）或者最小消耗小于所有图元都创建叶子节点的消耗（切割具有良好效果的情况）
//使用partition算法，按使用SAH找到最小消耗的切割位置，对桶进行排序，之后进行下一次遍历
//不然就直接生成节点，结束遍历
if (nPrimitives > maxPrimsInNode || minCost < leafCost) {
BVHPrimitiveInfo *pmid = std::partition(
&primitiveInfo[start], &primitiveInfo[end - 1] + 1,
[=](const BVHPrimitiveInfo &pi) {
int b = nBuckets *
centroidBounds.Offset(pi.centroid)[dim];
if (b == nBuckets) b = nBuckets - 1;
CHECK_GE(b, 0);
CHECK_LT(b, nBuckets);
return b <= minCostSplitBucket;
});
mid = pmid - &primitiveInfo[0];
} else {
// Create leaf _BVHBuildNode_
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;
}
}
break;
}
紧凑BVH树
使用SAH有两个缺点：
- 花费过多时间在使用SAH构建树上。
- 自上而下的BVH树的构建很难并行化。有一个解决方法就是构建若干个独立子树，不过这反过来限制了并行的可伸缩性。这也是GPU渲染所要面对的问题。
Linear bounding volume hierarchies (LBVHs) 就是为了解决这个问题而开发的。
LinearBVHNode大小为32字节，以满足32位内存对齐要求。
LBVH是基于莫顿码，讲原本的多维数据转换为排序过的一维的数据。
也就是将树中节点的相对位置按照规律排序：
如果用uint8来表示一个二叉树的2维状态：
$ y_3 x_3 y_2 x_2 y_1 x_1 y_0 x_0$
struct LinearBVHNode {
Bounds3f bounds;
union {
int primitivesOffset; // leaf
int secondChildOffset; // interior
};
uint16_t nPrimitives; // 0 -> interior node
uint8_t axis; // interior node: xyz
uint8_t pad[1]; // ensure 32 byte total size
};
hierarchical linear bounding volume hierarchy (HLBVH)
讲质心坐标转化为莫顿码，之后统计对应莫顿码（桶中）图元数量，之后以莫顿码作为index，将图元放入容器的对应位置，从而完成排序（使用了基数排序）。
ParallelFor([&](int i) {
// Initialize _mortonPrims[i]_ for _i_th primitive
//mortonScale=1024;
PBRT_CONSTEXPR int mortonBits = 10;
PBRT_CONSTEXPR int mortonScale = 1 << mortonBits;
mortonPrims[i].primitiveIndex = primitiveInfo[i].primitiveNumber;
//因为bounds.Offset返回的是[0,1]区间的百分比，所以为了转换成莫顿码需要再乘以一个比例系数
//PBRT使用int32来存储莫顿码，因为需要存储x、y、z三个维度，所以每个维度占用10个bit，所以比例系数为10，对于二进制莫顿码来说乘以10等于左移10个位，也就是2^10=1024
Vector3f centroidOffset = bounds.Offset(primitiveInfo[i].centroid);
//以边界盒的min为零点建立坐标系，使用质心位置*莫顿码缩放值，计算莫顿码
//EncodeMorton3，使用LeftShift3()分别计算x、y、z，之后再将y、z分别往左偏移1与2位
//LeftShift3()参看书中的图就可以明白了
mortonPrims[i].mortonCode = EncodeMorton3(centroidOffset * mortonScale);
}, primitiveInfo.size(), 512);
之后以索引对莫顿码进行排序（为了追求效率而没有选择std::sort）,这里如果不明白基数排序，就很难看懂。
for (int pass = 0; pass < nPasses; ++pass) {
//pass如果各个位都为1，着in为tempVector的引用，否则则为v
//每次循环out与in都进行互换，将上一次排序结果接着进行排序（第一次直接用外部未排序的容器，之后将每个pass都排序一遍，所有莫顿码即排序完成）
std::vector<MortonPrimitive> &in = (pass & 1) ? tempVector : *v;
std::vector<MortonPrimitive> &out = (pass & 1) ? *v : tempVector;
}
//一共64种可能的莫顿码
PBRT_CONSTEXPR int nBuckets = 1 << bitsPerPass;
int bucketCount[nBuckets] = {0};
PBRT_CONSTEXPR int bitMask = (1 << bitsPerPass) - 1;
for (const MortonPrimitive &mp : in) {
int bucket = (mp.mortonCode >> lowBit) & bitMask;
CHECK_GE(bucket, 0);
CHECK_LT(bucket, nBuckets);
//取得当前的pass的莫顿码，并且进行统计
++bucketCount[bucket];
}
//计算每个桶到第一个桶的莫顿码总量，也就是index的偏移量，毕竟当前位置的index+该位置桶内元素数量，就等于下一个桶到第一个桶的偏移量。
int outIndex[nBuckets];
outIndex[0] = 0;
for (int i = 1; i < nBuckets; ++i)
outIndex[i] = outIndex[i - 1] + bucketCount[i - 1];
//通过莫顿码讲图元节点放到对应的位置，从而完成排序。++是为了偏移放置下一个该桶元素的位置
for (const MortonPrimitive &mp : in) {
int bucket = (mp.mortonCode >> lowBit) & bitMask;
out[outIndex[bucket]++] = mp;
}
// 寻找每个小数中图元的间隔
std::vector<LBVHTreelet> treeletsToBuild;
for (int start = 0, end = 1; end <= (int)mortonPrims.size(); ++end) {
#ifdef PBRT_HAVE_BINARY_CONSTANTS
uint32_t mask = 0b00111111111111000000000000000000;
#else
uint32_t mask = 0x3ffc0000;
#endif
//遍历所有莫顿码高位不同的图元
if (end == (int)mortonPrims.size() ||
((mortonPrims[start].mortonCode & mask) !=
(mortonPrims[end].mortonCode & mask))) {
// Add entry to _treeletsToBuild_ for this treelet
int nPrimitives = end - start;
int maxBVHNodes = 2 * nPrimitives;
BVHBuildNode *nodes = arena.Alloc<BVHBuildNode>(maxBVHNodes, false);
treeletsToBuild.push_back({start, nPrimitives, nodes});
start = end;
}
}
之后对16个区块进行构建子树
emitLBVH
CHECK_GT(nPrimitives, 0);
//如果图元小于一定数量或者 则创建叶子节点
if (bitIndex == -1 || nPrimitives < maxPrimsInNode) {
(*totalNodes)++;
BVHBuildNode *node = buildNodes++;
Bounds3f bounds;
int firstPrimOffset = orderedPrimsOffset->fetch_add(nPrimitives);
for (int i = 0; i < nPrimitives; ++i) {
int primitiveIndex = mortonPrims[i].primitiveIndex;
orderedPrims[firstPrimOffset + i] = primitives[primitiveIndex];
bounds = Union(bounds, primitiveInfo[primitiveIndex].bounds);
}
node->InitLeaf(firstPrimOffset, nPrimitives, bounds);
return node;
} else {
// 从高位开始分割 （29-12）
int mask = 1 << bitIndex;
//如果所有图元都集中在分割的一边，则开始分割下一个子树
if ((mortonPrims[0].mortonCode & mask) ==
(mortonPrims[nPrimitives - 1].mortonCode & mask))
return emitLBVH(buildNodes, primitiveInfo, mortonPrims, nPrimitives,
totalNodes, orderedPrims, orderedPrimsOffset,
bitIndex - 1);
// Find LBVH split point for this dimension
int searchStart = 0, searchEnd = nPrimitives - 1;
while (searchStart + 1 != searchEnd) {
CHECK_NE(searchStart, searchEnd);
int mid = (searchStart + searchEnd) / 2;
//用二分法在这个轴线寻找分割点
if ((mortonPrims[searchStart].mortonCode & mask) ==
(mortonPrims[mid].mortonCode & mask))
searchStart = mid;
else {
CHECK_EQ(mortonPrims[mid].mortonCode & mask,
mortonPrims[searchEnd].mortonCode & mask);
searchEnd = mid;
}
}
int splitOffset = searchEnd;
CHECK_LE(splitOffset, nPrimitives - 1);
CHECK_NE(mortonPrims[splitOffset - 1].mortonCode & mask,
mortonPrims[splitOffset].mortonCode & mask);
//因为已经用莫顿码排序过，所以这里直接递归分割
(*totalNodes)++;
BVHBuildNode *node = buildNodes++;
BVHBuildNode *lbvh[2] = {
emitLBVH(buildNodes, primitiveInfo, mortonPrims, splitOffset,
totalNodes, orderedPrims, orderedPrimsOffset,
bitIndex - 1),
emitLBVH(buildNodes, primitiveInfo, &mortonPrims[splitOffset],
nPrimitives - splitOffset, totalNodes, orderedPrims,
orderedPrimsOffset, bitIndex - 1)};
int axis = bitIndex % 3;
node->InitInterior(axis, lbvh[0], lbvh[1]);
return node;
}
所有子树都创建后，buildUpperSAH将会构建所有子树的BVH。这里的操作和SAH差不多。之后就是flattenBVHTree，以深度优先顺序，对所有节点进行排序。
以下是便利过程：
while (true) {
const LinearBVHNode *node = &nodes[currentNodeIndex];
//是否与当前节点的边界盒相交
if (node->bounds.IntersectP(ray, invDir, dirIsNeg)) {
if (node->nPrimitives > 0) {
//如果是叶子节点，就对节点下的所有图元进行相交测试
for (int i = 0; i < node->nPrimitives; ++i)
if (primitives[node->primitivesOffset + i]->Intersect(
ray, isect))
hit = true;
//如果已经遍历了另外所需遍历的节点，则停止循环，不然设置下一个循环NodeIndex
if (toVisitOffset == 0) break;
currentNodeIndex = nodesToVisit[--toVisitOffset];
} else {
//如果是子树节点，先判断方向正负，是正就访问左节点，并且将右节点放入需要遍历的数组中，待之后的循环进行相交测试
if (dirIsNeg[node->axis]) {
nodesToVisit[toVisitOffset++] = currentNodeIndex + 1;
currentNodeIndex = node->secondChildOffset;
} else {
nodesToVisit[toVisitOffset++] = node->secondChildOffset;
currentNodeIndex = currentNodeIndex + 1;
}
}
} else {
if (toVisitOffset == 0) break;
currentNodeIndex = nodesToVisit[--toVisitOffset];
}
}
