Question Comment puis-je déterminer si un point 2D est dans un polygone?


J'essaye de créer un vite 2D point à l'intérieur de l'algorithme de polygone, à utiliser dans le test de hit (par ex. Polygon.contains(p:Point)). Des suggestions de techniques efficaces seraient appréciées.


426
2017-10-20 05:20


origine


Réponses:


Pour les graphiques, je préfère ne pas préférer les nombres entiers. De nombreux systèmes utilisent des entiers pour la peinture de l'interface utilisateur (les pixels sont ints après tout), mais macOS par exemple utilise float pour tout. macOS ne connait que des points et un point peut être traduit en un pixel, mais en fonction de la résolution du moniteur, il peut traduire quelque chose d'autre. Sur les écrans de rétine, un demi-point (0,5 / 0,5) est un pixel. Pourtant, je n'ai jamais remarqué que les interfaces utilisateur MacOS sont nettement plus lentes que les autres interfaces utilisateur. Après toutes les API 3D (OpenGL ou Direct3D) fonctionnent également avec des flotteurs et les bibliothèques graphiques modernes profitent très souvent de l'accélération GPU.

Maintenant, vous avez dit que la vitesse est votre principale préoccupation, d'accord, passons à la vitesse. Avant d'exécuter un algorithme sophistiqué, commencez par un test simple. Créé un Boîte de délimitation alignée sur l'axe autour de votre polygone. C'est très facile, rapide et peut déjà vous faire beaucoup de calculs. Comment ça marche? Itérer sur tous les points du polygone et trouver les valeurs min / max de X et Y.

Par exemple. vous avez les points (9/1), (4/3), (2/7), (8/2), (3/6). Cela signifie que Xmin est 2, Xmax est 9, Ymin est 1 et Ymax est 7. Un point en dehors du rectangle avec les deux arêtes (2/1) et (9/7) ne peut pas être 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!
}

Ceci est le premier test à exécuter pour n'importe quel point. Comme vous pouvez le voir, ce test est ultra rapide mais il est également très grossier. Pour gérer les points situés dans le rectangle englobant, nous avons besoin d'un algorithme plus sophistiqué. Il y a plusieurs façons de calculer cela. La méthode qui fonctionne dépend également du fait que le polygone peut avoir des trous ou sera toujours solide. Voici des exemples de solides (un convexe, un concave):

Polygon without hole

Et voici un avec un trou:

Polygon with hole

Le vert a un trou au milieu!

L'algorithme le plus simple, qui peut gérer les trois cas ci-dessus et est encore assez rapide, s'appelle Ray casting. L'idée de l'algorithme est assez simple: dessinez 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 hits est pair, c'est en dehors du polygone, si c'est impair, c'est à l'intérieur.

Demonstrating how the ray cuts through a polygon

le algorithme de numéro de bobinage serait une alternative, elle est plus précise pour les points très proches d'une ligne polygonale mais aussi beaucoup plus lente. Le lancer de rayons peut échouer pour des points trop proches d'un polygone en raison de la précision limitée des virgules flottantes et des problèmes d'arrondi, mais en réalité, cela ne pose guère de problème, car si un point est proche d'un côté, spectateur pour savoir s'il est déjà à l'intérieur ou à l'extérieur.

Vous avez toujours la boîte de délimitation ci-dessus, vous vous souvenez? 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 en dehors du polygone à coup sûr.

Mais comment ça e? Bien, e (en fait epsilon) donne la boîte englobante rembourrage. Comme je l'ai dit, le lancer de rayons échoue si nous commençons trop près d'une ligne polygonale. Depuis la zone de délimitation pourrait égaler le polygone (si le polygone est un rectangle aligné d'axe, la zone de délimitation est égal au polygone lui-même!), Nous avons besoin de rembourrage pour que cela soit sécuritaire, voilà tout. Quelle taille devriez-vous choisir e? Pas si gros. Cela dépend de l'échelle du système de coordonnées que vous utilisez pour dessiner. Si la largeur de votre pas de pixels est de 1,0, choisissez simplement 1.0 (pourtant, 0,1 aurait aussi bien fonctionné)

Maintenant que nous avons le rayon avec ses coordonnées de début et de fin, le problème change de "est le point dans le polygone" à "à quelle fréquence le rayon croise-t-il un côté polygonalPar conséquent, nous ne pouvons pas simplement travailler avec les points du polygone comme auparavant, nous avons maintenant besoin des côtés réels, un cô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 que le rayon est un vecteur et que chaque côté est un vecteur. Le rayon doit frapper de chaque côté exactement une fois ou jamais du tout. Il ne peut pas frapper le même côté deux fois. Deux lignes de l'espace 2D se croisent toujours une seule fois, à moins qu'elles ne soient parallèles, auquel cas elles ne se croisent jamais. Cependant, comme les vecteurs ont une longueur limitée, deux vecteurs peuvent ne pas être parallèles et ne se croisent jamais parce qu'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 testez-vous si deux vecteurs se croisent? Voici du 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 points d'extrémité du vecteur 1 (v1x1/v1y1 et v1x2/v1y2) et le vecteur 2 (v2x1/v2y1 et v2x2/v2y2). Donc, vous avez 2 vecteurs, 4 points, 8 coordonnées. YES et NO sont claires. YES augmente les intersections, NO ne fait rien.

Qu'en est-il de COLLINEAR? Cela signifie que les deux vecteurs se trouvent sur la même ligne infinie, en fonction de la position et de la longueur, qu'ils ne se croisent pas du tout ou qu'ils 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. Eh bien, ce cas est plutôt rare dans la pratique de toute façon en raison des erreurs d'arrondi à virgule flottante; un meilleur code ne teste probablement pas pour == 0.0f mais au lieu de quelque chose comme < epsilon, 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 la chose un peu en gardant les formes linéaires standard d'équations des côtés du polygone dans la mémoire, de sorte que vous n'avez pas à recalculer chaque fois que ceux-ci. Cela vous permettra d'économiser deux multiplications à virgule flottante et trois soustractions à virgule flottante sur chaque test en échange du stockage de trois valeurs à virgule flottante par côté polygone en mémoire. C'est un compromis typique entre la mémoire et le temps de calcul.

Last but not least: Si vous utilisez 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. Maintenant, nous allons OpenGL ou Direct3D peindre votre polygone (ou même tous vos polygones si vous voulez juste tester si le point est dans l'un d'eux, mais vous ne se soucient pas lequel) et remplir le polygone (s) avec un autre 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 récupération de mémoire O (1).

Bien sûr, cette méthode n'est utilisable que si votre surface de dessin ne doit pas être énorme. S'il ne peut pas entrer dans la mémoire GPU, cette méthode est plus lente que sur le processeur. Si cela devait être énorme et que votre GPU supporte les shaders modernes, vous pouvez toujours utiliser le GPU en implémentant le lancer de rayon montré ci-dessus en tant que shader GPU, ce qui est absolument possible. Pour un plus grand nombre de polygones ou un grand nombre de points à tester, cela sera payant, sachant que certains GPU pourront tester 64 à 256 points en parallèle. Notez cependant que le transfert de données du processeur vers le GPU et vers l'arrière est toujours coûteux, donc pour tester quelques points contre deux simples polygones, où les points ou les polygones sont dynamiques et changeront fréquemment, une approche GPU paiera rarement de.


634
2017-10-20 11:19



Je pense que le morceau de code suivant est la meilleure solution (tirée 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 dans le polygone. Que ce soit pour répéter le premier sommet à la fin a été discuté dans l'article mentionné ci-dessus.
  • vertx, verty: Tableaux contenant les coordonnées x et y des sommets du polygone.
  • testx, testy: Coordonnée X et Y du point de test.

Il est à la fois court et efficace et fonctionne à la fois pour les 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 polygonaux séparément.

L'idée derrière cela est assez simple. L'auteur le décrit comme suit:

Je cours un rayon semi-infini horizontalement (en augmentant x, y fixe) à partir du point de test, et compte combien d'arêtes il traverse. A chaque passage, le rayon change entre l'intérieur et l'extérieur. C'est ce qu'on appelle le théorème de la courbe de Jordan.

La variable c passe de 0 à 1 et 1 à 0 chaque fois que le rayon horizontal traverse une arête. Donc, fondamentalement, il reste à savoir si le nombre d'arêtes croisées est pair ou impair. 0 signifie pair et 1 signifie impair.


519
2018-05-27 16:08



Voici une version C # du réponse donnée par nirg, qui vient de ce professeur RPI. Notez que l'utilisation du code provenant de cette source RPI nécessite une attribution.

Un contrôle de boîte englobante a été ajouté en haut. Cependant, comme le souligne James Brown, le code principal est presque aussi rapide que la vérification du cadre de sélection, de sorte que la vérification du cadre de sélection peut ralentir l’opération dans le cas où la plupart des points vérifiés se trouvent dans le cadre . Vous pouvez donc quitter la boîte englobante ou utiliser une alternative pour précalculer les limites de vos polygones s'ils ne changent pas trop souvent de forme.

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

55
2018-05-06 04:14



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

39
2017-07-05 14:11



Calculer la somme orientée des angles entre le point p et chacun des sommets des polygones. Si l'angle total 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 parce qu'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 "atteindre" un sommet lors du calcul du nombre d'intersections.

EDIT: By The Way, cette méthode fonctionne avec les polygones concaves et convexes.

EDIT: j'ai récemment trouvé un ensemble Article Wikipedia sur le sujet.


25
2017-10-20 05:47



le 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 des angles est vraiment mauvaise comparée aux autres. Il est également intéressant de noter que des optimisations, telles que l’utilisation d’une grille de recherche pour subdiviser le polygone en secteurs «in» et «out», peuvent rendre le test incroyablement rapide, même sur des polygones de plus de 1000 côtés.

Quoi qu'il en soit, il est tôt mais mon vote va à la méthode des "passages à niveau", ce que je pense à peu près ce que décrit Mecki. Cependant, je l'ai trouvé le plus rapidement décrit et codifié par David Bourke. J'adore le fait qu'il n'y ait pas de véritable trigonométrie, et que cela fonctionne pour convexe et concave, et que le nombre de côtés augmente.

En passant, voici l'une des tables de performance de l'article d'Eric Haines pour l'intérêt, les tests sur 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

16
2017-12-29 04:28



Cette question est si intéressante. J'ai une autre idée pratique 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 l'angle formé par la cible et tous les deux points frontaliers sera de 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. C'est comme calculer le numéro d'enroulement.

Référez-vous à cette image pour avoir une compréhension de base de l'idée: enter image description here

Mon algorithme suppose que le sens des aiguilles d'une montre est la direction positive. Voici une contribution 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

15
2018-05-06 15:24



J'ai fait un travail là-dessus quand j'étais chercheur sous Michael Stonebraker - Vous savez, le professeur qui est venu avec Ingres, PostgreSQL, etc.

Nous nous sommes rendu compte que le moyen le plus rapide était de faire d'abord une boîte englobante parce que c'est SUPER rapide. Si c'est en dehors de la boîte englobante, c'est à l'extérieur. Sinon, vous faites le travail plus difficile ...

Si vous voulez un bon algorithme, regardez le code source du projet open source PostgreSQL pour le travail de géo ...

Je tiens à souligner que nous n'avons jamais eu d'idée sur le droit ou le gaucher (aussi exprimable comme un problème «intérieur» ou «extérieur» ...


METTRE À JOUR

Le lien de BKB fournissait un bon nombre d'algorithmes raisonnables. Je travaillais sur des problèmes de sciences de la Terre et j'avais donc besoin d'une solution qui fonctionne en latitude / longitude et qui a le problème particulier de la maîtrise des mains - la zone se situe-t-elle dans la zone plus petite ou la plus grande? La réponse est que la "direction" des vertiges a son importance - elle est soit gaucher, soit droitier et vous pouvez ainsi indiquer l'une ou l'autre des zones comme "à l'intérieur" d'un polygone donné. En tant que tel, mon travail a utilisé la solution trois énumérée sur cette page.

De plus, mon travail a utilisé des fonctions séparées pour les tests «en ligne».

... Depuis que quelqu'un a demandé: nous avons compris que les tests de boîte englobaient mieux quand le nombre de verticies dépassait un certain nombre - faites un test très rapide avant de faire le test le plus long si nécessaire ... Une boîte englobante est créée en prenant simplement la le plus grand x, le plus petit x, le plus grand y et le plus petit y et les assemblant pour faire quatre points d'une boîte ...

Une autre astuce pour ceux qui suivent: nous avons fait tous nos calculs plus sophistiqués et "gradation de lumière" dans un espace de grille, tous en points positifs sur un plan, puis re-projetés en longitude / latitude "réelle", évitant ainsi d'éventuelles erreurs de s'enrouler quand on a franchi la ligne 180 de longitude et lors de la manipulation des régions polaires. Travaillé super!


7
2017-10-20 05:44