Идентификация устойчивых двумерных фильтров Текст научной статьи по специальности «Математика»

ИДЕНТИФИКАЦИЯ УСТОЙЧИВЫХ ДВУМЕРНЫХ ФИЛЬТРОВ

Зимин Д.И., Фурсов В.А. Институт систем обработки изображений Самарский государственный аэрокосмический университет

Аннотация

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

Постановка задачи

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

1) выбор фрагментов, пригодных для определения характеристик искажений в канале регистрации изображения;

2) формирование из выбранных фрагментов тестовых образцов (компьютерное ретуширование изображений на фрагментах);

3) идентификация параметров искажающей системы и/или восстанавливающего фильтра по отобранным и тестовым (отретушированным) фрагментам;

4) выбор класса и уточнение параметров восстанавливающего фильтра по критериям вычислительной сложности и качества обработки;

5) выбор схемы декомпозиции исходного больше-форматного изображения;

6) обработка большеформатного изображения на многопроцессорной системе.

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

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

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

2 я (п1- т1, п2- т2)+

2 / (п1 - т1, п2 - т2) + п2 )

(1)

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

Поэтому задачу параметрической идентификации искажающей системы и/или восстанавливающего фильтра целесообразно решать в классе фильтров с бесконечной импульсной характеристикой (БИХ-фильтров), для которых связь отсчетов выходного Я (п1, п2) и входного / (п1, п2) изображений описывается соотношением:

Я (п1, П2 ) = - X

(т1 ,т2 )£0<

У ь

/ : т1 т

(т1,т2) £¡3/

где - а , Ь 2 - подлежащие определению коэффициенты фильтра, ^(п1, п2) - дискретный шум. При использовании БИХ-фильтров вида (1) параметры ат т , Ьт т даже в случае значительных

А т1,т2 у т1,т2

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

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

2 )_У ЬЩ

Н

)_ Р (Z1, 72) _

„т1 т2

2 V V

„т, т2

, V V

(2)

О (zl, 22) У а

устойчив, если многочлен О (, 72) не обращается в нуль ни при каких (7Х, 72) из области

Б = {(71, г): 1^1 < 1 п|2г\< 1} . (3)

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

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

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

Вначале рассмотрим простой случай, когда импульсная реакция фильтра:

Р ^ ) I Ьщт2 2?

H (z„ z2 ) =

g (z„ z2) x zm

(4)

обладает радиальной симметрией. В этом случае опорную область двумерного БИХ-фильтра обычно задают в виде квадрата с центральным отсчетом в точке пх, п2 (рис. 1).

Рис. 1 Опорная область БИХ фильтра Q

Рис. 2 Опорная область с системой окружностей

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

каждой отдельно взятой окружности одинакова, а координаты этих точек удовлетворяют условию mf + да2 = Const.

Пронумеруем окружности в порядке возрастания их радиусов: m = 0,1,2,... так, что m = 0 соответствует точке в центре опорной области. Теперь импульсную реакцию фильтра (4) можно записать в виде

H (z) = x hmZm.

(5)

Ясно, что ей будет соответствовать некоторая одномерная передаточная функция

Р (2 )_Х Ьт2т

H (Z ) =

G (z) X ^mZ^'

(6)

для которой анализ устойчивости существенно проще.

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

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

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

Идентификация устойчивого БИХ-фильтра

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

Пусть задана структура (порядки числителя и знаменателя) одномерной передаточной функции вида (6), которой по предположению может быть аппроксимирован двумерный БИХ-фильтр. Представив эту передаточную функцию по отрицательным степеням 2 :

H (z ) =

F (z ) = X

G (z) Xbmz~

(13)

можно записать разностное уравнение:

М М

§ (п)_-1 Ьт8 (п - т)+Х ат/(п - т). (14)

т_1 т _1

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

у _ Хс + #_ Хаа + ХьЬ +#. (15)

Здесь с _[а Ь] 2Мх1 вектор искомых параметров, составленный из двух Мх1-векторов а _[а0 аМ ] и Ь _[Ь0 ЬМ ], соответствующих

коэффициентам числителя и знаменателя передаточной функции (13).

- m

т

m=0

■ g (0) " ■ #(0) "

У = g(n) #(n)

_ g (N - M )_ #(N - M )

Xa =

g(0) g(1) ... g(m) ■

g (n) g (n) ... g (M + n)

g(N-M -1) g(N-M) ... g(N)

f (0) . f (M) " f (n) ... f (M + n)

f (N-M) ... f (N)

Пусть задан критерий качества оценок в виде функции потерь Q (с). Ставится задача по N* 1 - вектору y и N*2M - матрице X построить оценку с : Q (С) = min Q (с) (16)

Xb =

при условии

Ы < 1 + s, i = 1, M ,

(17)

где 21 - нули полинома О (2) (13) порядка М , а е -

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

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

1. Решение задачи (15) без учета условий (17).

2. Построение полинома О (2) с использованием

компонент Ь1,, _ 1,М вектора С, образующих вектор Ь.

3. Вычисление нулей 2, полинома О (2) и проверка условий (17):

|2,.| < 1 + е, , _ 1М . (18)

4. Если эти условия выполняются, задача решена, если нет - переход к следующему шагу.

5. Каждый неустойчивый полюс 2{ на комплексной плоскости Ъ заменяется ближайшим устойчивым 2* так, чтобы выполнялось условие (17).

6. С использованием набора устойчивых полюсов конструируется полином О (2), соответствующего

устойчивого фильтра с коэффициентами Ь*,, _ 1, М, образующими Мх 1 - вектор а* .

7. Решается задача оценивания вектора а, компоненты которого суть коэффициенты числителя передаточной функции (13). Соответствующее этой задаче векторно-матричное уравнение имеет вид

у - Х4Ь* _ Хаа + ц , (19)

где у - ХЬЬ - преобразованный Мх 1 - вектор выхода модели, а ц - Мх1 - вектор ошибок, учитывающих кроме прочего, ошибки аппроксимации связанные с ограничениями на устойчивость фильтра.

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

Пример реализации

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

Н (2)_ ао + а2 + а2222 , (20)

Ь0 + Ь12 + Ь2 2

т.е. с _ [а0 а1 а2 Ь0 Ь1 Ь2 ].

Для каждого положения опорной области 1-я компонента вектора у и элементы 1-й строки матрицы Х и вектора у вычислялись следующим образом:

у, _ § (0,0); Х,0 _ /(0,0);

/ (0,1) + / (1,0) + / (-1,0) + / (0,-1)

Xi,1 =-

4

X; 2 =

f (1,1)+ f (1,-1)+ f (-1,1) + f (-1,-1) .

4

X = g (0,1) + g (1,0)+ g (-1,0) + g (0,-1)

Xi,3 =

4

X;,2 =

g (1,1)+g (1,-1)+g (-1,1)+g (-1,-1)

4

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

С _ [-80058.52 -2.00 1.00 1.00 1.82 -0.82] .

Многочлен О (2), сформированный по полученным оценкам

О (2)_ 1 +1.822 - 0.8222 _ 0.

Корни соответствующего этому многочлену характеристического уравнения не удовлетворяют условию (17), поэтому строиться многочлен О (2) :

G (z )* = 1 + 1.19z-2.19z2, коэффициенты которого, составляют вектор a*.

. 14

f ■ »rm9V**n 1 '

_

Рис. 3 Искаженное изображение «Театр»

Рис. 4 Восстановленное изображение «Театр»

Рис. 5 Фрагменты для определения параметров фильтра

Затем решается уравнение (19), в результате чего получается вектор c* = [a* b* ]

с* =[178.43 -2.00 1.01 1 1.19 -2.19] (21) Фильтр, параметры которого являются компо-

* w

нентами вектора c , является устойчивым, поскольку полюса многочлена вида (20) лежат внутри единичной окружности.

На рисунке 4 показано восстановленное изображение.

Благодарности Работа выполнена при поддержке Министерства образования РФ, Администрации Самарской области и Американского фонда гражданских исследований и развития (CRDF Project SA-014-02) в рамках российско -американской программы "Фундаментальные исследования и высшее образование" («BRHE»), а также Российского фонда фундаментальных исследований (гранты №№ 03-01-00109, 04-07-90149, 04-07-96500).

Литература

1. Методы компьютерной обработки изображений / Под ред. Сойфера В.А., Москва, Физматлит, 2001.

2. Дроздов М.А, Зимин Д.И, Скуратов С.А, Попов С.Б, Фурсов В.А. Технология определения восстанавливающих фильтров и обработки больших изображений // Компьютерная оптика. № 25, 2003.

3. Зимин Д.И., Фурсов В.А., Распределенная обработка большеформатных изображений на многопроцессорных системах. // Тезисы международной конференции «Распределенные вычисления и Грид -технологии в науке и образовании», Дубна 29.07-2.08 2004 г.

4. Зимин Д.И., Фурсов В.А., Распределенная обработка большеформатных изображений на многопроцессорных системах. Тезисы третьего Международного научно-практического семинара «Высокопроизводительные параллельные вычисления на кластерных системах», Самара, 13.09-15ю09 2004 г.

5. Фурсов В.А. Идентификация моделей систем формирования изображений по малому числу наблюдений. Самара 1998.

6. Василенко Г.И., Тараторкин А.М., Восстановление изображений, М. Радио и связь, 1986.