Триангуляции точек: сколько в них треугольников и почему важны углы

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

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

Геодезист измерил высоту в четырёх точках: два конца лощины на высоте 0 и два склона хребта на высоте 100. Карта рельефа строится так: точки соединяют в треугольники, а внутри треугольника высоту считают линейно. Четыре точки - это два треугольника, и разрезать четырёхугольник можно по любой из двух диагоналей.

Ответь себе до чтения: зависит ли высота в точке около центра от того, какую диагональ провести? И если да - то насколько? Ответ: при одной диагонали получается 20 метров, при другой - 100. Одни и те же измерения дают то ли лощину, то ли хребет.

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

Натянуть палатку на четыре колышка разной высоты. Ткань ляжет либо складкой по одной диагонали, либо по другой - и в середине палатки окажется то впадина, то гребень. Колышки одни и те же, а форма поверхности определяется тем, какие колышки соединил шов.

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

Механизм: формула Эйлера и вектор углов

Почему число треугольников одно. Триангуляция - плоский граф: n вершин, E рёбер и F граней (T треугольников плюс одна внешняя грань). Формула Эйлера n − E + F = 2. Каждый треугольник ограничен тремя рёбрами, внешняя грань - h рёбрами оболочки, и каждое ребро граничит ровно с двумя гранями: 3T + h = 2E. Из двух равенств T = 2n − 2 − h и E = 3n − 3 − h. Число не зависит от того, как провели рёбра.

Сколько всего триангуляций. Очень много. Уже у выпуклого многоугольника с n вершинами их число - число Каталана C(n − 2): у восьмиугольника 132, у двадцатиугольника - почти полмиллиарда. Перебирать их, чтобы выбрать лучшую, бессмысленно.

Чем сравнивать. У триангуляции T треугольников и 3T углов. Выпишем углы по возрастанию - это вектор углов. Из двух триангуляций лучше та, у которой этот вектор лексикографически больше: прежде всего у неё больше наименьший угол. Именно это и сделает триангуляция Делоне - у неё вектор углов наибольший из всех возможных (в общем положении, когда Делоне одна).

Правило

Триангуляция n точек, из которых h на оболочке, всегда содержит 2n − 2 − h треугольников и 3n − 3 − h рёбер. Триангуляции различаются формой треугольников, и сравнивают их вектором углов, начиная с наименьшего. Тонкие треугольники портят интерполяцию: значение тянется из далёких точек.

Разбор: двадцать точек, тридцать треугольников

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

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]
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)
rebra = {}
for a, b, c in treug:
    for u, v in [[a, b], [b, c], [c, a]]:
        rebra[min(u, v) * 1000 + max(u, v)] = True
n = len(tochki)
h = len(H)
print(n, h, len(treug), 2 * n - 2 - h, len(rebra), 3 * n - 3 - h)

20 8 30 30 49 49

Двадцать точек, восемь на оболочке. Треугольников 30, рёбер 49 - ровно 2 · 20 − 2 − 8 и 3 · 20 − 3 − 8. Счёт подтверждает и само построение: веер даёт h − 2 = 6 треугольников, а каждая из 12 внутренних точек добавляет два - превращает один треугольник в три. 6 + 24 = 30.

Эта триангуляция корректна, но плоха: внутренняя точка разрезает тот треугольник, куда попала, и узкие клинья из веера остаются навсегда. Переделать её в хорошую - работа следующего урока.

Разбор: лощина или хребет

Четыре точки из «Как думаешь», две триангуляции. Высота внутри треугольника - сумма высот вершин с весами, равными долям площадей (барицентрическая интерполяция). Рядом печатаем синус наименьшего угла триангуляции.

import math

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

def vysota(tri, q):
    # линейная интерполяция высоты внутри треугольника: доли площадей
    a, b, c = tri
    s = povorot(a, b, c)
    wa = povorot(b, c, q) / s
    wb = povorot(c, a, q) / s
    wc = povorot(a, b, q) / s
    return wa * a[2] + wb * b[2] + wc * c[2]

def sin_min_ugla(tri):
    # синус наименьшего угла: наименьший угол меньше 60 градусов,
    # и меньший синус у меньшего из углов
    luchshij = 1.0
    for k in range(3):
        a = tri[k]
        b = tri[(k + 1) % 3]
        c = tri[(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

# два края долины (высота 0) и два склона хребта (высота 100)
A = [0, 5, 0]
B = [10, 0, 100]
C = [20, 5, 0]
D = [10, 10, 100]
q = [10, 4]
print("диагональ AC:", vysota([A, B, C], q),
      round(min(sin_min_ugla([A, B, C]), sin_min_ugla([A, C, D])), 3))
print("диагональ BD:", vysota([A, B, D], q),
      round(min(sin_min_ugla([A, B, D]), sin_min_ugla([B, C, D])), 3))

диагональ AC: 20.0 0.447
диагональ BD: 100.0 0.8

Точка q с диагональю AC оказывается в треугольнике ABC около долины: высота 20. С диагональю BD она лежит прямо на шве хребта: 100. Синус наименьшего угла - 0,447 (около 27°) против 0,8 (около 53°). Какая высота верна, по четырём точкам знать нельзя. Но толстые треугольники берут высоту у ближайших замеров, а клинья ABC и ACD растянуты на всю длину ромба и тянут ноль из точек A и C, до которых около 10 единиц, - а до B и D всего 4 и 6. Поэтому предпочитают разрез с большим наименьшим углом.

Синус здесь вместо самого угла - в движке нет арксинуса. Он годится для сравнения: наименьший угол треугольника меньше 60°, а на этом отрезке синус растёт вместе с углом.

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

Сетка для расчёта прочности с тонкими треугольниками. Метод конечных элементов на вытянутых треугольниках даёт плохо обусловленные системы: ошибки округления растут, решение «дрожит». Генераторы сеток прежде всего следят за наименьшим углом.

Рельеф по случайной триангуляции. Изолинии высот зигзагом, ложные хребты и ямы - как во втором разборе, где один шов превратил лощину в хребет.

Считать, что лучшую триангуляцию надо искать перебором. Их экспоненциально много, а нужная строится за n · log n.

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

Цифровые модели рельефа (TIN в геоинформационных системах) - триангуляция точек съёмки с интерполяцией по треугольникам. Сетки для расчётов теплопроводности и прочности. Трёхмерные модели в графике - тоже треугольники, и качество их формы видно на освещении. Во всех этих случаях начинают с триангуляции Делоне и потом её улучшают.

Теперь сам

1. 100 точек, из них 10 на выпуклой оболочке. Сколько треугольников и рёбер в любой их триангуляции?
2. Все n точек на выпуклой оболочке. Сколько треугольников?
3. Почему при сравнении триангуляций смотрят прежде всего на наименьший угол, а не на средний?