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

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

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

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

Van Herk / Gil-Werman (vHGW) является одним из самых быстрых алгоритмов для вычисления дилатации и эрозии бинарных и полутоновых изображений на последовательном процессоре.

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

Для распараллеливания вычислений использованы технологии OpenMP и CUDA, так как они в большей степени обладают высокой производительностью.

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

Алгоритм vHGW (van Herk Gil-Werman), описан ван Херком [10], Гилем и Верманом [6] в своих работах.

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

Для полутоновой морфологии дилатация вычисляется по формуле 1, а эрозия по формуле 2.

                      (1)

                                                                                                                              (2)

Для эрозии, вместо максимального значения нужно искать минимальное.

Рис. 1. vHGW алгоритм для горизонтального структурного элемента размера p

Fig. 1. vHGW algorithm for horizontal structural element of p size

          I.          Изображение разделяется на отрезки длины p с (p-1)/2 на каждой стороне, чтобы сформировать окно размером 2p-1, с центром в точке р-1, 2p-1, 3p-1, ...

       II.          Для каждого пикселя k = [0;p-1] в заданном окне w, массив префиксов R рассчитывается для пикселов слева от центра

a)                  ,

а массив суффиксов S рассчитывается для пикселов справа от центра (p-1) ... (2p-2),

b)                 

c)                  R [k] и S [k] объединяются вместе, чтобы вычислить максимум фильтра

     III.          Для каждого пикселя 
(p-1)/2<=j<p+(p-1)/2 в w (отрезок длины p) результат дилатации

В результате, мы не оцениваем первые p / 2 пикселей и, в худшем случае, последние 3p / 2 пикселей. Чтобы получить разумные результаты на всех пикселях в изображении, мы вкладываем реальное изображение в более крупное изображение увеличенных размеров, где добавлены пограничные пиксели, инициализированные 0 для дилатации и 255 для эрозии.

Алгоритм vHGW выполняет дилатацию с помощью SE размером р = 2n + 1 за O (n) времени (n = число пикселей изображения). При этом выполняет менее трех сравнений на пиксель, независимо от размера SE.

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

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

Поэтому, чтобы иметь возможность распараллелить программу с помощью OpenMP необходимо вынести некоторые вычисления в отдельные функции. После этого можно вставить директиву #pragmaompparallelforprivate(prefix, suffix) перед основным циклом алгоритма.

Технология OpenMP использует для распараллеливания модель «ветвление-слияние» (fork-joinparallelism), показанную на рисунке 2. Программа начинает выполняться в главном потоке, и когда встречается параллельная секция, создаются новые потоки, и нагрузка распределяется уже по нескольким потокам, после того как все потоки выполнились, снова начинает работать только один главный поток. С точки зрения программирования – это очень удобный подход, т.к. необходимо в последовательную часть программы вставить всего лишь одну или несколько директив [1].

Рис. 2. Модель «ветвление-слияние»

Fig. 2. Model "fork-joint parallelism"

Представленная реализация vHGW на CUDA использует следующие детали:

• Отдельный код для горизонтальной и вертикальной версии структурирующего элемента.

• Один поток вычисляет максимальное (минимальное) значения сразу для двух массивов S [K], R [K]

• Разделяемая память не используется, поскольку это ограничивает использование мультипроцессоров CUDA с увеличением размера изображений и структурирующих элементов.

При этом необходимо учесть следующие определения CUDA:

• CUDA мультипроцессор создаёт, управляет, планирует и выполняет потоки в группах из 32 параллельных нитей, называемых warp’ами [9]

• CUDA мультипроцессор содержит ряд блоков нитей, с максимумом 8 блоков нитей в CUDA CC версии 2.х и ниже.

• Размер блока должен быть кратным 
32 нитям (warp) с рекомендуемым минимумом 
64 (2 warp’а).

Если предположить, что возможно использовать только 32 потока (один warp), мы можем выполнить горизонтальную дилатацию на изображении с 1100 столбцов в строке и структурного элемента размером 11, с использованием более 32KB разделяемой памяти следующим образом:

• структурный элемент представляет собой горизонтальную линию 11 пикселов, поэтому р=11;

• каждая строка изображения составляет 1100 пикселей, поэтому 1100 / р = 100 окон;

• 2 массива в одном окне с помощью 2р значений и 2 нитей;

• каждый блок использует 32 потока;

• поэтому каждый блок требует 100*2р*32/2≈34K разделяемой памяти, которая равна 34KB для 8 битных целочисленных изображений или 136 КБ для 32 битных изображений с плавающей точкой;

• каждый мультипроцессор имеет 48kb разделяемой памяти для карт, поддерживающих CUDA CC 2,0 или выше [9].

71% разделяемой памяти мультипроцессора расходуется на изображение с 8 битами точности (для структурирующего элемента, размером 11 пикселей), а также превышает её в 1,8 раза для 32 битного изображения с плавающей точкой.

Распределение ресурсов, которое выделяет все 48kb разделяемой памяти на один блок с 32 нитями в блоке только 1 активный warp для этого мультипроцессора. Только с одним активным warp’ом из возможных 48 (CUDACC 2.0) [7] заполняемость (коэффициент использования нити) всего 2% от возможностей мультипроцессора, что значительно ограничивает потенциал производительности.

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

Для тестирования алгоритма морфологической обработки изображении vHGW, реализованного с помощью технологий OpenMP и CUDA, исходное изображение, представленное на рисунке 3, подавалось на вход 100 раз для более точного измерения времени его выполнения. Каждый проделанный эксперимент проводился с одним изображением, но разного разрешения. После выполнения обработки изображений, программа производила подсчет времени выполнения операций над изображением.

Рис. 3. Исходное изображение земной поверхности

Fig. 3. The original image of the Earth's surface

На рисунке 3 представлено исходное полутоновое изображение земной поверхности, снятое из космоса, а на рисунке 4 – результат эрозии этого изображения структурирующим элементом размера 3х3. Заметно, что изображение стало немного темнее, что вполне логично, учитывая то, что алгоритм находит локальный минимум интенсивности в заданной окрестности.

Рис. 4. Результат эрозии исходного изображения

Fig. 4. The result of erosion of the original image

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

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

Рис. 5. Исходное зашумленное изображение отпечатка пальца

Fig. 5. Original noisy fingerprint image

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

Fig. 6. The result of morphological filtering of a fingerprint image

Как видно из графика, представленного рисунке 7, при увеличении размера структурного элемента, время выполнения алгоритма не только не увеличивается, а даже уменьшается, как с использованием 1 потока, так и 12-ти потоков.

Рис. 7. Зависимость времени выполнения от размера структурирующего элемента для программы на 1 и 12 потоках

Fig. 7. Dependence of the execution time of the size of a structuring element for the application on 1 and 12 threads

На рисунке 8 видим, что с увеличением разрешения изображения программа с одним потоком работает гораздо медленнее, чем с 12-ю потоками.

Рис. 8. Зависимость времени выполнения от разрешения изображения для программы на 1 и 12 потоках

Fig. 8. Dependence of the execution time of the image resolution for the program on 1 and 12 threads

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

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

На рисунке 9. приводится график зависимости времени выполнения алгоритма от размера структурирующего элемента в трёх случаях:

- на 12 потоках CPU;

- с использованием стандартных морфологических операций на CUDA библиотеки OpenCV;

- с использованием алгоритма, рассматриваемого в данной работе на GPU.

Как видно из графика, представленного рисунке 9, алгоритм vHGW и его реализация на GPU гораздо эффективнее, чем алгоритмы, использующиеся в OpenCV.

Таким образом, можно сделать вывод, реализация алгоритма vHGW с использованием технологии NVIDIA CUDA позволяет получить ускорение приблизительно в 15 раз большее, по сравнению с использованием технологии OpenMP.

Рис. 9. Зависимость времени выполнения от размера структурирующего элемента

Fig. 9. Dependence of the execution time of the size of the structuring element

Рис. 10. Зависимость времени выполнения от размера изображения

Fig. 10. Dependence of the execution time of the image size

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

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