Матричное решение задачи оптимальной нелинейной фильтрации изображений Текст научной статьи по специальности «Физика»

Щербаков М.А., Сорокин С.В. МАТРИЧНОЕ РЕШЕНИЕ ЗАДАЧИ ОПТИМАЛЬНОЙ НЕЛИНЕЙНОЙ ФИЛЬТРАЦИИ ИЗОБРАЖЕНИЙ

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

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

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

тры; морфологические фильтры; фильтры, основанные на порядковых статистиках, и их разновидности: L-, R-, M-фильтры, медианные фильтры; расширенные фильтры Калмана; нейронные фильтры; полиномиальные фильтры.

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

тивно разрабатываются в настоящее время.

В общем случае цифровой полиномиальный фильтр размерности r и порядка M, определяется конечным дискретным рядом Вольтерра (функциональным полиномом) вида

ММ т

Хп)=/%+Х>’т(П) = /%+22---2Ат(П1’-"’Пт)ПХ(П_П<) ' (1)

m=1 m=1 ni nm i=1

где hm(n1, ..., пш) — многомерные импульсные характеристики (ядра) фильтра, зависящие от векторных аргументов n = [пц ... д*г]. Фильтры вида (1) часто называются также фильтрами (процессорами) Воль-

терра [2].

Выходной сигнал у(п) таких фильтров представляет собой сумму составляющих, характеризующих нелинейности различного порядка: первая ух(п) имеет вид линейной свертки, вторая У2(п) — квадратичной

свертки и т. д. Составляющая ym(n) фильтра, определяемая сверткой m-го порядка, является нелинейной относительно отсчетов входного сигнала, однако, остается линейной по отношению к коэффициентам фильтра.

При m = 1 ядро hm (ni) представляет собой обычную импульсную характеристику многомерного линейного

фильтра, в то время как при m = 2 , ..., M ядра hm(ni, ..., nm) можно рассматривать как импульсные характеристики высших порядков, характеризующие нелинейные свойства многомерных полиномиальных фильтров. Используя в качестве входного сигнала сумму пространственных 8-функций, импульсную характеристику hm(ni, ..., nm) m-го порядка можно интерпретировать аналогично одномерному случаю, рассматривая ее как составляющую реакции фильтра, обусловленную взаимодействием m пространственных импульсов.

Имеет место тесная взаимосвязь между многомерной линейной и полиномиальной фильтрацией. Как известно [1], линейная многомерная фильтрация сигнала u(ni, ..., ns) описывается многомерной линейной сверткой вида

да да

/...//. /. > . (2)

*1 = 0 is = 0

Как частные случаи, для s = 1 получаем линейную свертку для фильтрации одномерных сигналов, для s = 2 — линейную свертку, описывающую фильтрацию двухмерных полей, например, изображений и т. д.

Допустим теперь, что s-мерный входной сигнал является сепарабельной функцией, т. е. представим в виде произведения m сигналов меньшей размерности r = s/m. Используя векторные аргументы, запишем:

m

u (ni,...,nm ) = П - (nj ) . (3)

j=1

Для данного воздействия выражение (2) принимает вид

m

>1 >m j=1

Выделяя из выходного сигнала лишь диагональные блоки размерности г, т. е. полагая ni = n2 = ... = nm= n, получим:

m

Ут (П) = '= ■ ■ ■ = )ПХ(П “ ) •

>1 'm j=1

Это выражение есть не что иное как нелинейная свертка m-го порядка.

Таким образом, имеет место тесная взаимосвязь между многомерной линейной и полиномиальной фильтрацией, состоящая в следующем. Выходной сигнал нелинейного фильтра порядка m и размерности r может быть получен из реакции многомерного линейного фильтра (прототипа) размерности rm при сепарабельном воздействии вида (3) путем выделения из выходного сигнала данного фильтра лишь диагональных блоков размерности r. В частности, например, из шестимерного линейного фильтра в зависимости от представления s = 6 в виде произведения rm могут быть построены следующие нелинейные фильтры:

1. Нелинейный фильтр шестого порядка ( r = i, m = 6 )

6

u{nl,...,n6) = Y\x{nj) , j=1

У6(п) = У(пі,---,п6)

П=Пл = ... = Пї

2. Двухмерный нелинейный фильтр третьего порядка ( г = 2, т = 3 ) 3

и(пп,пи,и21,и22,«3!,П32) = Пх(п^,пі2) ,

у=1

у3(п) = у(пп,п12 ,п31,п32,п31,п32)

п = п, = П9 = 1Ь

3. Трехмерный нелинейный фильтр второго порядка ( г = 3, т = 2 ) 2

и(ип,«12,nlз,п21,п22,п23) =ППА,Иі,3) ,

у=1

У2(п) = у(пп>п12>пп>п21>п22>п2ъ)

п = п, = Пт

П1 п2

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

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

вивалентное матричное представление. Для перехода к матричной форме записи выражения (1) рассмотрение вектор входного сигнала

^=[*„(0) *„(1) *„(^-1)].

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

х(*) = X хп “ хп

г(* -1)

едем в

зектор-

Ут (П) = К Хп

г- (т)

(4)

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

Для двухмерных полиномиальных фильтров вместо векторной формы (4) можно использовать матричную. В этом случае опорная область фильтра определяется на двухмерной N х ^-решетке. Воспользуемся простым отображением для преобразования двухмерной индексации точек опорной области в одномерную. Для маски 3 х 3, например, это преобразование будет выглядеть так

(5)

"(0,0) (0,1) (0,2)' "0 1 2'

«2 = (1,0) (1,1) (1,2) 3 4 5

.(2,0) (2,1) (2,2). 6 7 8

Для такого представления опорной области ^2 входной сигнал может матрицы. Так например, для области (5) данная матрица принимает вид

быть записан в виде N х ^

0 *я хп (2)

Хп = 5? х(

6 ня ОО" х(

Здесь, в отличие от (4), входной сигнал не упорядочивается в вектор, а представляется в более естественном для двухмерной фильтрации матричном виде.

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

расположенные в соответствии с индексацией произведений отсчетов Хп (!±) •...-Хп (1т), содержащихся в матрице Хп(т). Тогда составляющая Ут (п) двухмерного полиномиального фильтра может быть представлена в следующей матричной форме:

Ут(п) = Ц Нга°Х«}, (6)

где О означает произведение Адамара (поэлементное произведение матриц), £{А} — сумму всех элементов матрицы А.

Следует заметить, что Ьт и Нт могут быть разбиты на группы одинаковых элементов. Действительно, без ограничения общности можно считать ядра Ьт(п1г ..., пт), т > 1, соответствующие коммутативным

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

мы можем оставить в (1) лишь члены, соответствующие различным сочетаниям (пх, ..., пт). При этом

общее количество уникальных членов для полиномиального нелинейного фильтра порядка М составит м

и.

^/-1т ___ґ-іМ

СЫ+т-1 = СЫ+М

(7)

где Ст — число сочетаний из п по т. В результате, многомерный полиномиальный фильтр, определяемый рядом (1), будет определяться Lмх 1 вектором коэффициентов

ьТ

ЬТ

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

Т

Хп =

п

т=0

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

Хп) = hrZn , (8)

линейной относительно вектора h, содержащего C^+м коэффициентов фильтра.

Используя матричное представление (8), задачу синтеза оптимального полиномиального фильтра общем случае можно сформулировать как задачу нахождения n х1 вектора коэффициентов hopt, минимизирующего некоторый функционал качества

F (h opt) = min F (hx (9)

hsG

hopt eG = {hER :g(h) = 0} , (10)

где F(h) — целевая функция, g(h) — векторная функция ограничений.

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

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

(ll)

где — область определения выходной реализации г-мерного сигнала

ЖГ = {П = (п^...,пг): 0 <п < К — 1; г = 1,...,г} ,

а | ^ | — количество элементов в ней.

Следует отличать от опорной области ^г фильтра. Согласно (7) размерность вектора Ь для поли-

номиального фильтра порядка М составит п = Ьм = С'М+м . Для изотропных фильтров Ьм будет определяться количеством классов эквивалентности и составит м

^м=1+21 ^т,

т=1

где |^т! — число элементов в классе эквивалентности для однородного фильтра т-го порядка.

Таблица х

Квадратичный изотропный фильтр с маской 3x3

m i Вектор h K Вектор х произведений входных отсчетов

0 0 h0 i i

l hi(0) 4 X0

l 2 hi(1) 4 Xi

З hi(4) i X4

4 h2(0, 0) 4 X0X0 + X2X2 + X6X6 + XeXe

5 h2(0, i) 16 2[X0 (Xi + xз) + X2 (Xi + X5) + X6 (X3 + X7) + Xg (X5 + X7)]

6 h2(0, 2) В 2(X2 + X6)(X0 + Xg )

7 h2(0, 4) В 2X4 (Xq + X2 + X4 + Xg)

В h2(0, 5) 16 2[X0 (X5 + X7) + X2 + X7) + X6 (Xi + X5) + Xg (Xi + xз)]

2 9 h2(0, В) 4 2(XQXg + X2X6)

l0 h2(i, i) 4 X1X1 + XзXз + X5X5 + X7X7

ll h2(1, З) В 2(Xi + X7 )(xз + X5)

l2 h2(i, 4) В 2X4 (Xi + Xз + X5 + X7)

із h2(1, 7) 4 2(XiX7 + XзX5)

і4 h2(4, 4) i X4X4

В частности, для изотропного квадратичного фильтра с маской 3 х 3 количество уникальных коэффициентов составит Ь2 = 1 + 3 + 11 = 15. Векторы Ь коэффициентов и х входных отсчетов для данного фильтра приведены в табл. 1 и соответствуют принятой нумерации элементов опорной области (5) ^2, где К - вес коэффициентов импульсных характеристик фильтра.

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

^0 = 0 2И1(І1) = 1 22к2(І1,І2) = 0 , (12)

*1 І1 І2

которые могут быть представлены в матричной форме АЬ = Ь . (13)

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

Г4 410000000000 0] Гг п А = , ЬТ = ' 01 .

(^^^^^^^^ 16 44884 1^^.

Целевая функция (11) может быть записана в матричной форме

F(h) = hT Rxh - 2hT rd

d Х "

(l4)

где Rx — автокорреляционная матрица произведений входных отсчетов RX “I т' I S ХпХп ;

Гсх — вектор взаимных корреляций между заданным выходным сигналом и произведениями отсчетов входного сигнала

г*х=т^л2 У(п)Хп

| А-1 пєа;

и га— среднеквадратическое значение заданного сигнала

ъ =^4 2 й2(п).

| Аг І пєА-

Таким образом, задача синтеза оптимального полиномиального фильтра сводится к минимизации квадратичной функции (14) на линейном подпространстве, определяемом условием (13). Решение этой задачи

хорошо известно [4] и достигается из любой начальной точки h[ нием

за один шаг в соответствии с выраже-

hopt = h + Z(ZTRzZ)-lZT(rdz — Rzh[0J) , (15)

где Z — матрица со столбцами, являющимися базисом нуль-пространства, образованного строками мат-

рицы А ограничений, т. е. удовлетворяющая условию AZ = 0.

а б в г

Рис. 1. Результаты обнаружения перепадов в зашумленном изображении: а — искаженное изображение,

используемое в качестве входного; б — эталонное изображение, содержащее идеальные перепады; в — результат обработки с помощью синтезированного оптимального фильтра; г — результат обработки с помощью оператора Собела

В качестве примера рассмотрим задачу синтеза квадратичного фильтра для обнаружения границ деталей изображения, устойчивого к воздействию шумов. Для этого воспользуемся синтезированным изображением размером 64 х 64, составленным из светлых и темных треугольников с уровнями яркости, равными соответственно 17 0 и 80. Данное изображение было искажено гауссовым шумом с дисперсией а2 = 4 00 и использовалось в качестве входного сигнала (рис. 1, а). Эталонное изображение, состоящее из выделенных перепадов, показано на рис. 1, б. Для его получения к исходному (неискаженному) изображению был применен известный оператор Собела с последующим сравнением с порогом.

При проектировании фильтров для обнаружения перепадов условия (12) сохранения уровня яркости не накладываются. Оптимальный по критерию (14) вектор h коэффициентов определялся из (15) при Z = I. На рис. 1, в, г приведены результаты обработки зашумленного изображения с помощью рассчитанного оптимального фильтра и оператора Собела, являющегося одним из наиболее эффективных среди известных детекторов перепада. Из сравнения полученных результатов видно, что синтезированный фильтр обладает

большей устойчивостью к шумам по сравнению с оператором Собела, что проявляется в существенно мень-

шем количестве ложных обнаружений.

ЛИТЕРАТУРА

1. Прэтт У. Цифровая обработка изображений. В 2 кн. — М.: Мир, 1982. — Кн. 1. — 312 с. Кн. 2. — 480 с.

2. Pitas I., Venetsanopoulos A. N. Nonlinear digital filters: principles and applications. —

Kluver Academic Publishers, 1990. — 391 p.

3. Щербаков М.А. Цифровая полиномиальная фильтрация: теория и приложение. — Пенза: Изд-во Пенз.

гос. техн. ун-та, 1997. — 246 с.

4. Гилл Ф., Мюррей У., Райт М. Практическая оптимизация. — М.: Мир, 1985. — 510 с.