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:


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?


1 Answer 1


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:

  1. dxx needs to be divided by L^2, not 1.

  2. dyy needs to be divided by L^2, not 1.

  3. The sign of dxy is incorrect. (This has no effect on the curvature formula, though.)

  4. The formulas for dyy and dxy are mixed up, as you note.

  5. A negative sign is missing from a term in the return value.

  6. 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.

  • Its always interesting to read your explicit elaborations in the morning.
    – Tomek
    Jun 11, 2012 at 5:57
  • @Tomek Now there is a diplomatic (=tactful and highly ambiguous) comment! :-)
    – whuber
    Jun 11, 2012 at 17:13
  • 1
    Thank you very much for such a complete answer! I'm going to report the formula errors since I'm now ensured that there is something to report. :)
    – Tapadi
    Jun 18, 2012 at 13:45
  • @whuber : I can confirm Tomek's reply that it is always interesting to read your comments on this forum, and I always learn something new from them!! Thank you for sharing your invaluable knowledge for free with us!! Would you mind if I ask just one more question: In any GIS application, when the curvature analysis of the terrain (raster) is performed, it is always the Gaussian curvature? Never the Mean curvature?
    – marco
    Sep 3, 2016 at 23:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.