Реконструкция смазанных и зашумленных изображений без использования граничных условий
ИКОНИКА – НАУКА ОБ ИЗОБРАЖЕНИИ
УДК 517.968: 517.983.54: 621.397.3
РЕКОНСТРУКЦИЯ СМАЗАННЫХ И ЗАШУМЛЕННЫХ ИЗОБРАЖЕНИЙ БЕЗ ИСПОЛЬЗОВАНИЯ ГРАНИЧНЫХ УСЛОВИЙ
© 2009 г. В. С. Сизиков, доктор техн. наук; М. В. Римских; Р. К. Мирджамолов Санкт-Петербургский государственный университет информационных технологий, механики и оптики, Санкт-Петербург
E-mail: sizikov2000@mail.ru, romefamily@yandex.ru
Предложен новый подход к проблеме учета интенсивностей, исходящих из точек вне заданных границ изображения. В ряде работ это учитывается путем введения так называемых граничных условий (анти-рефлективных и др.). В данной работе вместо них предлагается прием усечения, не требующий знания интенсивностей за границами изображения. На основе приема усечения изложены прямая и обратная задачи обработки изображения. Для решения обратной задачи реконструкции искаженного изображения использованы метод преобразования Фурье с регуляризацией Тихонова и метод квадратур (также с регуляризацией). Приведены численные иллюстрации, показавшие, что предложенная в работе методика заметно снижает погрешность реконструкции.
Коды OCIS: 100.0100, 100.2000.
Поступила в редакцию 03.12.2008.
Введение
В данной работе рассматривается задача реконструкции искаженных (смазанных или/и зашумленных) изображений [1, 2]. Рассматривается как прямая задача (моделирование искажения или экспериментальное получение реального искаженного изображения), так и обратная задача (реконструкция, восстановление изображения) [2–4].
Основные соотношения и уравнения, описывающие эти задачи в непрерывном виде, известны по работам [2–11]
x+Δ
(1/Δ) ò wy (ξ) dξ = gy (x) + δg,
x
(1)
¥
ò h(x - ξ)wy (ξ)dξ = gy (x) + δg,
-¥
(2)
¥¥
ò ò h(x - ξ,y - η)w(ξ,η)dξdη = g(x,y) + δg. (3)
-¥ -¥
Здесь Δ – смаз, h – функция рассеяния точки (ФРТ), обычно пространственно-инвариантная, w и g – распределение интенсивности по неискаженному и искаженному изображени-
ям соответственно, δg – помеха. В (1) и (2) ось x направлена вдоль смаза, а y играет роль параметра.
Интегральные уравнения (1) и (2) обычно используются в задаче смазывания, а (3) – в задаче дефокусирования, но часто, например, в работах [2, 5] уравнение (3) используется для решения обеих задач. Интегральные уравнения (2), (3) обычно решаются методами инверсной фильтрации, параметрической винеровской фильтрации, регуляризации Тихонова (прямой и итеративный вариант), простой итерации Фридмана, итераций Ландвебера, сопряженных градиентов (CG) и др. При численной реализации методов используются преобразование Фурье (ПФ), квадратуры, итерации и т. д.
Часто используются, например, в работах [9, 11–14] так называемые граничные условия (BCs): нулевые по Дирихле, периодические, копирующие границу, рефлективные (симметричные) относительно границы, анти-рефлективные относительно границы для учета интенсивностей вне заданных границ изображения [12, 13], когда функция w не является финитной, т. е. носитель изображения в действительности шире поля зрения – заданной области w и g. Непрерывным уравнениям (1)–(3) можно поставить в соответствие дискретные выражения (1)–(3) вида
38 “Оптический журнал”, 76, 5, 2009
g = Aw + δg,
(4)
где A – матрица, связанная с ФРТ и с граничными условиями: в случае нулевых условий матрица A есть блок теплицевых матриц, в случае периодических условий – блок циркулянтных матриц, в случае рефлективных относительно границы – сумма блоков теплицевых и ганкелевых матриц [12]. Заметим, что в некоторых случаях искомую функцию w можно считать финитной, например, когда рассматривается космический объект [12, Fig. 4]. Однако гораздо чаще искомая функция w не является таковой, как например, в работе [10, рис. 1], хотя авторы этой работы положили, что w финитна и, по существу, использовали нулевое граничное условие.
Основная цель данной работы – показать, что в случае, когда функция w не является финитной, т. е. w отлична от нуля в более широкой области, чем область поля зрения (FOV), можно адекватно решать как прямую, так и обратную задачу, не используя граничные условия. Данная работа продолжает работы [4, 7, 8, 15, 16].
Прием усечения изображения
О граничных условиях. Рассмотрим прямую задачу моделирования искажения изображения. Пусть нам дано исходное (неискаженное) изображение размером m×n пикселов, описываемое матрицей интенсивностей wm×n. Назовем область Dm×n заданной областью исходного изображения (FOV), а границу B этой области – границей исходного изображения.
В редких случаях (когда, например, исходное изображение является изображением космического объекта) функцию w можно считать финитной, т. е. равной нулю вне области D. В этом случае можно без потери информации моделировать искаженное изображение по формуле типа (4), положив размеры матрицы g, равными размерам матрицы w, и используя при этом граничное условие zero. Аналогично, при решении обратной задачи (реконструкции изображения) можно без потери информации полагать размеры матриц w и g одинаковыми, т. е. считать, что w = 0 вне области D.
Однако авторы данной работы, как и авторы работ [9, 11–14] рассматривают более типичный случай, когда вне области D интенсивности w, вообще говоря, отличны от нуля, но информация об этих интенсивностях отсутствует. В этом случае (при моделировании искажения, а также при решении обратной задачи) обычно используют
граничные условия (BCs), наиболее эффективным из которых считается anti-reflective, и которые являются различными вариантами экстраполяции значений w за границу В. Однако, во-первых, введение BCs является довольно искусственным приемом, а во-вторых, их использование порождает сложный математический аппарат. В данной работе предлагается более простое решение этой проблемы, которое будем называть приемом усечения изображения.
Усечение смазанного изображения. Рассмотрим простейший пример моделирования равномерного прямолинейного горизонтального смазывания изображения (прямая задача). Пусть нам дана матрица исходного изображения wm×n и необходимо смоделировать матрицу смазанного изображения g при смазе Δ и при нулевом угле смаза (горизонтальный смаз). Если формировать смазанное изображение такого же размера, как и исходное изображение, т. е. формировать g шириной n, то соотношение (4) более подробно запишется в виде (при замене интеграла в (1) конечной суммой по пикселам, без помехи)
ågj
(i)
=
Δ
1 +1
i+Δ k=i
wj
(k),
i =1, 2,…,n, j =1, 2,…, m,
(5)
где i (и k) – номер столбца, j – номер строки.
В (5) перед знаком суммы поставлен множитель
1/(Δ + 1) потому, что при Δ = 0 выражение (5)
должно переходить в gj(i) = wj(i). Заметим, что в выражении (1) перед знаком интеграла постав-
лен множитель 1/Δ, чтобы при Δ → 0 соотноше-
ние (1) переходило (без помехи) в wy(x) = gy(x) [4, с. 66].
Из (5) видно, что при i = n – Δ + 1, …, n требу-
ется знание значений wj (1), …, wj(n + Δ), однако значения wj(n + Δ)wj(n + Δ) неизвестны. Чтобы преодолеть этот недостаток информации, в об-
ширном цикле работ [11–14] предложены и ис-
пользуются граничные условия (было бы более
точно назвать их заграничными условиями):
от простейшего нулевого до анти-рефлектив-
ного, в которых делаются различные предпо-
ложения относительно недостающих значений
wj (n + 1), …, wj(n + Δ). Авторы статьи предполагают иное решение
данного вопроса, а именно, вместо соотношений
(5) использовать соотношения (рис. 1)
ågj
(i)
=
Δ
1 +1
i+Δ k=i
wj
(k),
i =1, 2,…, n - Δ, j =1, 2,…, m.
(6)
“Оптический журнал”, 76, 5, 2009
39
размыттие у g размыттие у g
i, k w шириной n
Dmxn (FOV)
j
(
m g шириной n –
B
w шириной n
n смаз
Рис. 1. Схема получения смазанного изображения с усечением.
g шириной n +
Рис. 2. Схема с введением размытых краев у изображения g.
Отличие выражения (6) от (5) состоит в том, что согласно (6) моделируется более узкое смазанное изображение. При этом требуются значения wj (1), …, wj(n), которые известны, и не нужно прибегать к такому приему как граничные условия. Правда, g получается несколько у′же, чем w (на Δ пикселов). Но, во-первых, практически эффект заужения невелик (обычно n ≈ 300÷600, а Δ ∼ 10), а, во-вторых, при решении обратной задачи, например, методом регуляризации Тихонова эффективно решается недоопределенная система линейных алгебраических уравнений (СЛАУ), вытекающая из выражения (6). Назовем вариант моделирования смазанного изображения на основе (6) приемом усечения искаженного изображения.
Схема для понижения эффекта Гиббса. Для сравнения рассмотрим схему с искусственным введением размытых краев у изображения (рис. 2)
ågj
(i)
=
Δ
1 +1
i+Δ k=i
qj
(k),
i =1, 2,…, n + Δ, j =1, 2,…, m,
(7)
где
qj (k) = îïïïíïìw0,j (k - Δ),
1£ k- Δ £n, èíà÷å.
(8)
Изображение g согласно выражению (7) получается шириной n + Δ с размытыми левым и правым краями. Заметим, что введение размытых краев не следует рассматривать как переход к финитности функции w. Введение размытых краев направлено на понижение возможных искажений типа “звона” (ложных волн, эффекта
Гиббса) [9, 13] при решении обратной задачи. Для их подавления можно воспользоваться m-функцией edgetaper [9] системы MatLab, также размывающей края изображения g, однако схема (7), как показало моделирование, является более эффективной.
Обратная задача. Рассмотрим задачу реконструкции смазанного изображения. Она сводится к решению одномерного интегрального уравнения Фредгольма I рода типа свертки (2) относительно wy(ξ) при каждом фиксированном значении y, играющем роль параметра. В уравнении (2), получаемом в результате преобразования уравнения (1) [3, 4, 6], ФРТ выражается формулой
= íïïîïïì10/,Δ,
h(x - ξ) =
x £ ξ £ x + Δ èëè - Δ £ x - ξ £ 0, èíà÷å.
(9)
Решение уравнения (2) методом ПФ с регуляризацией Тихонова имеет вид [4–8, 17]
¥
òwαy
(ξ)
=
1 2π
Wαy (ω)e-iωξdω,
-¥
где
(10)
Wαy
(ω)
=
|
H(-ω)Gy (ω) H(ω) |2 +αω2 p
,
(11)
¥¥
ò òH(ω) = h(x)eiωxdx, Gy (ω) = gy (x)eiωxdx,
-¥ -¥
a > 0 – параметр регуляризации, p ≥ 0 – порядок регуляризации. Далее метод ПФ с регуляри-
40 “Оптический журнал”, 76, 5, 2009
зацией в дискретном виде реализован в двух
вариантах вариант 1 – на основе использования g (6), в
этом случае матрица реконструированного изображения wα имеет размер m×(n – Δ);
вариант 2 – на основе использования g (7), в
этом случае матрица wα имеет размер m×(n + Δ). Предлагается в конце решения варианта 2 делать
усечение левого края матрицы wα на Δ пикселов, в результате чего ее размер станет равным m×n.
Рассмотрим также решение уравнения (2)
методом квадратур с регуляризацией Тихонова.
Запишем уравнение (2) в виде уравнения обще-
го типа
b
òAwy º h(x,ξ)wy (ξ)dξ = gy (x),
a
c £ x £ d,
(12)
где
h(x,ξ) = ïìïïíîï10/,Δ,
x £ ξ £ x + Δ, èíà÷å.
(13)
Здесь A – интегральный оператор. В дискретном виде, заменяя интеграл в уравнении (12) конечной суммой, получим для схемы с усечением (6) вариант 3
ån
Awj º h i kwj,k = gj,i,
k=1
i =1,…, n - Δ, j =1,…, m,
(14)
где
h
ik
=
îïíïìïï10/,(Δ
+1),
i£k£ èíà÷å.
i
+
Δ,
(15)
Для схемы с размытыми краями (7) получим
вариант 4
ån
Awj º h ikwj,k = gj,i,
k=1
i =1,…, n + Δ, j =1,…, m,
(16)
h ik
= îìïïïíï10/,(Δ
+1),
i-Δ£k èíà÷å.
£
i,
(17)
Рассмотрим для полноты также вариант 5 двойного усечения, когда в схеме на рис. 1 усекается как g, так и w до размера m×(n – Δ)
ån-Δ
Awj º hikwj,k = gj,i,
k=1
i =1,…, n - Δ, j =1,…, m.
(18)
В соотношениях (14), (16), (18) A – (ленточная) матрица при каждом фиксированном номере строки j, причем A – матрица размера (n – Δ)×n в (14), (n + Δ)×n в (16) и (n – Δ)×(n – Δ) в (18). При каждом фиксированном j имеется СЛАУ относительно искомого вектора w или w. При этом СЛАУ (14) является недоопределенной, СЛАУ (16) – переопределенной, а СЛАУ (18) – определенной. Чтобы реконструировать все изображение, нужно при каждом j (т. е. m раз) решить подобную СЛАУ.
Более кратко запишем СЛАУ (14), (16), (18) в виде
Aw = g.
(19)
Решение СЛАУ (19) методом регуляризации Тихонова [4, 8, 11, 13, 17]
wα = (αI+ ATA)–1АTg,
(20)
где I – единичная матрица, AT – транспониро-
ванная матрица. Отличительной особенностью
метода регуляризации Тихонова является то,
что он позволяет решить недоопределенную,
переопределенную и определенную СЛАУ, давая
приближение к точному решению, в качестве
которого выступает нормальное псевдореше-
ние (оно существует и является единственным)
[3, 17].
Прием усечения, изложенный применитель-
но к равномерному прямолинейному горизон-
тальному смазыванию изображения, может
быть распространен на смазывание под углом,
дефокусирование с различными ФРТ. Эти вари-
анты приема усечения изображения будут рас-
смотрены в последующих публикациях авторов
настоящей статьи.
Численные иллюстрации
В рамках системы программирования MatLab7 авторами был разработан ряд m-функций для моделирования прямой задачи смазывания изображения согласно выражениям (6)–(8) и решения обратной задачи реконструкции изображения согласно выражениям (10), (11) (варианты 1, 2) и (12)–(20) (варианты 3–5). В прямой задаче изображения смазывались и зашумлялись аддитивным гауссовым шумом. При решении обратной задачи параметр регуляризации α (11, 20) выбирался двумя способами: путем визуальной оценки реконструированного изображения wα, путем минимизации относительного среднеквадратического отклонения (СКО) изображения wα от исходного точного изображения w– [18]
“Оптический журнал”, 76, 5, 2009
41
σrel =
å åm n ëêé(wα )ji -wji úùû2
j=1i=1
åm ån w2ji
j=1i=1
(21)
(σreg можно вычислить лишь в модельной задаче, когда w– известно).
На рис. 3а дано исходное серое изображение cameraman.tif w размером 256×256 пикселов, а на рис. 3б – смазанное (Δ = 10) и зашумленное гауссовым шумом ||δg||/||g|| = 0,01 = 1% изображение с усечением g 256×246 (между белыми вертикальными линиями), а также изображение с размытыми краями g 256×266 (во всю ширину рис. 3б).
На рис. 4 приведены кривые зависимости σrel (α) при 1%-шуме и при p = 1 (11) для вариантов 1–5. Кривые построены в результате усреднения по 10 реализациям шумов. Видно, что кривые 3 и 4 имеют наиболее глубокие и широкие минимумы.
Для сравнения рассмотрен также вариант 6 – реконструкция изображения методом параметрической фильтрации Винера с использованием граничных условий [2, 9]. В варианте 6 использованы следующие m-функции системы MatLab7: fspecial, imfilter (с опциями circular (periodic) и symmetric (reflective)) и deconvwnr [9]. При этом в методе параметрической фильтрации Винера в качестве константы K [2, с. 392] использовано отношение шум/сигнал по мощности K = ||δg||2/||g||2 (параметр nsr функции deconvwnr). На рис. 4 пунктиром отмечены значения σrel для варианта 6c – варианта 6 с опцией circular и варианта 6s – варианта 6 с опцией symmetric.
(а)
n (б)
n– n
n+ Рис. 3. а – исходное изображение w. б – смазанное и зашумленное изображение g и g.
Значения αopt и σrel, усредненные по 10 реализациям шума (Δ = 10, p = 1)
Вариант
Отн. шум
1
2
3
||δg|| ||g||
4
5 6c 6s σrel σrel
lgαopt σrel(αopt)
lgαopt σrel(αopt)
lgαopt σrel(αopt)
lgαopt σrel(αopt)
lgαopt σrel(αopt)
0 0,001=0,1%
0,01=1% 0,03=3% 0,09=9%
–2,1 0,15 –2,1 0,15 –2,1 0,15 –2,0 0,17 –1,6 0,20
–2,7 –2,6 –2,5 –2,3 –1,8
0,15 0,15 0,15 0,17 0,21
≤–5 0,038
УДК 517.968: 517.983.54: 621.397.3
РЕКОНСТРУКЦИЯ СМАЗАННЫХ И ЗАШУМЛЕННЫХ ИЗОБРАЖЕНИЙ БЕЗ ИСПОЛЬЗОВАНИЯ ГРАНИЧНЫХ УСЛОВИЙ
© 2009 г. В. С. Сизиков, доктор техн. наук; М. В. Римских; Р. К. Мирджамолов Санкт-Петербургский государственный университет информационных технологий, механики и оптики, Санкт-Петербург
E-mail: sizikov2000@mail.ru, romefamily@yandex.ru
Предложен новый подход к проблеме учета интенсивностей, исходящих из точек вне заданных границ изображения. В ряде работ это учитывается путем введения так называемых граничных условий (анти-рефлективных и др.). В данной работе вместо них предлагается прием усечения, не требующий знания интенсивностей за границами изображения. На основе приема усечения изложены прямая и обратная задачи обработки изображения. Для решения обратной задачи реконструкции искаженного изображения использованы метод преобразования Фурье с регуляризацией Тихонова и метод квадратур (также с регуляризацией). Приведены численные иллюстрации, показавшие, что предложенная в работе методика заметно снижает погрешность реконструкции.
Коды OCIS: 100.0100, 100.2000.
Поступила в редакцию 03.12.2008.
Введение
В данной работе рассматривается задача реконструкции искаженных (смазанных или/и зашумленных) изображений [1, 2]. Рассматривается как прямая задача (моделирование искажения или экспериментальное получение реального искаженного изображения), так и обратная задача (реконструкция, восстановление изображения) [2–4].
Основные соотношения и уравнения, описывающие эти задачи в непрерывном виде, известны по работам [2–11]
x+Δ
(1/Δ) ò wy (ξ) dξ = gy (x) + δg,
x
(1)
¥
ò h(x - ξ)wy (ξ)dξ = gy (x) + δg,
-¥
(2)
¥¥
ò ò h(x - ξ,y - η)w(ξ,η)dξdη = g(x,y) + δg. (3)
-¥ -¥
Здесь Δ – смаз, h – функция рассеяния точки (ФРТ), обычно пространственно-инвариантная, w и g – распределение интенсивности по неискаженному и искаженному изображени-
ям соответственно, δg – помеха. В (1) и (2) ось x направлена вдоль смаза, а y играет роль параметра.
Интегральные уравнения (1) и (2) обычно используются в задаче смазывания, а (3) – в задаче дефокусирования, но часто, например, в работах [2, 5] уравнение (3) используется для решения обеих задач. Интегральные уравнения (2), (3) обычно решаются методами инверсной фильтрации, параметрической винеровской фильтрации, регуляризации Тихонова (прямой и итеративный вариант), простой итерации Фридмана, итераций Ландвебера, сопряженных градиентов (CG) и др. При численной реализации методов используются преобразование Фурье (ПФ), квадратуры, итерации и т. д.
Часто используются, например, в работах [9, 11–14] так называемые граничные условия (BCs): нулевые по Дирихле, периодические, копирующие границу, рефлективные (симметричные) относительно границы, анти-рефлективные относительно границы для учета интенсивностей вне заданных границ изображения [12, 13], когда функция w не является финитной, т. е. носитель изображения в действительности шире поля зрения – заданной области w и g. Непрерывным уравнениям (1)–(3) можно поставить в соответствие дискретные выражения (1)–(3) вида
38 “Оптический журнал”, 76, 5, 2009
g = Aw + δg,
(4)
где A – матрица, связанная с ФРТ и с граничными условиями: в случае нулевых условий матрица A есть блок теплицевых матриц, в случае периодических условий – блок циркулянтных матриц, в случае рефлективных относительно границы – сумма блоков теплицевых и ганкелевых матриц [12]. Заметим, что в некоторых случаях искомую функцию w можно считать финитной, например, когда рассматривается космический объект [12, Fig. 4]. Однако гораздо чаще искомая функция w не является таковой, как например, в работе [10, рис. 1], хотя авторы этой работы положили, что w финитна и, по существу, использовали нулевое граничное условие.
Основная цель данной работы – показать, что в случае, когда функция w не является финитной, т. е. w отлична от нуля в более широкой области, чем область поля зрения (FOV), можно адекватно решать как прямую, так и обратную задачу, не используя граничные условия. Данная работа продолжает работы [4, 7, 8, 15, 16].
Прием усечения изображения
О граничных условиях. Рассмотрим прямую задачу моделирования искажения изображения. Пусть нам дано исходное (неискаженное) изображение размером m×n пикселов, описываемое матрицей интенсивностей wm×n. Назовем область Dm×n заданной областью исходного изображения (FOV), а границу B этой области – границей исходного изображения.
В редких случаях (когда, например, исходное изображение является изображением космического объекта) функцию w можно считать финитной, т. е. равной нулю вне области D. В этом случае можно без потери информации моделировать искаженное изображение по формуле типа (4), положив размеры матрицы g, равными размерам матрицы w, и используя при этом граничное условие zero. Аналогично, при решении обратной задачи (реконструкции изображения) можно без потери информации полагать размеры матриц w и g одинаковыми, т. е. считать, что w = 0 вне области D.
Однако авторы данной работы, как и авторы работ [9, 11–14] рассматривают более типичный случай, когда вне области D интенсивности w, вообще говоря, отличны от нуля, но информация об этих интенсивностях отсутствует. В этом случае (при моделировании искажения, а также при решении обратной задачи) обычно используют
граничные условия (BCs), наиболее эффективным из которых считается anti-reflective, и которые являются различными вариантами экстраполяции значений w за границу В. Однако, во-первых, введение BCs является довольно искусственным приемом, а во-вторых, их использование порождает сложный математический аппарат. В данной работе предлагается более простое решение этой проблемы, которое будем называть приемом усечения изображения.
Усечение смазанного изображения. Рассмотрим простейший пример моделирования равномерного прямолинейного горизонтального смазывания изображения (прямая задача). Пусть нам дана матрица исходного изображения wm×n и необходимо смоделировать матрицу смазанного изображения g при смазе Δ и при нулевом угле смаза (горизонтальный смаз). Если формировать смазанное изображение такого же размера, как и исходное изображение, т. е. формировать g шириной n, то соотношение (4) более подробно запишется в виде (при замене интеграла в (1) конечной суммой по пикселам, без помехи)
ågj
(i)
=
Δ
1 +1
i+Δ k=i
wj
(k),
i =1, 2,…,n, j =1, 2,…, m,
(5)
где i (и k) – номер столбца, j – номер строки.
В (5) перед знаком суммы поставлен множитель
1/(Δ + 1) потому, что при Δ = 0 выражение (5)
должно переходить в gj(i) = wj(i). Заметим, что в выражении (1) перед знаком интеграла постав-
лен множитель 1/Δ, чтобы при Δ → 0 соотноше-
ние (1) переходило (без помехи) в wy(x) = gy(x) [4, с. 66].
Из (5) видно, что при i = n – Δ + 1, …, n требу-
ется знание значений wj (1), …, wj(n + Δ), однако значения wj(n + Δ)wj(n + Δ) неизвестны. Чтобы преодолеть этот недостаток информации, в об-
ширном цикле работ [11–14] предложены и ис-
пользуются граничные условия (было бы более
точно назвать их заграничными условиями):
от простейшего нулевого до анти-рефлектив-
ного, в которых делаются различные предпо-
ложения относительно недостающих значений
wj (n + 1), …, wj(n + Δ). Авторы статьи предполагают иное решение
данного вопроса, а именно, вместо соотношений
(5) использовать соотношения (рис. 1)
ågj
(i)
=
Δ
1 +1
i+Δ k=i
wj
(k),
i =1, 2,…, n - Δ, j =1, 2,…, m.
(6)
“Оптический журнал”, 76, 5, 2009
39
размыттие у g размыттие у g
i, k w шириной n
Dmxn (FOV)
j
(
m g шириной n –
B
w шириной n
n смаз
Рис. 1. Схема получения смазанного изображения с усечением.
g шириной n +
Рис. 2. Схема с введением размытых краев у изображения g.
Отличие выражения (6) от (5) состоит в том, что согласно (6) моделируется более узкое смазанное изображение. При этом требуются значения wj (1), …, wj(n), которые известны, и не нужно прибегать к такому приему как граничные условия. Правда, g получается несколько у′же, чем w (на Δ пикселов). Но, во-первых, практически эффект заужения невелик (обычно n ≈ 300÷600, а Δ ∼ 10), а, во-вторых, при решении обратной задачи, например, методом регуляризации Тихонова эффективно решается недоопределенная система линейных алгебраических уравнений (СЛАУ), вытекающая из выражения (6). Назовем вариант моделирования смазанного изображения на основе (6) приемом усечения искаженного изображения.
Схема для понижения эффекта Гиббса. Для сравнения рассмотрим схему с искусственным введением размытых краев у изображения (рис. 2)
ågj
(i)
=
Δ
1 +1
i+Δ k=i
qj
(k),
i =1, 2,…, n + Δ, j =1, 2,…, m,
(7)
где
qj (k) = îïïïíïìw0,j (k - Δ),
1£ k- Δ £n, èíà÷å.
(8)
Изображение g согласно выражению (7) получается шириной n + Δ с размытыми левым и правым краями. Заметим, что введение размытых краев не следует рассматривать как переход к финитности функции w. Введение размытых краев направлено на понижение возможных искажений типа “звона” (ложных волн, эффекта
Гиббса) [9, 13] при решении обратной задачи. Для их подавления можно воспользоваться m-функцией edgetaper [9] системы MatLab, также размывающей края изображения g, однако схема (7), как показало моделирование, является более эффективной.
Обратная задача. Рассмотрим задачу реконструкции смазанного изображения. Она сводится к решению одномерного интегрального уравнения Фредгольма I рода типа свертки (2) относительно wy(ξ) при каждом фиксированном значении y, играющем роль параметра. В уравнении (2), получаемом в результате преобразования уравнения (1) [3, 4, 6], ФРТ выражается формулой
= íïïîïïì10/,Δ,
h(x - ξ) =
x £ ξ £ x + Δ èëè - Δ £ x - ξ £ 0, èíà÷å.
(9)
Решение уравнения (2) методом ПФ с регуляризацией Тихонова имеет вид [4–8, 17]
¥
òwαy
(ξ)
=
1 2π
Wαy (ω)e-iωξdω,
-¥
где
(10)
Wαy
(ω)
=
|
H(-ω)Gy (ω) H(ω) |2 +αω2 p
,
(11)
¥¥
ò òH(ω) = h(x)eiωxdx, Gy (ω) = gy (x)eiωxdx,
-¥ -¥
a > 0 – параметр регуляризации, p ≥ 0 – порядок регуляризации. Далее метод ПФ с регуляри-
40 “Оптический журнал”, 76, 5, 2009
зацией в дискретном виде реализован в двух
вариантах вариант 1 – на основе использования g (6), в
этом случае матрица реконструированного изображения wα имеет размер m×(n – Δ);
вариант 2 – на основе использования g (7), в
этом случае матрица wα имеет размер m×(n + Δ). Предлагается в конце решения варианта 2 делать
усечение левого края матрицы wα на Δ пикселов, в результате чего ее размер станет равным m×n.
Рассмотрим также решение уравнения (2)
методом квадратур с регуляризацией Тихонова.
Запишем уравнение (2) в виде уравнения обще-
го типа
b
òAwy º h(x,ξ)wy (ξ)dξ = gy (x),
a
c £ x £ d,
(12)
где
h(x,ξ) = ïìïïíîï10/,Δ,
x £ ξ £ x + Δ, èíà÷å.
(13)
Здесь A – интегральный оператор. В дискретном виде, заменяя интеграл в уравнении (12) конечной суммой, получим для схемы с усечением (6) вариант 3
ån
Awj º h i kwj,k = gj,i,
k=1
i =1,…, n - Δ, j =1,…, m,
(14)
где
h
ik
=
îïíïìïï10/,(Δ
+1),
i£k£ èíà÷å.
i
+
Δ,
(15)
Для схемы с размытыми краями (7) получим
вариант 4
ån
Awj º h ikwj,k = gj,i,
k=1
i =1,…, n + Δ, j =1,…, m,
(16)
h ik
= îìïïïíï10/,(Δ
+1),
i-Δ£k èíà÷å.
£
i,
(17)
Рассмотрим для полноты также вариант 5 двойного усечения, когда в схеме на рис. 1 усекается как g, так и w до размера m×(n – Δ)
ån-Δ
Awj º hikwj,k = gj,i,
k=1
i =1,…, n - Δ, j =1,…, m.
(18)
В соотношениях (14), (16), (18) A – (ленточная) матрица при каждом фиксированном номере строки j, причем A – матрица размера (n – Δ)×n в (14), (n + Δ)×n в (16) и (n – Δ)×(n – Δ) в (18). При каждом фиксированном j имеется СЛАУ относительно искомого вектора w или w. При этом СЛАУ (14) является недоопределенной, СЛАУ (16) – переопределенной, а СЛАУ (18) – определенной. Чтобы реконструировать все изображение, нужно при каждом j (т. е. m раз) решить подобную СЛАУ.
Более кратко запишем СЛАУ (14), (16), (18) в виде
Aw = g.
(19)
Решение СЛАУ (19) методом регуляризации Тихонова [4, 8, 11, 13, 17]
wα = (αI+ ATA)–1АTg,
(20)
где I – единичная матрица, AT – транспониро-
ванная матрица. Отличительной особенностью
метода регуляризации Тихонова является то,
что он позволяет решить недоопределенную,
переопределенную и определенную СЛАУ, давая
приближение к точному решению, в качестве
которого выступает нормальное псевдореше-
ние (оно существует и является единственным)
[3, 17].
Прием усечения, изложенный применитель-
но к равномерному прямолинейному горизон-
тальному смазыванию изображения, может
быть распространен на смазывание под углом,
дефокусирование с различными ФРТ. Эти вари-
анты приема усечения изображения будут рас-
смотрены в последующих публикациях авторов
настоящей статьи.
Численные иллюстрации
В рамках системы программирования MatLab7 авторами был разработан ряд m-функций для моделирования прямой задачи смазывания изображения согласно выражениям (6)–(8) и решения обратной задачи реконструкции изображения согласно выражениям (10), (11) (варианты 1, 2) и (12)–(20) (варианты 3–5). В прямой задаче изображения смазывались и зашумлялись аддитивным гауссовым шумом. При решении обратной задачи параметр регуляризации α (11, 20) выбирался двумя способами: путем визуальной оценки реконструированного изображения wα, путем минимизации относительного среднеквадратического отклонения (СКО) изображения wα от исходного точного изображения w– [18]
“Оптический журнал”, 76, 5, 2009
41
σrel =
å åm n ëêé(wα )ji -wji úùû2
j=1i=1
åm ån w2ji
j=1i=1
(21)
(σreg можно вычислить лишь в модельной задаче, когда w– известно).
На рис. 3а дано исходное серое изображение cameraman.tif w размером 256×256 пикселов, а на рис. 3б – смазанное (Δ = 10) и зашумленное гауссовым шумом ||δg||/||g|| = 0,01 = 1% изображение с усечением g 256×246 (между белыми вертикальными линиями), а также изображение с размытыми краями g 256×266 (во всю ширину рис. 3б).
На рис. 4 приведены кривые зависимости σrel (α) при 1%-шуме и при p = 1 (11) для вариантов 1–5. Кривые построены в результате усреднения по 10 реализациям шумов. Видно, что кривые 3 и 4 имеют наиболее глубокие и широкие минимумы.
Для сравнения рассмотрен также вариант 6 – реконструкция изображения методом параметрической фильтрации Винера с использованием граничных условий [2, 9]. В варианте 6 использованы следующие m-функции системы MatLab7: fspecial, imfilter (с опциями circular (periodic) и symmetric (reflective)) и deconvwnr [9]. При этом в методе параметрической фильтрации Винера в качестве константы K [2, с. 392] использовано отношение шум/сигнал по мощности K = ||δg||2/||g||2 (параметр nsr функции deconvwnr). На рис. 4 пунктиром отмечены значения σrel для варианта 6c – варианта 6 с опцией circular и варианта 6s – варианта 6 с опцией symmetric.
(а)
n (б)
n– n
n+ Рис. 3. а – исходное изображение w. б – смазанное и зашумленное изображение g и g.
Значения αopt и σrel, усредненные по 10 реализациям шума (Δ = 10, p = 1)
Вариант
Отн. шум
1
2
3
||δg|| ||g||
4
5 6c 6s σrel σrel
lgαopt σrel(αopt)
lgαopt σrel(αopt)
lgαopt σrel(αopt)
lgαopt σrel(αopt)
lgαopt σrel(αopt)
0 0,001=0,1%
0,01=1% 0,03=3% 0,09=9%
–2,1 0,15 –2,1 0,15 –2,1 0,15 –2,0 0,17 –1,6 0,20
–2,7 –2,6 –2,5 –2,3 –1,8
0,15 0,15 0,15 0,17 0,21
≤–5 0,038