MathCAD: попадание точки в многоугольник

Собственно, subj. В реализации всего один цикл с числом шагов, равным числу вершин многоугольника, так что назвать её быстрой, думаю, не будет преувеличением.

Многоугольник не должен передаваться функции point_in_polygon как замкнутый, то есть, координаты первой и последней точек не совпадают ("замыкание" делает сам алгоритм). Служебная функция close_p делает "замыкание" только для отображения на графике, заодно видно, как в MathCAD показать многоугольник и точку на графике:

Mathcad: как вывести многоугольник и точку на график

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

Алгоритм основан на методе лучей (учёта пересечений) и может быть описан так:

p - многоугольник, n - число его вершин, a - точка, 
i,j - номера очередной и следующей вершина
eps - погрешность, с которой точка считается лежащей на стороне многоугольника
для i от 1 до n нц
 если p.x==a.x и p.y==a.y, то вернуть 0 {точка на вершине}
 если i<n, то j=i+1, иначе j=0
 если p[i].y==p[j].y, то continue {перейти к началу цикла; горизонтальная линия}
 если p[i].y>a.y и p[j].y>a.y, то continue {оба конца стороны выше точки}
 если p[i].y<a.y и p[j].y<a.y, то continue {оба конца стороны ниже точки}
 если max(p[i].y,p[j].y)==a.y, то нач
  с=с+1 {если верхняя  вершина лежит на одной горизонтали с точкой,}
  f=1   {считаем пересечение на стороне}
 кон
 иначе если min(p[i].y,p[j].y)==a.y, то continue {если нижняя, то нет}
 иначе нач
  считаем t - коэффициент пропорциональности деления стороны горизонтальным лучом из точки a
  считаем t2 = p[i].x+t*(p[j].x-p[i].x)
  если t>0 и t<1 и t2<=a.x, то нач
   c=c+1 {считаем пересечение}
   если |t2-a.x|<eps, то f=1 {если с точностью eps совпадает с точкой, точка считается на стороне}
  кон
 кон
кц
если c mod 2==0, вернуть -1 {точка не лежит в полигоне}
иначе если f==1, вернуть 0 {точка лежит на стороне полигона}
иначе вернуть 1 {точка принадлежит многоугольнику}

Выпуклым наш полигон быть не обязан, для выпуклых многоугольников полно тривиальных алгоритмов.

Тестовая точка a задана сначала как вектор-строка, но с ними MathCAD нормально не работает, так что в функции point_in_polygon объект a транспонируется, превращаясь в вектор-столбец из двух элементов.

Ниже реализация выведена на скриншот и прикреплена как файл MathCAD 15.

Принадлежность точки произвольному многоугольнику в MathCAD

 Скачать пример "Принадлежность точки многоугольнику", документ MathCAD 15 в архиве .zip (33 Кб)

P.S. В окончательной версии я отказался от параметра eps из-за найденных несоответствий между возвращаемыми значениями -1 и 0, так что 0 возвращается только если точка совпадает с одной из вершин многоугольника.

А вот показанные на Хабре алгоритмы работать не захотели никак :( Ну или я крив был.

28.11.2014, 10:00 [12924 просмотра]