Algorithme de point en polygone

J’ai vu l’algorithme ci-dessous pour vérifier si un point se trouve dans un polygone donné à partir de ce lien :

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

J’ai essayé cet algorithme et ça marche vraiment parfait. Mais malheureusement, je ne peux pas bien le comprendre après avoir passé du temps à essayer d’en avoir l’idée.

Donc, si quelqu’un est capable de comprendre cet algorithme, expliquez-le-moi un peu.

Je vous remercie.

L’algorithme est la diffusion de rayons vers la droite. A chaque itération de la boucle, le sharepoint test est vérifié sur l’un des bords du polygone. La première ligne du test if réussit si la coordonnée y du point est dans la scope du bord. La deuxième ligne vérifie si le sharepoint test est à gauche de la ligne (je pense – je n’ai aucun papier à remettre pour vérifier). Si cela est vrai, la ligne tracée à droite du sharepoint test traverse ce bord.

En inversant à plusieurs resockets la valeur de c , l’algorithme compte combien de fois la ligne droite croise le polygone. S’il traverse un nombre impair de fois, alors le point est à l’intérieur; si un nombre pair, le point est à l’extérieur.

Je serais préoccupé par a) la précision de l’arithmétique à virgule flottante, et b) les effets d’avoir un bord horizontal, ou un sharepoint test avec la même coordonnée y comme un sumt, cependant.

Chowlett est correct dans tous les sens, dans toutes les formes et dans toutes les formes. L’algorithme suppose que si votre point est sur la ligne du polygone, alors c’est en dehors – dans certains cas, c’est faux. Changer les deux opérateurs ‘>’ en ‘> =’ et changer ‘<' en '<=' résoudra ce problème.

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

Cela pourrait être aussi détaillé que possible pour expliquer l’algorithme de traçage des rayons dans le code actuel. Il n’est peut-être pas optimisé mais cela doit toujours venir après une compréhension complète du système.

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

Juste un mot sur l'optimisation: ce n'est pas vrai que le code le plus court (et / ou le plus souvent) est le plus rapide. Il est beaucoup plus rapide de lire et de stocker un élément d'un tableau et de l'utiliser (éventuellement) plusieurs fois dans l'exécution du bloc de code plutôt que d'accéder au tableau chaque fois qu'il est requirejs. Ceci est particulièrement important si le tableau est extrêmement grand. À mon avis, en stockant chaque terme d'un tableau dans une variable bien nommée, il est également plus facile d'évaluer son objective et de former ainsi un code beaucoup plus lisible. Juste mes deux cents...

L’algorithme est réduit aux éléments les plus nécessaires. Après avoir été développé et testé, tous les éléments inutiles ont été supprimés. En conséquence, vous ne pouvez pas le comprendre facilement, mais il fait le travail et aussi dans de très bonnes performances.


J’ai pris la liberté de le traduire en 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; } 

Cet algorithme fonctionne dans tout polygone fermé tant que les côtés du polygone ne se croisent pas. Triangle, pentagone, carré, même un élastique très sinueux, linéaire par morceaux, qui ne se croise pas.

1) Définissez votre polygone en tant que groupe de vecteurs dirigé. On entend par là que chaque côté du polygone est décrit par un vecteur allant du sumt et du sumt an + 1. Les vecteurs sont ainsi dirigés de telle sorte que la tête de l’un touche la queue de l’autre jusqu’à ce que le dernier vecteur touche la queue du premier.

2) Sélectionnez le point à tester à l’intérieur ou à l’extérieur du polygone.

3) Pour chaque vecteur Vn sur le périmètre du polygone, trouver le vecteur Dn qui commence au sharepoint test et se termine à la queue de Vn. Calculez le vecteur Cn défini comme DnXVn / DN * VN (X indique un produit croisé; * indique un produit scalaire). Appelez la magnitude de Cn par le nom Mn.

4) Ajouter tous les Mn et appeler cette quantité K.

5) Si K est zéro, le point est en dehors du polygone.

6) Si K n’est pas nul, le point est à l’intérieur du polygone.

Théoriquement, un point situé sur le bord du polygone produira un résultat indéfini.

La signification géomésortingque de K est l’angle total que la puce assise sur notre sharepoint test “a vu” la fourmi qui marchait au bord du polygone et marchait à gauche, moins l’angle tourné vers la droite. En circuit fermé, la fourmi se termine là où elle a commencé. En dehors du polygone, quel que soit son emplacement, la réponse est zéro.
À l’intérieur du polygone, quel que soit son emplacement, la réponse est “une fois autour du point”.


Cette méthode vérifie si le rayon du point (testx, testy) à O (0,0) coupe ou non les côtés du polygone.

Il y a une conclusion bien connue ici : si un rayon de 1 pointe et coupe les côtés d’un polygone pendant un temps impair, ce point appartiendra au polygone, sinon ce point sera en dehors du polygone.

Je pense que l’idée de base est de calculer des vecteurs à partir du point, un par arête du polygone. Si le vecteur traverse un bord, le point se trouve dans le polygone. Par polygones concaves, si elle croise un nombre impair d’arêtes, elle est également à l’intérieur (avertissement: même si on ne sait pas si cela fonctionne pour tous les polygones concaves).

Pour élargir la réponse de @chowlette où la deuxième ligne vérifie si le point est à gauche de la ligne, aucune dérivation n’est donnée mais c’est ce que j’ai élaboré: D’abord, il est utile d’imaginer 2 cas fondamentaux:

  • le point est à gauche de la ligne . / . / ou
  • le point est à droite de la ligne / .

Si notre propos était de tirer un rayon à l’horizontale, où irait-il sur le segment de ligne. Notre point est-il à gauche ou à droite? À l’intérieur ou à l’extérieur? Nous connaissons sa coordonnée y car elle est par définition la même que le point. Quelle serait la coordonnée x?

Prenez votre formule de ligne traditionnelle y = mx + b . m est la montée sur la course. Ici, nous essayons plutôt de trouver la coordonnée x du point sur ce segment de ligne qui a la même hauteur (y) que notre point .

Nous résolvons donc pour x: x = (y - b)/m . m est une élévation par rapport à la course, de sorte que celle-ci devient une surélévation ou (yj - yi)/(xj - xi) devient (xj - xi)/(yj - yi) . b est le décalage d’origine. Si nous supposons que yi est la base de notre système de coordonnées, b devient yi. Notre point testy est notre entrée, soustraire yi transforme la formule entière en un décalage de yi .

Nous avons maintenant (xj - xi)/(yj - yi) ou 1/m fois y ou (testy - yi) : (xj - xi)(testy - yi)/(yj - yi) mais testx n’est pas basé sur yi nous l’ajoutons donc pour comparer les deux (ou zéro testx aussi)

Voici une implémentation php de ceci:

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

Essai:

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

C’est l’algorithme que j’utilise, mais j’ai ajouté quelques astuces de prétraitement pour l’accélérer. Mes polygones ont ~ 1000 arêtes et elles ne changent pas, mais je dois vérifier si le curseur se trouve à l’intérieur d’un mouvement de souris.

Je divise fondamentalement la hauteur du rectangle englobant en intervalles de longueur égale et pour chacun de ces intervalles, je comstack la liste des arêtes qui se trouvent dans / se croisant avec elle.

Quand j’ai besoin de rechercher un point, je peux calculer – en O (1) le temps – quel intervalle il se trouve et je n’ai plus qu’à tester les arêtes qui se trouvent dans la liste des intervalles.

J’ai utilisé 256 intervalles et cela a réduit le nombre de bords que j’ai besoin de tester à 2-10 au lieu de ~ 1000.

J’ai modifié le code pour vérifier si le point se trouve dans un polygone, y compris le point situé sur une arête.

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

J’ai changé le code original pour le rendre un peu plus lisible (ceci utilise également Eigen). L’algorithme est identique.

 // 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::MasortingxX2d &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; }