2020-06-05计算机图形学 0次访问

Ray-tracing (4)—Acceleration

ps: Thumbnail 《月色真美》。

包围盒(Bounding Volumes)

对于复杂的物体，我们可以使用相对简单的物体把复杂的物体包起来，保证物体在简单物体的形状之内。(如上图所示复杂物体茶壶可以包含在简单物体中)
如果一个物体连包围盒都碰不到，那么更不可能碰到包围盒里面的物体。所以如果要测试是否击中物体，那么可以先测试是否计算了包围盒！

轴对齐包围盒(Axis-Aligned Bounding Box (AABB))

这是一种我们常用的包围盒，如上图所示分别沿着x,y,z的三个对面(pairs of slabs)来定一个包围盒。

光线和AABB求交

对于2维的情况:

上图中第一张图求得的是光线和垂直于x轴的两对面的进入和出去的时间。
第二张图求得的是光线和垂直于y轴的两对面的进入和出去的时间。
第三张图综合了前两张图，求得了最终光线进入包围盒和离开了包围盒的时间，光线进入时间等于x轴和y轴最晚的进入时间，光线出包围盒的时间等于x轴和y轴最早出去的时间。

对于3维情况:

光线进入包围盒只有当光线进入全部的三个对面(x,y,z)才算进入AABB。
只要光线离开任意一个对面就算离开了AABB。
因此可以得到计算公式t_enter = max{t_min},即所有对面中最晚进入的时间才是进入包围盒的时间;t_exit=min{t_max},即所有对面最早出去的时间就是光线离开包围盒的时间。

判断光线是否和AABB相交了

特殊情况1:如上图,当t_exit < 0 时，此时包围盒在光线源点后面，此时没有相交!

特殊情况2: 如上图，当t_exit >= 0 && t_enter < 0时，此时显然光线源点在包围盒内，这样显然是和包围盒相交的。

总结当t_enter < t_exit && t_exit >= 0 的时候就是有交点的情况。

交点的计算

由于我们使用的包围盒是沿着坐标轴的，所以计算交点十分的简单 (这里只是之前文章中推导过的光线和平面相交的特殊情况)，上图是垂直于x轴的对面，计算交点t的公式为 $$t = \frac{p\prime_x - o_x}{d_x}$$ 另外两对对面计算方式类似。在得到了t之后根据光线的表达式很容易就算出交点了。

Acceleration

下面将介绍的是如何利用上面的包围盒来加速光线求交——我们将使用加速结构。
目前已经诞生了很多种加速结构，从空间划分(Spatial Partitions)的角度来看有:

Oct-Tree：就是八叉树，每次将空间分为8个相等的部分，再递归的对子空间进行划分 ,当划分的子空间足够小或是空间中三角形面的数量很少的时候会停止划分。
KD-Tree：每次将每块空间划分为两部分，且划分依次沿着x-axis，y-axis，z-axis,依次递归划分，终止条件与八叉树类似。
BSP-Tree：KD-Tree类似，唯一不同的是划分不再沿着固定一轴，可以任意方向划分，缺点自然是划分的空间没有规则性，求交困难。

上面的几种方法我个人感觉KD-Tree比Oct-Tree划分更合理，比BSP-Tree更简单，是上面提到的这三种中较好的较平均的一个方法。(在这里没有赘述它们的具体的原理，这不是本文讨论的重点！)
但是即便如此它还是存在一些问题：

我们已经知道一个最小包围盒包含的所有三角形的信息是存在KD-Tree的叶子节点上的，因此我们需要求解AABB和哪些三角形求交。但问题是很难判定一个AABB和三角形的交集，即三角形和盒子的求交比较不好做。

一个物体可能会出现在不同的叶子节点中，会造成存储的冗余。

BVH

因此现在，更多的会使用BVH这种加速结构，它是通过物体来划分的(Object Partitions)。

BVH的步骤

首先找到一个包围盒。
递归的把包围盒中的物体划分成两部分。
每次要重新计算新划分出来的两部分物体的包围盒。
当划分出的两部分物体已经足够少的时候划分就结束了。同样的我们把物体的信息存储在叶子节点上，中间节点用作加速结构节点的判断。

从上面的步骤中可以看出，BVH每次把物体分成两堆，然后求它们的包围盒，区别于KD-Tree每次把包围盒分成两半再考虑包围盒和物体的相交情况。因此使用BVH不会涉及到包围盒和三角形求交的问题!
同时由于BVH每次是先分的三角形再对这堆三角形求解包围盒，因此每个三角形只会在一个包围盒中！因此之前的冗余的问题也解决了。
但是BVH也有不太好的地方:划分出的不同的包围盒之间可能会相交，但是其实只要划分策略不出本质性的错误，这一点问题不大。

划分两个子集的策略

每次沿着最长的一个轴 (x,y,z) 来进行划分，这样划分出的子集渐渐的会回到一个比较均匀的长度，即x,y,z三个方向的长度最终都会差不多。
经过上一步已经知道每次该往哪个轴劈一刀分成两个子集了，不过应该从哪里把它们分开呢? —— 找沿着划分轴的第 $\frac n2$ 个物体，然后劈开，这样保证左右两边的物体数量相同。 (这里的策略称之为 等量划分（SPLIT_EQUAL_COUNT）) 这里找第 $\frac n2$ 个物体可以用快速选择( 时间复杂度$O(n)$ )算法，来找出该物体。

注意: 场景中加入了新的物体，或者有物体动了，那么需要重新计算BVH，当然对于上面其他的加速结构也是这样，好像很耗时呢！

BVH Traversal

经过上面的一番说明，现在已经知道该如何去建立一个BVH了，但是该如何去计算和光线的求交呢？

首先我们要清楚BVH的数据结构:

其实BVH的数据结构就是一棵二叉树；
对于中间节点我们存储的是包围盒指向左右儿子结点的指针；
对于叶子节点我们存储的是包围盒和落在包围盒中的物体的信息。

光线和物体的求交

光线和物体求交伪代码

Intersect(Ray ray, BVH node) {
   if (ray misses node.bbox) return;

   if (node is a leaf node)
      test intersection with all objs;
      return closest intersection;
   hit1 = Intersect(ray, node.child1);
   hit2 = Intersect(ray, node.child2);
 
 
return the closer of hit1, hit2;
}

根据上述伪代码，可以知道这个过程是这样的:

如果光线没有和结点相交，则直接返回什么也不用计算;
如果光线和结点相交了，且相交的结点是叶子节点，那么遍历结点中所有的物体，返回光线最先打到的物体。
如果光线和结点相交，相交结点是中间节点，那么光线既可能和该节点的左孩子相交也可能和该节点的右孩子相交，因此要递归遍历左右孩子，然后返回光线最先打到的物体。

BVH核心代码

建立BVH >folded

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)
    //在把模型导入时，我们已经给每个三角形设置了一个包围盒，这样我们在之后求解大包围盒的时候可以很快的得到结果
    //union函数可以理解为合并包围盒
        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());
        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;
}

BVH和光线求交 >folded

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    // TODO Traverse the BVH to find intersection
    Intersection closestObj;
    std::array<int, 3> dirIsNeg = {(int)(ray.direction.x > 0),(int)(ray.direction.y > 0),(int)(ray.direction.z > 0)};
    if(node->bounds.IntersectP(ray,ray.direction_inv,dirIsNeg) == false) return closestObj;
    if(node->left == nullptr && node->right == nullptr){
        for(int i = 0;i < sizeof(node->object) / sizeof(Object);i++){
            Intersection inter = node->object[i].getIntersection(ray);
            if(inter.distance < closestObj.distance)
                closestObj = inter;
        }
        return closestObj;
    }

    Intersection hit1 = getIntersection(node->left,ray);
    Intersection hit2 = getIntersection(node->right,ray);

    return hit1.distance < hit2.distance ? hit1 : hit2;
}

SAH (Surface Area Heuristic)

在上述BVH的建树过程中，我们划分两个子集的策略是等量划分，当然除此之外的策略还有

中点划分（SPLIT_MIDDLE），即在先前选取的轴上找到最靠两侧的物体的包围盒，连接二者包围盒的重心得到一条线段，取其中点作为划分点，中点以左划分到左子树，中点以右划分到右子树。这个策略显然问题非常大。
表面积启发式算法划分（SPLIT_SAH, SAH=Surface Area Heuristic）。这是我们下面讨论的重点。

中点划分和等量划分存在的问题

某侧子树上可能会拥挤过多的物体，而另一侧子树上却太少。当然这个是中点划分才存在的问题。
包围盒之间互相“重叠”（overlapped）的情况。如果两棵子树对应的包围盒“重叠”的越多，那么一条射线穿过该区域时同时击中两子树的概率也就越大，这就意味着两棵子树都得进行相交测试。如果所有的子树都不得不进行相交测试，那便失去了“场景管理"本该具有的作用。这是两者都具有的问题。

而SAH可以较好的处理上述的两个问题，同时其性能较好，成本较低。该算法基于的理论是复杂度成本分析和概率论，是一个不错的策略!

SAH的思想是按某种方式进行划分之后可以使得计算光线相交测试等产生的开销最小。
首先最重要的开销就是光线穿过某区域与该区域内所有的图元做相交测试的花费，那么该如何计算呢?方法如下:
假设对于第i个物体，其做光线相交测试的时间复杂度为$t_{isect}(i)$代表光线和该图元进行相交计算所消耗的时间，那么总体的时间复杂度就为$$\sum_{i=1}^N t_{isect}(i)$$ $t_{isect}(i)$这里我们假设$t_{isect}(i)$对于所有的图元保持相同即假设光线和每个图元求交的时间一样，不妨简单的设置为1。因此此时的话可用某块区域中图元的数量代替花费。

SAH的总开销$$c(A,B) = t_{trav} + p_A\sum_{i=1}^{N_A} t_{isect}(a_i) + p_B\sum_{i=1}^{N_B} t_{isect}(b_i)$$
其中$t_{trav}$表示光线遍历内部结点和判断所穿越的子节点产生的开销。$p_A$和$p_B$表示光线穿越各个子节点的概率。$\sum_{i=1}^{N_A} t_{isect}(a_i)$和$p_B\sum_{i=1}^{N_B} t_{isect}(b_i)$表示将场景划分为A和B两个区域之后，射线分别和A和B两个区域求交所产生的花费。 SAH只是BVH的一种划分的方式，只不过在考虑从哪里划分开这点上有点讲究，确定划分的点之后每次依然还是把一块区域划分为两份

计算$p_A$和$p_B$
根据几何概率可计算$p_A$和$p_B$的概率值，该理念表明，针对包含另一个凸体B的凸体A，随机光线穿越B且穿越A的条件概率可以表示为表面积$S_A$和$S_B$的比值,即:$$p(A|B) = \frac{S_A}{S_B}$$
具体到如下图这种情况:

考虑一个父节点（C）下带有左（A）右（B）子节点（子树）这种情况。父节点C的表面积为S(C)，左子节点A的表面积为S(A)，右子节点B的表面积为S(B)。那么击中父节点下的左子节点A的概率为$$p(A|C) = \frac{S_A}{S_C}$$击中父节点下的右子节点B的概率为$$p(B|C) = \frac{S_B}{S_C}$$
注意$\frac{S_A}{S_C} + \frac{S_B}{S_C}$是有可能大于S(C)的，因为S(A)与S(B)可能会发生”重叠”,因此该值越小越好。

划分区域
在划分一个节点代表的空间区域时，可以通过不同的切法将该空间划分成两个子区域。切法可以平均等距地一刀一刀地切，也可以耍点小手段带点”智能“使其”自适应“。总之每次划分，都会得到两个子区域A和B。那么相应也会算出一个cost(A,B)来。比较所有划分方案下所得的cost(A,B)，取值最小的方案就是成本最低的划分方案，也作为划分该节点的最优方案。

因此基于上面的思想，我们可以这样做。首先把场景均匀的划分成k份，每一份我们称之为一个桶(如上图两条竖直的直线之间就是一个桶，注意桶的数量不要太多！)。然后依次以桶的边界作为划分的边界来进行划分，然后依次的求出以该边界划分之后产生的消耗，从中找出消耗最小的一条边界，它就是我们的划分方案，之后继续递归的进行划分。

伪代码

SAH伪代码

for each axis: x,y,z:
    /*初始化桶。沿着当前遍历的轴均匀的生成n个桶
    (相当于就是均匀的把坐标轴划分成n份)；*/
    initialize buckets;
    for each primitive p in node:
        /*计算图元该落到哪个桶内。只用看看图元中心的当前
        遍历轴的坐标在哪个桶的范围内即可。*/
        b = compute_bucket(p.centroid);
        /*计算桶内图元包围盒*/
        b.bbox.union(p.bbox);
        /*图元数量对应着光线求交的时间复杂度*/
        b.prim_count++;
    /*依次计算每种划分的花费*/
    for each of the B-1 possible partition evaluate SAH
/*按照产生最低花费的位置划分开或者如果达到数量已经比较少了
那么直接就可以装在叶子节点中不用再划分了*/
execute lowest cost partitioning found (or make node a leaf)

SAH核心代码

ps:本文中给出的核心代码中的英文注释有大部分是由我给出的(只是因为觉得在虚拟机上装输入法太麻烦了。。。)，因此注释的意思大概能够理解即可，不必深究。

使用SAH优化建树过程 >folded

//Buckets' num
const int nBuckets = 12;

// Build BVH WITH SAH
BVHBuildNode* BVHAccel::recursiveBuildWithSAH(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 {//patrition primitives using SAH
    	Bounds3 centroidBounds;
        for (int i = 0; i < objects.size(); ++i)
            centroidBounds =
                Union(centroidBounds, objects[i]->getBounds().Centroid());

        float cost[nBuckets];
        float minCost = std::numeric_limits<double>::max();
        int minCostSplit;
        //the pos to split primitives into two parts
        int splitPos = 0;
    	for(int dim = 0;dim < 3;dim++){//dim -> dimension
    		//Allocate BucketInfo for SAH partition buckets
	    	BucketInfo buckets[nBuckets];
	    	//Initialize BucketInfo for SAH partition buckets
	    	for(int i = 0;i < objects.size();i++){
	    		int b;
	    		switch(dim){
	    			case 0:
	    				b = nBuckets * centroidBounds.Offset(objects[i]->getBounds().Centroid()).x;
	    		 		break;
	    		 	case 1:
	    		 		b = nBuckets * centroidBounds.Offset(objects[i]->getBounds().Centroid()).y;
	    		 		break;
	    		 	case 2:
	    		 		b = nBuckets * centroidBounds.Offset(objects[i]->getBounds().Centroid()).z;
	    		}
	    		if(b == nBuckets) b = nBuckets - 1;
	    		buckets[b].count++;
	    		buckets[b].bounds = Union(buckets[b].bounds,objects[i]->getBounds());
	    	}

	    	//Compute costs for splitting after each bucket
    		for(int i = 0;i < nBuckets - 1;i++){
    			Bounds3 b0,b1;
    			int count0 = 0,count1 = 0;
    			//combined objects int 0~i buckets
    			for(int j = 0;j <= i;j++){
    				b0 = Union(b0,buckets[j].bounds);
    				count0 += buckets[j].count;
    			}
    			//combined objects int i+1~nBuckets-1 buckets;
    			//and now we split objects int 2 partions
    			for(int j = i + 1;j < nBuckets;j++){
    				b1 = Union(b1,buckets[j].bounds);
    				count1 += buckets[j].count;
    			}
    			//compute cost
    			//this algorithm include both "the average distribution" and "low overlapping"
    			//we just assume that ttrav = 1/8 and tisect = 1
    			cost[i] = 0.125f + (count0 * b0.SurfaceArea() + count1 * b1.SurfaceArea()) / centroidBounds.SurfaceArea();
    		}

    		//find bucket to split
    		for(int i = 1;i < nBuckets-1;i++){
    			if(cost[i] < minCost){
    				minCost = cost[i];
    				minCostSplit = i;
    				splitPos = 0;
    				for(int i = 0;i <= minCostSplit;i++){
			    		splitPos += buckets[i].count;
			    	}
    			}
    		}
    	}

    	auto beginning = objects.begin();
        auto middling = objects.begin() + splitPos;
        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;

}