¿Suavizar polígonos en el mapa de contorno?

51

Aquí hay un mapa de contorno para el que están disponibles todos los polígonos de niveles.

Deje que pregunte cómo suavizar los polígonos manteniendo todos los vértices conservados en sus ubicaciones exactas.

De hecho, el contorno se hace encima de los datos de una cuadrícula, puede sugerir que se suavicen los datos de la cuadrícula y, por lo tanto, el contorno resultante será más suave. Tenga en cuenta que esto no funciona como deseo, ya que la función de suavizado, como el filtro gaussiano, eliminará pequeños paquetes de datos y cambiará el rango de la tercera variable, por ejemplo, la altura que no está permitida en mi aplicación.

En realidad estoy buscando un fragmento de código (preferiblemente en Python ) que pueda suavizar polígonos 2D (cualquier tipo: convexo, cóncavo, auto-intersección, etc.) razonablemente indoloro (olvídate páginas de códigos) y precisa.

Para tu información, hay una función en ArcGIS que hace esto perfectamente, pero el uso de aplicaciones comerciales de terceros no es mi elección para esta pregunta.

1)

Scipy.interpolate:

¡Como ves, las splines resultantes (rojo) no son satisfactorias!

2)

Aquí está el resultado usando el código dado en aquí . ¡No está funcionando bien!

3)

Paramí,lamejorsolucióndeberíaseralgocomolasiguientefiguraenlaqueuncuadradosesuavizagradualmentecambiandosolounvalor.Esperounconceptosimilarparasuavizarcualquierformadepolígonos.

Satisfacer la condición de que la spline pase los puntos:

4)

Aquíestámiimplementacióndela"idea de Whuber" línea por línea en Python en sus datos. Posiblemente hay algunos errores ya que los resultados no son buenos.

K=2esundesastre,porloqueparak>=4.

5)

EliminéunpuntoenlaubicaciónproblemáticaylasplineresultanteahoraesidénticaaladeWhuber.Perosiguesiendounapregunta:¿porquéelmétodonofuncionaentodosloscasos?

6)

Un buen suavizado para los datos de Whuber puede ser el siguiente (elaborado por el software de gráficos vectoriales) en el que se ha agregado un punto adicional sin problemas (comparar con la actualización

4):

7)

VeaelresultadodelaversióndePythondelcódigodeWhuberparaalgunasformasicónicas:


Note que el método parece no funcionar para polilíneas. Para la polilínea de la esquina (contorno) el verde es lo que quiero pero se volvió rojo. Esto debe abordarse ya que los mapas de contorno son siempre polilíneas, aunque las polilíneas cerradas se pueden tratar como polígonos como en mis ejemplos. Tampoco es que el problema surgido en la actualización 4 aún no se haya resuelto.

8) [mi última]

Aquí está la solución final (¡no es perfecta!):

Recuerda que tendrás que hacer algo respecto al área señalada por las estrellas. Quizás haya un error en mi código o el método propuesto necesite un mayor desarrollo para considerar todas las situaciones y proporcionar los resultados deseados.

    
pregunta Developer 31.01.2018 - 07:16

3 respuestas

35

La mayoría de los métodos para dividir secuencias de números dividirán polígonos. El truco es hacer que las splines se "cierren" sin problemas en los puntos finales. Para hacer esto, "envuelve" los vértices alrededor de los extremos. Luego divida las coordenadas x e y por separado.

Aquí hay un ejemplo de trabajo en R . Utiliza el procedimiento por defecto cubic spline disponible en el paquete de estadísticas básicas. Para mayor control, sustituya casi cualquier procedimiento que prefiera: simplemente asegúrese de que divide a través de los números (es decir, los interpola) en lugar de simplemente utilizarlos como "puntos de control".

#
# 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, ]
}

Para ilustrar su uso, vamos a crear un polígono pequeño (pero complicado).

#
# 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 usando el código anterior. Para hacer que la spline sea más suave, aumente el número de vértices de 100; para que sea menos suave, disminuya el número de vértices.

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

Para ver los resultados, trazamos (a) el polígono original en rojo discontinuo, mostrando el espacio entre el primer y el último vértice (es decir, sin cerrar su polilínea límite); y (b) la spline en gris, una vez más mostrando su espacio. (Debido a que el espacio es muy pequeño, sus puntos finales están resaltados con puntos azules).

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)

    
respondido por el whuber 07.05.2012 - 18:48
2

Sé que es una publicación antigua, pero apareció en Google para algo que estaba buscando, así que pensé en publicar mi solución.

No veo esto como un ejercicio de ajuste de curva 2D, sino más bien un ejercicio 3D. Al considerar los datos como 3D, podemos asegurarnos de que las curvas nunca se crucen entre sí, y podemos usar la información de otros contornos para mejorar nuestra estimación de la actual.

El siguiente extracto de iPython utiliza la interpolación cúbica proporcionada por SciPy. Tenga en cuenta que los valores z que he trazado no son importantes, siempre que todos los contornos sean equidistantes en altura.

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

Los resultados aquí no parecen los mejores, pero con tan pocos puntos de control siguen siendo perfectamente válidos. Observe cómo se extrae la línea verde ajustada para seguir el contorno azul más amplio.

    
respondido por el Gilly 08.02.2016 - 13:12
1

Escribí casi exactamente el paquete que está buscando ... pero estaba en Perl, y fue hace más de una década: GD :: Polyline . Utilizó curvas de Bézier cúbicas 2D y "suavizaría" un polígono arbitrario o "polilínea" (mi nombre entonces para lo que ahora se denomina comúnmente "LineString").

El algoritmo fue de dos pasos: dados los puntos en el Polígono, agregue dos puntos de control Bezier entre cada punto; luego llame a un algoritmo simple para hacer una aproximación por partes de la spline.

La segunda parte es fácil; La primera parte fue un poco de arte. Aquí estaba la información: considere un "segmento de control" un vértice N: vN . El segmento de control fue de tres puntos co-lineales: [cNa, vN, cNb] . El punto central era el vértice. La pendiente de este control seg fue igual a la pendiente de Vertex N-1 a Vertex N + 1. La longitud de la porción izquierda de este segmento fue 1/3 de la longitud desde Vertex N-1 a Vertex N, y la longitud de la porción derecha de este segmento fue 1/3 de la longitud desde Vertex N hasta Vertex N + 1. / p>

Si la curva original era de cuatro vértices: [v1, v2, v3, v4] , cada vértice ahora obtiene un segmento de control de la forma: [c2a, v2, c2b] . Encadene esto de esta manera: [v1, c1b, c2a, v2, c2b, c3a, v3, c3b, c4a, v4] y mastíquelos cuatro a la vez como los cuatro puntos Bezier: [v1, c1b, c2a, v2] , luego [v2, c2b, c3a, v3] , y así sucesivamente. Debido a que [c2a, v2, c2b] era co-lineal, la curva resultante será suave en cada vértice.

Por lo tanto, esto también cumple con su requisito de parametrizar la "rigidez" de la curva: use un valor más pequeño que 1/3 para una curva "más estrecha", uno más grande para un ajuste "en bucle". En cualquier caso, la curva resultante siempre pasa a través de los puntos originales dados.

Esto resultó en una curva suave que "circunscribió" el Polígono original. También tenía alguna forma de "inscribir" una curva suave ... pero no veo eso en el código CPAN.

De todos modos, en este momento no tengo una versión disponible en Python, ni tengo ninguna cifra. PERO ... si / cuando lo haga a Python, me aseguraré de publicar aquí.

    
respondido por el Dan H 04.01.2016 - 14:34

Lea otras preguntas en las etiquetas