¿Suavizar polígonos en el mapa de contorno?

Question

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

Preguntemos ¿Cómo alisar los polígonos manteniendo todos los vértices en su ubicación exacta?

De hecho, el contorno se realiza sobre una cuadrícula, por lo que se 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 mi 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án permitidos en mi aplicación.

En realidad estoy buscando un trozo de código (preferiblemente en Python ) que puede hacer el suavizado de polígonos 2D (de cualquier tipo: convexo, cóncavo, auto-intersección, etc) razonablemente indoloro (olvídate de las páginas de códigos) y preciso.

Para su información, hay una función en ArcGIS que lo hace perfectamente, pero usar aplicaciones comerciales de terceros no es mi opción para esta cuestión.

---

1)

Scipy.interpolate:

Como ves, las splines resultantes (en rojo) no son satisfactorias.

2)

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

3)

Para mí la mejor solución debería ser algo como la siguiente figura en la que un cuadrado se está suavizando gradualmente cambiando sólo un valor. Espero un concepto similar para suavizar cualquier forma de polígonos.

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

4)

Aquí está mi implementación de la "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=2 es un desastre y lo mismo para k>=4.

5)

He eliminado un punto en el lugar problemático y la spline resultante es ahora idéntica a la de whuber. Pero sigue siendo una pregunta que por qué el método no funciona para todos los casos?

6)

Un buen alisado para los datos de Whuber puede ser el siguiente (dibujado con un software de gráficos vectoriales) en el que se ha añadido suavemente un punto extra (comparar con la actualización

4):

7)

Vea el resultado de la versión Python del código de Whuber para algunas formas icónicas:

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

8) [mi última]

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

Recuerda que tendrás que hacer algo en la zona señalada por las estrellas. Quizá 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.

Answer 1

La mayoría de los métodos para splinear secuencias de números splinearán polígonos. El truco consiste en hacer que los splines se "cierren" suavemente en los puntos finales. Para ello, "envuelva" los vértices alrededor de los extremos. Luego, spline las coordenadas x e y por separado.

Este es un ejemplo de trabajo en R . Utiliza el cubo por defecto spline procedimiento disponible en el paquete estadístico básico. Para un mayor control, sustituya casi cualquier procedimiento que prefiera: sólo asegúrese de que splines a través de los números (es decir, los interpola) en lugar de utilizarlos simplemente 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 pequeño (pero complicado) polígono.

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

Divídelo utilizando el código anterior. Para que la spline sea más suave, aumenta el número de vértices a partir de 100; para que sea menos suave, disminuye 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 hueco entre el primer y el último vértice (es decir, sin cerrar su polilínea límite); y (b) la spline en gris, mostrando una vez más su hueco. (Como el hueco es tan pequeño, sus puntos extremos se resaltan 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)

Answer 2

Sé que este es un post antiguo, pero apareció en Google para algo que estaba buscando, así que pensé en publicar mi solución.

No lo veo como un ejercicio de ajuste de curvas en 2D, sino en 3D. Al considerar los datos como 3D podemos asegurarnos de que las curvas nunca se cruzan entre sí, y podemos utilizar la información de otros contornos para mejorar nuestra estimación para el actual.

El siguiente extracto de iPython utiliza la interpolación cúbica proporcionada por SciPy. Tenga en cuenta que los valores de z que he trazado no son importantes, siempre y cuando todos los contornos son 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 son los mejores, pero con tan pocos puntos de control siguen siendo perfectamente válidos. Obsérvese cómo la línea verde ajustada se desplaza para seguir el contorno azul más amplio.

Answer 3

Escribí casi exactamente el paquete que estás buscando... pero era en Perl, y fue hace más de una década: GD::Polilínea . Utilizaba curvas de Bézier cúbicas en 2D, y "suavizaba" un polígono o "polilínea" arbitraria (mi nombre entonces para lo que ahora se llama comúnmente una "cadena de líneas").

El algoritmo constaba de dos pasos: dados los puntos del Polígono, añadir dos puntos de control Bezier entre cada punto; luego llamar a un algoritmo simple para hacer una aproximación a trozos de la spline.

La segunda parte es fácil; la primera fue un poco de arte. Esta es la idea: considere un "segmento de control" un Vértice N: vN . El segmento de control era de tres puntos colineales: [cNa, vN, cNb] . El punto central era el vértice. La pendiente de este control seg era igual a la pendiente del Vértice N-1 al Vértice N+1. La longitud de la porción izquierda de este segmento era 1/3 de la longitud del Vértice N-1 al Vértice N, y la longitud de la porción derecha de este segmento era 1/3 de la longitud del Vértice N al Vértice N+1.

Si la curva original tenía cuatro vértices: [v1, v2, v3, v4] entonces el cada vértice ahora obtener un segmento de control de la forma: [c2a, v2, c2b] . Encadena esto así: [v1, c1b, c2a, v2, c2b, c3a, v3, c3b, c4a, v4] y los pica de cuatro en cuatro como los cuatro puntos de Bézier: [v1, c1b, c2a, v2] entonces [v2, c2b, c3a, v3] y así sucesivamente. Porque [c2a, v2, c2b] fueran colineales, la curva resultante será suave en cada vértice.

Así que esto también cumple con su requisito de parametrizar la "estrechez" de la curva: utilice un valor menor que 1/3 para una curva "más ajustada", uno mayor para un ajuste "más bucle". En cualquier caso, la curva resultante siempre pasa por los puntos originales dados.

El resultado fue una curva suave que "circunscribía" 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 tampoco tengo cifras. PERO... si/cuando lo porte a Python, me aseguraré de publicarlo aquí.