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

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

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

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

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

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

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

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

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

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

x к + 1 = F k x k ⁺ w k , ⁽¹⁾

где x k+ 1 - изображение ВР, полученное из предыдущего изображения x k за счёт перемещения камеры и/или объекта в процессе получения изображений; F k – матрица сдвига размером L 2 х L ², характеризующая геометрические взаимные деформации (смещения) изображений; w _k – гауссовский шум с нулевым математическим ожиданием и ковариационной матрицей Q _k = g ^{q )2} I размером L 2 х L ².

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

У k = Hk xk + vk, (2) где yk - изображение с низким разрешением; xk -изображение с высоким разрешением, из которого производится yk ; Hk – матрица децимации размером M2 х L2, которая осуществляет преобразование xk в yk ; vk – вектор аддитивного гауссовского шу- ма с нулевым математическим ожиданием и матрицей ковариации Rk = Ск)2I размером M2 хM2.

Следует отметить, что матрицы H k , F k , R k и Q k , определяющие систему пространства состояний, считаются известными. В данной работе предполагается глобальный характер смещений между соседними изображениями, подразумевающий одинаковую величину смещений для всех пикселей одного изображения, что соответствует перемещению камеры относительно сцены.

Матрица сдвига F _k имеет структуру, обусловленную ядром фильтра для интерполяции изображения П _k = ЦП- ||, в результате свертки с которым исходного изображения получается смещённое изображение:

^F k =1 «1г

n

( k ) p , q ,

f ( k ) = J fr , t

r = ( j-1) L + i, t = (j + n) L + i + m, 0 < i + m < L, 0 < j + n < L;

0, r ^(j-1)L + i, t ^(j + n)L + i + m, 0 < i + m < L, 0 < j + n < L, где i = 1,...,L; j = 1,...,L ; p = 1,...,Nn ; q = 1,...,Nn ; m = p -(Nn+1) /2; n = q -(Nn+ 3) /2; матрица фильтра nk имеет размер Nпх Nn, Nn - нечётное число.

Матрица децимации H _k имеет похожую структуру. В случае если L / M = ц и ц - целое число, то:

H k =| ^ rk ’||;

⁶ k =16*^11;

f _А( k ) ⁶ p , q ’

U ( k ) _J "r , t

r = ( j-1) M + i, t = (j + n) L + i + m, 0 < i + m < L, 0 < j + n < L;

0, r ^( j -1) M + i, t ^ (j + n) L + i + m, 0 < i + m < L, 0 < j + n < L,

где i = 1,...,M ; j = 1,...,M ; p = 1,...,N6; q = 1,...,N6; m = p-(N6+1)/2; n = q-(N6+ 3)/2; i = ^(i-0.5)-ц"|; j = P( j - 0.5)-ц|. Оператор |"-"| означает округление до ближайшего целого в большую сторону. Ядро фильтра 6k размером N6 х N6 (N6 - нечётное число) характеризует функцию рассеяния точки (ФРТ) датчика камеры. В θk также следует учитывать весовые коэффициенты интерполяции соседних пикселей в случае, если ц не является целым.

Фильтр Калмана позволяет по полученной совокупности наблюдений y ^k = { y ₁,..., y _k } за k шагов получить оптимальные в среднеквадратичном оценки x _k | _k - 1 , x _k | _k вектора x _k . Оценка x _k k является изображением СР. Это оценка состояния объекта вида «точка в точке», полученная на основе апостериорной плотности распределения P( x _k | y ^k ). Оценка x _k k - 1 -это прогноз состояния объекта на один шаг вперёд, в данном случае – предсказание следующего изображения ВР. Основой получения прогноза является апостериорная плотность распределения P( x _k | _k - 1 | y ^{k - 1}). Переход от прогноза x _k k - 1 к оценке x ɶ _{k | k} осуществляется на основе матрицы ковариации ошибки P _k | _k - 1 = M [e _k £ T J , где ошибка £ _k = x _k - x _k k - 1 характеризует разницу между изображением ВР x _k и предсказанным изображением ВР x _k k - 1 .

В классической форме фильтра Калмана при поступлении очередного изображения НР y _k производится коррекция оценки x ɶ _{k | k} :

^x k | k = ^x k | k - 1 ⁺ W k ( ^y k ^{- H} k ^x k | k - 1 ) ,                  ⁽⁵⁾

W k = V k U k ¹, U k = H k P k | k - 1 H T + R k ,

V_k = P kk - 1 H T , P kk = P kk - 1 - W k V T .           (6)

Рекуррентное выражение для оценки xk+1k на k + 1-м шаге имеет вид xk+1k = Fxkik = Fxkik-1 + F Wk ( уk - Hkxkik-1). (7)

Рекуррентное выражение для матрицы P _{k + 1} k на k + 1-м шаге имеет вид

P k + 1 k = F k P kk F k ^T + Q k .                              (8)

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

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

Изображения делятся на B блоков по вертикали и горизонтали. Размер блока НР sbl х sbl (sbl • B = M), размер блока ВР sb = sbl -ц (sb • B = L). Из-за ограниченности размеров блоков изображения при использовании матрицы сдвига на краях блока смещённого изображения проявляется осцилляция яркости, что проиллюстрировано на рис. 1а. Аналогичный эффект возникает и при прореживании изображения.

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

Для преодоления данного эффекта требуется сшивка блоков изображений. Достичь этого можно за счёт расширения на k -м шаге блоков изображения ВР до размера s be х s be , где s be = s b + 2 A m , A m = ( N - 1) / 2 , N - ширина ядра фильтра, равная

N ₉ или N п . Для смещения блока используется расширенная матрица сдвига F ke ) размером s b х s be , а для прореживания блока применяется расширенная матрица H e ) размером s b х s be . Оценка x kk также расширяется:

X . =| H k ’|;

' ^ ( k | k )

^x ( j - 1 -A m ) L + i -A m

, i = A m + 1,

A m + L ,

x

(k) I i, j = 1

c ,

j = A m + 1,..., A m + L , i = 1,..., A m , j = 1,..., A m , и i = L + A m + 1,..., L + 2 A m ,

j = L + Am +1,..., L + 2A m, где с – некоторая константа.

Из расширенного изображения-оценки X _k извлекаются блоки G ^k размером s_be х s_be

G’ Pk =| Vq? ^{,( k l} k )l ^;

p , q ,( k | k ) = ( k | k )

® i , ^{j             ( p - 1)} s b ⁺ i ^{,( q - 1)} s b ^{+ j} , .

где i = 1,..., s be , j = 1,..., s be , p = 1,..., B , q = 1,..., B .

Каждый блок разворачивается в вектор g(pk,)q , по- сле чего к нему применяется операция сдвига:

( k + 1| k )       ( e ) ( k|k )

^x p , q k g p , q

где вектор 3x⁽ _p k + ^{1| k} )

является блоком оценки x _{k +} i| _k .

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

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

Рис. 2. Процесс блочного прогнозирования изображения

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

,( k )        ( e )„( k\k - 1)

z p , q      k g p , q    ’

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

sbe sbl

расширенное x _k|k

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

Расширенные матрицы Fke) и He) получаются за счёт добавления новых столбцов на места, соответствующие новым элементам в расширенном блоке g(kqk) или g^k+^k). Таким образом, значения пикселей на границах блоков x(k+*\к) и z(„k) вычисляются с помощью Р,q           Р,q элементов, содержащихся в перекрытиях блоков.

Формула (5) для коррекции блока имеет вид:

х( k\k )_=i( k\k -D_+W, (_v( ^k ) __z( ^k )

x p , q     x p , q ⁺ W k ⁽ y p , q z p , q ) ’                   ⁽¹³⁾

где x(kq)- развёрнутый в вектор блок изображения- оценки xk\k; x(kq 1 - вектор блока изображения-прогноза xk\k-1; y (kq - вектор блока обрабатываемо- го изображения НР yk ; z(pk,)q - вектор блока изобра- жения НР Zk ; p = 1,...,B , q = 1,...,B - индексы блока. Размер матрицы усиления Wk соответствует размеру блока и равен sb х sb.

В процессе прогнозирования и коррекции матрицы ковариации ошибки необходимо учесть, что расчёт блока оценки основан на применении расширенных матриц Fk(e) и H(ke) , в то время как ковариацион- ная матрица вычисляется на основе нерасширенных матриц Fk и Hk. Подобное несоответствие приводит к возникновению краевого эффекта на краях блоков изображения (рис. 1 б). Следовательно, требуется либо расчёт невязки в процессе вычислений, либо расширение матрицы Pk\k для вычислений вместе с матрицами Fke) и He), что приведёт к увеличению объёма вычислений. При реализации алгоритма блочной фильтрации могут применяться различные подходы к решению данной проблемы.

Первый вариант реализации алгоритма. Для построения прогноза матрицы Pk\k используется нерасширенная матрица Fk размером sb2 х sb2 . При этом необходим учёт использования расширенной матрицы сдвига Fke) для получения оценки xk+1\k. Для оценки разницы между воздействием матриц Fk(e) и Fk используется вектор диагональных элементов матрицы (k\k) r>(k\k)       r>(k\k) IT ковариации Pk\k p k\k = [ pn , p 22 ,...,p 2 2J , разме-’             ’                sb , sb ром sb2 х1 и расширенный вектор p(ke\k) размером sb2e х1 , полученный из вектора pk\k по формуле (9). В качестве константы c при расширении можно использовать среднее арифметическое вектора pk\k . Значение невязки выражается вектором qk :

q k = |^F k e ^kek ^- F k p k\k\ / о².                         ⁽¹⁴⁾

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

P k ₊₁ k = F k P kk F k + Q k + DQ k , ,,          ..                    L( k )          •

DQ k = dq ⁽ , ^кД, dqi, =       , =. j ;

"         "                  I 0, 1 * 1 .

Коррекция матрицы ковариации осуществляется похожим образом. Формируется вектор диагональных элементов матрицы ковариации    P _{k \ k - 1} ,

( k\k - 1) ( k\k - 1)        ( k\k - 1)11                     2

p k\k - 1 = [ ^p 11     , ^p 2 2    ,..., ^p -2 Л J , размером s b ^х1 и

,                   ,                      sb , sb расширенный вектор p(ke\k) -1 размером sb2e х1 . Различие заключается в том, что невязка выражается вектором rk размером sb х 1:

r k = |H ke ) p kek - 1 - H k p k\k -1| / ц ²                       (15)

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

P kk = P kk - 1 - V k U - V , V k = P kk - 1 H k ^T,

U k = H kVk + R k + DR k ,                    (16)

( k )

^DR k =1 d' S’ll- d-^ k *=p₀, ; ,;.^J;

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

~( k + 1| k ) _{= F}( e ) ( k | k ) _{= F} ~( k | k ) _F a( k )

x p , q        F k g p , q     F k x p , q ⁺ F k w p , q ,

,(k) =H(■)„(kk-l)=H =(k|k-1) „ -k z p, q H k g p, q      H k x p, q   + H k v p, q , где векторы w(,,q и vp,q содержат некоторые пиксели обрабатываемого изображения, а матрицы 1^k и Hk имеют следующую структуру:

F . =|| f , k '||;

_ ( k ) ₌

Jr , t

^, r = ( j-1) L+i, t = (j'+n)L+i+m,

i + m <  0или i + m >  L или j + n <  0или j + n >= L ;

0, r = ( j - 1 ) L + i , t = ( j’ + n ) L + i + m ,

0 <  i + m <  L , 0 <  j + n <  L ;

h k =1 | h !?| ;

0(pkq, r = (j -1) M+i, t = (j + n)L + i + m, i + m < 0 или i + m > L, j + n < 0 или j + n > L;

0, r = ( j - 1 ) M + i , t = ( j + n ) L + i + m ,

0 <  i + m <  L , 0 <  j + n <  L .

^_        _^_

hi,^k, ) =

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

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

P k + i| k = F k P kk F T + Q k = F k P k | k F k + F k P kk F T + Q k .

Коррекция осуществляется как в (16), но слагаемое DR k заменено другим членом:

U k = H k V k + R k =

= H k P kk - 1 H T + H k P k | k - 1 H T + R k .

Третий вариант реализации алгоритма . Альтернативой такому подходу является использование расширенных матриц H ⁽ _k ^{e )} и F _k ^{( e )} для обработки P _{k | k} и расширение самой матрицы ковариации ошибки во время данных операций. То есть на этапах прогноза и коррекции матрицы P _k | _k и P _k | _k , размером s b х s b подменяются своими расширенными версиями - матрицами P^ и Р^ размером s be ^х s be ^:

'X=1 pjk I p (e) _  „(k|k)(e)

Pk|k = pr, t pj), r = (i + Am - 1) sbe + Am +1,

С ,

„( k | k )( e ) _ J p r , t

0,

...,( i + A m ) s_be - A m , t = ( j + A m - 1) s_be +A m + 1, ...,( j + A m ) s_be - A m ;

r ^ ( i + A m - 1) s_be + A m + 1, ...,( i + A m ) s_be - A m , t ^ ( j + A m - 1) s_be + A m + 1, ...,( j + A m ) s_be -A m и r = t ;

r ^ (i + Am -1) sbe + Am +1, ...,(i + Am)sbe - Am, t ^ (j + Am -1) sbe + Am +1, ...,(j + Am)sbe-Am и r ^ t, где c - некоторая константа; i = 1,

...,

s b ; j = 1

, ..., s_b

. В

качестве c можно использовать среднее арифметическое вектора диагональных элементов P _k | _k . На этапе коррекции обрабатывается матрица P _k ⁽ _| ^e _k ⁾ _{- 1} :

( e )       ( e )      ( e )T                 ( e ) ( e )

V k = ^P k | k - 1 ^H k , ^U k = ^H k V k ^{+ R .}

( e )        ( e ) - 1      ( e )       ( e )           ( e ) ( e )T

W k = ^V k ^U k , ^P k | k = ^P k | k - 1 W k V k .

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

V^{( e} )   v² xv²  ТТ    v² xv²  W^{( e} )   v² xv²

^V k      s be ^X s bl , ^U k    s bl ^X s bl , W k      s be ^X s bl .

После чего матрицу P ^ ek ) следует преобразовать в P k\k . На этапе коррекции размер обработанного блока равен s_b ² _e X 1 .

( k | k )       ( k | k - 1)        ( e ) ( k )      ( k )

g p, q     g p, q +T 'k ( y p, q z p, q ) , где g^kqk) является вектором расширенного блока и нуждается в преобразовании в x(kqk) за счёт удаления элементов, соответствующих области.

При прогнозировании матрица ковариации ошибки опять расширяется до P _k ⁽ _| ^e _k ⁾ и используется вместе с расширенной матрицей сдвига F _k ^{( b )} :

( e ) ( e ) ( e )T

^P k + l|k = F k ^P k | k F k ^{+ Q} .

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

Необходимо также рассмотреть ещё один вариант организации обработки. Применение алгоритмов кал- мановского типа позволяет учитывать движение камеры при получении изображений непосредственно при задании матрицы децимации. В этом случае матрица сдвига считается единичной Fk=I, k=1,…,K. Учёт смещений в матрице Hk осуществляется за счёт модификации матрицы θk. Следует отметить, что в отличие от случая, когда Fk≠I, при котором сдвиги между изображениями являются смещениями между соседними k-м и k+1 -м кадрами, в данном случае смещения в матрице децимации являются смещениями обрабатываемого k -го кадра относительно некоторого опорного кадра из числа изображений НР в последовательности yk.

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

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

_ɶ                ( e ) _ɶ ( e )                ( e ) _ɶ ( e )

X k + 1| k = F k X k | k , Z k ^- H k X k | k - 1 ,

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

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

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

Размер блоков ВР ( s b X s b )	Первый вариант	Второй вариант	Третий вариант
4 х 4	18,2	18,2	21,2
8 х 8	5,2	5,3	6,2
16 X 16	2,0	2,0	2,4
32 х 32	3,0	4,2	4,5

X kk - X kk - 1 + W k ( Y k - Z k ) ,

где X ^ɶ _{k | k} , X ^{ɶ (} _k ^e _{| k} ⁾ , Y _k и Z _k - матрицы, каждый столбец которых является вектором-столбцом блока соответствующего изображения:

ɶ             ( k | k )

X k | k ^{- [ x} 1,1

,...,

( k | k )    ( k | k )

x B ,1 , x 1,2

,...

( k | k ) , x B ,2 ,

...,

( k | k ) x 1, B ,

...,

( k | k ) x B , B ] ,

V1 e ) -     k | k -1)

X k | k - 1 ^{- [ g} 1,1      ,

...

i k | k - 1)

, g B ,1 ,

...

i k | k - 1)

, g 1, B     ,

...

, g BB - ¹ ) ],

Y_k = [ У (1 Z k = [zu )

,...,

,...,

( k ) y B ,1

( k ) z B ,1

,...,

,...,

( k ) y 1, B

( k ) z 1, B

,...,

,...,

( k ) , y _{B , B} ] ,

( k ) z B , B ] .

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

O (sb • sb тируется

• B ² ) . Так как матрица ковариации коррек-однократно при поступлении наблюдения,

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

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

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

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

то порядок вычислений при коррекции и прогнозе P k\k составляет O ( s b ) . Эти данные относятся к первому и второму вариантам реализации алгоритма оптимальной фильтрации. В третьем варианте в процессе обработки P k\k используются расширенные матрицы, следовательно, потребуется порядка O ( s be ) операций. Кроме того, требуется проводить преобразования из P _k | _{k - 1} в P k ^ и обратно, в результате чего работа данного варианта алгоритма занимает наибольшее время. Размер блока ВР s_b не должен быть меньше Ц . Учитывая, что s be - s b + 2 A m , сложность алгоритмов зависит от ширины ядра размытия ФРТ

T L ²

=K -^7ZZ(xK't

TL t - 1 i - 1

x K ) ,

где x ⁱ _K ^{( t )} – пиксель изображения x ⁽ _K ^{t )} с номером i ; x ⁽ _K ^{t )} – последнее изображение ВР из последовательности, сгенерированной в ходе t -й реализации эксперимента; x ɶ ⁽ _K ^{t )} – восстановленное изображение; K =1,…,16.

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

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

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

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

На рис. 5 представлены результаты сравнения одного из предложенных вариантов построения алгоритмов (вариант 2) с известными алгоритмами сверхразрешения. Для сравнительного анализа были выбраны метод итеративных обратных проекций (ИОП) [1] и алгоритм максимального правдоподобия (МП) [2]. Первый из них основан на применении обратных проекций и является вариацией методологии, используемой в томографии. Второй алгоритм основан на методе максимального правдоподобия, являющемся широко распространённым способом оценки параметров статистической модели.

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

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

В данном случае происходит повышение разрешения в 4 раза по 16 изображениям НР и вычисляется средняя квадратичная ошибка между оригинальным и восстановленным изображениями, усредненная по T =10 реализациям в условиях шума с различными значениями СКО σ ⁽ _k ^{R )} .

На рис. 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 изображений НР с применением гауссовского шума ( o kR ⁾ = 0,05 ).

б)

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

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

Заключение

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

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