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:
dxx należy podzielić przez L ^ 2, a nie 1.
dyy należy podzielić przez L ^ 2, a nie 1.
Znak dxy jest niepoprawny. (Nie ma to jednak wpływu na wzór krzywizny.)
Jak zauważyłeś, formuły dla barwników i dxy są pomieszane.
Brak znaku ujemnego w terminie w wartości zwracanej.
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ść.