映射点与希尔伯特曲线
Posted
技术标签:
【中文标题】映射点与希尔伯特曲线【英文标题】：Mapping points to and from a Hilbert curve 【发布时间】：2021-09-22 15:57:19 【问题描述】：我一直在尝试为Hilbert curve 映射和逆映射编写函数。幸运的是上面有another SE post，并且接受的答案得到了高度评价，并且特色代码基于a paper in a peer-reviewed academic journal。
不幸的是，我玩弄了上面的代码并查看了论文，但我不确定如何让它工作。 似乎有问题的是我的代码向后绘制了 2 位二维希尔伯特曲线的后半部分。如果您在最后一列中绘制二维坐标，您会看到曲线的后半部分（位置 8 及以上）向后。
我不认为我被允许发布原始 C 代码，但下面的 C++ 版本只是经过轻微编辑。我的代码中有一些不同的地方。
C 代码对类型没有那么严格，所以我不得不使用std::bitset
除了前面提到的 SE 帖子中@PaulChernoch 提到的错误之外，下一个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-- ) // https://***.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 *** 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";
我注意到其他帖子的一些奇怪的事情：
除了上面提到的 SE 帖子中@PaulChernoch 提到的错误之外，下一个for 循环段错误也是如此。
没有人在谈论论文如何不提供单位间隔和单位立方体之间的映射，而是提供从整数到大立方体的映射，以及
我在这里没有看到任何关于论文用于无符号整数“转置”的奇怪表示。
如果你在最后一列画出二维坐标，你会看到曲线的后半部分（位置 8 及以上）向后。
【问题讨论】：
【参考方案1】：“除了@PaulChernoch 在上面提到的 SE 帖子中提到的错误之外，下一个 for 循环段错误也是如此。”实际上这是 my 代码中的一个错误——我是 having a hard time iterating over a container backwards。在我意识到还有其他 Python 包（例如 this 和 this）使用相同的底层 C 代码后，我开始审视自己。他们都没有抱怨另一个 for 循环。
其次，我上面的函数中也有一个错误会生成转置。第三，反函数，AxesToTranspose混淆了位数和维数。
All four correct functions（论文中提供的两个用于完成所有繁重的工作，另外两个用于在整数和“转置”之间进行转换）如下：
.
/**
* @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://***.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();
我也有一些单元测试here。
【讨论】：
以上是关于映射点与希尔伯特曲线的主要内容，如果未能解决你的问题，请参考以下文章