Вычислительная геометрия: примитивы, выпуклая оболочка, sweep line

Вычислительная геометрия — это единственный раздел алгоритмов, где вы можете написать формально доказанный алгоритм, прогнать его на реальных данных и получить бесконечный цикл. Не из-за ошибки в логике: логика будет верна. Из-за того, что double посчитал знак векторного произведения неправильно, и алгоритм узнал, что точка A левее точки B, а точка B левее точки A.

Это задаёт тон всей теме. Здесь два уровня понимания. Первый — комбинаторный: как устроены оболочки, развёртки, диаграммы Вороного. Второй — арифметический: почему эти конструкции разваливаются на вещественных числах и что с этим делают в CGAL, PostGIS и Google S2. Инженер, который знает только первый уровень, пишет код, падающий на продакшн-данных. Мы разберём оба.

Зачем это нужно на практике

Геометрия не «олимпиадная экзотика». Она живёт под капотом почти всего, что связано с пространством и физическим миром:

• Карты и GIS. «Какие магазины в этом полигоне доставки», «упростить контур области до 500 точек для зума 8», «пересекаются ли эти два района» — это point-in-polygon, Douglas–Peucker и булевы операции над полигонами.
• Игры и физика. Обнаружение столкновений: broad phase (пространственные индексы, sweep and prune) и narrow phase (GJK, SAT), где выпуклые оболочки — базовое представление тел.
• CAD/CAM и 3D-печать. Триангуляции, булевы операции над телами, offsetting контуров.
• Компьютерное зрение и ML. Non-maximum suppression — это IoU прямоугольников; выпуклая оболочка кластера, alpha-shapes, k-d деревья для kNN.
• Робототехника. Планирование пути в конфигурационном пространстве, диаграммы Вороного как «максимально безопасные» маршруты.
• Инфраструктура рекламы и вёрстки. Площадь объединения прямоугольников (viewability), overlap-детект в layout-движках.

Общее у всех этих задач одно: наивное решение квадратично или хуже, а данных миллионы. Вычислительная геометрия — это набор приёмов, как получить O(n log n).

Всё держится на одном предикате

Начнём с фундамента. Точка — пара чисел. Вектор — тоже пара. Два произведения:

• Скалярное dot(u, v) = ux·vx + uy·vy — про «сонаправленность», проекции, углы.
• Псевдоскалярное (векторное) произведение cross(u, v) = ux·vy − uy·vx — про ориентированную площадь параллелограмма и, главное, про сторону.

Практически всегда нам нужна трёхточечная форма:

cross(A, B, C) = (Bx − Ax)·(Cy − Ay) − (By − Ay)·(Cx − Ax)

Её знак отвечает на вопрос: «поворачиваем ли мы налево, идя по пути A → B → C?»

Это и есть предикат ориентации orient2d. Он один — и на нём стоит вся плоская геометрия:

Базовые примитивы на Python:

from math import hypot

Point = tuple[int, int]          # целые координаты — сознательный выбор, см. раздел о робастности

def sub(a: Point, b: Point) -> Point:
    return (a[0] - b[0], a[1] - b[1])

def dot(u: Point, v: Point) -> int:
    return u[0] * v[0] + u[1] * v[1]

def cross(o: Point, a: Point, b: Point) -> int:
    """Удвоенная ориентированная площадь треугольника OAB.
    > 0 — поворот налево (CCW), < 0 — направо (CW), == 0 — коллинеарны."""
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def sgn(x: int) -> int:
    return (x > 0) - (x < 0)

def dist2(a: Point, b: Point) -> int:
    """Квадрат расстояния. Держите расстояния в квадратах как можно дольше:
    так вы остаётесь в целых числах и не теряете точность на sqrt."""
    return (a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2

Правило номер один в геометрическом коде: не считать углы через atan2 и не считать расстояния через sqrt, если задачу решает знак cross или сравнение dist2. Тригонометрия вносит ошибку и заметно медленнее; сравнение целых — точное и мгновенное.

Площадь многоугольника: формула шнурования

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

def signed_area2(poly: list[Point]) -> int:
    """Удвоенная ориентированная площадь. > 0 — вершины идут против часовой стрелки.
    Работает и для невыпуклых простых многоугольников. O(n) времени, O(1) памяти."""
    s = 0
    n = len(poly)
    for i in range(n):
        x1, y1 = poly[i]
        x2, y2 = poly[(i + 1) % n]
        s += x1 * y2 - x2 * y1       # cross относительно начала координат
    return s

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

Полезное следствие — теорема Пика для многоугольника с целыми вершинами: S = I + B/2 − 1, где I — внутренние целые точки, B — целые точки на границе (B считается как сумма gcd(|dx|, |dy|) по рёбрам). Это мост к теории чисел, про которую — в https://courses.digitable.life/post/algorithms/13-number-theory/.

Пересечение отрезков

Классика из CLRS. Отрезки p1p2 и p3p4 пересекаются, если концы каждого лежат по разные стороны от другого — плюс аккуратная обработка коллинеарных случаев:

def on_segment(p: Point, q: Point, r: Point) -> bool:
    """Предполагая cross(p, q, r) == 0, лежит ли r внутри bounding box отрезка pq."""
    return (min(p[0], q[0]) <= r[0] <= max(p[0], q[0]) and
            min(p[1], q[1]) <= r[1] <= max(p[1], q[1]))

def segments_intersect(p1: Point, p2: Point, p3: Point, p4: Point) -> bool:
    d1 = sgn(cross(p3, p4, p1))
    d2 = sgn(cross(p3, p4, p2))
    d3 = sgn(cross(p1, p2, p3))
    d4 = sgn(cross(p1, p2, p4))

    # Собственное пересечение: концы строго по разные стороны в обе стороны.
    if d1 * d2 < 0 and d3 * d4 < 0:
        return True

    # Вырожденные случаи: конец одного отрезка лежит на другом.
    if d1 == 0 and on_segment(p3, p4, p1): return True
    if d2 == 0 and on_segment(p3, p4, p2): return True
    if d3 == 0 and on_segment(p1, p2, p3): return True
    if d4 == 0 and on_segment(p1, p2, p4): return True
    return False

Сложность O(1), и — важно — на целых координатах результат точен. Обратите внимание, сколько строк съедают вырожденные случаи: это типично. В геометрии основной объём кода и основной объём багов — не в «главной» ветке, а в касаниях, коллинеарностях и совпадающих точках.

Точка внутри многоугольника

Два подхода. Ray casting: пускаем луч из точки вправо и считаем пересечения с рёбрами; нечётное число — внутри. Winding number: суммируем ориентированные углы; ненулевая сумма — внутри. Для простых многоугольников они совпадают, для самопересекающихся дают разные (обе осмысленные) семантики — это ровно то же различие, что evenodd и nonzero в SVG и в правилах заливки графических движков.

def point_in_polygon(poly: list[Point], p: Point) -> bool:
    """Ray casting без деления: только целочисленные cross-предикаты.
    Точки строго на границе этот вариант не выделяет — проверяйте их отдельно.
    O(n) времени, O(1) памяти."""
    x, y = p
    inside = False
    n = len(poly)
    for i in range(n):
        a = poly[i]
        b = poly[(i + 1) % n]
        # Ребро пересекает горизонталь y (полуоткрытый интервал — так каждая
        # вершина учитывается ровно один раз и «луч через вершину» не ломает счёт).
        if (a[1] > y) != (b[1] > y):
            c = cross(a, b, p)
            if (c > 0) == (b[1] > a[1]):   # точка слева от восходящего ребра
                inside = not inside
    return inside

Ключевая деталь — полуоткрытый интервал (a.y > y) != (b.y > y). Именно он решает вечную проблему «луч попал точно в вершину»: вершина считается принадлежащей только одному из двух смежных рёбер. Наивная проверка min(a.y,b.y) <= y <= max(a.y,b.y) даёт двойной счёт и случайные ошибки на реальных данных.

Для выпуклого многоугольника есть более быстрый способ: O(log n) бинарным поиском по «вееру» треугольников от нулевой вершины. Приём тот же, что в https://courses.digitable.life/post/algorithms/03-searching-and-binary-search/: находим сектор, куда попадает точка, затем один предикат ориентации. Если запросов много, а полигон один — это разница между 10 мс и 10 мкс на запрос.

Робастность: почему код падает на реальных данных

Все алгоритмы ниже доказаны в предположении точной вещественной арифметики. Компьютер её не даёт. Последствия не «немного неточный ответ», а нарушение структурных инвариантов: выпуклая оболочка получается невыпуклой, sweep line теряет порядок в статусе, алгоритм зацикливается или падает по индексу.

Классический пример: orient2d для трёх почти коллинеарных точек с координатами порядка 1e7. Разность произведений даёт катастрофическую потерю значимости — знак получается случайным. А знак — это ответ на вопрос «слева или справа», на котором построена вся логика.

Практические правила:

1. Целые координаты — лучший выбор, если он возможен. Работайте в фиксированной точке: умножьте вход на 10^6 и округлите. Тогда cross точен, если промежуточные произведения влезают в int64 (координаты до ~2·10^9 в разности — безопасно). В Python целые произвольной точности, в C++/Go/Rust следите за переполнением: __int128 или проверка диапазона.
2. Один EPS на весь проект — это ошибка. Абсолютный 1e-9 бессмыслен для координат порядка 1e9 и слишком груб для нормализованных [0,1]. Если уж пользуетесь эпсилоном, делайте его относительным к масштабу входа.
3. Не сравнивайте точки на равенство по float. Дедупликация «почти совпадающих» точек должна быть явным шагом препроцессинга (снап к сетке), а не побочным эффектом сравнений.
4. Адаптивные предикаты. Джонатан Шевчук показал, как считать orient2d/incircle точно и при этом быстро: сначала дешёвая оценка с оценкой ошибки, и только если знак не определён — дорогое точное вычисление. На типичных данных дорогая ветка почти не срабатывает. Это стандарт де-факто: Robust Predicates.
5. Готовые ядра. CGAL даёт выбор ядра одной строкой: Exact_predicates_inexact_constructions_kernel (точные предикаты, приближённые конструкции — обычно то, что нужно) или ..._exact_constructions_kernel. GEOS/JTS используют похожие приёмы под PostGIS и Shapely.

Отдельная ловушка — конструкции против предикатов. Предикат возвращает знак и может быть сделан точным дёшево. Конструкция (точка пересечения двух отрезков) порождает новое число, которое почти всегда иррационально в исходном представлении, и дальнейшие предикаты уже с ним. В Bentley–Ottmann это главный источник проблем: точки пересечения идут обратно в очередь событий и участвуют в сравнениях.

Выпуклая оболочка

Определение. Выпуклая оболочка множества точек — минимальное выпуклое множество, их содержащее. Физическая аналогия: вбейте гвозди в точках и натяните резинку.

Зачем она нужна: оболочка — это компактное «внешнее описание» облака точек размера h ≤ n, и множество задач сводится к ней. Диаметр множества, ширина, минимальный охватывающий прямоугольник, проверка разделимости двух множеств, GJK-коллизии, первый шаг во многих задачах оптимизации.

Монотонная цепь Эндрю

Лучший практический алгоритм: сортируем по (x, y), строим отдельно нижнюю и верхнюю цепи стеком.

def convex_hull(points: list[Point]) -> list[Point]:
    """Монотонная цепь Эндрю. Возвращает вершины оболочки против часовой стрелки,
    без повторов и без точек на рёбрах.
    Время O(n log n) — доминирует сортировка. Память O(n)."""
    pts = sorted(set(points))
    if len(pts) <= 2:
        return pts

    def build(seq) -> list[Point]:
        chain: list[Point] = []
        for p in seq:
            # Пока последний поворот не строго левый — вершина не может быть на оболочке.
            while len(chain) >= 2 and cross(chain[-2], chain[-1], p) <= 0:
                chain.pop()
            chain.append(p)
        return chain

    lower = build(pts)
    upper = build(reversed(pts))
    return lower[:-1] + upper[:-1]   # крайние точки не дублируем

Почему O(n) после сортировки. Каждая точка добавляется в стек ровно один раз и выталкивается не более одного раза. Суммарное число операций pop не превосходит числа push, то есть n. Это учётный (амортизационный) аргумент — тот же приём, что для скользящего окна в https://courses.digitable.life/post/algorithms/04-two-pointers-sliding-window/ и для динамического массива; подробнее об этом стиле доказательств — https://courses.digitable.life/post/algorithms/01-analysis-and-proofs/.

Инвариант. После обработки префикса точек chain содержит нижнюю выпуклую оболочку этого префикса, и все повороты в ней строго левые.

Развилка <= 0 против < 0. Это не косметика:

• cross <= 0 выталкивает коллинеарные точки → в результате только «настоящие» углы. Обычно это то, что нужно.
• cross < 0 сохраняет точки на рёбрах оболочки. Нужно, например, если вы считаете количество точек на границе или строите оболочку для последующего сопоставления вершин. Осторожно: при таком варианте вырожденный вход (все точки на одной прямой) даёт «оболочку» с дублирующимися вершинами при обходе туда и обратно — обрабатывайте отдельно.

Обход Грэхема и семейство алгоритмов

Обход Грэхема (1972) — исторически первый O(n log n): выбрать самую нижнюю точку, отсортировать остальные по полярному углу вокруг неё, пройти стеком. Он делает то же, что монотонная цепь, но сортировка по углу требует либо atan2 (медленно и неточно), либо сравнения через cross (точно, но компаратор не является полным порядком при коллинеарных точках — источник тонких багов). Практический вывод: пишите монотонную цепь.

Если точки уже отсортированы (например, приходят потоком по времени), монотонная цепь работает за честные O(n) — так строят оболочку в оптимизации DP (техника «convex hull trick», см. https://courses.digitable.life/post/algorithms/07-dynamic-programming/).

В 3D картина меняется: оболочка имеет O(n) граней, а инкрементальный рандомизированный алгоритм даёт O(n log n) ожидаемо — это уже территория https://courses.digitable.life/post/algorithms/14-randomized-algorithms/. Практический инструмент — Qhull, который стоит под scipy.spatial.ConvexHull.

Что оболочка даёт дальше: вращающиеся штангенциркули

Диаметр множества (максимальное расстояние между двумя точками) достигается на вершинах оболочки. Наивно это O(h²), но пара «противоположных» точек движется монотонно — и весь обход становится O(h):

def hull_diameter_sq(hull: list[Point]) -> int:
    """Квадрат диаметра множества по его выпуклой оболочке (CCW, без коллинеарных).
    Rotating calipers: O(h) после построения оболочки."""
    n = len(hull)
    if n < 2:
        return 0
    if n == 2:
        return dist2(hull[0], hull[1])

    best = 0
    j = 1
    for i in range(n):
        ni = (i + 1) % n
        # Двигаем j, пока он «удаляется» от ребра (i, ni):
        # площадь треугольника растёт ровно до самой дальней точки.
        while cross(hull[i], hull[ni], hull[(j + 1) % n]) > cross(hull[i], hull[ni], hull[j]):
            j = (j + 1) % n
        best = max(best, dist2(hull[i], hull[j]), dist2(hull[ni], hull[j]))
    return best

Это ровно та же идея двух указателей, что и в линейных задачах: указатель j за весь цикл делает суммарно O(h) шагов, потому что никогда не идёт назад. Тем же приёмом решаются минимальный охватывающий прямоугольник (он всегда содержит ребро оболочки — теорема Фримена–Шапиры), ширина множества и расстояние между двумя выпуклыми полигонами.

Родственная задача — минимальная охватывающая окружность. Алгоритм Вельцля решает её за ожидаемое O(n) рандомизированным инкрементальным построением; это красивый пример того, как случайность убирает худший случай, — подробности в https://courses.digitable.life/post/algorithms/14-randomized-algorithms/.

Sweep line: главная парадигма

Идея развёртки (заметающей прямой) — самая продуктивная в вычислительной геометрии. Формулировка в одну фразу: сведите двумерную задачу к последовательности одномерных, двигая воображаемую вертикальную прямую слева направо.

Три обязательных компонента:

1. Очередь событий — точки по x, в которых «что-то меняется»: начало объекта, конец объекта, а иногда — новые события, порождаемые самим алгоритмом.
2. Статус развёртки — структура данных, описывающая пересечение прямой с объектами в текущий момент. Обычно упорядоченное множество (сбалансированное дерево), дерево отрезков или BIT.
3. Инвариант — «после обработки всех событий левее x статус корректен, а ответ для этой части уже учтён». Именно инвариант превращает развёртку из эвристики в доказуемый алгоритм.

Это схема алгоритма Бентли–Оттманна (1979) для поиска всех k пересечений n отрезков за O((n + k) log n) времени и O(n) памяти (в варианте Брауна). Наивный перебор пар — O(n²); при n = 100 000 и небольшом числе пересечений разница между «часы» и «секунды».

Ключевая лемма, на которой всё держится: два отрезка, пересекающиеся в точке p, непосредственно перед p являются соседями в порядке статуса. Поэтому достаточно проверять только пары соседей — O(1) проверок на событие вместо O(n).

Если нужен ответ «есть ли хоть одно пересечение», задача проще: алгоритм Шеймоса–Хоуи делает то же самое без событий-пересечений за O(n log n), останавливаясь на первой найденной паре. Именно это стоит за проверкой «валиден ли полигон» в GEOS/PostGIS (ST_IsValid).

Развёртка на практике: площадь объединения прямоугольников

Задача из реальной жизни (viewability рекламы, покрытие экрана элементами, покрытие области зонами доставки): даны n осепараллельных прямоугольников, найти площадь их объединения. Включения-исключения дают 2^n — нежизнеспособно. Развёртка:

def union_area(rects: list[tuple[int, int, int, int]]) -> int:
    """Площадь объединения осепараллельных прямоугольников (x1, y1, x2, y2).
    Простой вариант: сжатие координат по x + слияние интервалов в каждой полосе.
    Время O(n² log n), память O(n). Достаточно до ~5000 прямоугольников."""
    xs = sorted({x for (x1, _, x2, _) in rects for x in (x1, x2)})
    total = 0

    for i in range(len(xs) - 1):
        x_lo, x_hi = xs[i], xs[i + 1]
        width = x_hi - x_lo
        if width == 0:
            continue
        # Прямоугольники, накрывающие всю полосу целиком.
        active = sorted((y1, y2) for (a, y1, b, y2) in rects if a <= x_lo and b >= x_hi)

        covered, cur_end = 0, None
        for y1, y2 in active:                 # слияние отсортированных интервалов
            if cur_end is None or y1 > cur_end:
                covered += y2 - y1
                cur_end = y2
            elif y2 > cur_end:
                covered += y2 - cur_end
                cur_end = y2
        total += covered * width
    return total

Как довести до O(n log n): вместо пересборки интервалов в каждой полосе держать дерево отрезков с «счётчиком покрытия» — специализированную структуру, где в узле хранятся count (сколько раз интервал накрыт целиком) и covered_len (длина покрытой части). События — левые (+1) и правые (−1) стороны прямоугольников. Обновление и запрос корня — O(log n). Эта структура (её называют «дерево Кletka»/segment tree on coordinates) — классический спутник развёртки; аналогичная схема считает периметр объединения и объём объединения кубов.

Ближайшая пара точек

Каноническая задача «разделяй и властвуй» (см. https://courses.digitable.life/post/algorithms/05-recursion-and-divide-conquer/), которую часто решают и развёрткой. Наивно — O(n²). Классическое решение — O(n log n):

def closest_pair(points: list[Point]):
    """Ближайшая пара точек, divide & conquer.
    Возвращает (квадрат расстояния, пара). Время O(n log n), память O(n)."""
    pts = sorted(points)                       # по x, затем по y

    def better(a, b):
        return a if a[0] <= b[0] else b

    def rec(lo: int, hi: int):
        # Побочный эффект: после вызова pts[lo:hi] отсортирован по y.
        if hi - lo <= 3:
            best = (float('inf'), None)
            for i in range(lo, hi):
                for j in range(i + 1, hi):
                    best = better(best, (dist2(pts[i], pts[j]), (pts[i], pts[j])))
            pts[lo:hi] = sorted(pts[lo:hi], key=lambda p: p[1])
            return best

        mid = (lo + hi) // 2
        mid_x = pts[mid][0]
        best = better(rec(lo, mid), rec(mid, hi))

        # Слияние двух y-отсортированных половин — ровно шаг merge sort.
        left, right = pts[lo:mid], pts[mid:hi]
        merged, i, j = [], 0, 0
        while i < len(left) and j < len(right):
            if left[i][1] <= right[j][1]:
                merged.append(left[i]); i += 1
            else:
                merged.append(right[j]); j += 1
        merged.extend(left[i:]); merged.extend(right[j:])
        pts[lo:hi] = merged

        # Полоса шириной 2d вокруг разделяющей прямой.
        strip = [p for p in merged if (p[0] - mid_x) ** 2 < best[0]]
        for i in range(len(strip)):
            for j in range(i + 1, len(strip)):
                if (strip[j][1] - strip[i][1]) ** 2 >= best[0]:
                    break                      # дальше по y — только хуже
                best = better(best, (dist2(strip[i], strip[j]), (strip[i], strip[j])))
        return best

    return rec(0, len(pts))

Почему внутренний цикл константен. В полосе шириной 2d любые две точки одной половины находятся на расстоянии ≥ d. В прямоугольник 2d × d можно упаковать не более 8 таких точек (классическая оценка; аккуратный подсчёт даёт 6 для сравнений вперёд). Значит, break по y срабатывает после O(1) итераций, и слияние линейно. Рекуррентность T(n) = 2T(n/2) + O(n) даёт O(n log n) по мастер-теореме.

Частая ошибка реализации — сортировать полосу по y заново на каждом уровне. Это превращает алгоритм в O(n log² n). Приём выше (возвращать половины уже отсортированными по y и сливать их) — тот же трюк, что в merge sort из https://courses.digitable.life/post/algorithms/02-sorting/.

Триангуляция Делоне и диаграмма Вороного

Две конструкции, двойственные друг другу, и обе — рабочие лошадки прикладной геометрии.

Диаграмма Вороного множества точек (сайтов) разбивает плоскость на ячейки: ячейка сайта s — множество точек плоскости, для которых s ближайший сайт. Триангуляция Делоне — двойственный граф: сайты соединены ребром, если их ячейки соседние. Эквивалентное определение Делоне: триангуляция, в которой описанная окружность любого треугольника не содержит других точек внутри (empty circumcircle property).

Почему это важно:

• Делоне максимизирует минимальный угол среди всех триангуляций набора точек — то есть даёт «наименее вытянутые» треугольники. Для численных методов (FEM), интерполяции высот и генерации сеток это критично.
• Вороной даёт готовые ответы на «ближайший сосед», зоны обслуживания, зоны покрытия вышек, максимально удалённые от препятствий маршруты для робота.
• Минимальное остовное дерево евклидова множества точек — подграф Делоне. Это позволяет считать EMST за O(n log n) вместо O(n²) (см. https://courses.digitable.life/post/algorithms/10-mst-and-flows/).

Как это строят. Три подхода: развёртка Форчуна (O(n log n), изящный и сложный в реализации), инкрементальный рандомизированный алгоритм (ожидаемо O(n log n)), и самый концептуально красивый — подъём на параболоид: отобразите каждую точку (x, y) в (x, y, x² + y²) в трёхмерии. Тогда нижняя часть трёхмерной выпуклой оболочки этих точек, спроецированная обратно на плоскость, есть в точности триангуляция Делоне. Условие «точка внутри описанной окружности» превращается в «точка ниже плоскости» — то есть предикат incircle оказывается предикатом orient3d. Это объясняет, почему scipy.spatial.Delaunay — просто обёртка над Qhull.

Практически: scipy.spatial.Delaunay / Voronoi, CGAL Triangulation_2, JTS/GEOS DelaunayTriangulationBuilder, PostGIS ST_DelaunayTriangles и ST_VoronoiPolygons.

Пространственный поиск: как это работает в проде

Всё вышеперечисленное — алгоритмы «посчитать один раз». В проде чаще нужны запросы: «что попадает в этот прямоугольник», «кто из миллиона объектов ближайший». Здесь работают структуры данных.

Что важно понимать про каждую:

• k-d дерево. Делит пространство поочерёдно по осям. Прекрасно в 2–3 измерениях, быстро деградирует с ростом размерности («проклятие размерности»): при d ≳ 20 запрос kNN просматривает почти всё дерево, и линейный скан оказывается быстрее. Для ML-задач вместо него берут приближённый поиск (HNSW, IVF-PQ).
• R-дерево. Хранит минимальные ограничивающие прямоугольники (MBR) в сбалансированном дереве со страницами под размер блока диска. Это то, что реально стоит под CREATE INDEX ... USING GIST в PostgreSQL/PostGIS и под rtree в SQLite. Ключевой паттерн запроса — двухфазная фильтрация: индекс по MBR отсеивает кандидатов дёшево (&&), затем точный геометрический предикат (ST_Intersects) проверяет выживших. Один шаг стоит O(1) по индексу, второй — дорогой, но их мало.
• Геохеши и S2/H3. Кривая заполнения пространства (Гильберта или Z-order/Мортона) превращает 2D-координату в одномерный ключ так, что близкие точки чаще имеют общий префикс. Это позволяет использовать обычный отсортированный KV-store как пространственный индекс — именно так масштабируются гео-сервисы поверх Bigtable/Cassandra/DynamoDB. Библиотеки: Google S2 (сферические ячейки, корректная работа на сфере) и Uber H3 (гексагональная сетка, удобная для агрегаций и спроса/предложения). Важный нюанс: у геохешей есть краевой эффект — соседние точки по разные стороны границы ячейки имеют совсем разные префиксы, поэтому запрос «в радиусе R» обязан просматривать соседние ячейки, а не только свою.

Другие продакшн-инструменты

• Упрощение линий. Douglas–Peucker (O(n log n) в среднем, O(n²) худший) — стандарт для генерализации контуров под масштаб карты. Осторожно: он не гарантирует сохранение топологии, соседние полигоны могут «разъехаться» и породить щели. Топологически корректный вариант — упрощать общие рёбра один раз (TopoJSON).
• Булевы операции над полигонами. Отсечение по выпуклому окну — Sutherland–Hodgman (простой, O(n)), общий случай — алгоритм Ватти/Greiner–Hormann, промышленная реализация — Clipper2 (работает в целых координатах именно ради робастности) и GEOS OverlayNG.
• Broad phase в физике. Sweep and prune — та же развёртка: сортируем AABB по одной оси, поддерживаем список активных. На кадрах со слабым движением сортировка почти отсортированного массива стоит O(n) (insertion sort), что делает подход крайне дешёвым — практическая деталь в духе https://courses.digitable.life/post/algorithms/18-practical-optimization/.

Типичные ошибки

1. Считать в double там, где хватает целых. Самая частая и самая дорогая ошибка. Умножьте вход на масштабный коэффициент и работайте в целых.
2. atan2 для сортировки по углу. Медленно, неточно, ломает компаратор. Сортируйте по полуплоскости + cross.
3. Забыть про коллинеарные точки в оболочке. Решите заранее, нужны они или нет, и напишите это в комментарии рядом с <= 0.
4. Единый абсолютный EPS. Он либо слишком мал, либо слишком велик — почти никогда не «в самый раз» для реальных диапазонов координат.
5. Ray casting без полуоткрытого интервала. Луч, попавший в вершину, даёт двойной счёт и случайный ответ.
6. Пересортировка полосы в closest pair. Лишний log n, который никто не замечает, пока n не станет большим.
7. Игнорирование вырожденных входов. Все точки совпадают; все точки на одной прямой; полигон из двух вершин; отрезок нулевой длины. В геометрии это не экзотика, а типичные данные (особенно после квантования координат).
8. Забыть, что Земля — не плоскость. Для больших расстояний и полярных областей плоские формулы дают ошибки в километры. Используйте сферические/эллипсоидальные вычисления (ST_Distance с geography, S2) — и помните, что «выпуклая оболочка» на сфере определена иначе.
9. Полагаться на самопересекающиеся полигоны. Невалидная геометрия в PostGIS ведёт себя непредсказуемо: сначала ST_IsValid/ST_MakeValid, потом операции.
10. Строить оболочку там, где нужна вогнутая форма. Для «формы облака точек» выпуклая оболочка часто слишком груба — нужны alpha-shapes или concave hull (ST_ConcaveHull).

Мини-итог: сводка сложностей

Главное, что стоит унести:

• Один предикат. Знак cross — это 80% плоской геометрии. Научитесь видеть задачу как последовательность вопросов «слева или справа».
• Сортировка — вход в решение. Почти каждый эффективный геометрический алгоритм начинается с сортировки: по x, по углу, по времени события. O(n log n) в таблице выше почти всегда означает «сортировка плюс линейный проход».
• Развёртка — универсальный редуктор размерности. Если задача двумерная и наивное решение квадратично, первый вопрос: «что если двигать прямую слева направо и держать структуру данных на пересечении?»
• Робастность — не деталь реализации, а часть постановки задачи. Решите вопрос об арифметике до того, как напишете первую строку.

Источники

• T. Cormen, C. Leiserson, R. Rivest, C. Stein. Introduction to Algorithms, 4-е изд., глава «Computational Geometry» — MIT Press
• M. de Berg, O. Cheong, M. van Kreveld, M. Overmars. Computational Geometry: Algorithms and Applications, 3-е изд. — главный учебник по теме: Springer
• J. O’Rourke. Computational Geometry in C, 2-е изд. — с упором на реализацию и вырожденные случаи: Cambridge
• S. Skiena. The Algorithm Design Manual, раздел про геометрические задачи — algorist.com
• A. M. Andrew. Another efficient algorithm for convex hulls in two dimensions, 1979 — DOI
• J. L. Bentley, T. A. Ottmann. Algorithms for reporting and counting geometric intersections, 1979 — DOI
• S. Fortune. A sweepline algorithm for Voronoi diagrams, 1987 — DOI
• T. M. Chan. Optimal output-sensitive convex hull algorithms in two and three dimensions, 1996 — DOI
• E. Welzl. Smallest enclosing disks (balls and ellipsoids), 1991 — DOI
• J. R. Shewchuk. Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, 1997 — CMU
• D. Douglas, T. Peucker. Algorithms for the reduction of the number of points required to represent a digitized line, 1973 — DOI
• Документация CGAL, GEOS, Shapely, PostGIS, Qhull
• S2 Geometry и Uber H3 — гео-индексация в проде
• Разборы и эталонные реализации: cp-algorithms.com — геометрия

Смежные статьи трека: сортировка как предобработка — https://courses.digitable.life/post/algorithms/02-sorting/; бинарный поиск, в том числе по вееру и по ответу — https://courses.digitable.life/post/algorithms/03-searching-and-binary-search/; два указателя, из которых растут вращающиеся штангенциркули — https://courses.digitable.life/post/algorithms/04-two-pointers-sliding-window/; разделяй и властвуй для ближайшей пары — https://courses.digitable.life/post/algorithms/05-recursion-and-divide-conquer/; рандомизация в алгоритме Вельцля и инкрементальном Делоне — https://courses.digitable.life/post/algorithms/14-randomized-algorithms/; кэш и SIMD для геометрических ядер — https://courses.digitable.life/post/algorithms/18-practical-optimization/.

Что дальше

Мы несколько раз упирались в арифметику: точность предикатов, целочисленные координаты, теорема Пика, gcd для подсчёта точек на отрезке. Всё это — территория теории чисел, которая в алгоритмах отвечает не только за геометрию, но и за хеширование, криптографию и быстрые преобразования: Теория чисел и математика для алгоритмов и криптографии.