您需要 登录 才可以下载或查看，没有账号？注册
x
从零开始学图形学：通俗讲解简单加速结构——k-d tree, SAH, BVH 文/启思
整个系列的文章归档整理于 归档整理目录 中。由于笔者能力有限，欢迎指正。
引入上一篇文章中，我们讲了一种最简单的光线追踪方法：whitted-style ray tracing。其方法是从摄像机发射射线，在反射、折射的材质上反弹，在每个反弹点上计算能量，回溯时加起来做着色。之后解决了一个技术问题：如何判断射线与三角形相交。 现在，在能力上，我们已经可以用 whitted 方法渲染图像了，就如上一篇文章的代码所示。但是，考虑一下这个算法的时间复杂度：如果每条射线需要枚举所有的三角形做判断，那复杂度高达三角形个数*像素个数*平均反弹次数。 而现在想建模一个小物体的表面，往往要几千甚至几万个三角形，一个商业级产品，屏幕内甚至可能同时存在几亿个三角形。
虚幻五宣传demo，同屏几亿的三角形建模
这样的话，复杂度就显得太高了。有没有什么更加高效的方法？那肯定有的。 考虑一下，我们的目标是找到射线穿过的最近的三角形。尽管屏幕中三角形相当多，忽略“最近”这个条件，我们每次至多只需要检测一条射线所穿过的所有三角形即可。除非所有的三角形都整整齐齐排列成一排，跟糖葫芦一样被射线穿过，否则我们根本没必要去检测全部三角形。 换句话说，屏幕中大部分三角形，完全可以利用其在空间中的位置忽略掉，从而大大减少运行时间。 具体使用的数据结构，称其为加速结构。这里介绍两类加速结构：基于空间的、基于物体的。不过不论哪种方法，都没有完美的万金油，所以也说不上哪种更好。 基于空间的加速结构这种方法把物体所在的空间划分成若干组，再利用数据结构把它们串起来加速查询。原本的物体则被分配到组内。查询时，如果该组内所有物体都不可能被射线穿过，那么忽略。 一个简单的例子是 k-d tree，假设在 k 维空间有若干物体，建立这么一棵树： - 建立根节点，代表整个空间。指定当前节点为根节点。
- 选择第 i 维， ，取 k-1 维超平面 ，满足 ，即一个垂直于该轴的平面。将当前节点的空间划分成两部分，分别形成左右儿子，当前节点包含的所有物体按照空间位置被划分到两个儿子内，被该超平面覆盖的物体保留在当前节点或者复制两份放到子节点。
- 分别递归两个儿子，重复第二步，直到当前节点物体数量很少为止。
下图是 wiki 上的图，展示了物体是一个点的情况
当物体是一个点时 k-d tree方法具有很大的启发性。第二步中选择的维度 i 和超平面的唯一参数 M 是可以通过一些启发式方法得到的。 一般来说， i 是循环遍历 k 维的每一维（三维情况就是 x-axis, y-axis, z-axis 循环）。
对于 M 的选取，学问很大。比如上图中点的情况，M 一般选择所有点按该维度排序后的中位数。如果物体具有形形色色的形状（例如三角形，圆形），那么 M 可能就不能这么选了（当然，一般情况下也不会太差，可以选择中点之类来代替图形进而确定 M）。 在遍历的时候，只需要从根节点开始，判断当前射线是否穿过当前节点的空间，如果不穿过（或者不可能成为更优的答案），那么直接回溯；如果穿过，检查该节点内包含的所有物体是否与射线相交。
来源games101。随便画一画射线，就会发现总能跳过一些节点。
对于 k-d tree 的优点，首先它很好实现，二叉树的组织方式效率也较高；其次树中每个节点都是一个平行于坐标轴的包围盒，判断射线是否跟包围盒相交很简单，计算量很少。 缺点就是，对于 M 的选择，很有可能很多物体都与该超平面相交，会造成很多麻烦。对于这些物体，不论你是都放到子节点，还是保留在该节点，都不会是最优的，甚至复制还会造成内存占用问题。 基于物体的加速结构基于空间的加速结构有一个问题，那就是在空间上做划分对物体本身不友好，例如一个物体被超平面穿过。 基于物体的方法，直接对物体做类别划分，并组合成树形结构。保证一个物体只属于一个类别，但是每类的空间占用有可能重合。 Bounding Volume Hierarchy (BVH) 方法就是其中一种。按包围盒进行树形结构组织。
图源games101课件
建立的过程看图就明白了。与 k-d tree 相同的是，BVH 在建树时也需要做划分，满足如果射线不穿过父节点，那么肯定不穿过子节点。但是划分不是在空间结构上，而是对物体的集合做划分。这也就导致空间会重合。 查询跟 k-d tree 一样，就不赘述了。 Surface Area Heuristic （SAH）实际上是一种基于空间的划分，把它单列出来是怕行文太长。 简述一下基本思路。首先对 u 的子树的 cost 做了一个建模：
其中 C_u 表示 u 节点的花费，p_x 表示光线追踪时生成的射线穿过 x 节点的概率。C_{trav} 是一些额外花费，例如计算与包围盒是否相交之类的。可以看做与 u 有关的常数。 假设有 A 个物体划分给 x 子节点， B 个物体划分给 y 子节点，并且认为概率 p 与包围盒大小成正比，那么有如下近似：
其中 C_{isect} 是判断射线与物体是否相交的花费。 于是我们得到了优化目标：划分空间时，让这个式子最小。 对于每一个轴，将其等距离划分成 B 份，这样得到 B-1 个划分平面。对于第 i 个平面，该平面把三角形的中点集合划分成两份，则就按这两份计算花费，选取花费最小的平面。 三个轴都算一下，选最小的轴。
实验
渲染一个 Stanford bunny，其包含 69451 个三角形。 - 实现判断包围盒是否相交
- 实现利用 BVH 加速结构查询最近相交点
实现 SAH
解析
包围盒是否相交，只需要判断射线 进入与离开包围盒的时间即可。 对于包围盒的六个平面，互相平行的三对，只需要每一对挨个计算。假设对垂直于 x 轴那一对计算，实际上就是解 的 t 。具体实现如下： 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
auto t0 = (pMin - ray.origin) * invDir;
auto t1 = (pMax - ray.origin) * invDir;
auto v0 = Vector3f(std::min(t0.x, t1.x), std::min(t0.y, t1.y), std::min(t0.z, t1.z));
auto v1 = Vector3f(std::max(t0.x, t1.x), std::max(t0.y, t1.y), std::max(t0.z, t1.z));
float tenter = std::max(v0.x, std::max(v0.y, v0.z));
float texit = std::min(v1.x, std::min(v1.y, v1.z));
return tenter < texit && texit >= 0;
}
对于 BVH 的查询过程，先判断是否跟当前包围盒相交，然后判断是否是叶子，如果是叶子的话遍历所有三角形，否则递归，返回最小值。 Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
// TODO Traverse the BVH to find intersection
std::array<int, 3> dirIsNeg = {ray.direction.x>0, ray.direction.y>0, ray.direction.z>0};
if(!node->bounds.IntersectP(ray, ray.direction_inv, dirIsNeg)) {
return Intersection();
} else {
if(node->left && node->right) {
auto L = getIntersection(node->left, ray);
auto R = getIntersection(node->right, ray);
if(L.distance > R.distance) return R;
return L;
} else {
return node->object->getIntersection(ray);
}
}
}
一个简单的查询优化是，计算与包围盒是否相交时顺便返回 t_enter，如果一个子节点的递归答案已经早于另一个子节点的 t_enter 了，那么就不再递归。实测了一下没啥用……引入的额外计算反而有点慢。
对于 SAH ，查询和上面一样，主要是划分过程的实现。这里的 B 随便指定的 50，可以随意改。 if(splitMethod == SplitMethod::SAH) {
struct Data {
std::vector<Object*> lobjs, robjs;
double C;
};
auto deal = [](std::vector<Object*> objects, int axis_id) {
auto getValue = [](Object* obj, int id) {
if(id == 0) return obj->getBounds().Centroid().x;
if(id == 1) return obj->getBounds().Centroid().y;
if(id == 2) return obj->getBounds().Centroid().z;
assert(0);
};
int B = 50;
double maxV = -std::numeric_limits<double>::infinity();
double minV = std::numeric_limits<double>::infinity();
for(auto obj : objects) {
maxV = std::max(maxV, (double)getValue(obj, axis_id));
minV = std::min(minV, (double)getValue(obj, axis_id));
}
std::vector<std::vector<Object*>> buckets(B);
std::vector<Bounds3> lbd(B), rbd(B);
for(auto obj : objects) {
auto val = getValue(obj, axis_id);
int i = std::min((int)((val-minV)/(maxV-minV) * B), B-1);
buckets.push_back(obj);
lbd = Union(lbd, obj->getBounds());
}
for(int i = 0;i < B;i ++) rbd = lbd;
for(int i = 1;i < B;i ++) lbd = Union(lbd, lbd[i-1]);
for(int i = B-2;i >= 0;i --) rbd = Union(rbd, rbd[i+1]);
Bounds3 bd = lbd[B-1];
int lsz = 0, rsz = objects.size(), sz = objects.size();
double C = std::numeric_limits<double>::infinity();
int partition = -1;
for(int i = 0;i < B - 1;i ++) {
lsz += buckets.size();
rsz -= buckets.size();
if(lsz && rsz) {
double nowC = 1.0*lbd.area()/bd.area()*lsz
+ 1.0*rbd[i+1].area()/bd.area()*rsz;
if(nowC < C) {
C = nowC;
partition = i;
}
}
}
assert(partition != -1);
Data ret;
ret.C = C;
for(auto obj : objects) {
auto val = getValue(obj, axis_id);
int i = std::min((int)((val-minV)/(maxV-minV) * B), B-1);
if(i <= partition) ret.lobjs.push_back(obj);
else ret.robjs.push_back(obj);
}
return ret;
};
Data x_result = deal(objects, 0);
Data y_result = deal(objects, 1);
Data z_result = deal(objects, 2);
if(x_result.C > y_result.C) std::swap(x_result, y_result);
if(x_result.C > z_result.C) std::swap(x_result, z_result);
//std::cout << x_result.C << std::endl;
node->left = recursiveBuild(x_result.lobjs);
node->right = recursiveBuild(x_result.robjs);
node->bounds = Union(node->left->bounds, node->right->bounds);
}
对比一下，SAH 是 6s，BVH 是 7s，差别不大。
| 
评分
