Відновлення розфокусованих і змазаних зображень

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

Чому ж для усунення смаза та розфокусування практично нічого немає (unsharp mask не в рахунок) – може бути це в принципі неможливо? Насправді можливо – відповідний математичний апарат почав розроблятися приблизно 70 років тому, але, як і для багатьох інших алгоритмів обробки зображень, все це знайшло широке застосування тільки в недавній час. Ось, в якості демонстрації вау-ефекту, пара картинок:

Я не став використовувати замучену Лену, а знайшов свою фотку Венеції. Праве зображення чесно отримано з лівого, причому без використання хитрощів типу 48-бітового формату (в цьому випадку буде 100% відновлення вихідного зображення) – зліва самий звичайний PNG, розмитий штучно. Результат вражає… але на практиці не все так просто. Під катом детальний огляд теорії і практичні результати.
Обережно, багато картинок у форматі PNG!

Введення

Почнемо здалеку. Багато хто вважає, що розмиття необоротна операція і інформація безповоротно втрачається, оскільки кожен піксель перетворюється в пляму, все змішується, а при великому радіусі розмиття так і зовсім отримаємо однорідний колір по всьому зображенню. Це не зовсім так – вся інформація просто перерозподіляється по деякому закону і не може бути однозначно відновлено з деякими застереженнями. Виняток становить лише краю зображення шириною в радіус розмиття – там повноцінне відновлення неможливо.

Продемонструємо це «на пальцях», використовуючи невеликий приклад для одновимірного випадку – уявімо, що у нас є ряд з пікселів зі значеннями:
x1 | x2 | x3 | x4… – Початкове зображення

Після спотворення значення кожного пікселя підсумовується зі значенням лівого, тобто ‘ i = xi + xi-1. По ідеї, треба поділити на 2, але опустимо це для простоти. В результаті маємо розмите зображення зі значеннями пікселів:
x1 + x0 | x2 + x1 | x3 + x2 | x4 + x3… – Розмите зображення

Тепер будемо пробувати відновлювати, віднімемо послідовно по ланцюжку значення за схемою – з другого пікселя перший, із третього результат другого, четвертого результат третього і так далі, одержимо:
x1 + x0 | x2 — x0 | x3 + x0 | x4 — x0… – Відновлене зображення

У підсумку замість розмитого зображення отримали початкове зображення, до пікселям якого додана невідома константа x0 з чергується знаком. Це вже набагато краще – цю константу можна підібрати візуально, можна припустити, що вона приблизно дорівнює значенню x1, можна автоматично підібрати з таким критерієм, щоб значення сусідніх пікселів «скакали» як можна менше і т. д. Але все міняється, як тільки ми додаємо шум (які завжди є в реальних зображеннях). За описаною схемою на кожному кроці буде накопичуватися внесок шуму в загальну складову, що в підсумку може дати абсолютно неприйнятний результат, але, як ми переконалися, відновлення цілком реально навіть таким примітивним способом.

Модель процесу спотворення

А тепер перейдемо до більш формального і наукового опису цих процесів спотворення і відновлення. Будемо розглядати тільки напівтонові чорно-білі зображення в припущенні, що для обробки повнокольорового зображення досить повторити всі необхідні кроки для кожного з каналів кольорів RGB. Введемо наступні позначення:
f(x, y) – вихідне неспотворене зображення
h(x, y) – функція спотворює
n(x, y) – адитивний шум
g(x, y) – результат спотворення, тобто те, що ми спостерігаємо в результаті (змазане або розфокусований зображення)

Сформулюємо модель процесу спотворення наступним чином:
g(x, y) = h(x, y) * f(x, y) + n(x, y) (1)

Завдання відновлення спотвореного зображення полягає в знаходженні найкращого наближення f'(x, y) вихідного зображення. Розглянемо кожну складову більш докладно. З f(x, y) та g(x, y) все досить зрозуміло. А ось про функцію h(x, y) потрібно сказати пару слів – що ж вона з себе представляє? В процесі спотворення кожен піксель вихідного зображення перетворюється в пляму для випадку розфокусування і відрізок для випадку простого смаза. Або ж можна сказати навпаки, що кожен піксель спотвореного зображення «збирається» з пікселів деякій околиці вихідного зображення. Все це один на одного накладається і в результаті ми отримуємо спотворене зображення. Те, за яким законом розмазується або збирається один піксель і називається функцією спотворення. Інші синоніми – PSF (Point spread function, тобто функція розподілу точки), ядро спотворює оператора, kernel і інші. Розмірність цієї функції, як правило менше розмірності самого зображення – наприклад, у початковому розгляді прикладу «на пальцях» розмірність функції була 2, оскільки кожен піксель складався з двох.

Спотворюючі функції

Подивимося, як виглядають типові спотворюють функції. Тут і далі будемо використовувати, що вже став стандартним для таких цілей інструмент – Matlab, він містить у собі все необхідне для найрізноманітніших експериментів з обробкою зображень (і не тільки) і дозволяє зосередитися на самих алгоритмах, перекладаючи всю рутинну роботу бібліотеки функцій. Втім, за це доводиться розплачуватися продуктивністю. Отже, повернемося до PSF, ось приклади їх види:

PSF у разі розмиття по Гауссу функцією fspecial(‘gaussian’, 30, 8);

PSF у разі смаза фунцією fspecial(‘motion’, 40, 45);

Операція застосування спотворює функції до іншої функції (до зображення, в даному випадку) називається згорткою (convolution), тобто деяка область початкового зображення згортається в один піксель спотвореного зображення. Позначається через оператор «*», не плутати зі звичайним множенням! Математично для зображення f с розмірів M x N і спотворює функції h c розмірів m x n-це записується так:

(2)

Де a = (m — 1) / 2, b = (n – 1) / 2. Операція, зворотна згортку, називається деконволюцией (deconvolution) і рішення такої задачі досить нетривіально.

Модель шуму

Залишилося розглянути останнє доданок, що відповідає за шум, n(x, y) у формулі (1). Причини шуму у цифрових сенсорах можуть бути самими різними, але основні це – теплові коливання і темновые струми. На величину шуму також впливає ряд факторів, таких як значення ISO, тип матриці, розмір пікселя, температура, електромагнітні наведення та ін. У більшості випадків шум є Гауссовим (який задається двома параметрами – середнім і дисперсією), а також є адитивним, не корелює з зображенням і не залежить координат пікселів. Останні три припущення є дуже важливими для подальшої роботи.

Теорема про згортку

Повернемося тепер до первісної постановки задачі відновлення – нам необхідно якимось чином звернути згортку, при цьому не забуваючи про шум. З формули (2) видно, що f(x, y) з g(x, y) не так-то просто – якщо вирішувати, що називається, «в лоб», то вийде величезна система рівнянь. Але на допомогу приходить перетворення Фур’є, не будемо на ньому зупинятися, по цій темі вже було сказано чимало. Так от, є така теорема про згортку, яка говорить, що операція згортки в просторовій області еквівалентна звичайного множення в частотній області (причому поелементне множення, а не матричне). Відповідно, зворотна операція згортку еквівалентна поділу в частотній області, тобто це можна записати як:
(3)

Де H(u, v), F(u, v) – Фур’є-образи відповідних функцій. Значить процес спотворення з формули (1) можна переписати в частотній області:
(4)

Інверсна фільтрація
Тут же напрошується поділити цю рівність на H(u, v) і отримати таку оцінку F^(u, v) вихідного зображення:
(5)
Це називається інверсної фільтрацією, але на практиці майже ніколи не працює. Чому ж? Щоб відповісти на це питання, подивимося на останній доданок у формулі (5) – якщо функція H(u, v) приймає значення, близькі до нуля або нульові, то внесок цього доданку буде домінуючим. Це практично завжди зустрічається в реальних прикладах – для пояснення цього згадаємо як виглядає після спектр перетворення Фур’є.

Беремо початкове зображення,

перетворимо його в півтонове і, використовуючи Matlab, отримуємо спектр:

% Load image
I = imread('image_src.png');
figure(1); imshow(I); title('Початкове зображення'); 
% Convert image into grayscale
I = rgb2gray(I);
% Compute Fourier Transform and it center
fftRes = fftshift(fft2(I));
% Show result
figure(2); imshow(mat2gray(log(1+abs(fftRes)))); title('FFT - Амплітудний спектр (логарифмічна шкала)');
figure(3); imshow(mat2gray(angle(fftRes))); title('FFT - Фазовий спектр');

В результаті отримуємо дві компоненти: амплітудний і фазовий спектри. Про фазу, до речі, багато хто забуває. Зверніть увагу, що амплітудний спектр показаний у логарифмічній шкалі, оскільки його значення варіюються дуже сильно – на кілька порядків, в центрі максимальні значення (близько мільйона) і швидко зменшуються практично до нульових по мірі віддалення від центру. Саме з-за цього інверсна фільтрація буде працювати тільки при нульових або практично нульових значеннях шуму. Продемонструємо це на практиці з допомогою такого скрипта:

% Load image
I = im2double(imread('image_src.png')); 
figure(1); imshow(I); title('Початкове зображення');
% Blur image
Blurred = imfilter(I, PSF,'circular','conv' );
figure(2); imshow(Blurred); title('Розмите зображення');
% Add noise
noise_mean = 0;
noise_var = 0.0;
Blurred = imnoise(Blurred, 'gaussian', noise_mean, noise_var);
% Deconvolution
figure(3); imshow(deconvwnr(Blurred, PSF, 0)); title('Результат');

 noise_var = 0.0000001 noise_var = 0.000005

Добре видно, що додавання навіть дуже невеликого шуму призводить до значних перешкод, що сильно обмежує практичне застосування методу.

Існуючі підходи для деконволюції
Але є підходи, які враховують враховують наявність шуму на зображенні – один з найвідоміших і найбільш перше, це фільтр Вінера (Вінер). Він розглядає зображення і шум як випадкові процеси і знаходить таку оцінку f’ для неискаженного зображення f, щоб середньоквадратичне відхилення цих величин було мінімальним. Мінімум цього відхилення досягається на функцію в частотній області:
(6)
Цей результат був получине Вінером в 1942 році. Детальний висновок тут приводити не будемо, ті, хто цікавиться, можуть подивитися його тут . Функцією S тут позначаються енергетичні спектри шуму та вихідного зображення відповідно – оскільки, ці величини рідко бувають відомі, то дріб Sn / Sf замінюють на деяку константу K, яку можна приблизно охарактеризувати як співвідношення сигнал-шум.

Наступний метод, це «сглаживающая фільтрація методом найменших квадратів зі зв’язком», інші назви: «фільтрація по Тихонову», «Тихоновська регуляризация». Його ідея полягає у формулюванні задачі в матричному вигляді з подальшим рішенням відповідної задачі оптимізації. Це рішення записується у вигляді:
(7)
Де y – параметр регуляризації, а P(u, v) – Фур’є-перетворення оператора Лапласа (матриці 3 * 3).

Ще один цікавий підхід запропонували незалежно Ричардосн [Richardson, 1972] і Люсі [Lucy, 1974]. Метод так і називається «метод Люсі-Річардсона». Його особливість у тому, що він є нелінійним, на відміну від перших трьох – що потенційно може дати кращий результат. Друга особливість – метод є ітераційним, відповідно виникають труднощі з критерієм зупинки ітерацій. Основна ідея полягає у використанні методу максимальної правдоподібності для якого передбачається, що зображення підпорядковується розподілу Пуассона. Формули для обчислення досить прості, без використання перетворення Фур’є – все робиться в просторовій області:
(8)
Тут символом «*», як і раніше, позначається операція згортки. Цей метод широко використовується в програмах для обробки астрономічних фотографій – у них використання деконволюції (замість unsharp mask, як в фоторедакторі) є стандартом де-факто. В якості прикладу можна навести Astra Image, ось приклади деконволюції. Обчислювальна складність методу дуже велика – обробка середньої фотографії, залежно від кількості ітерацій, може знанимать багато години і навіть дні.

Останній розглянутий метод, а вірніше, ціле сімейство методів, які зараз активно розробляються і розвиваються – це сліпа деконволюция (blind deconvolution). У всіх попередніх методах передбачалося, що спотворює функція PSF точно відома, в реальності це не так, зазвичай PSF відома лише приблизно за характером видимих спотворень. Сліпа деконволюция як раз є спробою враховувати це. Принцип досить простий, якщо не заглиблюватися в деталі – вибирається перше наближення PSF, далі по одному з методів робиться деконволюция, після чого деяким критерієм визначається ступінь якості, на основі неї уточнюється функція PSF і ітерація повторюється до досягнення потрібного результату.

Практика

Тепер з теорією все – перейдемо до практики, почнемо з порівняння перерахованих методів на зображенні з штучним розмиванням і шумом.

% Load image
I = im2double(imread('image_src.png'));
figure(1); imshow(I); title('Початкове зображення');

% Blur image
PSF = fspecial('disk', 15);
Blurred = imfilter(I, PSF,'circular','conv' );

% Add noise
noise_mean = 0;
noise_var = 0.00001;
Blurred = imnoise(Blurred, 'gaussian', noise_mean, noise_var);
figure(2); imshow(Blurred); title('Розмите зображення');
estimated_nsr = noise_var / var(Blurred(:));

% Restore image
figure(3), imshow(deconvwnr(Blurred, PSF, estimated_nsr)), title('Wiener');
figure(4); imshow(deconvreg(Blurred, PSF)); title('Regul');
figure(5); imshow(deconvblind(Blurred, PSF, 100)); title('Blind');
figure(6); imshow(deconvlucy(Blurred, PSF, 100)); title('Lucy');

Результати:

Фільтр Вінера

Регуляризация по Тихонову

Фільтр Люсі-Річардсона

Сліпа деконволюция

Висновок
І в кінці першої частини трохи торкнемося приклади реальних зображень. До цього всі спотворення були штучними, що звичайно добре для обкатки і вивчення, але дуже цікаво подивитися, як все це буде працювати зі справжніми фотографіями. Ось один приклад такого зображення, знятого дзеркалкою Canon 500D з ручним уведенням фокуса:

Далі запускаємо нескладний скрипт:

% Load image
I = im2double(imread('IMG_REAL.PNG'));
figure(1); imshow(I); title('Початкове зображення');

%PSF
PSF = fspecial('disk', 8);
noise_mean = 0;
noise_var = 0.0001;
estimated_nsr = noise_var / var(I(:));

I = edgetaper(I, PSF);
figure(2); imshow(deconvwnr(I, PSF, estimated_nsr)); title('Результат');

І отримуємо наступний результат:

Як видно на зображенні з’явилися нові деталі, чіткість стала набагато вищою, правда з’явилися і перешкоди у вигляді «дзвону» на контрастних кордонах.

І приклад з реальним смазом — для його здійснення фотоапарат був встановлений на штатив, виставлена відносно довга витримка і рівномірним рухом у момент спрацьовування затвора був отриманий смаз:

Скрипт приблизно той же, тільки тип PSF тепер «motion»:

% Load image
I = im2double(imread('IMG_REAL_motion_blur.PNG'));
figure(1); imshow(I); title('Початкове зображення');

%PSF
PSF = fspecial('motion', 14, 0);
noise_mean = 0;
noise_var = 0.0001;
estimated_nsr = noise_var / var(I(:));

I = edgetaper(I, PSF);
figure(2); imshow(deconvwnr(I, PSF, estimated_nsr)); title('Результат');

Результат:

Якість, знову ж, помітно покращився — стали помітні рами на вікнах, машини. Артефакти вже інші, ніж у попередньому прикладі з расфокусировкой.

На цьому цікавому і закінчимо першу частину.
У другій частині я зупинюся на проблемах обробки реальних зображень — побудови PSF і їх оцінки, розгляну більш складні і просунуті техніки деконволюції, методи усунення дефектів типу дзвону, проведу огляд та порівняння існуючого та інше.

P. S. Не так давно була опублікована стаття на хабре про Виправлення змащених фотографій у новій версії Photoshop
Для тих, хто хоче погратися з схожою технологією усунення смаза (можливо, тієї самої, що буде використовуватися в фотошопі), можна за цим посиланням завантажити демо-версію програми, подивитися приклади відновлення, а також почитати про принцип роботи.

Література
Гонсалес Р., Вудс Р. Цифрова обробка зображень
Гонсалес Р., Вудс Р., Эддинс С. Цифрова обробка зображень в середовищі MATLAB