Algorithme pour trouver des points d'inflexion pour une polyligne

22

J'essaie de trouver les points d'inflexion, c'est-à-dire les points où les courbes d'une ligne commencent et se terminent. Si vous regardez l'image, la ligne verte peut être une route ou un ruisseau et les points noirs sont les points de départ et de fin des courbes. entrez la description de l'image ici

Quelles seraient les étapes de haut niveau pour automatiser la génération de ces points? J'ai un bureau ArcGIS et je suis assez maniable avec ArcObjects.

Devdatta Tengshe
la source
Les données source sont une polyligne composée de segments de ligne et vous souhaitez les approcher avec des courbes, ou contient-elle déjà des segments d'arc?
U2ros
Actuellement, il est composé de segments de ligne.
Devdatta Tengshe
1
L'illustration de cette question ressemble remarquablement à celle publiée sur esri.com/news/arcuser/0110/turning.html .
whuber
@whuber: Observation très astucieuse. C'était exactement la source de données que j'avais utilisée pour créer l'image.
Devdatta Tengshe

Réponses:

15

Lorsque la courbe est composée de segments de ligne, tous les points intérieurs de ces segments sont des points d'inflexion, ce qui n'est pas intéressant. Au lieu de cela, la courbe doit être considérée comme étant approximée par les sommets de ces segments. En splining une courbe deux fois différentiable par morceaux à travers ces segments, nous pouvons alors calculer la courbure. Un point d'inflexion proprement dit est alors un endroit où la courbure est nulle.

Dans l'exemple, il y a des étirements longs où la courbure est presque nulle. Cela suggère que les points indiqués devraient se rapprocher des extrémités de ces tronçons de régions à faible courbure.

Un algorithme efficace va donc spliner les sommets, calculer la courbure le long d'un ensemble dense de points intermédiaires, identifier les plages de courbure proche de zéro (en utilisant une estimation raisonnable de ce que signifie être "proche") et marquer les extrémités de ces plages .

Voici un Rcode de travail pour illustrer ces idées. Commençons par une chaîne de lignes exprimée comme une séquence de coordonnées:

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 les coordonnées x et y séparément pour obtenir un paramétrage de la courbe. (Le paramètre sera appelé time.)

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

Interpoler les splines pour le traçage et le calcul:

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

Nous avons besoin d'une fonction pour calculer la courbure d'une courbe paramétrée. Il doit estimer les première et deuxième dérivées de la spline. Avec de nombreuses splines (telles que des splines cubiques), il s'agit d'un calcul algébrique facile. Rfournit automatiquement les trois premiers dérivés. (Dans d'autres environnements, on peut vouloir calculer les dérivées numériquement.)

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

Je propose d' estimer un seuil de courbure nulle en termes d'étendue de la courbe. C'est au moins un bon point de départ; il doit être ajusté en fonction de la tortuosité de la courbe (c'est-à-dire augmentée pour les courbes plus longues). Il sera ensuite utilisé pour colorer les parcelles en fonction de la courbure.

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

Maintenant que les sommets ont été cannelés et la courbure calculée, il ne reste plus qu'à trouver les points d'inflexion . Pour les montrer, nous pouvons tracer les sommets, tracer la spline et y marquer les points d'inflexion.

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)

Terrain

Les points ouverts sont les sommets d'origine xyet les points noirs sont les points d'inflexion identifiés automatiquement avec cet algorithme. Étant donné que la courbure ne peut pas être calculée de manière fiable aux extrémités de la courbe, ces points ne sont pas spécialement marqués.

whuber
la source
Peut-être que la terminologie que j'ai utilisée était incorrecte. Ce que vous supposiez, c'est exactement ce que je voulais. Votre réponse semble prometteuse et je devrai travailler avec R pour traiter mon Shapefile.
Devdatta Tengshe
3

Vous pouvez utiliser l' outil Densifier . Dans ce cas, vous choisissez de densifier par angle, Ensuite, choisissez l'angle maximum accepté en ligne droite. Appliquez ensuite la ligne de résultat à l'outil Scinder la ligne aux sommets . Enfin, supprimez les lignes dont la longueur de forme est inférieure à la longueur de route minimale.

entrez la description de l'image ici

Dans cette image, nous voyons trois étapes:

1- Densifiez la ligne en utilisant l'angle. J'ai utilisé 10 degrés comme paramètre, et nous avons utilisé splitline. Dans l'image, la ligne courbe est dans sa phase initiale.

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

2- Sélectionnez les segments où la forme_longueur n'est pas redondante. Comme nous le voyons dans le tableau, je n'ai pas sélectionné ces longueurs redondantes. Ensuite, je les sélectionne dans une nouvelle classe d'entités.

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

3- Nous avons extrait les sommets situés dans les bords des lignes, qui sont des points d'inflexion.

arcpy.FeatureVerticesToPoints_management("line_split_selected" , "line_split_pnt" , "DANGLE")
geogeek
la source
J'ai les mêmes commentaires et questions concernant votre autre réponse: c'est une bonne idée, mais en même temps, il n'est pas clair qu'elle produira le résultat souhaité, ni comment choisir l'angle de seuil. Pourriez-vous fournir une illustration des résultats afin que les lecteurs puissent évaluer ce que fait réellement cette proposition? Fournir des exemples concrets est particulièrement important lors de la recommandation d'un logiciel ESRI dans le cadre d'une solution, car leurs algorithmes ne sont généralement pas documentés, ce qui rend impossible de savoir exactement ce qu'ils font.
whuber
pour être sûr que c'est une solution de travail, je dois la tester, mais je ne peux pas la tester, je manque les données, donc je suppose que les outils proposés par ESRI fonctionneront comme prévu, mais ces réponses doivent être testé davantage.
geogeek
nous pourrions leur
donner des
1
Voulez-vous que je les mette dans les commentaires, alors? BTW, si vous voulez tester les données, vous pouvez - pour commencer - utiliser les coordonnées que j'ai postées dans ma réponse, car elles sont proches de l'illustration de la question. Mais pourquoi ne pas simplement utiliser les données géographiques dont vous disposez?
whuber
2
oui vraiment cette solution fonctionne mieux pour extraire des lignes droites.
geogeek
1

Vous pouvez utiliser l' outil Généraliser qui a le décalage maximum par rapport à la ligne d'origine comme paramètre, afin que vous puissiez choisir le décalage qui correspond à votre cas.

entrez la description de l'image ici

Si nous nommons la ligne d'origine "line_cur" et la ligne généralisée "line_gen", nous pourrions couper "line_cur" par "line_gen". Le résultat sera le segment droit de "line_cur". Ensuite, nous pourrions nettoyer un segment très court en les supprimant avec une requête sql qui sélectionne la Shape_length supérieure à cette longueur de route minimale.

geogeek
la source
C'est une bonne idée. Cependant, il n'est pas clair à quel point cela fonctionnerait dans la pratique. Pourriez-vous peut-être montrer un exemple montrant les points d'inflexion trouvés?
whuber
j'ai fait un montage pour inclure une image, l'image explique comment cet outil peut faire coller une ligne avec des segments droits, nous devons donc faire un clip sur l'ancienne ligne, pour extraire uniquement les anciens segments de lignes droites
geogeek
quelque chose n'est pas clair je suis disponible pour répondre à vos questions?
geogeek
Je ne vois aucun point d'inflexion identifié dans l'illustration. Où seraient-ils exactement? Et comment choisir la tolérance de généralisation?
whuber
j'ai besoin de quelques données pour effectuer le test, mais je pense que nous devrions choisir la tolérance par expérimentation
geogeek