数据结构(1):简单八叉树学习--SimpleOctree

该代码定义了八叉树的数据结构，并在main.cpp中演示了“查询某一空间范围框内有哪些八叉树节点”这个demo

1 文件目录：

2 Vec3.h

3 OctreePoint.h

4 Stopwatch.h

5 Octree.h

6 main.cpp

7 比较八叉树查询和蛮力查询

1 文件目录：

其中README.md内容如下：

SimpleOctree
============A simple octree with good commenting for learning how octrees work. Blog post incoming with description, or read comments in Octree.hUsage
============
make && ./octree

Makefile文件内容：

default:g++ main.cpp -O3 -o octree
clean:@rm -f octree

LICENSE文件应该不需要，这里不展示了

2 Vec3.h

创建一个结构体容器Vec3，容器内包含三个float元素(三个元素组成的结构体/数组)，并定义该容器的各种运算符--operator

#ifndef Vec3_h_
#define Vec3_h_#include <cmath>struct Vec3;//声明要创建的数据结构
Vec3 operator*(float r, const Vec3& v);//声明要创建的函数，该函数在文件末尾，不在结构体内struct Vec3 {/*union共用体，也叫联合体，在一个“联合”内可以定义多种不同的数据类型， 一个被说明为该“联合”类型的变量中，允许装入该“联合”所定义的任何一种数据，这些数据共享同一段内存，以达到节省空间的目的。union变量所占用的内存长度等于最长的成员的内存长度。*/union {struct {float x,y,z;};float D[3];};Vec3() { } //初始化不赋值Vec3(float _x, float _y, float _z)//初始化并赋值:x(_x), y(_y), z(_z){ }float& operator[](unsigned int i) {//索引操作，会改变被索引对象return D[i];}const float& operator[](unsigned int i) const {//索引操作，不会改变被索引对象return D[i];}float maxComponent() const { //取最大值操作float r = x;if(y>r) r = y;if(z>r) r = z;return r;}float minComponent() const { //取最小值操作float r = x;if(y<r) r = y;if(z<r) r = z;return r;}Vec3 operator+(const Vec3& r) const { //元素相加操作return Vec3(x+r.x, y+r.y, z+r.z); }Vec3 operator-(const Vec3& r) const { //元素相减操作return Vec3(x-r.x, y-r.y, z-r.z); }Vec3 cmul(const Vec3& r) const { //元素相乘操作return Vec3(x*r.x, y*r.y, z*r.z);}Vec3 cdiv(const Vec3& r) const { //元素相除操作return Vec3(x/r.x, y/r.y, z/r.z);}Vec3 operator*(float r) const { //乘法操作return Vec3(x*r,y*r,z*r);}Vec3 operator/(float r) const { //除法操作return Vec3(x/r, y/r, z/r);}Vec3& operator+=(const Vec3& r) { //指针-元素相加操作x+=r.x;y+=r.y;z+=r.z;return *this;}Vec3& operator-=(const Vec3& r) { //指针-元素相减操作x-=r.x;y-=r.y;z-=r.z;return *this;}Vec3& operator*=(float r) { //指针-乘法操作x*=r; y*=r; z*=r;return *this;}// Inner/dot productfloat operator*(const Vec3& r) const { //对应元素平方和操作return x*r.x + y*r.y + z*r.z;}float norm() const {//对应元素平方和后开方操作return sqrtf(x*x+y*y+z*z);}float normSquared() const { //元素平方和操作return x*x + y*y + z*z;}// Cross productVec3 operator^(const Vec3& r) const { //元素交叉操作return Vec3(y * r.z - z * r.y, z * r.x - x * r.z, x * r.y - y * r.x);}Vec3 normalized() const { //归一化操作return *this / norm();}
};inline Vec3 operator*(float r, const Vec3& v) { //结构体数据乘以一个数值return Vec3(v.x*r, v.y*r, v.z*r);
}#endif

3 OctreePoint.h

声明八叉树中简单点OctreePoint的数据类型、初始化方式、以及OctreePoint位置position的两个操作函数。

#ifndef OctreePoint_H
#define OctreePoint_H#include "Vec3.h"// Simple point data type to insert into the tree.
// Have something with more interesting behavior inherit  from this in order to store other attributes in the tree.
// 要插入八叉树中简单点的数据类型。
// 继承一些具有更有趣行为的东西，以便在树中存储其他属性。
class OctreePoint {Vec3 position;
public:OctreePoint() { }//初始化OctreePointOctreePoint(const Vec3& position) : position(position) { } //初始化OctreePoint，并给OctreePoint的position赋值inline const Vec3& getPosition() const { return position; }  //返回/获取位置数据inline void setPosition(const Vec3& p) { position = p; }        //设置位置数据
};#endif

4 Stopwatch.h

定义一个用于计时的函数stopwatch()，时间单位为秒

#ifndef _WIN32
#include <sys/time.h>double stopwatch()
{struct timeval time;gettimeofday(&time, 0 );return 1.0 * time.tv_sec + time.tv_usec / (double)1e6; //单位为秒
}#else#include <windows.h>
double stopwatch()
{unsigned long long ticks;unsigned long long ticks_per_sec;QueryPerformanceFrequency( (LARGE_INTEGER *)&ticks_per_sec);QueryPerformanceCounter((LARGE_INTEGER *)&ticks);return ((float)ticks) / (float)ticks_per_sec;
}#endif

5 Octree.h

声明一个八叉树的类，基本组成为:当前节点的位置origin、当前节点对应空间范围halfDimension、当前节点的子节点children、当前节点的数据data。

然后递归定义了八叉树当前节点的八个子节点的八叉树；以及子节点序号查询函数和某一空间范围框内八叉树节点(当前节点、子节点、子子节点。。。。)的查询方法

#ifndef Octree_H
#define Octree_H#include <cstddef>
#include <vector>
#include "OctreePoint.h"namespace brandonpelfrey {/**!*定义一个八叉树的类，其组成包括:当前节点中心位置origin、当前节点对应空间尺寸halfDimension、当前节点的子节点children、当前节点的数据data*/class Octree {// Physical position/size. This implicitly defines the bounding box of this node   // 物体位置和大小，间接定义了当前节点框边界Vec3 origin;         //! The physical center of this node  该节点的物理中心Vec3 halfDimension;  //! Half the width/height/depth of this node 该节点的半宽/高/深//这棵树最多有八个子节点，分别可以额外存储一个点，不过在许多应用程序中，叶子将存储数据。// The tree has up to eight children and can additionally store // a point, though in many applications only, the leaves will store data.Octree *children[8]; //! Pointers to child octants   指向八叉树子节点的指针OctreePoint *data;   //! Data point to be stored at a node  一个指针--指向存储在一个八叉树节点中的数据/*Children follow a predictable pattern to make accesses simple.八叉树子节点遵循可预测的模式以使访问变得简单。Here, - means less than 'origin' in that dimension, + means greater than.这里，- 表示小于该维度(坐标)的“原点”，+ 表示大于。child:   0 1 2 3 4 5 6 7  //从下面xyz的+-表示中可以看出这个八个子节点的空间分布x:      - - - - + + + +y:      - - + + - - + +z:      - + - + - + - +*/public:Octree(const Vec3& origin, const Vec3& halfDimension) //初始化一个八叉树: origin(origin), halfDimension(halfDimension), data(NULL) {// Initially, there are no childrenfor(int i=0; i<8; ++i) children[i] = NULL;}Octree(const Octree& copy) //从另一个八叉数复制得到一个八叉树: origin(copy.origin), halfDimension(copy.halfDimension), data(copy.data) {}~Octree() {// Recursively destroy octants //递归删除子节点for(int i=0; i<8; ++i) delete children[i];}// 确定point在该八叉树的哪一个子节点空间内// Determine which octant of the tree would contain 'point'int getOctantContainingPoint(const Vec3& point) const {int oct = 0;if(point.x >= origin.x) oct |= 4; //"|="位或操作if(point.y >= origin.y) oct |= 2;if(point.z >= origin.z) oct |= 1;return oct;//这个值应该是0、1、2、3、4、5、6、7之中的一个}//检查是否是叶子节点(叶子节点--没有子节点的八叉树节点)bool isLeafNode() const {// This is correct, but overkill--这是正确的，但矫枉过正. See below./*for(int i=0; i<8; ++i)if(children[i] != NULL) return false;return true;*/// We are a leaf if we have no children. Since we either have none, or  all eight, it is sufficient to just check the first.// 如果没有子节点，就是一片叶子节点。 由于要么没有，要么全部都没有，所以只检查第一个就足够了。return children[0] == NULL;}//插入一个简单数据点，即将简单数据点插入到八叉树中void insert(OctreePoint* point) {// If this node doesn't have a data point yet assigned  and it is a leaf, then we're done!// 如果八叉树当前节点还没有分配一个数据点并且它是一个叶子节点，那么我们就完成了----把该简单数据点赋值到当前节点即可if(isLeafNode()) { //是否是叶子节点----没有子节点if(data==NULL) { //当是叶子节点并且没有被赋值时，将简单数据点赋值给它data = point;return;} else { // We're at a leaf, but there's already something here// We will split this node so that it has 8 child octants and then insert the old data that was here, along with this new data point // Save this data point that was here for a later re-insert// 是叶子节点并且已经被赋值时// 将分割这个叶子节点，使其有 8 个子节点，然后插入这里的旧数据以及这个新数据点// 保存此数据点，以便稍后重新插入OctreePoint *oldPoint = data; //保存叶子节点的数据data = NULL;// Split the current node and create new empty trees for each  child octant.// 拆分当前叶子节点，并为每个子节点创建新的空八叉树。for(int i=0; i<8; ++i) {// Compute new bounding box for this child 为子节点创建新的边界框Vec3 newOrigin = origin;//&运算符、三目运算符newOrigin.x += halfDimension.x * (i&4 ? .5f : -.5f);//i为4时newOrigin.y += halfDimension.y * (i&2 ? .5f : -.5f);//i为2、6时newOrigin.z += halfDimension.z * (i&1 ? .5f : -.5f);//i为1、3、5、7时children[i] = new Octree(newOrigin, halfDimension*.5f);}// Re-insert the old point, and insert this new point 重新插入旧点，并插入这个新点// (We wouldn't need to insert from the root, because we already know it's guaranteed to be in this section of the tree)//（我们不需要从根插入，因为我们已经知道它保证在八叉树的某一个子节点对应的空间内）// 这里对子节点进行插入“简单数据点”操作，观察children[]里面得到的子节点的序号children[getOctantContainingPoint(oldPoint->getPosition())]->insert(oldPoint);children[getOctantContainingPoint(point->getPosition())]->insert(point);}} else {// We are at an interior node. Insert recursively into the  appropriate child octant// 当八叉树有子节点时，根据需要插入的简单数据点得到子节点的序号，并将该简单数据点插入其中int octant = getOctantContainingPoint(point->getPosition());children[octant]->insert(point);}}// This is a really simple routine for querying the tree for points within a bounding box defined by min/max points (bmin, bmax)// 这是一个非常简单的例程，用于查询八叉树中（bmin，bmax）定义的边界框内的八叉树节点// All results are pushed into 'results'void getPointsInsideBox(const Vec3& bmin, const Vec3& bmax, std::vector<OctreePoint*>& results) {// If we're at a leaf node, just see if the current data point is inside the query bounding box// 如果我们在叶子节点，只需查看当前叶子节点是否在查询边界框内if(isLeafNode()) {if(data!=NULL) {const Vec3& p = data->getPosition();//当前叶子节点的数据if(p.x>bmax.x || p.y>bmax.y || p.z>bmax.z) return;if(p.x<bmin.x || p.y<bmin.y || p.z<bmin.z) return;results.push_back(data);}} else {// We're at an interior node of the tree. We will check to see if  the query bounding box lies outside the octants of this node.// 对于八叉树的内部节点----非叶子节点即有子节点的节点。 我们将检查这八个子节点是否位于查询边界框内for(int i=0; i<8; ++i) {// Compute the min/max corners of this child octant 计算这个某个子节点的最小/最大范围Vec3 cmax = children[i]->origin + children[i]->halfDimension;Vec3 cmin = children[i]->origin - children[i]->halfDimension;// If the query rectangle is outside the child's bounding box,  then continue// 如果查询矩形在子节点的边界框之外，则继续（其他子节点或向上一级节点搜索）if(cmax.x<bmin.x || cmax.y<bmin.y || cmax.z<bmin.z) continue;if(cmin.x>bmax.x || cmin.y>bmax.y || cmin.z>bmax.z) continue;// At this point, we've determined that this child is intersecting  the query bounding box// 至此，我们已经确定这个子节点在查询边界框内，进一步进行搜索(这个子节点是否有子子节点。。。。)children[i]->getPointsInsideBox(bmin,bmax,results);} }}};}
#endif

6 main.cpp

创建了一组随机点，比较使用八叉树查询和使用强制/蛮力查询所需要的时间；并指出，当在查询相对接近整个八叉树空间大小的情况下，八叉树实际上会比蛮力/强制查询个点慢一点！

#include <cstdlib>
#include <cstdio>
#include <vector>
#include <algorithm>#include "Octree.h"
#include "Stopwatch.h"using namespace brandonpelfrey;// Used for testing
std::vector<Vec3> points;
Octree *octree;
OctreePoint *octreePoints;
Vec3 qmin, qmax;float rand11() // Random number between [-1,1] //返回[-1,1]之间的随机数字
{ return -1.f + (2.f*rand()) * (1.f / RAND_MAX); }//创建随机向量，其中组成元素在[-1,1]范围内
Vec3 randVec3() // Random vector with components in the range [-1,1]
{ return Vec3(rand11(), rand11(), rand11()); }// Determine if 'point' is within the bounding box [bmin, bmax] //确定点是否在框里面
bool naivePointInBox(const Vec3& point, const Vec3& bmin, const Vec3& bmax) {return point.x >= bmin.x &&point.y >= bmin.y &&point.z >= bmin.z &&point.x <= bmax.x &&point.y <= bmax.y &&point.z <= bmax.z;
}void init() {// Create a new Octree centered at the origin with physical dimension 2x2x2// 创建一个中心在原点的八叉树，物理尺寸为2*2*2octree = new Octree(Vec3(0,0,0), Vec3(1,1,1));//初始化在原点附近，为得是接下来创建随机点直观些// Create a bunch of random points 创建一组随机点const int nPoints = 1 * 1000 * 1000;for(int i=0; i<nPoints; ++i) {points.push_back(randVec3());//随机点的组成元素值在[-1,1]之间}/*在使用多个输出函数连续进行多次输出时，有可能发现输出错误。因为下一个数据再上一个数据还没输出完毕，还在输出缓冲区中时，下一个printf就把另一个数据加入输出缓冲区，结果冲掉了原来的数据，出现输出错误。 在 prinf（）；后加上fflush(stdout); 强制马上输出，避免错误。*/printf("Created %ld points\n", points.size()); fflush(stdout);// Insert the points into the octree  在八叉树中插入点octreePoints = new OctreePoint[nPoints];for(int i=0; i<nPoints; ++i) {octreePoints[i].setPosition(points[i]);//设置八叉树简单数据点的位置octree->insert(octreePoints + i);//插入简单数据点}printf("Inserted points to octree\n"); fflush(stdout);// Create a very small query box. The smaller this box is the less work the octree will need to do. // 创建一个非常小的查询框。 这个框越小，八叉树需要做的工作就越少。// This may seem like it is exagerating the benefits, but often, we only need to know very nearby objects.// 这似乎夸大了好处，但通常，我们只需要知道非常近的物体。qmin = Vec3(-.05,-.05,-.05);qmax = Vec3(.05,.05,.05);// Remember: In the case where the query is relatively close to the size of the whole octree space, the octree will// actually be a good bit slower than brute forcing every point!// 请记住：在查询相对接近整个八叉树空间大小的情况下，八叉树实际上会比蛮力/强制查询个点慢一点！
}// Query using brute-force 使用强制查询框内的点
void testNaive() {double start = stopwatch();//用于计时std::vector<int> results;for(int i=0; i<points.size(); ++i) {if(naivePointInBox(points[i], qmin, qmax)) {results.push_back(i);}}double T = stopwatch() - start;printf("testNaive found %ld points in %.5f sec.\n", results.size(), T);
}// Query using Octree 使用八叉树查询
void testOctree() {double start = stopwatch();std::vector<OctreePoint*> results;octree->getPointsInsideBox(qmin, qmax, results);double T = stopwatch() - start;printf("testOctree found %ld points in %.5f sec.\n", results.size(), T);
}int main(int argc, char **argv) {init();testNaive();testOctree();return 0;
}

7 比较八叉树查询和蛮力查询

根据下面对比，发现查询框较小是八叉树查询很快；当查询框接近整个空间时，八叉树查询会慢很多；而蛮力查询所需时间一直很稳定，甚至当查询空间小于一般时，会降低

（1）查询框大小为

//main.cpp的63、64行qmin = Vec3(-.05,-.05,-.05);qmax = Vec3(.05,.05,.05);

编译并运行，得到八叉树查询用时：0.00003 sec.蛮力查询用时0.00711sec.

meng@meng:~/SimpleOctree$ make
g++ main.cpp -O3 -o octree
meng@meng:~/SimpleOctree$ ./octree
Created 1000000 points
Inserted points to octree
testNaive found 150 points in 0.00711 sec.
testOctree found 150 points in 0.00003 sec.

（2）查询框大小为

 qmin = Vec3(-.5,-.5,-.5);qmax = Vec3(.5,.5,.5);

删除可执行文件octree，重新编译

Created 1000000 points
Inserted points to octree
testNaive found 150 points in 0.00726 sec.
testOctree found 150 points in 0.00004 sec.

（3）查询框大小为

 qmin = Vec3(-.9,-.9,-.9);qmax = Vec3(.9,.9,.9);

meng@meng:~/SimpleOctree$ ./octree
Created 1000000 points
Inserted points to octree
testNaive found 729239 points in 0.00554 sec.
testOctree found 729239 points in 0.09651 sec.

更多对比自己验证啦

代码来源：https://github.com/brandonpelfrey/SimpleOctree

@meng