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

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

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

Аннотация

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

где ю=ехр{2пу^г}, т=0,1,—,N-1, справедливо со-

Введение

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

Есть хорошо известные программы обработки изображений, позволяющие выполнять фильтрацию (например, Adobe Photoshop). Однако такие программы, как правило, не предназначены для обработки большеформатных изображений, и, тем более, для фильтрации изображения по мере поступления информации. Так попытка простого «открытия» в Adobe Photoshop 220-мегабайтного изображения в формате TIFF требует около 300 Мб дискового пространства и выполняется около 9 минут на компьютере Pentium-III, 650МГц, 128Мб RAM.

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

Процедура восстановления и ее оптимизация

Вычисление свертки спектральным методом

Известно, что выходная последовательность фильтра представляет в математических терминах свертку исходной последовательности x(n) и импульсной характеристики фильтра h(n). В свою очередь, известный способ эффективного вычисления свертки состоит в использовании теоремы о свертке и быстрого преобразования Фурье (FFT) [1].

Пусть y (m )=(x *h) (m ) - дискретная круговая свертка двух N-периодических функций:

N-1

У(m)= ^ x(n)h (m-n). (1)

n=0

Тогда, как известно, для спектров дискретных преобразований Фурье

N-1

У(m)=X y (n)®mn , (2)

n=0

отношение

у(т)=х(т)/г(т) . (3)

Поэтому вычисление свертки (1) с помощью (2) и (3) может быть проведено по схеме приведенной на рис. 1, где ОДПФ - обратное дискретное преобразование Фурье:

N-1

У (n) = N-1 £ y (m)( m=0

(4)

Наличие быстрых алгоритмов вычисления преобразований (2) и (4) позволяет существенно снизить вычислительные затраты.

Рис. 1. Вычисление циклической свертки через БПФ

Очевидно, что суммарная сложность такого вычисления свертки определяется сложностью ККГ и равна

Ж (N )=2ЖРРТ (N )+Х , (5)

где ЖрТЩ - сложность преобразования Фурье, X - количество дополнительных операций для вычисления произведения спектров. Здесь предполагается, что сложность вычисления прямого и обратного преобразования Фурье одинакова, и спектр /г(т) импульсной характеристики фильтра Н(п) рассчитан заранее.

Секционированное вычисление свертки

Для задач линейной цифровой фильтрации сигналов, т.е. вычисления дискретной апериодической свертки

<Х)

У(п)= X И(к)X(п-к), (6)

к=-<»

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

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

Пусть импульсная характеристика (ядро свертки) имеет длину в М отсчетов и отлична от нуля на интервале [-М, М+] (М + М+1=М Тогда пределы суммирования в (6) становятся конечными:

Мт

у (п )= ^ к (к )х (п - к). (7)

к=-М _

При вычислении одномерной апериодической свертки вида (7) использование секционирования по методу перекрытия с накоплением заключается в следующем. Из входной последовательности выделяются блоки х(1) (к) размером в N отсчетов (К>М) с перекрытием в (М-1) точках (см. рис. 2):

х(/) (к)=х [к+/ ^-М+1)] при 0<к < N-1, /=0, ±1, ±2,— .

х(п)

X \ку-А

х(к)

х +1(к)

п

-+•

N

м + • ы%;

‘ •' у1 -1(к) •' /

(------4--------=---—--------(--------Г

м

У(п)

N

у(к)

М \

N

\М V. \4м *.

4 * « «

/ +1/7 Ч

У (к) ,

Рис. 2. Схема расположения блоков при секционировании одномерной свертки методом перекрытия с накоплением

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

U (М )=•

Ж (N )

(8)

у [к+/ (-М+1)=у(/)(к) ,

М + +1< к < N - М “-1.

Таким образом, на каждые (^М+1) отсчетов выходной последовательности вычисляется одна циклическая свертка, что требует выполнения Ж(Ы) арифметических операций (5). Считая входную последовательность бесконечной, определим число операций, приходящееся на один отсчет выходного сигнала. Для импульсной характеристики длиной М оно равно:

N - М+1

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

а) выбора наиболее быстрого алгоритма ДОП длины N, используемого при вычислении циклической свертки;

б) выбора оптимальной длины секции N. Результаты такой оптимизации представлены ниже. Данный подход легко обобщается на случай

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

У(пЬп2)=

М+ М+

= X X Ь (кь к2 )х(п1-kl,п2 -к2)

к=-М] к=-М2

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

п

на всей плоскости аргументов) разбивается на равные прямоугольные блоки размерами М]ХМ2 с шагами (N1 -Ы\ +1) и (М2 -М2 +1), где М1 = + М+ +1,

М2 = М2 + М2 +1 : х(/1,12) ((1, к2 )=

=X[к1 + /1 (( -М1 +1),к2 +12 ((2 -М2 + 1)]

при 0<А’1 <N1, 0<к2 <N2, /1, /2 = 0, ±1, ±2, К

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

формируются выходные блоки у(/1,/2 )(^1, к2) таких

же размеров. Далее из средних частей выходных блоков формируется искомый результат:

У [к1 +11 (^ -М1 +1) = к2 +12 (N2 -М2 +!)] =

=У(М2 )(к1, к2)

при М1 +1< к1 < N1 - М- -1, М2 +1< к2 < N2 - М2 -1 Аналог критерия (8), т.е. число арифметических операций, приходящееся на один отсчет выходного сигнала, в данном случае задается выражением:

Ж (Щ, N2 )

^ (М1, М2 )=т--------- .у ^---------- ,

11 2’ (N1 -М1 +1)(2 -М2 +1)

где Ж(^ ,^) - число операций при вычислении двумерной циклической свертки размерами ^х^.

Оптимизация спектрального алгоритма для вычисления свертки с секционированием

Наиболее популярным и известным алгоритмом вычисления ДПФ (2) до сих пор остается алгоритм быстрого преобразования Фурье (БПФ), предложенный в 1965г. Кули и Тьюки [1] и рассматриваемый большинством пользователей для последовательностей, длина (или период) которых равна целой степени двух: N=2'. Неблагоприятным следствием популярности БПФ-алгоритма Кули-Тьюки является широкое распространение мнения о том, что применять дискретное преобразование Фурье практично лишь при такой длине последовательности. В результате БПФ-алгоритмы стали диктовать параметры применяемых устройств вместо того, чтобы приложения диктовали выбор подходящего алгоритма БПФ [1].

В задачах, не предполагающих жестких аппаратурных требований к длине обрабатываемого сигнала, применение только традиционных алгоритмов БПФ Кули-Тьюки приводит к необходимости увеличивать период обрабатываемых сигналов - добавлять нулевые отсчеты до ближайшего «хорошего» целого числа вида N*=2'. Такие числа расположены в натуральном ряду весьма редко, что приводит к почти двукратному увеличению удельной вычислительной сложности для «плохих» N, равных, например, 2Г +1=N (см. рис.3 график 1).

На самом деле быстрые алгоритмы преобразования Фурье могут быть построены для большого

набора длин. В настоящей работе авторы производят выбор оптимальной структуры алгоритма длины N на основе, следующих БПФ-алгоритмов:

1) БПФ Кули-Тьюки с декомпозицией по основанию 2 [1], [2];

2) БПФ Кули-Тьюки с декомпозицией по основаниям 3 и 6 [4];

3) редукция Гуда-Томаса, для формирования алгоритмов промежуточных длин [1], [2];

4) специальные алгоритмы ДПФ коротких длин [2];

5) использование идеи совмещения для вещественных входных данных [1], [2].

На рис. 3 представлены графики удельной арифметической сложности вычисления спектра комплексной одномерной последовательности произвольной длины. Здесь N - реальный период; Ж*(№) -удельная арифметическая сложность (число арифметических операций на один отсчет входного сигнала) вычисления БПФ сигнала с увеличенным периодом ^>М На графике 1 приведена сложность вычисления спектра традиционным способом - на базе БПФ Кули-Тьюки по основанию 2.

Рис. 3. Удельная сложность вычисления одномерного спектра

На графике 2 представлена вычислительная сложность наилучшего алгоритма, выбранного из предложенного набора. При этом сложность преобразования уменьшается в среднем в 1,6 раза. Отметим, что для ряда длин удельная сложность вычисления спектра снижается более чем в два раза.

Рис. 4. Удельная сложность вычисления одномерной свертки

Рис. 4 показывает зависимость удельной сложности вычисления секционированной свертки длинного сигнала от размера импульсной характеристики. Здесь также первый график соответствует использованию традиционного БПФ Кули-Тьюки по основанию 2, а второй график показывает удельную сложность вычисления свертки при использовании предлагаемых алгоритмов. Сочетание секциониро-

вания и эффективных БА ДПФ приводит к снижению сложности вычисления свертки в 1,1-1,4 раза.

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

На рис. 5 приведена минимальная удельная сложность вычисления спектра двумерного массива объема NxN посредством преобразования Фурье объема N*xN*, где выбрано из предлагаемого

диапазона длин. График 1 построен для случая, когда используется посторочно-столбцовый способ преобразования при N*=2', график 2 показывает возможности построчно-столбцовых алгоритмов из расширенного набора длин, график 3 построен при использовании существенно двумерных алгоритмов.

у связаны с числами а и Ь соотношениями

Рис. 5. Удельная сложность вычисления двумерного спектра

Исследования показывают, что в двумерном случае расширение диапазона длин позволяет уменьшить вычислительную сложность в среднем на 63,8% (или до 92%, то есть практически в 2 раза в отдельных случаях). Переход к существенно двумерным алгоритмам позволяет увеличить этот выигрыш в среднем до 74% ( или до 97,5% в отдельных случаях).

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

Оптимальная реализация быстрого двумерного алгоритма ДПФ 216x216

Задача оптимизации спектрального алгоритма для вычисления свертки с секционированием была решена нами для использования в конкретной системе обработки большеформатных изображений с заданными размерами выходного блока не менее 170x170 пикселов и размерами импульсной характеристики 45x19 пикселов. В этих условиях оптимальным оказалось использование БПФ 216x216 (или 63х63), алгоритм которого предлагается ниже.

Представление комплексных чисел в виде г-кодов

Пусть

у=ехр)=2 (-1+*^) , 7=2 (-1-/л/з) .

Тогда для комплексного числа г, наряду с алгебраической формой представления г=а+Ьі, возможна и форма г=у5 +ур [4], где вещественные х и

3) а • р Т Х/зГ”-

(В частности, для вещественных г справедливы равенства г= -5= -р ).

Пару вещественных чисел (5, р), ассоциированную с представлением г в виде г=у«+ур , будем

называть у-кодом числа г. Арифметические действия над комплексными числами индуцируют следующие правила действия над кодами:

(5,у) + (р, V) = (5+и,р + V) , (9)

( 5, р )-(и, V ) = (( р - 5 )( V - и ) + SV• (р - 5 )( - и )-pv )

Таким образом, сложение комплексных чисел, представленных у-кодами, реализуется с помощью двух вещественных сложений, а умножение - с помощью трех вещественных умножений и трех сложений, то есть не сложнее, чем умножение и сложение комплексных чисел в обычной алгебраической форме. Особенно просто реализуется умножение на числа у и у , имеющие соответственно коды (1,0) и (0,1). Соотношение (9) в этом случае принимает вид: (1,0)-(и, V)=(и, и - V)

(0,1)-(и, v)=(v - и, -и ) ’

и не содержит нетривиальных умножений.

Отметим, что использование такого представления данных выгодно не только для алгоритмов преобразования Фурье по основанию 3 [4], но и для БПФ-алгоритмов по основанию 6, так как корень шестой степени из единицы равен -у:

Г 2пі 1 Г пі 1

ехр 1=ехр 17 |-у.

Описание алгоритма. Пусть N=216 и

х (Щ, п2) - двумерная входная комплексная последовательность, ю=ехр (2п*/N) , X (, М2) - двумерный дискретный спектр Фурье:

N-1 N-1

(10)

х

(М1, т2 )=Х X х (п1' п2 )

ют¥*1+т2п2 (Ц)

«1=0 П2 =0

Тогда редукция преобразований (11) к ДПФ массивов объема (NхN) осуществляется по формуле:

5 ь N6-1

X(гщ,М2)= X юаМ1+ М X х(6« +а,6«2 +Ь)х

а,Ь=0 «,«2=0

5

'= X Ба,ь(т1/т2)

а,Ь=0

(12)

/ 6\ М1п1+М2п2 т)

т.е. вычисление ДПФ объема 216x216 сводится к 36 ДПФ объема 36x36 и т.д. При этом уже на втором шаге, для вычисления ДПФ объема 6x6 при представлении степеней ю и входной последовательности у-кодами уже не требуется вещественных умножений.

Еще одним фактором, упрощающим алгоритм ДПФ, является то, что функции 8а,Ь(т1,т2) достаточно вычислять для пар (т1,т2), лежащих в относительно малой "фундаментальной области":

Л={ 0<т1, м2 <N-1} .

Действительно, если т=ехр(пУз) и (щ,М)єЛ ,

то

^а, Ь (т1 +а(-, М2 +Р^)=та“+ЬР^а, Ь (МЬ т2 ) .

Реализация предложенного алгоритма ДПФ предполагает выполнение следующей последовательности шагов.

Шаг 1. Формирование массива у-представ-лений входной двумерной последовательности размерами 216*216, а также констант алгоритма (степеней ю), если последние не сформированы заранее.

Шаги 2-4. Последовательное выполнение трех шагов вычислений согласно формуле редукции (12), и в результате - формирование двумерного дискретного Фурье-спектра, представленного в у-кодах.

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

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

Вычислительная сложность и структурная характеристика алгоритма. Анализ вычислительной сложности предложенного алгоритма ДПФ показывает, что он требует выполнения около 50 арифметических операций в расчете на один комплексный пиксел преобразуемого двумерного массива 216*216 и является наиболее простым в реализации среди всех известных алгоритмов ДПФ массива (Ш№) при 150 < N< 360.

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

Кроме того, существенным достоинством алгоритма является его простая структура.

Схема декомпозиции однородна, включает всего три шага. Вычисления (шаги 2-4) могут реализовываться как в “прямой” форме (от одноточечных преобразований к 6-точечным, 62-точечным и т.д.), так и с помощью обратной рекурсии.

Экспериментальное исследование

Разработанный алгоритм ДПФ 216x216 используется теперь в программном комплексе обработки и восстановления больших изображений. Время восстановления (фильтрации, свертки) полутонового изображения 15000x15000 пикселов (около 220Мб) на компьютере Pentium-III, 650МГц, 128Мб RAM составляет 5 минут. При замене эффективного алгоритма ДПФ 216x216 традиционным алгоритмом преобразования Фурье наиболее близкого размера 256x256 время обработки изображения возрастает до 8,5 минут, то есть в 1,7 раза.

Заключение

Таким образом, в работе представлена методика оптимизации восстановления большеформатных изображений по критерию вычислительной сложности. Показаны способы формирования эффективных алгоритмов одно- и двумерных БПФ в широком диапазоне длин. Разработанные алгоритмы эффективнее традиционного БПФ в 1,5-2 раза. Для программного комплекса обработки изображений при заданных параметрах разработан алгоритм БПФ 216x216 с представлением данных в у-кодах, являющийся наиболее эффективным в заданном диапазоне длин. Экспериментально подтверждена эффективность его использования.

Работа выполнена при финансовой поддержке Российского Фонда Фундаментальных Исследований, проект № 00-01-00600.

Литература

1. R.E. Blahut Fast algorithms for digital signal processing // Addison-Wesley Publishing Company, Inc., 1985.

2. H.G. Nussbaumer Fast Fourier transform and convolution algorithms // Springer-Verlag Berlin Heidelberg, 1981.

3. B. Aramberpola, P.J.W. Rayner Multidimensional fast Fourier transform algorithm // Electron. Lett. 1979. 15(3).

4. V.M. Chernov Algorithms, realizable in Hamilton-Eisenstein codes, for two-dimensional discrete orthogonal transforms // Problems Inform. Transmission, 1995. Vol. 31. N. 3. Р. 228-235.