总览
在之前的编程练习中,我们实现了基础的光线追踪算法,具体而言是光线传输、光线与三角形求交。我们采用了这样的方法寻找光线与场景的交点:遍历场景中的所有物体,判断光线是否与它相交。在场景中的物体数量不大时,该做法可以取得良好的结果,但当物体数量增多、模型变得更加复杂,该做法将会变得非常低效。因此,我们需要加速结构来加速求交过程。在本次练习中,我们重点关注物体划分算法 Bounding Volume Hierarchy (BVH)。本练习要求你实现 Ray-Bounding Volume 求交 与 BVH 查找。
首先,你需要从上一次编程练习中引用以下函数:
- Render() in Renderer.cpp: 将你的光线生成过程粘贴到此处,并且按照新框
架更新相应调用的格式。 - Triangle::getIntersection in Triangle.hpp: 将你的光线-三角形相交函数粘贴到此处,并且按照新框架更新相应相交信息的格式。
在本次编程练习中,你需要实现以下函数:
- IntersectP(const Ray& ray, const Vector3f& invDir,const std::array<int, 3>& dirIsNeg) in the Bounds3.hpp: 这个函数的作用是判断包围盒 BoundingBox 与光线是否相交,你需要按照课程介绍的算法实现求交过程。
- getIntersection(BVHBuildNode* node, const Ray ray) in BVH.cpp: 建立 BVH 之后,我们可以用它加速求交过程。该过程递归进行,你将在其中调用你实现的
Bounds3::IntersectP
.
开始实现
编译运行
基础代码只依赖于 CMake,下载基础代码后,执行下列命令,就可以编译这个项目:
mkdir build
cd ./build
cmake ..
make
在此之后,你就可以通过 ./Raytracing
来执行程序。
代码框架
我们修改了代码框架中的如下内容:
- Material.hpp: 我们从将材质参数拆分到了一个单独的类中,现在每个物体实
例都可以拥有自己的材质。 - Intersection.hpp: 这个数据结构包含了相交相关的信息。
- Ray.hpp: 光线类,包含一条光的源头、方向、传递时间
t
和范围range
. - Bounds3.hpp: 包围盒类,每个包围盒可由
pMin
和pMax
两点描述(请思考为什么)。Bounds3::Union
函数的作用是将两个包围盒并成更大的包围盒。与材质一样,场景中的每个物体实例都有自己的包围盒.- BVH.hpp: BVH 加速类。场景
scene
拥有一个BVHAccel
实例。从根节点开始,我们可以递归地从物体列表构造场景的 BVH.
代码框架解析
作业6的代码框架与作业5大体一致,主要有几处修改的地方。
一是将光线的源头、方向、传递时间 t
和范围 range
封装成了 Ray
光线类。因此,与光线有关的代码,需要有所修改。
struct Ray{
//Destination = origin + t*direction
Vector3f origin;
Vector3f direction, direction_inv;
double t;//transportation time,
double t_min, t_max;
Ray(const Vector3f& ori, const Vector3f& dir, const double _t = 0.0): origin(ori), direction(dir),t(_t) {
direction_inv = Vector3f(1./direction.x, 1./direction.y, 1./direction.z);
t_min = 0.0;
t_max = std::numeric_limits<double>::max();
}
Vector3f operator()(double t) const{return origin+direction*t;}
friend std::ostream &operator<<(std::ostream& os, const Ray& r){
os<<"[origin:="<<r.origin<<", direction="<<r.direction<<", time="<< r.t<<"]\n";
return os;
}
};
二是交点也被封装成了 Intersection
结构体,其主要包括,是否相交happened
、交点坐标coords
、交点法向量normal
、交点与相机距离 distance
,相交物体 obj
,材质 m
struct Intersection
{
Intersection(){
happened=false;
coords=Vector3f();
normal=Vector3f();
distance= std::numeric_limits<double>::max();
obj =nullptr;
m=nullptr;
}
bool happened;
Vector3f coords;
Vector3f normal;
double distance;
Object* obj;
Material* m;
};
#en
三是因为引入了BVH,所以在判断光线与物体是否有交点前,先对光线与包围框是否有交点进行判断。这里需要注意的是,空间中只有一个物体,而物体是由多个三角形网格组成。所以这里其实要构建两个 BVH,一个是针对空间中物体的 BVH,另一个是针对单一物体的组成网格的 BVH。运行代码,也可以看出,执行了两次 BVH 的构建。
那么在之后判断光线与物体中网格相交实际要执行三个步骤,第一步是判断光线是否与该物体包围框相交,如果相交,第二步判断是否与该物体中的三角形相交,而第二步由可以拆分成两个步骤,一是判断是否与三角形所在包围框相交,如果相交,二就是判断是否与三角形相交。
下面主要说明一下 BVH 的构建过程:
- 首先对于所有物体寻找一个包围框
- 递归地将物体的集合划分进两个子集中,然后重新计算子集的包围框
- 在合适的时候停止,比如一个包围框的物体个数为5个(代码里为1),在叶子节点存储物体信息
那么如何划分节点,一般地做法是选择一个维度去划分,一般是选最长的轴进行划分。在物体的中间位置划分节点 —— 每个子集有相同数量的物体,其对应的代码部分如下
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);
作业代码
Render
需要注意的地方在于,eyepos
的位置在
(
−
1
,
5
,
10
)
(-1, 5, 10)
(−1,5,10)。所以其实屏幕的中心坐标就不在
(
0
,
0
,
−
1
)
(0,0,-1)
(0,0,−1) 了,而是位于
(
−
1
,
5
,
9
)
(-1,5,9)
(−1,5,9)。所以其实只是相机源点发生了改变,而射线方向其实与作业5求法相同。
Vector3f dir = normalize(Vector3f(x, y, -1));
Ray ray(eye_pos, dir);
framebuffer[m++] = scene.castRay(ray, 0);
Triangle::getIntersection
因为将光线封装为 Ray
类,所以相应光线-三角形相交函数也需要有所修改:
inline Intersection Triangle::getIntersection(Ray ray)
{
Intersection inter;
if (dotProduct(ray.direction, normal) > 0)
return inter;
double u, v, t_tmp = 0;
Vector3f pvec = crossProduct(ray.direction, e2);
double det = dotProduct(e1, pvec);
if (fabs(det) < EPSILON)
return inter;
double det_inv = 1. / det;
Vector3f tvec = ray.origin - v0;
u = dotProduct(tvec, pvec) * det_inv;
if (u < 0 || u > 1)
return inter;
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
inter.happened = true;
inter.coords = ray.origin + t_tmp * ray.direction;
inter.normal = this->normal;
inter.distance = dotProduct(t_tmp * ray.direction, t_tmp * ray.direction);
inter.obj = this;
inter.m = this->m;
return inter;
}
IntersectP
这个函数的作用是判断包围盒 BoundingBox 与光线是否相交。其原理为:求光线与轴对齐包围盒交点 —— 光线进入所有平面内才算进入包围盒内,光线只要出了一个平面就算出了包围盒,所以
t
e
n
t
e
r
=
m
a
x
(
t
m
i
n
)
t_{enter}=max(t_{min})
tenter=max(tmin),
t
e
x
i
t
=
m
i
n
(
t
m
a
x
)
t_{exit}=min(t_{max})
texit=min(tmax) ,只要
t
e
n
t
e
r
<
t
e
x
i
t
&
&
t
e
x
i
t
>
=
0
t_{enter}<t_{exit} \&\& t_{exit}>=0
tenter<texit&&texit>=0 要求
t
e
x
i
t
>
=
0
t_{exit}>=0
texit>=0 是因为
t
e
x
i
t
<
0
t_{exit}<0
texit<0 说明包围盒在光线后面,没有交点
又因为包围盒是轴对齐的,所以其计算从向量计算简化为数值计算
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 t_enter;
float t_exit;
Vector3f t_enter_v3f = (pMin - ray.origin) * invDir;
Vector3f t_exit_v3f = (pMax - ray.origin) * invDir;
// for(int i=0; i<3; ++i) {
// if(!dirIsNeg[i])
// std::swap(t_enter_v3f[i], t_exit_v3f[i]);
// }
if(!dirIsNeg[0])
std::swap(t_enter_v3f.x, t_exit_v3f.x);
if(!dirIsNeg[1])
std::swap(t_enter_v3f.y, t_exit_v3f.y);
if(!dirIsNeg[2])
std::swap(t_enter_v3f.z, t_exit_v3f.z);
t_enter = std::max(t_enter_v3f.x, std::max(t_enter_v3f.y, t_enter_v3f.z));
t_exit = std::min(t_exit_v3f.x, std::min(t_exit_v3f.y, t_exit_v3f.z));
if (t_enter < t_exit && t_exit >= 0)
return true;
else
return false;
}
这里需要注意的是,如果一个轴光线的方向为负数,那么其离包围盒远端距离越小,所以这与光线方向为正,正好相反,所以需要 swap
一下
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 t_enter;
float t_exit;
Vector3f t_enter_v3f = (pMin - ray.origin) * invDir;
Vector3f t_exit_v3f = (pMax - ray.origin) * invDir;
// for(int i=0; i<3; ++i) {
// if(!dirIsNeg[i])
// std::swap(t_enter_v3f[i], t_exit_v3f[i]);
// }
if(!dirIsNeg[0])
std::swap(t_enter_v3f.x, t_exit_v3f.x);
if(!dirIsNeg[1])
std::swap(t_enter_v3f.y, t_exit_v3f.y);
if(!dirIsNeg[2])
std::swap(t_enter_v3f.z, t_exit_v3f.z);
t_enter = std::max(t_enter_v3f.x, std::max(t_enter_v3f.y, t_enter_v3f.z));
t_exit = std::min(t_exit_v3f.x, std::min(t_exit_v3f.y, t_exit_v3f.z));
if (t_enter < t_exit && t_exit >= 0)
return true;
else
return false;
}
getIntersection
其算法伪代码如下:
其思想就是如果光线与该节点包围盒没有交点,那么必定与包围盒内的物体没有交点。
如果有交点,如果当前节点是叶节点,那么就与包围盒内的物体判断是否有交点。如果不是叶节点,即中间节点,那么就递归其左右子树,取一个最近的交点返回。
Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
Intersection isect;
// TODO Traverse the BVH to find intersection
if(!node->bounds.IntersectP(ray, ray.direction_inv, std::array<int, 3> {ray.direction.x>0, ray.direction.y>0, ray.direction.z>0}))
return isect;
if(node->object != nullptr) // leaf node
return node->object->getIntersection(ray);
// internal node
Intersection isect_left, isect_right;
isect_left = getIntersection(node->left, ray);
isect_right = getIntersection(node->right, ray);
return isect_left.distance <= isect_right.distance ? isect_left : isect_right;
}
结果如下图所示
进阶代码
进阶代码要求实现 SAH(Surface Area Heuristic),其具体原理可以看这篇文章
简单来说如果图元分布不均匀,BVH 得到的划分结果会很差,比如下面的左图包围框有很大的重叠面积,那么射线与这两个包围框相交的概率肯定会比右图的大,计算就会增加。
SAH 就是考虑求交代价与包围盒面积,给出一个总的代价公式,我们寻找划分位置使得这个总代价最小。
这里的
S
(
A
)
S(A)
S(A) 表示左子树包围框面积,
S
(
B
)
S(B)
S(B) 表示右子树包围框面积,
S
(
C
)
S(C)
S(C) 表示当前节点包围框面积,
a
a
a表示左子树物体个数(认为求交代价为1),
b
b
b表示右子树物体个数,
0.125
0.125
0.125 为遍历划分代价。
其代码如下
BVHBuildNode* BVHAccel::recursiveBuild_SAH(std::vector<Object*> objects)
{
BVHBuildNode* node = new BVHBuildNode();
// std::cout<<objects.size()<<std::endl;
// Compute bounds of all primitives in SAH node
Bounds3 bounds;
for (int i = 0; i < objects.size(); ++i)
bounds = Union(bounds, objects[i]->getBounds());
if (objects.size() == 1) {
// Create leaf _SAHBuildNode_
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_SAH(std::vector{objects[0]});
node->right = recursiveBuild_SAH(std::vector{objects[1]});
node->bounds = Union(node->left->bounds, node->right->bounds);
return node;
}
else {
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;
}
int B_size = 10;
int minCostIdx = 1;
float minCost = std::numeric_limits<float>::infinity();
float SC = bounds.SurfaceArea();
// float SC = centroidBounds.SurfaceArea();
auto beginning = objects.begin();
auto ending = objects.end();
for(int i=1; i<B_size; ++i) {// search minCost and minCostIdx
auto middling = objects.begin() + std::max(1, (int)(objects.size() * i / B_size));
auto leftshapes = std::vector<Object*>(beginning, middling);
auto rightshapes = std::vector<Object*>(middling, ending);
assert(objects.size() == (leftshapes.size() + rightshapes.size()));
Bounds3 leftbounds, rightbounds;
for(int j=0; j<leftshapes.size(); ++j) {
leftbounds = Union(leftbounds, leftshapes[j]->getBounds());
// leftbounds = Union(leftbounds, leftshapes[j]->getBounds().Centroid());
}
for(int j=0; j<rightshapes.size(); ++j) {
rightbounds = Union(rightbounds, rightshapes[j]->getBounds());
// rightbounds = Union(rightbounds, rightshapes[j]->getBounds().Centroid());
}
float SA = leftbounds.SurfaceArea();
float SB = rightbounds.SurfaceArea();
float c_A_B = leftshapes.size() * SA / SC + rightshapes.size() * SB / SC + 0.125f;
if(c_A_B < minCost) {
minCost = c_A_B;
minCostIdx = i;
}
}
auto middling_minCost = objects.begin() + std::max(1, ((int)objects.size() * minCostIdx / B_size));
auto leftshapes_minCost = std::vector<Object*>(beginning, middling_minCost);
auto rightshapes_minCost = std::vector<Object*>(middling_minCost, ending);
assert(objects.size() == (leftshapes_minCost.size() + rightshapes_minCost.size()));
// std::cout<<leftshapes_minCost.size()<<" "<<rightshapes_minCost.size()<<std::endl;
node->left = recursiveBuild_SAH(leftshapes_minCost);
node->right = recursiveBuild_SAH(rightshapes_minCost);
node->bounds = Union(node->left->bounds, node->right->bounds);
}
return node;
}
这里考虑九个划分位置,即左边 1/10个物体,右边 9/10,左边 2/10,右边 8/10 …,然后面积我这里用的是整个包围框的面积,也有用质心包围框的面积,我感觉都可以,下面比较下两种方法的时间,可见 SAH 比 BVH 快了 1秒。