Все пересечения: события-точки и алгоритм Бентли-Оттмана

Если ты знаешь, как устроена очередь событий с точками пересечения, что происходит со статусом в такой точке и почему время (n + k) · log n, листай до «Теперь сам».

Как думаешь, когда это выгодно

Сорок отрезков - пар 780. Прошлый урок отвечал «есть ли пересечение» за n · log n, останавливаясь на первом. Теперь нужны все: например, чтобы построить наложение двух карт или найти все перекрёстки дорог.

Ответь себе до чтения: сколько проверок понадобится алгоритму, который находит все пересечения, если пересечений 11? А если 149? Всегда ли он лучше перебора 780 пар? Ответ: 59 проверок и 377 проверок. Чем больше пересечений, тем меньше выигрыш, - потому что каждое пересечение - отдельное событие со своей ценой.

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

Снова шоссе с полосами. Теперь машины могут обгонять друг друга - в момент обгона две машины в соседних полосах меняются местами. Регулировщик записывает каждый обгон. Обгон случается только между соседями, поэтому, заметив, что две соседние машины сближаются, он заранее отмечает в расписании, когда и где они поравняются. В этот момент он меняет их местами в своём списке и смотрит, к кому теперь подбираются новые соседи.

Расписание - очередь событий, список машин по полосам - статус, обгон - точка пересечения.

Механизм: очередь, статус и обмен соседей

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

Три вида событий.

левый конец s:     вставить s в статус; проверить s с соседями
правый конец s:    удалить s; проверить бывших соседей s друг с другом
точка s и t:       записать пересечение; поменять s и t местами;
                   проверить нижнего из них с соседом снизу,
                   верхнего - с соседом сверху

«Проверить» значит: если пара пересекается правее прямой и этой точки ещё нет в очереди - положить её в очередь. Одна и та же пара может стать соседями дважды - разойтись и снова сойтись, - поэтому найденные пары запоминаются.

Почему ничего не теряется. Тот же довод, что в прошлом уроке, но для каждого пересечения по очереди: перед точкой пересечения s и t обязательно становятся соседями на каком-то событии - и в этот момент точка кладётся в очередь.

Цена. Событий 2n + k. Каждое - извлечение из очереди и несколько операций со статусом и очередью; в сбалансированных деревьях это log n каждая (в очереди не больше 2n + k элементов, а log(n²) - это 2 log n). Итого (n + k) · log n. Нижняя граница для этой задачи - n · log n + k (в той же модели алгебраических решений, что у оболочки), и её достигают более сложные алгоритмы (Шазелл и Эдельсбруннер, 1992).

Правило

Все пересечения n отрезков: заметающая прямая с очередью событий по x и статусом по высоте. Концы отрезков - события заранее, точки пересечения - добавляются, когда пересекающиеся отрезки становятся соседями. В точке соседи меняются местами, и проверяются новые пары соседей. Время O((n + k) · log n), где k - число пересечений.

Разбор: четыре отрезка и журнал событий

Программа ниже - алгоритм целиком, на точной арифметике дробей. Она печатает журнал на четырёх отрезках, а потом проводит опыт, о котором следующий разбор. Учебные допущения: вертикальных отрезков нет, общих концов нет, никакие три отрезка не проходят через одну точку и никакие два события не делят x. Генератор nabor этого сам не гарантирует: набор с зерном 158 проверен отдельно, и у него все x событий различны. Статус и очередь - обычные списки; в настоящей реализации - сбалансированные деревья, а ключ пары min · 100000 + max годится, пока отрезков меньше ста тысяч.

import random

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

def ranshe(e, f):
    # событие: [x, знаменатель x, y, знаменатель y, вид, i, j]
    r = e[0] * f[1] - f[0] * e[1]
    if r != 0:
        return r < 0
    return e[2] * f[3] < f[2] * e[3]

def v_ochered(ochered, e):
    lo = 0
    hi = len(ochered)
    while lo < hi:
        mid = (lo + hi) // 2
        if ranshe(ochered[mid], e):
            lo = mid + 1
        else:
            hi = mid
    ochered.insert(lo, e)

def tochka(s, t):
    # точка пересечения отрезков s и t дробью или [] - если её нет
    a, b = s
    c, d = t
    if povorot(c, d, a) * povorot(c, d, b) >= 0:
        return []
    if povorot(a, b, c) * povorot(a, b, d) >= 0:
        return []
    r = [b[0] - a[0], b[1] - a[1]]
    q = [d[0] - c[0], d[1] - c[1]]
    zn = r[0] * q[1] - r[1] * q[0]
    ch = (c[0] - a[0]) * q[1] - (c[1] - a[1]) * q[0]
    if zn < 0:
        zn = -zn
        ch = -ch
    return [a[0] * zn + r[0] * ch, zn, a[1] * zn + r[1] * ch, zn]

def proverit(sost, u, v, x):
    # соседи status[u] (ниже) и status[v] (выше): пересекутся ли правее x
    otrezki, status, ochered, uzhe, schet = sost
    if u < 0 or v >= len(status):
        return
    s = status[u]
    t = status[v]
    schet[0] = schet[0] + 1
    p = tochka(otrezki[s], otrezki[t])
    klyuch = min(s, t) * 100000 + max(s, t)
    if p and p[0] * x[1] > x[0] * p[1] and klyuch not in uzhe:
        uzhe[klyuch] = True
        v_ochered(ochered, p + [1, s, t])

def bentli_ottman(otrezki, zhurnal):
    ochered = []
    for i in range(len(otrezki)):
        a, b = otrezki[i]
        v_ochered(ochered, [a[0], 1, a[1], 1, 0, i, -1])
        v_ochered(ochered, [b[0], 1, b[1], 1, 2, i, -1])
    status = []
    naideno = []
    uzhe = {}
    schet = [0]
    sost = [otrezki, status, ochered, uzhe, schet]
    while len(ochered) > 0:
        e = ochered.pop(0)
        x = [e[0], e[1]]
        if e[4] == 0:
            p = otrezki[e[5]][0]
            lo = 0
            hi = len(status)
            while lo < hi:
                mid = (lo + hi) // 2
                a, b = otrezki[status[mid]]
                if povorot(a, b, p) > 0:
                    lo = mid + 1
                else:
                    hi = mid
            status.insert(lo, e[5])
            proverit(sost, lo - 1, lo, x)
            proverit(sost, lo, lo + 1, x)
        elif e[4] == 2:
            k = status.index(e[5])
            status.pop(k)
            proverit(sost, k - 1, k, x)
        else:
            naideno.append([e[5], e[6]])
            k = status.index(e[5])
            status[k] = e[6]
            status[k + 1] = e[5]
            proverit(sost, k - 1, k, x)
            proverit(sost, k + 1, k + 2, x)
        if zhurnal:
            print(["начало", "точка", "конец"][e[4]], e[5:], status)
    return [naideno, schet[0]]

def nabor(n, L, M):
    # n отрезков на поле M × M, длина по x не больше L
    xs = list(range(M))
    random.shuffle(xs)
    otrezki = []
    for i in range(n):
        x1 = min(xs[2 * i], xs[2 * i + 1])
        x2 = max(xs[2 * i], xs[2 * i + 1])
        if x2 - x1 > L:
            x2 = x1 + (x2 - x1) % L + 1
        y1 = random.randint(0, M)
        otrezki.append([[x1, y1], [x2, random.randint(0, M)]])
    return otrezki

otrezki = [[[0, 2], [10, 8]], [[1, 7], [9, 1]], [[2, 4], [12, 5]],
           [[3, 9], [11, 6]]]
print(bentli_ottman(otrezki, True))

random.seed(158)
for L in [30, 400]:
    otrezki = nabor(40, L, 500)
    naideno, proverok = bentli_ottman(otrezki, False)
    n = len(otrezki)
    print(L, n, len(naideno), proverok, n * (n - 1) // 2)

начало [0, -1] [0]
начало [1, -1] [0, 1]
начало [2, -1] [0, 2, 1]
начало [3, -1] [0, 2, 1, 3]
точка [0, 2] [2, 0, 1, 3]
точка [0, 1] [2, 1, 0, 3]
точка [2, 1] [1, 2, 0, 3]
точка [0, 3] [1, 2, 3, 0]
конец [1, -1] [2, 3, 0]
конец [0, -1] [2, 3]
конец [3, -1] [2]
конец [2, -1] []
[[[0, 2], [0, 1], [2, 1], [0, 3]], 9]
30 40 11 59 780
400 40 149 377 780

Первые четыре строки - вставки; статус 0, 2, 1, 3 снизу вверх. Дальше идут четыре точки пересечения, и в каждой два соседа меняются местами: 0 и 2, потом 0 и 1, потом 2 и 1, потом 0 и 3. Обрати внимание на порядок: пересечения выдаются слева направо, а не в порядке номеров. Пересечение 0 и 3 нашлось, только когда 0 поднялся до соседства с 3 - после трёх обменов.

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

Разбор: когда пересечений мало и когда много

Последние две строки вывода - 40 отрезков на поле 500 × 500, сначала короткие (длина по x до 30), потом длинные (до 400). Столбцы: ограничение длины, число отрезков, число пересечений, число проверок, число пар для перебора.

Короткие отрезки: 11 пересечений, 59 проверок против 780 пар - в тринадцать раз меньше. Длинные: 149 пересечений, 377 проверок - уже только вдвое меньше пар. Проверок не больше 3n + 2k: вставка даёт две, удаление одну, точка две (на краях статуса соседа нет, и проверка не считается - поэтому у коротких отрезков всего 59). А пар n² / 2 всегда. Значит, выигрыш есть, пока k заметно меньше n².

Если пересекается почти всё со всем - k порядка n², - алгоритм делает n² · log n работы против n² у перебора и проигрывает. Его сила там, где пересечений мало: дорожки на плате, улицы на карте, границы участков - в реальных данных у отрезка обычно несколько соседей, а не все.

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

Точка пересечения в float. Два события с почти равными x встают в очередь в неверном порядке, обмен происходит не между соседями, и статус портится - дальше алгоритм пропускает пересечения или падает на поиске отрезка. Поэтому здесь числитель и знаменатель.

Без памяти о найденных парах. Пара, которая была соседями, разошлась и снова сошлась, попадёт в очередь дважды - и второй обмен вернёт их в неверный порядок.

Вырожденные случаи. Три отрезка через одну точку, общие концы, вертикальные отрезки - в учебнике де Берга и соавторов алгоритм обрабатывает их, упорядочивая статус «чуть ниже точки события» и переворачивая сразу весь блок отрезков, проходящих через точку. Учебная версия этого не делает.

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

Наложение карт. Слой районов и слой водоохранных зон надо совместить и получить участки «район × зона». Рёбра обоих слоёв - отрезки, их пересечения - новые вершины, и Бентли-Оттман находит их все. Чтобы из вершин и рёбер собрать грани, рёбра хранят в двусвязном списке рёбер: у каждого ребра ссылки на следующее ребро той же грани и на «близнеца» с другой стороны. Это стандартное представление разбиения плоскости в геобиблиотеках.

То же самое - булевы операции над многоугольниками (объединение, пересечение, разность) в графических редакторах и САПР: сначала все пересечения рёбер, потом обход граней.

Теперь сам

1. Статус снизу вверх a, b, c, d. Событие - пересечение b и c. Каким станет статус и какие пары проверить?
2. 1000 отрезков, 5000 пересечений. Во сколько раз примерно событий меньше, чем пар для перебора?
3. Почему точку пересечения добавляют в очередь, только если она правее прямой?