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)
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.
