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

В центре Алматы 7909 перекрёстков, районов в городе восемь. Каждому перекрёстку нужен свой район: кто чинит светофор, кто убирает снег. У программы есть только координаты перекрёстка и границы районов — ломаные из сотен точек. Ниже все перекрёстки раскрашены по районам, и каждый цвет поставил один и тот же тест: луч из точки на восток и подсчёт, сколько раз он пересёк границу. Коснитесь карты в любом месте.

Тест разложил по районам 7445 перекрёстков, ни один не попал в два района сразу. Ещё 459 стоят точно на границе — до метра, — и пять не попали ни в один район. Ни то ни другое не ошибка теста. Границы районов в OpenStreetMap во многих местах проведены прямо через перекрёстки, а на широких проспектах соседние районы делят улицу не по одной линии, и между их границами остаётся ничья полоса. Обо всём этом — ниже, а начнём с того, почему луч вообще работает.

Идея

Идите из точки по прямой в одну сторону, пока не выйдете за город. Каждое пересечение границы района меняет «внутри» на «снаружи» и обратно, а в конце пути вы точно снаружи. Значит, если пересечений было нечётное число, начинали вы внутри.

Проведите из точки луч и посчитайте, сколько раз он пересекает границу многоугольника. Нечётное число — точка внутри, чётное — снаружи. Сторона AB пересекает луч на восток из точки P, если ровно один её конец лежит выше P и точка пересечения восточнее P; второе проверяет знак векторного произведения, без деления.

Это правило чёт-нечет, или тест луча. Направление луча любое; на восток удобнее всего: пересечение со стороной ищется на прямой y = P_y. Сторона из A в B её пересекает, если один конец выше, а другой нет. А восточнее ли точки пересечение, говорит знак векторного произведения cross(A, B, P): если сторона идёт вверх, пересечение восточнее, когда P слева от неё, если вниз — когда справа. Ни дробей, ни деления: для целых координат ответ точный.

Почему это верно

За тестом стоит теорема Жордана (1887): простая замкнутая кривая делит плоскость на две части, внутреннюю и внешнюю. В общем виде её доказательство трудное, а для многоугольника его даёт сам луч. Двигайте точку, не задевая границу. Число пересечений её луча меняется, только когда луч проходит через вершину, и тогда меняется на 0 или на 2: две стороны при вершине либо сменяют одна другую, либо появляются на луче и исчезают с него вместе. Значит, чётность постоянна на любом куске плоскости, который можно обойти, не пересекая границу. Когда же точка пересекает сторону, её луч теряет или приобретает ровно одно пересечение, и чётность меняется. Поэтому точки делятся на два класса, а далеко от многоугольника пересечений ноль — там «снаружи». Освальд Веблен в 1905 году упрекал Жордана как раз в этом: тот принял теорему для многоугольников без доказательства.

Луч через вершину

Первая ловушка — луч, который проходит точно через вершину. Если считать сторону, когда y точки лежит в отрезке [min, max] её концов включительно, луч заденет обе стороны при вершине. Когда граница проходит через луч насквозь, это одно пересечение, а насчитано два, и ответ перевернётся. На карте городских районов это не редкость, а правило: координаты целые, границы проведены через перекрёстки, и лучи из соседних перекрёстков то и дело упираются в вершины. Включите в сцене «Наивный луч»: такой подсчёт ошибается на 355 перекрёстках из 7909 — каждом двадцать втором.

Лечится это одним условием: сторона считается, если ровно один её конец строго выше луча. Вершина на луче при этом как будто чуть ниже него. Если граница проходит через луч насквозь, одна сторона при вершине уходит вверх и посчитана, другая — вниз и не посчитана: одно пересечение. Если граница только касается луча, обе стороны по одну сторону от него: посчитаны обе или ни одной, чётность не меняется. Горизонтальная сторона на луче не считается никогда — её соседи справятся сами. Попробуйте на короне ниже: точка P стоит на высоте впадин между зубцами, а правый бок короны ломается как раз на этой высоте.

Точка на границе

Вторая ловушка — точка на самой стороне. Луч из неё начинается на границе, и чётность ничего не значит: ответ зависит от того, как правило округлит пересечение в самой точке. Поэтому «на границе» проверяют отдельно и первым делом: P лежит на стороне AB, если cross(A, B, P) = 0 и P внутри прямоугольника, натянутого на A и B. В Алматы таких перекрёстков 459: 222 лежат на границе одного района, 235 — на общей границе двух, и два — в углу, где сходятся три района. Все 459 — вершины границ: в OpenStreetMap границу района во многих местах провели прямо через узлы дорог.

Если соседние многоугольники делят стороны в точности, от точек на границе можно избавиться договорённостью. Так делают видеокарты: пиксель, центр которого лежит на общей стороне двух треугольников, рисует ровно один из них — в Direct3D это правило верхнего левого ребра. С районами Алматы так не выйдет: соседние районы делят широкие проспекты не по одной линии. На пересечении проспектов Абая и Абылай хана граница Алмалинского района проходит через перекрёсток с северной проезжей частью Абая, а граница Бостандыкского — через перекрёсток с южной, в пятнадцати метрах. Разделительная полоса между ними не принадлежит ни одному району. Пять перекрёстков, которые тест не отдал никому, стоят на таких полосах — выберите в сцене пример «Ничья полоса на проспекте Абая». Тест здесь прав, неточны данные.

Число оборотов

У луча есть второй способ считать. Будем различать, в какую сторону граница пересекает луч: снизу вверх — +1, сверху вниз — −1. Сумма — это число оборотов границы вокруг точки: сколько раз она обходит точку против часовой стрелки за вычетом обходов по часовой. Если стоять в точке и всё время смотреть на бегущего по границе, направление взгляда проходит через восток ровно тогда, когда бегущий пересекает луч, и поворачивается при этом против часовой стрелки, если тот бежит вверх. Каждый полный оборот проходит через восток один раз.

Для простого многоугольника число оборотов внутри — +1 или −1, смотря по направлению обхода, снаружи — 0, и ответ тот же, что у чётности: на всех перекрёстках Алматы, кроме лежащих на границе, оба теста согласны. Разница видна, когда граница пересекает сама себя. Выберите в сцене пентаграмму: из её центра луч пересекает две стороны, чётно — «снаружи», но обе пересекает снизу вверх, число оборотов 2 — «внутри». Отсюда два правила заливки, которые есть в SVG и у холста браузера: evenodd закрашивает точку при нечётном числе пересечений, nonzero — при ненулевом числе оборотов. По умолчанию и там и там действует nonzero. Петля в той же сцене обходит середину дважды, а правое ушко — по часовой стрелке: перетащите точку туда, и число оборотов станет −1.

Раньше число оборотов считали суммой углов, под которыми из точки видны стороны, через арктангенсы. Подсчёт пересечений со знаками даёт то же целое число без тригонометрии. Именно этот вариант разбирают Дэн Санди (2001) и Кай Хорманн с Александром Агатосом в статье того же года о точке в многоугольнике с самопересечениями.

Код

Три теста на маленьких фигурах: луч с полуоткрытым правилом и проверкой границы, наивный подсчёт и число оборотов. Корона — та же, что в сцене: наивный подсчёт ошибается для трёх точек на высоте 3, а для точки во впадине отвечает «внутри», хотя она на границе. Пентаграмма — для сравнения чётности и оборотов.

А это первая сцена в цифрах: все перекрёстки через все районы, с подсчётом того, сколько сторон пришлось проверить.

Много точек: полосы

Тест луча смотрит на все n сторон, хотя пересечь луч могут только те, что тянутся по высоте через P_y. Если точек много, стороны стоит разложить заранее. Разрежем район на горизонтальные полосы одинаковой высоты — столько полос, сколько сторон, — и запишем каждую сторону во все полосы, которые она задевает. Луч из точки лежит в одной полосе, и смотреть нужно только на её стороны: любая сторона, которая может его пересечь, там есть. В районах Алматы сторона — это десятки метров границы, и в полосе их обычно единицы. Листинг выше насчитал 16 429 057 проверок сторон без полос и 61 049 с ними — в 269 раз меньше. Подготовка — один проход по сторонам.

Перебор растёт с числом сторон линейно, полосы — куда медленнее: на десяти тысячах сторон они быстрее в сотни раз на Python и в десятки раз на Go. Слабое место — форма. Если многие стороны тянутся через весь многоугольник по высоте, как зубья гребёнки, они попадают в каждую полосу, и выигрыш пропадает. Гарантию O(log n) на точку при любой форме дают структуры для разбиения плоскости: полосы Добкина и Липтона (1976), которые режут плоскость через каждую вершину и ищут двоичным поиском дважды, и их потомки, экономящие память. Если же точки сами по себе — без многоугольника под рукой, — их сначала собирают в k-d дерево или R-дерево и спрашивают, какие из них попали в рамку района, а точный тест делают только для них.

Сколько это стоит

Одна точка — O(n) времени и O(1) памяти: одна проверка на сторону, две разности и два умножения в ней. Для простого многоугольника быстрее без подготовки нельзя: не посмотрев на какую-то сторону, нельзя знать, не проходит ли она между точкой и остальной границей. С полосами — O(n + S) на подготовку, где S — сколько раз стороны легли в полосы, и O(1 + k) на точку, где k — сторон в её полосе.

Точность — как у всей геометрии на векторных произведениях. В целых числах ответ точный: при координатах до 10⁹ разности до 2·10⁹, произведения до 4·10¹⁸, и всё помещается в int64. Поэтому перекрёстки и границы в этой статье переведены в целые метры. В числах с плавающей точкой точки у самой границы могут оказаться по любую её сторону — ловушка та же, что у поворота налево или направо. Если ответ «на границе» важен, координаты переводят в целые или считают с допуском и честно отвечают «слишком близко к границе».

Где работает

• Карты и доставка. В каком районе дом, попадает ли адрес в зону доставки, какой участок кадастра под точкой — тест луча после грубого отбора по рамке, как в листинге.
• Векторная графика. Правила заливки fill-rule в SVG и fill() холста — это чётность пересечений и число оборотов. Буква «О» — два контура, и дырка в ней остаётся дыркой по любому из правил, если внутренний контур обходится в обратную сторону: число оборотов там 1 − 1 = 0, пересечений два.
• Интерфейсы и игры. Попал ли клик в страну на карте или в фигуру на экране. В этой статье так выбирается район в обеих сценах.
• Растеризация. Видеокарта решает «внутри ли центр пикселя» для каждого треугольника, а правило верхнего левого ребра — что делать с пикселями на общей стороне.

История

В августе 1962 года журнал Communications of the ACM напечатал процедуру Моше Шимрата «Алгоритм 112: положение точки относительно многоугольника» — это самая ранняя публикация теста луча, которую обычно приводят. В декабре журнал опубликовал её проверку другим автором: тогда процедуры из раздела алгоритмов было принято перепроверять. С тех пор тест переоткрывали много раз, а спотыкались всё на тех же вершинах. В 1990 году в рассылке Ray Tracing News рассказывали о решении, которое при попадании луча в вершину предлагало просто выпустить другой луч. Марк Вандеветтеринг ответил, что проще договориться считать вершину на луче лежащей по одну сторону от него — и ответ выйдет верным. На вопрос, чья это идея, он написал, что её показал ему Эрик Хейнс, который на авторство тоже не претендует. Число оборотов без тригонометрии описал Дэн Санди в 2001 году; в том же году Кай Хорманн и Александр Агатос разобрали в журнале Computational Geometry точку в многоугольнике произвольной формы, с самопересечениями.

Задачи

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