О рекурсивном вычислении свертки изображения и двумерного неразделимого КИХ-фильтра Текст научной статьи по специальности «Математика»

ОБРАБОТКА ИЗОБРАЖЕНИЙ: МЕТОДЫ И ПРИКЛАДНЫЕ ЗАДАЧИ

О РЕКУРСИВНОМ ВЫЧИСЛЕНИИ СВЕРТКИ ИЗОБРАЖЕНИЯ И ДВУМЕРНОГО НЕРАЗДЕЛИМОГО КИХ-ФИЛЬТРА

Мясников В. В.

Институт систем обработки изображений РАН, Самарский государственный аэрокосмический университет

Аннотация

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

Введение

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

у(п1,п2) = ^ И(т1,т2)х( -т1,п2 -т2). (1)

( ,т2 )е В

Здесь у(п1,п2) - выходное изображение, В - конечная (ограниченная) область ненулевых отсчетов импульсной характеристики. Будем считать, что область В задана следующим образом:

В = {Ц, т2 )е[0, Ых ]х[0, Ы2 ]}. (2)

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

изображения и КИХ-фильтра одним из первых указал Л.П. Ярославский в работе [3]. Определенное обобщение и развитие этого подхода было представлено в работах [4, 5, 8]. Работы [6, 7, 9] использовали полиномиальные импульсные характеристики фильтра для рекурсивной реализации свертки. Авторская работа [9] дает описание метода построения алгоритма рекурсивного вычисления свертки изображения и неразделимого двумерного полиномиального КИХ-фильтра.

Настоящая работа представляет развитие метода, предложенного в [9], для случая, когда ИХ фильтра является произвольной неразделимой двумерной функцией с конечной областью определения.

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

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

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

Четвертый раздел дает алгоритм расчета параметров двумерной рекуррентной конечной ИХ для представления с ее помощью произвольной двумерной неразделимой ИХ линейного фильтра. Для этого используется аппроксимация данного двумерного линейного фильтра посредством вертикальных и горизонтальных двумерных рекуррентных соотношений. В качестве показателя качества аппроксимации

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

В пятом разделе приведена информация о поддержке данной работы со стороны грантов и спонсоров.

1. Рекурсивное вычисление свертки одномерного сигнала и КИХ-фильтра

Данный раздел содержит обзор некоторых известных результатов рекурсивной фильтрации для одномерного случая. Детали алгоритмов рекурсивного вычисления свертки сигнала и одномерного КИХ-фильтра могут быть найдены в работах [4, 5, 8].

Операция вычисления свертки сигнала и одномерного КИХ-фильтра может быть представлена в виде:

у(п) = ^Н(т)х(п - т).

(3)

т^В

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

Пусть ИХ фильтра удовлетворяет рекуррентному соотношению К-ой степени. Если данное утверждение не выполняется и ИХ не удовлетворяет рекуррентному соотношению, в этом случае всегда можно аппроксимировать ее с помощью такого рекуррентного соотношения. Поэтому предположим далее, что для импульсной характеристики фильтра выполняется следующее соотношение:

Н(т)=

Ьт, т = 0, К -1,

К

^Оакк{т - к), т = К,М -1,

к=1

0, т < 0, т > М.

(4)

Здесь величины ак (к = 1, К) являются коэффициентами рекуррентного соотношения, а величины Ьт (т = 0, К -1) - его начальными значениями [10]. Будем в дальнейшем называть импульсные характеристики, удовлетворяющие соотношению (4), одномерными рекуррентными импульсными характеристиками. Тогда выражение вычисления свертки (3) сигнала и фильтра с одномерной рекуррентной ИХ может быть представлено в рекурсивном виде [5, 8]:

У(п) = X а*У(п - к) + X И (т) х(п - т) -

к=1 т=0

X-1

-X И (т)х(п -М - т).

т=0

Функции И + (к) и И ~(к), использованные в данном выражении, могут быть вычислены заранее. Они задаются следующими выражениями:

т

И+ (т) = И (т) -Ха кИ(т - к),

И (т) = X акИ(М + т - к), т = 0, X -1.

к = т+1

Обозначим далее

X-1

у+ (п) = X И+ (т) х(п - т),

т=0 X-1

(5)

у (п) = XИ (т)х(п-т).

Тогда выражение, определяющее процедуру рекурсивного вычисления свертки, примет вид:

У(п) = XakУ(n - к) + у + (п)- у (п - М).

(6)

к=1

Алгоритм рекурсивного вычисления свертки, использующий выражения (5) - (6), требует и*(К) = 3К умножений и и+ (К) = 3К -1 сложений на один пиксель выходного сигнала. Будем в дальнейшем использовать для них следующие обозначения:

и+ (К) = 3К -1, и*(К) = 3К . (7)

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

Например, одной из наиболее простых рекуррентных импульсных характеристик является прямоугольная ИХ:

\Ъ, т е[0,М -1],

И (п ) =

0, т й[0,М-1].

(8)

Она удовлетворяет рекуррентному соотношению первого порядка (К = 1). Параметры рекуррентного соотношения для этой ИХ следующие: а1 = 1, Ъ0 = Ъ . Подставляя эти величины в выражение (5)-(6) можно легко получит следующую рекурсивную процедуру для вычисления свертки одномерного сигнала и фильтра с прямоугольным импульсным откликом:

у(п) = ау(п -1) + И(0)х(п) - И(М - 1)х(п -М).

В соответствии с выражением (7) эта процедура требует два сложения и три умножения на один отсчет выходного сигнала. Однако, принимая во внимание, что а1 = 1, И(т) = Ъ , окончательное выра-

жение для рекурсивного вычисления свертки принимает хорошо известное выражение:

у(п) = у(п -1) + Ъ(х(п)- х(п - М)). (9)

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

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

И(т1,т2) =

т1 = 0, К1 -1, т2

К

XatИ (т1 - к ,т2),

0, К2 -1,

т1

= К,,М, -1, т, = 0, К2-1,

(10)

"2

Xa,И(т1,тг -к),

--К 2, М 2 -1,

т = 0, М1 -1, т2 0, т ( [0,М1 -1] ог т2 ( [0,М2 -1].

Здесь величины а\ (к = 1, К1) задают коэффициенты вертикального рекуррентного соотношения, величины а^ (к = 1, К2) задают коэффициенты горизонтального рекуррентного соотношения, а величины Ът т , (т1 = 0, К1 - 1, т2 = 0, К2 -1) являются

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

И(т1 ,т2) =

_ 2

XагИ (т1,т2 - к),

т1 = 0, К1 -1, т2 = К2, М2 -1,

К1

XаkИ (т1 - к ,т2),

(11)

т1 = К1,М1 -1, т2 = 0, М2 -1,

к =0

т=0

к=1

к=1

к =1

к =1

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

Утверждение. Значения двумерной рекуррентной ИХ зависят только от начальных значений и коэффициентов вертикального и горизонтального рекуррентных соотношений. Значения двумерной рекуррентной ИХ не зависят от формы представления двумерного рекуррентного соотношения.

Доказательство приведенного утверждения заключается в очевидной проверке.

Для синтеза рекурсивного алгоритма вычисления свертки используется выражение (10). Для этого подставим горизонтальное рекуррентное соотношение, представленное в третьей строке выражения (10), в выражение для двумерной свертки (1). Это дает следующее полу-рекурсивное выражение, которое похоже на одномерное рекурсивное выражение (6):

— 2

У(«1,«2) = ^4akУ(п1,«2 -k) +

k=1

(12)

+ y+(«1> П2 )-У 2 («1 ' n2 - M2 )•

Данное выражение определяет основное соотношение в процедуре рекурсивного вычисления двумерной свертки по ее предшествующим значениям. К сожалению, величины у+(п1,п2),у-(п1,п2) определяются как двумерные свертки:

M1 -1 K2 -1

y+(«1, «2 )= X Х^+Ц, т2 )х(«1 - Щ, «2 - ^2),

^1=0 т2=0 -1 K2-1

(13)

У 21 (п1, п2 )= X X ^2 (т1, т2 )х(п1 - m1, п2 - т2).

Ш1 =0 т2=0

Легко проверить, что импульсные характеристики й+(т1, т2) и Й2"(т1, т2), используемые в этих выражениях, задаются следующими выражениями и могут быть рассчитаны заранее:

(т2 = 0, К2 -1, т1 = 0,М1 -1),

2

h+(m1;m2 ) = h(m1,т2 )-XakÄ(m1; m2 - k)

k=0

(14)

2

h-(m1,m2)= ^a^A(m1;M2 + m2 -k).

k=m2 +1

Окончательно, рекурсивное вычисление двумерной свертки (1) с использованием выражений (12) и (13) требует следующее количество арифметических операций:

и+ (,М2,К1,К2) = 2К2 ((М1 + 0,5)-1,

и* (М1,М2, К1, К2) = 2К2 (М1 + 0,5) . (15)

В полученных выражениях для вычислительной сложности рекурсивная часть (12) процедуры требует всего К2 +1 операций сложения и К2 операций умножения. А основная вычислительная сложность относится к вычислению значений двумерных свер-

ток (13) изображения и двумерных КИХ-фильтров h+(m1, m2 ) и h-(m1, m2 ), рассчитываемых на границах исходной импульсной характеристики. Следовательно, снижение вычислительной сложности этих операций приведет к снижению сложности рекурсивного алгоритма вычисления свертки в целом.

3. Рекурсивная процедура вычисления свертки изображения и двумерного неразделимого КИХ-фильтра

Утверждение. Если ИХ фильтр h(m1, т2 ) удовлетворяет вертикальному рекуррентному соотношению, заданному второй строкой выражения (10), тогда импульсные характеристики h+(m1, m2 ) и h-(m1, m2 ) также удовлетворяют этому рекуррентному соотношению.

Доказательство приведенного утверждения заключается в очевидной проверке.

С учетом приведенного утверждения следующие рекуррентные соотношения справедливы для импульсных характеристик h+(m1, m2 ) и h-(m1, m2 ):

(m1 ,m2 ) =

¿w m = 0, K1 -1, m2 = 0, K2 -1,

K1

Xakh# (m1 -k,m2 ),

(16)

m1 = K1,M1 -1, m2 = 0, K2 -1, 0, m1 ^ [0,M1 -1] or m2 £ [0, K2 -1] .

Здесь знак '#' обозначает либо знак «+», либо знак «-». Подставляя выражение (16) в выражение (13), получаем следующие два выражения для рекурсивного вычисления сверток (13):

(17)

У2+ ^ П2 ) =ХакУ2+ (п1 - k, п2) +

к =1

+ У2+1+ (П1, П2 ) - У2+Г (П1 -МР П2 ) ,

К,

У- (nl, п2) = X акУ- (п1 - k, п2) +

к=1

+ У2-1+ (nl, п2 ) - У-1- (п1 - M1, п2 ) .

Здесь величины У## (п1, п2) вычисляются как двумерные свертки, задаваемые следующими выражениями:

У## (П, П2) = (П1, П2) **х П П2) = К1 -1 К2-1 (18) = Х X т2 )х (п1 - ml, п2 - т2 ).

т1 =0 т2 = 0

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

Алгоритм

Шаг 1. Прямое вычисление четырех сверток (18): У## (п„ П2 ) = Й2#1# (п„ П2 ) * *х(п„ П2 ) .

k=1

ТТТяг 2 Рекурсивное вычисление значений У+ (ъ п2), У-(1, п2), используя выражение (17). Шаг 3. Рекурсивное вычисление очередного отсчета выходного изображения, используя выражение (12).

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

U + ((„ K 2 ) U * ( K 2 )

Шаг 1 4((K2 -1) 4KjK 2

Шаг 2 2(Kj + 2) 2Kj

Шаг 3 K 2 + 2 K 2

Итого 4KjK2 + 2Kj + K2 + 2 4KjK2 + 2Kj + K2

Таким образом, суммарное число всех арифметических операций (и + + и *) пропорционально величине 8К1К2 и не зависит от размеров области определения ИХ фильтра.

4. Представление двумерного неразделимого КИХ-фильтра в виде двумерного реккурентного соотношения Произвольная двумерная неразделимая нерекурсивная импульсная характеристика линейного фильтра может быть аппроксимирована двумерной рекурсивной ИХ. Получение коэффициентов вертикальных и горизонтальных рекуррентных соотношений связано с решением оптимизационной задачи. В качестве критерия задачи оптимизации может быть выбран критерий минимума среднеквадрати-ческой ошибки аппроксимации ИХ с помощью рекуррентной ИХ. К сожалению, такая задача оказывается существенно нелинейной относительно коэффициентов рекуррентных соотношений и не может быть решена в явном виде [5, 11]. Поэтому используем подход, предлагаемый в работе [11]. А именно, заменим критерий задачи оптимизации на другой, который позволяет записать решение задачи в явном виде. То есть, выберем в качестве показателя критерия величину линейного предсказания очередного отсчета ИХ посредством ее предшествующих отсчетов. Порядок предсказания должен соответствовать порядку рекуррентного соотношения. Тогда для вычисления коэффициентов рекуррентных соотношений в выражении (10) необходимо решить следующие две независимые задачи оптимизации:

В соответствие с результатами, представленными в работе [11] решение, получаемое для из-

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

Утверждение. Решения задач (19) дают оптимальное решение для первоначальной задачи оптимизации, если ИХ является рекуррентной и ее порядок равен порядку предсказания.

Утверждение. Решения задач (19) являются асимптотически оптимальным, то есть чем больше линейный размер ИХ, тем ближе получаемое решение к искомому решению первоначальной задачи оптимизации.

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

а1 = [я* ]-1 .

Нотации, использованные в этом выражении, обозначают следующие вектора и матрицы:

( 1 А а

а =

Bu =

b11

b

и

1Ki

и

b

и

bU =

f A b01

bE1

v boKi J

Здесь величины bfk вычисляются в соответствии с выражением:

M 2-1 M1 -1

bS = Z Zh(m1 - S'm2 )(m1 - k'm2 )•

7*2 —и Ш\— Л-J

M2 Mj-1 f Kj 2

=Z Z h(m1 ,m2 ЫаЖ - k,m2) ^ min, 4

m2 = 0 mj =Kj v k=1 J (19)

Mx M2-1 f K2 A 2

=Z Z h(m1 ,m2 Ы^Ц ,m2 - k) ^ min. 5

m!=0 m2 =K2 v k=1 J

Явное аналитическое выражение для ошибки аппроксимации также может быть легко получено.

Благодарности

Работа выполнена при поддержке российско-американской программы "Фундаментальные исследования и высшее образование" (BRHE) и гранта Президента РФ №НШ-1007.2003.01 .

Литература

1. Blahut R.E. Fast Algorithms for Digital Signal Processing // Reading, Wokingham: Addison-Wesley, 1985.

2. Dudgeon D.E., Mersereau R.M. Multidimensional Digital Signal Processing // Prentice-Hall, Inc., Cliffs, 1984.

3. Ярославский Л .П. О возможности параллельной и рекурсивной организации цифровых фильтров // Радиотехника, No. 3, 1984, стр. 87-91.

Сергеев В.В., Параллельно-рекурсивные КИХ-фильтры для обработки изображений // Компьютерная оптика, 1992. No.10-11, с.186-201. Chernov A.V. Fast Recursive Computation 1D and 2D Finite Convolution // Proceedings of 7th International Conference on Pattern Recognition and Image Analysis: New Information Technologies, St.Peterburg, Russia, 2004. P.1001-1004.

K, 1

6. Glumov N.I. , Myasnikov V.V. , Sergeyev V.V. Application of polynomial bases for image processing using sliding window // SPIE, Image Processing and Computer Optics, 1994. Vol.2363, P.40-49.

7. Glumov N.I. , Myasnikov V.V. , Sergeyev V.V. Parallel-Recursive Local Image Processing and Polynomial Bases // Proceedings of the Third IEEE International Conference on Electronics, Circuits, and Systems ICECS'96, Rodos, Greece, 1996. P.696-699

8. Myasnikov V.V. Methods for Designing Recursive FIR Filters // 7-th International Conference on Computer Vision and Graphics (ICCVG-2004), Warsaw, Poland, Springer, 2004. P.845-850.

9. Myasnikov V.V. Recursive algorithm of calculation the convolution of image and inseparable 2-D polynomial FIR-filter // Proc. of 7th Int. Conf. on Pattern Recognition and Image Analysis: New Information Technologies, St.Peterburg, Russia, 2004. P.327-330.

10. Lidl R. and Niederreiter H. Finite Fields // 2-nd Edition, Cambridge University Press, Cambridge, 1997.

11. Myasnikov V.V. On the solution of the recurrent equation used for the FIR-filter implementation // The IASTED International Conference on Signal and Image Processing (ACIT-SIP 2005), Novosibirsk, Russia, June 20-24, 2005 (printing).

12. Ifeachor E.C., Jervis B.W. Digital Signal Processing: A Practical Approach // Prentice-Hall, Inc., Cliffs, 2002.