Z-порядок и кривая Гильберта: плоскость в одном числе

Если ты перемешаешь биты двух координат в один номер, объяснишь, почему Z-порядок иногда прыгает через полсетки, и скажешь, чем кривая Гильберта за это платит, - листай до «Теперь сам».

Как думаешь, можно ли отсортировать точки плоскости

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

Прежде чем читать дальше, ответь себе: что будет, если отсортировать точки по x, а при равных - по y? Запрос узкой вертикальной полосы станет быстрым. Запрос небольшого квадрата - нет: в диапазон по x попадут все точки полосы от нижнего края карты до верхнего. Нужен порядок, в котором обе координаты равноправны.

На что это похоже

Библиотекарь расставляет книги, у каждой два признака: тема и год. Если сортировать по теме, а потом по году, все книги одного года разбросаны по всему залу. Можно иначе: сначала делить зал пополам по теме, каждую половину - пополам по году, каждую четверть - снова по теме, и так далее. Книга, близкая и по теме, и по году, окажется в том же углу зала.

Номер места при таком делении - последовательность ответов «левая или правая половина», в которой вопросы про тему и про год идут через один. В битах: биты одной координаты чередуются с битами другой.

Механизм: перемешивание битов и поворот четвертей

Z-порядок. Ось x направлена вправо, ось y - вниз, как строки на экране. Квадрат делится на четыре четверти, каждая четверть - снова на четыре. Номер клетки собирается из битов координат: бит i числа x идёт на место 2i, бит i числа y - на место 2i + 1.

x = 3 = 0 1 1                  y = 5 = 1 0 1
номер: y₂ x₂ y₁ x₁ y₀ x₀  =  1 0 0 1 1 1  =  39

Чтобы раздвинуть биты x по чётным местам, не нужен цикл по битам. Работает лестница масок из главы про счёт единиц, пройденная в обратную сторону: сначала раздвинуть байты, потом четвёрки, пары, одиночные биты.

В каждой строке x - результат предыдущей:

x                     = 0000 0000 0000 0000 abcd efgh ijkl mnop
(x | x << 8) & 00FF…  = 0000 0000 abcd efgh 0000 0000 ijkl mnop
(x | x << 4) & 0F0F…  = 0000 abcd 0000 efgh 0000 ijkl 0000 mnop
(x | x << 2) & 3333…  = 00ab 00cd 00ef 00gh 00ij 00kl 00mn 00op
(x | x << 1) & 5555…  = 0a0b 0c0d 0e0f 0g0h 0i0j 0k0l 0m0n 0o0p

Обход клеток в порядке номеров рисует букву Z в каждой четвёрке клеток - отсюда название. Его слабое место - переходы между четвертями: с правого нижнего угла одной четверти путь прыгает в левый верхний угол следующей.

Кривая Гильберта устраняет прыжки. Она тоже обходит четверти по очереди, но в каждой четверти рисунок повёрнут или отражён так, чтобы конец обхода одной четверти примыкал к началу следующей. Номер клетки по-прежнему читается парами битов - пара выбирает четверть, - но после каждой пары координаты приходится поворачивать. Простой перестановкой битов её не получить; нужен цикл по парам.

Правило

Z-номер клетки (x, y) - биты x на чётных местах, биты y на нечётных: razdvinut(x) | razdvinut(y) << 1. Клетки с общим началом номера длиной 2k битов лежат в одном квадрате со стороной 2^(m − k), где m - число битов координаты. У кривой Гильберта, кроме того, клетки с соседними номерами всегда соседние по стороне; у Z-порядка - нет. Обратное неверно для обеих: соседние клетки могут иметь далёкие номера.

Разбор: Z-номер и первые шестнадцать клеток кривой Гильберта

def razdvinut(x):
    x = (x | x << 8) & 0x00FF00FF
    x = (x | x << 4) & 0x0F0F0F0F
    x = (x | x << 2) & 0x33333333
    x = (x | x << 1) & 0x55555555
    return x

def z_kod(x, y):
    return razdvinut(x) | razdvinut(y) << 1

print(bin(razdvinut(0b1011)), z_kod(3, 5), bin(z_kod(3, 5)))
print(z_kod(0, 0), z_kod(1, 0), z_kod(0, 1), z_kod(1, 1), z_kod(2, 0))

0b1000101 39 0b100111
0 1 2 3 4

1011 раздвинулось в 1 0 0 0 1 0 1: между битами встали нули. Вторая строка - начало обхода: клетки (0,0), (1,0), (0,1), (1,1) получают номера 0-3 - это и есть буква Z, - а номер 4 достаётся клетке (2,0): из (1,1) путь ушёл по диагонали.

Для кривой Гильберта показано преобразование из номера в координаты; обратное устроено так же - циклом по парам битов. Здесь n - сторона сетки, степень двойки, d - номер клетки от 0 до n² − 1:

def gilbert(n, d):
    x = 0
    y = 0
    s = 1
    while s < n:
        rx = d >> 1 & 1
        ry = (d ^ rx) & 1
        if ry == 0:
            if rx == 1:
                x = s - 1 - x
                y = s - 1 - y
            t = x
            x = y
            y = t
        x = x + s * rx
        y = y + s * ry
        d = d >> 2
        s = s << 1
    return [x, y]

put = []
for d in range(16):
    put.append(gilbert(4, d))
print(put)

[[0, 0], [1, 0], [1, 1], [0, 1], [0, 2], [0, 3], [1, 3], [1, 2], [2, 2], [2, 3], [3, 3], [3, 2], [3, 1], [2, 1], [2, 0], [3, 0]]

Функция читает номер по два бита, с младших: rx и ry говорят, в какую четверть текущего квадрата попадает точка. Четверти с ry = 0 перед этим поворачиваются - координаты меняются местами, а в одной из них ещё и отражаются. s - сторона уже собранного квадрата, она удваивается на каждом круге.

Те же номера на сетке 4 на 4, y растёт вниз:

Z-порядок            кривая Гильберта
 0  1  4  5            0  1 14 15
 2  3  6  7            3  2 13 12
 8  9 12 13            4  7  8 11
10 11 14 15            5  6  9 10

У кривой Гильберта каждая следующая клетка отличается от предыдущей на единицу в одной координате: из (0,1) путь идёт в (0,2), а не прыгает, как у Z.

Разбор: длина пути на сетке 16 на 16, и чего не обещает ни одна кривая

Обойдём все 256 клеток в порядке номеров и сложим расстояния между соседними по номеру клетками - расстояние считается по клеткам, по горизонтали плюс по вертикали. Клетку по Z-номеру находят обратным действием: чётные биты номера - это x, нечётные - y.

порядок            длина пути    шагов длиннее одной клетки    самый длинный шаг
кривая Гильберта      255                  0                          1
Z-порядок             480                127                         16

У кривой Гильберта все 255 шагов - в соседнюю клетку. У Z-порядка половина шагов - скачки, и самый длинный, через середину сетки, - шестнадцать клеток. Для обхода, при котором данные читаются с диска или из памяти, это разница между «следующее рядом» и «следующее где-то».

За непрерывность кривая Гильберта платит ценой вычисления: Z-номер - четыре строки без цикла на координату, а в процессоре - одна команда PDEP, раскладывающая биты по местам, отмеченным маской; номер на кривой Гильберта - цикл по парам битов с ветвлениями.

Теперь граница, общая для обеих. Надёжное обещание одно: общее начало номера - общий квадрат. А вот «соседние клетки имеют близкие номера» выполняется часто, но не всегда, и всегда выполняться не может: на линии у клетки два соседа, на плоскости - четыре, и как ни укладывай плоскость в линию, у некоторых соседних клеток номера окажутся далеко. У кривой Гильберта на сетке 4 на 4 клетки (0,0) и (0,1) соседние, а номера у них - 0 и 3; клетки (1,0) и (2,0) - соседние с номерами 1 и 14. Поэтому запрос «всё в прямоугольнике» превращается не в один диапазон номеров, а в несколько, и хороший порядок - тот, при котором диапазонов мало.

Что ломается без этого

Первое: сортировка по одной координате. Индекс по x отвечает на запрос квадрата, читая всю вертикальную полосу. На карте мира это миллионы лишних точек на каждый запрос.

Второе: «соседние номера - значит рядом» для Z. На скачке между четвертями соседние номера лежат в разных концах сетки. Алгоритм, который полагается на близость соседей по Z-номеру, ошибается на каждом втором шаге.

Третье: «рядом - значит соседние номера». Не выполняется ни для какой кривой. Поиск ближайшей точки одним заглядыванием к соседям по номеру пропускает точки через границу четвертей.

Четвёртое: координата длиннее 16 битов. razdvinut рассчитан на 16 битов: старшие биты более длинной координаты не пропадают, а ложатся на чужие места - razdvinut(0x10000) и razdvinut(0x100) дают одно и то же, и две разные клетки получают один номер. Для 32-битных координат нужны 64-битные маски и пятая строка.

Зачем это в работе

Геохеш - строка вроде ucfv0 для точки на карте - это Z-номер широты и долготы, записанный буквами. Общее начало двух геохешей означает общий квадрат, и поиск «рядом» в обычной базе становится поиском по началу строки - с оговоркой: соседние точки по разные стороны границы квадрата имеют разные начала, поэтому смотрят ещё и восемь соседних квадратов. Библиотека S2, на которой работают карты и геозапросы в крупных базах данных, нумерует клетки земной поверхности кривой Гильберта - ради отсутствия скачков.

Аналитические базы упорядочивают строки таблицы по Z-номеру от нескольких столбцов сразу: запрос с условием на любой из них читает меньше файлов. Это та же идея равноправных координат, только координаты - дата, идентификатор клиента и регион.

В графических процессорах точки текстуры хранят в Z-порядке: соседние на картинке точки попадают в одну строку кеша. А обход матрицы по кривой Гильберта используют там, где важно, чтобы каждый следующий элемент был рядом с предыдущим, - от умножения больших матриц до растрового сглаживания изображений.

Теперь сам

1. Найди Z-номер клетки (5, 3).
2. Какой клетке соответствует Z-номер 27?
3. Раздвинь x = 0b0110 и выпиши результат в битах.
4. По списку из разбора: какой номер на кривой Гильберта у клетки (3, 3) и у её соседки (3, 2)? А у соседки (2, 3)?
5. Две точки имеют Z-номера с общими старшими шестью битами из восьми. В квадрате какого размера они обе лежат, если координаты четырёхбитные?