Algorytm znajdowania punktów przegięcia dla polilinii

22

Próbuję znaleźć punkty przegięcia, tj. Punkty, w których krzywe linii zaczynają się i kończą. Jeśli spojrzysz na obraz, zielona linia może być drogą lub strumieniem, a czarne punkty to punkty, w których krzywe zaczynają się i kończą. wprowadź opis zdjęcia tutaj

Jakie byłyby kroki na wysokim szczeblu w celu automatyzacji generowania tych punktów? Mam pulpit ArcGIS i jestem całkiem przydatny w ArcObjects.

Devdatta Tengshe
źródło
Dane źródłowe to polilinia złożona z segmentów linii i chcesz ją przybliżyć za pomocą krzywych, czy też ma już segmenty łukowe?
U2ros,
Obecnie składa się z segmentów linii.
Devdatta Tengshe,
1
Ilustracja w tym pytaniu wygląda zadziwiająco jak opublikowana na stronie esri.com/news/arcuser/0110/turning.html .
whuber
@whuber: Bardzo przebiegła obserwacja. To było dokładnie źródło danych, którego użyłem do stworzenia obrazu.
Devdatta Tengshe,

Odpowiedzi:

15

Gdy krzywa składa się z odcinków linii, wówczas wszystkie punkty wewnętrzne tych odcinków są punktami przegięcia, co nie jest interesujące. Zamiast tego należy uważać, że krzywa jest zbliżona do wierzchołków tych segmentów. Poprzez rozcięcie częściowo segmentowanej krzywej o dwóch różnicach między tymi segmentami, możemy następnie obliczyć krzywiznę. Mówiąc wprost, punkt przegięcia to miejsce, w którym krzywizna wynosi zero.

W tym przykładzie występują długie odcinki, w których krzywizna wynosi prawie zero. Sugeruje to, że wskazane punkty powinny przybliżać końce takich odcinków regionów o niskiej krzywiźnie.

Dlatego skuteczny algorytm splajtuje wierzchołki, oblicza krzywiznę wzdłuż gęstego zestawu punktów pośrednich, identyfikuje zakresy prawie zerowej krzywizny (wykorzystując rozsądne oszacowanie, co to znaczy być „blisko”), i zaznacza punkty końcowe tych zakresów .

Oto działający Rkod ilustrujący te pomysły. Zacznijmy od ciągu linii wyrażonego jako ciąg współrzędnych:

xy <- matrix(c(5,20, 3,18, 2,19, 1.5,16, 5.5,9, 4.5,8, 3.5,12, 2.5,11, 3.5,3, 
               2,3, 2,6, 0,6, 2.5,-4, 4,-5, 6.5,-2, 7.5,-2.5, 7.7,-3.5, 6.5,-8), ncol=2, byrow=TRUE)

Spline x i y współrzędne osobno do uzyskania parametryzacji krzywej. (Parametr zostanie wywołany time.)

n <- dim(xy)[1]
fx <- splinefun(1:n, xy[,1], method="natural")
fy <- splinefun(1:n, xy[,2], method="natural")

Interpoluj splajny na potrzeby kreślenia i obliczeń:

time <- seq(1,n,length.out=511)
uv <- sapply(time, function(t) c(fx(t), fy(t)))

Potrzebujemy funkcji do obliczenia krzywizny sparametryzowanej krzywej. Musi oszacować pierwszą i drugą pochodną splajnu. W przypadku wielu splajnów (takich jak splajny sześcienne) jest to łatwe obliczenie algebraiczne. Rzapewnia trzy pierwsze pochodne automatycznie. (W innych środowiskach można chcieć obliczać pochodne numerycznie).

curvature <- function(t, fx, fy) {
  # t is an argument to spline functions fx and fy.
  xp <- fx(t,1); yp <- fy(t,1)            # First derivatives
  xpp <- fx(t,2); ypp <- fy(t,2)          # Second derivatives
  v <- sqrt(xp^2 + yp^2)                  # Speed
  (xp*ypp - yp*xpp) / v^3                 # (Signed) curvature
  # (Left turns have positive curvature; right turns, negative.)
}

kappa <- abs(curvature(time, fx, fy))     # Absolute curvature of the data

Proponuję oszacować próg zerowej krzywizny pod względem zasięgu krzywej. To przynajmniej dobry punkt wyjścia; należy go wyregulować zgodnie z krętością krzywej (to znaczy zwiększyć dla dłuższych krzywych). Będzie to później wykorzystane do pokolorowania wykresów zgodnie z krzywizną.

curvature.zero <- 2*pi / max(range(xy[,1]), range(xy[,2])) # A small threshold
i.col <- 1 + floor(127 * curvature.zero/(curvature.zero + kappa)) 
palette(terrain.colors(max(i.col)))                        # Colors

Teraz, gdy wierzchołki zostały spline i obliczono krzywiznę, pozostaje tylko znaleźć punkty przegięcia . Aby je pokazać, możemy narysować wierzchołki, narysować splajn i zaznaczyć na nim punkty przegięcia.

plot(xy, asp=1, xlab="x",ylab="y", type="n")
tmp <- sapply(2:length(kappa), function(i) lines(rbind(uv[,i-1],uv[,i]), lwd=2, col=i.col[i]))
points(t(sapply(time[diff(kappa < curvature.zero/2) != 0], 
       function(t) c(fx(t), fy(t)))), pch=19, col="Black")
points(xy)

Wątek

Punkty otwarte to oryginalne wierzchołki, xya czarne punkty to punkty przegięcia automatycznie identyfikowane za pomocą tego algorytmu. Ponieważ krzywizny nie można wiarygodnie obliczyć w punktach końcowych krzywej, punkty te nie są specjalnie oznaczone.

Whuber
źródło
Może użyta przeze mnie terminologia była błędna. To, co zakładałeś, jest dokładnie tym, czego chciałem. Twoja odpowiedź wygląda obiecująco i będę musiał współpracować z R, aby przetworzyć mój plik kształtu.
Devdatta Tengshe
3

Możesz użyć narzędzia Densify . W tym przypadku wybierz zagęszczenie według kąta, następnie wybierz maksymalny kąt akceptowany w linii prostej. Następnie zastosuj do linii wyników do narzędzia Linia podziału w wierzchołkach . Na koniec usuń linie o kształcie_długości mniejszym niż minimalna długość drogi.

wprowadź opis zdjęcia tutaj

Na tym zdjęciu widzimy trzy kroki:

1- Zagęścić linię za pomocą kąta. Jako parametru użyłem 10 stopni i zastosowaliśmy linię podziału. Na zdjęciu zakrzywiona linia znajduje się w początkowej fazie.

arcpy.Densify_edit("line" , "ANGLE" , "","",10)
arcpy.SplitLine_management("line" , "line_split")

2 - Wybierz segmenty, w których długość kształtu nie jest redundantna. Jak widać w tabeli, nie wybrałem tych zbędnych długości. Następnie wybieram je do nowej klasy funkcji.

arcpy.Select_analysis("line_split" , "line_split_selected")

3- Wyodrębniliśmy wierzchołki znajdujące się na krawędziach linii, które są punktami przegięcia.

arcpy.FeatureVerticesToPoints_management("line_split_selected" , "line_split_pnt" , "DANGLE")
geogeek
źródło
Mam te same komentarze i pytania dotyczące twojej drugiej odpowiedzi: to fajny pomysł, ale jednocześnie nie jest jasne, czy przyniesie pożądany rezultat, ani jak wybrać kąt progowy. Czy możesz podać ilustrację wyników, aby czytelnicy mogli ocenić, co naprawdę robi ta propozycja? Podanie sprawdzonych przykładów jest szczególnie ważne, gdy zaleca się oprogramowanie ESRI jako część rozwiązania, ponieważ ich algorytmy zwykle nie są dokumentowane, co uniemożliwia dokładne określenie, co robią.
whuber
aby mieć pewność, że jest to działające rozwiązanie, muszę je przetestować, ale nie mogę go przetestować, brakuje mi danych, więc przypuszczam, że zaproponowane przez ESRI narzędzia będą działać zgodnie z oczekiwaniami, ale te odpowiedzi muszą być dalej testowanym.
geogeek
moglibyśmy nazwać je pomysłami, a nie odpowiedziami
geogeek
1
Czy chciałbyś, żebym je zamieścił w komentarzach? BTW, jeśli chcesz przetestować dane, możesz - na początek - użyć współrzędnych podanych w mojej odpowiedzi, ponieważ są one zbliżone do ilustracji w pytaniu. Ale dlaczego po prostu nie skorzystać z posiadanych danych geograficznych?
whuber
2
tak, naprawdę to rozwiązanie działa lepiej przy wydobywaniu prostych linii.
geogeek
1

Możesz użyć narzędzia Uogólnij , który ma maksymalne przesunięcie od oryginalnej linii jako parametr, dzięki czemu możesz wybrać przesunięcie, które pasuje do twojego przypadku.

wprowadź opis zdjęcia tutaj

Jeśli nazwiemy pierwotną linię „line_cur”, a uogólnioną „line_gen”, możemy przyciąć „line_cur” przez „line_gen”. Wynikiem będzie prosty segment „line_cur”. Następnie możemy wyczyścić bardzo krótki odcinek, usuwając go za pomocą zapytania SQL, które wybiera Shape_length większą niż minimalna długość drogi.

geogeek
źródło
To niezły pomysł. Nie jest jednak jasne, jak dobrze sprawdziłby się w praktyce. Czy mógłbyś pokazać przykład pokazujący znalezione punkty przegięcia?
whuber
dokonałem edycji, aby dołączyć zdjęcie, obraz wyjaśnia, jak to narzędzie może sprawić, że linia przylega do prostych odcinków, więc musimy zrobić klip do starej linii, aby wyodrębnić tylko stare odcinki prostych linii
geogeek
czy coś jest niejasne, czy mogę odpowiedzieć na twoje pytania?
geogeek
Nie widzę żadnych punktów przegięcia na ilustracji. Gdzie oni by dokładnie byli? A jak wybrać tolerancję dla uogólnienia?
whuber
potrzebuję trochę danych, aby wykonać test, ale myślę, że powinniśmy wybrać tolerancję poprzez eksperymenty
geogeek