Те, що ви хотіли знати про оптичний потік, але соромилися запитати

Оптичний потік (Optical flow) - технологія, що використовується в різних областях computer vision для визначення зрушень, сегментації, виділення об'єктів, компресії відео. Однак якщо ми захочемо його по-швидкому реалізувати в своєму проекті, прочитавши про нього на вікіпедії або де-небудь ще, то, швидше за все, дуже швидко натрапимо на те, що він працює дуже погано і збоїть при визначенні зрушень вже близько 1-2 пікселів (принаймні так було у мене). Тоді звернемося до готових реалізацій, наприклад, в OpenCV. Там він реалізований різними методами і абсолютно незрозуміло, ніж абревіатура PyrLK краще або гірше позначення Farneback або чого-небудь в цьому роді, та й доведеться забиратися зі змістом параметрів, яких в деяких реалізаціях дуже багато. Причому, що цікаво, ці алгоритми якось працюють, на відміну від того, що ми написали самі. У чому ж секрет?

Що ж таке оптичний потік

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

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

Є два варіанти розрахунку оптичного потоку: щільний (dense) і вибірковий (sparse). Sparse потік розраховує зсув окремих заданих точок (наприклад, точок, виділених деяким feature detector'ом), dense потік вважає зсув всіх точок зображення. Природно, вибірковий потік обчислюється швидше, однак для деяких алгоритмів різниця не така вже й велика, а для деяких завдань потрібно знаходження потоку в усіх точках зображення.

Для вироджених випадків можна застосовувати більш прості методи визначення зсуву. Зокрема, якщо всі точки зображення мають один і той же зсув (зображення зсунуто цілком), можна застосувати метод фазової кореляції: обчислити перетворення Фур'є для обох зображень, знайти згортку їх фаз і за нею визначити зсув (див. en.wikipedia.org/wiki/Phase_correlation). Також можна застосовувати поблочне порівняння (block matching): знаходити зсув, що мінімізує норму різниці зображень у вікні. У чистому вигляді такий алгоритм буде працювати довго і нестійко до поворотів та інших спотворень. Англійська вікіпедія зараховує перелічені алгоритми до різних варіантів обчислення оптичного потоку, але мені це здається не дуже коректним, так як ці алгоритми можуть бути застосовані і в інших цілях і не повністю вирішують дану задачу. Ми будемо називати оптичним потоком методи, засновані на локальних характеристиках зображень (те, що в англійській вікіпедії називається differential methods).

Стандартний підхід (метод Лукаса-Канаді)

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

Розглянемо математичну модель оптичного потоку, вважаючи, що у точки в результаті зміщення не змінилася інтенсивність.

Нехай - інтенсивність у певній точці (x, y) на першому зображенні (тобто у момент часу t). На другому зображенні ця точка зрушила на (dx, dy), при цьому минув час dt, тоді - це ми розклали за Тейлором функцію інтенсивності до першого члена (пізніше буде згадано, чому тільки до першого), тут - приватні похідні за координатами і часом, тобто по суті - зміна яскравості в точці (x, y) між двома кадрами.

Ми вважаємо, що у точки збереглася інтенсивність, значить

Отримуємо одне рівняння з двома невідомими (dx і dy), значить його недостатньо для вирішення, тобто тільки на цьому рівнянні далеко не поїдеш.

Найпростіше рішення проблеми - алгоритм Лукаса-Канаді. У нас же на зображенні об'єкти розміром більше 1 пікселя, значить, швидше за все, в околиці поточної точки у інших точок будуть приблизно такі ж зрушення. Тому ми візьмемо вікно навколо цієї точки і мінімізуємо (по МНК) в ньому сумарну похибку з ваговими коефіцієнтами, розподіленими по Гауссу, тобто так, щоб найбільшу вагу мали пікселі, що ближче всього знаходяться до досліджуваного. Після найпростіших перетворень, отримуємо вже систему з 2 рівнянь з 2 невідомими:

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

Чому не працює?

Якщо ми зупинимося на цьому етапі і реалізуємо цей алгоритм, то він буде успішно працювати. Але тільки якщо зсув між сусідніми зображеннями буде дуже маленький, близько 1 пікселя, і то не завжди. (Для аналізу якості генерувалися синтетичні послідовності з різним відносним зрушенням, причому це зрушення може виражатися нецілим числом пікселів, тоді результуюче зображення відповідним чином інтерполюється) Вже на зрушенні в 2 пікселі похибка буде велика, а якщо 3 і більше, то результат буде взагалі неадекватним. У чому ж справа?

Тут нам влаштувала підставу математика. Вона прищепила нам відчуття, що всі функції навколо безперервні і багато разів диференційовані. І взагалі нас в інституті привчили наближення функції в околиці точки записувати за допомогою формули Тейлора, і ми скрізь бездумно радісно користуємося цим. А тепер задумаємося, який фізичний сенс похідних в даному місці? Ми хочемо з їх допомогою визначити зміну значення функції в кінцевій околиці точки, а похідна дає уявлення про нескінченно малу околиці. Для розширення цієї околиці можна було б додати більш високий порядок похідних в розкладання Тейлора, але це призведе до нелінійностей в системі, від чого її стане істотно складніше вирішувати, а переваги будуть сумнівні, тим більше що на практиці ми маємо справу не з безперервними багаторазово диференційованими функціями, а з взагалі незрозуміло якими дискретними функціями. Тому логічніше буде шукати функцію g (x), для якої в нашому дискретному випадку якомога точніше виконується f (x) + g (x) = f (x + 1), f (x) + 2g (x) = f (x + 2), f (x) - g (x) = f (x-1), тощо. Таким чином, нам у цьому випадку потрібна не похідна, а деяка лінійна функція, яка найбільш близько лежить до точок вихідної функції. Прості математичні викладки призводять до вирішення, де. Якщо ми будували похідну по одній сусідній точці з кожного боку, то нам пощастило: у цьому випадку формула збігається з формулою наближеного обчислення похідних: g(x) = (f(x+1) – f(x-1)) / 2. Що характерно, в OpenCV при обчисленні оптичного потоку Лукаса-Канаді використовується саме така формула, до цього ми ще повернемося потім. А от якщо взяти більше точок, то формула вже стає зовсім не схожа на класичні різностні схеми для першої похідної.

Очевидно, якщо ми будуємо цю функцію, наприклад, по трьох навколишніх точках ліворуч і праворуч від вихідної, то вона ніяким чином не залежить від точок, розташованих далі, і, відповідно, при зрушенні більше трьох точок все одно у нас часто будуть виходити неадекватні результати. А ще, чим більше число точок, за якими ми будуємо цю функцію, тим більше середнє відхилення одержуваної лінії від використовуваних точок - знову ж таки через те, що у нас не лінійно міняються зображення, а чорт знає які. На практиці зрушення більше 2 пікселів вже дають неадекватно велику помилку, скільки б точок ми не взяли.

Іншим слабким місцем алгоритму є те, що ми знову ж таки маємо справу не з гладкими безперервними функціями, а з довільними, та ще й дискретними. Тому на деяких фрагментах зображення інтенсивність може «скакати» взагалі без явних закономірностей, наприклад на межах об'єктів, або через шуми. У цьому випадку ніяка функція g (x) не зможе достатньо точно описати зміни зображення в околиці точки. Щоб з цим поборотися (хоча б частково), пропонується вихідне зображення розмазати, причому корисно буде його розмазати досить сильно, тобто краще застосовувати навіть не всіма улюблений gaussian blur (усереднення з ваговими коефіцієнтами), а прямо таки box filter (рівномірне усереднення по вікну), та ще й кілька разів поспіль. Гладкість зображення для нас зараз важливіша, ніж деталізація.

Тим не менш, ці заходи так само не врятують нас від обмеження детектованого зрушення в 2-3 пікселя. І до речі, в OpenCV 1.0 була присутня така реалізація оптичного потоку, і працювала вона тільки в ідеальних умовах на дуже маленьких зрушеннях.

Що ж робити?

Отже, звичайний Лукас-Канаді добре визначає маленькі зрушення, такі, в рамках яких картинка схожа на своє лінійне наближення. Щоб з цим поборотися, скористаємося стандартним прийомом CV - multi-scaling'ом: побудуємо «піраміду» зображень різного масштабу (майже завжди береться масштабування в 2 рази по кожній осі, так простіше вважати) і пройдемо по них оптичним потоком від меншого зображення до більшого, тоді детектоване маленьке зрушення на маленькому зображенні буде відповідати великому зрушенню на великому зображенні. На найменшому зображенні ми виявляємо зсув не більше 1-2 пікселів, а переходячи від меншого масштабу до більшого, ми користуємося результатом з попереднього кроку і уточнюємо значення зсуву. Власне, в OpenCV його і реалізує функція calcOptFlowPyrLK. Використання цього пірамідального алгоритму дозволяє нам не заморочуватися обчисленням лінійної апроксимації за багатьма точками: простіше взяти більше рівнів піраміди, а на кожному рівні брати досить грубе наближення цієї функції. Тому в OpenCV і йде розрахунок всього по двох сусідніх точках. І тому стосовно цієї реалізації алгоритму наші умовиводи про перевагу апроксимуючої функції перед похідною виявилися марними: для такої кількості опорних точок похідна і є найкраща апроксимуюча функція.

А які ще бувають?

Цей алгоритм не є єдиним варіантом обчислення оптичного потоку. В OpenCV крім потоку Лукаса-Канаді є ще потік Farneback і Sim^ Flow, також часто посилаються на алгоритм Horn-Schunck.

Метод Horn-Schunck носить дещо більш глобальний характер, ніж метод Лукаса-Канаді. Він спирається на припущення про те, що на всьому зображенні оптичний потік буде досить гладким. Від того ж самого рівняння пропонується перейти до функціоналу, тобто додати вимогу на відсутність різкої зміни зрушень з ваговим коефіцієнтом порожніх. Мінімізація цього функціоналу призводить нас до системи з двох рівнянь:

У цих рівняннях лапласіан пропонують порахувати наближено: - різниця з середнім значенням. Отримуємо систему рівнянь, яку записуємо для кожного пікселя і вирішуємо загальну систему ітеративно:

У даному алгоритмі теж пропонують використовувати multi-scaling, причому рекомендують масштабувати зображення не в 2 рази, а з коефіцієнтом 0.65

Цей алгоритм був реалізований в перших версіях OpenCV, але надалі від нього відмовилися.

Farneback запропонував апроксимувати зміну інтенсивності в околиці за допомогою квадратичної форми: I = xAx + bx + c з симетричною матрицею A (по суті, розглядаючи розкладання по Тейлору до першого члена, ми брали лінійну апроксимацію I = bx + c, тобто зараз ми якраз вирішили підвищити точність наближення) Якщо зображення зрушилося в межах цієї околиці, то, підставляємо в квадратичне розкладання, розкриваємо дужки, отримуємо

.

Тепер ми можемо обчислити значення A, b, c на обох картинках, і тоді ця система стане надлишковою відносно d (особливо бентежить перше рівняння), і взагалі d можна отримати з другого рівняння: . Доводиться вдаватися до наступної апроксимації: . Позначимо ще для простоти, Тоді отримаємо просто.

Для компенсації шумів при обчисленні, знову звернемося до того припущення, що в околиці досліджуваної точки у всіх точок більш-менш однакове зрушення. Тому знову ж проінтегруємо похибку по вікну з гаусівськими ваговими коефіцієнтами w, і знайдемо вектор d, що мінімізує цю сумарну похибку. Тоді ми отримаємо оптимальне значення і відповідну мінімальну помилку. Тобто нам треба для кожної точки порахувати, усереднити по вікну, інвертувати матрицю і отримати результат. Відповідно ці твори можна порахувати для всієї картинки і використовувати заздалегідь розраховані значення для різних точок, тобто це якраз той випадок, коли має сенс рахувати dense потік.

Як завжди, у цього алгоритму є деяка кількість модифікацій і удосконалень, що в першу чергу дозволяють використовувати відому апріорну інформацію - задану початкову апроксимацію потоку - і, знову ж таки, multi-scaling.

В основі методу Sim^ Flow лежить наступна ідея: якщо ми все одно не вміємо визначати зсув більше ніж розмір вікна, за яким ми шукали похідні, то навіщо взагалі заморочуватися з обчисленням похідних? Давайте просто у вікні знайдемо найбільш схожу точку! А для вирішення неоднозначностей і для компенсації шумів врахуємо, що потік безперервний і в околиці даної точки всі точки мають майже однаковий зсув. А проблему з розміром вікна знову ж таки вирішимо за рахунок multi-scaling'a.

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

А якщо брати не інтенсивність?

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

Практика

Давайте випробуємо на практиці алгоритми, які нам пропонує OpenCV.

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

Детально синтаксис процедур описаний на сторінці з мануалом, тут я наведу витискання-переклад з моїми коментарями.

Класичний Лукас-Канаді реалізований з пірамідою в процедурі calcOpticalFlowPyrLK. Алгоритм розраховує sparse-потік, тобто для заданого набору точок на першому зображенні оцінює їх положення на другому. Вхідні параметри досить очевидні: два вхідних зображення, вхідний та вихідний набори точок, status - вихідний вектор, показує, чи знайдено успішно відповідну точку, err - вихідний вектор оцінених похибок відповідних точок, WinSize - розмір вікна, за яким відбувається гаусово усереднення, я брав 21х21 і працювало добре, maxLevel - кількість шарів у піраміді мінус один, тобто номер останнього шару, я брав 5, criteria - умова виходу з ітеративного процесу визначення зрушення (мінімізація похибки проводиться ітеративно) - цей параметр я залишав за замовчуванням, flags - додаткові прапори, наприклад можна використовувати початкове наближення потоку або вибрати метод оцінки похибки, minEigThreshold - порогове значення градієнта, нижче якого матриця вважається виродженою, я залишав за замовчуванням. Починаючи з OpenCV 2.4.1, можна використовувати заздалегідь обчислену піраміду вимасштабованих зображень.

Результат роботи - успішно і стабільно виявляються як малі, так і великі зрушення, стійкий до досить великих шумів, час роботи - близько 10 мс для 400 точок з 5-шаровою пірамідою (на core i7 950).

До речі, цей алгоритм реалізований так само на Gpu (CUDA), причому як dense, так і sparse версії.

Потік Farneback реалізується процедурою calcOpticalFlowFarneback, розраховується dense-потік, тобто зрушення кожної точки. Параметри: вхідні зображення, вихідний потік у форматі двоканальної матриці float'ів, pyr_scale визначає відношення масштабів між шарами піраміди, levels - кількість рівнів у піраміді, winsize - розмір вікна, за яким проводиться усереднення, iterations - кількість ітерацій на кожному рівні, poly_n - розмір полінома, за яким оцінюються значення A і b, poly_sigma - сигма гаусівського розмиття при згладжуванні похідних, рекомендовані значення параметрів вказані в мануалі, flags - додаткові прапорці, наприклад можна використовувати початкове наближення потоку або по-іншому усереднювати по вікну.

Цей алгоритм значно менш стабільний (за моїми спостереженнями), легше промахується на досить рівномірних картинках (мабуть, проблема у відсутності фільтрації невдалих точок), погано визначає великі зрушення. У мене відпрацьовував за 600 мс на зображенні 512х512.

Потік Sim^ Flow реалізує процедуру calcOpticalFlowSF (розраховується знову ж таки dense потік), і у неї є безліч загадкових параметрів без дефолтних значень, і взагалі на даний момент на сторінці інформація надана вельми лаконічно. Спробуємо розібратися. Перші 3 - вхідні зображення і вихідне двоканальне; layers - кількість шарів у піраміді, тобто скільки разів масштабуємо вихідне зображення; averaging_block_size - розмір вікна, в якому ми рахували функцію енергії пікселів; max_flow - максимальний зсув, який ми хочемо вміти визначати на кожному кроці, по суті він визначається розміром вікна (хоча не зовсім зрозуміло, чому він int). На цьому можна зупинитися, а можна задати ще кілька параметрів, сенс деяких з них від мене вислизає.

На сайті пропонують подивитися приклад його використання, в якому він запускається з такими параметрами: calcOpticalFlowSF(frame1, frame2, flow, 3, 2, 4, 4.1, 25.5, 18, 55.0, 25.5, 0.35, 18, 55.0, 25.5, 10);

У мене алгоритм працює значно повільніше інших, близько 9-12 секунд на картинку 512х512. Результат роботи здається більш правдоподібним, ніж Farneback, принаймні краще визначається зрушення на рівномірних картинках, помітно краще спрацьовує з великими зрушеннями.

Висновки

Якщо ви хочете використовувати десь оптичний потік, спочатку подумайте, чи потрібен він вам: часто можна обійтися більш простими методами. Братися реалізовувати потік самостійно варто тільки кілька разів подумавши: кожен алгоритм має безліч хитрощів, тонкощів і оптимізацій; що б ви не зробили, швидше за все, в OpenCV воно ж працює краще (природно, за умови, що воно там є). Тим більше що вони там щосили використовують логічні і хардварні оптимізації типу використання SSE інструкцій, багатопоточність, можливості обчислення з CUDA або OpenCL тощо. Якщо вам достатньо порахувати зсув деякого набору точок (тобто sparse потік), то можете сміливо використовувати функцію calcOpticalFlowPyrLK, воно працює добре, надійно і досить швидко. Для обчислення dense-потоку добре використовувати функцію calcOpticalFlowSF, але вона працює дуже повільно. Якщо швидкодія критична, то calcOpticalFlowFarneback, але треба ще впевнитися, що результати його роботи вас влаштують.

Література

docs.opencv.org/modules/video/doc/motion_analysis_and_object_tracking.html

Pyramidal Implementation of the Lucas Kanade Feature Tracker. Description of the algorithm — Jean-Yves Bouguet

Two-Frame Motion Estimation Based on Polynomial Expansion — Gunnar Farneback

SimpleFlow: A Non-iterative, Sublinear Optical Flow Algorithm — Michael Tao, Jiamin Bai, Pushmeet Kohli, and Sylvain Paris

Horn-Schunck Optical Flow with a Multi-Scale Strategy — Enric Meinhardt-Llopis, Javier Sanchez

en.wikipedia.org/wiki/Optical_flow