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

[updates-bot]

Обсуждение статьи "Полиномиальные преобразования - вычисления и практика"

http://gis-lab.info/qa/polynom-calc.html

[Максим Дубинин]

для убедительности добавлены сравнения коэффициентов расчитанных ERDAS и приводимым Excel-калькулятором.

24.02.2008
- добавлена реализация на delphi и листы рассчётов для MathCad.
- примеры реализации вынесены в отдельную статью:
http://gis-lab.info/qa/polynom-calc-examples.html

7.11.2008
- добавлена реализация на C#

[Desh]

Спасибо за статью, все подробно и ясно.
Однако, мне необходим случай когда количество точек привязки может быть любым (естественно больше минимального).
Написано, что нужный мне материал планируется в будущем, но времени в обрез.
Может у кого есть какие идеи...
Заранее спасибо.

[KolesovDmitry]

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

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

[Desh]

Я сделал немного иначе:

Рассматривается полиномиальное преобразование 2-й степени

у нас есть n пар уравнений
x' = a0+a1x+a2y+a3x^2+a4xy+a5y^2 + E1
y' = b0+b1x+b2y+b3x^2+b4xy+b5y^2 + E1
.................................................................

дальше нужно решить систему для x' (аналогично для y'), где
min(P) = E1^2 + E2^2 + E3^2 + ... + En^2 - минимально
дифференцируем и получаем коэффициенты A0, A1, A2, A3, A4, A5
а из уравнения для y' B0, B1, B2, B3, B4, B5

[UnderFelixAbove]

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

Имеется 10-30 точек - привязки, которые характеризуются пиксельными(местоположение на отсканированной карте) и географическими(долгота, широта) координатами. Эти точки предпологается использовать для привязки карты. Вопрос только в том, какие точки лучше для этого использовать. Предпологается использовать либо привязку по 3 точкам(афинные преобразования), либо по 6 (Полиномиальные преобразования). Максимальная из метрик карты местности до 400 км.

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

Если точка, положение которой мы хотим найти, зная ее географические координаты, лежит внутри области образованной точками-привязки, то погрешность подсчетов оказывается несколько меньше, если бы точка лежали вне этой области. При этом интересен факт, чем меньше расстояние от искомой точки до границ этой области, тем и погрешность меньше. Следовательно первоночальная идея о расположение точек-привязки по периметру была не совсем удачной(поправьте, пожалуйста, если я не прав).

Если сравнивать привязку по 3 точкам и по 6, то можно сделать общий вывод, что привязка по 6 точкам дает практически всегда небольшую и допустимую погрешность, чего нельзя сказать о привязке по 3 точкам. Конечно есть случаи, когда привязка по 3 точкам дает лучшие показания, чем по 6, но это наблюдается только тогда, когда область, образованная тремя точками-привязки образует треугольник по виду похожему на равносторонний, и искомая точка лежит внутри этой треугольной области. Понимаю, что ничего нового не сказал, это просто для полноты картины.

Из всего вышесказано можно сделать общий вывод:
1. Для более точного определения положения точки на карте необходимо выбирать точки-привязки, так, чтобы искомая точка была наиболее близка к ним приближена и по возможности лежала внутри области.
2. Привязку по 3 точкам выгоднее использовать, если искомая точка попала внутрь треугольной области.
3. В общем случае нужно использовать привязку по 6 точкам и стараться, чтобы искомая точка лежала внутри области, если этого сделать не получается, то стараться брать ближайшие 6 точек.

из п.3. сразу возникает нюанс. Если выбранные 6 ближайших контрольных точек будут лежать на одной прямой! Чтобы это избежать появляется идея заранее разбить наш базовый набор контрольных точек на группы по 6 точек, которые образуют внутри группы выпуклые оболочки. И уже, имея этот набор выпуклых оболочек, искать наиболее приближенную оболочку. Подход с выпуклыми оболочками исключит случай, когда контрольные точки лежат на одной прямой.

Поделитесь, пожалуйста, своими соображениями на этот счет!

P.S: Использование привязки более, чем по 6 точкам, например при порядке преобразования = 3(минимальное кол-во контрольных точек = 10) - является целесообразным для решения данной задачи? Если да, то подскажите где поискать математику по данному вопросу.
Заранее спасибо.

[KolesovDmitry]

...
К сожалению никогда не занимался привязкой карт, поэтому очень хочу услышать мнения людей, авторитетных в этом вопросе.

Я тоже не специалист, но поскольку спецы молчат, то пропробую прокоментировать

Имеется 10-30 точек - привязки, которые характеризуются пиксельными(местоположение на отсканированной карте) и географическими(долгота, широта) координатами. Эти точки предпологается использовать для привязки карты. Вопрос только в том, какие точки лучше для этого использовать. Предпологается использовать либо привязку по 3 точкам(афинные преобразования), либо по 6 (Полиномиальные преобразования).

Не думаю, что скажу вам что-то новое, но на всякий случай хочу все-таки расставить точки над i:

• -Действительно, линейное (аффинное) преобразование выполняется по трем точкам.
  -Но! Никто ведь не мешает взять контрольных точек не три а гораздо больше (обычно так и делается), а потом рассчитать такие параметры аффинного преобразования, чтобы погрешность привязки (в контрольных точках) была минимальной из всех возможных аффинных преобразований. Расчитать соответствующие параметры можно, к примеру, методом наименьших квадратов.
  -Соответственно, и для преобразования по полиному второго порядка можно брать не 6 точек, а хоть миллион, и рассчитать соответственно необходимые оптимальные коэффициенты преобразования тем же методом наименьших квадратов.

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

Не знаю, как у других, а у меня возник вопрос: почему? Как вы пришли к такому выводу -- расскажите немного подробнее.

Из всего вышесказано можно сделать общий вывод:
1. Для более точного определения положения точки на карте необходимо выбирать точки-привязки, так, чтобы искомая точка была наиболее близка к ним приближена и по возможности лежала внутри области.
2. Привязку по 3 точкам выгоднее использовать, если искомая точка попала внутрь треугольной области.
3. В общем случае нужно использовать привязку по 6 точкам и стараться, чтобы искомая точка лежала внутри области, если этого сделать не получается, то стараться брать ближайшие 6 точек.

Замечание по пункту 1: в практических задачах чаще всего заранее не известно, какая точка будет "искомой", поэтому карту привязывают так, чтобы были минимальны усредненные погрешности привязки (усредненные по всей карте). А для этого, при расчете параметров преобразования берут точек больше, чем требуется для решения системы уравнений (т.е. берут, к примеру, не 3 точки, а 3+n). Ну и используют, соответственно, не алгебраические методы расчетов параметров преобразований, а статистические.
Замечание по п.п. 2-3: эти пункты следуют из ваших предыдущих утверждений, но, прежде чем им следовать, хотелось бы посмотреть, как они получаются.

[PAV74]

Уважаемый KolesovDmitry.
Так совпало, что меня тоже интересуют вопросы заданные UnderFelixAbove, и он предвосхитил мои вопросы, которые я хотел задать в этом треде. Я новичок в этом деле, и первое что я попытался реализовать было почти тоже самое что написал UnderFelixAbove. Для начала я использовал как бы плавающее преобразование полиномом 2ой степени. Т.е. для преобразования координат, использовал кластер из 6 опорных точек, и по мере того как координаты менялись я добавлял новую ближнюю точку в него и удалял самую отдаленную. Таким образом я думал избавиться от скачков при изменении базиса преобразования. Но как показал эксперимент, на больших картах, это уже не работает, точнее погрешность вносимая изменением одной точки из 6 уже неприемлема для меня. Но это лирика. Пришлось копать в сторону регрессии методом МНК, но все мои поиски насчет того как это реализовать пока не увенчались успехом. Может у Вас есть какие нибудь ссылки по данной теме которыми Вы можете поделиться. Сейчас у меня по регрессии лишь понимание того, что по n опорным точкам должно получиться n уравнений плюс условие наименьших квадратов. Но вопрос в том как эту систему решить статистическим методом.

[KolesovDmitry]

...
Пришлось копать в сторону регрессии методом МНК, но все мои поиски насчет того как это реализовать пока не увенчались успехом. Может у Вас есть какие нибудь ссылки по данной теме которыми Вы можете поделиться. Сейчас у меня по регрессии лишь понимание того, что по n опорным точкам должно получиться n уравнений плюс условие наименьших квадратов. Но вопрос в том как эту систему решить статистическим методом.

Каких-то особых ссылок у меня нет, могу сказать лишь, что этот метод рассматривается в любом учебнике по мат.статистике. Могу дать на вскидку несколько таких ссылок:

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

Одна из первых попавшихся ссылок -- здесь подробнее, но тоже без лишней зауми. После этого по крайней мере будет понятно, откуда берется формула, приведенная в Википедии.

Если по пониманию или реализации метода будут вопросы -- спрашивайте, не стесняйтесь, если смогу -- постараюсь ответить.

[Максим Дубинин]

если вам, как мне, более привычна матричная форма, то вот тут я попытался разложить, что к чему, если число точек > минимально необходимого, т.о. если вы занимаетесь реализацией этого алгоритма, ваша программа должна уметь транспонировать, инвертировать и умножать матрицы
http://gis-lab.info/qa/polynom-calc-many.html

реализовать это с помощью R - очень просто, если нужно, могу подробно описать.

[PAV74]

Коллеги, спасибо за ссылки.
По этой проблеме все прояснилось. Алгоритм реализовал по статье Sim-а (http://gis-lab.info/qa/polynom-calc-many.html). В итоге получилось даже проще чем я предполагал. Ниже реализация на Delphi может кому еще сгодится.

Код: Выделить всё

Type
    // точка привязки
    TControlPoint = record
      E: Extended; // Easting -  X координата исходного пространства
      N: Extended; // Northing - Y координата исходного пространства
      X: Extended; // X координата целевого пространства (в которое осуществляется преобразование)
      Y: Extended; // Y координата целевого пространства 
    end;

    TControlPoints = array of TControlPoint;

    // параметры преобразования
    TTransformParams = record
      X: array of Extended;
      Y: array of Extended;
    end;

procedure Polynom2RMS_PrepareParams(const ControlPoints: TControlPoints; out Params: TTransformParams);
{GPS-Lab.info}
var
  PointsNumber: integer;
  XMT,XM,XMM,XMR: TReal2DArray;
  I,J:byte;
begin
PointsNumber:=Length(ControlPoints);
if (PointsNumber < 6) then Raise Exception.Create('мало точек привязки'); //. =>
{вычисление коэффициентов полинома преобразования E,N => X,Y}
//. заполнение матрицы XM
SetLength(XM,PointsNumber+1,7);
for I:=1 to PointsNumber do begin
  XM[I,1]:=1;
  XM[I,2]:=ControlPoints[I-1].E;
  XM[I,3]:=ControlPoints[I-1].N;
  XM[I,4]:=sqr(ControlPoints[I-1].E);
  XM[I,5]:=ControlPoints[I-1].E*ControlPoints[I-1].N;
  XM[I,6]:=sqr(ControlPoints[I-1].N);
  end;
//. заполнение транспонированной матрицы XMT
SetLength(XMT,7,PointsNumber+1);
for I:=1 to PointsNumber do
  for J:=1 to 6 do XMT[J,I]:=XM[I,J];
//. вычисление произведения матриц XMT и XM -> результат матрица XMM размером (6 x 6)
MultiplyMatrixes(6,PointsNumber,6, XMT,XM,{out} XMM);
//. вычисление обратной матрицы XMM
if (NOT Inverse(XMM,6)) then Raise Exception.Create('недействительные точки привязки'); //. =>
//. вычисление произведения матриц XMM и XMT -> результат матрица XMR размером (6 x PointsNumber)
MultiplyMatrixes(6,6,PointsNumber, XMM,XMT,{out} XMR);
//. умножение вектора с результирующей матрицей XMR
SetLength(Params.X,7); for I:=1 to 6 do Params.X[I]:=0;
for I:=1 to 6 do
  for J:=1 to PointsNumber do Params.X[I]:=Params.X[I]+XMR[I,J]*ControlPoints[J-1].X;
//. умножение вектора с результирующей матрицей XMR
SetLength(Params.Y,7); for I:=1 to 6 do Params.Y[I]:=0;
for I:=1 to 6 do
  for J:=1 to PointsNumber do Params.Y[I]:=Params.Y[I]+XMR[I,J]*ControlPoints[J-1].Y;
end;


procedure Polynom2RMS_Transform(const Params: TTransformParams; const E,N: Extended; out X,Y: Extended);
begin
X:=Params.X[1]+Params.X[2]*E+Params.X[3]*N+Params.X[4]*sqr(E)+Params.X[5]*E*N+Params.X[6]*sqr(N);
Y:=Params.Y[1]+Params.Y[2]*E+Params.Y[3]*N+Params.Y[4]*sqr(E)+Params.Y[5]*E*N+Params.Y[6]*sqr(N);
end;

Используется библиотека для работы с матрицами как в примере к статье "Полиномиальные преобразования - примеры реализации". Дополнительно к ней надо скачать модуль для произведения матриц из того же пакета ALGLIB (http://alglib.sources.ru/translator/dl/ ... delphi.zip).

Вопрос к Sim: Применительно к моей задаче нужно сделать преобразование большого участка карты, но с возможностью повышения точности преобразования в одной конкретной области, т.е. в окрестности какой либо точки привязки. Сейчас я реализовал это путем добавления нескольких экземпляров одной и той-же точки в набор точек преобразования. Таким образом я как бы повышаю вес данной привязки односительно других в наборе с использованием Вашей формулы. Это конечно тоже вариант, но при добавлении дубликата растет размерность матрицы. А нельзя ли "малой кровью" путем добавления вектора весовых коэффициентов точек привязки видоизменить формулу так, чтобы она считала по методу взвешенных наименьших квадратов. И еще маленкая ремарка там в статье неработает ссылка на примеры и опечатка в формулах, крайний справа вектор y1,y2,...(x1,x2,...) имеет размерность 1..6, наверно должно быть 1..n

[Максим Дубинин]

Про взвешенный МНК мне сложно сказать, но интересно, надо будет попробовать с этим поразбираться.

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

С другой стороны, мне кажется это от лукавого, если нужно точно привязать некоторый фрагмент карты, то лучше его вырезать и привязать. Ну а большую карту привязать отдельно-паралельно.

Спасибо за опечатки, поправил, заодно еще пяток нашел.

procedure Polynom2RMS_PrepareParams(const ControlPoints: TControlPoints; out Params: TTransformParams);
{GPS-Lab.info}

И у вас тоже опечаточка, в названии сайта

[UnderFelixAbove]

Извиняюсь за долгое отсутствие, уже все разобрались и все выяснилось.
KolesovDmitry спрашивал откуда берутся умозаключения по поводу

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

Это просто наблюдения.. математически никак не подкрепленные.

[ViTr]

Здравствуйте! Сразу оговорюсь, я пока плохо разбираюсь во всем этом. Переделал по аффинные преобразования, т.е убрал лишние коэффициенты.
Все не могу до конца разобраться в нюансах кода под делфи(. Как например передать координаты своих точек для привязки в процедуру map_calculate_ceeff;? И не могу разобраться как вобще происходит заполнение матрицы aa?
Если можете объясните вот эти строчки:
for i:=1 to 3 do
begin
aa[1,i]:=1;
aa[2,i]:=main.map.pnts[i-1].xg;
aa[3,i]:=main.map.pnts[i-1].yg;

[PAV74]

функция вычисления коэффициентов работает с глобальным объектом поэтому параметры в нее не передаются. Лучше передавать точки привязки через параметры процедуры, конечно, и на выходе иметь матрицу преобразования как результат. По поводу кода

>> Если можете объясните вот эти строчки:
for i:=1 to 3 do
begin
aa[1,i]:=1;
aa[2,i]:=main.map.pnts[i-1].xg;
aa[3,i]:=main.map.pnts[i-1].yg;

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