Jeden szybki i brudny sposób wykorzystuje rekurencyjny podział sferyczny . Zaczynając od triangulacji powierzchni ziemi, rekurencyjnie rozdzielaj każdy trójkąt od wierzchołka do środka jego najdłuższego boku. (Idealnie podzielisz trójkąt na dwie części o równej średnicy lub części o równej powierzchni, ale ponieważ wymagają one skomplikowanych obliczeń, po prostu podzielę boki dokładnie na pół: powoduje to, że różne trójkąty ostatecznie różnią się nieco wielkością, ale nie wydaje się to krytyczne dla tej aplikacji).
Oczywiście utrzymasz ten podział w strukturze danych, która pozwala szybko zidentyfikować trójkąt, w którym leży dowolny dowolny punkt. Drzewo binarne (oparte na wywołaniach rekurencyjnych) działa dobrze: za każdym razem, gdy trójkąt jest dzielony, drzewo jest dzielone w węźle tego trójkąta. Dane dotyczące płaszczyzny podziału są zachowywane, dzięki czemu można szybko określić, po której stronie płaszczyzny leży dowolny dowolny punkt: który określa, czy należy podróżować w lewo, czy w prawo po drzewie.
(Czy mówiłem o dzieleniu „płaszczyzny”? Tak - jeśli modelujemy powierzchnię ziemi jako kulę i używamy współrzędnych geocentrycznych (x, y, z), wówczas większość naszych obliczeń odbywa się w trzech wymiarach, gdzie boki trójkątów są kawałkami przecięcia sfery z płaszczyznami poprzez jej pochodzenie. Dzięki temu obliczenia są szybkie i łatwe).
Zilustruję to pokazując procedurę na jednej oktanie kuli; pozostałe siedem oktantów jest przetwarzanych w ten sam sposób. Taki oktan to trójkąt 90–90–90. W mojej grafice narysuję trójkąty euklidesowe obejmujące te same rogi: nie wyglądają zbyt dobrze, dopóki nie staną się małe, ale można je łatwo i szybko narysować. Oto trójkąt euklidesowy odpowiadający oktanowi: jest to początek procedury.
Ponieważ wszystkie boki mają równą długość, jedna jest wybierana losowo jako „najdłuższa” i dzieli się na:
Powtórz to dla każdego z nowych trójkątów:
Po n krokach będziemy mieli 2 ^ n trójkątów. Oto sytuacja po 10 krokach, pokazująca 1024 trójkąty w oktanie (i 8192 na kuli ogółem):
Jako kolejną ilustrację wygenerowałem losowy punkt w tej oktanie i przemierzyłem drzewo podziału, aż najdłuższy bok trójkąta osiągnie mniej niż 0,05 radianów. Trójkąty (kartezjańskie) są oznaczone kolorem punktu próbkowania na czerwono.
Nawiasem mówiąc, aby zawęzić położenie punktu do jednego stopnia szerokości geograficznej (w przybliżeniu), należy zauważyć, że jest to około 1/60 radianów, a więc obejmuje około (1/60) ^ 2 / (Pi / 2) = 1/6000 z powierzchnia całkowita. Ponieważ każdy poddział dzieli w przybliżeniu o połowę rozmiar trójkąta, około 13 do 14 pododdziałów oktanu załatwi sprawę. To niewiele obliczeń - jak zobaczymy poniżej - czyniąc efektywnym nie przechowywanie drzewa w ogóle, ale dokonywanie podziału w locie. Na początku należy zauważyć, w której ósemce znajduje się punkt - który jest określony przez znaki jego trzech współrzędnych, które można zapisać jako trzycyfrową liczbę dwójkową - i na każdym kroku chcesz pamiętać, czy punkt leży w lewej (0) lub prawej (1) trójkąta. To daje kolejną 14-cyfrową liczbę binarną. Możesz użyć tych kodów do grupowania dowolnych punktów.
(Zasadniczo, gdy dwa kody są bliskie jak rzeczywiste liczby binarne, odpowiednie punkty są bliskie; jednak punkty mogą nadal być bliskie i mieć wyjątkowo różne kody. Weźmy na przykład dwa punkty oddzielone równikiem o jeden metr: ich kody muszą się różnić nawet przed punktem binarnym, ponieważ są w różnych oktantach. Tego rodzaju rzeczy nie da się uniknąć przy żadnym stałym podziale przestrzeni.)
Użyłem Mathematica 8, aby to zaimplementować: możesz wziąć to tak, jak jest lub jako pseudokod do implementacji w twoim ulubionym środowisku programistycznym.
Określ, na której stronie płaszczyzny punkt 0-ab p leży:
side[p_, {a_, b_}] := If[Det[{p, a, b}] >= 0, left, right];
Udoskonal trójkąt abc na podstawie punktu p.
refine[p_, {a_, b_, c_}] := Block[{sides, x, y, z, m},
sides = Norm /@ {b - c, c - a, a - b} // N;
{x, y, z} = RotateLeft[{a, b, c}, First[Position[sides, Max[sides]]] - 1];
m = Normalize[Mean[{y, z}]];
If[side[p, {x, m}] === right, {y, m, x}, {x, m, z}]
]
Ostatnia cyfra została narysowana przez wyświetlenie oktanta, a ponadto przez renderowanie poniższej listy jako zestawu wielokątów:
p = Normalize@RandomReal[NormalDistribution[0, 1], 3] (* Random point *)
{a, b, c} = IdentityMatrix[3] . DiagonalMatrix[Sign[p]] // N (* First octant *)
NestWhileList[refine[p, #] &, {a, b, c}, Norm[#[[1]] - #[[2]]] >= 0.05 &, 1, 16]
NestWhileListwielokrotnie stosuje operację ( refine), dopóki spełniony jest warunek (trójkąt jest duży) lub do momentu osiągnięcia maksymalnej liczby operacji (16).
Aby wyświetlić pełną triangulację oktanu, zacząłem od pierwszego oktanu i powtórzyłem wyrafinowanie dziesięć razy. Zaczyna się to od niewielkiej modyfikacji refine:
split[{a_, b_, c_}] := Module[{sides, x, y, z, m},
sides = Norm /@ {b - c, c - a, a - b} // N;
{x, y, z} = RotateLeft[{a, b, c}, First[Position[sides, Max[sides]]] - 1];
m = Normalize[Mean[{y, z}]];
{{y, m, x}, {x, m, z}}
]
Różnica polega na tym, że splitzwraca obie połówki trójkąta wejściowego, a nie tę, w której leży dany punkt. Pełna triangulacja jest uzyskiwana przez iterację:
triangles = NestList[Flatten[split /@ #, 1] &, {IdentityMatrix[3] // N}, 10];
Aby to sprawdzić, obliczyłem miarę każdego trójkąta i spojrzałem na zakres. (Ten „rozmiar” jest proporcjonalny do figury w kształcie piramidy, znajdującej się w każdym trójkącie i środku kuli; w przypadku takich małych trójkątów rozmiar ten jest zasadniczo proporcjonalny do jego sferycznego obszaru).
Through[{Min, Max}[Map[Round[Det[#], 0.00001] &, triangles[[10]] // N, {1}]]]
{0,00523, 0,00739}
Zatem rozmiary różnią się w górę lub w dół o około 25% od ich średniej: wydaje się to rozsądne, aby osiągnąć w przybliżeniu jednolity sposób grupowania punktów.
Podczas skanowania tego kodu nie zauważysz żadnej trygonometrii : jedynym miejscem, w którym będzie ono potrzebne, będzie w ogóle konwersja w obie strony między współrzędnymi sferycznymi i kartezjańskimi. Kod nie rzutuje również powierzchni ziemi na żadną mapę, co pozwala uniknąć towarzyszących jej zniekształceń. W przeciwnym razie do wykonania całej pracy używa tylko uśredniania ( Mean), twierdzenia Pitagorasa ( Norm) i wyznacznika 3 na 3 ( Det). (Istnieje kilka prostych poleceń list-manipulacja, jak RotateLefti Flattenteż, wraz z poszukiwaniem najdłuższego boku każdego trójkąta).