Algoritmo para descobrir pontos de inflexão para uma polilinha

Estou tentando descobrir os pontos de inflexão, ou seja, os pontos em que as curvas de uma linha começam e terminam. Se você olhar para a imagem, a linha verde pode ser uma estrada ou um riacho, e os pontos pretos são os pontos onde as curvas começam e terminam.

Quais seriam as etapas de alto nível para automatizar a geração desses pontos? Eu tenho o ArcGIS desktop e sou bastante útil com o ArcObjects.

Quando a curva é composta por segmentos de linha, todos os pontos internos desses segmentos são pontos de inflexão, o que não é interessante. Em vez disso, a curva deve ser pensada como sendo aproximada pelos vértices desses segmentos. Ao dividir uma curva que pode ser diferenciada por duas vezes nesses segmentos, podemos calcular a curvatura. Um ponto de inflexão, a rigor, é então um lugar onde a curvatura é zero.

No exemplo, existem trechos longos, onde a curvatura é quase zero. Isso sugere que os pontos indicados devem aproximar as extremidades de tais trechos de regiões de baixa curvatura.

Um algoritmo eficaz, portanto, divide os vértices, calcula a curvatura ao longo de um conjunto denso de pontos intermediários, identifica faixas de curvatura quase zero (usando alguma estimativa razoável do que significa estar "próximo") e marca os pontos finais dessas faixas .

Aqui está o Rcódigo de trabalho para ilustrar essas idéias. Vamos começar com uma sequência de linhas expressa como uma sequência 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)

Divine as coordenadas x e y separadamente para obter uma parametrização da curva. (O parâmetro será chamado time.)

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

Interpole os splines para plotagem e computação:

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

Precisamos de uma função para calcular a curvatura de uma curva parametrizada. Ele precisa estimar a primeira e a segunda derivada do spline. Com muitos splines (como splines cúbicos), esse é um cálculo algébrico fácil. Rfornece os três primeiros derivados automaticamente. (Em outros ambientes, pode-se querer calcular os derivados numericamente.)

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

Proponho estimar um limiar para a curvatura zero em termos da extensão da curva. Este pelo menos é um bom ponto de partida; deve ser ajustado de acordo com a tortuosidade da curva (ou seja, aumentada para curvas mais longas). Mais tarde, isso será usado para colorir as plotagens de acordo com a 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

Agora que os vértices foram estriados e a curvatura calculada, resta apenas encontrar os pontos de inflexão . Para mostrá-los, podemos plotar os vértices, plotar o spline e marcar os pontos de inflexão nele.

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)

Os pontos em aberto são os vértices originais xye os pontos pretos são os pontos de inflexão identificados automaticamente com esse algoritmo. Como a curvatura não pode ser calculada com segurança nos pontos finais da curva, esses pontos não são especialmente marcados.

Você pode usar a ferramenta Densify . Nesse caso, você escolhe densificar por ângulo. Em seguida, escolha o ângulo máximo aceito em uma linha reta. Em seguida, aplique à linha de resultado na ferramenta Dividir linha nos vértices . Por fim, exclua as linhas com comprimento_forma menor ao comprimento mínimo da estrada.

Nesta foto, vemos três etapas:

1- Densifique a linha usando o ângulo. Eu usei 10 graus como parâmetro e usamos splitline. Na figura, a linha curva está em sua fase inicial.

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

2- Selecione os segmentos nos quais o comprimento da forma não é redundante. Como podemos ver na tabela, não selecionei esses comprimentos redundantes. Depois, seleciono-os em uma nova classe de recurso.

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

3- Extraímos os vértices localizados nas bordas das linhas, que são pontos de inflexão.

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

Você pode usar a ferramenta Generalizar , que possui o deslocamento máximo da linha original como parâmetro, para poder escolher o deslocamento adequado ao seu caso.

Se nomearmos a linha original "line_cur" e a generalizada "line_gen", poderíamos recortar "line_cur" por "line_gen". O resultado será o segmento direto de "line_cur". Em seguida, poderíamos limpar um segmento muito curto, excluindo-o com uma consulta sql que seleciona o Comprimento da forma maior que o comprimento mínimo da estrada.