将N维值映射到希尔伯特曲线上的点

我有一组巨大的N维点(数千万; N接近100).

我需要将这些点映射到单个维度,同时保留空间局部性.我想用希尔伯特空间填充曲线来做.

对于每个点,我想选择曲线上最近的点.该点的希尔伯特值(从曲线起点到拾取点的曲线长度)是我寻求的单维值.

计算不一定是即时的,但我希望它在不错的现代家用PC硬件上不会超过几个小时.

有关实施的建议吗？有没有可以帮助我的图书馆？(语言并不重要.)

我终于崩溃并取出了一些钱.AIP(美国物理学会)有一篇很好的短篇文章,其源代码在C中."John Hilling编写Hilbert曲线"(来自AIP Conf.Proc.707,381(2004))有一个代码为两个方向的映射.它适用于任意数量的维度> 1,不是递归的,不使用吞噬大量内存的状态转换查找表,并且主要使用位操作.因此它相当快并且具有良好的存储器占用空间.

如果您选择购买该文章,我在源代码中发现了一个错误.

以下代码行(在函数TransposetoAxes中找到)有错误:

for(i = n-1; i> = 0; i--)X [i] ^ = X [i-1];

校正是将大于或等于(> =)更改为大于(>).如果没有这种校正,当变量"i"变为零时,使用负索引访问X数组,从而导致程序失败.

我建议阅读这篇文章(长达七页,包括代码),因为它解释了算法是如何工作的,这是非常明显的.

我将他的代码翻译成C#供我自己使用.代码如下.Skilling执行转换,覆盖您传入的向量.我选择复制输入向量并返回一个新副本.此外,我将方法实现为扩展方法.

Skilling的代码将Hilbert索引表示为转置,存储为数组.我发现交错位并形成一个BigInteger更方便(在字典中更有用,更容易在循环中迭代等),但我优化了该操作,并且它与幻数,位操作等相反,并且代码很冗长,所以我省略了它.

namespace HilbertExtensions
{
    /// 

    /// 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.
    /// 
    /// 
    public static class HilbertCurveTransform
    {
        /// 

        /// 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.
        /// 
        /// The Hilbert index stored in transposed form.
        /// Number of bits per coordinate.
        /// Coordinate vector.
        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;
        }

        /// 

        /// 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.
        /// 
        /// Point in N-space.
        /// Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.
        /// The Hilbert distance (or index) as a transposed Hilbert index.
        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.

请参阅https://github.com/paulchernoch/HilbertTransformation

这里给出的从n-> 1和1-> n映射的算法 "使用希尔伯特空间填充曲线计算一维和n维值之间的映射"JK Lawder

如果您使用Google"SFC模块和Kademlia叠加层",您会发现一个声称将其用作系统一部分的群组.如果查看源代码,则可以提取相关函数.

我不清楚这将如何执行您想要的操作。考虑一下这种普通的3D情况：

001 ------ 101
 |\         |\
 | \        | \
 |  011 ------ 111
 |   |      |   |
 |   |      |   |
000 -|---- 100  |
  \  |       \  |
   \ |        \ |
    010 ------ 110

可以通过以下路径进行“ Hilbertized”：

001 -----> 101
  \          \
   \          \
    011        111
     ^          |
     |          |
000  |     100  |
  \  |       \  |
   \ |        \ V
    010        110

进入一维订单：

000 -> 010 -> 011 -> 001 -> 101 -> 111 -> 110 -> 100

这是讨厌的地方。考虑下面的线对和一维距离列表：

000 : 100 -> 7
010 : 110 -> 5
011 : 111 -> 3
001 : 101 -> 1

在所有情况下，左手值和右手值都是相同的3D距离（在第一位置为+/- 1），这似乎暗示了相似的“空间局部性”。但是通过任何尺寸顺序选择（在上面的示例中依次为y，z，z）进行线性化都会破坏该局部性。

换一种说法是，采用起点并按距起点的距离对其余点进行排序将提供截然不同的结果。以000开始为例，例如：

1D ordering : distance    3D ordering : distance
----------------------    ----------------------
        010 : 1           001,010,100 : 1
                          011,101,110 : sqrt(2)
                              111     : sqrt(3)
        011 : 2
        001 : 3
        101 : 4
        111 : 5
        110 : 6
        100 : 7

该效果随维数成倍增长（假设每个维具有相同的“大小”）。