Блочные алгоритмы обработки изображений на основе фильтра Калмана в задаче построения сверхразрешения Текст научной статьи по специальности «Компьютерные и информационные науки»

БЛОЧНЫЕ АЛГОРИТМЫ ОБРАБОТКИ ИЗОБРАЖЕНИИ НА ОСНОВЕ ФИЛЬТРА КАЛМАНА В ЗАДАЧЕ ПОСТРОЕНИЯ СВЕРХРАЗРЕШЕНИЯ

Сирота А.А., Иванков А.Ю.

Воронежский государственный университет

Аннотация

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

Ключевые слова: блочная обработка изображений, фильтр Калмана, сверхразрешение.

Введение

В настоящее время большое значение имеет развитие технологий обработки изображений и видеоданных. При этом одним из основных требований, предъявляемых к графическим данным, является высокое разрешение (ВР), способное обеспечить требуемый уровень детализации снимков и различение достаточно мелких деталей. Другим требованием является отсутствие на изображениях шумов и помех, искажающих информацию и мешающих работе алгоритмов распознавания образов. Однако не всегда складываются условия, необходимые для получения высококачественных графических данных. В силу естественных ограничений возможностей оборудования, применяемого для получения изображений, качество отдельных снимков не всегда может быть удовлетворительным. Одним из вариантов повышения качества изображений является применение методов сверхразрешения (СР) [1-6], позволяющих решить эту задачу с помощью алгоритмической обработки последовательностей изображений. Повышение качества в данном случае происходит за счёт использования нескольких кадров, содержащих схожую информацию. Это позволяет добиться эффекта восстановления изображения ВР из некоторого количества изображений низкого разрешения (НР). При этом изображения НР представляют собой различные отображения одной и той же сцены, смещённые на нецелое число пикселей. В этом случае новая информация, содержащаяся в каждом изображении НР, может быть использована для получения изображения ВР.

Существует много методов и алгоритмов построения сверхразрешения [1-6]. Одним из эффективных подходов с точки зрения получения оптимальных (в классе линейных) оценок является использование алгоритмов калмановского типа [4-10]. Однако общей проблемой при реализации таких алгоритмов в задаче сверхразрешения является большая размерность обрабатываемых векторов данных, получаемых после развёртки матрицы изображений (порядка 105... 106 элементов). Соответственно этому размерность матриц (например, матрицы ковариации ошибки), вычисляемых и обращаемых в процессе реализуемой обработки, становится огромной, что приводит к перерасходу используемых вычислительных ресурсов и возникновению случаев расходимости.

Перспективными подходами, которые могут обеспечить эффективное применение алгоритмов калманов-ского типа, являются методы распараллеливания вход-

ных данных и приёмы, позволяющие избежать обращения матрицы ковариации. К первым относятся матричные алгоритмы, заключающиеся в формировании блочной матрицы, содержащей в себе всю необходимую для данного этапа фильтрации информацию, и последующем её приведении к требуемому нижнему или верхнему треугольному виду [8, 9]. К последним относится техника скаляризации наблюдений [8], позволяющая избежать процедуры обращения.

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

Одним из способов сокращения объёма вычислений является блочная обработка изображений в процессе фильтрации. В [11] упомянуто, что для обработки изображений можно при незначительной потере в точности использовать сетки небольшого размера (блоки). В [9] приведена структура блочного фильтра Калмана, полученная за счёт организации векторной модели авторегрессии в блочной форме. Данная модель фильтра учитывает корреляцию соседних блоков, следующих друг за другом в пределах строки блоков, тем самым устраняя краевые эффекты на стыке блоков в пределах строки. Достигается это за счёт значительного расширения вектора состояний фильтра Калмана и, соответственно, ковариационной матрицы, что является серьёзным недостатком упомянутого метода, так как ведёт к увеличению объёма вычислений. Этот подход развит в [8], где приведена структура фильтра Калмана с сокращённым обновлением (reduced update Kalman filter). Данный алгоритм аппроксимирует пространство состояний в малом регионе с целью минимизации вычислений, связанных с обновлением ковариационной матрицы.

Применительно к задаче построения сверхразрешения изображений фильтр Калмана рассматривается в работах [3-6]. В [3, 4] представлены аппроксимации алгоритма калмановской фильтрации на основе мето-

да сопряжённых градиентов (Я-ББ) и наименьших средних квадратов (Я-ЬЫБ), в которых не используется обращение матрицы ковариации ошибки. Работы [5, 6] описывают аппроксимации фильтра Калмана: для случая стационарной, независимой от времени матрицы ковариации ошибки [5], и с применением аппроксимации ковариационной матрицы с помощью некоторых шаблонов, полученных опытным путём [6]. Несмотря на то, что обрабатываемые матрицы фильтра Калмана в указанных работах имеют разреженную структуру, позволяющую применять к матричным уравнениям различные итеративные методы решения, обработка всё равно занимает значительное время и требует больших объёмов памяти. В данной работе продемонстрированы варианты алгоритма построения сверхразрешения на основе калмановской фильтрации с использованием блочной обработки, позволяющей сократить объём вычислений. Если в процессе обработки не осуществлён учёт информации из соседних блоков изображений, на восстановленных изображениях возникают новые артефакты (краевые эффекты). Включение в обработку информации из соседних блоков предполагает расширение вектора состояний и ковариационной матрицы, что увеличивает объём вычислений. В настоящей работе предложены варианты решения данной проблемы.

Математическая модель фильтра Калмана в задаче построения сверхразрешения изображений

В качестве наблюдаемых данных выступает последовательность изображений НР у к, где к = 1,...,К . Изображения имеют размер МхМ пикселей. Для удобства последующего анализа каждое изображение представлено как одномерный вектор, полученный после развёртки исходной матрицы изображения по столбцам. В качестве оцениваемых данных рассматривается последовательность изображений хк с высоким разрешением, где каждое изображение имеет размер Ь х Ь , (Ь > М).

Последовательность изображений ВР хк удовлетворяет следующему уравнению:

(1)

Ч+1 = гк хк + «к,

где хк+1 - изображение ВР, полученное из предыдущего изображения хк за счёт перемещения камеры и/или объекта в процессе получения изображений; Fk - матрица сдвига размером Ь2 хЬ2, характеризующая геометрические взаимные деформации (смещения) изображений; «к - гауссовский шум с нулевым математическим ожиданием и ковариационной матрицей

Ок =оке)21 размером Ь2 х Ь2.

Каждое изображение НР связано с соответствующим изображением ВР соотношением

У к = нк хк + ^, (2)

где ук - изображение с низким разрешением; хк -изображение с высоким разрешением, из которого

производится ук; Нк - матрица децимации размером М 2 х Ь2 , которая осуществляет преобразование хк в ук; ук - вектор аддитивного гауссовского шума с нулевым математическим ожиданием и матрицей ковариации Як = окя)21 размером М2 хМ2.

Следует отметить, что матрицы Нк, Fk, Кк и Qk, определяющие систему пространства состояний, считаются известными. В данной работе предполагается глобальный характер смещений между соседними изображениями, подразумевающий одинаковую величину смещений для всех пикселей одного изображения, что соответствует перемещению камеры относительно сцены.

Матрица сдвига Fk имеет структуру, обусловленную ядром фильтра для интерполяции изображения

% = Цлкг ||, в результате свертки с которым исходного изображения получается смещённое изображение:

Fk =1 ^ 1;

Г (к) =„ ¿гЛ ~ 1

г = (у — 1)Ь+7, г = (у + п)Ь+7 + т, (3)

0 < 7 + т < Ь, 0 < у + п < Ь; 0, г ф( у — 1) Ь+1, г ф (у + п)Ь+7 + т, 0 < I + т < Ь, 0 < у + п < Ь, где 7 = 1,...,Ь ; у = 1,...,Ь ; р = 1,...,Ж^ ; д = 1,...,Ж^ ; т = р — (Ж^ +1)/2; п = д — (Ж^ + 3)/2; матрица

фильтра пк имеет размер Ж^ х Ж^, Ж^ - нечётное число.

Матрица децимации Нк имеет похожую структуру. В случае если Ь / М = | и | - целое число, то:

Нк =| И;

А = е(к)

"к = ер,д

(4)

"г,г

(к)

г = (у — 1)М +7, г = (7 + п)Ь + 7 + т, 0 < 7 + т < Ь, 0 < у + п < Ь; 0, г ф(у —1)М+7, г ф (у + п)Ь+7 + т, 0 < 7 + т < Ь, 0 < у + п < Ь, где 7 = 1,...,М ; у = 1,...,М ; р = 1,...,Же; д = 1,...,Же ; т = р—(Же+1)/2; п = д—(Же+ 3)/2; 7 = [(/ — 0.5)-|т] ; у = |(у — 0.5) |. Оператор [•] означает округление

до ближайшего целого в большую сторону. Ядро фильтра Ак размером Же х Же (Же - нечётное число) характеризует функцию рассеяния точки (ФРТ) датчика камеры. В 0к также следует учитывать весовые

рд

коэффициенты интерполяции соседних пикселей в случае, если т не является целым.

Фильтр Калмана позволяет по полученной совокупности наблюдений ук ={у1,...,ук} за к шагов получить оптимальные в среднеквадратичном оценки хк|к_1, хк|к вектора хк. Оценка Хк|к является изображением СР. Это оценка состояния объекта вида «точка в точке», полученная на основе апостериорной

плотности распределения Р(хк | ук). Оценка Xк|к_1 -

это прогноз состояния объекта на один шаг вперёд, в данном случае - предсказание следующего изображения ВР. Основой получения прогноза является апостериорная плотность распределения

Р(хк|к_1 | ук-1). Переход от прогноза Xк|к_1 к оценке Хк|к осуществляется на основе матрицы ковариации

ошибки Рк|к_1 = м [вк екк ] , где ошибка £к = Хк _ Xк|к_1 характеризует разницу между изображением ВР хк и предсказанным изображением ВР Хк|к-1.

В классической форме фильтра Калмана при поступлении очередного изображения НР ук производится коррекция оценки Xк|к :

хк|к = Хк|к_1 + Щк (ук Нк Хк|к_1 )

щ = Ук и-1, ик = НкРк|к_1Икт + Як, Ук = Рк|к_1Нк, Рк|к = Рк|к_1 _ ЩкУк. Рекуррентное выражение для оценки Х к + 1-м шаге имеет вид

к+1|к

(5)

(6)

на

Хк+1|к = гкХк|к = ркХк|к_1 + ркщк (ук _ НкХк|к_1). (7) Рекуррентное выражение для матрицы Рк+1|к на к + 1-м шаге имеет вид

Рк+1|к = ^ ^ + Ок.

(8)

Блочная обработка изображений с применением фильтра Калмана

Как уже упоминалось, главная проблема реализации метода калмановской фильтрации в классической форме - чрезмерно большой размер матриц, нуждающихся в обработке при получении сверхразрешения изображений. На практике это приводит к возникновению эффектов расходимости алгоритма фильтрации. Решением служит обработка изображений блоками с использованием специальных приёмов «сшивания» этих блоков в общее изображение.

Изображения делятся на В блоков по вертикали и горизонтали. Размер блока НР ¿ы х ¿ы (¿ы ■ В = М),

размер блока ВР ¿ь = ¿Ь1 • т (¿ь • В = Ь). Из-за ограниченности размеров блоков изображения при использовании матрицы сдвига на краях блока смещённого

изображения проявляется осцилляция яркости, что проиллюстрировано на рис. 1а. Аналогичный эффект возникает и при прореживании изображения.

Рис. 1. Проявление краевого эффекта: (а) при смещении изображения на 2 пикселя по вертикали и на 0,5 пикселя по горизонтали; (б) при повышении разрешения в 2 раза

Для преодоления данного эффекта требуется сшивка блоков изображений. Достичь этого можно за счёт расширения на к-м шаге блоков изображения ВР до размера ¿Ье х ¿Ье, где ¿Ье = ¿ь + 2Ат, Лот = (N _ 1) / 2 , N - ширина ядра фильтра, равная Ыв или N . Для смещения блока используется расширенная матрица сдвига Fke) размером ¿Ь х ¿1е, а для прореживания блока применяется расширенная матрица Нке) размером х ¿Ье. Оценка Хк|к также расширяется:

Хк = х(к/Ч ;

х(к) =

г (к|к)

, г = Ат +1,...,Ат + Ь,

(9)

(]_1_Ат)Ь+г_Ат'

] = Ат +1,..., Ат + Ь, с, г = 1,..., Ат, ] = 1,..., Ат, и г = Ь + Ат +1,..., Ь + 2Ат, 1 = Ь + Ат +1,..., Ь + 2Ат,

где с - некоторая константа.

Из расширенного изображения-оценки Хк извлекаются блоки Оразмером ¿Ье х¿Ье

Гк|к =1 |„Р,9,(к|к)|| .

'г, 1

»,9,(к|к) = х(к|к) 5? ■ — л'

'г, 1

где г = 1,...,

Ье

( Р_1>ь +г',(9_1>ь+1 '

1 = 1,..., 5ье , Р = 1,..., В, Я = 1,..., В

(10)

(к) Р,Я

по-

(11)

Каждый блок разворачивается в вектор g' сле чего к нему применяется операция сдвига:

Х(к+1|к) = ^е) Гк|к)

Р,Я к &>Р,Я '

где вектор ХР,+1|к) является блоком оценки Хк+1|к .

Процесс блочного прогнозирования иллюстрирует рис. 2.

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

!+2Дт расширенное %цк

Рис. 2. Процесс блочного прогнозирования изображения Процесс получения оценки НР показан на рис. 3. Требуется расширить оценку по формуле (9) и разбить её на блоки, как в (10). Каждый блок разворачивается в вектор —1), после чего подвергается децимации:

(к) = Н(е)„(к|к—1) к

(12)

"р,д к Ьр,д

где р = 1,...,В, д = 1,...,В. В результате получается развёрнутый в вектор блок изображения НР гк, с помощью которого можно оценить отклонение от ук.

М

Ь+2Кт расширенное хкк

Рис. 3. Процесс блочной децимации изображения

Расширенные матрицы Fk:e) и Нке) получаются за счёт добавления новых столбцов на места, соответствующие новым элементам в расширенном блоке ^р^)

или gРг+1|k). Таким образом, значения пикселей на границах блоков хр^+1|к) и г® вычисляются с помощью

элементов, содержащихся в перекрытиях блоков. Формула (5) для коррекции блока имеет вид:

(к|к) = х(к|к—1) + ' ( у(к) — г(к) ) р,д р,д к I * р,д р,д I

(13)

где х™)- развёрнутый в вектор блок изображения-с(рк'к—1) - вектор блока изображения-

оценки хк|к;

прогноза хк|к—1; ур^ - вектор блока обрабатываемого изображения НР ук; г(рк)д - вектор блока изображения НР Zk ; р = 1,...,В , д = 1,...,В - индексы блока. Размер матрицы усиления 'к соответствует размеру блока и равен х 'Ьэ1 .

В процессе прогнозирования и коррекции матрицы ковариации ошибки необходимо учесть, что расчёт блока оценки основан на применении расширенных матриц Fke) и Нке), в то время как ковариационная матрица вычисляется на основе нерасширенных матриц Fk и Нк. Подобное несоответствие приводит к возникновению краевого эффекта на краях блоков

изображения (рис. 1б). Следовательно, требуется либо расчёт невязки в процессе вычислений, либо расширение матрицы Рк|к для вычислений вместе с матрицами Fk:e) и Нкв), что приведёт к увеличению объёма вычислений. При реализации алгоритма блочной фильтрации могут применяться различные подходы к решению данной проблемы.

Первый вариант реализации алгоритма. Для построения прогноза матрицы Рк|к используется нерасширенная матрица Fk размером ^ х з^. При этом необходим учёт использования расширенной матрицы сдвига FJ(e) для получения оценки хк+1|к . Для оценки

разницы между воздействием матриц Fke) и Fk используется вектор диагональных элементов матрицы ковариации Р^ Рк|к = [А(?к),р?^,-,р^^, разме-

'ь ,8Ь '

ром х1 и расширенный вектор рк|к размером х1, полученный из вектора Рк|к по формуле (9). В качестве константы с при расширении можно использовать среднее арифметическое вектора рк|к . Значение невязки выражается вектором qk:

.2

qk = ркв)р$ — Fk Рк|к| / Ц2

(14)

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

Рк+1|к = Fk Рк|к FT + Ок + О0к,

»Ок=| т =Ц;у

Коррекция матрицы ковариации осуществляется похожим образом. Формируется вектор диагональ-

ных элементов матрицы ковариации Рк|к —1

р = [ р(к|к—1) р(к|к—1) р(к|к—1)]Т рязмером '2 х1 и рк|к—1 = [р11 , р2 2 р,2г2 ] , размером х1 и

'ь ,6'ь

расширенный вектор рк|к—1 размером х1. Различие заключается в том, что невязка выражается вектором гк размером х1:

Гк = Н^к—1—Нк рк|к—т2

и учитывается в ходе коррекции матрицы ковариации в форме диагональной матрицы »Як , диагональными элементами которой являются элементы вектора гк:

(15)

Рк|к = Рк|к—1 Уки— Ук , Ук = Рк|к—1Н,

ик = НкУк + Як + БЯк ,

Í(k)

к П Г Г , ' =

(16)

0,7 ф у.

Второй вариант реализации алгоритма. Использование расширенных блоков при прогнозировании и

у

ъ

У

ъ

коррекции оценки можно интерпретировать как добавку некоторых шумов к оценке, прогнозируемой и уточняемой на основе нерасширенных матриц Fk и Нк:

-(к+1|к) _ к(е)п(Цк) _

Р,Я

= F(e)g(klk) = FXI

= *к & Р,Я = 1Х

Х (к|к) к Х Р,Я

+ 1

(к)

к™ Р,д ■

z

(к) = Н(е)„(к|к _1) = к

Р,Я

к &Р,Я (к)

= НХ

(к|к_1)

к Р,Я

+ Н к V

к

ку Р,Я '

где векторы те р Я и Vря содержат некоторые пиксели обрабатываемого изображения, а матрицы 1к и IIк имеют следующую структуру:

1к =1 Л(кк II;

/•(к) = Jr,t ~

Чм, г=( 1 _1) ь+г, t = а+п)Ь+г+т, г+т < 0илш + т > Ьили 1+п < 0или у+п >= Ь; 0, г = (]_1)Ь+г, t = (у+п)Ь + г + да, 0 < г+да < Ь, 0 < ]+п < Ь;

Н =11ЕII;

^ =

С, г = (1 _ 1)М+г, t = (1 + п)Ь + г + т

г + т < 0или г + т > Ь, 1 + п < 0или 1 + п > Ь; 0, г = (]_1)М + г, t = (7 + п)Ь+Г+да, 0 < 7 + да < Ь, 0 < + п < Ь.

Матрицы 1к и Нк включают в расчёт пиксели, находящиеся за пределами блока и принадлежащие областям перекрытия расширенных блоков. Таким образом, при вычислении матрицы ковариации можно учесть эти шумы, если вычислить их матрицы ко-вариации. Так как это пиксели одного и того же изображения, то вычисление матриц ковариации шумов осуществляется на основе матриц Ркк и Рк+1к. Для вычисления используются матрицы сдвига и децимации 1к и Н к.

В данном случае прогноз матрицы ковариации ошибки выглядит следующим образом:

Рк+1|к = 1к Рк|к 1кТ + 0к = 1к Рк|к 1кТ+1 Рк|к ^ + Ок. Коррекция осуществляется как в (16), но слагаемое БЯк заменено другим членом:

ик = Нк Ук + Я к =

= Н к Рк|к _1Нк + Н к Рк|к_1Нк + Як. Третий вариант реализации алгоритма. Альтернативой такому подходу является использование

расширенных матриц Нке) и 1ке) для обработки Рк|к

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

коррекции матрицы Рк|к и Рк|к_1 размером ¿Ь х ¿Ь

подменяются своими расширенными версиями - мат-

рицами Р^ и ^ размером ¿Ье х ¿Ье :

Р =1 |„(к|к )|| Рк|к = \Рг, ] ,

(е)

Р(е) = р(к|к )(е) рк|к ~ pr,t

Р(к|к)(е) = РгЛ _

Р(к|к) и г, ] '

г = (г + Ат _ 1)^Ье + Ат +1, ...,(г + Ат)^Ье _Ат, t = (] + Ат _ 1)^Ье + Ат +1, ...,(1 + Ат)^Ье _Ат; с, г ф (г + Ат _ 1)^Ье + Ат +1, ...,(г + Ат)^Ье _Ат, t ф (1 + Ат _ 1)^Ье + Ат +1, ...,(1 + Ат)^Ье _Ат и г = ^

0, г ф (г + Ат _ 1)^Ье + Ат +1, ...,(г + Ат)^Ье _Ат, t ф (1 + Ат _ 1)^Ье + Ат +1, ...,(1 + Ат)^Ье _Ат и г ф ^

где с - некоторая константа; г' = 1,...,; ] = 1,...,. В качестве с можно использовать среднее арифметическое вектора диагональных элементов Рк|к . На этапе

коррекции обрабатывается матрица Рк(|ек)_1 :

У^е) = Рке^НГ, ик = Нке)уке)+Я.

Р(е) = Р(е) _ щ(е)у(е)Т к к к к|к к|к_1 к к

щке) = Уке) и_\

Размеры матриц, участвующих в вычислениях:

¿Ье '

¿2

¿Ь1 х ¿Ь1 , щ(е) - ¿Ье х ¿Ь1 .

(е)

После чего матрицу Р^ следует преобразовать в Ркк. На этапе коррекции размер обработанного блока равен ¿1е х1.

„(к|к) = „(к|к_1) + щ(е) ( у(к) _ (к) )

где ) является вектором расширенного блока и

нуждается в преобразовании в Х (Рк,|Як) за счёт удаления

элементов, соответствующих области.

При прогнозировании матрица ковариации ошибки опять расширяется до Рк(|ек) и используется вместе

7(Ь).

с расширенной матрицей сдвига

Р = 1(е) Р(е)1(е)т + О

Рк+1|к = *к Рк|к + О.

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

Необходимо также рассмотреть ещё один вариант организации обработки. Применение алгоритмов кал-

к

мановского типа позволяет учитывать движение камеры при получении изображений непосредственно при задании матрицы децимации. В этом случае матрица сдвига считается единичной Рк=Г, к=1,...,К. Учёт смещений в матрице Нк осуществляется за счёт модификации матрицы 0к. Следует отметить, что в отличие от случая, когда Р^^Г, при котором сдвиги между изображениями являются смещениями между соседними к-м и £+1 -м кадрами, в данном случае смещения в матрице децимации являются смещениями обрабатываемого к -го кадра относительно некоторого опорного кадра из числа изображений НР в последовательности ук.

Анализ вычислительной эффективности

Применение предложенных алгоритмов является наиболее эффективным в случае, если процесс получения изображений является стационарным, что подразумевает единство характеристик смещения и ФРТ датчика для всех точек изображения. В этом случае для обработки всех блоков изображения можно использовать одну и ту же матрицу ковариации ошибки, обновляя её однократно при поступлении очередного наблюдения НР. Уравнения (11)-(13) в данном случае можно представить как:

Ч+1|k

= F,w XX(e)

k Ak|k ' ^k

Z = H(e) XX e

*-k ^k|k-1 :

h\k = Xk|k-i + Wk (Yk - Zk) ,

He)

(17)

где XXк|к , X££ , Ук и Zk - матрицы, каждый столбец

которых является вектором-столбцом блока соответствующего изображения:

Х^ = [Х(к1к) ^(к|к) ^(к|к) -(к|к) - (к|к) -(к|к)-,

XX (e) = rg(k|k-1) Ak|k-1 = [g1,1 ,

*-k|k-1 = [Jk)

(k|k-1) g(k|k-1) ,B,1 '•••'61,B '

gBk|B-1)

у _ [у(к) у(к) у(к) у(к) ] ук _ [у1,1 ,...,уВ,1,..., у1,В ,...,уВ,В] ,

Z _ г,(к) „(к) „(к) „(к) ]

Прогнозирование блоков изображения потребует порядка О (• В2) операций, коррекция - порядка

О (^ • • В2). Так как матрица ковариации корректируется однократно при поступлении наблюдения, то порядок вычислений при коррекции и прогнозе Ркк

составляет О (). Эти данные относятся к первому и

второму вариантам реализации алгоритма оптимальной фильтрации. В третьем варианте в процессе обработки Ркк используются расширенные матрицы,

следовательно, потребуется порядка О () операций. Кроме того, требуется проводить преобразования из Рк|к _1 в РккЦ и обратно, в результате чего работа данного варианта алгоритма занимает наибольшее время. Размер блока ВР 8Ъ не должен быть меньше т. Учитывая, что 8Ъе _ 8Ъ + 2Аш, сложность алгоритмов зависит от ширины ядра размытия ФРТ

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

В табл. 1 показано время в секундах, затраченное на обработку 16 изображений НР размером 256 х 256 пикселей для повышения разрешения в 4 раза. В таблице приведены значения для блоков изображения ВР (размеры блоков НР в 4 раза меньше). Оценки проводились при использовании ПК с процессором Intel® Core™ i5-2500 (3,3 ГГц), 8 ГБ ОЗУ в среде Matlab R2013a (64 бит).

Таблица 1. Оценки времени (в секундах) работы различных вариантов реализации алгоритма

Размер блоков ВР (sb х sb) Первый вариант Второй вариант Третий вариант

4 х 4 18,2 18,2 21,2

8 х 8 5,2 5,3 6,2

16 х16 2,0 2,0 2,4

32 х 32 3,0 4,2 4,5

По результатам испытаний видно, что наилучшее быстродействие демонстрирует первый вариант, второй вариант работает медленнее, так как требует дополнительных матричных операций на этапе прогнозирования и коррекции. Третий вариант, как и предполагалось, имеет самое низкое быстродействие за счёт операций с матрицами большего размера. При малых размерах блоков (4*4, 8*8) все алгоритмы показывают снижение быстродействия вследствие операций по расширению блоков - чем больше блоков, тем больше операций приходится затратить на расширение. С другой стороны, при малом количестве блоков и, следовательно, при большом размере блоков (32*32) возрастает размер матриц сдвига, децимации и ковариации ошибки, что также сказывается на быстродействии.

Результаты работы алгоритмов

Далее приводятся результаты испытаний описанных алгоритмов. На рис. 4 представлен график средней квадратичной ошибки между оригинальным и восстановленным изображениями, усреднённой по 7=10 реализациям для различных длин серии изображений НР.

Ошибка рассчитывается по следующей формуле:

1 7 I2 . ,2

^ 11 (^ _ хКг)) ,

г _1 г_1

-к

где хКг) - пиксель изображения хК с номером ¡; хК -

последнее изображение ВР из последовательности, сгенерированной в ходе г-й реализации эксперимента;

X К - восстановленное изображение; К=1,.. .,16.

В ходе каждой реализации генерируется случайное поле, из которого производится последовательность из К изображений НР. На их основе происходит восстановление оригинального изображения. Для каждого К проводится 10 реализаций эксперимента. К изображениям добавлен гауссовский шум с величиной среднеквадратичного отклонения (СКО) а^ _ 0,05 (область значения пикселя в эксперименте: [0, 1]).

0.025

0.02

0.015

и 0.01

й и ю

I 0.005

■ I гггт

Вариант 1 Вариант 2 Вариант 3

0^

1 4 7 10 13 16

Число изображений НР в серии

Рис. 4. График зависимости ошибки восстановления от длины серии изображений НР

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

На рис. 5 представлены результаты сравнения одного из предложенных вариантов построения алгоритмов (вариант 2) с известными алгоритмами сверхразрешения. Для сравнительного анализа были выбраны метод итеративных обратных проекций (ИОП)

[1] и алгоритм максимального правдоподобия (МП)

[2]. Первый из них основан на применении обратных проекций и является вариацией методологии, используемой в томографии. Второй алгоритм основан на методе максимального правдоподобия, являющемся широко распространённым способом оценки параметров статистической модели.

0.05 т

0.04

0.03

я 0.02

й и ю

I 001

0

ф. Калмана

ИОП

МП

0 0.05 0.1 0.15

Среднеквадратичное отклонение шума

Рис. 5. График зависимости ошибки восстановления алгоритмов от величины СКО шума

В данном случае происходит повышение разрешения в 4 раза по 16 изображениям НР и вычисляется средняя квадратичная ошибка между оригинальным и восстановленным изображениями, усредненная по 7=10 реализациям в условиях шума с различными

значениями СКО .

На рис. 5 видно, что фильтр Калмана в блочной реализации демонстрирует лучшее качество восстановления при повышении величины СКО шума.

На рис. 6 представлены результаты повышения разрешения в 4 раза с помощью трёх описанных ме-

тодов обработки изображений. Для повышения разрешения используются 16 изображений.

Рис. 6. Повышение разрешения в 4 раза: (а) пример изображения НР; (б) результат работы первого варианта построения алгоритма; (в) результат работы второго варианта; (г) результат работы третьего варианта; (д) оригинальное изображение ВР

На рис. 7 представлены результаты работы алгоритмов ИОП и МП. Для сравнения показан результат работы фильтра Калмана (второй вариант построения). При реализации метода ИОП используется 15 итераций.

Рис. 7. Повышение разрешения в 4 раза: (а) результат применения фильтра Калмана; (б) результат применения метода ИОП; (в) результат применения метода МП; (г) оригинальное изображение ВР

Как видно из рис. 7, результаты восстановления с использованием указанных алгоритмов визуально практически неразличимы. В табл. 2 показано время в секундах, потребовавшееся на обработку 16 изображений размером от 32*32 до 128*128 пикселей для повышения разрешения в 4 раза для разных методов построения сверхразрешения.

Таблица 2. Быстродействие различных алгоритмов

Размер блока Фильтр Калмана ИОП МП

32 х 32 0,16 0,63 0,33

64 х 64 0,27 0,87 1,29

128х128 0,82 2,01 4,96

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

Метод МП требует наибольшего времени вычисления из-за необходимости обработки больших прореженных матриц и решения системы уравнений методом сопряжённых градиентов. Метод ИОП показывает промежуточные результаты с точки зрения быстродействия.

На рис. 8 продемонстрированы результаты работы предложенных вариантов алгоритма калмановской фильтрации в условиях шума. Для повышения разрешения изображения используются 16 изображений НР с применением гауссовского шума (о® = 0,05).

Рис. 8. Повышение разрешения в 4 раза в условиях шума: (а) пример изображения НР; (б) результат работы первого варианта построения алгоритма; (в) результат работы второго варианта; (г) результат работы третьего варианта

На рис. 8 видно, что в условиях шума предложенные варианты дают незначительно различающиеся результаты.

Заключение

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

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

Благодарности

Работа выполнена при поддержке гранта Минобр-науки РФ по программе «Развитие кооперации российских вузов и производственных предприятий» (Постановление Правительства № 218 от 09.04.2010 г. - 3 очередь, № 02.G25.31.0002).

Литература

1. Peleg, S. Improving resolution by image registration / S. Peleg, M. Irani // CVGIP: Graphical Models and Image Processing. -1991. - V. 53, N 3. - P. 231-239. - ISSN: 1049-9652.

2. Elad, M. Restoration of single super-resolution image from several blurred, noisy and down-sampled measured images / M. Elad, A. Feuer // Image Processing, IEEE Transactions on. - 1997. - V. 6, N 12. - P. 1646-1658. - ISSN 10577149.

3. Elad, M. Super-resolution reconstruction of continuous image sequences / M. Elad, A. Feuer // International Conference on Image Processing (ICIP 99). - 1999. - V. 3. - P. 459-463.

4. Elad, M. Super-resolution reconstruction of continuous image sequences: adaptive filtering approach / M. Elad, A. Feuer // Image Processing, IEEE Transactions on. - 1999. -V. 8, N. 3. - P. 387-395. - ISSN 1057-7149.

5. Newland, C. Time invariant steady-state Kalman filter for image super-resolution / C. Newland, D. Gray // Image and Vision Computing New Zealand 2005 (IVCNZ05), 2006. - P. 381-387.

6. Newland, C. Modified Kalman filtering for image super-resolution / C. Newland, D. Gray, D. Gibbins // Image and Vision Computing New Zealand 2006 (IVCNZ06), 2006. - P. 79-84.

7. Park, P. New square-root algorithms for Kalman filtering / P. Park, T. Kailath // IEEE Trans. Automat. Control. - 1995. - V. 40. - № 5. - P. 895-899. - ISSN 0018-9286.

8. Куликова, М.В. О скаляризованном вычислении функции правдоподобия в квадратно-корневых матричных алгоритмах фильтрации / М.В. Куликова // Автоматика и телемеханика. - 2009. - № 5. - С. 122-139. - ISSN 0005-2310.

9. Azimi-Sadjadi, M.R. Two-Dimensional Block Kalman Filtering For Image Restoration / M.R. Azimi-Sadjadi, P.W. Wong // IEEE Trans. on Acoustics, Speech and Signal Processing. -1987. - V. ASSP-35. - P. 1736-1749. - ISSN 0096-3518.

10. Patti, A.J. A new motion-compensated reduced-order model kalman filter for space-varying restoration of progressive and interlaced video / A.J. Patti, A.M. Tekalp, M.I. Sezan // Image Processing, IEEE Transactions on. - 1998. - V. 7, N. 4. - P. 543-554. - ISSN 1057-7149.

11. Васильев, К.К. Статистический анализ многомерных изображений / К.К. Васильев. - Ульяновск: УлГТУ, 2002. - 156 с.

References

1. Peleg, S. Improving Resolution by Image Registration / Pe-leg, S. Improving resolution by image registration / S. Peleg, M. Irani // CVGIP: Graphical Models and Image Processing. - 1991. - V. 53, N 3. - P. 231-239. - ISSN: 10499652.

2. Elad, M. Restoration of single super-resolution image from several blurred, noisy and down-sampled measured images / M. Elad, A. Feuer // Image Processing, IEEE Transactions

on. - 1997. - V. 6, N 12. - P. 1646-1658. - ISSN 10577149.

3. Elad, M. Super-resolution reconstruction of continuous image sequences / M. Elad, A. Feuer // International Conference on Image Processing (ICIP 99). - 1999. - V. 3. -P. 459-463.

4. Elad, M. Super-resolution reconstruction of continuous image sequences: adaptive filtering approach / M. Elad, A. Feuer // Image Processing, IEEE Transactions on. - 1999. -V. 8, N. 3. - P. 387-395. - ISSN 1057-7149.

5. Newland, C. Time invariant steady-state Kalman filter for image super-resolution / C. Newland, D. Gray // Image and Vision Computing New Zealand 2005 (IVCNZ05), 2006. - P. 381-387.

6. Newland, C. Modified Kalman filtering for image super-resolution / C. Newland, D. Gray, D. Gibbins // Image and Vision Computing New Zealand 2006 (IVCNZ06), 2006. - P. 79-84.

7. Park, P. New square-root algorithms for Kalman filtering / P. Park, T. Kailath // IEEE Trans. Automat. Control. - 1995. - V. 40. - № 5. - P. 895-899. - ISSN 0018-9286.

8. Kulikova, M.V. About scalarized calculating the likelihood function in the square-root matrix filtering algorithms / M.V. Kulikova // Automation and Remote Control. - 2009.

- V 5. - P. 122-139. - ISSN 0005-2310. - (In Russian).

9. Azimi-Sadjadi, M.R. Two-Dimensional Block Kalman Filtering For Image Restoration / M.R. Azimi-Sadjadi, P.W. Wong // IEEE Trans. on Acoustics, Speech and Signal Processing. - 1987. - V. ASSP-35. - P. 1736-1749. -ISSN 0096-3518.

10. Patti, A.J. A new motion-compensated reduced-order model kalman filter for space-varying restoration of progressive and interlaced video / A.J. Patti, A.M. Tekalp, M.I. Sezan // Image Processing, IEEE Transactions on. - 1998. - V. 7, N. 4. - P. 543-554. - ISSN 1057-7149.

11. Vasiljev, K.K Statistical analysis of multidimensional images / K.K. Vasiljev. - Ulyanovsk: UlSTU Publisher, 2002.

- (In Russian).

BLOCK ALGORITHMS OF IMAGE PROCESSING BASED ON KALMAN FILTER FOR SUPERRESOLUTION RECONSTRUCTION

A.A. Sirota, A.Y. Ivankov Voronezh State University Abstract

Different approaches to block image processing applied to the reconstruction of super-resolution images using the Kalman filter are analyzed. The results of the various approaches are represented. Key words: block based image processing, Kalman filter, super resolution.

Сведения об авторах Сирота Александр Анатольевич, 1954 года рождения, в 1976 году окончил Воронежский государственный университет по специальности «Радиофизика и электроника». Доктор технических наук (1995 год), профессор, заведует кафедрой технологий обработки и защиты информации Воронежского государственного университета. Область научных интересов: синтез и анализ систем сбора и обработки информации, методы и технологии компьютерного моделирования информационных процессов и систем, системный анализ в сфере информационной безопасности, компьютерная обработка изображений, нейронные сети и нейросетевые технологии в системах принятия решений. E-mail: sir@cs.vsu.ru .

https://sites.goodexoiTi/a/sc.vsu.m/tozi/sotmdniki-kafedrv/sotrudniki-kafedrv-tozi/sirota. Alexander Anatolievich Sirota (1954) graduated from Voronezh State University in 1976 majoring in Radiophysics and Electronics. Professor, Doctor of Technical Sciences (since 1995). Currently head of Information Processing and Security Technologies Chair at Voronezh State University. Research interests: analysis and design of information collection and processing systems, methods and techniques of information processes and systems computer modeling, system analysis in information security, digital image processing, neural networks and neural network technologies in decision-making systems.

Иванков Александр Юрьевич, 1990 года рождения, в 2012 году окончил Воронежский государственный университет по специальности 230400 «Информационные системы и технологии». Аспирант кафедры технологий обработки и защиты информации Воронежского государственного университета. Область научных интересов: обработка графических изображений, программирование.

E-mail: ivankovay@gmail.com .

Alexander Yurievich Ivankov (1990) graduated from Voronezh State University in 2012, majoring in Information Systems and Technologies. Currently a postgraduate student of Information Processing and Security Technologies chair at Voronezh State University. Research interests: computer graphics processing, programming.

Поступила в редакцию 3 октября 2013 г.