Реализации алгоритмов/Задача о принадлежности точки многоугольнику

В вычислительной геометрии известна задача об определении принадлежности точки многоугольнику. На плоскости даны многоугольник и точка. Требуется решить вопрос о принадлежности точки многоугольнику.

Многоугольник может быть как выпуклым, так и невыпуклым. Обычно предполагается, что многоугольник простой, т.е. без самопересечений, но задачу рассматривают и для не-простых многоугольников. В последнем случае разные способы определения принадлежности точки многоугольнику могут привести к разным результатам. Различают алгоритмы без предварительной обработки и алгоритмы с предварительной обработкой, в ходе которой создаются некоторые структуры данных, позволяющие в дальнейшем быстрее отвечать на множество запросов о принадлежности точек одному и тому же многоугольнику.

Алгоритм определяет точки границ многоугольника как точки, ему принадлежащие.

Описание править

Для того чтобы все результаты вычислений в программе могли быть представлены целочисленными переменными (манипулирование данными целого типа повышает быстродействие программы и является естественным для приложений компьютерной графики), вычисления и сравнения площадей треугольников заменяются вычислениями и сравнениями их удвоенных площадей. Тем самым исключается погрешность округления при программной реализации всего алгоритма, в целом.

Аргументами функции, реализующей проверку принадлежности данной точки данному многоугольнику произвольного вида, являются

указатель на массив пар целочисленных координат вершин многоугольника, а именно, на массив структур вида

    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 примет нечетное значение, то точка лежит внутри многоугольника. В противном случает точка лежит вне заданого многоугольника.

Замечание: В данной реализации алгоритма луч направлен вниз.