又是好久没更新博客了，最近在写算法分析与设计课程的期末作业，作业的题目随意，我就随兴写了烟花粒子的四叉树可视化程序和光追渲染器的八叉树求交优化。之前写的光追渲染器对每个三角网格模型的求交都是暴力遍历所有的三角形，对于三角形数量很多的模型来说效率非常低，所以我捡起了这个渲染器并为每个三角网格模型构建一颗八叉树加快射线与三角形的求交速度。还真别说，性能提升巨大。所以这篇博客本质上是一个期末作业。最后，新年快乐！

  分治算法因其巧妙的算法思想能够将算法的时间复杂度降低到一个非常低的程度，例如广为流传的二分搜索算法将线性搜索时间复杂度$O(n)$降低到了$O(log\ n)$，这是一个极其恐怖的性能提升。诸多的实际应用领域例如计算机图形学、计算机视觉都有着分治算法的身影。在这里我们将重点关注计算机图形学领域非常实用、好用、高效的分治算法——四叉树空间分割算法和八叉树空间分割算法，二分搜索算法本质上处理的是一维的数据，但在图形学领域我们通常面临的是二维或者三维的点和向量，而且在实际应用中这些数据的量都非常庞大（几百万个点的点云、高精度的网格模型等等），因此为了加速这些更高维数据的搜索，一些高效的数据结构算法被提出，其中四叉树和八叉树算法就是其中之一。（当然也有高维的二叉树，例如k-d树）

一、背景介绍

  我们先来看一下算法的应用背景，计算机图形学领域研究的主要内容就是关于二维三维空间的图形绘制、物理模拟、几何建模、电脑动画等等可视化课题，目前最为流行的模型表示方法就是采用三角形网格模型，即每个网格面采用一个三角形来表示，如下图1所示：

图1 三角网格模型
  这种表示方法就是采用了有限元的思想，利用很多个三角形去逼近一个曲面，三角形越多，网格精度越高，模型越接近目标物体。在一个大型的场景中，通常有很多个这样的物体模型，因此采用暴力的方法遍历每一个物体送入绘制管线、碰撞检测等等是非常不可行的，特别是对于游戏、虚拟现实等等实时性要求比较高的应用来说这是灾难性的做法，暴力搜索策略的碰撞检测将耗费大部分的时间。除了碰撞检测之外，还有射线与物体的求交运算，即发射一条射线，找到射线与物体相交的一点，这个在射击游戏中非常常见，此外在基于光线追踪的图形渲染技术中亦是如此。暴力、简单的做法就是将射线方程与场景中的每个三角形进行数学上的交点求取运算（其实就是解一个方程），遍历完所有的三角形，取一个最近的交点就是最终结果，算法复杂度为$O(n)$。这种方法不可取，复杂场景中的三角形数量几千万甚至上亿，在每一帧进行这样求交运算带来的后果就是帧率的急剧降低，用户看到的画面将非常卡顿，毫无流畅的游戏体验感。   事实上，可以很容易理解，射线并不会与场景中所有的三角形相交，暴力遍历的算法有$99.9\%$的计算量都在做无用功，因为在所有的三角形中只有一个三角形才是我们要寻找的会相交的、距离最最近的那个三角形。注意到这些，学者们提出了一些基于分治的空间分割算法思想，将二维空间、三维空间做一个划分，排除点那些不可能相交的三角形，加速整个搜索过程。在众多的空间分割算法中，四叉树和八叉树的算法是其中的一种简单、高效、好用的空间分割算法，其中四叉树对应的是二维的空间划分，而八叉树对应的是三维的空间划分，算法的思想并不难理解，本质上属于二分搜索的高维扩展。 四叉树和八叉树就是$2D$和$3D$的“二分法”，搜索过程与二叉树搜索也类似，二叉树中是将数组排序后存入二叉树中，从而在查找中实现时间复杂度为$log\ n$。而四叉树/八叉树是按平面/空间范围划分有序节点，将所有点/面片/网格模型放入所属节点中，达到类似于排序的结果，进而在搜索时可以快速排除掉那些不符合条件的 点/面片/网格模型。 ## 二、四叉树分割算法   四叉树或四元树也被称为Q树（QuadTree）。四叉树广泛应用于图像处理、空间数据索引、2D中的快速碰撞检测、存储稀疏数据等，而八叉树（Octree）主要应用于3D图形处理。 #### 1、四叉树分割算法——原理   四叉树索引的基本思想是将地理空间递归划分为不同层次的树结构。它将已知范围的空间等分成四个相等的子空间，如此递归下去，直至树的层次达到一定深度或者满足某种要求（例如数据对象数量少于一定的阈值）后停止分割。四叉树的结构比较简单，并且当空间数据对象分布比较均匀时，具有比较高的空间数据插入和查询效率，因此四叉树是GIS中常用的空间索引之一。常规四叉树的结构如图所示。
图2 一颗构建好的四叉树
  这里划分空间通常是一个轴向包围盒，这个包围盒可以用它的最低点和最高点来表示。总的来说，四叉树的定义是：它的每个节点下至多可以有四个子节点，通常把一部分二维空间细分为四个象限或区域并把该区域里的相关信息存入到四叉树节点中。 四叉树的每一个节点代表一个矩形区域，每一个矩形区域又可划分为四个小矩形区域，这四个小矩形区域作为四个子节点所代表的矩形区域。   这里的四叉树只有叶子节点才存储数据对象（如点、三角面片、网格模型等等），内部节点不存储数据对象，因而访问数据对象都要根据内部节点走到叶子节点去访问。一般点是没有大小的，因此数据对象是点的时候没有必要考虑跨越了多个区域的情况。而如果数据对象是面片、网格模型等有大小的时，四叉树的构建就需要小心一点。如图2所示，数据对象是一个矩形，每个矩形有一定的大小，图中的$2$、$4$、$8$这三个数据对象跨越了两个区域，为了防止遗漏，在构建四叉树的时候最好把这些跨越了多个区域的数据对象均放入它所涉及到的叶子节点上。一般情况下叶子节点不会直接存储原始的数据对象，而是将原始的数据对象用一个线性表存储，然后四叉树中的叶子节点存储的是数据对象在线性表中的索引，这是为了防止程序代码耦合度过高。   说了这么多，接下来我们以构建烟花粒子系统的四叉树为例展开相关的算法介绍和实现。烟花粒子系统中的数据对象是粒子点云，即点集，我们为这个系统构建一颗四叉树，并将其可视化出来。首先定义四叉树中的一个节点对象，代码如下所示。一个节点对象可能是内部节点，也可能是外部节点。对于内部节点，它应该有四个子节点的指针；对于叶子节点，它应该有一个存储数据对象的表。对于每一个节点，都应该有一个包围盒指定当前节点所覆盖的矩形空间范围，我们用矩形空间范围的最小点和最大点来表示。 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
struct TreeNode
{
	//! 子节点
	TreeNode *children[4];
	//! 包围盒
	glm::vec2 bMin, bMax;
	//! 叶节点的物体列表
	std::vector<glm::vec2> objects;
	//! 是否是叶节点
	bool isLeaf;

	TreeNode() :
		isLeaf(false), bMin(glm::vec2(0.0f)), bMax(glm::vec2(0.0f))
	{
		children[0] = children[1] = children[2] = children[3] = nullptr;
	}

	TreeNode(glm::vec2 min, glm::vec2 max) :
		isLeaf(false), bMin(min), bMax(max)
	{
		children[0] = children[1] = children[2] = children[3] = nullptr;
	}
};
#### 2、四叉树分割算法——构建   从一颗空树开始，给定粒子系统中所有的粒子点云的线性表，我们采用自顶向下的方式构建一颗四叉树。从根节点往下划分，直到到达给定的树的深度或者当前节点覆盖的数据对象少于一定的数量时不再往下划分。根据数据对象数量、当前节点的深度，自顶向下构建四叉树的过程大致可以分成三类： - 如果当前数据对象数量为零，代表当前节点覆盖的区域不包含任何数据对象，则返回空指针，表示当前节点是一个不包含任何对象的空指针节点； - 如果当前数据对象不为零且小于一定的数量，亦或者当前节点的深度达到最大深度，则不再往下划分，将当前的节点构建为叶子节点，并将数据对象添加到叶子节点的存储表中； - 出去上面的两种情况，剩下的情况就是当前节点为内部节点的情况，对于每一个内部节点，需要往下划分四个区域节点，我们根据当前区域的范围分割成上、下、左、右四个子区域，根据数据对象的位置将数据对象表分成四类，递归调用构建函数，返回子节点的指针，最后再返回当前内部节点的指针。   以上的步骤递归嵌套地执行，即可完成自顶向下的四叉树构建。详细的代码实现如下所示，这个函数的输入是当前节点的深度，当前节点的矩形范围以及数据对象列表，返回构建的节点指针。 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
TreeNode * QuadTree::recursiveBuild(unsigned int depth, glm::vec2 min, glm::vec2 max,
	const std::vector<glm::vec2>& objects)
{
	//! if there is no object at all, just return nullptr.
	if (objects.empty())
		return nullptr;

	//! if the number of objects is less than 10 or reach the maxDepth,
	//! just create the node as leaf and return it.
	if (objects.size() < 4 || depth == mMaxDepth)
	{
		TreeNode *cur = new TreeNode(min, max);
		for (auto &point : objects)
		{
			if (isContain(point, min, max))
				cur->objects.push_back(point);
		}
		cur->isLeaf = true;
		return cur;
	}

	//! otherwise just subdivied into four sub nodes.
	glm::vec2 center = (min + max) * 0.5f;
	float length = std::max(max.x - min.x, max.y - min.y);

	// ---------
	// | 3 | 2 |
	// ---------
	// | 0 | 1 |
	// ---------
	glm::vec2 subMin[4];
	glm::vec2 subMax[4];

	//! get the four subnodes' region.
	subMin[0] = min;
	subMax[0] = center;
	subMin[1] = center - glm::vec2(0.0f, length / 2);
	subMax[1] = center + glm::vec2(length / 2, 0.0f);
	subMin[2] = center;
	subMax[2] = max;
	subMin[3] = min + glm::vec2(0.0f, length / 2);
	subMax[3] = center + glm::vec2(0.0f, length / 2);

	//! subdivide the objects into four classes according to their positions.
	std::vector<glm::vec2> classes[4];
	for (auto &point : objects)
	{s
		if (isContain(point, subMin[0], subMax[0]))
			classes[0].push_back(point);
		else if (isContain(point, subMin[1], subMax[1]))
			classes[1].push_back(point);
		else if (isContain(point, subMin[2], subMax[2]))
			classes[2].push_back(point);
		else if (isContain(point, subMin[3], subMax[3]))
			classes[3].push_back(point);
	}

	//! allocate memory for current node.
	TreeNode *cur = new TreeNode(min, max);
	cur->children[0] = recursiveBuild(depth + 1, subMin[0], subMax[0], classes[0]);
	cur->children[1] = recursiveBuild(depth + 1, subMin[1], subMax[1], classes[1]);
	cur->children[2] = recursiveBuild(depth + 1, subMin[2], subMax[2], classes[2]);
	cur->children[3] = recursiveBuild(depth + 1, subMin[3], subMax[3], classes[3]);

	return cur;
}
  这里需要提一点的是上面代码中的$isContain$函数，它输入一个点和一个矩形包围盒范围，判断这个点是否在这个矩形包围盒之内，这个判断很简单，不再赘述。我们来看一下这种自顶向下构建方法的时间复杂度。四叉树的**每一层**都会处理$O(N)$个数据对象，若数据分布的比较均匀，则树高为$O(log_4 N)$，因此构建的时间复杂度为$O(N log_4 N)$，乍一看比暴力方法的时间复杂度还要高，但我们对每一个数据对象没有做很复杂的数据运算（如射线求交、碰撞检测等等），而且构建只需一次，对于那些检索密集型的应用来说非常划算。 #### 3、四叉树分割算法——销毁   我们采用四叉链表结构作为四叉树的实现方式，涉及到大量的动态指针，在销毁时需要手动释放这些指针占用堆空间。销毁四叉树并不难，直接采用**后序遍历**方法即可，后序遍历即先销毁子节点，将子节点指针置空，然后再销毁父节点，如此递归下去。 
1
2
3
4
5
6
7
8
9
10
11
12
13
void QuadTree::recursiveDestory(TreeNode *node)
{
	if (node == nullptr)
		return;

	recursiveDestory(node->children[0]);
	recursiveDestory(node->children[1]);
	recursiveDestory(node->children[2]);
	recursiveDestory(node->children[3]);

	delete node;
	node = nullptr;
}
#### 4、四叉树分割算法——插入   除了前面提到的给定数据对象列表直接自顶向下构建整颗四叉树，还有一种构建方法是动态插入法。即给一个数据对象，将其插入到四叉树的一个叶子节点上，可以称为逐个插入法。逐个插入法就是要找到输入的数据对象所在的四叉树叶子节点，然后将该数据对象添加到该叶子节点的数据对象列表上。插入过程要靠的情况稍微多一点，如果插入的叶子节点的列表长度超过了一定数量，则应该将当前的叶子节点分裂，这时它不再是叶子节点而是内部节点了，递归插入到该内部节点的叶子节点中。总的来说，插入过程分成以下几种情况： - 若输入的数据对象不在当前节点的覆盖范围之内，直接停止插入过程； - 如果走到了四叉树的最大深度，则一定不会再分裂，直接将该数据对象添加到节点的数据对象列表当中； - 如果当前节点是叶子节点，将数据对象添加到该节点的数据对象列表中。然后若数据对象表列长度超过了一定数量，则需要将该叶子节点分裂，往下递归划分，将数据对象列表分散到子节点中，自己不再是叶子节点； - 如果当前节点是内部节点，则需要确定当前的数据对象落到它的四个子节点中的哪一个，然后递归调用插入过程。   插入过程也不是很复杂，理清了思路就好，实现的代码如下所示，代码的输入是深度、当前节点指针、当前节点覆盖的范围以及要插入的数据对象： 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
bool QuadTree::recursiveInsert(unsigned int depth, TreeNode * node,
                               glm::vec2 min, glm::vec2 max, glm::vec2 object)
{
	if (!isContain(object, min, max))
		return false;

	glm::vec2 center = (max + min) * 0.5f;
	float length = std::max(max.x - min.x, max.y - min.y);
	//! get the four sub-nodes' region.
	glm::vec2 subMin[4];
	glm::vec2 subMax[4];
	subMin[0] = min;
	subMax[0] = center;
	subMin[1] = center - glm::vec2(0.0f, length / 2);
	subMax[1] = center + glm::vec2(length / 2, 0.0f);
	subMin[2] = center;
	subMax[2] = max;
	subMin[3] = min + glm::vec2(0.0f, length / 2);
	subMax[3] = center + glm::vec2(0.0f, length / 2);

	//! reach the max depth.
	if (depth == mMaxDepth)
	{
		node->objects.push_back(object);
		return true;
	}

	if (node->isLeaf)
	{
		node->objects.push_back(object);
		if (node->objects.size() > 4)
		{
			//! 超过四个就分裂，自己不再是叶子节点
			node->children[0] = new TreeNode(subMin[0], subMax[0]);
			node->children[1] = new TreeNode(subMin[1], subMax[1]);
			node->children[2] = new TreeNode(subMin[2], subMax[2]);
			node->children[3] = new TreeNode(subMin[3], subMax[3]);
			node->isLeaf = false;
			node->children[0]->isLeaf = true;
			node->children[1]->isLeaf = true;
			node->children[2]->isLeaf = true;
			node->children[3]->isLeaf = true;

			for (auto &point : node->objects)
			{
				if (isContain(point, subMin[0], subMax[0]))
					node->children[0]->objects.push_back(point);
				else if (isContain(point, subMin[1], subMax[1]))
					node->children[1]->objects.push_back(point);
				else if (isContain(point, subMin[2], subMax[2]))
					node->children[2]->objects.push_back(point);
				else if (isContain(point, subMin[3], subMax[3]))
					node->children[3]->objects.push_back(point);
			}
			std::vector<glm::vec2>().swap(node->objects);
		}
		return true;
	}

	//! 若为内部节点，往下深入搜索
	if (isContain(object, subMin[0], subMax[0]))
		return recursiveInsert(depth + 1, node->children[0], subMin[0], subMax[0], object);
	if (isContain(object, subMin[1], subMax[1]))
		return recursiveInsert(depth + 1, node->children[1], subMin[1], subMax[1], object);
	if (isContain(object, subMin[2], subMax[2]))
		return recursiveInsert(depth + 1, node->children[2], subMin[2], subMax[2], object);
	if (isContain(object, subMin[3], subMax[3]))
		return recursiveInsert(depth + 1, node->children[3], subMin[3], subMax[3], object);

	return false;
}
  插入一个数据对象的时间复杂就是四叉树的高度，理想情况下为$O(log_4 N)$。 #### 5、四叉树分割算法——删除   有插入就有删除，给定一个数据对象，我们要找到包含这个数据对象的节点，并将其从数据对象列表中删除。与插入的情况相反，删除对象的操作可能导致当前节点的数据对象列表为空，此时需要将当前节点从整颗四叉树中删除，因为它不再包含任何有效的数据。叶子节点的删除可能导致父节点的删除（全部子节点变为空指针），如此递归下去，因此需要仔细斟酌整个删除的过程。总的来说，删除过程需要考虑的情况有如下几种： - 当前节点为空或者当前节点覆盖的范围不包含输入的数据对象，直接终止删除过程； - 当前节点为叶子节点或者走到最大四叉树深度了，遍历当前节点的数据对象列表，找到需要删除对象的数组下标，如果没有则不执行删除操作，否则将其从数组中移除。如果移除一个对象之后的数据列表为空，则需要将其删除，返回一个标志告诉父节点将其删除； - 当前节点为内部节点，找到包含该数据对象的子节点，递归调用删除程序，如果程序返回一个删除标志，则释放该子节点的内存。如果内部节点的四个子节点指针均为空指针，则告诉该内部节点的父节点将其删除，如果递归下去。   将数据对象从数组中删除一个小小的技巧就是将数组尾部的内容覆盖到该数据对象所在的位置，然后将尾部的数据删除即可，避免数据的大规模移动。 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
bool QuadTree::recursiveRemove(unsigned int depth, TreeNode * node, glm::vec2 min, glm::vec2 max, glm::vec2 object)
{
	if (node == nullptr)
		return false;

	//! 确保在当前节点的范围之内
	if (!isContain(object, min, max))
		return false;

	//! 达到最大深度或者走到叶子节点了.
	if (depth == mMaxDepth || node->isLeaf)
	{
		//! 找到要删除元素的位置
		int index = -1;
		for (size_t i = 0; i < node->objects.size(); ++i)
		{
			if ((std::pow(node->objects[i].x - object.x, 2) + std::pow(node->objects[i].y - object.y, 2)) < 0.001f)
			{
				index = i;
				break;
			}
		}
		if (index != -1)
		{
			//! 与最后一个元素交换，方便删除
			node->objects[index] = node->objects.back();
			node->objects.pop_back();
		}
		//! 如果叶子节点空了，则告诉父节点需要将其删除掉.
		if (node->objects.empty())
			return true;
		return false;
	}

	//! 非叶子节点
	if (recursiveRemove(depth + 1, node->children[0], node->children[0]->bMin, node->children[0]->bMax, object))
	{
		delete node->children[0];
		node->children[0] = nullptr;
	}
	else if (recursiveRemove(depth + 1, node->children[1], node->children[1]->bMin, node->children[1]->bMax, object))
	{
		delete node->children[1];
		node->children[1] = nullptr;
		return false;
	}
	else if (recursiveRemove(depth + 1, node->children[2], node->children[2]->bMin, node->children[2]->bMax, object))
	{
		delete node->children[2];
		node->children[2] = nullptr;
	}
	else if (recursiveRemove(depth + 1, node->children[3], node->children[3]->bMin, node->children[3]->bMax, object))
	{
		delete node->children[3];
		node->children[3] = nullptr;
	}

	//! 父节点的全部子节点为空，自己也将被删除了
	if (node->children[0] == nullptr && node->children[1] == nullptr 
		&& node->children[2] == nullptr && node->children[3] == nullptr)
	{
		return true;
	}

	return false;
}
#### 6、四叉树分割算法——访问   最后，为了将四叉树的空间划分结果可视化出来，需要对整颗四叉树做一个遍历。树形数据结构遍历都是一些很基础的基本功，这里不再赘述，每遍历到一个节点，将该节点的矩形区域范围绘制出来。 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
void QuadTree::recursiveTraverse(TreeNode *node, glm::vec2 min, glm::vec2 max, std::vector<glm::vec2>& lines)
{
	glm::vec2 center = (max + min) * 0.5f;
	float length = std::max(max.x - min.x, max.y - min.y);
	glm::vec2 corners[4];
	corners[0] = min;
	corners[1] = min + glm::vec2(length, 0.0f);
	corners[2] = max;
	corners[3] = min + glm::vec2(0.0f, length);

	//! get the bounding box to draw.
	lines.push_back(corners[0]);
	lines.push_back(corners[1]);

	lines.push_back(corners[1]);
	lines.push_back(corners[2]);

	lines.push_back(corners[2]);
	lines.push_back(corners[3]);

	lines.push_back(corners[3]);
	lines.push_back(corners[0]);

	if (node == nullptr || node->isLeaf)
		return;

	//! get the four sub-nodes' region.
	glm::vec2 subMin[4];
	glm::vec2 subMax[4];
	subMin[0] = min;
	subMax[0] = center;
	subMin[1] = center - glm::vec2(0.0f, length / 2);
	subMax[1] = center + glm::vec2(length / 2, 0.0f);
	subMin[2] = center;
	subMax[2] = max;
	subMin[3] = min + glm::vec2(0.0f, length / 2);
	subMax[3] = center + glm::vec2(0.0f, length / 2);

	recursiveTraverse(node->children[0], subMin[0], subMax[0], lines);
	recursiveTraverse(node->children[1], subMin[1], subMax[1], lines);
	recursiveTraverse(node->children[2], subMin[2], subMax[2], lines);
	recursiveTraverse(node->children[3], subMin[3], subMax[3], lines);
}
#### 7、四叉树可视化实验结果   采用OpenGL图形API将烟花粒子系统和它的四叉树结构可视化出来如下图所示，烟花粒子系统的实现和OpenGL渲染API的使用超出了本课程的范围，因此不再赘述。这里构建的四叉树数据对象是粒子点云，下面一系列图片中左边的图展示是烟花粒子效果，右边的图则展示了构建出来的粒子点云四叉树的可视化结果（绿色部分）。可以看到空间分割的情况符合我们的预期，在粒子点密集的区域四叉树往下分割了不少层，而点云稀疏甚至没有的部分不会往下继续分割，这个就是二维分割的四叉树结构，当我们要搜索某一个区域周围的粒子点时就可以通过四叉树迅速地进行邻域查询。
  关于四叉树的可视化部分就到这里，下面来看看三维八叉树在计算机图形学光线追踪渲染技术中的实际应用和性能提升。 ## 三、利用八叉树算法加速射线求交   四叉树和八叉树从本质上来讲没有区别，算法的思想都是一样的，只不过四叉树是针对二维的情况，而八叉树针对的是三维的情况。在三维空间均匀划分的是一个长方体，每一个长方体分割成八个相同的子长方体，因而三维情况的空间分割树有八个子节点，被称为八叉树。Wiki的八叉树定义为： 八叉树（Octree）是一种用于描述三维空间的树状数据结构。八叉树的每个节点表示一个正方体的体积元素，每个节点有八个子节点，这八个子节点所表示的体积元素加在一起就等于父节点的体积。一般中心点作为节点的分叉中心。
图3 三维空间的八叉树

  想象一个立方体，我们最少可以切成多少个相同等分的小立方体？答案就是8个。再想象 我们有一个房间，房间里某个角落藏着一枚金币，我们想很快的把金币找出来。 我们可以把房间当成一个立方体，先切成八个小立方体，然后排除掉没有放任何东西的小立方体，再把有可能藏金币的小立方体继续切八等份…. 如此下去，平均在$Log_8$(房间内的所有物品数)的时间内就可找到金币。八叉树就是用在3D空间中的场景管理，可以加速我们的物体搜索、碰撞检测、射线求交、邻域搜索等等空间查找操作。

  一般情况下，八叉树的构建过程如下：

  • (1).设定最大递归深度;
  • (2).找出场景的最大尺寸，并以此尺寸建立第一个立方体;
  • (3).依序将单位元元素丢入能被包含且没有子节点的立方体;
  • (4).若没有达到最大递归深度，就进行细分八等份，再将该立方体所装的单位元元素全部分担给八个子立方体;
  • (5).若发现子立方体所分配到的单位元元素数量不为零且跟父立方体是一样的，则该子立方体停止细分，因为跟据空间分割理论，细分的空间所得到的分配必定较少，若是一样数目，则再怎么切数目还是一样，会造成无穷切割的情形；
  • (6).重复步骤(3)，直到达到最大递归深度。

  光线追踪是一个递归的过程。发射一束光线到场景，求出光线和几何图形间最近的交点，如果该交点的材质是反射性或折射性的，可以在该交点向反射方向或折射方向继续追踪，如此递归下去，直到设定的最大递归深度或者射线追踪到光源处（或者背景色），如此便计算处一个像素的着色值。光线追踪算法涉及到大量的光线与几何体的求交点运算，因而与渲染效率息息相关。复杂物体大多采用三角形网格表示，因而进一步来说是射线与三角形的求交运算。下面我将为每个三角网格模型构建一颗八叉树，加速射线与三角网格模型的求交。

1、八叉树分割算法——构建与销毁

  在开始构建之前，我们首先要弄清楚一点的就是在这里数据对象是一个个三角形而非一个个点，三角形是由一定的大小的，它有可能刚好处在分割平面上，导致三角形跨越了多个子节点。这时为了防止遗漏，应该把三角形都存储到由交集的节点中。同时，每个叶子节点中并不会直接存储三角形图元数据，而是存储三角形面片的索引。定义八叉树节点如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
struct Node
{
	//! it's leaf or not.
	bool m_isLeaf;
	//! the eight children.
	Node *m_children[8];
	//! extent.
	AABB m_boundingBox;
	//! objects' id list.
	std::vector<std::tuple<unsigned int, unsigned int, unsigned int>> m_objectIds;

	Node() : m_isLeaf(true), m_boundingBox(AABB(Vector3D(0.0f, 0.0f, 0.0f), Vector3D(0.0f, 0.0f, 0.0f)))
	{
		m_children[0] = m_children[1] = m_children[2] = m_children[3] = nullptr;
		m_children[4] = m_children[5] = m_children[6] = m_children[7] = nullptr;
	}

	Node(Vector3D min, Vector3D max, bool isLeaf) : m_isLeaf(isLeaf),
		m_boundingBox(AABB(min, max))
	{
		m_children[0] = m_children[1] = m_children[2] = m_children[3] = nullptr;
		m_children[4] = m_children[5] = m_children[6] = m_children[7] = nullptr;
	}

};

  上面的$m_objectIds$s是一个三元组的vector，三元组的每个元素存储一个三角形面片的三个顶点索引。然后输入模型的三角形面片列表，我们自顶向下构建八叉树。构建的过程其实跟四叉树的构建过程差别不大，区别在于子区域的划分：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
Node *Octree::recursiveBuild(unsigned int depth, Vector3D min, Vector3D max,
	const MeshHitable * target, std::vector<unsigned int> ids)
{
	//! reach the max depth or there are just a few objects, we dont' need to go further.
	if (depth == m_maxDepth || ids.size() <= 10 * 3)
	{
		Node *cur = new Node(min, max, true);
		const auto &vertices = target->m_vertices;
		for (size_t i = 0; i < ids.size(); i += 3)
		{
			if (isContain(cur->m_boundingBox, vertices[ids[i + 0]].m_position,
				vertices[ids[i + 1]].m_position, vertices[ids[i + 2]].m_position))
				cur->m_objectIds.push_back(std::make_tuple(ids[i + 0], ids[i + 1], ids[i + 2]));
		}
		return cur;
	}

	//! otherwise, just divide into 8 sub-nodes.
	Node *cur = new Node(min, max, false);
	std::vector<unsigned int> subIds[8];
	const auto &vertices = target->m_vertices;
	std::vector<AABB> subRegions = cur->m_boundingBox.getEightSubAABB();
	for (size_t i = 0; i < ids.size(); i += 3)
	{
		const auto &p1 = vertices[ids[i + 0]].m_position;
		const auto &p2 = vertices[ids[i + 1]].m_position;
		const auto &p3 = vertices[ids[i + 2]].m_position;
		for (size_t j = 0; j < 8; ++j)
		{
			if (isContain(subRegions[j], p1, p2, p3))
			{
				subIds[j].push_back(ids[i + 0]);
				subIds[j].push_back(ids[i + 1]);
				subIds[j].push_back(ids[i + 2]);
			}
		}
	}

	cur->m_children[0] = recursiveBuild(depth + 1, subRegions[0].getMin(), subRegions[0].getMax(), target, subIds[0]);
	cur->m_children[1] = recursiveBuild(depth + 1, subRegions[1].getMin(), subRegions[1].getMax(), target, subIds[1]);
	cur->m_children[2] = recursiveBuild(depth + 1, subRegions[2].getMin(), subRegions[2].getMax(), target, subIds[2]);
	cur->m_children[3] = recursiveBuild(depth + 1, subRegions[3].getMin(), subRegions[3].getMax(), target, subIds[3]);
	cur->m_children[4] = recursiveBuild(depth + 1, subRegions[4].getMin(), subRegions[4].getMax(), target, subIds[4]);
	cur->m_children[5] = recursiveBuild(depth + 1, subRegions[5].getMin(), subRegions[5].getMax(), target, subIds[5]);
	cur->m_children[6] = recursiveBuild(depth + 1, subRegions[6].getMin(), subRegions[6].getMax(), target, subIds[6]);
	cur->m_children[7] = recursiveBuild(depth + 1, subRegions[7].getMin(), subRegions[7].getMax(), target, subIds[7]);

	return cur;
}

  构建过程的时间复杂度为$Nlog_8\ N$。同样的，八叉树销毁过程采用后序遍历的方式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
void Octree::recursiveDestory(Node * node)
{
	//! backward traverse for deletion.
	if (node == nullptr)
		return;

	recursiveDestory(node->m_children[0]);
	recursiveDestory(node->m_children[1]);
	recursiveDestory(node->m_children[2]);
	recursiveDestory(node->m_children[3]);
	recursiveDestory(node->m_children[4]);
	recursiveDestory(node->m_children[5]);
	recursiveDestory(node->m_children[6]);
	recursiveDestory(node->m_children[7]);

	delete node;
	node = nullptr;
}

  这里有一点需要特别注意的就是需要判断一个立方体区域是否与三角形相交或包含，我采用的方法是首先判断三角形的三个点中是否至少有一个在立方体内部，如果是则返回真。而如果三角形的三个顶点均不在立方体区域内部，但仍有可能三角形与立方体区域存在交集，此时可以通过依次判断三角形的三条边是否与立方体相交。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
bool Octree::isContain(const AABB & box, const Vector3D & p1, const Vector3D & p2, const Vector3D & p3)
{
	//! at least one point of the triangle is inside the box, return true.
	if (box.isInside(p1) || box.isInside(p2) || box.isInside(p3))
		return true;

	//! otherwise testing for each edge.
	Ray edge1(p1, p2 - p1);
	Ray edge2(p2, p3 - p2);
	Ray edge3(p3, p1 - p3);
	float length1 = (p2 - p1).getLength();
	float length2 = (p3 - p2).getLength();
	float length3 = (p1 - p3).getLength();
	if (box.hit(edge1, 0.0f, length1) || box.hit(edge2, 0.0f, length2) || box.hit(edge3, 0.0f, length3))
		return true;

	//! if there is no intersection at all, just return false.
	return false;
}

2、八叉树分割算法——访问

  构建了八叉树之后，八叉树的优势就体现出来了，在光线追踪过程每发射一条射线我们就要检测该射线与场景中的物体是否发生了碰撞，如果发生了碰撞则需要将碰撞的点计算出来。有了八叉树结构，我们首先检测射线是否与当前节点的区域发生了碰撞，如果不发生碰撞则直接退出；而如果发生了碰撞且当前节点是叶子节点的话，则遍历叶子节点存储的三角形，计算射线与三角形的交点，取距离射线源最近的结果；而当前节点是内部节点的话，则递归调用遍历程序，判断射线是否与该区域的八个子区域碰撞，有碰撞才会继续往下走，无碰撞则不会往下深入。

  这个算法的实现并不复杂，实现代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
bool Octree::recursiveTraveling(Node * node, const Ray & ray, float & t_min, float & t_max,
	std::function<bool(unsigned int i1, unsigned int i2, unsigned int i3, float &t)> func)
{
	//! just return false if node is nullptr.
	if (node == nullptr)
		return false;

	//! if the ray doesn't intersect with the aabb box, just return false.
	if (!node->m_boundingBox.hit(ray, t_min, t_max))
		return false;

	//! if it's a leaf, travels every objects of this node.
	if (node->m_isLeaf)
	{
		bool ret = false;
		for (size_t i = 0; i < node->m_objectIds.size(); ++i)
		{
			const auto &face = node->m_objectIds[i];
			ret |= func(std::get<0>(face), std::get<1>(face), std::get<2>(face), t_max);
		}
		return ret;
	}

	//! otherwise, divide into 8 sub-nodes.
	bool ret = false;
	std::vector<AABB> subRegions = node->m_boundingBox.getEightSubAABB();
	if (subRegions[0].hit(ray, t_min, t_max))
		ret |= recursiveTraveling(node->m_children[0], ray, t_min, t_max, func);
	if (subRegions[1].hit(ray, t_min, t_max))
		ret |= recursiveTraveling(node->m_children[1], ray, t_min, t_max, func);
	if (subRegions[2].hit(ray, t_min, t_max))
		ret |= recursiveTraveling(node->m_children[2], ray, t_min, t_max, func);
	if (subRegions[3].hit(ray, t_min, t_max))
		ret |= recursiveTraveling(node->m_children[3], ray, t_min, t_max, func);
	if (subRegions[4].hit(ray, t_min, t_max))
		ret |= recursiveTraveling(node->m_children[4], ray, t_min, t_max, func);
	if (subRegions[5].hit(ray, t_min, t_max))
		ret |= recursiveTraveling(node->m_children[5], ray, t_min, t_max, func);
	if (subRegions[6].hit(ray, t_min, t_max))
		ret |= recursiveTraveling(node->m_children[6], ray, t_min, t_max, func);
	if (subRegions[7].hit(ray, t_min, t_max))
		ret |= recursiveTraveling(node->m_children[7], ray, t_min, t_max, func);

	return ret;
}

  可以看到，通过八叉树的结构我们可以避免很多无用功，大大加速整个求交的过程，有了八叉树后射线与三角形求交的时间复杂度就变成了$O(log_8 N)$，其中$N$是输入规模，相比于原来的$O(N)$复杂度提升非常明显。下面就通过实验比较有无八叉树的算法效率。至于光线追踪方面的算法，与本门课的主题关系比较小，而且涉及的内容实在太多，因此不再赘述。

3、八叉树加速实验结果

  我们将对比两个场景渲染的实验，两个场中均涉及到复杂的三角网格模型。控制的变量有每个像素发射的采样光线个数、有无八叉树加速。一般情况下，采样光线数量越多越好，数量过少会导致渲染出来的图片噪声过多，质量太低。场景中三角网格比较多的模型如下图5所示：

图5 从左到右三角形面片的数量依次为100000个、69666个和5022个
  下面给出了渲染绘制出来的两个场景，从左到右、从上到下的光线采样数量依次为16、64、128和512，可以看到低采样数时渲染出来的图片噪声明显，质量很低。因此为了渲染出高质量的图片，采样数不能太低，这意味着求交数量也将大大增加。我们记下面的两个场景分别为dragonSquare和dragonBox。下面的图片可以看到利用八叉树并没有损失光线追踪的渲染质量。
图6 场景一：dragonSquare

图7 场景二：dragonBox

  下面的表1给出了两个场景在有无八叉树时的渲染时间对比，表格中的s是second的缩写，即单位为秒。表格的最后一列有两个“>3小时“是因为等不下去了，就不再等了，花费的时间实在太长了。通过表格可以看到，有了八叉树的加速，时间效率提升极其明显，特别是采样数量越多的时候，提升倍数越多，完全是几个数量级的提升。这是因为采样数量呈指数增值时，射线与三角形的求交数量也呈指数增长。

表1 渲染时间比较

scene 16次采样 64次采样 128次采样 512次采样
dragonSquare（无八叉树） 427.776s 2195.918s 8715.48s >3小时
dragonSquare（有八叉树） 17.824s 70.382s 147.72s 614.562s
dragonBox（无八叉树） 380.944s 1818.491s 6548.502 >3小时
dragonBox（有八叉树） 13.136s 58.661s 105.642s 422.652s

四、总结

  分治算法思想是一个非常有用的算法思想，我们围绕分治算法这一主题展开其在计算机图形学中应用——四叉树和八叉树空间分割算法实践。抛开诸多细节不说，四叉树和八叉树在本质上与二分查找类似。与二分查找的排序预处理一样，四叉树和八叉树都需要在最开始构建，构建的复杂度与排序的时间复杂度一致，对于那些查找密集型的应用场景来说非常划算（提升的效率极其恐怖）。当然，四叉树和八叉树并不完美的，它们也有缺点，一个比较大的缺点就是它要求空间物体对象分布比较均匀，如果分布不均匀，那么构建出来的树将不会很平衡，这将影响后续的访问效率。针对这一问题，后续的学者已经有比较多的算法变种，限于篇幅这里不再详细展开。

 上一页

Computer Graphics Ray Tracer

 评论


博客内容遵循 署名-非商业性使用-相同方式共享 4.0 国际 (CC BY-NC-SA 4.0) 协议

本站使用 Material X 作为主题 ， 总访问量为 次 。
载入天数...载入时分秒...