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

УДК 621.391.837.1:004.932.4:519.712.2

М.А. ПОПОВ*, С.А. СТАНКЕВИЧ*, С.В. ШКЛЯР*

АЛГОРИТМ ПОВЫШЕНИЯ РАЗРЕШЕНИЯ СУБПИКСЕЛЬНО СМЕЩЕННЫХ ИЗОБРАЖЕНИЙ

Научный Центр аэрокосмических исследований Земли Института геологических наук НАН Украины, Киев, Украина

Киевский национальный университет имени Тараса Шевченко, Киев, Украина

Анотація. Викладено математичну модель і алгоритм відновлення зображення підвищеної розрізненості шляхом спільної обробки вхідних зображень низької розрізненості, субпіксельно зміщених одно відносно одного. Досягається підвищення розрізненості, близьке до формального, в результаті передискретизації.

Ключові слова: підвищення просторової розрізненості, субпіксельно зміщені зображення, переди-скретизація, регуляризація, згортка в частотній області.

Аннотация. Изложены математическая модель и алгоритм восстановления изображения повышенного разрешения путём совместной обработки исходных изображений низкого разрешения, субпиксельно смещённых относительно друг друга. Достигается повышение разрешения, близкое к формальному, в результате передискретизации.

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

Abstract. Mathematical model and algorithm of enhanced resolution image reconstruction by joint processing of input low-resolution of subpixel displaced images one relative to other are presented. Resolution enhancement, which is close to the formal one as oversampling result is achieved.

Keywords: spatial resolution enhancement, subpixel displaced images, oversampling, regularization, frequency domain convolution.

1. Введение

Основным средством регистрации изображений в современных оптико-электронных системах являются многоэлементные твердотельные полупроводниковые фотоприёмники. Во многих отраслях народного хозяйства, науки и техники, таких как системы слежения и безопасности, навигации и управления, дистанционного зондирования Земли и астрономических объектов, основным критическим показателем качества цифровых изображений выступает обеспечиваемое пространственное разрешение [1]. Существенный прогресс в полупроводниковых технологиях, в частности, в производстве многоэлементных - линей-ковых и матричных - фотоприёмников позволяет обеспечивать необходимое пространственное разрешение в большинстве приложений путём увеличения количества фотоприёмных элементов (фотодетекторов) в изображающей матрице [2]. Однако в некоторых случаях требуемое разрешение выходит за пределы возможностей существующих технологий изготовления матричных фотоприёмников. К таким областям относятся астрономическая и детальная аэрокосмическая съёмка, предельная микроскопия, прецизионная видеоспектрометрия и др. Размеры фотодетекторов могут уменьшаться только до определённого, фундаментально физически ограниченного предела, определяемого длиной волны регистрируемого излучения. Кроме того, увеличение количества фотодетекторов в фотоприёмной матрице приводит к значительному усложнению конструкции и производства и резко повышает стоимость съёмочной аппаратуры. Выход из такого положения состоит в последовательной регистрации нескольких субпиксельно смещённых изображений и формировании из них единого изображения повышенного разрешения [3]. В зарубежной научной литературе данный метод получил название “сверхразрешение” (superresolution) [4].

© Попов М.А., Станкевич С.А., Шкляр С.В., 2015

ISSN 1028-9763. Математичні машини і системи, 2015, № 1

29

Сверхразрешение заключается в итерационном получении двух или нескольких изображений низкого разрешения одной и той же сцены, смещённых относительно друг друга на определённую долю геометрического размера пиксела. Для восстановления значений изображения в дискретных частях пикселов применяется программная обработка набора исходных изображений низкого разрешения [5]. Ниже будет рассмотрен алгоритм повышения разрешения двух матричных изображений, смещённых на произвольные доли пиксела одно относительно другого как по горизонтали, так и по вертикали. Алгоритм не претерпевает принципиальных изменений и в случае использования нескольких субпиксельно смещённых изображений [6].

2. Состояние вопроса

Достаточно полный обзор известных методов субпиксельного повышения разрешения цифровых изображений сделан в работах [7] и [8].

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

Основные этапы обработки:

1) определение смещения между изображениями и субпиксельная привязка изображений к некоторой референтной субпиксельной сетке;

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

3) обработка передискретизированного изображения: интерполяция и повышение чёткости.

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

Не все методы повышения разрешения изображения строго удовлетворяют этой схеме. Например, в методе maximum a posteriori (MAP) чёткость повышается уже на этапе создания передискретизированного изображения. Иногда используется привязка исходных изображений низкого изображения к передискретизированному изображению. При этом алгоритм получается итерационным. Существуют методы, требующие на начальном этапе только грубого (с точностью до пиксела) совмещения исходных изображений. Точная оценка смещения получается одновременно с повышением разрешения. Это методы, работающие в частотной области, и методы, сводящие задачу повышения разрешения к задаче total least squares (TLS) или structured total least squares (STLS) [10].

В зависимости от того, используется ли при обработке изображения преобразование Фурье, методы делят на работающие в пространственной области (spatial domain methods) и работающие в частотной области (frequency domain methods). Такие действия, как оценка смещения одного изображения относительно другого, могут быть произведены (в случае “поступательного” смещения или в более общем случае аффинного преобразования) в частотной области. Повышение чёткости и устранение артефактов (на частоте Найквиста) также могут быть произведены в частотной области. Обычно считается, что проводить обработку в частотной области можно только в случае поступательных смещений. Хотя в [11] предложен итерационный метод error-energy reduction (EER), позволяющий учитывать любые смещения изображений (предполагается, что сами смещения были найдены раньше) и на каждой итерации использующий прямое и обратное преобразования Фурье. Полностью алгоритм, включающий в себя оценку поступательного смещения по методу максимизации корреляции, работающему в частотной области, и повышение разрешения по методу ERR, приведен в [12].

30

ISSN 1028-9763. Математичні машини і системи, 2015, № 1

Методы projecting onto convex sets (POCS) [13] и MAP [14], работающие в пространственной области, используют регуляризацию и поэтому требуют некоторого предположения об изображении высокого разрешения.

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

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

В методе MAP уравнения заменены на приближенные равенства и ищется решение задачи следующего вида [15]:

где первое слагаемое отвечает за регуляризацию (в вероятностном подходе - за априорные предположения об изображении повышенного разрешения), а второе слагаемое - сумма квадратов невязок - в вероятностном подходе соответствует правдоподобию. Решение задачи минимизации находится в явном виде. Не обязательно использовать евклидову норму; например, в работе [16] была использована А1-норма.

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

Пространственная частота дискретизации увеличивается в два раза по каждому из направлений. Количество пикселов в передискретизированном изображении будет в 4 раза больше, чем в исходном изображении.

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

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

где x - развернутое в вектор изображение повышенного разрешения, у - развернутые в

один вектор два исходных изображения, A - матрица линейного оператора. Система (1) состоит из уравнений вида

j| x( x у) s( x - х0г, у - у0г) dxdy = у.

к=1

3. Метод

Ax » y, Ax = у + e,

(1)

(2)

a1 xij + a 2 x

'2й- i+1, j

+ a x.

'3xi,j+1 + a4xi+1,j+1 ~ yiJ ,

(3)

ISSN 1028-9763. Математичні машини і системи, 2015, № 1

31

где a, ..a4 обозначают относительные размеры площади пикселов в изображении повышенного разрешения, a +a2 + a3 + a4 = 1. В каждой строке матрицы A четыре ненулевых элемента, а именно а1, ..., a4.

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

X = E-1A(AT E-1A + Ee)-1(y - A [і x) + і ^ (4)

где mx и Sx - соответственно априорные среднее и ковариационная матрица вектора x, Se - ковариационная матрица вектора ошибок.

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

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

4. Преобразование в окне

Так как обработка всего изображения описанным способом очень ресурсоёмка, то применялась обработка изображения в окне. Для каждого пиксела изображения повышенного разрешения строилось свое окно размером не менее 23*23 так, чтоб пиксел находился в центре окна (конечно, если такое возможно, то есть, если пиксел находится достаточно далеко от границы изображения). Проводилось повышение разрешения изображения в окне, у изображения повышенного разрешения вычислялся и использовался только центральный пиксел.

Целью использования преобразования в окне является уменьшение вычислительных затрат. Если обрабатывается все изображение размером m X n без использования окна, то количество операций имеет порядок O(m3 Xn3). При использовании окон размером dX d количество операций имеет порядок O(d6 + mnd2) .

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

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

5. Алгоритм

Алгоритм повышения разрешения использует следующие входные данные: два изображения низкого разрешения, смещённые одно относительно другого; величины смещений второго изображения относительно первого; шумовое изображение (например, полученное как разность двух изображений, полученных со смещением на целое число пикселов и потом совмещённых).

32

ISSN 1028-9763. Математичні машини і системи, 2015, № 1

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

Шаг 1. Вычисление автоковариационных матриц исходного и шумового изображений - оценка матриц Sx и Se.

Шаг 2. Регуляризированное обращение матрицы A по формуле

B = S-1A(AT S-1A + Хе)-1.

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

Шаг 3. Вычисление оценок значения каждого пиксела высокого разрешения. По сути, восстановление изображения высокого разрешения в окне проводится по формуле

x = B(y-д) + д. (5)

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

6. Повышение разрешения в частотной области

Такое же самое изображение повышенного разрешения (исключая пикселы, близкие к границе изображения) можно получить при помощи преобразования Фурье. Сначала объединяются два исходных изображения в одно, в результате чего формируется передискрети-зированное изображение (рис. 1). Пикселы изображения повышенного разрешения, лежащие в центре “своего” окна, можно получить свёртывая передискретизированное изображение с шаблоном, значениями которого являются элементы средней строки матрицы B , формула (5). Этот шаблон различается в зависимости от свойств чётности координат пиксела, подлежащего оценке. Всего необходимо использовать четыре версии шаблона. Можно вычислить передаточную функцию свёртки.

Покажем, как вычислить значения пикселов повышенного разрешения, имеющие чётные координаты. Другие пикселы вычисляются аналогично, но при их вычислении нужно внимательно совместить исходные изображения. Пусть окно имеет размер (4dl + 3) х (4dl + 3) . Тогда значение пиксела, который в изображении повышенного разрешения имеет чётные координаты, можно вычислить по формуле

ij +d

x

211,2 J1

1= x

J1 +d

x

bt IJHyV

І —ij—dj J —Jj—dj

ij +d J1 +d

+ x x

І —ij—dj J —Jj —dj

b^i.A—jJ

(6)

где коэффициенты bfJ составляют центральную строку матрицы B. Для одновременного

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

(x2I2J) = Y(1) *b(1) + Y(2)*b

(1)

К2)

л(2)

(7)

где (x212J) - изображение, состоящее из чётных пикселов восстановленного изображения повышенного разрешения, b(1) и b(2)

- матрицы, составленные из элементов центральной строки матрицы B. Здесь знак “ * ” обозначает двумерную свёртку с сохранением числа элементов первого множителя. После преобразования Фурье

(х2Ш)Л = Y(1) • b(1) + Y(2) • b(2). (8)

ISSN 1028-9763. Математичні машини і системи, 2015, № 1

33

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

Два исходных изображения низкого Пересемпелированное изображение. Пик-разрешения селы, обозначенные •, взяты из первого

изображения. Заштрихованные пикселы взяты из второго изображения

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

пиксел по каждой из координат

Преобразование Фурье передискретизированного изображения

У;;*ГІ,(гУ д 2) = У(1,(гу д 2)exp

(д1+д 2)ani

• ^

2

+ 'Y(2)(d1, д 2)exp

—(д1+д 2)ani

• ^

2

(9)

а передаточная функция свёртки с шаблоном

( -(д1+д2)ani ^

ЬЄГГ1)(д1, д 2) = Ь(1)(д1, д 2) exp

2

+ Ь(2)(д^ д 2)exp

(д1+д 2)ani

2

(10)

здесь a - размер пиксела изображения низкого разрешения, i - мнимая единица. Знаки перед J и J2 зависят от направления взаимного смещения изображений. Приведенные формулы справедливы, если обе координаты вектора смещения второго изображения относительно первого положительные. Можно показать, что свёртка Yeqnter1) * ЬЄіПєгГ) имеет

такие же пикселы с чётными координатами, как и сумма свёрток в правой части (7), но

также имеет и ненулевые пикселы с нечётными координатами.

Точки изображения с чётными координатами можно получить таким способом.

Шаг 1. Найти матрицу B, из элементов центральной строки составить матрицы Ь(1)

и Ь(2), по формуле (9) найти передаточную функцию.

Шаг 2. Составить передискретизированное изображение Y(inter1); в точках, где изо-

eq

бражение не задано, записать нули.

Шаг 3. Свернуть передискретизированное изображение с шаблоном. Для этого провести преобразование Фурье Yeqnter1), умножить на передаточную функцию b(lqnterl) (д, д2 ),

выполнить обратное преобразование Фурье. Получим

X(res.conv)

Y(inter1) jb(interl) xeq ueq

(12)

Шаг 4. У полученного изображения использовать только пикселы с чётными координатами:

34

ISSN 1028-9763. Математичні машини і системи, 2015, № 1

x _ Y(res. conv)

A2I,2J Л21,2J

(13)

В статье [17] описан похожий способ повышения разрешения. Обоснование способа проводится в частотной области. В случае, если окно достаточно большое, взаимное смещение изображений известно точно и составляет ровно половину диагонали пиксела, ФРТ сенсора постоянная в “квадратной” чувствительной зоне, то изображения, восстановленные нашим методом и методом small-kernel Wiener filter, описанным в статье [18], должны быть похожими.

7. Эксперимент

Экспериментальная проверка изложенного алгоритма проводилась на тестовых субпиксельно смещённых инфракрасных изображениях мир Фуко, полученных на прецизионном измерительном стенде, разработанном специалистами Казённого предприятия (КП) специального приборостроения “Арсенал” [19] (рис. 2).

а б

. - ь

* J :

■ I \ V

Рис. 2. Фрагменты тестового цифрового изображения инфракрасной миры, полученного на специализированном измерительном стенде КП “Арсенал”: а, б - исходные субпиксельно смещённые изображения низкого разрешения, в - изображение повышенного разрешения, восстановленное

при помощи описанного алгоритма

в

Статистическая оценка пространственного разрешения, выполненная непосредственно по цифровому изображению рис. 2, демонстрирует повышение в 1,87 раза, что достаточно близко к формальному двукратному повышению в результате передискретизации.

8. Возможные улучшения

Следует разработать способ оценивания матрицы S x без использованного предположения об изображенном объекте.

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

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

9. Выводы

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

ISSN 1028-9763. Математичні машини і системи, 2015, № 1

35

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

СПИСОК ЛИТЕРАТУРЫ

1. Торшина И.П. Оценка пространственного разрешения оптико-электронной системы с помощью компьютерного моделирования / И.П. Торшина, Ю.Г. Якушенков // Ползуновский альманах. -2011. - № 1. - С. 5 - 7.

2. Титов В.С. Направления развития методов, алгоритмов и аппаратных средств повышения качества изображений оптико-электронных систем / В.С. Титов, М.И. Труфанов // Известия ВУЗов. Приборостроение. - 2013. - Т. 56, № 6. - С. 7 - 10.

3. Subpixel image acquisition for detailed aerospace imaging / V.Y. Belenok, V.G. Burachek, V.I. Zatser-kovny [et al.] // Proc. of the Eighth International Conference on Digital Technologies (DT’2011). - Zilina: University of Zilina, 2011. - P. 190 - 193.

4. Vandewalle P. Superresolution images reconstructed from aliased images / P. Vandewalle,

5. Susstrunk, M. Vetterli // Proc. of the SPIE. - 2003. - Vol. 5150. - P. 1398 - 1405.

5. Станкевич С.А. Підвищення просторової розрізненості аерознімання на основі субпіксельної реєстрації зображень / С.А. Станкевич, С.В. Шкляр, М.С. Лубський // Зб. наук. праць Державного науково-дослідного інституту авіації. - Київ: ДНДІА, 2013. - Вип. 9 (16). - С. 110 - 117.

6. Stankevich S.A. Satellite imagery resolution enhancement using subpixel frames acquisition / S.A. Stankevich, S.V. Shklyar, V.M. Tyagur // Journal of Information, Control and Management Systems. - 2013. - Vol. 11, N 2. - P. 135 - 144.

7. Boreman S. Spatial resolution enhancement of low-resolution image sequences: A comprehensive review with directions for future research / S. Boreman, R. Stevenson // Laboratory for Image and Signal Analysis (LISA) Technical Report. - Notre Dame: University of Notre Dame, 1998. - 64 p.

8. Super-Resolution Imaging / Ed. by P. Milanfar. - Boca Raton: CRC Press, 2010. - 490 p.

9. Young S.S. Signal Processing and Performance Analysis for Imaging Systems / Young S.S., Driggers R.G., Jacobs E.L. - Norwood: Artech House, 2008. - 304 p.

10. Lemmerling P. The structured total least-squares approach for non-linearly structured matrices / P. Lemmerling, S. Van Huffel, B. De Moor // Numerical Linear Algebra with Applications. - 2002. -Vol. 9, N 4. - P. 321 - 332.

11. Gerchberg R. Super-resolution through error energy reduction / R. Gerchberg // Journal of Modern Optics. - 1974. - Vol. 21, N 9. - P. 709 - 720.

12. Young S.S. Super-resolution image reconstruction from a sequence of aliased imagery / S.S. Young, R.G. Driggers // Proc. of the SPIE. - 2005. - Vol. 5784. - P. 114 - 114.

13. Stark H. High-resolution image recovery from image-plane arrays using convex projections / H. Stark, P. Oskoui // Journal of the Optical Society of America. - 1989. - Vol. 6, N11. - P. 1715 - 1726.

14. Schultz R.R. A Bayesian approach to image expansion for improved definition / R.R. Schultz,

R. L. Stevenson // IEEE Transactions on Image Processing. - 1994. - Vol. 3, N 3. - P. 233 - 242.

15. Lee S.L. Regularized adaptive high-resolution image reconstruction considering inaccurate subpixel registration / E.S. Lee, M.G. Kang // IEEE Transactions on Image Processing. - 2003. - Vol. 12, N 7. -P. 826 - 837.

16. Fast and robust multi-frame super-resolution / S. Farsiu, D. Robinson, M. Elad [et al.] // IEEE Transactions on Image Processing. - 2004. - Vol. 13, N 10. - P. 1327 - 1344.

17. Shi J. Small-kernel superresolution methods for microscanning imaging systems / J. Shi,

S. E. Reichenbach, J.D. Howe // Applied Optics. - 2006. - Vol. 45, N 6. - P. 1203 - 1214.

18. Reichenbach S.E. Small convolution kernels for high-fidelity restoration / S.E. Reichenbach, S.K. Park // IEEE Transactions on Signal Processing. - 1991. - Vol. 39, N 10. - P. 2263 - 2274.

19. Повышение пространственного разрешения путём субпиксельной обработки изображений / М.А. Попов, С.А. Станкевич, В.М. Тягур [та ін.] // Матеріали Восьмої Міжнар. наук.-техн. конф. “Проблеми телекомунікацій-2014” (ПТ-2014). - Київ: ІТС НТУУ “КПІ”, 2014. - С. 57 - 60.

Стаття надійшла до редакції 12.06.2014

36

ISSN 1028-9763. Математичні машини і системи, 2015, № 1