Skip to main content

Трекинг точек. Lucas-Kanade

Введение

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

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

Описание задачи

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



Пример трекинга точки

Особые точки

Точечная особенность изображения m — это такая точка изображения, окрестность которой o(m) можно отличить от окрестности любой другой точки изображения o(n) в некоторой другой окрестности особой точки. o2(m)Слежение за точечными особенностями сцены (Point feature tracking)

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

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

Области применения

Области применения алгоритмов трекинга точек имеют разнообразные направления:

  • Слежение за объектами, получение информации об их положении в пространстве;
  • Стабилизация видеопоследовательности на каком-либо объекте;
  • Отслеживание перемещения объектов при недостатке информации, полученной из изображения (например, при перекрытиях);
  • Алгоритмы сжатия и улучшения качества видеопотока;
  • Калибровка камеры.
  • 3D реконструкция
  • Сегментирование
  • Распознавание жестов
  • Генерация панорам

Формальная постановка задачи

Пусть у нас есть последовательность изображений, обозначим её I\left(x,t\right)
Изображения представляются как некоторая дискретная двумерная функция (массив) интенсивностей пикселей в каждый момент времени t.
Рассмотрим точку u = [u_x, u_y] , принадлежащую I(x,t). Цель алгоритма — найти такую точку v = u + d на изображении I\left(v,t+1\right), что u и v «похожи». Степень похожести можно определять по разному. В общем случае вычисляется некоторый дескриптор точки, и точки сравниваются по некоторой метрике. Дескриптор и метрика выбираются в зависимости от алгоритма.

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

Существующие методы

Все работы опираются на работы Lukasa и Kanade 1981 года. Позже формулировки были изменены. Следующие алгоритмы написаны на ее основе с учетом различных афинных преобразований движения и изменения освещения.

Lucas-Kanade

В этом алгоритме движение рассчитывается самым простым образом и не учитывает возможные изменения структуры региона вокруг отслеживаемой точки.

Метод Лукаса и Канаде подразумевает, что смещение между двумя кадрами мало и одинаково в пределах соседей некоторой точки p по определению.
I_x(q_1) V_x + I_y (q_1) V_y = -I_t(q_1)
I_x(q_2) V_x + I_y (q_2) V_y = -I_t(q_2)
\vdots
I_x(q_n) V_x + I_y (q_n) V_y = -I_t(q_n)

где q_1,q_2,\dots,q_n — точки внутри окна поиска, и I_x(q_i),I_y(q_i),I_t(q_i) частные производные изображения I по x, y и времени t, вычисляемые в точке q_i в текущий момент времени.

Эти выражения можно переписать в матричной форме A v = b, где
A = \begin{bmatrix}  I_x(q_1) & I_y(q_1) \\[10pt]  I_x(q_2) & I_y(q_2) \\[10pt]  \vdots  & \vdots  \\[10pt]  I_x(q_n) & I_y(q_n)  \end{bmatrix},  \quad\quad  v =  \begin{bmatrix}  V_x\\[10pt]  V_y  \end{bmatrix},  \quad\quad  b =  \begin{bmatrix}  -I_t(q_1)\\ [10pt]  -I_t(q_2)\\ [10pt]  \vdots \\[10pt]  -I_t(q_n)  \end{bmatrix}

Эта система уравнений обычно сильно избыточна, поэтому обычно ее решают методом наименьших квадратов.
A^T A v=A^T b или 	\mathrm{v}=(A^T A)^{-1}A^T b

\begin{bmatrix}  V_x\\[10pt]  V_y  \end{bmatrix}  =  \begin{bmatrix}  \sum_i I_x(q_i)^2      & \sum_i I_x(q_i)I_y(q_i) \\[10pt]  \sum_i I_x(q_i)I_y(q_i) & \sum_i I_y(q_i)^2  \end{bmatrix}^{-1}  \begin{bmatrix}  -\sum_i I_x(q_i)I_t(q_i) \\[10pt]  -\sum_i I_y(q_i)I_t(q_i)  \end{bmatrix}

суммы берутся от i=1 до n.

Tomasi-Kanade

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

latex path not specified.:   \left\{  \begin{array}{lcl}  d_0 = 0 \\[10pt]  d_{k+1} = d_k + C^{-1}\sum_W  \left [ \left ( I \left ( x , t \right ) —  I \left ( x + d_k , t+1 \right ) \right ) \right ]  \end{array}  \right.

Shi-Tomasi-Kanade

В этом алгоритме учитываются возможные аффинные искажения. Модель смещения принимает вид  Ax + d, где A — аффинное искажение размером 2х2, а d — смещение, размером 2х1. Для вычисления смещения ищутся параметры, минимизирующие

  \epsilon = \int\limits_{x\epsilon\Omega}W\left(x\right)\left[J\left(Ax+d\right)-I\left(x\right)\right] dx,

где \Omega — окно поиска, W(x) — весовая функция в окне.

Для поиска минимума выражение дифференцируется и приравнивается к нулю. Потом производится разложение в ряд Тейлора:

  J\left(Ax+d\right) = J\left(x\right) + g^t\left(u\right).

Такое разложение дает нам систему Tz = a, где z = [ d_{xx} d_{yx} d_{xy} d_{yy} d_{x} d_{y}]^T

Тогда вектор ощибки можно записать, как
  a = \int\limits_{x\epsilon\Omega}W\left(x\right)\left[I\left(x\right)-J\left(x\right)\right]  \begin{bmatrix}  xg_x \\[10pt]  xg_y \\[10pt]  yg_x \\[10pt]  yg_y \\[10pt]  g_x \\[10pt]  g_y \\[10pt]  \end{bmatrix}  dx

Матрицу T можно разложить как   T = \int\limits_{\Omega}W\left(x\right)  \begin{bmatrix}  U & V \\[10pt]  V^T & Z  \end{bmatrix}

  U =  \begin{bmatrix}  x^2g^2_x & x^2g_xg_y & xyg^2_x & xyg_xg_y \\[10pt]   x^2g_xg_y & x^2g^2_y & xyg_xg_y &  xyg^2_y  \\[10pt]   xyg^2_x & xyg_xg_y & y^2g^2_x & y^2g_xg_y  \\[10pt]   xyg_xg_y & xyg^2_y  & y^2g_xg_y & y^2g^2_y    \\[10pt]   \end{bmatrix}

  V^T =  \begin{bmatrix}  xg^2_x & xg_xg_y & yg^2_x & yg_xg_y  \\[10pt]   xg_xg_y & xg^2_y & yg_xg_y &  yg^2_y  \end{bmatrix}

   Z =  \begin{bmatrix}  g^2_x & g_xg_y \\[10pt]   g_xg_y & g^2_y  \end{bmatrix}

Выбранный для реализации / разработанный алгоритм

Для задания был реализовал пирамидальный итеративный алгоритмы Lukas-Kanade трекера.
Мы будем искать такую точку, которая минимизирует взвешенную квадратичную модель ограничений на производные первого порядка в некотором окне поиска \Omega.

  \sum_{ x \epsilon \Omega} W^2 \left(x\right) \left [ \nabla I \left (x \right )\times d + I_t(x) \right ]^2

где W(x) — Функция весов в окне, отражающая важность близких и дальних соседей.
I_t — полуразность первого и второго кадра
\nabla I — градиент функции I, приближаемый сверткой по ядру \frac{1}{2} \begin{pmatrix}  0 & 0 & 0 \\  -1 & 0 & 1 \\0 & 0 & 0  \end{pmatrix}.
Решение имеет вид A^T W^2 A v = A^T W^2
Где для n точек из \Omega на левом кадре

A = \left[\nabla I \left ( x_1 \right ) \ldots \nabla I \left ( x_n \right ) \right ]^T

W = diag \left [ W \left( x_1 \right ) \ldots W \left (x_n \right) \right ]

 b =  -\left [ I_t \left ( x1 \right ) \ldots I_t \left ( x_n \right ) \right]^T

 v=\left[A^TW^2A\right]^{-1} A^TW^2b

Решается приближенно, когда матрица A^TW^2A невырождена, так как это матрица размером 2\times2

  A^TW^2A =  \begin{pmatrix}  \sum\limits_{ x \epsilon \Omega}W^2\left(x\right)I_x^2\left(x\right) & \sum\limits_{ x \epsilon \Omega}W^2\left(x\right)I_x\left(x\right)I_y\left(x\right)  \\   \sum\limits_{ x \epsilon \Omega}W^2\left(x\right)I_x\left(x\right)I_y\left(x\right) & \sum\limits_{ x \epsilon \Omega}W^2\left(x\right)I_x^2\left(x\right)  \end{pmatrix}.

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

Основными ограничениями алгоритма являются его предположения — постоянная яркость движущихся точек и одинаковое направление движения в пределах окна поиска. Также алгоритм плохо срабатывает на слаботекстурированных областях (когда матрица близка к вырожденной)

Эксперименты

Алгоритм был протестирован на реальных данных, а так же на стандартных последовательностях. Было произведено несколько вариантов экспериментов:

  • Построение плотного оптического потока ( для каждого пикселя );
  • Построение промежуточного кадра для видеопоследовательности;
  • Трекинг точки в GUI.

В первом случае алгоритм показал себя довольно неплохо ( Смотрите изображения справа ). Оптический поток получается рваным, потому что не накладывается никаких ограничений. В связи с этим различные точки исходных изображений могут сходится в одну на итоговом, в связи с чем картинка получается рваной. В целом такой алгоритм с настройками на низкое качество можно использовать как хорошое начальное приближение для других методов вычисления оптического потока. Правда, такие вычисления очень медленные для работы на CPU.

Разница построенного и целевого изображения
Построенное изображение

Целевой кадр
HSV карта

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

Первый кадр
Второй кадр (Оригинальный)
Второй кадр (построенный)
Третий кадр
HSV карта
Разница построенного и целевого кадра

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

Так же планирую модифицировать код для соответствия базе Middlebury и протестировать алгоритм на их базе.

Программная реализация

Реализация алгоритма была написана на C++ (Qt) с использованием OpenMP. Qt использовался для загрузки изображений. Также был реализован простой GUI для тестирования алгоритма на gif анимациях.

В целом логика алгоритма очень хорошо поддается распараллеливанию, поэтому я начал переписывать ее с использованием технологии CUDA.

Псевдокод

Вход — два изображения I,J. пусть u — точка на I. Цель — найти соответсвующую точку на J
Предобработка — построить пирамидальное представление изображений I и J

  \{I^L\}, L=0,\dots,L_m

  \{J^L\}, L=0,\dots,L_m

Инициализация предположения
  g^{L_m} = \left[ g_x^{L_m} \quad g_y^{L_m}\right]^T = \left[ 0 \quad 0\right]^T

For L=L_m down to 0 with step -1
Положение точки u на изображении I^L: u^L=[p_x\quad p_y]^T=\frac{u}{2^L}
Производная I^L по x: I_x\left(x,y\right)=\frac{I^L\left(x+1,y\right)-I^L\left(x-1,y\right)}{2}
Производная I^L по y: I_x\left(x,y\right)=\frac{I^L\left(x,y+1\right)-I^L\left(x,y+1\right)}{2}
Матрица градиентов: latex path not specified.: G = \sum\limits_{x=p_x — w_x}^{p_x + w_x}\sum\limits_{y=p_y — w_y}^{p_y + w_y}  \begin{bmatrix}  I^2_x\left(x,y\right) & I_{x}\left(x,y\right)I_{y}\left(x,y\right) \\  I_{x}\left(x,y\right)I_{y}\left(x,y\right) & I^2_y\left(x,y\right)  \end{bmatrix}
Инициализация итеративного L-K: v^0 = \left[0\quad0\right]^T
for k=1 to K with step of 1 (или до сходимости)
Разница изображений: \delta I_k \left(x,y\right)=I^L\left(x,y\right)-J^L\left(x+g^L_x+v^{k-1}_x,y+g^L_y+v^{k-1}_y\right)
Вектор несоответствия изображений: latex path not specified.:  b^k = \sum\limits_{x=p_x — w_x}^{p_x + w_x}\sum\limits_{y=p_y — w_y}^{p_y + w_y}  \begin{bmatrix}  \delta I_k \left(x,y\right)I_{x}\left(x,y\right)\\  \delta I_k \left(x,y\right)I_{y}\left(x,y\right)  \end{bmatrix}
Сдвиг: n^k = G^{-1}b_k
Предположение для следующей итерации: v^k=v^{k-1}+b^k
End of for-loop on k
Итоговый сдвиг на уровне L: d^L = v^K
Предположение для следующего уровня L-1: g^{L-1}= \left[g^{L_1}_x\quad g^{L-1}_y\right]^T = 2\left(g^L + d^L\right)
End of for-loop on L
Итоговый сдвиг: d  = g^0 + d^0

Позиция точки на J: v = u + d

Решение: соответствующая точка v на изображении J.

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

Планы на будущее

  1. Модифицировать алгоритм для использования информации о цвете;
  2. Реализовать отбор хороших точек для отслеживания;
  3. Переписать код с использованием OpenCV или других популярных фреймворков;
  4. Модифицировать код для соотвествия базе Middlebury и протестировать алгоритм на их базе.
  5. Оптимизировать код;
  6. Перенести часть или всю логику алгоритма на GPU;
  7. Реализовать работу с видеопоследовательностью (многокадровый трекинг);

Ссылки

  1. http://cgm.computergraphics.ru/content/view/54
  2. http://www.cs.cmu.edu/~efros/courses/AP06/Papers/li-siggraph-05.pdf – раздел 4.1 (Bi-directional feature tracking)
  3. http://www.robots.ox.ac.uk/~abm/docs/buchanan06interactive.pdf
  4. Реализация алгоритма
  5. простой GUI