compress/wavelet/remap.hpp

#pragma once

#include <climits>
#include <cstdint>
#include <utility>

namespace remap {
// there are some more efficient ways in
// http://graphics.stanford.edu/~seander/bithacks.html#InterleaveTableObvious
// but I was too lazy to understand them

// some templates dealing with sizes
// clang-format off
template <typename T> struct unsigned_size {};
template <> struct unsigned_size<uint8_t>  { using twice = uint16_t;      /* 4 bits */      };
template <> struct unsigned_size<uint16_t> { using twice = uint32_t; using half = uint8_t;  };
template <> struct unsigned_size<uint32_t> { using twice = uint64_t; using half = uint16_t; };
template <> struct unsigned_size<uint64_t> {      /* 128 bits */     using half = uint32_t; };
// clang-format on

// shortcuts
template <typename T>
using twice = typename unsigned_size<T>::twice;
template <typename T>
using half = typename unsigned_size<T>::half;

// converts x,y coordinates to z order
template <typename T>
twice<T> to_z_order(const T xin, const T yin) {
	const twice<T> x = xin;
	const twice<T> y = yin;
	twice<T> z = 0;
	for (int i = 0; i < sizeof(T) * CHAR_BIT; i++) {
		z |= (x & (1U << i)) << i | (y & (1U << i)) << (i + 1);
	}
	return z;
}

// converts z to x,y
template <typename T>
std::pair<half<T>, half<T>> from_z_order(const T z) {
	half<T> x = 0;
	half<T> y = 0;
	for (int i = 0; i < sizeof(T) * CHAR_BIT; i += 2) {
		x |= (z & (1U << i)) >> (i / 2);
		y |= (z & (1U << (i + 1))) >> (i / 2 + 1);
	}
	return {x, y};
}


// from wikipedia
// rotate/flip a quadrant appropriately
template <typename T>
void hilbert_rot(T n, T * x, T * y, T rx, T ry) {
	if (ry == 0) {
		if (rx == 1) {
			*x = n - 1 - *x;
			*y = n - 1 - *y;
		}

		// Swap x and y
		int t = *x;
		*x = *y;
		*y = t;
	}
}

template <typename T>
twice<T> to_hilbert(twice<T> n, T xin, T yin) {
	twice<T> x = xin;
	twice<T> y = yin;
	twice<T> d = 0;
	for (twice<T> s = n / 2; s > 0; s /= 2) {
		twice<T> rx = (x & s) > 0;
		twice<T> ry = (y & s) > 0;
		d += s * s * ((3 * rx) ^ ry);
		hilbert_rot(s, &x, &y, rx, ry);
	}
	return d;
}

// convert d to (x,y)
template <typename T>
std::pair<half<T>, half<T>> from_hilbert(T n, T d) {
	T t = d;
	T x = 0;
	T y = 0;
	for (T s = 1; s < n; s *= 2) {
		T rx = 1 & (t / 2);
		T ry = 1 & (t ^ rx);
		hilbert_rot(s, &x, &y, rx, ry);
		x += s * rx;
		y += s * ry;
		t /= 4;
	}
	return {half<T>(x), half<T>(y)};
}

}