Где я нахожусь: локализация точки полосами и трапециями

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

Как думаешь, как ответить быстро

Карта города разбита на районы, границы - отрезки. Каждую секунду приходят тысячи вызовов такси с координатами, и для каждого надо назвать район. Проверять точку в каждом многоугольнике лучом из главы о примитивах - это O(n) на запрос, где n - число рёбер всех районов.

Ответь себе до чтения: если карта не меняется, а запросов миллионы, сколько можно заплатить один раз за подготовку, чтобы каждый запрос стоил логарифм? И чем за это расплачиваются? Ответ урока: двумя двоичными поисками, а расплачиваются памятью. На «лестнице» из 60 отрезков структура хранит 1860 записей.

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

Библиотека, где книги разложены по стеллажам, а на каждом стеллаже - по полкам. Чтобы найти книгу, не перебирают все книги: сначала по каталогу находят стеллаж, потом на нём - полку. Два поиска вместо перебора. Но если одна огромная книга занимает полку на всех стеллажах сразу, её карточку приходится класть в каталог каждого стеллажа - и каталог разбухает.

Полосы - стеллажи, отрезки внутри полосы - полки, а длинный отрезок, пересекающий много полос, - та огромная книга.

Механизм: полосы и упорядоченные отрезки

Разбиение задано отрезками-границами, которые не пересекаются (кроме общих концов). Область, где лежит точка, определяется отрезком прямо под ней: он и область над ним - одно и то же, если у каждого отрезка записано, какая область над ним (массив «номер отрезка → район»). Точка ровно на границе - договор: в программе ниже «выше» строгое, и такая точка относится к району под границей, как будто лежит чуть ниже её. Если под точкой отрезков нет вовсе, она во внешней области; обычно карту обносят рамкой, и тогда снизу отрезок есть всегда.

Метод полос (Добкин и Липтон, 1976). Через каждый конец отрезка проводится вертикаль. Между двумя соседними вертикалями - полоса, в которой нет ни одного конца. Значит, каждый отрезок либо пересекает полосу целиком, либо не касается её, и отрезки внутри полосы не пересекаются - их можно упорядочить снизу вверх.

запрос q:
  1. двоичный поиск по x вертикалей - полоса, где q        O(log n)
  2. двоичный поиск по отрезкам полосы: поворот            O(log n)
     «q выше отрезка?» - последний отрезок ниже q

Почему память квадратична. Полос до 2n, и в каждой может оказаться порядка n отрезков - если длинные отрезки тянутся над множеством концов коротких. Тогда записей порядка n².

Трапецоидная карта (Сейдел, 1991; Малмулей, 1990) чинит память. Из каждого конца отрезка проводят вертикаль не через всю плоскость, а только до ближайших отрезков сверху и снизу. Плоскость режется на трапеции, и их не больше 3n + 1: левая стенка каждой трапеции проходит через конец отрезка; левый конец задаёт не больше двух трапеций, правый - одну, и ещё одна самая левая. Высоты отрезков в полосе программа сравнивает как дроби с положительными знаменателями - перекрёстным умножением. Отрезки вставляют по одному в случайном порядке, обновляя карту и граф поиска - ориентированный граф проверок «левее / правее конца» и «выше / ниже отрезка». Ожидаемая память O(n), ожидаемое время запроса O(log n), построение O(n · log n). Случайный порядок здесь не украшение: для любого фиксированного порядка найдётся вход, где граф поиска глубокий, а при случайном порядке вставки плохой граф получается с малой вероятностью - для любого входа.

Правило

Локализация: подготовить разбиение один раз, отвечать на запрос за логарифм. Метод полос - вертикали через все концы, в полосе отрезки по высоте, запрос - два двоичных поиска за O(log n), память до O(n²). Трапецоидная карта со случайным порядком вставки - ожидаемо O(n) памяти и O(log n) на запрос.

Разбор: карта и лестница

Программа строит полосы и отвечает на 200 случайных запросов, сверяя каждый ответ с перебором всех отрезков. Концы отрезков стоят на чётных x, запросы - на нечётных: так запрос не попадает на вертикаль, где отрезок начинается или кончается, и спора «принадлежит ли он полосе» не возникает.

def povorot(a, b, c):
    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])

def vyshe_v(s, t, x2):
    # выше ли отрезок s отрезка t в точке x = x2 / 2 (точно, в целых)
    a, b = s
    c, d = t
    ys = a[1] * 2 * (b[0] - a[0]) + (x2 - 2 * a[0]) * (b[1] - a[1])
    zs = 2 * (b[0] - a[0])
    yt = c[1] * 2 * (d[0] - c[0]) + (x2 - 2 * c[0]) * (d[1] - c[1])
    zt = 2 * (d[0] - c[0])
    return ys * zt > yt * zs

def postroit(otrezki):
    xs = []
    for s in otrezki:
        xs.append(s[0][0])
        xs.append(s[1][0])
    xs = sorted(set(xs))
    polosy = []
    for k in range(len(xs) - 1):
        x2 = xs[k] + xs[k + 1]
        v_polose = []
        for i in range(len(otrezki)):
            a, b = otrezki[i]
            if a[0] <= xs[k] and b[0] >= xs[k + 1]:
                # вставка по высоте в середине полосы
                j = len(v_polose)
                v_polose.append(i)
                while j > 0:
                    if not vyshe_v(otrezki[v_polose[j - 1]], [a, b], x2):
                        break
                    v_polose[j] = v_polose[j - 1]
                    j = j - 1
                v_polose[j] = i
        polosy.append(v_polose)
    return [xs, polosy]

def nayti(struktura, otrezki, q, schet):
    # номер отрезка прямо под точкой q или -1
    xs, polosy = struktura
    lo = 0
    hi = len(xs) - 1
    if q[0] <= xs[0] or q[0] >= xs[-1]:
        return -1
    while hi - lo > 1:
        schet[0] = schet[0] + 1
        mid = (lo + hi) // 2
        if xs[mid] <= q[0]:
            lo = mid
        else:
            hi = mid
    sp = polosy[lo]
    a = 0
    b = len(sp)
    while a < b:
        schet[0] = schet[0] + 1
        m = (a + b) // 2
        s = otrezki[sp[m]]
        if povorot(s[0], s[1], q) > 0:
            a = m + 1
        else:
            b = m
    if a == 0:
        return -1
    return sp[a - 1]
def perebor(otrezki, q):
    luchshij = -1
    for i in range(len(otrezki)):
        a, b = otrezki[i]
        if a[0] < q[0] < b[0] and povorot(a, b, q) > 0:
            if luchshij == -1:
                luchshij = i
            elif vyshe_v([a, b], otrezki[luchshij], 2 * q[0]):
                luchshij = i
    return luchshij

import random
random.seed(6)
# 1) карта из коротких отрезков: каждый в своей клетке сетки 6 × 6
karta = []
for i in range(6):
    for j in range(6):
        # концы на чётных x, запросы - на нечётных:
        # запрос никогда не встаёт на границу полосы
        x1 = 100 * i + 2 * random.randint(1, 15)
        x2 = 100 * i + 2 * random.randint(30, 49)
        y1 = 100 * j + random.randint(1, 99)
        karta.append([[x1, y1], [x2, 100 * j + random.randint(1, 99)]])
# 2) лестница: 30 длинных отрезков и 30 коротких под ними
lestnica = []
for k in range(30):
    lestnica.append([[0, 100 + 10 * k], [1000, 101 + 10 * k]])
for k in range(30):
    lestnica.append([[10 + 30 * k, 0], [20 + 30 * k, 5]])
for otrezki in [karta, lestnica]:
    s = postroit(otrezki)
    pamyat = 0
    for p in s[1]:
        pamyat = pamyat + len(p)
    sovpalo = 0
    schet = [0]
    for t in range(200):
        q = [2 * random.randint(0, 499) + 1, random.randint(1, 599)]
        if nayti(s, otrezki, q, schet) == perebor(otrezki, q):
            sovpalo = sovpalo + 1
    print(len(otrezki), len(s[1]), pamyat, sovpalo, schet[0] / 200)

36 58 174 200 4.69
60 61 1860 200 10.765

Столбцы: число отрезков, число полос, записей во всех полосах, совпавших с перебором ответов, шагов поиска на запрос. Карта из 36 коротких отрезков: 58 полос и 174 записи - по три отрезка на полосу, память почти линейна. Запрос - около пяти шагов поиска.

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

Вторая строка - «лестница»: 30 длинных отрезков через весь экран и 30 коротких под ними. Концы коротких нарезали 61 полосу, и в каждой лежат все 30 длинных отрезков: 1860 записей на 60 отрезков. Удвой число ступенек - записей станет вчетверо больше: это n². Шагов на запрос около одиннадцати - log₂ 61 + log₂ 31 ≈ 6 + 5, логарифм держится, разбухает только память.

На такой лестнице трапецоидная карта хранит не больше 3 · 60 + 1 = 181 трапеции: вертикаль из конца короткого отрезка упирается в первый же длинный отрезок над ним и не идёт выше. Именно поэтому в геоинформационных системах полосы на практике не строят, а используют трапецоидные карты или деревья с ограничивающими прямоугольниками.

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

Запрос на вертикали полосы. Отрезок, начавшийся ровно в x запроса, одна функция считает принадлежащим полосе, другая - нет. Нужен договор: полуоткрытые полосы, как полуоткрытые интервалы в главе о заметающей прямой.

Порядок отрезков в полосе по высоте в float. Два близких отрезка могут поменяться местами, и запрос вернёт не тот район. Сравнение в середине полосы делается в целых - удвоением координаты середины.

Фиксированный порядок вставки в трапецоидную карту. Например, по длине или по x. Для такого порядка существует вход, где граф поиска вырождается в цепочку, и запрос стоит n.

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

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

Теперь сам

1. Отрезков 1000, и все они короткие и не перекрываются по x. Сколько примерно записей в полосах?
2. Почему для запроса хватает отрезка прямо под точкой, а не всего многоугольника района?
3. Зачем трапецоидной карте случайный порядок вставки, если результат - сама карта - от порядка не зависит?