Obliczanie sferycznego centroidu wielokąta

Chciałbym ogólny sposób obliczania centroidów dla wielokątów na kuli.

Jak dotąd najlepszym źródłem informacji online jest:

Narzędzia do grafiki i kształtów autorstwa Jeffa Jennessa.

Opisana tam metoda sugeruje rozkład wielokąta na wiele trójkątów sferycznych i obliczenie średniej sferycznych centroidów trójkątów, ważonej przez obszar trójkąta sferycznego.

Wiem, że istnieje kilka sposobów definiowania sferycznego centroidu wielokąta, ale szukam czegoś analogicznego do następujących definicji punktów i polilinii:

• Punkty : średnia arytmetyczna wektorów kartezjańskich reprezentujących punkty.
• Polilinie : średnia ważona wektorów kartezjańskich reprezentujących punkty środkowe każdego segmentu linii, ważona przez (sferyczną) długość każdego segmentu.

Wydaje się, że rozsądną kontynuacją jest zdefiniowanie centroidów wielokątów jako średniej ważonej rozkładu trójkątnego, ważonej według powierzchni.

Moje pytanie brzmi, czy metoda z powyższego odniesienia będzie działać niezależnie od zastosowanego rozkładu trójkąta. W szczególności wspomina o rozkładzie na trójkąty względem dowolnego punktu, nawet zewnętrznego względem wielokąta, tak że niektóre trójkąty będą miały ujemne obszary, które mają ujemny ciężar.

Powiązane: Jak znaleźć środek geometrii obiektu?

Nie będzie działać konsekwentnie, nawet jeśli wykonasz wszystkie triangulacje w stosunku do jednego, stałego punktu. Problem polega na tym, że obliczenia sferyczne i euklidesowe są mieszane bez uwzględnienia ich znaczenia.

Jednym ze sposobów, aby uczynić to oczywistym, jest rozważenie dość ekstremalnego trójkąta, takiego jak prawie połowa półkuli. Na przykład, zaczynając od (lon, lat) = (-179, 0), biegnij wzdłuż równika do (0, 0), następnie do bieguna północnego w (0, 90), a następnie z powrotem do początku o (- 179, 0). Jest to trójkąt 90-179-90 obejmujący większość północnej połowy zachodniej półkuli. Problem polega na tym, że jego punkty końcowe (pokazane na rysunku jako białe kropki) leżą praktycznie w płaszczyźnie: jedna znajduje się na słupie, a pozostałe dwa są prawie po przeciwnych stronach. Tak więc ich średnia, rzutowana z powrotem na kulę (czerwona kropka), znajduje się prawie na biegunie - ale to jest tak daleko od jakiegokolwiek rozsądnego centrum, jak można uzyskać:

Jako inny przykład, triangulujmy wielokąt reprezentujący górną półkulę w stosunku do jej centrum, Bieguna Północnego. Zawsze podzielimy półkulę zachodnią na dwie równe połowy, z których każda ma trójkąt 90–90–90 (unikając w ten sposób problemów z dużymi trójkątami obejmującymi półkulę). Półkula wschodnia zostanie jednak podzielona na n równych pół-lun. Wierzchołki lune k ( k = 1, 2, ..., n ) mają współrzędne (lon, lat)

((k-1) * 180/n, 0),  (k * 180/n, 0),  (k * 180/n, 90).

Ten rysunek pokazuje konfigurację dla k = 8. Czerwone kropki to pojedyncze „centra” trójkąta obliczone zgodnie z dokumentem „Narzędzia do grafiki i kształtów”, s. 65–67.

Na podstawie obliczeń stwierdzam, że przy k = 2 środek ważony powierzchniowo rzeczywiście znajduje się na biegunie północnym (co byłoby wskazane w rozważaniach dotyczących symetrii), ale wraz ze wzrostem n wynik szybko przesuwa się na półkulę zachodnią, a na granica, zbliża się do szerokości geograficznej 89.556 stopni wzdłuż długości -90 stopni. Jest to około 50 kilometrów na południe od samego bieguna północnego.

Wprawdzie błąd +/- 50 kilometrów dla wielokąta rozciągającego się na 20 000 kilometrów jest niewielki; całkowita ilość dowolnych zmian spowodowanych różnymi triangulacjami w tym przypadku wynosi tylko 0,5%. Oczywiście względne błędy mogą być dowolnie duże przez włączenie ujemnych trójkątów (wystarczy dodać i odjąć kilka naprawdę dużych trójkątów w stosunku do małego trójkąta). Niezależnie od tego, każdy, kto stara się wykonywać obliczenia sferyczne, najwyraźniej stara się uniknąć błędów projekcji, więc poszukuje wysokiej dokładności. Ta metoda triangulacji nie może być zalecana.

Dobrym pomysłem jest wyliczenie właściwości, które powinna mieć centroid wielokąta. Oto moje kryteria:

(a) Jest to właściwość wnętrza wielokąta (zamiast wierzchołków lub krawędzi). Zatem podzielenie krawędzi na dwie części poprzez wstawienie dodatkowego wierzchołka nie powinno zmieniać położenia środka ciężkości. Zauważ, że definicja środka ciężkości według Jennessa nie spełnia tego kryterium, ponieważ położenie środka ciężkości będzie zależeć od podziału wielokąta na trójkąty.

(b) Lekkie poruszenie kształtu wielokąta powinno nieznacznie przesunąć środek ciężkości. Konieczne jest tutaj nałożenie ograniczenia na ogólny zasięg wielokąta (np. Na pojedynczą półkulę). Bez tego ograniczenia łatwo jest skonstruować przypadki, w których centroid nagle przesunie się na przeciwną stronę ziemi z lekkim ruchem wierzchołka. Ten warunek wyklucza metody wymagające, aby środek ciężkości leżał wewnątrz wielokąta.

(c) Powinien sprowadzać się do płaskiej definicji centroidu dla małych wielokątów.

Oto dwa podejścia spełniające te kryteria:

(1) Obliczyć środek ciężkości wielokąta elipsoidalnego w trzech wymiarach i rzutować z powrotem na powierzchnię elipsoidy (wzdłuż linii normalnej do elipsoidy). Duża zaleta: centroid można obliczyć, dzieląc wielokąt na prostsze kształty.

(2) Środek ciężkości jest punktem o minimalnej odległości geodezyjnej RMS do wszystkich punktów we wnętrzu wielokąta. Patrz Buss i Fillmore, „Średnie sferyczne i zastosowania do splajnów sferycznych i interpolacji”, Transakcje ACM na grafice 20 , 95–126 (2001). Duża zaleta: wynikowy punkt nie zależy od sposobu osadzenia powierzchni w R. ³ .

Niestety, żadna z tych definicji nie jest łatwa do wprowadzenia w życie. jednak pierwszą metodę można przeprowadzić po prostu dla kuli. Najlepszym do użycia „elementarnym” obszarem jest czworobok ograniczony krawędzią wielokąta, dwa południki przez punkty końcowe krawędzi i równik. Wynik dla całego wielokąta wymaga zsumowania wkładów na krawędziach. (Jeśli wielokąt otacza biegun), należy podjąć dodatkowe kroki.

Załóżmy, że punktami końcowymi krawędzi są (φ ₁ , λ ₁ ) i (φ ₂ , λ ₂ ). Niech azymuty krawędzi i punktów końcowych przez α ₁ i α ₂ . Zakładając, że promień kuli wynosi 1, powierzchnia czworoboku jest równa

  A = α ₂ - α ₁
      = 2 tan ⁻¹ [tan ½ (λ ₂ - λ ₁ ) sin ½ (φ ₂ + φ ₁ ) / cos ½ (φ ₂ + φ ₁ )]

(Ta formuła dla obszaru, ze względu na Bessela, jest znacznie lepiej zachowana numerycznie niż powszechnie stosowana formuła L'Huiliera dla powierzchni trójkąta.)

Składniki środka ciężkości dla tego czworoboku podano przez

  2 ⟨ x ⟩ = φ ₂ sin (λ ₂ - λ ₀ ) - φ ₁ sin (λ ₁ - λ ₀ )   2 ⟨ Y ⟩ = cos α ₀ (σ ₂ - σ ₁ ) - (φ ₂ cos (λ ₂ - λ ₀ ) - φ ₁ cos (λ ₁ - λ ₀ ))   2 ⟨ z ⟩ = (λ ₂ - λ ₁ ) - sin α ₀ (σ ₁ )

₂ - σ

gdzie σ ₂ - σ ₁ jest długością krawędzi i λ ₀ i α ₀ są długość i azymut krawędzi, gdzie przecina równik, a x i y osie są usytuowane tak, że przejście równik wynosi x = 1, y = 0. ( z jest oczywiście osią przechodzącą przez biegun).