Zakładając serię punktów w przestrzeni 2d, które nie przecinają się same, jaka jest skuteczna metoda wyznaczania pola powierzchni powstałego wielokąta?

Na marginesie, to nie jest praca domowa i nie szukam kodu. Szukam opisu, który mógłbym wykorzystać do wdrożenia własnej metody. Mam pomysły na wyciągnięcie sekwencji trójkątów z listy punktów, ale wiem, że jest kilka przypadków krawędzi dotyczących wielokątów wypukłych i wklęsłych, których prawdopodobnie nie złapię.

Oto standardowa metoda AFAIK. Zasadniczo zsumuj iloczyny poprzeczne wokół każdego wierzchołka. Znacznie prostsze niż triangulacja.

Kod Pythona, biorąc pod uwagę wielokąt reprezentowany jako lista współrzędnych wierzchołków (x, y), niejawnie zawijający się od ostatniego wierzchołka do pierwszego:

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 komentuje: Warto wspomnieć, dlaczego ten algorytm działa: Jest to aplikacja twierdzenia Greena dla funkcji −y i x; dokładnie tak, jak działa planymetr . Dokładniej:

Wzór powyżej =
integral_over_perimeter(-y dx + x dy) =
integral_over_area((-(-dy)/dy+dx/dx) dy dx) =
2 Area

Produkt krzyżowy to klasyk.

Jeśli masz do wykonania milion takich obliczeń, wypróbuj następującą zoptymalizowaną wersję, która wymaga o połowę mniej mnożeń:

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;

Dla przejrzystości używam indeksu tablicy. Bardziej efektywne jest używanie wskaźników. Chociaż dobrzy kompilatorzy zrobią to za Ciebie.

Zakłada się, że wielokąt jest „zamknięty”, co oznacza, że ​​kopiujesz pierwszy punkt jako punkt z indeksem dolnym N. Zakłada się również, że wielokąt ma parzystą liczbę punktów. Dołącz dodatkową kopię pierwszego punktu, jeśli N nie jest parzyste.

Algorytm uzyskuje się poprzez rozwinięcie i połączenie dwóch kolejnych iteracji klasycznego algorytmu krzyżowego.

Nie jestem pewien, jak oba algorytmy porównują się pod względem dokładności numerycznej. Mam wrażenie, że powyższy algorytm jest lepszy od klasycznego, ponieważ mnożenie ma tendencję do przywracania utraty precyzji odejmowania. W przypadku ograniczenia do używania pływaków, tak jak w przypadku GPU, może to mieć znaczący wpływ.

EDYCJA: „Obszar trójkątów i wielokątów 2D i 3D” opisuje jeszcze wydajniejszą metodę

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

Ta strona pokazuje, że formuła

można uprościć do:

Jeśli napiszesz kilka terminów i pogrupujesz je według wspólnych czynników xi, nie jest trudno zauważyć równość.

Końcowe sumowanie jest bardziej wydajne, ponieważ wymaga tylko nmnożenia zamiast 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

Nauczyłem się tego uproszczenia od Joe Kingtona tutaj .

Jeśli masz NumPy, ta wersja jest szybsza (dla wszystkich oprócz bardzo małych tablic):

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

Zestaw punktów bez żadnych innych ograniczeń niekoniecznie musi jednoznacznie definiować wielokąt.

Więc najpierw musisz zdecydować, jaki wielokąt zbudować z tych punktów - może wypukły kadłub? http://en.wikipedia.org/wiki/Convex_hull

Następnie wykonaj triangulację i oblicz obszar. http://www.mathopenref.com/polygonirregulararea.html

Aby rozszerzyć obszary trójkątów i trójkątów sumujących, działają one, jeśli zdarzy się, że masz wypukły wielokąt LUB wybierzesz punkt, który nie generuje linii do każdego innego punktu przecinającego wielokąt.

W przypadku ogólnego nieprzecinającego się wielokąta należy zsumować iloczyn poprzeczny wektorów (punkt odniesienia, punkt a), (punkt odniesienia, punkt b), gdzie aib są „obok” siebie.

Zakładając, że masz listę punktów, które definiują wielokąt w kolejności (kolejność, w której punkty i i i + 1 tworzą linię wielokąta):

Suma (iloczyn poprzeczny ((punkt 0, punkt i), (punkt 0, punkt i + 1)) dla i = 1 do n - 1.

Weź wielkość tego iloczynu krzyżowego i masz pole powierzchni.

Pozwoli to obsłużyć wklęsłe wielokąty bez martwienia się o wybranie dobrego punktu odniesienia; dowolne trzy punkty, które generują trójkąt, który nie znajduje się wewnątrz wielokąta, będą miały iloczyn poprzeczny wskazujący w przeciwnym kierunku niż dowolny trójkąt wewnątrz wielokąta, więc obszary zostaną poprawnie zsumowane.

Aby obliczyć pole wielokąta

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

Lub wykonaj całkę konturową. Twierdzenie Stokesa pozwala wyrazić całkę powierzchniową jako całkę konturową. Trochę kwadratury Gaussa i Bob jest twoim wujem.

rozwiązanie niezależne od języka:

PODANE: wielokąt może ZAWSZE składać się z n-2 trójkątów, które nie zachodzą na siebie (n = liczba punktów LUB boków). 1 trójkąt = wielokąt trójstronny = 1 trójkąt; 1 kwadrat = wielokąt czteroboczny = 2 trójkąty; etc ad mdłości QED

dlatego też wielokąt można zmniejszyć przez „ucinanie” trójkątów, a całkowita powierzchnia będzie sumą powierzchni tych trójkątów. wypróbuj to za pomocą kartki papieru i nożyczek, najlepiej wizualizować proces przed wykonaniem.

jeśli weźmiesz dowolne 3 kolejne punkty na ścieżce wielokąta i utworzysz trójkąt z tymi punktami, będziesz mieć jeden i tylko jeden z trzech możliwych scenariuszy:

1. wynikowy trójkąt jest całkowicie wewnątrz oryginalnego wielokąta
2. wynikowy trójkąt jest całkowicie poza oryginalnym wielokątem
3. wynikowy trójkąt jest częściowo zawarty w oryginalnym wielokącie

interesują nas tylko przypadki, które mieszczą się w pierwszej opcji (całkowicie zawarte).

za każdym razem, gdy znajdujemy jeden z nich, odcinamy go, obliczamy jego powierzchnię (łatwy peasy, nie wyjaśniam tutaj wzoru) i tworzymy nowy wielokąt z jednym bokiem mniej (odpowiednik wielokąta z odciętym trójkątem). dopóki nie zostanie nam tylko jeden trójkąt.

jak zaimplementować to programowo:

utwórz tablicę (kolejnych) punktów, które reprezentują ścieżkę WOKÓŁ wielokąta. zacznij od punktu 0. uruchom tablicę tworząc trójkąty (po jednym naraz) z punktów x, x + 1 i x + 2. przekształć każdy trójkąt z kształtu w obszar i przetnij go z obszarem utworzonym z wielokąta. JEŻELI wynikowe przecięcie jest identyczne z oryginalnym trójkątem, wówczas wspomniany trójkąt jest całkowicie zawarty w wielokącie i można go odciąć. usuń x + 1 z tablicy i zacznij ponownie od x = 0. w przeciwnym razie (jeśli trójkąt jest na zewnątrz [częściowo lub całkowicie] wielokąta), przejdź do następnego punktu x + 1 w tablicy.

dodatkowo, jeśli chcesz zintegrować się z mapowaniem i zaczynasz od punktów geograficznych, najpierw musisz dokonać konwersji z punktów geograficznych na punkty ekranowe. wymaga to ustalenia modelu i wzoru na kształt ziemi (chociaż myślimy o ziemi jako kuli, w rzeczywistości jest to nieregularny jajo (kształt jajka) z wgnieceniami). istnieje wiele modeli, aby uzyskać więcej informacji na wiki. Ważną kwestią jest to, czy uznasz ten obszar za płaszczyznę, czy za zakrzywiony. Ogólnie rzecz biorąc, „małe” obszary, w których punkty są oddalone od siebie do kilku kilometrów, nie będą generować znaczącego błędu, jeśli uznamy je za płaskie, a nie wypukłe.

Jednym ze sposobów byłoby rozłożenie wielokąta na trójkąty , obliczenie powierzchni trójkątów i potraktowanie sumy jako powierzchni wielokąta.

1. Ustaw punkt bazowy (najbardziej wypukły punkt). To będzie twój punkt obrotu trójkątów.
2. Oblicz punkt znajdujący się najbardziej po lewej stronie (dowolny), inny niż punkt bazowy.
3. Oblicz drugi skrajny lewy punkt, aby ukończyć swój trójkąt.
4. Zapisz ten triangulowany obszar.
5. W każdej iteracji przesuń o jeden punkt w prawo.
6. Zsumuj obszary podzielone na trójkąty

Lepsze niż sumowanie trójkątów jest sumowanie trapezów w przestrzeni kartezjańskiej:

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

Implementację formuły Shoelace można było wykonać w Numpy. Zakładając te wierzchołki:

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

Aby znaleźć pole, możemy zdefiniować następującą funkcję:

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

I uzyskiwanie wyników:

print PolyArea(x,y)
# 0.26353377782163534

Unikanie pętli sprawia, że ​​ta funkcja jest ~ 50X szybsza niż 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

Uwaga: napisałem tę odpowiedź na inne pytanie , wspominam o tym tutaj, aby mieć pełną listę rozwiązań.

Wolałabym po prostu zacząć odcinać trójkąty. Nie widzę, jak cokolwiek innego mogłoby uniknąć okropnego owłosienia.

Weź trzy kolejne punkty, które składają się na wielokąt. Upewnij się, że kąt jest mniejszy niż 180. Masz teraz nowy trójkąt, którego obliczenie nie powinno stanowić problemu, usuń środkowy punkt z listy punktów wielokąta. Powtarzaj, aż pozostaną tylko trzy punkty.

C sposób na zrobienie tego:

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

Kod Pythona

Jak opisano tutaj: http://www.wikihow.com/Calculate-the-Area-of-a-Polygon

Z pandami

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

Podam kilka prostych funkcji do obliczania pola powierzchni wielokąta 2D. Działa to zarówno w przypadku wielokątów wypukłych, jak i wklęsłych. po prostu dzielimy wielokąt na wiele pod-trójkątów.

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