Как определить, находится ли 2D-точка внутри многоугольника?

Я пытаюсь создать быструю двумерную точку внутри алгоритма многоугольника, для использования при хит-тестировании (например, Polygon.contains(p:Point) ). Предложения по эффективным методам будут оценены.

30 Solutions collect form web for “Как определить, находится ли 2D-точка внутри многоугольника?”

Для графики я предпочитаю не использовать целые числа. Многие системы используют целые числа для рисования UI (пикселы – это ints), но macOS, например, использует float для всего. macOS знает только точки и точка может перевести на один пиксель, но в зависимости от разрешения монитора он может перевести что-то еще. На сетчатых экранах половина точки (0,5 / 0,5) – это пиксель. Тем не менее, я никогда не замечал, что пользовательские интерфейсы MacOS значительно медленнее, чем другие пользовательские интерфейсы. Ведь 3D-API (OpenGL или Direct3D) также работает с поплавками, а современные графические библиотеки очень часто используют ускорение GPU.

Теперь вы сказали, что скорость – это ваша главная забота, хорошо, давайте идти на скорость. Прежде чем запускать какой-либо сложный алгоритм, сначала выполните простой тест. Создайте ось с выравниванием по оси вокруг вашего многоугольника. Это очень просто, быстро и уже может вас сэкономить много расчетов. Как это работает? Перейдем по всем точкам многоугольника и найдем значения min / max для X и Y.

Например, у вас есть точки (9/1), (4/3), (2/7), (8/2), (3/6) . Это означает, что Xmin равно 2, Xmax равно 9, Ymin равно 1, а Ymax равно 7. Точка вне прямоугольника с двумя ребрами (2/1) и (9/7) не может находиться внутри многоугольника.

 // p is your point, px is the x coord, py is the y coord if (px < Xmin || px > Xmax || py < Ymin || py > Ymax) { // Definitely not within the polygon! } 

Это первый тест для любой точки. Как вы можете видеть, этот тест очень быстрый, но он также очень грубый. Чтобы обрабатывать точки, находящиеся внутри ограничивающего прямоугольника, нам нужен более сложный алгоритм. Есть несколько способов, как это можно вычислить. Какой метод работает также зависит от того, может ли многоугольник иметь отверстия или всегда будет твердым. Вот примеры твердых (один выпуклый, один вогнутый):

Многоугольник без отверстия

И вот один с отверстием:

Многоугольник с отверстием

У зеленого есть отверстие посередине!

Самый простой алгоритм, который может обрабатывать все три случая выше и все еще довольно быстро, называется лучей . Идея алгоритма довольно проста: нарисуйте виртуальный луч из любой точки вне полигона до своей точки и подсчитайте, как часто он попадает на сторону многоугольника. Если количество обращений равно, это вне полигона, если оно странно, оно внутри.

Демонстрация того, как луч прорезает многоугольник

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

Помните, у вас все еще есть ограничивающая рамка сверху? Просто выберите точку за пределами frameworks и используйте ее в качестве отправной точки для своего луча. Например, точка (Xmin - e/py) находится за пределами многоугольника.

Но что такое e ? Ну, e (фактически epsilon) дает ограничивающее поле некоторую отступу . Как я уже сказал, трассировка лучей терпит неудачу, если мы начнем слишком близко к линии многоугольника. Поскольку ограничивающая рамка может равняться многоугольнику (если многоугольник представляет собой прямоугольник, выровненный по оси, ограничивающий прямоугольник равен самому многоугольнику!), Нам нужно некоторое дополнение, чтобы сделать это безопасным, вот и все. Насколько большой вы должны выбрать e ? Не слишком большой. Это зависит от масштаба системы координат, который вы используете для рисования. Если ваша ширина шага пикселя равна 1,0, то просто выберите 1.0 (еще 0,1 работал бы также)

Теперь, когда у нас есть луч с его начальными и конечными координатами, проблема сдвигается от « находится в пределах полигона » к « как часто луч пересекает сторону многоугольника ». Поэтому мы не можем просто работать с точками полигона по-прежнему, теперь нам нужны фактические стороны. Сторона всегда определяется двумя точками.

 side 1: (X1/Y1)-(X2/Y2) side 2: (X2/Y2)-(X3/Y3) side 3: (X3/Y3)-(X4/Y4) : 

Вам нужно протестировать луч со всех сторон. Рассмотрим луч как вектор, а каждая сторона – вектор. Луч должен ударить каждую сторону ровно один раз или никогда вообще. Он не может ударить по одной и той же стороне дважды. Две строки в 2D-пространстве всегда будут пересекаться ровно один раз, если они не параллельны, и в этом случае они никогда не пересекаются. Однако, поскольку векторы имеют ограниченную длину, два вектора могут быть не параллельными и все же никогда не пересекаться, потому что они слишком коротки, чтобы встречаться друг с другом.

 // Test the ray against all sides int intersections = 0; for (side = 0; side < numberOfSides; side++) { // Test if current side intersects with ray. // If yes, intersections++; } if ((intersections & 1) == 1) { // Inside of polygon } else { // Outside of polygon } 

До сих пор так хорошо, но как вы проверяете, пересекаются ли два вектора? Вот код C (не протестирован), который должен сделать трюк:

 #define NO 0 #define YES 1 #define COLLINEAR 2 int areIntersecting( float v1x1, float v1y1, float v1x2, float v1y2, float v2x1, float v2y1, float v2x2, float v2y2 ) { float d1, d2; float a1, a2, b1, b2, c1, c2; // Convert vector 1 to a line (line 1) of infinite length. // We want the line in linear equation standard form: A*x + B*y + C = 0 // See: http://en.wikipedia.org/wiki/Linear_equation a1 = v1y2 - v1y1; b1 = v1x1 - v1x2; c1 = (v1x2 * v1y1) - (v1x1 * v1y2); // Every point (x,y), that solves the equation above, is on the line, // every point that does not solve it, is not. The equation will have a // positive result if it is on one side of the line and a negative one // if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector // 2 into the equation above. d1 = (a1 * v2x1) + (b1 * v2y1) + c1; d2 = (a1 * v2x2) + (b1 * v2y2) + c1; // If d1 and d2 both have the same sign, they are both on the same side // of our line 1 and in that case no intersection is possible. Careful, // 0 is a special case, that's why we don't test ">=" and "< =", // but "<" and ">". if (d1 > 0 && d2 > 0) return NO; if (d1 < 0 && d2 < 0) return NO; // The fact that vector 2 intersected the infinite line 1 above doesn't // mean it also intersects the vector 1. Vector 1 is only a subset of that // infinite line 1, so it may have intersected that line before the vector // started or after it ended. To know for sure, we have to repeat the // the same test the other way round. We start by calculating the // infinite line 2 in linear equation standard form. a2 = v2y2 - v2y1; b2 = v2x1 - v2x2; c2 = (v2x2 * v2y1) - (v2x1 * v2y2); // Calculate d1 and d2 again, this time using points of vector 1. d1 = (a2 * v1x1) + (b2 * v1y1) + c2; d2 = (a2 * v1x2) + (b2 * v1y2) + c2; // Again, if both have the same sign (and neither one is 0), // no intersection is possible. if (d1 > 0 && d2 > 0) return NO; if (d1 < 0 && d2 < 0) return NO; // If we get here, only two possibilities are left. Either the two // vectors intersect in exactly one point or they are collinear, which // means they intersect in any number of points from zero to infinite. if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR; // If they are not collinear, they must intersect in exactly one point. return YES; } 

v1x1/v1y1 значениями являются две конечные точки вектора 1 ( v1x1/v1y1 и v1x2/v1y2 ) и вектор 2 ( v2x1/v2y1 и v2x2/v2y2 ). Таким образом, у вас есть 2 вектора, 4 точки, 8 координат. YES и NO ясны. YES увеличивает пересечения, NO ничего не делает.

Как насчет COLLINEAR? Это означает, что оба вектора лежат на одной и той же бесконечной линии, в зависимости от положения и длины, они не пересекаются вообще или пересекаются в бесконечном числе точек. Я не совсем уверен, как справиться с этим делом, я бы не стал считать его пересечением в любом случае. Ну, этот случай на практике довольно редок в любом случае из-за ошибок округления с плавающей запятой; лучший код, вероятно, не будет проверять на == 0.0f а вместо этого на что-то вроде < epsilon , где epsilon - довольно небольшое число.

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

И последнее, но не менее важное: если вы можете использовать 3D-оборудование для решения проблемы, есть интересная альтернатива. Просто позвольте GPU делать всю работу за вас. Создайте поверхность для рисования, которая выключена. Заполните его полностью черным цветом. Теперь пусть OpenGL или Direct3D нарисуют ваш полигон (или даже все ваши полигоны, если вы просто хотите проверить, находится ли точка в любом из них, но вам все равно, какой) и заполнить многоугольник (ы) другим цвет, например белый. Чтобы проверить, находится ли точка внутри полигона, получите цвет этой точки с поверхности чертежа. Это всего лишь выборка памяти O (1).

Конечно, этот метод можно использовать только в том случае, если ваша поверхность рисования не должна быть огромной. Если он не может вписаться в память графического процессора, этот метод работает медленнее, чем при работе с ЦП. Если это должно быть огромным, и ваш GPU поддерживает современные шейдеры, вы все равно можете использовать графический процессор, реализуя лучевое кастинг, показанный выше, как графический шейдер, что абсолютно возможно. Для большего количества полигонов или большого количества точек для тестирования это окупится, подумайте, что некоторые графические процессоры смогут протестировать от 64 до 256 точек параллельно. Обратите внимание, однако, что передача данных с CPU на GPU и обратно всегда стоит дорого, поэтому для простого тестирования пары точек против нескольких простых полигонов, где точки или полигоны являются динамическими и будут часто меняться, подход GPU редко будет платить выкл.

Я думаю, что следующий fragment кода – лучшее решение (взято отсюда ):

 int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy) { int i, j, c = 0; for (i = 0, j = nvert-1; i < nvert; j = i++) { if ( ((verty[i]>testy) != (verty[j]>testy)) && (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) ) c = !c; } return c; } 

аргументы

  • nvert : Число вершин в многоугольнике. Повторение первой вершины в конце было обсуждено в статье, упомянутой выше.
  • vertx, verty : Массивы, содержащие x- и y-координаты вершин многоугольника.
  • testx, testy : X- и y-координаты контрольной точки.

Он как короткий, так и эффективный и работает как для выпуклых, так и вогнутых полигонов. Как было предложено ранее, вы должны сначала проверить ограничивающий прямоугольник и обрабатывать многоугольные отверстия отдельно.

Идея этого довольно проста. Автор описывает это следующим образом:

Я запускаю полубесконечный луч горизонтально (увеличивая x, фиксированный y) из тестовой точки и подсчитывая, сколько ребер пересекает. На каждом пересечении луч переключается между внутренней и внешней. Это называется теоремой о кривых Жордана.

Переменная c переключается с 0 на 1 и 1 на 0 каждый раз, когда горизонтальный луч пересекает любое ребро. Поэтому в основном он отслеживает, четное или нечетное число пересекаемых ребер. 0 означает, что четный и 1 означает нечетный.

Вот версия ответа C #, предоставленная nirg , которая исходит от этого профессора RPI . Обратите внимание, что использование кода из этого источника RPI требует атрибуции.

В верхней части была добавлена ​​проверка frameworks. Тем не менее, как отмечает Джеймс Браун, основной код почти так же быстро, как и ограничивающий блок, сам проверяет, поэтому проверка ограничивающего блока может фактически замедлить общую операцию в том случае, если большинство точек, которые вы проверяете, находятся в ограничивающей рамке , Таким образом, вы можете оставить контрольную коробку, или альтернативой было бы прекомпотировать ограничивающие прямоугольники ваших полигонов, если они не меняют форму слишком часто.

 public bool IsPointInPolygon( Point p, Point[] polygon ) { double minX = polygon[ 0 ].X; double maxX = polygon[ 0 ].X; double minY = polygon[ 0 ].Y; double maxY = polygon[ 0 ].Y; for ( int i = 1 ; i < polygon.Length ; i++ ) { Point q = polygon[ i ]; minX = Math.Min( qX, minX ); maxX = Math.Max( qX, maxX ); minY = Math.Min( qY, minY ); maxY = Math.Max( qY, maxY ); } if ( pX < minX || pX > maxX || pY < minY || pY > maxY ) { return false; } // http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html bool inside = false; for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ ) { if ( ( polygon[ i ].Y > pY ) != ( polygon[ j ].Y > pY ) && pX < ( polygon[ j ].X - polygon[ i ].X ) * ( pY - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X ) { inside = !inside; } } return inside; } 

Вот вариант варианта ответа М. Каца на основе подхода Нирга:

 function pointIsInPoly(p, polygon) { var isInside = false; var minX = polygon[0].x, maxX = polygon[0].x; var minY = polygon[0].y, maxY = polygon[0].y; for (var n = 1; n < polygon.length; n++) { var q = polygon[n]; minX = Math.min(qx, minX); maxX = Math.max(qx, maxX); minY = Math.min(qy, minY); maxY = Math.max(qy, maxY); } if (px < minX || px > maxX || py < minY || py > maxY) { return false; } var i = 0, j = polygon.length - 1; for (i, j; i < polygon.length; j = i++) { if ( (polygon[i].y > py) != (polygon[j].y > py) && px < (polygon[j].x - polygon[i].x) * (py - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) { isInside = !isInside; } } return isInside; } 

Вычислить ориентированную сумму углов между точкой p и каждой из вершин многоугольника. Если общий ориентированный угол равен 360 gradleусам, точка внутри. Если сумма равна 0, точка находится вне.

Мне нравится этот метод лучше, потому что он более надежный и менее зависимый от числовой точности.

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

EDIT: By The Way, этот метод работает с вогнутыми и выпуклыми многоугольниками.

EDIT: Недавно я нашел целую статью в Википедии по этой теме.

Статья Эрика Хейнса, процитированная bobobobo, действительно превосходна. Особенно интересны таблицы, сравнивающие производительность алгоритмов; метод суммирования углов действительно плох по сравнению с другими. Интересно также, что оптимизация, например, использование сетки поиска для дальнейшего разбиения полигона на «in» и «out», может сделать тест невероятно быстрым даже на многоугольниках с> 1000 сторонами.

Во всяком случае, это ранние дни, но мой голос переходит к методу «пересечений», что в значительной степени отражает мнение Мекки. Однако я нашел его наиболее сугубо описанным и кодифицированным Дэвидом Бурком . Мне нравится, что нет никакой реальной тригонометрии, и она работает на выпуклую и вогнутую, и она работает достаточно хорошо, как увеличивается количество сторон.

Кстати, вот одна из таблиц производительности из статьи Эрика Хейнса для интереса, тестирование на случайных полигонах.

  number of edges per polygon 3 4 10 100 1000 MacMartin 2.9 3.2 5.9 50.6 485 Crossings 3.1 3.4 6.8 60.0 624 Triangle Fan+edge sort 1.1 1.8 6.5 77.6 787 Triangle Fan 1.2 2.1 7.3 85.4 865 Barycentric 2.1 3.8 13.8 160.7 1665 Angle Summation 56.2 70.4 153.6 1403.8 14693 Grid (100x100) 1.5 1.5 1.6 2.1 9.8 Grid (20x20) 1.7 1.7 1.9 5.7 42.2 Bins (100) 1.8 1.9 2.7 15.1 117 Bins (20) 2.1 2.2 3.7 26.3 278 

Этот вопрос так интересен. У меня есть другая работающая идея, отличная от других ответов этого сообщения.

Идея состоит в том, чтобы использовать сумму углов, чтобы решить, находится ли цель внутри или снаружи. Если цель находится внутри области, сумма формы угла по цели и каждые две граничные точки будет равна 360. Если цель находится снаружи, сумма не будет равна 360. Угол имеет направление. Если угол идет назад, угол является отрицательным. Это точно так же, как вычисление числа обмотки .

Обратитесь к этому изображению, чтобы получить базовое представление об этой идее: введите описание изображения здесь

Мой алгоритм предполагает, что по часовой стрелке есть положительное направление. Вот потенциальный ввод:

 [[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]] 

Ниже приведен код python, который реализует идею:

 def isInside(self, border, target): degree = 0 for i in range(len(border) - 1): a = border[i] b = border[i + 1] # calculate distance of vector A = getDistance(a[0], a[1], b[0], b[1]); B = getDistance(target[0], target[1], a[0], a[1]) C = getDistance(target[0], target[1], b[0], b[1]) # calculate direction of vector ta_x = a[0] - target[0] ta_y = a[1] - target[1] tb_x = b[0] - target[0] tb_y = b[1] - target[1] cross = tb_y * ta_x - tb_x * ta_y clockwise = cross < 0 # calculate sum of angles if(clockwise): degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C))) else: degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C))) if(abs(round(degree) - 360) <= 3): return True return False по def isInside(self, border, target): degree = 0 for i in range(len(border) - 1): a = border[i] b = border[i + 1] # calculate distance of vector A = getDistance(a[0], a[1], b[0], b[1]); B = getDistance(target[0], target[1], a[0], a[1]) C = getDistance(target[0], target[1], b[0], b[1]) # calculate direction of vector ta_x = a[0] - target[0] ta_y = a[1] - target[1] tb_x = b[0] - target[0] tb_y = b[1] - target[1] cross = tb_y * ta_x - tb_x * ta_y clockwise = cross < 0 # calculate sum of angles if(clockwise): degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C))) else: degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C))) if(abs(round(degree) - 360) <= 3): return True return False 

Я немного поработал над этим, когда я был исследователем под руководством Майкла Стоунбрейкера – вы знаете, профессор, который придумал Ingres , PostgreSQL и т. Д.

Мы поняли, что самым быстрым способом было сначала сделать ограничительную рамку, потому что это СУПЕР быстро. Если он находится за пределами frameworks, он находится снаружи. В противном случае вы делаете более трудную работу …

Если вам нужен отличный алгоритм, посмотрите на исходный код PostgreSQL с открытым исходным кодом для геолокации …

Я хочу отметить, что у нас никогда не было никакого представления о правильном и левостороннем (также выражаемом как проблема «внутри» и «снаружи» …


ОБНОВИТЬ

Ссылка BKB обеспечила большое количество разумных алгоритмов. Я работал над проблемами науки о Земле и поэтому нуждался в решении, которое работает в широте / долготе, и у него есть своеобразная проблема ручности – это область внутри меньшей площади или большей площади? Ответ заключается в том, что «направление» вершин имеет значение – оно либо левое, либо правое, и таким образом вы можете указать любую область как «внутри» любого заданного полигона. Таким образом, моя работа использовала решение, перечисленное на этой странице.

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

… Так как кто-то спросил: мы выяснили, что тесты с ограничивающей коробкой лучше всего, когда количество вершин выходит за frameworks некоторого числа – сделайте очень быстрый тест, прежде чем выполнять более длительный тест, если это необходимо … Ограничивающий прямоугольник создается, просто принимая наибольший x, наименьший x, наибольший y и наименьший y и складывающий их вместе, чтобы сделать четыре точки коробки …

Еще один совет для следующих: мы выполнили все наши более сложные и «световые тускнеющие» вычисления в пространстве сетки все в положительных точках на плоскости, а затем снова проецировались обратно в «настоящую» долготу / широту, что позволило избежать возможных ошибок обертывание, когда одна перекрестная линия 180 долготы и при обработке полярных областей. Отлично!

На самом деле нравится решение, отправленное Nirg и отредактированное bobobobo. Я просто сделал его javascript дружественным и немного более parsingчивым для моего использования:

 function insidePoly(poly, pointx, pointy) { var i, j; var inside = false; for (i = 0, j = poly.length - 1; i < poly.length; j = i++) { if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside; } return inside; } 

Быстрая версия ответа по nirg :

 extension CGPoint { func isInsidePolygon(vertices: [CGPoint]) -> Bool { guard !vertices.isEmpty else { return false } var j = vertices.last!, c = false for i in vertices { let a = (iy > y) != (jy > y) let b = (x < (jx - ix) * (y - iy) / (jy - iy) + ix) if a && b { c = !c } j = i } return c } } 

Ответ Дэвида Сегонда в значительной степени является стандартным общим ответом, а Ричард Т – наиболее распространенная оптимизация, хотя некоторые другие. Другие сильные оптимизации основаны на менее общих решениях. Например, если вы собираетесь проверять один и тот же полигон с большим количеством точек, триангуляция полигона может ускорить работу, так как существует множество очень быстрых алгоритмов поиска TIN. Другим является то, что многоугольник и точки находятся на ограниченной плоскости при низком разрешении, например, на экране, вы можете рисовать многоугольник в буфер отображения отображаемой памяти в заданном цвете и проверять цвет данного пикселя, чтобы увидеть, в многоугольниках.

Как и многие оптимизации, они основаны на конкретных, а не общих случаях, и дают преимущества, основанные на амортизированном времени, а не на одном использовании.

Работая в этом поле, я нашел полезную геометрию вычислений Джозефа О’Руркса в C ‘ISBN 0-521-44034-3.

Тривиальное решение состояло бы в том, чтобы разделить многоугольник на треугольники и проверить тестовые треугольники, как описано здесь

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

Я понимаю, что это старо, но вот алгоритм кастинга лучей, реализованный в Cocoa, в случае, если кто-то заинтересован. Не уверен, что это самый эффективный способ сделать что-то, но это может помочь кому-то.

 - (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point { NSBezierPath *currentPath = [path bezierPathByFlatteningPath]; BOOL result; float aggregateX = 0; //I use these to calculate the centroid of the shape float aggregateY = 0; NSPoint firstPoint[1]; [currentPath elementAtIndex:0 associatedPoints:firstPoint]; float olderX = firstPoint[0].x; float olderY = firstPoint[0].y; NSPoint interPoint; int noOfIntersections = 0; for (int n = 0; n < [currentPath elementCount]; n++) { NSPoint points[1]; [currentPath elementAtIndex:n associatedPoints:points]; aggregateX += points[0].x; aggregateY += points[0].y; } for (int n = 0; n < [currentPath elementCount]; n++) { NSPoint points[1]; [currentPath elementAtIndex:n associatedPoints:points]; //line equations in Ax + By = C form float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y; float _B_FOO = point.x - (aggregateX/[currentPath elementCount]); float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y); float _A_BAR = olderY - points[0].y; float _B_BAR = points[0].x - olderX; float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY); float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO); if (det != 0) { //intersection points with the edges float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det; float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det; interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint); if (olderX <= points[0].x) { //doesn't matter in which direction the ray goes, so I send it right-ward. if ((interPoint.x >= olderX && interPoint.x < = points[0].x) && (interPoint.x > point.x)) { noOfIntersections++; } } else { if ((interPoint.x >= points[0].x && interPoint.x < = olderX) && (interPoint.x > point.x)) { noOfIntersections++; } } } olderX = points[0].x; olderY = points[0].y; } if (noOfIntersections % 2 == 0) { result = FALSE; } else { result = TRUE; } return result; } 

Obj-C версия ответа nirg с образцом метода для тестирования точек. Ответ Нирга сработал хорошо для меня.

 - (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test { NSUInteger nvert = [vertices count]; NSInteger i, j, c = 0; CGPoint verti, vertj; for (i = 0, j = nvert-1; i < nvert; j = i++) { verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue]; vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue]; if (( (verti.y > test.y) != (vertj.y > test.y) ) && ( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) ) c = !c; } return (c ? YES : NO); } - (void)testPoint { NSArray *polygonVertices = [NSArray arrayWithObjects: [NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)], [NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)], [NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)], [NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)], [NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)], [NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)], nil ]; CGPoint tappedPoint = CGPointMake(23.0, 70.0); if ([self isPointInPolygon:polygonVertices point:tappedPoint]) { NSLog(@"YES"); } else { NSLog(@"NO"); } } 

sample polygon

I too thought 360 summing only worked for convex polygons but this isn’t true.

This site has a nice diagram showing exactly this, and a good explanation on hit testing: Gamasutra – Crashing into the New Year: Collision Detection

C# version of nirg’s answer is here: I’ll just share the code. It may save someone some time.

 public static bool IsPointInPolygon(IList polygon, Point testPoint) { bool result = false; int j = polygon.Count() - 1; for (int i = 0; i < polygon.Count(); i++) { if (polygon[i].Y < testPoint.Y && polygon[j].Y >= testPoint.Y || polygon[j].Y < testPoint.Y && polygon[i].Y >= testPoint.Y) { if (polygon[i].X + (testPoint.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) * (polygon[j].X - polygon[i].X) < testPoint.X) { result = !result; } } j = i; } return result; } 

Java Version:

 public class Geocode { private float latitude; private float longitude; public Geocode() { } public Geocode(float latitude, float longitude) { this.latitude = latitude; this.longitude = longitude; } public float getLatitude() { return latitude; } public void setLatitude(float latitude) { this.latitude = latitude; } public float getLongitude() { return longitude; } public void setLongitude(float longitude) { this.longitude = longitude; } } public class GeoPolygon { private ArrayList points; public GeoPolygon() { this.points = new ArrayList(); } public GeoPolygon(ArrayList points) { this.points = points; } public GeoPolygon add(Geocode geo) { points.add(geo); return this; } public boolean inside(Geocode geo) { int i, j; boolean c = false; for (i = 0, j = points.size() - 1; i < points.size(); j = i++) { if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) && (geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude())) c = !c; } return c; } } 

.Net port:

  static void Main(string[] args) { Console.Write("Hola"); List vertx = new List(); List verty = new List(); int i, j, c = 0; vertx.Add(1); vertx.Add(2); vertx.Add(1); vertx.Add(4); vertx.Add(4); vertx.Add(1); verty.Add(1); verty.Add(2); verty.Add(4); verty.Add(4); verty.Add(1); verty.Add(1); int nvert = 6; //Vértices del poligono double testx = 2; double testy = 5; for (i = 0, j = nvert - 1; i < nvert; j = i++) { if (((verty[i] > testy) != (verty[j] > testy)) && (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i])) c = 1; } } 

There is nothing more beutiful than an inductive definition of a problem. For the sake of completeness here you have a version in prolog which might also clarify the thoughs behind ray casting :

Based on the simulation of simplicity algorithm in http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

Some helper predicates:

 exor(A,B):- \+A,B;A,\+B. in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)). inside(false). inside(_,[_|[]]). inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) + X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]). get_line(_,_,[]). get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]). 

The equation of a line given 2 points A and B (Line(A,B)) is:

  (YB-YA) Y - YA = ------- * (X - XA) (XB-YB) 

It is important that the direction of rotation for the line is setted to clock-wise for boundaries and anti-clock-wise for holes. We are going to check whether the point (X,Y), ie the tested point is at the left half-plane of our line (it is a matter of taste, it could also be the right side, but also the direction of boundaries lines has to be changed in that case), this is to project the ray from the point to the right (or left) and acknowledge the intersection with the line. We have chosen to project the ray in the horizontal direction (again it is a matter of taste, it could also be done in vertical with similar restrictions), so we have:

  (XB-XA) X < ------- * (Y - YA) + XA (YB-YA) 

Now we need to know if the point is at the left (or right) side of the line segment only, not the entire plane, so we need to restrict the search only to this segment, but this is easy since to be inside the segment only one point in the line can be higher than Y in the vertical axis. As this is a stronger restriction it needs to be the first to check, so we take first only those lines meeting this requirement and then check its possition. By the Jordan Curve theorem any ray projected to a polygon must intersect at an even number of lines. So we are done, we will throw the ray to the right and then everytime it intersects a line, toggle its state. However in our implementation we are goint to check the lenght of the bag of solutions meeting the given restrictions and decide the innership upon it. for each line in the polygon this have to be done.

 is_left_half_plane(_,[],[],_). is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] = [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA)); is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test). in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon), in_range(Y,YA,YB). all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line), in_y_range_at_poly(Coordinate,Line,Polygon), Lines). traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, < ), IntersectingLines), length(IntersectingLines, Count). % This is the entry point predicate inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside). 

VBA VERSION:

Note: Remember that if your polygon is an area within a map that Latitude/Longitude are Y/X values as opposed to X/Y (Latitude = Y, Longitude = X) due to from what I understand are historical implications from way back when Longitude was not a measurement.

CLASS MODULE: CPoint

 Private pXValue As Double Private pYValue As Double '''''X Value Property''''' Public Property Get X() As Double X = pXValue End Property Public Property Let X(Value As Double) pXValue = Value End Property '''''Y Value Property''''' Public Property Get Y() As Double Y = pYValue End Property Public Property Let Y(Value As Double) pYValue = Value End Property 

MODULE:

 Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean Dim i As Integer Dim j As Integer Dim q As Object Dim minX As Double Dim maxX As Double Dim minY As Double Dim maxY As Double minX = polygon(0).X maxX = polygon(0).X minY = polygon(0).Y maxY = polygon(0).Y For i = 1 To UBound(polygon) Set q = polygon(i) minX = vbMin(qX, minX) maxX = vbMax(qX, maxX) minY = vbMin(qY, minY) maxY = vbMax(qY, maxY) Next i If pX < minX Or pX > maxX Or pY < minY Or pY > maxY Then isPointInPolygon = False Exit Function End If ' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html isPointInPolygon = False i = 0 j = UBound(polygon) Do While i < UBound(polygon) + 1 If (polygon(i).Y > pY) Then If (polygon(j).Y < pY) Then If pX < (polygon(j).X - polygon(i).X) * (pY - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then isPointInPolygon = True Exit Function End If End If ElseIf (polygon(i).Y < pY) Then If (polygon(j).Y > pY) Then If pX < (polygon(j).X - polygon(i).X) * (pY - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then isPointInPolygon = True Exit Function End If End If End If j = i i = i + 1 Loop End Function Function vbMax(n1, n2) As Double vbMax = IIf(n1 > n2, n1, n2) End Function Function vbMin(n1, n2) As Double vbMin = IIf(n1 > n2, n2, n1) End Function Sub TestPointInPolygon() Dim i As Integer Dim InPolygon As Boolean ' MARKER Object Dim p As CPoint Set p = New CPoint pX =  pY =  ' POLYGON OBJECT Dim polygon() As CPoint ReDim polygon() 'Amount of vertices in polygon - 1 For i = 0 To  'Same value as above Set polygon(i) = New CPoint polygon(i).X =  'Source a list of values that can be looped through polgyon(i).Y =  'Source a list of values that can be looped through Next i InPolygon = isPointInPolygon(p, polygon) MsgBox InPolygon End Sub 

Heres a point in polygon test in C that isn’t using ray-casting. And it can work for overlapping areas (self intersections), see the use_holes argument.

 /* math lib (defined below) */ static float dot_v2v2(const float a[2], const float b[2]); static float angle_signed_v2v2(const float v1[2], const float v2[2]); static void copy_v2_v2(float r[2], const float a[2]); /* intersection function */ bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr, const bool use_holes) { /* we do the angle rule, define that all added angles should be about zero or (2 * PI) */ float angletot = 0.0; float fp1[2], fp2[2]; unsigned int i; const float *p1, *p2; p1 = verts[nr - 1]; /* first vector */ fp1[0] = p1[0] - pt[0]; fp1[1] = p1[1] - pt[1]; for (i = 0; i < nr; i++) { p2 = verts[i]; /* second vector */ fp2[0] = p2[0] - pt[0]; fp2[1] = p2[1] - pt[1]; /* dot and angle and cross */ angletot += angle_signed_v2v2(fp1, fp2); /* circulate */ copy_v2_v2(fp1, fp2); p1 = p2; } angletot = fabsf(angletot); if (use_holes) { const float nested = floorf((angletot / (float)(M_PI * 2.0)) + 0.00001f); angletot -= nested * (float)(M_PI * 2.0); return (angletot > 4.0f) != ((int)nested % 2); } else { return (angletot > 4.0f); } } /* math lib */ static float dot_v2v2(const float a[2], const float b[2]) { return a[0] * b[0] + a[1] * b[1]; } static float angle_signed_v2v2(const float v1[2], const float v2[2]) { const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]); return atan2f(perp_dot, dot_v2v2(v1, v2)); } static void copy_v2_v2(float r[2], const float a[2]) { r[0] = a[0]; r[1] = a[1]; } 

Note: this is one of the less optimal methods since it includes a lot of calls to atan2f , but it may be of interest to developers reading this thread (in my tests its ~23x slower then using the line intersection method).

For Detecting hit on Polygon we need to test two things:

  1. If Point is inside polygon area. (can be accomplished by Ray-Casting Algorithm)
  2. If Point is on the polygon border(can be accomplished by same algorithm which is used for point detection on polyline(line)).

To deal with the following special cases in Ray casting algorithm :

  1. The ray overlaps one of the polygon’s side.
  2. The point is inside of the polygon and the ray passes through a vertex of the polygon.
  3. The point is outside of the polygon and the ray just touches one of the polygon’s angle.

Check Determining Whether A Point Is Inside A Complex Polygon . The article provides an easy way to resolve them so there will be no special treatment required for the above cases.

You can do this by checking if the area formed by connecting the desired point to the vertices of your polygon matches the area of the polygon itself.

Or you could check if the sum of the inner angles from your point to each pair of two consecutive polygon vertices to your check point sums to 360, but I have the feeling that the first option is quicker because it doesn’t involve divisions nor calculations of inverse of trigonometric functions.

I don’t know what happens if your polygon has a hole inside it but it seems to me that the main idea can be adapted to this situation

You can as well post the question in a math community. I bet they have one million ways of doing that

If you are looking for a java-script library there’s a javascript google maps v3 extension for the Polygon class to detect whether or not a point resides within it.

 var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3); var isWithinPolygon = polygon.containsLatLng(40, -90); 

Google Extention Github

When using qt (Qt 4.3+), one can use QPolygon’s function containsPoint

The answer depends on if you have the simple or complex polygons. Simple polygons must not have any line segment intersections. So they can have the holes but lines can’t cross each other. Complex regions can have the line intersections – so they can have the overlapping regions, or regions that touch each other just by a single point.

For simple polygons the best algorithm is Ray casting (Crossing number) algorithm. For complex polygons, this algorithm doesn’t detect points that are inside the overlapping regions. So for complex polygons you have to use Winding number algorithm.

Here is an excellent article with C implementation of both algorithms. I tried them and they work well.

http://geomalgorithms.com/a03-_inclusion.html

Scala version of solution by nirg (assumes bounding rectangle pre-check is done separately):

 def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = { val length = polygon.length @tailrec def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = { if (i == length) tracker else { val intersects = (polygon(i).y > py) != (polygon(j).y > py) && px < (polygon(j).x - polygon(i).x) * (py - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x oddIntersections(i + 1, i, if (intersects) !tracker else tracker) } } oddIntersections(0, length - 1, tracker = false) } 

I’ve made a Python implementation of nirg’s c++ code :

Inputs

  • bounding_points: nodes that make up the polygon.
  • bounding_box_positions: candidate points to filter. (In my implementation created from the bounding box.

    (The inputs are lists of tuples in the format: [(xcord, ycord), ...] )

Возвращает

  • All the points that are inside the polygon.
 def polygon_ray_casting(self, bounding_points, bounding_box_positions): # Arrays containing the x- and y-coordinates of the polygon's vertices. vertx = [point[0] for point in bounding_points] verty = [point[1] for point in bounding_points] # Number of vertices in the polygon nvert = len(bounding_points) # Points that are inside points_inside = [] # For every candidate position within the bounding box for idx, pos in enumerate(bounding_box_positions): testx, testy = (pos[0], pos[1]) c = 0 for i in range(0, nvert): j = i - 1 if i != 0 else nvert - 1 if( ((verty[i] > testy ) != (verty[j] > testy)) and (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ): c += 1 # If odd, that means that we are inside the polygon if c % 2 == 1: points_inside.append(pos) return points_inside 

Again, the idea is taken from here

This only works for convex shapes, but Minkowski Portal Refinement, and GJK are also great options for testing if a point is in a polygon. You use minkowski subtraction to subtract the point from the polygon, then run those algorithms to see if the polygon contains the origin.

Also, interestingly, you can describe your shapes a bit more implicitly using support functions which take a direction vector as input and spit out the farthest point along that vector. This allows you to describe any convex shape.. curved, made out of polygons, or mixed. You can also do operations to combine the results of simple support functions to make more complex shapes.

More info: http://xenocollide.snethen.com/mpr2d.html

Also, game programming gems 7 talks about how to do this in 3d (:

Давайте будем гением компьютера.