¿Comprender el filtro de curvatura del análisis de terreno ráster QGIS?

He leído el código fuente de varios filtros ráster QGis-1.7.4 que calculan la pendiente, el aspecto y la curvatura.

Hay una fórmula en el filtro que calcula la curvatura total que me desconcierta.

El archivo fuente está en la versión actual de QGis, con la siguiente ruta:

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

El propósito de este filtro es calcular la curvatura total de la superficie en una ventana de nueve celdas. El código de función es el siguiente:

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;
}

Estoy de acuerdo con la fórmula "dxx" y con la expresión de retorno. Pero creo que las fórmulas "dyy" y "dxy" están invertidas: esto hace que el resultado total sea asimétrico con respecto a las dimensiones x e y.

¿Me estoy perdiendo algo o debo reemplazar las expresiones de derivadas dobles 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

¿Podría decirme su opinión sobre estas fórmulas, si son incorrectas como pensaba o si estoy equivocado? En este último caso, ¿sabe por qué las fórmulas tienen que ser asimétricas con respecto a x e y?

Tus conjeturas son correctas. Verificar la simetría es una excelente idea: la curvatura (gaussiana) es una propiedad intrínseca de una superficie. Por lo tanto, girar una cuadrícula no debería cambiarla. Sin embargo, las rotaciones introducen un error de discretización, excepto las rotaciones por múltiplos de 90 grados. Por lo tanto, cualquier rotación de este tipo debería preservar la curvatura.

Podemos entender lo que está sucediendo al capitalizar la primera idea del cálculo diferencial: los derivados son límites de cocientes de diferencia. Eso es todo lo que realmente necesitamos saber.

dxxse supone que es una aproximación discreta para la segunda derivada parcial en la dirección x. Esta aproximación particular (de las muchas posibles) se calcula muestreando la superficie a lo largo de un transecto horizontal a través de la celda. Al ubicar la celda central en la fila 2 y la columna 2, escrita (2,2), el transecto pasa a través de las celdas en (1,2), (2,2) y (3,2).

A lo largo de este transecto, las primeras derivadas se aproximan por sus cocientes de diferencia, (* x32- * x22) / L y (* x22- * x12) / L donde L es la distancia (común) entre celdas (evidentemente igual a cellSizeAvg). Los segundos derivados se obtienen por los cocientes de diferencia de estos, produciendo

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

¡Observe la división por L ^ 2!

Del mismo modo, dyyse supone que es una aproximación discreta para la segunda derivada parcial en la dirección y. El transecto es vertical, pasando a través de las celdas en (2,1), (2,2) y (2,3). La fórmula se verá igual que para dxxpero con los subíndices transpuestos. Esa sería la tercera fórmula en la pregunta, pero aún necesita dividir por L ^ 2.

La segunda derivada parcial mixta dxy, se puede estimar separando las diferencias entre dos celdas. Por ejemplo, la primera derivada con respecto a x en la celda (2,3) (¡la celda central superior, no la celda central!) Puede estimarse restando el valor a su izquierda, * x13, del valor a su derecha, * x33, y dividiendo por la distancia entre esas celdas, 2L. La primera derivada con respecto a x en la celda (2,1) (la celda central inferior) se estima mediante (* x31 - * x11) / (2L). Su diferencia, dividida por 2L, estima la mezcla parcial, dando

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

No estoy realmente seguro de qué se entiende por curvatura "total", pero probablemente se pretende que sea la curvatura gaussiana (que es el producto de las curvaturas principales). Según Meek & Walton 2000 , Ecuación 2.4, la curvatura gaussiana se obtiene dividiendo dxx * dyy - dxy ^ 2 (¡observe el signo menos! - esto es un determinante ) por el cuadrado de la norma del gradiente de la superficie. Por lo tanto, el valor de retorno citado en la pregunta no es una curvatura, pero parece una expresión parcial desordenada para la curvatura gaussiana.

Encontramos, entonces, seis errores en el código , la mayoría de ellos críticos:

1. dxx debe dividirse entre L ^ 2, no 1.

2. dyy necesita ser dividido por L ^ 2, no 1.

3. El signo de dxy es incorrecto. (Sin embargo, esto no tiene ningún efecto en la fórmula de curvatura).

4. Las fórmulas para dyy y dxy se mezclan, como se observa.

5. Falta un signo negativo de un término en el valor de retorno.

6. En realidad, no calcula una curvatura, sino solo el numerador de una expresión racional para la curvatura.

Como una comprobación muy simple, verifiquemos que la fórmula modificada devuelva valores razonables para ubicaciones horizontales en superficies cuadráticas. Tomando tal ubicación como el origen del sistema de coordenadas, y tomando su elevación como a altura cero, todas esas superficies tienen ecuaciones de la forma

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

para constante a, by c. Con el cuadrado central en las coordenadas (0,0), el que está a su izquierda tiene coordenadas (-L, 0), etc. Las nueve elevaciones son

*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 donde, por la 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.

La curvatura se estima como 2a * 2c - (2b) ^ 2 = 4 (ac - b ^ 2). (El denominador en la fórmula de Meek & Walton es uno en este caso). ¿Tiene sentido? Pruebe algunos valores simples de a, byc:

• a = c = 1, b = 0. Este es un paraboloide circular; su curvatura gaussiana debe ser positiva. El valor de 4 (ac-b ^ 2) de hecho es positivo (igual a 4).

• a = c = 0, b = 1. Este es un hiperboloide de una hoja - una silla de montar - el ejemplo estándar de una superficie de curvatura negativa . Efectivamente, 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = -1. Esta es otra ecuación del hiperboloide de una hoja (rotada 45 grados). Una vez más, 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = 0. Esta es una superficie plana doblada en forma parabólica. Ahora, 4 (ac-b ^ 2) = 0: la curvatura gaussiana cero detecta correctamente la planitud de esta superficie.

Si prueba el código en la pregunta de estos ejemplos, encontrará que siempre obtiene un valor erróneo.