Entendendo o filtro de curvatura da análise de terreno raster QGIS?

Eu li o código fonte de vários filtros raster QGis-1.7.4, computando inclinação, aspecto e curvatura.

Existe uma fórmula no filtro que calcula a curvatura total que me intriga.

O arquivo de origem está na versão atual do QGis, com o seguinte caminho:

qgis-1.7.4 / src / analysis / raster / qgstotalcurvaturefilter.cpp

O objetivo desse filtro é calcular a curvatura total da superfície em uma janela de nove células. O código da função é o seguinte:

float QgsTotalCurvatureFilter::processNineCellWindow( 
   float* x11, float* x21, float* x31, 
   float* x12, float* x22, float* x32, 
   float* x13, float* x23, float* x33 ) {

  ... some code deleted ...

  double dxx = ( *x32 - 2 * *x22 + *x12 ) / ( 1 );
  double dyy = ( -*x11 + *x31 + *x13 - *x33 ) / ( 4 * cellSizeAvg*cellSizeAvg );
  double dxy = ( *x21 - 2 * *x22 + *x23 ) / ( 1 );

  return dxx*dxx + 2*dxy*dxy + dyy*dyy;
}

Eu estou bem com a fórmula "dxx" e com a expressão de retorno. Mas acho que as fórmulas "dyy" e "dxy" são invertidas: isso torna o resultado total assimétrico em relação às dimensões xey.

Estou faltando alguma coisa ou devo substituir as expressões derivadas duplas por:

  double dxx = ( *x32 - 2 * *x22 + *x12 ) / ( 1 ); // unchanged
  // inversion of the two following:
  double dxy = ( -*x11 + *x31 + *x13 - *x33 ) / ( 4 * cellSizeAvg*cellSizeAvg );
  double dyy = ( *x21 - 2 * *x22 + *x23 ) / ( 1 );
  return dxx*dxx + 2*dxy*dxy + dyy*dyy; // unchanged

Você poderia me dizer sua opinião sobre essas fórmulas, se estiverem incorretas como eu pensava ou se estiver errado? Neste último caso, você sabe por que as fórmulas precisam ser assimétricas em relação a x e y?

Suas suposições estão corretas. Verificar a simetria é uma excelente idéia: a curvatura (gaussiana) é uma propriedade intrínseca de uma superfície. Portanto, girar uma grade não deve alterá-la. No entanto, as rotações introduzem erros de discretização - exceto as rotações em múltiplos de 90 graus. Portanto, qualquer rotação desse tipo deve preservar a curvatura.

Podemos entender o que está acontecendo capitalizando a primeira idéia do cálculo diferencial: derivadas são limites de quocientes de diferença. É tudo o que realmente precisamos saber.

dxxdeve ser uma aproximação discreta para a segunda derivada parcial na direção x. Essa aproximação específica (dentre as muitas possíveis) é calculada amostrando a superfície ao longo de um transecto horizontal através da célula. Localizando a célula central na linha 2 e na coluna 2, escrita (2,2), o transecto passa pelas células em (1,2), (2,2) e (3,2).

Ao longo desse transecto, as primeiras derivadas são aproximadas por seus quocientes de diferença (* x32- * x22) / L e (* x22- * x12) / L, onde L é a distância (comum) entre as células (evidentemente igual a cellSizeAvg). Os segundos derivativos são obtidos pelos quocientes de diferença destes, produzindo

dxx = ((*x32-*x22)/L - (*x22-*x12)/L)/L
    = (*x32 - 2 * *x22 + *x12) / L^2.

Observe a divisão por L ^ 2!

Da mesma forma, dyyé suposto ser uma aproximação discreta para a segunda derivada parcial na direção y. O transecto é vertical, passando pelas células em (2,1), (2,2) e (2,3). A fórmula terá a mesma aparência de dxxmas com os subscritos transpostos. Essa seria a terceira fórmula da pergunta - mas você ainda precisa dividir por L ^ 2.

A segunda derivada parcial mista,, dxypode ser estimada separando-se duas células. Por exemplo, a primeira derivada com relação a x na célula (2,3) (a célula do meio superior, não a célula central!) Pode ser estimada subtraindo o valor à sua esquerda, * x13, do valor à direita, * x33 e dividindo pela distância entre essas células, 2L. A primeira derivada em relação a x na célula (2,1) (a célula do meio do fundo) é estimada por (* x31 - * x11) / (2L). Sua diferença, dividida por 2L, estima a parcial mista, dando

dxy = ((*x33 - *x13)/(2L) - (*x31 - *x11)/(2L))/(2L)
    = (*x33 - *x13 - *x31 + *x11) / (4 L^2).

Não tenho muita certeza do que se entende por curvatura "total", mas provavelmente deve ser a curvatura gaussiana (que é o produto das principais curvaturas). De acordo com Meek & Walton 2000 , Equação 2.4, a curvatura gaussiana é obtida dividindo-se dxx * dyy-dxy ^ 2 (observe o sinal de menos! - este é um determinante ) pelo quadrado da norma do gradiente da superfície. Assim, o valor de retorno citado na pergunta não é exatamente uma curvatura, mas parece uma expressão parcial confusa da curvatura gaussiana.

Encontramos, então, seis erros no código , a maioria deles críticos:

1. dxx precisa ser dividido por L ^ 2, não 1.

2. dyy precisa ser dividido por L ^ 2, não 1.

3. O sinal do dxy está incorreto. (Porém, isso não afeta a fórmula da curvatura.)

4. As fórmulas para dyy e dxy são misturadas, como você observa.

5. Um sinal negativo está ausente de um termo no valor de retorno.

6. Na verdade, ele não calcula uma curvatura, mas apenas o numerador de uma expressão racional para a curvatura.

Como uma verificação muito simples, vamos verificar se a fórmula modificada retorna valores razoáveis ​​para locais horizontais em superfícies quadráticas. Tomando esse local como a origem do sistema de coordenadas e levando a sua elevação a uma altura zero, todas essas superfícies têm equações da forma

elevation = a*x^2 + 2b*x*y + c*y^2.

para constante a, bec. Com o quadrado central nas coordenadas (0,0), o quadrado à esquerda possui coordenadas (-L, 0) etc. As nove elevações são

*x13 *x23 *x33     (a-2b+c)L^2, (c)L^2, (a+2b+c)L^2
*x12 *x22 *x32  =  (a)L^2,      0,      (a)L^2
*x11 *x21 *x31     (a+2b+c)L^2, (c)L^2, (a-2b+c)L^2

De onde, pela fórmula modificada,

dxx = (a*L^2 - 2*0 + a*L^2) / L^2
    = 2a;

dxy = ((a+2b+c)L^2 - (a-2b+c)L^2 - (a-2b+c)L^2 + (a+2b+c)L^2)/(4L^2)
    = 2b;

dyy = ... [computed as in dxx] ... = 2c.

A curvatura é estimada como 2a * 2c - (2b) ^ 2 = 4 (ac - b ^ 2). (O denominador na fórmula de Meek & Walton é um neste caso.) Isso faz sentido? Tente alguns valores simples de a, bec:

• a = c = 1, b = 0. Este é um parabolóide circular; sua curvatura gaussiana deve ser positiva. O valor de 4 (ac-b ^ 2) é de fato positivo (igual a 4).

• a = c = 0, b = 1. Este é um hiperboloide de uma folha - uma sela - o exemplo padrão de uma superfície de curvatura negativa . Com certeza, 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = -1. Essa é outra equação do hiperboloide de uma folha (girada 45 graus). Mais uma vez, 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = 0. Esta é uma superfície plana dobrada em uma forma parabólica. Agora, 4 (ac-b ^ 2) = 0: a curvatura Gaussiana zero detecta corretamente o nivelamento dessa superfície.

Se você tentar o código da pergunta nesses exemplos, descobrirá que ele sempre recebe um valor incorreto.