XZordering 的C++实现

网上搜了很多XZ-Ordering实现没搜到,通过 sfcurve-master和geomesa-geomesa-3.2.2 得scala代码改编而来, 环境为C++, vs2015, 理论

#include "stdafx.h"
#include <vector>
#include <list>
#include <queue>
#include <map>
#include <algorithm>
#include <iostream>
#include <deque>
#include <math.h>
#include <time.h>
#include <functional>  
#include <memory>
#include <set>
#include <map>
#include <queue>
#include <string>
#include <memory>
#include <functional>
#include <math.h>
using namespace std;
typedef unsigned long long ull;

struct IndexRange {
	ull lower;
	ull upper;
	bool contained;
	IndexRange()
	{

	}
	IndexRange(ull a, ull b, bool c)
		:lower(a), upper(b), contained(c)
	{

	}
	IndexRange(const IndexRange& o)
	{
		lower = o.lower;
		upper = o.upper;
		contained = o.contained;
	}

	int operator ()(const IndexRange& o)const
	{
		if (lower != o.lower) return lower > o.lower ? 1 : -1;
		if (upper != o.upper) return upper > o.upper ? 1 : -1;

		return 0;
	}
};
/// 用于 std::sort 对结构体按指定成员比较
struct {
	bool operator()(const IndexRange &a, const IndexRange &b) const
	{
		if (a.lower != b.lower) return a.lower > b.lower ? 1 : -1;
		if (a.upper != b.upper) return a.upper > b.upper ? 1 : -1;

		return 0;
	}
}CompA;

#pragma region XZOrdering
struct XZOrdering
{
	double xLo = 0;
	double xHi = 0;
	double yLo = 0;
	double yHi = 0;

	double xSize = 0;
	double  ySize = 0;
	short DefaultPrecision = 12;
	double	 LogPointFive = log(0.5);
	short g = 0;
	XZOrdering(short sg, double xmin,  double ymin, double xmax, double ymax)
	{
		g = sg;
		xLo = xmin;
		xHi = xmax;
		yLo = ymin;
		yHi = ymax;
		xSize = xHi - xLo;
		ySize = yHi - yLo;
	}
	bool predicate(double min, double max, double w2)
	{
		return max <= (floor(min / w2) * w2) + (2 * w2);
	};
	ull index(double xmin, double ymin, double xmax, double ymax, bool lenient = false)
	{
		// normalize inputs to [0,1]
		double nxmin, nymin, nxmax, nymax;
		normalize(xmin, ymin, xmax, ymax, lenient, nxmin, nymin, nxmax, nymax);

		// calculate the length of the sequence code (section 4.1 of XZ-Ordering paper)

		double maxDim = max(nxmax - nxmin, nymax - nymin);

		// l1 (el-one) is a bit confusing to read, but corresponds with the paper's definitions
		int l1 = floor(log(maxDim) / LogPointFive);

		// the length will either be (l1) or (l1 + 1)
		ull length = 0;
		if (l1 >= g)
		{
			length = g;
		}
		else
		{
			double w2 = pow(0.5, l1 + 1);
			if (predicate(nxmin, nxmax, w2) && predicate(nymin, nymax, w2))
				length = l1 + 1;
			else
				length = l1;
		}

		return sequenceCode(nxmin, nymin, length);
	}

	ull  sequenceCode(double x, double y, int length) {
		double xmin = 0.0;
		double ymin = 0.0;
		double xmax = 1.0;
		double ymax = 1.0;

		ull cs = 0;

		int i = 0;
		while (i < length) {
			double xCenter = (xmin + xmax) / 2.0;
			double yCenter = (ymin + ymax) / 2.0;
			bool t1 = x < xCenter;
			bool t2 = y < yCenter;

			if (t1 == true && t2 == true)
			{
				cs += 1L; xmax = xCenter; ymax = yCenter;
			}
			else if (t1 == false && t2 == true)
			{
				cs += 1L + 1L * ((ull)pow(4, g - i) - 1L) / 3L; xmin = xCenter; ymax = yCenter;
			}
			else if (t1 == true && t2 == false) {
				cs += 1L + 2L * ((ull)pow(4, g - i) - 1L) / 3L; xmax = xCenter; ymin = yCenter;
			}
			else  if (t1 == false && t2 == false)
			{
				cs += 1L + 3L * ((ull)pow(4, g - i) - 1L) / 3L; xmin = xCenter; ymin = yCenter;
			}
			i += 1;
		}

		return cs;
	}

	void sequenceInterval(double x, double y, short length, bool partial, ull& omin, ull&omax)
	{
		ull min = sequenceCode(x, y, length);
		ull max = min;
		if (partial)
		{
			max = min;
		}
		else
		{
			max = min + ((ull)pow(4, g - length + 1) - 1L) / 3L;
		}
		omin = min;
		omax = max;
	}
	void  normalize(double xmin, double ymin, double xmax, double ymax, bool lenient,
		double &dxmin, double&dymin, double &dxmax, double& dymax)
	{
		dxmin = (xmin - xLo) / xSize;
		dymin = (ymin - yLo) / ySize;
		dxmax = (xmax - xLo) / xSize;
		dymax = (ymax - yLo) / ySize;

	}

	struct QueryWindow
	{
		double xmin; double xmax; double ymin; double ymax;
		QueryWindow()
		{

		}
		QueryWindow(double dxmin, double dymin , double dxmax, double dymax) :
			xmin(dxmin), xmax(dxmax), ymin(dymin), ymax(dymax)
		{

		}

		QueryWindow(const QueryWindow& o)
		{
			xmin = o.xmin;
			ymin = o.ymin;
			xmax = o.xmax;
			ymax = o.ymax;
		}
	};
	struct XElement
	{
		double xmin; double xmax; double ymin; double ymax; double length;

		double xext = 0;
		double yext = 0;
		XElement() {}
		XElement(double dxmin, double dymin, double dxmax, double dymax, double dlength)
			:xmin(dxmin), xmax(dxmax), ymin(dymin), ymax(dymax), length(dlength)
		{
			xext = xmax + length;
			yext = ymax + length;
		}
		XElement(const XElement& o)
		{
			xmin = o.xmin;
			xmax = o.xmax;
			ymin = o.ymin;
			ymax = o.ymax;
			length = o.length;
			xext = o.xext;
			yext = o.yext;
		}
		XElement clone()
		{
			return XElement(xmin, ymin, xmax, ymax, length);
		}
		bool isContained(QueryWindow window)
		{
			return window.xmin <= xmin && window.ymin <= ymin && window.xmax >= xext && window.ymax >= yext;
		}
		bool overlaps(QueryWindow window)
		{
			return 	window.xmax >= xmin && window.ymax >= ymin && window.xmin <= xext && window.ymin <= yext;
		}
		std::vector<XElement*> children() {
			double xCenter = (xmin + xmax) / 2.0;
			double yCenter = (ymin + ymax) / 2.0;
			double len = length / 2.0;
			std::vector<XElement*> ret;
			XElement* c0 = new XElement(xmin,  ymin, xmax, ymax, length);
			c0->xmax = xCenter;
			c0->ymax = yCenter;
			c0->length = len;
			XElement* c1 = new XElement(xmin,  ymin, xmax, ymax, length);
			c1->xmin = xCenter; c1->ymax = yCenter; c1->length = len;
			XElement* c2 = new XElement(xmin, ymin, xmax, ymax, length);
			c2->xmax = xCenter; c2->ymin = yCenter; c2->length = len;
			XElement* c3 = new XElement(xmin,  ymin, xmax, ymax, length);
			c3->xmin = xCenter; c3->ymin = yCenter; c3->length = len;
			ret.push_back(c0);
			ret.push_back(c1);
			ret.push_back(c2);
			ret.push_back(c3);
			return ret;
		}
	};
	// the initial level of quads
	std::vector<XElement*> LevelOneElements = XElement(0.0, 0.0, 1.0, 1.0, 1.0).children();
	XElement* LevelTerminator = new XElement(-1.0, -1.0, -1.0, -1.0, 0);
	std::map<short, XZOrdering*> cache;
	XZOrdering(short g)
	{
		auto sfc = cache.find(g);
		if (sfc == cache.end())
		{
			XZOrdering *sfct = new XZOrdering(g, -180.0, -90, 180.0, 90);
			cache[g] = sfct;
		}
	}

	list<IndexRange> ranges(vector<QueryWindow>& query, int rangeStop = INT_MAX)
	{
		// stores our results - initial size of 100 in general saves us some re-allocation
		vector<IndexRange> ranges;
		//ranges;
		// values remaining to process - initial size of 100 in general saves us some re-allocation
		deque<XElement*> remaining;
		//remaining.resize(100);

		auto isContained = [&query](XElement& quad)
		{
			int i = 0;
			while (i < query.size())
			{
				if (quad.isContained(query[i]))
				{
					return true;
				}
				i += 1;
			};
			return false;
		};


		// checks if a quad overlaps the search space
		auto isOverlapped = [&query](XElement & quad)
		{
			int i = 0;
			while (i < query.size())
			{
				if (quad.overlaps(query[i]))
				{
					return true;
				}
				i += 1;
			}
			return false;
		};

		XZOrdering* ths = this;
		auto checkValue = [&ranges, &remaining, isContained, isOverlapped, ths](XElement * quad, short level)
		{
			ull  min, max;
			if (isContained(*quad))
			{
				ths->sequenceInterval(quad->xmin, quad->ymin, level, false, min, max);
				ranges.emplace_back(min, max, true);
			}
			else if (isOverlapped(*quad))
			{
				ths->sequenceInterval(quad->xmin, quad->ymin, level, true, min, max);
				ranges.emplace_back(min, max, false);
				auto cds = quad->children();
				for (int i = 0; i < cds.size(); i++)
					remaining.push_back(cds[i]);

			}
		};


		for (int i = 0; i < LevelOneElements.size(); i++)
			remaining.push_back(LevelOneElements[i]);
		remaining.push_back(LevelTerminator);
		// level of recursion
		short level = 1;

		while (level < g && !remaining.empty() && ranges.size() < rangeStop)
		{
			auto next = remaining.front();
			remaining.pop_front();
			if (next == (LevelTerminator))
			{
				if (!remaining.empty())
				{
					level = (level + 1);
					remaining.push_back(LevelTerminator);
				}
			}
			else {
				checkValue(next, level);
			}
		};

		while (!remaining.empty()) {
			auto quad = remaining.front();
			remaining.pop_front();
			if (quad == LevelTerminator)
			{
				level = (level + 1);
			}
			else {
				ull min, max;
				sequenceInterval(quad->xmin, quad->ymin, level, false, min, max);
				ranges.emplace_back(min, max, false);
			}
		}

		std::sort(ranges.begin(), ranges.end(), CompA);


		auto current = ranges[0];
		list<IndexRange> result;
		int i = 1;
		while (i < ranges.size()) {
			auto range = ranges[i];
			if (range.lower <= current.upper + 1)
			{
				current = IndexRange(current.lower, max(current.upper, range.upper), current.contained && range.contained);
			}
			else {
				// append the last range and set the current range for future merging
				result.push_back(current);
				current = range;
			}
			i += 1;
		}
		result.push_back(current);
		return result;
	}
};



int testXZ()
{
	XZOrdering f(31, -180, -90, 180, 90);
	ull g = f.index(78, 15, 112, 45);

	XZOrdering sfc(12, -180,  -90, 180, 90);
	ull poly = sfc.index(10, 10, 12, 12);
	ull polyright = sfc.index(12.1, 10, 14, 12);
	ull polyleft = sfc.index(8, 10, 10, 12);
	ull polytop = sfc.index(10, 13, 12, 15);
	ull polybottom = sfc.index(10, 8, 12, 10);
	ull polybottomright = sfc.index(12, 8, 14, 10);

	std::vector<XZOrdering::QueryWindow> containing_overlapping = {
		XZOrdering::QueryWindow(9.0, 9.0, 13.0, 13.0),
		XZOrdering::QueryWindow(-180.0, -90.0, 180.0, 90.0),
		XZOrdering::QueryWindow(0.0, 0.0, 180.0, 90.0),
		XZOrdering::QueryWindow(0.0, 0.0, 20.0, 20.0) ,
		XZOrdering::QueryWindow(11.0, 11.0, 13.0, 13.0),
		XZOrdering::QueryWindow(9.0, 9.0, 11.0, 11.0),
		XZOrdering::QueryWindow(10.5, 10.5, 11.5, 11.5),
		XZOrdering::QueryWindow(11.0, 11.0, 11.0, 11.0)
	};
	// note: in general, some disjoint ranges will match due to false positives
	std::vector<XZOrdering::QueryWindow>  disjoint = {
		XZOrdering::QueryWindow(-180.0, -90.0, 8.0, 8.0),
		XZOrdering::QueryWindow(0.0, 0.0, 8.0, 8.0),
		XZOrdering::QueryWindow(9.0, 9.0, 9.5, 9.5),
		XZOrdering::QueryWindow(20.0, 20.0, 180.0, 90.0) };
	for each(auto item  in containing_overlapping)
	{
		std::vector<XZOrdering::QueryWindow> q;
		q.push_back(item);
		auto ranges = sfc.ranges(q);
		for each(auto itemr  in ranges)
		{
			if (itemr.lower <= poly && itemr.upper <= poly)
			{
				printf("$bbox -  match");
			}
			else
			{
				printf("$bbox - no match");
			}
		}

	}
	//forall(disjoint) {
	//	bbox = >
	//		val ranges = sfc.ranges(Seq(bbox)).map(r = > (r.lower, r.upper))
	//		val matches = ranges.exists(r = > r._1 <= poly && r._2 >= poly)
	//		if (matches) {
	//			logger.warn(s"$bbox - invalid match")
	//		}
	//	matches must beFalse

	return 0;
}

int main()
{
	testXZ();
	return 0;
}