Как определить, находится ли 2D-точка внутри многоугольника?
Я пытаюсь создать быструю 2D точку внутри алгоритма многоугольника, для использования в тестировании попаданий (например, Polygon.contains(p:Point)
). Предложения для эффективных методов будут оценены.
43 ответа
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)
}
Барицентрические координаты можно использовать для описания положения точки внутри треугольника. Они представляют собой отношения расстояний между точкой и вершинами треугольника.
В треугольнике с вершинами A, B и C любая точка P внутри треугольника может быть выражена с использованием барицентрических координат (u, v, w), где u, v и w представляют собой отношения расстояний между P и вершинами. А, Б и С соответственно. Координаты должны удовлетворять условию u + v + w = 1.
Стоит отметить, что барицентрические координаты применимы к любому многоугольнику, независимо от его формы.
Для вогнутого многоугольника барицентрические координаты точки внутри многоугольника по-прежнему определяются отношениями расстояний между точкой и вершинами. Единственное отличие состоит в том, что некоторые барицентрические координаты могут быть отрицательными, что указывает на то, что точка находится вне многоугольника. В таких случаях сумма координат не будет равна 1.
Аналогично, для выпуклого многоугольника барицентрические координаты точки внутри многоугольника всегда будут неотрицательными и в сумме будут равны 1.
Является ли форма вогнутой или выпуклой, не влияет на вычисление барицентрических координат. Ключевым моментом является обеспечение правильного расчета координат на основе относительных расстояний между точкой и вершинами многоугольника.
Ниже вы можете увидеть реализацию Matlab для формы с четырьмя вершинами.
function isInside = isPointInsideQuadrilateral(x, y, x1, y1, x2, y2, x3, y3, x4, y4)
% Calculate barycentric coordinates
denominator1 = (y2 - y1)*(x4 - x1) + (x1 - x2)*(y4 - y1);
alpha1 = ((y2 - y1)*(x - x1) + (x1 - x2)*(y - y1)) / denominator1;
beta1 = ((y1 - y4)*(x - x1) + (x4 - x1)*(y - y1)) / denominator1;
gamma1 = 1 - alpha1 - beta1;
denominator2 = (y3 - y2)*(x4 - x2) + (x2 - x3)*(y4 - y2);
alpha2 = ((y3 - y2)*(x - x2) + (x2 - x3)*(y - y2)) / denominator2;
beta2 = ((y2 - y4)*(x - x2) + (x4 - x2)*(y - y2)) / denominator2;
gamma2 = 1 - alpha2 - beta2;
isInside = (alpha1 >= 0 && beta1 >= 0 && gamma1 >= 0) || (alpha2 >= 0 && beta2 >= 0 && gamma2 >= 0);
end
Этот метод довольно быстрый и имеет высокую точность.
Вот golang-версия ответа @nirg (вдохновлена кодом 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
}
Для решения следующих особых случаев в алгоритме приведения лучей:
- Луч перекрывает одну из сторон многоугольника.
- Точка находится внутри многоугольника, и луч проходит через вершину многоугольника.
- Точка находится за пределами многоугольника, и луч просто касается одного из углов многоугольника.
Проверьте, находится ли точка внутри сложного многоугольника. В статье представлен простой способ их решения, поэтому для указанных выше случаев не требуется никакого специального лечения.
При использовании qt (Qt 4.3+) можно использовать функцию QPolygon containsPoint
Ответ зависит от того, есть ли у вас простые или сложные полигоны. Простые многоугольники не должны иметь пересечений отрезков. Таким образом, они могут иметь отверстия, но линии не могут пересекаться друг с другом. Сложные области могут иметь пересечения линий, поэтому они могут иметь перекрывающиеся области или области, которые касаются друг друга только одной точкой.
Для простых полигонов лучшим алгоритмом является алгоритм преобразования лучей (число пересечений). Для сложных полигонов этот алгоритм не обнаруживает точки, которые находятся внутри перекрывающихся областей. Так что для сложных полигонов вы должны использовать алгоритм числа обмоток.
Вот отличная статья с реализацией C обоих алгоритмов. Я попробовал их, и они хорошо работают.
Вы можете сделать это, проверив, совпадает ли область, образованная путем соединения желаемой точки с вершинами вашего многоугольника, с областью самого многоугольника.
Или вы можете проверить, что сумма внутренних углов от вашей точки до каждой пары из двух последовательных вершин многоугольника до вашей контрольной точки равна 360, но у меня есть ощущение, что первый вариант быстрее, потому что он не включает деления или вычисления инверсии тригонометрических функций.
Я не знаю, что произойдет, если в вашем многоугольнике будет дыра, но мне кажется, что основная идея может быть адаптирована к этой ситуации.
Вы также можете опубликовать вопрос в математическом сообществе. Могу поспорить, у них есть миллион способов сделать это
Как и в ответе Дэвида Сегонда, я использую подход суммирования углов, полученный из моего алгоритма рисования вогнутых многоугольников . Он основан на суммировании приблизительных углов подтреугольников вокруг точки для получения веса. Вес вокруг
1.0
означает, что точка находится внутри треугольника, вес вокруг
0.0
означает снаружи, вес вокруг
-1.0
это то, что происходит, когда внутри многоугольника, но в обратном порядке (например, с одной из половин четырехугольника в форме галстука-бабочки) и весом
NAN
если точно по краю. Причина, по которой это не медленно, заключается в том, что углы вообще не нужно точно оценивать. Отверстия можно обрабатывать, рассматривая их как отдельные полигоны и вычитая веса.
typedef struct { double x, y; } xy_t;
xy_t sub_xy(xy_t a, xy_t b)
{
a.x -= b.x;
a.y -= b.y;
return a;
}
double calc_sharp_subtriangle_pixel_weight(xy_t p0, xy_t p1)
{
xy_t rot, r0, r1;
double weight;
// Rotate points (unnormalised)
rot = sub_xy(p1, p0);
r0.x = rot.x*p0.y - rot.y*p0.x;
r0.y = rot.x*p0.x + rot.y*p0.y;
r1.y = rot.x*p1.x + rot.y*p1.y;
// Calc weight
weight = subtriangle_angle_approx(r1.y, r0.x) - subtriangle_angle_approx(r0.y, r0.x);
return weight;
}
double calc_sharp_polygon_pixel_weight(xy_t p, xy_t *corner, int corner_count)
{
int i;
xy_t p0, p1;
double weight = 0.;
p0 = sub_xy(corner[corner_count-1], p);
for (i=0; i < corner_count; i++)
{
// Transform corner coordinates
p1 = sub_xy(corner[i], p);
// Calculate weight for each subtriangle
weight += calc_sharp_subtriangle_pixel_weight(p0, p1);
p0 = p1;
}
return weight;
}
Таким образом, для каждого сегмента многоугольника формируется подтреугольник с оцениваемой точкой, затем каждый подтреугольник поворачивается, чтобы оценить его приблизительные углы и добавить к весу.
Звонки в
subtriangle_angle_approx(y, x)
можно заменить на
atan2(y, x) / (2.*pi)
, однако очень грубое приближение будет достаточно точным:
double subtriangle_angle_approx(double y, double x)
{
double angle, d;
int obtuse;
if (x == 0.)
return NAN;
obtuse = fabs(y) > fabs(x);
if (obtuse)
swap_double(&y, &x);
// Core of the approximation, a very loosely approximate atan(y/x) / (2.*pi) over ]-1 , 1[
d = y / x;
angle = 0.13185 * d;
if (obtuse)
angle = sign(d)*0.25 - angle;
return angle;
}
Кажется, это работает в R (извиняюсь за уродство, хотел бы увидеть лучшую версию!).
pnpoly <- function(nvert,vertx,verty,testx,testy){
c <- FALSE
j <- nvert
for (i in 1:nvert){
if( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i])*(testy-verty[i])/(verty[j]-verty[i])+vertx[i]))
{c <- !c}
j <- i}
return(c)}
from typing import Iterable
def pnpoly(verts, x, y):
#check if x and/or y is iterable
xit, yit = isinstance(x, Iterable), isinstance(y, Iterable)
#if not iterable, make an iterable of length 1
X = x if xit else (x, )
Y = y if yit else (y, )
#store verts length as a range to juggle j
r = range(len(verts))
#final results if x or y is iterable
results = []
#traverse x and y coordinates
for xp in X:
for yp in Y:
c = 0 #reset c at every new position
for i in r:
j = r[i-1] #set j to position before i
#store a few arguments to shorten the if statement
yneq = (verts[i][1] > yp) != (verts[j][1] > yp)
xofs, yofs = (verts[j][0] - verts[i][0]), (verts[j][1] - verts[i][1])
#if we have crossed a line, increment c
if (yneq and (xp < xofs * (yp - verts[i][1]) / yofs + verts[i][0])):
c += 1
#if c is odd store the coordinates
if c%2:
results.append((xp, yp))
#return either coordinates or a bool, depending if x or y was an iterable
return results if (xit or yit) else bool(c%2)
Эта версия Python универсальна. Вы можете ввести одно значение x и одно значение y для получения результата True / False или использовать
range
для
x
а также
y
чтобы пройти всю сетку точек. Если используются диапазоны
list
пар x / y для всех
True
баллы возвращается. В
vertices
аргумент ожидает двумерного
Iterable
пар x / y, таких как:
[(x1,y1), (x2,y2), ...]
пример использования:
vertices = [(25,25), (75,25), (75,75), (25,75)]
pnpoly(vertices, 50, 50) #True
pnpoly(vertices, range(100), range(100)) #[(25,25), (25,26), (25,27), ...]
На самом деле, даже это сработает.
pnpoly(vertices, 50, range(100)) #check 0 to 99 y at x of 50
pnpoly(vertices, range(100), 50) #check 0 to 99 x at y of 50
Для обнаружения попадания в полигон нам нужно протестировать две вещи:
- Если точка находится внутри полигона. (может быть выполнено с помощью алгоритма преобразования лучей)
- Если точка находится на границе многоугольника (может быть выполнено тем же алгоритмом, который используется для обнаружения точки на полилинии (линии)).
Это работает только для выпуклых фигур, но Minkowski Portal Refinement и GJK также являются отличными вариантами для проверки, находится ли точка в многоугольнике. Вы используете вычитание Минковского, чтобы вычесть точку из многоугольника, затем запустите эти алгоритмы, чтобы увидеть, содержит ли многоугольник начало координат.
Кроме того, что интересно, вы можете описывать свои фигуры немного более неявно, используя вспомогательные функции, которые принимают вектор направления в качестве входных данных и выплевывают самую дальнюю точку вдоль этого вектора. Это позволяет вам описать любую выпуклую форму... изогнутую, изготовленную из многоугольников или смешанную. Вы также можете выполнить операции, чтобы объединить результаты простых функций поддержки для создания более сложных фигур.
Дополнительная информация: http://xenocollide.snethen.com/mpr2d.html
Также в игровом программировании Gems 7 рассказывается, как это сделать в 3d (:
Вот Rust-версия ответа @nirg (версия javascript Филиппа Ленссена). Я даю этот ответ, потому что я получаю много помощи от этого сайта, и я перевожу версию javascript на ржавчину для исключения и надеюсь, что может помочь кому-то, последняя причина в том, что в моя работа, я переведу этот код как wasm, чтобы улучшить производительность моего холста, это для начала. мой плохой английский ааа ..., простите меня `
pub struct Point {
x: f32,
y: f32,
}
pub fn point_is_in_poly(pt: Point, polygon: &Vec<Point>) -> bool {
let mut is_inside = false;
let max_x = polygon.iter().map(|pt| pt.x).reduce(f32::max).unwrap();
let min_x = polygon.iter().map(|pt| pt.x).reduce(f32::min).unwrap();
let max_y = polygon.iter().map(|pt| pt.y).reduce(f32::max).unwrap();
let min_y = polygon.iter().map(|pt| pt.y).reduce(f32::min).unwrap();
if pt.x < min_x || pt.x > max_x || pt.y < min_y || pt.y > max_y {
return is_inside;
}
let len = polygon.len();
let mut j = len - 1;
for i in 0..len {
let y_i_value = polygon[i].y > pt.y;
let y_j_value = polygon[j].y > pt.y;
let last_check = (polygon[j].x - polygon[i].x) * (pt.y - polygon[i].y)
/ (polygon[j].y - polygon[i].y)
+ polygon[i].x;
if y_i_value != y_j_value && pt.x < last_check {
is_inside = !is_inside;
}
j = i;
}
is_inside
}
let pt = Point {
x: 1266.753,
y: 97.655,
};
let polygon = vec![
Point {
x: 725.278,
y: 203.586,
},
Point {
x: 486.831,
y: 441.931,
},
Point {
x: 905.77,
y: 445.241,
},
Point {
x: 1026.649,
y: 201.931,
},
];
let pt1 = Point {
x: 725.278,
y: 203.586,
};
let pt2 = Point {
x: 872.652,
y: 321.103,
};
println!("{}", point_is_in_poly(pt, &polygon));// false
println!("{}", point_is_in_poly(pt1, &polygon)); // true
println!("{}", point_is_in_poly(pt2, &polygon));// true
`