Делоне и Вороной: двойственность, параболоид и что лежит внутри Делоне

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

Как думаешь, какие рёбра нужны дереву

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

Ответь себе до чтения: может быть, многие пары можно выбросить заранее, не рискуя потерять ребро дерева? Какие? Ответ урока: достаточно 49 рёбер Делоне - дерево получится то же. В общем случае 3n − 3 − h (h - точек на оболочке) вместо n²/2.

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

Карта стран и карта столиц. Карта Вороного - это страны: у каждой столицы своя территория. Делоне - это сеть дорог между столицами соседних стран: если у двух стран есть общая граница, между их столицами проложена дорога. По одной карте другую восстановить можно сразу: границы знают, кто сосед, а дороги знают, где проходит граница - посередине и поперёк.

Механизм: двойственность и подъём на параболоид

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

Параболоид. Поднимем каждую точку (x, y) в (x, y, x² + y²). Окружность на плоскости x² + y² = αx + βy + γ поднимается в сечение параболоида плоскостью z = αx + βy + γ, и точка внутри окружности поднимается под эту плоскость. Значит, «окружность через p, q, r пуста» означает: плоскость через поднятые p, q, r оставляет все остальные поднятые точки над собой - это грань нижней выпуклой оболочки. Триангуляция Делоне - проекция нижней оболочки поднятых точек. Отсюда два следствия: определитель «в окружности» - это поворот в трёх измерениях, а Делоне в d измерениях - это оболочка в d + 1. Так её и строит scipy.spatial.Delaunay: через оболочку в библиотеке Qhull.

Что лежит в Делоне. Ближайший сосед q точки p: окружность с диаметром pq пуста - иначе внутри нашлась бы точка ближе q. Значит, pq - ребро Делоне. Ребро минимального остовного дерева pq: если бы в круге с диаметром pq лежала точка r, то pr и rq обе короче pq, и замена pq на одно из них дала бы дерево легче. Значит, круг пуст, и pq - ребро Делоне. При попарно разных расстояниях цепочка вложений: граф ближайших соседей ⊆ остовное дерево ⊆ граф Габриеля (рёбра с пустым кругом на диаметре) ⊆ Делоне.

Правило

Делоне двойственна Вороному: рёбра Делоне соединяют соседние ячейки, треугольники соответствуют вершинам. Делоне - проекция нижней оболочки точек, поднятых на z = x² + y², а «в окружности» - поворот в трёх измерениях. Делоне содержит граф ближайших соседей и минимальное остовное дерево, поэтому дерево строится Крускалом по 3n − 3 − h её рёбрам за O(n · log n).

Разбор: определитель окружности - это поворот в пространстве

Проверим тождество на двухстах случайных четвёрках точек: определитель «в окружности» и трёхмерный поворот поднятых точек.

import random

def v_kruge(a, b, c, d):
    m = []
    for p in [a, b, c]:
        dx = p[0] - d[0]
        dy = p[1] - d[1]
        m.append([dx, dy, dx * dx + dy * dy])
    return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]))

def podnyat(p):
    # точка плоскости - на параболоид z = x² + y²
    return [p[0], p[1], p[0] * p[0] + p[1] * p[1]]

def povorot3(a, b, c, d):
    # знак: по какую сторону от плоскости abc лежит d
    m = []
    for p in [a, b, c]:
        m.append([p[0] - d[0], p[1] - d[1], p[2] - d[2]])
    return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]))

random.seed(9)
sovpalo = 0
vnutri = 0
for t in range(200):
    ch = []
    for k in range(4):
        ch.append([random.randint(-50, 50), random.randint(-50, 50)])
    a, b, c, d = ch
    x = v_kruge(a, b, c, d)
    y = povorot3(podnyat(a), podnyat(b), podnyat(c), podnyat(d))
    if x == y:
        sovpalo = sovpalo + 1
    if x > 0:
        vnutri = vnutri + 1
print(sovpalo, vnutri)
print(v_kruge([0, 5], [10, 0], [20, 5], [10, 10]))
print(povorot3(podnyat([0, 5]), podnyat([10, 0]), podnyat([20, 5]),
               podnyat([10, 10])))

200 90
15000
15000

Двести из двухсот - значения совпали до единицы. Это не совпадение знаков, а равенство чисел: в столбце определителя «в окружности» стоит (ax − dx)² + (ay − dy)², а у поднятых точек - (ax² + ay²) − (dx² + dy²). Разница между ними - −2dx(ax − dx) − 2dy(ay − dy), то есть сумма первых двух столбцов с коэффициентами, и определитель от этого не меняется. Девяносто положительных значений - те четвёрки, где d внутри окружности при обходе a, b, c против часовой стрелки или снаружи при обходе по часовой. Последние две строки - ромб из прошлого урока: те же 15 000.

Разбор: остовное дерево по рёбрам Делоне

Двадцать точек из прошлых уроков, Делоне переворотами, затем Крускал дважды: по рёбрам Делоне и по всем парам.

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 obolochka(tochki):
    t = []
    for i in range(len(tochki)):
        t.append([tochki[i][0], tochki[i][1], i])
    t.sort()
    niz = []
    for p in t:
        while len(niz) >= 2 and povorot(niz[-2], niz[-1], p) <= 0:
            niz.pop()
        niz.append(p)
    verh = []
    for p in t[::-1]:
        while len(verh) >= 2 and povorot(verh[-2], verh[-1], p) <= 0:
            verh.pop()
        verh.append(p)
    otvet = []
    for p in niz[:-1] + verh[:-1]:
        otvet.append(p[2])
    return otvet

def triangulirovat(tochki):
    H = obolochka(tochki)
    treug = []
    for k in range(1, len(H) - 1):
        treug.append([H[0], H[k], H[k + 1]])
    for i in range(len(tochki)):
        if i in H:
            continue
        p = tochki[i]
        for k in range(len(treug)):
            a, b, c = treug[k]
            if povorot(tochki[a], tochki[b], p) > 0:
                if povorot(tochki[b], tochki[c], p) > 0:
                    if povorot(tochki[c], tochki[a], p) > 0:
                        treug[k] = [a, b, i]
                        treug.append([b, c, i])
                        treug.append([c, a, i])
                        break
    return [H, treug]
def v_kruge(a, b, c, d):
    m = []
    for p in [a, b, c]:
        dx = p[0] - d[0]
        dy = p[1] - d[1]
        m.append([dx, dy, dx * dx + dy * dy])
    return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]))

def perevernut_odno(tochki, treug):
    rebro = {}
    for k in range(len(treug)):
        a, b, c = treug[k]
        rebro[a * 10000 + b] = k
        rebro[b * 10000 + c] = k
        rebro[c * 10000 + a] = k
    for k in range(len(treug)):
        for s in range(3):
            u = treug[k][s]
            v = treug[k][(s + 1) % 3]
            x = treug[k][(s + 2) % 3]
            klyuch = v * 10000 + u
            if klyuch not in rebro:
                continue
            m = rebro[klyuch]
            w = treug[m][0] + treug[m][1] + treug[m][2] - u - v
            if v_kruge(tochki[u], tochki[v], tochki[x], tochki[w]) > 0:
                treug[k] = [u, w, x]
                treug[m] = [w, v, x]
                return True
    return False

def delone_perevorotami(tochki, treug):
    perevorotov = 0
    while perevernut_odno(tochki, treug):
        perevorotov = perevorotov + 1
    return perevorotov
def kv(p, q):
    return (p[0] - q[0]) ** 2 + (p[1] - q[1]) ** 2

def kraskal(n, rebra, tochki):
    # дерево Крускала: рёбра по длине, лес предков против циклов;
    # сравниваем квадраты длин - дерево то же, вес - сумма квадратов
    spisok = []
    for u, v in rebra:
        spisok.append([kv(tochki[u], tochki[v]), u, v])
    spisok.sort()
    roditel = list(range(n))
    ves = 0
    vybrano = []
    for d, u0, v0 in spisok:
        u = u0
        v = v0
        while roditel[u] != u:
            u = roditel[u]
        while roditel[v] != v:
            v = roditel[v]
        if u != v:
            roditel[u] = v
            ves = ves + d
            vybrano.append([u0, v0])
    return [ves, vybrano]

random.seed(4)
tochki = []
while len(tochki) < 20:
    p = [random.randint(0, 100), random.randint(0, 100)]
    if p not in tochki:
        tochki.append(p)
n = len(tochki)
H, treug = triangulirovat(tochki)
delone_perevorotami(tochki, treug)
rebra_d = {}
for a, b, c in treug:
    for u, v in [[a, b], [b, c], [c, a]]:
        rebra_d[min(u, v) * 1000 + max(u, v)] = [min(u, v), max(u, v)]
vse = []
for u in range(n):
    for v in range(u + 1, n):
        vse.append([u, v])
print(len(rebra_d), len(vse))
po_delone = kraskal(n, list(rebra_d.values()), tochki)
polnoe = kraskal(n, vse, tochki)
print(po_delone[0], len(po_delone[1]), po_delone == polnoe)
# ближайший сосед каждой точки - ребро Делоне?
v_delone = 0
for u in range(n):
    luchshij = -1
    for v in range(n):
        if v != u:
            d = kv(tochki[u], tochki[v])
            if luchshij == -1 or d < d_min:
                luchshij = v
                d_min = d
    if min(u, luchshij) * 1000 + max(u, luchshij) in rebra_d:
        v_delone = v_delone + 1
print(v_delone, n)

49 190
5704 19 True
20 20

Рёбер Делоне 49, пар 190. Деревья совпали ребро в ребро (True - списки выбранных рёбер равны): 19 рёбер, сумма квадратов длин 5704. Крускал сравнивает квадраты, дерево от этого не меняется, потому что квадрат возрастает вместе с длиной. И у каждой из двадцати точек ребро до ближайшего соседа оказалось ребром Делоне.

Сколько это экономит. Для миллиона точек пар 5 · 10¹¹, рёбер Делоне около 3 · 10⁶. Делоне строится за n · log n, Крускал по её рёбрам - тоже n · log n. Граница приёма: он только для евклидовых расстояний на плоскости (или в малой размерности). Если длина ребра - цена по дорогам, вложение в Делоне не выполняется.

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

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

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

Перенос на дорожную сеть. На графе дорог Делоне не поможет: расстояние другое, и вложение держится только на геометрии кругов.

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

Минимальное дерево и кластеризация: разрезать самые длинные рёбра остовного дерева - это одиночная связь в иерархической кластеризации, и Делоне делает её быстрой. Восстановление поверхности по облаку точек и объёмные сетки. И подъём на параболоид - общий приём: задачу про окружности на плоскости он превращает в задачу про плоскости в пространстве, где работают уже известные алгоритмы оболочки.

Теперь сам

1. Сколько вершин у диаграммы Вороного 100 точек, если их Делоне имеет 190 треугольников?
2. Почему ребро между точкой и её ближайшим соседом - всегда ребро Делоне?
3. Поднимем точки (0, 0), (2, 0), (0, 2) и (1, 1) на параболоид. Где окажется (1, 1) относительно плоскости через три первые - выше или ниже - и что это значит?