网上搜了很多Z3-Ordering实现没搜到,通过 sfcurve-master和geomesa-geomesa-3.2.2 得scala代码改编而来, 环境为C++, vs2015, 理论上windows和Linux都可以用. 不依赖任何库, 这项自身理解和翻译断断续续进行, 最近终于有一点进展, 本次放出Z3, 待全部实现完毕将直接挂出
详细实现代码如下:
#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>
#include <array>
#include <tuple>
using namespace std;
typedef unsigned long long ull;
#pragma region baseinterface
struct IndexRange {
ull lower;
ull upper;
bool contained;
IndexRange()
{
}
IndexRange(ull a, ull b, bool c = true)
: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 false;
}
}CompA;
//也不清楚scala的seq能对应std哪个容器,先用定义后面直接改这个定义即可
using SeqIndexRange = std::list<IndexRange>;
using dblbox = std::array<double, 4>;
using SeqBoxs = std::vector<std::array<double, 4>>;
using SeqTimes = std::vector<std::array<ull, 2>>;
struct XZSFC {
short DefaultPrecision = 12;
double LogPointFive = log(0.5);
};
struct ZRange {
ull min = 0, max = 0, mid = 0, length = 0;
ZRange(ull umin, ull umax)
{
min = umin;
max = umax;
mid = (max + min) >> 1;
length = max - min + 1;
}
ZRange(const ZRange& o)
{
min = o.min;
max = o.max;
mid = o.mid;
length = o.length;
}
bool contains(ull bits)
{
return bits >= min && bits <= max;
}
bool contains(const ZRange& r)
{
return contains(r.min) && contains(r.max);
}
bool overlaps(const ZRange& r) {
return contains(r.min) || contains(r.max);
}
};
using ArrayZRange = std::vector<ZRange>;
class SpaceTimeFillingCurve {
int FullPrecision = 64;
int maxRanges = -1;
public:
virtual ull index(double x, double y, ull t, bool lenient = false) = 0;
virtual void invert(ull z, double& x, double & y, ull&t) = 0;
virtual SeqIndexRange ranges(SeqBoxs seqboxs, SeqTimes seqtimes, int precision, int maxRanges) = 0;
virtual SeqIndexRange ranges(double x1, double x2, double y1, double y2, ull t1, ull t2)
{
SeqBoxs d;
d.resize(1);
d[0][0] = x1; d[0][1] = y1; d[0][2] = x2; d[0][3] = y1;
SeqTimes tseq;
tseq.resize(1);
tseq[0][0] = (t1); tseq[0][1] = (t2);
ranges(d, tseq, FullPrecision, maxRanges);
}
virtual SeqIndexRange ranges(double x1, double x2, double y1, double y2, ull t1, ull t2, int precision)
{
SeqBoxs d;
d.resize(1);
d[0][0] = x1; d[0][1] = y1; d[0][2] = x2; d[0][3] = y1;
SeqTimes tseq;
tseq.resize(1);
tseq[0][0] = (t1); tseq[0][1] = (t2);
ranges(d, tseq, precision, maxRanges);
}
virtual SeqIndexRange ranges(double x1, double x2, double y1, double y2, ull t1, ull t2, int precision, int maxRanges)
{
SeqBoxs d;
d.resize(1);
d[0][0] = x1; d[0][1] = y1; d[0][2] = x2; d[0][3] = y1;
SeqTimes tseq;
tseq.resize(1);
tseq[0][0] = (t1); tseq[0][1] = (t2);
ranges(d, tseq, precision, maxRanges);
}
};
class SpaceFillingCurve {
int FullPrecision = 64;
int maxRanges = -1;
virtual ull index(double x, double y, bool lenient = false) = 0;
virtual void invert(ull z, double& x, double & y) = 0;
virtual SeqIndexRange ranges(SeqBoxs seqboxs, int precision, int maxRanges) = 0;
virtual SeqIndexRange ranges(double x1, double x2, double y1, double y2)
{
SeqBoxs d;
d.resize(1);
d[0][0] = x1; d[0][1] = y1; d[0][2] = x2; d[0][3] = y1;
ranges(d, FullPrecision, maxRanges);
}
virtual SeqIndexRange ranges(double x1, double x2, double y1, double y2, int precision)
{
SeqBoxs d;
d.resize(1);
d[0][0] = x1; d[0][1] = y1; d[0][2] = x2; d[0][3] = y1;
ranges(d, precision, maxRanges);
}
virtual SeqIndexRange ranges(double x1, double x2, double y1, double y2, int precision, int maxRanges)
{
SeqBoxs d;
d.resize(1);
d[0][0] = x1; d[0][1] = y1; d[0][2] = x2; d[0][3] = y1;
ranges(d, precision, maxRanges);
}
};
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;
}
};
struct BitNormalizedDimension {
double min, max;
int maxIndex;
double normalizer;
double denormalizer;
BitNormalizedDimension()
{
}
BitNormalizedDimension(const BitNormalizedDimension & o)
{
min = o.min;
max = o.max;
maxIndex = o.maxIndex;
normalizer = o.denormalizer;
denormalizer = o.denormalizer;
}
BitNormalizedDimension(double dmin, double dmax, int precision)
{
if (!(precision > 0 && precision < 32))
return;
min = dmin;
max = dmax;
ull bins = pow(2, precision);// .toLong1 << precision;
normalizer = bins / (max - min);
denormalizer = (max - min) / bins;
maxIndex = (bins - 1);
}
//case class NormalizedLat(precision : Int) extends BitNormalizedDimension(-90d, 90d, precision)
//case class NormalizedLon(precision : Int) extends BitNormalizedDimension(-180d, 180d, precision)
//case class NormalizedTime(precision : Int, override val max : Double) extends BitNormalizedDimension(0d, max, precision)
int normalize(double x) {
if (x >= max)
return maxIndex;
else
return floor((x - min) * normalizer);
}
double denormalize(int x) {
if (x >= maxIndex)
{
return min + (maxIndex + 0.5) * denormalizer;
}
else
{
return min + (x + 0.5) * denormalizer;
}
}
};
#pragma endregion
#pragma region ZNOrdering
struct ZPrefix
{
ull prefix;
int precision;
ZPrefix()
{
}
ZPrefix(ull llprefix, int iprecision)
{
prefix = llprefix;
precision = iprecision;
}
ZPrefix(ZPrefix&o)
{
prefix = o.prefix;
precision = o.precision;
}
};
/**
* N-dimensional z-curve base class
*/
struct ZN {
// number of bits used to store each dimension
int BitsPerDimension;
// number of dimensions
int Dimensions;
// max value for this z object - can be used to mask another long using &
ull MaxMask;
// total bits used - usually bitsPerDim * dims
int TotalBits;
// number of quadrants in our quad/oct tree - has to be lazy to instantiate correctly
double Quadrants;
int DefaultRecurse = 7;
ZN(int Dimensions)
{
Quadrants = pow(2, Dimensions);
}
/**
* Insert (Dimensions - 1) zeros between each bit to create a zvalue from a single dimension.
* Only the first BitsPerDimension can be considered.
*
* @param value value to split
* @return
*/
virtual ull split(ull value) = 0;
/**
* Combine every (Dimensions - 1) bits to re-create a single dimension. Opposite of split.
*
* @param z value to combine
* @return
*/
virtual int combine(ull z) = 0;
/**
* Is the value contained in the range. Considers user-space.
*
* @param range range
* @param value value to be tested
* @return
*/
virtual bool contains(ZRange &range, ull value) = 0;
/**
* Is the value contained in the range. Considers user-space.
*
* @param range range
* @param value value to be tested
* @return
*/
bool contains(ZRange& range, ZRange value)
{
return contains(range, value.min) && contains(range, value.max);
}
/**
* Does the value overlap with the range. Considers user-space.
*
* @param range range
* @param value value to be tested
* @return
*/
virtual bool overlaps(ZRange& range, ZRange& value) = 0;
/**
* Returns (litmax, bigmin) for the given range and point
*
* @param p point
* @param rmin minimum value
* @param rmax maximum value
* @return (litmax, bigmin)
*/
void zdivide(ull p, ull rmin, ull rmax, ull& minz, ull&maxz)
{
//: (Long, Long) = ZN.zdiv(load, Dimensions)(p, rmin, rmax)
}
// virtual SeqIndexRange ranges(SeqBoxs seqboxs, int precision, int maxRanges) = 0;
virtual SeqIndexRange zranges(ArrayZRange & zbounds)
{
return zranges(zbounds, 64, -1, DefaultRecurse);
}
virtual SeqIndexRange zranges(ArrayZRange & zbounds, int precision)
{
return zranges(zbounds, precision, -1, DefaultRecurse);
}
virtual SeqIndexRange zranges(ArrayZRange & zbounds, int precision, int maxRanges)
{
return zranges(zbounds, precision, maxRanges, DefaultRecurse);
}
/**
* Calculates ranges in index space that match any of the input bounds. Uses breadth-first searching to
* allow a limit on the number of ranges returned.
*
* To improve performance, the following decisions have been made:
* uses loops instead of foreach/maps
* uses java queues instead of scala queues
* allocates initial sequences of decent size
* sorts once at the end before merging
*
* @param zbounds search space
* @param precision precision to consider, in bits (max 64)
* @param maxRanges loose cap on the number of ranges to return. A higher number of ranges will have less
* false positives, but require more processing.
* @param maxRecurse max levels of recursion to apply before stopping
* @return ranges covering the search space
*/
struct ZPrefix
{
ull prefix;
int precision;
};
virtual SeqIndexRange zranges(ArrayZRange & zbounds, int precision = 64, int maxRanges = -1, int maxRecurse = 7)
{
// stores our results - initial size of 100 in general saves us some re-allocation
//val ranges = new java.util.ArrayList[IndexRange](100)
std::list<IndexRange> ranges;
ranges.resize(100);
// values remaining to process - initial size of 100 in general saves us some re-allocation
//val remaining = new java.util.ArrayDeque[(Long, Long)](100)
std::deque<std::pair<ull, ull>> remaining;
remaining.resize(100);
// calculate the common prefix in the z-values - we start processing with the first diff
ZPrefix tmpZPrefix = longestCommonPrefix(zbounds);
ull commonPrefix = tmpZPrefix.prefix;
int commonBits = tmpZPrefix.precision;
int offset = 64 - commonBits;
// checks if a range is contained in the search space
auto isContained = [&zbounds, this](ZRange& range)
{
int i = 0;
while (i < zbounds.size())
{
if (contains(zbounds[i], range))
return true;
i++;
}
return false;
};
auto isOverlapped = [&zbounds, this](ZRange& range)
{
int i = 0;
while (i < zbounds.size())
{
if (overlaps(zbounds[i], range))
return true;
i++;
}
return false;
};
auto checkValue = [&offset, &isContained, &isOverlapped, &precision, &ranges, &remaining](ull prefix, ull quadrant) ->unsigned int {
ull min = prefix | (quadrant << offset);
ull max = min | (1L << offset) - 1;
ZRange quadrantRange = ZRange(min, max);
if (isContained(quadrantRange) || offset < 64 - precision) {
// whole range matches, happy day
ranges.push_back(IndexRange(min, max, true));
}
else if (isOverlapped(quadrantRange)) {
// some portion of this range is excluded
// queue up each sub-range for processing
remaining.push_back(std::make_pair(min, max));
}
};
std::pair<ull, ull> LevelTerminator = std::make_pair(-1, -1);
auto bottomOut = [this, &remaining, &ranges, LevelTerminator]() ->unsigned int {
do {
std::pair<ull, ull> minMax = remaining.front();
auto min2 = std::get<0>(LevelTerminator);
auto max2 = std::get<1>(LevelTerminator);
if (minMax.first == min2 && minMax.second == max2)
{
ranges.emplace_back(minMax.first, minMax.second, false);
}
remaining.pop_front();
} while (!remaining.empty());
};
// initial level - we just check the single quadrant
checkValue(commonPrefix, 0);
remaining.push_back(LevelTerminator);
offset -= Dimensions;
// level of recursion
int level = 0;
int rangeStop = maxRanges == -1 ? INT_MAX : maxRanges;
int recurseStop = (maxRecurse == 7) ? (DefaultRecurse) : maxRecurse;
do {
auto next = remaining.front();
remaining.pop_front();
if (next.first == LevelTerminator.first && next.second == LevelTerminator.second)
{
if (!remaining.empty()) {
level += 1;
offset -= Dimensions;
if (level >= recurseStop || offset < 0) {
bottomOut();
}
else {
remaining.push_back(LevelTerminator);
}
}
}
else {
auto prefix = next.first;
ull quadrant = 0L;
while (quadrant < Quadrants) {
checkValue(prefix, quadrant);
quadrant += 1;
};
// subtract one from remaining.size to account for the LevelTerminator
if (ranges.size() + remaining.size() - 1 >= rangeStop) {
bottomOut();
}
}
} while (!remaining.empty());
// we've got all our ranges - now reduce them down by merging overlapping values
//std::sort(ranges.begin(), ranges.end());
ranges.sort(CompA);
auto current = ranges.front(); // note: should always be at least one range
SeqIndexRange result;
int i = 1;
while (i < ranges.size()) {
auto range = ranges.front();
if (range.lower <= current.upper + 1) {
// merge the two ranges
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;
}
// append the last range - there will always be one left that wasn't added
result.push_back(current);
return result;
}
/**
* Cuts Z-Range in two and trims based on user space, can be used to perform augmented binary search
*
* @param xd: division point
* @param inRange: is xd in query range
*/
static SeqIndexRange g_empty;
SeqIndexRange cut(ZRange &r, ull xd, bool inRange) {
SeqIndexRange ret;
if (r.min == r.max) {
return g_empty;
}
else if (inRange) {
if (xd == r.min) { // degenerate case, two nodes min has already been counted
ret.emplace_back(r.max, r.max);
}
else if (xd == r.max) { // degenerate case, two nodes max has already been counted
ret.emplace_back(r.min, r.min);
}
else {
ret.emplace_back(r.min, xd - 1);
ret.emplace_back(xd + 1, r.max);
}
}
else {
ull litmax0, bigmin = 0;
std::tuple<ull, ull> tup = zdiv(Dimensions, xd, r.min, r.max);
ret.emplace_back(r.min, std::get<0>(tup));
ret.emplace_back(std::get<1>(tup), r.max);
}
}
/**
* Calculates the longest common binary prefix between two z longs
*
* @return (common prefix, number of bits in common)
*/
//这里不是指针是不确定参数大小
ZPrefix longestCommonPrefix(ArrayZRange & values) {
if (values.size() < 1)
return ZPrefix();
std::vector<ull> llu;
for (auto& var : values)
{
llu.push_back(var.min);
llu.push_back(var.max);
}
int bitShift = TotalBits - Dimensions;
ull head = llu[0] >> bitShift;
int tail = llu.size() - 1;
for (int i = llu.size() - 1; i > 0; i--)
{
ull v = llu[i] >> bitShift;
if ((v == head) && bitShift > -1)
{
bitShift -= Dimensions;
head = llu[0] >> bitShift;
}
}
bitShift += Dimensions; // increment back to the last valid value
ZPrefix ret;
ret.prefix = (llu[0] & (((ull)~((ull)0)) << bitShift));
ret.precision = 64 - bitShift;
return ret;
}
// case class ZPrefix(prefix : Long, precision : Int) // precision in bits
/** Loads either 1000... or 0111... into starting at given bit index of a given dimension */
ull load(ull target, ull p, int bits, int dim)
{
auto mask = ~(split(MaxMask >> (BitsPerDimension - bits)) << dim);
auto wiped = target & mask;
return wiped | (split(p) << dim);
}
/**
* Implements the the algorithm defined in: Tropf paper to find:
* LITMAX: maximum z-index in query range smaller than current point, xd
* BIGMIN: minimum z-index in query range greater than current point, xd
*
* @param load: function that knows how to load bits into appropraite dimension of a z-index
* @param xd: z-index that is outside of the query range
* @param rmin: minimum z-index of the query range, inclusive
* @param rmax: maximum z-index of the query range, inclusive
* @return (LITMAX, BIGMIN)
*/
private:
//std::tuple<ull, ull> LevelTerminator = std::make_tuple(-1,-1);
std::tuple<ull, ull> zdiv(int dims, ull xd, ull rmin, ull rmax)
{
//require(rmin < rmax, "min ($rmin) must be less than max $(rmax)")
ull zmin = rmin;
ull zmax = rmax;
ull bigmin = 0L;
ull litmax = 0L;
//int dims = 2;
auto bit = [&](ull x, int idx) { return (int)((x & (1ULL << idx)) >> idx); };
auto over = [&](ull bits) {return 1L << (bits - 1); };
auto under = [&](ull bits) { return (1L << (bits - 1)) - 1; };
auto i = 64;
while (i > 0) {
i -= 1;
auto bits = i / dims + 1;
auto dim = i % dims;
auto a = bit(xd, i);
auto b = bit(zmin, i);
auto c = bit(zmax, i);
if (a == 0 && b == 0 && c == 1)
{
zmax = load(zmax, under(bits), bits, dim);
bigmin = load(zmin, over(bits), bits, dim);
}
else if (a == 0 && b == 1 && c == 1)
{
bigmin = zmin;
return std::make_tuple(litmax, bigmin);
}
else if (a == 1 && b == 0 && c == 0)
{
litmax = zmax;
return std::make_tuple(litmax, bigmin);
}
else if (a == 1 && b == 0 && c == 1)
{
litmax = load(zmax, under(bits), bits, dim);
zmin = load(zmin, over(bits), bits, dim);
}
else
{
continue;
}
};
return std::make_tuple(litmax, bigmin);
}
};
#pragma endregion
#pragma region Z3Ordering =Z3SFC.scala
enum class TimePeriod
{
eDay,
eWeek,
eMonth,
eYear,
};
struct BinnedTime
{
//case TimePeriod.Day = > ChronoUnit.DAYS.getDuration.toMillis
//case TimePeriod.Week = > ChronoUnit.WEEKS.getDuration.toMillis / 1000L
//case TimePeriod.Month = > (ChronoUnit.DAYS.getDuration.toMillis / 1000L) * 31L
// based on 365 days + 1 leap day, with a fudge factor of 10 minutes to account for leap seconds added each year
//case TimePeriod.Year = > (ChronoUnit.DAYS.getDuration.toMinutes * 366L) + 10L
//86400000
//604800
//2678400
//527050
static ull maxOffset(TimePeriod period)
{
switch (period)
{
case TimePeriod::eDay:
return 86400000;
case TimePeriod::eWeek:
return 604800;
case TimePeriod::eMonth:
return 2678400;
case TimePeriod::eYear:
return 527050;
default:
break;
}
}
};
struct Z3 :public ZN {
int Dimensions = 2;
int BitsPerDimension = 31;
int TotalBits = 62;
ull MaxMask = 0x1fffffL;
ull z = 0;
int d0 = 0, d1 = 0, d2 = 0;
int DefaultRecurse = 7;
Z3(int dim = 3) :ZN(dim)
{
Dimensions = dim;
}
Z3(ull z) :ZN(3)
{
d0 = combine(z);
d1 = combine(z >> 1);
d2 = combine(z >> 2);
}
Z3(int x, int y, int t) :ZN(3)
{
ull z = (split(x) | split(y) << 1 | split(t) << 2);
d0 = combine(z);
d1 = combine(z >> 1);
d2 = combine(z >> 2);
}
//def < (other: Z3) = z < other.z
//def >(other: Z3) = z > other.z
//def >= (other: Z3) = z >= other.z
//def <= (other: Z3) = z <= other.z
//def + (offset: Long) = new Z3(z + offset)
//def - (offset: Long) = new Z3(z - offset)
//def == (other: Z3) = other.z == z
//def decode : (Int, Int, Int) = (d0, d1, d2)
int dim(int i)
{
if (i == 0) return d0;
else if (i == 1)return d1;
else if (i == 2)return d2;
else {
return 0;//throw new IllegalArgumentException(s"Invalid dimension $i - valid dimensions are 0,1,2")
}
}
bool inRange(Z3 &rmin, Z3 &rmax)
{
int x = d0; int y = d1; int z = d2;
return x >= rmin.d0 &&
x <= rmax.d0 &&
y >= rmin.d1 &&
y <= rmax.d1 &&
z >= rmin.d2 &&
z <= rmax.d2;
}
Z3 mid(Z3 &p)
{
int x = d0; int y = d1; int z = d2;
int px = p.d0; int py = p.d1; int pz = p.d2;
return Z3((x + px) >> 1, (y + py) >> 1, (z + pz) >> 1);
}
//toString 方法未实现
//def bitsToString = f"(${z.toBinaryString.toLong}%016d)(${d0.toBinaryString.toLong}%08d,${d1.toBinaryString.toLong}%08d,${d2.toBinaryString.toLong}%08d)"
//override def toString = f"$z $decode"
//def unapply(z : Z3) : Option[(Int, Int, Int)] = Some(z.decode)
/** insert 00 between every bit in value. Only first 21 bits can be considered. */
virtual ull split(ull value)
{
auto x = value & MaxMask;
x = (x | x << 32) & 0x1f00000000ffffL;
x = (x | x << 16) & 0x1f0000ff0000ffL;
x = (x | x << 8) & 0x100f00f00f00f00fL;
x = (x | x << 4) & 0x10c30c30c30c30c3L;
return (x | x << 2) & 0x1249249249249249L;
}
/** combine every third bit to form a value. Maximum value is 21 bits. */
virtual int combine(ull value)
{
auto x = z & 0x1249249249249249L;
x = (x ^ (x >> 2)) & 0x10c30c30c30c30c3L;
x = (x ^ (x >> 4)) & 0x100f00f00f00f00fL;
x = (x ^ (x >> 8)) & 0x1f0000ff0000ffL;
x = (x ^ (x >> 16)) & 0x1f00000000ffffL;
x = (x ^ (x >> 32)) & MaxMask;
return (int)x;// x.toInt
}
virtual bool contains(ZRange& range, ull value)
{
Z3 ztmp(value);
int x = ztmp.d0, y = ztmp.d1, z = ztmp.d2;
return x >= Z3(range.min).d0 && x <= Z3(range.max).d0 &&
y >= Z3(range.min).d1 && y <= Z3(range.max).d1 &&
z >= Z3(range.min).d2 && z <= Z3(range.max).d2;
}
virtual bool overlaps(int a1, int a2, int b1, int b2)
{
return max(a1, b1) <= min(a2, b2);
}
virtual bool overlaps(ZRange &range, ZRange& value)
{
return overlaps(Z3(range.min).d0, Z3(range.max).d0, Z3(value.min).d0, Z3(value.max).d0) &&
overlaps(Z3(range.min).d1, Z3(range.max).d1, Z3(value.min).d1, Z3(value.max).d1) &&
overlaps(Z3(range.min).d2, Z3(range.max).d2, Z3(value.min).d2, Z3(value.max).d2);
}
};
/**
* Z3 space filling curve
*
* @param period time period used to bin results
* @param precision bits used per dimension - note all precisions must sum to less than 64
*/
class Z3SpaceTimeFillingCurve :public SpaceTimeFillingCurve {
BitNormalizedDimension lat;// (-90, 90, precision);
BitNormalizedDimension lon;// (-180, 180, precision);
BitNormalizedDimension time;// (0, maxtime, precision);
public:
Z3SpaceTimeFillingCurve(TimePeriod period, int precision = 21)
{
if (precision <= 0 || precision >= 22)
{
printf("error Z3SpaceTimeFillingCurve Precision (bits) per dimension must be in [1,21]");
return;
}
lat = BitNormalizedDimension(-90, 90, precision);
lon = BitNormalizedDimension(-180, 180, precision);
double maxtime = BinnedTime::maxOffset(period);
lon = BitNormalizedDimension(0, maxtime, precision);
SeqTimes wholePeriod;
wholePeriod.resize(1);
wholePeriod[0][0] = time.min;
wholePeriod[0][1] = time.max;
}
virtual ull index(double x, double y, ull t, bool lenient = false)
{
if (x >= lon.min && x <= lon.max && y >= lat.min && y <= lat.max && t >= time.min && t <= time.max)
{
return Z3(lon.normalize(x), lat.normalize(y), time.normalize(t)).z;
}
else
{
return 0;
}
//catch {
//case _: IllegalArgumentException if lenient = > lenientIndex(x, y, t)
}
ull lenientIndex(double x, double y, ull t)
{
double bx = 0; if (x < lon.min) { bx = lon.min; }
else if (x > lon.max) { bx = lon.max; }
else { bx = x; }
double by = 0; if (y < lat.min) { by = lat.min; }
else if (y > lat.max) { by = lat.max; }
else { by = y; }
double bt = 0; if (t < time.min) { bt = time.min; }
else if (t > time.max) { bt = time.max; }
else { bt = t; }
return Z3(lon.normalize(bx), lat.normalize(by), time.normalize(bt)).z;
}
virtual void invert(ull z, double& x, double & y, ull&t)
{
x = lon.denormalize(Z3(z).d0);
y = lat.denormalize(Z3(z).d1);
t = time.denormalize(Z3(z).d2);
}
virtual SeqIndexRange ranges(SeqBoxs seqboxs, SeqTimes seqtimes, int precision, int maxRanges)
{
ArrayZRange az;
for (int i = 0; i < seqboxs.size(); i++)
{
ull min = index(seqboxs[i][0], seqboxs[i][1], seqtimes[i][0]);
ull max = index(seqboxs[i][2], seqboxs[i][3], seqtimes[i][1]);
az.emplace_back(min, max);
}
Z3 zn(3);
return zn.zranges(az, precision, maxRanges, 7);
}
};
#pragma endregion
