Пустая окружность и переворот рёбер: как получить Делоне из любой триангуляции

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

Как думаешь, как починить триангуляцию

В прошлом уроке веер и вставка точек дали корректную, но плохую триангуляцию: синус наименьшего угла 0,012 - это меньше градуса. Переделать её с нуля? Перебрать все триангуляции? Их экспоненциально много.

Ответь себе до чтения: что, если смотреть только на пары соседних треугольников - на четырёхугольники - и в каждом выбирать лучшую диагональ? Хватит ли таких местных исправлений, чтобы получить лучшую триангуляцию вообще? Не будет ли одно исправление портить другое бесконечно? Ответ урока: хватит, и не будет. На двадцати точках это 24 переворота.

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

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

Складка - нелегальное ребро, расправление - переворот, «ровнее» - лучше вектор углов.

Механизм: определитель, четырёхугольник и вектор углов

Точка внутри окружности. Для треугольника a, b, c против часовой стрелки и точки d:

| ax − dx   ay − dy   (ax − dx)² + (ay − dy)² |
| bx − dx   by − dy   (bx − dx)² + (by − dy)² |
| cx − dx   cy − dy   (cx − dx)² + (cy − dy)² |
определитель > 0   ⇔   d строго внутри окружности abc

Откуда здесь квадраты. Окружность - это множество точек, у которых x² + y² выражается линейно через x и y: x² + y² = αx + βy + γ. Добавив третью координату z = x² + y², точки поднимаются на параболоид, окружность превращается в плоское сечение, а «внутри окружности» - в «ниже плоскости». Определитель - это поворот в трёх измерениях, знак «выше или ниже плоскости через три поднятые точки». Подробнее об этом подъёме - следующий урок.

Нелегальное ребро и переворот. Ребро uv делят треугольники uvx и vuw. Если w внутри окружности uvx, ребро нелегально. Тогда четырёхугольник u, w, v, x выпуклый, и его можно разрезать другой диагональю xw - получатся треугольники uwx и wvx. После переворота новое ребро легально: это свойство вписанных углов - если w внутри окружности uvx, то v снаружи окружности uwx, а u снаружи окружности wvx.

Почему перевороты кончаются. Переворот нелегального ребра заменяет шесть углов двух треугольников так, что наименьший из них строго растёт. Значит, вектор углов всей триангуляции (все её углы по возрастанию, из прошлого урока) лексикографически растёт с каждым переворотом. Триангуляций конечное число, и ни одна не может повториться - процесс останавливается. Остановиться он может только там, где нелегальных рёбер нет. А триангуляция, у которой все рёбра легальны, - это и есть Делоне: пустота каждой окружности следует из легальности каждого ребра. Довод: пусть точка p внутри окружности треугольника t. Идём от t к p по треугольникам, которые пересекает отрезок. У легального ребра часть окружности за ним лежит внутри окружности соседа, поэтому p остаётся внутри окружности каждого следующего треугольника - вплоть до соседа самой p. А там p - противоположная вершина, и ребро нелегально. Отсюда же теорема: Делоне максимизирует наименьший угол - любую другую триангуляцию переворотами можно только улучшить до неё.

Почему O(n²). Конечность ещё не даёт границы. Её даёт другой довод: подними точки на параболоид (x, y, x² + y²). Переворот нелегального ребра только опускает поверхность из поднятых треугольников, поэтому выброшенное ребро не возвращается. Рёбер не больше n(n − 1)/2 - и переворотов не больше.

Единственность. Если никакие четыре точки не лежат на одной окружности, триангуляция Делоне одна. Если лежат - как у вершин квадрата, - годится любая диагональ, и выбирают по договору.

Правило

Триангуляция Делоне - все описанные окружности пусты. Ребро нелегально, если противоположная вершина соседнего треугольника внутри окружности; переворот заменяет его другой диагональю. Перевороты в любом порядке кончаются за O(n²) и дают Делоне - триангуляцию с наибольшим наименьшим углом. Предикат «в окружности» - определитель 3 × 3 в целых.

Разбор: ромб из прошлого урока

def v_kruge(a, b, c, d):
    # > 0: d строго внутри окружности через a, b, c (a, b, c против
    # часовой стрелки); определитель в целых - ответ точный
    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):
    # найти одно нелегальное ребро и перевернуть; True, если нашлось
    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
A = [0, 5]
B = [10, 0]
C = [20, 5]
D = [10, 10]
print(v_kruge(A, B, C, D), v_kruge(C, D, A, B))
tochki = [A, B, C, D]
treug = [[0, 1, 2], [0, 2, 3]]
print(delone_perevorotami(tochki, treug), treug)
print(v_kruge(A, B, D, C), v_kruge(B, C, D, A))

15000 15000
1 [[2, 3, 1], [3, 0, 1]]
-15000 -15000

Первая строка: D внутри окружности ABC, и, симметрично, B внутри окружности CDA - ребро AC нелегально с обеих сторон. Перевороты сделали ровно один шаг: треугольники стали C, D, B и D, A, B - ребро BD. Последняя строка: C вне окружности ABD, A вне окружности BCD. Делоне выбрала разрез с углами 53°, как и предсказывал прошлый урок.

Словарь rebro хранит направленное ребро «из a в b» одним числом a * 10000 + b - это верно, пока номеров меньше 10 000. У соседа то же ребро идёт обратно, «из v в u», отсюда ключ v * 10000 + u. Третья вершина соседа - сумма его вершин минус u и v.

Числа большие: в определителе квадраты разностей. При координатах до тысячи - до 10¹⁴. Целым Python это не страшно, но во float (точен до 2⁵³ ≈ 9 · 10¹⁵) запас уже не огромный.

Разбор: двадцать точек и 24 переворота

Триангуляция из прошлого урока - с тем же допущением: точку ровно на ребре строгие проверки пропустят, но здесь таких нет, треугольников 30 - и счётчик плохих пар «треугольник - точка внутри его окружности».

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
import math

def sin_min(tochki, treug):
    luchshij = 1.0
    for t in treug:
        for k in range(3):
            a = tochki[t[k]]
            b = tochki[t[(k + 1) % 3]]
            c = tochki[t[(k + 2) % 3]]
            ab = math.sqrt((b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2)
            ac = math.sqrt((c[0] - a[0]) ** 2 + (c[1] - a[1]) ** 2)
            luchshij = min(luchshij, abs(povorot(a, b, c)) / (ab * ac))
    return luchshij

def plohih(tochki, treug):
    # сколько пар «треугольник - точка внутри его окружности»
    schet = 0
    for a, b, c in treug:
        for i in range(len(tochki)):
            if i == a or i == b or i == c:
                continue
            if v_kruge(tochki[a], tochki[b], tochki[c], tochki[i]) > 0:
                schet = schet + 1
    return schet

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)
H, treug = triangulirovat(tochki)
print(plohih(tochki, treug), round(sin_min(tochki, treug), 3))
k = delone_perevorotami(tochki, treug)
print(plohih(tochki, treug), round(sin_min(tochki, treug), 3), k)

149 0.012
0 0.057 24

До переворотов 149 плохих пар и наименьший угол меньше градуса (синус 0,012). После 24 переворотов плохих пар ноль - окружности пусты, это Делоне. Синус наименьшего угла вырос до 0,057, около 3,3°.

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

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

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

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

Один проход вместо «пока есть нелегальные». Переворот делает нелегальными соседние рёбра, и один проход по всем рёбрам оставляет их нетронутыми.

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

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

Делоне - стандартная основа сеток: рельеф (TIN), конечные элементы, интерполяция по разбросанным замерам. Перевороты - не только учебный алгоритм: они же чинят триангуляцию после сдвига точек - не надо строить заново, достаточно перевернуть рёбра, ставшие нелегальными. scipy.spatial.Delaunay, CGAL, библиотека Triangle Шевчука считают тот же предикат - с точной или адаптивной арифметикой.

Теперь сам

1. Квадрат (0, 0), (4, 0), (4, 4), (0, 4), разрезанный диагональю. Легальна ли она? Сколько у квадрата триангуляций Делоне?
2. Треугольник (0, 0), (6, 0), (0, 6) и точка (3, 3). Внутри ли она окружности?
3. Почему переворот, а не просто «найти тупой угол и что-то с ним сделать»?