Czy rozumiesz filtr krzywizny analizy terenu rastrowego QGIS?

Przeczytałem kod źródłowy kilku filtrów rastrowych QGis-1.7.4 obliczających nachylenie, aspekt i krzywiznę.

Filtr zawiera formułę obliczającą całkowitą krzywiznę, która mnie zastanawia.

Plik źródłowy znajduje się w bieżącej wersji QGis, z następującą ścieżką:

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

Celem tego filtra jest obliczenie całkowitej krzywizny powierzchni w oknie dziewięciu komórek. Kod funkcji jest następujący:

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

Nie przeszkadza mi formuła „dxx” i wyrażenie zwrotne. Ale myślę, że formuły „dyy” i „dxy” są odwrócone: to sprawia, że ​​całkowity wynik jest asymetryczny w odniesieniu do wymiarów xiy.

Czy coś mi brakuje, czy powinienem zastąpić wyrażenia podwójnych pochodnych przez:

  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

Czy możesz mi powiedzieć swoją opinię na temat tych formuł, jeśli są one nieprawidłowe, jak myślałem, lub jeśli się mylę? Czy w tym ostatnim przypadku wiesz, dlaczego formuły muszą być asymetryczne względem xiy?

Twoje przypuszczenia są prawidłowe. Sprawdzanie symetrii to doskonały pomysł: krzywizna (gaussowska) jest nieodłączną właściwością powierzchni. Dlatego obracanie siatki nie powinno jej zmieniać. Jednak obroty wprowadzają błąd dyskretyzacji - z wyjątkiem obrotów o wielokrotności 90 stopni. Dlatego każdy taki obrót powinien zachować krzywiznę.

Możemy zrozumieć, co się dzieje , wykorzystując pierwszą ideę rachunku różniczkowego: pochodne są granicami ilorazów różnic. To wszystko, co naprawdę musimy wiedzieć.

dxxma być dyskretnym przybliżeniem drugiej pochodnej cząstkowej w kierunku x. To szczególne przybliżenie (spośród wielu możliwych) jest obliczane przez próbkowanie powierzchni wzdłuż poziomego transektu przez komórkę. Lokalizując centralną komórkę w wierszu 2 i kolumnie 2, zapisaną (2,2), transekty przechodzą przez komórki w (1,2), (2,2) i (3,2).

Wzdłuż tego transetu pierwsze pochodne są aproksymowane przez ich iloraz różnic, (* x32- * x22) / L i (* x22- * x12) / L, gdzie L jest (wspólną) odległością między komórkami (ewidentnie równą cellSizeAvg). Drugie pochodne otrzymuje się przez iloraz różnic tych, uzyskując

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

Zwróć uwagę na podział według L ^ 2!

Podobnie dyyjest powinien być dyskretne przybliżenie drugiego częściowego pochodnej w kierunku y. Transekt jest pionowy, przechodząc przez komórki w (2,1), (2,2) i (2,3). Formuła będzie wyglądać tak samo jak ta, dxxale z transponowanymi indeksami dolnymi. To byłaby trzecia formuła w pytaniu - ale nadal musisz podzielić przez L ^ 2.

Mieszaną drugą pochodną cząstkową dxymożna oszacować, biorąc różnice między dwiema komórkami. Np. Pierwszą pochodną względem x w komórce (2,3) (górna środkowa komórka, a nie centralna komórka!) Można oszacować, odejmując wartość po jej lewej stronie, * x13, od wartości po prawej stronie, * x33 i dzieląc przez odległość między tymi komórkami, 2L. Pierwsza pochodna w odniesieniu do x w komórce (2,1) (dolna środkowa komórka) jest szacowana przez (* x31 - * x11) / (2L). Ich różnica, podzielona przez 2L, szacuje mieszane częściowe, dając

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

Nie jestem do końca pewien, co należy rozumieć przez „całkowitą” krzywiznę, ale prawdopodobnie ma to być krzywizna Gaussa (która jest produktem głównych krzywizn). Według Meek & Walton 2000 , równanie 2.4, krzywiznę Gaussa uzyskuje się dzieląc dxx * dyy - dxy ^ 2 (zauważ znak minus! - jest to wyznacznik ) przez kwadrat normy gradientu powierzchni. Zatem wartość zwrotna cytowana w pytaniu nie jest do końca krzywizną, ale wygląda na zepsute wyrażenie częściowe dla krzywizny Gaussa.

Mamy zatem sześć błędów w kodzie , z których większość ma krytyczne znaczenie:

1. dxx należy podzielić przez L ^ 2, a nie 1.

2. dyy należy podzielić przez L ^ 2, a nie 1.

3. Znak dxy jest niepoprawny. (Nie ma to jednak wpływu na wzór krzywizny.)

4. Jak zauważyłeś, formuły dla barwników i dxy są pomieszane.

5. Brak znaku ujemnego w terminie w wartości zwracanej.

6. W rzeczywistości nie oblicza krzywizny, a jedynie licznik racjonalnego wyrażenia krzywizny.

Jako bardzo proste sprawdzenie, sprawdźmy, czy zmodyfikowana formuła zwraca rozsądne wartości dla poziomych lokalizacji na powierzchniach kwadratowych. Przyjmując takie położenie za początek układu współrzędnych i przyjmując, że jego rzędna ma być na zerowej wysokości, wszystkie takie powierzchnie mają równania kształtu

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

dla stałej a, b i c. Z centralnym kwadratem o współrzędnych (0,0), ten po jego lewej stronie ma współrzędne (-L, 0) itp. Dziewięć rzędnych to

*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

Skąd, według zmodyfikowanej formuły,

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.

Krzywizna jest szacowana jako 2a * 2c - (2b) ^ 2 = 4 (ac - b ^ 2). (W tym przypadku mianownik w formule Meek & Walton jest jeden.) Czy to ma sens? Wypróbuj kilka prostych wartości a, b i c:

• a = c = 1, b = 0. Jest to okrągły paraboloid; jego krzywizna Gaussa powinna być dodatnia. Wartość 4 (ac-b ^ 2) rzeczywiście jest dodatnia (równa 4).

• a = c = 0, b = 1. Jest to hiperboloid jednego arkusza - siodła - standardowy przykład powierzchni krzywizny ujemnej . Jasne, 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = -1. To kolejne równanie hiperboloidu jednego arkusza (obróconego o 45 stopni). Jeszcze raz 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = 0. Jest to płaska powierzchnia złożona w kształt paraboliczny. Teraz 4 (ac-b ^ 2) = 0: zerowa krzywizna Gaussa poprawnie wykrywa płaskość tej powierzchni.

Jeśli wypróbujesz kod w pytaniu na tych przykładach, przekonasz się, że zawsze otrzymuje on błędną wartość.