地理空間索引實作：z 曲線、希爾伯特曲線、四叉樹, 最鄰近幾何特征查詢、範圍查詢

我的GIS/CS學習筆記：https://github.com/yunwei37/ZJU-CS-GIS-ClassNotes

<一個浙江大學大學生的計算機、地理資訊科學知識庫 >

詳細代碼可在其中檢視

空間索引

在談論空間索引之前，我們必須了解資料索引的概念：索引是為了提高資料集的檢索效率。打個比喻，一本書的目錄就是這本書的内容的“索引”，我們檢視感興趣的内容前，通過檢視書的目錄去快速查找對應的内容，而不是一字一句地找我們感興趣的内容；就像這樣，事先建構的索引可以有效地加速查詢的速度。

然而，和一般的資料相比，有效地查詢地理空間資料是相當大的挑戰，因為資料是二維的（有時候甚至更高），不能用如傳統的B+樹這樣标準的索引技術來加速查詢位置相關的資料。是以，我們就引入了空間索引技術來解決這個問題。

空間索引的定義：

依據空間實體的位置和形狀或空間實體之間的某種空間關系，按一定順序排列的一種資料結構，其中包含空間實體的概要資訊，如對象的辨別，最小邊界矩形及指向空間實體資料的指針

常見的空間索引技術有網格索引、四叉樹索引，空間填充曲線索引，以及最用于地理空間資料庫的R樹索引以及相關變體等等。

網格索引

網格索引的基本思想是将研究區域用橫豎線劃分大小相等和不等的網格，每個網格可視為一個桶(bucket)，建構時記錄落入每一個網格區域内的空間實體編号。

進行空間查詢時，先計算出查詢對象所在網格，再在該網格中快速查

詢所選空間實體

網格索引優點：簡單，易于實作，具有良好的可擴充性；
網格索引缺點：網格大小影響網格索引檢索性能

理想的情況下，網格大小是使網格索引記錄不至于過多，同時每個網格内的要素個數的均值與最大值盡可能地少。如要獲得較好的網格劃分，可以根據使用者的多次試驗來獲得經驗最佳值，也可以通過建立地理要素的大小和空間分布等特征值來定量确定網格大小。

網格索引的實作這裡暫時沒有涉及。

空間填充曲線索引

常用的空間索引曲線有z曲線、希爾伯特曲線，其目的是在空間網格的基礎上降低空間次元，以便于在順序讀取的磁盤上存取資訊。

空間填充曲線(space-filling curve)是一條連續曲線，自身沒有任何交叉；通過通路所有單元格來填充包含均勻網格的四邊形。

z曲線

地理空間索引實作：z 曲線、希爾伯特曲線、四叉樹, 最鄰近幾何特征查詢、範圍查詢
希爾伯特曲線

地理空間索引實作：z 曲線、希爾伯特曲線、四叉樹, 最鄰近幾何特征查詢、範圍查詢

Z曲線和Hilbert曲線共同特點：

填充曲線值臨近的網格，其空間位置通常也相對臨近；
任何一種空間排列都不能完全保證二維資料空間關系的維護（編号相鄰，空間位置可能很遠）

不同點：

Hilbert曲線的資料聚集特性更優，Z曲線的資料聚集特性較差
Hilbert曲線的映射過程較複雜，Z曲線的映射過程較簡單

z曲線實作：

Z-curve曲線的二維坐标與Z值的互相轉換：

基于bit-shuffling思想，實作二維坐标(coor)與Z值(value)的互相轉換。order為Z-Curve的階數，order為4時，平面網格大小為(2^4, 2^4)，即16x16。

/*
 * zorder - Calculate the index number (Z-value) of the point (represented by coordinate) in the Z-curve with given order
 */
void zorder(int order, int& value, int coor[2])
{
	// Calculate the z-order by bit shuffling
	value = 0;
	for (int i = 0; i < 2; ++i) {
		for (int j = 0; j < order; ++j) {
			// Task 1.1 zorder，修改以下代碼
			int mask = 1<<j;
			// Check whether the value in the position is 1
			if (coor[1-i] & mask)
				// Do bit shuffling
				value |= 1 << (i + j*2);
		}
	}
}

/*
 * izorder - Calculate the coordinate of the point with given Z-value and order
 */
void izorder(int order, int value, int coor[2])
{
	// Initialize the coordinate to zeros
	for (int i = 0; i < 2; ++i)
		coor[i] = 0;
	// Task 1.2 izoder
	// Write your code here
	for (int i = 0; i < 2; ++i) {
		for (int j = 0; j < order; ++j) {
			// Task 1.1 zorder，修改以下代碼
			int mask = 1 << j*2+i;
			// Check whether the value in the position is 1
			if (value & mask)
				// Do bit shuffling
				coor[1-i] |= 1 << j;
		}
	}
}

希爾伯特曲線實作

Hilbert Curve的二維坐标與H值的互相轉換：

基于linear-quadtrees遞歸填充思想，實作二維坐标(coor)到H值(value)的互相轉換。order為Hilbert Curve的階數，order為4時，平面網格大小為(2^4, 2^4)，即16x16。

/*
 * horder - Calculate the index number (H-value) of the point (represented by coordinate) in Hilbert curve with given order
 */
void horder(int order, int& value, int coor[2])
{
	// order = 1
	int num = int(pow(2, 1));
	int * hcurve = new int[num * num];
    hcurve[0] = 1;
    hcurve[1] = 2;
    hcurve[2] = 0;
    hcurve[3] = 3;

    for(int i = 2; i <= order; ++i) {
		int add = (int)pow(2, 2 * i - 2);		// the number of values in order - 1
		int blockLen = (int)pow(2, i - 1);
		int * temp = hcurve;

        num = int(pow(2, i));
		hcurve = new int[num * num];

		for (int j = 0; j < blockLen; ++j) {
			for (int k = 0; k < blockLen; ++k) {
				// Task 2.1 horder，修改以下四行代碼
				hcurve[j * num + k] = temp[j * blockLen + k] + add;
				hcurve[j * num + k + blockLen] = temp[j * blockLen + k] + add*2;
				hcurve[(j + blockLen) * num + k] = temp[(blockLen - k-1) * blockLen + (blockLen - j-1)] + 0;
				hcurve[(j + blockLen) * num + k + blockLen] = temp[k * blockLen + j] + add * 3;
			}
		}

		delete temp;
    }

	// Task 2.1 horder，修改以下一行代碼
	value = hcurve[num*(num - coor[1] - 1)+coor[0]];
	delete[] hcurve;
}

/*
 * ihorder - Calculate the coordinate of the point with given H-value and order
 */
void ihorder(int order, int value, int coor[2])
{
	// order = 1
	int num = int(pow(2, 1));
	int * hcurve = new int[num * num];
	hcurve[0] = 1;
	hcurve[1] = 2;
	hcurve[2] = 0;
	hcurve[3] = 3;
	
	for (int i = 2; i <= order; ++i) {
		int add = (int)pow(2, 2 * i - 2);		// the number of values in order - 1
		int blockLen = (int)pow(2, i - 1);
		int* temp = hcurve;

		num = int(pow(2, i));
		hcurve = new int[num * num];

		for (int j = 0; j < blockLen; ++j) {
			for (int k = 0; k < blockLen; ++k) {
				// Task 2.1 horder，修改以下四行代碼
				hcurve[j * num + k] = temp[j * blockLen + k] + add;
				hcurve[j * num + k + blockLen] = temp[j * blockLen + k] + add * 2;
				hcurve[(j + blockLen) * num + k] = temp[(blockLen - k - 1) * blockLen + (blockLen - j - 1)] + 0;
				hcurve[(j + blockLen) * num + k + blockLen] = temp[k * blockLen + j] + add * 3;
			}
		}

		delete temp;
	}
	for (int i = 0; i < 2; ++i)
		coor[i] = 0;
	int i = 0;
	for (; i < num * num; ++i) {
		if (value == hcurve[i])
			break;
	}
	coor[0] = i % num;
	coor[1] = num - 1 - i / num;
	delete[] hcurve;
}

四叉樹索引

四叉樹索引是在網格索引的思想基礎上，為了實作要素真正被網格分割，同時保證桶内要素不超過一個量而提出的一種空間索引方法。

構造方法：

首先将整個資料空間分割成為四個相等的矩陣，分别對應西北(NW)，東北(NE)，西南(SW)，東南(SE)四個象限；
若每個象限内包含的要素不超過給定的桶量則停止，否則對超過桶量的矩形再按照同樣的方法進行劃分，直到桶量滿足要求或者不再減少為止，最終形成一顆有層次的四叉樹。

地理空間索引實作：z 曲線、希爾伯特曲線、四叉樹, 最鄰近幾何特征查詢、範圍查詢
四叉樹優缺點：

與網格索引相比，四叉樹在一定程度上實作了地理要素真正被網格分割，保證了桶内要素不超過某個量，提高了檢索效率；
對于海量資料，四叉樹的深度會很深，影響查詢效率
可擴充性不如網格索引：當擴大區域時，需要重新劃分空間區域，重建四叉樹，當增加或删除一個對象，可能導緻深度加一或減一，葉節點也有可能重新定位。

四叉樹索引建構：

四叉樹建立輸入一組幾何特征，将節點分裂為四個子節點，每個特征加到包圍盒重疊的子節點中（即一個特征可能在多個節點中），删除目前節點的幾何特征記錄（即所有特征隻存儲在葉節點中），如果子節點的幾何特征個數大于capacity，遞歸分裂子節點。

/*
 * QuadNode
 */
void QuadNode::split(size_t capacity)
{
	for (int i = 0; i < 4; ++i) {
		delete []nodes[i];
		nodes[i] = NULL;
	}
	
	if (features.size() > capacity) {
		double midx = (bbox.getMinX() + bbox.getMaxX()) / 2.0;
		double midy = (bbox.getMinY() + bbox.getMaxY()) / 2.0;
		Envelope e0(bbox.getMinX(), midx, midy, bbox.getMaxY());
		nodes[0] = new QuadNode(e0);
		Envelope e1(midx, bbox.getMaxX(), midy, bbox.getMaxY());
		nodes[1] = new QuadNode(e1);
		Envelope e2(midx, bbox.getMaxX(), bbox.getMinY(), midy);
		nodes[2] = new QuadNode(e2);
		Envelope e3(bbox.getMinX(), midx, bbox.getMinY(), midy);
		nodes[3] = new QuadNode(e3);

		for (Feature f : features) {
			for (int i = 0; i < 4; ++i) {
				if (nodes[i]->getEnvelope().intersect(f.getEnvelope())) {
					nodes[i]->add(f);
				}
			}
		}

		for (int i = 0; i < 4; ++i) {
			nodes[i]->split(capacity);
		}
		features.clear();
	}
}

/*
 * QuadTree
 */
bool QuadTree::constructQuadTree(vector<Feature>& features)
{
	if (features.empty())
		return false;

	// Task 5.1 construction
	// Write your code here

	Envelope e(DBL_MAX, -DBL_MAX, DBL_MAX, -DBL_MAX);
	for (Feature f : features) {
		e = e.unionEnvelope(f.getEnvelope());
	}
	root = new QuadNode(e);

	root->add(features);

	root->split(capacity);

	bbox = e; // 注意此行代碼需要更新為features的包圍盒，或根節點的包圍盒
	//bbox = Envelope(-74.1, -73.8, 40.6, 40.8);
	return true;
}

效果：

四叉樹上的最鄰近幾何特征查詢、範圍查詢：

最鄰近幾何特征查詢：

最鄰近幾何特征查詢(K-NN)輸入查詢點(x, y)，傳回與該點最鄰近的幾何特征，存儲在feature。首先，通過pointInLeafNode查詢點(x, y)所在的葉節點，計算查詢點(x, y)與該葉節點内的幾何特征包圍盒的最大距離的最小值minDist，即通過包圍盒而非原始幾何加速最小距離計算；然後，構造查詢區域 (x – minDist, x + minDist, y – minDist, y + minDist)，查詢幾何特征的包圍盒與該區域相交的幾何特征（filter），再查詢與查詢點(x, y)距離最近的幾何特征（refine）。

範圍查詢：

地理空間索引實作：z 曲線、希爾伯特曲線、四叉樹, 最鄰近幾何特征查詢、範圍查詢

空間索引

網格索引

空間填充曲線索引

z曲線實作：

希爾伯特曲線實作

四叉樹索引

四叉樹索引建構：

四叉樹上的最鄰近幾何特征查詢、範圍查詢：

具體代碼可在我的github項目中檢視：

繼續閱讀

swmm與lisflood-fp源碼如何一起編譯 CMake指令

CSMA/CD1． CSMA/CD的概述2． CSMA 的工作原理3． CSMA/CD控制規程及特點4． CSMA/CD協定5． CSMA/CD的優點6．結束語

Windows下VS開發環境環境安裝工程項目設定關于Debug和Release的提示

一文看懂字元串的加減乘除

極大似然法(ML)與最大期望法(EM)

C++ 第十五周報告1--《冒泡法排序》

C++實作簡單順序表

筆試面試題目：滑動視窗(二)

資料結構與算法（27）——排序（二）

Dijkstra--簡易版（最短路徑）

C經典書籍筆記——C陷阱與缺陷②(文法陷阱之優先級)一、錯誤案列二、優先級規律

GitHub連夜封殺！這份阿裡 10W 字内部 Java 字面試手冊到底有多強？

線性表之順序表的實作

C++判斷素數、求最大公約數代碼判斷一個數是否為素數求兩個數的最大公約數

SequoiaDB巨杉資料庫C++驅動概述

hdu7108哈希