Что даёт оболочка: диаметр калиперами и точка внутри за логарифм

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

Как думаешь, сколько пар проверять

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

Ответь себе до чтения: сколько пар на самом деле надо проверить, если сначала построить оболочку? Её вершин окажется 45. Пар вершин - 990. А вращающиеся калиперы обойдутся двумястами двадцатью поворотами - около пяти на вершину.

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

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

Калиперы - буквально «штангенциркуль» по-английски. Алгоритм вращает две параллельные прямые вокруг многоугольника, перескакивая от ребра к ребру.

Механизм: опорные прямые и монотонность

Почему диаметр на оболочке. Возьми самую далёкую пару p, q. Если бы p лежала внутри оболочки, её можно было бы сдвинуть от q, оставаясь в оболочке, и расстояние выросло бы. А в выпуклом многоугольнике самая далёкая от q точка - вершина: расстояние до q растёт к краю и достигает максимума в углу. Значит, обе точки диаметра - вершины.

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

Почему пар мало. Обходим рёбра по порядку против часовой стрелки. Самая далёкая от ребра вершина при этом тоже движется только вперёд: на выпуклом многоугольнике высота над ребром растёт до максимума и потом падает. Поэтому указатель j на дальнюю вершину никогда не возвращается назад, и за полный обход он проходит круг с небольшим - меньше 2h шагов. Каждая проверка - два поворота, проверок - шаги указателя плюс по одной останавливающей на ребро, итого меньше 6h поворотов - линейно. Почему больший поворот значит большую высоту: |поворот(a, b, c)| - удвоенная площадь треугольника, то есть |ab|, умноженная на расстояние от c до прямой; при одном ребре он растёт вместе с расстоянием. Это тот же приём двух указателей, что в задачах на отсортированных массивах.

Точка внутри за логарифм. Проведём из вершины H[0] лучи ко всем остальным вершинам - многоугольник разрежется на веер треугольников. Лучи упорядочены по углу, и поворот относительно луча H[0] → H[k] меняет знак ровно один раз при движении по k. Значит, нужный треугольник ищется двоичным поиском, а потом одна проверка: по какую сторону от ребра H[k] → H[k+1] лежит точка.

Правило

Диаметр множества точек - пара вершин его оболочки; вращающиеся калиперы находят его за O(h) после построения оболочки: для каждого ребра по порядку двигают вперёд указатель на самую далёкую вершину. Точка внутри выпуклого многоугольника проверяется двоичным поиском по вееру из одной вершины за O(log h). Обе оценки держатся на выпуклости - для невыпуклого многоугольника они неверны.

Разбор: двести датчиков

Датчики стоят в кольце радиусов от 950 до 1000. Сравним калиперы с перебором пар вершин оболочки и пар всех точек.

import math
import random
schet = [0]

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

def dal(a, b):
    return (b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2

def endryu(tochki):
    t = sorted(tochki)
    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)
    return niz[:-1] + verh[:-1]

def diametr(H):
    h = len(H)
    luchshee = 0
    j = 1
    for i in range(h):
        a = H[i]
        b = H[(i + 1) % h]
        # двигаем j, пока точка уходит дальше от прямой ab
        while povorot(a, b, H[(j + 1) % h]) > povorot(a, b, H[j]):
            j = (j + 1) % h
        luchshee = max(luchshee, dal(a, H[j]), dal(b, H[j]))
    return luchshee

random.seed(3)
tochki = []
for i in range(200):
    ugol = random.uniform(0, 2 * math.pi)
    r = random.uniform(950, 1000)
    tochki.append([int(r * math.cos(ugol)), int(r * math.sin(ugol))])
H = endryu(tochki)
schet[0] = 0
d = diametr(H)
print(len(tochki), len(H), d, schet[0])
perebor = 0
for i in range(len(H)):
    for j in range(i + 1, len(H)):
        perebor = max(perebor, dal(H[i], H[j]))
vse = 0
for i in range(len(tochki)):
    for j in range(i + 1, len(tochki)):
        vse = max(vse, dal(tochki[i], tochki[j]))
print(perebor, vse, len(H) * (len(H) - 1) // 2)

200 45 3983994 220
3983994 3983994 990

Оболочка из 45 вершин, квадрат диаметра 3 983 994 (около 1996 - почти два радиуса). Калиперы потратили 220 поворотов. Перебор пар вершин оболочки дал то же число через 990 сравнений, перебор всех точек - через 19 900. Расстояния везде квадратами: корень не нужен, чтобы найти наибольшее, и квадраты целых точны.

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

Разбор: внутри ли точка, за восемь поворотов

На той же оболочке проверим 304 точки двумя способами: перебором рёбер («слева от каждого ребра») и двоичным поиском по вееру.

import math
import random
schet = [0]

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

def endryu(tochki):
    t = sorted(tochki)
    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)
    return niz[:-1] + verh[:-1]

def vnutri_perebor(H, p):
    for i in range(len(H)):
        if povorot(H[i], H[(i + 1) % len(H)], p) < 0:
            return False
    return True

def vnutri_poisk(H, p):
    # веер треугольников из H[0]; граница считается внутренней
    h = len(H)
    if povorot(H[0], H[1], p) < 0 or povorot(H[0], H[h - 1], p) > 0:
        return False
    lo = 1
    hi = h - 1
    while hi - lo > 1:
        mid = (lo + hi) // 2
        if povorot(H[0], H[mid], p) >= 0:
            lo = mid
        else:
            hi = mid
    return povorot(H[lo], H[lo + 1], p) >= 0

random.seed(3)
tochki = []
for i in range(200):
    ugol = random.uniform(0, 2 * math.pi)
    r = random.uniform(950, 1000)
    tochki.append([int(r * math.cos(ugol)), int(r * math.sin(ugol))])
H = endryu(tochki)
zaprosy = [[0, 0], [2000, 0], H[5], [975, 0]]
for i in range(300):
    x = random.randint(-1000, 1000)
    zaprosy.append([x, random.randint(-1000, 1000)])
soglasny = 0
vnutri = 0
schet_perebor = 0
schet_poisk = 0
for p in zaprosy:
    schet[0] = 0
    a = vnutri_perebor(H, p)
    schet_perebor = schet_perebor + schet[0]
    schet[0] = 0
    b = vnutri_poisk(H, p)
    schet_poisk = schet_poisk + schet[0]
    if a == b:
        soglasny = soglasny + 1
    if b:
        vnutri = vnutri + 1
print(len(H), len(zaprosy), soglasny, vnutri)
print(schet_perebor, schet_poisk)
print(vnutri_poisk(H, [0, 0]), vnutri_poisk(H, [2000, 0]))
print(vnutri_poisk(H, H[5]))

45 304 304 231
11801 2486
True False
True

Ответы совпали на всех 304 точках, внутри оказались 231. Перебор потратил 11 801 поворот - почти 39 на запрос (внешние точки отсекаются раньше, внутренние проверяются всеми 45 рёбрами). Двоичный поиск - 2486, около восьми: log₂ 45 ≈ 5,5 шагов поиска плюс две проверки на краях веера и одна на ребре. Вершина H[5] засчитана внутренней - так договорились в коде: граница внутри.

Граница применимости: оба приёма требуют выпуклости. Вырожденная оболочка-отрезок (h = 2) калиперам не мешает - все повороты нули, и диаметр - сам отрезок, - а веер из двух вершин треугольников не образует, и такую оболочку проверяют отдельно. В букве Г лучи из одной вершины выходят за фигуру, и двоичный поиск по вееру отвечает чепуху - там работает только луч и чётность из прошлой главы.

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

Калиперы на невыпуклом контуре. Указатель на дальнюю вершину перестаёт двигаться монотонно, и алгоритм пропускает настоящий диаметр. Сначала оболочка, потом калиперы.

Диаметр через корни. sqrt для сравнения расстояний не нужен, а в float он добавляет шум, и две почти равные пары сравниваются неверно. Сравнивай квадраты - на целых они точны.

Указатель j с нуля для каждого ребра. Это снова h²: весь выигрыш калиперов - в том, что j продолжает с того места, где остановился.

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

Самая далёкая пара - размер кластера на карте, худший случай длины кабеля, разброс отпечатков. Калиперы решают и родственные задачи тем же обходом: ширину фигуры (наименьшее расстояние между параллельными опорными прямыми - пролезет ли деталь в щель) и наименьший охватывающий прямоугольник любой ориентации - одна его сторона обязательно лежит на ребре оболочки, это нужно для упаковки и для рамок объектов в компьютерном зрении. Точка в выпуклом многоугольнике за логарифм - проверка «внутри ли поле зрения», «внутри ли выпуклой зоны», когда запросов миллионы, а многоугольник один.

Теперь сам

1. Оболочка - квадрат (0, 0), (4, 0), (4, 4), (0, 4). Какие антиподальные пары найдут калиперы и каков квадрат диаметра?
2. Сколько поворотов примерно потратит двоичный поиск по вееру на оболочке из 1000 вершин?
3. Почему самая далёкая от точки q точка выпуклого многоугольника - вершина, а не середина ребра?