地理空间索引实现：z 曲线、希尔伯特曲线、四叉树, 最邻近几何特征查询、范围查询

依据空间实体的位置和形状或空间实体之间的某种空间关系，按一定顺序排列的一种数据结构，其中包含空间实体的概要信息，如对象的标识，最小边界矩形及指向空间实体数据的指针

/*
 * 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;
		}
	}
}

/*
 * 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;
}

/*
 * 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;
}