将点映射到希尔伯特曲线和从希尔伯特曲线映射点

Question

我一直在尝试为 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.

编写一个函数

不幸的是，我玩弄了上面的代码并查看了论文，但我不确定如何让它工作。 似乎有问题的是我的代码向后绘制了 2 位二维希尔伯特曲线的后半部分。如果你在最后一列中绘制二维坐标，你将看到曲线的后半部分（位置 8 及以上）向后。

我不认为我被允许 post 原始的 C 代码，但下面的 C++ 版本只是轻微编辑。我的代码中有几处不同。

C 代码对类型不那么严格，所以我不得不使用 std::bitset
除了前面提到的SEpost中@PaulChernoch提到的bug，接下来的for循环段错误也是如此。
论文奇怪地表示一维坐标。他们称其为数字的“转置”。我编写了一个从常规整数生成“转置”的函数。

关于此算法的另一件事：它不会生成单位间隔和单位超立方体之间的映射。相反，它扩展了问题并在间隔和立方体之间映射，单位为 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)

这是代码（在 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-- ) // 
        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 Whosebug 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";
    }


}

关于另一个我注意到的一些奇怪的事情post:

除了上述SEpost中@PaulChernoch提到的bug外，下一个for循环段错误也是
没有人在谈论论文如何不提供单位区间和单位立方体之间的映射，而是提供从整数到大立方体的映射，并且
我在这里没有看到任何关于论文用于无符号整数“转置”的奇怪表示的提及。

如果您在最后一列中绘制二维坐标，您会看到曲线的后半部分（位置 8 及以上）向后。

Answer 1

“除了@PaulChernoch 在上述 SE post 中提到的错误之外，还有下一个 for 循环段错误。”实际上这是 my 代码中的一个错误——我是 . I started looking at myself after I realized there were other Python packages (e.g. this and this) 使用相同的底层 C 代码。他们都没有抱怨另一个 for 循环。

其次，我上面的函数中也存在一个生成转置的错误。第三，反函数，AxesToTranspose混淆了位数和维数。

All four correct functions（论文中提供的两个完成所有繁重的工作，另外两个用于在整数和“Transpose”之间进行转换）如下：

.

/**
 * @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-- ) // 
        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();

}

我也有一些单元测试here。

将点映射到希尔伯特曲线和从希尔伯特曲线映射点

Mapping points to and from a Hilbert curve

c++

hilbert-curve