Widzę MerseyViking zaleciła QuadTree . Chciałem zasugerować to samo i aby to wyjaśnić, oto kod i przykład. Kod jest zapisany, Rale powinien być łatwy do przeniesienia, powiedzmy, do Pythona.
Pomysł jest niezwykle prosty: podziel punkty w przybliżeniu na pół w kierunku x, a następnie rekurencyjnie podziel dwie połówki wzdłuż kierunku y, zmieniając kierunki na każdym poziomie, aż nie będzie już konieczne dzielenie.
Ponieważ celem jest ukrycie rzeczywistych lokalizacji punktów, przydatne jest wprowadzenie losowości do podziałów . Jednym szybkim i prostym sposobem na dokonanie tego jest podzielenie na zestaw kwantyli małej losowej kwoty od 50%. W ten sposób (a) jest mało prawdopodobne, aby wartości podziału były zbieżne ze współrzędnymi danych, tak że punkty będą jednoznacznie wpadać do kwadrantów utworzonych przez podział, oraz (b) współrzędnych punktów nie będzie można dokładnie odtworzyć z kwadratu.
Ponieważ intencją jest utrzymanie minimalnej kliczby węzłów w każdym liściu czterokąta, implementujemy ograniczoną formę czterokąta. Będzie obsługiwał (1) punkty grupowania w grupy posiadające pomiędzy ki 2 * k-1 każdy element oraz (2) mapowanie kwadrantów.
Ten Rkod tworzy drzewo węzłów i liści końcowych, rozróżniając je według klasy. Etykietowanie klas przyspiesza przetwarzanie końcowe, takie jak drukowanie, przedstawione poniżej. Kod używa wartości liczbowych dla identyfikatorów. Działa to do głębokości 52 w drzewie (przy użyciu podwójnych; jeśli używane są długie liczby całkowite bez znaku, maksymalna głębokość wynosi 32). W przypadku głębszych drzew (które są bardzo mało prawdopodobne w żadnej aplikacji, ponieważ w grę wchodzą co najmniej k* 2 ^ 52 punkty), identyfikatory musiałyby być łańcuchami.
quadtree <- function(xy, k=1) {
d = dim(xy)[2]
quad <- function(xy, i, id=1) {
if (length(xy) < 2*k*d) {
rv = list(id=id, value=xy)
class(rv) <- "quadtree.leaf"
}
else {
q0 <- (1 + runif(1,min=-1/2,max=1/2)/dim(xy)[1])/2 # Random quantile near the median
x0 <- quantile(xy[,i], q0)
j <- i %% d + 1 # (Works for octrees, too...)
rv <- list(index=i, threshold=x0,
lower=quad(xy[xy[,i] <= x0, ], j, id*2),
upper=quad(xy[xy[,i] > x0, ], j, id*2+1))
class(rv) <- "quadtree"
}
return(rv)
}
quad(xy, 1)
}
Należy zauważyć, że rekurencyjna konstrukcja tego algorytmu typu „dziel i rządź” (i w konsekwencji większość algorytmów przetwarzania końcowego) oznacza, że wymagany czas to O (m), a użycie pamięci RAM to O (n), gdzie mjest liczba komórki i njest liczbą punktów. mjest proporcjonalne don podzielonej przez minimalną punktów na komórkę,k. Jest to przydatne do szacowania czasów obliczeń. Na przykład, jeśli podzielenie n = 10 ^ 6 punktów na komórki 50-99 punktów (k = 50) zajmuje 13 sekund, m = 10 ^ 6/50 = 20000. Jeśli zamiast tego chcesz podzielić na 5-9 punktów na komórkę (k = 5), m jest 10 razy większy, więc czas wzrasta do około 130 sekund. (Ponieważ proces dzielenia zestawu współrzędnych wokół ich środków przyspiesza, gdy komórki stają się mniejsze, faktyczny czas wyniósł zaledwie 90 sekund.) Aby przejść do k = 1 punkt na komórkę, zajmie to około sześć razy dłużej jeszcze dziewięć minut i możemy spodziewać się, że kod będzie trochę szybszy.
Zanim przejdziemy dalej, wygenerujmy interesujące dane o nieregularnych odstępach i stwórzmy ograniczony quadree (czas, który upłynął 0,29 sekundy):
Oto kod do tworzenia tych wykresów. Wykorzystuje Rpolimorfizm: points.quadtreebędzie wywoływany za każdym razem, gdy pointsfunkcja zostanie zastosowana na przykład do quadtreeobiektu. Siła tego jest widoczna w ekstremalnej prostocie funkcji pokolorowania punktów zgodnie z ich identyfikatorem skupienia:
points.quadtree <- function(q, ...) {
points(q$lower, ...); points(q$upper, ...)
}
points.quadtree.leaf <- function(q, ...) {
points(q$value, col=hsv(q$id), ...)
}
Wykreślenie samej siatki jest nieco trudniejsze, ponieważ wymaga wielokrotnego przycinania progów użytych do partycjonowania quadtree, ale to samo rekurencyjne podejście jest proste i eleganckie. W razie potrzeby użyj wariantu, aby zbudować wielokątne reprezentacje ćwiartek.
lines.quadtree <- function(q, xylim, ...) {
i <- q$index
j <- 3 - q$index
clip <- function(xylim.clip, i, upper) {
if (upper) xylim.clip[1, i] <- max(q$threshold, xylim.clip[1,i]) else
xylim.clip[2,i] <- min(q$threshold, xylim.clip[2,i])
xylim.clip
}
if(q$threshold > xylim[1,i]) lines(q$lower, clip(xylim, i, FALSE), ...)
if(q$threshold < xylim[2,i]) lines(q$upper, clip(xylim, i, TRUE), ...)
xlim <- xylim[, j]
xy <- cbind(c(q$threshold, q$threshold), xlim)
lines(xy[, order(i:j)], ...)
}
lines.quadtree.leaf <- function(q, xylim, ...) {} # Nothing to do at leaves!
Jako inny przykład wygenerowałem 1 000 000 punktów i podzieliłem je na grupy po 5-9 sztuk. Czas wyniósł 91,7 sekundy.
n <- 25000 # Points per cluster
n.centers <- 40 # Number of cluster centers
sd <- 1/2 # Standard deviation of each cluster
set.seed(17)
centers <- matrix(runif(n.centers*2, min=c(-90, 30), max=c(-75, 40)), ncol=2, byrow=TRUE)
xy <- matrix(apply(centers, 1, function(x) rnorm(n*2, mean=x, sd=sd)), ncol=2, byrow=TRUE)
k <- 5
system.time(qt <- quadtree(xy, k))
#
# Set up to map the full extent of the quadtree.
#
xylim <- cbind(x=c(min(xy[,1]), max(xy[,1])), y=c(min(xy[,2]), max(xy[,2])))
plot(xylim, type="n", xlab="x", ylab="y", main="Quadtree")
#
# This is all the code needed for the plot!
#
lines(qt, xylim, col="Gray")
points(qt, pch=".")
Jako przykład interakcji z GIS , wypiszmy wszystkie komórki z czterema drzewami jako plik kształtu wielokąta za pomocą shapefilesbiblioteki. Kod emuluje procedury obcinania lines.quadtree, ale tym razem musi wygenerować opisy wektorowe komórek. Są one wyprowadzane jako ramki danych do użytku z shapefilesbiblioteką.
cell <- function(q, xylim, ...) {
if (class(q)=="quadtree") f <- cell.quadtree else f <- cell.quadtree.leaf
f(q, xylim, ...)
}
cell.quadtree <- function(q, xylim, ...) {
i <- q$index
j <- 3 - q$index
clip <- function(xylim.clip, i, upper) {
if (upper) xylim.clip[1, i] <- max(q$threshold, xylim.clip[1,i]) else
xylim.clip[2,i] <- min(q$threshold, xylim.clip[2,i])
xylim.clip
}
d <- data.frame(id=NULL, x=NULL, y=NULL)
if(q$threshold > xylim[1,i]) d <- cell(q$lower, clip(xylim, i, FALSE), ...)
if(q$threshold < xylim[2,i]) d <- rbind(d, cell(q$upper, clip(xylim, i, TRUE), ...))
d
}
cell.quadtree.leaf <- function(q, xylim) {
data.frame(id = q$id,
x = c(xylim[1,1], xylim[2,1], xylim[2,1], xylim[1,1], xylim[1,1]),
y = c(xylim[1,2], xylim[1,2], xylim[2,2], xylim[2,2], xylim[1,2]))
}
Same punkty można odczytać bezpośrednio, używając read.shplub importując plik danych o współrzędnych (x, y).
Przykład zastosowania:
qt <- quadtree(xy, k)
xylim <- cbind(x=c(min(xy[,1]), max(xy[,1])), y=c(min(xy[,2]), max(xy[,2])))
polys <- cell(qt, xylim)
polys.attr <- data.frame(id=unique(polys$id))
library(shapefiles)
polys.shapefile <- convert.to.shapefile(polys, polys.attr, "id", 5)
write.shapefile(polys.shapefile, "f:/temp/quadtree", arcgis=TRUE)
(Użyj tutaj dowolnego pożądanego zasięgu, aby xylimprzejść do podregionu lub rozszerzyć mapowanie na większy region; ten kod domyślnie określa zakres punktów).
Już samo to wystarczy: przestrzenne połączenie tych wielokątów z pierwotnymi punktami zidentyfikuje skupiska. Po zidentyfikowaniu operacje „podsumowania” w bazie danych wygenerują statystyki podsumowujące punkty w każdej komórce.
