Последовательность операций по фильтрации шумов на искаженных изображениях
ИКОНИКА – НАУКА ОБ ИЗОБРАЖЕНИИ
УДК 621.397.3
ПОСЛЕДОВАТЕЛЬНОСТЬ ОПЕРАЦИЙ ПРИ ФИЛЬТРАЦИИ ШУМОВ НА ИСКАЖЕННЫХ ИЗОБРАЖЕНИЯХ
© 2013 г.
В. С. Сизиков*, доктор техн. наук; Р. А. Экземпляров**, аспирант ** Санкт-Петербургский национальный исследовательский университет информационных ** технологий, механики и оптики, Санкт-Петербург
** Санкт-Петербургский государственный политехнический университет, Санкт-Петербург
** Е-mail: sizikov2000@mail.ru, rexe@yandex.ru
Рассмотрены вопросы, связанные с зашумлением смазанных или/и дефокусированных изображений. Проанализирована последовательность выполнения фильтрации шумов на таких изображениях – до устранения смазывания/дефокусирования или после него. Введены понятия предшествующей и последующей фильтрации шумов. Устранение смазывания/дефокусирования ряда изображений выполнялось методами параметрической фильтрации Винера и регуляризации Тихонова, а фильтрация шумов – методами адаптивной фильтрации Винера и медианной фильтрации.
Ключевые слова: зашумленное смазанное или дефокусированное изображение, предшествующая и последующая фильтрация шумов, устранение смазывания и дефокусирования.
Коды OCIS: 030.4280, 100.0100
Поступила в редакцию 26.06.2012
Постановка задачи
Искаженные (смазанные, дефокусированные) изображения обычно зашумлены [1–6]. Цель данной работы – анализ следующих вопросов: способны ли такие устойчивые методы, как метод параметрической фильтрации Винера, метод регуляризации Тихонова и др. [1–3, 5–10] не только устранять (в той или иной мере) смазывание или дефокусирование изображения, но и фильтровать шумы; в какой степени способствуют фильтрации шумов такие специальные методы, как метод адаптивной фильтрации Винера, метод медианной фильтрации и др. [1–6]; какова эффективность перечисленных методов в зависимости от типа шума (гауссов, импульсный и др. [1–6]); какова наиболее эффективная последовательность (очередность) фильтрации шума (до или после устранения смаза/дефокусирования) [11] – основной вопрос и др. Данные вопросы рассмотрены путем проведения экспериментальной проверки с получением количественных оценок качества восстановленных изображений при различных типах шумов и очередности их фильтрации на репрезентативных выбор-
ках. Первые три вопроса в той или иной мере рассмотрены в публикациях [1–13] и при рассмотрении 4-го (основного) вопроса результаты этих публикаций будут учтены.
Качественный анализ
Сначала обсудим вопрос о фильтрации шумов качественно. Если (равномерное и прямолинейное) смазывание изображения объекта произошло из-за сдвига камеры (например, сдвига фотоаппарата или рассогласования вращения телескопа и небесной сферы), а между объектом и камерой имелся импульсный шум (пыль, капли влаги и т. п.), то будут смазанными и объект, и шум (все точки объекта и шумовые точки превратятся в штрихи, рис. 1б). В этом случае логично сначала устранить смазывание объекта и шума некоторым методом (параметрической фильтрации Винера, регуляризации Тихонова и др.). Смазанный импульсный шум должен превратиться практически в точки. Затем нужно выполнить последующую фильтрацию шума (адаптивной фильтрацией Винера, медианной фильтрацией и др.). Заметим, что в данном случае по шумовым
“Оптический журнал”, 80, 1, 2013
39
(а) (б)
(в) (г)
Рис. 1. а – точное изображение w–, б – изображение, зашумленное импульсным шумом и смазанное с шумом, с размытием краев, в – устранение смазывания методом параметрической фильтрации Винера (K = 10–4) и методом квадратур с регуляризацией Тихонова (a = 10–6), г – последующая медианная фильтрация изображения (в) с маской 3×3, srel = 0,105.
томобиля [12] и др.), а в среде между объектом и неподвижной камерой имелся шум (туман, пыль, капли, снежинки и др.). В результате на изображении объект будет смазанным, а шум – неподвижным. В этом случае нужно также в ыполнить сначала предшествующую фильтрацию шума, а затем устранить смазывание.
Наконец, если изображение объекта несмазано и недефокусировано, а импульсный шум смазан (дождь, падающие снежинки и др.), то целесообразно выполнить лишь фильтрацию шума, например, медианным фильтром.
Могут существовать и другие ситуации с зашумлением, смазыванием или/и дефокусированием [1–3, 5, 6], причем шум может быть различных типов (импульсный, гауссов, рав номерный, мультипликативный и др.), изображение (а также шум) – черно-белым и цветным, а искажение – более сложным, чем смазывание или дефокусирование [2, 13, 14].
(а) (б)
штрихам можно оценить смаз D (и его направление), т. е. шум способствует определению функции рассеяния точки (ФРТ). Если же происходил не сдвиг камеры, а неправильная установка фокуса, то будут дефокусированными и объект, и шум (все точки превратятся в размытые пятна в виде ФРТ, рис. 2б). В этом случае нужно сначала устранить дефокусирование объекта и шума, а затем выполнить последующую фильтрацию шума.
Если изображение объекта смазалось за счет сдвига камеры (или дефокусировалось), а импульсный шум возник в результате сбоя в работе сенсоров матрицы ПЗС, то изображение объекта будет смазанным (или дефокусированным), а импульсный шум представлен в виде так называемых “битых пикселов” (рис. 3а и 4а). В этом случае естественно сначала выполнить предшествующую фильтрацию шума, а затем устранить смазывание (или дефокусирование) изображения объекта. Аналогичную последовательность операций следует соблюсти и в случае, когда изображение смазалось за счет движения объекта (самолета [5, 6, 12], ав-
40
(в) (г)
Рис. 2. а – точное изображение Луны w–, б – изображение, зашумленное импульсным шумом и дефокусированное вместе с шумом, в – изображение (б), рефокусированное методом двумерного ПФ с регуляризацией Тихонова, a = 10–8 (изображение и импульсный шум восстановились), г – последующая медианная фильтрация изображения (в) с маской 3×3, srel = 0,0259.
“Оптический журнал”, 80, 1, 2013
(а) (б)
(а) (б)
(в) (г)
(в) (г)
(д) (е)
Рис. 3. а – смазанное изображение Сатурна с четырьмя спутниками и с наложенным несмазанным импульсным шумом, б – обработка изображения (а) методом параметрической фильтрации Винера, K = 10–1, srel = 0,192, в – предшествующая медианная фильтрация изображения (а) с маской 3×3, г – параметрическая фильтрация Винера изображения (в) (K = 10–3, srel = 0,050).
Заметим, что последовательность фильтрации шума в вышеописанных ситуациях может казаться естественной и очевидной, однако приведенные ниже результаты численных экспериментов говорят о том, что отмеченный порядок не всегда является оптимальным.
Основные математические соотношения
Как прямая задача моделирования искажений изображений, так и обратная задача восстановления изображений описывается соотношениями [5–10]
¥
ò h(x - ξ)wy (ξ)dξ = gy (x) + dg,
-¥
¥¥
ò ò h(x -ξ, y -h)w(ξ, h)dξdh = g(x, y) + dg,
-¥ -¥
(1) (2)
где h – ФРТ, аппаратная функция, обычно пространственно-инвариантная (разностная),
Рис. 4. а – смазанное изображение и несма-
занный импульсный шум, б – изображение,
аналогичное (а), но с размытием краев, в –
обработка изображения (а) методом пара-
метрической фильтрации Винера, K = 10–1,
srel = 0,253, г – предшествующая медианная фильтрация изображения (а) с маской
3×3, srel = 0,243, д – параметрическая трация Винера изображения (г), K =
филь10–2,5,
srel = 0,139, e – восстановление методом квадратур с регуляризацией Тихонова изображе-
ния, близкого (г), a = 10–3, srel = 0,133.
w и g – распределения интенсивности по неискаженному и искаженному изображениям соответственно, dg – помеха (шум). В выражении (1) ось x направлена вдоль смаза, а y – перпендикулярно смазу (играет роль параметра).
Если рассматривается обратная задача, то соотношения (1) и (2) являются интегральными уравнениями (ИУ) Фредгольма I рода. Их решение есть некорректная задача [5, 6, 8, 15, 16], поэтому будем использовать для их решения устойчивые методы – методы регуляризации Тихонова и параметрической фильтрации
“Оптический журнал”, 80, 1, 2013
41
Винера [2, 3, 5–10, 15, 16], а решения другими методами (Люси–Ричардсона, “слепой” деконволюции, итераций и проч.), полученные в других работах [1–3, 5, 6, 13], также будем учитывать.
Решение ИУ типа свертки (1) методом одномерного преобразования Фурье (ПФ) с регуляризацией Тихонова имеет вид [5–10, 12, 15–17]
òway (ξ) =
¥
1 2p
-¥
H(-w)Gy (w) H(w) 2 + aw2 p
exp(-iwξ)dw,
где
¥
H(w) = ò h(x)exp(iwx)dx,
-¥
¥
òGy (w) = gy (x)exp(iwx)dx
-¥
(3) (4)
– спектры ПФ, a > 0 – параметр регуляризации, p ≥ 0 – порядок регуляризации. Для решения ИУ (1) согласно соотношениям (3), (4) авторами разработана в рамках системы программирования MatLab7 m-функция desmearingf.m [5, 6].
ИУ (1) может быть решено также методом квадратур с регуляризацией Тихонова [5, 6, 8, 9, 16] (исследования показали [5, 6, 9], что такой метод дает обычно более точное восстановление изображений, чем метод ПФ с регуляризацией). Согласно этому методу, ИУ (1) при каждом значении y приводится к системе линейных алгебраических уравнений (СЛАУ) вида Awy = gy, где A – матрица, wy и gy – векторы, а решение СЛАУ методом регуляризации Тихонова выглядит как (при каждом y)
way = (aI + AT A)-1 AT gy,
(5)
где I – единичная матрица, а AT – транспони рованная матрица (m-функция desmearingq.m). Решения (3)–(5) были использованы для задачи восстановления смазанных изображений (так называемый вариант 2 решения (3) и вариант 4 решения (5) [5, 6, 9]).
Что же касается задачи восстановления дефокусированных изображений, то она описывается двумерным ИУ (2) и ее решение методом двумерного ПФ с регуляризацией Тихонова представляется в виде (ср. (3), (4)) [5–8, 10, 15–17]
ò òwa (ξ, h) =
1 4p2
¥¥
Wa (w1, w2)´
-¥ -¥
´exp[-i(w1ξ + w2h)]dw1dw2,
(6)
где
Wa (w1, w2) =
H*(w1, H(w1, w2)
w2)G(w1, w2) 2 + a(w12 + w22
)
p
,
(7)
H(w1, w2) =
¥¥
= ò ò h(x, y)expi(w1x + w2y)dxdy, -¥ -¥
(8)
G(w1, w2) =
¥¥
= ò ò g(x, y)expi(w1x + w2y)dxdy -¥ -¥
(9)
(m-функция refocusingT.m), а решение методом параметрической фильтрации Винера принимает вид (ср. (6)–(9)) [2, 3, 5–8, 10]
ò ò=
1 4p2
¥¥ -¥ -¥
WK
(w1,
wK (ξ, h) =
w2 ) exp[-i(w1ξ
+
(10)
w2h)]dw1dw2,
где
WK
(w1,
w2 )
=
H*(w1, w2)G(w1, w2) H(w1, w2) 2 + K
,
(11)
K ≥ 0 – некоторая константа (параметр, дающий оценку отношению шум/сигнал по мощности). Метод параметрической фильтрации Винера реализован в системе MatLab7 в виде m-функции deconvwnr.m [3].
Заметим, что метод параметрической фильтрации Винера и m-функция deconvwnr.m часто используются [2, 3] не только для восстановления дефокусированных, но и смазанных изображений, причем в обоих случаях в двумерном виде (8)–(11).
В данной работе рассмотрены два типа шумов – гауссов и биполярный импульсный (типа “salt & pepper”), а также специальные фильтры для фильтрации шумов – среднеарифметический фильтр, метод адаптивной фильтрации Винера и медианный фильтр [1–4]. Что же касается других типов шумов (равномерный, мультипликативный) и методов их фильтрации (адаптивный медианный фильтр, ранговый фильтр и др.), то результаты их использования в различных публикациях [1–4] также учтены в данной работе.
В заключение данного пункта остановимся на вопросе оценки обработки изображений. К ачественно оценивать искаженные и восстановленные изображения будем визуально. А количественно степень искажения и, особенно, погрешность восстановления изображения
42 “Оптический журнал”, 80, 1, 2013
будем характеризовать относительным среднеквадратическим отклонением – СКО изображения w~ от точного изображения w– [5, 6, 9, 10]
srel =
w -w L2 || w ||L2
=
å å( )m n wji -wji 2
j=1 i=1
.
å åm n wj2i
j=1 i=1
(12)
Под w~ будем подразумевать искаженное
(смазанное, дефокусированное, зашумленное) изображение w~, а также изображение, восста-
новленное методом параметрической фильтра-
ции Винера wK или методом регуляризации Тихонова wa. Для вычисления srel согласно формуле (12) авторами разработана m-функция
rmsd.m [5, 6]. Будем обозначать через Kopt и aopt значения K и a, дающие минимум srel.
Вместо формулы (12) можно использовать
выражения
srel = || w -w ||L1 || w ||L1 , s = || w -w ||L2 , s = || w -w ||L1 ,
( )PSNR =10lg wmaxmn || w -w ||2L2
и др. [5, 6]. Однако в данной работе будем поль-
зоваться выражением (12) – типа (безразмерного) отношения шум/сигнал.1
1 1. Приведенные пользовать лишь
ввыслруачжаеен, икяогддлаяиsзrоeбl риаsжемноиженw–о
исиз-
вестно, т. е. при обработке модельных изображений
(фантомов), что и было сделано в данной работе.
Специальный же, весьма сложный вопрос об оценке
погрешности восстановления изображений без знания w– в результате решения некорректной задачи
[12, 15–17] в данной работе не рассматривается.
2. Можно использовать также некоторые способы выбора параметра регуляризации a и параметра K (без знания w–): способ невязки [8, 15–17] (нужно знать погрешность ||dg|| измерения искаженного изображения g), обобщенный принцип невязки [16] (нужно знать ||dg|| и погрешность знания ФРТ ||dh||), способ псевдо обратного оператора [16] и др. Однако эти способы часто дают завышенные по сравнению с aopt и Kopt значения a и K, а значит, заглаженные решения wa и wK. Кроме того, они требуют знания погрешностей ||dg|| и || dh||, которые обычно известны неточно. Поэтому использование способов выбора a и K отложим на дальнейшее.
3. В данной работе будем считать, что параметры искажений изображений (величина смаза D, угол смаза q, уровень и тип шума и проч.) известны точно. Вопросы, связанные с неточностями параметров искажений, специально рассмотрены в различных публикациях ([2, 5, 6, 13, 17, 18] и др.). Мы же, полагая, что w– и параметры искажений известны точно, будем тем самым определять потенциальные возможности рассматриваемой в данной работе методики.
Численные результаты
Ниже приведены результаты обработки изо-
бражений нескольких объектов. Для каждо-
го объекта приведено точное (неискаженное)
изображение. Это сделано специально, чтобы
оценить количественно качество обработки
(посредством соотношения (12)). Затем вноси-
лись искажения, выполненные искусственно
(программным путем). Что касается обработ-
ки изображений, когда истинное изображение
неизвестно, а искажения имеют естественную
природу, то такая обработка была также вы-
полнена [5, 6] и в заключительных выводах
у чтена.
Смазанное изображение и несмазанный
шум. На рис. 4а представлено черно-белое изо-
бражение girl.jpg 256×256, смазанное (D = 20,
q = 27°) и зашумленное несмазанным импульс-
ным шумом (доля зашумленных пикселов
d = 0,02), а на рис. 4б – аналогичное изобра-
жение, но с размытием краев (для уменьшения
эффекта Гиббса [5, 6, 9, 10]). Далее выполня-
лось устранение смазывания методами параме-
трической фильтрации Винера при различных
K и квадратур с регуляризацией Тихонова при
различных a, а также фильтрация шума меди-
анным фильтром и адаптивным фильтром Ви-
нера с различными масками, причем фильтра-
ция шума выполнялась как предшествующая
(pre), так и последующая (post).
На рис. 5 представлены полученные зави-
симости относительных погрешностей восста-
новления изображения srel(K) (под w~ в форму-
ле (12) подразумевается нию (10)) и srel(a) (под w~
wK согласно выражев формуле (12) подраз-
умевается wa согласно выражению (5)) для слу-
чая, когда фильтрация импульсного шума вы-
полнялась медианным фильтром (с помощью
m-функции medfilt2.m [3, 4]). Точное изображение w– представлено на рис. 1а. По миниму-
мам кривых srel(K) и srel(a) определялись Kopt и aopt, а также srel(Kopt) и srel(aopt) – значения, при которых изображение восстанавливается
с наименьшей погрешностью согласно коли-
чественному критерию (по визуальному кри-
терию “наилучшему” восстановлению могут
соответствовать несколько другие K и a).
Из приведенного рисунка видно, что значе-
ния погрешностей srel(Kopt) и srel(aopt) заметно зависят от последовательности выполнения
фильтрации шума (предшествующая или по-
следующая фильтрация), а также от размера
маски медианного фильтра. В данном примере
“Оптический журнал”, 80, 1, 2013
43
srel(K), srel(a)
1
3Wpost 0,8
1W 3Tpost
1T
9preW 9Wpost 0,6
0,4 3preW 9Tpost
0,2 3preT9preT
0 –6 –5 –4 –3 –2 –1 0
lgK, lga
Рис. 5. Погрешность восстановления изображения на рис. 4а и 4б: srel(K) – методом параметрической фильтрации Винера (W) и srel(a) – методом квадратур с регуляризацией Тихонова (T) с использованием предшествующей (pre) и последующей (post) фильтрации импульсного шума медианным фильтром с масками 1×1, 3×3 и 9×9.
предшествующая фильтрация импульсного шума при умеренной маске 3×3 ведет к погрешности srel, примерно в 1,5 раза меньшей, чем последующая фильтрация, что согласуется с вышеприведенным качественным анализом. Представленные на рис. 5 зависимости показывают, что метод регуляризации Тихонова дает чуть меньшую погрешность srel, чем метод параметрической фильтрации Винера.
Заметим, что при маске 1×1 медианная фильтрация посуществу не выполняется, поэтому кривая 1W на рис. 5 аналогична кри-
вым 1preW и 1Wpost, а 1T аналогична 1preT и 1Tpost.
На рис. 4в приведен результат параметрической фильтрации Винера после медианной фильтрации шума с маской 1×1. Он показывает, что параметрическая фильтрация Винера не только устраняет смазывание изображения, но и подавляет шум (K = 10–1, srel = 0,253), однако делает это лишь частично (метод квадратур с регуляризацией Тихонова дает аналогичный результат). Значительно более точный результат получается, когда выполняется предшествующая медианная фильтрация с бо′льшим размером маски (рис. 4г при маске 3×3), а затем устраняется смазывание (рис. 4д, е).
Данный пример был промоделирован также при использовании адаптивного фильтра Винера (с помощью m-функции wiener2.m [3, 4]) применительно к фильтрации как импульсного, так и гауссова шумов. В табл. 1 приведены результаты для импульсного шума, а в табл. 2 – для гауссова шума.
Данные, приведенные в табл. 1, показывают (в дополнение к рис. 4 и 5), что если шум импульсный, то медианный фильтр фильтрует его заметно лучше, чем адаптивный винеровский (а также среднеарифметический) фильтр: погрешность s примерно в 1,8–1,9 раза меньше в случае предшествующей фильтрации и умеренной маски 3×3, в 1,5 раза меньше в случае маски 5×5 (результат не приведен в табл. 1) и в 1,2 раза меньше в случае маски 9×9.
В табл. 2 приведены аналогичные результаты для случая гауссова 1%-го шума. Они показывают, что погрешность фильтрации s гауссова шума адаптивным фильтром Винера чуть меньше (примерно в 1,1 раза в данном
Таблица 1. Значения K = Kopt, s = srel(Kopt), a = aopt и s = srel(aopt)
Маска 1×1
Маска 3×3
Маска 9×9
Смазанное изображение girl.jpg
(D = 20, q = 27°), несмазанный Метод пар. Метод кв. Метод пар. Метод кв. Метод пар. Метод кв.
импульсный шум (d = 0,02) фильтрации с регуляр. фильтрации с регуляр. фильтрации с регуляр.
Винера Тихонова Винера Тихонова Винера Тихонова
Медианный предшествующая K = 10–1
фильтр
фильтрация
s = 0,253
a = 10–1 s = 0,252
K = 10–2,5 s = 0,139
a = 10–3 s = 0,133
K = 10–2 s = 0,208
a = 10–2,5 s = 0,187
последующая фильтрация
K = 10–1 s = 0,253
a = 10–1 s = 0,252
K = 10–1,5 s = 0,207
a = 10–1,5 s = 0,214
K = 10–2 s = 0,201
a = 10–2 s = 0,204
Адаптивный предшествующая
фильтр
фильтрация
Винера
последующая
фильтрация
K = 10–1 s = 0,253
K = 10–1 s = 0,253
a = 10–1 s = 0,253
a = 10–1 s = 0,253
K = 10–1 s = 0,250
K = 10–1 s = 0,237
a = 10–1 s = 0,251
a = 10–1 s = 0,244
K = 10–1,5 s = 0,236
K = 10–1,5 s = 0,221
a = 10–1,5 s = 0,236
a = 10–1,5 s = 0,227
44 “Оптический журнал”, 80, 1, 2013
Таблица 2. Значения K = Kopt, s = srel(Kopt), a = aopt и s = srel(aopt)
Маска 1×1
Маска 3×3
Маска 9×9
Смазанное изображение girl.jpg
(D = 20, q = 27°), несмазанный Метод пар. Метод кв. Метод пар. Метод кв. Метод пар. Метод кв.
1%-ый гауссов шум
фильтрации с регуляр. фильтрации с регуляр. фильтрации с регуляр.
Винера Тихонова Винера Тихонова Винера Тихонова
Медианный предшествующая
фильтр
фильтрация
K = 10–2,5 s = 0,126
a = 10–3 s = 0,117
K = 10–2,5 s = 0,140
a = 10–3 s = 0,129
K = 10–2 s = 0,207
a = 10–2,5 s = 0,185
последующая фильтрация
K = 10–1 s = 0,126
a = 10–3 s = 0,117
K = 10–3 s = 0,120
a = 10–3,6 s = 0,114
K = 10–3,2 s = 0,175
a = 10–3,4 s = 0,167
Адаптив- предшествующая
ный фильтр фильтрация
Винера
последующая
фильтрация
K = 10–1 s = 0,126
K = 10–1 s = 0,126
a = 10–3 s = 0,117
a = 10–3 s = 0,117
K = 10–3 s = 0,129
K = 10–3,5 s = 0,109
a = 10–3,4 s = 0,117
a = 10–3,7 s = 0,105
K = 10–2,3 s = 0,194
K = 10–3,6 s = 0,145
a = 10–2,8 s = 0,177
a = 10–3,8 s = 0,139
прим ере), чем медианным фильтром. Это подтверждается и расчетами в случае гауссова 5%-го шума.
Однако табл. 1 и 2 демонстрируют и неожиданный результат: если шум несмазанный импульсный и используется адаптивный фильтр Винера (см. табл. 1), то последующая фильтрация шума дает чуть меньшую погрешность s, чем предшествующая фильтрация (примерно в 1,05 раз для маски 3×3 в данном примере). Еще значительнее этот эффект в случае несмазанного гауссова шума и использования как медианного фильтра, так и адаптивного фильтра Винера (примерно в 1,1–1,15 раз в случае маски 3×3, табл. 2).
Смазанное изображение и смазанный шум. На рис. 1 представлены результаты обработки изображения girl.jpg в случае, когда изображение зашумлено импульсным шумом и смазано вместе с шумом (рис. 1б). Для устранения смазывания изображения использованы методы параметрической фильтрации Винера, а также квадратур с регуляризацией Тихонова и при каждом значении параметров K и a выполнена последующая медианная фильтрация с масками 1×1, 3×3, 9×9 и др. В результате построены кривые типа рис. 5 и получены следующие значения K = Kopt, s = srel(Kopt), a = aopt и s = srel(aopt): при маске 1×1 – K = 10–2,6, s = 0,148, a = 10–3, s = 0,141; при маске 3×3 – K = 10–4, s = 0,106, a = 10–6, s = 0,105; при маске 9×9 – K = 10–3,7, s = = 0,165, a = 10–6, s = 0,165. На рис. 1в приведено и зображение, полученное после устранения смазывания методами параметрической фильтрации Винера и квадратур с регуляриза-
цией Тихонова (результаты для обоих методов практически одинаковы), а на рис. 1г – изображение, полученное в результате последующей медианной фильтрации шума с маской 3×3. Рис. 1в демонстрирует хорошее восстановление изображения объекта и шума, а рис. 1г – окончательное восстановление объекта практически без смазывания и шума.
Дефокусированное изображение и дефокусированный шум. На рис. 2а представлено точное изображение Луны w– (файл moon. tif 455×325), а на рис. 2б – ее дефокусированное изображение (ФРТ задана в виде однородного диска радиуса r = 10 пкс), зашумленное импульсным шумом с долей зашумления d = 0,02 (шум также дефокусирован и поэтому на рис. 2б невидим). Было устранено дефокусирование методом двумерного ПФ с регуляризацией Тихонова согласно выражениям (6)–(9), а также методом параметрической фильтрации Винера согласно соотношениям (8)–(11) для ряда значений параметров a и K и при каждом a и K выполнена последующая фильтрация шума медианным фильтром с масками 1×1, 3×3, 9×9 и др.
На рис. 6 представлены полученные зависимости srel(K) и srel(a) – относительные погрешности восстановления изображения Луны. По минимумам этих зависимостей были определены оптимальные значения K = = Kopt, s = srel(Kopt), a = aopt и s = srel(aopt): при маске 1×1 – K = 10–2,9, s = 0,0794, a = 10–0,5, s = 0,0977; при маске 3×3 – K = = 10–3,5, s = 0,0691, a = 10–8, s = 0,0259; при маске 9×9 – K = 10–4, s = 0,0640, a = 10–8, s = 0,0549.
“Оптический журнал”, 80, 1, 2013
45
srel(K), srel(a)
0,5
0,4
0,3
3Wpost
1W
0,2 1T
9Wpost
0,1 9Tpost 3Tpost
0 –8 –7 –6 –5 –4 –3 –2 –1 0
lgK, lga
Рис. 6. Погрешности восстановления изображения Луны (см. рис. 2). Обозначения такие же, как на рис. 5.
На рис. 2в представлено изображение Луны, рефокусированное методом двумерного ПФ с регуляризацией Тихонова при a = 10–8, s = 0,2156 (видно, что изображение Луны и импульсный шум восстановились), а на рис. 2г – изображение, полученное в результате последующей фильтрации импульсного шума медианным фильтром с маской 3×3 (s = 0,0259). Восстановление изображения методом параметрической фильтрации Винера заметно хуже (погрешность s = 0,0691, т. е. в 2,7 раз больше, чем методом Тихонова).
Смазанное изображение с мелкими деталями и несмазанный шум. Рассмотрим изображение, представленное на рис. 3а, аналогичное изображению на рис. 4а, поскольку оба изображения – это смазанные изображения с наложенным несмазанным импульсным шумом. Однако изображение на рис. 3а отличается наличием мелких деталей, а именно, четырех маленьких спутников Сатурна, слабо видимых на фоне шума. Изображение saturn.png 1500×1200×3 – цветное, преобразованное в серое 1500×1200, D = 40 пкс, q = 35°. На изображение наложен несмазанный импульсный шум с долей зашумления d = 0,02. Смазывание изображения Сатурна со спутниками произошло из-за рассогласования вращения телескопа и небесной сферы (за время довольно продолжительной экспозиции). А импульсный шум, как и на рис. 4а, возник в результате выхода из строя ряда сенсоров матрицы ПЗС или/и из-за электромагнитных помех. Другими словами, импульсный шум
в данном примере – это аппаратурный шум, но не шум в среде.
Была выполнена фильтрация шума медианным фильтром (а также адаптивным фильтром Винера) и устранено смазывание методами параметрической фильтрации Винера при различных K, а также квадратур с регуляризацией Тихонова при различных a, причем фильтрация шума выполнялась как предшествующая, так и последующая. Результаты получились качественно весьма похожими на рис. 5 и табл. 1, однако количественно относительные погрешности s заметно меньше.
На рис. 3б приведены результаты параметрической фильтрации Винера и регуляризации Тихонова без выполнения специальной (медианной или адаптивной винеровской) фильтрации шума при K = Kopt = a = aopt = = 10–1 и s = srel(Kopt) = srel(aopt) = 0,192. Видим, как и на рис. 4в, что параметрическая фильтрация Винера и регуляризация Тихонова устранили смазывание и подавили шум, но лишь частично. Более качественный результат получается при добавлении предшествующей медианной фильтрации (с маской 3×3 или большей). На рис. 3в приведен результат предшествующей медианной фильтрации изображения, представленного на рис. 3а, с маской 3×3 (видим, что импульсный шум отфильтровался, а спутники в виде штрихов сохранились). На рис. 3г – результат параметрической фильтрации Винера и регуляризации Тихонова по отношению к изображению на рис. 3в. Погрешности s получились значительно меньшими, чем на рис. 4д, е, 5 и в табл. 1, а именно, K = a = 10–3, s = 0,050, т. е. примерно в 2,7 раз меньше, чем на рис. 5 и в табл. 1. Рис. 3г показывает, что удалось отфильтровать шум и устранить смазывание не только крупного объекта – Сатурна, но и мелких объектов – его спутников.
Для повышения репрезентативности исследования были обработаны также другие изображения: черно-белые cameraman.tif (портретное изображение), textp.bmp, textn.bmp (тексты), plane.bmp (быстролетящая цель – самолет), а также цветные flower.bmp, pears. png, peppers.png, onion.png, football.jpg (предметы), auto.jpg (подвижный объект), pair.jpg, rehoitus.jpg, причем как модельные (фантомы), так и реальные (натурные) изображения [5–12], в том числе, томографические [19].
На основе этих результатов, а также результатов других публикаций ([1–4, 13, 14, 17–19]
46 “Оптический журнал”, 80, 1, 2013
и др.) можно сделать следующие выводы в отношение задачи фильтрации шумов на искаженных (смазанных или дефокусированных) изображениях.
Выводы
1. Если устранять смазывание/дефокусирование, а также фильтровать шум каким-то одним методом (параметрической фильтрации Винера, регуляризации Тихонова, максимального правдоподобия Люси–Ричардсона, “слепой” деконволюции и др. [1–3, 5, 6]), то этого может оказаться недостаточно для получения удовлетворительного результата (см. рис. 3б и 4в). Более точная обработка может получиться, если указанные методы использовать в основном для устранения смазывания/ дефокусирования, а для фильтрации шума дополнительно использовать такие специальные методы, как медианный фильтр, адаптивный медианный фильтр, среднеарифметический фильтр, ранговая фильтрация, адаптивный фильтр Винера и др. [1–4] (см. рис. 1г, 2г, 3г, 4д, е).
2. Погрешность обработки s зависит от последовательности выполнения фильтрации шума на смазанном/дефокусированном изображении – до устранения смазывания/дефокусирования (предшествующая фильтрация) или после него (последующая фильтрация). Например, если на смазанное изображение наложен несмазанный (аппаратурный) импульсный шум (рис. 3а, 4а), то нужно сначала выполнить предшествующую фильтрацию шума медианным фильтром (рис. 3в, 4г), а потом устранить смазывание (рис. 3г, 4д, е). Если же эти операции поменять местами (смазывание –
последующая фильтрация), то это может привести к заметному повышению погрешности s (см. табл. 1, медианный фильтр, маска 3×3).
3. Если импульсный несмазанный шум фильтровать адаптивным фильтром Винера, то последующая фильтрация шума дает чуть меньшую погрешность s, чем предшествующая фильтрация (табл. 1). Такой же результат получается для гауссова несмазанного шума (для любого фильтра, табл. 2). Это – несколько неожиданный результат, не согласующийся с вышеприведенным качественным анализом. Он связан, во-первых, с различием свойств импульсного и гауссова шумов: биполярный импульсный шум типа “salt & pepper” имеет интенсивности лишь 0 и 255 (импульсы), присутствует лишь в d-ой доле (d ∈ [0, 1]) точек изображения и замещает интенсивности на изображении, а гауссов шум имеет различные интенсивности, присутствует в каждой точке изображения и добавляется к нему. Во-вторых, медианная и винеровская адаптивная фильтрации имеют также различия: в медианной практически исключаются импульсы, а в винеровской адаптивной импульсы с большими значениями интенсивности сильно влияют на результат фильтрации.
4. Импульсный шум заметно лучше фильтруется медианным фильтром, чем адаптивным Винера или среднеарифметическим, а гауссов шум, наоборот, лучше фильтруется адаптивным фильтром Винера, чем медианным фильтром (см. табл. 1).
Эти выводы подтвердились обработкой иных изображений: cameraman.tif, textp.bmp, textn.bmp, plane.bmp, flower.bmp, pears.png, peppers.png, onion.png, football.jpg, auto.jpg, pair.jpg, rehoitus.jpg и др.
* * * * *
ЛИТЕРАТУРА
1. Lim J.S. Two-dimensional signal and image processing. New Jersey: Prentice Hall PTR, 1990. 694 p.
2. Гонсалес Р., Вудс Р. Цифровая обработка изображений. М.: Техносфера, 2006. 1072 с.
3. Гонсалес Р., Вудс Р., Эддинс С. Цифровая обработка изображений в среде MATLAB. М.: Техносфера, 2006. 616 с.
4. Дьяконов В., Абраменкова И. MATLAB. Обработка сигналов и изображений. СПб.: Питер, 2002. 608 с.
5. Сизиков В.С. Обратные прикладные задачи и MatLab. СПб.: Лань, 2011. 256 с.
6. Сизиков В.С. Интегральные уравнения и MatLab в задачах томографии, иконики и спектроскопии. Saarbrücken: LAP, 2011. 252 с.
“Оптический журнал”, 80, 1, 2013
47
7. Сизиков В.С., Белов И.А. Реконструкция смазанных и дефокусированных изображений методом регуляризации // Оптический журнал. 2000. Т. 67. № 4. С. 60–63.
8. Сизиков В.С. Математические методы обработки результатов измерений. СПб.: Политехника, 2001. 240 с. 9. Сизиков В.С., Римских М.В., Мирджамолов Р.К. Реконструкция смазанных и зашумленных изображений
без использования граничных условий // Оптический журнал. 2009. Т. 76. № 5. С. 38–46. 10. Сизиков В.С. Прием “усечение–размытие–поворот” для восстановления искаженных изображений // Оп-
тический журнал. 2011. Т. 78. № 5. С. 18–26. 11. Экземпляров Р.А., Сизиков В.С. Устранение смазывания и дефокусирования изображений с предваритель-
ной фильтрацией шумов в рамках системы MatLab // Мат. международ. конф. “XXXIX Неделя науки СПбГПУ”, ч. XIII (ИМОП) / Под ред. Рудского А.И. СПб.: СПбГПУ, 2010. С. 223–225. 12. Petrov Yu.P., Sizikov V.S. Well-Posed, Ill-Posed, and intermediate problems with applications. Leiden–Boston: VSP, 2005. 234 p. 13. Горшков А.В. Улучшение разрешения изображений при обработке данных физического эксперимента и нахождение неизвестной аппаратной функции по программам пакета REIMAGE // Приборы и техника эксперимента. 1995. № 2. С. 68–78. 14. Donatelli M., Estatico C., Martinelli A., Serra-Capizzano S. Improved image deblurring with anti-reflective boundary conditions and re-blurring // Inverse Problems. 2006. V. 22. P. 2035–2053. 15. Engl H.W., Hanke M., Neubauer A. Regularization of inverse problems. Dordr echt: Kluwer, 1996. 328 p. 16. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. Киев: Наук. думка, 1986. 544 с. 17. Воскобойников Ю.Е., Литасов В.А. Устойчивый алгоритм восстановления изображения при неточно заданной аппаратной функции // Автометрия. 2006. Т. 42. № 6. С. 3–15. 18. Остриков В.Н., Плахотников О.В. Идентификация функции рассеяния точки канала наблюдения по калибрующему изображению посредством метода наименьших квадратов // Оптический журнал. 2006. Т. 73. № 2. С. 26–30. 19. Захаров Д.Д., Сизиков В.С., Шемплинер В.В., Щекотин Д.С. Новые способы устранения артефактов на томографических и иных изображениях // Научно-техн. вестник ИТМО. 2006. В. 32. С. 138–143.
48 “Оптический журнал”, 80, 1, 2013
УДК 621.397.3
ПОСЛЕДОВАТЕЛЬНОСТЬ ОПЕРАЦИЙ ПРИ ФИЛЬТРАЦИИ ШУМОВ НА ИСКАЖЕННЫХ ИЗОБРАЖЕНИЯХ
© 2013 г.
В. С. Сизиков*, доктор техн. наук; Р. А. Экземпляров**, аспирант ** Санкт-Петербургский национальный исследовательский университет информационных ** технологий, механики и оптики, Санкт-Петербург
** Санкт-Петербургский государственный политехнический университет, Санкт-Петербург
** Е-mail: sizikov2000@mail.ru, rexe@yandex.ru
Рассмотрены вопросы, связанные с зашумлением смазанных или/и дефокусированных изображений. Проанализирована последовательность выполнения фильтрации шумов на таких изображениях – до устранения смазывания/дефокусирования или после него. Введены понятия предшествующей и последующей фильтрации шумов. Устранение смазывания/дефокусирования ряда изображений выполнялось методами параметрической фильтрации Винера и регуляризации Тихонова, а фильтрация шумов – методами адаптивной фильтрации Винера и медианной фильтрации.
Ключевые слова: зашумленное смазанное или дефокусированное изображение, предшествующая и последующая фильтрация шумов, устранение смазывания и дефокусирования.
Коды OCIS: 030.4280, 100.0100
Поступила в редакцию 26.06.2012
Постановка задачи
Искаженные (смазанные, дефокусированные) изображения обычно зашумлены [1–6]. Цель данной работы – анализ следующих вопросов: способны ли такие устойчивые методы, как метод параметрической фильтрации Винера, метод регуляризации Тихонова и др. [1–3, 5–10] не только устранять (в той или иной мере) смазывание или дефокусирование изображения, но и фильтровать шумы; в какой степени способствуют фильтрации шумов такие специальные методы, как метод адаптивной фильтрации Винера, метод медианной фильтрации и др. [1–6]; какова эффективность перечисленных методов в зависимости от типа шума (гауссов, импульсный и др. [1–6]); какова наиболее эффективная последовательность (очередность) фильтрации шума (до или после устранения смаза/дефокусирования) [11] – основной вопрос и др. Данные вопросы рассмотрены путем проведения экспериментальной проверки с получением количественных оценок качества восстановленных изображений при различных типах шумов и очередности их фильтрации на репрезентативных выбор-
ках. Первые три вопроса в той или иной мере рассмотрены в публикациях [1–13] и при рассмотрении 4-го (основного) вопроса результаты этих публикаций будут учтены.
Качественный анализ
Сначала обсудим вопрос о фильтрации шумов качественно. Если (равномерное и прямолинейное) смазывание изображения объекта произошло из-за сдвига камеры (например, сдвига фотоаппарата или рассогласования вращения телескопа и небесной сферы), а между объектом и камерой имелся импульсный шум (пыль, капли влаги и т. п.), то будут смазанными и объект, и шум (все точки объекта и шумовые точки превратятся в штрихи, рис. 1б). В этом случае логично сначала устранить смазывание объекта и шума некоторым методом (параметрической фильтрации Винера, регуляризации Тихонова и др.). Смазанный импульсный шум должен превратиться практически в точки. Затем нужно выполнить последующую фильтрацию шума (адаптивной фильтрацией Винера, медианной фильтрацией и др.). Заметим, что в данном случае по шумовым
“Оптический журнал”, 80, 1, 2013
39
(а) (б)
(в) (г)
Рис. 1. а – точное изображение w–, б – изображение, зашумленное импульсным шумом и смазанное с шумом, с размытием краев, в – устранение смазывания методом параметрической фильтрации Винера (K = 10–4) и методом квадратур с регуляризацией Тихонова (a = 10–6), г – последующая медианная фильтрация изображения (в) с маской 3×3, srel = 0,105.
томобиля [12] и др.), а в среде между объектом и неподвижной камерой имелся шум (туман, пыль, капли, снежинки и др.). В результате на изображении объект будет смазанным, а шум – неподвижным. В этом случае нужно также в ыполнить сначала предшествующую фильтрацию шума, а затем устранить смазывание.
Наконец, если изображение объекта несмазано и недефокусировано, а импульсный шум смазан (дождь, падающие снежинки и др.), то целесообразно выполнить лишь фильтрацию шума, например, медианным фильтром.
Могут существовать и другие ситуации с зашумлением, смазыванием или/и дефокусированием [1–3, 5, 6], причем шум может быть различных типов (импульсный, гауссов, рав номерный, мультипликативный и др.), изображение (а также шум) – черно-белым и цветным, а искажение – более сложным, чем смазывание или дефокусирование [2, 13, 14].
(а) (б)
штрихам можно оценить смаз D (и его направление), т. е. шум способствует определению функции рассеяния точки (ФРТ). Если же происходил не сдвиг камеры, а неправильная установка фокуса, то будут дефокусированными и объект, и шум (все точки превратятся в размытые пятна в виде ФРТ, рис. 2б). В этом случае нужно сначала устранить дефокусирование объекта и шума, а затем выполнить последующую фильтрацию шума.
Если изображение объекта смазалось за счет сдвига камеры (или дефокусировалось), а импульсный шум возник в результате сбоя в работе сенсоров матрицы ПЗС, то изображение объекта будет смазанным (или дефокусированным), а импульсный шум представлен в виде так называемых “битых пикселов” (рис. 3а и 4а). В этом случае естественно сначала выполнить предшествующую фильтрацию шума, а затем устранить смазывание (или дефокусирование) изображения объекта. Аналогичную последовательность операций следует соблюсти и в случае, когда изображение смазалось за счет движения объекта (самолета [5, 6, 12], ав-
40
(в) (г)
Рис. 2. а – точное изображение Луны w–, б – изображение, зашумленное импульсным шумом и дефокусированное вместе с шумом, в – изображение (б), рефокусированное методом двумерного ПФ с регуляризацией Тихонова, a = 10–8 (изображение и импульсный шум восстановились), г – последующая медианная фильтрация изображения (в) с маской 3×3, srel = 0,0259.
“Оптический журнал”, 80, 1, 2013
(а) (б)
(а) (б)
(в) (г)
(в) (г)
(д) (е)
Рис. 3. а – смазанное изображение Сатурна с четырьмя спутниками и с наложенным несмазанным импульсным шумом, б – обработка изображения (а) методом параметрической фильтрации Винера, K = 10–1, srel = 0,192, в – предшествующая медианная фильтрация изображения (а) с маской 3×3, г – параметрическая фильтрация Винера изображения (в) (K = 10–3, srel = 0,050).
Заметим, что последовательность фильтрации шума в вышеописанных ситуациях может казаться естественной и очевидной, однако приведенные ниже результаты численных экспериментов говорят о том, что отмеченный порядок не всегда является оптимальным.
Основные математические соотношения
Как прямая задача моделирования искажений изображений, так и обратная задача восстановления изображений описывается соотношениями [5–10]
¥
ò h(x - ξ)wy (ξ)dξ = gy (x) + dg,
-¥
¥¥
ò ò h(x -ξ, y -h)w(ξ, h)dξdh = g(x, y) + dg,
-¥ -¥
(1) (2)
где h – ФРТ, аппаратная функция, обычно пространственно-инвариантная (разностная),
Рис. 4. а – смазанное изображение и несма-
занный импульсный шум, б – изображение,
аналогичное (а), но с размытием краев, в –
обработка изображения (а) методом пара-
метрической фильтрации Винера, K = 10–1,
srel = 0,253, г – предшествующая медианная фильтрация изображения (а) с маской
3×3, srel = 0,243, д – параметрическая трация Винера изображения (г), K =
филь10–2,5,
srel = 0,139, e – восстановление методом квадратур с регуляризацией Тихонова изображе-
ния, близкого (г), a = 10–3, srel = 0,133.
w и g – распределения интенсивности по неискаженному и искаженному изображениям соответственно, dg – помеха (шум). В выражении (1) ось x направлена вдоль смаза, а y – перпендикулярно смазу (играет роль параметра).
Если рассматривается обратная задача, то соотношения (1) и (2) являются интегральными уравнениями (ИУ) Фредгольма I рода. Их решение есть некорректная задача [5, 6, 8, 15, 16], поэтому будем использовать для их решения устойчивые методы – методы регуляризации Тихонова и параметрической фильтрации
“Оптический журнал”, 80, 1, 2013
41
Винера [2, 3, 5–10, 15, 16], а решения другими методами (Люси–Ричардсона, “слепой” деконволюции, итераций и проч.), полученные в других работах [1–3, 5, 6, 13], также будем учитывать.
Решение ИУ типа свертки (1) методом одномерного преобразования Фурье (ПФ) с регуляризацией Тихонова имеет вид [5–10, 12, 15–17]
òway (ξ) =
¥
1 2p
-¥
H(-w)Gy (w) H(w) 2 + aw2 p
exp(-iwξ)dw,
где
¥
H(w) = ò h(x)exp(iwx)dx,
-¥
¥
òGy (w) = gy (x)exp(iwx)dx
-¥
(3) (4)
– спектры ПФ, a > 0 – параметр регуляризации, p ≥ 0 – порядок регуляризации. Для решения ИУ (1) согласно соотношениям (3), (4) авторами разработана в рамках системы программирования MatLab7 m-функция desmearingf.m [5, 6].
ИУ (1) может быть решено также методом квадратур с регуляризацией Тихонова [5, 6, 8, 9, 16] (исследования показали [5, 6, 9], что такой метод дает обычно более точное восстановление изображений, чем метод ПФ с регуляризацией). Согласно этому методу, ИУ (1) при каждом значении y приводится к системе линейных алгебраических уравнений (СЛАУ) вида Awy = gy, где A – матрица, wy и gy – векторы, а решение СЛАУ методом регуляризации Тихонова выглядит как (при каждом y)
way = (aI + AT A)-1 AT gy,
(5)
где I – единичная матрица, а AT – транспони рованная матрица (m-функция desmearingq.m). Решения (3)–(5) были использованы для задачи восстановления смазанных изображений (так называемый вариант 2 решения (3) и вариант 4 решения (5) [5, 6, 9]).
Что же касается задачи восстановления дефокусированных изображений, то она описывается двумерным ИУ (2) и ее решение методом двумерного ПФ с регуляризацией Тихонова представляется в виде (ср. (3), (4)) [5–8, 10, 15–17]
ò òwa (ξ, h) =
1 4p2
¥¥
Wa (w1, w2)´
-¥ -¥
´exp[-i(w1ξ + w2h)]dw1dw2,
(6)
где
Wa (w1, w2) =
H*(w1, H(w1, w2)
w2)G(w1, w2) 2 + a(w12 + w22
)
p
,
(7)
H(w1, w2) =
¥¥
= ò ò h(x, y)expi(w1x + w2y)dxdy, -¥ -¥
(8)
G(w1, w2) =
¥¥
= ò ò g(x, y)expi(w1x + w2y)dxdy -¥ -¥
(9)
(m-функция refocusingT.m), а решение методом параметрической фильтрации Винера принимает вид (ср. (6)–(9)) [2, 3, 5–8, 10]
ò ò=
1 4p2
¥¥ -¥ -¥
WK
(w1,
wK (ξ, h) =
w2 ) exp[-i(w1ξ
+
(10)
w2h)]dw1dw2,
где
WK
(w1,
w2 )
=
H*(w1, w2)G(w1, w2) H(w1, w2) 2 + K
,
(11)
K ≥ 0 – некоторая константа (параметр, дающий оценку отношению шум/сигнал по мощности). Метод параметрической фильтрации Винера реализован в системе MatLab7 в виде m-функции deconvwnr.m [3].
Заметим, что метод параметрической фильтрации Винера и m-функция deconvwnr.m часто используются [2, 3] не только для восстановления дефокусированных, но и смазанных изображений, причем в обоих случаях в двумерном виде (8)–(11).
В данной работе рассмотрены два типа шумов – гауссов и биполярный импульсный (типа “salt & pepper”), а также специальные фильтры для фильтрации шумов – среднеарифметический фильтр, метод адаптивной фильтрации Винера и медианный фильтр [1–4]. Что же касается других типов шумов (равномерный, мультипликативный) и методов их фильтрации (адаптивный медианный фильтр, ранговый фильтр и др.), то результаты их использования в различных публикациях [1–4] также учтены в данной работе.
В заключение данного пункта остановимся на вопросе оценки обработки изображений. К ачественно оценивать искаженные и восстановленные изображения будем визуально. А количественно степень искажения и, особенно, погрешность восстановления изображения
42 “Оптический журнал”, 80, 1, 2013
будем характеризовать относительным среднеквадратическим отклонением – СКО изображения w~ от точного изображения w– [5, 6, 9, 10]
srel =
w -w L2 || w ||L2
=
å å( )m n wji -wji 2
j=1 i=1
.
å åm n wj2i
j=1 i=1
(12)
Под w~ будем подразумевать искаженное
(смазанное, дефокусированное, зашумленное) изображение w~, а также изображение, восста-
новленное методом параметрической фильтра-
ции Винера wK или методом регуляризации Тихонова wa. Для вычисления srel согласно формуле (12) авторами разработана m-функция
rmsd.m [5, 6]. Будем обозначать через Kopt и aopt значения K и a, дающие минимум srel.
Вместо формулы (12) можно использовать
выражения
srel = || w -w ||L1 || w ||L1 , s = || w -w ||L2 , s = || w -w ||L1 ,
( )PSNR =10lg wmaxmn || w -w ||2L2
и др. [5, 6]. Однако в данной работе будем поль-
зоваться выражением (12) – типа (безразмерного) отношения шум/сигнал.1
1 1. Приведенные пользовать лишь
ввыслруачжаеен, икяогддлаяиsзrоeбl риаsжемноиженw–о
исиз-
вестно, т. е. при обработке модельных изображений
(фантомов), что и было сделано в данной работе.
Специальный же, весьма сложный вопрос об оценке
погрешности восстановления изображений без знания w– в результате решения некорректной задачи
[12, 15–17] в данной работе не рассматривается.
2. Можно использовать также некоторые способы выбора параметра регуляризации a и параметра K (без знания w–): способ невязки [8, 15–17] (нужно знать погрешность ||dg|| измерения искаженного изображения g), обобщенный принцип невязки [16] (нужно знать ||dg|| и погрешность знания ФРТ ||dh||), способ псевдо обратного оператора [16] и др. Однако эти способы часто дают завышенные по сравнению с aopt и Kopt значения a и K, а значит, заглаженные решения wa и wK. Кроме того, они требуют знания погрешностей ||dg|| и || dh||, которые обычно известны неточно. Поэтому использование способов выбора a и K отложим на дальнейшее.
3. В данной работе будем считать, что параметры искажений изображений (величина смаза D, угол смаза q, уровень и тип шума и проч.) известны точно. Вопросы, связанные с неточностями параметров искажений, специально рассмотрены в различных публикациях ([2, 5, 6, 13, 17, 18] и др.). Мы же, полагая, что w– и параметры искажений известны точно, будем тем самым определять потенциальные возможности рассматриваемой в данной работе методики.
Численные результаты
Ниже приведены результаты обработки изо-
бражений нескольких объектов. Для каждо-
го объекта приведено точное (неискаженное)
изображение. Это сделано специально, чтобы
оценить количественно качество обработки
(посредством соотношения (12)). Затем вноси-
лись искажения, выполненные искусственно
(программным путем). Что касается обработ-
ки изображений, когда истинное изображение
неизвестно, а искажения имеют естественную
природу, то такая обработка была также вы-
полнена [5, 6] и в заключительных выводах
у чтена.
Смазанное изображение и несмазанный
шум. На рис. 4а представлено черно-белое изо-
бражение girl.jpg 256×256, смазанное (D = 20,
q = 27°) и зашумленное несмазанным импульс-
ным шумом (доля зашумленных пикселов
d = 0,02), а на рис. 4б – аналогичное изобра-
жение, но с размытием краев (для уменьшения
эффекта Гиббса [5, 6, 9, 10]). Далее выполня-
лось устранение смазывания методами параме-
трической фильтрации Винера при различных
K и квадратур с регуляризацией Тихонова при
различных a, а также фильтрация шума меди-
анным фильтром и адаптивным фильтром Ви-
нера с различными масками, причем фильтра-
ция шума выполнялась как предшествующая
(pre), так и последующая (post).
На рис. 5 представлены полученные зави-
симости относительных погрешностей восста-
новления изображения srel(K) (под w~ в форму-
ле (12) подразумевается нию (10)) и srel(a) (под w~
wK согласно выражев формуле (12) подраз-
умевается wa согласно выражению (5)) для слу-
чая, когда фильтрация импульсного шума вы-
полнялась медианным фильтром (с помощью
m-функции medfilt2.m [3, 4]). Точное изображение w– представлено на рис. 1а. По миниму-
мам кривых srel(K) и srel(a) определялись Kopt и aopt, а также srel(Kopt) и srel(aopt) – значения, при которых изображение восстанавливается
с наименьшей погрешностью согласно коли-
чественному критерию (по визуальному кри-
терию “наилучшему” восстановлению могут
соответствовать несколько другие K и a).
Из приведенного рисунка видно, что значе-
ния погрешностей srel(Kopt) и srel(aopt) заметно зависят от последовательности выполнения
фильтрации шума (предшествующая или по-
следующая фильтрация), а также от размера
маски медианного фильтра. В данном примере
“Оптический журнал”, 80, 1, 2013
43
srel(K), srel(a)
1
3Wpost 0,8
1W 3Tpost
1T
9preW 9Wpost 0,6
0,4 3preW 9Tpost
0,2 3preT9preT
0 –6 –5 –4 –3 –2 –1 0
lgK, lga
Рис. 5. Погрешность восстановления изображения на рис. 4а и 4б: srel(K) – методом параметрической фильтрации Винера (W) и srel(a) – методом квадратур с регуляризацией Тихонова (T) с использованием предшествующей (pre) и последующей (post) фильтрации импульсного шума медианным фильтром с масками 1×1, 3×3 и 9×9.
предшествующая фильтрация импульсного шума при умеренной маске 3×3 ведет к погрешности srel, примерно в 1,5 раза меньшей, чем последующая фильтрация, что согласуется с вышеприведенным качественным анализом. Представленные на рис. 5 зависимости показывают, что метод регуляризации Тихонова дает чуть меньшую погрешность srel, чем метод параметрической фильтрации Винера.
Заметим, что при маске 1×1 медианная фильтрация посуществу не выполняется, поэтому кривая 1W на рис. 5 аналогична кри-
вым 1preW и 1Wpost, а 1T аналогична 1preT и 1Tpost.
На рис. 4в приведен результат параметрической фильтрации Винера после медианной фильтрации шума с маской 1×1. Он показывает, что параметрическая фильтрация Винера не только устраняет смазывание изображения, но и подавляет шум (K = 10–1, srel = 0,253), однако делает это лишь частично (метод квадратур с регуляризацией Тихонова дает аналогичный результат). Значительно более точный результат получается, когда выполняется предшествующая медианная фильтрация с бо′льшим размером маски (рис. 4г при маске 3×3), а затем устраняется смазывание (рис. 4д, е).
Данный пример был промоделирован также при использовании адаптивного фильтра Винера (с помощью m-функции wiener2.m [3, 4]) применительно к фильтрации как импульсного, так и гауссова шумов. В табл. 1 приведены результаты для импульсного шума, а в табл. 2 – для гауссова шума.
Данные, приведенные в табл. 1, показывают (в дополнение к рис. 4 и 5), что если шум импульсный, то медианный фильтр фильтрует его заметно лучше, чем адаптивный винеровский (а также среднеарифметический) фильтр: погрешность s примерно в 1,8–1,9 раза меньше в случае предшествующей фильтрации и умеренной маски 3×3, в 1,5 раза меньше в случае маски 5×5 (результат не приведен в табл. 1) и в 1,2 раза меньше в случае маски 9×9.
В табл. 2 приведены аналогичные результаты для случая гауссова 1%-го шума. Они показывают, что погрешность фильтрации s гауссова шума адаптивным фильтром Винера чуть меньше (примерно в 1,1 раза в данном
Таблица 1. Значения K = Kopt, s = srel(Kopt), a = aopt и s = srel(aopt)
Маска 1×1
Маска 3×3
Маска 9×9
Смазанное изображение girl.jpg
(D = 20, q = 27°), несмазанный Метод пар. Метод кв. Метод пар. Метод кв. Метод пар. Метод кв.
импульсный шум (d = 0,02) фильтрации с регуляр. фильтрации с регуляр. фильтрации с регуляр.
Винера Тихонова Винера Тихонова Винера Тихонова
Медианный предшествующая K = 10–1
фильтр
фильтрация
s = 0,253
a = 10–1 s = 0,252
K = 10–2,5 s = 0,139
a = 10–3 s = 0,133
K = 10–2 s = 0,208
a = 10–2,5 s = 0,187
последующая фильтрация
K = 10–1 s = 0,253
a = 10–1 s = 0,252
K = 10–1,5 s = 0,207
a = 10–1,5 s = 0,214
K = 10–2 s = 0,201
a = 10–2 s = 0,204
Адаптивный предшествующая
фильтр
фильтрация
Винера
последующая
фильтрация
K = 10–1 s = 0,253
K = 10–1 s = 0,253
a = 10–1 s = 0,253
a = 10–1 s = 0,253
K = 10–1 s = 0,250
K = 10–1 s = 0,237
a = 10–1 s = 0,251
a = 10–1 s = 0,244
K = 10–1,5 s = 0,236
K = 10–1,5 s = 0,221
a = 10–1,5 s = 0,236
a = 10–1,5 s = 0,227
44 “Оптический журнал”, 80, 1, 2013
Таблица 2. Значения K = Kopt, s = srel(Kopt), a = aopt и s = srel(aopt)
Маска 1×1
Маска 3×3
Маска 9×9
Смазанное изображение girl.jpg
(D = 20, q = 27°), несмазанный Метод пар. Метод кв. Метод пар. Метод кв. Метод пар. Метод кв.
1%-ый гауссов шум
фильтрации с регуляр. фильтрации с регуляр. фильтрации с регуляр.
Винера Тихонова Винера Тихонова Винера Тихонова
Медианный предшествующая
фильтр
фильтрация
K = 10–2,5 s = 0,126
a = 10–3 s = 0,117
K = 10–2,5 s = 0,140
a = 10–3 s = 0,129
K = 10–2 s = 0,207
a = 10–2,5 s = 0,185
последующая фильтрация
K = 10–1 s = 0,126
a = 10–3 s = 0,117
K = 10–3 s = 0,120
a = 10–3,6 s = 0,114
K = 10–3,2 s = 0,175
a = 10–3,4 s = 0,167
Адаптив- предшествующая
ный фильтр фильтрация
Винера
последующая
фильтрация
K = 10–1 s = 0,126
K = 10–1 s = 0,126
a = 10–3 s = 0,117
a = 10–3 s = 0,117
K = 10–3 s = 0,129
K = 10–3,5 s = 0,109
a = 10–3,4 s = 0,117
a = 10–3,7 s = 0,105
K = 10–2,3 s = 0,194
K = 10–3,6 s = 0,145
a = 10–2,8 s = 0,177
a = 10–3,8 s = 0,139
прим ере), чем медианным фильтром. Это подтверждается и расчетами в случае гауссова 5%-го шума.
Однако табл. 1 и 2 демонстрируют и неожиданный результат: если шум несмазанный импульсный и используется адаптивный фильтр Винера (см. табл. 1), то последующая фильтрация шума дает чуть меньшую погрешность s, чем предшествующая фильтрация (примерно в 1,05 раз для маски 3×3 в данном примере). Еще значительнее этот эффект в случае несмазанного гауссова шума и использования как медианного фильтра, так и адаптивного фильтра Винера (примерно в 1,1–1,15 раз в случае маски 3×3, табл. 2).
Смазанное изображение и смазанный шум. На рис. 1 представлены результаты обработки изображения girl.jpg в случае, когда изображение зашумлено импульсным шумом и смазано вместе с шумом (рис. 1б). Для устранения смазывания изображения использованы методы параметрической фильтрации Винера, а также квадратур с регуляризацией Тихонова и при каждом значении параметров K и a выполнена последующая медианная фильтрация с масками 1×1, 3×3, 9×9 и др. В результате построены кривые типа рис. 5 и получены следующие значения K = Kopt, s = srel(Kopt), a = aopt и s = srel(aopt): при маске 1×1 – K = 10–2,6, s = 0,148, a = 10–3, s = 0,141; при маске 3×3 – K = 10–4, s = 0,106, a = 10–6, s = 0,105; при маске 9×9 – K = 10–3,7, s = = 0,165, a = 10–6, s = 0,165. На рис. 1в приведено и зображение, полученное после устранения смазывания методами параметрической фильтрации Винера и квадратур с регуляриза-
цией Тихонова (результаты для обоих методов практически одинаковы), а на рис. 1г – изображение, полученное в результате последующей медианной фильтрации шума с маской 3×3. Рис. 1в демонстрирует хорошее восстановление изображения объекта и шума, а рис. 1г – окончательное восстановление объекта практически без смазывания и шума.
Дефокусированное изображение и дефокусированный шум. На рис. 2а представлено точное изображение Луны w– (файл moon. tif 455×325), а на рис. 2б – ее дефокусированное изображение (ФРТ задана в виде однородного диска радиуса r = 10 пкс), зашумленное импульсным шумом с долей зашумления d = 0,02 (шум также дефокусирован и поэтому на рис. 2б невидим). Было устранено дефокусирование методом двумерного ПФ с регуляризацией Тихонова согласно выражениям (6)–(9), а также методом параметрической фильтрации Винера согласно соотношениям (8)–(11) для ряда значений параметров a и K и при каждом a и K выполнена последующая фильтрация шума медианным фильтром с масками 1×1, 3×3, 9×9 и др.
На рис. 6 представлены полученные зависимости srel(K) и srel(a) – относительные погрешности восстановления изображения Луны. По минимумам этих зависимостей были определены оптимальные значения K = = Kopt, s = srel(Kopt), a = aopt и s = srel(aopt): при маске 1×1 – K = 10–2,9, s = 0,0794, a = 10–0,5, s = 0,0977; при маске 3×3 – K = = 10–3,5, s = 0,0691, a = 10–8, s = 0,0259; при маске 9×9 – K = 10–4, s = 0,0640, a = 10–8, s = 0,0549.
“Оптический журнал”, 80, 1, 2013
45
srel(K), srel(a)
0,5
0,4
0,3
3Wpost
1W
0,2 1T
9Wpost
0,1 9Tpost 3Tpost
0 –8 –7 –6 –5 –4 –3 –2 –1 0
lgK, lga
Рис. 6. Погрешности восстановления изображения Луны (см. рис. 2). Обозначения такие же, как на рис. 5.
На рис. 2в представлено изображение Луны, рефокусированное методом двумерного ПФ с регуляризацией Тихонова при a = 10–8, s = 0,2156 (видно, что изображение Луны и импульсный шум восстановились), а на рис. 2г – изображение, полученное в результате последующей фильтрации импульсного шума медианным фильтром с маской 3×3 (s = 0,0259). Восстановление изображения методом параметрической фильтрации Винера заметно хуже (погрешность s = 0,0691, т. е. в 2,7 раз больше, чем методом Тихонова).
Смазанное изображение с мелкими деталями и несмазанный шум. Рассмотрим изображение, представленное на рис. 3а, аналогичное изображению на рис. 4а, поскольку оба изображения – это смазанные изображения с наложенным несмазанным импульсным шумом. Однако изображение на рис. 3а отличается наличием мелких деталей, а именно, четырех маленьких спутников Сатурна, слабо видимых на фоне шума. Изображение saturn.png 1500×1200×3 – цветное, преобразованное в серое 1500×1200, D = 40 пкс, q = 35°. На изображение наложен несмазанный импульсный шум с долей зашумления d = 0,02. Смазывание изображения Сатурна со спутниками произошло из-за рассогласования вращения телескопа и небесной сферы (за время довольно продолжительной экспозиции). А импульсный шум, как и на рис. 4а, возник в результате выхода из строя ряда сенсоров матрицы ПЗС или/и из-за электромагнитных помех. Другими словами, импульсный шум
в данном примере – это аппаратурный шум, но не шум в среде.
Была выполнена фильтрация шума медианным фильтром (а также адаптивным фильтром Винера) и устранено смазывание методами параметрической фильтрации Винера при различных K, а также квадратур с регуляризацией Тихонова при различных a, причем фильтрация шума выполнялась как предшествующая, так и последующая. Результаты получились качественно весьма похожими на рис. 5 и табл. 1, однако количественно относительные погрешности s заметно меньше.
На рис. 3б приведены результаты параметрической фильтрации Винера и регуляризации Тихонова без выполнения специальной (медианной или адаптивной винеровской) фильтрации шума при K = Kopt = a = aopt = = 10–1 и s = srel(Kopt) = srel(aopt) = 0,192. Видим, как и на рис. 4в, что параметрическая фильтрация Винера и регуляризация Тихонова устранили смазывание и подавили шум, но лишь частично. Более качественный результат получается при добавлении предшествующей медианной фильтрации (с маской 3×3 или большей). На рис. 3в приведен результат предшествующей медианной фильтрации изображения, представленного на рис. 3а, с маской 3×3 (видим, что импульсный шум отфильтровался, а спутники в виде штрихов сохранились). На рис. 3г – результат параметрической фильтрации Винера и регуляризации Тихонова по отношению к изображению на рис. 3в. Погрешности s получились значительно меньшими, чем на рис. 4д, е, 5 и в табл. 1, а именно, K = a = 10–3, s = 0,050, т. е. примерно в 2,7 раз меньше, чем на рис. 5 и в табл. 1. Рис. 3г показывает, что удалось отфильтровать шум и устранить смазывание не только крупного объекта – Сатурна, но и мелких объектов – его спутников.
Для повышения репрезентативности исследования были обработаны также другие изображения: черно-белые cameraman.tif (портретное изображение), textp.bmp, textn.bmp (тексты), plane.bmp (быстролетящая цель – самолет), а также цветные flower.bmp, pears. png, peppers.png, onion.png, football.jpg (предметы), auto.jpg (подвижный объект), pair.jpg, rehoitus.jpg, причем как модельные (фантомы), так и реальные (натурные) изображения [5–12], в том числе, томографические [19].
На основе этих результатов, а также результатов других публикаций ([1–4, 13, 14, 17–19]
46 “Оптический журнал”, 80, 1, 2013
и др.) можно сделать следующие выводы в отношение задачи фильтрации шумов на искаженных (смазанных или дефокусированных) изображениях.
Выводы
1. Если устранять смазывание/дефокусирование, а также фильтровать шум каким-то одним методом (параметрической фильтрации Винера, регуляризации Тихонова, максимального правдоподобия Люси–Ричардсона, “слепой” деконволюции и др. [1–3, 5, 6]), то этого может оказаться недостаточно для получения удовлетворительного результата (см. рис. 3б и 4в). Более точная обработка может получиться, если указанные методы использовать в основном для устранения смазывания/ дефокусирования, а для фильтрации шума дополнительно использовать такие специальные методы, как медианный фильтр, адаптивный медианный фильтр, среднеарифметический фильтр, ранговая фильтрация, адаптивный фильтр Винера и др. [1–4] (см. рис. 1г, 2г, 3г, 4д, е).
2. Погрешность обработки s зависит от последовательности выполнения фильтрации шума на смазанном/дефокусированном изображении – до устранения смазывания/дефокусирования (предшествующая фильтрация) или после него (последующая фильтрация). Например, если на смазанное изображение наложен несмазанный (аппаратурный) импульсный шум (рис. 3а, 4а), то нужно сначала выполнить предшествующую фильтрацию шума медианным фильтром (рис. 3в, 4г), а потом устранить смазывание (рис. 3г, 4д, е). Если же эти операции поменять местами (смазывание –
последующая фильтрация), то это может привести к заметному повышению погрешности s (см. табл. 1, медианный фильтр, маска 3×3).
3. Если импульсный несмазанный шум фильтровать адаптивным фильтром Винера, то последующая фильтрация шума дает чуть меньшую погрешность s, чем предшествующая фильтрация (табл. 1). Такой же результат получается для гауссова несмазанного шума (для любого фильтра, табл. 2). Это – несколько неожиданный результат, не согласующийся с вышеприведенным качественным анализом. Он связан, во-первых, с различием свойств импульсного и гауссова шумов: биполярный импульсный шум типа “salt & pepper” имеет интенсивности лишь 0 и 255 (импульсы), присутствует лишь в d-ой доле (d ∈ [0, 1]) точек изображения и замещает интенсивности на изображении, а гауссов шум имеет различные интенсивности, присутствует в каждой точке изображения и добавляется к нему. Во-вторых, медианная и винеровская адаптивная фильтрации имеют также различия: в медианной практически исключаются импульсы, а в винеровской адаптивной импульсы с большими значениями интенсивности сильно влияют на результат фильтрации.
4. Импульсный шум заметно лучше фильтруется медианным фильтром, чем адаптивным Винера или среднеарифметическим, а гауссов шум, наоборот, лучше фильтруется адаптивным фильтром Винера, чем медианным фильтром (см. табл. 1).
Эти выводы подтвердились обработкой иных изображений: cameraman.tif, textp.bmp, textn.bmp, plane.bmp, flower.bmp, pears.png, peppers.png, onion.png, football.jpg, auto.jpg, pair.jpg, rehoitus.jpg и др.
* * * * *
ЛИТЕРАТУРА
1. Lim J.S. Two-dimensional signal and image processing. New Jersey: Prentice Hall PTR, 1990. 694 p.
2. Гонсалес Р., Вудс Р. Цифровая обработка изображений. М.: Техносфера, 2006. 1072 с.
3. Гонсалес Р., Вудс Р., Эддинс С. Цифровая обработка изображений в среде MATLAB. М.: Техносфера, 2006. 616 с.
4. Дьяконов В., Абраменкова И. MATLAB. Обработка сигналов и изображений. СПб.: Питер, 2002. 608 с.
5. Сизиков В.С. Обратные прикладные задачи и MatLab. СПб.: Лань, 2011. 256 с.
6. Сизиков В.С. Интегральные уравнения и MatLab в задачах томографии, иконики и спектроскопии. Saarbrücken: LAP, 2011. 252 с.
“Оптический журнал”, 80, 1, 2013
47
7. Сизиков В.С., Белов И.А. Реконструкция смазанных и дефокусированных изображений методом регуляризации // Оптический журнал. 2000. Т. 67. № 4. С. 60–63.
8. Сизиков В.С. Математические методы обработки результатов измерений. СПб.: Политехника, 2001. 240 с. 9. Сизиков В.С., Римских М.В., Мирджамолов Р.К. Реконструкция смазанных и зашумленных изображений
без использования граничных условий // Оптический журнал. 2009. Т. 76. № 5. С. 38–46. 10. Сизиков В.С. Прием “усечение–размытие–поворот” для восстановления искаженных изображений // Оп-
тический журнал. 2011. Т. 78. № 5. С. 18–26. 11. Экземпляров Р.А., Сизиков В.С. Устранение смазывания и дефокусирования изображений с предваритель-
ной фильтрацией шумов в рамках системы MatLab // Мат. международ. конф. “XXXIX Неделя науки СПбГПУ”, ч. XIII (ИМОП) / Под ред. Рудского А.И. СПб.: СПбГПУ, 2010. С. 223–225. 12. Petrov Yu.P., Sizikov V.S. Well-Posed, Ill-Posed, and intermediate problems with applications. Leiden–Boston: VSP, 2005. 234 p. 13. Горшков А.В. Улучшение разрешения изображений при обработке данных физического эксперимента и нахождение неизвестной аппаратной функции по программам пакета REIMAGE // Приборы и техника эксперимента. 1995. № 2. С. 68–78. 14. Donatelli M., Estatico C., Martinelli A., Serra-Capizzano S. Improved image deblurring with anti-reflective boundary conditions and re-blurring // Inverse Problems. 2006. V. 22. P. 2035–2053. 15. Engl H.W., Hanke M., Neubauer A. Regularization of inverse problems. Dordr echt: Kluwer, 1996. 328 p. 16. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. Киев: Наук. думка, 1986. 544 с. 17. Воскобойников Ю.Е., Литасов В.А. Устойчивый алгоритм восстановления изображения при неточно заданной аппаратной функции // Автометрия. 2006. Т. 42. № 6. С. 3–15. 18. Остриков В.Н., Плахотников О.В. Идентификация функции рассеяния точки канала наблюдения по калибрующему изображению посредством метода наименьших квадратов // Оптический журнал. 2006. Т. 73. № 2. С. 26–30. 19. Захаров Д.Д., Сизиков В.С., Шемплинер В.В., Щекотин Д.С. Новые способы устранения артефактов на томографических и иных изображениях // Научно-техн. вестник ИТМО. 2006. В. 32. С. 138–143.
48 “Оптический журнал”, 80, 1, 2013