Фильтр Гаусса на стероидах: подход на точность вычислений / forpes.ru

Главная
Фильтр Гаусса на стероидах: подход на точность вычислений

Фильтр Гаусса на стероидах: подход на точность вычислений +9

18.02.2025 15:14

SmartEngines 0 1300 Источник

Hello, world! Перед вами вторая часть хабростатьи Smart Engines, посвященной быстрой фильтрации изображений. Да-да, создавая топовый продукт по распознаванию документов, нам приходится разбираться в методах обработки изображений на экспертном уровне (иначе не получилось бы распознать изображение паспорта за 150 мс на мобильном телефоне). В первой части мы начали разбирать аппроксимации гауссовского фильтра, которые стали темой нашей недавней публикации в научном журнале MDPI Applied Sciences [1]. О том, как работает оригинальный фильтр Гаусса, уже упоминалось, сейчас только напомним о его использовании всюду, где возникает обработка изображений: от редактирования фотографий на смартфоне – для размытия фона за объектом в режиме "портрет", до анализа рентгеновских снимков – чтобы убрать шум и улучшить читаемость изображения.

В прошлый раз мы рассказали, что для аппроксимации фильтра Гаусса можно использовать КИХ- или БИХ-фильтры, и успели познакомить читателей с тремя вариантами КИХ-фильтров. Но не спешите закрывать этот текст, если эти термины вам не знакомы! Достаточно понимать следующее: КИХ-фильтр формирует выходной сигнал y[n] , опираясь только на входные значения x[n] , тогда как БИХ-фильтр, помимо x[n] , задействует также предыдущие значения y[k], k < n .

По сути, все рассмотренные в первой части методы приближали гауссиану композициями полиномов 0-го или 1-го порядков, за счет чего были очень быстрыми. Сейчас мы рассмотрим более сложные варианты: полиномы 2-го порядка, БИХ-фильтры, а также способ вычисления свертки через преобразование Фурье. Мы сгруппировали их вместе по принципу скорости — они работают чуть дольше, но итоговое приближение получается на порядок точнее.

Кусочно-параболическая аппроксимация

Данный метод по своему устройству относится к КИХ-фильтрам. Его идея заключается в том, чтобы приближать гауссиану участками парабол. Мы разбили отрезок $[-3\sigma;\, 3\sigma]$ на три подотрезка p_k, k = 1, 2, 3, так, что $\cup p_k = [-3\sigma; 3\sigma]$ , и приближали части гауссианы на каждом подотрезке куском параболы (Рис. 1). В наших экспериментах мы делили отрезок в точках $\pm \sigma$ .

Рис. 1. Аппроксимация гауссианы () кусочками парабол. — Рис. 1. Аппроксимация гауссианы ( $\sigma=10$ ) кусочками парабол.

Рис. 1. Аппроксимация гауссианы ( $\sigma=10$ ) кусочками парабол.

Выбрав три различные точки на каждом подотрезке p_k , мы можем решить систему из трех уравнений:

$\Big\{ a_k \, x_i^2 + b_k \, x_i + c_k = g_\sigma[x_i], \;\; i = 1, 2, 3, \;\; x_i \in p_k,$

и определить коэффициенты a_k, b_k, c_k соответствующей параболы.

Теперь задача свелась к вычислению сверток с константой, и x^2 . Вычисление "влоб" неоптимально, поэтому мы используем технику быстрого вычисления. Запишем выражение для результата свертки с параболой $h[n] = a_k \, n^2 + b_k \, n + c_k$ :

$y[i] = \sum_{j=i-M}^{i+M} \big(a_k (i-j)^2 + b_k (i-j) + c_k\big)\,x[j]=\\ = \; a_k \sum_{j=i-M}^{i+M} (i-j)^2\,x[j]\; + \; b_k \sum_{j=i-M}^{i+M}(i-j)\,x[j]\; + \; c_k \sum_{j=i-M}^{i+M}x[j] =\\ = \; a_k \, y^{F2}[i] + b_k \, y^{F1}[i] + c_k \, y^C[i].$

На каждом из трех подотрезков $p_k,\, k = 1, 2, 3,$ мы аппроксимируем гауссиану кусочком параболы, так что каждый подсегмент соответствует ограниченному диапазону. Проделав ряд математических преобразований (которые можно посмотреть в [1]), мы получили

$y^C[i + 1] = y^C[i] - x[i-M] + x[i+M+1],\\ y^F[i+1] = y^F[i] + y^C[i] - \big(M + 1\big) \cdot x[i-M] - M \cdot x[i+M+1],\\ y^{F2}[i+1] = y^{F2}[i] + 2\cdot y^{F}[i] + y^C[i] - (M + 1)^2 \cdot x[i-M] - M^2 \cdot x[i+M+1].$

Так, результат свертки для одного пикселя можно вычислить за 11 сложений и 8 умножений (за исключением крайнего левого пикселя). По сравнению с наивным вычислением это значительно ускоряет обработку.

Рекурсивные фильтры

Задача проектирования рекурсивных фильтров заключается в определении вещественных коэффициентов и b_k передаточной функции вида

$H(z) = \frac{\sum_{k=0}^{m-1} b_k \, z^{-k}}{1 + \sum_{k=1}^{n} a_k \, z^{-k}}.$

Она описывает рекурсивный фильтр порядка, который представляется в виде разностного уравнения

$y[i] = \sum_{k=0}^{m-1} b_k \, x[i - k] - \sum_{k=1}^{n} a_k \, y[i - k].$

Коэффициенты подбираются так, чтобы передаточная функция H(z) наилучшим образом приближала, в соответствии с заданным критерием ошибки, передаточную функцию дискретного гауссовского фильтра для заданного значения $\sigma_0$ . Затем коэффициенты масштабируются на остальные значения $\sigma$ .

1) Фильтр Дериша

В работе 1993 года [2] Рашид Дериш предложил рекурсивные фильтры для аппроксимации гауссианы, её первой и второй производных. Передаточная функция фильтра представлялась в виде суммы каузальной и антикаузальной частей:

$H(z) = H^{+}(z) + H^{-}(z).$

Для случая гауссовской фильтрации автор сконструировал фильтр 4-го порядка. С целью снижения вычислительной сложности, в [1] мы рассмотрели фильтры 1-го и 2-го порядков, строящиеся по аналогичной схеме. Здесь для краткости изложения мы остановимся на одном варианте.

В случае фильтра 2-го порядка каузальная часть $H^{+}$ и антикаузальная часть $H^{-}$ имеют вид

$H^{+}(z) = \frac{b_0^{+} + b_1^{+} \, z^{-1}}{1 + a_1\, z^{-1} + a_2 \,z^{-2}}, \qquad H^{-}(z) = \frac{b_1^{-} z + b_2^{-} \,z^2}{1 + a_1 \,z + a_2 \,z^2}.$

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

$y[k] = y^{+}[k] + y^{-}[k],$

где

$y^{+}[k] = b_{0}^{+}\, x[k] + b_{1}^{+}\, x[k-1] - a_{1}^{+}\, y^{+}[k-1] - a_{2}^{+}\, y^{+}[k-2],\\ y^{-}[k] = b_{1}^{-}\, x[k+1] + b_{2}^{-}\, x[k+2] - a_{1}^{-}\, y^{-}[k+1] - a_{2}^{-}\, y^{-}[k+2].$

Несложно видеть, что вычисление такого фильтра требует 7 сложений и 8 умножений на элемент. Импульсная характеристика берется в виде

$h[n] = \Big(\gamma_1\, \cos\left(\frac{\omega\, n}{\sigma}\right) + \gamma_2\, \sin\left(\frac{\omega\, n}{\sigma}\right)\Big) \cdot e^{- \frac{b\, n}{\sigma}}, \quad \gamma_1,\, \gamma_2,\, \omega,\, b \in \mathbb{R}.$

Оптимальные значения $\gamma_1, \gamma_2, \omega$ и ищутся из условия

$\sum_{n=0}^{\infty} (g_\sigma[n] - h[n])^2 \to \min$

и подставляются в выражение для коэффициентов фильтра. Пример полученной аппроксимации представлен на Рис. 2.

Рис. 2. Аппроксимация гауссианы () фильтром Дериша 2-го порядка. — Рис. 2. Аппроксимация гауссианы ( $\sigma= 10$ ) фильтром Дериша 2-го порядка.

2) Фильтр VYV

Этот фильтр был предложен в 1998 году группой из трех авторов [3] Лукасом ван Влиетом, Ианом Янгом и Питом Вербиком. Мы будем сокращенно называть его VYV-фильтром. Если у Дериша передаточная функция системы состояла из суммы каузальной и антикаузальной частей, то для VYV-фильтра — это их произведение:

$H(z) = H^{+}(z) \cdot H^{-}(z).$

В [1] мы описали фильтры 2-го и 3-го порядков, построенные по авторской схеме. Для фильтра 1-го порядка точность получилась достаточно низкой (около 8e-2), поэтому мы не стали включать его в текст. Здесь мы рассмотрим фильтр порядка 2. Каузальная и антикаузальная части задаются через

$H^{+}(z) = \frac{\alpha}{1 + b_1 \,z^{-1} + b_2 \,z^{-2}},\quad H^{-}(z) = \frac{\alpha}{1 + b_1 \,z + b_2 \,z^2},$

и выражения для фильтра имеют следующий вид:

$y^+[k] = \alpha\, x[k] - b_1\, y^+[k-1] - b_2\, y^+[k-2],\\y^-[k] = \alpha\, y^+[k] - b_1\, y^-[k+1] - b_2\, y^-[k+2].$

Вычислительная стоимость – 4 сложения и 6 умножений на элемент. Стоит обратить внимание на иное устройство VYV-фильтра в сравнении с фильтром Дериша: первый при вычислении y^-[k] задействует значения y^+[k] . По этой причине, например, не получится вычислять y^+[k] и y^-[k] параллельно для фиксированной строки изображения. Коэффициенты получились путем минимизации функционала

$L^2 = \sqrt{\frac{1}{2 \pi} \int_{-\pi}^\pi \big(H\left(\omega ; \sigma_0\right)-G\left(\omega ; \sigma_0\right)\big)^2 \mathrm{d} \omega},$

в котором $G_{\sigma_0}$ — передаточная функция гауссианы. На Рис. 3 можно видеть, что полученный результат дает хорошее приближение.

Рис. 3. Аппроксимация передаточной функции гауссианы () фильтром VYV 2-го порядка. — Рис. 3. Аппроксимация передаточной функции гауссианы ( $\sigma=10$ ) фильтром VYV 2-го порядка.

Вычисление свертки через преобразование Фурье

Вычисление свёртки с гауссовским ядром можно ускорить, если перейти из пространственной области в частотную через дискретное преобразование Фурье (ДПФ) и использовать способ его быстрого вычисления. В частотной области свёртка преобразуется в покомпонентное умножение спектров сигналов:

$F(y) = F(g_\sigma * x) = F(g_\sigma) \cdot F(x).$

Поэтому, сначала вычисляют преобразование Фурье F(x) исходного сигнала, затем умножают его на спектр гауссовского ядра $g_\sigma$ , который имеет аналитическое представление

$F(g_\sigma) = e^{-\frac{2\pi^2 \sigma^2 k^2}{N^2}},$

и после этого преобразуют результат в пространственную область при помощи обратного преобразования Фурье:

$y = F^{-1} \big( F(g_\sigma) \cdot F(x) \big).$

Важно учитывать, что оба сигнала должны иметь одинаковую длину. Если это не так, более короткий сигнал следует дополнить нулями (zero-padding).

С вычислительной точки зрения, сложность свёртки одномерного сигнала с использованием быстрого преобразования Фурье составляет $O(N \log N)$ для всего сигнала. Или $O(\log N)$ в пересчете на один элемент (или пиксель).

Сводка по точности и сложности

Далее мы дополнили таблицу 1 “сложность-точность”, которую начали заполнять в прошлый раз. Мы добавили значения среднеквадратичной ошибки и сложности вычислений на пиксель (с учетом того, что для приближения гауссовского фильтра одномерным фильтром нужно выполнить 2 прохода в горизонтальном и вертикальном направлениях). Точность добавленных методов на порядок выше, чем у ранее рассмотренных вариантов. В то же время методы требуют больше арифметических операций.

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

Метод	# +	# ·	Ошибка	Доп. память
Фильтр Гаусса	2(K-1)	2K	зависит от K	О(1)
Stack blur	12	0	9.35e−6	О(1)
Bell blur	28	0	4.40e−6	О(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})
Кусочно-параболическая	22	16	1.92e-7	импл. зависимо
ДПФ	O(log w+log h)	O(log w+log h)	–	O(1)
Фильтр Дериша (2 пор.)	14	16	1.80e-7	O(1)
Фильтр VYV (2 пор.)	8	12	1.39e-7	O(1)

Тестирование на устройствах: больше методов, больше результатов

Наше тестирование заключалось в запуске запрограммированных на C++ методов на процессорах Intel Core i9-9900KF (с архитектурой x86_64) и ARM Cortex-A73 (с архитектурой ARMv8). Для свертки через ДПФ мы использовали реализацию из библиотеки fftw3 [4]. Тестировали 2 типа входных данных: 8-битные беззнаковые целые (uint8) и 32-битные с плавающей точкой (real32). Для целочисленных реализаций мы делали квантование весов. Из рассмотренных сегодня методов такой трюк проделывался для кусочно-параболической аппроксимации. Для БИХ-фильтров квантованные реализации оказались неудачными, т.к. при обходе строки изображения накапливались слишком большие значения, которые не помещались в используемый тип аккумулятора. Аппроксимации были распараллелены с помощью SIMD: на x86_64 использовались SSE-интринсики, на ARMv8 — NEON. Задействованы 128-битные регистры, хранящие 4 числа с плавающей точкой или 16 целых. Каждый метод тестировался на 1000 случайных изображениях $1024 \times 1024$ , время усреднялось по серии запусков. В таблице 2 приведены численные результаты.

Таблица 2. Среднее время работы рассмотренных аппроксимаций ( $\sigma=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
Кусочно-параболическая (real32)	20.0	5.8	47.9	34.8
Кусочно-параболическая (uint8)	15.0	4.9	65.2	29.8
ДПФ	10.7	–	86.1	–
Фильтр Дериша (2 пор.)	12.5	4.2	38.2	22.7
Фильтр VYV (2 пор.)	12.2	3.4	36.0	18.9

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

Заключение

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

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

Источники:

E. O. Rybakova, E. E. Limonova, D. P. Nikolaev, “Fast Gaussian Filter Approximations Comparison on SIMD Computing Platforms,” Appl. Sci. 2024, 14, 4664.
R. Deriche, “Recursively implementing the Gaussian and its derivatives,” Technical Report 1893, INRIA, Unit´e de Recherche Sophia-Antipolis, 1993.
L. J. van Vliet, I. T. Young, P. W. Verbeek, “Recursive Gaussian derivative filters,” in Proc. ICPR, Brisbane, QLD, Australia, 1998, pp. 509–514, vol. 1.
https://www.fftw.org/