Su conjetura es correcta. La comprobación de la simetría es una excelente idea: (Gauss) la curvatura es una propiedad intrínseca de una superficie. Por lo tanto, la rotación de un grid no debe cambiar. Sin embargo, las rotaciones de introducir error de discretización, excepto las rotaciones por los múltiplos de 90 grados. Por lo tanto, cualquier rotación debe preservar la curvatura.
Podemos entender lo que está pasando , aprovechando la primera idea de cálculo diferencial: los derivados son los límites de la diferencia de cocientes. Eso es todo lo que realmente necesita saber.
dxx supone una aproximación discreta de la segunda derivada parcial en la dirección x. Esta particular aproximación (de los muchos posibles) se calcula mediante el muestreo de la superficie a lo largo de un eje horizontal transecto a través de la célula. La localización de la central celda de la fila 2 y la columna 2, escrito (2,2), el eje que pasa a través de las células en (1,2), (2,2), (3,2).
A lo largo de este transecto, el primero de los derivados se aproximan por su diferencia de cocientes, (*x32-*x22)/L y (*x22-*x12)/L, donde L es la (común) de la distancia entre las células (evidentemente igual a cellSizeAvg). El segundo de los derivados se obtienen por la diferencia de los cocientes de estas, dando
dxx = ((*x32-*x22)/L - (*x22-*x12)/L)/L
= (*x32 - 2 * *x22 + *x12) / L^2.
Aviso de la división por L^2!
Del mismo modo, dyy está supuesta a ser una aproximación discreta de la segunda derivada parcial en la dirección del eje y. El corte es vertical, pasando a través de las células en (2,1), (2,2), (2,3). La fórmula tendrá el mismo aspecto que el de dxx pero con los subíndices transpuesto. Que sería la tercera fórmula, en la pregunta--pero usted todavía tiene que dividir por L^2.
La mezcla de la segunda derivada parcial, dxy, se puede estimar diferencias de dos células separadas. E. g., la primera derivada con respecto a x en la celda (2,3) (la parte superior central de la célula, no la celda central!) se puede calcular restando el valor a su izquierda, *x13, a partir del valor en su derecho, *x33, y dividiendo por la distancia entre esas células, 2L. La primera derivada con respecto a x en la celda (2,1) (la parte central de la célula) se calcula por (*x31 - *x11)/(2L). Su diferencia, dividida por 2L, las estimaciones de la mezcla parcial, dando
dxy = ((*x33 - *x13)/(2L) - (*x31 - *x11)/(2L))/(2L)
= (*x33 - *x13 - *x31 + *x11) / (4 L^2).
No estoy realmente seguro de lo que se entiende por "total", la curvatura, pero es probablemente la intención de ser la curvatura de Gauss (que es el producto de las curvaturas principales). Según Manso Y Walton 2000, la Ecuación 2.4, la Gaussiana de la curvatura se obtiene dividiendo dxx * dyy - dxy^2 (observe el signo menos!--este es un factor 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 no se ve como un mal parcial de la expresión de la curvatura de Gauss.
Nos encontramos, entonces, seis errores en el código, la mayoría de ellas críticas:
dxx necesita ser dividida por L^2, no 1.
dyy necesita ser dividida por L^2, no 1.
El signo de dxy es incorrecta. (Esto no tiene ningún efecto sobre la curvatura de la fórmula, sin embargo.)
Las fórmulas para dyy y dxy se mezclan, como nota.
Un signo negativo es la falta de un término en el valor de retorno.
En realidad no calcular la curvatura, pero sólo el numerador de una expresión racional de la curvatura.
Como una forma muy sencilla de comprobar, vamos a comprobar que la modificación de la fórmula devuelve valores razonables para ubicaciones horizontales en superficies cuadráticas. Tomando un lugar para ser el origen del sistema de coordenadas, y tomando su elevación a estar en el cero de altura, todas las superficies tienen ecuaciones de la forma
elevation = a*x^2 + 2b*x*y + c*y^2.
para las constantes a, b, y c de. Con la plaza central en las coordenadas (0,0), el que a su izquierda tiene coordenadas (-L,0), etc. Los 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 modificación de la fórmula,
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 el Manso Y Walton fórmula uno en este caso). ¿Esto tiene sentido? Pruebe algunos de los simples valores de a, b, y c:
a = c = 1, b = 0. Este es un paraboloide circular; su la curvatura Gaussiana debe ser positiva. El valor de 4(ac-b^2) en efecto es positivo (igual a 4).
a = c = 0, b = 1. Este es un hyperboloid de una hoja-una silla-el estándar ejemplo de una superficie de negativa de la curvatura. Bastante seguro, 4(ac-b^2) = -4.
a = 1, b = 0, c = -1. Esta es otra ecuación de la hyperboloid de una hoja (girado 45 grados). Una vez más, 4(ac-b^2) = -4.
a = 1, b = 0, c = 0. Esta es una superficie plana doblado en una forma parabólica. Ahora, 4(ac-b^2) = 0: el cero de la curvatura Gaussiana detecta correctamente la planitud de la superficie de este.
Si prueba el código de la pregunta en estos ejemplos, usted encontrará que siempre se obtiene un valor erróneo.
