28.01.2025 13:39
SmartEngines 21 2600 Источник

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

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

В мае 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}{2\pi\sigma^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}{2\pi \sigma^2} \sum_{m=-M}^M\sum_{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]=\\=2\pi \sigma^2\sum_{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, на отрезке [-3\sigma, 3\sigma] (Рис. 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 точек, равномерно взятых на интервале [-3\sigma_0, 3\sigma_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.

Комментарии (21)


  1. Sdima1357
    28.01.2025 16:33
    #27846292

    1. В такой статье обойтись без упоминания https://en.m.wikipedia.org/wiki/Central_limit_theorem ? Повторная свертка прямоугольников (бегущее среднее) очень быстро сходится к гауссиану.

    1. Когда ошибка 10^-6, картинки показывать не нужно, они одинаковые


    1. Refridgerator
      28.01.2025 16:33
      #27848486

      Делал такой трюк для полосового фильтра с линейной ФЧХ - 4 скользящих средних с пропорциями 5,6,7,8 дают подавление боковых лепестков в 65 дБ. 6, 7, 8, 9, 10 - 80 дБ.


  1. adeshere
    28.01.2025 16:33
    #27847164

    Глянул на годы выхода публикаций. Для меня это просто шок-контент ;-) Я был уверен, что все это придумано и описано лет на 50 раньше. Мы у себя в программе используем подобные трюки где-то с 1989г, когда у нас встал вопрос о скорости вычислений. Нам тогда и в голову не приходило, что в этом может быть какая-то новизна.. Единственное существенное отличие - у нас в сигналах присутствует много Nan-значений, которые приходится обрабатывать особо. Из-за них алгоритмы получаются совсем не такими красивыми, как у классиков :-( Ну и еще у нас временные ряды, то есть задачка во всех случаях исключительно одномерная.


    1. Sdima1357
      28.01.2025 16:33
      #27847470

      Нет тут особой новизны. У меня тоже лет 20 назад были примерно такие же результаты как в статье( для feature detection DOG и похожих ) но на заметно более слабых компьютерах. Точные цифры не помню, надо посмотреть реализацию и замеры для точных цифр, но было быстрее чем open cv раз в пять Тут важно правильно заботиться о попадании в кеш


      1. adeshere
        28.01.2025 16:33
        #27847562

        Тут важно правильно заботиться о попадании в кеш

        На 8086? Или на 80286?

        Я тогда такого слова вообще не знал ;-)

        Хотя, в начале 1990-х уже как бы появился 486 (с кэшем L1), но мы жили и работали в горной экспедиции, куда новая техника подъезжала с опозданием лет эдак на 10. Да и в институте (куда я вернулся после распада СССР) в 90-е было не сильно лучше. Приходилось выкручиваться на старом железе за счет алгоритмов. Так как как раз в это время наблюдения массово переводились на минутный опрос, и мегабайтные ряды у нас стали нормой.


        1. Sdima1357
          28.01.2025 16:33
          #27848306

          Я уже про 20 лет назад говорю. Это пентиум 4 и дальше. А кеши появились уже у некоторых 80386 , это где то 1986 год


          1. adenis78
            28.01.2025 16:33
            #27849458

            Вопрос с кэшем оказался сложнее чем показался вначале :)


        1. adenis78
          28.01.2025 16:33
          #27849654

          пока проверял, когда появились контроллеры кэша для 386го и как кэш был устроен, истекло время редактирования.
          Вот про микросхему 82385, контроллер совмещнного кэша программ-данных для intel 80386 https://theretroweb.com/chips/3095
          К ней требовались еще микросхемы быстрой статической памяти для хранения кэшируемых данных.
          Первые экземпляры = 1987 год.


  1. c0mmandor
    28.01.2025 16:33
    #27847364

    в таблице производительность в наносекундах на элемент или в миллисекундах на все изображение? явно не сказано, а кажется сомнительным что фильтр гаусса с сигмой 10 для FP32 изображения 1024х1024 можно посчитать на x86 процессоре за 81.3 нс


    1. Sdima1357
      28.01.2025 16:33
      #27847446

      На точку конечно. 4 мегабайта данных только прочитать из памяти нужно порядка миллисекунды или сотен микросекунд (10-89 гигабайт в секунду)


    1. Jijiki
      28.01.2025 16:33
      #27849730

      1024x1024 уже в памяти - когда я жму на применить фильтр к картинке например в Гимп на ползунки X/Y фильтр применяется на всю картинку почти мгновенно, тут скорее всего время без загрузки картинки 1024х1024 и в примере чб(поидее это 1024х1024х2хsizeof(...) байт ) и еще процессор i9 в тесте написано по крайней мере (мне просто стало интересно я запустил Гимп создал чб 1024х1024(32FP) накалякал чтото там чтоб был чб и применил фильтр, плюс в Гимп может модный фильтр, тут пишут что способ Stack Blur разновидност Гаусса быстрее)


      1. c0mmandor
        28.01.2025 16:33
        #27849996

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


        1. Jijiki
          28.01.2025 16:33
          #27850190

          допустим 1024 на 1024 на 2 канала частота 4 ГГц доступ к памяти 3МГЦ, допустим выборка окна 4 на 4 где 4 на 4 в лучшем случае на 2 канала в 1 регистр тогда за такт выполнится 1 к операций поидее ну допустим при 4 ГГЦ это 4 тысячи, 65536 / 4000 (просто предположил)


          1. c0mmandor
            28.01.2025 16:33
            #27850476

            тут как чего не допускай, но за 81 наносекунду не может быть обработано 1048576 элементов (это 1024*1024), т.к. это около 77 фемто секунд на элемент получится (если арифметика меня не подвела или я ее) 81e-9 / (1024*1024). это ж процессор на терра герцах работать должен, какой бы он i9 ни был, но это пока рановато


            1. Jijiki
              28.01.2025 16:33
              #27850510

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


              1. c0mmandor
                28.01.2025 16:33
                #27850582

                да нет, тут как ни крути, оценка в 81 наносекунду для всего изображения не физична. просто в таблице стоит добавить что это nsec / pixel тогда получается 84 миллисекунды на все изображение, вот в это поверю.


    1. SmartEngines Автор
      28.01.2025 16:33
      #27851308

      Вы правы, наносекунды получатся, если оценивать время на пиксель изображения. В Таблице 2 приведено время обработки всего изображения, поэтому правильная единица измерения - миллисекунды. Спасибо, исправили в статье.


      1. c0mmandor
        28.01.2025 16:33
        #27851518

        вам спасибо, статья получилась интересная и познавательная


  1. Andy_U
    28.01.2025 16:33
    #27849824

    А со скоростью вычисления свертки с помощью преобразования Фурье вы свои свои алгоритмы не сравнивали?


    1. SmartEngines Автор
      28.01.2025 16:33
      #27851296

      Да, это сравнение мы тоже делали. Оно есть в статье, а здесь мы напишем про него в следующей части нашего материала. Спойлер: оказалось не быстрее.


  1. lgorSL
    28.01.2025 16:33
    #27851500

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

    Главный подвох в том, что количество фотонов от пикселя яркостью в 127 и 255 отличается совсем не в два раза, а примерно в четыре. Если быть точным, в (255/127)^2.1 раз