Документ взят из кэша поисковой машины. Адрес оригинального документа : http://imaging.cs.msu.ru/pub/2011.Book.Krylov_Nasonov.Resampling.ru.pdf Дата изменения: Fri Aug 19 13:57:55 2011 Дата индексирования: Mon Oct 1 19:50:08 2012 Кодировка: Windows-1251 Поисковые слова: принцип гаусса

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ имени М. В. ЛОМОНОСОВА Факультет вычислительной математики и кибернетики

А.С. Крылов, А.В. Насонов

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

Учебное пособие

Москва 2011


УДК 517.6 ББК 22.193 М18

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

Разгулин А.В., профессор, Баяковский Ю.М., доцент,
Крылов А.С., Насонов А.В.

д.ф.-м.н. к.ф.-м.н.

Компьютерное повышение разрешения изображений с использовани-

М18

ем методов математической физики:

Учебное пособие. М.: Издательский отдел факультета ВМК МГУ имени М.В. Ломоносова (лицензия ИД N 05899 от 24.09.2001 г.); МАКС Пресс, 2011. 72 с. ISBN 000-0-00000-000-0
Учебное пособие посвящено использованию регуляризирующих методов решения некорректных задач для повышения разрешения изображений и суперразрешения. Приводится обзор существующих методов повышения разрешения изображений. Строится математическая модель получения изображений для цифровых камер. На основе этой модели ставится задача повышения разрешения изображений, строятся и анализируются регуляризирующие методы для е? решения. Описаны базовые понятия, алгоритмы и численные методы минимизации используемых регуляризирующих функционалов. Пособие рассчитано на студентов старших курсов и аспирантов, специализирующихся в области математической физики и математических методов обработки изображений.

УДК 517.6 ББК 22.193

Рукопись подготовлена при поддержке ФЦП ?Научные и научно-педагогические кадры инновационной России? на 20092013 годы
КРЫЛОВ Андрей Серджевич, НАСОНОВ Андрей Владимирович КОМПЬЮТЕРНОЕ ИЗОБРАЖЕНИЙ С ПОВЫШЕНИЕ РАЗРЕШЕНИЯ МЕТОДОВ

ИСПОЛЬЗОВАНИЕМ ФИЗИКИ

МАТЕМАТИЧЕСКОЙ Учебное пособие

Издательский отдел Факультета вычислительной математики и кибернетики МГУ имени М.В. Ломоносова Лицензия ИД B 05899 от 24.09.01 г. 119992, ГСП-2, Москва, Ленинские горы, МГУ имени М.В. Ломоносова, 2-й учебный корпус Напечатано с готового оригинал-макета в издательстве ООО ?МАКС Пресс? Лицензия ИД N 00510 от 01.12.99 г. Подписано к печати 00.00.0000 г. Формат 60х90 1/16 Усл.печ.л. 4,5 Тираж 100 экз. Заказ 000. 119992, ГСП-2, Москва, Ленинские горы, МГУ имени М.В. Ломоносова, 2-й учебный корпус, 627 к. Тел. (495) 939-3890, 939-3891, Тел./Факс 939-3891.

ISBN 000-0-00000-000-0

c Крылов А.С., Насонов А.В., 2011 c Факультет вычислительной математики и кибернетики МГУ имени М.В. Ломоносова, 2011


3

Оглавление

Введение

5

1

Обзор методов повышения разрешения изображений

11

1.1 1.2
2

Линейные методы повышения разрешения изображений . . Нелинейные методы повышения разрешения изображений

11 15
20

Постановки задач повышения разрешения изображений

2.1 2.2 2.3
3

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

20 22 25

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

3.1 3.2 3.3

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

4

Численные методы, используемые для решения поставленных задач 56

4.1

Субградиентный метод минимизации регуляризирующих функционалов ......................... 56


4

4.2

Метод минимизации квадратичных функционалов на множестве функций с ограниченной полной вариацией 4.2.1 4.2.2 .... 62 62 63 66 Метод условного градиента . . . . . . . . . . . . . . Применение метода условного градиента для задачи подавления эффекта Гиббса . . . . . . . . . . . . . .

4.3

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


5

Введение

Быстрый прогресс в области компьютерной техники позволяет применять для решения задач обработки и анализа изображений вс? более современные математические методы. Более того, многие достижения в области обработки и анализа изображений во многом связаны именно с применением математических методов. Математические методы, ранее возникшие и развившиеся при решении задач в различных прикладных областях, в первую очередь в физике и химии, являются ценным источником создания новых мультимедийных алгоритмов. Во многом, это относится и к методам математической физики, в том числе и к методам решения некорректных задач математической физики, которые в данном пособии составляют базу построения эффективных алгоритмов компьютерного повышения разрешения изображений. Задача повышения разрешения изображений является важной для широкого класса практических приложений, таких как обработка и анализ медицинских изображений, обработка аэрокосмических снимков, обработка данных видеонаблюдения, трансляция видеопотока низкого разрешения на современных широкоформатных дисплеях и ряда других задач. Рост производительности компьютеров в настоящее время позволяет в реальном времени использовать сложные итерационные методы повышения разрешения изображений, а также решать задачу
ния суперразреше-

[1, 2]. Суперразрешение позволяет сразу по нескольким различным

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


6

го качества по сравнению с применением повышения разрешения каждого из изображений низкого разрешения по отдельности. Задача повышения разрешения изображений обладает особенностями, не позволяющими эффективно применять общие методы интерполяции функций для е? решения, поэтому для повышения разрешения изображений требуется разработка специальных методов. Среди таких особенностей можно выделить следующие: 1. Наличие специфичной априорной информации о содержании и структуре изображения. Примером такой информации является информация о спектре непрерывного изображения: если изображение удовлетворяет условиям теоремы Котельникова [3], то возможно его однозначное восстановление по дискретному изображению низкого разрешения. Другими примерами являются предположения о сохранении полной вариации изображения при повышении разрешения, предположения о структуре контуров объектов. Для различных классов изображений эта априорная информация различается, поэтому невозможно разработать универсальный метод повышения разрешения изображений для произвольных изображений. 2. Значимость субъективной оценки результата повышения разрешения. Например, в случае отображения видеопотока низкого разрешения на экранах высокой ч?ткости, визуальное качество изображений имеет решающее значение. При этом имеет значение как отсутствие артефактов, так и правдоподобность получаемого результата (отсутствие искажений, таких как исчезновение мелких деталей или появление новых деталей). Наиболее распростран?нными артефактами, возникающими при повышении разрешения изображений, являются артефакты, связанные с искажением высокочастотной информации: (ступенчатость контуров) и
эффект размытия, алиасинг эффект Гиббса.

В задачах обработки изобложного оконтури-

ражений эффект Гиббса проявляется как эффект


7
вания,

возникающий при недостатке информации о высоких частотах

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

Размытие

Алиасинг (ступенчатость контуров)

Ложное оконтуривание (эффект Гиббса)

Рис. 1: Примеры артефактов, возникающих при повышении разрешения изображений.

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

, а коэффициент отношения шага
коэффициентом

крупной сетки к коэффициенту шага мелкой сетки
увеличения изображений

.

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


8

ражений и повышении разрешения каждого из этих изображений по отдельности. При этом часто производится переход в другое цветовое пространство, например, YUV, где первая компонента представляет из себя значение интенсивности, а оставшиеся две компоненты определяют цвет. В связи с тем, что чувствительность человеческого восприятия к яркости выше, чем к цвету, такой подход позволяет понизить вычислительную сложность при повышении разрешения цветных изображений, применяя для повышения разрешения цветовых компонент быстрые алгоритмы. Существуют алгоритмы, которые работают с компонентами пикселя цветного изображения как с единым целым. Например, в задачах выделения контуров используется цветовой градиент [4], позволяющий находить контуры, которые нельзя найти, используя только яркостную компоненту изображения. Однако в задачах повышения разрешения изображения такие алгоритмы не получили популярности из-за высокой сложности и отсутствия существенного повышения субъективного качества изображения. В настоящее время алгоритмы интерполяции цветных изображений применяются в основном для демозаикинга интерполяции байеровских шаблонов [5]. В большинстве устройств получения цветных цифровых изображений используются матрицы, состоящие из различных фотоэлементов, чувствительных к свету с определенной длиной волны. Используется три типа элементов: чувствительных к красному, зеленому и синему цветам. Эти три типа элементов расположены в виде мозаики, называемой обычно байеровским шаблоном. Задача демозаикинга состоит в получении полноцветного изображения по его байеровскому шаблону. Цифровое изображение в градациях серого

v = {vi,j },

i = 0, 1, . . . , Nx , j = 0, 1, . . . , Ny ,

представляет собой конечномерную двумерную таблицу размера (Nx +


9

1) Ч (Ny + 1), элементам которой, называемым

пикселями,

присвоено зна-

чение из некоторого конечного множества значений, в качестве которого обычно используется множество целых чисел на отрезке от 0 до 255. Данное значение называется интенсивностью. В случае цветного изображения значением пикселя является вектор из тр?х значений, соответствующих цветовым компонентам. Удобно использовать представление изображения в виде вещественной сеточной функции u(ihx , j hy ) = vi,j , заданной на двухмерной равномерной сетке

hx

,hy

= {(ihx , j hy ) : 0 i Nx , 0 j Ny }, hx > 0, hy > 0.

Для того, чтобы избавиться от задания граничных условий в явном виде для каждого конкретного метода повышения разрешения изображений, удобно произвести продолжение сеточной функции u с конечномерной сетки hx
,hy

на бесконечномерную



hx ,hy

= {(ihx , j hy ) : i, j Z}, hx > 0, hy > 0.

(1)

При этом значения функции u совпадают со значениями пикселей изображения u для 0 i Nx , 0 j Ny , а значения за пределами этого прямоугольника вычисляются пут?м ч?тного (можно также использовать и другие способы, например неч?тное продолжение) продолжения относительно его сторон. Ч?тное продолжение для изображения размера


10

(Nx + 1) Ч (Ny + 1) выглядит следующим образом: u i,j при 0 i < Nx , 0 j < Ny , u2N -i,j при Nx i < 2Nx , 0 j < Ny , x ui+2Nx p,j +2Ny q = u i,2Ny -j при 0 i < Nx , Ny j < 2Ny , u2N -i,2N -j при Nx i < 2Nx , Ny j < 2Ny , x y
где p, q Z.

(2)

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

u(x, y ) =
i,j =-

ui,j (x - ihx ) (y - ihy ),

(3)

где (t) дельта-функция:
+

f (t) (t - t0 )dt = f (t0 ).
-

В данном учебном пособии рассматривается задача повышения разрешения и суперразрешения однокомпонентных изображений. Эта задача ставится как обратная задача для задачи понижения разрешения изображений. Она является некорректной, для е? решения используются методы, основанные на регуляризации Тихонова [6]. Также рассматривается задача подавления эффекта Гиббса после повышения разрешения. Описаны численные методы для решения поставленных задач.


11

Глава 1
Обзор методов повышения разрешения изображений

1.1

Линейные методы повышения разрешения изображений

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


f (x0 , y0 ) =

u(x0 - x, y0 - y )K
- - + +

xy , hx hy

dxdy =
(1.1)

=
i=- j =-

ui,j K

x0 y0 - i, -j , hx hy

где f (x, y ) интерполированное изображение. На данном этапе не накладывается никаких ограничений на значения пикселей и размер исходного изображения. Как правило, используется ядро вида K (x, y ) =

K (x)K (y ). Это позволяет существенно повысить скорость работы алгоритма, так как в этом случае интерполяция разбивается на последовательную одномерную интерполяцию сначала по одной оси, затем по


12

другой:

+

fi (y ) =
j =- +

ui,j K fi (y )K
i=-

y -j , hy x -i . hx

(1.2)

f (x, y ) =

При этом на ядро K (t) накладываются следующие условия: 1. Условие интерполяции

K (0) = 1, K (n) = 0, n = +1, +2, . . . .

Это условие нужно для выполнения равенства значений исходного и интерполированного изображений

f (ihx , j hy ) = ui,j ,
в узлах сетки hx
,hy

i, j Z

.

2. Финитность ядра

K (t) = 0,

|t| > p > 0.

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

p.
3. Условие нормировки коэффициентов
+

K (t + i) = 1,
i=-

0 t < 1.

Данное условие следует из того, что интерполяция изображения, все пиксели которого имеют одинаковое значение ui,j = C , должна давать константную функцию f (x, y ) = C . В работах [810] собраны воедино, обобщены и подробно рассмотре-


13

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

1 для - 1/2 t < 1/2, K (t) = 0 иначе.
Это самый быстрый метод интерполяции. В случае ресамплинга с целым коэффициентом увеличения он представляет собой простое повторение пикселей изображения. Данный метод обладает серь?зным недостатком: получаемая функция f (x, y ) является разрывной, а на изображении ярко выражен эффект ступенчатости. 2. Метод интерполяции первого порядка:

1 - |t| для - 1 < t < 1, K (t) = 0 иначе.
В двумерном случае этот метод называют
ей билинейной интерполяци-

. Получаемая с помощью этого метода функция f (x, y ) является непре-

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

):

(a + 2)|t|3 - (a + 3)|t|2 + 1 для |t| 1, K (t) = a|t|3 - 5a|t|2 + 8a|t| - 4t для 1 < |t| < 2, 0 иначе.

(1.3)

Этот метод наиболее часто применяется для повышения разрешения


14

изображений из-за малого значения p (p = 2) и хорошего баланса между артефактами размытия, алиасинга и оконтуривания кра?в. 4. ?Идеальная? интерполяция, основанная на использовании теоремы Котельникова [3]:

K (t) = sinc t =

sin t . t

(1.4)

Теорема Котельникова гарантирует восстановление непрерывного сигнала по дискретному при условии, что в спектре непрерывного сигнале не было частот выше
1 2 max(hx ,hy )

. На практике применение этого метода за-

труднительно из-за того, что ядро (1.4) не является финитным (p = ), а спектр естественных изображений не является ограниченным. При интерполяции изображения с помощью данного метода, на изображении высокого разрешения наблюдается ярко выраженный эффект ложного оконтуривания, причиной которого является эффект Гиббса. 5. Метод Ланцоша, представляющий собой приближение идеальной интерполяции финитным ядром:

Kp (t) =

sinc(t) ћ sinc 0

t p

для |t| < p, иначе.

(1.5)

Параметр p N зада?т размер ядра. На практике обычно используют

p = 2 и p = 3.
6. Метод интерполяции с использованием фильтра Гаусса:

1 K (t) = e 2

-

t2 2 2

.

Данный метод интерполяции не удовлетворяет первому условию, тем не менее, он обладает рядом полезных свойств и широко используется в обработке изображений. Например, он используется для вычисления производных любых порядков в произвольных точках на изображении. Ядро фильтра не является финитным, тем не менее, оно быстро убывает при росте t, поэтому на практике обычно принимают K (t) = 0 при


15

|t| > t , где значение t выбирается таким образом, чтобы ошибка не
превышала заданного значения. При использовании t = 3 значение ошибки существенно ниже разницы между отображаемыми градациями интенсивности современных мониторов, поэтому на практике чаще всего используют именно это значение. Для эффективного вычисления св?ртки с фильтром Гаусса используется его приближение с помощью рекуррентных фильтров [11]. Среди класса линейных методов невозможно выбрать наилучший метод. Любой линейный метод представляет собой баланс между тремя типами артефактов: размытия, алиасинга и эффекта Гиббса (см. во введении, рис. 1). Подавление одного из артефактов приводит к усилению других артефактов.
1.2 Нелинейные методы повышения разрешения изображений

Добиться более качественных результатов можно с помощью нелинейных методов, в которых функция усреднения K (x, y ) зада?тся отдельно для каждого интерполируемого пикселя и, вообще говоря, зависит от значений пикселей интерполируемого изображения. Примером нелинейных методов является класс градиентных методов [12, 13]. В основе градиентных алгоритмов лежит тот факт, что интерполяция вдоль кра?в (контуров) деталей изображения да?т лучшие результаты, чем обычная линейная интерполяция. В целом результат, получаемый градиентными методами, получается близок к результату бикубической интерполяции, но эффект алиасинга оказывается практически полностью подавлен. Один из способов реализации такой интерполяции основан на использовании функции Гаусса с переменными радиусами по разным направлениям:
- 1 K (x , y , x0 , y0 ) = e 2 x y
x2 2 2x

-

y2 2 2y

.


16

y
1.2

1

0.8

0.6

x


K(x, y)

0.4 3 2 1 0 0 -1 -0.2 -3 -2 -1 0 1
x

0.2

-2 2 -3 3
y

Рис. 1.1: Пример ядра K (x, y ) в точке контура при интерполяции с помощью градиентных методов.

Здесь ядро K (x , y , x0 , y0 ) задано в декартовой системе координат

Ox y c центром в интерполируемом пикселе O = (x0 , y0 ) и осью Ox , совпадающую с направлением градиента исходного изображения в (x0 , y0 ). Радиус y фиксирован, а радиус x выбирается на основе анализа модуля градиента изображения: чем больше модуль градиента, тем меньше значение x . Пример ядра при интерполяции в точке границы привед?н на рис. 1.1. Ещ? один алгоритм, основанный на использовании градиента это метод WADI (Warped Distance) [14]. В н?м значение интерполируемого пикселя представляет собой взвешенную сумму значений четыр?х ближайших пикселей, прич?м веса выбираются в зависимости от расстояния до этих пикселей и модуля производной в этих пикселях. Чем больше производная, тем меньше весовые коэффициенты. Алгоритм NEDI (New Edge-Directed Interpolation) увеличивает изображение в 2 раза, используя предположение о самоподобии фрагментов изображения высокого разрешения и соответствующих им фрагментов изображения низкого разрешения [15, 16]. Значения интерполируемых пикселей представляют собой взвешенную сумму четыр?х соседних пикселей, при этом веса вычисляются в предположении, что исходное изображение было получено с теми же весами пут?м увеличения уменьшен-


17

ного в 2 раза исходного изображения. NEDI хорошо справляется с контурами объектов (визуально лучше, чем градиентные методы), но плохо обрабатывает участки изображения, для которых предположение о самоподобии неверно. К методу NEDI близок класс фрактальных алгоритмов, в основе которых лежит свойство самоподобия целых блоков изображения. Принцип работы алгоритмов фрактального ресамплинга аналогичен алгоритмам фрактального сжатия и основан на методах фрактального кодирования. Различные методы фрактального кодирования, в том числе использующие итерации сжимающих отображений (Iterated Function System, IFS), рассматриваются в статьях [1720]. Широкий класс методов представляет задачу повышения разрешения изображений как обратную задачу для уравнения

Az = u,

(1.6)

где u исходное изображение низкого разрешения, z искомое изображения высокого разрешения, A оператор понижения разрешения. Таким образом, задача формулируется в виде: построить такое изображение высокого разрешения, которое после уменьшения даст известное изображение низкого разрешения. Эта постановка задачи повышения разрешения как обратной задачи позволяет использовать для ее решения методы теории решения обратных задач. Эта теория была создана А.Н.Тихоновым [6] и ориентирована, в первую очередь, на решение некорректных задач математической физики. Однако, за последнее время, регуляризирующие методы решения обратных задач математической физики нашли широкое применение в различных областях прикладной математики. Созданный математический аппарат для использования дополнительной информации о решении оказался эффективным и в области обработки изображений. Обратная задача, определяемая уравнением (1.6), является некоррект-


18

но поставленной. Для е? решения применяются следующие методы, основанные на использовании априорной информации об изображении: 1. Регуляризиризующие методы [21]. В общем случае задаются два функционала: функционал соответствия изображения высокого разрешения изображению низкого разрешения Id (z , u), например, невязка

Id (z , u) = Az - u

2 2

и функционал соответствия изображения высокого разрешения априорной информации Is (z ) стабилизатор [22]. Для нахождения изображения высокого разрешения производится минимизация регуляризирующего функционала

z = arg min (Id (z , u) + Is (z )) ,
z

где параметр регуляризации контролирует баланс между соответствием изображению низкого разрешения и априорной информации. Также используются следующие регуляризирующие методы [2325]

zC = arg min Id (z , u),
z M
C

MC = {z |Is (z ) C }

и

z = arg min Is (z ),
z M


M = {z |Id (z , u) }.

2. Метод обратной проекции ошибки [26, 27], заключающийся в итерационной минимизации невязки Az - u :

zk

+1

= zk + U (Azk - u),

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

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


19

результат с хорошим визуальным качеством. 3. Метод проекции на выпуклые множества, заключающийся в проектировании приближения z0 на множество изображений M , где M выпуклый компакт [28, 29]:

zR = prM z0 = arg min z - z0 2 . 2
z M


20

Глава 2
Постановки задач повышения разрешения изображений

2.1

Модель получения цифровых изображений

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







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

В цифровой камере изображение проецируется на матрицу, состоящую из светочувствительных элементов (сенсоров), каждый из которых соответствует пикселю на изображении. Значение пикселя ui,j равно сум-


21

марной интенсивности непрерывного светового потока

f (x, y ) F = L(R2 ),
попавшего на сенсор:

ui,j =
x,y

f (x, y )Ki,j (x - xi , y - yi )dxdy ,

(2.1)

где (xi , yi ) координаты центра пикселя, Ki,j (x, y ) функция усреднения в координатах относительно центра пикселя, называемая также
point spread function

(PSF), или

функцией распределения точки

, L(R2 )

банахово пространство вещественных функций, определ?нных на R2 с некоторой нормой

ћ

F

. (1).

Предполагается, что все сенсоры одинаково усредняют световой поток

Ki,j (x, y ) = K (x, y ) и расположены в узлах равномерной сетки hx

,hy

Дополнительно накладывается условие центральной симметричности на

K (x, y ). Этим условиям удовлетворяет большинство современных камер.
При данных условиях модель (2.1) принимает вид:

ui,j = [f K ](xi , yj ),
где f K операция св?ртки
+ +

(2.2)

[f (x, y ) K (x, y )](x0 , y0 ) =
- -

f (x, y )K (x0 - x, y0 - y )dxdy .

Для прямоугольных сенсоров функция распределения точки принимает вид

1 при - K (x, y ) = 0 иначе.

hx 2

x<

hx 2

,-

hy 2

y<

hy 2

,

(2.3)

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


22

мерная функция Гаусса

K (x, y ) = G
где
радиусы

x

,y

x2 y2 1 exp - 2 - 2 , (x, y ) = 2 x y 2 x 2 y

2 2 x и y (квадратные корни дисперсий x , y ) выбираются

исходя из размеров сенсоров hx и hy . Для сенсоров с одинаковыми размерами по вертикали и горизонтали

h = hx = hy используются одинаковые значения радиусов = x = y ,
и функция распределения точки записывается в виде

G (x, y ) =

1 x2 + y exp - 2 2 2 2

2

.

(2.4)

Значение радиуса бер?тся пропорционально размеру сенсоров =

0 h. Параметр 0 определяется конструкцией камеры, обычно используются значения 0 из интервала [0.4, 0.5]. При обработке реальных изображений в подавляющем большинстве случаев нет информации о камере, поэтому используются универсальные параметры, хорошо приближающие параметры камеры. При этом предполагается, что камеры имеют квадратные сенсоры hx = hy = h, а функция распределения точки описывается в виде (2.4). Построение методов повышения разрешения изображений производится исходя из этих параметров.
2.2 Постановка задачи ресамплинга изображений

Обозначим пространство (h ) банахово пространство вещественных сеточных функций, заданных на бесконечномерной равномерной сетке с шагом h

h = {(xi , yj ) : xi = ih, yj = j h,
с нормой

i, j Z}

ћ

(h )

.


23
Построение непрерывного изображения

Согласно введ?нной в предыдущем разделе математической модели (2.2), задача построения непрерывного изображения f L(R2 ) по изображению u (h ) формулируется в виде

ui,j = [f G0 h ](xi , yj ) = [Dh H0 h f ]i,j ,
где H
0

h

оператор св?ртки с функцией Гаусса G0 h ,

Dh : L(R2 ) (h )
оператор дискретизации:

[Dh f ]i,j = f (ih, j h).
Таким образом, в непрерывном случае задача повышения разрешения изображений ставится в виде обратной задачи для задачи получения дискретного изображения с камеры:

Ah f = Dh H0 h f = u,
Ресамплинг изображений

f L(R2 ), u (h ).

(2.5)

В реальных задачах требуется не построение непрерывного изображения

f , а переход на сетку с меньшим шагом ресамплинг изображений.
Задача понижения разрешения изображений ставится в следующем виде: по имеющемуся изображению z (h/s ), полученного с помощью камеры с размером сенсоров

h/s,

построить

изображение

u (h ) для камеры с размером сенсоров h, где s коэффициент
масштабирования, s > 1.


24

Согласно математической модели (2.2):

zi,j = [f G

0

h/s

]

ih j h s, s

,

ui,j = [f G0 h ](ih, j h).
Используя свойство Ga Gb = G записана в виде
a2 +b
2

, вторая формула может быть

ui,j = [f G

0

h/s

G h
0

1-1/s2

]

ih j h s, s

.

(2.6)

Для вычисления ui,j необходимо иметь непрерывное изображение

f G

0

h/s

, но известным является только дискретное изображение zi,j . В
0

этом случае используется -представление z (3) вместо f G ^ в качестве аппроксимации:

h/s

в (2.6)

ui,j = [z G h ^
0

1-1/s2

]

ih j h s, s

.

При этом оста?тся лишь зависимость от s:
i,j Z

zi,j G h
0

1-1/s
2

2

(i0 h -

ih s

, j0 h -

jh s

) =

ui

0

,j0

=

i,j Z i,j Z

G h
0 0

1-1/s

(i0 h -

ih s

, j0 h -

jh s

)

zi,j G G

0

1-1/s2

i (i0 - s , j0 - j ) s

=

i,j Z i,j Z

s2

-1

(i0 s - i, j0 s - j ) (si0 - i, sj0 - j )

=

zi,j G G
0


0

s2 -1

=



i,j Z

s2 -1

(si0 - i, sj0 - j )

.

Если коэффициент масштабирования s целый, то все узлы грубой сетки h являются узлами исходной сетки
h/s

. В этом случае для нахож-


25

дения u достаточно вычислить св?ртку:

ui0

,j0

=
i,j Z

zi,j G

0


0

s2 -1

(i0 s - i, j0 s - j ) =
s2 -1

= [z G

s2 -1 ]i0 s,j0 s

= [H0

z ]i

0

s,j0 s

.

Если коэффициент s не является целым, то для вычисления значений пикселей в узлах сетки h используются методы интерполяции, например, билинейная интерполяция. Таким образом, задача понижения разрешения изображений в s раз формулируется в виде:

u = As z = Ds H


0

s2 -1

z,

(2.7)

где Ds это оператор перехода на сетку с s раз большим шагом, действующий из пространства (h/s ) в пространство (h ):

[Ds z ]i,j = zsi,sj .
Оператор As является линейным оператором. Задача повышения разрешения изображений ставится в виде обратной задачи для задачи понижения разрешения:

As z = u,
2.3

u (h ), z (h/s ).

(2.8)

Постановка задачи суперразрешения

Суперразрешение

это построение изображения высокого разрешения

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


26

Рис. 2.2: Иллюстрация для задачи суперразрешения: восстановление изображения с вдвое большим разрешениям по четыр?м изображениям низкого разрешения.

2h имеют прямоугольную функцию распределения точки (2.3), а изображение было снято 4 раза: первый раз (u1 ) без сдвига, второй (u2 ) со сдвигом на h по горизонтали, третий (u3 ) со сдвигом на h по вертикали, четв?ртый (u4 ) со сдвигом на h по горизонтали и по вертикали. Требуется построить изображение z без сдвига с вдвое большим разрешением. Значения пикселей первого изображения низкого разрешения можно выразить через пиксели z :

1 u1 = (z2i,2j + z2i+1,2j + z2 i,j 4

i,2j +1

+ z2i+1,

2j +1

),

Для второго изображения справедлива формула

1 u2 = (z2i i,j 4

+1,2j

+ z2i+2,2j + z2

i+1,2j +1

+ z2i+2

,2j +1

).

Аналогично выписываются соотношения для пикселей третьего и чет-


27

в?ртого изображений. Объединяя эти соотношения, получим:

1 (zi,j + zi+1,j + zi,j 4
где

+1

+ zi+1,j +1 ) = vi,j ,

v2i,2j = u1 , i,j v2i+1 v2i, v2i+1
,2j

= u2 , i,j = u3 , i,j = u4 . i,j

2j +1

,2j +1

Таким образом, рассмотренная задача свелась к задаче обращения св?ртки изображения z с прямоугольным фильтром 2 Ч 2. Привед?м общую постановку задачи суперразрешения, использующую описанную выше модель понижения разрешения изображений. Пусть f L(R2 ) исходное непрерывное изображение. С течением времени объекты на изображении смещаются, поэтому камера получает на вход изображения

f k L(R2 ), k = 1, 2, . . . , K,
полученные из f в результате применения операторов движения Tk . Оператор движения
k k Tk = (Tx , Ty ) k k определяется двумя функциями Tx (x, y ) и Ty (x, y ), задающими соответ-

ствие между координатами (x, y ) на k -м изображении и координатами
k k (x + Tx (x, y ), y + Ty (x, y )) на исходном изображении: k k [Tk f ](x, y ) = f (x + Tx (x, y ), y + Ty (x, y )).

(2.9)

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

uk (h ), k = 1, 2, . . . , K,


28

прич?м

uk = Dh H0 h Tk f .

(2.10)

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

Ak f = Dh H0 h Tk f = uk , h

k = 1, 2, . . . , K,

f L(R2 ), uk (h ),
где операторы Dh и H фильтром Гаусса.
0

(2.11)

h

это введ?нные ранее для задачи повышения

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

uk (h ),

k = 1, 2, . . . , K

изображения низкого разрешения, полученные с помощью камеры с размером сенсоров h, z (h/s ) реконструируемое изображение высокого разрешения, полученное с помощью камеры с размером сенсоров

h/s, s коэффициент масштабирования, s > 1.
Согласно математической модели (2.2):

zi,j = [f G

0

h/s

](ih/s, j h/s) k = 1, 2, . . . , K.

uk = [(Tk f ) G0 h ](ih, j h), i,j

Пут?м замены второй формулы на следующее приближение:

uk = [(Tk (f G i,j

0

h/s

)) G h
0

1-1/s

2

](ih, j h),

k = 1, 2, . . . , K.

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


29

носительно z :

Ak z = Ds H s
где операторы Ds и H

0


0

s2 -1 Tk

z = uk ,

k = 1, 2, . . . , K,

(2.12)

s2 -1

это введ?нные ранее для задачи ресам-

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


30

Глава 3
Построение и анализ регуляризирующих методов решения задачи повышения разрешения изображений и суперразрешения

Для нахождения решения задачи повышения разрешения изображений в непрерывной постановке (2.5) требуется построение обратного оператора для оператора Ah = Dh H0 h . Обратный оператор для оператора св?ртки с фильтром Гаусса H0 h является неограниченным, поэтому решение обратной задачи (2.5) является неустойчивым. В дискретной постановке (2.8) задача обращения св?ртки является плохо обусловленной. Помимо неустойчивости, отсутствует единственность решения, так как операторы Dh , Ds являются линейными операторами с

det Dh = 0,

det Ds = 0.

По изображению низкого разрешения можно построить несколько изображений высокого разрешения, прич?м все эти варианты будут допустимыми (см. рис.3.1). В связи с этим не требуется единственность решения задачи повышения разрешения изображений, однако требуется компактность множества решений. Вышеописанное верно и для задачи суперразрешения (2.11, 2.12). Но


31

Рис. 3.1: Пример неоднозначной интерполяции ступенчатого контура на изображении.

в отличие от задачи повышения разрешения изображений, всегда имеющей решение, решение задачи суперразрешения может отсутствовать в случае, если системы уравнений (2.11), (2.12) являются несовместными. Таким образом, задачи повышения разрешения изображений и суперразрешения являются некорректно поставленными. Для нахождения решения задач повышения разрешения изображений и суперразрешения используются регуляризирующие алгоритмы, базирующиеся на методе регуляризации А.Н. Тихонова [6]. Этот метод был создан в 60-х годах ХХ столетия как метод решения обратных задач математической физики, но в настоящее время является крайне актуальным и востребованным и при решении задач обработки изображений.
3.1 Построение регуляризирующих методов решения задачи повышения разрешения изображений и суперразрешения

В данном разделе производится построение регуляризирующих методов для задачи нахождения решения операторного уравнения

Af = u

(3.1)


32

относительно f , где f F, u U . F и U банаховы пространства с нормами

ћ

F

и

ћ

U

соответственно и метриками

F (f1 , f2 ) = f1 - f

2F

, ,

U (u1 , u2 ) = u1 - u2

U

A линейный непрерывный оператор, отображающий F на U .
Введ?м обозначения:

u0 точное значение правой части уравнения (3.1); u приближ?нное значение правой части уравнения (3.1), прич?м U (u0 , u ) ; A-1 u = {f F : Af = u} множество решений уравнения для
правой части u;

A-1 u0 множество точных решений уравнения при u = u0 ; A-1 u множество приближ?нных решений уравнения при u = u .
В случае единственности решения задача (3.1) называется устойчивой [6], если для любого > 0 существует > 0 такое, что для любых допустимых значений правой части u0 и u таких, что U (u , u0 ) выполняется

F (A-1 u , A-1 u0 ) < .
Для случая неединственности решения используются понятия
расстояния бета-

и

бета-сходимости

для оценки отклонения множества ре-

шений A-1 u для приближ?нной правой части от множества точных решений A-1 u0 [30].
Определение 1. Бета-расстоянием

между множествами M1 и M2 в

метрическом пространстве называется величина

(M1 , M2 ) = sup (x, M2 ),
xM
1

где

(x, M ) = inf (x, y ).
y M


33
Определение 2.

Последовательность множеств {Mn } в метрическом
бета-сходящейся

пространстве называется

к M при n , если

n

lim (Mn , M ) = 0,

где

(Mn , M ) = sup (xn , M ) = sup
xn M
n

xn M

n

xM

inf (xn , x) .

Отметим, что если последовательность {Mn } бета-сходится к некоторому множеству M , то она является бета-сходящейся к любому множеству M , содержащему M .
Утверждение 1.

Пусть множества Mn являются подмножествами

компакта M0 в банаховом пространстве, а множество MP множество всех предельных точек всевозможных последовательностей {xn }, составленных из элементов множеств Mn : xn Mn . Тогда последовательность

{Mn } является бета-сходящейся к MP .
Доказательство.

Докажем от противного. Пусть последовательность

{Mn } не является бета-сходящейся к MP . Тогда последовательность {yn },
где

yn = sup (xn , MP )
xn M
n

не сходится к нулю при n , а значит, существует > 0, для которого можно выделить подпоследовательность {ynk } такую, что

y

nk

=
x

sup (xnk , MP ) > .
nk

M

nk

В каждом из множеств Mnk возьм?м элемент xnk Mnk такой, что

(xnk , MP ) .
Так как M0 компакт, то из последовательности {xnk } можно выде-


34

лить сходящуюся подпоследовательность {x
m

nk

m

} x, прич?м

lim (xnkm , MP ) = (x, MP ) ,

а значит, предельная точка x не принадлежит множеству MP . Но по условию последовательность {xnkm } является подпоследовательностью некоторой последовательности, составленной из элементов множеств Mn , и все е? предельные точки, включая x, принадлежат MP . Полученное противоречие показывает, что последовательность {Mn } является бетасходящейся к MP . Утверждение доказано. В терминах бета-сходимости задача (3.1) является корректно поставленной, если для любого u0 из множества допустимых значений правой части множество решений A-1 u0 непусто и устойчиво к малым изменениям правой части, то есть для любого > 0 существует > 0 такое, что

(A-1 u , A-1 u0 ) = sup F (f , A-1 u0 ) <
f A
-1

u



для любых u таких, что U (u0 , u ) . Предположим, что вместо точного значения правой части u = u0 уравнения (3.1) известно его приближ?нное значение u . Множество точных решений уравнения (3.1) не является устойчивым к малым изменениям

u, поэтому нельзя использовать множество A-1 u в качестве приближ?нного решения уравнения (3.1). Альтернативой этому является использование метода регуляризации Тихонова, заключающегося в рассмотрении задачи, близкой к (3.1), решение которой обладает устойчивостью и бетасходится к точному решению уравнения (3.1) при 0. Задачу нахождения приближ?нного решения для (3.1) можно поставить одним из следующих способов [24], где [f ] стабилизатор, описанный ниже: I. Нахождение inf ( Af - u
U

+ [f ]) , > 0.
U

II. Нахождение inf [f ] для Af - u

, > 0.


35

III. Нахождение inf Af - u Iа. Нахождение inf

U

для [f ] C, C > 0.

Также используется следующая постановка:

Af - u

2 U

+ [f ] , > 0.

Для задачи (3.1) рассмотрим построение регуляризирующего метода в постановке II, а затем применим теоремы эквивалентности для перехода к постановкам I, Iа и III. Для уравнения (3.1) с неединственным решением используется понятие
регуляризирующего оператора

f = R(u , ),
который обладает свойствами:

u U, f F,

1) существует такое число 1 > 0, что оператор R(u , ) опред?лен для всякого , 0 1 и любого u U такого, что

U (u , u0 ) ;
2) для всякого > 0 существует 0 = 0 ( , u ) 1 такое, что из неравенства

U (u , u0 ) 0
следует неравенство

(f , A-1 u0 ) = sup F (f , A-1 u0 ) .
f f


Используем понятие
затора

стабилизирующего функционала

, или

стабили-

неотрицательного функционала [f ], определ?нного на всюду

плотном в F подмножестве F1 множества F и удовлетворяющего условиям: 1) элементы f A-1 u0 принадлежат его области определения; 2) для всякого числа d > 0 множество F компактно в F . Определим функционал R(u , ). Пусть для некоторого u существует
1,d

= {f F1 : [f ] d}


36

непустое множество элементов f , которые минимизируют [f ] на множестве

Q = {f : f F1 , U (Af , u ) }.
Множество f можно рассматривать как применение к правой части уравнения (3.1) некоторого оператора R, зависящего от параметра , т.е.

f = R(u , ) = arg min [f ] = {f Q : [f ] = min [f ]}.
f Q f Q

(3.2)

Это множество является подмножеством компактного в F множества

F

1,d

при d = min [f ].
f Q

Предполагается, что такое множество существует для любых 0 и

u U . Тогда справедлива
Теорема 1

. R(u , ) (3.2) является регуляризирующим оператором

для уравнения (3.1). Доказательство теоремы аналогично доказательству соответствующей теоремы из [6] со следующими отличиями: 1) Вместо сходимости R(u , ) к единственному точному решению уравнения Az = u0 при 0 происходит -сходимость к множеству точных решений A-1 u0 . 2) Множества значений R(u , ) являются подмножествами компактного в F множества

F

1,0

= {f F1 : [f ]

f0 A

min [f0 ]}, -1
u0

таким образом, R(u , ) бета-сходится к компактному подмножеству множества точных решений A-1 u0

F1,0 .

Рассмотренный подход, заключающийся в минимизации стабилизатора [f ] на множестве приближ?нных решений U (Af , u ) < (постановка II), неудобен с вычислительной точки зрения. Для задачи повышения разрешения изображений эффективнее ставить задачу регуляризации в


37

виде минимизации регуляризирующего функционала (постановка I):

R(u , ) = arg min ( Af - u
f F
1

U

+ [f ]) ,

(3.3)

где

параметр регуляризации

.

Для случая единственности решения уравнения (3.1) эквивалентность регуляризирующих методов в вариантах постановки I, II и III достигается при выполнении следующих условий [24]: 1) Пространства F и U нормированные, метрика порождена нормой пространства

U (u1 , u2 ) = u1 - u2

U

,

область определения D(A) оператора A линейное многообразие. 2) [f ] выпуклый функционал с областью определения D(A). 3) Для любого f D(A) функция () = [f ] является строго возрастающей функцией параметра 0. Для случая неединственности решения уравнения (3.1) используется видоизмен?нное третье условие и накладывается дополнительное условие на регуляризирующий функционал: 3) Для любого f D(A), f = 0, [f ] = 0 функция () = [f ] является строго возрастающей функцией параметра 0. 4) Если множество f F решение задачи регуляризации в постановке I, то для любых f1 , f2 f выполняется

Af1 - u

U

= Af2 - u

U

.

В этом случае эквивалентность регуляризирующих методов в постановках I и II достигается при выполнении условий теоремы:
Теорема 2

. Пусть A линейный оператор с областью определения

F1 (F1 = F ) и областью значений AF1 U , F и U нормированные


38

пространства, выполнено неравенство

u - u0

U

< ,

0<<

1 u0 2

U

, u0 = 0.

Если f решение II, прич?м [f ] > 0, то существует положительное число > 0 такое, что F минимизирует функционал

M [f ; u , ] = Af - u
прич?м Af - u II.
U

U

+ [f ],

= . И наоборот, если f минимизирует функционал
U

I для некоторого > 0 и = Af - u

, [f ] > 0, то f решение

Доказательство теоремы аналогично доказательству соответствующей теоремы из [24].
Замечание

. От условия 4 можно избавиться, если рассматривать эк-

вивалентность постановок Ia и II. В этом случае на множестве f выполняется

Af - u

2 U

+ [f ] const.

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

Af1 - u

U

= Af2 - u

U

выполняется для любых f1 , f2 f . Таким образом, в случае выполнения условия 4, регуляризирующие методы для задачи повышения разрешения изображений в постановках I и Ia будут эквиваленты. В этом случае используется постановка, наиболее удобная для минимизации конкретного регуляризирующего функционала.


39

Аналогично формулируется и доказывается теорема эквивалентности регуляризирующих методов в постановках I, Ia и III.
Применение метода регуляризации к задаче повышения разрешения изображений

Для задачи повышения разрешения изображений (2.5)

Ah f = Dh H0 h f = u,
определим оператор

f F, u U

f = Rh [u; ] = arg min h2 Ah f - u
f F

U

+ [f ]

(3.4)

Здесь мы добавили нормировочный коэффициент h2 к невязке, чтобы избавиться от зависимости от h. Достаточным условием того, что Rh [u; ] является регуляризирующим оператором для задачи повышения разрешения изображений, является выполнение условий теорем 1 и 2. В частности, теорема 1 требует, чтобы [f ] являлся стабилизатором, для которого выполнено условие компактности множества {f : [f ] d} для всех d > 0. Это условие не выполняется в пространстве F = L(R2 ) для обычно используемых в задачах повышения разрешения изображений стабилизаторов и норм пространств F и U . В связи с этим на пространства F и U накладываются ограничения. Размеры матрицы сенсоров камеры, с помощью которой регистрируются реальные изображения низкого разрешения, конечны. Получаемые с помощью камеры изображения u определены на конечномерной сетке

h = {(ih, j h) : 0 i Nx , 0 j Ny }.
Значения пикселей изображения u (h ) с индексами (i, j ), лежащими за пределами прямоугольника [0, Nx ] Ч [0, Ny ], вычисляются с по-


40

мощью ч?тного продолжения (2). Для пространства изображений u, полученных из u (h ) пут?м продолжения на сетку h используется обозначение (h ). В этом случае мы вводим предположение о том, что изображение высокого разрешения f L(R2 ), соответствующее изображению низкого разрешения из (h ), также является ч?тным продолжением изображения

f



L([0, Rx ] Ч [0, Ry ]),

заданного

на

прямоугольнике

[0, Rx ] Ч [0, Ry ], Rx = hNx , Ry = hNy : f (x + 2Rx p, y + 2Ry q ) = f f = f f (x, y ) (2Rx - x, y ) (x, 2Ry - y )
при 0 x Rx , 0 y Ry , при Rx x < 2Rx , 0 y < Ry , при 0 x < Rx , Ry y < 2Ry ,

(2Rx - x, 2Ry - y ) при Rx x < 2Rx , Ry y < 2Ry , p, q Z.

Пространство

ч?тных

продолжений

изображений

из

L([0, Rx ] Ч [0, Ry ]) на R2 обозначается как L(R2 ).
Значения норм элементов в пространствах (h ) и L(R2 ) определяются через значения норм соответствующих им элементов из пространств

(h ) и L([0, Rx ] Ч [0, Ry ]): u f
(h ) L(R2 )

=u =f

(h )

, .

L([0,Rx ]Ч[0,Ry ])

Использование пространств F = L(R2 ) и U = (h ) вместо L(R2 ) и

(h ) для задачи повышения разрешения изображений (2.5) не накладывает ограничений на е? практическое применение. В данных простран-


41

ствах метод минимизации (3.4) принимает вид

f = arg min h2 Ah f - u
f F 2

U

+ [f ] ,

(3.5)

F = L(R ), U = (h ).
В разделе 3.2 мы опишем выбор стабилизаторов и норм пространств

F , U , при которых метод (3.5) является регуляризирующим для задачи
повышения разрешения изображений. Для нахождения приближ?нного решения задачи повышения разрешения изображений в дискретном случае (2.8) используется следующий метод:

z = arg min s2 As z - u
z Z

U

+ [z ] ,

(3.6)

Z = (h/s ), U = (h ).
Применение метода регуляризации к задаче суперразрешения

Для задачи суперразрешения основной проблемой являются ошибки вычисления векторов движения, прич?м решения систем уравнений (2.11), (2.12) являются неустойчивыми к малым изменениям векторов движения. Также нельзя априори определить, являются ли системы уравнений (2.11), (2.12) совместными. Приближ?нное решение задачи суперразрешения ищется в виде решения задачи минимизации
K

f = arg min
f F k =1

h2 k A f -u Kh

U

+ [f ] .

(3.7)

в непрерывной постановке и
K

z = arg min
z Z k =1

s2 k A z-u Ks

U

+ [z ] .

(3.8)

в дискретной постановке [3133]. Методы, основанные на минимизации (3.7) и (3.8) не являются регуля-


42

ризирующими, тем не менее, их практическое применение да?т хорошие субъективные результаты.
3.2 Оптимизация параметров регуляризирующих методов

Рассомтрим выбор стабилизатора и норм пространств F , Z , U для задач повышения разрешения изображений и суперразрешения и проверим выполнение условий того, что методы (3.5) и (3.6) являются регуляризирующими при данном выборе стабилизаторов и норм пространств.
Выбор норм пространств и стабилизатора

Для задач повышения разрешения изображений мы используем нормы

L1 и L2 : f
1

Rx Ry

=
0 0 Rx Ry

|f (x, y )|dxdy ,

f
и их дискретные аналоги

2 2

=
0 0

|f (x, y )|2 dxdy

N

x

Ny

z

1

=
i=0 j =0 N
x

|zi,j |,
N
y

z

2 2

=
i=0 j =0

|zi,j |2 .

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


43

построении регуляризирующих методов повышения разрешения изображений мы рассматриваем следующие варианты выбора стабилизатора: 1) Квадратичный стабилизатор f 2) Функционал полной вариации. 3) Функционал билатеральной полной вариации, представляющий собой адаптацию функционала полной вариации для задачи повышения разрешения изображений в дискретном случае.
Квадратичный стабилизатор

2 2

.

Квадратичный стабилизатор [6] являются одними из наиболее часто используемых в методах регуляризации стабилизаторов:
Rx R
y

[f ] = f

2 2

=
0 0

2f 2f + x2 y 2

2

dxdy .
2 2

Уравнение Эйлера для задачи минимизации (3.5) при [f ] = f записывается в виде:

(h2 A A + )f = A u,

(3.9)

где сопряж?нный к A оператор A определяется из условия (Af , u) =

(f , A u), для любых f F, u U . В силу линейности оператора A
сопряж?нный к нему оператор A будет линейным, а значит и оператор

(h2 A A + ) из (3.9) также будет линейным.
Это означает, что при выборе квадратичного стабилизатора и квадратичной нормы результату повышения разрешения изображений присущи артефакты линейных методов, а именно при больших возникает эффект размытия, а при малых появляется эффект Гиббса, проявляющийся в виде ореолов возле контуров объектов.


44
Функционал полной вариации

Важнейшей частью изображений являются контуры (edges). Алгоритмы повышения разрешения изображений не должны удалять, либо добавлять контуры. Так, в линейных методах ресамплинга контуры теряются при размытии изображения, а новые контуры добавляются вследствие ложного оконтуривания и алиасинга. Для контроля этих артефактов можно использовать функционал полной вариации [34]. Полная вариация дифференцируемой функции f (x) L(R2 ) определяется в виде
Rx R
y

f
где

V

=
0 0

- | f (x, y )|dxdy ,

- | f (x, y )| =

f (x, y ) x

2

+

f (x, y ) y

2

1/2

.

Для недифференцируемых функций полная вариация вычисляется в виде
Rx R
y

f

V

= lim

h0 0 0

f (x + h, y ) - f (x, y ) h
2 1/2

2

+

f (x, y + h) - f (x, y ) + h

dxdy .

В дискретном случае для z (h ) полная вариация определяется как
N
x

Ny

z

V

=h
i=0 j =0

(zi

+1,j

- zi,j )2 + (zi,j

+1

- zi,j )2

1/2

.

С вычислительной точки зрения более эффективно использовать определение полной вариации в виде
N
x

N

y

z

TV

=h
i=0 j =0

(|zi

+1,j

- zi,j | + |zi,j

+1

- zi,j |) .

(3.10)


45

Для него справедливо

z

V

z

TV





2z

V

,

так как a2 + b2 (|a| + |b|)2 2(a2 + b2 ). Важнейшим свойством полной вариации является е? взаимосвязь с осцилляциями Гиббса, возникающими при действии идеального фильтра высоких частот [34]. Пусть f = f g , где g ядро идеального частотного фильтра со значениями преобразования Фурье

g (w) = 1[- ^
Если f (x) L2 (R2 ), то lim то f
V

, ]

1, |w| , = 0, |w| > .
2

+

f -f

= 0. Однако если f

V

< +,

= + для любого 0 < < . Таким образом, для предотвра-

щения возникновения эффекта Гиббса можно наложить дополнительное условие на ограниченность полной вариации.
Теорема 3.

Метод, основанный на минимизации (3.5) является реV

гуляризирующим для стабилизатора [f ] = f

в пространствах F =
1

L1 (R2 ), U =

1

(h ) пространствах L(R2 ) и (h ) с нормами L1 и

соответственно.
Доказательство.

Для доказательства теоремы достаточно показать,

что множество

MC = {f F : f

V

C, Ah f - u

U

}

является компактным в F для любого C > 0. Из ограниченности полной вариации следует ограниченность разности между наибольшим

sup f =

sup
[0,Rx ]Ч[0,Ry ]

f (x, y )


46

и наименьшим

inf f =
значениями функции f :

[0,Rx ]Ч[0,Ry ]

inf

f (x, y )

sup f - inf f C.
Оператор Ah не увеличивает разность между наибольшим inf f и наименьшим sup f значениями функции f , поэтому

inf f [Ah f ]i,j sup f ,
откуда следует

(i, j ) h ,

0 [Ah f ]i,j - inf f sup f - inf f C,

(i, j ) h .

При вычислении нормы неравенство сохраняется

Ah f - inf f

U

=
i,j

([Ah f ]i,j - inf f )

(sup f - inf f )(Rx + 1)(Ry + 1) C (Rx + 1)(Ry + 1).
С другой стороны

| inf f |(Rx + 1)(Ry + 1) = inf f
откуда следует, что

U

Ah f - inf f

U

+ Ah f

U

,

| inf f | C +
Из условия Ah f - u образом
U

Ah f U . (Rx + 1)(Ry + 1) u
U

получается Ah f

+ . Таким

| inf f | C +

u U + = C1 . (Rx + 1)(Ry + 1)

Из sup f - inf f C получается неравенство для sup f :

| sup f | C1 + C.


47

Отсюда следует, что значения функции f ограничены константой

C2 = max(C1 , C1 + C ) = C1 + C,
значит множество MC является ограниченным. Доказательство компактности MC производится с помощью теоремы Арцела [35], для чего строится конечная -сеть множества MC . Множество [0, Rx ] Ч [0, Ry ] разбивается на N Ч N прямоугольников

P P

i,j

=

N ,j

P P

i,N

N ,N

(i - 1)Rx , N (N - 1)Rx = N (i - 1)Rx , = N (N - 1)Rx = N

iRx N , Rx iRx N , Rx

(j - 1)Ry j Ry , N N (j - 1)Ry j Ry , Ч N N (N - 1)Ry Ч , Ry N (N - 1)Ry Ч , Ry N Ч

, , , ,

1 i, j < N , 1 j < N, 1 i < N,

и записывается функционал полной вариации в виде суммы полных вариаций на этих прямоугольниках:
N N

f

V

=
i=1 j =1

f

V ,Pi,j

.

Рассматривается кусочно-постоянная функция

g (x, y ) =

(x ,y )P

inf

f (x , y ),
i,j

для (x, y ) Pi,j .
V ,Pi,j

Для не? справедливо f - g

V ,Pi,j

=f

на каждом отдельном


48

прямоугольнике Pi,j . Для f - g
Rx R
y

F

на вс?м множестве выполняется:

f -g
N

F

=
0 N 0

|f (x, y ) - g (x, y )|dxdy =

=
i=1 j =1 P N N
i,j

|f (x, y ) - g (x, y )|dxdy =

=
i=1 j =1 P
i,j

f (x, y ) - inf f (x, y ) dxdy
P
i,j

Rx Ry N2 Rx Ry N2
N

N

N

sup f (x, y ) - inf f (x, y )
i=1 j =1 N P
i,j

P

i,j

f
i=1 j =1

V ,Pi,j

=

Rx Ry f N2

V

=

Rx Ry C . N2

Значение N выбирается так, чтобы

Rx Ry C . N2 2
Далее определим конечное множество значений

C2 m , m = 0, +1, . . . + M , M
которые могут принимать кусочно-постоянные функции. При выборе M таком, чтобы

Rx Ry C2 , 2M + 1 2

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

Таким образом, нами построена -сеть из функций, покрывающих MC . Значит, множество MC является компактным, а метод минимизации (3.2) со стабилизатором [f ] = f на.
Следствие 1.

V

регуляризирующим. Теорема доказа-

Теорема справедлива и для пространств F = L2 (R2 ),


49

U= f
TV

2

(h ) с квадратичной нормой, а также для стабилизатора [f ] =

.

При использовании в качестве стабилизатора функционала полной вариации удобно использовать регуляризирующий метод в постановке III:

fC = arg min
f
V

C

Ah f - u

U

.

(3.11)

В отличие от постановок I и II, где регуляризирующие параметры и

априори неизвестны, есть возможность оценить значение C вследствие
связи полной вариации с суммарной длиной линий уровня функции f . Обозначим

t = {(x, y ) R : f (x, y ) > t}.
Если f (x, y ) непрерывна, то граница t множества t является линией уровня всех точек (x, y ) таких, что f (x, y ) = t. Пусть H 1 (t ) длина t . Справедлива теорема [36], что если f
+ V

< +, тогда

f

V

=
-

H 1 (t )dt.

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

C= u

V

.


50
Функционал билатеральной полной вариации

В дискретном случае использование (3.10) в качестве полной вариации да?т лишь приближ?нную информацию о полной вариации изображения. Для повышения точности используется функционал билатеральной полной вариации [1]:
p

z

BT V

=h
s,t=-p



s,t i,j

|zi+

s,j +t

- zi,j | ,

вычисляющий взвешенную сумму разностных производных на изображении z . Параметр p ограничивает количество направлений. Чем больше значение p, тем выше точность метода, но и выше его вычислительная сложность. Параметр s,t зада?т веса для каждого из направлений. В [1] используются значения
s,t

= 0.8|

s|+|t|

.

Анализ выбора стабилизатора и нормы для задачи повышения разрешения изображений

Регуляризирующий метод повышения разрешения изображений был протестирован на базе из 124 фотографий, включающих в себя изображения природы и зданий [37]. Для оценки качества метода эталонные изображения из базы были уменьшены в некоторое число раз, используя модель (2.7), затем увеличены с помощью описанного регуляризирующего метода до исходного размера и сравнены c эталонными изображениями с помощью метрики структурного подобия S S I M [38]:

S S I M (u, v ) =

(2чu чv + c1 )(2uv + c2 ) , 2 2 (ч2 + ч2 + c1 )(u + v + c2 ) u v

(3.12)

2 2 где чu и чv средние значения u и v соответственно, u и v диспер-

сии,

uv

ковариация u и v . Значения c1 = (0, 01L)2 , c2 = (0, 03L)2 , где

L динамический диапазон интенсивностей. Для компьютерных изображений с 8-битным значением интенсивности используется L = 255. Данная метрика хорошо коррелирует с восприятием человека [38]. Чем


51

больше значение метрики, чем ближе изображения друг к другу. При равенстве изображений значение S S I M равно 1. Был провед?н экспериментальный анализ выбора наилучшей комбинации нормы пространств Z , U , F и стабилизатора для увеличения изображений в 2, 3 и в 4 раза с использованием данной метрики. Для каждого изображения значение параметра регуляризации выбиралось таким, чтобы значение S S I M было максимальным. В качестве критерия оценки качества результатов для конкретного выбора нормы и стабилизатора использовалось среднее значение S S I M по всей базе при наилучшем выборе . Также анализ был провед?н для случая зашумл?нных изображений низкого разрешения. Был использован шум с нормальным распределением со стандартным отклонением 10.
``` ``` ``` Норма Стабилизатор ```` `

Увеличение

в 2 раза

в 3 раза
2

в 4 раза
2

ћ

1

ћ

ћ

1

ћ

ћ

1

ћ

2

z z TV z B T V (p = 1) z B T V (p = 2) z B T V (p = 3) Среднее

2 2

0,9849 0,9849 0,9850 0,9839 0,9833 0,9844

0,9852 0,9844 0,9845 0,9829 0,9820 0,9838

0,9732 0,9731 0,9734 0,9736 0,9740 0,9735

0,9738 0,9731 0,9736 0,9737 0,9739 0,9736

0,9652 0,9654 0,9658 0,9658 0,9659 0,9656

0,9653 0,9649 0,9654 0,9655 0,9655 0,9653

Таблица 3.1: Средние значения S S I M при использовании различных норм пространств и стабилизаторов при повышении разрешения незашумл?нных изображений в 2, 3 и в 4 раза. Увеличение в 2 раза в 3 раза
2

``` ``` Норма `` Стабилизатор ````` `

в 4 раза
2

ћ

1

ћ

ћ

1

ћ

ћ

1

ћ

2

z 2 2 z TV z B T V (p = 1) z B T V (p = 2) z B T V (p = 3) Среднее

0,9614 0,9692 0,9687 0,9683 0,9683 0,9672

0,9704 0,9732 0,9732 0,9711 0,9687 0,9713

0,9454 0,9538 0,9556 0,9548 0,9546 0,9528

0,9565 0,9603 0,9607 0,9600 0,9596 0,9594

0,9262 0,9451 0,9458 0,9448 0,9444 0,9413

0,9449 0,9492 0,9499 0,9495 0,9491 0,9485

Таблица 3.2: Средние значения S S I M при использовании различных норм пространств и стабилизаторов при повышении разрешения зашумл?нных изображений в 2, 3 и в 4 раза.


52

В таблицах 3.1 и 3.2 приведены средние значения S S I M для рассмотренных норм пространств и стабилизаторов для различных коэффициентов масштабирования. На основе полученных результатов можно сделать следующие наблюдения: 1. Для случая незашумл?нных изображений нормы

ћ

1

и

ћ

2

дают

практически одинаковые значения S S I M с незначительным преимуществом ћ 1 , а для зашумл?нных наилучшие результаты показывает норма

ћ 2 . Отметим, что были рассмотрены и другие варианты выбора нормы ћ
q

, q > 1, q = 2, но они не дали удовлетворительных результатов. ћ
2

2. Для коэффициента масштабирования s = 2 и нормы

наилучV

шие результаты достигаются при использовании стабилизаторов z

,

z

TV

,z

BT V

(p = 1). ћ
2

3. Для коэффициентов масштабирования s = 3, 4 и нормы

все

рассмотренные стабилизаторы дают примерно одинаковые результаты. С вычислительной точки зрения эффективнее использовать стабилизаторы с наименьшее число операций. Такими стабилизаторами являются

z

TV

иz

BT V

(p = 1).

В результате обобщения этих наблюдений формулируется вывод: для повышения разрешения изображений с помощью разработанного регуляризирующего метода предпочтительнее использование нормы стабилизатора z
BT V

ћ

2

и

при p = 1.

Анализ выбора стабилизатора и нормы для задачи суперразрешения

Аналогичный анализ был произвед?н нами для задачи суперразрешения на изображениях из той же базы. Для каждого эталонного изображения высокого разрешения z было сгенерировано 32 изображения низкого разрешения uk в соответствии с моделью (2.10) с операторами движения Tk , состоящими из случайного сдвига в пределах +5 пикселей по горизонтали и вертикали и случайного поворота в пределах +10 градусов. Основной проблемой в задаче суперразрешения являются ошибки вы-


53

числения операторов движения Tk , поэтому после моделирования изображений низкого разрешения к операторам Tk добавлялась ошибка. Для четверти операторов Tk ошибка представляла собой равномерно распредел?нное случайное смещение каждого пикселя в интервале [-n0 , n0 ], для остальных операторов ошибка находилась в интервале [-n0 /4, n0 /4]. Тестирование методов проводилось как в условиях недостатка данных, когда число неизвестных в уравнении (2.12) больше числа уравнений, то есть при s2 > K , так и в условиях избытка данных s2 < K . Увеличение производилось в 2, 3 и 4 раза при числе изображений низкого разрешения s2 /2, s2 и 2s2 . Рассматривалось три уровня шума: n0 =

0.1, 0.5, 1.
Полученные результаты оказались близки к результатам для задачи повышения разрешения одиночного изображения.
3.3 Постобработка результатов повышения разрешения и суперразрешения

В предыдущем разделе была показана связь функционала полной вариации с эффектом Гиббса. Появление эффекта Гиббса приводит к значительному увеличению значения полной вариации. В связи с этим ставится задачу подавления эффекта Гиббса на изображении z в виде задачи проектирования изображения z на множество изображений с ограниченной полной вариацией

zC = prC z = arg

min
z: z
TV

C

z - z 2,

(3.13)

где параметр C отвечает за силу подавления эффекта Гиббса. При подавлении эффекта Гиббса после повышения разрешения изображений параметр C в (3.13) может быть оцен?н через изображение низкого разрешения u [39, 40]. Полная вариация равна сумме длин линий уровня [34]. Так как алгоритм повышения разрешения изображений не


54

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

C= u

TV

.

Так как изображение может состоять из достаточно разнородных участков, глобальное подавление полной вариации может привести к тому, что в одних участках помимо эффекта Гиббса могут быть подавлены мелкие детали, тогда как в других участках эффект Гиббса не будет полностью подавлен. Для того, чтобы учесть локальные особенности изображений, можно использовать блочный подход [41]. В этом случае изображение

z разбивается на перекрывающиеся блоки, подавление эффекта Гиббса
производится независимо для каждого блока, при этом в качестве параметра C используется значение полной вариации соответствующих блоков на изображении u. Затем из блоков восстанавливается изображение высокого разрешения. Пример подавления эффекта Гиббса после повышения разрешения изображений привед?н на рис. 3.2.


55

Изображение низкого разрешения T V = 8, 6515

Увеличение изображения в 2 раза T V = 17, 24, RE = 8, S S I M = 0, 9779

Глобальный метод T V = 8, 927, RE = -0, 7, S S I M = 0, 9827

Блочный метод T V = 8, 736, RE = -1, 5, S S

IM

= 0, 9812

Рис. 3.2: Пример подавления эффекта Гиббса после повышения разрешения изображений.


56

Глава 4
Численные методы, используемые для решения поставленных задач

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

Для нахождения решения задач повышения разрешения изображений (3.6), суперразрешения (3.8), подавления эффекта Гиббса (3.13) требуется произвести минимизацию функционала вида

z = arg min J (z ).
z Z

В нашем случае используется функционал вида J (z ) = Az - u +

[z ].
При использовании регуляризирующих методов со стабилизатором

[z ] в виде полной вариации или билатеральной полной вариации регуляризирующий функционал не является дифференцируемым. Для того, чтобы иметь возможность использовать градиентные методы недифференцируемых функционалов, используется понятие субградиента [42].
Определение.

Пусть функционал J (z ) определ?н на множестве Z из
субградиентом

E n . Вектор c = c(z0 ) E n называется

функционала J (z )


57

в точке z0 , если условие

J (z ) J (z0 ) + (c(z0 ), z - z0 )
выполняется для любого z E n . Множество всех субградиентов функционала в точке z0 называется
субдифференциалом

и обозначается через J (z0 ).

Субградиент обладает следующими свойствами [42]: 1. Если J (z ) выпуклый функционал, а Z выпуклое множество, то субдифференциал является непустым для любого z0 Z . 2. Если функционал J (z ) является дифференцируемым в точке

z0 Z , то его градиент J (z0 ) является субградиентом, а если z0 внутренняя точка Z , то его субдифференциал содержит только один вектор, равный градиенту J (z0 ). Для минимизации J (z ) используется субградиентный метод [43]

z
где g
(k )

(k +1)

=z

(k )

- k g ( k ) ,
(k )

один из субградиентов функционала J (z ) в точке z = z

.

В [43] рассматриваются следующие стратегии выбора шага k и для них доказывается сходимость субградиентного метода: 1. Постоянный коэффициент шага k = , либо постоянная длина шага k =
g
(k )

. В этом случае сходимость является -оптимальной [43],
(k )

то есть начиная с некоторой итерации все приближ?нные решения z от . 2. Последовательность k , удовлетворяющая условиям
2 k k =1

находятся в -окрестности точного решения z , при этом зависит только

k 0,

< ,
k =1

k = .

3. Несуммируемая убывающая последовательность коэффициентов ша-


58

га



k 0,

k

lim k = ,
k =1

k = ,

либо несуммируемая убывающая последовательность длин шага


k 0,

k

lim k = ,
k =1

k = ,

k k = (k) . g
Скорость сходимости промежуточного решения z может быть оценена с помощью формулы
(k )

субградиентного

метода к точному решению z при наилучшем выборе k = (R/G)/ k



z
где R = z
(0)

(k )

-z



RG , k

(4.1)

- z , g

(k )

G.

Вычисление субградиента регуляризирующего функционала

Рассмотрим задачу минимизации функционала вида

J (z ) = Az - u

n Z

+ [z ].

Для функции (x) = |x| субградиент при x = 0 полагается равным нулю: (x)|x
=0

= 0. В этом случае получаются следующие субградиенn Z

ты для невязки Az - u

:

( Az - u 1 ) = A sign(Az - u), = Az - u 1 ћ A sign(Az - u), 1 ( Az - u 2 ) = A (Az - u), Az - u 2 Az - u
2 1



Az - u

2 2

= 2A (Az - u),

где sign u оператор, действующий из U в U и представляющий из себя


59

покомпонентное применение функции sign:

[sign u]i,j

1, ui,j > 0, = 0, ui,j = 0, -1, ui,j < 0.

Для стабилизатора z 2 : 2

z

2 2

= 2 z .
TV

Субградиенты стабилизаторов z зуя их представление в виде

иz

BT V

можно найти, исполь-

z z

TV

= h ( Sx z - z
p

1

+ Sy z - z 1 ) ,
1

BT V

=h
s,t=-p

st s,t Sx Sy z - z

,

s t где Sx и Sy операторы сдвига изображения по горизонтали и по вер-

тикали на s и t пикселей соответственно.
t s Для операторов Sx и Sy субградиент выглядит следующим образом s Sx z - z t Sy z - z - s = Sx s sign(Sx z - z ), - t = Sy t sign(Sy z - z ). s2 -1

1 1

Для операторов Ah = Dh H0 h , As = Ds H0

и сопряж?нные

операторы также могут быть выписаны в явном виде.
Оператор H = H как в непрерывном, так и в дискретном случае. В


60

самом деле,

(H f , g ) = (x - 1 exp - = 2 2 1 (x - = exp - 2 2 1 = f (x , y ) exp 2 2 1 exp = f (x , y ) 2 2 =

[H f ](x, y )g (x, y )dxdy = f (x , y )dx dy g (x, y )dxdy =

x )2 + (y - y )2 2 2 x )2 + (y - y )2 2 2 (x - x)2 + ( - 2 2 (x - x)2 + ( - 2 2

f (x , y )g (x, y )dx dy dxdy = y - y )2 y - y )2 g (x, y )dx dy dxdy = g (x, y )dxdy dx dy =

f (x , y )[H g ](x , y )dx dy = (f , H g ).

Аналогично доказывается, что = .
Далее строится v = Dh u:

(Dh f , u) =
i,j

f (ih, j h)ui,j =

f (x, y )
i,j

ui,j (x - ih) (y - j h)dxdy = (f , Dh u),

откуда следует, что

v (x, y ) =
i,j

ui,j (x - ih) (y - j h).

В дискретном случае для целого s
[Ds u]i,j

ui/s,j /s , если i и j кратны s, = 0, иначе.

Практическое применение субградиентного метода для задачи ресамплинга изображений

Субградиентный метод имеет низкую скорость сходимости (4.1). При программной реализации субградиентного метода особое внимание необ-


61

ходимо уделить выбору начального приближения z циентов k .

(0)

и выбору коэффи-

При повышении разрешения изображений основной задачей является не нахождение точного решения задачи минимизации, а нахождение приближ?нного решения приемлемого качества за минимальное количество итераций. Анализ, провед?нный в [21], показал, что для повышения разрешения изображения можно использовать шаг вида

k =

N k , g (k) 2

k = 0 q k ,

(4.2)

где N число пикселей изображения. Основным параметром здесь является q , отвечающий за скорость уменьшения длины шага. Субградиентный метод останавливается, когда k < . При этом количество итераций является фиксированным и зависит только от 0 , q и

. С ростом q точность метода должна повышаться. Однако анализ зависимости метрики значение
SSI M SSI M

от q показал, что, начиная с некоторого q = q0 ,

переста?т повышаться. Это связано с тем, что точное

решение задачи минимизации отличается от эталонного изображения. Поэтому для получения изображения наилучшего качества достаточно использовать значение q = q0 . Анализ выбора зависимости результата от начального приближения показал, что чем больше значение
SSI M

для начального приближения

z (0) , тем больше значение

SSI M

для результата субградиентного метода.

Однако с ростом q зависимость от начального приближения снижается. В общем случае нельзя задать параметры 0 , q и начального приближения z
(0)

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


62
4.2 Метод минимизации квадратичных функционалов на множестве функций с ограниченной полной вариацией

Рассмотрим задачу минимизации квадратичного функционала

F (z ) = Az - u 2 , 2
где A линейный непрерывный оператор с det A = 0, на множестве функций с ограниченной полной вариацией

MC = {z : z

TV

C },

используемую в задачах подавления эффекта Гиббса (3.13). Для нахождения е? решения используется метод минимизации, основанный на методе условного градиента [44].
4.2.1 Метод условного градиента

Метод условного градиента применяется для минимизации функционала

F (z ) на множестве M выпуклом многограннике с вершинами {Tj }.
В методе условного градиента минимизирующая z ная z
(k ) (k )

и вспомогатель-

последовательности строятся согласно следующему алгоритму:
(0)

1. Выбирается некоторое начальное приближение z 2. Пусть известно приближение z ное приближение z
(k ) (k )

M.

, k = 0, 1, . . .. Тогда вспомогатель-

полагается равным решению задачи:
(k )

(F (z

), z

(k )

) = min(F (z
z M

(k )

), z ).

Так как функционал (F (z выпуклым, то z
(k )

(k )

), z ) является линейным, а множество M

следует искать на границе множества M . Для слу-

чая выпуклого многогранника достаточно перебрать все его вершины

Tn .


63

3. Пусть значение z

(k )

найдено. Далее z

(k +1)

строится по формуле
(k )

z

(k +1)

=z

(k )

+ k (z

(k )

-z

),

где k [0, 1] является решением задачи минимизации

F (z
4.2.2

(k +1)

) = F (z

(k )

+ k (z

(k )

-z

(k )

)) = min F (z
[0,1]

(k )

+ (z

(k )

-z

(k )

)).

Применение метода условного градиента для задачи подавления эффекта Гиббса

Метод условного градиента неприменим в явном виде к задаче минимизации

z1 = arg min F (z ),
z M
C

(4.3)

так как множество MC не является ограниченным.
Одномерный случай

Для того, чтобы метод условного градиента был применим, рассматривается задача минимизации

z2 = arg min F (z ),
z M

(4.4)

где M некоторый выпуклый многогранник, и строится такое множество M MC , при котором задачи (4.3) и (4.4) будут эквивалентны. Решения задач (4.3) и (4.4) существуют и единственны в силу строгой выпуклости минимизируемого функционала и выпуклости множеств M и MC . Для доказательства эквивалентности достаточно показать, что

z1 = z2 . Это равенство выполняется в случае, если z1 M .
Используя рассуждения, аналогичные рассуждениям в теореме 3 из раздела 3.2, можно показать, что значение Az1
2

ограничено некоторой
-1

константой C1 , зависящей от u и C . Обратный оператор A

существу-

ет и ограничен, так как det A = 0, а пространство конечномерно. Отсюда


64

следует, что решение (4.3) ограничено:

z1
и для него справедливо

2

C2 = C1 A

-1

|z1,i | B ,

i = 1, 2, . . . , n,

где B константа, зависящая от A, u и C . Таким образом, получено, что z1 принадлежит ограниченному множеству

Mz = {z Z : max |zi | B }.
i

Выпуклый многогранник M , такой, что Mz M MC , может быть построен в явном виде. В [44] рассмотрено построение вершин много0 гранника для множества MC , состоящего из элементов z Z , удовле-

творяющих условию

|z2 - z1 | + |z3 - z2 | + . . . + |zn - zn-1 | C, zn = 0.
Для данного множества вершины T имеют вид
(j )

, j = -(n-1), . . . , -1, 1, . . . , n-1

T

(j ) i

C, i j, = 0, i > j, = -Ti ,
(j )

j = 1, 2, . . . , n - 1, j = 1, 2, . . . , n - 1.

Ti

(-j )

0 От условия zn = 0 можно избавиться, используя множество Mz , со-

стоящее из элементов z Z , удовлетворяющих условию

z1 = z2 = . . . = zn , |zn | B .
0 Множество M = MC Mz состоит из всех элементов z , для которых


65

z

TV

C и |zn | B , и представляет собой призму с вершинами C - B , i j, (j ) Ti = j = 1, 2, . . . , n - 1, -B , i > j, C + B , i j, (j +n-1) Ti = j = 1, 2, . . . , n - 1, B , i > j, Ti
(-j )

= -T

(j ) i

,

j = 1, 2, . . . , 2n - 2.

Очевидно, что M MC и Mz M , а значит, задача (4.4) эквивалентна задаче (4.3). Таким образом, нахождение решения задачи (4.3) сведено к нахождению решения задачи (4.4), для которой метод условного градиента применим.
Двухмерный случай

В двухмерном случае построение многогранника является достаточно затруднительным и неэффективным с вычислительной точки зрения. Метод условного градиента применяется только для задачи подавления эффекта Гиббса после повышения разрешения изображения, в которой изображение разбивается на перекрывающиеся блоки небольшого размера. При использовании метода условного градиента вместо квадратных блоков используются горизонтальные и вертикальные блоки с высотой (шириной) в один пиксель. Пусть QH (z , u) оператор подавления эффекта Гиббса на изображении высокого разрешения z при известном изображении низкого разрешения u с помощью метода условного градиента с горизонтальными блоками, QV (z , u) с вертикальными блоками. Тогда результатом подавления эффекта Гиббса мы считаем изображение [41]

z=

1 (QH (QV (z , u), u) + QV (QH (z , u), u)) . 2


66

4.3

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

До сих пор не была затронута задача выбора параметра регуляризации

для задач повышения разрешения изображений (3.5), (3.6).
Это связано с тем, что в общем случае параметр регуляризации выбрать нельзя, так как значение ошибки приближ?нно заданного изображения низкого разрешения u неизвестно. Для каждого изображения параметр регуляризации, при котором будут наилучшими значения метрик оценки качества изображений, будет различаться. Тем не менее, можно осуществить выбор параметра регуляризации, при котором качество результата повышения разрешения изображений будет близко к наилучшему. Используются две стратегии выбора параметра регуляризации: предварительное вычисление параметра регуляризации на основе изображений из базы эталонных изображений и адаптивное вычисление параметра регуляризации на основе анализа самого изображения.
Предварительное вычисление параметра регуляризации

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

между эта-

лонным изображением высокого разрешения и результатом повышения

, будет


67

близко для всех изображений из одного класса. Таким образом, достаточно вычислить наилучшее значение параметра регуляризации для базы изображений определ?нного класса, а затем использовать полученное значение при повышении разрешения всех изображений этот класса. В таблицах 4.1 и 4.2 приведены средние геометрические значения наилучших параметров регуляризации при использовании различных норм пространств и стабилизаторов при повышении разрешения изображений в 2, 3 и в 4 раза, вычисленные для базы, состоящей из 124 фотографий природы и зданий.
``` ``` Норма `` Стабилизатор ````` `

Увеличение

в 2 раза

в 3 раза
2

в 4 раза
2

ћ

1

ћ

ћ

1

ћ

ћ

1

ћ

2

z z z

B B B

z 2 2 z TV T V (p = 1) T V (p = 2) T V (p = 3)

0,00328 0,04936 0,04720 0,01124 0,01060

0,01576 0,58840 0,28412 0,11164 0,04492

0,01215 0,21240 0,10665 0,03636 0,02205

0,07164 0,60741 0,37071 0,08343 0,04968

0,00624 0,26688 0,13152 0,03024 0,01856

0,01696 0,35920 0,23600 0,06928 0,03696

Таблица 4.1: Средние геометрические значения наилучших параметров регуляризации при использовании различных норм пространств и стабилизаторов при повышении разрешения незашумл?нных изображений в 2, 3 и в 4 раза. Увеличение в 2 раза в 3 раза
2

``` ``` Норма `` Стабилизатор ````` `

в 4 раза
2

ћ

1

ћ

ћ

1

ћ

ћ

1

ћ

2

z z z

B B B

z 2 2 z TV T V (p = 1) T V (p = 2) T V (p = 3)

0,01556 0,66348 0,35892 0,09628 0,05260

0,48404 6,75636 3,88072 1,32696 0,78588

0,03645 0,73980 0,50553 0,12663 0,05805

1,92033 10,24659 5,24529 1,39761 0,68310

0,02912 1,48672 0,79424 0,17952 0,08016

2,15136 14,57408 6,79312 1,80144 0,84352

Таблица 4.2: Средние геометрические значения наилучших параметров регуляризации при использовании различных норм пространств и стабилизаторов при повышении разрешения зашумл?нных изображений в 2, 3 и в 4 раза.

Адаптивное вычисление параметра регуляризации

Альтернативой предварительного вычисления параметра регуляризации является вычисление параметра регуляризации на основе анализа самого


68

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


69

Литература

1. Fast and robust multi-frame super-resolution / S. Farsiu, D. Robinson, M. Elad, P. Milanfar // 2.
IEEE Transactions on Image Processing

. 2004. Vol. 13, no. 10.

Pp. 13271344.
Pham T. Q., van Vliet L. J., Schutte K.

Robust super-resolution by minimizing a . 2008.

gaussian-weighted 3.

2

error norm //

Journal of Physics: Conference Series

Vol. 124, no. 1. Pp. 1237.
Котельников В. А.

О пропускной способности эфира и проволоки в электросвя-

зи // Всесоюзный энергетический комитет. Материалы к I Всесоюзному съезду по вопросам технической реконструкции дела связи и развития слаботочной промышленности. М.: Управление связи РККА, 1933. С. 119. 4. 5.
Zenzo S. D.

A note on the gradient of a multi-image //

Journal on Computer Vision,

Graphics, and Image Processing

. 1986. Vol. 33, no. 1. Pp. 116125. Color plane interpolation using
Image Processing

Gunturk B. K., Altunbasak Y., Mersereau R. M.

alternating pro jections // Vol. 11. Pp. 9971013. 6. 7. 8. 9. 10. 11.

IEEE

Transactions

on

. 2002.

Тихонов А. Н., Арсенин В. Я.

Методы решения некорректных задач. М.: НауImage interpolation and resampling // Handbook of

ка, 1979. С. 288.
Thevenaz P., Blu T., Unser M.

medical imaging. Academic Press, Inc., 2000. Pp. 393420.
Blu T., Thevenaz P., Unser M. on Image Processing

Linear interpolation revitalized // Interpolation revisited //

IEEE Transactions

. 2004. Vol. 13, no. 5. Pp. 710719.
IEEE Transactions on

Thevenaz P., Blu T., Unser M. Medical Imaging Turkowski K.

. 2000. Vol. 19, no. 7. Pp. 739758.

Filters for common resampling tasks // Graphics gems. Academic Recursive implementation of the gaussian lter //

Press, Inc., 1990. Pp. 147165.
Young I. T., van Vliet L. J. Processing Signal

. 1995. no. 44. Pp. 139151.


70

12.

Lee Y. J., Yoon J.

Nonlinear image upsampling method based on radial basis function . 2010. Vol. 19, no. 10.

interpolation // Pp. 26822692. 13.

IEEE Transactions on Image Processing

Sun J., Xu Z., Shum H.-Y.

Image super-resolution using gradient prole prior //

IEEE

Conference on Computer Vision and Pattern Recognition (CVPR2008)

. 2008.

Pp. 18. 14. High accuracy WADI image interpolation with local gradient features / S. Yuan, M. Abe, A. Taguchi, M. Kawamata // Proceedings of 2005 International Symposium on Intelligent Signal Processing and Communication Systems. 2005. Pp. 8588. 15.
Li X., Orchard M. T. Processing

New edge-directed interpolation //

IEEE Transactions on Image

. 2001. Vol. 10, no. 10. Pp. 15211527. Content-adaptive video up-scaling for high-

16.

Leitao J. A., Zhao M., de Haan G.

denition displays //
'2003

Proceedings of Image and Video Communications and Processing

. 2003. Vol. 5022. Pp. 612622. Resampling and reconstruction with fractal interpolation . 1998. Vol. 5, no. 9. Pp. 228230.
IEEE

17.

Price J. R., Hayes M. H.

functions // 18.

IEEE signal processing letters

Wohlberg B., de Jager G.

A review of the fractal image coding literature // . 1999. Vol. 8, no. 12. Pp. 17161729.

Transactions on Image Processing

19.

Polidori E., Jean-Luc D.

Zooming using iterated function system //

NATO ASI Conf.

Fractal Image Encoding and Analysis

. 1995. Vol. 5. Pp. 111123.

20. Resolution enhancement of images using fractal coding / M. Gharavi-Alkhansari, R. Denardo, Y. Tenda, T. S. Huang //
Processing Visual Communications and Image

. 1997. Vol. 3024. Pp. 10891100. Image interpolation by super-resolution //

21.

Lukin A. S., Krylov A. S., Nasonov A. V. Proceedings of GraphiCon'2006

. 2006. Pp. 239242.

22.

Aly H. A., Dubois E.

Image up-sampling using total-variation regularization with a
IEEE Transactions on Image Processing

new observation model // no. 10. Pp. 16471659. 23.
Иванов В. К.

. 2005. Vol. 14,

О линейных некорректных задачах //

ДАН СССР

. 1962. Т.

145. С. 270272. 24.
Васин В. В.

О связи некоторых вариационных методов приближенного решения
Математические заметки

некорректных задач // 272.

. 1970. Т. 7, 3. С. 265


71

25. 26.

Rudin L., Osher S. J., Fatemi E.

Nonlinear total variation based noise removal

algorithms //
Irani M.,

Physica D S.

. 1992. no. 60. Pp. 259268.
Proceedings of

Peleg

Super resolution from image sequences //

International Conference on Pattern Recognition (ICPR'90)

. 1990. Vol. 2.

Pp. 115120. 27. 28. 29.
Irani M., Peleg S.

Improving resolution by image registration //

CVGIP: Graphical

models and image processing Ratakonda K., Ahuja N.

. 1991. Vol. 53, no. 3. Pp. 231239.
International

POCS based adaptive image magnication //

Conference on Image Processing

. 1998. Vol. 3, no. 3. Pp. 203207. POCS-based restoration of space-varying . 1994. Vol. 3, no. 4.

Ozkan M. N., Tekalp A. M., Sezan M. I.

blurred images // Pp. 450454. 30. 31.
Танана В. П.

IEEE Transactions on Image Processing

Методы решения операторных уравнений. М.: Наука, 1981. Применение метода суперразрешения
Системы высокой

С. 156.
Насонов А. В., Крылов А. С., Ушмаев О. С.

для биометрических задач распознавания лиц в видеопотоке //
доступности

. 2009. Т. 1. С. 2634. Face image super-resolution from video . 2008.

32.

Krylov A. S., Nasonov A. V., Sorokin D. V.

data with non-uniform illumination // Pp. 150155. 33.

Proceedings of GraphiCon'2008

Насонов А. В., Крылов А. С., Ушмаев О. С.

Развитие методов повышения ка.

чества изображений лиц в видеопотоке // 2009. Т. 3, 1. С. 1928. 34. 35. 36. 37.
Малла С.

Информатика и е? применения

Вэйвлеты в обработке сигналов. М.: Мир, 2005. С. 672. Элементы теории функций и функционального

Колмогоров А. Н., Фомин С. В.

анализа. М.: Наука, 1972. С. 496.
Ziemer W. P.

Weakly Dierentiable Functions. New York: Springer-Verlag, 1989. Edge-preserving nonlinear iterative image

Krylov A. S., Lukin A. S., Nasonov A. V.

resampling method //
(ICIP'09)

Proceedings of International Conference on Image Processing

. 2009. Pp. 385388. .

38. Image quality assessment: from error visibility to structural similarity / Z. Wang, A. Bovik, H. Sheikh, E. Simoncelli // 2004. Vol. 13, no. 4. Pp. 600612.
IEEE Transactions on Image Processing


72

39.

Nasonov A. V., Krylov A. S., Lukin A. S.

Post-processing by total variation quasi. 2007.

solution method for image interpolation // Pp. 178181. 40.
Насонов А. В.

Proceedings of GraphiCon'2007

Программное повышение разрешения и подавление эффекта ГиббСборник аннотаций работ 6-й Курчатовской молод?ж-

са на изображениях //
ной научной школы

. 2008. С. 153.
A. V.

41.

Krylov

A.

S.,

Nasonov

Adaptive total variation deringing method for

image interpolation //
(ICIP'08)

Proceedings of International Conference on Image Processing

. 2008. Pp. 26082611. Численные методы решения экстремальных задач. М.: Наука,

42.

Васильев Ф. П.

1988. С. 552. 43.
Boyd S., Xiao L., Mutapcic A.

Subgradient methods. Stanford University: Lecture

notes of EE392o, 2003. 44. Регуляризирующие алгоритмы и априорная информация / А. Н. Тихонов, А. В. Гончарский, В. В. Степанов, А. Г. Ягола. М.: Наука, 1983. С. 200.