今天依然在放假中，在此将以前在学校写的四叉树的东西拿出来和大家分享。
四叉树索引的基本思想是将地理空间递归划分为不同层次的树结构。它将已知范围的空间等分成四个相等的子空间，如此递归下去，直至树的层次达到一定深度或者满足某种要求后停止分割。四叉树的结构比较简单，并且当空间数据对象分布比较均匀时，具有比较高的空间数据插入和查询效率，因此四叉树是gis中常用的空间索引之一。常规四叉树的结构如图所示，地理空间对象都存储在叶子节点上，中间节点以及根节点不存储地理空间对象。
四叉树示意图
四叉树对于区域查询，效率比较高。但如果空间对象分布不均匀，随着地理空间对象的不断插入，四叉树的层次会不断地加深，将形成一棵严重不平衡的四叉树，那么每次查询的深度将大大的增多，从而导致查询效率的急剧下降。
本节将介绍一种改进的四叉树索引结构。四叉树结构是自顶向下逐步划分的一种树状的层次结构。传统的四叉树索引存在着以下几个缺点：
(1)空间实体只能存储在叶子节点中，中间节点以及根节点不能存储空间实体信息，随着空间对象的不断插入，最终会导致四叉树树的层次比较深，在进行空间数据窗口查询的时候效率会比较低下。
(2)同一个地理实体在四叉树的分裂过程中极有可能存储在多个节点中，这样就导致了索引存储空间的浪费。
(3)由于地理空间对象可能分布不均衡，这样会导致常规四叉树生成一棵极为不平衡的树，这样也会造成树结构的不平衡以及存储空间的浪费。
相应的改进方法，将地理实体信息存储在完全包含它的最小矩形节点中，不存储在它的父节点中，每个地理实体只在树中存储一次，避免存储空间的浪费。首先生成满四叉树，避免在地理实体插入时需要重新分配内存，加快插入的速度，最后将空的节点所占内存空间释放掉。改进后的四叉树结构如下图所示。四叉树的深度一般取经验值4-7之间为最佳。
图改进的四叉树结构
为了维护空间索引与对存储在文件或数据库中的空间数据的一致性，作者设计了如下的数据结构支持四叉树的操作。
(1)四分区域标识
分别定义了一个平面区域的四个子区域索引号，右上为第一象限0，左上为第二象限1，左下为第三象限2，右下为第四象限3。
typedef enum
{
ur = 0,// ur第一象限
ul = 1, // ul为第二象限
ll = 2, // ll为第三象限
lr = 3 // lr为第四象限
}quadrantenum;
(2)空间对象数据结构
空间对象数据结构是对地理空间对象的近似，在空间索引中，相当一部分都是采用mbr作为近似。
/*空间对象mbr信息*/
typedef struct shpmbrinfo
{
int nid; //空间对象id号
maprect box; //空间对象mbr范围坐标
}shpmbrinfo;
nid是空间对象的标识号，box是空间对象的最小外包矩形（mbr）。
(3)四叉树节点数据结构
四叉树节点是四叉树结构的主要组成部分，主要用于存储空间对象的标识号和mbr，也是四叉树算法操作的主要部分。
/*四叉树节点类型结构*/
typedef struct quadnode
{
maprect box; //节点所代表的矩形区域
int nshpcount; //节点所包含的所有空间对象个数
shpmbrinfo* pshapeobj; //空间对象指针数组
int nchildcount; //子节点个数
quadnode *children[4]; //指向节点的四个孩子
}quadnode;
box是代表四叉树对应区域的最小外包矩形，上一层的节点的最小外包矩形包含下一层最小外包矩形区域；nshpcount代表本节点包含的空间对象的个数；pshapeobj代表指向空间对象存储地址的首地址，同一个节点的空间对象在内存中连续存储；nchildcount代表节点拥有的子节点的数目；children是指向孩子节点指针的数组。
上述理论部分都都讲的差不多了，下面就贴上我的c语言实现版本代码。
头文件如下：
#ifndef __quadtree_h_59cae94a_e937_42ad_aa27_794e467715bb__
#define __quadtree_h_59cae94a_e937_42ad_aa27_794e467715bb__
/* 一个矩形区域的象限划分：:
ul(1) | ur(0)
----------|-----------
ll(2) | lr(3)
以下对该象限类型的枚举
*/
typedef enum
{
ur = 0,
ul = 1,
ll = 2,
lr = 3
}quadrantenum;
/*空间对象mbr信息*/
typedef struct shpmbrinfo
{
int nid; //空间对象id号
maprect box; //空间对象mbr范围坐标
}shpmbrinfo;
/* 四叉树节点类型结构 */
typedef struct quadnode
{
maprect box; //节点所代表的矩形区域
int nshpcount; //节点所包含的所有空间对象个数
shpmbrinfo* pshapeobj; //空间对象指针数组
int nchildcount; //子节点个数
quadnode *children[4]; //指向节点的四个孩子
}quadnode;
/* 四叉树类型结构 */
typedef struct quadtree_t
{
quadnode *root;
int depth; // 四叉树的深度
}quadtree;
//初始化四叉树节点
quadnode *initquadnode();
//层次创建四叉树方法（满四叉树）
void createquadtree(int depth,geolayer *polayer,quadtree* pquadtree);
//创建各个分支
void createquadbranch(int depth,maprect &rect,quadnode** node);
//构建四叉树空间索引
void buildquadtree(geolayer*polayer,quadtree* pquadtree);
//四叉树索引查询(矩形查询)
void searchquadtree(quadnode* node,maprect &queryrect,vector& itemsearched);
//四叉树索引查询(矩形查询)并行查询
void searchquadtreepara(vector resnodes,maprect &queryrect,vector& itemsearched);
//四叉树的查询（点查询）
void ptsearchqtree(quadnode* node,double cx,double cy,vector& itemsearched);
//将指定的空间对象插入到四叉树中
void insert(long key,maprect &itemrect,quadnode* pnode);
//将指定的空间对象插入到四叉树中
void insertquad(long key,maprect &itemrect,quadnode* pnode);
//将指定的空间对象插入到四叉树中
void insertquad2(long key,maprect &itemrect,quadnode* pnode);
//判断一个节点是否是叶子节点
bool isquadleaf(quadnode* node);
//删除多余的节点
bool delfalsenode(quadnode* node);
//四叉树遍历(所有要素)
void traversalquadtree(quadnode* quadtree,vector& resvec);
//四叉树遍历（所有节点）
void traversalquadtree(quadnode* quadtree,vector& arrnode);
//释放树的内存空间
void releasequadtree(quadnode** quadtree);
//计算四叉树所占的字节的大小
long calbytequadtree(quadnode* quadtree,long& nsize);
#endif 源文件如下：
#include "quadtree.h"
quadnode *initquadnode()
{
quadnode *node = new quadnode;
node->box.maxx = 0;
node->box.maxy = 0;
node->box.minx = 0;
node->box.miny = 0;
for (int i = 0; i < 4; i )
{
node->children[i] = null;
}
node->nchildcount = 0;
node->nshpcount = 0;
node->pshapeobj = null;
return node;
}
void createquadtree(int depth,geolayer *polayer,quadtree* pquadtree)
{
pquadtree->depth = depth;
geoenvelope env; //整个图层的mbr
polayer->getextent(&env);
maprect rect;
rect.minx = env.minx;
rect.miny = env.miny;
rect.maxx = env.maxx;
rect.maxy = env.maxy;
//创建各个分支
createquadbranch(depth,rect,&(pquadtree->root));
int ncount = polayer->getfeaturecount();
geofeature **pfeatureclass = new geofeature*[ncount];
for (int i = 0; i < polayer->getfeaturecount(); i )
{
pfeatureclass[i] = polayer->getfeature(i);
}
//插入各个要素
geoenvelope envobj; //空间对象的mbr
//#pragma omp parallel for
for (int i = 0; i < ncount; i )
{
pfeatureclass[i]->getgeometry()->getenvelope(&envobj);
rect.minx = envobj.minx;
rect.miny = envobj.miny;
rect.maxx = envobj.maxx;
rect.maxy = envobj.maxy;
insertquad(i,rect,pquadtree->root);
}
//delfalsenode(pquadtree->root);
}
void createquadbranch(int depth,maprect &rect,quadnode** node)
{
if (depth != 0)
{
*node = initquadnode(); //创建树根
quadnode *pnode = *node;
pnode->box = rect;
pnode->nchildcount = 4;
maprect boxs[4];
pnode->box.split(boxs,boxs 1,boxs 2,boxs 3);
for (int i = 0; i < 4; i )
{
//创建四个节点并插入相应的mbr
pnode->children[i] = initquadnode();
pnode->children[i]->box = boxs[i];
createquadbranch(depth-1,boxs[i],&(pnode->children[i]));
}
}
}
void buildquadtree(geolayer *polayer,quadtree* pquadtree)
{
assert(polayer);
geoenvelope env; //整个图层的mbr
polayer->getextent(&env);
pquadtree->root = initquadnode();
quadnode* rootnode = pquadtree->root;
rootnode->box.minx = env.minx;
rootnode->box.miny = env.miny;
rootnode->box.maxx = env.maxx;
rootnode->box.maxy = env.maxy;
//设置树的深度( 根据等比数列的求和公式)
//pquadtree->depth = log(polayer->getfeaturecount()*3/8.0 1)/log(4.0);
int ncount = polayer->getfeaturecount();
maprect rect;
geoenvelope envobj; //空间对象的mbr
for (int i = 0; i < ncount; i )
{
polayer->getfeature(i)->getgeometry()->getenvelope(&envobj);
rect.minx = envobj.minx;
rect.miny = envobj.miny;
rect.maxx = envobj.maxx;
rect.maxy = envobj.maxy;
insertquad2(i,rect,rootnode);
}
delfalsenode(pquadtree->root);
}
void searchquadtree(quadnode* node,maprect &queryrect,vector& itemsearched)
{
assert(node);
//int corenum = omp_get_num_procs();
//vector * presarr = new vector[corenum];
if (null != node)
{
for (int i = 0; i < node->nshpcount; i )
{
if (queryrect.contains(node->pshapeobj[i].box)
|| queryrect.intersects(node->pshapeobj[i].box))
{
itemsearched.push_back(node->pshapeobj[i].nid);
}
}
//并行搜索四个孩子节点
/*#pragma omp parallel sections
{
#pragma omp section
if ((node->children[0] != null) &&
(node->children[0]->box.contains(queryrect)
|| node->children[0]->box.intersects(queryrect)))
{
int tid = omp_get_thread_num();
searchquadtree(node->children[0],queryrect,presarr[tid]);
}
#pragma omp section
if ((node->children[1] != null) &&
(node->children[1]->box.contains(queryrect)
|| node->children[1]->box.intersects(queryrect)))
{
int tid = omp_get_thread_num();
searchquadtree(node->children[1],queryrect,presarr[tid]);
}
#pragma omp section
if ((node->children[2] != null) &&
(node->children[2]->box.contains(queryrect)
|| node->children[2]->box.intersects(queryrect)))
{
int tid = omp_get_thread_num();
searchquadtree(node->children[2],queryrect,presarr[tid]);
}
#pragma omp section
if ((node->children[3] != null) &&
(node->children[3]->box.contains(queryrect)
|| node->children[3]->box.intersects(queryrect)))
{
int tid = omp_get_thread_num();
searchquadtree(node->children[3],queryrect,presarr[tid]);
}
}*/
for (int i = 0; i < 4; i )
{
if ((node->children[i] != null) &&
(node->children[i]->box.contains(queryrect)
|| node->children[i]->box.intersects(queryrect)))
{
searchquadtree(node->children[i],queryrect,itemsearched);
//node = node->children[i]; //非递归
}
}
}
/*for (int i = 0 ; i < corenum; i )
{
itemsearched.insert(itemsearched.end(),presarr[i].begin(),presarr[i].end());
}*/
}
void searchquadtreepara(vector resnodes,maprect &queryrect,vector& itemsearched)
{
int corenum = omp_get_num_procs();
omp_set_num_threads(corenum);
vector* searcharrs = new vector[corenum];
for (int i = 0; i < corenum; i )
{
searcharrs[i].clear();
}
#pragma omp parallel for
for (int i = 0; i < resnodes.size(); i )
{
int tid = omp_get_thread_num();
for (int j = 0; j < resnodes[i]->nshpcount; j )
{
if (queryrect.contains(resnodes[i]->pshapeobj[j].box)
|| queryrect.intersects(resnodes[i]->pshapeobj[j].box))
{
searcharrs[tid].push_back(resnodes[i]->pshapeobj[j].nid);
}
}
}
for (int i = 0; i < corenum; i )
{
itemsearched.insert(itemsearched.end(),
searcharrs[i].begin(),searcharrs[i].end());
}
delete [] searcharrs;
searcharrs = null;
}
void ptsearchqtree(quadnode* node,double cx,double cy,vector& itemsearched)
{
assert(node);
if (node->nshpcount >0) //节点
{
for (int i = 0; i < node->nshpcount; i )
{
if (node->pshapeobj[i].box.ispointinrect(cx,cy))
{
itemsearched.push_back(node->pshapeobj[i].nid);
}
}
}
else if (node->nchildcount >0) //节点
{
for (int i = 0; i < 4; i )
{
if (node->children[i]->box.ispointinrect(cx,cy))
{
ptsearchqtree(node->children[i],cx,cy,itemsearched);
}
}
}
//找出重复元素的位置
sort(itemsearched.begin(),itemsearched.end()); //先排序,默认升序
vector::iterator unique_iter =
unique(itemsearched.begin(),itemsearched.end());
itemsearched.erase(unique_iter,itemsearched.end());
}
void insert(long key, maprect &itemrect,quadnode* pnode)
{
quadnode *node = pnode; //保留根节点副本
shpmbrinfo pshpinfo;
//节点有孩子
if (0 < node->nchildcount)
{
for (int i = 0; i < 4; i )
{
//如果包含或相交，则将节点插入到此节点
if (node->children[i]->box.contains(itemrect)
|| node->children[i]->box.intersects(itemrect))
{
//node = node->children[i];
insert(key,itemrect,node->children[i]);
}
}
}
//如果当前节点存在一个子节点时
else if (1 == node->nshpcount)
{
maprect boxs[4];
node->box.split(boxs,boxs 1,boxs 2,boxs 3);
//创建四个节点并插入相应的mbr
node->children[ur] = initquadnode();
node->children[ul] = initquadnode();
node->children[ll] = initquadnode();
node->children[lr] = initquadnode();
node->children[ur]->box = boxs[0];
node->children[ul]->box = boxs[1];
node->children[ll]->box = boxs[2];
node->children[lr]->box = boxs[3];
node->nchildcount = 4;
for (int i = 0; i < 4; i )
{
//将当前节点中的要素移动到相应的子节点中
for (int j = 0; j < node->nshpcount; j )
{
if (node->children[i]->box.contains(node->pshapeobj[j].box)
|| node->children[i]->box.intersects(node->pshapeobj[j].box))
{
node->children[i]->nshpcount = 1;
node->children[i]->pshapeobj =
(shpmbrinfo*)malloc(node->children[i]->nshpcount*sizeof(shpmbrinfo));
memcpy(node->children[i]->pshapeobj,&(node->pshapeobj[j]),sizeof(shpmbrinfo));
free(node->pshapeobj);
node->pshapeobj = null;
node->nshpcount = 0;
}
}
}
for (int i = 0; i < 4; i )
{
//如果包含或相交，则将节点插入到此节点
if (node->children[i]->box.contains(itemrect)
|| node->children[i]->box.intersects(itemrect))
{
if (node->children[i]->nshpcount == 0) //如果之前没有节点
{
node->children[i]->nshpcount = 1;
node->pshapeobj =
(shpmbrinfo*)malloc(sizeof(shpmbrinfo)*node->children[i]->nshpcount);
}
else if (node->children[i]->nshpcount > 0)
{
node->children[i]->nshpcount = 1;
node->children[i]->pshapeobj =
(shpmbrinfo *)realloc(node->children[i]->pshapeobj,
sizeof(shpmbrinfo)*node->children[i]->nshpcount);
}
pshpinfo.box = itemrect;
pshpinfo.nid = key;
memcpy(node->children[i]->pshapeobj,
&pshpinfo,sizeof(shpmbrinfo));
}
}
}
//当前节点没有空间对象
else if (0 == node->nshpcount)
{
node->nshpcount = 1;
node->pshapeobj =
(shpmbrinfo*)malloc(sizeof(shpmbrinfo)*node->nshpcount);
pshpinfo.box = itemrect;
pshpinfo.nid = key;
memcpy(node->pshapeobj,&pshpinfo,sizeof(shpmbrinfo));
}
}
void insertquad(long key,maprect &itemrect,quadnode* pnode)
{
assert(pnode != null);
if (!isquadleaf(pnode)) //非叶子节点
{
int ncorver = 0; //跨越的子节点个数
int iindex = -1; //被哪个子节点完全包含的索引号
for (int i = 0; i < 4; i )
{
if (pnode->children[i]->box.contains(itemrect)
&& pnode->box.contains(itemrect))
{
ncorver = 1;
iindex = i;
}
}
//如果被某一个子节点包含，则进入该子节点
if (/*pnode->box.contains(itemrect) ||
pnode->box.intersects(itemrect)*/1 <= ncorver)
{
insertquad(key,itemrect,pnode->children[iindex]);
}
//如果跨越了多个子节点，直接放在这个节点中
else if (ncorver == 0)
{
if (pnode->nshpcount == 0) //如果之前没有节点
{
pnode->nshpcount = 1;
pnode->pshapeobj =
(shpmbrinfo*)malloc(sizeof(shpmbrinfo)*pnode->nshpcount);
}
else
{
pnode->nshpcount = 1;
pnode->pshapeobj =
(shpmbrinfo *)realloc(pnode->pshapeobj,sizeof(shpmbrinfo)*pnode->nshpcount);
}
shpmbrinfo pshpinfo;
pshpinfo.box = itemrect;
pshpinfo.nid = key;
memcpy(pnode->pshapeobj pnode->nshpcount-1,&pshpinfo,sizeof(shpmbrinfo));
}
}
//如果是叶子节点，直接放进去
else if (isquadleaf(pnode))
{
if (pnode->nshpcount == 0) //如果之前没有节点
{
pnode->nshpcount = 1;
pnode->pshapeobj =
(shpmbrinfo*)malloc(sizeof(shpmbrinfo)*pnode->nshpcount);
}
else
{
pnode->nshpcount = 1;
pnode->pshapeobj =
(shpmbrinfo *)realloc(pnode->pshapeobj,sizeof(shpmbrinfo)*pnode->nshpcount);
}
shpmbrinfo pshpinfo;
pshpinfo.box = itemrect;
pshpinfo.nid = key;
memcpy(pnode->pshapeobj pnode->nshpcount-1,&pshpinfo,sizeof(shpmbrinfo));
}
}
void insertquad2(long key,maprect &itemrect,quadnode* pnode)
{
quadnode *node = pnode; //保留根节点副本
shpmbrinfo pshpinfo;
//节点有孩子
if (0 < node->nchildcount)
{
for (int i = 0; i < 4; i )
{
//如果包含或相交，则将节点插入到此节点
if (node->children[i]->box.contains(itemrect)
|| node->children[i]->box.intersects(itemrect))
{
//node = node->children[i];
insert(key,itemrect,node->children[i]);
}
}
}
//如果当前节点存在一个子节点时
else if (0 == node->nchildcount)
{
maprect boxs[4];
node->box.split(boxs,boxs 1,boxs 2,boxs 3);
int cnt = -1;
for (int i = 0; i < 4; i )
{
//如果包含或相交，则将节点插入到此节点
if (boxs[i].contains(itemrect))
{
cnt = i;
}
}
//如果有一个矩形包含此对象，则创建四个孩子节点
if (cnt > -1)
{
for (int i = 0; i < 4; i )
{
//创建四个节点并插入相应的mbr
node->children[i] = initquadnode();
node->children[i]->box = boxs[i];
}
node->nchildcount = 4;
insertquad2(key,itemrect,node->children[cnt]); //递归
}
//如果都不包含，则直接将对象插入此节点
if (cnt == -1)
{
if (node->nshpcount == 0) //如果之前没有节点
{
node->nshpcount = 1;
node->pshapeobj =
(shpmbrinfo*)malloc(sizeof(shpmbrinfo)*node->nshpcount);
}
else if (node->nshpcount > 0)
{
node->nshpcount = 1;
node->pshapeobj =
(shpmbrinfo *)realloc(node->pshapeobj,
sizeof(shpmbrinfo)*node->nshpcount);
}
pshpinfo.box = itemrect;
pshpinfo.nid = key;
memcpy(node->pshapeobj,
&pshpinfo,sizeof(shpmbrinfo));
}
}
//当前节点没有空间对象
/*else if (0 == node->nshpcount)
{
node->nshpcount = 1;
node->pshapeobj =
(shpmbrinfo*)malloc(sizeof(shpmbrinfo)*node->nshpcount);
pshpinfo.box = itemrect;
pshpinfo.nid = key;
memcpy(node->pshapeobj,&pshpinfo,sizeof(shpmbrinfo));
}*/
}
bool isquadleaf(quadnode* node)
{
if (null == node)
{
return 1;
}
for (int i = 0; i < 4; i )
{
if (node->children[i] != null)
{
return 0;
}
}
return 1;
}
bool delfalsenode(quadnode* node)
{
//如果没有子节点且没有要素
if (node->nchildcount ==0 && node->nshpcount == 0)
{
releasequadtree(&node);
}
//如果有子节点
else if (node->nchildcount > 0)
{
for (int i = 0; i < 4; i )
{
delfalsenode(node->children[i]);
}
}
return 1;
}
void traversalquadtree(quadnode* quadtree,vector& resvec)
{
quadnode *node = quadtree;
int i = 0;
if (null != node)
{
//将本节点中的空间对象存储数组中
for (i = 0; i < node->nshpcount; i )
{
resvec.push_back((node->pshapeobj i)->nid);
}
//遍历孩子节点
for (i = 0; i < node->nchildcount; i )
{
if (node->children[i] != null)
{
traversalquadtree(node->children[i],resvec);
}
}
}
}
void traversalquadtree(quadnode* quadtree,vector& arrnode)
{
deque nodequeue;
if (quadtree != null)
{
nodequeue.push_back(quadtree);
while (!nodequeue.empty())
{
quadnode* queuehead = nodequeue.at(0); //取队列头结点
arrnode.push_back(queuehead);
nodequeue.pop_front();
for (int i = 0; i < 4; i )
{
if (queuehead->children[i] != null)
{
nodequeue.push_back(queuehead->children[i]);
}
}
}
}
}
void releasequadtree(quadnode** quadtree)
{
int i = 0;
quadnode* node = *quadtree;
if (null == node)
{
return;
}
else
{
for (i = 0; i < 4; i )
{
releasequadtree(&node->children[i]);
}
free(node);
node = null;
}
node = null;
}
long calbytequadtree(quadnode* quadtree,long& nsize)
{
if (quadtree != null)
{
nsize = sizeof(quadnode) quadtree->nchildcount*sizeof(shpmbrinfo);
for (int i = 0; i < 4; i )
{
if (quadtree->children[i] != null)
{
nsize = calbytequadtree(quadtree->children[i],nsize);
}
}
}
return 1;
} 代码有点长，有需要的朋友可以借鉴并自己优化。