Grundlegendes zum Krümmungsfilter der QGIS-Raster-Terrain-Analyse?

Ich habe den Quellcode mehrerer QGis-1.7.4-Rasterfilter gelesen, die Steigung, Aspekt und Krümmung berechnen.

Der Filter enthält eine Formel zur Berechnung der Gesamtkrümmung, die mich verwundert.

Die Quelldatei befindet sich in der aktuellen Version von QGis mit dem folgenden Pfad:

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

Der Zweck dieses Filters besteht darin, die Gesamtkrümmung der Oberfläche in einem Fenster mit neun Zellen zu berechnen. Der Funktionscode lautet wie folgt:

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

Ich bin mit der "dxx" Formel und mit dem Rückkehrausdruck okay. Ich denke jedoch, dass die "dyy" - und die "dxy" -Formel invertiert sind: Dies macht das Gesamtergebnis in Bezug auf x- und y-Dimensionen asymmetrisch.

Vermisse ich etwas oder soll ich die Ausdrücke für doppelte Ableitungen ersetzen durch:

  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

Könnten Sie mir Ihre Meinung zu diesen Formeln sagen, ob sie falsch sind, wie ich dachte, oder ob ich falsch liege? Wissen Sie im letzten Fall, warum die Formeln in Bezug auf x und y asymmetrisch sein müssen?

Ihre Vermutungen sind richtig. Die Prüfung auf Symmetrie ist eine hervorragende Idee: (Gaußsche) Krümmung ist eine intrinsische Eigenschaft einer Oberfläche. Das Drehen eines Gitters sollte es also nicht ändern. Drehungen führen jedoch zu Diskretisierungsfehlern - mit Ausnahme von Drehungen um ein Vielfaches von 90 Grad. Daher sollte eine solche Drehung die Krümmung bewahren.

Wir können verstehen, was vor sich geht, wenn wir die allererste Idee der Differentialrechnung nutzen: Derivate sind Grenzen von Differenzquotienten. Das ist alles was wir wirklich wissen müssen.

dxxsoll eine diskrete Näherung für die zweite partielle Ableitung in x-Richtung sein. Diese bestimmte Annäherung (aus den vielen möglichen) wird berechnet, indem die Oberfläche entlang eines horizontalen Durchgangs durch die Zelle abgetastet wird. Durch Lokalisieren der zentralen Zelle in Zeile 2 und Spalte 2, geschrieben (2,2), durchläuft der Transekt die Zellen in (1,2), (2,2) und (3,2).

Entlang dieses Transektors werden die ersten Ableitungen durch ihre Differenzquotienten (* x32- * x22) / L und (* x22- * x12) / L approximiert, wobei L der (gemeinsame) Abstand zwischen Zellen ist (offensichtlich gleich cellSizeAvg). Die zweiten Derivate werden durch die Differenzquotienten dieser erhalten, was ergibt

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

Beachten Sie die Division durch L ^ 2!

In ähnlicher Weise dyywird angenommen , eine diskrete Näherung für die zweite partielle Ableitung in Y-Richtung sein. Der Schnitt verläuft vertikal durch die Zellen in (2,1), (2,2) und (2,3). Die Formel sieht genauso aus wie für, dxxnur dass die Indizes transponiert sind. Das wäre die dritte Formel in der Frage - aber Sie müssen immer noch durch L ^ 2 dividieren.

Die gemischte zweite partielle Ableitung dxykann geschätzt werden, indem die Differenzen zweier Zellen auseinandergenommen werden. Zum Beispiel kann die erste Ableitung in Bezug auf x in Zelle (2,3) (die obere mittlere Zelle, nicht die zentrale Zelle!) Geschätzt werden, indem der Wert * x13 links von dem Wert * rechts davon subtrahiert wird. x33 und dividiert durch den Abstand zwischen diesen Zellen, 2L. Die erste Ableitung in Bezug auf x bei Zelle (2,1) (die untere mittlere Zelle) wird durch (* x31 - * x11) / (2L) geschätzt. Ihre Differenz, dividiert durch 2L, schätzt das gemischte partielle Geben

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

Ich bin mir nicht sicher, was unter "vollständiger" Krümmung zu verstehen ist, aber es ist wahrscheinlich die Gaußsche Krümmung (die das Produkt der Hauptkrümmungen ist). Nach Meek & Walton 2000 , Gleichung 2.4, erhält man die Gaußsche Krümmung, indem man dxx * dyy - dxy ^ 2 (beachten Sie das Minuszeichen! - dies ist eine Determinante ) durch das Quadrat der Norm des Gradienten der Oberfläche dividiert . Der in der Frage angegebene Rückgabewert ist also keine richtige Krümmung, sondern es sieht aus wie ein verzerrter Teilausdruck für die Gaußsche Krümmung.

Wir finden also sechs Fehler im Code , von denen die meisten kritisch sind:

1. dxx muss durch L ^ 2 geteilt werden, nicht durch 1.

2. dyy muss durch L ^ 2 geteilt werden, nicht durch 1.

3. Das Vorzeichen von dxy ist falsch. (Dies hat jedoch keine Auswirkung auf die Krümmungsformel.)

4. Wie Sie bemerken, sind die Formeln für dyy und dxy gemischt.

5. Ein negatives Vorzeichen fehlt in einem Term im Rückgabewert.

6. Es berechnet eigentlich keine Krümmung, sondern nur den Zähler eines rationalen Ausdrucks für die Krümmung.

Vergewissern Sie sich zur Vereinfachung, dass die geänderte Formel sinnvolle Werte für horizontale Positionen auf quadratischen Flächen zurückgibt. Wenn ein solcher Ort als Ursprung des Koordinatensystems und seine Höhe als Höhe Null angenommen wird, haben alle diese Flächen Formgleichungen

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

für Konstante a, b und c. Wenn sich das zentrale Quadrat bei den Koordinaten (0,0) befindet, hat das linke die Koordinaten (-L, 0) usw. Die neun Höhen sind

*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

Woher, durch die modifizierte Formel,

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.

Die Krümmung wird mit 2a * 2c - (2b) ^ 2 = 4 (ac - b ^ 2) geschätzt. (Der Nenner in der Meek & Walton-Formel ist in diesem Fall einer.) Ist das sinnvoll? Probieren Sie einige einfache Werte von a, b und c aus:

• a = c = 1, b = 0. Dies ist ein kreisförmiges Paraboloid; seine Gaußsche Krümmung sollte positiv sein. Der Wert von 4 (ac-b ^ 2) ist in der Tat positiv (gleich 4).

• a = c = 0, b = 1. Dies ist ein Hyperboloid eines Blattes - ein Sattel - das Standardbeispiel für eine Oberfläche mit negativer Krümmung. Sicher genug, 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = -1. Dies ist eine andere Gleichung des Hyperboloids eines Blatts (um 45 Grad gedreht). Nochmals 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = 0. Dies ist eine flache Oberfläche, die in eine parabolische Form gefaltet ist. Jetzt ist 4 (ac-b ^ 2) = 0: Die Null-Gaußsche Krümmung erkennt die Ebenheit dieser Oberfläche korrekt.

Wenn Sie den Code in der Frage in diesen Beispielen ausprobieren, werden Sie feststellen, dass er immer einen falschen Wert erhält.