Point in Polygon Algorithm

Vi que el algoritmo de abajo funciona para verificar si un punto está en un polígono dado desde este enlace :

int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy) { int i, j, c = 0; for (i = 0, j = nvert-1; i testy) != (verty[j]>testy)) && (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) ) c = !c; } return c; } 

Intenté este algoritmo y realmente funciona perfecto. Pero lamentablemente no puedo entenderlo bien después de pasar algún tiempo tratando de entenderlo.

Entonces, si alguien puede entender este algoritmo, explíquemelo un poco.

Gracias.

El algoritmo está lanzando rayos a la derecha. Cada iteración del ciclo, el punto de prueba se compara con uno de los bordes del polígono. La primera línea de la prueba if tiene éxito si el punto y-coord está dentro del scope del borde. La segunda línea verifica si el punto de prueba está a la izquierda de la línea (creo, no tengo ningún papel a mano para verificar). Si eso es cierto, la línea dibujada hacia la derecha desde el punto de prueba cruza ese borde.

Al invertir repetidamente el valor de c , el algoritmo cuenta cuántas veces la línea a la derecha cruza el polígono. Si cruza un número impar de veces, entonces el punto está adentro; si es un número par, el punto está afuera.

Me preocuparía a) la precisión de la aritmética de coma flotante, yb) los efectos de tener un borde horizontal, o un punto de prueba con la misma y-coord como un vértice, sin embargo.

Chowlett tiene razón en todos los sentidos, formas y formas. El algoritmo supone que si su punto está en la línea del polígono, entonces está afuera; en algunos casos, esto es falso. Cambiar los dos operadores ‘>’ a ‘> =’ y cambiar ‘<' a '<=' lo arreglará.

 bool PointInPolygon(Point point, Polygon polygon) { vector 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; } 

Esto podría ser tan detallado como podría ser para explicar el algoritmo de trazado de rayos en el código real. Puede que no se optimice, pero eso siempre debe venir después de una comprensión completa del sistema.

  //method to check if a Coordinate is located in a polygon public boolean checkIsInPolygon(ArrayList 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; } 

Solo una palabra sobre la optimización: no es cierto que el código más corto (y / o el más corto) sea el más rápido implementado. Es un proceso mucho más rápido leer y almacenar un elemento de una matriz y usarlo (posiblemente) muchas veces en la ejecución del bloque de código que acceder a la matriz cada vez que se requiere. Esto es especialmente significativo si la matriz es extremadamente grande. En mi opinión, al almacenar cada término de una matriz en una variable bien nombrada, también es más fácil evaluar su propósito y así formar un código mucho más legible. Solo mis dos centavos ...

El algoritmo se reduce a los elementos más necesarios. Después de que fue desarrollado y probado, todas las cosas innecesarias han sido eliminadas. Como resultado, no puede entenderlo fácilmente, pero cumple su función y también en muy buen rendimiento.


Me tomé la libertad de traducirlo a 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; } 

Este algoritmo funciona en cualquier polígono cerrado siempre que los lados del polígono no se crucen. Triángulo, pentágono, cuadrado, incluso una goma elástica lineal muy discontinua que no se cruza.

1) Defina su polígono como un grupo dirigido de vectores. Con esto se quiere decir que cada lado del polígono está descrito por un vector que va del vértice un al vértice un + 1. Los vectores están dirigidos de tal manera que la cabeza de uno toca la cola de la siguiente hasta que el último vector toca la cola de la primera.

2) Seleccione el punto para probar dentro o fuera del polígono.

3) Para cada vector Vn a lo largo del perímetro del polígono, encuentre el vector Dn que comienza en el punto de prueba y termina en la cola de Vn. Calcule el vector Cn definido como DnXVn / DN * VN (X indica producto cruzado; * indica producto punto). Llame la magnitud de Cn por el nombre Mn.

4) Agregue todo Mn y llame a esta cantidad K.

5) Si K es cero, el punto está fuera del polígono.

6) Si K no es cero, el punto está dentro del polígono.

Teóricamente, un punto que se encuentra EN EL borde del polígono producirá un resultado indefinido.

El significado geométrico de K es el ángulo total que la pulga que está sentada en nuestro punto de prueba “vio” a la ant caminando en el borde del polígono caminando hacia la izquierda menos el ángulo caminó hacia la derecha. En un circuito cerrado, la ant termina donde comenzó. Fuera del polígono, independientemente de la ubicación, la respuesta es cero.
Dentro del polígono, independientemente de su ubicación, la respuesta es “una vez alrededor del punto”.


Este método verifica si el rayo del punto (testx, irritable) a O (0,0) corta los lados del polígono o no.

Aquí hay una conclusión bien conocida: si un rayo de 1 punto corta los lados de un polígono por un tiempo impar, ese punto pertenecerá al polígono, de lo contrario ese punto estará fuera del polígono.

Creo que la idea básica es calcular vectores desde el punto, uno por borde del polígono. Si el vector cruza un borde, entonces el punto está dentro del polígono. Por polígonos cóncavos si cruza un número impar de bordes, también está dentro (descargo de responsabilidad: aunque no estoy seguro de si funciona para todos los polígonos cóncavos).

Para expandir la respuesta de @ chowlette donde la segunda línea verifica si el punto está a la izquierda de la línea, no se proporciona ninguna derivación, pero esto es lo que resolví: primero ayuda a imaginar 2 casos básicos:

  • el punto queda de la línea . / . / o
  • el punto es correcto de la línea / .

Si nuestro punto fuera disparar un rayo horizontalmente, ¿dónde golpearía el segmento de línea? ¿Está nuestro punto a la izquierda o a la derecha de eso? Por dentro o por fuera? Sabemos su coordenada y porque es, por definición, lo mismo que el punto. ¿Cuál sería la coordenada x?

Tome su fórmula de línea tradicional y = mx + b . m es el aumento sobre la carrera. Aquí, en cambio, estamos tratando de encontrar la coordenada x del punto en ese segmento de línea que tiene la misma altura (y) que nuestro punto .

Entonces resolvemos para x: x = (y - b)/m . m es la subida sobre la carrera, por lo que se convierte en carrera sobre la elevación o (yj - yi)/(xj - xi) convierte en (xj - xi)/(yj - yi) . b es el desplazamiento desde el origen. Si asumimos yi como la base de nuestro sistema de coordenadas, b se convierte en yi. Nuestro punto testy es nuestra entrada, restando yi convierte toda la fórmula en un desplazamiento de yi .

Ahora tenemos (xj - xi)/(yj - yi) o 1/m veces y o (testy - yi) : (xj - xi)(testy - yi)/(yj - yi) pero testx no está basado en yi así que lo agregamos de nuevo para comparar los dos (o cero testx también)

Aquí hay una implementación php de esto:

 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; } 

Prueba de funcionamiento:

  pointInPolygon($point_to_find); echo $isPointInPolygon; var_dump($isPointInPolygon); 

Este es el algoritmo que uso, pero agregué un poco de trucos de preprocesamiento para acelerarlo. Mis polígonos tienen ~ 1000 bordes y no cambian, pero necesito buscar si el cursor está dentro de uno en cada movimiento del mouse.

Básicamente dividí la altura del rectángulo delimitador en intervalos de igual longitud y para cada uno de estos intervalos compilo la lista de bordes que se encuentran dentro / se cruzan con él.

Cuando necesito buscar un punto, puedo calcular, en O (1) tiempo, en qué intervalo se encuentra y luego solo necesito probar los bordes que están en la lista del intervalo.

Usé 256 intervalos y esto redujo la cantidad de aristas que necesito probar a 2-10 en lugar de ~ 1000.

Modifiqué el código para verificar si el punto está en un polígono, incluido el punto en un borde.

 bool point_in_polygon_check_edge(const vec& v, vec 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& pi = polygon[i); const vec& 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; } 

Cambié el código original para hacerlo un poco más legible (también esto usa Eigen). El algoritmo es idéntico.

 // 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; }