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

Технические науки

УДК 681.323

А. П. Бахрушин

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

В статье рассматриваются вопросы цифровой фильтрации зашумленных изображений. Описывается общие принципы построения цифровых фильтров с конечной импульсной характеристикой (КИХ-фильтры). Отмечается, что одно из существенных преимуществ КИХ-фильтров перед фильтрами с бесконечной импульсной характеристикой (БИХ-фильтров) заключается в возможности синтеза и практической реализации КИХ-фильтров с вещественными частотными характеристиками. Предлагается новый метод фильтрации зашумленных изображений на основе медианной и инверсной фильтрации. Экспериментальные результаты подтвердили эффективность предложенного метода.

Ключевые слова: цифровая фильтрация, КИХ-фильтры, медианные фильтры, инверсные фильтры.

Alexander P. Bahrushin. APPLICATION OF SPECTRAL ANALYSIS METHODS IN THE PROCESSING NOISY IMAGES(Far Eastern State Academy for Social and Humanity Studies).

This paper examines the problems of noised images digital filtering. The general ideas of digital filters building with final impulse response filters (FIR filters) are described. It is noted that one of the FIR filters advantages in front of filters with infinite impulse response (IIR filters) is the possibility of FIR filters synthesis and practical realization with real frequency responses. A new method of noised images digital filtering on base of median and inverse filtering is proposed. The experimental results have confirmed the efficiency of the proposed method.

Keywords: digital filtering, FIR filters, median filters, inverse filters.

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

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

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

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

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

Следует заметить, что математические методы, применяемые при синтезе цифровых фильтров, в значительной степени зависят от того, конечна длительность импульсной характеристики фильтра или бесконечна. Особенностью фильтра с конечной импульсной характеристикой, получившего название КИХ-фильтра, является то, что он не имеет полюсов, а имеет только нули. Так как импульсный отклик КИХ-фильтра содержит лишь ко-

нечное число ненулевых отсчетов, то он абсолютно суммируем, и, следовательно, КИХ-фильтры всегда устойчивы.

Отметим также, что одно из существенных преимуществ КИХ-фильтров перед фильтрами с бесконечной импульсной характеристикой (БИХ-фильтров) заключается в возможности синтеза и практической реализации КИХ-фильтров с чисто вещественными частотными характеристиками. Такие фильтры получили название фильтров с нулевой фазой. В частотной области условие нулевой фазы можно определить следующим образом:

Н(р) = Н* (р).

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

Рассмотрим прежде всего круг проблем, возникающих при реализации одномерного КИХ-фильтра с линейной фазовой характеристикой. на основе

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

N —1

к(I) = ^Н(и)ехр(і(2ж/N)иї) , I=0,...,N—1,

и = 0

Согласно данному выражению, для получения импульсной характеристики И(I) достаточно частотную ось разбить на N равноотстоящих точек и вычислить обратное ДПФ. Однако в действительности эта прямая процедура не представляет особого практического интереса из-за возникновения явления Гиббса, которое выражается в виде постоянного выброса и пульса-

ций амплитудно-частотной характеристики (АЧХ) фильтра. Эти пульсации обычно описываются с помощью функции sin х / х и некоторых интегралов от этой функции.

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

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

Запишем z - преобразование последовательности отсчетов импульсной характеристики h (l):

N -1

H( z) = X h (l) z -' = X

l = 0

n = 0

1 N- 1 2n

— У H

N„4-0

z =

l -1

N

l = 0

exp(i—u) z 1

N

H (u) 1 - exp (i 2 n u) z

N

- N

i í  2 n X -1 1 - exp ( i — u) z

N

u = 0

H

1 - exp

где г — комплексная переменная;

Н (г) — функция комплексной переменной.

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

произведем замену переменной z частотой о:

-і

на комплексную экспоненту exp(-ІО) с

H (о) =

N

У

u = 0

H

1-exp(-io)exp( i—u) N

N -1

u = 0

z

exp

( .NgI (. Nа Л ( . N а ЛІ

і * I exp * — exp і П I

V 2 ) V 2 ) V 2 )_

N

I

H

exp

if w

G її '

-------u

2 2

-j

V v

exp

if w

G її '

------u

22

г

VV

exp

)

(а її '

— г---------u

V V 2 2 У)

exp(— zN—11 JV_ jH (u)exp | — і—u I sin[(N g)/2]

V 2 _2 ^_____________V N )_____________

N

sin

—

Формула (1) обеспечивает возможность рекурсивного синтеза КИХ-фильтра [1] и соответствует последовательному соединению фильтров, причем первый из них представляет собой нерекурсивный фильтр с передаточной функцией:

н (z) =

1 — z

N

Данная функция имеет N нулей на единичной окружности в точках:

zu = exp

u = 0,..., N—1.

Разностное уравнение данного фильтра можно представить в виде:

5 (0 5 (I - N

y(l)=-

N

N

где ^ (I) — некоторая дискретная последовательность.

Как известно, такой фильтр получил название гребенчатого фильтра. Второй фильтр фактически представляет собой совокупность N параллельно соединенных рекурсивных КИХ-фильтров первого порядка:

н (z) = H (u )/

ґ

1— exp

.2жпІ — ,

і----I z

N )

имеющих полюс в точке exp(/(2я/N)) и .

Рассматривая соответствующую составляющую с полюсом z=exp (-i(2я/N)) u, получаем фильтр второго порядка с действительными коэффициентами.

Разностное уравнение такого фильтра имеет вид [2]:

X

u = 0

u = 0

у (І)=2 А (и) \ соб в (и) я(1) - СОБ

в (и)-—и

ж

^ (І - т + 2соб

2п ) —г

N

у

причем Н(и)=Л(и) ехр(/'0(и)), где как Л(и), так и 0(и) удовлетворяют условиям симметрии. Очевидно, что выполнение данного условия необходимо для того, чтобы импульсная характеристика была действительной (как это обычно и бывает на практике).

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

I, = «и

и=0,

где N может быть как четным, так и нечетным.

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

Запишем частотную характеристику фильтра в виде модуля и аргумента:

Н(и)=| Н(и)|ехр(і р (и)) ,

где | н(и)\ — амплитудно-частотная характеристика фильтра,

р(и) — фазо-частотная характеристика фильтра.

Тогда условия симметрии при нечетном N можно представить в следующем виде:

Н

р(и) =

при и = 0,...,(Ж-1)/2, и = (N-1)/2,..., N-1.

Подставив выражение (3) в (2) и произведя ряд преобразований, получаем частотную характеристику в виде:

н (в)=|нМб,П( ^+/2

N б,п(ю/2) N и=1

Н (и)|ехр

. 2п ( N-1)! ( • п ) •

- і — I ехр - і—и I БІП

_ N 1 2 1 N )

N а 2 .

БІп

+

1

N-1

+ - I

N u = (N - l)/2

\H (u) | exp

2п íлт \N

- i — (N - u)-N2

1] С • n V

exp -1—u I sin

1 N )

N a 2

sin

Перепишем полученное выражение, сделав замену m=N -u:

H wsiníN^/2) + INI/2

N sin (a /2) u= i

\H (tt)|exp

. 2n С N-1 ЇІ С • n V

-і—ы I exp -1—ы I sin

_ N 1 2 JJ 1 N )

N a 2

sin

-+

(Nd)/2 |H (N - m)| +1

exp

. 2n С N-1 Л]

-1—m 2 JJ exp

L n

n

N

u I Sin

m = 1

N

sin

H(0) sin (Na/2) (Nv1)/2 |H(u) I . Na

---------------------------+ I J-------------1 sin-----

N sin (a/2) u=1 N N

sin

Тогда с помощью тождества sin(a + un) = sin(a-n) = (- l)u sin a преобразуем частотную характеристику к классической форме, удобной для использования в технике синтеза цифровых фильтров:

jj, X |H(0)| sin(Na /2) (N^)/2 |H(u)| . N a

H (a) =----------------------+ I -------sin —

N sin (a/2) u 1 N N

sin [N (a /2 - nu / N)] sin N (a/2-n u / N)]'

sin

Если N четно, то прежде всего нужно обеспечить условие симметрии выборок импульсной характеристики, для выполнения которого необходимо начало координат разместить между отсчетами N /2 и N /2+1.

Тогда:

(р(ы ) =

2n СN-1

--ы I --

N I 2

2í(n - ы)СN-1

N К \ 2 0

при ы = 0, .. ., N /2 -1;

при ы = N/2 + 1, . . ., N-1; при N/2.

+

Заметим, что для фильтра с линейной фазой при четном N должно выполняться условие Н() = 0 при =л, чем и объясняется то, что

(р(N/2) = 0 и |И(N/2)| = 0 .

С учетом полученных выражений частотную характеристику (6) для случая четного N можно записать в следующем виде:

. |H(0)| sin (Na /2) sin [N(a /2 - nu /N)]

H (a) =------------------------------------------------;-1-;-+

N sin

sin [n (a/2-nu / N)] N^-1 |H (ы)\ . N a

+---- —------------------------------— + I J-1 sin .

sin (a/2-nu/N) ы = 1 N N

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

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

Таким образом, необходимо так подобрать число выборок в переходной полосе и их значения, чтобы ошибка аппроксимации частотной характеристики была минимальной. Решению этой проблемы посвящен ряд интересных публикаций. Так, например, в работе [3] предлагается использовать метод наименьшего спуска, а в работе [44] рассматривается более эффективный подход, основанный на сведении задачи аппроксимации к задаче линейного программирования.

Рассмотрим один из возможных подходов к решению проблемы повышения резкости зашумленных изображений на основе синтеза двумерного цифрового фильтра. С этой целью представим уравнение формирования изображения в виде [4]:

ад ад

gT

- ад - ад

где gT (x,y) — сформированное изображение;

h (x - xj, y - y1) — импульсная характеристика системы передачи изображения (системная функция);

^ (xj, y) — функция распределения яркости объекта.

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

С учетом аддитивного шума модель процесса формирования изображения (8) можно переписать в виде:

Ев (х У) =№

| |И(х - хр У - У!) 5 (%!, У!)йхх йух

+ п

где ^ (х, у) — фактически записанное изображение;

Ш — характеристика записывающего процесса; п (х, у) — помехи.

В приведенной модели предполагается, что шум п (х, у) не зависит от записанного изображения. В этом случае задача восстановления изображения, например, его резкости, сводится к определению исходного изображения ^ (хх, у) на основе зашумленного изображения ^ (х, у).

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

У— [&

где ** - двумерная свертка.

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

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

да — да

Ш 1 можно применять к слагаемым по отдельности:

1

что равносильно замене шума п(х,у) на новый шумовой процесс Ш_1 [п(х, у)].

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

gF

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

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

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

(IА х, к А у)« Ш

N - 1 N - 1

I I к [

р = 0 ц = 0

+ п (IА х, к А у),

где А х и А у - шаги дискретизации вдоль координатных осей I и к, соответственно.

Если выполняется предположение, сделанное при выводе соотношения (12), то без потери общности можно допустить, что если Ах = Ау = 1, то:

N - 1 N - 1

gF

р = 0 ^ = 0

Заметим, что операция дискретной свертки в выражении (13) имеет аналог в пространстве дискретного преобразования Фурье (ДПФ), и поэтому имеем:

0(ы, V) = Н (и, V) ^ (и, V) + N (и, V),

где и и V — номера пространственных частот вдоль координатных осей, причем и, V = 0,1, . . ., N -1.

Для определения системной функции Н (и, V) необходимо воспользоваться прямым преобразованием Фурье:

N - 1 Л - 1 / 2Ж

Н (и, V) = I I к(I, к) ехр [----------------(1и + kv)].

к = 0 I = 0 N

N - 1 N - 1

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

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

где 8п (и, V) и 5^ (и, V) — энергетические спектры шума и сигнала соответственно.

Как и в случае инверсной фильтрации, обработка изображения производится на основе двумерных ДПФ, а в результате выполнения обратного двумерного преобразование Фурье от (15) формируется скорректированное изображение. Из (14) следует, что если энергия шума в широком диапазоне частот мала, то винеровский фильтр практически переходит в обратный фильтр. Если же энергия шума в высокочастотной части спектра достаточно велика, то значение передаточной функции в этом диапазоне частот падает и, следовательно, фильтр перестает усиливать как высокочастотные составляющие шума, так и составляющие полезного сигнала. Очевидно, что изображение в этом случае неизбежно станет менее резким, а возможно, и вообще потеряет мелкие детали.

Н

Тогда:

£ (и, V) = Н^ (и, V) GF (и, V) .

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

где 5 — значение яркости I - го пиксела.

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

центральных элементов вектора ц , если N четно.

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

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

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

Для количественной оценки качества фильтрации зашумленного изображения использовался показатель PSNR (Peak Signal to Noise Ratio):

где M — число уровней квантования исходного изображения;

S — исходное изображение;

Si — зашумленное изображение;

S 2 — отфильтрованное изображение;

m и к — число пикселов исходного изображения по вертикали и горизонтали, соответственно;

В ходе эксперимента использовались тестовые изображения в виде “Лены”, “Барборы”, размера 512х512пикселов при числе уровней квантования 256.

Соотношение сигнал /помеха принималось равным 3.

Данные, полученные во время экспериментальных исследований, показали, что для зашумленных изображений показатель PSNR составлял величину порядка 20—30, в то время как для отфильтрованных изображений порядка 50—60. Сопоставление результатов исследования свидетельствуют о высокой эффективности предложенного метода цифровой фильтрации зашумленных изображений.

Литература

1. Rabiner I.R. and Shafer R.W. Recursive and nonrecursive realization of digital filters designed by frequency sampling techniques/ IEEE Trans. Audio Electroacoustics AU-20, 1, 1971, pp 38—56.

2. Cappellini V., Constantinides A.G. and Emitani U. Digital filters and their applications/ Academic Press, London, 1978. — 397 p.

3. Rabiner L.R., Gold D. and NcGoneral C.A. An approach to the approximation problem for nonrecursive digital filters. IEEE Trans. Audio Electroacoustics 18, 2, 1970, pp 65—96.

4. Оппенгейман Э. Применение цифровой обработки сигналов. М.: Мир, 1980. 552 с.

5. Hellstrom C.W. Image restoration by the method of least squares. Opt. Soc. Amer., 1967, pp. 297-303.

6. Хуанг Т. Быстрые алгоритмы в обработке изображений. М.: Радио и связь,

г

\

, n = 1,2,

1984. 271 с.