Про систему блюра – Laplace Blur

Можна блюрить Лапласом замість Гаусса, у скільки разів це швидше, і варто того втрата 1/32 точності.

(Laplace Blur — Пропонована оригінальна назва алгоритму)

Сьогодні мій внутрішній демосценер штовхнув мене і змусив таки написати статтю, яку потрібно було написати вже півроку тому. Як любитель на дозвіллі розробляти оригінальні алгоритми ефектів, хотів би запропонувати громадськості алгоритм «майже гаусиан блюра», що відрізняється застосуванням виключно швидких процесорних інструкцій (зрушень та масок), а тому доступний до реалізації аж до мікроконтролерів (надзвичайно швидкий в обмеженому оточенні).

За сформованою у мене традиції написання статей на хабр, я буду наводити приклади на JS, як мовою максимально популярному, і хочете вірте хочете ні, дуже зручному для цілей швидкого прототипування алгоритмів. До того ж, можливість реалізувати це ефективно на JS прийшла разом з типізованими масивами. На моєму не дуже потужний ноут фулскрин зображення обробляється зі швидкістю 30fps (багатопоточність воркеров задіяна не була).

Дисклеймер для крутих математиківВідразу ж скажу, що я знімаю перед вами капелюха, тому що сам вважаю себе не досить підкованим фундаментальної математики. Однак я завжди керуюся загальним духом фундаментального підходу. Тому перед тим як хейтить мій, деякою мірою «наглядова» підхід до апроксимації, озаботьтесь підрахунком бітової складності алгоритму, який, як ви думаєте, може бути отриманий методами класичної апроксимації поліномами. Я вгадав? Ними ви хотіли швидко апроксимувати? З урахуванням того, що вони зажадають плаваючої арифметики, вони будуть значно повільніше, ніж один бітовий зсув, до якого я наведу пояснення в підсумку. Одним словом, не поспішайте з теоретичним фундаменталізмом, і не забувайте про контекст в якому я вирішую завдання.

Даний опис присутній тут скоріше для того, щоб пояснити хід моїх думок і здогадів, які привели мене до результату. Для тих кому буде цікаво:

Оригінальна Функція Гауса:

g(x) = a * e ** ( — ((x-b)**2) / c), де
а — амплітуда (якщо колір у нас вісім біт на канал, то вона = 256)
e — постійна Ейлера ~2.7
b — зсув графіка по х (нам не потрібен = 0)
з — параметр впливає на ширину графіка пов’язаний з нею як ~w/2.35

Наша приватна функція (мінус показника ступеня прибраний з заміною множення діленням):

g(x) = 256 / e ** ( x*x / c)

Читайте також  Трагікомедія в NaN актах: як розробники зробили гру на JS і випустили її в Steam

Так розпочнеться брудне аппроксимационное дійство:

Зауважимо, що параметр c — дуже близький до полушіріна і поставимо 8 (це пов’язане з тим, за скільки кроків можна зрушити по одному 8 біт каналу).

Також грубо замінимо на e 2, зазначивши однак, що це впливає більше на кривизну «дзвони» ніж на його межі. Насправді впливає на 2/e раз, але сюрприз в тому, що ця помилка компенсує внесену при визначенні параметра с, так що граничні умови у нас все ще в порядку, а похибка з’являється лише злегка некоректному «нормальному розподілі», для графічних алгоритмів це буде позначатися на динаміці градієнтних переходів кольорів, але це майже неможливо помітити оком.

Отже, тепер наша функція така:

gg(x) = 256 / 2 ** ( x*x / 8) або gg(x) = 2 ** (8 — x*x / 8)
Зауважимо, що показник ступеня (x*x/8) має ту ж область значення [0-8], що і функція меншого порядку abs(x), тому остання є кандидатом на заміну. Швидко перевіримо гіпотезу, глянувши як змінюється графік при ній gg(x) = 256/ (2 ** abs(x)):

GaussBlur vs LaplasBlur:

Здається, відхилення занадто великі, до того ж функція, втративши гладкість, тепер має пік. Але стривайте.
По-перше, не будемо забувати, що гладкість одержуваних при розмитті градієнтів залежить не від функції щільності ймовірності (якою є функція Гауса), а від її інтеграла функції розподілу. На той момент я цього факту не знав, але насправді, провівши «згубну» апроксимацію щодо функції щільності ймовірності (Гауса), функція розподілу при цьому залишилася досить схожою.

Було:

Стало:

Пруф, знятий з готового алгоритму, збігається:

(забігаючи вперед скажу, що помилка розмиття мого алгоритму щодо Gausian x5 склала всього 3%)
Отже, ми перейшли набагато ближче до функції розподілу Лапласа. Хто б міг подумати, але їм можна мити зображення на 97% не гірше.

Пруф, різниці гаусиан блюра х5 і «Лаплас блюра» х7:

(це не чорне зображення! Можете повивчати в редакторі)

Допущення цього перетворення дозволило перейти до ідеї отримання значення ітеративної фільтрацією, до якого я планував звести спочатку.

Перед розповіддю конкретного алгоритму буде чесно, якщо я забіжу вперед і відразу опишу його єдиний недолік (хоча реалізацію можна виправити, з втратою швидкості). Але цей алгоритм реалізований з використанням сдвиговой арифметики, і ступеня 2ки є його обмеженням. Так оригінальний зроблений для розмиття х7 (що в тестах ближче всього співвідноситися з Gausian x5). Пов’язано це обмеження реалізації з тим, що при восьмибитном кольорі, зрушуючи за крок значення в накопичувачі фільтра на один біт, всяке вплив від точки закінчується максимум за 8 кроків. Я також реалізував трохи більш повільну версію через пропорції і додаткові додавання, яка реалізує швидке поділ на 1.5 (отримавши в підсумку радіус х15). Але з подальшим застосуванням такого підходу похибка зростає, а швидкість падає, що не дозволяє його так використовувати. З іншого боку, варто зауважити, що х15 вже достатньо, що б не сильно помічати різницю, отриманий результат з оригіналу або із даунсемплированого изоображения. Так що метод цілком придатний, якщо потрібна надзвичайна швидкість в обмеженому оточенні.

Отже, ядро алгоритму просте, здійснюється чотири однотипних проходи:

1. Половина значення накопичувача t (спочатку рівне нулю) складається з половиною значення чергового пікселя, результат присвоюється нього ж. Продовжуємо так до кінця рядка зображення. Для всіх рядків.

Читайте також  Повнолітня журналістика: від Росії до Кремля

По завершенню першого проходу зображення виявляється змазано в одну сторону.

2. Другим проходом робимо те ж саме у зворотний бік для всіх ліній.
Отримуємо зображення, повністю розмите по горизонталі.

3-4. Тепер робимо те ж саме по вертикалі.
Готово!

Спочатку я використовував двухпроходний алгоритм з реалізацією зворотного розмиття через стек, але він складний для розуміння, не граціозний, і до того ж на нинішніх архітектурах виявився повільніше. Можливо, однопрохідний алгоритм буде швидше на мікроконтролерах, плюсом так само стане можливість виводити результат прогресивно.
Поточний четырехпроходный спосіб реалізації я підглянув на хабре у попереднього гуру за алгоритмами розмиття. habr.com/post/151157 Користуючись нагодою, висловлюю йому свою солідарність та глибоку вдячність.

Але на цьому хакі не закінчилися. Тепер з приводу того, як обчислювати всі три канали кольорів за одну інструкцію процесора! Справа в тому, що використаний в якості ділення на два бітовий зсув дозволяє дуже добре контролювати позиції біт результату. Єдиною проблемою стає те, що молодші розряди каналів сповзають в сусідні старші, але можна просто очистити їх, ніж виправити проблему, з деякою втратою точності. А згідно описаної формулою фільтра, складання половини значення накопичувача з половиною черговий значення комірки (за умови обнулення съехавших розрядів) ніколи не призводить до переповнення, так що турбуватися за нього не варто. І формула фільтра для одночасного обчислення всіх розрядів ставати такою:

buf32[i] = t = ( ( ( t >> 1 ) & 0x7F7F7F ) + ( ( buf32[i] >> 1 ) & 0x7F7F7F );

Однак потрібно ще одне доповнення: досвідченим шляхом було з’ясовано, що втрата точності в даній формулі надто значна, візуально відчутно стрибає яскравість картинки. Стало зрозуміло, що теряемый біт потрібно округляти до цілого, а не відкидати. Простим способом зробити це в цілочисельної арифметики є збільшення половини дільника перед поділом. Дільник у нас двійка, значить потрібно додавати одиницю, у всі розряди, — константу 0x010101. Але з будь-яким складанням потрібно побоюватися отримати переповнення. Так що ми не можемо використовувати таку корекцію для обчислення половинного значення черговий осередку. (Якщо там білий колір, ми отримаємо переповнення, тому її ми коригувати не будемо). Але виявилося, що основну помилку вносило багаторазове поділ накопичувача, який ми як раз скорегувати можемо. Тому що, насправді, навіть з такою корекцією, значення в накопичувачі не підніметься вище 254. Зате при складанні з 0x010101 переповнення гарантовано не буде. І формула фільтра з корекцією набуває вигляду:

Читайте також  Несподівана зустріч. Глава 18 [остання глава, вихід на краудфандінг]

buf32[i] = t = ( ( ( ( 0x010101 + t ) >> 1 ) & 0x7F7F7F ) + ( ( buf32[i] >> 1 ) & 0x7F7F7F );

Насправді формула досить добре виконує корекцію, так що при повторному багаторазовому застосуванні даного алгоритму до зображення, артефакти починають виднітися тільки на другому десятку проходів. (не факт що повторення гаусиан блюра не дасть подібних артефактів).

Додатково, є чудова властивість, при безлічі проходів. (Є заслугою не мого алгоритму, а самої «нормальності» нормального розподілу). Вже на другому проході Лапласа блюра функція щільності ймовірності (якщо я правильно все прикинув) буде виглядати якось так:

Що, погодьтеся, дуже близько до гаусиана.

Досвідченим шляхом я знайшов, що використовувати модифікації з великим радіусом допустимо в парі, т. к. описане вище властивість компенсує помилки, якщо останній прохід є більш точним (найточніший це описаний тут алгоритм х7 блюра).

демо: impfromliga.github.io/LaplaceBlur
реп: github.com/impfromliga/LaplaceBlur
кодпен: codepen.io/impfromliga/pen/brXZjJ

Для крутих математиків:

Що було б цікаво дізнатися, наскільки коректно використовувати фільтр сепарабельно, не впевнений виходить там симетрична картина розподілу. Хоча на око неоднорідностей і не видно.

Степан Лютий

Обожнюю технології в сучасному світі. Хоча частенько і замислююся над тим, як далеко вони нас заведуть. Не те, щоб я прям і знаюся на ядрах, пікселях, коллайдерах і інших парсеках. Просто приходжу в захват від того, що може в творчому пориві вигадати людський розум.

You may also like...

Залишити відповідь

Ваша e-mail адреса не оприлюднюватиметься. Обов’язкові поля позначені *