Box-фильтры с большим размером окна эффективно считаются с помощью интегрального изображения.

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

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

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

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

По первой ссылке какой-то хренов плагин хочет. На кой оно мне надо?

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

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

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

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

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

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

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

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

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

#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;
}

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

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

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

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

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

То что в примере - всего лишь масштабирование синуса на ширину окна.

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

там только конечно просто <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

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

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

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

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

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

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

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

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

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

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

Я на яваскрипте пишу

please, make me unsee it!

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

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

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

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

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

Спасибо!

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Разве в JS кэш вообще важен? Проходи массив блоками.

Вот:

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

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

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