Comprendre le filtre de courbure de l'analyse de terrain raster QGIS?

J'ai lu le code source de plusieurs filtres raster QGis-1.7.4 calculant la pente, l'aspect et la courbure.

Il y a une formule dans le filtre calculant la courbure totale qui me laisse perplexe.

Le fichier source se trouve dans la version actuelle de QGis, avec le chemin suivant:

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

Le but de ce filtre est de calculer la courbure totale de la surface dans une fenêtre à neuf cellules. Le code de fonction est le suivant:

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

Je suis d'accord avec la formule "dxx" et avec l'expression de retour. Mais je pense que les formules "dyy" et "dxy" sont inversées: cela rend le résultat global asymétrique par rapport aux dimensions x et y.

Suis-je en train de manquer quelque chose ou dois-je remplacer les expressions dérivées doubles par:

  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

Pourriez-vous me donner votre avis sur ces formules, si elles sont incorrectes comme je le pensais ou si je me trompe? Dans ce dernier cas, savez-vous pourquoi les formules doivent être asymétriques par rapport à x et y?

Vos hypothèses sont correctes. Vérifier la symétrie est une excellente idée: la courbure (gaussienne) est une propriété intrinsèque d'une surface. Ainsi, la rotation d'une grille ne doit pas la modifier. Cependant, les rotations introduisent une erreur de discrétisation - à l'exception des rotations par multiples de 90 degrés. Par conséquent, une telle rotation doit conserver la courbure.

Nous pouvons comprendre ce qui se passe en capitalisant sur la toute première idée du calcul différentiel: les dérivés sont des limites de quotients de différence. C'est tout ce que nous devons vraiment savoir.

dxxest censé être une approximation discrète pour la dérivée seconde partielle dans la direction x. Cette approximation particulière (parmi les nombreuses possibles) est calculée en échantillonnant la surface le long d'un transect horizontal à travers la cellule. En localisant la cellule centrale à la ligne 2 et à la colonne 2, écrit (2,2), le transect passe à travers les cellules en (1,2), (2,2) et (3,2).

Le long de ce transect, les premières dérivées sont approximées par leurs quotients de différence, (* x32- * x22) / L et (* x22- * x12) / L où L est la distance (commune) entre les cellules (évidemment égale à cellSizeAvg). Les dérivées secondes sont obtenues par les quotients de différence de ceux-ci, donnant

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

Remarquez la division par L ^ 2!

De même, dyyest supposé être une approximation discrète pour la dérivée seconde partielle dans la direction y. Le transect est vertical, traversant les cellules en (2,1), (2,2) et (2,3). La formule sera la même que celle de dxxmais avec les indices transposés. Ce serait la troisième formule de la question - mais vous devez toujours diviser par L ^ 2.

La deuxième dérivée partielle mixte,, dxypeut être estimée en séparant les différences de deux cellules. Par exemple, la première dérivée par rapport à x à la cellule (2,3) (la cellule centrale supérieure, pas la cellule centrale!) Peut être estimée en soustrayant la valeur à sa gauche, * x13, de la valeur à sa droite, * x33, et en divisant par la distance entre ces cellules, 2L. La dérivée première par rapport à x à la cellule (2,1) (la cellule centrale du bas) est estimée par (* x31 - * x11) / (2L). Leur différence, divisée par 2L, estime le partiel mixte, donnant

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

Je ne suis pas vraiment sûr de ce que l'on entend par courbure "totale", mais il est probablement destiné à être la courbure gaussienne (qui est le produit des courbures principales). Selon Meek & Walton 2000 , équation 2.4, la courbure gaussienne est obtenue en divisant dxx * dyy - dxy ^ 2 (remarquez le signe moins! - c'est un déterminant ) par le carré de la norme du gradient de la surface. Ainsi, la valeur de retour citée dans la question n'est pas tout à fait une courbure, mais elle ressemble à une expression partielle foirée pour la courbure gaussienne.

On retrouve donc six erreurs dans le code , la plupart critiques:

1. dxx doit être divisé par L ^ 2, pas 1.

2. dyy doit être divisé par L ^ 2, pas 1.

3. Le signe de dxy est incorrect. (Cela n'a cependant aucun effet sur la formule de courbure.)

4. Les formules pour le colorant et le dxy sont mélangées, comme vous le constatez.

5. Un signe négatif manque dans un terme de la valeur de retour.

6. Il ne calcule pas réellement une courbure, mais seulement le numérateur d'une expression rationnelle pour la courbure.

Comme vérification très simple, vérifions que la formule modifiée renvoie des valeurs raisonnables pour les emplacements horizontaux sur les surfaces quadratiques. En prenant un tel emplacement pour être à l'origine du système de coordonnées, et en prenant son élévation pour être à une hauteur nulle, toutes ces surfaces ont des équations de la forme

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

pour les constantes a, b et c. Avec le carré central aux coordonnées (0,0), celui à sa gauche a les coordonnées (-L, 0), etc. Les neuf élévations sont

*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

D'où, par la formule modifiée,

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 courbure est estimée à 2a * 2c - (2b) ^ 2 = 4 (ac - b ^ 2). (Le dénominateur de la formule de Meek & Walton en est un dans ce cas.) Cela a-t-il un sens? Essayez quelques valeurs simples de a, b et c:

• a = c = 1, b = 0. Il s'agit d'un paraboloïde circulaire; sa courbure gaussienne devrait être positive. La valeur de 4 (ac-b ^ 2) est en effet positive (égale à 4).

• a = c = 0, b = 1. Il s'agit d'un hyperboloïde d'une feuille - une selle - l'exemple standard d'une surface de courbure négative . Effectivement, 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = -1. Ceci est une autre équation de l'hyperboloïde d'une feuille (tournée de 45 degrés). Encore une fois, 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = 0. Il s'agit d'une surface plane pliée en forme parabolique. Maintenant, 4 (ac-b ^ 2) = 0: la courbure gaussienne nulle détecte correctement la planéité de cette surface.

Si vous essayez le code de la question sur ces exemples, vous constaterez qu'il obtient toujours une valeur erronée.