Быстрое размытие по Гауссу

    Фильтр размытия по гауссу (широко известный “gaussian blur” в фотошопе) достаточно часто применяется сам по себе или как часть других алгоритмов обработки изображений. Далее будет описан метод, позволяющий получать размытие со скоростью, не зависящей от радиуса размытия, используя фильтры с бесконечной импульсной характеристикой.

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

    Итак, пусть исходное изображение будет задано яркостью x(m,n). Гаусово размытие с радиусом r считается по формуле

    image

    Пределы суммирования по u и v можно выбирать как плюс минус несколько сигм, т.е. радиусов r, что даёт сложность алгоритма порядка O(r2) операций на пиксель. Для больших r и многомегапиксельных изображений не слишком здорово, правда?

    Первое ускорение даёт свойство сепарабельности гауссова размытия. То есть, можно провести фильтрацию по оси x для каждой строки, полученное изображение отфильтровать по y по каждому столбцу и получить тот же результат со сложностью O(r) операций на пиксель. Уже лучше. Это свойство мы тоже будем использовать, поэтому дальше все рассуждения будут для одномерного случая, где нужно получить y(n) имея x(n).

    В этом нам помогут фильтры с бесконечной импульсной характеристикой. Идея фильтра такова: значения y(n) рекуррентно рассчитываются по формуле:

    image

    Где ak и bi – некоторые предрасчитанные коэффициенты, а y(n) и x(n) при n<0 полагаются равными нулю.

    Часть формулы, зависящая от bi сводится к простой свёртке с конечным ядром. Поскольку мы хотим, чтобы фильтр считал побыстрее, то наш фильтр будем искать для случая P=0, то есть рассматривая один единственный коэффициент b0.

    На всякий случай, напомню, что любой линейный фильтр (а БИХ фильтр также является линейным) полностью характеризуется своим откликом на дельта функцию, равную 1 в точке 0, и 0 во всех остальных точках. Часть, зависящая от коэффициентов a при отклике на такую функцию либо разойдётся и уйдёт в бесконечность (этого хотелось бы избежать), либо даст нам красивый убывающий хвост. Например, вот такой:



    Похоже на гаусс? Ну да, что-то уже есть, но как-то оно кособоко выглядит. Поэтому идея алгоритма будет такой: мы отфильтруем в одну сторону (циклом от 0 до nmax), а потом результат в обратную (от nmax до 0), с теми же коэффициентами. Математически должна получиться строго симметричная кривая (поэтому не важно, фильтровать сначала туда потом обратно, или наоборот). Несколько забегая вперёд, вот что получится если профильтровать это ещё в обратном порядке:



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

    Итак, осталось рассчитать коэффициенты фильтра a. Далее все формулы пойдут для случая Q=3.

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

    image

    Или для нашего случая Q=3, P=1:

    image

    Поскольку в числителе и знаменателе находятся многочлены, их можно разложить на множители. Числитель, в общем случае, даст нам нули (в нашем случае нулей нет), а знаменатель – полюса:

    image

    в нашем случае 1/Z0, 1/Z1, 1/Z2. Заметим, что полюса однозначно определяют коэффициенты ak нашего фильтра (все БИХ фильтры, обычно, ищутся именно через полюса, а не напрямую). Кстати, получить коэффициенты по полюсам проще, чем наоборот (для этого нужно перемножить полиномы, а наоборот – решить степенное уравнение).

    Самый главный теоретический постулат, который нам пригодиться: если комплексные полюса лежат внутри единичного круга, т.е. по модулю меньше единицы, (или, что тоже самое, обратные значения Z0, Z1, Z2 лежат вне единичного круга), то фильтр будет устойчивым, то результат не убежит в бесконечность.

    Отметим также, что назвав произвольные комплексно сопряженные Z0, Z1 и вещественное Z2 (опять же, все по модулю больше единицы) мы можем построить по ним фильтр. В произведении всё превратится в строго вещественные коэффициенты a.

    Коэффициент b0 работает как “коэффициент громкости” фильтра.

    Итак, задача свелась к определению трёх коэффициентов Z0, Z1, Z2. Поскольку, повторюсь, два из них комплексно сопряженные, а третий вещественные, то нужно найти вещественную и мнимую часть от первого, и вещественную от третьего. То есть, три вещественных числа: real(Z0), im(Z0), Z2

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

    Дальше я пошёл немного другим путём, чем в работе «Recursive Gaussian Derivative Filters». Там они распространяли случай с r=2 на все другие, а я определил лучшие коэффициенты для всех радиусов от 2 до 2048, с экспоненциальным шагом. Я написал оптимизационный алгоритм, который ищет наиболее близкую кривую по методу минимизации модуля максимальной разницы. Дополнительным условием было сохранение фильтром полной энергии, т.е. чтобы функция x(n)=const переходила саму в себя, что даёт условие

    b0=1-(a1+a2+a3)

    Я перепробовал разные оптимизационные алгоритмы, и лучшие выдавал немного измененный генетический алгоритм. (Пожалуй, про оптимизацию можно написать отдельную заметку).

    Желающие могут посмотреть результаты на google spreadsheet.

    Видно, что для случая r=2 результаты отличаются от данных в работе. Почему так, я сказать не могу, при этом по моим расчетам мои коэффициенты дают меньшую, процентов на 40, погрешность.

    Дальше я воспользовался не буквой а духом статьи, и стал искать данные как

    real(Z0) = cos(W(r)/r)*eA(r)/r
    im(Z0) = sin(W(r)/r)*eA(r)/r
    Z2 = eB(r)/r

    Т.е. дальше нужно просто подобрать три похожие функции, каждая из которых стремится в пределе в константе. Я искал функции в виде отношения

    (k3 r3 + k2 r2 + k1 r + k0)/r3

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

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

    int gaussCoef(double sigma, double a[3], double *b0)
    {
    	double sigma_inv_4;
    
    	sigma_inv_4 = sigma*sigma; sigma_inv_4 = 1.0/(sigma_inv_4*sigma_inv_4);
    
    	double coef_A = sigma_inv_4*(sigma*(sigma*(sigma*1.1442707+0.0130625)-0.7500910)+0.2546730);
    	double coef_W = sigma_inv_4*(sigma*(sigma*(sigma*1.3642870+0.0088755)-0.3255340)+0.3016210);
    	double coef_B = sigma_inv_4*(sigma*(sigma*(sigma*1.2397166-0.0001644)-0.6363580)-0.0536068);
    
    	double z0_abs   = exp(coef_A);
    
    	double z0_real  = z0_abs*cos(coef_W);
    	double z0_im    = z0_abs*sin(coef_W);
    	double z2       = exp(coef_B);
    
    	double z0_abs_2 = z0_abs*z0_abs;
    	
    	a[2] = 1.0/(z2*z0_abs_2);
    	a[0] = (z0_abs_2+2*z0_real*z2)*a[2];
    	a[1] = -(2*z0_real+z2)*a[2];
    
    	*b0 = 1.0 - a[0] - a[1] - a[2];
    
    	return 0;
    };


    Всё! Теперь, нужно написать сам код. Расчёт, конечно, нужно производить во float, но современные компьютеры считают на числа с плавающей запятой (особенно с sse) довольно быстро. Программисты из Интел, кстати, озаботились оптимизацией Gauss-IIR фильтров под векторные инструкции процессора уже написали целую статью. Правда, там считают немножко другим методом, но основные способы оптимизации описаны хорошо.

    В конце можно дать пример того, что получилось:



    Картинка практически не отличается от «честной». Впрочем, если открыть её в фотошопе и поизучать внимательно, разницу найти можно.

    P.S. Это правда мой первый пост на Хабре, я мог что-то пропустить. Если будут вопросы, я отвечу в комментариях.
    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 43

      +1
      Поправьте меня меня если я не прав
      Так как фильтр Гауса является сепарабельным то количество операций на точку (при реализации горизонтальный +вертикальный проходы) не будет превышать O(r*2) (умножить) ценой 2x памяти.
        +1
        Во-первых, O® это тоже самое, что O(r*2). Это — для расчёта FIR фильтром.

        В описанном методе идёт O(1) для изображения фиксированого размера.

        Память нужна на источник, промежуточный буфер хранения, и приёмник. Если источник не нужен, или совпадает с приёмником, можно писать туда же.

        Если очень большое изображение и очень дорога память, можно размывать кусочками, но ценой лишних вычислений.
          0
          Про O(..) само собой — я просто хотел явно указать что не O(r^2)
          Память нужна ровно 2x если не нужно исходное изображение.
            0
            Ну да, но если считать просто свёртку лобовым методом, то память тоже нужна 2x
              0
              Про сепарабельность еще:
              В общем данный способ всегда используют в гейм деве для размытия.
              Гоняются данные с нужным кернелом между 2умя текстурами.
                0
                Некоторое время назад мы тоже были озабочены похожей оптимизацией для размытия в играх — но из за того, что алгоритм требует чтения и записи данных из одного и того же источника нам это не подошло.
                  0
                  А почему нельзя читать и записывать из одного и того же источника?
                    0
                    Результат при чтении и записи данных в одну и туже текстуру не определен.

                    Коротко о причине:
                    Из за кэшов на чтение/запись и затраты на их синхронизацию.
                      0
                      Это что за архитектура? Насколько я понял, рассчёт на видеокартах идёт?

                      И да, почему нельзя записывать в цикле y[i-3], а все промежуточные y держать в регистре каком-нибудь?
                        0
                        На GPU.
                        На GPU эти подпрограммы (шейдера) выполняются параллельно для каждого пиксела без какой либо связи с другим пикселами — именно по этому на регистрах не получится так как единственная связь между вызовами программы для каждого пикселя это источники данных (текстуры)

                          0
                          Ну я немного знаком с GPU, и по прежнему не понимаю, в чём проблема.

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

                          То, что вам нужно получать y как результат — ну не храните его в памяти. На момент расчёта y4 вам уже не нужен будет y0, и его можно писать в память.

                          Ну да, вам придётся отфильтровать текстуру 4 раза. Но сбросить весь кэш можно же между проходами, верно?
                            0
                            Когда GPU вызывает пиксельный шейдер для координаты (0,0) то результат этого пиксельного шейдера записывается только в (0,0) и никуда иначе.
                            Когда для (0,1) это происходит — то никаких регистров из (0,0) в (0,1) не переходит
                          0
                          Более того — никто не гарантирует что пиксель (0,1) будет рендерится сразу после (0,0)
                          Так что даже если бы поддерживалась возможность чтения записи из одного источника данных то ничего не работало бы
                            0
                            Ааааа, я понял.

                            Вам надо выдавать значения яркости пикселя по текстурным координатам.

                            Т… е имеем x,y нужно выдать f(x,y), но x и y — это куда уж нам повезло что луч туда попал.

                            Ну тогда да. Но такие вещи решаются прерасчётом и мипмэппингом.
                0
                В общем — при повторном чтении уже нашел про сепарабельность фильтра в статье — так что мой первый комментарии можно игнорить.
              0
              Не так хорошо дружу с математикой, как хотелось бы. Приведенная формула для y(n), насколько я понял, для одномерного случая. Как быть для двухмерного?
              Не могли бы вы привести сам код размытия, с помощью которого получилась финальная картинка?
                +1
                Для двумерного — сначала фильтруем по строкам, потом результат по столбцам. Ссылку на код на сайте интела я тоже дал.

                Двумерную картинку я вообще получил имитацией метода в excel, а потом перевёл из таблицы в текстовый pbm, вот так. Ну лень писать было :-)

                Ну вот код из оптимизатора, например, который по IIR фильтру фильтрует туда-обратно:

                int CSignal::setAsConvIIR(CSignal &signal, CFilterIIR &filter)
                {
                	if (signal.nOf!=nOf)
                		return -1;
                
                	int orderOfFilter = filter.orderOf();
                
                	flpointt* memTmp = (flpointt*) _alloca(sizeof(flpointt)*nOf);
                	
                	for(int i=0; i<nOf; i++)
                	{
                		memTmp[i]=0;
                	}
                	
                	for(int i=offs; i<nOf; i++)
                	{
                		memTmp[i]=signal.s[i]*filter.b[0];
                
                		for(int j=1; j<filter.p; j++)
                		{
                			memTmp[i]+=signal.s[i-j]*filter.b[j];
                		}
                
                		for(int j=0; j<filter.q; j++)
                		{
                			memTmp[i]+=memTmp[i-j-1]*filter.a[j];
                		}
                	}
                
                	for(int i=nOf-orderOfFilter-1; i<nOf; i++)
                	{
                		s[i]=0;
                	}
                
                	for(int i=nOf-orderOfFilter-1; i>=0; i--)
                	{
                		s[i]=memTmp[i]*filter.b[0];
                
                		for(int j=1; j<filter.p; j++)
                		{
                			s[i]+=memTmp[i+j]*filter.b[j];
                		}
                
                		for(int j=0; j<filter.q; j++)
                		{
                			s[i]+=s[i+j+1]*filter.a[j];
                		}
                	}
                
                	return 0;
                }
                
                  0
                  Я думаю что если автор приведет пример того как работает (что это такое в данном случает) сепарабельность то станет понятно почему из 2d все вдруг стало 1d
                    +1
                    Берём изображение1.

                    Для каждого столбца: считаем строку как 1d.

                    Получаем изображение2. Обрабатываем уже его:

                    Для каждой строки: считаем столбец как 1d.

                    Получаем изображение3. Это и есть результат.

                    Я даже не знаю, как попроще объяснить…
                      0
                      Я всегда объясняю это так:

                      Если сначала в фотошопе сделать размытие по гаусу только вертикальное а потом размыть получившуюся картинку по вертикали то это будет тоже самое если применить фильтр гауса сразу (все с тем же размером кернела)
                  +2
                  Рядом «Честную»-то выложить бы не помешало!
                    0
                    Они похожи, на взгляд не отличить.
                      0

                      И тем не менее

                    0
                    Респект! Побольше бы таких статей на хабре.
                      0
                      Подобный метод применяется в детекторе границ Deriche. Довольно любопытный способ.
                        0
                        Кстати, в публикации, которую вы привели, есть ссылка на работу R. Deriche.
                        0
                        Напишите, пожалуйста, «гаусса» с большой буквы.
                        Спасибо.
                          +1
                          Одна из самых быстрых честных реализаций — через преобразование Фурье.
                          1. Делаем прямое преобразование изображения и фильтра (не только гауссового, но и любого другого)
                          2. Поэлементно перемножаем полученные спектры
                          3. С результатом делаем обратное преобразование Фурье
                          4. PROFIT

                          Сложность преобразования Фурье — O(nlog(n)), где n — количество пикселов в изображении
                          Время работы алгоритма не зависит от размера фильтра — 1 пиксел или же 1000.
                          Время работы — около 0.5 секунд для 12 мегапиксельной картинки на процессоре i7.
                            0
                            Предложенный алгоритм имеет сложность O(n) где n количество пикселей на изображении. Это меньше n log(n)

                            Кстати, если расскажите про алгоритм Блюстейна для БПФ произвольного размера, то тоже будет хорошо :-) Я, например, не очень понял.
                              0
                              Насколько я понял, предложенный алгоритм не относится к «честным», т.е. дает очень похожую, но не идентичную картинку. Это уже не гауссово размытие, а имитация гауссового размытия. В ряде областей это очень критично, например когда идет сложная обработка или восстановление изображений с участием гауссового фильтра.
                                0
                                Если честно, я не могу представить метода где нужна погрешность в расчёте гаусианы больше, чем в третьем — четвёртом знаке после запятой :-)
                              0
                              Маленькая ремарка: преобразование Фурье гауссовского распределения — тоже гауссовское, так что реально только надо преобразовывать саму картинку.
                                0
                                Ага, но еще нужно сообразить, какое у гауссианы будет \sigma в спектральной области, причем с учетом того, какой вариант записи DFT используется и в 95% случаев сделать это правильно не получается. Из тех, кому удалось, 4% отваливается на необходимости доопределить исходную картинку нулями, предварительно центрировав. Ну а необходимость определить ДПФ гауссианы с учетом хитрого порядка отсчетов при двумерном FFT доводит процент отвалившихся до 99.(9).
                              0
                              Отличнейшая статья! Очень пригодилось на практике. Благодаря ей удалось разобраться в методе довольно быстро. Но после реализации возникли вопросы:

                              А чему равна sigma при подсчёте коэффициентов? Она должна быть равна радиусу или половине радиуса? У меня получилось, что размытие вдвое сильнее, так как у нас два прохода.

                              Кроме того, во второй формуле указан знак минус, но на практике надо +.

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

                              Возможно, я что-то делаю неверно? Похоже мне не удаётся верно реализовать разворот между проходами на правой границе. Глазом это не заметно, но как0то не правильно это.
                                0
                                Впрочем, в процентном отношении около 0.5%, что не так уж много… видимо поэтому на глаз и не видно
                                  0
                                  Выложите excel файл, сейчас по виду понять сложно
                                    0
                                    Здесь уже немного другие входные данные, чем на скриншоте: максимум 255 в точке 10

                                    iir-gaussian.ods
                                      0
                                      Вы мне результат показали, я вижу что он кривоват.

                                      А вы мне покажите, как вы его получили.
                                        0
                                        Идея реализовать всё это в экселе частично помогла.

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

                                        1iir-gaussian2.ods

                                        Небольшой комментарий по реализации в экселе. В таблице есть три столбца. Первый соответствует входным данным, второй (G) первому проходу, а столбец J — второму и финальному проходу (справа налево). На первом проходе первые три точки считаются не так как все. На втором проходе последние три считаются не как все.

                                        Внизу и на втором листе есть график.
                                          0
                                          И да, всё это для радиуса 3.0 (соответственно, сигма была взята 1.5)
                                            0
                                            Ясно, что если бы расчёт был бы честным, то там нулей было бы куда больше. И если бы максимум был бы 1, а не 255, то относительная погрешность ещё бы выросла. Хотя сейчас, после исправления, разница между левой и правой частью составляет 0.01% от максимума. Возможно, это и есть предел точности для метода?
                                              0
                                              Во-первых, я не понял что значит «первые три отсчёта считаем не как все». Там формула та же, просто нужно ставить вместо отсчётов с отрицательными координатами нули.

                                              Во-вторых, вы взяли до радиуса 3.0 сигм?

                                              В целом, я думаю что поставив больше экспериментов (увеличив дальность расчётов, например) сами ответите на все вопросы :-) Математически всё должно получаться симметрично, зачит вопрос только в ошибках округления и насколько далеко считаем.
                                      0
                                      В функции gaussCoef есть переменная z0_im, значение которой не используется. Так задумано или это ошибка?

                                      Only users with full accounts can post comments. Log in, please.