Понимание фильтра кривизны растрового анализа QGIS?

Я прочитал исходный код нескольких растровых фильтров QGis-1.7.4, вычисляющих наклон, аспект и кривизну.

В фильтре есть формула для вычисления полной кривизны, которая меня озадачивает.

Исходный файл находится в текущей версии QGis со следующим путем:

QGIS-1.7.4 / SRC / анализ / растровой / qgstotalcurvaturefilter.cpp

Целью этого фильтра является вычисление полной кривизны поверхности в окне из девяти ячеек. Код функции выглядит следующим образом:

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

Я в порядке с формулой "dxx" и с выражением возврата. Но я думаю, что формулы «dyy» и «dxy» инвертированы: это делает общий результат асимметричным относительно размеров x и y.

Я что-то упустил или я должен заменить выражения двойных производных на:

  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

Не могли бы вы высказать свое мнение об этих формулах, если они неверны, как я думал, или я ошибаюсь? В последнем случае, знаете ли вы, почему формулы должны быть асимметричными относительно x и y?

Ваши предположения верны. Проверка на симметрию - отличная идея: кривизна (гауссова) является внутренним свойством поверхности. Таким образом, вращение сетки не должно ее менять. Однако при поворотах возникает ошибка дискретизации - за исключением поворотов, кратных 90 градусам. Следовательно, любое такое вращение должно сохранять кривизну.

Мы можем понять, что происходит , извлекая выгоду из самой первой идеи дифференциального исчисления: производные - это пределы разностей. Это все, что нам действительно нужно знать.

dxxпредполагается, что это дискретное приближение для второй частной производной в направлении х. Это конкретное приближение (из множества возможных) вычисляется путем выборки поверхности вдоль горизонтального разреза через ячейку. Расположив центральную ячейку в строке 2 и столбце 2, написанную (2,2), трансект проходит через ячейки в точках (1,2), (2,2) и (3,2).

Вдоль этого разреза первые производные аппроксимируются их разностными коэффициентами, (* x32- * x22) / L и (* x22- * x12) / L, где L - (общее) расстояние между клетками (очевидно, равное cellSizeAvg). Вторые производные получаются по разности их, что дает

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

Обратите внимание на деление на L ^ 2!

Точно так же, dyyкак предполагается , будет дискретным приближением для второй частной производной в у-направлении. Трансект является вертикальным, проходя через клетки в (2,1), (2,2) и (2,3). Формула будет выглядеть так же, как для, dxxно с транспонированными индексами. Это будет третья формула в вопросе - но вам все равно нужно разделить на L ^ 2.

Смешанную вторую частную производную dxyможно оценить, взяв различия между двумя ячейками. Например, первая производная по x в ячейке (2,3) (верхняя средняя ячейка, а не центральная ячейка!) Может быть оценена путем вычитания значения слева * x13 из значения справа * х33, и деление на расстояние между этими ячейками, 2л. Первая производная по x в ячейке (2,1) (нижняя средняя ячейка) оценивается как (* x31 - * x11) / (2L). Их разница, разделенная на 2 л, оценивает смешанную частичную, давая

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

Я не совсем уверен, что подразумевается под «общей» кривизной, но, вероятно, она предназначена для гауссовой кривизны (которая является продуктом основных кривизн). Согласно Meek & Walton 2000 , уравнение 2.4, кривизна Гаусса получается делением dxx * dyy - dxy ^ 2 (обратите внимание на знак минус! - это определитель ) на квадрат нормы градиента поверхности. Таким образом, возвращаемое значение, указанное в вопросе, не совсем кривизна, но оно выглядит как испорченное частичное выражение для гауссовой кривизны.

Итак, мы находим в коде шесть ошибок , большинство из которых критические:

1. dxx нужно разделить на L ^ 2, а не 1.

2. dyy нужно разделить на L ^ 2, а не 1.

3. Знак dxy неверен. (Это никак не влияет на формулу кривизны.)

4. Как вы заметили, формулы для dyy и dxy перепутаны.

5. Отрицательный знак отсутствует в термине возвращаемого значения.

6. На самом деле он вычисляет не кривизну, а только числитель рационального выражения для кривизны.

В качестве очень простой проверки давайте проверим, что измененная формула возвращает разумные значения для горизонтальных положений на квадратичных поверхностях. Принимая такое местоположение как начало системы координат, и принимая его высоту равной нулю, все такие поверхности имеют уравнения вида

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

для константы а, б и в. С центральным квадратом в координатах (0,0), один слева имеет координаты (-L, 0) и т. Д.

*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

Откуда по модифицированной формуле

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.

Кривизна оценивается как 2a * 2c - (2b) ^ 2 = 4 (ac - b ^ 2). (Знаменатель в формуле Мика и Уолтона в этом случае равен единице.) Имеет ли это смысл? Попробуйте несколько простых значений a, b и c:

• a = c = 1, b = 0. Это круглый параболоид; его гауссова кривизна должна быть положительной. Значение 4 (ac-b ^ 2) действительно является положительным (равно 4).

• a = c = 0, b = 1. Это гиперболоид одного листа - седло - стандартный пример поверхности отрицательной кривизны. Конечно же, 4 (ac-b ^ 2) = -4.

• а = 1, б = 0, с = -1. Это еще одно уравнение гиперболоида одного листа (повернутого на 45 градусов). Еще раз, 4 (ac-b ^ 2) = -4.

• a = 1, b = 0, c = 0. Это плоская поверхность, сложенная в параболическую форму. Теперь 4 (ac-b ^ 2) = 0: нулевая гауссова кривизна правильно определяет плоскостность этой поверхности.

Если вы попробуете код в вопросе на этих примерах, вы обнаружите, что он всегда получает ошибочное значение.