LINUX.ORG.RU

Помогите разобраться что за фильтр и с аппроксимациями gaussian filter

 gaussian blur,


1

1

tl;dr

Посоны, я школу закончил давно, поэтому в математике теперь как кутузов - понимаю только очень крупным планом и неспешно. Помогите разобраться с наиболее быстрой «размывкой картинок» для больших радиусов фильтров.


ВОПРОС 1.

Вот тут http://incubator.quasimondo.com/processing/fast_blur_deluxe.php (там и на JS ссылка есть) перец заявил «StackBlur», который «типа лучше бокса и ближе к гауссу». А все остальные этот алгоритм юзают и не задумываются.

Но чем больше я смотрю на егоную математику, тем больше недоумеваю. Такое впечатление, что чувак решил сделать tent фильтр, и с какого-то хрена решил:

а) что это быстрее, чем 2 прохода box-фильтра
б) что tent-фильтр можно считать по частям (по вертикали и горизонтали) как прямоугольный и гауссовский

Кто-нибудь может проверить, он реально такую ересь слепил, или я чего-то по его коду не понимаю?


ВОПРОС 2.

Вот тут http://www.csse.uwa.edu.au/~pk/research/pkpapers/FastGaussianSmoothing.pdf описывается, что гауссовский фильтр можно считать быстрее, если крутить не один box-фильтр много раз, а сделать пару box-фильтров разной ширины.

Глубже я не копал, т.к. побоялся что взорвется мозг. Кто-нибудь может сказать по простому, будет ли профит по сравнению с 4-проходным box-фильтром на радиусах от 1 до 3?

Мне надо всего лишь картинки блюрить бысто, чтобы делать unsharp mask и local contrast enhancement.

★★★★★

Ответ на: комментарий от yoghurt

Я знаю как считаются бокс-фильтры, чай не совсем дебил. Но мне нужен гауссовский. И хочется понять, налажал чувак в StackBlur или нет.

Vit ★★★★★
() автор топика
Ответ на: комментарий от yoghurt

Я далек от обработки изображений, но тащем то любой линейный фильтр (в т.ч. и гаусса) ЕМНИП имеет подобную эффективную реализацию.

AIv ★★★★★
()

БПФ сейчас вполне шустро делаются. Я в реальном времени по десятку кадров в секунду обрабатывал: медианный фильтр + лапласиан гауссианы + поиск и выделение объектов + определение их центра тяжести по полученной маске.

Что тут медленного-то? А если нужно здоровенную картинку хитро как-то обработать (т.е. чтобы был профит от копирования туда-сюда кучи данных), то можно на GPU забульбенить.

Eddy_Em ☆☆☆☆☆
()
Ответ на: комментарий от Eddy_Em

Мамма миа! Он свертку циклами фигачит! И даже не попытался openMP туда воткнуть... Мде, хорош оптимизатор!!!

Eddy_Em ☆☆☆☆☆
()
Ответ на: комментарий от AIv

Примеры в студию. Как апроксимировать гаусса произвольного радиуса быстрее чем 4-6 проходами бокса.

Vit ★★★★★
() автор топика
Ответ на: комментарий от Vit

А почему ты БПФ не хочешь использовать? Рисуешь гаусса нужного размера, получаешь его фурье-образ, умножаешь на фурье-образ изображения, делаешь обратное БПФ → вуаля! И никакого извращения. И быстро.

Eddy_Em ☆☆☆☆☆
()
Ответ на: комментарий от Vit

Может у Eddy_Em есть ссылкив печатном виде. Я просто для линейного фильтра с косинусом делал такую реализацию. В 2D случае чуток сложнее но тоже реализуемо, лет дцать назад у нас студент на эту тему диплом защищал (делал таким образом вейвлеты). Там сложность О(N)

AIv ★★★★★
()
Ответ на: комментарий от AIv

Я до конца не уверен, но в 2D случае для линейного фильтра скорей всего расщепление работает (т.е. сначала по одной координате прогнали, потом по другой).

ЗЫ Если конечно ядро как произвдение двух 1D функций.

AIv ★★★★★
()
Последнее исправление: AIv (всего исправлений: 1)
Ответ на: комментарий от Vit

В общем в 1D там следующее значение отфильтрованного изображения считается как (предыдущее_значение - значение_с_заднего_фронта*A)*B + значение_с_переднего_фронта*C.

Если расписать отфильтрованное значение как сумму, вычесть из него следующее отфильтрованное значение и перегруппировать члены это видно.

AIv ★★★★★
()
Ответ на: комментарий от AIv

Или чуток посложнее. Математику эту вспоминать щас лень, вот мой древний код для линейного фильтра с косинусом. Для Гаусса по идее должно быть что то похожее (ЕМНИП это подходит для любого линейного фильтра):

#include <iostream.h>
#include <math.h>
#include <list>
using namespace std;

int main(int NArg,char** Arg){
        int l=1;
        if(NArg>1) l=atoi(Arg[1]);
        else{
                char buf[1024];
                cin.getline(buf,1024);
                while(!!cin){
                        cout<<buf<<endl;
                        cin.getline(buf,1024);
                }
                return 0;
        }

        list<double> history;
        double x,y,x0,y0;
        cin>>x0>>y0>>x>>y;
        history.push_back(y0);

        double 
                delta=x-x0,
                L=(.5+l)*delta,
                C=.5*delta/L,
                cosPiL=cos(M_PI*delta/L),
                sinPiL=sin(M_PI*delta/L),
                Acos=L*sinPiL/(M_PI*delta),
                Asin=L*(1.-cosPiL)/(M_PI*delta),
                S=y0, 
                Scos=(-sinPiL*Asin-cosPiL*Acos)*y0, 
                Ssin=(-cosPiL*Asin+sinPiL*Acos)*y0, 
                addcos, addsin;
        y0=0.;

        while(!!cin){
                history.push_back(y);
                S+=y-y0;
                addcos=Scos-Acos*(y-y0);
                addsin=Ssin-Asin*(y-y0);
                Scos=cosPiL*addcos+sinPiL*addsin;
                Ssin=cosPiL*addsin-sinPiL*addcos;
                if(history.size()>=2*l+1){
                        cout<<x-L+.5*delta<<" "<<(S+Scos)*C<<endl;
                        y0=history.front();
                        history.pop_front();
                };
                cin>>x>>y;
        }
        return 0;
}

AIv ★★★★★
()
Ответ на: комментарий от AIv

Там же надо сумму весов учитывать. Вычитать и добавлять с краев прокатит только когда все веса одинаковые, то есть box filter. Ну можно еще для ядро елкой (1 2 3 4 3 2 1), если интегральные суммы заранее посчитать. Но подсчет интегральных сумм займет столько же, сколько второй проход бокса. К тому же елка

Где почитать можно про умножание на эти магические коэффициент из вашей формулы?

Vit ★★★★★
() автор топика
Ответ на: комментарий от AIv

Я тут статьи пытался осилить по оптимизацим фильтров, в одной было написано прямым текстом - tent на раскладывается.

Vit ★★★★★
() автор топика
Ответ на: комментарий от Vit

Любой линейный фильтр наск я помню так реализуется. Коэффициенты как раз от весов и зависят. Я же привел пример с ядром в виде косинуса - он довльно похож на Гаусса и работает;-)

AIv ★★★★★
()

Марио сознался. Сказал что действительно пирамиду считает и что сам не уверен что пирамидка поддается разбивке, но так как визуальный эффект хороший, то особо не заморачивался.

Vit ★★★★★
() автор топика
Ответ на: комментарий от Vit

В примере линейный фильтр с косинусом в качестве ядра. Масштабирование там тоже есть конечно:-)

AIv ★★★★★
()
Ответ на: комментарий от Vit

там только конечно просто <iostream> и <stdlib.h> не хватает.

$ python -c 'for i in range(43): print i*.1, 1 if i==21 else 0' | ./linef-cos 10 
1 0
1.1 0.000707305
1.2 0.00487562
1.3 0.0128419
1.4 0.0238982
1.5 0.0370623
1.6 0.0511644
1.7 0.0649514
1.8 0.0771984
1.9 0.0868172
2 0.092953
2.1 0.0950607
2.2 0.092953
2.3 0.0868172
2.4 0.0771984
2.5 0.0649514
2.6 0.0511644
2.7 0.0370623
2.8 0.0238982
2.9 0.0128419
3 0.00487562
3.1 0.000707305
3.2 -5.3686e-17

AIv ★★★★★
()
Ответ на: комментарий от AIv

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

Ну то есть код-то такой написать можно, только он по-моему будет фигню выдавать :)

Vit ★★★★★
() автор топика

Это все хорошо, но ты должен понимать, что на практике правильный паттерн доступа к памяти часто важнее теоретической сложности. Если размер фильра не очень большой(до 10 пикселей), то просто считай двумя проходами, если большой, то используй box-blur, или как это называется infinite impulse response approximations.

В плане производительности очень поможет одновременно с горизонтальным проходом транспонировать изображение, чтобы второй вертикальный проход был на самом деле горизонтальным, так ты получишь выигрышь в векторном считывании памяти. Здесь все описано https://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementa... .

В реальности, все конечно используют БПФ, тем более для гауссианы фурье-образ можно считать аналитически на лету, не теребя лишний раз память.

maggotroot
()
Ответ на: комментарий от Vit

Ну вот я же привел пример кода, и он выдает ровно то что нужно же. И да, берет только крайние точки.

AIv ★★★★★
()
Ответ на: комментарий от maggotroot

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

За ссылку спасибо большое, очень годная.

Если не сложно, объясните пожалуйста «на пальцах» два момента:

- зачем делается 2 прохода, справа налево и слева направо?
- как в таких фильтрах поступают с краями картинки, когда точек перестает хватать?

Vit ★★★★★
() автор топика
Ответ на: комментарий от Vit

- зачем делается 2 прохода, справа налево и слева направо?

Таким проходом можно аппроксимровать ядро определённое только при x >= 0. Если, например, двигаться только слева на паво, то каждое R[x] значение результата будет зависеть от U[y] y <= x исходных значений. Чтобы аппроксимировать всю неограниченную гауссиану на (-Inf, +Inf) приходится делать два прохода каждый из которых аппроксимирует по половине ядра, при x >=0 и x <= 0.

как в таких фильтрах поступают с краями картинки, когда точек перестает хватать?

Так же как и везде, т.е. зависит от «вкуса», ~ продолжают краевым значением, зеркально продолжают значения за границы, различные экспраполяции приграничных значений.

mashina ★★★★★
()
Последнее исправление: mashina (всего исправлений: 1)
Ответ на: комментарий от mashina

Ну и ещё есть проблема симметричности если пытаться делать за один проход.

mashina ★★★★★
()
Ответ на: комментарий от Vit

В ява скрипте есть хоть какая-то векторизация? Чтобы быстро умножить 2 массива поэлементно, например как x.*y в матлабе? Если есть то точно делай фурье, если нету, то скорее всего тоже фурье будет лучше :) Всякие аппроксимации очень хороши только за счет умного паттерна доступа к памяти.

ПС. Можно же как-то webgl заюзать?

maggotroot
()
Ответ на: комментарий от maggotroot

Если есть то точно делай фурье, если нету, то скорее всего тоже фурье будет лучше :)

Чтобы оно было точно лучше придётся изрядно потрахать мозги, т.к. простой blur через преобразование Фурье не эквивалентен аналогичному через приближения: ~ Фурье неявно задаёт схему продления картинки за краями методом бесконечного тиражирования изображения с каждой из сторон. Это плохой способ продолжения картинки за границы т.к. делает разрывы с соответствующими артефактам (~ слева у границы 0, справа максимум - результат будет содержать нежелательные затемнение слева и осветление справа). Решение проблемы с граничными условиями для Фурье сложнее.

Сложность Фурье выше FSM, O(N) vs примерно O(NlogN). Есть смысл использовать Фурье если не удаётся придумать FSM для ядра. Например, если n-мерное ядро нельзя разделить на n одномерных. Но GB отлично реализуется на FSM.

mashina ★★★★★
()
Ответ на: комментарий от mashina

Про сложность ты, очевидно, прав, но если есть способ быстрого поэлементного умножения, то фурье будет быстрее, так как для FIR придется писать поэлементный цикл на JS, что будет работать дико медленно.

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

maggotroot
()
Ответ на: комментарий от maggotroot

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

WebGL заюзать можно, но эта штука есть далеко не везде, поэтому тема выделена в отдельную задачу . И тут даже денег за нее предлагал - не по моим знаниям в подобное вписываться (мне не лень, просто еще другие проекты веду).

Vit ★★★★★
() автор топика
Ответ на: комментарий от mashina

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

Vit ★★★★★
() автор топика
Ответ на: комментарий от Vit

Для примера - у ресайзера входные данные RGBA лежат в Uint8Array. Я пробовал сделать Uint32Array, который мапится на тот же буфер, чтобы по 4 байта за раз выгребать. Стало хуже. Вот так и живем :)

Vit ★★★★★
() автор топика
Ответ на: комментарий от Vit

Я мало что знаю про JS, и извини, что увожу тему, но ведь вроде в CSS есть разные фильтры для картинок? Вроде можно с помощью то ли CSS то ли с SVG наложить 2 картинки и получить результат? Да там видимо 8-битная арифметика, но может можно как-то справиться?

maggotroot
()
Ответ на: комментарий от maggotroot

CSS «write only». Фильтры наложить можно, но извлечь результат обратно не получится. Реально либо полностью делать на webgl с векторами, либо на яваскрипте без векторов. Сейчас ходят слухи, что добавят SIMD, но это не универсально и пока только слухи. Универсальный вариант делать все равно прийдется. А всякие изыски - это уже на потом.

Для меня твоя ссылка на интел полезна тем, что теоретически позволит добиться такой же скорости, что двойной box blur. Теоретически - потому что проходов по памяти столько же, но я не знаю, на сколько яваскриптовая плавучка всё просадит. Проверять надо.

Но в принципе жить можно. Практика показывает, что хорошо заточенный JS уступает С всего в полтора раза. Конечно, если не считать SSE и с некоторыми ограничениями по арифметике.

Vit ★★★★★
() автор топика
Ответ на: комментарий от maggotroot

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

Ядро в данном случае это такая же картинка по размеру. Костылить с заграничной областью, конечно, можно, но на неё накладываются точно такие же условия - т.е. не должно быть нигде разрывов в т.ч. и на границах между продолжениями. В 2D это уже сделать не просто. Но ещё нужно думать о прочих геометрических искажениях которые будет накладывать цикличность и продолжение после обратного преобразования. В общем случае совсем избавиться от «цикличных» искажений от DFT нельзя.

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

Это является проблемой всегда когда пытаются делать операции с нециклическими свёртками через DFT.

mashina ★★★★★
()
Последнее исправление: mashina (всего исправлений: 1)

Что думают всякие немецкие демосценеры

Почитай тут: http://fgiesen.wordpress.com/2012/07/30/fast-blurs-1/ http://fgiesen.wordpress.com/2012/08/01/fast-blurs-2/

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

nikitos ★★★
()
Ответ на: Что думают всякие немецкие демосценеры от nikitos

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

Vit ★★★★★
() автор топика
Ответ на: комментарий от mashina

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

maggotroot
()
Ответ на: комментарий от maggotroot

А есть где-нибудь описание простых принципов, как эффективно использовать кеш? Ну помимо транспарирования картинки на вертикальном проходе.

У меня картинка в типизированном массиве. На практике это обычный кусок памяти, только обращение не по указателю а по индексу от начала. Из операций:

- либо дергаем элементы в пределах строки (как в интеловской статье)
- либо дергаем элементы в пределах строки + дергаем из таблицы с предрасчитанным фильтром (для ресайза, ланцеш, общая длина = строка*3)
- еще наверное JS будет отписывать в память локальные переменные, которые в регистры не влезут.

Про L1 и L2 чего-то там «краями ухов» слышал, но деталей не представляю.

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

Vit ★★★★★
() автор топика
Ответ на: комментарий от Vit

Общие рекомендации типа этой прочитал.

Интересуют нюансы под конкретную задачу - стоит ли нарезать большие картинки на прямоугольники и т.п.

Vit ★★★★★
() автор топика
Ответ на: комментарий от maggotroot

Вот:

http://nodeca.github.io/glur/

Подбил комрада, с которым делали порт zlib на яваскрипт, зафигачить имплементацию примерно по интеловской статье. Правда он более оптимальную апроксимацию IIR откопал.

Еще не оптимизировали скорость до конца, но вцелом годно. Мне для unsharp mask нужно блюрить только по одному каналу (светимость), так что скорость точно проблемой не станет.

Vit ★★★★★
() автор топика
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.