Lissage des polygones dans la carte de contour?

52

Voici une carte de contour pour laquelle tous les polygones de niveaux sont disponibles.

Laissez-vous demander comment lisser les polygones en conservant tous les sommets à leur emplacement exact?

En effet, le contour étant créé au-dessus d’une donnée de grille, vous pouvez alors suggérer de lisser les données de la grille afin que le contour résultant soit plus lisse. Notez que cela ne fonctionne pas comme je le souhaite car la fonction de lissage, telle que le filtre gaussien, supprimera les petits paquets de données et modifiera la plage de la troisième variable, par exemple la hauteur non autorisée dans mon application.

En fait, je cherche un morceau de code (de préférence en Python ) capable de lisser des polygones 2D (tout type: convexe, concave, s'entrecroisant, etc.) raisonnablement indolore (oublier les pages de codes) et précis.

Pour votre information, il existe une fonction dans ArcGIS qui le fait parfaitement, mais l’utilisation d’applications commerciales tierces n’est pas mon choix pour cette question.

entrez la description de l'image ici


1)

Scipy.interpolate:

entrez la description de l'image ici

Comme vous le voyez, les splines obtenues (en rouge) ne sont pas satisfaisantes!

2)

Voici le résultat en utilisant le code donné ici . Cela ne fonctionne pas bien!

entrez la description de l'image ici

3)

Pour moi, la meilleure solution devrait être quelque chose comme la figure suivante dans laquelle un carré est lissé progressivement en ne modifiant qu'une seule valeur. J'espère un concept similaire pour lisser toute forme de polygone.

entrez la description de l'image ici

Satisfaire à la condition que spline passe les points:

entrez la description de l'image ici

4)

Voici mon implémentation de "l'idée de whuber" ligne par ligne en Python sur ses données. Il y a peut-être des bugs car les résultats ne sont pas bons.

entrez la description de l'image ici

K = 2 est un désastre et donc pour k> = 4.

5)

J'ai supprimé un point de l'emplacement problématique et la spline résultante est maintenant identique à celle de whuber. Mais il reste une question que pourquoi la méthode ne fonctionne pas pour tous les cas?

entrez la description de l'image ici

6)

Un bon lissage des données de whuber peut être le suivant (dessiné par un logiciel de graphisme vectoriel) dans lequel un point supplémentaire a été ajouté sans à-coups (comparer avec update

4):

entrez la description de l'image ici

sept)

Voir le résultat de la version Python du code de whuber pour certaines formes emblématiques:

entrez la description de l'image ici
Notez que la méthode semble ne pas fonctionner pour les polylignes. Pour la polyligne d’angle (contour), le vert est ce que je veux, mais j’ai eu le rouge. Ceci doit être traité car les cartes de contour sont toujours des polylignes, bien que les polylignes fermées puissent être traitées comme des polygones, comme dans mes exemples. En outre, le problème apparu dans la mise à jour 4 n’a pas encore été résolu.

8) [mon dernier]

Voici la solution finale (pas parfaite!):

entrez la description de l'image ici

N'oubliez pas que vous devrez faire quelque chose à propos de la zone indiquée par les étoiles. Il y a peut-être un bogue dans mon code ou la méthode proposée nécessite un développement supplémentaire pour prendre en compte toutes les situations et fournir les résultats souhaités.

développeur
la source
Comment générez-vous des contours «polygonaux»? ne seraient-ils pas toujours des lignes, puisqu'un contour coupant le bord d'un DEM ne se fermerait jamais sur lui-même?
pistachionut
J'ai utilisé la fonction v.generalize dans GRASS pour lisser les courbes de niveau avec des résultats corrects, bien que cela puisse prendre un certain temps pour les cartes avec des contours très denses.
pistachionut
@pistachionut Vous pouvez considérer que les niveaux de contour sont des polylignes. Je cherche du code pur au premier stade. Si non disponible, allumer le paquet pour Python.
Développeur
Peut-être regardez-vous sur scipy.org/Cookbook/Interpolation parce que vous avez l’impression de vouloir spline
PolyGeo
1
La courbe de @Pablo Bézier dans votre lien fonctionne bien pour les polylignes. Whuber fonctionne presque bien pour les polygones. Ainsi, ensemble, ils pourraient répondre à la question. Merci beaucoup de partager vos connaissances gratuitement.
Développeur

Réponses:

37

La plupart des méthodes pour spline des suites de nombres spline des polygones. L'astuce consiste à faire en sorte que les splines se "referment" en douceur aux extrémités. Pour ce faire, "envelopper" les sommets autour des extrémités. Ensuite, séparez les coordonnées x et y séparément.

Voici un exemple de travail dans R. Il utilise la splineprocédure cubique par défaut disponible dans le package de statistiques de base. Pour plus de contrôle, remplacer presque toute procédure que vous préférez: assurez - vous qu'il splines par les chiffres (qui est, interpole eux) plutôt que de simplement les utiliser comme des « points de contrôle ».

#
# Splining a polygon.
#
#   The rows of 'xy' give coordinates of the boundary vertices, in order.
#   'vertices' is the number of spline vertices to create.
#              (Not all are used: some are clipped from the ends.)
#   'k' is the number of points to wrap around the ends to obtain
#       a smooth periodic spline.
#
#   Returns an array of points. 
# 
spline.poly <- function(xy, vertices, k=3, ...) {
    # Assert: xy is an n by 2 matrix with n >= k.

    # Wrap k vertices around each end.
    n <- dim(xy)[1]
    if (k >= 1) {
        data <- rbind(xy[(n-k+1):n,], xy, xy[1:k, ])
    } else {
        data <- xy
    }

    # Spline the x and y coordinates.
    data.spline <- spline(1:(n+2*k), data[,1], n=vertices, ...)
    x <- data.spline$x
    x1 <- data.spline$y
    x2 <- spline(1:(n+2*k), data[,2], n=vertices, ...)$y

    # Retain only the middle part.
    cbind(x1, x2)[k < x & x <= n+k, ]
}

Pour illustrer son utilisation, créons un petit polygone (mais compliqué).

#
# Example polygon, randomly generated.
#
set.seed(17)
n.vertices <- 10
theta <- (runif(n.vertices) + 1:n.vertices - 1) * 2 * pi / n.vertices
r <- rgamma(n.vertices, shape=3)
xy <- cbind(cos(theta) * r, sin(theta) * r)

Spline le en utilisant le code précédent. Pour rendre la spline plus lisse, augmentez le nombre de sommets de 100; pour le rendre moins lisse, diminuez le nombre de sommets.

s <- spline.poly(xy, 100, k=3)

Pour voir les résultats, on trace (a) le polygone d'origine en pointillé rouge, montrant l'écart entre le premier et le dernier sommet (c'est-à-dire ne ferme pas sa polyligne limite); et (b) la spline en gris, montrant une fois de plus son écart. (L'écart étant si petit, ses extrémités sont mises en évidence par des points bleus.)

plot(s, type="l", lwd=2, col="Gray")
lines(xy, col="Red", lty=2, lwd=2)
points(xy, col="Red", pch=19)
points(s, cex=0.8)
points(s[c(1,dim(s)[1]),], col="Blue", pch=19)

Polygone spliné

whuber
la source
5
Bonne réponse. Existe-t-il un moyen de garantir que les contours ne se croisent pas à la suite du lissage?
Kirk Kuykendall
C'est une bonne question, @Kirk. Je ne suis au courant d'aucune méthode pour garantir le non-franchissement de cette forme de lissage. (En fait, je ne vois même pas comment garantir que la polyligne lissée ne s'entrecroise pas. Ce n'est toutefois pas un gros problème pour la plupart des contours.) Pour ce faire, vous devez revenir à l'original. DEM et utilise plutôt une meilleure méthode pour calculer les contours en premier lieu. (Il existe de meilleures méthodes - elles sont connues depuis longtemps - mais autant que je sache, certains des SIG les plus populaires ne les utilisent pas.)
whuber
Premièrement, je travaille toujours à implémenter votre réponse en Python, mais sans succès. Deuxièmement, quel sera le résultat si vous appliquez votre méthode sur un carré? Vous pouvez vous référer à ceux que j'ai dessinés dans la question.
Développeur
1
J'ai accepté ceci comme réponse car cela donne une bonne solution. Même si ce n’est pas une solution parfaite, mais que cela m’a donné quelques idées, espérons que je trouverai une solution qui satisfasse aux points que j’ai mentionnés ci-dessus dans ma question et mes commentaires. Vous pouvez également considérer les commentaires de whuber pour la question [QC], il y a de bons trucs là. Enfin, je devrais dire que la traduction en python est presque simple après l'installation du joli paquet Scipy. Considérez également le commentaire de Pablo dans QC comme solution possible pour les polylignes, c'est-à-dire les courbes de Bézier. Bonne chance à tous.
Développeur
1
En voyant vos réponses, je regrette de ne pas avoir bien pris soin de mes calculs !!!
vinayan
2

Je sais que ceci est un ancien message, mais il est apparu sur Google pour quelque chose que je cherchais et j'ai donc pensé poster ma solution.

Je ne vois pas cela comme un exercice d’ajustement de courbe en 2D, mais plutôt en 3D. En considérant les données en 3D, nous pouvons nous assurer que les courbes ne se croisent jamais et que nous pouvons utiliser les informations provenant d'autres contours pour améliorer notre estimation pour le graphique actuel.

L'extrait iPython suivant utilise une interpolation cubique fournie par SciPy. Notez que les valeurs z que j'ai tracées ne sont pas importantes, tant que tous les contours sont équidistants en hauteur.

In [1]: %pylab inline
        pylab.rcParams['figure.figsize'] = (10, 10)
        Populating the interactive namespace from numpy and matplotlib

In [2]: import scipy.interpolate as si

        xs = np.array([0.0, 0.0, 4.5, 4.5,
                       0.3, 1.5, 2.3, 3.8, 3.7, 2.3,
                       1.5, 2.2, 2.8, 2.2,
                       2.1, 2.2, 2.3])
        ys = np.array([0.0, 3.0, 3.0, 0.0,
                       1.1, 2.3, 2.5, 2.3, 1.1, 0.5,
                       1.1, 2.1, 1.1, 0.8,
                       1.1, 1.3, 1.1])
        zs = np.array([0,   0,   0,   0,
                       1,   1,   1,   1,   1,   1,
                       2,   2,   2,   2,
                       3,   3,   3])
        pts = np.array([xs, ys]).transpose()

        # set up a grid for us to resample onto
        nx, ny = (100, 100)
        xrange = np.linspace(np.min(xs[zs!=0])-0.1, np.max(xs[zs!=0])+0.1, nx)
        yrange = np.linspace(np.min(ys[zs!=0])-0.1, np.max(ys[zs!=0])+0.1, ny)
        xv, yv = np.meshgrid(xrange, yrange)
        ptv = np.array([xv, yv]).transpose()

        # interpolate over the grid
        out = si.griddata(pts, zs, ptv, method='cubic').transpose()

        def close(vals):
            return np.concatenate((vals, [vals[0]]))

        # plot the results
        levels = [1, 2, 3]
        plt.plot(close(xs[zs==1]), close(ys[zs==1]))
        plt.plot(close(xs[zs==2]), close(ys[zs==2]))
        plt.plot(close(xs[zs==3]), close(ys[zs==3]))
        plt.contour(xrange, yrange, out, levels)
        plt.show()

Résultat Interpolé Cubique

Les résultats obtenus ne sont pas optimaux, mais avec si peu de points de contrôle, ils restent parfaitement valables. Notez comment la ligne ajustée verte est tirée pour suivre le contour bleu plus large.

Gilly
la source
Les courbes lisses ajustées doivent rester aussi proches que possible du polygone / polyligne original.
Développeur
1

J'ai écrit presque exactement le package que vous recherchez ... mais c'était en Perl et l'était il y a plus de dix ans: GD :: Polyline . Il utilisait des courbes de Bézier cubiques 2D et "lisserait" un polygone ou une "polyligne" arbitraire (mon nom s’appelait alors à ce que l’on appelle communément un "LineString").

L'algorithme comportait deux étapes: compte tenu des points du polygone, ajoutez deux points de contrôle de Bézier entre chaque point; appelez ensuite un algorithme simple pour faire une approximation par morceaux de la spline.

La deuxième partie est facile. la première partie était un peu d'art. Voici l'aperçu était: envisager un « segment de contrôle » un Vertex N: vN. Le segment de contrôle était de trois points de co-linéaire: [cNa, vN, cNb]. Le point central était le sommet. La pente de ce segment de contrôle était égale à la pente du sommet N-1 au sommet N + 1. La longueur de la partie gauche de ce segment était égale à 1/3 de la longueur allant du sommet N-1 au sommet N, et la longueur de la partie droite de ce segment était égale à 1/3 de la longueur du sommet N au sommet N + 1.

Si la courbe d' origine était de quatre sommets: [v1, v2, v3, v4]alors chaque sommet maintenant obtenir un segment de contrôle de la forme: [c2a, v2, c2b]. Enchaînez-les comme ceci: [v1, c1b, c2a, v2, c2b, c3a, v3, c3b, c4a, v4]et broyez-les quatre à la fois comme les quatre points de Bézier:, [v1, c1b, c2a, v2]puis [v2, c2b, c3a, v3], et ainsi de suite. Parce que [c2a, v2, c2b]co-linéaires, la courbe résultante sera lisse à chaque sommet.

Donc, cela répond également à votre exigence de paramétrer le "serrage" de la courbe: utilisez une valeur inférieure à 1/3 pour une courbe "plus étroite", une valeur plus grande pour un ajustement "plus bouclé". Dans les deux cas, la courbe résultante passe toujours par les points donnés d'origine.

Cela a abouti à une courbe lisse qui "circonscrit" le polygone d'origine. J'avais aussi un moyen "d'inscrire" une courbe lisse ... mais je ne le vois pas dans le code CPAN.

Quoi qu'il en soit, je n'ai pas encore de version disponible en Python ni de chiffres. MAIS ... si / quand je porte ceci en Python, je serai sûr de poster ici.

Dan H
la source
Impossible d'évaluer le code Perl, ajoutez des graphiques pour démontrer son fonctionnement, si possible.
Développeur