Podsumowanie
Ta odpowiedź umieszcza pytanie w szerszym kontekście, opisuje efektywny algorytm mający zastosowanie do reprezentacji plików kształtów funkcji (jako „wektorów” lub „linii” punktów), pokazuje niektóre przykłady jego zastosowania i podaje kod roboczy do używania lub przenoszenia do środowisko GIS.
tło
To jest przykład rozszerzenia morfologicznego. Zasadniczo dylatacja „rozciąga” punkty regionu na ich sąsiedztwa; zbiór punktów, w których kończą, to „dylatacja”. Aplikacje w GIS są liczne: modelowanie rozprzestrzeniania się ognia, przemieszczania się cywilizacji, rozprzestrzeniania się roślin i wiele innych.
Matematycznie iw bardzo wielkiej (ale użytecznej) ogólności dylatacja rozkłada zbiór punktów w rozmaitości Riemanniana (np. Płaszczyzna, kula lub elipsoida). Rozłożenie jest określone przez podzbiór wiązki stycznej w tych punktach. Oznacza to, że w każdym z punktów podany jest zestaw wektorów (kierunków i odległości) (nazywam to „sąsiedztwem”); każdy z tych wektorów opisuje ścieżkę geodezyjną rozpoczynającą się w punkcie bazowym. Punkt bazowy jest „rozciągnięty” na końce wszystkich tych ścieżek. (Aby uzyskać znacznie bardziej ograniczoną definicję „dylatacji”, która jest tradycyjnie stosowana w przetwarzaniu obrazu, zobacz artykuł w Wikipedii . Funkcja rozprzestrzeniania jest znana jako mapa wykładnicza w geometrii różnicowej).
„Buforowanie” cechy jest jednym z najprostszych przykładów takiej dylatacji: dysk o stałym promieniu (promień bufora) jest tworzony (przynajmniej koncepcyjnie) wokół każdego punktu cechy. Połączenie tych dysków jest buforem.
To pytanie wymaga obliczenia nieco bardziej skomplikowanej dylatacji, w której rozprzestrzenianie jest dozwolone tylko w określonym zakresie kątów (to znaczy w sektorze kołowym). Ma to sens tylko w przypadku elementów, które nie obejmują wyraźnie zakrzywionej powierzchni (takich jak małe elementy na kuli lub elipsoidzie lub jakiekolwiek elementy w płaszczyźnie). Kiedy pracujemy w płaszczyźnie, sensowne jest również zorientowanie wszystkich sektorów w tym samym kierunku. (Gdybyśmy jednak modelowali rozprzestrzenianie się ognia przez wiatr, chcielibyśmy, aby sektory były zorientowane na wiatr, a ich rozmiary również mogą się różnić w zależności od prędkości wiatru: to jedna z motywacji do ogólnej definicji dylatacji, którą podałem. ) (Na zakrzywionych powierzchniach, takich jak elipsoida, ogólnie niemożliwe jest zorientowanie wszystkich sektorów w „tym samym” kierunku).
W następujących okolicznościach dylatacja jest stosunkowo łatwa do obliczenia:
Obiekt znajduje się w płaszczyźnie (czyli rozszerzamy mapę obiektu i mam nadzieję, że mapa jest dość dokładna).
Dylatacja będzie stała : rozprzestrzenianie się w każdym punkcie cechy nastąpi w przystających dzielnicach o tej samej orientacji.
To wspólne sąsiedztwo jest wypukłe. Wypukłość znacznie upraszcza i przyspiesza obliczenia.
To pytanie mieści się w takich wyspecjalizowanych okolicznościach: prosi o dylatację dowolnych wielokątów przez sektory kołowe, których początki (środki dysków, z których pochodzą) znajdują się w punktach bazowych. Pod warunkiem, że sektory te nie przekraczają 180 stopni, będą wypukłe. (Większe sektory można zawsze podzielić na pół na dwa wypukłe sektory; połączenie dwóch mniejszych rozszerzeń da pożądany rezultat.)
Realizacja
Ponieważ wykonujemy euklidesowe obliczeń - robi rozprzestrzenia się w samolocie - możemy rozszerzać punkt jedynie poprzez tłumaczenia sąsiedztwo dylatacje do tego punktu. (Aby to zrobić, sąsiedztwo potrzebuje pochodzeniaktóry będzie odpowiadał punktowi bazowemu. Na przykład źródłem sektorów w tym pytaniu jest środek koła, z którego zostały utworzone. To pochodzenie dzieje się na granicy sektora. W standardowej operacji buforowania GIS sąsiedztwo jest kołem, którego początkiem jest środek; teraz pochodzenie leży we wnętrzu koła. Wybór źródła nie jest wielką sprawą obliczeniową, ponieważ zmiana pochodzenia jedynie przesuwa całą dylatację, ale może to być duża sprawa w zakresie modelowania zjawisk naturalnych. sectorFunkcja w kodzie poniżej ilustruje, w jaki sposób można określić pochodzenie).
Dylatacja odcinka linii może być trudna, ale dla wypukłego sąsiedztwa możemy stworzyć dylatację jako połączenie dylatacji dwóch punktów końcowych wraz ze starannie wybranym równoległobokiem. (W interesie przestrzeni kosmicznej nie zatrzymam się, aby udowodnić takie twierdzenia matematyczne, ale zachęcę czytelników do wypróbowania własnych dowodów, ponieważ jest to wnikliwe ćwiczenie.) Oto ilustracja wykorzystująca trzy sektory (pokazane na różowo). Mają promienie jednostkowe, a ich kąty podano w tytułach. Sam segment linii ma długość 2, jest poziomy i jest pokazany na czarno:
Równolegogramy można znaleźć, lokalizując różowe punkty, które są jak najdalej od segmentu tylko w kierunku pionowym . Daje to dwa dolne punkty i dwa górne punkty wzdłuż linii równoległych do segmentu. Musimy tylko połączyć cztery punkty w równoległobok (pokazany na niebiesko). Zwróć uwagę, po prawej, jak to ma sens, nawet jeśli sam sektor jest tylko segmentem liniowym (a nie prawdziwym wielokątem): tam każdy punkt na tym segmencie został przetłumaczony w kierunku 171 stopni na wschód od północy na odległość w zakresie od 0 do 1. Zestawem tych punktów końcowych jest pokazany równoległobok. Szczegóły tego obliczenia pojawiają się w bufferfunkcji zdefiniowanej dilate.edgesw poniższym kodzie.
Aby rozszerzyć polilinię , tworzymy połączenie rozszerzeń punktów i segmentów, które ją tworzą. Dwie ostatnie linie dilate.edgeswykonują tę pętlę.
Dylatacja wielokąta wymaga uwzględnienia jego wnętrza wraz z rozszerzeniem jego granicy. (To twierdzenie przyjmuje pewne założenia dotyczące sąsiedztwa dylatacji. Jednym z nich jest to, że wszystkie sąsiedztwa zawierają punkt (0,0), co gwarantuje, że wielokąt jest objęty dylatacją. W przypadku dzielnic zmiennych zakłada również, że dylatacja dowolnego wnętrza punkt wielokąta nie będzie rozciągał się poza rozszerzeniem punktów granicznych. Dotyczy to stałych sąsiedztw.)
Spójrzmy na kilka przykładów tego, jak to działa, najpierw z nonagonem (wybranym w celu ujawnienia szczegółów), a następnie z kółkiem (wybranym w celu dopasowania do ilustracji w pytaniu). W przykładach nadal będą używane te same trzy dzielnice, ale zmniejszone do promienia 1/3.
Na tej figurze wnętrze wielokąta jest szare, rozszerzenia punktowe (sektory) są różowe, a rozszerzenia krawędzi (równoległoboki) są niebieskie.
„Koło” to tak naprawdę zaledwie 60-gram, ale ładnie zbliża się do koła.
Wydajność
Gdy element podstawowy jest reprezentowany przez N punktów, a sąsiedztwo dylatacji przez M punktów, algorytm ten wymaga wysiłku O (N M) . Następnie należy uprościć bałagan wierzchołków i krawędzi w złączu, co może wymagać wysiłku O (N M log (N M)): o to należy poprosić GIS; nie powinniśmy tego programować.
Wysiłek obliczeniowy można poprawić do O (M + N) dla wypukłych elementów bazowych (ponieważ można dowiedzieć się, jak poruszać się po nowej granicy, odpowiednio łącząc listy wierzchołków opisujących granice dwóch oryginalnych kształtów). Nie wymagałoby to również późniejszego czyszczenia.
Gdy sąsiedztwo dylatacji powoli zmienia rozmiar i / lub orientację w miarę przesuwania się wokół elementu podstawowego, rozszerzenie dylatacji krawędzi można ściśle przybliżyć od wypukłego kadłuba połączenia dylatacji jego punktów końcowych. Jeśli dwa sąsiedztwa dylatacyjne mają punkty M1 i M2, można to znaleźć przy wysiłku O (M1 + M2) przy użyciu algorytmu opisanego w Shamos i Preparata, Geometria obliczeniowa . Dlatego pozwalając K = M1 + M2 + ... + M (N) być całkowitą liczbą wierzchołków w sąsiedztwie dylatacji N, możemy obliczyć dylatację w czasie O (K * log (K)).
Dlaczego mielibyśmy zmierzyć się z takim uogólnieniem, jeśli wszystko, czego chcemy, to prosty bufor? W przypadku dużych obiektów na ziemi sąsiedztwo dylatacji (takie jak dysk), które w rzeczywistości ma stały rozmiar, może mieć różną wielkość na mapie, na której wykonywane są te obliczenia. W ten sposób uzyskujemy sposób wykonywania obliczeń, które są dokładne dla elipsoidy, przy jednoczesnym korzystaniu ze wszystkich zalet geometrii euklidesowej.
Kod
Przykłady zostały opracowane przy użyciu tego Rprototypu, który można łatwo przenieść na ulubiony język (Python, C ++ itp.). Pod względem struktury przypomina analizę opisaną w tej odpowiedzi i dlatego nie wymaga osobnego wyjaśnienia. Komentarze wyjaśniają niektóre szczegóły.
(Ciekawe może być to, że obliczenia trygonometryczne są używane tylko do tworzenia przykładowych cech - które są regularnymi wielokątami - i sektorów. Żadna część obliczeń dylatacyjnych nie wymaga żadnej trygonometrii).
#
# Dilate the vertices of a polygon/polyline by a shape.
#
dilate.points <- function(p, q) {
# Translate a copy of `q` to each vertex of `p`, resulting in a list of polygons.
pieces <- apply(p, 1, function(x) list(t(t(q)+x)))
lapply(pieces, function(z) z[[1]]) # Convert to a list of matrices
}
#
# Dilate the edges of a polygon/polyline `p` by a shape `q`.
# `p` must have at least two rows.
#
dilate.edges <- function(p, q) {
i <- matrix(c(0,-1,1,0), 2, 2) # 90 degree rotation
e <- apply(rbind(p, p[1,]), 2, diff) # Direction vectors of the edges
# Dilate a single edge from `x` to `x+v` into a parallelogram
# bounded by parts of the dilation shape that are at extreme distances
# from the edge.
buffer <- function(x, v) {
y <- q %*% i %*% v # Signed distances orthogonal to the edge
k <- which.min(y) # Find smallest distance, then the largest *after* it
l <- (which.max(c(y[-(1:k)], y[1:k])) + k-1) %% length(y)[1] + 1
list(rbind(x+q[k,], x+v+q[k,], x+v+q[l,], x+q[l,])) # A parallelogram
}
# Apply `buffer` to every edge.
quads <- apply(cbind(p, e), 1, function(x) buffer(x[1:2], x[3:4]))
lapply(quads, function(z) z[[1]]) # Convert to a list of matrices
}
#----------------------- (This ends the dilation code.) --------------------------#
#
# Display a polygon and its point and edge dilations.
# NB: In practice we would submit the polygon, its point dilations, and edge
# dilations to the GIS to create and simplify their union, producing a single
# polygon. We keep the three parts separate here in order to illustrate how
# that polygon is constructed.
#
display <- function(p, d.points, d.edges, ...) {
# Create a plotting region covering the extent of the dilated figure.
x <- c(p[,1], unlist(lapply(c(d.points, d.edges), function(x) x[,1])))
y <- c(p[,2], unlist(lapply(c(d.points, d.edges), function(x) x[,2])))
plot(c(min(x),max(x)), c(min(y),max(y)), type="n", asp=1, xlab="x", ylab="y", ...)
# The polygon itself.
polygon(p, density=-1, col="#00000040")
# The dilated points and edges.
plot.list <- function(l, c) lapply(l, function(p)
polygon(p, density=-1, col=c, border="#00000040"))
plot.list(d.points, "#ff000020")
plot.list(d.edges, "#0000ff20")
invisible(NULL) # Doesn't return anything
}
#
# Create a sector of a circle.
# `n` is the number of vertices to use for approximating its outer arc.
#
sector <- function(radius, arg1, arg2, n=1, origin=c(0,0)) {
t(cbind(origin, radius*sapply(seq(arg1, arg2, length.out=n),
function(a) c(cos(a), sin(a)))))
}
#
# Create a polygon represented as an array of rows.
#
n.vertices <- 60 # Inscribes an `n.vertices`-gon in the unit circle.
angles <- seq(2*pi, 0, length.out=n.vertices+1)
angles <- angles[-(n.vertices+1)]
polygon.the <- cbind(cos(angles), sin(angles))
if (n.vertices==1) polygon.the <- rbind(polygon.the, polygon.the)
#
# Dilate the polygon in various ways to illustrate.
#
system.time({
radius <- 1/3
par(mfrow=c(1,3))
q <- sector(radius, pi/12, 2*pi/3, n=120)
d.points <- dilate.points(polygon.the, q)
d.edges <- dilate.edges(polygon.the, q)
display(polygon.the, d.points, d.edges, main="-30 to 75 degrees")
q <- sector(radius, pi/3, 4*pi/3, n=180)
d.points <- dilate.points(polygon.the, q)
d.edges <- dilate.edges(polygon.the, q)
display(polygon.the, d.points, d.edges, main="-150 to 30 degrees")
q <- sector(radius, -9/20*pi, -9/20*pi)
d.points <- dilate.points(polygon.the, q)
d.edges <- dilate.edges(polygon.the, q)
display(polygon.the, d.points, d.edges, main="171 degrees")
})
Czas obliczeniowy dla tego przykładu (z ostatniej cyfry), przy N = 60 i M = 121 (po lewej), M = 181 (po środku) i M = 2 (po prawej), wynosił kwadrans. Jednak większość z nich dotyczyła wyświetlacza. Zazwyczaj ten Rkod będzie obsługiwał około N M = 1,5 miliona na sekundę (zajmuje to tylko 0,002 sekundy, aby wykonać wszystkie przedstawione przykładowe obliczenia). Niemniej jednak pojawienie się produktu M N oznacza rozszerzenie wielu postaci lub skomplikowanych postaci w szczegółowym sąsiedztwie, może zająć sporo czasu, więc uważaj! Porównać czas dla mniejszych problemów przed rozwiązaniem dużego. W takich okolicznościach można spojrzeć na rozwiązanie oparte na rastrze (które jest znacznie łatwiejsze do wdrożenia, wymagając zasadniczo tylko jednego obliczenia sąsiedztwa).
