I have been trying to write a function for the Hilbert curve map and inverse map. Fortunately there was another SE post on it, and the accepted answer was highly upvoted, and featured code based on a paper in a peer-reviewed academic journal.
Unfortunately, I played around with the code above and looked at the paper, and I'm not sure how to get this to work. What appears to be broken is that my code drawing the second half of a 2-bit 2-dimensional Hilbert Curve backwards. If you draw out the 2-d coordinates in the last column, you'll see the second half of the curve (position 8 and on) backwards.
I don't think I'm allowed to post the original C code, but the C++ version below is only lightly edited. A few things that are different in my code.
C code is less strict on types, so I had to use
std::bitsetIn addition to the bug mentioned by @PaulChernoch in the aforementioned SE post, the next
forloop segfaults, too.The paper represents one-dimensional coordinates weirdly. They call it the number's "Transpose." I wrote a function that produces a "Transpose" from a regular integer.
Another thing about this algorithm: it doesn't produce a map between unit intervals and unit hypercubes. Rather, it stretches the problem out and maps between intervals and cubes with unit spacing.
NB: HTranspose is the representation of H they use in the paper
H, HTranspose, mappedCoordinates
------------------------------------
0: (00, 00), (0, 0)
1: (00, 01), (1, 0)
2: (01, 00), (1, 1)
3: (01, 01), (0, 1)
4: (00, 10), (0, 2)
5: (00, 11), (0, 3)
6: (01, 10), (1, 3)
7: (01, 11), (1, 2)
8: (10, 00), (3, 2)
9: (10, 01), (3, 3)
10: (11, 00), (2, 3)
11: (11, 01), (2, 2)
12: (10, 10), (2, 1)
13: (10, 11), (3, 1)
14: (11, 10), (3, 0)
15: (11, 11), (2, 0)
Here's the code (in c++).
#include <array>
#include <bitset>
#include <iostream>
#include <cmath>
namespace hilbert {
/// The Hilbert index is expressed as an array of transposed bits.
///
/// Example: 5 bits for each of n=3 coordinates.
/// 15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored
/// as its Transpose ^
/// X[0] = A D G J M X[2]| 7
/// X[1] = B E H K N <-------> | /X[1]
/// X[2] = C F I L O axes |/
/// high low 0------> X[0]
template<size_t num_bits, size_t num_dims>
std::array<std::bitset<num_bits>,num_dims> TransposeToAxes(std::array<std::bitset<num_bits>,num_dims> X)
{
using coord_t = std::bitset<num_bits>;
using coords_t = std::array<coord_t, num_dims>;
coord_t N = 2 << (num_bits-1);
// Gray decode by H ^ (H/2)
coord_t t = X[num_dims-1] >> 1;
for(size_t i = num_dims-1; i > 0; i-- ) // https://stackoverflow.com/a/10384110
X[i] ^= X[i-1];
X[0] ^= t;
// Undo excess work
for( coord_t Q = 2; Q != N; Q <<= 1 ) {
coord_t P = Q.to_ulong() - 1;
for( size_t i = num_dims-1; i > 0 ; i-- ){ // did the same stackoverflow thing
if( (X[i] & Q).any() )
X[0] ^= P;
else{
t = (X[0]^X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
}
}
return X;
}
template<size_t num_bits, size_t num_dims>
std::array<std::bitset<num_bits>,num_dims> AxesToTranspose(std::array<std::bitset<num_bits>, num_dims> X)
{
using coord_t = std::bitset<num_bits>;
using coords_t = std::array<coord_t, num_dims>;
coord_t M = 1 << (num_bits-1);
// Inverse undo
for( coord_t Q = M; Q.to_ulong() > 1; Q >>= 1 ) {
coord_t P = Q.to_ulong() - 1;
for(size_t i = 0; i < num_bits; i++ ){
if( (X[i] & Q).any() )
X[0] ^= P;
else{
coord_t t = (X[0]^X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
}
} // exchange
// Gray encode
for( size_t i = 1; i < num_bits; i++ )
X[i] ^= X[i-1];
coord_t t = 0;
for( coord_t Q = M; Q.to_ulong() > 1; Q >>= 1 ){
if( (X[num_dims-1] & Q).any() )
t ^= Q.to_ulong()-1;
}
for( size_t i = 0; i < num_bits; i++ )
X[i] ^= t;
return X;
}
template<size_t num_bits, size_t num_dims>
std::array<std::bitset<num_bits>,num_dims> makeHTranspose(unsigned int H)
{
using coord_t = std::bitset<num_bits>;
using coords_t = std::array<coord_t, num_dims>;
using big_coord_t = std::bitset<num_bits*num_dims>;
big_coord_t Hb = H;
coords_t X;
for(size_t dim = 0; dim < num_dims; ++dim){
coord_t tmp;
unsigned c = num_dims - 1;
for(size_t rbit = dim; rbit < num_bits*num_dims; rbit += num_dims){
tmp[c] =Hb[num_bits*num_dims - 1 - rbit];
c--;
}
X[dim] = tmp;
}
return X;
}
} //namespace hilbert
int main()
{
constexpr unsigned nb = 2;
constexpr unsigned nd = 2;
using coord_t = std::bitset<nb>;
using coords_t = std::array<coord_t,nd>;
std::cout << "NB: HTranspose is the representation of H they use in the paper\n";
std::cout << "H, HTranspose, mappedCoordinates \n";
std::cout << "------------------------------------\n";
for(unsigned H = 0; H < pow(2,nb*nd); ++H){
// H with the representation they use in the paper
coords_t weirdH = hilbert::makeHTranspose<nb,nd>(H);
std::cout << H << ": ("
<< weirdH[0] << ", "
<< weirdH[1] << "), ("
<< hilbert::TransposeToAxes<nb,nd>(weirdH)[0].to_ulong() << ", "
<< hilbert::TransposeToAxes<nb,nd>(weirdH)[1].to_ulong() << ")\n";
}
}
Some strange things I noticed about the other post:
In addition to the bug mentioned by @PaulChernoch in the above mentioned SE post, the next
forloop segfaults, too.Nobody is talking about how the paper doesn't provide a mapping between the unit interval and the unit cubes, but rather a mapping from integers to big cubes, and
I don't see any mention here about the weird representation the paper uses for unsigned integers "Transpose".
If you draw out the 2-d coordinates in the last column, you'll see the second half of the curve (position 8 and on) backwards.
