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];
}
}