GIS-LAB

Географические информационные системы и дистанционное зондирование

Полиномиальные преобразования - примеры реализации для большого количества точек

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

Обсудить в форуме Комментариев — 16

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

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

Оглавление

  1. Математика
  2. Тестовый набор данных
  3. Пример пересчета в Excel
  4. Пример пересчета в R
  5. Листы рассчётов для MathCad
  6. Процедура для пересчета на Delphi
  7. Процедура для пересчета на C#

1Математика

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

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

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

Точка

x

y

x'

y'

1

84

-36

557125

5479747

2

110

-583

564345

5376737

3

1038

-435

646175

5421503

4

539

-694

603773

5363472

5

831

-352

626858

5433468

6

633

-219

607905

5455043

7

500

-300

595703

5437671

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

система уравнений

3Пример пересчета в 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.

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

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

x = c(84,110,1038,539,831,633,500)
y = c(-36,-583,-435,-694,-352,-219,-300)
x2 = c(557125,564345,646175,603773, 626858, 607905, 595703)
y2 = c(5479747,5376737,5421503,5363472,5433468,5455043, 5437671)

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

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

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

an  = solve(t(mat)%*%mat)%*%t(mat)%*%x2
bn  = solve(t(mat)%*%mat)%*%t(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)

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

[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 лишь на доли метра. Что и требовалось доказать.

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

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

6Процедура для пересчета на 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;

7Процедура для пересчета на 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();
        }
    }
}

 

Обсудить в форуме Комментариев — 16

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

Последнее обновление: September 09 2021

Donate

Новое на форумеRSS

© GIS-Lab и авторы, 2002-2021. При использовании материалов сайта, ссылка на GIS-Lab и авторов обязательна. Содержание материалов - ответственность авторов. (подробнее).