Wie kann ich feststellen, ob ein 2D-Punkt innerhalb eines Polygons liegt?

Question

Ich versuche, eine schnell 2D-Punkt-in-Polygon-Algorithmus zur Verwendung bei der Trefferprüfung (z. B. Polygon.contains(p:Point) ). Für Vorschläge für wirksame Techniken wären wir dankbar.

Answer 1

Bei Grafiken würde ich eher keine ganzen Zahlen bevorzugen. Viele Systeme verwenden Integer für die UI-Malerei (Pixel sind schließlich Ints), aber macOS, zum Beispiel, verwendet Float für alles. macOS kennt nur Punkte und ein Punkt kann mit einem Pixel übersetzt werden, aber je nach Bildschirmauflösung kann er mit etwas anderem übersetzt werden. Auf Retina-Bildschirmen ist ein halber Punkt (0,5/0,5) ein Pixel. Trotzdem ist mir nie aufgefallen, dass die macOS-Benutzeroberfläche wesentlich langsamer ist als andere Benutzeroberflächen. Schließlich arbeiten 3D-APIs (OpenGL oder Direct3D) auch mit Fließkommazahlen und moderne Grafikbibliotheken nutzen sehr oft die Vorteile der GPU-Beschleunigung.

Sie sagten, die Geschwindigkeit sei Ihr Hauptanliegen, okay, dann nehmen wir die Geschwindigkeit. Bevor Sie einen ausgeklügelten Algorithmus anwenden, führen Sie zunächst einen einfachen Test durch. Erstellen Sie eine axial ausgerichteter Begrenzungsrahmen um Ihr Polygon herum. Das ist sehr einfach, schnell und kann Ihnen bereits eine Menge an Berechnungen ersparen. Wie funktioniert das? Iterieren Sie über alle Punkte des Polygons und finden Sie die Min/Max-Werte von X und Y.

Sie haben z.B. die Punkte (9/1), (4/3), (2/7), (8/2), (3/6) . Das bedeutet: Xmin ist 2, Xmax ist 9, Ymin ist 1 und Ymax ist 7. Ein Punkt außerhalb des Rechtecks mit den beiden Kanten (2/1) und (9/7) kann nicht innerhalb des Polygons liegen.

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

Dies ist der erste Test, der für jeden Punkt durchgeführt wird. Wie Sie sehen können, ist dieser Test extrem schnell, aber auch sehr grob. Um Punkte zu behandeln, die innerhalb des begrenzenden Rechtecks liegen, brauchen wir einen ausgefeilteren Algorithmus. Es gibt mehrere Möglichkeiten, wie dies berechnet werden kann. Welche Methode am besten funktioniert, hängt auch davon ab, ob das Polygon Löcher haben kann oder ob es immer massiv ist. Hier sind Beispiele für solide Polygone (eines konvex, eines konkav):

Und hier ist einer mit einem Loch:

Der grüne hat in der Mitte ein Loch!

Der einfachste Algorithmus, der alle drei oben genannten Fälle bewältigen kann und trotzdem recht schnell ist, heißt Strahlenguss . Die Idee des Algorithmus ist ziemlich einfach: Zeichne einen virtuellen Strahl von irgendwo außerhalb des Polygons zu deinem Punkt und zähle, wie oft er eine Seite des Polygons trifft. Wenn die Anzahl der Treffer gerade ist, liegt er außerhalb des Polygons, wenn sie ungerade ist, liegt er innerhalb.

El Wicklungszahl-Algorithmus wäre eine Alternative, sie ist genauer für Punkte, die sehr nahe an einer Polygonlinie liegen, aber sie ist auch viel langsamer. Ray Casting kann bei Punkten, die zu nahe an einer Polygonkante liegen, wegen der begrenzten Fließkommagenauigkeit und Rundungsproblemen scheitern, aber in der Realität ist das kaum ein Problem, denn wenn ein Punkt so nahe an einer Kante liegt, ist es für einen Betrachter oft visuell nicht einmal möglich zu erkennen, ob er bereits innerhalb oder noch außerhalb liegt.

Sie haben immer noch die Bounding Box von oben, erinnern Sie sich? Wählen Sie einfach einen Punkt außerhalb des Begrenzungsrahmens und verwenden Sie ihn als Startpunkt für Ihren Strahl. Z.B. der Punkt (Xmin - e/p.y) liegt mit Sicherheit außerhalb des Polygons.

Aber was ist e ? Nun ja, e (eigentlich Epsilon) gibt dem Begrenzungsrahmen etwas Polsterung . Wie ich bereits sagte, schlägt das Raytracing fehl, wenn wir zu nahe an einer Polygonlinie beginnen. Da das Begrenzungsrechteck dem Polygon entsprechen könnte (wenn das Polygon ein achsenausgerichtetes Rechteck ist, ist das Begrenzungsrechteck gleich dem Polygon selbst!), brauchen wir etwas Polsterung, um dies sicher zu machen, das ist alles. Wie groß sollten Sie wählen e ? Nicht zu groß. Das hängt von der Skala des Koordinatensystems ab, das Sie zum Zeichnen verwenden. Wenn Ihre Pixelschrittweite 1,0 ist, dann wählen Sie einfach 1,0 (0,1 hätte aber auch funktioniert)

Da wir nun den Strahl mit seinen Start- und Endkoordinaten haben, verschiebt sich das Problem von " ist der Punkt innerhalb des Polygons " zu " wie oft schneidet der Strahl eine Polygonseite ". Deshalb können wir nicht nur mit den Punkten des Polygons arbeiten, sondern brauchen jetzt die tatsächlichen Seiten. Eine Seite ist immer durch zwei Punkte definiert.

side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:

Sie müssen den Strahl von allen Seiten testen. Betrachte den Strahl als einen Vektor und jede Seite als einen Vektor. Der Strahl muss jede Seite genau einmal oder überhaupt nicht treffen. Er kann nicht zweimal auf dieselbe Seite treffen. Zwei Linien im 2D-Raum schneiden sich immer genau einmal, es sei denn, sie sind parallel, in diesem Fall schneiden sie sich nie. Da Vektoren jedoch eine begrenzte Länge haben, kann es sein, dass zwei Vektoren nicht parallel sind und sich trotzdem nicht schneiden, weil sie zu kurz sind, um sich jemals zu treffen.

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

So weit, so gut, aber wie prüft man, ob sich zwei Vektoren kreuzen? Hier ist etwas C-Code (nicht getestet), der den Zweck erfüllen sollte:

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

Die Eingabewerte sind die zwei Endpunkte des Vektors 1 ( v1x1/v1y1 y v1x2/v1y2 ) und Vektor 2 ( v2x1/v2y1 y v2x2/v2y2 ). Sie haben also 2 Vektoren, 4 Punkte und 8 Koordinaten. YES y NO sind klar. YES erhöht die Zahl der Kreuzungen, NO tut nichts.

Was ist mit COLLINEAR? Das bedeutet, dass beide Vektoren auf der gleichen unendlichen Linie liegen, je nach Position und Länge schneiden sie sich entweder gar nicht oder in unendlich vielen Punkten. Ich bin mir nicht ganz sicher, wie ich diesen Fall behandeln soll, ich würde ihn so oder so nicht als Schnittpunkt betrachten. Nun, dieser Fall ist in der Praxis wegen der Fließkomma-Rundungsfehler ohnehin eher selten; besserer Code würde wahrscheinlich nicht auf == 0.0f sondern stattdessen für etwas wie < epsilon , wobei epsilon eine eher kleine Zahl ist.

Wenn Sie eine größere Anzahl von Punkten testen müssen, können Sie das Ganze sicherlich etwas beschleunigen, indem Sie die Standardformen der linearen Gleichungen für die Polygonseiten im Speicher halten, so dass Sie diese nicht jedes Mal neu berechnen müssen. Das erspart Ihnen bei jedem Test zwei Fließkommamultiplikationen und drei Fließkommasubtraktionen im Austausch für die Speicherung von drei Fließkommawerten pro Polygonseite im Speicher. Das ist ein typischer Kompromiss zwischen Speicher und Rechenzeit.

Und nicht zuletzt: Wenn Sie 3D-Hardware verwenden können, um das Problem zu lösen, gibt es eine interessante Alternative. Lassen Sie einfach die GPU die ganze Arbeit für Sie erledigen. Erstellen Sie eine Malfläche, die sich außerhalb des Bildschirms befindet. Füllen Sie sie vollständig mit der Farbe Schwarz. Lassen Sie nun OpenGL oder Direct3D Ihr Polygon malen (oder sogar alle Ihre Polygone, wenn Sie nur testen wollen, ob der Punkt innerhalb eines von ihnen liegt, aber es ist Ihnen egal, welches) und füllen Sie das/die Polygon(e) mit einer anderen Farbe, z. B. Weiß. Um zu prüfen, ob ein Punkt innerhalb des Polygons liegt, ermitteln Sie die Farbe dieses Punktes aus der Zeichenfläche. Dies ist nur ein O(1)-Speicherabruf.

Natürlich ist diese Methode nur dann brauchbar, wenn die Zeichenfläche nicht riesig sein muss. Wenn sie nicht in den GPU-Speicher passt, ist diese Methode langsamer als die Ausführung auf der CPU. Wenn sie sehr groß sein muss und Ihr Grafikprozessor moderne Shader unterstützt, können Sie den Grafikprozessor trotzdem nutzen, indem Sie das oben gezeigte Raycasting als GPU-Shader implementieren, was durchaus möglich ist. Bei einer größeren Anzahl von Polygonen oder einer großen Anzahl von Punkten, die getestet werden müssen, wird sich dies auszahlen, da einige GPUs in der Lage sind, 64 bis 256 Punkte parallel zu testen. Beachten Sie jedoch, dass die Übertragung von Daten von der CPU zur GPU und zurück immer kostspielig ist, so dass sich ein GPU-Ansatz für das Testen von nur ein paar Punkten gegen ein paar einfache Polygone, bei denen entweder die Punkte oder die Polygone dynamisch sind und sich häufig ändern werden, nur selten lohnt.

Answer 2

Hier ist eine C#-Version der Antwort von nirg die von dieser RPI-Professor . Beachten Sie, dass die Verwendung des Codes aus dieser RPI-Quelle eine Namensnennung erfordert.

Oben wurde ein Kontrollkästchen für Begrenzungsrahmen hinzugefügt. Wie James Brown jedoch feststellt, ist der Hauptcode fast so schnell wie die Bounding-Box-Prüfung selbst, so dass die Bounding-Box-Prüfung den Gesamtvorgang verlangsamen kann, wenn die meisten der zu prüfenden Punkte innerhalb der Bounding-Box liegen. Sie könnten also die Bounding-Box-Prüfung weglassen oder alternativ die Bounding-Boxen Ihrer Polygone vorberechnen, wenn sie ihre Form nicht allzu oft ändern.

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

    // https://wrf.ecse.rpi.edu/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;
}

Answer 3

Hier ist eine JavaScript-Variante der Antwort von M. Katz, die auf dem Ansatz von Nirg basiert:

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

Answer 4

Berechnen Sie die orientierte Winkelsumme zwischen dem Punkt p und jedem der Polygonscheitelpunkte. Wenn die orientierte Winkelsumme 360 Grad beträgt, ist der Punkt innen. Wenn die Summe 0 ist, ist der Punkt außen.

Mir gefällt diese Methode besser, weil sie robuster ist und weniger von der numerischen Präzision abhängt.

Methoden, die die Gleichmäßigkeit der Anzahl der Schnittpunkte berechnen, sind begrenzt, weil man während der Berechnung der Anzahl der Schnittpunkte einen Scheitelpunkt "treffen" kann.

EDIT: Diese Methode funktioniert übrigens auch mit konkaven und konvexen Polygonen.

EDIT: Ich habe kürzlich eine ganze Wikipedia-Artikel zu diesem Thema.

Answer 5

Diese Frage ist sehr interessant. Ich habe eine weitere praktikable Idee, die sich von anderen Antworten auf diesen Beitrag unterscheidet. Die Idee ist, die Summe der Winkel zu verwenden, um zu entscheiden, ob das Ziel innerhalb oder außerhalb liegt. Besser bekannt als Windungszahl .

x sei der Zielpunkt. Das Feld [0, 1, .... n] sei die Gesamtheit der Punkte des Gebietes. Verbinden Sie den Zielpunkt mit jedem Grenzpunkt durch eine Linie. Wenn der Zielpunkt innerhalb der Fläche liegt. Die Summe aller Winkel wird 360 Grad betragen. Wenn nicht, sind die Winkel kleiner als 360.

Die folgende Abbildung vermittelt ein grundlegendes Verständnis der Idee:

Mein Algorithmus geht davon aus, dass der Uhrzeigersinn die positive Richtung ist. Hier ist eine mögliche Eingabe:

[[-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]]

Im Folgenden finden Sie den Python-Code, der diese Idee umsetzt:

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