I have read the source code of several QGis-1.7.4 raster filters computing slope, aspect and curvature.
There is a formula in the filter computing the total curvature that puzzles me.
The source file is in the current version of QGis, with the following path:
qgis-1.7.4/src/analysis/raster/qgstotalcurvaturefilter.cpp
The purpose of this filter is to compute the total curvature of the surface in a nine-cells window. The function code is as follow:
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;
}
I'm ok with the "dxx" formula, and with the return expression. But I think that the "dyy" and the "dxy" formulas are inverted: this makes the total result asymmetric regarding x and y dimensions.
Am I missing something or should I replace the double derivatives expressions by:
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
Could you tell me your opinion on these formulas, if they are incorrect as I thought or if I am wrong? In this last case, do you know why the formulas have to be asymmetric regarding x and y?
Your surmises are correct. Checking for symmetry is an excellent idea: (Gaussian) curvature is an intrinsic property of a surface. Thus, rotating a grid should not change it. However, rotations introduce discretization error--except rotations by multiples of 90 degrees. Therefore, any such rotation should preserve the curvature.
We can understand what's going on by capitalizing on the very first idea of differential calculus: derivatives are limits of difference quotients. That's all we really need to know.
dxx is supposed to be a discrete approximation for the second partial derivative in the x-direction. This particular approximation (out of the many possible) is computed by sampling the surface along a horizontal transect through the cell. Locating the central cell at row 2 and column 2, written (2,2), the transect passes through the cells at (1,2), (2,2), and (3,2).
Along this transect, the first derivatives are approximated by their difference quotients, (*x32-*x22)/L and (*x22-*x12)/L where L is the (common) distance between cells (evidently equal to cellSizeAvg). The second derivatives are obtained by the difference quotients of these, yielding
dxx = ((*x32-*x22)/L - (*x22-*x12)/L)/L
= (*x32 - 2 * *x22 + *x12) / L^2.
Notice the division by L^2!
Similarly, dyy is supposed to be a discrete approximation for the second partial derivative in the y-direction. The transect is vertical, passing through cells at (2,1), (2,2), and (2,3). The formula will look the same as that for dxx but with the subscripts transposed. That would be the third formula in the question--but you still need to divide by L^2.
The mixed second partial derivative, dxy, can be estimated by taking differences two cells apart. E.g., the first derivative with respect to x at cell (2,3) (the top middle cell, not the central cell!) can be estimated by subtracting the value to its left, *x13, from the value at its right, *x33, and dividing by the distance between those cells, 2L. The first derivative with respect to x at cell (2,1) (the bottom middle cell) is estimated by (*x31 - *x11)/(2L). Their difference, divided by 2L, estimates the mixed partial, giving
dxy = ((*x33 - *x13)/(2L) - (*x31 - *x11)/(2L))/(2L)
= (*x33 - *x13 - *x31 + *x11) / (4 L^2).
I'm not really sure what is meant by "total" curvature, but it's probably intended to be the Gaussian curvature (which is the product of the principal curvatures). According to Meek & Walton 2000, Equation 2.4, the Gaussian curvature is obtained by dividing dxx * dyy - dxy^2 (notice the minus sign!--this is a determinant) by the square of the norm of the gradient of the surface. Thus, the return value quoted in the question is not quite a curvature, but it does look like a messed-up partial expression for the Gaussian curvature.
We find, then, six errors in the code, most of them critical:
dxx needs to be divided by L^2, not 1.
dyy needs to be divided by L^2, not 1.
The sign of dxy is incorrect. (This has no effect on the curvature formula, though.)
The formulas for dyy and dxy are mixed up, as you note.
A negative sign is missing from a term in the return value.
It does not actually compute a curvature, but only the numerator of a rational expression for the curvature.
As a very simple check, let's verify that the modified formula returns reasonable values for horizontal locations on quadratic surfaces. Taking such a location to be the origin of the coordinate system, and taking its elevation to be at zero height, all such surfaces have equations of the form
elevation = a*x^2 + 2b*x*y + c*y^2.
for constant a, b, and c. With the central square at coordinates (0,0), the one to its left has coordinates (-L,0), etc. The nine elevations are
*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
Whence, by the modified formula,
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.
The curvature is estimated as 2a * 2c - (2b)^2 = 4(ac - b^2). (The denominator in the Meek & Walton formula is one in this case.) Does this make sense? Try some simple values of a, b, and c:
a = c = 1, b = 0. This is a circular paraboloid; its Gaussian curvature should be positive. The value of 4(ac-b^2) indeed is positive (equal to 4).
a = c = 0, b = 1. This is a hyperboloid of one sheet--a saddle--the standard example of a surface of negative curvature. Sure enough, 4(ac-b^2) = -4.
a = 1, b = 0, c = -1. This is another equation of the hyperboloid of one sheet (rotated by 45 degrees). Once again, 4(ac-b^2) = -4.
a = 1, b = 0, c = 0. This is a flat surface folded into a parabolic shape. Now, 4(ac-b^2) = 0: the zero Gaussian curvature correctly detects the flatness of this surface.
If you try out the code in the question on these examples, you will find it always gets an erroneous value.
