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

Вычислительные технологии

Том 11, Специальный выпуск, 2006

ОПТИМИЗАЦИЯ НЕКОТОРЫХ АЛГОРИТМОВ ОБРАБОТКИ ИЗОБРАЖЕНИЙ С ИСПОЛЬЗОВАНИЕМ УЗЛОВ МИНИМАЛЬНЫХ КУБАТУРНЫХ ФОРМУЛ*

В. Б. Кашкин, О. В. Колямкин, М.В. Носков Красноярский государственный технический университет, Россия

e-mail: noskov@fivt.krasn.ru

New algorithms for the image processing (in particular, the algorithm of inverse filtration) that use nodes of the minimal cubature formulae for trigonometric polynomials are described in the paper.

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

изображения — N2 и числа отсчетов (узлов кубатурной формулы) — N2/2 [2, 3],

N2

тов для стандартного ДПФ на два подмножества, каждое из которых является набором узлов минимальной кубатурной формулы точности N — 1 для двумерного случая, а именно формулой с центрированной решеткой (такое разбиение носит в данной статье название "двойной центрированной решетки"). Для двойной центрированной решетки в статье предложены алгоритмы некоторых задач обработки изображений, в том числе алгоритмы фильтрации и инверсной фильтрации. Применение этих алгоритмов дает выигрыш по времени обработки изображений примерно на 15-20% за счет уменьшения числа операций, Конечно, в данном случае ДПФ применяется не ко всему растровому изображению, а по отдельности — к подмножествам пикселов, на которые разбивается все их множество,

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-01-00823) и Министерства образования и науки РПН (грант № 2.1.2.559 н).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

1. Применение решетчатых кубатурных формул

к двумерному дискретному преобразованию Фурье

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

/(х,у) = £аа,в£ |ав,в| < ~

— это ее ряд Фурье, При дискретном преобразовании Фурье устанавливается взаимно-

/(х, у)

(х^у^), 1 ^ з ^ N и множеством чисел Ааа,ра, 1 ^ 5 ^ N, такое, что тригонометрический полином

n

Т(х,у) = £ Л^ае2™(аах+вау)

5 = 1

/(х, у)

Т(х{'),у{')) = /(х(Л, у(Л), 1 < з < N, (1)

при этом число 1 ^ 5 ^ Nесть приближенное значение коэффициента Фурье:

Г 1 n

а«.,А = / Пх, ¿х<1у « - £ >+/^>) = ^^

N

[0,1)

2 3 = 1

Рассмотрим кубатурную формулу вида

1 1 N

I[f] = j J f (x,y) dxdy « £ Cjf (x(j),y(j)) = Q[f]. (2)

0 0 j=1

Тригонометрическим мономом будем называть функцию вида

(x,y) = e2ni(ax+ey),

число |а| + |вмонома фа,в(x, y); будем говорить, что кубатурная формула (2)

d

I [Фа,в] = Q [Фа,в] (3)

для всех мономов степени не выше d и существует хотя бы один моном степени d + 1, для которого равенство (3) не выполняется. Если d = 2m + 1, то минимальное число узлов кубатурной формулы (2) тригонометрической степени d равно N = 2(m + 1)2 [5],

Кубатурные формулы тригонометрической степени d = 2m + 1 с таким числом узлов называются минимальными. Более подробная информация о минимальных кубатурных

формулах дана в [4]. В настоящей статье мы будем использовать минимальную кубатур-ную формулу

[0,1)2

/ (х, у) dxdy ~

№

■§1 §2 \ ¡81 1 8 2 1

у / +/ , —+

2(к +1)2 ^ ' \к + 1' к + 1/ ' \к +1 2(к + IV к +1 2(к + 1)

Конфигурация узлов этой формулы называется центрированнной решеткой. Легко видеть, что сдвиг множества узлов на 1/2(к + 1) вдоль любой из осей дает с учетом периодичности функции такое множество точек, что в объединении с множеством узлов формулы (4) мы получим все множество точек вида

г,] =0,1, 2,..., к,

использующееся для ДПФ в двумерном случае с числом отсчетов 4(к + 1)2,

1

2. Преобразование центрированной решетки

Возьмем обычную прямоугольную решетку с шагом по каждой из осей, вдвое большим, чем в исходном изображении. Такая решетка имеет узлов в четыре раза меньше, чем пикселов в изображении. Наложим на нее такую же решетку, но сдвинутую вниз и вправо на один пиксел. Новая решетка, являющаяся объединением двух предыдущих, содержит уже половину пикселов изображения, что дает достаточную точность, поскольку узлы располагаются на изображении равномерно, подобно одноцветным клеткам на шахматной доске. При этом узлы обладают следующим свойством: если одна из координат узла (в пикселах) является четным числом, то и вторая — четное число, аналогично и в нечетном случае. Это позволяет нам перейти от двумерной матрицы к двум векторам:

= ШМ+1) + 1 2(к+1-(к\М)М) + 1\ 1к+1 1 \ 2М ' 2М ) '

где 0 ^ к ^ М2 — 1, изображение имеет разрешение 2М х 2М, а символ обозначает целочисленное деление.

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

Т(к\М + 1, к +1 — (к\М)М) = /к+1

1 V 2М ' 2М ) 1 ' Тн(к\М + 1, к +1 — (к\М)М) = ,

где 0 < к < И2 - 1,0 < з,1 < И.

Ясно, что половина пикселов при обратном преобразовании не заполняется. Этот недостаток можно частично устранить таким же образом, как это предлагалось в [3]. Мы не будем здесь рассматривать эти приемы отдельно, так как наша цель — восстановление изображения и его обработка, в частности инверсная фильтрация, без "пустых" пикселов. Отметим, что в случае выбора узлов из центрированной решетки фильтрация осуществляется следующим образом: заранее известная нам передаточная функция фильтра К (и; у) разделяется на два одномерных вектора по указанным выше формулам, т, е.

г^нечет ту Кк+1 = К

/ 2(к\М + 1) + 1 2(к + 1 — (к\М)М) + 1\ V ш ' 2М )

где 0 < к < И2 - 1.

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

3. Использование двойной центрированной решетки в обработке изображений

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

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

Значения функции /(х, у) разделяются по двум совокупностям узлов — /1 (к) и /2(к) со спектрами ¿1(и) и ¿2(и). Передаточная функция инверсного фильтра имеет следующий вид:

НЕ (и) = Н (и)-1,

где Н(и) — передаточная функция при расфокусировке, которая считается известной.

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

*1>в = ¿1 (и)Н (и) *2>в = ¿2 (и)Н (и)

■ВД Н{иУ

(У) Н(иУ

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

=

*2,в =

^(ц)Р(ц,цс)

Н(и) Р2(и)Р(и,ис) Н(и)

где ис — граничная частота; Р(и, ис) — характеристическая функция интервала радиуса ис с центром в и.

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

Н (и) =

1

а/1 + а\и\

Для сглаживания изображения использовалась передаточная функция следующего вида:

Н (и) =

а|и|

у/ттщ

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

Нв(и) = л/1 + а\и\ Р{и, глс), и|

/ х \/1 + Ъ и\ ,

а| и|

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

Результаты применения предложенных алгоритмов можно видеть на рис. 1-3. На рис. 1,а представлено исходное изображение, а на рис. 2 — действительная и мнимая части спектра этого изображения, вычисленные по двойной центрированной решетке.

Рис. 1.

Рис. 2.

Рис. 3.

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

4. Оценка выигрыша, получаемого при использовании алгоритма с выбором узлов по двойной центрированной решетке

Сравним число операций при вычислении спектра классическим алгоритмом с применением быстрого преобразования по Кули — Тьюки и при вычислении спектра алгоритмом, основанным на использовании двойной центрированной решетки, также с применением быстрого преобразования по Кули — Тьюки. В таблице I — классический алгоритм. II — алгоритм с выбором узлов по двойной центрированной решетке; при этом дополнительно затрачивается 524 288 операций (преимущественно сложений) на заполнение решетки. Обычно в расчете числа операций, требуемых алгоритмом, учитывают только число

Количество операций, затраченных на вычисление спектра изображения

Алгоритм Количество умножений Общее количество операций

I 6 024192 14416 896

II 4 849 644 12 255 256

умножений, так как умножение является гораздо более трудоемким с точки зрения машинных ресурсов действием, нежели сложение. Деление рассматривается как один из случаев умножения, вычитание — сложения.

Классический алгоритм. Число операций для вычисления спектра одной строки или столбца изображения размерностью N х N элементов с применением быстрого преобразования по Кули — Тьюки равно N log2 N. Проводя преобразование по каждой из строк изображения, получаем N2 log2 N операций. Повторяя этот процесс по каждому из столбцов, мы выполняем еще N2 log2 N операций, получая в итоге 2N2 log2 N. Таково общее число умножений при вычислении спектра классическим алгоритмом с применением быстрого преобразования по Кули — Тьюки,

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

iV2

Т

log2

iV2

Т

iV2

Т

(log2 N - 1)

(5)

операций. Проводя четыре прохода, чтобы обработать все узлы изображения, получаем итоговое количество умножений — 2^2(^2 N — 1),

Таким образом, отношение числа операций, требуемых алгоритмом, основанным на использовании двойной центрированной решетки, к числу операций, требуемых классическим алгоритмом, равно

2iV2(log2 N — I) 2N2 log2 N

1-

1

log2 N

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

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

Список литературы

[1] Кашкин В.Б., Носков М.В., Осипов H.H. Вариант дискретного преобразования Фурье с узлами на параллелепипедальных сетках // Журн. вычисл. математики и мат. физики. 2001. Т. 41, № 3. С. 355-359.

[2] Носков М.В., Кашкин В.В., Осипов H.H. и др. Численные эксперименты по применению кубатурных формул высокой тригонометрической точности к дискретным преобразованиям Фурье // Кубатурные формулы и их приложения: Матер. VI Междунар. семинара-совещания. Уфа: ИМВЦ УфНЦ РАН, БГПУ, 2001. С. 76-79.

[3] Kashkin V.B., Noskov M.V., Osipov N.N. Application of latticed cubature formulas to 2D discrete Fourier transform // Pattern Recognition and Image Analysis. 2002. Vol. 12, N 4. P. 397-399.

[4] Кашкин В.Б., Носков M.B., Осипов H.H. Применение минимальных кубатурных формул к двумерному дискретному преобразованию Фурье // Кубатурные формулы и их приложения: Матер. V Междунар. семинара-совещания. Красноярск: КГТУ, 2000. С. 84-89.

[5] Носков М.В., Шмид X. Кубатурные формулы высокой тригонометрической точности // Журн. вычисл. математики и мат. физики 2004. Т. 44, № 5. С. 729-738.

Поступила в редакцию 15 сентября 2006 г.