Реализации алгоритмов/Задача о принадлежности точки многоугольнику
В вычислительной геометрии известна задача об определении принадлежности точки многоугольнику. На плоскости даны многоугольник и точка. Требуется решить вопрос о принадлежности точки многоугольнику.
Многоугольник может быть как выпуклым, так и невыпуклым. Обычно предполагается, что многоугольник простой, т.е. без самопересечений, но задачу рассматривают и для не-простых многоугольников. В последнем случае разные способы определения принадлежности точки многоугольнику могут привести к разным результатам. Различают алгоритмы без предварительной обработки и алгоритмы с предварительной обработкой, в ходе которой создаются некоторые структуры данных, позволяющие в дальнейшем быстрее отвечать на множество запросов о принадлежности точек одному и тому же многоугольнику.
Алгоритм определяет точки границ многоугольника как точки, ему принадлежащие.
Описание
[править]Для того чтобы все результаты вычислений в программе могли быть представлены целочисленными переменными (манипулирование данными целого типа повышает быстродействие программы и является естественным для приложений компьютерной графики), вычисления и сравнения площадей треугольников заменяются вычислениями и сравнениями их удвоенных площадей. Тем самым исключается погрешность округления при программной реализации всего алгоритма, в целом.
Аргументами функции, реализующей проверку принадлежности данной точки данному многоугольнику произвольного вида, являются
- указатель на массив пар целочисленных координат вершин многоугольника, а именно, на массив структур вида
struct Point {
int x;
int y;
};
- число вершин многоугольника;
- целочисленное значение координаты X заданной точки;
- целочисленное значение координаты Y заданной точки.
Функция возвращает 1, если точка принадлежит многоугольнику, иначе — 0.
Функция имеет следующий вид.
int IsPointInsidePolygon (Point *p, int Number, int x, int y)
{
int i1, i2, n, N, S, S1, S2, S3, flag;
N = Number;
for (n=0; n<N; n++)
{
flag = 0;
i1 = n < N-1 ? n + 1 : 0;
while (flag == 0)
{
i2 = i1 + 1;
if (i2 >= N)
i2 = 0;
if (i2 == (n < N-1 ? n + 1 : 0))
break;
S = abs (p[i1].x * (p[i2].y - p[n ].y) +
p[i2].x * (p[n ].y - p[i1].y) +
p[n].x * (p[i1].y - p[i2].y));
S1 = abs (p[i1].x * (p[i2].y - y) +
p[i2].x * (y - p[i1].y) +
x * (p[i1].y - p[i2].y));
S2 = abs (p[n ].x * (p[i2].y - y) +
p[i2].x * (y - p[n ].y) +
x * (p[n ].y - p[i2].y));
S3 = abs (p[i1].x * (p[n ].y - y) +
p[n ].x * (y - p[i1].y) +
x * (p[i1].y - p[n ].y));
if (S == S1 + S2 + S3)
{
flag = 1;
break;
}
i1 = i1 + 1;
if (i1 >= N)
i1 = 0;
break;
}
if (flag == 0)
break;
}
return flag;
}
Очень быстрый алгоритм
[править]В основе алгоритма лежит идея подсчёта количества пересечений луча, исходящего из данной точки в направлении горизонтальной оси, со сторонами многоугольника. Если оно чётное, точка не принадлежит многоугольнику. В данном алгоритме луч направлен влево.
bool pnpoly(int npol, float * xp, float * yp, float x, float y)
{
bool c = false;
for (int i = 0, j = npol - 1; i < npol; j = i++)
{
if ((((yp[i] <= y) && (y < yp[j])) || ((yp[j] <= y) && (y < yp[i]))) &&
(((yp[j] - yp[i]) != 0) && (x > ((xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))))
c = !c;
}
return c;
}
Замечание: Так как умножение быстрее деления, условие можно записать так:
int pnpoly(int npol, float * xp, float * yp, float x, float y)
{
int c = 0;
for (int i = 0, j = npol - 1; i < npol; j = i++)
{
if ((
(yp[i] < yp[j]) && (yp[i] <= y) && (y <= yp[j]) &&
((yp[j] - yp[i]) * (x - xp[i]) > (xp[j] - xp[i]) * (y - yp[i]))
) || (
(yp[i] > yp[j]) && (yp[j] <= y) && (y <= yp[i]) &&
((yp[j] - yp[i]) * (x - xp[i]) < (xp[j] - xp[i]) * (y - yp[i]))
))
c = !c;
}
return c;
}
Однако, стоит заметить, что данный алгоритм не эквивалентен предыдущему, поэтому его использование может привести к неправильным результатам.
Perl
[править] my $x = -40; my $y = -60; # Проверяемая точка
my @xp = (-73,-33,7,-33); # Массив X-координат полигона
my @yp = (-85,-126,-85,-45); # Массив Y-координат полигона
&InPoly(\@xp,\@yp,$x,$y);
sub InPoly()
{
my($xp, $yp, $x, $y) = @_;
my $npol = @{$xp};
my $j = $npol - 1;
my $c = 0;
for(my $i = 0; $i < $npol;$i++) {
if ((((@{$yp}[$i]<=$y) && ($y<@{$yp}[$j])) || ((@{$yp}[$j]<=$y) && ($y<@{$yp}[$i]))) &&
($x > (@{$xp}[$j] - @{$xp}[$i]) * ($y - @{$yp}[$i]) / (@{$yp}[$j] - @{$yp}[$i]) + @{$xp}[$i]))
{
$c = !$c
}
$j = $i;
}
return $c;
}
Delphi (Object Pascal)
[править] type
tPolygon = array of tPoint; //tPoint - это запись, с двумя полями, x и y
...
function IsMouseInPoly(x,y: integer; myP: tPolygon): boolean; //x и y - это координаты мыши
var //myP - массив с вершинами полигона
i,j,npol: integer;
inPoly: boolean;
begin
npol:=length(myP)-1;
j:=npol;
inPoly:=false;
for i:=0 to npol do
begin
if ((((myP[i].y<=y) and (y<myP[j].y)) or ((myP[j].y<=y) and (y<myP[i].y))) and
(x>(myP[j].x-myP[i].x)*(y-myP[i].y) / (myP[j].y-myP[i].y)+myP[i].x))
then inPoly:=not inPoly;
j:=i;
end;
result:=inPoly;
end;
JavaScript
[править] var x = -40;
var y = -60;
var xp = new Array(-73,-33,7,-33); // Массив X-координат полигона
var yp = new Array(-85,-126,-85,-45); // Массив Y-координат полигона
function inPoly(x,y){
var npol = xp.length;
var j = npol - 1;
var c = 0;
for (var i = 0; i < npol;i++){
if ((((yp[i]<=y) && (y<yp[j])) || ((yp[j]<=y) && (y<yp[i]))) &&
(x > (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i])) {
c = !c
}
j = i;
}
return c;
}
inPoly(x,y);
Python 3
[править]На Python программа несколько отличается от других языков в сторону компактности из-за особенностей адресации элементов массива. Не нужны дополнительные переменные. Не работает с многоугольниками вогнутого типа.
def inPolygon(x, y, xp, yp):
c=0
for i in range(len(xp)):
if (((yp[i]<=y and y<yp[i-1]) or (yp[i-1]<=y and y<yp[i])) and
(x > (xp[i-1] - xp[i]) * (y - yp[i]) / (yp[i-1] - yp[i]) + xp[i])): c = 1 - c
return c
print( inPolygon(100, 0, (-100, 100, 100, -100), (100, 100, -100, -100)))
Быстрый алгоритм для случая, когда луч пересекает одну или несколько вершин
[править]Функция Cross определяет, пересекает ли луч j-ое ребро многоугольника:
bool Cross(int j)
{
int first = j;
int second = j == n - 1 ? 0 : j + 1;
double y = (xh - points[first].x) * (points[second].y - points[first].y) / (points[second].x - points[first].x) + points[first].y;
double minimal = min(points[first].x, points[second].x);
double maximal = max(points[first].x, points[second].x);
return (points[first].x != points[second].x) && (yh >= y) && (xh > minimal) && (xh <= maximal);
}
Фрагмент основной программы:
...
int count = 0;
for (int i = 0; i < n; i++)
{
count += Cross(i);
}
...
Если переменная count примет нечетное значение, то точка лежит внутри многоугольника. В противном случает точка лежит вне заданого многоугольника.
Замечание: В данной реализации алгоритма луч направлен вниз.