CC BY
69
ЯЛГЙОяО/
Уфа : УГАТУ. 2011______________________________ __________________________________________________Т. 15, №5(45). С. 137-141
МАТЕМАТИЧЕСКОЕ И ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ...
УДК 004.7+621.372.542
А. В. Никоноров, В. А. Фурсов, П. Ю. Якимов
МАССИВНО-МНОГОПОТОЧНАЯ РЕАЛИЗАЦИЯ ДВУМЕРНЫХ БИХ-ФИЛЬТРОВ
Рассматривается информационная технология реализации двумерного фильтра с бесконечной импульсной характеристикой (БИХ-фильтра). Для обеспечения физической реализуемости БИХ-фильтр строится в виде параллельного соединения 4 фильтров одного квадранта. Параметры каждого фильтра одного квадранта определяются путем идентификации по тестовым фрагментам, формируемым из исходного искаженного изображения. После определения параметров фильтра осуществляется анализ его устойчивости. Для обработки крупноформатного изображения с использованием построенных указанным способом четырех фильтров одного квадранта предлагается схема эффективной массивно-многопоточной реализации БИХ-фильтра в системах на основе графических процессоров (GPU). Обработка изображений; восстановление изображений; БИХ-фипьтры; устойчивость БИХ-фипьтров; массивно-многопоточное программирование', CUDA, GPU
В работе [1] рассматривалась информационная технология определения характеристик и обработки изображений с использованием двумерных фильтров с бесконечной импульсной характеристикой (БИХ-фильтров). В частности, для решения задачи идентификации параметров фильтра в указанной работе предложено использовать малые тестовые фрагменты, которые формируются из искаженного изображения с использованием априорной информации
о геометрической форме регистрируемых объектов. При этом двумерный БИХ-фильтр строится в виде параллельного соединения физически реализуемых фильтров с опорной областью в виде квадранта, обеспечивающих компенсацию сильных искажений с использованием опорной области небольших размеров.
Вопросы повышения производительности обработки одномерными БИХ-фильтрами за счет использования параллельных вычислительных ресурсов достаточно хорошо исследованы. В работе [2] рассматривалась реализация одномерных БИХ-фильтров в многопоточных вычислительных средах. В работе [3] приведены результаты СРи-реализации БИХ фильтра. В ряде случаев возможна декомпозиция обработки изображений, в результате которой задачи сводятся к последовательному применению одномерных фильтров [4]. Однако в общем случае такую декомпозицию провести нельзя.
В настоящей работе рассматривается достаточно универсальная массивно-многопоточная реализация двумерных БИХ-фильтров, с использованием графических процессоров и технологии С1ГОА, обеспечивающая значительное
Контактная информация: adminfiUncdk.com Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии-2011»
увеличение производительности по сравнению не только с СРи-процессорами, но и с обычными С1ГОА-процедурами, использующими неоп-тимизированные коды.
1. ТЕХНОЛОГИЯ ФОРМИРОВАНИЯ И ИДЕНТИФИКАЦИИ ДВУМЕРНОГО БИХ-ФИЛЬТРА
В общем случае выражение, определяющее способ вычисления ВЫХОДНОГО отсчета у(Пх, л2) двумерного БИХ-фильтра, имеет вид [5]:
у{п
Г1 гг
“X (1)
(*!*2 * }\П2 )
Нетрудно заметить, что в данном случае требование рекурсивной вычислимости не выполняется (вычисляемый выходной отсчет на рис. 1, а обозначен кружком в центре опорной области). Для физической реализуемости двумерного БИХ-фильтра он представляется в виде параллельного соединения БИХ-фильтров одного квадранта [6] (рис. 1, 6-0).
а б в
• • о • • л
г • д ,
. . . О * *
Рис. 1. Исходная маска 5x5 (а) и соответствующие ей маски и допустимые направления рекурсии (б-д)
Для заданной опорной маски в виде квадранта и некоторой пары тестовых фрагментов, один из которых сформирован путем компьютерного ретуширования, записывается система линейных уравнений [1]:
У = Эф + % , (2)
где, в соответствии с (1),
>1 (« ,п2) *1 (г1 ,Г2 ) • • Уі (к1 ,к2) •
У = У2 («1 ,П2 ) , э = *2 (г1 ,Г2 ) • • У2 (к1 ,к2 ) •
_УN («1 ,П2 )_ 1 X N • УN (к1 ,к2) •_
% =
Х1 («1 -«2 ) Х2 («1 -«2 )
Х N («1- «2 )
N - число отсчетов на тестовых фрагментах, по которым сформирована система (2), ф = [а(п1 -г1, п2- г2),...Ь(п1 - г1, п2 - г2),...]г - искомый вектор параметров фильтра.
Для обеспечения возможности идентификации фильтра высокого порядка в работе [1] матрицу S в (5) предложено составлять из блочных матриц, каждая из которых формируется по простым фрагментам, на каждом из которых функция яркости изменяется в одном направлении, но эти направления для разных фрагментов различны. Примеры таких фрагментов приведены в работе [1].
Для вычисления параметров фильтра по данным системы (2) может использоваться простейшая оценка метода наименьших квадратов (МНК):
ф = [эт э]"18 т У,
(3)
либо более устойчивая к грубым ошибкам типа сбоев оценка метода наименьших модулей (МНМ-оценка).
2. АНАЛИЗ УСТОЙЧИВОСТИ ФИЛЬТРА
Анализ устойчивости фильтра необходим вследствие того, что получающиеся в результате идентификации фильтры одного квадранта могут оказаться неустойчивыми, что может привести к искажениям обработанного изображения [6]. Поэтому в качестве промежуточного этапа технологии предлагается осуществлять анализ устойчивости полученного фильтра.
Передаточная функция БИХ-фильтра, который мы идентифицируем, выглядит следующим образом:
н* (^ *2) = -^
N1 N2
ЕЕ а ( ^ Г2 ) г
г =0 г =0
N1 N2
(4)
ЕЕь (к к2) 41 к2
к =0 к2 =0
где N1, N2 - это размеры входной маски (будем считать, что выходная маска имеет такие же размеры, что и входная).
Нам необходимо исследовать поведение функции, находящейся в знаменателе. Для удобства обозначим:
N1 N2
В ( *1, *2 ) = ЕЕ Ь ( k1, к2 ) 2\к
к =0 к =0
В качестве основы для проведения анализа устойчивости целесообразно опираться на теорему Шэнкса [7]. В соответствии с этой теоремой для каждой точки Ь = \z2\ единичной окружности \z2\ = 1 на комплексной плоскости z2 для всех \г1\ > 1 должно выполняться условие:
В{2Х, Ь) ф 0,
а для каждой точки а = \^1 \ единичной окружности \^\ = 1 на комплексной плоскости ^ для всех \г2\ > 1 должно выполняться условие:
В(а, z2) Ф 0.
Для проверки выполнения первого из указанных условий строится годограф корней характеристического уравнения
В(^, Ь) = 0
на плоскости ^ при условии, что параметр Ь = = \z2\ «пробегает» все точки на единичной окружности в плоскости z2. Если этот годограф пересекает единичную окружность на плоскости ^1, это свидетельствует о неустойчивости фильтра. Устойчивость возможна только в случае, если пересечения годографов с единичной окружностью отсутствуют и все они оказываются внутри круга единичного радиуса плоскости .
Аналогично годограф корней характеристического уравнения
В(а, z2) = 0
на плоскости z2, образованный при условии, что параметр а = \^\ «пробегает» все точки на единичной окружности в плоскости ^1, должен находиться внутри круга единичного радиуса на плоскости z2.
При проведении анализа устойчивости следует иметь в виду, что числа а и Ь в общем случае являются комплексными, поэтому при вычислении корней следует приравнивать нулю как действительную, так и мнимую часть полинома в левой части характеристического уравнения.
На рис. 2, б и 2, в показано поведение годографа устойчивого фильтра первого квадранта, идентифицированного с использованием миры, показанной на рис. 2, а.
б
Рис. 2. Поведение годографа устойчивого фильтра первого квадранта: а - мира, использованная для идентификации параметров БИХ-фильтра; б - годограф корней характеристического уравнения B(z1, Ь) = 0; в - годограф корней характеристического уравнения B(a, z2) = 0
Из рис. 2 можно заключить, что исследуемый БИХ-фильтр, полученный идентификацией из тестовых фрагментов, является устойчивым. На рисунке хорошо видно, как годографы корней уравнения B(z1, Z2) хорошо вписываются в единичную окружность.
3. МАССИВНО-МНОГОПОТОЧНАЯ РЕАЛИЗАЦИЯ БИХ-ФИЛЬТРА В СРи СИСТЕМЕ
Поквадрантная реализация БИХ-фильтра [1] позволяет естественным образом провести распределение вычислений на четыре потока. В настоящей работе предлагается массивномногопоточных алгоритм реализации БИХ-фильтра. Такая реализация в первую очередь ориентирована на использование в ОРОРИ системах, однако также может быть использована в современных многоядерных СРи системах.
В работе [9] была предложена процедура эффективной реализации рекуррентной обработки изображения локальным окном. С небольшими изменениями эта процедура может
быть использована для решения задачи двумерной КИХ фильтрации.
В настоящей работе предлагается эффективная процедура реализации квадрантного БИХ-фильтра. БИХ-фильтр можно представить в виде последовательного соединения двух фильтров-компонент. Первый КИХ компонент зависит только от входного сигнала и соответствует числителю формулы (4). Второй зависит только от выхода и соответствует знаменателю в формуле (4). Будем называть его рекурсивным компонентом фильтра.
Проанализируем реализацию рекурсивного компонента. Будем рассматривать один квадрант, (х>0, у>0), для остальных трех квадрантов процедура выполняется аналогично. Передаточная функция для такого компонента имеет вид:
1
I
i=0,j=0
b0,0 1.
Схема в терминах
вычислений описывается далее потоков (thread), блоков (block) и сетки (grid). Эти термины стандартны для архитектур CUDA и OpenCL [S, 10]. Предлагаемая схема содержит два уровня параллелизма - на уровне блоков и на уровне потоков.
Каждый блок выполняет расчет для одного столбца данных. Для единообразия полагаем, что номер блока равен номеру столбца данных, т. е. блок с индексом g должен N раз выполнить следующее суммирование:
yk,.
I Mk-,,g-j, k=0,1,...N.
(З)
i=0, j=0
Так как вычисления выполняются рекурсивно, блоку для вычисления значения в некоторой точке (к, і) необходимо, чтобы блоки с номерами меньше j выполнили вычисления для не менее чем і строк. Таким образом, если все блоки начинают вычисления одновременно, блок с индексом Ь должен дождаться завершения вычислений предыдущими блоками.
Потоки каждого блока выполняют суммирование (5) параллельно, по схеме каскадного суммирования или редукции [8, 10]. В таком суммировании должно быть задействовано
22]1оє? т[ ~
2 потоков и тогда количество операций выполняемых потоками параллельно составит
°е, = 2]1о§2 ^, (6)
где под ]т[ понимается наименьшее целое, большее т.
Тогда количество вычислений для каждого блока составит:
а
в
Оь, = 2 N ]1о§2 т[. (7)
Так как количество (6) для всех блоков одинаковое и блок с индексом g должен ожидать завершения вычисления для ^-1) блоков, то время ожидания для блока составит примерно
0ЬК = 2(g - 1)]1о§2 т[. (8)
Если предположить, что N блоков начинают вычисления одновременно, то для последнего блока, с учетом (7) и (8), количество операций составит
О = (4N - 2)]1о§2 т[. (9)
Количество операций, требуемое для последовательной реализации алгоритма, составляет Ох = ^т2. Таким образом, теоретическая оценка ускорения параллельной реализации позволяет предположить значительное ускорение. Однако модель многопоточности современных ОРИ накладывает некоторые дополнительные условия на реализацию предлагаемого алгоритма.
В общем случае, если количество блоков равно О, то алгоритм выполняется полосами ]МО[ раз, а суммарная сложность должна составить:
О =^ / О[(4N - 2)] 1о§2 т[.
Зависимость (7) имеет квадратичный характер по N, что должно сказываться на производительности при О << N.
Возможность увеличения количества блоков О зависит от возможностей ОРИ. В настоящей работе проведено исследование производительности при увеличении О для видеоадаптера ОБ 9600. Время обработки для одной полосы и для всего изображения для N = 512, т = 5 и числе потоков 16 для различных значений О приведено на рис. 3.
□ Одна полоса Н Все изображение
160
140
120
о 2 100
80
60
40
20
0
К
К
К
й
й
16 32 64 128 256 512
Рис. 3. Зависимость времени обработки от количества блоков
Увеличение количества блоков приводит к увеличению накладных расходов планировщика потоков. Однако планировщик ОРИ ис-
ключает тупиковые ситуации, и увеличение числа блоков позволяет увеличить производительность.
4. ПРИМЕР ОБРАБОТКИ
Экспериментальные исследования качества предложенной технологии проводились в следующей последовательности:
• автоматизированный поиск тестовых фрагментов;
• ретушь найденных тестовых фрагментов;
• идентификация параметров БИХ-фильтров по тестовым фрагментам;
• проверка на устойчивость;
• обработка изображений с помощью полученных квадрантных БИХ-фильтров.
В качестве тестового изображения использована тестовая мира. В результате все идентифицированные фильтры получились устойчивыми, что позволило получить значительное улучшение обработанных изображений.
На рис. 4 приведен пример обработки изображений с использованием описанной технологии.
а б
Рис. 4. Результат обработки изображения миры: а - размытое изображение, б - обработанное изображение
ЗАКЛЮЧЕНИЕ
Рассмотренная технология обработки изображений с использованием БИХ-фильтров является существенным развитием технологии, описанной в работе [1]. Включение в качестве одного из этапов технологии проверки устойчивости полученного в результате идентификации БИХ-фильтра обеспечивает повышение качества и существенное повышение надежности обработки крупноформатных изображений.
Наиболее важным результатом работы является рекуррентная реализация БИХ-фильтра на ОРИ-процессорах. В данном случае путем орги-
нальной организации процедур обработки удалось преодолеть традиционно существовавшее мнение, что эффективная реализация рекурсивных процедур на GPU-процессорах невозможна. Также интересным для дальнейших исследований является анализ эффективности параллельной реализации БИХ-фильтров в RDMA системах.
СПИСОК ЛИТЕРАТУРЫ
1. Милюткин М. Г., Никоноров А. В., Фурсов В. А. Параллельная реализация двумерных БИХ-фильтров в распределенной системе обработки изображений // Труды международной конференции ПаВТ 2010. С. 268-275.
2. Sung W., Mitra S. K. Effiicient Multi-Proccessor Implementation of Recursive Digital Filter // ICASSP, 1986.
3. Mccool M. D. Signal Processing and General-Purpose Computing and GPUs // Signal Processing Magazine, IEEE. 2008. P. 109-114.
4. Мурызин С. А., Сергеев В. В., Фролова Л. Г. Исследование эффективности двумерных параллельно-рекурсивных КИХ-фильтров // Компьютерная оптика. М.: МЦНТИ, 1992. Вып. 12. С. 65-71.
5. Методы компьютерной обработки изображений / Под ред. В. А. Сойфера. М.: Физматлит, 2001.
6. Зимин Д. И., Фурсов В. А. Построение устойчивых алгоритмов обработки изображений путем аппроксимации фильтров с бесконечной импульсной характеристикой // Компьютерная оптика. 2005. № 28. C. 124-127.
7. Dudgeon D. E. , Mersereau R. M. Multidimensional digital signal processing. Prentice-Hall, Inc., Englewood Cliffs, 1984.
8. Боресков А. В., Харламов А. А. Основы работы с технологией CUDA. ДМК Пресс, 2010. 230 c.
9. Memory Optimization for Recurrent CUDA Image Processing / S. Bibikov [et al.] // PRIA-2010 Proceedings. 2010. P. 176-179.
10. NVIDIA CUDA C Best Practices Guide // Santa Clara, 2010. 75 p.
ОБ АВТОРАХ
Никоноров Артем Владимирович, доцент кафедры общей информатики Самарск. гос. аэро-космическ. ун-та. Дипл. инженер (Самарск. гос. аэрокосмическ. ун-т, 2002). Иссл. в обл. распознавания образов и анализа изображений, идентификации систем, параллельных и распределенных вычислений, вычислений с использованием графических процессоров.
Фурсов Владимир Алексеевич, проф., зав. той же каф. Д-р техн. наук. Иссл. в обл. теории и м е тодов оценивания по малому числу измерений, методов обработки и распознавания изображений, построения параллельных алгоритмов обработки и распознавания изображений, реализуемых с использованием многопроцессорных вычислительных систем.
Якимов Павел Юрьевич, магистрант той же каф. Дипл. бакалавр (Самарск. гос. аэрокосми-ческ. ун-т, 2009). Иссл. в обл. распознавания образов и анализа изображений, параллельных и распределенных вычислений, вычислений с помощью графических процессоров.
