Как я могу определить, находится ли 2D-точка внутри многоугольника?
https://stackoverflow.com/questions/217578
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Я пытаюсь создать быстро Алгоритм 2D point inside polygon для использования при тестировании попаданий (например, Polygon.contains(p:Point)).Мы были бы признательны за предложения по эффективным методам.
Решение
Для графики я бы предпочел не использовать целые числа.Многие системы используют целые числа для отображения пользовательского интерфейса (в конце концов, пиксели - это целые числа), но macOS, например, использует float для всего.macOS знает только точки, и точка может быть преобразована в один пиксель, но в зависимости от разрешения монитора она может быть преобразована во что-то другое.На экранах retina половина точки (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)
:
Вам нужно протестировать луч со всех сторон.Рассмотрим луч как вектор, а каждую сторону - как вектор.Луч должен попасть в каждую сторону ровно по одному разу или вообще никогда.Он не может попасть в одну и ту же сторону дважды.Две прямые в двумерном пространстве всегда будут пересекаться ровно один раз, если только они не параллельны, и в этом случае они никогда не пересекаются.Однако, поскольку векторы имеют ограниченную длину, два вектора могут быть не параллельны и все равно никогда не пересекаться, потому что они слишком короткие, чтобы когда-либо встретиться друг с другом.
// 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-оборудование для решения этой проблемы, есть интересная альтернатива.Просто позвольте графическому процессору сделать всю работу за вас.Создайте окрашиваемую поверхность, которая находится за пределами экрана.Полностью заполните его черным цветом.Теперь позвольте OpenGL или Direct3D нарисовать ваш полигон (или даже все ваши полигоны, если вы просто хотите проверить, находится ли точка внутри какого-либо из них, но вам все равно, в каком именно) и залейте полигон (полигоны) другим цветом, напримерБелый.Чтобы проверить, находится ли точка внутри многоугольника, получите цвет этой точки с поверхности рисования.Это всего лишь выборка из памяти O (1).
Конечно, этот метод применим только в том случае, если ваша поверхность для рисования не обязательно должна быть огромной.Если он не может поместиться в память графического процессора, этот метод работает медленнее, чем выполнение его на центральном процессоре.Если он должен быть огромным, а ваш графический процессор поддерживает современные шейдеры, вы все равно можете использовать графический процессор, реализовав приведенное выше преобразование лучей в виде графического шейдера, что абсолютно возможно.При большем количестве полигонов или большом количестве тестируемых точек это окупится, учитывая, что некоторые графические процессоры смогут тестировать от 64 до 256 точек параллельно.Однако обратите внимание, что передача данных с CPU на GPU и обратно всегда обходится дорого, поэтому для простого тестирования пары точек на паре простых полигонов, где либо точки, либо полигоны динамичны и будут часто меняться, подход 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;
}
Аргументы
- нверт:Количество вершин в многоугольнике.Вопрос о том, повторять ли первую вершину в конце, обсуждался в статье, упомянутой выше.
- вершина, вершина:Массивы, содержащие x- и y-координаты вершин многоугольника.
- тестх, вспыльчивый:X- и y -координаты контрольной точки.
Он короткий и эффективный и работает как для выпуклых, так и для вогнутых многоугольников.Как предлагалось ранее, сначала вам следует проверить ограничивающий прямоугольник и обработать полигональные отверстия отдельно.
Идея, стоящая за этим, довольно проста.Автор описывает это следующим образом:
Я провожу полубесконечный луч по горизонтали (увеличивающий x, фиксированный y) от контрольной точки и подсчитываю, сколько ребер он пересекает.При каждом пересечении луч переключается между внутренним и внешним.Это называется теоремой о кривой Жордана.
Переменная c переключается с 0 на 1 и с 1 на 0 каждый раз, когда горизонтальный луч пересекает какое-либо ребро.Таким образом, в основном это отслеживание того, является ли количество пересеченных ребер четным или нечетным.0 означает четный, а 1 - нечетный.
Вот версия C # для ответ, данный nirg, который исходит из этот профессор 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, точка находится за пределами.
Мне больше нравится этот метод, потому что он более надежный и менее зависит от точности чисел.
Методы, которые вычисляют равномерность числа пересечений, ограничены, потому что вы можете «поразить» вершину во время вычисления количества пересечений.
РЕДАКТИРОВАТЬ: Кстати, этот метод работает с вогнутыми и выпуклыми многоугольниками.
РЕДАКТИРОВАТЬ. Недавно я обнаружил целую статью Википедии по этой теме.
статья Эрика Хейнса , цитируемая Бобобобо, действительно превосходна. Особенно интересны таблицы, сравнивающие производительность алгоритмов; метод суммирования углов действительно плох по сравнению с другими. Также интересно то, что оптимизации, такие как использование поисковой сетки для дальнейшего разбиения полигона на & Quot; in & Quot; и " 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
Swift-версия ответа по 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
}
}
Я проделал кое-какую работу над этим еще тогда, когда был исследователем в Майкл Стоунбрейкер - вы знаете, профессор, который придумал Ингрес, PostgreSQL, и т.д.
Мы поняли, что самый быстрый способ - это сначала сделать ограничивающую рамку, потому что это ОЧЕНЬ быстро.Если это находится за пределами ограничивающей рамки, значит, это снаружи.В противном случае вы выполняете более тяжелую работу...
Если вам нужен отличный алгоритм, обратитесь к исходному коду проекта PostgreSQL с открытым исходным кодом для работы с geo...
Я хочу отметить, что мы так и не получили никакого представления о соотношении правши и левши (это также можно выразить как проблему "внутри" и "снаружи"...
Обновить
Ссылка BKB предоставила большое количество разумных алгоритмов.Я работал над проблемами науки о Земле, и поэтому мне нужно было решение, которое работает по широте / долготе, и у него есть специфическая проблема управляемости - находится ли область внутри меньшей области или большей области?Ответ заключается в том, что имеет значение "направление" вершин - оно либо левостороннее, либо правостороннее, и таким образом вы можете указать любую область как "внутреннюю" любого заданного многоугольника.Таким образом, в моей работе использовалось третье решение, перечисленное на этой странице.
Кроме того, в моей работе использовались отдельные функции для тестов "на линии".
...С тех пор, как кто -то спросил:мы выяснили, что тесты с ограничивающими рамками были лучшими, когда количество вершин превышало некоторое число - проведите очень быстрый тест, прежде чем выполнять более длительный тест, если необходимо...Ограничивающий прямоугольник создается простым взятием наибольшего x, наименьшего x, наибольшего y и наименьшего y и соединением их вместе, чтобы получились четыре точки прямоугольника...
Еще один совет для тех, кто последует за мной:мы выполнили все наши более сложные вычисления с "уменьшением освещенности" в сеточном пространстве, все в положительных точках на плоскости, а затем перепроецировали обратно на "реальную" долготу / широту, таким образом избегая возможных ошибок обтекания при пересечении линии 180 долготы и при обработке полярных регионов.Сработало отлично!
Очень нравится решение, опубликованное 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;
}
Ответ Дэвида Сегонда в значительной степени является стандартным общим ответом, а Ричард Т - наиболее распространенная оптимизация, хотя есть и другие. Другие сильные оптимизации основаны на менее общих решениях. Например, если вы собираетесь проверять один и тот же многоугольник с множеством точек, триангуляция многоугольника может значительно ускорить процесс, поскольку существует ряд очень быстрых алгоритмов поиска TIN. Другой вариант: если многоугольник и точки находятся на ограниченной плоскости при низком разрешении, например, на экране, вы можете нарисовать многоугольник в отображаемом в память буфере отображения заданного цвета и проверить цвет данного пикселя, чтобы увидеть, лежит ли он в многоугольниках.
Как и многие другие оптимизации, они основаны на конкретных, а не на общих случаях, и приносят выгоды, основанные на амортизированном времени, а не на единичном использовании.
Работая в этой области, я обнаружил, что вычислительная геометрия Джозефа О'Руркса в C 'ISBN 0-521-44034-3 очень помогает.
Тривиальным решением было бы разделить многоугольник на треугольники и выполнить проверку треугольников, как описано здесь
Если ваш полигон CONVEX , возможно, есть более подходящий подход. Посмотрите на многоугольник как на коллекцию бесконечных линий. Каждая строка делит пространство на две части. для каждой точки легко сказать, находится ли она на одной стороне или другой стороне линии. Если точка находится на одной стороне всех линий, она находится внутри многоугольника.
Я понимаю, что это старый, но вот алгоритм наложения лучей, реализованный в Какао, на случай, если кому-то будет интересно. Не уверен, что это самый эффективный способ сделать что-то, но он может помочь кому-то.
- (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");
}
}
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;
}
Нет ничего более красивого, чем индуктивное определение проблемы. Для полноты здесь у вас есть версия в прологе, которая может также прояснить мысли, стоящие за приведением лучей :
На основе алгоритма симуляции простоты в 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]).
Уравнение прямой, заданной двумя точками A и B (линия (A, B)):
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
Важно, чтобы направление вращения линии установлен по часовой стрелке для границ и против часовой стрелки для дырок. Мы собираемся проверить, находится ли точка (X, Y), то есть проверенная точка слева полуплоскость нашей линии (это вопрос вкуса, это также может быть правая сторона, но и направление линий границ должно быть изменено в в этом случае), это для проецирования луча из точки вправо (или влево) и подтвердите пересечение с линией. Мы выбрали проект луч в горизонтальном направлении (опять же, дело вкуса, это также может быть сделано по вертикали с аналогичными ограничениями), поэтому мы имеем:
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
Теперь нам нужно знать, находится ли точка слева (или справа) от только отрезок, а не вся плоскость, поэтому нам нужно ограничить поиск только этим сегментом, но это легко, так как быть внутри сегмента, только одна точка в линии может быть выше чем Y по вертикальной оси. Поскольку это более сильное ограничение, должен быть первым, чтобы проверить, поэтому мы берем сначала только те строки выполнение этого требования, а затем проверить его возможность. Иордан Теорема кривой любой луч, проецируемый на многоугольник, должен пересекаться на четное количество строк. Итак, мы сделали, мы бросим луч к право, а затем каждый раз, когда он пересекает линию, переключать свое состояние. Однако в нашей реализации мы собираемся проверить длину пакет решений, отвечающих заданным ограничениям, и решите участие в этом. для каждой линии в многоугольнике это должно быть сделано.
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).
Версия 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;
}
}
.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;
}
}
ВЕРСИЯ VBA:
Примечание:Помните, что если ваш полигон является областью на карте, то Широта / Долгота - это значения Y / X, в отличие от X / Y (Широта = Y, Долгота = X), поскольку, насколько я понимаю, это исторические последствия с давних времен, когда Долгота не была измерением.
МОДУЛЬ КЛАССА:Точка
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
Я создал реализацию Python для у нирга c++ код:
Входные данные
- граничные точки: узлы, из которых состоит полигон.
ограничение_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
Опять же, идея взята из здесь
Вот точка в тесте полигонов в 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 раза медленнее, чем при использовании метода пересечения линии ).
Для обнаружения попадания в полигон нам нужно проверить две вещи:
<Ол>Для решения следующих особых случаев в Алгоритм приведения лучей :
<Ол>Проверьте определение того, находится ли точка внутри сложного многоугольника . В статье представлен простой способ их решения, поэтому для указанных выше случаев не требуется особого подхода.
Это можно сделать, проверив, совпадает ли область, образованная путем соединения нужной точки с вершинами вашего многоугольника, с областью самого многоугольника.
Или вы можете проверить, равна ли сумма внутренних углов от вашей точки каждой паре двух последовательных вершин многоугольника до вашей контрольной точки сумме 360, но у меня есть ощущение, что первый вариант быстрее, потому что он не включает деления, ни вычисления обратных тригонометрических функций.
Я не знаю, что произойдет, если у вашего многоугольника есть дыра внутри, но мне кажется, что основная идея может быть адаптирована к этой ситуации
Вы также можете опубликовать вопрос в математическом сообществе. Могу поспорить, у них есть миллион способов сделать это
Если вы ищете библиотеку java-script, для класса Polygon есть расширение javascript google maps v3, чтобы определить, находится ли в нем точка.
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
При использовании qt (Qt 4.3+) можно использовать функцию QPolygon containsPoint
Ответ зависит от того, есть ли у вас простые или сложные полигоны. Простые многоугольники не должны иметь пересечений отрезков. Таким образом, они могут иметь отверстия, но линии не могут пересекаться друг с другом. Сложные области могут иметь пересечения линий, поэтому они могут иметь перекрывающиеся области или области, которые соприкасаются друг с другом только на одну точку.
Для простых полигонов лучший алгоритм - это алгоритм преобразования лучей (число пересечений). Для сложных полигонов этот алгоритм не обнаруживает точки, которые находятся внутри перекрывающихся областей. Так что для сложных полигонов вы должны использовать алгоритм числа обмоток.
Вот отличная статья с реализацией C обоих алгоритмов. Я попробовал их, и они хорошо работают.
Scala-версия решения по nirg (предполагается, что предварительная проверка ограничивающего прямоугольника выполняется отдельно):
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 > 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
oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
}
}
oddIntersections(0, length - 1, tracker = false)
}
Вот версия ответа @nirg на golang (вдохновленная кодом C # от @@m-katz)
func isPointInPolygon(polygon []point, testp point) bool {
minX := polygon[0].X
maxX := polygon[0].X
minY := polygon[0].Y
maxY := polygon[0].Y
for _, p := range polygon {
minX = min(p.X, minX)
maxX = max(p.X, maxX)
minY = min(p.Y, minY)
maxY = max(p.Y, maxY)
}
if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
return false
}
inside := false
j := len(polygon) - 1
for i := 0; i < len(polygon); i++ {
if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
inside = !inside
}
j = i
}
return inside
}
Удивлен никто раньше, но для прагматиков, которым требуется база данных: MongoDB имеет отличную поддержку запросов Geo, включая этот.
То, что вы ищете, это:
db.neighborhoods.findOne ({geometry: {$ geoIntersects: {$ geometry: { тип: " точка " ;, координаты: [" долгота " ;, " широта " ]}}} })
Neighborhoods - это коллекция, в которой хранится один или несколько полигонов в стандартном формате GeoJson. Если запрос возвращает ноль, он не пересекается, в противном случае это так.
Очень хорошо задокументировано здесь:
https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/
Производительность для более чем 6000 точек, классифицированных в сетке из 330 нерегулярных полигонов, была меньше одной минуты без какой-либо оптимизации вообще, включая время для обновления документов соответствующим полигоном.
