games101 作業6 詳解SAH
可參考資料：https://www.pbr-book.org/3ed-2018/Primitives_and_Intersection_Acceleration/Bounding_Volume_Hierarchies
程式碼分析
作業6 的程式碼框架在作業五的基礎上進行了進一步地豐富與改進
一開始還是建立場景 新增物體以及燈光 這次的物體只有一個用objloader讀取的兔子的3d模型 但是在將該模型轉換為meshtriangle資料結構的時候，也為該模型中的各個三角形片元建立BVH（Bounding volume hierarchies）的加速結構 這裡雖然也為場景中的各個物體也構建了BVH 但是目前場景中只有一個物體
MeshTriangle bunny(Utils::PathFromAsset("model/bunnyAssignment6/bunny.obj"));
scene.Add(&bunny);
scene.Add(std::make_unique<Light>(Vector3f(-20, 70, 20), 1));
scene.Add(std::make_unique<Light>(Vector3f(20, 70, 20), 1));
scene.buildBVH();
MeshTriangle(const std::string& filename){
.......
std::vector<Object*> ptrs;
for (auto& tri : triangles)
ptrs.push_back(&tri);
bvh = new BVHAccel(ptrs);
}
BVHAccel::BVHAccel(std::vector<Object*> p, int maxPrimsInNode,
SplitMethod splitMethod)
: maxPrimsInNode(std::min(255, maxPrimsInNode)), splitMethod(splitMethod),
primitives(std::move(p))
{
time_t start, stop;
time(&start);
if (primitives.empty())
return;
root = recursiveBuild(primitives);
time(&stop);
double diff = difftime(stop, start);
int hrs = (int)diff / 3600;
int mins = ((int)diff / 60) - (hrs * 60);
int secs = (int)diff - (hrs * 3600) - (mins * 60);
printf(
"\rBVH Generation complete: \nTime Taken: %i hrs, %i mins, %i secs\n\n",
hrs, mins, secs);
}
而BVH的構建過程就是一顆樹的遞迴構建的過程 劃分BVH節點的方法 就是在最長的維度上 對所有的三角形進行排序 找出位於中間的這個三角形的位置 左邊的就是左節點 右邊就是右節點，劃分結束的終止條件就是隻有一個物體
對應課件的ppt：
pbrt圖示：
BVHBuildNode* BVHAccel::recursiveBuild(std::vector<Object*> objects)
{
BVHBuildNode* node = new BVHBuildNode();
// Compute bounds of all primitives in BVH node
Bounds3 bounds;
for (int i = 0; i < objects.size(); ++i)
bounds = Union(bounds, objects[i]->getBounds());
if (objects.size() == 1) {
// Create leaf _BVHBuildNode_
node->bounds = objects[0]->getBounds();
node->object = objects[0];
node->left = nullptr;
node->right = nullptr;
return node;
}
else if (objects.size() == 2) {
node->left = recursiveBuild(std::vector{objects[0]});
node->right = recursiveBuild(std::vector{objects[1]});
node->bounds = Union(node->left->bounds, node->right->bounds);
return node;
}
else {
//用每個三角形對應的box的中點 拼成一個box 然後找出這個box最長的維度 在這個維度上進行劃分
Bounds3 centroidBounds;
for (int i = 0; i < objects.size(); ++i)
centroidBounds =
Union(centroidBounds, objects[i]->getBounds().Centroid());
int dim = centroidBounds.maxExtent();
switch (dim) {
case 0:
std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
return f1->getBounds().Centroid().x <
f2->getBounds().Centroid().x;
});
break;
case 1:
std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
return f1->getBounds().Centroid().y <
f2->getBounds().Centroid().y;
});
break;
case 2:
std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
return f1->getBounds().Centroid().z <
f2->getBounds().Centroid().z;
});
break;
}
auto beginning = objects.begin();
auto middling = objects.begin() + (objects.size() / 2);
auto ending = objects.end();
auto leftshapes = std::vector<Object*>(beginning, middling);
auto rightshapes = std::vector<Object*>(middling, ending);
assert(objects.size() == (leftshapes.size() + rightshapes.size()));
node->left = recursiveBuild(leftshapes);
node->right = recursiveBuild(rightshapes);
node->bounds = Union(node->left->bounds, node->right->bounds);
}
return node;
}
之後就是光線追蹤castRay 求交 然後進行著色計算 和作業5類似
但是這次求交要利用我們構建的BVH加速結構進行加速 這也是我們要完成的部分
分析與解決
與BVH求交
這裡主要是對構建的BVH加速結構進行一個遍歷 和節點的boundingbox求交 如果沒有交點就返回 如果有 則檢視該節點是不是葉子節點 不是葉子節點 則繼續遍歷該節點的子節點 如果是葉子節點 則對其中的物體進行相交測試 找出最近的交點返回：
Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
// TODO Traverse the BVH to find intersection
Intersection hit1;
Intersection hit2;
std::array<int, 3> dirIsNeg{ 0,0,0 };
for (int i = 0; i < 3; i++)
{
if (ray.direction[i] < 0) {
dirIsNeg[i] = 1;
}
}
if (!node->bounds.IntersectP(ray, ray.direction_inv, std::move(dirIsNeg)))
{
return Intersection{};
}
if (node->left == nullptr && node->right == nullptr) {
return node->object->getIntersection(ray);
}
hit1 = getIntersection(node->left, ray);
hit2 = getIntersection(node->right, ray);
return hit1.distance > hit2.distance ? hit2 : hit1;
}
boundingbox求交
boundingbox是個長方體 並且與各個面與xyz軸平行 我們只需要求出和每組平面 的tenter 與 texit 並從enter中找出最大值 exit中找出最小值 這樣可以防止錯誤的將交點找到長方體的延展面上 同時也要考慮dir光線為負的情況 這樣的話會和座標值較小的平面先相交：
inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
const std::array<int, 3>& dirIsNeg) const
{
// invDir: ray direction(x,y,z), invDir=(1.0/x,1.0/y,1.0/z), use this because Multiply is faster that Division
// dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x>0),int(y>0),int(z>0)], use this to simplify your logic
// TODO test if ray bound intersects
float txmin = (pMin - ray.origin).x * invDir.x;
float txmax = (pMax - ray.origin).x * invDir.x;
float tymin = (pMin - ray.origin).y * invDir.y;
float tymax = (pMax - ray.origin).y * invDir.y;
float tzmin = (pMin - ray.origin).z * invDir.z;
float tzmax = (pMax - ray.origin).z * invDir.z;
//如果方向某個維度是負的 會和pMax對應維度上的面上先相交
if (dirIsNeg[0]) {
std::swap(txmin, txmax);
}
if (dirIsNeg[1]) {
std::swap(tymin, tymax);
}
if (dirIsNeg[2]) {
std::swap(tzmin, tzmax);
}
float tenter = std::max(std::max(txmin, tymin), tzmin);
float texit = std::min(std::min(txmax, tymax), tzmax);
return tenter < texit && texit >= 0;
}
三角形求交
和作業5中的實現類似
這裡需要返回交點的屬性 比如材質 法線等等：
inline Intersection Triangle::getIntersection(Ray ray)
{
Intersection inter;
//如果法線與光線異側 說明光線是從背側射入
if (dotProduct(ray.direction, normal) > 0)
return inter;
double u, v, t_tmp = 0;
// D E2 = S1
Vector3f pvec = crossProduct(ray.direction, e2);
// S1 . E1
double det = dotProduct(e1, pvec);
if (fabs(det) < EPSILON)
return inter;
double det_inv = 1. / det;
//S
Vector3f tvec = ray.origin - v0;
u = dotProduct(tvec, pvec) * det_inv;
if (u < 0 || u > 1)
return inter;
//S2
Vector3f qvec = crossProduct(tvec, e1);
v = dotProduct(ray.direction, qvec) * det_inv;
if (v < 0 || u + v > 1)
return inter;
t_tmp = dotProduct(e2, qvec) * det_inv;
// TODO find ray triangle intersection
if (t_tmp < 0) {
return inter;
}
inter.happened = true;
inter.coords = ray(t_tmp);
inter.normal = normal;
inter.distance = t_tmp;
inter.obj = this;
inter.m = this->m;
return inter;
}
SAH
SAH 全稱是Surface Area Heuristic 是在構建BVH過程中 進行劃分的一種方法 之前我們也提到了一種簡單的劃分方法 之前的方法在下圖所示的情況下:
會選擇圖中藍線這樣的劃分 最終boundingbox會有部分重疊 從而降低光線求交的效率
SAH提出了一個函式來量化地估計劃分割槽域的代價 以將一個boundingbox劃分為AB兩個為例：
其中t-trav是遍歷樹節點的代價 這裡固定為一個常數0.125 pbrt的解釋是我們將光線求交的代價設定為1 遍歷效率大概是求交的1/8 重要的是它們之間的相對關係 而不是數值
t-isect代表著求交的代價 具體設定為每個區域內圖元的數目
pa pb代表著光線穿過子節點的機率 使用A的boundingbox面積 與 總boundingbox面積的比值來確定
具體流程概括如下：
1.首先我們和之前的劃分演算法一樣 找出最長的維度
2.這裡我們定義一定數目的bucket 在整體的boundingbox上進行劃分 計算出每個物體在哪個bucket的索引
3.每個bucket都代表著一種劃分方法 我們遍歷計算各方法的代價 然後找出最小代價對應的bucket索引：
4.藉助最小索引將物體分為兩堆
程式碼：參考https://github.com/mmp/pbrt-v3/blob/master/src/accelerators/bvh.cpp
BVHBuildNode* BVHAccel::SAHrecursiveBuild(std::vector<Object*> objects)
{
BVHBuildNode* node = new BVHBuildNode();
// Compute bounds of all primitives in BVH node
Bounds3 bounds;
for (int i = 0; i < objects.size(); ++i)
bounds = Union(bounds, objects[i]->getBounds());
//假如物體數目過少 就直接使用原來的劃分方法 沒必要
if (objects.size() <= 2) {
node = recursiveBuild(objects);
return node;
}
else {
Bounds3 centroidBounds;
for (int i = 0; i < objects.size(); ++i)
centroidBounds =
Union(centroidBounds, objects[i]->getBounds().Centroid());
int dim = centroidBounds.maxExtent();
//計算索引
constexpr int nBuckets = 8;
BucketInfo buckets[nBuckets];
for (int i = 0; i < objects.size(); ++i) {
int b = nBuckets *
centroidBounds.Offset(
objects[i]->getBounds().Centroid())[dim];
if (b == nBuckets) b = nBuckets - 1;
buckets[b].count++;
buckets[b].bounds =
Union(buckets[b].bounds, objects[i]->getBounds());
}
//計算每個劃分的代價
float cost[nBuckets - 1];
for (int i = 0; i < nBuckets - 1; ++i) {
Bounds3 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] = 0.125f +
(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 = objects.size();
int mid = 0;
//如果包含的物體過多 或者縮小劃分代價 則進行劃分
if (objects.size() > maxPrimsInNode || minCost < leafCost) {
for (int i = 0; i < objects.size(); ++i) {
int b = nBuckets *
centroidBounds.Offset(
objects[i]->getBounds().Centroid())[dim];
if (b == nBuckets) b = nBuckets - 1;
// 滿足條件的元素放到陣列前半部分
if (b <= minCostSplitBucket) {
std::swap(objects[i], objects[mid]);
mid++;
}
}
}
auto beginning = objects.begin();
auto middling = objects.begin() + mid;
auto ending = objects.end();
auto leftshapes = std::vector<Object*>(beginning, middling);
auto rightshapes = std::vector<Object*>(middling, ending);
assert(objects.size() == (leftshapes.size() + rightshapes.size()));
node->left = SAHrecursiveBuild(leftshapes);
node->right = SAHrecursiveBuild(rightshapes);
node->bounds = Union(node->left->bounds, node->right->bounds);
}
return node;
}
結果展示
BVH:
BVH-SAH:
