Многоугольник: площадь шнурком, направление обхода и точка внутри

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

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

Участок в форме буквы Г: вершины (0, 0), (6, 0), (6, 2), (2, 2), (2, 5), (0, 5). Площадь ты сосчитаешь в уме: прямоугольник 6 × 2 и сверху 2 × 3, всего 18. А как это посчитает программа, которой дали только список вершин? Многоугольник может быть каким угодно невыпуклым, и резать его на прямоугольники программа не умеет.

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

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

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

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

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

Формула шнурка. Вершины P₀, P₁, …, Pₙ₋₁, после последней снова идёт P₀. Удвоенная площадь:

2S  =  P₀ × P₁  +  P₁ × P₂  +  …  +  Pₙ₋₁ × P₀
       где  P × Q  =  Px · Qy − Py · Qx

Почему именно так. Pᵢ × Pᵢ₊₁ - это поворот трёх точек O, Pᵢ, Pᵢ₊₁, где O - начало координат, то есть удвоенная ориентированная площадь треугольника O Pᵢ Pᵢ₊₁. Получается веер треугольников из O на каждую сторону. Сторона, которую обход проходит против часовой стрелки, если смотреть из O, даёт плюс, сторона, пройденная по часовой, - минус, и всё, что лежит вне многоугольника, покрыто плюсом и минусом поровну. Если O лежит внутри выпуклого многоугольника, вычитать нечего - сумма просто складывается из треугольников, режущих фигуру. От выбора O ответ не зависит: сдвиг всех вершин не меняет площадь.

Знак суммы - направление обхода: против часовой стрелки площадь положительна, по часовой - отрицательна, ровно как у поворота трёх точек.

Луч. Из точки p проводим луч вправо, параллельно оси x. Ребро ab пересекает луч, если один конец ребра выше p, а другой - нет, и точка пересечения лежит правее p. Второе проверяется поворотом без деления: для ребра, идущего вверх, p должна быть слева от него, а для ребра, идущего вниз, - справа: направление ребра перевернулось, и вместе с ним «лево» и «право».

Правило полуоткрытого ребра. «Выше» - строго: вершина на уровне луча считается лежащей НИЖЕ него. Тогда через вершину, где луч проходит насквозь, засчитывается ровно одно ребро из двух, а вершина, которой луч касается, даёт ноль или два - чётное число, и ответ не меняется.

Правило

Удвоенная площадь многоугольника - сумма Pᵢ × Pᵢ₊₁ по всем сторонам, знак - направление обхода (плюс - против часовой стрелки). Точка внутри, если луч из неё пересекает нечётное число рёбер; ребро засчитывается, когда ровно один его конец СТРОГО выше точки и пересечение правее точки. Точки границы проверяются отдельно - лучом их не различить.

Разбор: буква Г

def udv_ploshchad(poly):
    s = 0
    n = len(poly)
    for i in range(n):
        a = poly[i]
        b = poly[(i + 1) % n]
        s = s + a[0] * b[1] - a[1] * b[0]
    return s

bukva_g = [[0, 0], [6, 0], [6, 2], [2, 2], [2, 5], [0, 5]]
s = udv_ploshchad(bukva_g)
print(s, s / 2)
naoborot = bukva_g[::-1]
print(udv_ploshchad(naoborot))
sdvig = []
for p in bukva_g:
    sdvig.append([p[0] + 100, p[1] - 40])
print(udv_ploshchad(sdvig))

36 18.0
-36
36

Удвоенная площадь 36, площадь 18 - как у прямоугольников, хотя программа о них ничего не знала. Вершины перечислены против часовой стрелки, поэтому знак плюс. Тот же контур в обратном порядке даёт −36: площадь та же, обход другой. Сдвиг на (100, −40) не изменил ничего - начало координат в формуле лишь опорная точка веера, и куда её поставить, всё равно.

Модуль нужен, когда важна только площадь, знак - когда важно направление. Например, форматы карт требуют определённого обхода: в GeoJSON (RFC 7946) внешний контур перечисляется против часовой стрелки, а дыры - по часовой; в формате Shapefile наоборот. Знак формулы шнурка - самый дешёвый способ это проверить и развернуть контур.

Разбор: ромб, вершина на луче и формула Пика

vnutri - луч с правилом полуоткрытого ребра. naivno - то, что пишут первым: ребро засчитывается, если y точки между y концов включительно. na_granice - отдельная проверка границы.

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

def vnutri(poly, p):
    schet = 0
    n = len(poly)
    for i in range(n):
        a = poly[i]
        b = poly[(i + 1) % n]
        if (a[1] > p[1]) != (b[1] > p[1]):
            z = povorot(a, b, p)
            if (b[1] > a[1] and z > 0) or (b[1] < a[1] and z < 0):
                schet = schet + 1
    return schet % 2 == 1

def naivno(poly, p):
    schet = 0
    n = len(poly)
    for i in range(n):
        a = poly[i]
        b = poly[(i + 1) % n]
        lo = min(a[1], b[1])
        hi = max(a[1], b[1])
        if a[1] != b[1] and lo <= p[1] <= hi:
            x = a[0] + (p[1] - a[1]) * (b[0] - a[0]) / (b[1] - a[1])
            if x > p[0]:
                schet = schet + 1
    return schet % 2 == 1

def na_granice(poly, p):
    n = len(poly)
    for i in range(n):
        a = poly[i]
        b = poly[(i + 1) % n]
        if povorot(a, b, p) == 0:
            if min(a[0], b[0]) <= p[0] <= max(a[0], b[0]):
                if min(a[1], b[1]) <= p[1] <= max(a[1], b[1]):
                    return True
    return False

romb = [[0, 2], [2, 0], [4, 2], [2, 4]]
for p in [[1, 2], [2, 3], [5, 2], [3, 3]]:
    print(p, vnutri(romb, p), naivno(romb, p), na_granice(romb, p))

bukva_g = [[0, 0], [6, 0], [6, 2], [2, 2], [2, 5], [0, 5]]
vnutr = 0
gran = 0
for x in range(7):
    for y in range(6):
        if na_granice(bukva_g, [x, y]):
            gran = gran + 1
        elif vnutri(bukva_g, [x, y]):
            vnutr = vnutr + 1
print(vnutr, gran, vnutr + gran / 2 - 1)

[1, 2] True False False
[2, 3] True True False
[5, 2] False False False
[3, 3] False False True
8 22 18.0

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

Последняя строка - проверка формулы Пика: у простого многоугольника без дыр с вершинами в целых точках площадь равна I + B/2 − 1, где I - число целых точек строго внутри, B - на границе. Перебор всех узлов сетки нашёл 8 внутренних и 22 граничных: 8 + 11 − 1 = 18, та же площадь, что дал шнурок. Граничные точки можно было и не перебирать: на стороне от (x₁, y₁) до (x₂, y₂) целых точек, не считая одного конца, ровно НОД(|x₂ − x₁|, |y₂ − y₁|): если g - этот НОД, то шаг (dx/g, dy/g) - наименьший целый шаг вдоль стороны, и таких шагов ровно g. Для буквы Г это 6 + 2 + 4 + 3 + 2 + 5 = 22. Почему формула Пика именно такая: узел внутри «отвечает» за единицу площади, узел на границе - только за половину (половина его окрестности снаружи), а углы в вершинах в сумме недобирают ровно одну целую - сумма внешних углов многоугольника составляет полный оборот.

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

Луч через вершину. Ошибка не случайная: точки на одной высоте с вершинами в реальных данных встречаются постоянно - сетки, выровненные по координатам здания, кварталы с прямоугольной застройкой. Наивный счёт ошибается ровно на них.

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

Самопересекающийся контур. У «бабочки» (0, 0), (4, 4), (4, 0), (0, 4) формула шнурка даёт 0: два треугольника обходятся в разные стороны и взаимно гасятся. У пятиконечной звезды, нарисованной одним росчерком, расходятся уже правила «внутри»: луч из центра пересекает два ребра - по чётности центр снаружи, а контур обматывает центр дважды - по правилу «обмотка не ноль» он внутри. В SVG это видно глазами: атрибут fill-rule со значением evenodd оставляет в центре звезды дыру, nonzero - заливает.

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

Сервис доставки получает адрес и должен понять, в какую зону он попадает, - это точка в многоугольнике на каждый заказ, ST_Contains в PostGIS. Карта района прислана в GeoJSON, а база ждёт другой обход - знак площади это ловит до того, как построятся вывернутые полигоны. А если надо сосчитать целые точки внутри многоугольника - узлы сетки, пиксели с целыми центрами, - формула Пика даёт ответ по площади и границе, без перебора.

Теперь сам

1. Площадь четырёхугольника (1, 0), (5, 1), (4, 4), (0, 3) по формуле шнурка. В какую сторону перечислены вершины?
2. Для квадрата (0, 0), (4, 0), (4, 4), (0, 4) и точки (2, 0) что скажет vnutri? А для (2, 4)?
3. Треугольник (0, 0), (4, 0), (0, 4). Сколько целых точек строго внутри?