任务
首先,你需要从上一次编程练习中引用以下函数:
实现
本次作业需要自己实现的两个函数比较简单,难得在BVH树的建立的理解和SAH。
Render
void Renderer::Render(const Scene& scene)
{
std::vector<Vector3f> framebuffer(scene.width * scene.height);
float scale = tan(deg2rad(scene.fov * 0.5));
float imageAspectRatio = scene.width / (float)scene.height;
Vector3f eye_pos(-1, 5, 10);
int m = 0;
for (uint32_t j = 0; j < scene.height; ++j) {
for (uint32_t i = 0; i < scene.width; ++i) {
// generate primary ray direction
float x = (2 * (i + 0.5) / (float)scene.width - 1) *
imageAspectRatio * scale;
float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;
Vector3f dir = normalize(Vector3f(x, y, -1));
Ray ray = Ray(eye_pos,dir);
framebuffer[m++] = scene.castRay(ray, 0);
}
UpdateProgress(j / (float)scene.height);
}
UpdateProgress(1.f);
// save framebuffer to file
FILE* fp = fopen("binary.ppm", "wb");
(void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);
for (auto i = 0; i < scene.height * scene.width; ++i) {
static unsigned char color[3];
color[0] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].x));
color[1] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].y));
color[2] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].z));
fwrite(color, 1, 3, fp);
}
fclose(fp);
}
Triangle::getIntersection
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
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 = m;
return inter;
}
IntersectP
实现原理是课上讲的方法,对每一个轴都计算一次tmin和tmax,tmin是光线进入包围盒某一个轴的时间,tmax是光线离开包围盒某一个轴的时间,计算完三个轴的时间后,对t_enter是max{tmin},t_exit是min(tmax),只有当t_enter<t_exit并且t_exit >=0时才算光线经过包围盒。
这里有关dirIsNeg的数组,用来判断光线方向的。有什么作用呢?在这里使用这个数组提前判断方向从而可以简化计算。以x轴为例,图里的的光源是在左边,x0<x1,那么假设光源在右边,那么只对x的值求出来是t_enter > t_exit的,明显有问题,设置了dirIsNeg数组后,在x轴上是1,说明光线的x分量是负方向的,因此交换t_enter和t_exit才能保证t_enter<t_exit。为什么能这么交换呢?因为这里只要求线是否经过包围盒。
inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
const std::array<int, 3>& dirIsNeg) const
{
double x_min = (pMin[0] - ray.origin[0])*invDir[0];
double x_max = (pMax[0] - ray.origin[0])*invDir[0];
double y_min = (pMin[1] - ray.origin[1])*invDir[1];
double y_max = (pMax[1] - ray.origin[1])*invDir[1];
double z_min = (pMin[2] - ray.origin[2])*invDir[2];
double z_max = (pMax[2] - ray.origin[2])*invDir[2];
if(!dirIsNeg[0]){
double temp = x_min;
x_min = x_max;
x_max = temp;
}
if(!dirIsNeg[1]){
double temp = y_min;
y_min = y_max;
y_max = temp;
}
if(!dirIsNeg[2]){
double temp = z_min;
z_min = z_max;
z_max = temp;
}
double t_min = std::fmax(x_min, std::fmax(y_min,z_min) );
double t_max = std::fmin(x_max, std::fmin(y_max,z_max) );
if(t_min < t_max && t_max >= 0)return true;
return false;
}
getIntersection
这里的实现更加简单,就是简单的二叉树的遍历,如果光线与包围盒有交点,那么就寻找这个包围盒的子包围盒,直到找到距离最近且有交点的包围盒。然后在包围盒里遍历三角形是其他的函数实现的。
Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
std::array<int, 3> dirIsNeg;
dirIsNeg[0] = (ray.direction[0]>0);
dirIsNeg[1] = (ray.direction[1]>0);
dirIsNeg[2] = (ray.direction[2]>0);
Intersection inter;
if( node->bounds.IntersectP(ray,ray.direction_inv,dirIsNeg) == false )return inter;
if( node->left ==nullptr && node->right == nullptr ){
return node->object->getIntersection(ray);
}
Intersection l = getIntersection( node->left,ray );
Intersection r = getIntersection( node->right,ray);
if(l.distance < r.distance){
return l;
}else {
return r;
}
return inter;
}
结果
提高
本作业框架的BVH树的建立
先遍历传入的全部对象,创建一个包住全部对象的包围盒,找长度最大的轴,找到后将整个空间一分为二 ,包围盒的中心点为依据,按最长轴左右排序,将一个整体包围盒B0按对象的数量,均分为B1和B2,然后对B1和B2也进行这样的操作,直到所有的单个对象都处于叶子节点。
//递归构建BVH树
BVHBuildNode* BVHAccel::recursiveBuild(std::vector<Object*> objects)
{
BVHBuildNode* node = new BVHBuildNode();
//先遍历传入的全部的对象,创建一个最小的嫩能包住全部对象的包围盒
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 {
//如果有三个以上的对象,那么需要进行递归创建,先遍历全部对象创建一个把全部对象包住的包围盒
Bounds3 centroidBounds;
for (int i = 0; i < objects.size(); ++i)
centroidBounds =
Union(centroidBounds, objects[i]->getBounds().Centroid());
//查询最长轴是哪个轴,0是x,1是y,2是z
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()));
//根据划分出来的两个新的对象数组,递归创建BVG树
node->left = recursiveBuild(leftshapes);
node->right = recursiveBuild(rightshapes);
node->bounds = Union(node->left->bounds, node->right->bounds);
}
return node;
}
SAH方法建立BVH树
SAH的步骤
-
分区选择:选择一个维度(如x轴、y轴或z轴)作为分区轴,并将该轴上的范围划分为多个候选分区点。
-
成本计算:对于每个候选分区点,计算将其作为分区边界时的SAH成本。这涉及估算射线穿过每个子区域的可能性以及子区域内图元的相交测试成本。
-
分区决策:选择SAH成本最低的分区点作为最终的分区边界,并将图元分配到相应的子区域中。
-
递归构建:对每个子区域重复上述步骤,直到达到所有的图元都在BVH树的叶子上
这里使用简化的方式,认为求交的时间为1,因为是比较大小,Ctrav是固定开销,不影响大小的比较,因此直接当成0,那么总时间可以求出来是T = Pa*Na + Pb*Nb,概率通过面积求比值得出,该作业框架里包围盒的面积可以直接调用Bounds3::SurfaceArea()方法求得。
// 递归构建BVH树
BVHBuildNode *BVHAccel::recursiveBuild(std::vector<Object *> objects)
{
BVHBuildNode *node = new BVHBuildNode();
// 先遍历传入的全部的对象,创建一个最小的嫩能包住全部对象的包围盒
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
{
// 如果有三个以上的对象,那么需要进行递归创建,先遍历全部对象创建一个把全部对象包住的包围盒
Bounds3 centroidBounds;
for (int i = 0; i < objects.size(); ++i)
centroidBounds =
Union(centroidBounds, objects[i]->getBounds().Centroid());
// 查询最长轴是哪个轴,0是x,1是y,2是z
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;
}
//使用SAH方法划分BVH树
float min_time = std::numeric_limits<float>::max();
int index = 0;
for (int i = 0; i < objects.size(); i++)
{
auto middling = objects.begin() + i;
auto leftobjects = std::vector<Object *>(objects.begin(), middling);
auto rightobjects = std::vector<Object *>(middling, objects.end());
float time_left, time_right, time;
auto left_bound = Bounds3();
auto right_bound = Bounds3();
for (int j = 0; j < leftobjects.size(); j++)
{
left_bound = Union(left_bound, leftobjects[j]->getBounds());
}
for (int k = 0; k < rightobjects.size(); k++)
{
right_bound = Union(right_bound, rightobjects[k]->getBounds());
}
auto big_bound = Union(left_bound, right_bound);
float Sbig = big_bound.SurfaceArea();
float Sleft = left_bound.SurfaceArea();
float Sright = right_bound.SurfaceArea();
time_left = Sleft / Sbig * leftobjects.size();
time_right = Sright / Sbig * rightobjects.size();
time = time_left + time_right;
if (time < min_time)
{
min_time = time;
index = i;
}
}
auto beginning = objects.begin();
auto middling = objects.begin() + index;
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;
}
开销的变化
该图是使用原来作业框架自带的方法构建BVH树,因为虚拟机性能调的比较高,因此0秒
该图是使用SAH方法划分BVH树,可以看到开销大了很多
该图是使用默认划分的BVH树加速光线追踪所消耗的时间
该图是使用SAH方法划分BVH树加速光线追踪所消耗的时间
可以很明显的看到,使用SAH方法构建树的时间开销很大,但是使用其划分出来的BVH树加速光线追踪会比默认的快,但是由于作业中给的场景不复杂,因此光线追踪所消耗的时间变化不太能感觉出来。