Разработка мультимедийных приложений с использованием библиотек OpenCV и IPP

Отслеживание движения и алгоритмы сопровождения ключевых точек

Разбить на страницы
Показывать лекцию целиком

Введение

Одни существа идут, другие — следуют за ними.

Презентацию к лекции Вы можете скачать здесь.

Движение – смещение одних объектов относительно других. В задачах компьютерного зрения выделяется несколько принципиально разных случаев движения [1]:

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

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

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

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

    В настоящее время выделяется несколько основных групп задач, связанных с понятием движения на наборе последовательных кадров [9]:

  • Определение движения или областей движения. Задача рассматривается в условиях неподвижной камеры и относительно постоянного фона. Наиболее простая задача, т.к. ее решение не требует анализа полученной информации о движении. Подобные задачи, как правило, возникают при разработке охранных систем слежения, в которых необходимо обнаружить несанкционированный доступ на территорию закрытого объекта, где, например, исключено движение в определенное время суток.
  • Поиск движущихся объектов . Данная задача является более сложной по сравнению с предыдущей группой задач, т.к. необходимо не просто определить область, в которой происходит движение, но и выделить движущиеся объекты. Один из возможных подходов к решению основан на определении областей движения и последующем распознавании объектов внутри полученных областей. Задача усложняется, когда требуется построить траекторию движения объекта – восстановить положения объекта на последовательном наборе изображений. В этом случае выполняется либо сопровождение областей движения и распознавание объектов в каждой области полученной последовательности, либо поиск новых объектов на сцене и их дальнейшее сопровождение. Поиск и построение траектории движущихся объектов представляет интерес для многих систем видеонаблюдения, т.к. решение данной задачи обеспечивает возможность анализа поведения объектов (покупателей в торговых центрах, автомобилей на дороге и т.п.).
  • Определение трехмерных свойств объекта из набора его последовательных двумерных изображений . Трехмерные свойства объекта, как правило, необходимы для восстановления трехмерной модели объекта, а как следствие, реконструкции трехмерной сцены. Задача актуальна для приложений компьютерной графики, связанных с созданием реалистичных сцен.
  • В настоящем разделе остановимся на некоторых наиболее известных методах решения задач первой и второй группы.

    1. Определение областей движения

    1.1. Постановка задачи выделения областей движения

    Задача выделения областей движения на видео – одна из классических задач компьютерного зрения. Входом указанной задачи является последовательность кадров $$I_1,I_2, ... ,I_N$$ некоторого видео V:

    $$I_k= \lbrace I_k(x,y),0 \leqslant x < width,\, 0 \leqslant y < height \rbrace ,k=\overline{1,N},$$

    где width – ширина кадра, height – высота кадра, а $$I_k(x,y)$$ в общем случае представляет собой вектор фиксированной размерности В случае полутонового изображения $$I_k(x,y)$$ – вектор размерности 1, который представляет интенсивность пикселя (x,y) . Если исходное изображение $$I_k$$ задано в формате RGB, то $$I_k(x,y)$$ – вектор размерности 3, каждая компонента которого определяет интенсивность пикселя по соответствующему каналу. .

    Решением настоящей задачи является совокупность областей изображения для каждого кадра видео, в которых происходит движение одного или нескольких объектов. Таким образом, в результате обработки видео необходимо сформировать набор бинарных изображений, в которых белые пиксели (интенсивность 255) соответствуют пикселям, принадлежащим движущимся объектам, а черные (интенсивность 0) – пикселям фона (4.2).

    $$M_k(x,y)=\begin{cases} 255,\text{(x,y)-пиксель объекта}\\ 0,\text{(x,y)-пиксель фона} \end{cases},k=\overline{1,N}$$

    1.2. Вычитание фона

    Наиболее простой подход к решению данной задачи состоит в том, чтобы использовать механизм вычитания фона из кадра видео (background subtraction) [6 – 9]. Процедура вычитания предполагает, что для данного видео построена модель фона (4.3):

    $$F= \lbrace F(x,y),0 \leqslant x < width,\, 0 \leqslant y < height \rbrace,$$

    также возможно существует механизм обновления модели фона с течением времени. Для одноканального (в оттенках серого) изображения, т.е. когда $$I_k(x,y),F(x,y) \in \lbrace 0,...,255 \rbrace,k=\overline{1,N},$$ процедуру вычитания можно разбить на два этапа:

  • Вычитание фонового изображения из текущего кадра видео. Данный шаг включает в себя попиксельное вычитание интенсивностей кадра видео и фонового изображения (4.4). $$D_k(x,y)=abs(I_k(x,y)-F(x,y)),k=\overline{1,N}$$
  • Отбор пикселей, принадлежащих фону и объекту, – построение бинарного изображения (маски). Считается, что пиксель принадлежит объекту и имеет белый цвет в маске, если разность интенсивности фона и текущего кадра для данного пикселя превышает некоторое пороговое значение в противном случае, принимается, что пиксель принадлежит фону (4.5). $$M_k(x,y)=\begin{cases} 255,D_k(x,y)\geqslant \tau\\ 0,D_k(x,y)< \tau \end{cases},k=\overline{1,N}$$
  • Дополнительно к указанным операциям с целью повышения качества поиска может выполняться, например, фильтрация кадров исходного потока видеоданных, либо фильтрация бинарного видео, также могут применяться морфологические операции к полученному отсечению с целью удаления шумов [8]. Если имеется цветное изображение, то его всегда можно преобразовать в оттенки серого.

    Качество определения положения движущихся областей посредством вычитания фона во многом зависит от качества построенной модели фона. Множество всех техник вычитания фона подразделяется на две группы в зависимости от механизма построения фонового изображения:

    1. Нерекурсивные. Нерекурсивные методы обновляют модель фона для текущего кадра на основании информации об интенсивностях пикселей некоторого набора предшествующих моделей фона [6] (или кадров) и текущего кадра. К наиболее распространенным нерекурсивным методам относятся следующие методы:

  • Метод вычитания текущего и предыдущего кадра . Согласно данному методу считается, что для кадра $$I_k$$ модель фона $$F_k$$ совпадает с предыдущим кадром, т.е. $$I_k(x,y)$$. Тогда на первом этапе алгоритма вычитания фона вычисляется разница пары последовательно идущих кадров (4.6). $$D_k(x,y)=abs(I_k(x,y)-F(x,y)) = \\ = abs(I_k(x,y)-I_{k-1}(x,y)),k=\overline{2,N} $$
  • Метод усреднения определенного количества предшествующих кадров . Обозначим количество кадров, по которым будет выполняться построение модели фона, как s. Тогда для кадра $$I_k$$ модель фона $$F_k$$ определяется в соответствии с формулой (4.7). $$F_k(x,y)=\frac 1 s \sum^{s-1}_{j=0} {I_{k-j}(x,y)}$$
  • Метод определения медианы фиксированного количества предшествующих кадров . Пусть s – число кадров, на основании которых будет обновляться модель фона. Тогда модель фона вычисляется согласно формуле (4.8). $$F_k(x,y)=med_{j=\overline{0,s-1}} \lbrace I_{k-j}(x,y) \rbrace$$
  • Основное преимущество нерекурсивных методов – простота реализации и скорость обновления моделей фона при переходе от кадра к кадру. При этом заметим, что качество работы методов данной группы в значительной степени зависит от скорости объектов. Медленно перемещающиеся объекты, как правило, обнаруживаются плохо. Более того, приведенные методы не дают качественный результат при изменении света на сцене, либо при наличии динамического фона (листва деревьев, струящаяся вода и т.п.). Чтобы сгладить влияние указанных эффектов, обновленная модель фона для кадра $$I_k$$ представляется выпуклой оболочкой модели фона $$F_{k-1}$$ с текущим изображением (4.9). Такая процедура называется $$\alpha$$-смешиванием.

    $$F_k(x,y)=\alpha I_k (x,y)+(1+\alpha)F_{k-1}(x,y)$$

    2. Рекурсивные. Рекурсивные методы для обновления модели фона используют информацию об интенсивностях пикселей только текущего кадра. К методам данной группы относятся гистограммный метод, метод представления модели фона смесью Гауссовых распределений (Gaussian mixture model) [2, 3, 9], метод "шифровальной" книги (codebook) [10, 11], метод извлечения визуального фона (Visual Background Extractor, ViBe) [12]. Рассмотрим некоторые из перечисленных методов.

  • Гистограммный метод. Идея гистограммного метода состоит в том, что все цветовое пространство разбивается на отдельные бины (для полутонового изображения пространство представляет отрезок изменения интенсивности, в случае цветного изображения – трехмерный куб). Для каждого изображения в последовательности выполняется построение гистограммы. Осуществляется проход по всем пикселям изображения, в зависимости от того, какая интенсивность/цвет наблюдается в пикселе, увеличивается на единицу величина соответствующего бина гистограммы. Принимается, что пиксели, составляющие некоторый бин, принадлежат фону, если величина данного бина меньше фиксированного порогового значения, в противном случае, считается, что они принадлежат объекту. Основная проблема применения гистограммного метода состоит в необходимости использования дополнительной памяти, а также в выполнении большого количества операций обращения к памяти в процессе реализации.
  • Смесь Гауссовых распределений. При построении фона с использованием данного метода считается, что для любого пикселя $$(x_0,y_0)$$ изображения $$I_k$$ известна история изменения его интенсивности/цвета на всех предшествующих кадрах $$\lbrace X_1,X_2,...,X_k \rbrace = \lbrace I_j (x_0,y_0),j=\overline{1,k} \rbrace$$ . Тогда вероятность того, что наблюдается значение $$X_k$$, может быть представлена смесью из Гауссовых распределений (4.10).
  • $$P(x_k)=\sum^s_{j=1} {\omega^k_jN(x_k\lvert \mu^k_j,\Sigma^k_j)},$$

    где $$\omega^k_j$$ – вес j-ого распределения Гаусса для кадра с номером $$k,\mu^k_j$$, – математическое ожидание, $$\Sigma^k_j)$$ – среднеквадратичное отклонение, $$N(x_k\lvert \mu^k_j,\Sigma^k_j)$$ – функция плотности нормального распределения (4.11).

    $$N(x_k\lvert \mu^k_j,\Sigma^k_j)=\frac 1 {(2\pi)^{\frac D 2}\lvert \Sigma^k_j \rvert ^{\frac 1 2}}e^{-\frac 1 2 (x_k-\mu^k_j)^T(\Sigma^k_j)^{-1}(x_k-\mu^k_j)}$$

    Предполагается, что компоненты цвета независимы и имеют одинаковое среднеквадратичное отклонение. Поэтому матрица ковариации имеет вид $$\Sigma^k_j=(\sigma^k_j)^2E$$, где E – единичная матрица.

    Указанное предположение позволяет снизить вычислительную трудоемкость метода за счет отсутствия необходимости вычислять матрицу, обратную к матрице ковариации $$\Sigma^k_j$$ в (4.11). Таким образом, задано распределение наблюдаемых значений цвета для каждого пикселя. Новое значение будет представляться одной из основных компонент построенной смеси Гауссовых распределений и использоваться для обновления параметров модели. Распределения сортируются в порядке уменьшения величины $$r^k_j=\frac {\omega^k_j} {\sigma^k_j}$$ . Такая сортировка предполагает, что пиксель фона отвечает распределению с большим весом и малой дисперсией.

    Принимается, что первые $$B^k$$ распределений, удовлетворяющих условию (4.12), соответствуют распределению цвета фоновых пикселей.

    $$B^k=argmin_b \lbrace \sum^b_{j=1} {\omega^k_j>T}\rbrace,$$

    где T – некоторое пороговое значение, параметр модели. Когда приходит очередной кадр $$I_{k+1}$$, для каждого пикселя изображения выполняется тест, который позволяет определить с использованием расстояния Махаланобиса, какому распределению соответствует полученное значение (4.13).

    $$\sqrt {(x_{k+1}-\mu^k_j)^T(\sigma^k_j)^{-1}(x_{k+1}-\mu^k_j)}< 2,5\sigma^k_j$$
  • Если нашлось соответствующее распределение Гаусса, то в зависимости от того, определяет ли оно распределение фоновых пикселей (входит в группу из $$B^k$$ распределений) или нет, текущий пиксель классифицируется как фоновый, либо как принадлежащий объекту.
  • Если не обнаружилось ни одного распределения, удовлетворяющего условию (4.13), то считается, что пиксель принадлежит объекту.
  • На основании такого правила формируется двумерная маска.

    Чтобы обработать следующий кадр, необходимо обновить параметры распределений: математическое ожидание $$\mu^k_j$$ и среднеквадратичное отклонение $$\sigma^k_j$$ . В зависимости от того, нашлось ли соответствующее распределение для цвета текущего пикселя, обновление выполняется по-разному.

  • Соответствие обнаружено. Тогда весовые коэффициенты, составляющие смесь Гауссовых распределений, которым соответствует $$X_{k+1}$$ и параметры распределений пересчитываются согласно формулам (4.14) – (4.16). $$\omega^{k+1}_j=(1-\alpha)\omega^{k}_j+\alpha,$$ $$\mu^{k+1}_j=(1-\rho)\mu^{k}_j+\rho X_{k+1},$$ $$(\sigma^{k+1}_j)^2=(1-\rho)(\sigma^{k}_j)^2+\rho(x_{k+1}-\mu^{k+1}_j)(x_{k+1}-\mu^{k+1}_j)^T,$$ где $$\alpha$$ – заданная константа, $$\rho=\alpha N(x_k\lvert \mu^k_j,\Sigma^k_j).$$. Для всех распределений, которым $$X_{k+1}$$ не соответствует, параметры не изменяются, только пересчитываются коэффициенты $$\omega^{k+1}_j$$ согласно (4.17). $$\omega^{k+1}_j=(1-\alpha)\omega^{k}_j$$
  • Соответствие не найдено. В данном случае крайнее (в смысле введенного отношения порядка) распределение Гаусса замещается распределением с новыми параметрами. Математическое ожидание выбирается равным текущему значению цвета пикселя $$\mu^{k+1}_s=X_{k+1}$$ дисперсия $$(\sigma^{k+1}_j)^2$$ максимально возможной, а вес $$\omega^{k+1}_s$$ минимально допустимым.
  • В завершении отметим, что количество распределений определяется сложностью фона и имеющимися вычислительными мощностями (в [2] предлагается использовать значение в пределах от 3 до 5). Начальная инициализация параметров распределений может выполняться с использованием метода k-средних [2], либо EM-алгоритма (Expectation Maximization) [4, 32].

    На рис.4.1 показан пример применения метода вычитания фона с представлением модели фона смесью Гауссовых распределений.

    (рис 4.1) Пример применения метода вычитания фона с использованием метода представления фона смесью Гауссовых распределений

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

    Метод извлечения визуального фона (Visual Background Extractor , ViBe) [12]. В соответствии с данным методом модель фона на кадре с номером $$ k$$ представляется набором множеств $$M^k(p)=\lbrace v_1,v_2,...,v_N\rbrace$$ для всех пикселей $$p=(x,y)$$ где $$v_i$$ – интенсивность/цвет пикселя (в общем случае вектор).

    Для классификации пикселя в цветовом пространстве строится сфера $$S_R(v(x))$$ радиуса R и определяется количество векторов множества $$M(p)$$ которые попадают вовнутрь построенной сферы, $$K=\lvert S_R(v(p)) \cap M(p) \rvert$$. Если $$K>T_{min}$$ , где $$T_{min}$$ – фиксированное пороговое значение, то принимается, что пиксель принадлежит фону, в противном случае, объекту.

    На начальном этапе необходимо выполнить инициализацию множеств $$M^0(p)$$ для всех пикселей следующим образом:

    $$M^0(p)=\lbrace v^0(y),y \in N_G(p) \rbrace,$$

    где $$N_G(p)$$ – окрестность пикселя размера 3x3 (9 клеток, включая текущий пиксель), y выбирается $$N$$ раз случайным образом. Обновление модели фона для кадра $$I_k$$ выполняется в два шага:

  • Если p проклассифицирован как пиксель фона, то из множества $$M^k(p)$$ случайно выбирается компонента, которая заменяется значением $$v(p)$$.
  • Случайным образом выбирается один соседний пиксель из окрестности $$N_G(p)$$ , для которого выполняется предыдущий шаг.
  • Множество методов построения моделей фона не ограничивается набором методов, представленных в этом разделе. В данном направлении ведутся активные исследования до настоящего момента. Поэтому здесь представлены лишь принципиально разные подходы к решению задачи построения фоновых моделей.

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

    Еще один простейший способ оценить движение на нескольких изображениях – перебрать все возможные варианты смещений изображений или их отдельных фрагментов (translational alignment, [5, гл. 8.1]). Для этого первоначально необходимо выбрать метрику для оценки степени сходства фрагментов. Как следствие, исходная задача определения движения может быть сведена к минимизации функции ошибки по всем возможным направлениям смещения (рис.4.2).

    (рис 4.2) Пример перебора вариантов по четырем принципиальным направлениям

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

  • 1. Квадратичная функция ошибки (4.19). $$E(u,v)=\sum_i {(I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i))^2},$$ где $$(u,v)$$ – вектор смещения изображения или фрагмента.
  • 2. Грубые оценочные функции ошибки. Общий вид функций данной группы определяется формулой (4.20). При этом $$\rho (e)$$ – функция, которая имеет меньшую степень роста по сравнению квадратичной функцией. $$E(u,v)=\sum_i {\rho (I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i))},$$ Примером грубой функции ошибки может служить сумма абсолютных разностей (4.21). $$E(u,v)=\sum_i {\lvvert I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i) \rvert},$$
  • 3. Взвешенная квадратичная функция ошибки (4.22). Весовые коэффициенты – дискретные функции, которые принимают нулевые значения за пределами изображения. $$E(u,v)=\sum_i {\omega_{k-1} (x_i,y_i)\omega_k(x_i+u,y_i+v) \cdot (I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i))^2},$$
  • 4. Кросс-корреляция пары выровненных изображений (4.23). Заметим, что при выборе кросс-корреляции в качестве функции ошибки выполняется не минимизация, а максимизация по всем возможным направлениям смещения. $$E(u,v)=\sum_i {I_{k-1} (x_i,y_i) \cdot I_k(x_i+u,y_i+v)$$
  • В общем случае если задан набор цветных изображений, то всегда можно выполнить их конвертирование в полутоновые, либо ввести дополнительную сумму по числу каналов, и использовать одну из предложенных функций ошибки.

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

    Еще один способ ускорения вычисления при реализации переборной схемы – использование подхода, основанного на вычислении быстрого преобразования Фурье (БПФ) [30]. Метод опирается на тот факт, что модули образов Фурье исходного и смещенного сигналов связаны уравнением (4.24).

    $$F(I_k(x_j+u,y_j+v))=F(I_k(x_i,y_i))e^{-i(\begin{bmatrix} u \\ v \end{bmatrix},\overline{\omega})}=I_k(\overline{\omega})e^{-i(\begin{bmatrix} u \\ v \end{bmatrix},\overline{\omega})},$$

    где $$\overline{\omega}=\begin{bmatrix} \omega_x \\ \omega_y \end{bmatrix}$$ – вектор частот двумерного БПФ, $$I_k(\overline{\omega})=F(I_k(x_i,y_i))$$ – образ Фурье сигнала $$I_k(x_i,y_i)$$. Допустим, что в качестве функции ошибки выбрана кросс-корреляция, вычислим БПФ от этой функции (4.25). В результате получаем преобразование Фурье от свертки функций. Исходя из свойств преобразования, образ Фурье от полученной свертки функций – произведение образа $$I_{k-1}(\overline{\omega})$$ первой функции на комплексно- сопряженный образ $$I^*_k(\overline{\omega})$$ второй функции.

    $$F(E(u,v))=F \Biggl( \sum_i {I_{k-1} (x_i,y_i) \cdot I_k(x_i+u,y_i+v) \Biggr) = F(I_{k-1}(u,v)+I_k(u,v))=I_{k-1}(\overline{\omega})I^*_k(\overline{\omega})$$

    Отсюда эффективность оценивания кросс-корреляции с использованием преобразования Фурье очевидна. Достаточно вычислить образы Фурье для исходной пары изображений с помощью БПФ (сложность составляет $$O(NM\,log\, NM))$$ для изображения размера $$N \times M$$ , перемножить полученные образы (за время, не превышающее $$O(N^2M^2))$$) , затем вычислить обратное преобразование Фурье. В результате получится матрица значений кросс-корреляционной функции, из которой необходимо выбрать максимальное значение. Приведенный пример не является единственным примером ускорения вычислений функции ошибки с помощью БПФ [5].

    Для некоторых практических приложений, в частности, для стабилизации видео необходима субпиксельная точность определения смещения – смещение в дробных координатах. Один из наиболее известных подходов состоит в том, чтобы использовать процедуру уточнения направления смещения. Впервые была введена при разработке алгоритма сопровождения Лукаса-Канаде, о котором речь пойдет в одном из следующих разделов. Процедура основана на применении метода градиентного спуска для минимизации квадратичной функции ошибки в смещенной точке – функции энергии (4.26), которая приближается разложением в ряд Тейлора.

    $$E(u+\Delta u,v+\Delta v)=\sum_i {(I_k(x_i+u+\Delta u,y_i+v+\Delta v)-I_{k-1} (x_i,y_i))^2} \approx$$ $$\approx \sum_i {\Biggl( I_k(x_i+u,y_i+v)+\nabla I_k(x_i+u,y_i+v)}\begin{bmatrix} \Delta u \\ \Delta v \end{bmatrix}-I_{k-1} (x_i,y_i) \Biggr) ^2=$$ $$= \sum_i {(\nabla I_k(x_i+u,y_i+v)\begin{bmatrix} \Delta u \\ \Delta v \end{bmatrix}+e_i)^2}=$$ $$\sum_i {(J_k(x_i+u,y_i+v)\begin{bmatrix} \Delta u \\ \Delta v \end{bmatrix}+e_i)^2} \rightarrow min_{(u,v)},$$ $$e_i=I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i),$$ $$J_k(x_i+u,y_i+v)=\nabla I_k(x_i+u,y_i+v)=[\frac {\deltaI_k} {\delta x},\frac {\deltaI_k} {\delta y}]\rvert_{(x_i+u,y_i+v)}$$

    Полученная задача минимизации (4.27) сводится к решению системы линейных алгебраических уравнений вида (4.28). Доказательство данного факта можно найти в [31].

    $$A\begin{bmatrix} \Delta u \\ \Delta v \end{bmatrix}=b$$

    где $$A=\sum_i {J^T_k(x_i+u,y_i+v)J_k(x_i+u,y_i+v)},b=-\sum_i {e_iJ^T_k(x_i+u,y_i+v)}.$$ .

    Эффективность решения задачи достигается за счет того, что для текущего изображения градиенты $$J_k(x_i+u,y_i+v)$$ в смещенной точке можно приближенно заменить градиентами предыдущего изображения $$J_{k-1}(x_i,y_i)$$ в исходной точке, т.е. $$J_k(x_i+u,y_i+v)\approx J_{k-1}(x_i,y_i)$$.

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

    1.4. Параметрические модели движения

    Многие прикладные задачи, такие, как склеивание изображений в панораму, стабилизация видео, требуют построения более сложных моделей движения, т.к. движение не описывается только линейным сдвигом. Поэтому рассматриваются пространственные поля смещений, и выполняется построение параметрических моделей движения (parametric motion, [5, гл. 8.2]). В параметрических моделях вместо постоянного вектора смещения$$[u,v]^T$$ рассматривается параметризованное поле смещений $$x'(x,y;p)$$, где p-вектор параметров. Приведем некоторые примеры параметризованных моделей:

  • Линейный сдвиг (4.29). $$x'(x,y;p)=\begin{bmatrix} E p \\ {0^T} 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix},$$ где $$p=\begin{bmatrix} {p_1} \\ {p_2} \end{bmatrix}$$ – вектор сдвига.
  • Поворот со сдвигом (4.30). $$x'(x,y;p)=\begin{bmatrix} {cos\,\phi} {-sin\,\phi} {p_1} \\ {sin\,\phi} {cos\,\phi} {p_2} \\ 0 0 1 \end{bmatrix}\begin{bmatrix} x \\ y \\ 1 \end{bmatrix},$$ где $$R=\begin{bmatrix} {cos\,\phi} {-sin\,\phi} \\ {sin\,\phi} {cos\,\phi} \end{bmatrix}$$ – ортонормированная матрица поворота, т.е. $$RR^T=1$$ и $$det\,R=1$$ , $$p=(p_1,p_2,\phi)$$ – набор параметров модели движения.
  • Преобразование подобия (4.31). $$x'(x,y;p)=\begin{bmatrix} a -b {p_1} \\ b a {p_2} \\ 0 0 1 \end{bmatrix}\begin{bmatrix} x \\ y \\ 1 \end{bmatrix},$$ где $$sR=\begin{bmatrix} a -b \\ b a \end{bmatrix}$$ – матрица поворота, умноженная на коэффициент масштабирования $$s,\,p=(p_1,p_2,a,b)$$ – набор параметров модели движения.
  • Аффинное преобразование (4.32). $$x'(x,y;p)=\begin{bmatrix} a_{00} a_{01} a_{02} \\ a_{10} a_{11} a_{12} \end{bmatrix}\begin{bmatrix} x \\ y \\ 1 \end{bmatrix}$$
  • Перспективная проекция или гомография (4.33). $$x'(x,y;p)=\begin{bmatrix} h_{11} h_{12} h_{13} \\ h_{21} h_{22} h_{23} \\ h_{31} h_{32} h_{33} \end{bmatrix}\begin{bmatrix} x \\ y \\ 1 \end{bmatrix}=H\overline{x}$$
  • Аналогично процедуре уточнения направления, описанной в разделе 1.3, можно выполнить определение направления движения для любой параметрической модели (4.35), расписав квадратичную функцию ошибки (4.34) с помощью разложения в ряд Тейлора по вектору параметров и выполнив ее минимизацию.

    $$E(p+\Delta p)=$$ $$\sum_i {(I_k(x'(x_i,y_i;p+\Delta p))-I_{k-1}(x_i,y_i))^2}\Approx$$ $$\sum_i {(I_k(x'_i) +I_k(x'_i)\Delta p-I_{k-1}(x_i,y_i))^2}=$$ $$\sum_i {(J_k(x'_i)\Delta p+e_i)^2} \rightarrow min_{\Delta p},$$ $$e_i=I_k(x'_i)-I_{k-1}(x_i,y_i),$$ $$J_k(x'_i)=\frac {\delta I_k} {\delta p}=\nabla I_k(x'_i)\frac {\delta x'_i} {\delta p}\rvert_{(x_i,y_i)}.$$

    1.5. Многоуровневое движение

    Во многих случаях визуальное движение вызвано смещением небольшого количества объектов, находящихся на разной глубине изображения [5]. Поэтому движение пикселей можно описать более эффективно, если сгруппировать их в слои [14], и отслеживать многоуровневое движение (layered motion, [5, гл. 8.5]) построенных слоев, например, с помощью параметрических моделей [15].

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

    1. Как представить слой? Слой определяется набором из трех карт (матриц):

  • матрица интенсивности слоя (в компьютерной графике текстурная карта);
  • $$\alpha$$карта, которая определяет прозрачность слоя в каждом точке изображения;
  • карта скоростей описывает изменение положения точек с течением времени.
  • Ниже (рис.4.3) показан пример представления кадров в виде набора слоев. Изображения в строке (а) соответствуют фоновому слою, в строке (b) – слою, содержащему объект (рука), в (с) – последовательность изображений, восстановленная на основании выделенных слоев.

    (рис 4.3) Пример представления изображения в виде набора слоев [15]

    Предположим, что $$E_0(x,y)$$ – матрица интенсивности фонового слоя, $$E_1(x,y)$$– слоя, содержащего объект, $$\alpha_1(x,y)$$ – $$\alpha$$-карта слоя $$E_1(x,y)$$. Поскольку слои перекрываются, то для представленного примера восстановить можно исходное изображение согласно (4.36). $$I_1(x,y)=E_0(x,y)(1-\alpha_1(x,y))+E_1(x,y)\alpha_1(x,y)$$

    Отметим, что данную процедуру можно распространить на случай большего количества слоев (4.37). $$I_k(x,y)=I_{k-1}(x,y)(1-\alpha_k(x,y))+E_k(x,y)\alpha_k(x,y)$$

    2. Как разбить множество пикселей на слои и определить движение каждого слоя? Процедура разбиения на слои состоит из нескольких шагов:

  • Оценивание движения с использованием оптического потока и аффинной модели для набора неперекрывающихся блоков (рис.4.4, а), полученных в результате разбиения исходного изображения. Данная операция выполняется независимо для каждого кадра из подмножества изображений.
  • Кластеризация полученных оценок с помощью метода k-средних (рис.4.4, b). Кластеризация также осуществляется независимо для каждого изображения. Каждый кластер определяет сегмент, отвечающий некоторому слою изображения.
  • Применение медианной фильтрации для получения смешанных слоев, устойчивых к незначительным изменениям интенсивности, а также для выявления перекрытий между слоями (рис.4.4, c).
  • (рис 4.4) Пример выделения слоев [14]

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

    1.6. Вычисление оптического потока

    Другой распространенный подход к решению задачи детектирования областей движения – вычисление оптического потока (optical flow) [5, 9, 13]. Оптический поток позволяет определить смещение каждой точки кадра, т.е. построить поле скоростей.

    Применение данного подхода требует выполнения двух основных условий:

  • яркость каждой точки объекта не изменяется с течением времени;
  • ближайшие точки, принадлежащие одному объекту, в плоскости изображения двигаются с похожей скоростью.
  • Рассмотрим схему работы метода вычисления оптического потока. Допустим, что имеется непрерывно изменяющееся во времени изображение. Обозначим через $$I_k(x,y)$$ яркость пикселя с координатами $$(x,y)$$ в момент времени k. Считая динамическое изображение функцией положения и времени, яркость пикселя в следующий момент времени можно выразить с помощью разложения в ряд Тейлора:

    $$I_{k+\Delta x}(x+\Delta x,y+\Delta y)=I_k(x,y)+\frac {\delta I_k} {\delta x}\Delta x+\frac {\delta I_k} {\delta y}\Delta y+\frac {\delta I_k} {\delta k}\Delta k+O(\delta^2)$$

    где $$\frac {\delta I_k} {\delta x}\Delta x,\frac {\delta I_k} {\delta y}\Delta y,\frac {\delta I_k} {\delta k}\Delta k$$ – частные производные Под частными производными в силу дискретности представления изображения понимается центральный разностный оператор. функции $$I_k$$. Предполагается, что за небольшой промежуток времени пиксель смещается незначительно, т.е. $$(\Delta x, \Delta y) \rightarrow 0$$ при $$\Delta k \rightarrow 0$$, значит можно считать, что яркость почти не изменяется, поэтому остаточным членом ряда можно пренебречь:

    $$I_{k+\Delta x}(x+\Delta x,y+\Delta y)=I_k(x,y)$$ $$-\frac {\delta I_k} {\delta k}= \frac {\delta I_k} {\delta x}\frac {\delta x} {\delta k}+\frac {\delta I_k} {\delta y}\frac {\delta y} {\delta k}$$

    Уравнение (4.39) называется уравнением оптического потока . Дальнейшая цель – построить вектор скорости пикселя $$c=(\frac {\delta x} {\delta k},\frac {\delta y} {\delta k})=(u,v)$$ . Для полноты постановки задачи вводится условие гладкости изменения скорости. В результате исходная задача сводится к задаче минимизации квадратичной ошибки $$(\frac {\delta I_k} {\delta x}u+\frac {\delta I_k} {\delta y}v+\frac {\delta I_k} {\delta k}k)^2$$ при наличии ограничений в виде равенств $$u^2_x+u^2_y=0$$ и $$v^2_x+v^2_y=0$$ (условия гладкости). Получаем задачу математического программирования. Процедура минимизации применяется к каждому пикселю текущего изображения, в результате чего обеспечивается построение поля векторов смещения всех пикселей.

    2. Методы сопровождения объектов

    2.1. Постановка задачи сопровождения объектов

    Сопровождение (трекинг) движущихся объектов – это один из составляющих компонентов многих систем реального времени таких, как системы слежения, анализа видео и других. Входными данными любого алгоритма сопровождения является последовательность изображений (кадров видео) $$I_1,I_2,...,I_N$$(4.1) с нарастающим объемом информации, которую необходимо обрабатывать и анализировать.

    Задача сопровождения состоит в том, чтобы построить траектории движения целевых объектов на входной последовательности кадров. Допустим, что положение объекта на изображении с номером обозначается $$P_k$$ . Тогда траекторией движения объекта называется последовательность его положений $$P_s,P_{s+1},...,P_{s+l-1}$$ , где s – номер первого кадра, на котором был обнаружен объект, l – количество кадров последовательности, где наблюдается объект. Заметим, что в зависимости от метода сопровождения положение объекта может определяться по- разному (координаты и размер сторон окаймляющего прямоугольника, координаты центра масс контура и т.п.).

    2.2. Классификация методов сопровождения

    Существует несколько категорий методов сопровождения объектов [16]:

  • Методы сопровождения особых точек (point tracking). В таких методах принимается, что положение объекта определяется расположением набора характерных точек. Один и тот же объект на последовательных кадрах представляется наборами соответствующих пар точек. Данная группа методов разделяется на две подгруппы:

    Детерминистские методы [17] используют качественные эвристики движения (небольшое изменение скорости, неизменность расстояния в трехмерном пространстве между парой точек, принадлежащих объекту), по существу задача сводится к минимизации функции соответствия наборов точек. Методы, основанные на вычислении плотного и разреженного оптического потока, а также методы сопоставления (matching) дескрипторов ключевых точек являются типичными представителями детерминистских методов.

    Вероятностные методы используют подход, основанный на понятии пространства состояний. Считается, что движущийся объект имеет определенное внутреннее состояние, которое измеряется на каждом кадре. Чтобы оценить следующее состояние объекта, требуется максимально обобщить полученные измерения, т.е. определить новое состояние при условии, что получен набор измерений для состояний на предыдущих кадрах. Типичным примерами таких методов являются методы на базе фильтра Кальмана [6, 18 – 20] и фильтра частиц (particle filter) [21, 22, 33].

  • Методы сопровождения компонент (kernel tracking). Под компонентой понимается форма объекта. В простейшем случае компонента может быть представлена шаблоном прямоугольной или овальной формы, в более сложных – трехмерной моделью объекта, спроецированной на плоскость изображения. Как правило, методы данной группы применяются, если движение определяется обычным смещением, поворотом или аффинным преобразованием. Трекинг компонент – итеративная процедура локализации, основанная на максимизации некоторого критерия подобия. На практике реализуется с использованием сдвига среднего (mean shift) [8, 23] и его непрерывной модификации (Continuous Adaptive Mean Shift, CAM Shift) [8, 24].
  • Методы сопровождения силуэта (silhouette tracking). Силуэт может быть задан контуром, либо набором связанных простых геометрических примитивов. Задача трекеров силуэта состоит в том, чтобы на каждом кадре определить область, в которой находится объект, с использованием модели его силуэта, построенной на основании предшествующих кадров. Можно выделить две группы методов:
  • Методы сопоставления и сопровождения фрагментов изображения, содержащих объект. При сопоставлении естественным образом должна быть введена мера сходства пары областей, ограниченных контуром. На практике используется расстояние Хаусдорфа, норма $$L_2$$ , также значение оценочной функции может вычисляться с использованием нескольких мер, включая, например, кросс-корреляцию, расстояние Бхачатария. Определение положения фрагмента на следующем кадре последовательности выполняется посредством вычисления оптического потока (раздел 4.6) для внутренних точек области.
  • Методы сопровождения контура. Позволяют прогнозировать положение контура на следующем кадре. Первый подход состоит в использовании моделей пространства состояний (по типу фильтра Кальмана), второй – в минимизации функции энергии контура с использованием прямых техник таких, как градиентный спуск (раздел 4.3).
  • 2.3. Методы сопоставления ключевых точек

    В настоящее время детерминистские методы сопровождения особых точек чаще остальных используются на практике. Необходимым условием использования этих методов является наличие выделенных точек, которые наилучшим образом представляют объект. Выделение особых точек выполняется с использованием специальных детекторов и дескрипторов Обзор существующих детекторов и дескрипторов приведен в предыдущей лекции. . Если детектор находит положение особой точки на изображении, то дескриптор дополнительно строит вектор признаков, характерных для полученной точки.

    Схема сопровождения особых точек может быть представлена в виде последовательности действий [34]:

  • Поиск особых точек на предыдущем (рис.4.5, слева) и текущем (рис.4.5, справа) кадрах посредством выбранного детектора (LoG, DoG, SURF, SIFT, ORB и др.).
  • Вычисление дескрипторов для полученного набора точек (SURF, SIFT, GLOH, DAISY, BRIEF и т.д.) – n-мерных векторов-описателей.
  • Сопоставление (matching) дескрипторов, полученных на текущем и предыдущем кадре. По существу необходимо построить преобразование точек n-мерного пространства дескрипторов. Прямой перебор дескрипторов предыдущего изображения и поиск наиболее близких дескрипторов текущего не является эффективным из- за большого количества ложных соответствий (рис.4.6). Близость дескрипторов определяется посредством вычисления некоторой метрики (норма $$L_2$$ , кросс-корреляция и т.п.).
  • (рис 4.5) Пример исходного изображения объекта (слева) и объекта внутри некоторой сцены (справа) (рис 4.6) Набор соответствий до применения алгоритма RANSAC

    После сопоставления для отсечения ложных срабатываний (рис.4.7) применяется алгоритм RANSAC (RANdom SAmple Consensus, [6, 7]).

    (рис 4.7) Набор соответствий после применения алгоритма RANSAC

    RANSAC – это общий метод, который используется для оценки параметров модели на основании случайных выборок. При сопоставлении модель представляет собой матрицу преобразования (гомография). На входе алгоритма имеется два множества дескрипторов, полученных на предыдущем и текущем изображении. Схема работы RANSAC состоит из многократного повторения трех этапов:

  • Выбор точек и построение параметров модели. Из входных множеств дескрипторов выбираются случайным образом без повторений наборы фиксированного размера. На основании полученных наборов строится матрица преобразования.
  • Проверка построенной модели. Для каждого дескриптора предыдущего кадра находится проекция на текущем кадре и выполняется поиск наиболее близкого дескриптора из множества дескрипторов текущего кадра. Дескриптор помечается как выброс, если расстояние между проекцией и соответствующим дескриптором текущего изображения больше некоторого порога.
  • Замещение модели. После проверки всех точек проверяется, является ли построенная модель лучшей среди набора предшествующих моделей.
  • В результате применения RANSAC строится наилучшая матрица гомографии. Вычислив перспективную проекцию набора дескрипторов предыдущего кадра, достаточно выполнить проход по всем соответствиям, полученным в процессе перебора, и проверить, является ли соответствующий дескриптор текущего кадра достаточно близким к проекции дескриптора предыдущего кадра. Если не является, то пара отбрасывается.

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

    2.4. Методы сопровождения объектов, основанные на вычислении оптического потока

    С каждым пикселем изображения можно связать некоторый вектор скорости, который определяет, какое расстояние "прошел" пиксель в течение временного промежутка между предыдущим и текущим кадрами. Такая конструкция, построенная для каждого пикселя изображения, представляет собой плотный оптический поток. Одним из наиболее известных алгоритмов сопровождения, основанных на вычислении плотного оптического потока, является метод Хорна (Horn-Schunk method) [8, 25]. По существу данный метод предполагает решение системы дифференциальных уравнений, полученных в задаче математического программирования (раздел 1.6), с помощью итерационных методов для построения компонент вектора скорости в каждой точке изображения.

    Другой класс алгоритмов, использующих плотный оптический поток, – алгоритмы сопоставления блоков (block matching algorithms, [8]). Идея состоит в том, что предыдущее и текущее изображения разбиваются на блоки, как правило, квадратные и перекрывающиеся, а затем определяется движение этих блоков посредством сопоставления. Поскольку алгоритмы сопоставления работают с блоками, то изображение поля скоростей обычно имеет меньшее разрешение по сравнению с исходным изображением.

    Алгоритм Лукаса-Канаде (Lucas-Kanade) основан на вычислении разреженного оптического потока, т.е. на построении векторного поля скоростей для выделенного набора точек. Задача слежения без учета аффинных искажений с помощью данного алгоритма сводится к поиску оптического потока в особых точках [26]. Впоследствии появились модификации данного алгоритма Томаши-Канаде (Tomasi-Kanade) и Ши- Томаши-Канаде (Shi-Tomasi-Kanade). Трекер Ши-Томаши-Канаде впервые учитывает аффинные искажения окрестных точек. Также существует модификация метода для случая переменного освещения [27].

    В литературе можно обнаружить различные модификации перечисленных методов. В [28] предлагается сравнение производительности различных алгоритмов сопровождения с использованием оптического потока.

    2.5. Применение метода "сдвиг среднего" к решению задачи сопровождения

    Несколько в стороне стоят методы сопровождения, основанные на использовании техники "сдвиг среднего" (mean shift) [8, 23, 24]. Идея методов состоит в том, что для каждой особой точки (в общем случае, для каждого объекта) выбирается окно поиска, вычисляется центр масс распределения интенсивностей (гистограммы). Соответственно центр окна смещается в центр масс, который представляет собой положение точки на текущем кадре. Определение положения точки на последующих кадрах сводится к применению очередного шага метода "сдвига среднего". Метод останавливается, когда центр масс перестает смещаться.

    2.6. Применение фильтра Кальмана для оценки положения объекта

    Задачу сопровождения можно рассматривать как хорошо изученную проблему теории управления, которая состоит в том, чтобы оценить состояние системы на основании последовательности зашумленных измерений [9]. Формально имеется модель объекта, которая наблюдается в зашумленном пространстве сцены. Обозначим через x – модель объекта, z – наблюдение объекта. В общем случае x,z – вектора признаков, которые вполне могут иметь разную размерность. В процессе сопровождения объекта модели и наблюдения могут использоваться двумя способами:

  • Последовательность наблюдений $$z_1,z_2,...$$ применяется для уточнения базовой модели x . Заметим, что модель с течением времени также может изменяться, поэтому каждое новое наблюдение $$z_k$$ может давать новую оценку модели $$x_k$$.
  • Построенная оценка модели $$x_k$$ используется для предсказания модели $$x_{k+1}$$ и наблюдения $$z_{k+1}$$.
  • В результате получаем механизм обратной связи: наблюдаем $$z_k$$, оцениваем $$x_k$$, предсказываем $$x_{k+1}$$ и $$z_{k+1}$$ , обновляем $$x_{k+1}$$ на основании $$z_{k+1}$$. Такой механизм лежит в основе методов сопровождения, основанных на фильтрах Кальмана и фильтрах частиц. В данном разделе остановимся на применении фильтра Кальмана.

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

    $$x_{k+1}=F_kx_k+w_k,$$ $$z_k=H_kx_k+v_k,$$

    где $$F_k$$– матрица преобразования состояния системы; $$w_k$$ – "белый" шум с нормальным распределением $$N(0,Q_k)$$ с математическим ожиданием, равным 0, и матрицей ковариации $$Q_k;H_k$$– матрица связи модели и наблюдения; $$v_k$$ – "белый" шум с нормальным распределением $$N(0,R_k)$$. Заметим, что по определению матрицы ковариации $$Q_k=M(w_kw^T_k),$$, $$(Q_k)_{ij}=M(w^i_kw^j_k),$$ . Если X,Y – две случайные величины, определенные на одном и том же вероятностном пространстве, тогда ковариация численно равна математическому ожиданию от произведения случайных величин X-MX и Y-MY , где MX,MY – математическое ожидание случайных величин X,Y , т.о. cov(x,Y)=M((x-MX)(Y-MY)).

    Приведем пример модели для случая равномерного движения [8]. Состояние системы определяется положением и скоростью точки, тогда состояние может быть представлено вектором вида (4.42).

    $$x_k=\begin{bmatrix} x \\ y \\ v_x \\ v_y \end{bmatrix}_k$$

    Матрица преобразования состояния $$F_k=F$$ и, исходя из физических соображений, может быть записана согласно (4.43).

    $$F=\begin{bmatrix} 1 0 dt 0\\ 0 1 0 dt\\ 0 0 1 0 \\ 0 0 0 1 \end{bmatrix}$$

    Поскольку при сопровождении камера регистрирует только координаты положения объекта, то $$z_k=\begin{bmatrix} z_x \\ z_y \end{bmatrix}_k$$ Поэтому матрицу связи наблюдения и предсказания $$H_k=H$$ можно записать в соответствии с (4.44).

    $$H=\begin{bmatrix} 1 0 \\ 0 1 \\ 0 0 \\ 0 0 \end{bmatrix}$$

    Т.к. в реальных условиях нельзя наблюдать равномерное движение, то вводится Гауссов шум с матрицей ковариации $$Q_k$$ . Матрица ковариации $$R_k$$ строится, исходя из того, насколько точно выполняются измерения положения объекта.

    Теперь необходимо получить обобщенные уравнения для обновления состояния и модели. Идея состоит в том, что сначала строится априорная оценка $$\tilde{x}^-_k$$ состояния согласно (4.40):

    $$\tilde{x}^-_k=F_{k-1}x_{k-1}+w_{k-1}.$$

    Затем с использованием построенной оценки вычисляется наблюдение

    $$z_k=H_k\tilde{x}^-_k+v_k.$$

    Также можно определить апостериорную оценку $$\tilde{x}^+$$ состояния, которая строится после наблюдения. Введем ошибки $$e^-_k$$(4.47) и $$e^+_k$$(4.48), связанные с каждой оценкой состояния.

    $$e^-_k=x_k-\tilde{x}^-_k$$ $$e^+_k=x_k-\tilde{x}^+_k$$

    Фильтр Кальмана оперирует разностью $$z_k-H_k\tilde{x}^-_k$$ , в которую вносит вклад ошибка $$e^-_k$$ и случайный шум $$v_k$$ . В идеальном случае шум отсутствует и оценка состояния идеальна, поэтому указанная разность обращается в ноль. Задача состоит в том, чтобы построить матричный коэффициент Кальмана $$K_k$$ с целью обновления апостериорной оценки согласно (4.49). Отсюда, если $$K_k$$ известен, то известен закон получения обновленной модели $$x_k$$ .

    $$\tilde{x}^+_k=\tilde{x}^-_k+K_k(z_k-H_k\tilde{x}^-_k)$$

    Уравнения (4.47) – (4.49) позволяют получить выражение связи априорной и апостериорной ошибки (4.50).

    $$e^+_k=x_k-\tilde{x}^+_k=x_k-((E-K_kH_k)\tilde{x}^-_k-K_kz_k)=x_k-(E-K_kH_k)\tilde{x}^-_k+K_k(H_kx_k+v_k)=$$ $$=(E-K_kH_k)e^-_k+K_kv_k$$

    Из определения матрицы ковариации пары случайных величин следуют выражения (4.51).

    $$P^-_k=M(e^-_ke^{-T}_k),\,\,\,\,P^+_k=M(e^+_ke^{+^T}_k),\,\,\,\,R_k=M(v_k,v^T_k)$$

    Вследствие независимости ошибок $$e^-_k$$ и $$e^+_k$$ получаем нулевую ковариацию для величин $$e^-_k$$ и $$v_k$$ (4.52).

    $$M(e^-_kv^T_k)=M(v_ke^{-T}_k)=0$$

    Проделав несложные преобразования (4.53) (умножение (4.50) на $$e^{+^T}_k$$ и определение математического ожидания от полученного выражения), можно перейти к уравнению с матрицами ковариации (4.54).

    $$e^+_ke^{+^T}_k=(E-K_kH_k)e^-_ke^{+^T}+K_kv_ke^{+^T}$$ $$=(E-K_kH_k)e^-_k((E-K_kH_k)e^-_k+K_kv_k)^T$$ $$+K_kv_k((E-K_kH_k)e^-_k+K_kv_k)^T$$ $$=(E-K_kH_k)e^-_ke^{-T}_k(E-K_kH_k)^T+(E-K_kH_k)e^-_kv^T_kK^T_k$$ $$+K_kv_ke^{-T}_k(E-K_kH_k)^T+K_kv_kv^T_kK^T_k$$ $$P^+_k=(E-K_kH_k)P^-_k(E-K_kH_k)^T+K_kR_kK^T_k$$

    Отсюда получаем, что матрицу $$K_k$$ необходимо выбрать так, чтобы сумма диагональных элементов (след) матрицы $$P^+_k$$ был минимален.

    $$trace(p^+_k) \rightarrow min_{K_k}$$

    Решение данной задачи определяется посредством дифференцирования функции $$trace(p^+_k)$$ по $$K_k$$ (4.56) Результат дифференцирования получается благодаря применению свойства $$\frac {\delta} {\delta A} (trace(ABA^T))=2AB$$ при условии, что матрица B является симметричной. .

    $$-2(E-K_kH_k)P^-_kH^T_k+2K_kR_k=0$$

    Как следствие, получаем формулу (4.57) для вычисления матрицы $$K_k$$.

    $$K_k=P^-_kH^T_k(R_k+H_kP^-_kH^T_k)^{-1},$$

    где

    $$P^-_k=F_kP^+_{k-1}F^T_k+Q_{k-1}.$$

    Заметим, что посредством несложных математических выкладок можно вывести формулу (4.59) для вычисления

    $$P^+_k=(E-K_kH_k)P^-_k$$

    Подводя итог, более четко выделим схему работы фильтра:

  • Предсказание. Предполагает вычисление априорной оценки состояния $$\tilde{x}^-_k$$ (4.45) и наблюдения $$z_k$$ (4.46).
  • Коррекция. Включает определение матричного коэффициента Кальмана $$K_k$$ (4.57) и построение апостериорной оценки состояния $$\tilde{x}^+_k$$ согласно (4.49).
  • 2.7. Применение фильтра частиц к задаче сопровождения

    Для некоторых практических задач, чтобы получить более точную оценку состояния системы, необходимо уйти от предположения, что шум имеет Гауссово распределение [5]. В этом случае вводится понятие мультимодального распределения Мультимодальным распределением называется распределение, имеющее несколько мод или локальных максимумов. Мультимодальное распределение зачастую представляется смесью нескольких распределений. шума, а для моделирования подобных систем используются фильтры частиц. Фильтры частиц являются более общим подходом к решению задачи сопровождения с применением вероятностных методов [9, 21, 33, 35].

    Алгоритм воспроизведения условной плотности (CONditional DENSity propAGATION, CONDENSATION) [21, 35] – базовый алгоритм фильтрации частиц, на основании которого строится большинство алгоритмов данной группы, применяемых в компьютерном зрении. Поэтому остановимся более детально на рассмотрении схемы работы именно этого алгоритма.

    Предполагая, что система может находиться в состояниях $$X_t={x_1,x_2,...,x_t}$$, в момент времени получим плотность распределения вероятности. Так же, как и для фильтра Кальмана, последовательность наблюдений будем обозначать $$Z_t={z_1,z_2,...,z_t}$$. Наряду с этим введем предположение о том, что состояние $$x_{t}$$ зависит только от предыдущего состояния $$x_{t-1}$$ – условие Марковской цепи. Таким образом, получаем систему с независимым набором наблюдений. Техники фильтрации частиц представляют распределение вероятности в виде коллекции взвешенных выборок – частиц, появление которых регулируется посредством введения весов. Тогда множество $$S_t$$ (4.60) определяет функцию плотности вероятности для состояния $$x_t$$ при заданном наборе наблюдений $$Z_t$$ . $$S_t$$ задает приближенное распределение $$p(x_t \lvert Z_t)$$.

    $$S_t=\lbrace (s^t_i,\pi^t_i),i=\overline{1,N}, \sum_{i=1}^N {\pi^t_i=1} \rbrace$$

    Задача состоит в том, чтобы построить метод восстановления множества $$S_{t}$$ на основании $$S_{t-1}$$. Формально алгоритм можно представить в виде последовательности этапов [9, 21, 35]:

  • Пусть коллекция взвешенных выборок в момент времени $$t-1$$ построена (4.61). $$S_{t-1}=\lbrace (s^{t-1}_i,\pi^{t-1}_i),i=\overline{1,N}, \sum_{i=1}^N {\pi^{t-1}_i=1} \rbrace$$ Дополнительно вычислим интегральные веса согласно (4.62). $$c_i=c_{i-1}+\pi^{t-1}_i,i=\overline{1,N},c_0=0$$
  • Определим n-ый экземпляр выборки $$S_t$$ . Для этого случайным образом выберем число r из отрезка [0,1] и вычислим $$j=arg\,min_i{c_i>r}$$. Отсюда получаем текущую оценку состояния $$s^{t-1}_j$$ .
  • Выполним предсказание следующего состояния. Предсказание (4.63) выполняется аналогично фильтру Кальмана (4.40), разница лишь в том, что нет ограничений, связанных с линейностью системы и видом распределения шума. $$s^{t}_n=F_{t-1}s^{t-1}_j+w_{t-1}$$
  • Выполним коррекцию. Используя текущее наблюдение и его распределение, необходимо установить вес полученного экземпляра согласно (4.64). $$\pi^t_n=p(z_t \lvert x_t=s^t_n)$$
  • Построим множество частиц , повторив шаги 2 – 4 -раз.
  • Нормализуем последовательность весов $$\pi^t_i$$ так, чтобы $$\sum^N_{i=1} {\pi^t_i}=1$$ .
  • Вычислим наилучшую оценку для состояния $$x_t$$, например, как линейную свертку полученного набора экземпляров выборки (4.65). Таким образом, фактически определим некоторую среднюю частицу. $$x_t \sum^N_{i=1} {\pi^t_is^t_i}$$
  • Описанный процесс можно проинтерпретировать графически (рис.4.8) с использованием понятия частицы. На входе итерации алгоритма имеется множество частиц $$\lbrace (s^{t-1}_i,\pi ^{t-1}_i) \rbrace$$ (верхний уровень диаграммы). В результате N-кратного случайного выбора частиц из $$S_{t-1}$$ получается некоторый набор экземпляров (2-ой уровень сверху). Применение шага предсказания приводит к формированию множества оценочных состояний частиц (3-ий уровень), затем для каждой оценки выполняется коррекция на основании имеющихся наблюдений. Как следствие, создается множество частиц $$\lbrace (s^{t}_i,\pi ^{t}_i) \rbrace$$ в следующий момент времени (последний уровень диаграммы).

    (рис 4.8) Итерация алгоритма воспроизведения условной плотности (CONDENSATION, [35])

    3. Контрольные вопросы

  • В каких случаях схема вычитания фона в задаче определения движения работает неэффективно?
  • Докажите, что формула (4.25) справедлива, т.е. образ Фурье от корреляции пары функций равен произведению образу Фурье первой функции на комплексно-сопряженный образ второй функции.
  • Основное свойство аффинного преобразования?
  • Приведите примеры задач, в которых используется модель перспективной проекции (гомография).
  • Докажите справедливость формул (4.58) и (4.59).
  • Получите формулы вычисления матричного коэффициента Кальмана для случая равномерного движения.
  • Приведите схему построения фильтра Кальмана в случае одномерного состояния и наблюдения при некоторых фиксированных значениях $$F_k=F=const$$ и $$H_k=H=const$$.
  • Страницы:

    Введение

    Одни существа идут, другие — следуют за ними.

    Презентацию к лекции Вы можете скачать здесь.

    Движение – смещение одних объектов относительно других. В задачах компьютерного зрения выделяется несколько принципиально разных случаев движения [1]:

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

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

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

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

    В настоящее время выделяется несколько основных групп задач, связанных с понятием движения на наборе последовательных кадров [9]:

  • Определение движения или областей движения. Задача рассматривается в условиях неподвижной камеры и относительно постоянного фона. Наиболее простая задача, т.к. ее решение не требует анализа полученной информации о движении. Подобные задачи, как правило, возникают при разработке охранных систем слежения, в которых необходимо обнаружить несанкционированный доступ на территорию закрытого объекта, где, например, исключено движение в определенное время суток.
  • Поиск движущихся объектов . Данная задача является более сложной по сравнению с предыдущей группой задач, т.к. необходимо не просто определить область, в которой происходит движение, но и выделить движущиеся объекты. Один из возможных подходов к решению основан на определении областей движения и последующем распознавании объектов внутри полученных областей. Задача усложняется, когда требуется построить траекторию движения объекта – восстановить положения объекта на последовательном наборе изображений. В этом случае выполняется либо сопровождение областей движения и распознавание объектов в каждой области полученной последовательности, либо поиск новых объектов на сцене и их дальнейшее сопровождение. Поиск и построение траектории движущихся объектов представляет интерес для многих систем видеонаблюдения, т.к. решение данной задачи обеспечивает возможность анализа поведения объектов (покупателей в торговых центрах, автомобилей на дороге и т.п.).
  • Определение трехмерных свойств объекта из набора его последовательных двумерных изображений . Трехмерные свойства объекта, как правило, необходимы для восстановления трехмерной модели объекта, а как следствие, реконструкции трехмерной сцены. Задача актуальна для приложений компьютерной графики, связанных с созданием реалистичных сцен.
  • В настоящем разделе остановимся на некоторых наиболее известных методах решения задач первой и второй группы.

    1. Определение областей движения

    1.1. Постановка задачи выделения областей движения

    Задача выделения областей движения на видео – одна из классических задач компьютерного зрения. Входом указанной задачи является последовательность кадров $$I_1,I_2, ... ,I_N$$ некоторого видео V:

    $$I_k= \lbrace I_k(x,y),0 \leqslant x < width,\, 0 \leqslant y < height \rbrace ,k=\overline{1,N},$$

    где width – ширина кадра, height – высота кадра, а $$I_k(x,y)$$ в общем случае представляет собой вектор фиксированной размерности В случае полутонового изображения $$I_k(x,y)$$ – вектор размерности 1, который представляет интенсивность пикселя (x,y) . Если исходное изображение $$I_k$$ задано в формате RGB, то $$I_k(x,y)$$ – вектор размерности 3, каждая компонента которого определяет интенсивность пикселя по соответствующему каналу. .

    Решением настоящей задачи является совокупность областей изображения для каждого кадра видео, в которых происходит движение одного или нескольких объектов. Таким образом, в результате обработки видео необходимо сформировать набор бинарных изображений, в которых белые пиксели (интенсивность 255) соответствуют пикселям, принадлежащим движущимся объектам, а черные (интенсивность 0) – пикселям фона (4.2).

    $$M_k(x,y)=\begin{cases} 255,\text{(x,y)-пиксель объекта}\\ 0,\text{(x,y)-пиксель фона} \end{cases},k=\overline{1,N}$$

    1.2. Вычитание фона

    Наиболее простой подход к решению данной задачи состоит в том, чтобы использовать механизм вычитания фона из кадра видео (background subtraction) [6 – 9]. Процедура вычитания предполагает, что для данного видео построена модель фона (4.3):

    $$F= \lbrace F(x,y),0 \leqslant x < width,\, 0 \leqslant y < height \rbrace,$$

    также возможно существует механизм обновления модели фона с течением времени. Для одноканального (в оттенках серого) изображения, т.е. когда $$I_k(x,y),F(x,y) \in \lbrace 0,...,255 \rbrace,k=\overline{1,N},$$ процедуру вычитания можно разбить на два этапа:

  • Вычитание фонового изображения из текущего кадра видео. Данный шаг включает в себя попиксельное вычитание интенсивностей кадра видео и фонового изображения (4.4). $$D_k(x,y)=abs(I_k(x,y)-F(x,y)),k=\overline{1,N}$$
  • Отбор пикселей, принадлежащих фону и объекту, – построение бинарного изображения (маски). Считается, что пиксель принадлежит объекту и имеет белый цвет в маске, если разность интенсивности фона и текущего кадра для данного пикселя превышает некоторое пороговое значение в противном случае, принимается, что пиксель принадлежит фону (4.5). $$M_k(x,y)=\begin{cases} 255,D_k(x,y)\geqslant \tau\\ 0,D_k(x,y)< \tau \end{cases},k=\overline{1,N}$$
  • Дополнительно к указанным операциям с целью повышения качества поиска может выполняться, например, фильтрация кадров исходного потока видеоданных, либо фильтрация бинарного видео, также могут применяться морфологические операции к полученному отсечению с целью удаления шумов [8]. Если имеется цветное изображение, то его всегда можно преобразовать в оттенки серого.

    Качество определения положения движущихся областей посредством вычитания фона во многом зависит от качества построенной модели фона. Множество всех техник вычитания фона подразделяется на две группы в зависимости от механизма построения фонового изображения:

    1. Нерекурсивные. Нерекурсивные методы обновляют модель фона для текущего кадра на основании информации об интенсивностях пикселей некоторого набора предшествующих моделей фона [6] (или кадров) и текущего кадра. К наиболее распространенным нерекурсивным методам относятся следующие методы:

  • Метод вычитания текущего и предыдущего кадра . Согласно данному методу считается, что для кадра $$I_k$$ модель фона $$F_k$$ совпадает с предыдущим кадром, т.е. $$I_k(x,y)$$. Тогда на первом этапе алгоритма вычитания фона вычисляется разница пары последовательно идущих кадров (4.6). $$D_k(x,y)=abs(I_k(x,y)-F(x,y)) = \\ = abs(I_k(x,y)-I_{k-1}(x,y)),k=\overline{2,N} $$
  • Метод усреднения определенного количества предшествующих кадров . Обозначим количество кадров, по которым будет выполняться построение модели фона, как s. Тогда для кадра $$I_k$$ модель фона $$F_k$$ определяется в соответствии с формулой (4.7). $$F_k(x,y)=\frac 1 s \sum^{s-1}_{j=0} {I_{k-j}(x,y)}$$
  • Метод определения медианы фиксированного количества предшествующих кадров . Пусть s – число кадров, на основании которых будет обновляться модель фона. Тогда модель фона вычисляется согласно формуле (4.8). $$F_k(x,y)=med_{j=\overline{0,s-1}} \lbrace I_{k-j}(x,y) \rbrace$$
  • Основное преимущество нерекурсивных методов – простота реализации и скорость обновления моделей фона при переходе от кадра к кадру. При этом заметим, что качество работы методов данной группы в значительной степени зависит от скорости объектов. Медленно перемещающиеся объекты, как правило, обнаруживаются плохо. Более того, приведенные методы не дают качественный результат при изменении света на сцене, либо при наличии динамического фона (листва деревьев, струящаяся вода и т.п.). Чтобы сгладить влияние указанных эффектов, обновленная модель фона для кадра $$I_k$$ представляется выпуклой оболочкой модели фона $$F_{k-1}$$ с текущим изображением (4.9). Такая процедура называется $$\alpha$$-смешиванием.

    $$F_k(x,y)=\alpha I_k (x,y)+(1+\alpha)F_{k-1}(x,y)$$

    2. Рекурсивные. Рекурсивные методы для обновления модели фона используют информацию об интенсивностях пикселей только текущего кадра. К методам данной группы относятся гистограммный метод, метод представления модели фона смесью Гауссовых распределений (Gaussian mixture model) [2, 3, 9], метод "шифровальной" книги (codebook) [10, 11], метод извлечения визуального фона (Visual Background Extractor, ViBe) [12]. Рассмотрим некоторые из перечисленных методов.

  • Гистограммный метод. Идея гистограммного метода состоит в том, что все цветовое пространство разбивается на отдельные бины (для полутонового изображения пространство представляет отрезок изменения интенсивности, в случае цветного изображения – трехмерный куб). Для каждого изображения в последовательности выполняется построение гистограммы. Осуществляется проход по всем пикселям изображения, в зависимости от того, какая интенсивность/цвет наблюдается в пикселе, увеличивается на единицу величина соответствующего бина гистограммы. Принимается, что пиксели, составляющие некоторый бин, принадлежат фону, если величина данного бина меньше фиксированного порогового значения, в противном случае, считается, что они принадлежат объекту. Основная проблема применения гистограммного метода состоит в необходимости использования дополнительной памяти, а также в выполнении большого количества операций обращения к памяти в процессе реализации.
  • Смесь Гауссовых распределений. При построении фона с использованием данного метода считается, что для любого пикселя $$(x_0,y_0)$$ изображения $$I_k$$ известна история изменения его интенсивности/цвета на всех предшествующих кадрах $$\lbrace X_1,X_2,...,X_k \rbrace = \lbrace I_j (x_0,y_0),j=\overline{1,k} \rbrace$$ . Тогда вероятность того, что наблюдается значение $$X_k$$, может быть представлена смесью из Гауссовых распределений (4.10).
  • $$P(x_k)=\sum^s_{j=1} {\omega^k_jN(x_k\lvert \mu^k_j,\Sigma^k_j)},$$

    где $$\omega^k_j$$ – вес j-ого распределения Гаусса для кадра с номером $$k,\mu^k_j$$, – математическое ожидание, $$\Sigma^k_j)$$ – среднеквадратичное отклонение, $$N(x_k\lvert \mu^k_j,\Sigma^k_j)$$ – функция плотности нормального распределения (4.11).

    $$N(x_k\lvert \mu^k_j,\Sigma^k_j)=\frac 1 {(2\pi)^{\frac D 2}\lvert \Sigma^k_j \rvert ^{\frac 1 2}}e^{-\frac 1 2 (x_k-\mu^k_j)^T(\Sigma^k_j)^{-1}(x_k-\mu^k_j)}$$

    Предполагается, что компоненты цвета независимы и имеют одинаковое среднеквадратичное отклонение. Поэтому матрица ковариации имеет вид $$\Sigma^k_j=(\sigma^k_j)^2E$$, где E – единичная матрица.

    Указанное предположение позволяет снизить вычислительную трудоемкость метода за счет отсутствия необходимости вычислять матрицу, обратную к матрице ковариации $$\Sigma^k_j$$ в (4.11). Таким образом, задано распределение наблюдаемых значений цвета для каждого пикселя. Новое значение будет представляться одной из основных компонент построенной смеси Гауссовых распределений и использоваться для обновления параметров модели. Распределения сортируются в порядке уменьшения величины $$r^k_j=\frac {\omega^k_j} {\sigma^k_j}$$ . Такая сортировка предполагает, что пиксель фона отвечает распределению с большим весом и малой дисперсией.

    Принимается, что первые $$B^k$$ распределений, удовлетворяющих условию (4.12), соответствуют распределению цвета фоновых пикселей.

    $$B^k=argmin_b \lbrace \sum^b_{j=1} {\omega^k_j>T}\rbrace,$$

    где T – некоторое пороговое значение, параметр модели. Когда приходит очередной кадр $$I_{k+1}$$, для каждого пикселя изображения выполняется тест, который позволяет определить с использованием расстояния Махаланобиса, какому распределению соответствует полученное значение (4.13).

    $$\sqrt {(x_{k+1}-\mu^k_j)^T(\sigma^k_j)^{-1}(x_{k+1}-\mu^k_j)}< 2,5\sigma^k_j$$
  • Если нашлось соответствующее распределение Гаусса, то в зависимости от того, определяет ли оно распределение фоновых пикселей (входит в группу из $$B^k$$ распределений) или нет, текущий пиксель классифицируется как фоновый, либо как принадлежащий объекту.
  • Если не обнаружилось ни одного распределения, удовлетворяющего условию (4.13), то считается, что пиксель принадлежит объекту.
  • На основании такого правила формируется двумерная маска.

    Чтобы обработать следующий кадр, необходимо обновить параметры распределений: математическое ожидание $$\mu^k_j$$ и среднеквадратичное отклонение $$\sigma^k_j$$ . В зависимости от того, нашлось ли соответствующее распределение для цвета текущего пикселя, обновление выполняется по-разному.

  • Соответствие обнаружено. Тогда весовые коэффициенты, составляющие смесь Гауссовых распределений, которым соответствует $$X_{k+1}$$ и параметры распределений пересчитываются согласно формулам (4.14) – (4.16). $$\omega^{k+1}_j=(1-\alpha)\omega^{k}_j+\alpha,$$ $$\mu^{k+1}_j=(1-\rho)\mu^{k}_j+\rho X_{k+1},$$ $$(\sigma^{k+1}_j)^2=(1-\rho)(\sigma^{k}_j)^2+\rho(x_{k+1}-\mu^{k+1}_j)(x_{k+1}-\mu^{k+1}_j)^T,$$ где $$\alpha$$ – заданная константа, $$\rho=\alpha N(x_k\lvert \mu^k_j,\Sigma^k_j).$$. Для всех распределений, которым $$X_{k+1}$$ не соответствует, параметры не изменяются, только пересчитываются коэффициенты $$\omega^{k+1}_j$$ согласно (4.17). $$\omega^{k+1}_j=(1-\alpha)\omega^{k}_j$$
  • Соответствие не найдено. В данном случае крайнее (в смысле введенного отношения порядка) распределение Гаусса замещается распределением с новыми параметрами. Математическое ожидание выбирается равным текущему значению цвета пикселя $$\mu^{k+1}_s=X_{k+1}$$ дисперсия $$(\sigma^{k+1}_j)^2$$ максимально возможной, а вес $$\omega^{k+1}_s$$ минимально допустимым.
  • В завершении отметим, что количество распределений определяется сложностью фона и имеющимися вычислительными мощностями (в [2] предлагается использовать значение в пределах от 3 до 5). Начальная инициализация параметров распределений может выполняться с использованием метода k-средних [2], либо EM-алгоритма (Expectation Maximization) [4, 32].

    На рис.4.1 показан пример применения метода вычитания фона с представлением модели фона смесью Гауссовых распределений.

    (рис 4.1) Пример применения метода вычитания фона с использованием метода представления фона смесью Гауссовых распределений

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

    Метод извлечения визуального фона (Visual Background Extractor , ViBe) [12]. В соответствии с данным методом модель фона на кадре с номером $$ k$$ представляется набором множеств $$M^k(p)=\lbrace v_1,v_2,...,v_N\rbrace$$ для всех пикселей $$p=(x,y)$$ где $$v_i$$ – интенсивность/цвет пикселя (в общем случае вектор).

    Для классификации пикселя в цветовом пространстве строится сфера $$S_R(v(x))$$ радиуса R и определяется количество векторов множества $$M(p)$$ которые попадают вовнутрь построенной сферы, $$K=\lvert S_R(v(p)) \cap M(p) \rvert$$. Если $$K>T_{min}$$ , где $$T_{min}$$ – фиксированное пороговое значение, то принимается, что пиксель принадлежит фону, в противном случае, объекту.

    На начальном этапе необходимо выполнить инициализацию множеств $$M^0(p)$$ для всех пикселей следующим образом:

    $$M^0(p)=\lbrace v^0(y),y \in N_G(p) \rbrace,$$

    где $$N_G(p)$$ – окрестность пикселя размера 3x3 (9 клеток, включая текущий пиксель), y выбирается $$N$$ раз случайным образом. Обновление модели фона для кадра $$I_k$$ выполняется в два шага:

  • Если p проклассифицирован как пиксель фона, то из множества $$M^k(p)$$ случайно выбирается компонента, которая заменяется значением $$v(p)$$.
  • Случайным образом выбирается один соседний пиксель из окрестности $$N_G(p)$$ , для которого выполняется предыдущий шаг.
  • Множество методов построения моделей фона не ограничивается набором методов, представленных в этом разделе. В данном направлении ведутся активные исследования до настоящего момента. Поэтому здесь представлены лишь принципиально разные подходы к решению задачи построения фоновых моделей.

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

    Еще один простейший способ оценить движение на нескольких изображениях – перебрать все возможные варианты смещений изображений или их отдельных фрагментов (translational alignment, [5, гл. 8.1]). Для этого первоначально необходимо выбрать метрику для оценки степени сходства фрагментов. Как следствие, исходная задача определения движения может быть сведена к минимизации функции ошибки по всем возможным направлениям смещения (рис.4.2).

    (рис 4.2) Пример перебора вариантов по четырем принципиальным направлениям

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

  • 1. Квадратичная функция ошибки (4.19). $$E(u,v)=\sum_i {(I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i))^2},$$ где $$(u,v)$$ – вектор смещения изображения или фрагмента.
  • 2. Грубые оценочные функции ошибки. Общий вид функций данной группы определяется формулой (4.20). При этом $$\rho (e)$$ – функция, которая имеет меньшую степень роста по сравнению квадратичной функцией. $$E(u,v)=\sum_i {\rho (I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i))},$$ Примером грубой функции ошибки может служить сумма абсолютных разностей (4.21). $$E(u,v)=\sum_i {\lvvert I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i) \rvert},$$
  • 3. Взвешенная квадратичная функция ошибки (4.22). Весовые коэффициенты – дискретные функции, которые принимают нулевые значения за пределами изображения. $$E(u,v)=\sum_i {\omega_{k-1} (x_i,y_i)\omega_k(x_i+u,y_i+v) \cdot (I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i))^2},$$
  • 4. Кросс-корреляция пары выровненных изображений (4.23). Заметим, что при выборе кросс-корреляции в качестве функции ошибки выполняется не минимизация, а максимизация по всем возможным направлениям смещения. $$E(u,v)=\sum_i {I_{k-1} (x_i,y_i) \cdot I_k(x_i+u,y_i+v)$$
  • В общем случае если задан набор цветных изображений, то всегда можно выполнить их конвертирование в полутоновые, либо ввести дополнительную сумму по числу каналов, и использовать одну из предложенных функций ошибки.

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

    Еще один способ ускорения вычисления при реализации переборной схемы – использование подхода, основанного на вычислении быстрого преобразования Фурье (БПФ) [30]. Метод опирается на тот факт, что модули образов Фурье исходного и смещенного сигналов связаны уравнением (4.24).

    $$F(I_k(x_j+u,y_j+v))=F(I_k(x_i,y_i))e^{-i(\begin{bmatrix} u \\ v \end{bmatrix},\overline{\omega})}=I_k(\overline{\omega})e^{-i(\begin{bmatrix} u \\ v \end{bmatrix},\overline{\omega})},$$

    где $$\overline{\omega}=\begin{bmatrix} \omega_x \\ \omega_y \end{bmatrix}$$ – вектор частот двумерного БПФ, $$I_k(\overline{\omega})=F(I_k(x_i,y_i))$$ – образ Фурье сигнала $$I_k(x_i,y_i)$$. Допустим, что в качестве функции ошибки выбрана кросс-корреляция, вычислим БПФ от этой функции (4.25). В результате получаем преобразование Фурье от свертки функций. Исходя из свойств преобразования, образ Фурье от полученной свертки функций – произведение образа $$I_{k-1}(\overline{\omega})$$ первой функции на комплексно- сопряженный образ $$I^*_k(\overline{\omega})$$ второй функции.

    $$F(E(u,v))=F \Biggl( \sum_i {I_{k-1} (x_i,y_i) \cdot I_k(x_i+u,y_i+v) \Biggr) = F(I_{k-1}(u,v)+I_k(u,v))=I_{k-1}(\overline{\omega})I^*_k(\overline{\omega})$$

    Отсюда эффективность оценивания кросс-корреляции с использованием преобразования Фурье очевидна. Достаточно вычислить образы Фурье для исходной пары изображений с помощью БПФ (сложность составляет $$O(NM\,log\, NM))$$ для изображения размера $$N \times M$$ , перемножить полученные образы (за время, не превышающее $$O(N^2M^2))$$) , затем вычислить обратное преобразование Фурье. В результате получится матрица значений кросс-корреляционной функции, из которой необходимо выбрать максимальное значение. Приведенный пример не является единственным примером ускорения вычислений функции ошибки с помощью БПФ [5].

    Для некоторых практических приложений, в частности, для стабилизации видео необходима субпиксельная точность определения смещения – смещение в дробных координатах. Один из наиболее известных подходов состоит в том, чтобы использовать процедуру уточнения направления смещения. Впервые была введена при разработке алгоритма сопровождения Лукаса-Канаде, о котором речь пойдет в одном из следующих разделов. Процедура основана на применении метода градиентного спуска для минимизации квадратичной функции ошибки в смещенной точке – функции энергии (4.26), которая приближается разложением в ряд Тейлора.

    $$E(u+\Delta u,v+\Delta v)=\sum_i {(I_k(x_i+u+\Delta u,y_i+v+\Delta v)-I_{k-1} (x_i,y_i))^2} \approx$$ $$\approx \sum_i {\Biggl( I_k(x_i+u,y_i+v)+\nabla I_k(x_i+u,y_i+v)}\begin{bmatrix} \Delta u \\ \Delta v \end{bmatrix}-I_{k-1} (x_i,y_i) \Biggr) ^2=$$ $$= \sum_i {(\nabla I_k(x_i+u,y_i+v)\begin{bmatrix} \Delta u \\ \Delta v \end{bmatrix}+e_i)^2}=$$ $$\sum_i {(J_k(x_i+u,y_i+v)\begin{bmatrix} \Delta u \\ \Delta v \end{bmatrix}+e_i)^2} \rightarrow min_{(u,v)},$$ $$e_i=I_k(x_i+u,y_i+v)-I_{k-1}(x_i,y_i),$$ $$J_k(x_i+u,y_i+v)=\nabla I_k(x_i+u,y_i+v)=[\frac {\deltaI_k} {\delta x},\frac {\deltaI_k} {\delta y}]\rvert_{(x_i+u,y_i+v)}$$

    Полученная задача минимизации (4.27) сводится к решению системы линейных алгебраических уравнений вида (4.28). Доказательство данного факта можно найти в [31].

    $$A\begin{bmatrix} \Delta u \\ \Delta v \end{bmatrix}=b$$

    где $$A=\sum_i {J^T_k(x_i+u,y_i+v)J_k(x_i+u,y_i+v)},b=-\sum_i {e_iJ^T_k(x_i+u,y_i+v)}.$$ .

    Эффективность решения задачи достигается за счет того, что для текущего изображения градиенты $$J_k(x_i+u,y_i+v)$$ в смещенной точке можно приближенно заменить градиентами предыдущего изображения $$J_{k-1}(x_i,y_i)$$ в исходной точке, т.е. $$J_k(x_i+u,y_i+v)\approx J_{k-1}(x_i,y_i)$$.

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

    1.4. Параметрические модели движения

    Многие прикладные задачи, такие, как склеивание изображений в панораму, стабилизация видео, требуют построения более сложных моделей движения, т.к. движение не описывается только линейным сдвигом. Поэтому рассматриваются пространственные поля смещений, и выполняется построение параметрических моделей движения (parametric motion, [5, гл. 8.2]). В параметрических моделях вместо постоянного вектора смещения$$[u,v]^T$$ рассматривается параметризованное поле смещений $$x'(x,y;p)$$, где p-вектор параметров. Приведем некоторые примеры параметризованных моделей:

  • Линейный сдвиг (4.29). $$x'(x,y;p)=\begin{bmatrix} E p \\ {0^T} 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix},$$ где $$p=\begin{bmatrix} {p_1} \\ {p_2} \end{bmatrix}$$ – вектор сдвига.
  • Поворот со сдвигом (4.30). $$x'(x,y;p)=\begin{bmatrix} {cos\,\phi} {-sin\,\phi} {p_1} \\ {sin\,\phi} {cos\,\phi} {p_2} \\ 0 0 1 \end{bmatrix}\begin{bmatrix} x \\ y \\ 1 \end{bmatrix},$$ где $$R=\begin{bmatrix} {cos\,\phi} {-sin\,\phi} \\ {sin\,\phi} {cos\,\phi} \end{bmatrix}$$ – ортонормированная матрица поворота, т.е. $$RR^T=1$$ и $$det\,R=1$$ , $$p=(p_1,p_2,\phi)$$ – набор параметров модели движения.
  • Преобразование подобия (4.31). $$x'(x,y;p)=\begin{bmatrix} a -b {p_1} \\ b a {p_2} \\ 0 0 1 \end{bmatrix}\begin{bmatrix} x \\ y \\ 1 \end{bmatrix},$$ где $$sR=\begin{bmatrix} a -b \\ b a \end{bmatrix}$$ – матрица поворота, умноженная на коэффициент масштабирования $$s,\,p=(p_1,p_2,a,b)$$ – набор параметров модели движения.
  • Аффинное преобразование (4.32). $$x'(x,y;p)=\begin{bmatrix} a_{00} a_{01} a_{02} \\ a_{10} a_{11} a_{12} \end{bmatrix}\begin{bmatrix} x \\ y \\ 1 \end{bmatrix}$$
  • Перспективная проекция или гомография (4.33). $$x'(x,y;p)=\begin{bmatrix} h_{11} h_{12} h_{13} \\ h_{21} h_{22} h_{23} \\ h_{31} h_{32} h_{33} \end{bmatrix}\begin{bmatrix} x \\ y \\ 1 \end{bmatrix}=H\overline{x}$$
  • Аналогично процедуре уточнения направления, описанной в разделе 1.3, можно выполнить определение направления движения для любой параметрической модели (4.35), расписав квадратичную функцию ошибки (4.34) с помощью разложения в ряд Тейлора по вектору параметров и выполнив ее минимизацию.

    $$E(p+\Delta p)=$$ $$\sum_i {(I_k(x'(x_i,y_i;p+\Delta p))-I_{k-1}(x_i,y_i))^2}\Approx$$ $$\sum_i {(I_k(x'_i) +I_k(x'_i)\Delta p-I_{k-1}(x_i,y_i))^2}=$$ $$\sum_i {(J_k(x'_i)\Delta p+e_i)^2} \rightarrow min_{\Delta p},$$ $$e_i=I_k(x'_i)-I_{k-1}(x_i,y_i),$$ $$J_k(x'_i)=\frac {\delta I_k} {\delta p}=\nabla I_k(x'_i)\frac {\delta x'_i} {\delta p}\rvert_{(x_i,y_i)}.$$

    1.5. Многоуровневое движение

    Во многих случаях визуальное движение вызвано смещением небольшого количества объектов, находящихся на разной глубине изображения [5]. Поэтому движение пикселей можно описать более эффективно, если сгруппировать их в слои [14], и отслеживать многоуровневое движение (layered motion, [5, гл. 8.5]) построенных слоев, например, с помощью параметрических моделей [15].

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

    1. Как представить слой? Слой определяется набором из трех карт (матриц):

  • матрица интенсивности слоя (в компьютерной графике текстурная карта);
  • $$\alpha$$карта, которая определяет прозрачность слоя в каждом точке изображения;
  • карта скоростей описывает изменение положения точек с течением времени.
  • Ниже (рис.4.3) показан пример представления кадров в виде набора слоев. Изображения в строке (а) соответствуют фоновому слою, в строке (b) – слою, содержащему объект (рука), в (с) – последовательность изображений, восстановленная на основании выделенных слоев.

    (рис 4.3) Пример представления изображения в виде набора слоев [15]

    Предположим, что $$E_0(x,y)$$ – матрица интенсивности фонового слоя, $$E_1(x,y)$$– слоя, содержащего объект, $$\alpha_1(x,y)$$ – $$\alpha$$-карта слоя $$E_1(x,y)$$. Поскольку слои перекрываются, то для представленного примера восстановить можно исходное изображение согласно (4.36). $$I_1(x,y)=E_0(x,y)(1-\alpha_1(x,y))+E_1(x,y)\alpha_1(x,y)$$

    Отметим, что данную процедуру можно распространить на случай большего количества слоев (4.37). $$I_k(x,y)=I_{k-1}(x,y)(1-\alpha_k(x,y))+E_k(x,y)\alpha_k(x,y)$$

    2. Как разбить множество пикселей на слои и определить движение каждого слоя? Процедура разбиения на слои состоит из нескольких шагов:

  • Оценивание движения с использованием оптического потока и аффинной модели для набора неперекрывающихся блоков (рис.4.4, а), полученных в результате разбиения исходного изображения. Данная операция выполняется независимо для каждого кадра из подмножества изображений.
  • Кластеризация полученных оценок с помощью метода k-средних (рис.4.4, b). Кластеризация также осуществляется независимо для каждого изображения. Каждый кластер определяет сегмент, отвечающий некоторому слою изображения.
  • Применение медианной фильтрации для получения смешанных слоев, устойчивых к незначительным изменениям интенсивности, а также для выявления перекрытий между слоями (рис.4.4, c).
  • (рис 4.4) Пример выделения слоев [14]

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

    1.6. Вычисление оптического потока

    Другой распространенный подход к решению задачи детектирования областей движения – вычисление оптического потока (optical flow) [5, 9, 13]. Оптический поток позволяет определить смещение каждой точки кадра, т.е. построить поле скоростей.

    Применение данного подхода требует выполнения двух основных условий:

  • яркость каждой точки объекта не изменяется с течением времени;
  • ближайшие точки, принадлежащие одному объекту, в плоскости изображения двигаются с похожей скоростью.
  • Рассмотрим схему работы метода вычисления оптического потока. Допустим, что имеется непрерывно изменяющееся во времени изображение. Обозначим через $$I_k(x,y)$$ яркость пикселя с координатами $$(x,y)$$ в момент времени k. Считая динамическое изображение функцией положения и времени, яркость пикселя в следующий момент времени можно выразить с помощью разложения в ряд Тейлора:

    $$I_{k+\Delta x}(x+\Delta x,y+\Delta y)=I_k(x,y)+\frac {\delta I_k} {\delta x}\Delta x+\frac {\delta I_k} {\delta y}\Delta y+\frac {\delta I_k} {\delta k}\Delta k+O(\delta^2)$$

    где $$\frac {\delta I_k} {\delta x}\Delta x,\frac {\delta I_k} {\delta y}\Delta y,\frac {\delta I_k} {\delta k}\Delta k$$ – частные производные Под частными производными в силу дискретности представления изображения понимается центральный разностный оператор. функции $$I_k$$. Предполагается, что за небольшой промежуток времени пиксель смещается незначительно, т.е. $$(\Delta x, \Delta y) \rightarrow 0$$ при $$\Delta k \rightarrow 0$$, значит можно считать, что яркость почти не изменяется, поэтому остаточным членом ряда можно пренебречь:

    $$I_{k+\Delta x}(x+\Delta x,y+\Delta y)=I_k(x,y)$$ $$-\frac {\delta I_k} {\delta k}= \frac {\delta I_k} {\delta x}\frac {\delta x} {\delta k}+\frac {\delta I_k} {\delta y}\frac {\delta y} {\delta k}$$

    Уравнение (4.39) называется уравнением оптического потока . Дальнейшая цель – построить вектор скорости пикселя $$c=(\frac {\delta x} {\delta k},\frac {\delta y} {\delta k})=(u,v)$$ . Для полноты постановки задачи вводится условие гладкости изменения скорости. В результате исходная задача сводится к задаче минимизации квадратичной ошибки $$(\frac {\delta I_k} {\delta x}u+\frac {\delta I_k} {\delta y}v+\frac {\delta I_k} {\delta k}k)^2$$ при наличии ограничений в виде равенств $$u^2_x+u^2_y=0$$ и $$v^2_x+v^2_y=0$$ (условия гладкости). Получаем задачу математического программирования. Процедура минимизации применяется к каждому пикселю текущего изображения, в результате чего обеспечивается построение поля векторов смещения всех пикселей.

    2. Методы сопровождения объектов

    2.1. Постановка задачи сопровождения объектов

    Сопровождение (трекинг) движущихся объектов – это один из составляющих компонентов многих систем реального времени таких, как системы слежения, анализа видео и других. Входными данными любого алгоритма сопровождения является последовательность изображений (кадров видео) $$I_1,I_2,...,I_N$$(4.1) с нарастающим объемом информации, которую необходимо обрабатывать и анализировать.

    Задача сопровождения состоит в том, чтобы построить траектории движения целевых объектов на входной последовательности кадров. Допустим, что положение объекта на изображении с номером обозначается $$P_k$$ . Тогда траекторией движения объекта называется последовательность его положений $$P_s,P_{s+1},...,P_{s+l-1}$$ , где s – номер первого кадра, на котором был обнаружен объект, l – количество кадров последовательности, где наблюдается объект. Заметим, что в зависимости от метода сопровождения положение объекта может определяться по- разному (координаты и размер сторон окаймляющего прямоугольника, координаты центра масс контура и т.п.).

    2.2. Классификация методов сопровождения

    Существует несколько категорий методов сопровождения объектов [16]:

  • Методы сопровождения особых точек (point tracking). В таких методах принимается, что положение объекта определяется расположением набора характерных точек. Один и тот же объект на последовательных кадрах представляется наборами соответствующих пар точек. Данная группа методов разделяется на две подгруппы:

    Детерминистские методы [17] используют качественные эвристики движения (небольшое изменение скорости, неизменность расстояния в трехмерном пространстве между парой точек, принадлежащих объекту), по существу задача сводится к минимизации функции соответствия наборов точек. Методы, основанные на вычислении плотного и разреженного оптического потока, а также методы сопоставления (matching) дескрипторов ключевых точек являются типичными представителями детерминистских методов.

    Вероятностные методы используют подход, основанный на понятии пространства состояний. Считается, что движущийся объект имеет определенное внутреннее состояние, которое измеряется на каждом кадре. Чтобы оценить следующее состояние объекта, требуется максимально обобщить полученные измерения, т.е. определить новое состояние при условии, что получен набор измерений для состояний на предыдущих кадрах. Типичным примерами таких методов являются методы на базе фильтра Кальмана [6, 18 – 20] и фильтра частиц (particle filter) [21, 22, 33].

  • Методы сопровождения компонент (kernel tracking). Под компонентой понимается форма объекта. В простейшем случае компонента может быть представлена шаблоном прямоугольной или овальной формы, в более сложных – трехмерной моделью объекта, спроецированной на плоскость изображения. Как правило, методы данной группы применяются, если движение определяется обычным смещением, поворотом или аффинным преобразованием. Трекинг компонент – итеративная процедура локализации, основанная на максимизации некоторого критерия подобия. На практике реализуется с использованием сдвига среднего (mean shift) [8, 23] и его непрерывной модификации (Continuous Adaptive Mean Shift, CAM Shift) [8, 24].
  • Методы сопровождения силуэта (silhouette tracking). Силуэт может быть задан контуром, либо набором связанных простых геометрических примитивов. Задача трекеров силуэта состоит в том, чтобы на каждом кадре определить область, в которой находится объект, с использованием модели его силуэта, построенной на основании предшествующих кадров. Можно выделить две группы методов:
  • Методы сопоставления и сопровождения фрагментов изображения, содержащих объект. При сопоставлении естественным образом должна быть введена мера сходства пары областей, ограниченных контуром. На практике используется расстояние Хаусдорфа, норма $$L_2$$ , также значение оценочной функции может вычисляться с использованием нескольких мер, включая, например, кросс-корреляцию, расстояние Бхачатария. Определение положения фрагмента на следующем кадре последовательности выполняется посредством вычисления оптического потока (раздел 4.6) для внутренних точек области.
  • Методы сопровождения контура. Позволяют прогнозировать положение контура на следующем кадре. Первый подход состоит в использовании моделей пространства состояний (по типу фильтра Кальмана), второй – в минимизации функции энергии контура с использованием прямых техник таких, как градиентный спуск (раздел 4.3).
  • 2.3. Методы сопоставления ключевых точек

    В настоящее время детерминистские методы сопровождения особых точек чаще остальных используются на практике. Необходимым условием использования этих методов является наличие выделенных точек, которые наилучшим образом представляют объект. Выделение особых точек выполняется с использованием специальных детекторов и дескрипторов Обзор существующих детекторов и дескрипторов приведен в предыдущей лекции. . Если детектор находит положение особой точки на изображении, то дескриптор дополнительно строит вектор признаков, характерных для полученной точки.

    Схема сопровождения особых точек может быть представлена в виде последовательности действий [34]:

  • Поиск особых точек на предыдущем (рис.4.5, слева) и текущем (рис.4.5, справа) кадрах посредством выбранного детектора (LoG, DoG, SURF, SIFT, ORB и др.).
  • Вычисление дескрипторов для полученного набора точек (SURF, SIFT, GLOH, DAISY, BRIEF и т.д.) – n-мерных векторов-описателей.
  • Сопоставление (matching) дескрипторов, полученных на текущем и предыдущем кадре. По существу необходимо построить преобразование точек n-мерного пространства дескрипторов. Прямой перебор дескрипторов предыдущего изображения и поиск наиболее близких дескрипторов текущего не является эффективным из- за большого количества ложных соответствий (рис.4.6). Близость дескрипторов определяется посредством вычисления некоторой метрики (норма $$L_2$$ , кросс-корреляция и т.п.).
  • (рис 4.5) Пример исходного изображения объекта (слева) и объекта внутри некоторой сцены (справа) (рис 4.6) Набор соответствий до применения алгоритма RANSAC

    После сопоставления для отсечения ложных срабатываний (рис.4.7) применяется алгоритм RANSAC (RANdom SAmple Consensus, [6, 7]).

    (рис 4.7) Набор соответствий после применения алгоритма RANSAC

    RANSAC – это общий метод, который используется для оценки параметров модели на основании случайных выборок. При сопоставлении модель представляет собой матрицу преобразования (гомография). На входе алгоритма имеется два множества дескрипторов, полученных на предыдущем и текущем изображении. Схема работы RANSAC состоит из многократного повторения трех этапов:

  • Выбор точек и построение параметров модели. Из входных множеств дескрипторов выбираются случайным образом без повторений наборы фиксированного размера. На основании полученных наборов строится матрица преобразования.
  • Проверка построенной модели. Для каждого дескриптора предыдущего кадра находится проекция на текущем кадре и выполняется поиск наиболее близкого дескриптора из множества дескрипторов текущего кадра. Дескриптор помечается как выброс, если расстояние между проекцией и соответствующим дескриптором текущего изображения больше некоторого порога.
  • Замещение модели. После проверки всех точек проверяется, является ли построенная модель лучшей среди набора предшествующих моделей.
  • В результате применения RANSAC строится наилучшая матрица гомографии. Вычислив перспективную проекцию набора дескрипторов предыдущего кадра, достаточно выполнить проход по всем соответствиям, полученным в процессе перебора, и проверить, является ли соответствующий дескриптор текущего кадра достаточно близким к проекции дескриптора предыдущего кадра. Если не является, то пара отбрасывается.

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

    2.4. Методы сопровождения объектов, основанные на вычислении оптического потока

    С каждым пикселем изображения можно связать некоторый вектор скорости, который определяет, какое расстояние "прошел" пиксель в течение временного промежутка между предыдущим и текущим кадрами. Такая конструкция, построенная для каждого пикселя изображения, представляет собой плотный оптический поток. Одним из наиболее известных алгоритмов сопровождения, основанных на вычислении плотного оптического потока, является метод Хорна (Horn-Schunk method) [8, 25]. По существу данный метод предполагает решение системы дифференциальных уравнений, полученных в задаче математического программирования (раздел 1.6), с помощью итерационных методов для построения компонент вектора скорости в каждой точке изображения.

    Другой класс алгоритмов, использующих плотный оптический поток, – алгоритмы сопоставления блоков (block matching algorithms, [8]). Идея состоит в том, что предыдущее и текущее изображения разбиваются на блоки, как правило, квадратные и перекрывающиеся, а затем определяется движение этих блоков посредством сопоставления. Поскольку алгоритмы сопоставления работают с блоками, то изображение поля скоростей обычно имеет меньшее разрешение по сравнению с исходным изображением.

    Алгоритм Лукаса-Канаде (Lucas-Kanade) основан на вычислении разреженного оптического потока, т.е. на построении векторного поля скоростей для выделенного набора точек. Задача слежения без учета аффинных искажений с помощью данного алгоритма сводится к поиску оптического потока в особых точках [26]. Впоследствии появились модификации данного алгоритма Томаши-Канаде (Tomasi-Kanade) и Ши- Томаши-Канаде (Shi-Tomasi-Kanade). Трекер Ши-Томаши-Канаде впервые учитывает аффинные искажения окрестных точек. Также существует модификация метода для случая переменного освещения [27].

    В литературе можно обнаружить различные модификации перечисленных методов. В [28] предлагается сравнение производительности различных алгоритмов сопровождения с использованием оптического потока.

    2.5. Применение метода "сдвиг среднего" к решению задачи сопровождения

    Несколько в стороне стоят методы сопровождения, основанные на использовании техники "сдвиг среднего" (mean shift) [8, 23, 24]. Идея методов состоит в том, что для каждой особой точки (в общем случае, для каждого объекта) выбирается окно поиска, вычисляется центр масс распределения интенсивностей (гистограммы). Соответственно центр окна смещается в центр масс, который представляет собой положение точки на текущем кадре. Определение положения точки на последующих кадрах сводится к применению очередного шага метода "сдвига среднего". Метод останавливается, когда центр масс перестает смещаться.

    2.6. Применение фильтра Кальмана для оценки положения объекта

    Задачу сопровождения можно рассматривать как хорошо изученную проблему теории управления, которая состоит в том, чтобы оценить состояние системы на основании последовательности зашумленных измерений [9]. Формально имеется модель объекта, которая наблюдается в зашумленном пространстве сцены. Обозначим через x – модель объекта, z – наблюдение объекта. В общем случае x,z – вектора признаков, которые вполне могут иметь разную размерность. В процессе сопровождения объекта модели и наблюдения могут использоваться двумя способами:

  • Последовательность наблюдений $$z_1,z_2,...$$ применяется для уточнения базовой модели x . Заметим, что модель с течением времени также может изменяться, поэтому каждое новое наблюдение $$z_k$$ может давать новую оценку модели $$x_k$$.
  • Построенная оценка модели $$x_k$$ используется для предсказания модели $$x_{k+1}$$ и наблюдения $$z_{k+1}$$.
  • В результате получаем механизм обратной связи: наблюдаем $$z_k$$, оцениваем $$x_k$$, предсказываем $$x_{k+1}$$ и $$z_{k+1}$$ , обновляем $$x_{k+1}$$ на основании $$z_{k+1}$$. Такой механизм лежит в основе методов сопровождения, основанных на фильтрах Кальмана и фильтрах частиц. В данном разделе остановимся на применении фильтра Кальмана.

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

    $$x_{k+1}=F_kx_k+w_k,$$ $$z_k=H_kx_k+v_k,$$

    где $$F_k$$– матрица преобразования состояния системы; $$w_k$$ – "белый" шум с нормальным распределением $$N(0,Q_k)$$ с математическим ожиданием, равным 0, и матрицей ковариации $$Q_k;H_k$$– матрица связи модели и наблюдения; $$v_k$$ – "белый" шум с нормальным распределением $$N(0,R_k)$$. Заметим, что по определению матрицы ковариации $$Q_k=M(w_kw^T_k),$$, $$(Q_k)_{ij}=M(w^i_kw^j_k),$$ . Если X,Y – две случайные величины, определенные на одном и том же вероятностном пространстве, тогда ковариация численно равна математическому ожиданию от произведения случайных величин X-MX и Y-MY , где MX,MY – математическое ожидание случайных величин X,Y , т.о. cov(x,Y)=M((x-MX)(Y-MY)).

    Приведем пример модели для случая равномерного движения [8]. Состояние системы определяется положением и скоростью точки, тогда состояние может быть представлено вектором вида (4.42).

    $$x_k=\begin{bmatrix} x \\ y \\ v_x \\ v_y \end{bmatrix}_k$$

    Матрица преобразования состояния $$F_k=F$$ и, исходя из физических соображений, может быть записана согласно (4.43).

    $$F=\begin{bmatrix} 1 0 dt 0\\ 0 1 0 dt\\ 0 0 1 0 \\ 0 0 0 1 \end{bmatrix}$$

    Поскольку при сопровождении камера регистрирует только координаты положения объекта, то $$z_k=\begin{bmatrix} z_x \\ z_y \end{bmatrix}_k$$ Поэтому матрицу связи наблюдения и предсказания $$H_k=H$$ можно записать в соответствии с (4.44).

    $$H=\begin{bmatrix} 1 0 \\ 0 1 \\ 0 0 \\ 0 0 \end{bmatrix}$$

    Т.к. в реальных условиях нельзя наблюдать равномерное движение, то вводится Гауссов шум с матрицей ковариации $$Q_k$$ . Матрица ковариации $$R_k$$ строится, исходя из того, насколько точно выполняются измерения положения объекта.

    Теперь необходимо получить обобщенные уравнения для обновления состояния и модели. Идея состоит в том, что сначала строится априорная оценка $$\tilde{x}^-_k$$ состояния согласно (4.40):

    $$\tilde{x}^-_k=F_{k-1}x_{k-1}+w_{k-1}.$$

    Затем с использованием построенной оценки вычисляется наблюдение

    $$z_k=H_k\tilde{x}^-_k+v_k.$$

    Также можно определить апостериорную оценку $$\tilde{x}^+$$ состояния, которая строится после наблюдения. Введем ошибки $$e^-_k$$(4.47) и $$e^+_k$$(4.48), связанные с каждой оценкой состояния.

    $$e^-_k=x_k-\tilde{x}^-_k$$ $$e^+_k=x_k-\tilde{x}^+_k$$

    Фильтр Кальмана оперирует разностью $$z_k-H_k\tilde{x}^-_k$$ , в которую вносит вклад ошибка $$e^-_k$$ и случайный шум $$v_k$$ . В идеальном случае шум отсутствует и оценка состояния идеальна, поэтому указанная разность обращается в ноль. Задача состоит в том, чтобы построить матричный коэффициент Кальмана $$K_k$$ с целью обновления апостериорной оценки согласно (4.49). Отсюда, если $$K_k$$ известен, то известен закон получения обновленной модели $$x_k$$ .

    $$\tilde{x}^+_k=\tilde{x}^-_k+K_k(z_k-H_k\tilde{x}^-_k)$$

    Уравнения (4.47) – (4.49) позволяют получить выражение связи априорной и апостериорной ошибки (4.50).

    $$e^+_k=x_k-\tilde{x}^+_k=x_k-((E-K_kH_k)\tilde{x}^-_k-K_kz_k)=x_k-(E-K_kH_k)\tilde{x}^-_k+K_k(H_kx_k+v_k)=$$ $$=(E-K_kH_k)e^-_k+K_kv_k$$

    Из определения матрицы ковариации пары случайных величин следуют выражения (4.51).

    $$P^-_k=M(e^-_ke^{-T}_k),\,\,\,\,P^+_k=M(e^+_ke^{+^T}_k),\,\,\,\,R_k=M(v_k,v^T_k)$$

    Вследствие независимости ошибок $$e^-_k$$ и $$e^+_k$$ получаем нулевую ковариацию для величин $$e^-_k$$ и $$v_k$$ (4.52).

    $$M(e^-_kv^T_k)=M(v_ke^{-T}_k)=0$$

    Проделав несложные преобразования (4.53) (умножение (4.50) на $$e^{+^T}_k$$ и определение математического ожидания от полученного выражения), можно перейти к уравнению с матрицами ковариации (4.54).

    $$e^+_ke^{+^T}_k=(E-K_kH_k)e^-_ke^{+^T}+K_kv_ke^{+^T}$$ $$=(E-K_kH_k)e^-_k((E-K_kH_k)e^-_k+K_kv_k)^T$$ $$+K_kv_k((E-K_kH_k)e^-_k+K_kv_k)^T$$ $$=(E-K_kH_k)e^-_ke^{-T}_k(E-K_kH_k)^T+(E-K_kH_k)e^-_kv^T_kK^T_k$$ $$+K_kv_ke^{-T}_k(E-K_kH_k)^T+K_kv_kv^T_kK^T_k$$ $$P^+_k=(E-K_kH_k)P^-_k(E-K_kH_k)^T+K_kR_kK^T_k$$

    Отсюда получаем, что матрицу $$K_k$$ необходимо выбрать так, чтобы сумма диагональных элементов (след) матрицы $$P^+_k$$ был минимален.

    $$trace(p^+_k) \rightarrow min_{K_k}$$

    Решение данной задачи определяется посредством дифференцирования функции $$trace(p^+_k)$$ по $$K_k$$ (4.56) Результат дифференцирования получается благодаря применению свойства $$\frac {\delta} {\delta A} (trace(ABA^T))=2AB$$ при условии, что матрица B является симметричной. .

    $$-2(E-K_kH_k)P^-_kH^T_k+2K_kR_k=0$$

    Как следствие, получаем формулу (4.57) для вычисления матрицы $$K_k$$.

    $$K_k=P^-_kH^T_k(R_k+H_kP^-_kH^T_k)^{-1},$$

    где

    $$P^-_k=F_kP^+_{k-1}F^T_k+Q_{k-1}.$$

    Заметим, что посредством несложных математических выкладок можно вывести формулу (4.59) для вычисления

    $$P^+_k=(E-K_kH_k)P^-_k$$

    Подводя итог, более четко выделим схему работы фильтра:

  • Предсказание. Предполагает вычисление априорной оценки состояния $$\tilde{x}^-_k$$ (4.45) и наблюдения $$z_k$$ (4.46).
  • Коррекция. Включает определение матричного коэффициента Кальмана $$K_k$$ (4.57) и построение апостериорной оценки состояния $$\tilde{x}^+_k$$ согласно (4.49).
  • 2.7. Применение фильтра частиц к задаче сопровождения

    Для некоторых практических задач, чтобы получить более точную оценку состояния системы, необходимо уйти от предположения, что шум имеет Гауссово распределение [5]. В этом случае вводится понятие мультимодального распределения Мультимодальным распределением называется распределение, имеющее несколько мод или локальных максимумов. Мультимодальное распределение зачастую представляется смесью нескольких распределений. шума, а для моделирования подобных систем используются фильтры частиц. Фильтры частиц являются более общим подходом к решению задачи сопровождения с применением вероятностных методов [9, 21, 33, 35].

    Алгоритм воспроизведения условной плотности (CONditional DENSity propAGATION, CONDENSATION) [21, 35] – базовый алгоритм фильтрации частиц, на основании которого строится большинство алгоритмов данной группы, применяемых в компьютерном зрении. Поэтому остановимся более детально на рассмотрении схемы работы именно этого алгоритма.

    Предполагая, что система может находиться в состояниях $$X_t={x_1,x_2,...,x_t}$$, в момент времени получим плотность распределения вероятности. Так же, как и для фильтра Кальмана, последовательность наблюдений будем обозначать $$Z_t={z_1,z_2,...,z_t}$$. Наряду с этим введем предположение о том, что состояние $$x_{t}$$ зависит только от предыдущего состояния $$x_{t-1}$$ – условие Марковской цепи. Таким образом, получаем систему с независимым набором наблюдений. Техники фильтрации частиц представляют распределение вероятности в виде коллекции взвешенных выборок – частиц, появление которых регулируется посредством введения весов. Тогда множество $$S_t$$ (4.60) определяет функцию плотности вероятности для состояния $$x_t$$ при заданном наборе наблюдений $$Z_t$$ . $$S_t$$ задает приближенное распределение $$p(x_t \lvert Z_t)$$.

    $$S_t=\lbrace (s^t_i,\pi^t_i),i=\overline{1,N}, \sum_{i=1}^N {\pi^t_i=1} \rbrace$$

    Задача состоит в том, чтобы построить метод восстановления множества $$S_{t}$$ на основании $$S_{t-1}$$. Формально алгоритм можно представить в виде последовательности этапов [9, 21, 35]:

  • Пусть коллекция взвешенных выборок в момент времени $$t-1$$ построена (4.61). $$S_{t-1}=\lbrace (s^{t-1}_i,\pi^{t-1}_i),i=\overline{1,N}, \sum_{i=1}^N {\pi^{t-1}_i=1} \rbrace$$ Дополнительно вычислим интегральные веса согласно (4.62). $$c_i=c_{i-1}+\pi^{t-1}_i,i=\overline{1,N},c_0=0$$
  • Определим n-ый экземпляр выборки $$S_t$$ . Для этого случайным образом выберем число r из отрезка [0,1] и вычислим $$j=arg\,min_i{c_i>r}$$. Отсюда получаем текущую оценку состояния $$s^{t-1}_j$$ .
  • Выполним предсказание следующего состояния. Предсказание (4.63) выполняется аналогично фильтру Кальмана (4.40), разница лишь в том, что нет ограничений, связанных с линейностью системы и видом распределения шума. $$s^{t}_n=F_{t-1}s^{t-1}_j+w_{t-1}$$
  • Выполним коррекцию. Используя текущее наблюдение и его распределение, необходимо установить вес полученного экземпляра согласно (4.64). $$\pi^t_n=p(z_t \lvert x_t=s^t_n)$$
  • Построим множество частиц , повторив шаги 2 – 4 -раз.
  • Нормализуем последовательность весов $$\pi^t_i$$ так, чтобы $$\sum^N_{i=1} {\pi^t_i}=1$$ .
  • Вычислим наилучшую оценку для состояния $$x_t$$, например, как линейную свертку полученного набора экземпляров выборки (4.65). Таким образом, фактически определим некоторую среднюю частицу. $$x_t \sum^N_{i=1} {\pi^t_is^t_i}$$
  • Описанный процесс можно проинтерпретировать графически (рис.4.8) с использованием понятия частицы. На входе итерации алгоритма имеется множество частиц $$\lbrace (s^{t-1}_i,\pi ^{t-1}_i) \rbrace$$ (верхний уровень диаграммы). В результате N-кратного случайного выбора частиц из $$S_{t-1}$$ получается некоторый набор экземпляров (2-ой уровень сверху). Применение шага предсказания приводит к формированию множества оценочных состояний частиц (3-ий уровень), затем для каждой оценки выполняется коррекция на основании имеющихся наблюдений. Как следствие, создается множество частиц $$\lbrace (s^{t}_i,\pi ^{t}_i) \rbrace$$ в следующий момент времени (последний уровень диаграммы).

    (рис 4.8) Итерация алгоритма воспроизведения условной плотности (CONDENSATION, [35])

    3. Контрольные вопросы

  • В каких случаях схема вычитания фона в задаче определения движения работает неэффективно?
  • Докажите, что формула (4.25) справедлива, т.е. образ Фурье от корреляции пары функций равен произведению образу Фурье первой функции на комплексно-сопряженный образ второй функции.
  • Основное свойство аффинного преобразования?
  • Приведите примеры задач, в которых используется модель перспективной проекции (гомография).
  • Докажите справедливость формул (4.58) и (4.59).
  • Получите формулы вычисления матричного коэффициента Кальмана для случая равномерного движения.
  • Приведите схему построения фильтра Кальмана в случае одномерного состояния и наблюдения при некоторых фиксированных значениях $$F_k=F=const$$ и $$H_k=H=const$$.
  • Вернуться к учебному плану