Suavizando polígonos no mapa de contorno?

52

Aqui está um mapa de contorno para o qual todos os polígonos dos níveis estão disponíveis.

Vamos perguntar como suavizar os polígonos, mantendo todos os vértices preservados em seus locais exatos?

Na verdade, o contorno é feito sobre os dados da grade, você pode sugerir suavizar os dados da grade e, portanto, o contorno resultante será mais suave. Observe que isso não está funcionando como meu desejo, pois a função de suavização, como o filtro Gaussian, removerá pequenos pacotes de dados e alterará o alcance da terceira variável, por exemplo, a altura que não é permitida no meu aplicativo.

Na verdade, estou procurando um pedaço de código (de preferência em Python ) que possa suavizar os polígonos 2D (qualquer tipo: convexo, côncavo, auto-interceptável etc.) razoavelmente indolor (esqueça as páginas de códigos) e preciso.

Para sua informação, existe uma função no ArcGIS que faz isso perfeitamente, mas o uso de aplicativos comerciais de terceiros não é minha escolha para esta pergunta.

insira a descrição da imagem aqui


1)

Scipy.interpolate:

insira a descrição da imagem aqui

Como você vê, as splines resultantes (vermelhas) não são satisfatórias!

2)

Aqui está o resultado usando o código fornecido aqui . Não está funcionando bem!

insira a descrição da imagem aqui

3)

Para mim, a melhor solução deve ser algo como a figura a seguir, na qual um quadrado está sendo suavizado gradualmente, alterando apenas um valor. Espero um conceito semelhante para suavizar qualquer forma de polígonos.

insira a descrição da imagem aqui

Satisfazer a condição que o spline passa nos pontos:

insira a descrição da imagem aqui

4)

Aqui está minha implementação da "idéia do whuber's" linha por linha no Python em seus dados. Possivelmente existem alguns erros, pois os resultados não são bons.

insira a descrição da imagem aqui

K = 2 é um desastre e, portanto, para k> = 4.

5)

Eu removi um ponto no local problemático e o spline resultante agora é idêntico ao do whuber's. Mas ainda é uma pergunta por que o método não funciona para todos os casos?

insira a descrição da imagem aqui

6)

Uma boa suavização para os dados do whuber pode ser a seguinte (desenhada pelo software de gráficos vetoriais), na qual um ponto extra foi adicionado sem problemas (compare com a atualização

4):

insira a descrição da imagem aqui

7)

Veja o resultado da versão Python do código do whuber para obter algumas formas icônicas:

insira a descrição da imagem aqui
Observe que o método parece não funcionar para polilinhas. Para o canto polilinha (contorno), verde é o que eu quero, mas fiquei vermelho. Isso precisa ser resolvido, pois os mapas de contorno são sempre polilinhas, embora as polilinhas fechadas possam ser tratadas como polígonos, como nos meus exemplos. Também não é que o problema surgido na atualização 4 ainda não tenha sido resolvido.

8) [meu último]

Aqui está a solução final (não perfeita!):

insira a descrição da imagem aqui

Lembre-se de que você terá que fazer algo sobre a área apontada por estrelas. Talvez haja um erro no meu código ou o método proposto precise de mais desenvolvimento para considerar todas as situações e fornecer os resultados desejados.

desenvolvedor
fonte
como você está gerando contornos de 'polígono'? elas não seriam sempre linhas, uma vez que um contorno que cruza a borda de um DEM nunca se fecha?
Pistachionut
Eu usei a função v.generalize no GRASS para suavizar as linhas de contorno com resultados decentes, embora possa levar algum tempo para mapas com contornos muito densos.
Pistachionut
@ pistachionut Você pode considerar que os níveis de contorno são polilinhas. Estou procurando um código puro no primeiro estágio. Se não estiver disponível, leve o pacote para Python.
Desenvolvedor
Talvez veja scipy.org/Cookbook/Interpolation porque parece que você deseja spline
PolyGeo
11
A curva @Pablo Bezier no seu link funciona bem para polilinhas. whuber's funciona quase bem para polígonos. Para que juntos pudessem resolver a questão. Muito obrigado por compartilhar seu conhecimento gratuitamente.
Desenvolvedor

Respostas:

37

A maioria dos métodos para spline seqüências de números spline polígonos. O truque é fazer com que os splines "fechem" suavemente nos pontos finais. Para fazer isso, "envolva" os vértices nas extremidades. Em seguida, divine as coordenadas x e y separadamente.

Aqui está um exemplo de trabalho em R. Ele usa o splineprocedimento cúbico padrão disponível no pacote de estatísticas básicas. Para obter mais controle, substituto quase qualquer procedimento que você prefere: apenas certifique-se estrias através dos números (isto é, interpola eles) em vez de simplesmente usá-los como "pontos de controle."

#
# 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 seu uso, vamos criar um polígono pequeno (mas 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)

Divine-o usando o código anterior. Para tornar o spline mais suave, aumente o número de vértices de 100; para torná-lo menos suave, diminua o número de vértices.

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

Para ver os resultados, plotamos (a) o polígono original em vermelho tracejado, mostrando o espaço entre o primeiro e o último vértices (isto é, não fechando sua polilinha de limite); e (b) o spline em cinza, mais uma vez mostrando seu espaço. (Como a diferença é muito pequena, seus pontos finais são destacados com pontos azuis).

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)

Polígono estriado

whuber
fonte
5
Boa resposta. Existe alguma maneira de garantir que os contornos não acabem cruzando como resultado da suavização?
Kirk Kuykendall
Essa é uma boa pergunta, @Kirk. Não conheço nenhum método para garantir a não-passagem dessa forma de suavização. (Na verdade, eu nem vejo como garantir que a polilinha suavizada não se intercepte. Esse não é um grande problema para a maioria dos contornos.) Para fazer isso, você precisaria retornar ao original. DEM e, em vez disso, use um método melhor para calcular os contornos em primeiro lugar. (Não são melhores métodos - eles são conhecidos há muito tempo -. Mas AFAIK alguns dos mais populares GISes não usá-los)
whuber
Primeiro, ainda estou trabalhando para implementar sua resposta em Python, mas sem êxito. Segundo, qual será o resultado se você aplicar seu método em um quadrado? Você pode se referir àqueles que desenhei na pergunta.
Desenvolvedor
11
Aceitei isso como resposta, pois fornece uma boa solução. Embora não seja perfeito, mas me deu algumas idéias, espero encontrar uma solução que satisfaça os pontos mencionados acima na minha pergunta e nos comentários. Você também pode considerar os comentários do whuber para a pergunta [CQ], existem bons truques por lá. Por fim, devo dizer que a tradução para python é quase direta, com o adorável pacote Scipy instalado. Considere também o comentário de Pablo no CQ como uma possível solução para polilinhas, isto é, curvas de Bezier. Boa sorte a todos.
Desenvolvedor
11
vendo suas respostas, lamento não cuidar bem da minha matemática !!!
Vinayan #
2

Sei que este é um post antigo, mas apareceu no Google para algo que eu estava procurando, então pensei em publicar minha solução.

Não vejo isso como um exercício de ajuste de curva 2D, mas sim em 3D. Ao considerar os dados em 3D, podemos garantir que as curvas nunca se cruzem e podemos usar informações de outros contornos para melhorar nossa estimativa da atual.

O seguinte extrato do iPython usa interpolação cúbica fornecida pelo SciPy. Observe que os valores z que plotei não são importantes, desde que todos os contornos tenham altura eqüidistante.

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

Resultado interpolado cúbico

Os resultados aqui não parecem os melhores, mas com tão poucos pontos de controle, eles ainda são perfeitamente válidos. Observe como a linha verde ajustada é puxada para seguir o contorno azul mais amplo.

Gilly
fonte
As curvas suaves ajustadas devem permanecer o mais próximo possível do polígono / polilinha original.
Desenvolvedor
1

Eu escrevi quase exatamente o pacote que você está procurando ... mas estava no Perl e há mais de uma década: GD :: Polyline . Ele usava curvas cúbicas de Bezier em 2D e "suavizava" um polígono ou "polilinha" arbitrária (meu nome para o que agora é chamado de "LineString").

O algoritmo foi composto de duas etapas: dados os pontos no polígono, adicione dois pontos de controle de Bezier entre cada ponto; em seguida, chame um algoritmo simples para fazer uma aproximação por partes do spline.

A segunda parte é fácil; a primeira parte foi um pouco de arte. Aqui foi o insight: considerar um "segmento de controle" a Vertex N: vN. O segmento de controlo era três pontos co-lineares: [cNa, vN, cNb]. O ponto central era o vértice. A inclinação deste controle seg era igual à inclinação do Vertex N-1 ao Vertex N + 1. O comprimento da porção esquerda deste segmento era 1/3 do comprimento do vértice N-1 ao vértice N, e o comprimento da porção direita desse segmento era 1/3 do comprimento do vértice N ao vértice N + 1.

Se a curva original quatro vértices: [v1, v2, v3, v4]então a cada vértice agora obter um segmento da forma de controle: [c2a, v2, c2b]. Junte-os assim: [v1, c1b, c2a, v2, c2b, c3a, v3, c3b, c4a, v4]e mastigue-os quatro de cada vez, como os quatro pontos de Bezier:, [v1, c1b, c2a, v2]então [v2, c2b, c3a, v3], e assim por diante. Por [c2a, v2, c2b]serem co-lineares, a curva resultante será suave em cada vértice.

Portanto, isso também atende aos seus requisitos para parametrizar o "aperto" da curva: use um valor menor que 1/3 para uma curva "mais apertada", e um valor maior para um ajuste "loopier". Em qualquer um dos casos, a curva resultante sempre passa pelos pontos dados originais.

Isso resultou em uma curva suave que "circunscreveu" o polígono original. Eu também tinha uma maneira de "inscrever" uma curva suave ... mas não vejo isso no código CPAN.

De qualquer forma, no momento não tenho uma versão disponível em Python, nem tenho números. MAS ... se / quando eu portar isso para Python, não deixe de postar aqui.

Dan H
fonte
Não é possível avaliar o código Perl, adicione gráficos para demonstrar como estava funcionando, se possível.
Desenvolvedor