Отображение N-мерного значения в точку на кривой Гильберта

Ответ 1

Я, наконец, сломался и вырубил деньги. AIP (Американский физический институт) имеет красивую короткую статью с исходным кодом на C. "Программирование кривой Гильберта" Джона Скиллинга (из AIP Conf. Proc. 707, 381 (2004)) имеет приложение с кодом для отображения в обоих направлениях. Он работает для любого количества измерений > 1, не является рекурсивным, не использует таблицы поиска состояния, которые поглощают огромные объемы памяти и в основном используют битовые операции. Таким образом, он достаточно быстр и имеет хорошую память.

Если вы решили приобрести статью, я обнаружил ошибку в исходном коде.

Следующая строка кода (найденная в функции TransposetoAxes) имеет ошибку:

для (i = n-1; я >= 0; я -) X [i] ^ = X [i-1];

Поправка заключается в изменении большего или равного ( > =) на большее, чем ( > ). Без этой коррекции массив X обращается с использованием отрицательного индекса, когда переменная "i" становится нулевой, что приводит к сбою программы.

Я рекомендую прочитать статью (которая длится семь страниц, включая код), поскольку она объясняет, как работает алгоритм, что далеко не очевидно.

Я перевел свой код на С# для собственного использования. Код следует. Skilling выполняет преобразование на месте, перезаписывая вектор, который вы проходите. Я решил сделать клон входного вектора и вернуть новую копию. Кроме того, я реализовал методы как методы расширения.

Skilling code представляет индекс Гильберта как транспозицию, хранящуюся как массив. Мне удобнее чередоваться с битами и формировать один 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 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 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 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;
        }

    }
}

Я отправил рабочий код в С# в github.

См. https://github.com/paulchernoch/HilbertTransformation

Ответ 2

Алгоритм отображения из n- > 1 и 1- > n, приведенный здесь "Вычисление сопоставлений между одним и n-мерными значениями с использованием критической заполняющей гильбертовой кривой" J K Lawder

Если вы Google для "SFC module и Kademlia overlay", вы найдете группу, которая заявляет, что использует ее как часть своей системы. Если вы просмотрите источник, вы, вероятно, можете извлечь соответствующую функцию.

Ответ 3

Мне не ясно, как это будет делать то, что вы хотите. Рассмотрим этот трехмерный трехмерный случай:

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

который может быть "гильбертизирован" по следующему пути:

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

в 1-й порядок:

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

Вот неприятный бит. Рассмотрим список пар и 1D расстояния ниже:

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

Этот эффект растет экспоненциально с количеством измерений (при условии, что каждый размер имеет тот же "размер" ).

Ответ 4

Еще одна возможность заключалась бы в создании kd-tree ваших данных, а затем на обход дерева по порядку заказ. Построение kd-дерева требует только того, чтобы у вас был хороший медианный алгоритм поиска, которого много.

Ответ 5

Я потратил немного времени на перевод кода Павла Чернова на Java и его очистку. Возможно, в моем коде есть ошибка, особенно потому, что у меня нет доступа к тому документу, из которого он был изначально. Тем не менее, он передает те модульные тесты, которые я смог написать. Он ниже.

Обратите внимание, что я оценил как Z-Order, так и кривые Гильберта для пространственного индексирования на довольно больших наборах данных. Должен сказать, что Z-Order обеспечивает гораздо лучшее качество. Но не стесняйтесь попробовать сами.

    /**
     * Convert the Hilbert index into an N-dimensional point expressed as a vector of uints.
     *
     * Note: In Skilling paper, this function is named TransposetoAxes.
     * @param transposedIndex The Hilbert index stored in transposed form.
     * @param bits Number of bits per coordinate.
     * @return Point in N-space.
     */
    static long[] HilbertAxes(final long[] transposedIndex, final int bits) {
        final long[] result = transposedIndex.clone();
        final int dims = result.length;
        grayDecode(result, dims);
        undoExcessWork(result, dims, bits);
        return result;
    }

    static void grayDecode(final long[] result, final int dims) {
        final long swap = result[dims - 1] >>> 1;
        // Corrected error in Skilling paper on the following line. The appendix had i >= 0 leading to negative array index.
        for (int i = dims - 1; i > 0; i--)
            result[i] ^= result[i - 1];
        result[0] ^= swap;
    }

    static void undoExcessWork(final long[] result, final int dims, final int bits) {
        for (long bit = 2, n = 1; n != bits; bit <<= 1, ++n) {
            final long mask = bit - 1;
            for (int i = dims - 1; i >= 0; i--)
                if ((result[i] & bit) != 0)
                    result[0] ^= mask; // invert
                else
                    swapBits(result, mask, i);
        }
    }

    /**
     * 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 paper, this function is called AxestoTranspose.
     * @param hilbertAxes Point in N-space.
     * @param bits Depth of the Hilbert curve. If bits is one, this is the top-level Hilbert curve.
     * @return The Hilbert distance (or index) as a transposed Hilbert index.
     */
    static long[] HilbertIndexTransposed(final long[] hilbertAxes, final int bits) {
        final long[] result = hilbertAxes.clone();
        final int dims = hilbertAxes.length;
        final long maxBit = 1L << (bits - 1);
        inverseUndo(result, dims, maxBit);
        grayEncode(result, dims, maxBit);
        return result;
    }

    static void inverseUndo(final long[] result, final int dims, final long maxBit) {
        for (long bit = maxBit; bit != 0; bit >>>= 1) {
            final long mask = bit - 1;
            for (int i = 0; i < dims; i++)
                if ((result[i] & bit) != 0)
                    result[0] ^= mask; // invert
                else
                    swapBits(result, mask, i);
        } // exchange
    }

    static void grayEncode(final long[] result, final int dims, final long maxBit) {
        for (int i = 1; i < dims; i++)
            result[i] ^= result[i - 1];
        long mask = 0;
        for (long bit = maxBit; bit != 0; bit >>>= 1)
            if ((result[dims - 1] & bit) != 0)
                mask ^= bit - 1;
        for (int i = 0; i < dims; i++)
            result[i] ^= mask;
    }

    static void swapBits(final long[] array, final long mask, final int index) {
        final long swap = (array[0] ^ array[index]) & mask;
        array[0] ^= swap;
        array[index] ^= swap;
    }

Ответ 6

Я не вижу, как вы можете использовать кривую Гильберта в одном измерении.

Если вы заинтересованы в сопоставлении точек с меньшей размерностью при сохранении расстояний (с минимальной ошибкой), вы можете посмотреть алгоритмы "Многомерное масштабирование".

Имитированный отжиг - это один подход.

Изменить: Спасибо за комментарий. Теперь я понимаю, что вы подразумеваете под подходом Гильберта Кривой. Тем не менее, это сложная проблема, и, учитывая N = 100 и 10 миллионов точек данных, я не думаю, что какой-либо подход будет хорошо сохранять местность и работать в разумные сроки. Я не думаю, что здесь будут работать kd-деревья.

Если поиск полного заказа не важен для вас, вы можете посмотреть на хэширование на основе местоположения и другие приближенные схемы ближайших соседей. Иерархическое многомерное масштабирование с ведрами точек для уменьшения размера ввода может дать вам хороший порядок, но опять же сомнительно в таком высоком размере.