Mapping points to and from a Hilbert curve

Question

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.

1. C code is less strict on types, so I had to use std::bitset

2. In addition to the bug mentioned by @PaulChernoch in the aforementioned SE post, the next for loop segfaults, too.

3. 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:

1. In addition to the bug mentioned by @PaulChernoch in the above mentioned SE post, the next for loop segfaults, too.

2. 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

3. 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.

Answer 1

"In addition to the bug mentioned by @PaulChernoch in the above mentioned SE post, the next for loop segfaults, too." Actually it was a bug in my code--I was having a hard time iterating over a container backwards. I started looking at myself after I realized there were other Python packages (e.g. this and this) that use the same underlying C code. Neither of them were complaining about the other for loop.

Second, there was a bug in my function above that generates the transpose, too. And third, the inverse function, AxesToTranspose confused number of bits with number of dimensions.

All four correct functions (the two provided in the paper that do all the heavy lifting, and two more for converting between integers and "Transpose"s) are as follows:

.

/**
 * @brief converts an integer in a transpose form to a position on the Hilbert Curve.
 * Code is based off of John Skilling , "Programming the Hilbert curve", 
 * AIP Conference Proceedings 707, 381-387 (2004) https://doi.org/10.1063/1.1751381
 * @file resamplers.h
 * @tparam num_bits how "accurate/fine/squiggly" you want the Hilbert curve
 * @tparam num_dims the number of dimensions the curve is in
 * @param X an unsigned integer in a "Transpose" form.
 * @return a position on the hilbert curve 
 */
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>;


    // Gray decode by H ^ (H/2)
    coord_t t = X[num_dims-1] >> 1;
    for(int i = num_dims-1; i > 0; i-- ) // https://stackoverflow.com/a/10384110
        X[i] ^= X[i-1];
    X[0] ^= t;

    // Undo excess work
    coord_t N = 2 << (num_bits-1);
    for( coord_t Q = 2; Q != N; Q <<= 1 ) {
        coord_t P = Q.to_ulong() - 1;
        for( int i = num_dims - 1; i >= 0; i--){
            if( (X[i] & Q).any() ){ // invert low bits of X[0]
                X[0] ^= P;
            } else{ // exchange low bits of X[i] and X[0]
                t = (X[0]^X[i]) & P;
                X[0] ^= t;
                X[i] ^= t;
            }
        }
    }

    return X;
}


/**
 * @brief converts a position on the Hilbert curve into an integer in a "transpose" form.
 * Code is based off of John Skilling , "Programming the Hilbert curve", 
 * AIP Conference Proceedings 707, 381-387 (2004) https://doi.org/10.1063/1.1751381
 * @file resamplers.h
 * @tparam num_bits how "accurate/fine/squiggly" you want the Hilbert curve
 * @tparam num_dims the number of dimensions the curve is in
 * @param X a position on the hilbert curve (each dimension coordinate is in base 2)
 * @return a position on the real line (in a "Transpose" form)
 */
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>;

    // Inverse undo
    coord_t M = 1 << (num_bits-1);
    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_dims; 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_dims; 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_dims; i++ )
        X[i] ^= t;

    return X;
}


/**
 * @brief converts an integer on the positive integers into its "Transpose" representation..
 * This code supplements the above two functions that are 
 *  based off of John Skilling , "Programming the Hilbert curve", 
 * AIP Conference Proceedings 707, 381-387 (2004) https://doi.org/10.1063/1.1751381
 * @file resamplers.h
 * @tparam num_bits how "accurate/fine/squiggly" you want the Hilbert curve
 * @tparam num_dims the number of dimensions the curve is in
 * @param H a position on the hilbert curve (0,1,..,2^(num_dims * num_bits) )
 * @return a position on the real line (in a "Transpose" form)
 */
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 dim_coord_tmp;
        unsigned start_bit = num_bits*num_dims-1-dim;
        unsigned int c = num_bits - 1;
        for(int bit = start_bit; bit >= 0; bit -= num_dims){
            dim_coord_tmp[c] = Hb[bit];
            c--;
        }
        X[dim] = dim_coord_tmp;
    }
    return X;
}


/**
 * @brief converts an integer in its "Transpose" representation into a positive integer..
 * This code supplements two functions above that are 
 *  based off of John Skilling , "Programming the Hilbert curve", 
 * AIP Conference Proceedings 707, 381-387 (2004) https://doi.org/10.1063/1.1751381
 * @file resamplers.h
 * @tparam num_bits how "accurate/fine/squiggly" you want the Hilbert curve
 * @tparam num_dims the number of dimensions the curve is in
 * @param Htrans a position on the real line (in a "Transpose" form)  
 * @return a position on the hilbert curve (0,1,..,2^(num_dims * num_bits) )
 */
template<size_t num_bits, size_t num_dims>
unsigned int makeH(std::array<std::bitset<num_bits>,num_dims> Htrans)
{
    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 H;
    unsigned int which_dim = 0;
    unsigned which_bit;
    for(int i = num_bits*num_dims - 1; i >= 0; i--){
        which_bit = i / num_dims;
        H[i] = Htrans[which_dim][which_bit];
        which_dim = (which_dim + 1) % num_dims;
    }
    return H.to_ulong();

}

I have some unit tests here as well.