Реализации алгоритмов/Задача о принадлежности точки многоугольнику
В вычислительной геометрии известна задача об определении принадлежности точки многоугольнику. На плоскости даны многоугольник и точка. Требуется решить вопрос о принадлежности точки многоугольнику.
Многоугольник может быть как выпуклым, так и невыпуклым. Обычно предполагается, что многоугольник простой, т.е. без самопересечений, но задачу рассматривают и для не-простых многоугольников. В последнем случае разные способы определения принадлежности точки многоугольнику могут привести к разным результатам. Различают алгоритмы без предварительной обработки и алгоритмы с предварительной обработкой, в ходе которой создаются некоторые структуры данных, позволяющие в дальнейшем быстрее отвечать на множество запросов о принадлежности точек одному и тому же многоугольнику.
Алгоритм определяет точки границ многоугольника как точки, ему принадлежащие.
Описание
правитьДля того чтобы все результаты вычислений в программе могли быть представлены целочисленными переменными (манипулирование данными целого типа повышает быстродействие программы и является естественным для приложений компьютерной графики), вычисления и сравнения площадей треугольников заменяются вычислениями и сравнениями их удвоенных площадей. Тем самым исключается погрешность округления при программной реализации всего алгоритма, в целом.
Аргументами функции, реализующей проверку принадлежности данной точки данному многоугольнику произвольного вида, являются
- указатель на массив пар целочисленных координат вершин многоугольника, а именно, на массив структур вида
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 примет нечетное значение, то точка лежит внутри многоугольника. В противном случает точка лежит вне заданого многоугольника.
Замечание: В данной реализации алгоритма луч направлен вниз.