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

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

2012-08-25T01:43:42Z

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

← Предыдущая версия		Версия от 03:43, 25 августа 2012
Строка 255:		Строка 255:
	Для расчетов также понадобится дополнительная библиотека для осуществления операций с матрицами. Библиотеку LU-декомпозиции для решения системы, а так же пример программы можно скачать [http://gis-lab.info/other/polynom-cs.rar здесь]. Прислал Deshchenko Sergey.		Для расчетов также понадобится дополнительная библиотека для осуществления операций с матрицами. Библиотеку LU-декомпозиции для решения системы, а так же пример программы можно скачать [http://gis-lab.info/other/polynom-cs.rar здесь]. Прислал Deshchenko Sergey.

			<syntaxhighlight lang="cpp">
	using System;		using System;

Строка 359:		Строка 359:
	}		}
	}		}
			</syntaxhighlight>

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

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

2012-08-25T01:41:32Z

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

← Предыдущая версия		Версия от 03:41, 25 августа 2012
Строка 205:		Строка 205:

	Для расчетов также понадобится дополнительная [http://gis-lab.info/other/polynom-calc-inv-delphi.zip библиотека] для осуществления операций с матрицами.		Для расчетов также понадобится дополнительная [http://gis-lab.info/other/polynom-calc-inv-delphi.zip библиотека] для осуществления операций с матрицами.
			<syntaxhighlight lang="delphi">
	map:record		map:record
	pnts:array of record // набор точек для привязки карты		pnts:array of record // набор точек для привязки карты
Строка 247:		Строка 247:
	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;		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;		end;
			</syntaxhighlight>

	[http://gis-lab.info/other/polynom-calc-mathcad.rar ]		[http://gis-lab.info/other/polynom-calc-mathcad.rar ]

Bolotoved: /* Пример пересчета в R */

2012-08-25T01:39:37Z

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

← Предыдущая версия		Версия от 03:39, 25 августа 2012
Строка 147:		Строка 147:

	Зададим исходные данные:		Зададим исходные данные:
			<syntaxhighlight lang="rsplus">
	x = c(83.786,109.929,1038.000,539.107,831.036,632.786)~~<br />~~ y = c(-36.107,-582.929,-434.786,-694.036,-352.000,-219.107)~~<br />~~ x2 = c(557124.596,564344.898,646174.994,603772.500, 626857.5, 607905.000) ~~<br />~~ y2 = c(5479746.857,5376737.207,5421503.083,5363472.000,5433468.000,5455042.500)		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)
			</syntaxhighlight>

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

			<syntaxhighlight lang="rsplus">
	mat = matrix( c(1, 1, 1, 1, 1, 1, x, y, x^2, x*y, y^2), nrow = 6, ncol = 6)		mat = matrix( c(1, 1, 1, 1, 1, 1, x, y, x^2, x*y, y^2), nrow = 6, ncol = 6)
			</syntaxhighlight>

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

	an = solve(mat, x2)~~<br />~~ bn = solve(mat, y2)		<syntaxhighlight lang="rsplus">
			an = solve(mat, x2)
			bn = solve(mat, y2)
			</syntaxhighlight>

	Или обратную задачу (x,y <- x2,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)~~<br />~~ aninv = solve(matinv, x)~~<br />~~ bninv = solve(matinv, y)		<syntaxhighlight lang="rsplus">
			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~~<br />~~ [1] 5.485053e+06 1.809379e+01 1.888970e+02 -2.225621e-04 -3.448929e-04 -6.190561e-04		<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
			</syntaxhighlight>

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

	testpoint = c(500, -300)		<syntaxhighlight lang="rsplus">
	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		testpoint = c(500, -300)
	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		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
			</syntaxhighlight>

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

	596703.3		<syntaxhighlight lang="rsplus">
	5437371.0		596703.3
			5437371.0
			</syntaxhighlight>

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

Bolotoved: /* Пример пересчета в Excel */

2012-08-25T01:32:24Z

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

← Предыдущая версия		Версия от 03:32, 25 августа 2012
Строка 103:		Строка 103:
	Для быстрой проверки можно воспользоваться [http://gis-lab.info/other/polynom-calc.xls данной таблицей в MS Excel], которая использует вышеприведенные расчеты для вычисления трансформации и предсказания новых координат. С помощью этой таблицы, используем тестовый набор данных и сравним коэффициенты преобразования с полученными в ERDAS IMAGINE:		Для быстрой проверки можно воспользоваться [http://gis-lab.info/other/polynom-calc.xls данной таблицей в MS Excel], которая использует вышеприведенные расчеты для вычисления трансформации и предсказания новых координат. С помощью этой таблицы, используем тестовый набор данных и сравним коэффициенты преобразования с полученными в ERDAS IMAGINE:

	{\| ~~width~~="~~60%~~" align="center"		{\| class="wikitable" align="center"
	\|		\|
	Коэффициент		Коэффициент

Bolotoved: /* Тестовый набор данных */

2012-08-25T01:32:02Z

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

← Предыдущая версия		Версия от 03:32, 25 августа 2012
Строка 14:		Строка 14:
	Для дальнейших пересчетов, зададим тестовый набор данных, на котором будем демонстрировать корректность вычислений. Данный пример иллюстрирует преобразование из локальной системы координат в прямоугольную (спроектированную), но на месте конечных координат могут быть и географического координаты долгота/широта.		Для дальнейших пересчетов, зададим тестовый набор данных, на котором будем демонстрировать корректность вычислений. Данный пример иллюстрирует преобразование из локальной системы координат в прямоугольную (спроектированную), но на месте конечных координат могут быть и географического координаты долгота/широта.

	{\| ~~width~~="~~60%~~" align="center"		{\| class="wikitable" align="center"
	\|		\|
	Точка		Точка
Строка 97:		Строка 97:
	Так же, используя тот же инструмент ERDAS IMAGINE, посмотрим вычисленные им коэффициенты преобразования. Примечание: по какой-то причине, ERDAS IMAGINE показыват коэффициенты не прямой, а обратной задачи, поэтому при вычислении например в Excel, нужно будет поменять местами x,y и x',y'. Итак, для нашего тестового набора данных коэффициенты расчитанные ERDAS IMAGINE следующие:		Так же, используя тот же инструмент ERDAS IMAGINE, посмотрим вычисленные им коэффициенты преобразования. Примечание: по какой-то причине, ERDAS IMAGINE показыват коэффициенты не прямой, а обратной задачи, поэтому при вычислении например в Excel, нужно будет поменять местами x,y и x',y'. Итак, для нашего тестового набора данных коэффициенты расчитанные ERDAS IMAGINE следующие:

	[[Image:polynom-calc-09.gif\|425px\|система уравнений]]		[[Image:polynom-calc-09.gif\|425px\|система уравнений\|center]]

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

Bolotoved: Новая страница: «{{Статья|Опубликована|polynom-calc-examples}} {{Аннотация|Математические выкладки решения задачи, п…»

2012-08-25T01:30:50Z

Новая страница: «{{Статья|Опубликована|polynom-calc-examples}} {{Аннотация|Математические выкладки решения задачи, п…»

Новая страница

{{Статья|Опубликована|polynom-calc-examples}}
{{Аннотация|Математические выкладки решения задачи, применяемой при привязке данных}}

При операции географической привязки данных, то есть перевода данных из локальной системы координат в географическую или прямоугольную, сам пересчет обычно происходит "за сценой" и его особенности часто понятны пользователю только интуитивно. Эта статья использует [http://gis-lab.info/qa/polynom-calc.html математические выкладки] алгоритма и показывает их реализацию, с помощью различных инструментов. Основная цель - быстрая интеграция кода в свое программное обеспечение и просто лучшее понимание процесса привязки.

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

==Математика==

Математические выкладки описывающие аналитическое решение задачи приводятся [http://gis-lab.info/qa/polynom-calc.html в отдельной статье].

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

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

{| width="60%" align="center"
|
Точка
|
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 следующие:

[[Image:polynom-calc-09.gif|425px|система уравнений]]

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

Для быстрой проверки можно воспользоваться [http://gis-lab.info/other/polynom-calc.xls данной таблицей в MS Excel], которая использует вышеприведенные расчеты для вычисления трансформации и предсказания новых координат. С помощью этой таблицы, используем тестовый набор данных и сравним коэффициенты преобразования с полученными в ERDAS IMAGINE:

{| width="60%" align="center"
|
Коэффициент
|
a
|
b
|-
|
0
| - 15071.583923079
| width="234" | -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)

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

[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-го порядков ([http://gis-lab.info/other/polynom-calc-mathcad.rar скачать]). Прислал Александр Г.

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

Для расчетов также понадобится дополнительная [http://gis-lab.info/other/polynom-calc-inv-delphi.zip библиотека] для осуществления операций с матрицами.

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;

[http://gis-lab.info/other/polynom-calc-mathcad.rar ]

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

Для расчетов также понадобится дополнительная библиотека для осуществления операций с матрицами. Библиотеку LU-декомпозиции для решения системы, а так же пример программы можно скачать [http://gis-lab.info/other/polynom-cs.rar здесь]. Прислал 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();
}
}
}

==Ссылки по теме==
*[http://gis-lab.info/qa/polynom.html Полиномиальные преобразования]
*[http://gis-lab.info/qa/rmse.html Среднеквадратичная ошибка (RMSE)]
*[http://gis-lab.info/qa/polynom-calc.html Полиномиальные преобразования - математика]