Как определить, находится ли 2D-точка внутри многоугольника?
Я пытаюсь создать быструю 2D точку внутри алгоритма многоугольника, для использования в тестировании попаданий (например, Polygon.contains(p:Point)). Предложения для эффективных методов будут оценены.
43 ответа
Для графики я бы предпочел не предпочитать целые числа. Многие системы используют целые числа для рисования пользовательского интерфейса (в конце концов, пиксели - это целые числа), но, например, macOS использует float для всего. macOS знает только точки, и точка может переводиться в один пиксель, но в зависимости от разрешения монитора она может переводиться во что-то другое. На экранах сетчатки половина точки (0,5/0,5) составляет пиксель. Тем не менее, я никогда не замечал, что пользовательские интерфейсы macOS значительно медленнее, чем другие пользовательские интерфейсы. Ведь 3D API (OpenGL или Direct3D) также работают с плавающей точкой, а современные графические библиотеки очень часто используют преимущества ускорения GPU.
Теперь вы сказали, что скорость - ваша главная задача, хорошо, давайте перейдем к скорости. Прежде чем запускать какой-либо сложный алгоритм, сначала выполните простой тест. Создайте выровненный по оси ограничивающий прямоугольник вокруг вашего многоугольника. Это очень легко, быстро и уже может обезопасить вас от многих вычислений. Как это работает? Переберите все точки многоугольника и найдите минимальные / максимальные значения 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, p.x is the x coord, p.y is the y coord
if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
// Definitely not within the polygon!
}
Это первый тест для любой точки. Как видите, этот тест очень быстрый, но он также очень грубый. Для обработки точек, которые находятся внутри ограничивающего прямоугольника, нам нужен более сложный алгоритм. Есть несколько способов, как это можно рассчитать. Какой метод работает, зависит также от того, может ли многоугольник иметь отверстия или всегда будет сплошным. Вот примеры твердых (одна выпуклая, одна вогнутая):
А вот один с дыркой:
Зеленый имеет отверстие в середине!
Самый простой алгоритм, который может обрабатывать все три приведенных выше случая и который все еще довольно быстр, называется приведением лучей. Идея алгоритма довольно проста: нарисуйте виртуальный луч из любой точки за пределами многоугольника к своей точке и посчитайте, как часто он попадает на сторону многоугольника. Если число совпадений четное, оно находится за пределами многоугольника, если оно нечетное, оно внутри.
Алгоритм числа обмоток был бы альтернативой, он более точен для точек, находящихся очень близко к линии многоугольника, но он также намного медленнее. Приведение лучей может быть неудачным для точек, расположенных слишком близко к стороне многоугольника, из-за ограниченной точности с плавающей точкой и проблем с округлением, но в действительности это вряд ли проблема, так как если точка находится близко к стороне, это часто визуально даже невозможно для зритель распознает, находится ли он уже внутри или еще снаружи.
У вас все еще есть ограничительная рамка выше, помните? Просто выберите точку за пределами ограничительной рамки и используйте ее в качестве отправной точки для вашего луча. Например, точка (Xmin - e/p.y) вне полигона точно.
Но что это 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;
}
Входными значениями являются две конечные точки вектора 1 (v1x1/v1y1 а также v1x2/v1y2) и вектор 2 (v2x1/v2y1 а также v2x2/v2y2). Итак, у вас есть 2 вектора, 4 точки, 8 координат. YES а также NO ясно. YES увеличивает пересечения, NO ничего не делает.
Как насчет Коллинеар? Это означает, что оба вектора лежат на одной бесконечной линии, в зависимости от положения и длины, они вообще не пересекаются или пересекаются в бесконечном количестве точек. Я не совсем уверен, как справиться с этим делом, я бы не посчитал это пересечением в любом случае. Ну, во всяком случае, этот случай довольно редко встречается на практике из-за ошибок округления с плавающей точкой; лучший код, вероятно, не будет проверять == 0.0f но вместо этого для чего-то вроде < epsilon где эпсилон - довольно небольшое число.
Если вам нужно протестировать большее количество точек, вы, безусловно, можете немного ускорить все это, сохранив в памяти стандартные формы сторон многоугольника в памяти, поэтому вам не придется каждый раз пересчитывать их. Это сэкономит вам два умножения с плавающей запятой и три вычитания с плавающей запятой при каждом тесте в обмен на сохранение в памяти трех значений с плавающей запятой на каждую сторону многоугольника. Это типичное соотношение памяти и времени вычислений.
И последнее, но не менее важное: если вы можете использовать 3D-оборудование для решения проблемы, есть интересная альтернатива. Просто пусть GPU сделает всю работу за вас. Создайте поверхность для рисования вне экрана. Заполните его полностью черным цветом. Теперь позвольте OpenGL или Direct3D нарисовать ваш многоугольник (или даже все ваши многоугольники, если вы просто хотите проверить, находится ли точка в каком-либо из них, но вам не важно, какая из них), и заполните многоугольник (и) другим цвет, например белый. Чтобы проверить, находится ли точка внутри многоугольника, получите цвет этой точки с поверхности рисования. Это всего лишь O(1) выборка из памяти.
Конечно, этот метод можно использовать только в том случае, если ваша поверхность для рисования не обязательно должна быть огромной. Если он не может поместиться в память графического процессора, этот метод медленнее, чем в CPU. Если он должен быть огромным, а ваш графический процессор поддерживает современные шейдеры, вы все равно можете использовать графический процессор, реализовав приведенное выше приведение лучей в качестве графического шейдера, что абсолютно возможно. Для большего количества полигонов или большого количества точек для тестирования это окупится, учитывая, что некоторые графические процессоры смогут тестировать от 64 до 256 точек параллельно. Тем не менее, обратите внимание, что передача данных из CPU в GPU и обратно всегда обходится дорого, поэтому для простого тестирования пары точек на пару простых многоугольников, где либо точки, либо многоугольники являются динамическими и будут часто меняться, подход с использованием графического процессора редко платит выкл.
Я думаю, что следующий фрагмент кода - лучшее решение (взято отсюда):
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 означает нечетное.
Вот версия ответа Nirg на C# от этого профессора RPI. Обратите внимание, что использование кода из этого источника RPI требует указания авторства.
Ограничительный флажок был добавлен вверху. Однако, как отмечает Джеймс Браун, основной код работает почти так же быстро, как и сам флажок ограничивающего прямоугольника, поэтому флажок ограничивающего прямоугольника может на самом деле замедлить всю операцию в случае, если большинство проверяемых точек находятся внутри ограничивающего прямоугольника., Таким образом, вы можете оставить ограничивающий флажок, или альтернативой будет предварительный расчет ограничивающих прямоугольников ваших полигонов, если они не слишком часто меняют форму.
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( q.X, minX );
maxX = Math.Max( q.X, maxX );
minY = Math.Min( q.Y, minY );
maxY = Math.Max( q.Y, maxY );
}
if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > 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 > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
{
inside = !inside;
}
}
return inside;
}
Вот вариант ответа М. Каца на JavaScript, основанный на подходе Нирга:
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(q.x, minX);
maxX = Math.max(q.x, maxX);
minY = Math.min(q.y, minY);
maxY = Math.max(q.y, maxY);
}
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
return false;
}
var i = 0, j = polygon.length - 1;
for (i, j; i < polygon.length; j = i++) {
if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
isInside = !isInside;
}
}
return isInside;
}
Вычислить ориентированную сумму углов между точкой p и каждым из вершин многоугольника. Если общий угол ориентирования составляет 360 градусов, точка находится внутри. Если сумма равна 0, точка находится за пределами.
Мне больше нравится этот метод, потому что он более надежный и менее зависит от точности чисел.
Методы, которые вычисляют равномерность количества пересечений, ограничены, потому что вы можете "поразить" вершину во время вычисления количества пересечений.
РЕДАКТИРОВАТЬ: Кстати, этот метод работает с вогнутыми и выпуклыми полигонами.
РЕДАКТИРОВАТЬ: Я недавно нашел целую статью в Википедии по этой теме.
Этот вопрос так интересен. У меня есть другая работоспособная идея, отличная от других ответов этого поста.
Идея состоит в том, чтобы использовать сумму углов, чтобы решить, находится ли цель внутри или снаружи. Если цель находится внутри области, сумма угла, образованного целью и каждыми двумя граничными точками, будет 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
Статья Эрика Хейнса, на которую ссылается Бобобо, действительно превосходна. Особенно интересны таблицы, сравнивающие производительность алгоритмов; метод суммирования углов действительно плох по сравнению с другими. Также интересно то, что оптимизация, например использование поисковой сетки для дальнейшего разделения полигона на секторы "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
Большинство ответов на этот вопрос не подходят для всех угловых случаев. Некоторые тонкие угловые случаи, как показано ниже:
Это версия javascript с хорошо обработанными угловыми случаями.
/** Get relationship between a point and a polygon using ray-casting algorithm
* @param {{x:number, y:number}} P: point to check
* @param {{x:number, y:number}[]} polygon: the polygon
* @returns -1: outside, 0: on edge, 1: inside
*/
function relationPP(P, polygon) {
const between = (p, a, b) => p >= a && p <= b || p <= a && p >= b
let inside = false
for (let i = polygon.length-1, j = 0; j < polygon.length; i = j, j++) {
const A = polygon[i]
const B = polygon[j]
// corner cases
if (P.x == A.x && P.y == A.y || P.x == B.x && P.y == B.y) return 0
if (A.y == B.y && P.y == A.y && between(P.x, A.x, B.x)) return 0
if (between(P.y, A.y, B.y)) { // if P inside the vertical range
// filter out "ray pass vertex" problem by treating the line a little lower
if (P.y == A.y && B.y >= A.y || P.y == B.y && A.y >= B.y) continue
// calc cross product `PA X PB`, P lays on left side of AB if c > 0
const c = (A.x - P.x) * (B.y - P.y) - (B.x - P.x) * (A.y - P.y)
if (c == 0) return 0
if ((A.y < B.y) == (c > 0)) inside = !inside
}
}
return inside? 1 : -1
}
Очень нравится решение, опубликованное Nirg и отредактированное bobobobo. Я просто сделал его дружественным к JavaScript и немного более разборчивым для моего использования:
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 = (i.y > y) != (j.y > y)
let b = (x < (j.x - i.x) * (y - i.y) / (j.y - i.y) + i.x)
if a && b { c = !c }
j = i
}
return c
}
}
Я работал над этим, когда был исследователем при Майкле Стоунбрейкере - вы знаете, профессор, который придумал Ingres, PostgreSQL и т. Д.
Мы поняли, что самым быстрым способом было сначала создать ограничивающую рамку, потому что это СУПЕР быстро. Если это вне ограничительной рамки, то это снаружи. В противном случае, вы делаете тяжелую работу...
Если вы хотите отличный алгоритм, обратитесь к исходному коду PostgreSQL с открытым исходным кодом для гео работы...
Я хочу отметить, что мы никогда не понимали право против левши (также выражается как проблема "внутри" и "снаружи"...
ОБНОВИТЬ
Ссылка БКБ предоставила большое количество разумных алгоритмов. Я работал над проблемами наук о Земле и поэтому нуждался в решении, которое работает по широте / долготе, и у него есть особая проблема с ручностью - область внутри меньшей области или большей области? Ответ заключается в том, что "направление" вершин имеет значение - оно либо левостороннее, либо правостороннее, и таким образом вы можете указать любую область как "внутри" любого данного многоугольника. Таким образом, моя работа использовала решение три, перечисленное на этой странице.
Кроме того, в моей работе использовались отдельные функции для тестов "на линии".
... Так как кто-то спросил: мы выяснили, что тесты ограничивающего прямоугольника были наилучшими, когда число вершин превысило некоторое число - проведите очень быстрый тест, прежде чем выполнять более длинный тест, если это необходимо... Ограничивающий прямоугольник создается простым взятием наибольший x, наименьший x, наибольший y и наименьший y и объединение их в четыре точки в квадрате...
Еще один совет для тех, кто следует: мы выполнили все наши более сложные и "светорегуляторные" вычисления в пространстве сетки, все в положительных точках на плоскости, а затем повторно проецировали обратно в "реальную" долготу / широту, избегая, таким образом, возможных ошибок обтекание, когда одна пересекает линию 180 долготы и при работе с полярными областями. Работал отлично!
Тривиальным решением было бы разделить многоугольник на треугольники и нажать "проверить треугольники", как описано здесь
Если у вас многоугольник CONVEX, то может быть лучше. Посмотрите на многоугольник как на коллекцию бесконечных линий. Каждая строка делит пространство на две части. для каждой точки легко сказать, находится ли она на одной стороне или другой стороне линии. Если точка находится на одной стороне всех линий, то она находится внутри многоугольника.
Версия Java:
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<Geocode> points;
public GeoPolygon() {
this.points = new ArrayList<Geocode>();
}
public GeoPolygon(ArrayList<Geocode> 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;
}
}
Ответ Дэвида Сегонда в значительной степени является стандартным общим ответом, а Ричард Т - наиболее распространенная оптимизация, хотя есть и другие. Другие сильные оптимизации основаны на менее общих решениях. Например, если вы собираетесь проверять один и тот же многоугольник с множеством точек, триангуляция многоугольника может значительно ускорить процесс, поскольку существует ряд очень быстрых алгоритмов поиска TIN. Другой вариант: если многоугольник и точки находятся на ограниченной плоскости при низком разрешении, например, на экране, вы можете нарисовать многоугольник в отображаемом в память буфере отображения заданного цвета и проверить цвет данного пикселя, чтобы увидеть, лежит ли он в многоугольниках.
Как и многие оптимизации, они основаны на конкретных, а не общих случаях, и приносят выгоды, основанные на амортизированном времени, а не на единичном использовании.
Работая в этой области, я обнаружил, что вычислительная геометрия Джозефа О'Руркса в C' ISBN 0-521-44034-3 очень помогает.
Я сделал реализацию Python nirg C++ nirg:
входные
- bounding_points: узлы, которые составляют многоугольник.
bounding_box_positions: точки-кандидаты для фильтрации. (В моей реализации создано из ограничительной рамки.
(Входные данные представляют собой списки кортежей в формате:
[(xcord, ycord), ...])
Возвращает
- Все точки, которые находятся внутри многоугольника.
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
Опять же, идея взята отсюда
Я понимаю, что это старый, но вот алгоритм наложения лучей, реализованный в Какао, на случай, если кому-то будет интересно. Не уверен, что это самый эффективный способ сделать что-то, но он может помочь кому-то.
- (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;
}
Нет ничего более прекрасного, чем индуктивное определение проблемы. Для полноты здесь у вас есть версия в прологе, которая может также прояснить соображения, лежащие в основе приведения лучей:
На основе алгоритма симуляции простоты в http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
Некоторые предикаты помощника:
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]).
Уравнение прямой с учетом 2 точек A и B (Линия (A,B)) имеет вид:
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
Важно, чтобы направление вращения линии было установлено по часовой стрелке для границ и против часовой стрелки для отверстий. Мы собираемся проверить, находится ли точка (X,Y), то есть проверенная точка, в левой полуплоскости нашей линии (это вопрос вкуса, это также может быть правая сторона, но также и направление границ линии должны быть изменены в этом случае), это проецирует луч из точки вправо (или влево) и подтверждает пересечение с линией. Мы решили проецировать луч в горизонтальном направлении (опять же, это вопрос вкуса, это также может быть сделано по вертикали с аналогичными ограничениями), поэтому мы имеем:
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
Теперь нам нужно знать, находится ли точка на левой (или правой) стороне только сегмента линии, а не на всей плоскости, поэтому нам нужно ограничить поиск только этим сегментом, но это легко, поскольку находиться внутри сегмента только одна точка на линии может быть выше, чем Y на вертикальной оси. 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).
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");
}
}
VBA VERSION:
Примечание. Помните, что если полигон является областью на карте, то широта / долгота являются значениями Y/X, а не X/Y (широта = Y, долгота = X) из-за того, что, как я понимаю, являются историческими последствиями, когда Долгота не была измерением.
МОДУЛЬ КЛАССА: 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
МОДУЛЬ:
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(q.X, minX)
maxX = vbMax(q.X, maxX)
minY = vbMin(q.Y, minY)
maxY = vbMax(q.Y, maxY)
Next i
If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > 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 > p.Y) Then
If (polygon(j).Y < p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - 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 < p.Y) Then
If (polygon(j).Y > p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - 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
p.X = <ENTER X VALUE HERE>
p.Y = <ENTER Y VALUE HERE>
' POLYGON OBJECT
Dim polygon() As CPoint
ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1
For i = 0 To <ENTER VALUE HERE> 'Same value as above
Set polygon(i) = New CPoint
polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through
polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through
Next i
InPolygon = isPointInPolygon(p, polygon)
MsgBox InPolygon
End Sub
Еще одна реализация numpyic, которая, на мой взгляд, является самой лаконичной из всех ответов на данный момент.
Например, предположим, что у нас есть многоугольник с впадинами многоугольника, который выглядит следующим образом:
2D-координаты вершин большого многоугольника:
[[139, 483], [227, 792], [482, 849], [523, 670], [352, 330]]
Координаты квадратной полости:
[[248, 518], [336, 510], [341, 614], [250, 620]]
Координаты полого треугольника:
[[416, 531], [505, 517], [495, 616]]
Скажем, мы хотим проверить две точки и находятся ли они в красной области (за исключением краев). Если мы найдем две точки, это будет выглядеть так:
Очевидно,
[296, 557] не находится внутри области чтения, тогда как
[422, 730] делает.
Мое решение основано на алгоритме числа оборотов . Ниже мой 4-строчный код Python, использующий только
numpy:
def detect(points, *polygons):
import numpy as np
endpoint1 = np.r_[tuple(np.roll(p, 1, 0) for p in polygons)][:, None] - points
endpoint2 = np.r_[polygons][:, None] - points
p1, p2 = np.cross(endpoint1, endpoint2), np.einsum('...i,...i', endpoint1, endpoint2)
return ~((p1.sum(0) < 0) ^ (abs(np.arctan2(p1, p2).sum(0)) > np.pi) | ((p1 == 0) & (p2 <= 0)).any(0))
Чтобы протестировать реализацию:
points = [[296, 557], [422, 730]]
polygon1 = [[139, 483], [227, 792], [482, 849], [523, 670], [352, 330]]
polygon2 = [[248, 518], [336, 510], [341, 614], [250, 620]]
polygon3 = [[416, 531], [505, 517], [495, 616]]
print(detect(points, polygon1, polygon2, polygon3))
Выход:
[False True]
Я тоже думал, что суммирование 360 работает только для выпуклых многоугольников, но это не так.
На этом сайте есть хорошая диаграмма, показывающая именно это, и хорошее объяснение тестирования на удар: Gamasutra - столкновение с Новым годом: обнаружение столкновений
C# версия ответа nirg здесь: я просто поделюсь кодом. Это может сэкономить кому-то время.
public static bool IsPointInPolygon(IList<Point> 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;
}
Удивило, что никто не поднял этот вопрос раньше, кроме тех прагматиков, которым нужна база данных: MongoDB имеет отличную поддержку запросов Geo, включая этот.
То, что вы ищете, это:
db.neighborhoods.findOne ({geometry: {$ geoIntersects: {$ geometry: {тип: "Точка", координаты: [ "долгота", "широта" ] } } } })
Neighborhoods это коллекция, которая хранит один или несколько полигонов в стандартном формате GeoJson. Если запрос возвращает значение NULL, он не пересекается, иначе это так.
Очень хорошо задокументировано здесь: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/
Производительность для более чем 6000 точек, классифицированных в сетке из 330 нерегулярных многоугольников, была меньше одной минуты без какой-либо оптимизации вообще, включая время для обновления документов соответствующим полигоном.
Порт.Net:
static void Main(string[] args)
{
Console.Write("Hola");
List<double> vertx = new List<double>();
List<double> verty = new List<double>();
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;
}
}
Это предположительно немного менее оптимизированная версия кода C отсюда, который был получен с этой страницы.
Моя версия на C++ использует std::vector<std::pair<double, double>>и два двойных символа x и y. Логика должна быть точно такой же, как и в исходном коде C, но мне легче читать мой. Я не могу говорить о спектакле.
bool point_in_poly(std::vector<std::pair<double, double>>& verts, double point_x, double point_y)
{
bool in_poly = false;
auto num_verts = verts.size();
for (int i = 0, j = num_verts - 1; i < num_verts; j = i++) {
double x1 = verts[i].first;
double y1 = verts[i].second;
double x2 = verts[j].first;
double y2 = verts[j].second;
if (((y1 > point_y) != (y2 > point_y)) &&
(point_x < (x2 - x1) * (point_y - y1) / (y2 - y1) + x1))
in_poly = !in_poly;
}
return in_poly;
}
Исходный код C
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;
}
Вот алгоритм быстрее, чем любой другой алгоритм для большинства случаев.
Он новый и элегантный. Мы тратим
Вместо трассировки лучей или углов вы можете получить значительно более быстрые результаты для нескольких проверок одного и того же полигона, используя таблицу сканирующего луча. Мы должны предварительно построить таблицу активных границ сканирующего луча, что и делает большая часть кода.
Мы вычисляем сканирующий луч и активные ребра для этой позиции в направлении Y. Мы составляем список точек, отсортированных по их y-компоненте, и проходим по этому списку для двух событий. Start-Y и End-Y, мы отслеживаем активные ребра по мере обработки списка. Мы обрабатываем события по порядку, и для каждого сканирующего луча мы записываем значение y события и активные ребра в каждом событии (события — начало-y и конец-y), но мы записываем их только тогда, когда наше событие-y отличается от последний раз (поэтому все в точке события обрабатывается до того, как мы отметим это в нашей таблице).
Получаем нашу таблицу:
- []
- р6р5, р6р7
- р6р5, р6р7, р2р3, р2р1
- р6р7, р5р4, р2р3, р3р1
- p7p8, p5p4, p2p3, p2p1
- p7p8, p5p4, p3p4, p2p1
- р7р8, р2р1
- р7р8, р1р0
- р8р0, р1р0
- []
Фактический код, выполняющий эту работу, находится всего в паре строк после построения этой таблицы.
- Примечание: код здесь использует комплексные значения в качестве точек. Так
def point_in_scantable(actives_table, events, xi, point):
beam = bisect(events, point.imag) - 1 # Binary search in sorted array.
actives_at_y = actives_table[beam]
total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
return bool(total % 2)
Мы бинарно ищем события, чтобы найти actives_at_y для определенного значения. Для всех активов в этом y мы вычисляем значение x сегментов в конкретном y нашей точки. Мы добавляем 1 каждый раз, когда этот x-перехват больше, чем x-компонента нашей точки. Затем мы модифицируем эту сумму на 2. (Это правило заполнения чет-нечет, вы можете легко адаптировать его к любому другому правилу заполнения).
Полный код:
from bisect import bisect
def build_edge_list(polygon):
edge_list = []
for i in range(1, len(polygon)):
if (polygon[i].imag, polygon[i].real) < (polygon[i - 1].imag, polygon[i - 1].real):
edge_list.append((polygon[i], i))
edge_list.append((polygon[i - 1], ~i))
else:
edge_list.append((polygon[i], ~i))
edge_list.append((polygon[i - 1], i))
def sort_key(e):
return e[0].imag, e[0].real, ~e[1]
edge_list.sort(key=sort_key)
return edge_list
def build_scanbeam(edge_list):
actives_table = []
events = []
y = -float("inf")
actives = []
for pt, index in edge_list:
if y != pt.imag:
actives_table.append(list(actives))
events.append(y)
if index >= 0:
actives.append(index)
else:
actives.remove(~index)
y = pt.imag
return actives_table, events
def point_in_polygon(polygon, point):
def x_intercept(e, y):
pt0 = polygon[e-1]
pt1 = polygon[e]
if pt1.real - pt0.real == 0:
return pt0.real
m = (pt1.imag - pt0.imag) / (pt1.real - pt0.real)
b = pt0.imag - (m * pt0.real)
return (y - b) / m
edge_list = build_edge_list(polygon)
actives_table, events = build_scanbeam(edge_list)
try:
if len(point):
return [point_in_scantable(actives_table, events, x_intercept, p) for p in point]
except TypeError:
return point_in_scantable(actives_table, events, x_intercept, point)
def point_in_scantable(actives_table, events, xi, point):
beam = bisect(events, point.imag) - 1 # Binary search in sorted array.
actives_at_y = actives_table[beam]
total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
return bool(total % 2)
Если мы проигнорируем, то
Есть несколько трюков с производительностью, например, мы можем снизить временную сложность до
Кроме того, поскольку это в основном становится таблицей поиска и некоторыми минимальными вычислениями для x-перехвата, это гораздо проще сделать с помощью графического процессора. Вам больше не нужно зацикливаться на многоугольнике. Таким образом, вы можете на самом деле массово вычислять эти точки с помощью таких вещей, как numpy, где вы получаете повышение производительности для выполнения всех вычислений одновременно.
Если вы используете SDK Google Map и хотите проверить, находится ли точка внутри многоугольника, вы можете попробовать использовать GMSGeometryContainsLocation. Отлично работает!! Вот как это работает
if GMSGeometryContainsLocation(point, polygon, true) {
print("Inside this polygon.")
} else {
print("outside this polygon")
}
Вот ссылка: https://developers.google.com/maps/documentation/ios-sdk/reference/group___geometry_utils
Вот точка в тесте полигонов в C, которая не использует приведение лучей. И это может работать для перекрывающихся областей (самопересечения), см. use_holes аргумент.
/* 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];
}
Примечание: это один из менее оптимальных методов, поскольку он включает в себя множество вызовов atan2f, но это может представлять интерес для разработчиков, читающих этот поток (в моих тестах он примерно в 23 раза медленнее, чем при использовании метода пересечения линий).
Если вы ищете библиотеку java-script, есть класс javascript google maps v3 для класса Polygon, чтобы определить, находится ли в нем точка.
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
Для полноты, вот lua-реализация алгоритма, предоставленного nirg и обсуждаемого Меки:
function pnpoly(area, test)
local inside = false
local tx, ty = table.unpack(test)
local j = #area
for i=1, #area do
local vxi, vyi = table.unpack(area[i])
local vxj, vyj = table.unpack(area[j])
if (vyi > ty) ~= (vyj > ty)
and tx < (vxj - vxi)*(ty - vyi)/(vyj - vyi) + vxi
then
inside = not inside
end
j = i
end
return inside
end
Переменная areaпредставляет собой таблицу точек, которые, в свою очередь, хранятся в виде 2D-таблиц. Пример:
> A = {{2, 1}, {1, 2}, {15, 3}, {3, 4}, {5, 3}, {4, 1.5}}
> T = {2, 1.1}
> pnpoly(A, T)
true
Ссылка на GitHub Gist.