Я наконец сломался и выложил немного денег. В AIP (Американский институт физики) есть хорошая короткая статья с исходным кодом на языке C. «Программирование кривой Гильберта» Джона Скиллинга (из AIP Conf. Proc. 707, 381 (2004)), а также приложение с кодом для отображения в обоих направлениях. Он работает для любого числа измерений> 1, не является рекурсивным, не использует таблицы поиска с изменением состояния, которые поглощают огромные объемы памяти, и в основном использует битовые операции. Таким образом, он достаточно быстрый и имеет хороший объем памяти.
Если вы решили приобрести статью, я обнаружил ошибку в исходном коде.
Следующая строка кода (находится в функции TransposetoAxes) имеет ошибку:
для (i = n-1; i> = 0; i--) X [i] ^ = X [i-1];
Корректировка заключается в изменении значения больше или равно (> =) на значение больше (>). Без этого исправления доступ к массиву X осуществляется с использованием отрицательного индекса, когда переменная «i» становится равной нулю, что приводит к сбою программы.
Я рекомендую прочитать статью (которая состоит из семи страниц, включая код), поскольку она объясняет, как работает алгоритм, что далеко не очевидно.
Я перевел его код на C # для собственного использования. Код следует. Skilling выполняет преобразование на месте, перезаписывая вектор, который вы передаете. Я решил сделать клон входного вектора и вернуть новую копию. Также я реализовал методы как методы расширения.
Код Скиллинга представляет индекс Гильберта в виде транспонирования, сохраненного в виде массива. Я считаю более удобным чередовать биты и формировать один BigInteger (более полезный в словарях, легче перебирать циклы и т. Д.), Но я оптимизировал эту операцию и ее инверсию с помощью магических чисел, битовых операций и т. код длинный, поэтому я его опустил.
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.
См. https://github.com/paulchernoch/HilbertTransformation