将 N 维值映射到希尔伯特曲线上的一个点 [英] Mapping N-dimensional value to a point on Hilbert curve
问题描述
我有大量的 N 维点(数千万；N 接近 100).
I have a huge set of N-dimensional points (tens of millions; N is close to 100).
我需要将这些点映射到一个维度，同时保留空间局部性.我想使用 Hilbert 空间填充曲线来做到这一点.
I need to map these points to a single dimension while preserving spatial locality. I want to use Hilbert space-filling curve to do it.
对于每个点，我想选择曲线上最近的点.该点的希尔伯特值(从曲线起点到选取点的曲线长度)是我寻求的单维值.
For each point I want to pick the closest point on the curve. The Hilbert value of the point (curve length from the start of curve to the picked point) is the single dimension value I seek.
计算不一定是即时的，但我预计不会超过几个小时在体面的现代家用 PC 硬件上.
Computation does not have to be instant, but I expect it to be no more than several hours on decent modern home PC hardware.
对实施有什么建议吗?有没有可以帮助我的图书馆?(语言并不重要.)
Any suggestions on implementation? Are there any libraries that would help me? (Language does not matter much.)
推荐答案
我终于崩溃了，掏出了一些钱.AIP(美国物理学会)有一篇不错的短文，源代码为 C.John Skilling 的Programming the Hilbert curve"(来自 AIP Conf. Proc. 707, 381 (2004))有一个附录，其中包含以下代码:两个方向的映射.它适用于任意数量的大于 1 的维度，不是递归的，不使用吞噬大量内存的状态转换查找表，并且主要使用位操作.因此它相当快并且具有良好的内存占用.
I finally broke down and shelled out some money. The AIP (American Institute of Physics) has a nice, short article with source code in C. "Programming the Hilbert curve" by John Skilling (from the AIP Conf. Proc. 707, 381 (2004)) has an appendix with code for mappings in both directions. It works for any number of dimensions > 1, is not recursive, does not use state-transition lookup tables that gobble up huge amounts of memory, and mostly uses bit operations. Thus it is reasonably fast and has a good memory footprint.
如果您选择购买文章，我发现源代码中有错误.
If you choose to purchase the article, I discovered an error in the source code.
以下代码行(在函数 TransposetoAxes 中找到)有错误:
The following line of code (found in the function TransposetoAxes) has the error:
for( i = n-1; i >= 0; i-- ) X[i] ^= X[i-1];
for( i = n-1; i >= 0; i-- ) X[i] ^= X[i-1];
修正是将大于等于(>=)改为大于(>).如果没有这种更正，当变量i"变为零时，X 数组将使用负索引访问，从而导致程序失败.
The correction is to change the greater than or equal (>=) to a greater than (>). Without this correction, the X array is accessed using a negative index when the variable "i" becomes zero, causing the program to fail.
我推荐阅读这篇文章(它有七页长，包括代码)，因为它解释了算法的工作原理，这远非显而易见.
I recommend reading the article (which is seven pages long, including the code), as it explains how the algorithm works, which is far from obvious.
我将他的代码翻译成 C# 供我自己使用.代码如下.Skilling 执行原地转换，覆盖您传入的向量.我选择克隆输入向量并返回一个新副本.此外，我将这些方法实现为扩展方法.
I translated his code into C# for my own use. The code follows. Skilling performs the transformation in place, overwriting the vector that you pass in. I chose to make a clone of the input vector and return a new copy. Also, I implemented the methods as extension methods.
Skilling 的代码将希尔伯特索引表示为转置，存储为数组.我发现交错位并形成单个 BigInteger 更方便(在字典中更有用，更容易在循环中迭代等)，但我优化了该操作及其逆运算，并使用幻数、位运算等，以及代码很长，所以我省略了它.
Skilling's code represents the Hilbert index as a transpose, stored as an array. I find it more convenient to interleave the bits and form a single BigInteger (more useful in Dictionaries, easier to iterate over in loops, etc), but I optimized that operation and its inverse with magic numbers, bit operations and the like, and the code is lengthy, so I have omitted it.
namespace HilbertExtensions
{
/// <summary>
/// Convert between Hilbert index and N-dimensional points.
///
/// 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]
///
/// NOTE: This algorithm is derived from work done by John Skilling and published in "Programming the Hilbert curve".
/// (c) 2004 American Institute of Physics.
///
/// </summary>
public static class HilbertCurveTransform
{
/// <summary>
/// Convert the Hilbert index into an N-dimensional point expressed as a vector of uints.
///
/// Note: In Skilling's paper, this function is named TransposetoAxes.
/// </summary>
/// <param name="transposedIndex">The Hilbert index stored in transposed form.</param>
/// <param name="bits">Number of bits per coordinate.</param>
/// <returns>Coordinate vector.</returns>
public static uint[] HilbertAxes(this uint[] transposedIndex, int bits)
{
var X = (uint[])transposedIndex.Clone();
int n = X.Length; // n: Number of dimensions
uint N = 2U << (bits - 1), P, Q, t;
int i;
// Gray decode by H ^ (H/2)
t = X[n - 1] >> 1;
// Corrected error in Skilling's paper on the following line. The appendix had i >= 0 leading to negative array index.
for (i = n - 1; i > 0; i--)
X[i] ^= X[i - 1];
X[0] ^= t;
// Undo excess work
for (Q = 2; Q != N; Q <<= 1)
{
P = Q - 1;
for (i = n - 1; i >= 0; i--)
if ((X[i] & Q) != 0U)
X[0] ^= P; // invert
else
{
t = (X[0] ^ X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
} // exchange
return X;
}
/// <summary>
/// Given the axes (coordinates) of a point in N-Dimensional space, find the distance to that point along the Hilbert curve.
/// That distance will be transposed; broken into pieces and distributed into an array.
///
/// The number of dimensions is the length of the hilbertAxes array.
///
/// Note: In Skilling's paper, this function is called AxestoTranspose.
/// </summary>
/// <param name="hilbertAxes">Point in N-space.</param>
/// <param name="bits">Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.</param>
/// <returns>The Hilbert distance (or index) as a transposed Hilbert index.</returns>
public static uint[] HilbertIndexTransposed(this uint[] hilbertAxes, int bits)
{
var X = (uint[])hilbertAxes.Clone();
var n = hilbertAxes.Length; // n: Number of dimensions
uint M = 1U << (bits - 1), P, Q, t;
int i;
// Inverse undo
for (Q = M; Q > 1; Q >>= 1)
{
P = Q - 1;
for (i = 0; i < n; i++)
if ((X[i] & Q) != 0)
X[0] ^= P; // invert
else
{
t = (X[0] ^ X[i]) & P;
X[0] ^= t;
X[i] ^= t;
}
} // exchange
// Gray encode
for (i = 1; i < n; i++)
X[i] ^= X[i - 1];
t = 0;
for (Q = M; Q > 1; Q >>= 1)
if ((X[n - 1] & Q)!=0)
t ^= Q - 1;
for (i = 0; i < n; i++)
X[i] ^= t;
return X;
}
}
}
我已将 C# 中的工作代码发布到 github.
I have posted working code in C# to github.
参见 https://github.com/paulchernoch/HilbertTransformation
更新:我刚刚在 crates.io 上发布了(2019 年秋季)一个名为hilbert"的 Rust crate.它还使用斯基林算法.请参阅 https://crates.io/crates/hilbert
UPDATE: I just published (Fall 2019) a Rust crate on crates.io called "hilbert". It also uses Skilling's algorithm. See https://crates.io/crates/hilbert
这篇关于将 N 维值映射到希尔伯特曲线上的一个点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
