Jak zmierzyć podobieństwo obiektów SpatialLines

9

Stworzyłem dwa SpatialLinesobiekty w R: postać.

Te obiekty zostały utworzone w ten sposób:

library(sp)
xy <- cbind(x,y)
xy.sp = sp::SpatialPoints(xy)
spl1 <- sp::SpatialLines(list(Lines(Line(xy.sp), ID="a")))

Teraz chcę jakoś dojść do wniosku, że jest to ta sama linia obrócona i odwrócona, a różnica między nimi jest równa 0 (tzn. Kształt jest równy).

Aby to zrobić, można użyć maptoolspaczki i obrócić linię nr 1, np .:

spl180 <- maptools::elide(spl1, rotate=180)

Każdą obróconą linię należy następnie sprawdzić w stosunku do linii nr 2 za pomocą rgeospakietu, np .:

hdist <- rgeos::gDistance(spl180, spl2, byid=FALSE, hausdorff=TRUE)

Jest to jednak tak drogi pod względem obliczeniowym sposób dopasowywania SpatialLinesobiektów, zwłaszcza jeśli liczba obiektów wynosi około 1000.

Czy jest jakiś sprytny sposób na wykonanie tej pracy?

PS Ponadto wyżej opisane podejście nie gwarantuje wszystkich możliwych obrotów i przerzutów.

P.S2. Jeśli linia nr 1 jest pomniejszona w stosunku do linii nr 2, różnica między liniami nr 1 i nr 2 musi być nadal równa 0.

AKTUALIZACJA:

wprowadź opis zdjęcia tutaj

Klausos Klausos
źródło

Odpowiedzi:

9

Każda skuteczna metoda prawdziwie ogólnego zastosowania ustandaryzuje reprezentacje kształtów , aby nie ulegały zmianie po obrocie, translacji, odbiciu lub trywialnych zmianach reprezentacji wewnętrznej.

Jednym ze sposobów na to jest wyszczególnienie każdego połączonego kształtu jako naprzemiennej sekwencji długości krawędzi i (podpisanych) kątów, zaczynając od jednego końca. (Kształt powinien być „czysty” w tym sensie, że nie ma krawędzi zerowych ani prostych kątów.) Aby ten niezmiennik był odbijany, zaneguj wszystkie kąty, jeśli pierwszy niezerowy jest ujemny.

(Ponieważ jakakolwiek połączona polilinia n wierzchołków będzie miała n- 1 krawędzi oddzielonych n- kątami 2, w Rponiższym kodzie uznałem za wygodne użycie struktury danych składającej się z dwóch tablic, jednej dla długości krawędzi, $lengthsa drugiej dla kąty $angles,. Segment linii w ogóle nie będzie miał kątów, dlatego ważne jest, aby obsługiwać tablice o zerowej długości w takiej strukturze danych.)

Takie przedstawienia można uporządkować leksykograficznie. Należy uwzględnić pewne błędy zmiennoprzecinkowe nagromadzone podczas procesu normalizacji. Elegancka procedura oszacowałaby te błędy w funkcji pierwotnych współrzędnych. W poniższym rozwiązaniu zastosowano prostszą metodę, w której dwie długości są uważane za równe, gdy różnią się względnie bardzo małą ilością . Kąty mogą się różnić tylko bardzo małą wartością bezwzględną.

Aby stały się niezmiennikami przy odwróceniu orientacji leżącej u podstaw, wybierz najwcześniejsze leksykograficznie przedstawienie między polilinią a jej odwróceniem.

Aby obsługiwać polilinie wieloczęściowe, ustaw ich komponenty w porządku leksykograficznym.

Aby znaleźć klasy równoważności w transformacjach euklidesowych ,

  • Twórz standardowe reprezentacje kształtów.

  • Wykonaj leksykograficzny typ znormalizowanych reprezentacji.

  • Przejdź przez posortowaną kolejność, aby zidentyfikować sekwencje równych reprezentacji.

Czas obliczania jest proporcjonalny do O (n * log (n) * N), gdzie n to liczba cech, a N to największa liczba wierzchołków w dowolnej funkcji. To jest wydajne.

Prawdopodobnie warto wspomnieć mimochodem, że wstępne grupowanie oparte na łatwo obliczalnych niezmiennych właściwościach geometrycznych, takich jak długość polilinii, środek i momenty wokół tego środka, często można zastosować w celu usprawnienia całego procesu. Wystarczy znaleźć podgrupy przystających funkcji w każdej takiej wstępnej grupie. Podana tutaj pełna metoda byłaby potrzebna dla kształtów, które w przeciwnym razie byłyby tak niezwykle podobne, że takie proste niezmienniki wciąż by ich nie rozróżniały. Proste funkcje zbudowane z danych rastrowych mogą mieć na przykład takie cechy. Ponieważ jednak podane tutaj rozwiązanie jest tak wydajne, że jeśli ktoś podejmie wysiłek jego wdrożenia, może samo w sobie działać dobrze.


Przykład

Rysunek po lewej stronie pokazuje pięć polilinii plus 15 innych, które zostały uzyskane z nich przez losowe tłumaczenie, obrót, odbicie i odwrócenie orientacji wewnętrznej (która nie jest widoczna). Prawa figura koloruje je zgodnie z klasą równoważności euklidesowej: wszystkie figury tego samego koloru są zgodne; różne kolory nie są zgodne.

Postać

Rkod następuje. Kiedy dane wejściowe zostały zaktualizowane do 500 kształtów, 500 dodatkowych (przystających) kształtów, ze średnią 100 wierzchołków na kształt, czas wykonania na tym komputerze wynosił 3 sekundy.

Ten kod jest niekompletny: ponieważ Rnie ma rodzimego sortowania leksykograficznego i nie miałem ochoty kodować go od zera, po prostu wykonuję sortowanie według pierwszej współrzędnej każdego znormalizowanego kształtu. Będzie to dobrze w przypadku tworzonych tutaj losowych kształtów, ale do prac produkcyjnych należy wdrożyć pełny rodzaj leksykograficzny. Ta funkcja order.shapemiałaby wpływ tylko na tę funkcję . Jego dane wejściowe to lista znormalizowanych kształtów, sa dane wyjściowe to sekwencja indeksów, sktóre by je posortowały.

#
# Create random shapes.
#
n.shapes <- 5      # Unique shapes, up to congruence
n.shapes.new <- 15 # Additional congruent shapes to generate
p.mean <- 5        # Expected number of vertices per shape
set.seed(17)       # Create a reproducible starting point
shape.random <- function(n) matrix(rnorm(2*n), nrow=2, ncol=n)
shapes <- lapply(2+rpois(n.shapes, p.mean-2), shape.random)
#
# Randomly move them around.
#
move.random <- function(xy) {
  a <- runif(1, 0, 2*pi)
  reflection <- sign(runif(1, -1, 1))
  translation <- runif(2, -8, 8)
  m <- matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2) %*%
    matrix(c(reflection, 0, 0, 1), 2, 2)
  m <- m %*% xy + translation
  if (runif(1, -1, 0) < 0) m <- m[ ,dim(m)[2]:1]
  return (m)
}
i <- sample(length(shapes), n.shapes.new, replace=TRUE)
shapes <- c(shapes, lapply(i, function(j) move.random(shapes[[j]])))
#
# Plot the shapes.
#
range.shapes <- c(min(sapply(shapes, min)), max(sapply(shapes, max)))
palette(gray.colors(length(shapes)))
par(mfrow=c(1,2))
plot(range.shapes, range.shapes, type="n",asp=1, bty="n", xlab="", ylab="")
invisible(lapply(1:length(shapes), function(i) lines(t(shapes[[i]]), col=i, lwd=2)))
#
# Standardize the shape description.
#
standardize <- function(xy) {
  n <- dim(xy)[2]
  vectors <- xy[ ,-1, drop=FALSE] - xy[ ,-n, drop=FALSE]
  lengths <- sqrt(colSums(vectors^2))
  if (which.min(lengths - rev(lengths))*2 < n) {
    lengths <- rev(lengths)
    vectors <- vectors[, (n-1):1]
  }
  if (n > 2) {
    vectors <- vectors / rbind(lengths, lengths)
    perps <- rbind(-vectors[2, ], vectors[1, ])
    angles <- sapply(1:(n-2), function(i) {
      cosine <- sum(vectors[, i+1] * vectors[, i])
      sine <- sum(perps[, i+1] * vectors[, i])
      atan2(sine, cosine)
    })
    i <- min(which(angles != 0))
    angles <- sign(angles[i]) * angles
  } else angles <- numeric(0)
  list(lengths=lengths, angles=angles)
}
shapes.std <- lapply(shapes, standardize)
#
# Sort lexicographically.  (Not implemented: see the text.)
#
order.shape <- function(s) {
  order(sapply(s, function(s) s$lengths[1]))
}
i <- order.shape(shapes.std)
#
# Group.
#
equal.shape <- function(s.0, s.1) {
  same.length <- function(a,b) abs(a-b) <= (a+b) * 1e-8
  same.angle <- function(a,b) min(abs(a-b), abs(a-b)-2*pi) < 1e-11
  r <- function(u) {
    a <- u$angles
    if (length(a) > 0) {
      a <- rev(u$angles)
      i <- min(which(a != 0))
      a <- sign(a[i]) * a
    }
    list(lengths=rev(u$lengths), angles=a)
  }
  e <- function(u, v) {
    if (length(u$lengths) != length(v$lengths)) return (FALSE)
    all(mapply(same.length, u$lengths, v$lengths)) &&
      all(mapply(same.angle, u$angles, v$angles))
    }
  e(s.0, s.1) || e(r(s.0), s.1)
}
g <- rep(1, length(shapes.std))
for (j in 2:length(i)) {
  i.0 <- i[j-1]
  i.1 <- i[j]
  if (equal.shape(shapes.std[[i.0]], shapes.std[[i.1]])) 
    g[j] <- g[j-1] else g[j] <- g[j-1]+1
}
palette(rainbow(max(g)))
plot(range.shapes, range.shapes, type="n",asp=1, bty="n", xlab="", ylab="")
invisible(lapply(1:length(i), function(j) lines(t(shapes[[i[j]]]), col=g[j], lwd=2)))
Whuber
źródło
Gdy do grupy transformacji zalicza się dowolne rozszerzenia (lub „izotetyki”), klasy równoważności są klasami zgodności o geometrii afinicznej . Z tą komplikacją łatwo sobie poradzić: na przykład standaryzuj wszystkie polilinie, aby uzyskać całkowitą długość jednostki.
whuber
Dziękuję bardzo. Tylko jedno pytanie: czy kształty powinny być reprezentowane jako SpatialLines czy SpatialPolygons?
Klausos Klausos
Wieloboki tworzą kolejną komplikację: ich granice nie mają określonych punktów końcowych. Istnieje wiele sposobów radzenia sobie z tym, na przykład standaryzacja reprezentacji na początek (powiedzmy) wierzchołka, który sortuje najpierw w porządku leksykograficznym xy i postępuje w kierunku przeciwnym do ruchu wskazówek zegara wokół wielokąta. (Połączony wielokąt z topologicznie „czystym” będzie miał tylko jeden taki wierzchołek.) To, czy kształt zostanie uznany za wielokąt lub polilinię, zależy od rodzaju reprezentowanej przez niego cechy: nie ma żadnego wewnętrznego sposobu, aby powiedzieć o zamkniętej liście punktów, czy jest to przeznaczony do polilinii lub wielokąta.
whuber
Przepraszam za proste pytanie, ale powinienem je zadać, aby zrozumieć twój przykład. Twój obiekt shapes.std ma zarówno $ długości, jak i $ kąty. Jeśli jednak uruchomię ten kod na moich danych xy (np. [1,] 3093,5 -2987.8 [2,] 3072.7 -2991.0 itd.), Nie oszacuje kątów, ani nie narysuje kształtów. Jeśli uruchomię wykres (kształty [[1]]), wtedy wyraźnie zobaczę swoją polilinię. Jak więc zapisać polilinie w R, aby móc przetestować kod na moich danych?
Klausos Klausos
Zacząłem od tej samej struktury danych, którą zrobiłeś: tablicy współrzędnych (x, y). Moje tablice umieszczają te współrzędne w kolumnach (tak jakbyś używał rbind(x,y)zamiast cbind(x,y)). To wszystko, czego potrzebujesz: spbiblioteka nie jest używana. Jeśli chcesz śledzić, co dzieje się w szczegółach, proponujemy rozpocząć się z, powiedzmy n.shapes <- 2, n.shapes.new <- 3oraz p.mean <- 1. Następnie shapes, shapes.stditp są wystarczająco małe, aby być łatwo kontrolowane. Eleganckim - i „właściwym” sposobem rozwiązania tego wszystkiego byłoby stworzenie klasy znormalizowanych reprezentacji funkcji.
whuber
1

Prosisz o wiele z dowolną rotacją i dylatacją! Nie jestem pewien, jak przydatna byłaby odległość Hausdorffa, ale sprawdź to. Moje podejście polegałoby na zmniejszeniu liczby spraw do sprawdzenia za pomocą tanich danych. Na przykład można pominąć drogie porównania, jeśli długość dwóch linii nie jest liczbą całkowitą ( przy założeniu skalowania liczb całkowitych / stopniowanych ). Możesz również sprawdzić, czy obszar ramki granicznej lub ich wypukłe obszary kadłuba są w dobrym stosunku. Jestem pewien, że istnieje wiele tanich testów, które można wykonać względem środka ciężkości, takich jak odległości lub kąty od początku / końca.

Tylko wtedy, jeśli wykryjesz skalowanie, cofnij je i wykonaj naprawdę drogie kontrole.

Wyjaśnienie: Nie znam używanych pakietów. Przez stosunek liczb całkowitych rozumiałem, że powinieneś podzielić obie odległości, sprawdzić, czy wynikiem jest liczba całkowita, jeśli nie, odwrócić tę wartość (możliwe, że wybrano niewłaściwą kolejność) i ponownie sprawdzić. Jeśli dostaniesz liczbę całkowitą lub wystarczająco blisko, możesz wywnioskować, że być może miało miejsce skalowanie. Lub mogą to być dwa różne kształty.

Jeśli chodzi o obwiednię, prawdopodobnie masz przeciwne punkty prostokąta, który ją reprezentuje, więc wydobycie z nich obszaru jest prostą arytmetyką. Zasada porównania współczynników jest taka sama, tyle że wynik byłby podniesiony do kwadratu. Nie przejmuj się wypukłymi kadłubami, jeśli nie możesz ich dobrze wyciągnąć z tego pakietu R, to był tylko pomysł (prawdopodobnie i tak nie wystarczająco tani).

lynxlynxlynx
źródło
Wielkie dzięki. Czy możesz wyjaśnić, jak wykryć, czy długość dwóch linii nie jest liczbą całkowitą? Bardzo doceniam też, jeśli możesz podać przykład sprawdzenia, czy „obszar ramki granicznej lub obszary wypukłego kadłuba są w dobrym stosunku”
Klausos Klausos
Na przykład, jeśli wyodrębnię przestrzenną ramkę graniczną z danych przestrzennych, wtedy otrzymam tylko dwa punkty: spl <- sp :: SpatialLines (lista (Linie (Linia (xy.sp), ID = i))) b <- bbox ( spl)
Klausos Klausos
Rozszerzono główny post.
lynxlynxlynx
„Jeśli dostaniesz liczbę całkowitą lub wystarczająco blisko, możesz wywnioskować, że być może miało miejsce skalowanie”. Czy użytkownik nie mógł zastosować skali około 1,4?
Germán Carrillo
Jasne, ale moje założenie stało się jasne, zwłaszcza w przypadku późniejszych zmian. Wyobrażałem sobie powiększanie w stylu mapy internetowej, gdzie jedna jest ładnie ograniczona.
lynxlynxlynx
1

Dobrym sposobem na porównanie tych polilinii byłoby poleganie na reprezentacji jako sekwencji (odległości, kątów skrętu) na każdym wierzchołku: w przypadku linii złożonej z punktów P1, P2, ..., PNtaka sekwencja byłaby:

(odległość (P1P2), kąt (P1, P2, P3), odległość (P2P3), ..., kąt (P (N-2), P (N-1), PN), odległość (P (N-1) ) PN)).

Zgodnie z twoimi wymaganiami, dwie linie są równe wtedy i tylko wtedy, gdy odpowiadające im sekwencje są takie same (modulo kolejności i kierunku kąta). Porównywanie sekwencji numerów jest banalne.

Obliczając każdą sekwencję polilinii tylko raz i, jak sugeruje lynxlynxlynx, testując podobieństwo sekwencji tylko dla polilinii o takich samych trywialnych cechach (długość, liczba wierzchołków ...), obliczenia powinny być naprawdę szybkie!

Julien
źródło
To jest właściwy pomysł. Aby jednak zadziałało, należy rozwiązać wiele szczegółów, takich jak radzenie sobie z odbiciami, orientacja wewnętrzna, możliwość wielu połączonych komponentów i błąd zaokrąglenia zmiennoprzecinkowego. Są one omówione w rozwiązaniu, które podałem.
whuber
Tak, opisałem tylko główny pomysł. Twoja odpowiedź jest znacznie bardziej kompletna (jak często :-)
Julien