Фильтр Гаусса на стероидах: секреты ускорения вычислений. smartengines.. smartengines. аппроксимация.. smartengines. аппроксимация. Блог компании Smart Engines.. smartengines. аппроксимация. Блог компании Smart Engines. гауссовский фильтр.. smartengines. аппроксимация. Блог компании Smart Engines. гауссовский фильтр. математика.. smartengines. аппроксимация. Блог компании Smart Engines. гауссовский фильтр. математика. Обработка изображений.. smartengines. аппроксимация. Блог компании Smart Engines. гауссовский фильтр. математика. Обработка изображений. ускорение.

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

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

В мае 2024 года группа исследователей Smart Engines опубликовала статью [1] в журнале MDPI Applied Sciences, в которой проводилось сравнение нескольких методов, предложенных для ускоренного приближенного вычисления гауссовской фильтрации. Сегодня мы расскажем о том, как она устроена математически и познакомим читателей с алгоритмами, которые помогают значительно быстрее посчитать ее приближенные версии.

Как работает фильтр Гаусса?

Основная задача гауссовской фильтрации — сглаживать изображение, уменьшая шум и устраняя ненужные детали, которые могут мешать анализу. Примеры до/после применения фильтра представлены на Рис. 1.

Рис. 1. Результаты применения гауссовской фильтрации.

Рис. 1. Результаты применения гауссовской фильтрации.

С математической точки зрения обработка изображения фильтром Гаусса есть свертка изображения с гауссовским ядром. На практике рассматриваются ядра конечного размера. Каждый пиксель заменяется взвешенным средним его окружения внутри окна размера K times K, приложенного в рассматриваемом пикселе. Чем ближе соседние пиксели к текущему, тем больший “вклад” они дают в ответ. Рис. 2 иллюстрирует описанный процесс.

Рис. 2. Свертка изображения с ядром 5 x 5.

Рис. 2. Свертка изображения с ядром 5 x 5.

Формально гауссовское ядро, далее обозначаемое g_sigma, определяется уравнением

g_sigma [m, n]=frac{1}{2pisigma^2} e^{-frac{1}{2}left(frac{m^2 + n^2}{sigma^2}right)},

где sigma– стандартное отклонение. Процедура гауссовой фильтрации имеет вид:

y[i, j]=sum_{m=-M}^M sum_{n=-M}^M g_sigma[m, n] cdot x[i-m, j-n]

Тут K^2умножений и K^2-1сложений на пиксель, а x[i, j] обозначает значение (i, j)-го пикселя исходного изображения, y[i, j] – значение(i, j)-го пикселя результирующего изображения, а K=2M + 1. Согласно данному выражению, зависимость числа арифметических операций, необходимых для вычисления ответа, от размера окна квадратичная.

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

y[i, j]=frac{1}{2pi sigma^2} sum_{m=-M}^Msum_{n=-M}^M e^{-frac{1}{2}left( frac{m^2+n^2}{sigma^2}right)}x[i-m,j-n]=\=sum_{m=-M}^M frac{1}{sqrt{2 pi}sigma}e^{-frac{1}{2}left(frac{m^2}{sigma^2}right)} cdot sum_{n=-M}^M frac{1}{sqrt{2 pi}sigma}e^{-frac{1}{2}left(frac{n^2}{sigma^2}right)} x[i-m, j-n]=\=2pi sigma^2sum_{m=-M}^M g_sigma[m] cdot sum_{n=-M}^M g_sigma[n] cdot x[i-m, j-n].

Так, свертка с двумерным гауссовым ядром может быть разложена на пару сверток с одномерными гауссовыми ядрами: вертикальную и горизонтальную. Отметим, что такое разложение само по себе не является приближением и по-прежнему дает точный результат. Присмотревшись, видим, что вычислительные затраты сократились до 2K умножений и 2(K-1) сложений на пиксель.

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

Подходы к приближенному вычислению

Существует два основных типа аппроксимаций: фильтры с конечной импульсной характеристикой (КИХ) и фильтры с бесконечной импульсной характеристикой (БИХ). 

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

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

y[i]=sum_{n=i-M}^{i+M}h[i-n] cdot x[n].

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

y[n]=b_0 x[n]+b_1 x[n-1] + dots+b_P x[n-P] - a_1y[n-1]-dots-a_Qy[n-Q].

 Коэффициенты b_i, i=0, dots, P, и a_j, j=1, dots, Q, необходимо настраивать.

С БИХ-фильтрами мы разберемся в следующий раз, а сейчас речь пойдет о нескольких аппроксимациях первого типа. Здесь мы постарались объяснить схему работы “на пальцах”. Заинтересовавшиеся могут найти строгое математическое изложение в [1].

Stack blur и Bell blur

Алгоритм Stack blur был предложен Марио Клингеманном [2]. Используется метод “скользящего окна”. Значения внутри окна выбраны так, чтобы крайние пиксели имели единичный вес, а затем веса увеличивались на 1 по мере продвижения к центральному пикселю – его вес будет наибольшим. Во время исполнения алгоритма окно перемещается вдоль строки или столбца (при горизонтальном или вертикальном проходах соответственно). При каждом наложении производится умножение элементов окна на соответствующие элементы изображения и последующее суммирование – пока это просто определение КИХ-фильтра. Весь смысл заключается в способе быстрого вычисления.

Обратимся к Рис. 3, который иллюстрирует то, как работает Stack blur в горизонтальном направлении с окном длины 5. Допустим, что мы посчитали ответ y[i] для пикселя x[i] со значением c и хотим посчитать ответ y[i+1] для пикселя x[i+1] со значением d. Для этого необходимо всего лишь добавить к y[i] сумму “входящих” элементов (отмеченных “+1”) и вычесть сумму “уходящих” элементов (отмеченных “-1”), т. е. затратить 1 сложение и 1 вычитание. На соседних итерациях алгоритма “входящие” элементы отличаются на пару крайних пикселей – чтобы пересчитать сумму “входящих” элементов, нужно выполнить 1 сложение и 1 вычитание. Аналогично для уходящих. Итого 6 аддитивных операций.

Рис. 3. Принцип работы метода Stack blur.

Рис. 3. Принцип работы метода Stack blur.

При таком алгоритме одномерное гауссовское ядро gσ аппроксимируется кусочно-линейной функцией

S_sigma[m]=r-|m|+1=begin{cases}r-m+1, & if~~~m>0,\ r+m+1,& otherwise, end{cases}

что на графике выглядит как уголок, близкий к гауссиане (Рис. 4).

Рис. 4. Аппроксимация гауссианы () методами Stack blur и Bell blur.

Рис. 4. Аппроксимация гауссианы (sigma=10) методами Stack blur и Bell blur.

Метод Stack blur можно усложнить, если заменить суммы входящих и уходящих пикселей на взвешенные суммы. Эта идея была описана в [3] и названа Quadratic Stack blur, или Bell blur, т. к. в качестве весов выбраны те же, что используются в алгоритме Stack blur, а полученное кусочно-постоянное приближение гауссианы визуально имеет форму колокола  (Рис. 4). Принцип работы аналогичен методу Stack blur и схематично представлен на Рис.5.

Рис. 5. Принцип работы метода Bell blur.

Рис. 5. Принцип работы метода Bell blur.

Итого для получения ответа y[i+1] нужно к ответу y[i] прибавить Stack blur по “входящим” пикселям и вычесть Stack blur по “уходящим” пикселям – всего 14 аддитивных операций. Математическое выражение для получающейся кусочно-постоянной функции-аппроксиматора гауссианы здесь мы не приводим.

Running sums

Еще одну кусочно-постоянную аппроксимацию гауссианы и способ быстрого вычисления свертки с ней предложили Эльбохер и Верман [4]. В отличие от метода Bell blur, в котором значения весов являются целыми числами, метод Running sums оперирует вещественными весами. Авторы определяют функцию из k ступенек с высотами c_i in mathbb{R}, i=1, dots, k, на отрезке [-3sigma, 3sigma] (Рис. 6). Высота и длина ступеней определяются путем минимизации среднеквадратичной ошибки.

Рис. 6. Аппроксимация гауссианы () для .

Рис. 6. Аппроксимация гауссианы (sigma=5) для k=3, 4.

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

I[i]=sum_{n=0}^i x[n],

так что в i-том элементе массива содержится сумма значений всех пикселей в строке/столбце от нулевого до i-го. Для того чтобы понять, как это упрощает вычисления, обратимся к Рис. 7. Он иллюстрирует идею вычисления для случая k=4.

Рис. 7. Идея метода Running sums для .

Рис. 7. Идея метода Running sums для k=4.

Запишем выражение для результата фильтрации в пикселе x[i] со значением e и переупорядочим в нем слагаемые:

c_1(d+e+f)+c_2(c+g)+c_3(b+h)+c_4(a+i)=\=(c_1-c_2)(d+e+f)+\+(c_2-c_3)(c+d+e+f+g)+\+(c_3-c_4)(b+c+d+e+f+g+h)+\+c_4(a+b+c+d+e+f+g+h+i).

Левая часть данного тождества соответствует суммированию по вертикальным цветным прямоугольникам, правая – по горизонтальным штрихованным. Обозначив w_i=c_i - c_{i + 1}, i=1, dots, k (c_{k+1}=0) и выразив второй множитель внутри каждого слагаемого через элементы массива I, можно переписать результат следующим образом:

w_1(I[i+1]-I[i-2])+\+w_2(I[i+2]-I[i-3])+\+w_3(I[i+3]-I[i-4])+\+w_4(I[i+4]-I[i-5]).

То есть для получения результата необходимо затратить 4 вычитания, 4 сложения (1 идет из заполнения массива I) и 4 умножения. В случае произвольного k метод требует 2k аддитивных и k мультипликативных операций на пиксель. Заметим также, что метод требует использование дополнительной памяти для хранения кумулятивной суммы по строке или столбцу – O(max {h, w}), где h и w – высота и ширина изображения.

Что по точности и сложности?

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

MSE=frac{1}{S}sum_{n=1}^S(g_{sigma_0}[n] - h[n])^2,

для S=500 точек, равномерно взятых на интервале [-3sigma_0, 3sigma_0] с sigma_0 = 10. Результаты представлены в Таблице 1.

Таблица 1. Сравнение аппроксимаций гауссовского фильтра для изображения размера h times w по количеству аддитивных (# +) и мультипликативных (# ·) операций на пиксель, ошибке и дополнительной памяти. K – размер ядра.

Метод

# +

# ·

Ошибка

Доп. память

Фильтр Гаусса

2(K − 1)

2K

зависит от K

O(1)

Stack blur

12

0

9.35e−6

O(1)

Bell blur

28

0

4.40e−6

O(1)

Running sums, k = 3

12

6

1.43e−5

O(max{h, w})

Running sums, k = 4

16

8

7.91e−6

O(max{h, w})

Тестирование на устройствах

Мы запрограммировали методы на языке C++ и опробовали их на процессорах Intel Core i9-9900KF (с архитектурой x86_64) и ARM Cortex-A73 (с архитектурой ARMv8). Тестировали 2 типа входных данных: 8-битные беззнаковые целые (uint8) и 32-битные с плавающей точкой (real32). Фильтр Гаусса и метод Running sums используют вещественные коэффициенты – для целочисленных реализаций мы делали квантование весов. Также аппроксимации были распараллелены по данным при помощи SIMD расширений. На x86_64 применялись интринсики семейства SSE, на ARMv8 – интринсики семейства NEON. Использовались 128-битные регистры, способные хранить 4 значения с плавающей точкой одинарной точности или 16 беззнаковых целых значений. Каждый метод запускался на 1000 случайно инициализированных изображениях размером 1024 x 1024, время работы усреднялось по серии этих запусков. Численные результаты приведены в Таблице 2.

Таблица 2. Среднее время работы рассмотренных аппроксимаций (sigma_0=10).

Метод

Время работы, нс (x86_64)

Время работы, нс (ARMv8)

без SIMD

с SIMD

без SIMD

с SIMD

Фильтр Гаусса (real32)

81.3

281.5

Фильтр Гаусса (uint8)

12.0

128.0

Stack blur (real32)

6.0

3.5

20.4

17.7

Stack blur (uint8)

3.8

1.6

15.5

14.7

Bell blur (real32)

6.7

4.3

31.6

26.8

Bell blur (uint8)

5.7

2.6

22.5

22.0

Running sums, k = 3 (real32)

8.3

5.6

40.7

23.2

Running sums, k = 3 (uint8)

7.0

2.7

32.9

14.1

Running sums, k = 4 (real32)

9.7

5.9

46.8

25.0

Running sums, k = 4 (uint8)

8.3

3.0

39.4

15.7

Метод Stack Blur в реализации над целыми числами оказался самым быстрым на x86_64. На ARMv8 его обогнал метод Running sums в 3-ступенчатом варианте. В сравнении с настоящей гауссовской фильтрацией эти методы оказались быстрее в 50 и 20 раз соответственно.

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

Рис. 7. Результат применения гауссовской фильтрации в точном и приближенном вариантах: (а) исходное изображение; (б) фильтр Гаусса (real32); (в) фильтр Гаусса (uint8); (г) Stack blur (real32); (д) Stack blur (uint8); (е) Bell blur (real32); (ж) Bell blur (uint8); (з) Running sums,  (real32); (и) Running sums,  (uint8); (к) Running sums,  (real32); (л) Running sums,  (uint8).

Рис. 7. Результат применения гауссовской фильтрации в точном и приближенном вариантах: (а) исходное изображение; (б) фильтр Гаусса (real32); (в) фильтр Гаусса (uint8); (г) Stack blur (real32); (д) Stack blur (uint8); (е) Bell blur (real32); (ж) Bell blur (uint8); (з) Running sums, k=3 (real32); (и) Running sums, k=3 (uint8); (к) Running sums, k=4 (real32); (л) Running sums, k=4 (uint8).

Заключение

Сегодня мы рассмотрели три способа приближенного вычисления гауссовского фильтра. Метод Stack blur показал себя как самый быстрый на платформе x86_64, благодаря минимальным вычислительным затратам. В то же время метод Running sums в 3-ступенчатом варианте продемонстрировал высокую производительность на ARMv8, что делает его отличным выбором для мобильных устройств и встраиваемых систем. Метод Bell blur обеспечивает более точное приближение гауссовского ядра при незначительном увеличении сложности расчетов.

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

Источники:

  1. Rybakova, E. O.; Limonova, E. E.; Nikolaev, D. P. “Fast Gaussian Filter Approximations Comparison on SIMD Computing Platforms,” Appl. Sci. 2024, 14, 4664.

  2. https://underdestruction.com/2004/02/25/stackblur-2004/

  3. https://observablehq.com/@jobleonard/mario-klingemans-stackblur

  4. Elboher, E.; Werman, M. “Efficient and accurate Gaussian image filtering using running sums,” Proc. ISDA, 2012, стр. 897–902.

Автор: SmartEngines

Источник

Rambler's Top100