Точка в алгоритме многоугольника
Я видел, что приведенный ниже алгоритм работает, чтобы проверить, находится ли точка в данном многоугольнике по этой ссылке:
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
Я попробовал этот алгоритм, и он на самом деле работает просто отлично. Но, к сожалению, я не могу понять это хорошо, потратив некоторое время, пытаясь понять это.
Поэтому, если кто-то может понять этот алгоритм, пожалуйста, объясните мне его немного.
Спасибо.
13 ответов
Алгоритм направляет луч вправо. На каждой итерации цикла контрольная точка проверяется по одному из ребер многоугольника. Первая строка if-теста завершается успешно, если y-координата точки находится в области видимости края. Вторая строка проверяет, находится ли контрольная точка слева от линии (я думаю - у меня нет под рукой макулатуры). Если это так, линия, проведенная вправо от контрольной точки, пересекает этот край.
Путем многократного инвертирования значения cалгоритм подсчитывает, сколько раз правая линия пересекает многоугольник. Если он пересекает нечетное количество раз, то точка находится внутри; если четное число, точка находится за пределами.
У меня были бы проблемы с а) точностью арифметики с плавающей точкой и б) эффектом наличия горизонтального ребра или контрольной точки с той же координатой y, что и у вершины.
Чоулетт прав во всех отношениях, формах и формах. Алгоритм предполагает, что если ваша точка находится на линии многоугольника, то это за пределами - в некоторых случаях это неверно. Изменение двух операторов ">" на ">=" и изменение "<" на "<=" исправят это.
bool PointInPolygon(Point point, Polygon polygon) {
vector<Point> points = polygon.getPoints();
int i, j, nvert = points.size();
bool c = false;
for(i = 0, j = nvert - 1; i < nvert; j = i++) {
if( ( (points[i].y >= point.y ) != (points[j].y >= point.y) ) &&
(point.x <= (points[j].x - points[i].x) * (point.y - points[i].y) / (points[j].y - points[i].y) + points[i].x)
)
c = !c;
}
return c;
}
Я изменил исходный код, чтобы сделать его немного более читабельным (также он использует Eigen). Алгоритм идентичен.
// This uses the ray-casting algorithm to decide whether the point is inside
// the given polygon. See https://en.wikipedia.org/wiki/Point_in_polygon#Ray_casting_algorithm
bool pnpoly(const Eigen::MatrixX2d &poly, float x, float y)
{
// If we never cross any lines we're inside.
bool inside = false;
// Loop through all the edges.
for (int i = 0; i < poly.rows(); ++i)
{
// i is the index of the first vertex, j is the next one.
// The original code uses a too-clever trick for this.
int j = (i + 1) % poly.rows();
// The vertices of the edge we are checking.
double xp0 = poly(i, 0);
double yp0 = poly(i, 1);
double xp1 = poly(j, 0);
double yp1 = poly(j, 1);
// Check whether the edge intersects a line from (-inf,y) to (x,y).
// First check if the line crosses the horizontal line at y in either direction.
if ((yp0 <= y) && (yp1 > y) || (yp1 <= y) && (yp0 > y))
{
// If so, get the point where it crosses that line. This is a simple solution
// to a linear equation. Note that we can't get a division by zero here -
// if yp1 == yp0 then the above if be false.
double cross = (xp1 - xp0) * (y - yp0) / (yp1 - yp0) + xp0;
// Finally check if it crosses to the left of our test point. You could equally
// do right and it should give the same result.
if (cross < x)
inside = !inside;
}
}
return inside;
}
Это может быть как можно более подробно для объяснения алгоритма трассировки лучей в реальном коде. Это не может быть оптимизировано, но это всегда должно происходить после полного понимания системы.
//method to check if a Coordinate is located in a polygon
public boolean checkIsInPolygon(ArrayList<Coordinate> poly){
//this method uses the ray tracing algorithm to determine if the point is in the polygon
int nPoints=poly.size();
int j=-999;
int i=-999;
boolean locatedInPolygon=false;
for(i=0;i<(nPoints);i++){
//repeat loop for all sets of points
if(i==(nPoints-1)){
//if i is the last vertex, let j be the first vertex
j= 0;
}else{
//for all-else, let j=(i+1)th vertex
j=i+1;
}
float vertY_i= (float)poly.get(i).getY();
float vertX_i= (float)poly.get(i).getX();
float vertY_j= (float)poly.get(j).getY();
float vertX_j= (float)poly.get(j).getX();
float testX = (float)this.getX();
float testY = (float)this.getY();
// following statement checks if testPoint.Y is below Y-coord of i-th vertex
boolean belowLowY=vertY_i>testY;
// following statement checks if testPoint.Y is below Y-coord of i+1-th vertex
boolean belowHighY=vertY_j>testY;
/* following statement is true if testPoint.Y satisfies either (only one is possible)
-->(i).Y < testPoint.Y < (i+1).Y OR
-->(i).Y > testPoint.Y > (i+1).Y
(Note)
Both of the conditions indicate that a point is located within the edges of the Y-th coordinate
of the (i)-th and the (i+1)- th vertices of the polygon. If neither of the above
conditions is satisfied, then it is assured that a semi-infinite horizontal line draw
to the right from the testpoint will NOT cross the line that connects vertices i and i+1
of the polygon
*/
boolean withinYsEdges= belowLowY != belowHighY;
if( withinYsEdges){
// this is the slope of the line that connects vertices i and i+1 of the polygon
float slopeOfLine = ( vertX_j-vertX_i )/ (vertY_j-vertY_i) ;
// this looks up the x-coord of a point lying on the above line, given its y-coord
float pointOnLine = ( slopeOfLine* (testY - vertY_i) )+vertX_i;
//checks to see if x-coord of testPoint is smaller than the point on the line with the same y-coord
boolean isLeftToLine= testX < pointOnLine;
if(isLeftToLine){
//this statement changes true to false (and vice-versa)
locatedInPolygon= !locatedInPolygon;
}//end if (isLeftToLine)
}//end if (withinYsEdges
}
return locatedInPolygon;
}
Несколько слов об оптимизации: неверно, что самый короткий (и / или кратчайший) код реализован быстрее всего. Это гораздо более быстрый процесс чтения и сохранения элемента из массива и его использования (возможно) много раз в ходе выполнения блока кода, чем для доступа к массиву каждый раз, когда это требуется. Это особенно важно, если массив очень большой. На мой взгляд, сохраняя каждый член массива в хорошо названной переменной, также легче оценить его назначение и, таким образом, сформировать гораздо более читаемый код. Просто мои два цента...
Этот алгоритм работает в любом замкнутом многоугольнике, пока стороны многоугольника не пересекаются. Треугольник, пятиугольник, квадрат, даже очень изогнутая кусочно-линейная резинка, которая не пересекает себя.
1) Определите свой многоугольник как направленную группу векторов. Под этим подразумевается, что каждая сторона многоугольника описывается вектором, который идет от вершины an к вершине an+1. Векторы направлены так, что голова одного касается хвоста следующего, пока последний вектор не коснется хвоста первого.
2) Выберите точку для тестирования внутри или снаружи многоугольника.
3) Для каждого вектора Vn по периметру многоугольника найдите вектор Dn, который начинается в контрольной точке и заканчивается в хвосте Vn. Вычислите вектор Cn, определенный как DnXVn/DN*VN (X обозначает перекрестное произведение; * обозначает точечное произведение). Назовите величину Cn по имени Mn.
4) Добавьте все Mn и назовите это количество K.
5) Если K равно нулю, точка находится за пределами многоугольника.
6) Если K не равно нулю, точка находится внутри многоугольника.
Теоретически, точка, лежащая на краю многоугольника, даст неопределенный результат.
Геометрическое значение K - это общий угол, по которому блоха, сидящая на нашей контрольной точке, "видела" муравья, идущего по краю многоугольника, идущего влево, минус угол, пройденный вправо. В замкнутом контуре муравей заканчивается там, где он начался. За пределами многоугольника, независимо от местоположения, ответ ноль.
Внутри многоугольника, независимо от местоположения, ответ "один раз вокруг точки".
Алгоритм урезан до самых необходимых элементов. После того, как он был разработан и протестирован, все ненужные вещи были удалены. В результате вы не можете выполнить это легко, но это делает работу, а также в очень хорошем исполнении.
Я взял на себя смелость перевести его на ActionScript-3:
// not optimized yet (nvert could be left out)
public static function pnpoly(nvert: int, vertx: Array, verty: Array, x: Number, y: Number): Boolean
{
var i: int, j: int;
var c: Boolean = false;
for (i = 0, j = nvert - 1; i < nvert; j = i++)
{
if (((verty[i] > y) != (verty[j] > y)) && (x < (vertx[j] - vertx[i]) * (y - verty[i]) / (verty[j] - verty[i]) + vertx[i]))
c = !c;
}
return c;
}
Чтобы развернуть ответ @chowlette, где вторая строка проверяет, находится ли точка слева от линии, не дается деривация, но это то, что я разработал: во-первых, это помогает представить 2 основных случая:
- точка слева от линии
. /или же - точка находится справа от линии
/ .
Если бы наша точка была направлена на луч горизонтально, куда бы он попал на отрезок линии. Наша точка слева или справа? Внутри или снаружи? Мы знаем его координату y, потому что она по определению совпадает с точкой. Какой будет координата х?
Возьмите свою традиционную формулу линии y = mx + b, м - подъем над трассой. Здесь вместо этого мы пытаемся найти координату x точки на этом отрезке, который имеет ту же высоту (y), что и наша точка.
Итак, мы решаем для х: x = (y - b)/m, m это подъем над бегом, так что это становится бегом над подъемом или (yj - yi)/(xj - xi) становится (xj - xi)/(yj - yi), b - смещение от начала координат. Если мы предположим, yi в качестве базы для нашей системы координат b становится yi. Наша точка зрения testy это наш вклад, вычитая yi превращает всю формулу в смещение от yi,
Теперь у нас есть (xj - xi)/(yj - yi) или же 1/m раз у или (testy - yi): (xj - xi)(testy - yi)/(yj - yi) но testx не основан на yi поэтому мы добавим его обратно, чтобы сравнить два (или ноль testx)
Этот метод проверяет, разрезал ли луч от точки (testx, testy) до O (0,0) стороны многоугольника или нет.
Здесь есть общеизвестный вывод: если луч из 1 точки обрезает стороны многоугольника нечетное время, эта точка будет принадлежать многоугольнику, в противном случае эта точка будет находиться за пределами многоугольника.
Вот реализация php этого:
<?php
class Point2D {
public $x;
public $y;
function __construct($x, $y) {
$this->x = $x;
$this->y = $y;
}
function x() {
return $this->x;
}
function y() {
return $this->y;
}
}
class Point {
protected $vertices;
function __construct($vertices) {
$this->vertices = $vertices;
}
//Determines if the specified point is within the polygon.
function pointInPolygon($point) {
/* @var $point Point2D */
$poly_vertices = $this->vertices;
$num_of_vertices = count($poly_vertices);
$edge_error = 1.192092896e-07;
$r = false;
for ($i = 0, $j = $num_of_vertices - 1; $i < $num_of_vertices; $j = $i++) {
/* @var $current_vertex_i Point2D */
/* @var $current_vertex_j Point2D */
$current_vertex_i = $poly_vertices[$i];
$current_vertex_j = $poly_vertices[$j];
if (abs($current_vertex_i->y - $current_vertex_j->y) <= $edge_error && abs($current_vertex_j->y - $point->y) <= $edge_error && ($current_vertex_i->x >= $point->x) != ($current_vertex_j->x >= $point->x)) {
return true;
}
if ($current_vertex_i->y > $point->y != $current_vertex_j->y > $point->y) {
$c = ($current_vertex_j->x - $current_vertex_i->x) * ($point->y - $current_vertex_i->y) / ($current_vertex_j->y - $current_vertex_i->y) + $current_vertex_i->x;
if (abs($point->x - $c) <= $edge_error) {
return true;
}
if ($point->x < $c) {
$r = !$r;
}
}
}
return $r;
}
Тестовый забег:
<?php
$vertices = array();
array_push($vertices, new Point2D(120, 40));
array_push($vertices, new Point2D(260, 40));
array_push($vertices, new Point2D(45, 170));
array_push($vertices, new Point2D(335, 170));
array_push($vertices, new Point2D(120, 300));
array_push($vertices, new Point2D(260, 300));
$Point = new Point($vertices);
$point_to_find = new Point2D(190, 170);
$isPointInPolygon = $Point->pointInPolygon($point_to_find);
echo $isPointInPolygon;
var_dump($isPointInPolygon);
Это алгоритм, который я использую, но я добавил немного хитрости предварительной обработки, чтобы ускорить его. Мои полигоны имеют ~1000 ребер, и они не меняются, но мне нужно посмотреть, находится ли курсор внутри одного при каждом движении мыши.
Я в основном делю высоту ограничивающего прямоугольника на интервалы равной длины, и для каждого из этих интервалов я составляю список ребер, которые лежат внутри / пересекаются с ним.
Когда мне нужно найти точку, я могу вычислить - за O(1) время - в каком интервале она находится, и тогда мне нужно только проверить те ребра, которые есть в списке интервала.
Я использовал 256 интервалов, и это уменьшило количество ребер, которые мне нужно проверить, до 2-10 вместо ~1000.
Я думаю, что основная идея состоит в том, чтобы вычислять векторы из точки, по одному на ребро многоугольника. Если вектор пересекает одно ребро, то точка находится внутри многоугольника. Вогнутыми полигонами, если он пересекает нечетное число ребер, он тоже внутри (отказ от ответственности: хотя и не уверен, работает ли он для всех вогнутых полигонов).
В ответах выше используется деление с плавающей запятой для вычисления пересечения края и бесконечного луча, указывающего прямо из точки. Но точное расположение перекрестка на самом деле не нужно. Достаточно знать, находится ли пересечение слева или справа от ребра. Перекрестное произведение дает эту информацию, поэтому арифметика с плавающей запятой не требуется.
Часто на практике не имеет значения, каков будет результат, если точка находится на ребре. Пример такого «быстрого и грязного» целочисленного кода:
#include <vector>
struct Point
{
using T = int; // ..or float or double
T x, y;
};
/** Checks if a point is inside a simple polygon.
* The vertices are given as a coordinate vector.
* Note that this concise version does not give consistent
* answers for points which are on an edge.
*/
bool const pointInPolygon(const Point& p, const std::vector<Point>& vertex)
{
if (vertex.size() < 3)
return false;
bool isInside = false;
auto prev = vertex.size() - 1u;
for (decltype(prev) next = 0; next < vertex.size(); prev = next++)
if ((vertex[prev].y > p.y) != (vertex[next].y > p.y)) {
// Now an infinitely long ray stretching right from the point
// may intersect the edge between vertices prev and next.
auto diffVertexX = vertex[next].x - vertex[prev].x;
auto diffVertexY = vertex[next].y - vertex[prev].y;
auto diffPointX = p.x - vertex[prev].x;
auto diffPointY = p.y - vertex[prev].y;
// Use the sign of the cross product of the differences
// to determine the side on which the point lies compared to the edge.
// Calculate the two parts of the cross product separately.
auto crossProductXY = diffVertexX * diffPointY;
auto crossProductYX = diffVertexY * diffPointX;
// Switch the sign of the cross product depending on whether the edge
// is pointing upwards or downwards.
if (diffVertexY < 0)
isInside = isInside ^ (crossProductYX > crossProductXY);
else
isInside = isInside ^ (crossProductXY > crossProductYX);
}
return isInside;
}
Я изменил код, чтобы проверить, находится ли точка в многоугольнике, включая точку на ребре.
bool point_in_polygon_check_edge(const vec<double, 2>& v, vec<double, 2> polygon[], int point_count, double edge_error = 1.192092896e-07f)
{
const static int x = 0;
const static int y = 1;
int i, j;
bool r = false;
for (i = 0, j = point_count - 1; i < point_count; j = i++)
{
const vec<double, 2>& pi = polygon[i);
const vec<double, 2>& pj = polygon[j];
if (fabs(pi[y] - pj[y]) <= edge_error && fabs(pj[y] - v[y]) <= edge_error && (pi[x] >= v[x]) != (pj[x] >= v[x]))
{
return true;
}
if ((pi[y] > v[y]) != (pj[y] > v[y]))
{
double c = (pj[x] - pi[x]) * (v[y] - pi[y]) / (pj[y] - pi[y]) + pi[x];
if (fabs(v[x] - c) <= edge_error)
{
return true;
}
if (v[x] < c)
{
r = !r;
}
}
}
return r;
}