Полиномиальные преобразования - примеры реализации

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/polynom-calc-examples.html

Математические выкладки решения задачи, применяемой при привязке данных

При операции географической привязки данных, то есть перевода данных из локальной системы координат в географическую или прямоугольную, сам пересчет обычно происходит "за сценой" и его особенности часто понятны пользователю только интуитивно. Эта статья использует математические выкладки алгоритма и показывает их реализацию, с помощью различных инструментов. Основная цель - быстрая интеграция кода в свое программное обеспечение и просто лучшее понимание процесса привязки.

Так как одно из наиболее часто используемых преобразований при привязке - полиномиальное преобразование 2-й степени, мы иллюстрируем наши рассчёты на его примере. Вычисления для аффинного преобразования (оно же полиномиальное преобразование 1-й степени) выполняются аналогичным образом, с меньшим количеством коэффициентов.

Математика

Математические выкладки описывающие аналитическое решение задачи приводятся в отдельной статье.

Тестовый набор данных

Для дальнейших пересчетов, зададим тестовый набор данных, на котором будем демонстрировать корректность вычислений. Данный пример иллюстрирует преобразование из локальной системы координат в прямоугольную (спроектированную), но на месте конечных координат могут быть и географического координаты долгота/широта.

Точка	x	y	x'	y'
1	83.786	-36.107	557124.596	5479746.857
2	109.929	-582.929	564344.898	5376737.207
3	1038.000	-434.786	646174.994	5421503.083
4	539.107	-694.036	603772.500	5363472.000
5	831.036	-352.000	626857.500	5433468.000
6	632.786	-219.107	607905.000	5455042.500

Сделаем нашей тестовой точкой точку с координатами: x = 500 и y = -300. Используя возможности ERDAS IMAGINE получим его вариант предсказания: 596703.345492103, 5437370.8467262

Так же, используя тот же инструмент ERDAS IMAGINE, посмотрим вычисленные им коэффициенты преобразования. Примечание: по какой-то причине, ERDAS IMAGINE показыват коэффициенты не прямой, а обратной задачи, поэтому при вычислении например в Excel, нужно будет поменять местами x,y и x',y'. Итак, для нашего тестового набора данных коэффициенты расчитанные ERDAS IMAGINE следующие:

Пример пересчета в Excel

Для быстрой проверки можно воспользоваться данной таблицей в MS Excel, которая использует вышеприведенные расчеты для вычисления трансформации и предсказания новых координат. С помощью этой таблицы, используем тестовый набор данных и сравним коэффициенты преобразования с полученными в ERDAS IMAGINE:

Коэффициент	a	b
0	- 15071.583923079	-24679.2116161184
1 (X)	0.01203458721463060	-0.0017086683749
2 (Y)	0.00267941522144888	0.0040310417433
3 (X2)	-0.00000000027645868	0.0000000001492
4 (XY)	-0.00000000011992843	0.0000000000882
5 (Y2)	-0.00000000019248837	0.0000000001062

Результат - не отличается от результата ERDAS IMAGINE.

Пример пересчета в R

Зададим исходные данные:

x = c(83.786,109.929,1038.000,539.107,831.036,632.786) 
y = c(-36.107,-582.929,-434.786,-694.036,-352.000,-219.107)
x2 = c(557124.596,564344.898,646174.994,603772.500, 626857.5, 607905.000) 
y2 = c(5479746.857,5376737.207,5421503.083,5363472.000,5433468.000,5455042.500)

Построим матрицу:

mat = matrix( c(1, 1, 1, 1, 1, 1, x, y, x^2, x*y, y^2), nrow = 6, ncol = 6)

И решим прямую задачу (x,y -> x2,y2) для нахождения коэффициентов a и b:

an = solve(mat, x2)
bn = solve(mat, y2)

Или обратную задачу (x,y <- x2,y2), для этого нам также понадобится переопределить матрицу:

matinv <- matrix( c(1, 1, 1, 1, 1, 1, x2, y2, x2^2, x2*y2, y2^2), nrow = 6, ncol = 6)
aninv = solve(matinv, x)
bninv = solve(matinv, y)

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

<syntaxhighlight lang="rsplus">
[1] 5.493158e+05 8.948879e+01 -8.564444e+00 2.322324e-04 2.727497e-04 6.293799e-04
[1] 5.485053e+06 1.809379e+01 1.888970e+02 -2.225621e-04 -3.448929e-04 -6.190561e-04

Создадим тестовую точку и проверим результат:

testpoint = c(500, -300)
xpred = an[1] + an[2]* testpoint[1] + an[3]* testpoint[2] + an[4]* testpoint[1]^2 + an[5]* testpoint[1]*testpoint[2] + an[6]*testpoint[2]^2
ypred = bn[1] + bn[2]* testpoint[1] + bn[3]* testpoint[2] + bn[4]* testpoint[1]^2 + bn[5]* testpoint[1]*testpoint[2] + bn[6]*testpoint[2]^2

Наш результат:

596703.3
5437371.0

Что отличается от результатов ERDAS лишь на доли метра. Что и требовалось доказать.

Листы рассчётов для MathCad

Файлы представляют из себя лист рассчётов для среды MathCad (версия 11 и выше), где на примере показано как можно трансформировать координаты с использованием полиномиальных преобразований 1-го и 2-го порядков (скачать). Прислал Александр Г.

Процедура для пересчета на Delphi

Для расчетов также понадобится дополнительная библиотека для осуществления операций с матрицами.

 map:record
      pnts:array of record                    // набор точек для привязки карты
      xr,yr,xg,yg:extended;                 // xr,yr - растровые координаты; xg,yg - географические координаты
      end;
      coeff_x:array of extended;
      coeff_y:array of extended;
      end;
 
 function getx(x,y:extended):extended;var m:word;
 begin
     if length(main.map.coeff_x)>0 then getx:=main.map.coeff_x[1]+main.map.coeff_x[2]*x+main.map.coeff_x[3]*y+main.map.coeff_x[4]*x*x+main.map.coeff_x[5]*x*y+main.map.coeff_x[6]*y*y else getx:=-1;
 end;
 function gety(x,y:extended):extended;
 begin
     if length(main.map.coeff_y)>0 then gety:=main.map.coeff_y[1]+main.map.coeff_y[2]*x+main.map.coeff_y[3]*y+main.map.coeff_y[4]*x*x+main.map.coeff_y[5]*x*y+main.map.coeff_y[6]*y*y else gety:=-1;
 end;
 
 procedure map_calculate_ceeff;var aa:TReal2DArray;var i,j:byte;
 begin
 // вычисление коэффициентов полинома
 //--------------------------------------------------------
 // заполнение обратной матрицы
  setlength(aa,7,7);
  for i:=1 to 6 do begin
    aa[1,i]:=1;
    aa[2,i]:=main.map.pnts[i-1].xg;
    aa[3,i]:=main.map.pnts[i-1].yg;
    aa[4,i]:=main.map.pnts[i-1].xg*main.map.pnts[i-1].xg;
    aa[5,i]:=main.map.pnts[i-1].xg*main.map.pnts[i-1].yg;
    aa[6,i]:=main.map.pnts[i-1].yg*main.map.pnts[i-1].yg;
  end;
  // вычисление обратной матрицы
  if Inverse(aa,6)=false then begin main.show_err('Ошибка. Необходим другой набор контрольных точек.');exit;end;
 
  // умножение с обратной матрицей X
  setlength(main.map.coeff_x,7);for j:=1 to 6 do main.map.coeff_x[j]:=0;
  for i:=1 to 6 do for j:=1 to 6 do main.map.coeff_x[i]:=main.map.coeff_x[i]+aa[j,i]*main.map.pnts[j-1].xr;
  // умножение с обратной матрицей Y
  setlength(main.map.coeff_y,7);for j:=1 to 6 do main.map.coeff_y[j]:=0;
  for i:=1 to 6 do for j:=1 to 6 do main.map.coeff_y[i]:=main.map.coeff_y[i]+aa[j,i]*main.map.pnts[j-1].yr;
 end;

[1]

Процедура для пересчета на C#

Для расчетов также понадобится дополнительная библиотека для осуществления операций с матрицами. Библиотеку LU-декомпозиции для решения системы, а так же пример программы можно скачать здесь. Прислал Deshchenko Sergey.

 
 using System;
 
 namespace Transformation
 {
     class Program
     {
         static void Main()
         {
             double[] x1 = { 83.786, 109.929, 1038.000, 539.107, 831.036, 632.786 };
             double[] y1 = { -36.107, -582.929, -434.786, -694.036, -352.000, -219.107 };
             double[] x2 = { 557124.596, 564344.898, 646174.994, 603772.500, 626857.500, 607905.000 };
             double[] y2 = { 5479746.857, 5376737.207, 5421503.083, 5363472.000, 5433468.000, 5455042.500 };
             int n = 6;
             int p = 6;
 
             /*************************************************************************
             Входные параметры:
             m   -   Матрица системы.
                 Массив с нумерацией элементов [1..N, 1..N].
             bForX   -   Правая часть для X.
             bForY   -   Правая часть для Y.
                 Массив с нумерацией элементов [1..N, 1..N].
             n   -   Размерность системы.
             p   -   Количество точек.
             a   -   Решение системы с X.
             b   -   Решение системы с Y.
             *************************************************************************/
 
             double[,] m = new double[n + 1, n + 1];
             double[] bForX = new double[n + 1];
             double[] bForY = new double[n + 1];
             double[] a = new double[n + 1];
             double[] b = new double[n + 1];
 
             for (int i = 0; i < p; i++)
             {
 
                 m[1, 1] = n;
                 m[1, 2] = m[2, 1] += x1[i];
                 m[1, 3] = m[3, 1] += y1[i];
                 m[1, 4] = m[4, 1] = m[2, 2] += Math.Pow(x1[i], 2);
                 m[1, 5] = m[5, 1] = m[2, 3] = m[3, 2] += x1[i] * y1[i];
                 m[1, 6] = m[6, 1] = m[3, 3] += Math.Pow(y1[i], 2);
                 m[2, 4] = m[4, 2] += Math.Pow(x1[i], 3);
                 m[2, 5] = m[5, 2] = m[3, 4] = m[4, 3] += Math.Pow(x1[i], 2) * y1[i];
                 m[2, 6] = m[6, 2] = m[3, 5] = m[5, 3] += x1[i] * Math.Pow(y1[i], 2);
                 m[3, 6] = m[6, 3] += Math.Pow(y1[i], 3);
                 m[4, 4] += Math.Pow(x1[i], 4);
                 m[4, 5] = m[5, 4] += Math.Pow(x1[i], 3) * y1[i];
                 m[4, 6] = m[6, 4] = m[5, 5] += Math.Pow(x1[i], 2) * Math.Pow(y1[i], 2);
                 m[5, 6] = m[6, 5] += x1[i] * Math.Pow(y1[i], 3);
                 m[6, 6] += Math.Pow(y1[i], 4);
 
                 bForX[1] += x2[i];
                 bForX[2] += x1[i] * x2[i];
                 bForX[3] += y1[i] * x2[i];
                 bForX[4] += Math.Pow(x1[i], 2) * x2[i];
                 bForX[5] += x1[i] * y1[i] * x2[i];
                 bForX[6] += Math.Pow(y1[i], 2) * x2[i];
                 bForY[1] += y2[i];
                 bForY[2] += x1[i] * y2[i];
                 bForY[3] += y1[i] * y2[i];
                 bForY[4] += Math.Pow(x1[i], 2) * y2[i];
                 bForY[5] += x1[i] * y1[i] * y2[i];
                 bForY[6] += Math.Pow(y1[i], 2) * y2[i];
             }
             {
                 System.Console.Write("Error! Degenerate matrix A!");
                 System.Console.WriteLine();
                 System.Console.ReadKey();
                 return;
             }
             if (!linsolve.solvesystem(m, bForY, n, ref b))
             {
                 System.Console.Write("Error! Degenerate matrix A!");
                 System.Console.WriteLine();
                 System.Console.ReadKey();
                 return;
             }
             double testPointX1 = 500;
             double testPointY1 = -300;
             double testPointX2 = a[1] + a[2] * testPointX1 + a[3] * testPointY1 + a[4] * Math.Pow(testPointX1, 2) + a[5] * testPointX1 * testPointY1 + a[6] * Math.Pow(testPointY1, 2);
             double testPointY2 = b[1] + b[2] * testPointX1 + b[3] * testPointY1 + b[4] * Math.Pow(testPointX1, 2) + b[5] * testPointX1 * testPointY1 + b[6] * Math.Pow(testPointY1, 2);
             System.Console.Write("X = " + testPointX2);
             System.Console.WriteLine();
             System.Console.Write("Y = " + testPointY2);
             System.Console.WriteLine();
             for (int i = 0; i < n; i++)
             {
                 System.Console.WriteLine();
                 System.Console.Write("a[" + i + "] = " + a[i + 1]);
             }
             System.Console.WriteLine();
             for (int i = 0; i < n; i++)
             {
                 System.Console.WriteLine();
                 System.Console.Write("b[" + i + "] = " + b[i + 1]);
 
             }
             System.Console.ReadKey();
         }
     }
 }

Ссылки по теме

• Полиномиальные преобразования
• Среднеквадратичная ошибка (RMSE)
• Полиномиальные преобразования - математика