Wie berechne ich die Fläche eines 2D-Polygons?

81

Was ist eine effiziente Methode zur Bestimmung der Fläche des resultierenden Polygons, wenn eine Reihe von Punkten im 2D-Raum angenommen wird, die sich nicht selbst schneiden?

Nebenbei bemerkt, dies sind keine Hausaufgaben und ich suche keinen Code. Ich suche nach einer Beschreibung, mit der ich meine eigene Methode implementieren kann. Ich habe meine Ideen, eine Folge von Dreiecken aus der Liste der Punkte zu ziehen, aber ich weiß, dass es eine Reihe von Randfällen in Bezug auf konvexe und konkave Polygone gibt, die ich wahrscheinlich nicht fangen werde.

Jason Z.
quelle
6
Der Begriff "Oberfläche" ist etwas irreführend. Was Sie zu wollen scheinen, ist nur der (reguläre) Bereich. In 3D ist die Oberfläche die Fläche der Außenfläche, daher wäre die natürliche 2D-Verallgemeinerung dieses Konzepts die Länge des Umfangs des Polygons, was eindeutig nicht das ist, wonach Sie suchen.
Batty
def area (polygon): return abs (numpy.cross (polygon, numpy.roll (polygon, -1, 0)). sum () / 2)
iouvxz

Antworten:

110

Hier ist die Standardmethode AFAIK. Summieren Sie grundsätzlich die Kreuzprodukte um jeden Scheitelpunkt. Viel einfacher als Triangulation.

Python-Code mit einem Polygon, das als Liste von (x, y) Scheitelpunktkoordinaten dargestellt wird und implizit vom letzten zum ersten Scheitelpunkt umläuft:

def area(p):
    return 0.5 * abs(sum(x0*y1 - x1*y0
                         for ((x0, y0), (x1, y1)) in segments(p)))

def segments(p):
    return zip(p, p[1:] + [p[0]])

David Lehavi kommentiert: Es ist erwähnenswert, warum dieser Algorithmus funktioniert: Es ist eine Anwendung des Greenschen Theorems für die Funktionen −y und x; genau so, wie ein Planimeter funktioniert. Genauer:

Formel oben =
integral_over_perimeter(-y dx + x dy) =
integral_over_area((-(-dy)/dy+dx/dx) dy dx) =
2 Area

Darius Bacon
quelle
7
Es ist erwähnenswert, warum dieser Algorithmus funktioniert: Es ist eine Anwendung des Greenschen Theorems für die Funktionen -y und x; genau so, wie ein Planimeter funktioniert. Genauer gesagt: Formel oben = Integral_Permieter (-y dx + x dy) = Integralbereich ((- (- dy) / dy + dx / dx) Dydyx = 2 Bereich
David Lehavi
6
Der Link im Beitrag ist tot. Hat jemand einen anderen?
Yakov
1
Die verknüpfte Diskussion auf der Mailingliste [email protected] ist für mich nicht verfügbar. Ich habe die Nachricht aus dem Google Cache kopiert: gist.github.com/1200393
Andrew 7ндрей Листочкин
2
@ perfectm1ng Richtungswechsel würden das Vorzeichen in der Summe umdrehen, aber abs()das Vorzeichen entfernen.
Darius Bacon
3
Einschränkungen: Diese Methode liefert die falsche Antwort für sich selbst schneidende Polygone, bei denen eine Seite die andere kreuzt, wie rechts gezeigt. Es funktioniert jedoch korrekt für Dreiecke, regelmäßige und unregelmäßige Polygone, konvexe oder konkave Polygone. ( mathopenref.com/coordpolygonarea.html )
OneWorld
14

Das Kreuzprodukt ist ein Klassiker.

Wenn Sie zig Millionen solcher Berechnungen durchführen müssen, versuchen Sie die folgende optimierte Version, die halb weniger Multiplikationen erfordert:

area = 0;
for( i = 0; i < N; i += 2 )
   area += x[i+1]*(y[i+2]-y[i]) + y[i+1]*(x[i]-x[i+2]);
area /= 2;

Ich benutze Array-Index zur Klarheit. Es ist effizienter, Zeiger zu verwenden. Gute Compiler erledigen das für Sie.

Es wird angenommen, dass das Polygon "geschlossen" ist, was bedeutet, dass Sie den ersten Punkt als Punkt mit dem Index N kopieren. Es wird auch angenommen, dass das Polygon eine gerade Anzahl von Punkten hat. Fügen Sie eine zusätzliche Kopie des ersten Punkts hinzu, wenn N nicht gerade ist.

Der Algorithmus wird erhalten, indem zwei aufeinanderfolgende Iterationen des klassischen produktübergreifenden Algorithmus abgewickelt und kombiniert werden.

Ich bin mir nicht so sicher, wie sich die beiden Algorithmen hinsichtlich der numerischen Genauigkeit vergleichen lassen. Mein Eindruck ist, dass der obige Algorithmus besser ist als der klassische, da die Multiplikation dazu neigt, den Genauigkeitsverlust der Subtraktion wiederherzustellen. Wenn die Verwendung von Floats wie bei der GPU eingeschränkt ist, kann dies einen erheblichen Unterschied bewirken.

EDIT: "Bereich der Dreiecke und Polygone 2D & 3D" beschreibt eine noch effizientere Methode

// "close" polygon
x[N] = x[0];
x[N+1] = x[1];
y[N] = y[0];
y[N+1] = y[1];

// compute area
area = 0;
for( size_t i = 1; i <= N; ++i )
  area += x[i]*( y[i+1] - y[i-1] );
area /= 2;
chmike
quelle
1
Ich kann mir nicht vorstellen, dass das zweite Code-Snippet funktioniert. Es ist ziemlich offensichtlich, dass seine Fläche umso größer ist, je weiter das Polygon auf der X-Achse liegt.
Cygon
1
Es ist eine korrekte mathematische Neuordnung des oben beschriebenen Algorithmus, wobei einige Multiplikationen eingespart werden. Sie haben Recht, aber die durch andere Scheitelpunkte definierten Bereiche werden subtrahiert. Dies kann jedoch tatsächlich zu einer Verschlechterung der Präzision führen.
Chmike
2
Was Sie übersehen haben, ist, dass die Addition aufgrund der y-Subtraktion immer einige negative Terme hat. Betrachten Sie eine beliebige polygonale 2D-Form und vergleichen Sie die y-Werte aufeinanderfolgender Scheitelpunkte. Sie werden sehen, dass eine Subtraktion einen negativen und eine positive Wert ergibt.
Chmike
2
In der Tat ist dieser letzte Absatz das, worum ich mich nicht kümmern konnte! Mit i <= N funktioniert es. Vielen Dank für Ihre Geduld, ich nehme alles zurück :)
Cygon
1
Nebenbei bemerkt, der vom Algorithmus zurückgegebene Bereich ist "vorzeichenbehaftet" (negativ oder positiv, basierend auf der Reihenfolge der Punkte). Wenn Sie also immer einen positiven Bereich wünschen, verwenden Sie einfach den absoluten Wert.
NightElfik
11

Diese Seite zeigt, dass die Formel

Geben Sie hier die Bildbeschreibung ein

kann vereinfacht werden zu:

Geben Sie hier die Bildbeschreibung ein

Wenn Sie einige Begriffe aufschreiben und sie nach gemeinsamen Faktoren von gruppieren xi, ist die Gleichheit nicht schwer zu erkennen.

Die endgültige Summierung ist effizienter, da nstattdessen nur Multiplikationen erforderlich sind 2n.

def area(x, y):
    return abs(sum(x[i] * (y[i + 1] - y[i - 1]) for i in xrange(-1, len(x) - 1))) / 2.0

Ich habe gelernt , diese Vereinfachung von Joe Kington, hier .


Wenn Sie NumPy haben, ist diese Version schneller (für alle außer sehr kleinen Arrays):

def area_np(x, y):        
    x = np.asanyarray(x)
    y = np.asanyarray(y)
    n = len(x)
    shift_up = np.arange(-n+1, 1)
    shift_down = np.arange(-1, n-1)    
    return (x * (y.take(shift_up) - y.take(shift_down))).sum() / 2.0
unutbu
quelle
1
Danke für die NumPy-Version.
Physikmichael
4

Um die Triangulations- und Summen-Dreiecksbereiche zu erweitern, funktionieren diese, wenn Sie zufällig ein konvexes Polygon haben ODER wenn Sie einen Punkt auswählen, der keine Linien zu jedem anderen Punkt erzeugt, der das Polygon schneidet.

Für ein allgemeines nicht schneidendes Polygon müssen Sie das Kreuzprodukt der Vektoren (Referenzpunkt, Punkt a) (Referenzpunkt, Punkt b) summieren, wobei a und b "nebeneinander" liegen.

Angenommen, Sie haben eine Liste von Punkten, die das Polygon in der Reihenfolge definieren (die Reihenfolge ist, dass die Punkte i und i + 1 eine Linie des Polygons bilden):

Summe (Kreuzprodukt ((Punkt 0, Punkt i), (Punkt 0, Punkt i + 1)) für i = 1 bis n - 1.

Nehmen Sie die Größe dieses Kreuzprodukts und Sie haben die Oberfläche.

Dies behandelt konkave Polygone, ohne sich um die Auswahl eines guten Referenzpunkts kümmern zu müssen. Alle drei Punkte, die ein Dreieck erzeugen, das sich nicht innerhalb des Polygons befindet, haben ein Kreuzprodukt, das in die entgegengesetzte Richtung eines Dreiecks zeigt, das sich innerhalb des Polygons befindet, sodass die Flächen korrekt summiert werden.

MSN
quelle
3

Berechnen der Fläche des Polygons

http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=geometry1#polygon_area

int cross(vct a,vct b,vct c)
{
    vct ab,bc;
    ab=b-a;
    bc=c-b;
    return ab.x*bc.y-ab.y*bc.x;
}    
double area(vct p[],int n)
{ 
    int ar=0;
    for(i=1;i+1<n;i++)
    {
        vct a=p[i]-p[0];
        vct b=p[i+1]-p[0];
        area+=cross(a,b);
    }
    return abs(area/2.0);
}    
Steve
quelle
Dies ist eine 3 Jahre alte Frage mit 34 positiven Stimmen für die akzeptierte Antwort. Sagen Sie uns, wie Ihre Antwort besser ist als alle anderen bereits veröffentlichten Antworten.
Mark Taylor
3
Es ist ein Beispiel in c und nicht in Python. Nicht besser, aber schön, es in verschiedenen Sprachen zu haben
underdoeg
2

Oder machen Sie ein Konturintegral. Mit dem Stokes-Theorem können Sie ein Flächenintegral als Konturintegral ausdrücken. Eine kleine Gauß-Quadratur und Bob ist dein Onkel.

Duffymo
quelle
2

Sprachunabhängige Lösung:

GEGEBEN: Ein Polygon kann IMMER aus n-2 Dreiecken bestehen, die sich nicht überlappen (n = Anzahl der Punkte ODER Seiten). 1 Dreieck = 3-seitiges Polygon = 1 Dreieck; 1 Quadrat = 4-seitiges Polygon = 2 Dreiecke; etc ad nauseam QED

Daher kann ein Polygon durch "Abhacken" von Dreiecken reduziert werden, und die Gesamtfläche ist die Summe der Flächen dieser Dreiecke. Versuchen Sie es mit einem Stück Papier und einer Schere. Am besten visualisieren Sie den Vorgang, bevor Sie fortfahren.

Wenn Sie 3 aufeinanderfolgende Punkte in einem Polygonpfad nehmen und mit diesen Punkten ein Dreieck erstellen, haben Sie nur eines von drei möglichen Szenarien:

  1. Das resultierende Dreieck befindet sich vollständig im ursprünglichen Polygon
  2. Das resultierende Dreieck liegt vollständig außerhalb des ursprünglichen Polygons
  3. Das resultierende Dreieck ist teilweise im ursprünglichen Polygon enthalten

Wir sind nur an Fällen interessiert, die in die erste Option fallen (vollständig enthalten).

Jedes Mal, wenn wir eines davon finden, hacken wir es ab, berechnen seine Fläche (einfach, hier wird die Formel nicht erklärt) und erstellen ein neues Polygon mit einer Seite weniger (entspricht dem Polygon mit diesem abgeschnittenen Dreieck). bis wir nur noch ein Dreieck haben.

wie man dies programmatisch umsetzt:

Erstellen Sie ein Array von (aufeinanderfolgenden) Punkten, die den Pfad um das Polygon darstellen. Beginnen Sie bei Punkt 0. Führen Sie das Array aus den Punkten x, x + 1 und x + 2 nacheinander aus. Transformieren Sie jedes Dreieck von einer Form in eine Fläche und schneiden Sie es mit einer aus Polygon erstellten Fläche. WENN der resultierende Schnittpunkt mit dem ursprünglichen Dreieck identisch ist, ist das Dreieck vollständig im Polygon enthalten und kann abgeschnitten werden. Entfernen Sie x + 1 aus dem Array und beginnen Sie erneut mit x = 0. Andernfalls (wenn sich das Dreieck außerhalb des [teilweise oder vollständig] Polygons befindet) zum nächsten Punkt x + 1 im Array wechseln.

Wenn Sie in die Kartierung integrieren möchten und von Geopunkten ausgehen, müssen Sie außerdem ZUERST von Geopunkten in Bildschirmpunkte konvertieren. Dies erfordert die Entscheidung über eine Modellierung und Formel für die Erdform (obwohl wir die Erde eher als Kugel betrachten, handelt es sich tatsächlich um ein unregelmäßiges Ovoid (Eierform) mit Dellen). Es gibt viele Modelle, für weitere Informationen Wiki. Eine wichtige Frage ist, ob Sie den Bereich als eben oder gekrümmt betrachten. Im Allgemeinen erzeugen "kleine" Bereiche, in denen die Punkte bis zu einigen km voneinander entfernt sind, keinen signifikanten Fehler, wenn sie planar und nicht konvex betrachtet werden.

user836725
quelle
1
  1. Legen Sie einen Basispunkt fest (den konvexesten Punkt). Dies ist Ihr Drehpunkt der Dreiecke.
  2. Berechnen Sie den am weitesten links liegenden Punkt (beliebig) außer Ihrem Basispunkt.
  3. Berechnen Sie den Punkt ganz links, um Ihr Dreieck zu vervollständigen.
  4. Speichern Sie diesen triangulierten Bereich.
  5. Verschieben Sie jede Iteration um einen Punkt nach rechts.
  6. Summiere die triangulierten Bereiche
Joe Phillips
quelle
Stellen Sie sicher, dass Sie den Dreiecksbereich negieren, wenn sich der nächste Punkt "rückwärts" bewegt.
rekursiv
1

Besser als das Summieren von Dreiecken ist das Summieren von Trapezoiden im kartesischen Raum:

area = 0;
for (i = 0; i < n; i++) {
  i1 = (i + 1) % n;
  area += (vertex[i].y + vertex[i1].y) * (vertex[i1].x - vertex[i].x) / 2.0;
}
David Hanak
quelle
1

Die Implementierung der Schnürsenkelformel könnte in Numpy erfolgen. Angenommen, diese Eckpunkte:

import numpy as np
x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)

Wir können die folgende Funktion definieren, um einen Bereich zu finden:

def PolyArea(x,y):
    return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))

Und Ergebnisse erzielen:

print PolyArea(x,y)
# 0.26353377782163534

Durch das Vermeiden von Schleifen ist diese Funktion ~ 50-mal schneller als PolygonArea:

%timeit PolyArea(x,y)
# 10000 loops, best of 3: 42 µs per loop
%timeit PolygonArea(zip(x,y))
# 100 loops, best of 3: 2.09 ms per loop

Hinweis: Ich habe diese Antwort auf eine andere Frage geschrieben . Ich erwähne dies hier nur, um eine vollständige Liste der Lösungen zu erhalten.

Mahdi
quelle
0

Meine Neigung wäre, einfach Dreiecke abzuschneiden. Ich sehe nicht ein, wie irgendetwas anderes es vermeiden könnte, schrecklich haarig zu sein.

Nehmen Sie drei aufeinanderfolgende Punkte, aus denen das Polygon besteht. Stellen Sie sicher, dass der Winkel kleiner als 180 ist. Sie haben jetzt ein neues Dreieck, dessen Berechnung kein Problem darstellen sollte. Löschen Sie den Mittelpunkt aus der Punkteliste des Polygons. Wiederholen, bis nur noch drei Punkte übrig sind.

Loren Pechtel
quelle
Der haarige Teil dabei ist, dass Sie ein Problem haben, wenn Ihre drei aufeinander folgenden Punkte ein Dreieck außerhalb oder teilweise außerhalb des Polygons definieren.
Richard
@ Richard: Deshalb ist die Qualifikation um 180 Grad. Wenn Sie ein Dreieck außerhalb des Polygons abschneiden, erhalten Sie zu viele Grad.
Loren Pechtel
Möglicherweise müssen Sie besser beschreiben, wie Sie den Winkel finden. In der Ebenengeometrie gibt es keine Möglichkeit, 3 Punkte als Teil eines Dreiecks zu haben und einen Winkel oder eine Kombination von Winkeln über 180 Grad zu haben - die Prüfung scheint bedeutungslos zu sein.
Richard
@ Richard: Auf Ihrem Polygon haben Sie den Winkel jeder Kreuzung. Wenn das relevante Dreieck außerhalb des Polygons liegen würde, wäre der Winkel zwischen den beiden Segmenten größer als 180 Grad.
Loren Pechtel
Sie meinen, der Innenwinkel der beiden benachbarten Kantensegmente wäre größer als 180 Grad.
Richard
0

C Art und Weise, das zu tun:

float areaForPoly(const int numVerts, const Point *verts)
{
    Point v2;
    float area = 0.0f;

    for (int i = 0; i<numVerts; i++){
        v2 = verts[(i + 1) % numVerts];
        area += verts[i].x*v2.y - verts[i].y*v2.x;
    }

    return area / 2.0f;
}
Joseph
quelle
0

Python-Code

Wie hier beschrieben: http://www.wikihow.com/Calculate-the-Area-of-a-Polygon

Mit Pandas

import pandas as pd

df = pd.DataFrame({'x': [10, 20, 20, 30, 20, 10, 0], 'y': [-10, -10, -10, 0, 10, 30, 20]})
df = df.append(df.loc[0])

first_product = (df['x'].shift(1) * df['y']).fillna(0).sum()
second_product = (df['y'].shift(1) * df['x']).fillna(0).sum()

(first_product - second_product) / 2
600
Firelynx
quelle
0

Ich werde ein paar einfache Funktionen zur Berechnung der Fläche eines 2D-Polygons geben. Dies funktioniert sowohl für konvexe als auch für konkave Polygone. Wir teilen das Polygon einfach in viele Unterdreiecke.

//don't forget to include cmath for abs function
struct Point{
  double x;
  double y;
}
// cross_product
double cp(Point a, Point b){ //returns cross product
  return a.x*b.y-a.y*b.x;
}

double area(Point * vertices, int n){  //n is number of sides
  double sum=0.0;
  for(i=0; i<n; i++){
    sum+=cp(vertices[i], vertices[(i+1)%n]); //%n is for last triangle
  }
  return abs(sum)/2.0;
}
abe312
quelle
cpnimmt zwei Argumente, aber Sie nennen es mit einem.