轉自原文 四叉樹空間索引原理及其實現
四叉樹索引的基本思想是將地理空間遞歸划分為不同層次的樹結構。它將已知范圍的空間等分成四個相等的子空間，如此遞歸下去，直至樹的層次達到一定深度或者滿足某種要求后停止分割。四叉樹的結構比較簡單，並且當空間數據對象分布比較均勻時，具有比較高的空間數據插入和查詢效率，因此四叉樹是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<int>& ItemSearched); //四叉樹索引查詢(矩形查詢)並行查詢 void SearchQuadTreePara(vector<QuadNode*> resNodes,MapRect &queryRect,vector<int>& ItemSearched); //四叉樹的查詢（點查詢） void PtSearchQTree(QuadNode* node,double cx,double cy,vector<int>& 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<int>& resVec); //四叉樹遍歷（所有節點） void TraversalQuadTree(QuadNode* quadTree,vector<QuadNode*>& 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<int>& ItemSearched) { assert(node); //int coreNum = omp_get_num_procs(); //vector<int> * pResArr = new vector<int>[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<QuadNode*> resNodes,MapRect &queryRect,vector<int>& ItemSearched) { int coreNum = omp_get_num_procs(); omp_set_num_threads(coreNum); vector<int>* searchArrs = new vector<int>[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<int>& 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<int>::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<int>& 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<QuadNode*>& arrNode) { deque<QuadNode*> 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; }