Программная технология определения пространственых перемещений объектов по данным дистанционного зондирования земли Текст научной статьи по специальности «Математика»

УДК 528.852

А.А. Бучнев, В.П. Пяткин

ИВМиМГ СО РАН, Новосибирск

ПРОГРАММНАЯ ТЕХНОЛОГИЯ ОПРЕДЕЛЕНИЯ ПРОСТРАНСТВЕНЫХ ПЕРЕМЕЩЕНИЙ ОБЪЕКТОВ ПО ДАННЫМ ДИСТАНЦИОННОГО ЗОНДИРОВАНИЯ ЗЕМЛИ

А.А. Buchnev, V.P. Pyatkin ICM&MG SB RAS

prospect Akademika Lavrentjeva, 6, Novosibirsk, 630090, Russia

SOFTWARE DETERMINATION OF OBJECTS SPATIAL MOTION BY THE EARTH REMOTE SENSING DATA

The method based on coefficient maximums discovering of cross correlation is one of the ways of object spatial motion detection on satellite images occurring at different time. Here, the correlation is used as the means of equivalents discovering of object target which is the part of the current image, on the next image of images series.

In the work presented, to increase the results authenticity of discovering new positions, the targets are transformed, i.e. adjusted to scale, rotation and translation.

The problem solution consists of the following steps: 1) square objects-targets searching on the current image based on the results of contrast maximum or entropy maximum; 2) the searching of the positions of obtained targets on the next image, based on the achievement of maximal value of cross correlation coefficient; 3) building of vector fields of spatial objects motions in accordance with obtained positions.

While scanning of the field of search the target is transformed, i.e. scaled and rotated. Then the transformed target is successively carried to each admissible position of the searching field on the next image.

For multispectral images, the procedure of searching of targets and their equivalents is applied to data of each spectral range.

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

Согласно [1] одним из способов определения пространственных перемещений объектов по разновременным спутниковым изображениям является метод, основанный на нахождении максимумов коэффициента взаимной корреляции. В [2] аналогичный подход рассматривается в качестве метода распознавания образов, известного как корреляционное сопоставление. В обоих случаях корреляция используется как средство поиска эквивалентов объекта-эталона, представленного в виде изображения w(x, у) размерами j х к , на изображении f(x, у) размерами м х N ; предполагается, что J < м и к < N. Коэффициент взаимной корреляции

X X [f(Х + ^ У + 1) “ fm (X, У )][^) - ^ ]

К X, У ) = -^---------------------------------. (1)

°м,° Г

Здесь wm - среднее значение пикселов в эталоне w, /т - среднее значение элементов изображения / в области, покрываемой эталоном. Знаменатель в (1) является произведением стандартного отклонения пикселов эталона w на

стандартное отклонение аг пикселов изображения / в области, покрываемой эталоном.

В общем случае размеры и ориентация объекта-эталона, найденного на очередном изображении серии изображений, на последующем изображении могут иметь другие значения. Нормированное относительно изменений амплитуд выражение (1) не является таковым относительно поворота или изменений размеров (масштабирования) эталона. В [2] отмечается, что «Нормировка относительно размеров связана с пространственным масштабированием, что само по себе связано с весьма трудоемкими вычислениями. Нормировка относительно поворота является еще более трудной задачей...». Иногда для решения задачи сопоставления (совмещения) изображений с учетом масштабирования и поворота эталона используются менее точные по сравнению с (1), но более простые с вычислительной точки зрения критерии (например, в работе [3] и ряде других работ используется сумма абсолютных значений разностей соответствующих компонент изображений).

Согласно [4] поиск позиций найденных эталонов (определение смещений) на следующем изображении серии может быть реализован одним из трех методов: определением максимума коэффициента взаимной

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

В представляемой работе определение смещений эталонов производится на основе определения максимума коэффициента взаимной корреляции в пространственной области в соответствии с формулой (1). При этом эталон может подвергаться преобразованию, состоящему из масштабирования, поворота и переноса. Ниже приведен соответствующий алгоритм и представлены результаты вычислительных экспериментов на изображениях, полученных с КА «METEOSAT-8». Эти результаты свидетельствуют как о необходимости учета масштабирования и поворота эталона, так и о приемлемом времени соответствующих вычислений.

Решение задачи состоит из следующих основных шагов:

поиск на текущем изображении Ж квадратных объектов-эталонов, основанный на достижении максимума контраста либо максимума энтропии; центр эталона совпадает с центром квадрата;

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

построение векторных полей пространственных перемещений объектов в соответствии с найденными позициями.

Поиск эталонов основан на методике, предложенной EUMETSAT [4]. Согласно этой методике предусматриваются эталоны двух типов: главные и вторичные. Позиции главных эталонов совпадают с концами векторов перемещений эталонов из предыдущего изображения (для первого изображения серии главных эталонов нет).

Поиск вторичных эталонов размера TargSize производится в узлах сетки с размером GridSize. Размер квадратной области с центром в узлах сетки для поиска эталонов задается параметром TargSearch. Допустимое минимальное расстояние между эталонами управляется параметром Targ Dist. Оптимальной позицией для эталона внутри области поиска считается та, в которой достигается максимум значения управляющего параметра Par - контраста либо энтропии. При поиске эталонов внутри области размера TargSearch используются локальные средние и стандартные отклонения, вычисленные по окрестности 3*3. Контраст определяется как разность между максимальным и минимальным значениями локальных средних. Кроме перечисленных параметров, определяющих «физические» характеристики, при поиске эталонов используются параметры, характеризующие «изменчивость» изображения внутри области, покрываемой эталоном: MinStDev - минимальное значение локального стандартного отклонения, NumGrSD - минимальное количество пикселов со стандартным отклонением большим MinStDev.

Для многоспектральных изображений процедура поиска эталонов применяется к данным каждого спектрального диапазона.

Смещение определяется для каждого из KTargs найденных эталонов. Поиск новой позиции эталона производится внутри квадратной области размера Search Size. Центр области поиска совпадает с исходной позицией эталона. Новой позицией эталона считается позиция, в которой достигается максимальное значение коэффициента взаимной корреляции Corr. В процессе сканирования области поиска эталон подвергается преобразованиям масштабирования и поворота. Масштаб и угол поворота эталона определяются следующими параметрами: ScaleMin - минимальное значение масштаба эталона; ScaleMax - максимальное значение масштаба эталона; ScaleDelta - приращение масштаба эталона; AngleBeg - начальный угол поворота эталона в градусах; Angle End - конечный угол поворота эталона в градусах; Angle Delta - приращение угла поворота эталона в градусах.

Если s, Scale _ Min < s < Scale _ Max - текущий масштаб эталона, p, Angle_Beg <ф< Angle_End - текущий угол поворота эталона, то при применении преобразований масштабирования и поворота координаты (x,y)

пикселов эталона переходят в координаты (х ,у ) выходного образа, равные (х , у) = (л, у)Т, где Т - матрица преобразования:

T = s ■

(2)

Sinp cos p)

Известно (см., например, [1]), что применение прямого преобразования с матрицей (2) к объектам с дискретными координатами приводит к тому, что пиксел с координатами (x',у) может лежать между пикселами выходного образа. Использование техники назначения результату преобразования ближайшего соседа пиксела с координатами (х', у ) не является решением вопроса из-за возможного появления пустот в выходном образе. Простейшим (в идейном смысле) методом является распределение значения входного пиксела между несколькими выходными пикселами. В этом методе пикселы рассматриваются как квадраты, и часть площади входного пиксела, которая покрывает выходной пиксел, рассматривается в качестве весового коэффициента. Ясно, что в реализации техника подобного рода является достаточно громоздкой. Более простой альтернативой является использование обратного преобразования: если из каких-то дополнительных соображений известны координаты (x , у ) пикселов выходного образа, то соответствующие им координаты (x,y) пикселов эталона равны (х , у ) =(х , y)T~l.

Такая схема позволяет избежать появления пустот в выходном образе. Но здесь, также как и в прямом преобразовании, пиксел эталона, координаты которого получены с использованием матрицы T -1, может лежать между других пикселов. Следовательно, возникает проблема интерполяции пикселов эталона. Вопросы использования интерполяции будут описаны ниже.

Приведем более подробное изложение реализованной вычислительной схемы. Применяя к координатам вершин квадратного эталона с центром в начале координат прямое преобразование (2) для текущих значений s и p, получим координаты (округленные до ближайшего соседа) преобразованного квадрата. Далее к полученному квадрату применяется процедура, аналогичная методу сканирующих строк в заполнении многоугольников [5]. Первым шагом является сортировка вершин квадрата по возрастанию координаты у. Результатом ее является набор координат

(Х^ уmin X (Х2 , у2 X (Х3 , у3 X (Х4 , у,тах ) , где ут» < у2 < уз < у,тах . Таким обрa3Oм,

преобразованный квадрат занимает r = у^ - утп +1 строк. Для каждой сканирующей строки у ,1 < s < r , определяются (округленные до ближайших целых) координаты точек ее пересечения с ребрами квадрата (xbs, у),(xl, у s) . В строке у преобразованный квадрат занимает ^ = xI - xb +1 пикселов. Для каждого из ^ пикселов с координатами (x‘s, у), xb < xi < xI , с использованием матрицы T 1 определяются (округленные до ближайших целых) соответствующие координаты исходного квадрата. Таким образом, преобразованный квадрат полностью

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

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

Преобразованный ранее квадрат прибавлением к координатам пикселов значений (хп, уп) переносится в центр эталона на текущем изображении Ж,

определяя тем самым преобразованный эталон. В массив ргх_1 заносятся соответствующие пикселы w изображения Ж, подсчитывается их сумма и в конце определяется среднее значение wm . Затем элементы массива ргх_1 меняются на w' = w - wИ!, суммируются квадраты величин w' и определяется стандартное отклонение а.

Далее на последующем изображении F для каждой допустимой позиции (х7,у7)внутри области размера $еагск_$12в с центром в (хи,уп) выполняются

следующие действия. Преобразованный ранее квадрат прибавлением к координатам (х\, у) значений (х7, у7) переносится на изображение ^ В массив ргх_2 заносятся соответствующие пикселы / изображения ^ подсчитывается их сумма и в конце определяется среднее значение /т. Затем элементы массива ргх_2 меняются на 7 = 7 - /т , суммируются квадраты величин / и определяется стандартное отклонение а 1 . Производится свертка массивов ргх_1 и ргх_2, которая делится на произведение аа. Тем самым определяется значение коэффициента корреляции у(х7, у7) в позиции (х7, у7) . Если значение у(х7, у7) больше текущего значения коэффициента корреляции у(, то производится присвоение у{ =у(х7,у7) и запоминаются координаты (х7, у7) позиции нового максимума коэффициента корреляции. Для многоспектральных изображений процедура применяется к данным каждого спектрального диапазона.

По окончанию обработки всех К_Тат%8 эталонов производится переход к новой паре значений масштаба 5 и угла поворота (р, определяющих новую пару матриц преобразований Т и Т 1.

В дальнейшем из найденных Ырв8 позиций могут быть выделены Ырв8_С позиций, для которых значения коэффициентов Согг не меньше порогового значения ТИгв8_Согг, а длина векторов смещений больше значения параметра Мт_Ьвп$Н. Именно эти Ыро8_С позиций (по умолчанию Ыро8_С=Ыро8) участвуют в последующих операциях: отображение векторов смещений и определение новых эталонов.

В табл. 1 приведены результаты вычислительных экспериментов, выполненных на персональном компьютере с процессором AMD А1:Ыоп(1:т) ХР 3200+. В качестве исходного изображения Ж взято изображение с геостационарного КА «METEOSAT-8» в оптическом

диапазоне (длина волны 620 нм), полученное 15.03.2006 г. в 11 часов 30 мин. В качестве следующего изображения F взято изображение того же канала через 15 минут. Поиск эталонов производился со следующими параметрами: Та^_812в = 15, Grid_Size = 32, Та^_8вагсИ = 40, Таг£_В1^ = 30, Min_St_Dev = 30, ^ит^г_8& = 110. Выделено К_Та^=3461 эталонов. При определении коэффициентов корреляции значение Scale_Delta=0.1. Значение параметра Angle_Delta = 3 градуса. В таблице собраны данные о количестве позиций Ыро8 с ненулевыми векторами смещений, минимальном и максимальном значениях коэффициента корреляции, времени выполнения программы, а также значений Ыро8_С для ТИ^ Согг 0.5, 0.6, 0.7, 0.8, 0.9 при длине векторов смещений Min_Length>3.

Таблица 1

8еа1е_ Мт Беа1е_ Мах Ап§1е_ Мт Ап§1е_ Мах КрОБ Мт_ Согг Мах Согг Ите, Бес ТЬгеБ Согг

0.5 0.6 0.7 0.8 0.9

1 1 0 0 2878 0.292 0.995 2 1173 1149 1105 966 616

0.8 1.2 -45 45 2903 0.413 0.995 62 1221 1201 1163 1037 668

0.8 1.2 -90 90 2904 0.466 0.995 124 1227 1208 1166 1038 668

0.8 1.2 0 357 2904 0.466 0.995 240 1232 1214 1171 1044 672

Анализ данных таблицы показывает, что наиболее существенные изменения по сравнению с первой строкой, относящейся к тождественному преобразованию, содержатся во второй строке (диапазон изменения угла поворота эталона [-45,45] градусов). Последующие строки, связанные с более широкими диапазонами изменения угла поворота, содержат лишь незначительные уточнения. Вероятно, это связано с малым интервалом времени между изображениями.

Рассмотрим теперь вопросы интерполяции пикселов эталонов. Согласно [1] интерполяция необходима, т. к. преобразованная дискретная сетка точек входного изображения, вообще говоря, не совпадает больше с дискретной сеткой выходного изображения и наоборот.

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

изображение и затем выполнить дискретизацию на новой дискретной сетке.

Известно [1, 2], что восстановление непрерывной функции из дискретных отсчетов может рассматриваться как операция свертки

Рг (Х) = ХР(ХШ,п ЖХ - Хт,п X (3)

т,п

где непрерывная интерполирующая маска И является двумерной sinc функцией

= 81п(ЛХ1/ Ау1) 8Ш(я2/ Ау 2) яЫ / Аы ях2 / Ау2

Передаточной функцией для функции рассеяния точки (4) является прямоугольная функция шириной 2к~ = 1 / Ахж :

к~(К) = П(к~ /2,к~ /2), где К = 2к№Ат№, ^ = 1,2. (5)

Заключение об интерполирующей природе ядра свертки (4) может быть получено из следующих соображений. Интерполированные значения в уравнении (3) в точках дискретной сетки Х должны воспроизводить

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

Г1, т = 0, п = 0

Ж(Хт,п) = 1 Л . (6)

[ иначе 0 47

Интерполирующая маска (4) удовлетворяет этому условию. Тот факт, что интерполяция является операцией свертки и вследствие этого может быть описана передаточной функцией в пространстве Фурье, дает удобный инструмент для оценки ошибок интерполяции. Прямоугольная форма передаточной функции для идеальной интерполирующей маски просто означает, что все частоты внутри интервала возможных частот \к^\ < 1 /(2Ах№)

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

Идеальная интерполирующая функция в (3) разделима. Следовательно, необходимо решить только одномерную задачу интерполяции. Как только она решена, мы можем решить и двумерную задачу (и даже задачу интерполяции большей размерности).

Классическим подходом к интерполяции является линейная интерполяция [1]. Интерполируемая точка лежит на прямой, соединяющей соседние точки дискретной сетки. Для целей симметрии примем координаты двух соседних точек равными - 1/2 и 1/2. Это дает уравнение интерполяции

Р(х) = Р1/2 *Р-1/2 + (р1/2 - р—1 /2)х для |х| < 1/2. (7)

Сравнивая уравнение (7) с (3), можно увидеть, что непрерывная интерполирующая маска для линейной интерполяции имеет вид

7 / Ч Г1 - Ы, Ы < 1

Мх) = Г 11 11 . (8)

[ иначе 0

Это треугольная функция. Соответствующая передаточная функция является квадратом функции 8тс:

к) = ^2(як /2) (9)

1 ( ) (як/2)2 . (9)

Сравнение (9) с передаточной функцией (5) для идеальной интерполяции показывает, что линейная интерполяция вводит два искажения:

1. В то время, как низкие частоты (и особенно значение к = 0) интерполируются правильно, высокие частоты слабо уменьшаются по амплитуде, приводя к некоторому сглаживанию. При к = 1 передаточная функция уменьшается примерно на 40 %: И~(1) = (2/^)2 « 0.4. Так как И~(£) не равно нулю при частотах к > 1, появляются ложные высокие частоты.

2. Возвращаясь в уравнении (7) к исходным координатам 0 и 1 соседних точек, получим уравнение линейной интерполяции в виде

Р(х) = (1 - х)Ро + ХР1 . (10)

Пусть теперь х, у - некоторые координаты внутри единичного квадрата с вершинами (0,0) и (1,1) (фактически о < х, у < 1 - дробные части координат X, У , полученных при выполнении геометрических преобразований эталона: X = їпі;(X) + х, У = їпі;(У) + у).

Пусть Р = {р00,р10,р01,рп} - набор пикселов эталона в вершинах единичного квадрата (верхний левый, верхний правый, нижний левый, нижний правый пикселы соответственно). Тогда значение р(х,у) пиксела, которое будет участвовать в вычислении коэффициента корреляции, при линейной интерполяции по двум координатам (билинейной интерполяции) определяется двумерной сверткой

Р( x, У) =(1, х)

( 1 0 Л Р 0 0 Р01 ^ (1 -10 (11

V- 1 1 у V Р10 Р11 .у V 0 V у у

(11)

Развернув это матричное уравнение, получим классическую формулу билинейной интерполяции [2]:

р(х, у) = (1 - х)(1 - у)р00 + (1 - Х)УР0! + х(1 - у)Рю + ХУР 11. (12)

Обобщая линейную интерполяцию, мы можем использовать полином степени N с N + 1 неизвестными коэффициентами , который должен проходить через N + 1 точек:

N

Р( Х ) = Х апХ

(13)

п=0

Из условия интерполяции в точках дискретной сетки р(Хя) = Р„ получается система из N + 1 линейных уравнений относительно N+1 неизвестных коэффициентов аи . В частности, при N = 2 (квадратичная

интерполяция) получается следующая система уравнений для получения коэффициентов полинома:

Р 0 (1 -1 П 0

Р1 = 1 0 0 а

VР2 у 1 1у V а 2 У

ґаЛ

с решением

а

Vа 2 у

( 0 2

-1 0

1 - 2

0 0( Р0Л

Р1

V Р 2 у

(14)

Соответствующая биквадратная интерполяция (интерполяция по 9 пикселам):

1

2 \ Т

р(х, у) = - (1, х, х )ЫРЫ (1, у, у )

(15)

где М - матрица коэффициентов, равная

М

2 0Л 0 1

V

1 - 2 1

, Р - матрица

пикселов, равная Р

Р-1,-1 Ро,-1 V Р1,-1

Р-1,0

Ро,о

Р1,0

Р-1,1 Р0,1 Р1,1 у

Т - символ транспонирования.

Аналогичные соображения применимы для N = 3 - бикубической интерполяции - интерполяции по 16 пикселам. В [1] построен соответствующий интерполяционный многочлен, равный 1

2304

(1, *, *2, *3) МРМ1 (1, ус, у2, уЗ) т

(16)

где *с = * - 0.5, ус = у - 0.5

М - матрица коэффициентов, равная М =

Г- 3 27 27 - 3^

2 - 54 54 - 2

12 -12 -12 12 ^ 8 24 - 24 8 у

Р - матрица пикселов, равная Р =

-2,-2

-1,-2

Р

Р

-2,-1

-1,-1

Р

Р-

-2,0

1,0

Р

Р-

-2,1

1,1

Р0,-

Р

0,-1

Р

0,0

Р0,1

V Р1,-2 Р1,-1 Р1,0 Р1,1 у

Т - символ транспонирования.

Заметим, что для целей симметрии в этом случае координаты дискретной сетки были сдвинуты на -1/2 (начало координат находится в центре сетки 4*4 пиксела), что позволяет более корректно учитывать влияние соседних пикселов.

С увеличением степени интерполяционного многочлена передаточная функция лучше приближается к прямоугольной функции (см. [1]). Сходимость, однако, является слабой.

Анализ данных вычислительных экспериментов позволяет сделать следующие выводы. Прежде всего, при слишком маленьких значениях масштаба эталона (например, 0.5 - уменьшение в два раза) получаются слишком большие значения минимальных коэффициентов корреляции. Объясняется это тем, что уменьшение масштаба соответствует уменьшению пространственного разрешения, что, согласно свойствам преобразования Фурье, приводит к уменьшению ширины спектра изображения эталона, т. е. высокочастотные составляющие эталона подавляются. Но алгоритм поиска эталонов, основанный на поиске областей с большой дисперсией, как раз предполагает наличие в изображении эталона высокочастотных элементов. С другой стороны, слишком большие значения масштаба, кроме существенных затрат времени вычислений, не приводят к каким-либо положительным результатам. Для масштабов эталона представляются оптимальными значения в диапазоне 0.9 - 1.1. Результаты являются гораздо более устойчивыми к

2

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

В табл. 2 приведены данные, полученные с использованием геометрических преобразований эталонов с параметрами 8саЫ_Мт=0.9, 8са1е_Мах=1.1, 8са\е_Веиа=0.1, Angle_Min=-45, Angle_Max=45, Angle_Delta=3 и различными типами интерполяций: 0 - тождественное преобразование, 1 - округление до ближайшего целого, 2 - билинейная, 3 -биквадратная и 4 - бикубическая интерполяция соответственно.

Результаты с использованием билинейной интерполяции, на первый взгляд, являются более лучшими, чем, например, для бикубической интерполяции. Вероятно, это объясняется отмеченными ранее недостатками билинейной интерполяции (появление в спектре ложных высоких частот). Учитывая небольшую разницу в затратах времени вычислений, более предпочтительной является бикубическая интерполяция.

Таблица 2

Interpol Npos Min Corr Max Corr Time, s Thres Corr

0.5 0.7 0.8 0.9

0 3479 0.353 0.998 3 3351 1939 1675 1139

1 3507 0.450 0.999 301 3438 2045 1790 1210

2 3515 0.491 0.998 307 3467 2044 1828 1282

3 3524 0.490 0.998 309 3468 2037 1813 1263

4 3505 0.500 0.998 310 3453 2047 1816 1248

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

Работа выполнена частично при финансовой поддержке Российского фонда фундаментальных исследований (проект № 07-07-00085).

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Bernd Jahne. Digital Image Processing. Springer-Verlag, Berlin, Heidelberg, 2005. -607 pp.

2. Гонсалес Р., Цифровая обработка изображений / Р. Гонсалес, Р. Вудс. / Пер. с англ. под ред. П.А. Чочиа. Москва: Техносфера, 2005. - 1072 с.

3. Забелин В.А.Идентификация контрольных точек в аэрокосмических изображениях. / В.А. Забелин, В.П. Пяткин / В сб. “Математические и технические проблемы обработки изображений”. ВЦ СО АН СССР, Новосибирск, 1980, с. 50 - 59.

4. MSG Meteorological Products Extraction Facility. Algorithm Specification Document. - Doc. No. EUM/MSG/SPE/022. Issue 2.6. 1 June 2004.

5. Фоли Дж. Основы интерактивной машинной графики. Т. 2. / Дж. Фоли, А. вэн Дэм./ Москва, Мир, 1985. 368 с.

© А.А. Бучнев, В.П. Пяткин, 2008