Comment puis-je déterminer si un point 2D se trouve dans un polygone ?
https://stackoverflow.com/questions/217578
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
J'essaie de créer un rapide Algorithme de point 2D à l'intérieur d'un polygone, à utiliser dans les tests de détection (par ex. Polygon.contains(p:Point)).Des suggestions de techniques efficaces seraient appréciées.
La solution
Pour les graphiques, je préférerais ne pas préférer les entiers. De nombreux systèmes utilisent des entiers pour peindre l'interface utilisateur (les pixels sont après tout ints), mais macOS, par exemple, utilise float pour tout. macOS ne connaît que les points et un point peut être converti en un pixel, mais selon la résolution de l'écran, il peut se traduire par autre chose. Sur les écrans de rétine, un demi-point (0,5 / 0,5) est un pixel. Malgré tout, je n'ai jamais remarqué que les interfaces utilisateur macOS sont nettement plus lentes que les autres. Toutes les API 3D (OpenGL ou Direct3D) fonctionnent également avec des flottants et les bibliothèques graphiques modernes tirent très souvent parti de l’accélération GPU.
Maintenant, vous avez dit que la vitesse était votre principale préoccupation. Très bien, passons à la vitesse. Avant d’exécuter un algorithme sophistiqué, effectuez d’abord un test simple. Créez un cadre de sélection aligné sur l’axe autour de votre polygone. Ceci est très facile, rapide et peut déjà vous protéger de nombreux calculs. Comment ça marche? Parcourez tous les points du polygone et recherchez les valeurs min / max de X et Y.
E.g. vous avez les points (9/1), (4/3), (2/7), (8/2), (3/6). Cela signifie que Xmin est égal à 2, Xmax à 9, Ymin à 1 et Ymax à 7. Un point situé à l'extérieur du rectangle avec les deux arêtes (2/1) et (9/7) ne peut pas figurer dans le polygone.
// p is your point, p.x is the x coord, p.y is the y coord
if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
// Definitely not within the polygon!
}
Il s'agit du premier test à exécuter pour un point quelconque. Comme vous pouvez le constater, ce test est ultra rapide mais il est également très grossier. Pour gérer les points situés dans le rectangle de délimitation, nous avons besoin d'un algorithme plus sophistiqué. Cela peut être calculé de différentes manières. Quelle méthode fonctionne également si le polygone peut avoir des trous ou sera toujours solide. Voici des exemples de solides (un convexe, un concave):
Et en voici un avec un trou:
Le vert a un trou au milieu!
L’algorithme le plus simple, capable de traiter les trois cas ci-dessus et toujours rapide, est appelé lancer de rayons . L'idée de l'algorithme est assez simple: tracez un rayon virtuel de n'importe où en dehors du polygone jusqu'à votre point et comptez la fréquence à laquelle il frappe un côté du polygone. Si le nombre de résultats est pair, c'est en dehors du polygone, s'il est impair, c'est à l'intérieur.
L'algorithme de nombre de tours serait une alternative, il est plus précis pour les points très proches d'une ligne de polygone, mais il est également beaucoup plus lent. La diffusion de rayons peut échouer pour des points situés trop près d'un polygone en raison de problèmes de précision de la virgule flottante et d'arrondi limités, mais en réalité, ce n'est pas un problème spectateur à reconnaître s’il est déjà à l’intérieur ou encore à l’extérieur.
Vous avez toujours le cadre de sélection ci-dessus, vous vous rappelez? Il suffit de choisir un point en dehors du cadre de sélection et de l’utiliser comme point de départ pour votre rayon. Par exemple. le point (Xmin - e/p.y) est bien en dehors du polygone.
Mais qu'est-ce que e? Eh bien, v1x1/v1y1 (en fait, epsilon) donne à la boîte de sélection un remplissage . Comme je l'ai dit, le lancer de rayon échoue si nous commençons trop près d'une ligne polygonale. Puisque le cadre de sélection peut être égal au polygone (si le polygone est un rectangle aligné par axe, le cadre de sélection est égal au polygone lui-même!), Nous avons besoin d'un peu de remplissage pour sécuriser ce point, c'est tout. Quelle taille devriez-vous choisir v1x2/v1y2? Pas si gros. Cela dépend de l'échelle du système de coordonnées que vous utilisez pour dessiner. Si votre largeur de pas en pixels est de 1,0, choisissez simplement 1,0 (mais 0,1 aurait également fonctionné)
Maintenant que nous avons le rayon avec ses coordonnées de début et de fin, le problème passe de & "; est le point dans le polygone &"; à & "; à quelle fréquence le rayon intersecte un côté du polygone &"; Par conséquent, nous ne pouvons pas simplement utiliser les points polygonaux comme auparavant. Nous avons maintenant besoin des côtés réels. UNEcôté est toujours défini par deux points.
side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:
Vous devez tester le rayon contre tous les côtés. Considérez le rayon comme un vecteur et chaque côté comme un vecteur. Le rayon doit frapper chaque côté exactement une fois ou jamais du tout. Il ne peut pas frapper deux fois le même côté. Deux lignes d'un espace 2D se coupent toujours exactement une fois, sauf si elles sont parallèles, auquel cas elles ne se croisent jamais. Cependant, comme les vecteurs ont une longueur limitée, il est possible que deux vecteurs ne soient pas parallèles et ne se croisent jamais car ils sont trop courts pour se rencontrer.
// Test the ray against all sides
int intersections = 0;
for (side = 0; side < numberOfSides; side++) {
// Test if current side intersects with ray.
// If yes, intersections++;
}
if ((intersections & 1) == 1) {
// Inside of polygon
} else {
// Outside of polygon
}
Jusqu'ici tout va bien, mais comment tester si deux vecteurs se croisent? Voici un code C (non testé), qui devrait faire l'affaire:
#define NO 0
#define YES 1
#define COLLINEAR 2
int areIntersecting(
float v1x1, float v1y1, float v1x2, float v1y2,
float v2x1, float v2y1, float v2x2, float v2y2
) {
float d1, d2;
float a1, a2, b1, b2, c1, c2;
// Convert vector 1 to a line (line 1) of infinite length.
// We want the line in linear equation standard form: A*x + B*y + C = 0
// See: http://en.wikipedia.org/wiki/Linear_equation
a1 = v1y2 - v1y1;
b1 = v1x1 - v1x2;
c1 = (v1x2 * v1y1) - (v1x1 * v1y2);
// Every point (x,y), that solves the equation above, is on the line,
// every point that does not solve it, is not. The equation will have a
// positive result if it is on one side of the line and a negative one
// if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector
// 2 into the equation above.
d1 = (a1 * v2x1) + (b1 * v2y1) + c1;
d2 = (a1 * v2x2) + (b1 * v2y2) + c1;
// If d1 and d2 both have the same sign, they are both on the same side
// of our line 1 and in that case no intersection is possible. Careful,
// 0 is a special case, that's why we don't test ">=" and "<=",
// but "<" and ">".
if (d1 > 0 && d2 > 0) return NO;
if (d1 < 0 && d2 < 0) return NO;
// The fact that vector 2 intersected the infinite line 1 above doesn't
// mean it also intersects the vector 1. Vector 1 is only a subset of that
// infinite line 1, so it may have intersected that line before the vector
// started or after it ended. To know for sure, we have to repeat the
// the same test the other way round. We start by calculating the
// infinite line 2 in linear equation standard form.
a2 = v2y2 - v2y1;
b2 = v2x1 - v2x2;
c2 = (v2x2 * v2y1) - (v2x1 * v2y2);
// Calculate d1 and d2 again, this time using points of vector 1.
d1 = (a2 * v1x1) + (b2 * v1y1) + c2;
d2 = (a2 * v1x2) + (b2 * v1y2) + c2;
// Again, if both have the same sign (and neither one is 0),
// no intersection is possible.
if (d1 > 0 && d2 > 0) return NO;
if (d1 < 0 && d2 < 0) return NO;
// If we get here, only two possibilities are left. Either the two
// vectors intersect in exactly one point or they are collinear, which
// means they intersect in any number of points from zero to infinite.
if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR;
// If they are not collinear, they must intersect in exactly one point.
return YES;
}
Les valeurs d'entrée sont les deux extrémités du vecteur 1 (v2x1/v2y1 et v2x2/v2y2) et du vecteur 2 (YES et NO). Donc, vous avez 2 vecteurs, 4 points, 8 coordonnées. == 0.0f et < epsilon sont clairs. <=> augmente les intersections, <=> ne fait rien.
Et COLLINEAR? Cela signifie que les deux vecteurs se trouvent sur la même ligne infinie. Selon leur position et leur longueur, ils ne se croisent pas du tout ou se croisent en un nombre infini de points. Je ne suis pas absolument sûr de savoir comment gérer ce cas, je ne le considérerais pas comme une intersection de toute façon. De toute façon, ce cas est plutôt rare dans la pratique en raison d’erreurs d’arrondi en virgule flottante; Un meilleur code ne serait probablement pas testé pour <=> mais plutôt pour quelque chose comme <=>, où epsilon est un nombre plutôt petit.
Si vous avez besoin de tester un plus grand nombre de points, vous pouvez certainement accélérer un peu le processus en conservant en mémoire les formes standard des équations linéaires des côtés des polygones, afin que vous n'ayez pas à les recalculer à chaque fois. Cela vous permettra d'économiser deux multiplications en virgule flottante et trois soustractions en virgule flottante à chaque test en échange du stockage de trois valeurs en virgule flottante par côté de polygone en mémoire. C’est un compromis typique entre mémoire et temps de calcul.
Dernier point mais non le moindre: si vous pouvez utiliser du matériel 3D pour résoudre le problème, il existe une alternative intéressante. Laissez simplement le GPU faire tout le travail pour vous. Créez une surface de peinture hors écran. Remplissez-le complètement avec la couleur noire. Laissez maintenant OpenGL ou Direct3D peindre votre polygone (ou même l’ensemble de vos polygones si vous souhaitez simplement tester si le point se trouve dans l’un d’entre eux, mais vous ne vous en souciez pas) et remplissez le ou les polygones avec un couleur, par exemple blanc. Pour vérifier si un point se trouve dans le polygone, obtenez la couleur de ce point à partir de la surface de dessin. Ceci est juste une extraction de mémoire O (1).
Bien sûr, cette méthode n’est utilisable que si votre surface de dessin n’est pas nécessairement immense. Si elle ne peut pas tenir dans la mémoire du processeur graphique, cette méthode est plus lente que sur le processeur. Si cela devait être énorme et que votre GPU prenne en charge les shaders modernes, vous pouvez toujours utiliser le GPU en implémentant la projection de rayons ci-dessus comme un GPU shader, ce qui est tout à fait possible. Pour un plus grand nombre de polygones ou un grand nombre de points à tester, cela rapportera, considérez que certains GPU pourront tester 64 à 256 points en parallèle. Notez toutefois que le transfert de données d'unité centrale à l'autre du processeur graphique est toujours coûteux. Par conséquent, pour tester quelques points sur quelques polygones simples, où les points ou les polygones sont dynamiques et changent fréquemment, une approche GPU sera rarement rentable. off.
Autres conseils
Je pense que le code suivant est la meilleure solution (extrait de ici ):
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;
}
Arguments
- Nvert : nombre de sommets du polygone. La répétition du premier sommet à la fin a été discutée dans l'article susmentionné.
- vertx, verty : tableaux contenant les coordonnées x et y des sommets du polygone.
- testx, testy : coordonnées X et y du point de test.
Court et efficace, il convient aux polygones convexes et concaves. Comme suggéré précédemment, vous devez d'abord vérifier le rectangle de délimitation et traiter les trous de polygone séparément.
L’idée derrière ceci est assez simple. L'auteur le décrit comme suit:
Je lance un rayon semi-infini horizontalement (en augmentant x, y fixe) à partir du point de test et compte le nombre d'arêtes qu'il traverse. A chaque croisement, le rayon bascule entre intérieur et extérieur. C'est ce qu'on appelle le théorème de la courbe de Jordan.
La variable c bascule de 0 à 1 et de 1 à 0 chaque fois que le rayon horizontal traverse un bord. Donc, fondamentalement, il faut savoir si le nombre d'arêtes croisées est pair ou impair. 0 signifie pair et 1 signifie impair.
Voici une version C # de la réponse donnée par nirg , qui provient de ce professeur RPI . Notez que l'utilisation du code de cette source RPI nécessite une attribution.
Un contrôle de boîte de sélection a été ajouté en haut. Cependant, comme le souligne James Brown, le code principal est presque aussi rapide que le contrôle de la boîte de sélection, de sorte que la vérification du cadre de sélection peut réellement ralentir le fonctionnement global, dans le cas où la plupart des points que vous vérifiez se trouvent à l'intérieur du cadre de sélection. . Vous pouvez donc laisser la case à cocher cochée ou vous pouvez également pré-calculer les boîtes de sélection de vos polygones si elles ne changent pas de forme trop souvent.
public bool IsPointInPolygon( Point p, Point[] polygon )
{
double minX = polygon[ 0 ].X;
double maxX = polygon[ 0 ].X;
double minY = polygon[ 0 ].Y;
double maxY = polygon[ 0 ].Y;
for ( int i = 1 ; i < polygon.Length ; i++ )
{
Point q = polygon[ i ];
minX = Math.Min( q.X, minX );
maxX = Math.Max( q.X, maxX );
minY = Math.Min( q.Y, minY );
maxY = Math.Max( q.Y, maxY );
}
if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > maxY )
{
return false;
}
// http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
bool inside = false;
for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ )
{
if ( ( polygon[ i ].Y > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
{
inside = !inside;
}
}
return inside;
}
Voici une variante JavaScript de la réponse de M. Katz basée sur l'approche de Nirg:
function pointIsInPoly(p, polygon) {
var isInside = false;
var minX = polygon[0].x, maxX = polygon[0].x;
var minY = polygon[0].y, maxY = polygon[0].y;
for (var n = 1; n < polygon.length; n++) {
var q = polygon[n];
minX = Math.min(q.x, minX);
maxX = Math.max(q.x, maxX);
minY = Math.min(q.y, minY);
maxY = Math.max(q.y, maxY);
}
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
return false;
}
var i = 0, j = polygon.length - 1;
for (i, j; i < polygon.length; j = i++) {
if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
isInside = !isInside;
}
}
return isInside;
}
Calculez la somme orientée des angles entre le point p et chacun des sommets du polygone.Si l’angle total orienté est de 360 degrés, le point est à l’intérieur.Si le total est 0, le point est à l'extérieur.
J'aime mieux cette méthode car elle est plus robuste et moins dépendante de la précision numérique.
Les méthodes qui calculent la régularité du nombre d'intersections sont limitées car vous pouvez « toucher » un sommet pendant le calcul du nombre d'intersections.
MODIFIER:À propos, cette méthode fonctionne avec des polygones concaves et convexes.
MODIFIER:J'ai récemment trouvé un tout Article Wikipédia sur le sujet.
L'article d'Eric Haines cité par bobobobo est vraiment excellent. Les tableaux comparant les performances des algorithmes sont particulièrement intéressants. la méthode de sommation d'angle est vraiment mauvaise comparée aux autres. Il est également intéressant de noter que les optimisations telles que l’utilisation d’une grille de recherche pour subdiviser davantage le polygone en & "; Dans &"; et " out " les secteurs peuvent rendre le test incroyablement rapide même sur des polygones avec > 1000 côtés.
Quoi qu’il en soit, c’est le début, mais mon vote va aux " croisements " méthode, qui est à peu près ce que Mecki décrit je pense. Cependant, je l’ai trouvé très succinctement décrit et codifié par David Bourke . J'aime le fait qu'il n'y ait pas de réelle trigonométrie requise, et cela fonctionne pour les convexes et les concaves, et il fonctionne assez bien lorsque le nombre de côtés augmente.
À propos, voici l'une des tables de performances de l'article d'intérêt d'Eric Haines, qui teste des polygones aléatoires.
number of edges per polygon
3 4 10 100 1000
MacMartin 2.9 3.2 5.9 50.6 485
Crossings 3.1 3.4 6.8 60.0 624
Triangle Fan+edge sort 1.1 1.8 6.5 77.6 787
Triangle Fan 1.2 2.1 7.3 85.4 865
Barycentric 2.1 3.8 13.8 160.7 1665
Angle Summation 56.2 70.4 153.6 1403.8 14693
Grid (100x100) 1.5 1.5 1.6 2.1 9.8
Grid (20x20) 1.7 1.7 1.9 5.7 42.2
Bins (100) 1.8 1.9 2.7 15.1 117
Bins (20) 2.1 2.2 3.7 26.3 278
Cette question est tellement intéressante. J'ai une autre idée réalisable différente des autres réponses de cet article.
L’idée est d’utiliser la somme des angles pour décider si la cible est à l’intérieur ou à l’extérieur. Si la cible est à l'intérieur d'une zone, la somme de la forme de l'angle par la cible et tous les deux points de bordure sera 360. Si la cible est à l'extérieur, la somme ne sera pas 360. L'angle a une direction. Si l'angle recule, l'angle est négatif. Cela revient à calculer le numéro sinueux .
Reportez-vous à cette image pour obtenir une compréhension de base de l'idée:
Mon algorithme suppose que le sens des aiguilles d’une montre est le sens des aiguilles d’une montre. Voici une entrée potentielle:
[[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]]
Voici le code python qui implémente l'idée:
def isInside(self, border, target):
degree = 0
for i in range(len(border) - 1):
a = border[i]
b = border[i + 1]
# calculate distance of vector
A = getDistance(a[0], a[1], b[0], b[1]);
B = getDistance(target[0], target[1], a[0], a[1])
C = getDistance(target[0], target[1], b[0], b[1])
# calculate direction of vector
ta_x = a[0] - target[0]
ta_y = a[1] - target[1]
tb_x = b[0] - target[0]
tb_y = b[1] - target[1]
cross = tb_y * ta_x - tb_x * ta_y
clockwise = cross < 0
# calculate sum of angles
if(clockwise):
degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
else:
degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
if(abs(round(degree) - 360) <= 3):
return True
return False
Version rapide du réponse par nirg :
extension CGPoint {
func isInsidePolygon(vertices: [CGPoint]) -> Bool {
guard !vertices.isEmpty else { return false }
var j = vertices.last!, c = false
for i in vertices {
let a = (i.y > y) != (j.y > y)
let b = (x < (j.x - i.x) * (y - i.y) / (j.y - i.y) + i.x)
if a && b { c = !c }
j = i
}
return c
}
}
J'ai travaillé sur ce sujet à l'époque où j'étais chercheur sous Michael Stonebraker - vous savoir, le professeur qui a proposé Ingres , PostgreSQL , etc.,
Nous avons compris que le moyen le plus rapide était de commencer par créer un cadre de sélection, car c Si c'est en dehors de la boîte englobante, c'est en dehors. Sinon, vous faites le travail le plus dur ...
Si vous voulez un excellent algorithme, regardez le code source du projet Open Source PostgreSQL pour le travail géographique ...
Je tiens à signaler que nous n’avons jamais eu une idée précise du droit vs droitier (également exprimable en tant que & "; intérieur &" contre & "extérieur &"; problème .. .
MISE À JOUR
Le lien de BKB fournissait un bon nombre d’algorithmes raisonnables. Je travaillais sur des problèmes liés aux sciences de la Terre et j’avais donc besoin d’une solution qui fonctionne en latitude / longitude et qui présente le problème particulier de l’adhérence: la zone est-elle située dans la zone la plus petite ou dans la plus grande? La réponse est que la & Quot; direction & Quot; des vertices est important - que ce soit pour gaucher ou droitier - de cette façon, vous pouvez indiquer l’une des zones comme & "à l’intérieur de &"; un polygone donné. En tant que tel, mon travail utilisait la solution trois énumérée sur cette page.
De plus, mon travail utilisait des fonctions distinctes pour & "; sur la ligne &"; tests.
... Puisque quelqu'un a demandé: nous avons compris que les tests de la boîte englobante étaient préférables lorsque le nombre de vertices dépassait un certain nombre - faites un test très rapide avant de faire le test plus long si nécessaire ... Une boîte englobante est créée par en prenant simplement le plus grand x, le plus petit x, le plus grand et le plus petit y et en les assemblant pour former quatre points d'une boîte ...
Autre astuce pour ceux qui suivent: nous avons mis au point notre plus sophistiqué et & "obscurcissement de la lumière &"; calculer dans un espace de grille, le tout en points positifs sur un plan, puis re-projeté dans & "Real &"; longitude / latitude, évitant ainsi les erreurs possibles d’embobinage lorsqu’on croise la ligne 180 de longitude et lors de la manipulation de régions polaires. A bien fonctionné!
J'aime beaucoup la solution publiée par Nirg et modifiée par bobobobo. Je viens de le rendre javascript amical et un peu plus lisible pour mon utilisation:
function insidePoly(poly, pointx, pointy) {
var i, j;
var inside = false;
for (i = 0, j = poly.length - 1; i < poly.length; j = i++) {
if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside;
}
return inside;
}
La réponse de David Segond est à peu près la réponse générale standard, et celle de Richard T est l’optimisation la plus courante, bien qu’il en existe d’autres. D'autres optimisations fortes reposent sur des solutions moins générales. Par exemple, si vous allez vérifier le même polygone avec beaucoup de points, sa triangulation peut considérablement accélérer les choses car il existe un certain nombre d'algorithmes de recherche TIN très rapides. Une autre est que si le polygone et les points sont sur un plan limité à basse résolution, disons un affichage à l'écran, vous pouvez peindre le polygone sur une mémoire tampon mappée en mémoire dans une couleur donnée et vérifier la couleur d'un pixel donné pour voir si elle se trouve dans les polygones.
Comme de nombreuses optimisations, celles-ci sont basées sur des cas spécifiques plutôt que sur des cas généraux et génèrent des avantages basés sur le temps amorti plutôt que sur une utilisation unique.
Travaillant dans ce domaine, j’ai trouvé la "Géométrie de calcul en C" de ISBN 0-521-44034-3 de Joeseph O'Rourkes très utile.
La solution la plus simple serait de diviser le polygone en triangles et de tester les triangles, comme expliqué ici
Si votre polygone est CONVEX , il pourrait y avoir une meilleure approche. Regardez le polygone comme une collection de lignes infinies. Chaque ligne divisant l'espace en deux. pour chaque point, il est facile de dire si c'est d'un côté ou de l'autre côté de la ligne. Si un point se trouve du même côté de toutes les lignes, il se trouve à l'intérieur du polygone.
Je réalise que c'est vieux, mais voici un algorithme de casting de rayons implémenté dans Cocoa, au cas où quelqu'un serait intéressé. Pas sûr que ce soit le moyen le plus efficace de faire les choses, mais cela peut aider quelqu'un.
- (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point
{
NSBezierPath *currentPath = [path bezierPathByFlatteningPath];
BOOL result;
float aggregateX = 0; //I use these to calculate the centroid of the shape
float aggregateY = 0;
NSPoint firstPoint[1];
[currentPath elementAtIndex:0 associatedPoints:firstPoint];
float olderX = firstPoint[0].x;
float olderY = firstPoint[0].y;
NSPoint interPoint;
int noOfIntersections = 0;
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
aggregateX += points[0].x;
aggregateY += points[0].y;
}
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
//line equations in Ax + By = C form
float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y;
float _B_FOO = point.x - (aggregateX/[currentPath elementCount]);
float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y);
float _A_BAR = olderY - points[0].y;
float _B_BAR = points[0].x - olderX;
float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY);
float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO);
if (det != 0) {
//intersection points with the edges
float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det;
float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det;
interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint);
if (olderX <= points[0].x) {
//doesn't matter in which direction the ray goes, so I send it right-ward.
if ((interPoint.x >= olderX && interPoint.x <= points[0].x) && (interPoint.x > point.x)) {
noOfIntersections++;
}
} else {
if ((interPoint.x >= points[0].x && interPoint.x <= olderX) && (interPoint.x > point.x)) {
noOfIntersections++;
}
}
}
olderX = points[0].x;
olderY = points[0].y;
}
if (noOfIntersections % 2 == 0) {
result = FALSE;
} else {
result = TRUE;
}
return result;
}
Version obj-C de la réponse de nirg avec exemple de méthode pour tester des points. La réponse de Nirg a bien fonctionné pour moi.
- (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test {
NSUInteger nvert = [vertices count];
NSInteger i, j, c = 0;
CGPoint verti, vertj;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue];
vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue];
if (( (verti.y > test.y) != (vertj.y > test.y) ) &&
( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) )
c = !c;
}
return (c ? YES : NO);
}
- (void)testPoint {
NSArray *polygonVertices = [NSArray arrayWithObjects:
[NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)],
[NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)],
[NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)],
[NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)],
nil
];
CGPoint tappedPoint = CGPointMake(23.0, 70.0);
if ([self isPointInPolygon:polygonVertices point:tappedPoint]) {
NSLog(@"YES");
} else {
NSLog(@"NO");
}
}
La version C # de la réponse de nirg est ici: je partagerai simplement le code. Cela peut faire gagner du temps à quelqu'un.
public static bool IsPointInPolygon(IList<Point> polygon, Point testPoint) {
bool result = false;
int j = polygon.Count() - 1;
for (int i = 0; i < polygon.Count(); i++) {
if (polygon[i].Y < testPoint.Y && polygon[j].Y >= testPoint.Y || polygon[j].Y < testPoint.Y && polygon[i].Y >= testPoint.Y) {
if (polygon[i].X + (testPoint.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) * (polygon[j].X - polygon[i].X) < testPoint.X) {
result = !result;
}
}
j = i;
}
return result;
}
Il n’ya rien de plus beau qu'une définition inductive d’un problème. Par souci d’exhaustivité, vous avez ici une version dans prolog qui pourrait également clarifier les points cachés derrière la diffusion de rayons :
Basé sur l'algorithme de simulation de la simplicité dans http: // www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
Quelques prédicats auxiliaires:
exor(A,B):- \+A,B;A,\+B.
in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)).
inside(false).
inside(_,[_|[]]).
inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) + X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]).
get_line(_,_,[]).
get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]).
L’équation d’une ligne à 2 points A et B (ligne (A, B)) est la suivante:
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
Il est important que le sens de rotation de la ligne soit réglé dans le sens des aiguilles d'une montre pour les limites et dans le sens inverse des aiguilles d'une montre pour les trous. Nous allons vérifier si le point (X, Y), c’est-à-dire que le point testé est à gauche demi-plan de notre ligne (c'est une question de goût, ça pourrait aussi être du côté droit, mais aussi la direction des lignes de démarcation doit être modifiée en dans ce cas), c’est projeter le rayon du point vers la droite (ou la gauche) et reconnaître l'intersection avec la ligne. Nous avons choisi de projeter le rayon dans le sens horizontal (encore une fois, c'est une question de goût, cela pourrait également se faire verticalement avec des restrictions similaires), nous avons donc:
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
Maintenant, nous devons savoir si le point est à gauche (ou à droite) de le segment de ligne seulement, pas l’ensemble de l’avion, nous devons donc restreindre la recherche uniquement à ce segment, mais cela est facile car être à l'intérieur du segment, un seul point dans la ligne peut être plus élevé que Y dans l'axe vertical. Comme il s'agit d'une restriction plus forte, doit être le premier à vérifier, donc nous prenons d'abord seulement les lignes répondre à cette exigence et ensuite vérifier sa possition. Par la Jordanie Théorème de courbe, tout rayon projeté sur un polygone doit se croiser à un angle nombre pair de lignes. Donc, nous avons fini, nous allons jeter le rayon au à droite et à chaque fois qu’il coupe une ligne, change son état. Cependant, dans notre implémentation, nous allons vérifier la longueur de le sac de solutions répondant aux restrictions données et décide de la intendance sur elle. cela doit être fait pour chaque ligne du polygone.
is_left_half_plane(_,[],[],_).
is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] = [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA));
is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test).
in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon), in_range(Y,YA,YB).
all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line), in_y_range_at_poly(Coordinate,Line,Polygon), Lines).
traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, <), IntersectingLines), length(IntersectingLines, Count).
% This is the entry point predicate
inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside).
Version Java:
public class Geocode {
private float latitude;
private float longitude;
public Geocode() {
}
public Geocode(float latitude, float longitude) {
this.latitude = latitude;
this.longitude = longitude;
}
public float getLatitude() {
return latitude;
}
public void setLatitude(float latitude) {
this.latitude = latitude;
}
public float getLongitude() {
return longitude;
}
public void setLongitude(float longitude) {
this.longitude = longitude;
}
}
public class GeoPolygon {
private ArrayList<Geocode> points;
public GeoPolygon() {
this.points = new ArrayList<Geocode>();
}
public GeoPolygon(ArrayList<Geocode> points) {
this.points = points;
}
public GeoPolygon add(Geocode geo) {
points.add(geo);
return this;
}
public boolean inside(Geocode geo) {
int i, j;
boolean c = false;
for (i = 0, j = points.size() - 1; i < points.size(); j = i++) {
if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) &&
(geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude()))
c = !c;
}
return c;
}
}
Port .Net:
static void Main(string[] args)
{
Console.Write("Hola");
List<double> vertx = new List<double>();
List<double> verty = new List<double>();
int i, j, c = 0;
vertx.Add(1);
vertx.Add(2);
vertx.Add(1);
vertx.Add(4);
vertx.Add(4);
vertx.Add(1);
verty.Add(1);
verty.Add(2);
verty.Add(4);
verty.Add(4);
verty.Add(1);
verty.Add(1);
int nvert = 6; //Vértices del poligono
double testx = 2;
double testy = 5;
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 = 1;
}
}
VERSION VBA:
Remarque: rappelez-vous que si votre polygone est une zone de la carte, Latitude / Longitude sont des valeurs Y / X et non X / Y (Latitude = Y, Longitude = X) en raison de ce que je comprends comme étant des implications historiques provenant de retour lorsque la longitude n’était pas une mesure.
MODULE DE CLASSE: CPoint
Private pXValue As Double
Private pYValue As Double
'''''X Value Property'''''
Public Property Get X() As Double
X = pXValue
End Property
Public Property Let X(Value As Double)
pXValue = Value
End Property
'''''Y Value Property'''''
Public Property Get Y() As Double
Y = pYValue
End Property
Public Property Let Y(Value As Double)
pYValue = Value
End Property
MODULE:
Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean
Dim i As Integer
Dim j As Integer
Dim q As Object
Dim minX As Double
Dim maxX As Double
Dim minY As Double
Dim maxY As Double
minX = polygon(0).X
maxX = polygon(0).X
minY = polygon(0).Y
maxY = polygon(0).Y
For i = 1 To UBound(polygon)
Set q = polygon(i)
minX = vbMin(q.X, minX)
maxX = vbMax(q.X, maxX)
minY = vbMin(q.Y, minY)
maxY = vbMax(q.Y, maxY)
Next i
If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > maxY Then
isPointInPolygon = False
Exit Function
End If
' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
isPointInPolygon = False
i = 0
j = UBound(polygon)
Do While i < UBound(polygon) + 1
If (polygon(i).Y > p.Y) Then
If (polygon(j).Y < p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
ElseIf (polygon(i).Y < p.Y) Then
If (polygon(j).Y > p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
End If
j = i
i = i + 1
Loop
End Function
Function vbMax(n1, n2) As Double
vbMax = IIf(n1 > n2, n1, n2)
End Function
Function vbMin(n1, n2) As Double
vbMin = IIf(n1 > n2, n2, n1)
End Function
Sub TestPointInPolygon()
Dim i As Integer
Dim InPolygon As Boolean
' MARKER Object
Dim p As CPoint
Set p = New CPoint
p.X = <ENTER X VALUE HERE>
p.Y = <ENTER Y VALUE HERE>
' POLYGON OBJECT
Dim polygon() As CPoint
ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1
For i = 0 To <ENTER VALUE HERE> 'Same value as above
Set polygon(i) = New CPoint
polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through
polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through
Next i
InPolygon = isPointInPolygon(p, polygon)
MsgBox InPolygon
End Sub
J'ai réalisé une implémentation Python de de nirg's c ++ code :
Entrées
- bounding_points: les nœuds constituant le polygone.
-
bounding_box_positions: candidats à filtrer. (Dans mon implémentation créée à partir du cadre de sélection.
(Les entrées sont des listes de n-uplets au format:
[(xcord, ycord), ...])
Retourne
- Tous les points situés à l'intérieur du polygone.
def polygon_ray_casting(self, bounding_points, bounding_box_positions):
# Arrays containing the x- and y-coordinates of the polygon's vertices.
vertx = [point[0] for point in bounding_points]
verty = [point[1] for point in bounding_points]
# Number of vertices in the polygon
nvert = len(bounding_points)
# Points that are inside
points_inside = []
# For every candidate position within the bounding box
for idx, pos in enumerate(bounding_box_positions):
testx, testy = (pos[0], pos[1])
c = 0
for i in range(0, nvert):
j = i - 1 if i != 0 else nvert - 1
if( ((verty[i] > testy ) != (verty[j] > testy)) and
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ):
c += 1
# If odd, that means that we are inside the polygon
if c % 2 == 1:
points_inside.append(pos)
return points_inside
Encore une fois, l’idée vient de ici
Voici un point du test polygonal en C qui n’utilise pas le lancer de rayons. Et cela peut fonctionner pour les zones qui se chevauchent (auto-intersections), voir l'argument use_holes.
/* math lib (defined below) */
static float dot_v2v2(const float a[2], const float b[2]);
static float angle_signed_v2v2(const float v1[2], const float v2[2]);
static void copy_v2_v2(float r[2], const float a[2]);
/* intersection function */
bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr,
const bool use_holes)
{
/* we do the angle rule, define that all added angles should be about zero or (2 * PI) */
float angletot = 0.0;
float fp1[2], fp2[2];
unsigned int i;
const float *p1, *p2;
p1 = verts[nr - 1];
/* first vector */
fp1[0] = p1[0] - pt[0];
fp1[1] = p1[1] - pt[1];
for (i = 0; i < nr; i++) {
p2 = verts[i];
/* second vector */
fp2[0] = p2[0] - pt[0];
fp2[1] = p2[1] - pt[1];
/* dot and angle and cross */
angletot += angle_signed_v2v2(fp1, fp2);
/* circulate */
copy_v2_v2(fp1, fp2);
p1 = p2;
}
angletot = fabsf(angletot);
if (use_holes) {
const float nested = floorf((angletot / (float)(M_PI * 2.0)) + 0.00001f);
angletot -= nested * (float)(M_PI * 2.0);
return (angletot > 4.0f) != ((int)nested % 2);
}
else {
return (angletot > 4.0f);
}
}
/* math lib */
static float dot_v2v2(const float a[2], const float b[2])
{
return a[0] * b[0] + a[1] * b[1];
}
static float angle_signed_v2v2(const float v1[2], const float v2[2])
{
const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]);
return atan2f(perp_dot, dot_v2v2(v1, v2));
}
static void copy_v2_v2(float r[2], const float a[2])
{
r[0] = a[0];
r[1] = a[1];
}
Remarque: c’est l’une des méthodes les moins optimales car elle inclut de nombreux appels à atan2f, mais elle peut également intéresser les développeurs qui lisent ce fil (dans mes tests, il ralentit environ 23 fois puis utilise la méthode d’intersection de lignes. ).
Pour détecter un polygone, nous devons tester deux choses:
- Si Point est dans la surface du polygone. (peut être accompli par l'algorithme Ray-Casting)
- Si Point se trouve sur la frontière du polygone (peut être réalisé avec le même algorithme que celui utilisé pour la détection de point sur une polyligne (ligne)).
Pour traiter les cas spéciaux suivants dans Algorithme de diffusion de rayon :
- Le rayon recouvre l'un des côtés du polygone.
- Le point est à l'intérieur du polygone et le rayon passe par un sommet du polygone.
- Le point est en dehors du polygone et le rayon ne fait que toucher l'un des angles du polygone.
Vérifiez Déterminer si un point est à l'intérieur d'un polygone complexe . Cet article fournit un moyen simple de les résoudre afin qu'aucun traitement spécial ne soit nécessaire pour les cas ci-dessus.
Vous pouvez le faire en vérifiant si l'aire formée en reliant le point souhaité aux sommets de votre polygone correspond à l'aire du polygone lui-même.
Vous pouvez également vérifier si la somme des angles intérieurs de votre point à chaque paire de deux sommets de polygone consécutifs à votre point de contrôle correspond à 360, mais j’ai le sentiment que la première option est plus rapide car elle n’implique pas divisions ni calculs d'inverse de fonctions trigonométriques.
Je ne sais pas ce qui se passe si votre polygone a un trou à l'intérieur, mais il me semble que l'idée principale peut être adaptée à cette situation
Vous pouvez également poster la question dans une communauté mathématique. Je parie qu'ils ont un million de façons de le faire
Si vous recherchez une bibliothèque de scripts java, une extension javascript de Google Maps v3 permet à la classe Polygon de détecter si un point y réside ou non.
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
Lors de l'utilisation de qt (Qt 4.3+) , on peut utiliser la fonction QPolygon containsPoint
La réponse dépend de si vous avez des polygones simples ou complexes. Les polygones simples ne doivent avoir aucune intersection de segment de droite. Ainsi, ils peuvent avoir des trous mais les lignes ne peuvent pas se croiser. Les régions complexes peuvent avoir des intersections de lignes - elles peuvent donc avoir les régions qui se chevauchent ou les régions qui se touchent simplement par un seul point.
Pour les polygones simples, le meilleur algorithme est l’algorithme Ray casting (Crossing number). Pour les polygones complexes, cet algorithme ne détecte pas les points situés dans les régions qui se chevauchent. Donc, pour les polygones complexes, vous devez utiliser l’algorithme du nombre de Winding.
Voici un excellent article avec une implémentation en C des deux algorithmes. Je les ai essayées et elles fonctionnent bien.
Version Scala de la solution par nirg (en supposant que le contrôle préliminaire du rectangle est effectué séparément):
def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = {
val length = polygon.length
@tailrec
def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = {
if (i == length)
tracker
else {
val intersects = (polygon(i).y > p.y) != (polygon(j).y > p.y) && p.x < (polygon(j).x - polygon(i).x) * (p.y - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x
oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
}
}
oddIntersections(0, length - 1, tracker = false)
}
Voici la version golang de @nirg answer (inspirée du code C # de @@ m-katz)
func isPointInPolygon(polygon []point, testp point) bool {
minX := polygon[0].X
maxX := polygon[0].X
minY := polygon[0].Y
maxY := polygon[0].Y
for _, p := range polygon {
minX = min(p.X, minX)
maxX = max(p.X, maxX)
minY = min(p.Y, minY)
maxY = max(p.Y, maxY)
}
if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
return false
}
inside := false
j := len(polygon) - 1
for i := 0; i < len(polygon); i++ {
if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
inside = !inside
}
j = i
}
return inside
}
Surpris, personne n’avait évoqué cette question plus tôt, mais pour les pragmatiques exigeant une base de données: MongoDB offre un excellent support pour les requêtes Geo, y compris celle-ci.
Ce que vous recherchez, c'est:
db.neighborhoods.findOne ({geometry: {$ geoIntersects: {$ geometry: { tapez: & "Point &"; coordonnées: [& "longitude &"; & "latitude &"; ]}}} })
Neighborhoods est la collection qui stocke un ou plusieurs polygones au format GeoJson standard. Si la requête renvoie la valeur null, elle n'est pas intersectée, sinon.
Très bien documenté ici: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/ <
La performance pour plus de 6 000 points classés dans une grille de 330 polygones irréguliers était inférieure à une minute, sans aucune optimisation et incluant le temps de mise à jour des documents avec leur polygone respectif.
