Algoritmo para encontrar puntos de inflexión para una polilínea

Estoy tratando de averiguar los puntos de inflexión, es decir, los puntos donde comienzan y terminan las curvas en una línea. Si observa la imagen, la línea verde puede ser una carretera o una corriente, y los puntos negros son los puntos donde comienzan y terminan las curvas.

¿Cuáles serían los pasos de alto nivel para automatizar la generación de estos puntos? Tengo el escritorio ArcGIS y soy bastante útil con ArcObjects.

Cuando la curva se compone de segmentos de línea, todos los puntos interiores de esos segmentos son puntos de inflexión, lo que no es interesante. En cambio, la curva debe considerarse aproximada por los vértices de esos segmentos. Al dividir una curva por partes dos veces diferenciable a través de esos segmentos, podemos calcular la curvatura. Un punto de inflexión, estrictamente hablando, es un lugar donde la curvatura es cero.

En el ejemplo, hay tramos largos donde la curvatura es casi cero. Esto sugiere que los puntos indicados deberían aproximarse a los extremos de tales tramos de regiones de baja curvatura.

Por lo tanto, un algoritmo eficaz dividirá los vértices, calculará la curvatura a lo largo de un conjunto denso de puntos intermedios, identificará rangos de curvatura cercana a cero (utilizando una estimación razonable de lo que significa estar "cerca") y marcará los puntos finales de esos rangos .

Aquí hay un Rcódigo de trabajo para ilustrar estas ideas. Comencemos con una cadena de línea expresada como una secuencia de coordenadas:

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 las coordenadas x e y por separado para lograr una parametrización de la curva. (Se llamará al parámetro time).

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

Interpolar las splines para el trazado y el cálculo:

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

Necesitamos una función para calcular la curvatura de una curva parametrizada. Necesita estimar la primera y segunda derivada de la spline. Con muchas splines (como splines cúbicas), este es un cálculo algebraico fácil. Rproporciona los primeros tres derivados automáticamente. (En otros entornos, uno podría querer calcular las derivadas numéricamente).

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

Propongo estimar un umbral para la curvatura cero en términos de la extensión de la curva. Este al menos es un buen punto de partida; debe ajustarse de acuerdo con la tortuosidad de la curva (es decir, aumentada para curvas más largas). Esto luego se usará para colorear las parcelas de acuerdo con la curvatura.

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

Ahora que los vértices se han dividido y se ha calculado la curvatura, solo queda encontrar los puntos de inflexión . Para mostrarlos, podemos trazar los vértices, trazar la spline y marcar los puntos de inflexión en ella.

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)

Los puntos abiertos son los vértices originales xyy los puntos negros son los puntos de inflexión identificados automáticamente con este algoritmo. Debido a que la curvatura no puede calcularse de manera confiable en los puntos finales de la curva, esos puntos no están especialmente marcados.

Puede usar la herramienta Densificar . Para este caso, elige densificar por ángulo. Luego, elija el ángulo máximo aceptado en una línea recta. Luego aplique a la línea de resultado a la herramienta Dividir línea en vértices . Finalmente, elimine las líneas que tienen shape_length más pequeñas que la longitud mínima del camino.

En esta imagen, vemos tres pasos:

1- Densificar la línea usando el ángulo. He usado 10 grados como parámetro, y usamos splitline. En la imagen, la línea curva está en su fase inicial.

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

2- Seleccione los segmentos donde tener shape_length no sea redundante. Como vemos en la tabla, no he seleccionado esas longitudes redundantes. Luego, los selecciono en una nueva clase de entidad.

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

3- Hemos extraído los vértices ubicados en los bordes de las líneas, que son puntos de inflexión.

arcpy.FeatureVerticesToPoints_management("line_split_selected" , "line_split_pnt" , "DANGLE")

Puede utilizar la herramienta Generalizar que tiene el desplazamiento máximo de la línea original como parámetro, por lo que puede elegir el desplazamiento que se ajuste a su caso.

Si nombramos la línea original "line_cur" y la generalizada "line_gen", podríamos recortar "line_cur" por "line_gen". El resultado será el segmento recto de "line_cur". Luego podríamos limpiar algunos segmentos muy cortos eliminándolos con una consulta sql que selecciona Shape_length mayor que la longitud mínima del camino.