Локализация точки в выпуклом многоугольнике

  • Tutorial
Листая страницы хаба «Алгоритмы», наткнулся на топик, посвященный решению задачи локализации точки в многоугольнике: задан многоугольник (замкнутая ломаная линия без самопересечений), требуется определить — находится ли заданная точка A внутри этого многоугольника или нет. В одном из последних комментариев к топику было высказано недоумение, какое отношение такая чисто математическая задача имеет к теории алгоритмов. Имеет-имеет, причем самое непосредственное. Задача локализации является классической задачей вычислительной геометрии (не путать с компьютерной графикой). В качестве разминки предлагается взглянуть на картинку справа, на которой изображен многоугольник типа кривой Пеано (источник [1]), и попытаться ответить на вопрос — красная точка ты видишь суслика? и я не вижу, а он есть! находится внутри или снаружи многоугольника? А ниже мы (исключительно в образовательных целях) рассмотрим простую вариацию данной задачи, когда заданный многоугольник является выпуклым.


Понятно, что если под многоугольником понимается треугольник или четырехугольник, то никаких алгоритмов в полном смысле этого слова не возникает, а возникают три-четыре формулы, которые и надо запрограммировать. Однако, совсем другая история начинается, если число вершин многоугольника оказывается большим. Например, миллион или сто миллионов вершин. Такая задача кажется уже вполне себе алгоритмической. Наивный алгоритм (запуск луча из заданной точки и подсчет числа его пересечений со сторонами многоугольника) является линейным по числу вершин многоугольника n, т.к. мы должны проверить пересечение луча со всеми сторонами многоугольника. Соответственно возникает вопрос, а можно ли сделать это быстрее, т.е. имеется ли сублинейный алгоритм решения задачи локализации в данной постановке? Правильный ответ на этот вопрос — да, задача может быть решена за время O(log2n). Подробности решения в общем случае, когда многоугольник не обязан быть выпуклым, можно найти в замечательной книжке [2]. А ниже мы рассмотрим, как можно применить алгоритм двоичного поиска к задаче локализации точки в выпуклом многоугольнике.

Небольшой ликбез


Пусть заданы три точки A, B и C. Предположим, что мы смотрим из точки A в точку B. Где при этом окажется точка C — справа или слева относительно направления нашего взгляда? Ответить на этот вопрос можно с помощью векторного произведения векторов a=AB и b=BC, точнее с помощью z-координаты такого произведения, которая вычисляется по простой формуле z=axby-aybx. Если z>0, то искомый поворот левый, если z<0 — то правый. Если же монета упала на ребро z=0, то три точки лежат на одной прямой (говорят, что векторы a и b коллинеарны).

Код на питоне, вычисляющий направление поворота, элементарен (точки представляются списками длины 2):
def rotate(A,B,C):
  return (B[0]-A[0])*(C[1]-B[1])-(B[1]-A[1])*(C[0]-B[0])

Имея такую полезную функцию, можно делать много добрых дел. Например, можно определять, пересекаются ли два заданных отрезка AB и CD? Это, кстати, еще одна базовая операция вычислительной геометрии, и нам она тоже скоро понадобится (правда, всего один раз). Идея простая — два отрезка пересекаются тогда и только тогда, когда концы одного отрезка лежат по разные стороны другого и наоборот.
Чтобы проверить, располагаются ли точки C и D по разные стороны относительно отрезка (вектора) AB, надо посмотреть на направления поворотов rotate(A,B,C) и rotate(A,B,D). Если знаки этих выражений различны, то прямая AB пересекает отрезок CD (причем во внутренней точке). Знаки двух чисел различны, в том и только в том случае, когда их произведение отрицательно. Следовательно, критерием пересечения отрезков AB и CD является отрицательность двух выражений rotate(A,B,C)*rotate(A,B,D) и rotate(C,D,A)*rotate(C,D,B). Если заменить строгую отрицательность на неположительность, то сможем отлавливать пересечение и в концевых точках отрезков. Нам же потребуется некоторый промежуточный вариант, а именно, одно произведение будет нестрого сравниваться с нулем, второе — строго:
def intersect(A,B,C,D): 
  return rotate(A,B,C)*rotate(A,B,D)<=0 and rotate(C,D,A)*rotate(C,D,B)<0

Причиной такого странного выбора знаков неравенств является особенность будущего применения данной функции (в частности, точки A, B и C будут различными вершинами выпуклого многоугольника, а точка D — точкой, которую мы пытаемся локализовать).

Двоичный поиск


Теперь переходим собственно к локализации. Заданы выпуклый многоугольник P, состоящий из n вершин и точка A. При этом предполагается, что вершины в P пронумерованы против часовой стрелки (говорят, что направление обхода по периметру P положительно).

Идея алгоритма: возьмем первую вершину многоугольника P0 и попытаемся определить, в какой сегмент (угол) PiP0Pi+1 попадает точка A. Для начала проверим, что A попадает в сегмент Pn-1P0P1, если это не так, то A гарантированно лежит вне многоугольника.

Код:
def pointloc(P,A):
  n = len(P)
  if rotate(P[0],P[1],A)<0 or rotate(P[0],P[n-1],A)>0:
    return False

Дальше все как в классическом двоичном поиске — полагаем p=1, r=n-1 (границы текущего сегмента), вычисляем среднюю вершину q=(p+r)/2, смотрим, где находится A относительно вектора P0Pq, если слева, то заменяем r на q, если справа, то заменяем p на q.

Продолжаем этот процесс, пока не окажется, что r-p=1:
  p, r = 1, n-1
  while r-p>1:
    q = (p+r)/2
    if rotate(P[0],P[q],A)<0: r = q
    else: p = q

Теперь осталось проверить, пересекаются ли отрезки P0A и PpPr?
Если пересекаются, то точка A лежит вне, если не пересекаются, то внутри.

Код:
  return not intersect(P[0],A,P[p],P[r])

Вот собственно и все…

Итоги


Нетрудно убедиться, что сложность алгоритма такая же, как и у классического алгоритма двоичного поиска — O(log2n). И все было бы замечательно, если бы этот алгоритм не был заточен под выпуклые многоугольники. Справедливости ради, рассмотренный нами алгоритм работает и с некоторыми невыпуклыми многоугольниками, но весьма специфического вида — в таких многоугольниках все диагонали, ведущие из нулевой вершины должны лежать внутри многоугольника:

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

PS: Задача с картинкой вверху топика легко решается в любом графическом редакторе: выбираем инструмент заливка (fill), берем цвет отличный от фона, кликаем по краю картинки (вне многоугольника) и вуаля!

Литература


  1. P. Prusinkiewicz, A. Lindenmayer, The algorithmic beauty of plants, 1996
  2. M. de Berg, O. Cheong, M. van Kreveld, M. Overmars, Computational Geometry. Algorithms and Applications, 2008
  3. Ф. Препарата, М. Шеймос, Вычислительная геометрия, 1989
Share post

Comments 45

    0
    Пока задачи такой не вставало, но за статью спасибо. Вдруг когда-нибудь понадобится.
      +11
      не помню где, не помню когда, слышал что такая задача решается следующим образом:
      из точки провести прямую и посчитать количество пересечений с границей многоугольника.
      если количество пересечений четное — точка вне многоугольника.
        +2
        провести луч (вероятно вы хотели написать)
          +2
          о, да, конечно луч
            0
            Все правильно, нужно провести луч и посчитать количество пересечений, единственное, что может быть сложным — аккуратно отследить и правильно обработать все возможные случаи пересечения, а их, на сколько я помню, 5 штук.
              0
              Если вершины многоугольника заданы целочисленными координатами (или заданы по какой-то сетке), то достаточно сместить на пол-клетки исходную точку, и проводить луч уже от такой новой точки. Тогда не придется рассматривать различные случаи когда луч проходит по одной или нескольким граням многоугольника. Здесь есть еще пара нюансов в виде того, что такой проход надо делать два раза (X+0.5, X-0.5), но в целом имплементировать и тестить гораздо проще.
                0
                Не всегда все задано так, как вы говорите.
          +1
          Такая реализация — наиболее простая.
            +4
            Например, вы могли (бы) услышать про такой алгоритм из поста выше — цитирую: «запуск луча из заданной точки и подсчет числа его пересечений со сторонами многоугольника», там же чуть ниже объясняется основной недостаток этого алгоритма — его линейность…
              +1
              кроме недостатков, наверное, у такого метода есть и достоинства, ведь так?
              иначе выглядит безальтернативно, и значит тема могла бы быть раскрыта лучше.
                0
                ну если проще, было бы здорово обзор возможных подходов в качестве вступления.
                тогда можно было бы использовать статью как начальную точку отсчета при выборе инструмента
                  +2
                  Было бы здорово, я согласен, но это требует существенно большего времени на подготовку и написание. Основной же тезис данного поста — чисто математическая задача при больших количествах параметров становится объектом теории алгоритмов.
                  +1
                  Конечно, достоинства есть — это универсальность и простота. Ситуация примерно такая, как с поиском числа в массиве. Есть универсальный метод последовательного поиска — работает медленно (на больших массивах), но всегда, а есть двоичный поиск — работает быстро, но только на упорядоченных массивах.
                    +1
                    На правах оффтопа -> а еще есть хеш-тейблы, которые работают еще быстрее, с высокой вероятностью за одну итерацию, но в отличие от двоичного поиска, ищут только совпадения, не могут найти ближайшее меньшее/большее.
                    0
                    Основное достоинство — метод с лучом работает для произвольного многоугольника, а описанный в статье — для выпуклого.
                  0
                  Если многоугольник задан обходом вершин в известном направлении, то достаточно найти ближайшую к началу луча сторону, с которой этот луч пересекается. и посмотреть, в какую сторону она идет — вверх или вниз.
              • UFO just landed and posted this here
                  +3
                  Все-таки многоугольник — это ломаная линия, поэтому он обычно задается именно последовательностью точек, а не набором отрезков. Максимум плохого здесь — эта линия может быть неверно ориентирована (по часовой), но эта проблема легко решается не более чем за линейное время.
                    0
                    Вообще говоря можно и без бинпоиска за O(N). Для _выпуклого_ многоугольника надо просто проверить, чтобы точка лежала по одинаковую сторону от каждой стороны (то бишь, N векторных произведений посчитать).

                    Луч уместно применять, только когда многоугольник невыпуклый, имхо, потому что там без чисел с плавающей точкой не обойтись (?).
                      0
                      Нет, можно обойтись целочисленной арифметикой (при условии, что координаты являются целыми), используя функцию типа intersect, определенную выше. Главная проблема там — учет всяких граничных состояний, но тоже все решается целочисленно.
                    0
                    А не проще ли подсчитать сумму углов? Если она равна 360° то точка внутри, если 0°, то снаружи.
                      0
                      Можно и так, но это 1) опять-таки будет линейно, 2) требует применения тригонометрических функций, которые в отличие от арифметики вычисляются очень и очень долго.
                        0
                        В книге Никулина Е. А. «Компьютерная геометрия и алгоритмы машинной графики» приведен угловой тест с использованием октантов, отпадает необходимость использования тригонометрических функций.
                      +5
                      Странно, что в начале задачи дается задача в общем виде, потом говорится, что задача в общем виде решается за O(log2n), потом дается частное решение в случае выпуклого полигона. Нас где-то обманули.

                      В общем случае, для такого алгоритма требуется предварительная обработка полигона, которая производится за O(n), но только один раз.
                        0
                        Нас где-то обманули
                        если только чуть-чуть :)
                        0
                        можно ещё сделать «заливку» от заданной точки.
                        и проверить цвет точки, которая заведомо вне многоугольника лежит )
                          0
                          Заливка может производится если есть массив пикселов. Полигон же может быть с любыми сторонами, даже меньше одного пиксела в длине.
                            0
                            хм. в таком разрезе — да.
                            спасибо за поправку
                          +2
                          С меньше или равно и нулём — я бы осторожнее. Если арифметика — не целочисленная, то из-за ошибок округления всё может поехать к чОртовой матери.

                          Юзайте сравнения с эпсилон!
                            0
                            Моя оплошность — не указал явно, что все координаты предполагаются целыми. Для вещественных чисел, конечно, все будет не так шоколадно…
                            +1
                            Выпухлые полигоны — зверь редкий, но очень вкусный.
                            Половина современной хитрой математики определена именно что для них.
                            Представим что у меня есть произвольно заданный полигон без самопересечений(которые вообще можно отдетектить за nlogn)
                            Что же делать?
                            И тут нам на помощь приходят монотоны
                            За O(n) вы находите сегменты линии монотонные по Y( y либо растет либо убывает, в общем не меняет знак производной)
                            Далее можно отработать по старинке — пустить луч из точки и считать пересечения.
                            Но монотоны это отсортированная последовательность, и нужный сегмент для тестирования можно найти за O(logn) — поиском делением пополам.
                            И из всего набора отрезков проверять только нужные.
                            Можно отработать местами еще хитрее, но о них не буду рассказывать пока не прочитаю что там за ссылкой в самом топике скрыто…
                              0
                              Интересно… Однако, вопрос: какие будут монотоны вот для такого многоугольника (и сколько их будет)?
                                +1
                                хренова тут будет.
                                Выкручиваемся двумя способами — либо через теже самые монотоны или BSP бьем на конвексы, и работаем как в топике и сказано.
                                Либо поднимаем «spatial bullets» они же grid — но построение такой структуры относительно дорогое.
                                Вот для данного примера — оно смысла вообще не имеет — быстрее 50 раз проверить в лоб чем построить этот справочник и гонять по чему почти что за O(1)
                                  0
                                  Ага, вот и у меня такое же чувство :) Кажется, что это плохой пример для всех быстрых методов…
                              0
                              Хорошее оформление статьи, очень наглядно.
                                0
                                ЭЭЭ.

                                Достаточно провести любую прямую через точку А.
                                Если кол-во пересечений луча с многоугольником нечетно — то внутри.
                                Если кол-во пересечений луча с многоугольником 0 или четно — то снаружи.
                                  0
                                  Простите, а как вы будете искать точки пересечения, например?
                                    0
                                    А как у вас задан многоугольник?
                                    В виде функции, системы функций? или последовательности координат вершин многоугольника? :)

                                    Для второго случая выполним пересечение двух прямых (одна из них наша прямая, а вторая — образавана из двух вершин многоугольника), в результате, если прямые не параллельны, получим какую-то точку, которую надо проверить на принадлежность обоим отрезкам; для этого достаточно проверить, что эта точка принадлежит обоим отрезкам в проекции на ось X и на ось Y.

                                    Это самый универсальный способ.
                                    Перед ним можно сделать проверки, чтоб отсеять явно не пересекающие стороны многоугольника.
                                      0
                                      O(n) — чего и требуется избежать
                                  0
                                  ->в таких многоугольниках все диагонали, ведущие из нулевой вершины должны лежать внутри многоугольника:
                                  Проще говоря, точки должны быть не только отсортированы по углу, еще и порядок обхода этих точек должен совпадать с этой сортировкой.
                                  Кстати, этот метод преподается(вался) в политехе. Так же, как и всякие интерполяции сплайнами, методы наименьших квадратов и т.п.
                                    +1
                                    Если у нас один (невыпуклый) многоугольник и много запросов принадлежности точки, то задача решается за O(log(N)) на запрос плюс O(N^2*log(N)) на предвычисления (и O(N^2) памяти): режем плоскость на горизонтальные полосы линиями, проходящими через вершины, и для каждой полосы перечисляем стороны многоугольника (слева направо), которые эту полосу пересекают. Задача решается двумя бинарными поисками: по координате y находим полосу, а потом — по (x,y) находим положение точки в этой полосе.
                                      0
                                      bin, кстати надо будет его попробовать.
                                      Только почему предвычисления N^2*log(N)
                                        +1
                                        Я не очень думал над алгоритмом предвычисления. Первое, что пришло в голову — сортируем y-координаты (N*log(N) операций), получаем границы полос. Для каждой полосы (их N+1) ищем точки пересечения «средней линии» со всеми сторонами (общее время всего O(N^2), и в худшем случае — когда наш многоугольник закручен в спираль — получается O(N^2) точек пересечения), а потом сортируем эти точки для каждой полосы (общее время в худшем случае — как раз O(N^2*log(N)).
                                        Вероятно, можно воспользоваться тем, что при переходе от одной полосы к другой последовательность сторон меняется только в одном месте. При использовании простых вставок время получится O(N^2). Интересно, можно ли довести его на O(N*log(N)), и какой должна быть структура представления многоугольника, укладывающаяся в такую память.
                                          0
                                          Если K запросов на принадлежность точки известны заранее, то ответить на них можно за O(K*log(K)+N*log(N)). А заодно получается алгоритм, позволяющий за O(N*log(N)) проверить, есть ли у замкнутой ломаной самопересечения. Память в обоих случаях линейна.
                                      0
                                      Я для каждого многоугольника вычисляю «коробку», т.е. две пары максимальных-минимальных координат. Если искомая точка не лежит между ними, то вычислять многоугольник не имеет смысла.

                                      Only users with full accounts can post comments. Log in, please.