CC BY
13УДК 681.518.3
САМОНАСТРАИВАЮЩИЙСЯ АЛГОРИТМ ВЫСОКОТОЧНОГО СОВМЕЩЕНИЯ АФФИННО-ПРЕОБРАЗОВАННЫХ ИЗОБРАЖЕНИЙ1
П.К. Кузнецов, Е.Ю. Чекотило, Б.В. Мартемьянов
Самарский государственный технический университет 443100, г. Самара, ул. Молодогвардейская, 244
В статье предложен самонастраивающийся алгоритм высокоточного совмещения аффинно-преобразованных изображений, которые получены оптико-электронными системами при съемке подвижных объектов. Алгоритм основан на применении метода функционализации изображения для получения оценок параметров аффинного преобразования изображений.
Ключевые слова: яркостный объект, оптико-электронные системы, аффинное преобразование, высокоточное совмещение изображений, метод функционализации.
Введение
Высокоточная оценка параметров движения и смещения положения изображения подвижного объекта в последовательности кадров его изображений, получаемых с помощью оптоэлектронных систем наблюдения, является одной из актуальных задач в различных технологиях визуального контроля [1, 5]. Например, высокоточное совмещение изображений используется при создании управляющих систем роботизации производств, систем навигации подвижных транспортных платформ, в частности, летательных и космических аппаратов, а также при решении более специфических задач, например, устранения “смаза” изображений, получаемых с борта летательных или космических аппаратов при нештатных режимах съемки.
Известные точные методы определения параметров движения изображений [1], основанные на вычислении взаимно корреляционных функций совмещаемых изображений, обладают большой вычислительной сложностью и слабо защищены от ложных отождествлений фрагментов изображений (ложное опознавание). Вследствие этого они мало приспособлены для реализации в реальном времени, особенно при идентификации параметров сложного, «аффинного», движения. Другая группа методов аналогичного назначения основана на вычислении производных по времени и пространству от функции плотности изображения. Это так называемые «градиентные методы» [2-5]. Они обладают высоким быстродействием, но имеют значительные методические погрешности, связанные с дифференцированием изображений, которые, вообще говоря, не дифференцируемы. Еще один недостаток этой группы методов состоит в том, что они работоспособны лишь при малых смещениях изображения от кадра к кадру.
В данной работе предлагается использовать компенсационный вариант метода функционализации [6-9] для получения оценок параметров взаимного смещения аф-финно - преобразованных изображений, расположенных в двух последовательных кадрах анализируемой серии изображений. Метод функционализации позволяет без
1 Работа выполнена при финансовой поддержке целевой программы «Научные и научнопедагогические кадры инновационной России на 2009-2013 гг.» (проект П2612).
Павел Константинович Кузнецов - д.т.н., профессор.
Борис Викторович Мартемьянов - к.т.н., доцент.
Елена Юрьевна Чекотило - к.т.н.
применения поисковых алгоритмов и вычисления сверток с высокой точностью определять параметры смещения как протяженных, так и малоразмерных объектов со случайным распределением яркостей, как при малых, так и при больших смещениях. Метод является обобщением градиентного метода и основан на использовании соотношения (ФС - уравнения), функционально связывающего параметры движения наблюдаемого объекта с изменениями измеримых характеристик изображения [6]. Под «измеримыми» понимаются характеристики, значения которых определены на подобластях изображения ненулевой меры. Здесь этот метод используется для определения параметров движения изображения, обладающего тремя степенями свободы (движение смещения по осям декартовой системы координат, лежащей в плоскости изображений оптической системы, и вращение относительно нормали к этой плоскости). При таком движении трансформация изображения от кадра к кадру, получаемого оптико-электронной системой, описывается аффинным преобразованием без растяжения.
Постановка задачи
Возьмем две плоскости Р1 и Р2 изображений с заданными на них декартовыми системами координат (СК) и совместим их так, чтобы совпали как начала двух СК, так и их одноименные оси. В результате, на двух плоскостях изображений будем иметь единую СК, которую обозначим как ОХУ. Будем считать, что на плоскостях Р1 и Р2 заданы изображения Е1(г) и Е2 (Г), соответственно, функциями распределения плотности (освещенности) изображений:
Ев (г) е Я1, г = (х,у)т, 0 < Е(г) < М, (1)
где в е{1,2} ; ()т - знак транспонирования; х, у е Я1 - координаты в СК ОХУ.
Условия непрерывности и дифференцируемости на функции плотности изображений не накладываются.
На плоскости Р1 расположим некоторым образом область анализа (ОблА) ¥Э, которую будем называть эталонной. Фрагмент изображения Е1(г) из области ¥Э будем трактовать как эталонное изображение. Обозначим его как Ег(г) (Ег(г) = Е1(Я: г е¥ * )).
Считается, что на плоскости Р2 существует ОблА ¥к, которая содержит фрагмент Ек (г) изображения Е2 (г), являющийся результатом аффинного преобразования изображения Ег (г) (Ек (г) = Е2(г : г е ¥к)).
Необходимо на плоскости Р2 найти область ¥к по заданной области ¥Э.
Будем считать, что области ¥Э и ¥к имеют прямоугольную форму. Стороны области ¥Э сориентируем параллельно осям системы координат ОХУ. Ограничения, накладываемые здесь на форму и ориентацию области анализа ¥Э, не существенны, сделаны для упрощения изложения и легко могут быть сняты.
Зададим группу П = {л1} аффинных преобразований ж1 (без растяжения) области ¥Э в произвольную область ¥! плоскости Р2:
¥! = ж1 (¥Э ) = {г :(г = Л1- (х, у, 1)т ) л (х, у) е ¥Э } , (2)
"1 0 Л(І)" cos ^3(І) - эт ^3(І) 0" соэ р3 (І) - Рз(І) Л(І)"
А1 = 0 1 Л(І) X sin ^Р3(І) соэ ^Р3(І) 0 = ИП ^з(І) соэ р3 (І) _Р2 (І)
0 0 1 0 0 1 0 0 1
- матрица I -го аффинного преобразования в однородной системе координат;
^93(1) е Я1 - угол поворота области ¥Э относительно начала СК;
д(1), ^92(1) е Я1 - сдвиги области ¥Э по осям ОХ и ОУ, соответственно.
Обозначим как ^(1) = (д(1), ^2(1), ^Р3(1))Т - вектор параметров аффинного преобразования ж1.
Матрица Л1 предполагает преобразования вектора из области ¥Э в последовательности: поворот, затем сдвиг.
Требуется на основе анализа эталонного изображения Ег(г) и изображения Е2(г) найти такое преобразование жк е П , что
¥к = жк (¥Э ). (3)
Тогда р(К) = (р1(К), р2(К), р3(К ))Т - вектор параметров искомого аффинного преобразования жк эталонной области ¥Э в область ¥к
Итерационная процедура определения параметров аффинного
преобразования
Для решения поставленной задачи используем следующую итерационную процедуру вычисления последовательности векторов параметров преобразования (2):
Р(0) = Д»
р (І +1) = ^ (І) + Л-А^ (I), (4)
1 = 0,1, 2,..., т;
где р0 - начальный вектор параметров преобразования (2);
Л - матрица коэффициентов обратной связи:
Л = diag {Яі = 1...3. (5)
На каждой итерации процедуры (4), как и во всех известных способах совмещения изображений, основанных на использовании «градиентного метода» [4], вычисляют поправки А^(І) к текущей оценке ^З(І) вектора параметров аффинного преобразования. Итерации процедуры определения параметров аффинного преобразования заканчиваются, когда при некотором I = т выполняется условие малости величины приращения А^(!):
(таХ { \АРі (І)|, \АР2 (і)| Ь Є* ) Л (| Р3 (І)| ^ Єа), где в^, ва> 0 - заранее заданные числа.
Поправки А^(І) вычисляются на основе анализа изображений в эталонной об-
На первой итерации алгоритма, исходя из априорных оценок ожидаемого смещения, задают начальное значение р0 вектора р(0) параметров аффинного преоб-
Итерационная процедура приводит к получению оценки искомого вектора ^(К) с субпиксельной точностью.
Варианты итерационной процедуры в различных литературных источниках отличаются значениями коэффициентов обратной связи Л і и способом вычисления поправки А^(І).
Основным отличием процедуры (4) от известных процедур является применение для вычисления указанной поправки А^(І) метода функционализации [6-7]. Здесь метод функционализации применяется следующим образом.
На каждой I -ой итерации, включая первую, используя текущий вектор ^З(І), по
заданным вершинам прямоугольной области ¥Э с помощью преобразования ж1 вычисляют вершины прямоугольной области ¥І. Обе области (¥Э и ¥І) покрывают идентичной системой окон анализа (ОА) {р“ : а є {Э,I}і = 1 N}, і - но-
Э I
мер ОА) и для каждой одноименной по і пары ОА (Рі , Рі ) формируют систему уравнений функциональной связи вида [8]:
где К* (г) - непрерывная, дифференцируемая почти всюду по всем своим аргументам и равная нулю на границе ОА О1* функция веса;
ласти ¥Э и в последовательности областей ¥1, перемещаемых на каждой итерации.
разования ж . При этом должно выполняться условие:
¥0 п¥к ^0.
(7)
В1 • Ар (I) = АФ1, где АФ1 = (ФіІ - ФЭ), фа (а є {I, Э}) - функционал вида:
(8)
фа =Цка (г) • £а(г) Л , ра
(9)
В1 = Й ] - матрица коэффициентов, элементы которой определяют из выраже-
ний:
Ьік = 0,5 • (Ьік + Ьік ) ;
Ьа = \\\к'ха (?)соэР3 + куа (гМпР3. Еа (г;
где К'х* иК'у* - частные производные функции К* (г) по х и у, соответственно;
Е* (г) - фрагмент изображения в ОА О* .
Система (8) в общем случае является переопределенной и может быть решена, например, методом квазиобращения.
Приведенный алгоритм требует настройки параметров Л г и размеров ОА О* по условиям сходимости процедуры.
В [10] найдена аналитическая оценка сходимости итерационного метода на примере тестового моночастотного изображения вида
Е (г) = sin( Л2ат) + Л1ат) (11)
отдельно для плоско-параллельного движения (р = 0) и вращения в плоскости изображения (р1, р2 = 0).
Оказалось, что достаточным условием сходимости итерационной процедуры для случая плоско-параллельного движения является выполнение следующих условий, накладываемых на начальное рассогласование положения (перенос) объекта в анализируемых областях и значения коэффициентов X г (г = 1,2 ) в (5):
А. р2 е[-^], | Л21 < 4®2С082 р (21)@ , (12)
3т 3т 2
Условием сходимости итерационного процесса для случая чистого вращения при изменении угла поворота р3 в пределах р3 е [0, п/2] является:
ад^3(1)
Ср3(1)
(13)
где Ар (I) - приращение оценки угла поворота на итерации с номером I.
Таким образом, процесс сходится, если начальное рассогласование фрагментов
¥Э и ¥к по каждой из координат смещения не превышает трети периода Т = —
ю
тестового изображения: \рх\, \р2\ < Т/3 . Этот факт положен в основу процесса самонастройки алгоритма совмещения изображений.
В общем случае изображение по своей контрастно-частотной характеристике не является моночастотным. Получение условий сходимости итерационного процесса для случая, когда Фурье-спектр изображения полигармонический, является в общем случае достаточно сложной задачей.
Здесь предлагается упрощенная, эвристическая, методика настройки параметров алгоритма совмещения изображений для случая, когда в его Фурье-спектре существуют доминирующие частоты, т.е. частоты, на которых спектр имеет явно выраженные максимумы. Эти доминирующие частоты можно выделить с использованием полосовых фильтров и выбрать те из них, на которых условия сходимости (12) и (13) процесса совмещения изображений выполняются.
Заметим, что функционал (9), сам по себе, является фильтром с ядром К (г). Поэтому, настраивая параметры ядра и ОА, можно выделить частоту, для которой условия сходимости (12) и (13) выполняются. Для простоты реализации используем
«пирамидальную» функцию веса, которая (для одиночного ОА с центром в начале локальной системы координат ОА) в границах ОА имеет вид:
V(х, у) (у| < \х1у / 1Х (к(х, у)= (1 - \х\11Х)),
V(х, у) (у| > |х1у / Іх (к(х, у) = (1 -1у\/ 1у )),
где (2Іх ,2Іу ) - размеры ОА по осям х и у, соответственно.
По аналогии нетрудно записать функцию веса для функционала с произвольным расположением на изображении окна, на котором функционал (9) определен.
Необходимая фильтрация обеспечивается за счет подбора параметров (Іх, І ) весовой функции.
На рис. 1. представлена амплитудно-частотная характеристика фильтра с «пирамидальной» функцией веса.
В реальных наблюдениях функция распределения облученности изображения имеет случайный характер. Поэтому, несмотря на то, что измерения производятся на конкретных реализациях изображения, корректнее использовать не частотноконтрастный спектр изображения, а спектральную плотность мощности (энергетический спектр) видеосигнала (рис. 2):
ТО ТО / ч
рв(^х,^) = І |в(x,у) • е‘^хХ+ШуУиЫу, (15)
-ТО -ТО
где В(х,у) - автокорреляционная функция изображения.
Р и с. 1. - Амплитудно-частотная „ „
. Р и с. 2. - Спектральная плотность
характеристика фильтра _ _
„ , „ мощности изображения типа «Г оры»
с «пирамидальной» функцией веса
Эти выводы легли в основу алгоритма настройки параметров процесса измерения, входящего в состав итерационной процедуры совмещения изображений.
Алгоритм настройки параметров процесса совмещения изображений
Исходными данными являются:
- функция спектральной плотности изображения ¥в (^х, wy);
- априорная оценка начального смещения изображения ро фрагмента ¥к
относительно фрагмента ¥Э;
- предельно допустимое значение соотношения сигнал/шум - у > 0.
1. Определение допустимых доминирующих частот
1.1. Определяем доминирующие частоты Е = {мхУ,МуУ)|, (V = 1...иу). Доминирующие частоты - это частоты, являющиеся экстремумами функции спектральной плотности ¥В (мх, Му), для которых выполняется условие:
сигнал/шум > у.
1.2. В соответствии с достаточным признаком сходимости итерационной процедуры (12) по текущей оценке смещения изображения р(1) вычисляем предельно допустимую по условиям сходимости частоту спектральной плотности сигнала:
м (I) = Л^ (I), ^ (I)): ^ (I) = 1п V м (I) = 1п
пр V / 1 \ хпр V ' ? у пр V /) хпр V ' з г-^1 (I) упр^' 3 (I)
1.3. Формируем множество пар допустимых доминирующих частот:
Е* С Е , Е*: {^х ^, Му *^: ™х ^ ™хпр ® V Му ^ < Мупр С)}
2. Настройка параметров процесса совмещения
2.1. Для каждой пары допустимых доминирующих частот Е * выбираем типоразмер (/хV х 1уV ) окна анализа, максимизирующий амплитудно-частотную характеристику Ж(мх, му, 1х, 1у) фильтра с функцией веса \К(г)] г на этих частотах:
тах Ж К Му *% 1Х, 1у ) = Ж * К *, ^ С, ).
[х, [у
Каждый типоразмер окна анализа используем далее для образования покрытия фрагмента изображения. В результате получаем несколько покрытий (система слоев покрытия) одних и тех же участков изображения окнами разного типоразмера. Слои покрытия образуем последовательно, начиная со слоя, соответствующего окну наиболее крупного размера. Для текущей системы покрытий определяем вектор - столбец В1 коэффициентов системы уравнений функциональной связи (8). Для каждого типоразмера окна анализа (/,/ х ] выбираем критическое (максимальное) значение
коэффициента обратной связи X *р , отвечающее достаточному признаку сходимости (12). Назначаем значение коэффициента X3 = 10 • X3кр .
3. Слои покрытия ранжируем. Правило ранжирования не жестко. Например, слои могут ранжироваться по площади окон анализа или по отношению доминирующих частот, им соответствующим. Ранжируем и условия останова итерационного процесса (6). Вычисляем для каждого слоя новое значение оценки смещения р(1 +1), начиная с типоразмера наивысшего ранга, соответствующего самым низким доминирующим частотам. Полученные оценки параметров обратного оператора используем в качестве начальных данных для уточнения смещений по слою следую-
щего ранга и так до исчерпания всех слоев покрытия на I - ой итерации. Повторяем пункты 1 - 2. Процесс останавливаем при выполнении условия (6) для покрытия самого высокого ранга. В результате получается двойной итерационный процесс вычисления параметров аффинного преобразования. Внутренние итерации последовательно смещают области ¥1, приближая изображение в очередной из этих областей к изображению в эталонной области ¥Э. Внешние итерации перебирают слои покрытия фрагментов ¥1 и ¥Э изображений системами ОА \р] ] «от грубого к точному».
Для апробации предложенной итерационной процедуры и определения погрешности получаемых оценок было применено компьютерное моделирование. В качестве изображений использовались фрагменты типа «Река», «Г оры» космических снимков поверхности Земли формата 1000*1000 пикселей и по 10 бит/пиксель.
Экспериментально установлено, что параметры аффинного преобразования вычисляются с заданной погрешностью не более чем за 8 - 9 итераций. Оценка параметра рз углового смещения изображения в экспериментах имела постоянное смещение порядка 1% как при простом повороте (без смещения), так и при полном аффинном движении (рис. 3.). Этот факт можно объяснить влиянием дискретного представления видеосигнала.
0,06
0,05
0,04
0,03
0,02
0,01
0
а) составляющие переноса p, p
5 6 7 8 9
б) вращение p2 итерации
Р и с. 3. - Погрешность в пикселях вычисления оценок переноса Sj , 82 и вращения 83 (в относительных единицах) аффинно -преобразованного изображения при начальном
рассогласовании pj, Р2 = 5 пикс., p3 = 50
Проведенные численные эксперименты показали работоспособность метода и его существенное преимущество по быстродействию в сравнении с другими известными методами.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Котов В.В. Оценка смещения фона при многокадровой обработке в системах технического зрения/ В.В. Котов // Мехатроника, автоматизация, управление. - 2005. - №3. - С. 26-28.
2. B.K.P. Horn, B.G. Schunck, Determining optical flow, Artif. Intell. 17 (1981) 185-203.
3. B. D. Lucas and T. Kanade. An iterative image registration technique with an application in stereo vision. In Seventh International JointConference on Artificial Intelligence, pages 674-679, Vancouver, 1981.
4. Y. Keller, A. Averbuch, Global parametric image alignment via high-order approximation. Computer Vision and Image Understanding 109 (2008) 244-259
5. Современные информационные технологии в задачах навигации и наведения беспилотных маневренных летательных аппаратов / Под ред. М.Н. Красильщикова, Г.Г. Себрякова. - М.: ФИЗМАТ-ЛИТ, 2009. - 556 с. - ISBN 978-5-9221-1168-3.
4
6. Абакумов А.М. и др. А.С. № 753244 G01 C21/00/ (СССР). Способ измерения параметров движения протяженного объекта со случайным распределением яркостей. Приоритет 14.02.1979г. Опубликовано 20.09.2006г. /Абакумов А.М., Бережной И.А., Денкевиц В.А., Кузнецов П.К., Ляпидов В.С., Мишин В.Ю., Семавин В.И.
7. Кузнецов П.К., Семавин В.И. Метод определения параметров движения яркостного поля//Известия ВУЗов. Приборостроение. - 1990. - №6. - С. 26 - 30.
8. Кузнецов П.К., Семавин В.И., Мишин В.Ю., Владимиров М.В. Метод функциональных преобразований в задаче определения скорости движения яркостных полей// Вестник Самарского технического университета. - 1994. - №1. - С. 66 - 76.
9. Кузнецов П.К., Мартемьянов Б.В., Семавин В.И., Чекотило Е.Ю. Метод определения вектора скорости движения подстилающей поверхности// Вестник Самарского государственного технического университета. Серия Технические науки. - 2008. - № 2(22). - С. 96 - 110.
10. Кузнецов П.К., Чекотило Е.Ю., Мартемьянов Б.В. Исследование сходимости итерационной процедуры определения параметров движения изображений методом функционализации// Вестник Самарского государственного технического университета. Серия Технические науки. - 2010. -№2(26). - С. 80 - 85.
Статья поступила в редакцию 20 января 2011 г.
UDC 681.518.3
SELF-ADAPTIVE ALGORITHM OF HIGH-PRECISION SUPERPOSITION OF THE AFFINE-TRANSFORMED IMAGE
P.K. Kuznesov, E. Yu. Chekotilo, B. V. Martemyanov
Samara State Technical University
244, Molodogvardeyskaya st., Samara, 443100
This article presents a self-adaptive algorithm of high-precision superposition of the affine-transformed images. Methodology is based on using method of functionalization of parameters of geometrical mismatch of the affine-transformed images of mobile objects.
Keywords: brightness object, optical-electron system, affine transformation, high-precision superposition of images, functionalization method.
Pavel K. Kuznesov - Doctor of Technical Sciences, Professor.
Elena Yu. Chekotilo - Candidate of Technical Sciences.
Boris V. Martemyanov - Candidate of Technical Sciences, Associate professor.
