САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
КАФЕДРА КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ И
МНОГОПРОЦЕССОРНЫХ СИСТЕМ
МАГИСТЕРСКАЯ ДИССЕРТАЦИЯ
Тема: «Оценка и визуализация оптического
потока в задачах обработки изображений»
Направление:
02.04.02 – ФИиИТ
Магистерская программа: ВМ.5502.2014 – Вычислительные технологии
Выполнил студент гр. 633
Смирнов Константин Валерьевич
Научный руководитель,
к. т. н., доцент
В. М. Гришкин
Санкт-Петербург – 2016
Содержание
Введение
Глава 1
4
Обзор существующих методов вычисления
оптического потока
6
§1
Алгоритм Лукаса-Канаде . . . . . . . . . . . . . . . . . . . .
6
§2
Алгоритм Horn-Schunck . . . . . . . . . . . . . . . . . . . . .
14
§3
Обобщения для предположений постоянства производных .
19
§4
Детали дискретизации . . . . . . . . . . . . . . . . . . . . . .
22
Глава 2
Алгоритм построения поля скоростей
24
§1
Смешаный алгоритм Лукаса-Канаде и Horn Schunck . . . . .
24
§2
Модификация для выделения границ . . . . . . . . . . . . .
28
§3
Система линейных уравнений . . . . . . . . . . . . . . . . . .
30
§4
Численные методы решения систем линейных уравнений . .
32
§5
Сходимость блочных методов . . . . . . . . . . . . . . . . . .
33
Глава 3
Программная реализация вычисления оптического
потока
35
§1
Общая схема реализации вычислительной системы . . . . .
35
§2
Средства реализации . . . . . . . . . . . . . . . . . . . . . .
37
§3
Описание жизненного цикла задачи . . . . . . . . . . . . . .
39
Глава 4
Примеры и анализ работы алгоритмов
42
§1
Пример: движение прямоугольника . . . . . . . . . . . . . .
42
§2
Пример: сдвиг изображения . . . . . . . . . . . . . . . . . . .
47
§3
Пример: сдвиг карты . . . . . . . . . . . . . . . . . . . . . .
50
§4
Анализ работы алгоритмов . . . . . . . . . . . . . . . . . . .
53
Заключение
58
2
Введение
В настоящее время широкий интерес вызывают задачи распознавания [1-9]. В данной работе рассматривается задача распознавания движения на изображении или на видео. Современные информационные технологии позволяют ускорить методы, которые применяются для решения данной задачи, например, за счет параллельной обработки различных частей
видео или изображений из одной временной последовательности.
Целью данной работы является построение прототипа распределенной системы для распознавания движения на видео или последовательности изображений.
В качестве основы можно рассмотреть, семейство методов, основанных на представлении движения на изображении в форме поля скоростей
или оптического потока [4-5, 10]. Поле скоростей (оптический поток) - одна из форм представления движения объектов на изображении, при котором каждой точке изображения ставится в соответствие вектор скорости,
отвечающий соответствующей точке на объекте. При постановке задачи
распознавания движения в форме поиска оптического потока существует
несколько методов для ее решения.
1. Фазовая корреляция — основан на инверсия нормализованного перекрестного спектра.
2. Блочные методы — минимизация суммы квадратов или суммы модулей разностей.
3. Группа дифференциальных методов оценки оптического потока,
основанных на частных производных:
3.1. Алгоритм Лукаса-Канаде — рассматриваются части изображения
и аффинная модель движения.
3.2. Horn–Schunck — минимизация функционала, выражающего предположения постоянства цвета объектов на изображении и гладкости полу4
чаемого векторного поля.
3.3. Buxton–Buxton — основан на модели движения границ объектов
в последовательности изображений.
4. Общие вариационные методы — модификации метода Horn-Schunck,
использующие другие ограничения на данные и другие ограничения на
гладкость.
5. Дискретные методы оптимизации — поисковое пространство квантуется, затем каждому пикселю изображения ставится в соответствие метка таким образом, чтобы расстояние между последовательными кадрами
было минимальным. Оптимальное решение часто ищется с помощью алгоритмов нахождения минимального разреза и максимального потока в
графе, линейного программирования или belief propagation.
В этой работе будут рассматриваться дифференциальные методы оценки поля скоростей, в его основе лежит вычисление частных производных по
горизонтальному и вертикальному направлениям изображения. Как будет
показано далее, одних только производных недостаточно чтобы определить
смещения точек объектов на изображении. Поэтому существует множество
модификаций названых алгоритмов, заточенных под определенные частные задачи.
Для построения прототипа распределенной системы для распознавания движения на видео или последовательности изображений в работе
решаются следующие задачи:
1.Анализ существующих методов оценки поля скоростей, выбор основного алгоритма
2.Создание распределенной компьютерной реализации алгоритма вычисления поля скоростей на основе облачных технологий.
3.Анализ результатов.
5
Глава 1 Обзор существующих методов вычисления
оптического потока
§1 Алгоритм Лукаса-Канаде
В данной главе рассматриваются несколько дифференциальных методов оценки поля скоростей. Начнем с рассмотрения самого простого из
них алгоритма Лукаса-Канаде.
Будем считать, что рассматриваются черно-белые изображения, в
этом случае цвет точки изображения определяется одним числом, которое
будем далее называть яркостью.
В основе алгоритма Лукаса-Канаде лежит одно очень важное предположение: будем считать, что значения яркостей пикселей переходят из
одного кадра в следующий без изменений. Таким образом, мы делаем допущение, что пиксели, относящиеся к одному и тому же объекту, могут
сместиться в какую либо сторону, но их значение останется неизменным.
Конечно, это предположение слабо связано с реальностью, так как от кадра к кадру могут меняться глобальные условия освещения и освещенность
самого движущегося объекта. Это допущение в дальнейшем приводит к
некоторым проблемам, но, не смотря на это, оно достаточно хорошо работает на практике.
Рис. 1: Пример линейных изображений.
На математическом языке это допущение можно записать так: I(x, y, t) =
I(x + ux , y + uy , t + 1). Где I — это функция яркости пикселей от положения
6
Рис. 2: Линейные изображения как функции.
на кадре и времени. Другими словами x и y — это вещественные координаты точки в плоскости кадра, ux и uy — это смещение, а t — это номер
кадра в последовательности. Условимся, что между двумя соседними кадрами проходит единичный отрезок времени.
Для начала рассмотрим одномерный случай. Представим себе два
одномерных кадра 1 пиксель в высоту и 20 пикселей в ширину. На втором
кадре изображение немного смещено вправо. Именно это смещение мы и
хотим найти. Для этого представим эти же кадры в виде функций. На входе
позиция пикселя, на выходе — его интенсивность. В таком представление
искомое смещение d видно еще более наглядно. В соответствии с нашим
предположением, f2 это просто смещенная f1 , то есть можем сказать, что
f2 = I(x − d).
Обратим внимание, что f1 и f2 при желании можно записать и в
общем виде: f1 (x) = I(x, y, t), f2 (x) = I(x, y, t) где y и t зафиксированы и
равны нулю.
Для каждой координаты нам известны значения f1 и f2 в этой точке,
кроме того мы можем вычислить их производные. Свяжем известные значения со смещением d. Для этого запишем разложение в ряд Тейлора для
f1 (x − d):
7
f1 (x − d) = f1 (x) − df10 (x) + O(d2 f100 )
Сделаем второе важное предположение: Предположим, что f1 (x −
d) достаточно хорошо аппроксимируется первой производной. Сделав это
предположение, отбросим всё что после первой производной:
f1 (x − d) = f1 (x) − df10 (x)
Насколько это корректно? В общем-то, не очень, тут мы теряем в
точности, если только наша функция/изображение не строго линейна, как
в нашем искусственном примере. Зато это существенно упрощает метод,
а для достижения требуемой точности можно сделать последовательное
приближение, которое мы рассмотрим позже.
Мы почти у цели. Смещение d — это наша искомая величина, поэтому
надо что-то сделать с f1 (x−d). Как мы условились ранее, f2 (x) = f1 (x−d),
поэтому просто перепишем:
f2 (x) = f1 (x) − df10 (x)
То есть:
d=
f1 (x) − f2 (x)
f10 (x)
Теперь перейдем от одномерного случая к двумерному. Запишем разложение в ряд Тейлора для I(x+ux , y +uy , t) и сразу отбросим все старшие
производные. Вместо первой производной появляется градиент:
I(x + ux , y + uy , t) = I(x, y, t) + u∇I(x, y, t)
Где u = (ux , uy )T — вектор смещения. В соответствии со сделанным
допущением I(x, y, t) = I(x + ux , y + uy , t + 1). Обратите внимание, что это
8
выражение эквивалентно I(x + ux , y + uy , t) = I(x, y, t + 1). Это то, что нам
нужно. Перепишем:
I(x, y, t + 1) = I(x, y, t) + u∇I(x, y, t)
I(x, y, t) − I(x, y, t + 1) + u∇I(x, y, t) = 0
Поскольку между двумя кадрами проходит единичный интервал времени, то можно сказать, что I(x, y, t) − I(x, y, t + 1) есть не что иное, как
производная по времени. Перепишем:
∂I(x, y, t)
+ u∇I(x, y, t) = 0
∂t
Перепишем ещё раз, раскрыв градиент:
∂I(x, y, t)
∂I(x, y, t)
∂I(x, y, t)
+ ux
+ uy
=0
∂t
∂x
∂y
Мы получили уравнение, которое говорит нам о том, что сумма частных производных должны быть равна нулю. Проблема только в том, что
уравнение у нас одно, а неизвестных в нем два: ux и uy . На этом моменте
начинается полет фантазии и разнообразие подходов.
Сделаем третье предположение: Предположим, что соседние пиксели смещаются на одинаковое расстояние. Возьмем фрагмент изображения,
скажем 5 на 5 пикселей, и условимся, что для каждого из 25 пикселей ux и
uy равны. Тогда вместо одного уравнения мы получим сразу 25 уравнений!
Очевидно, что в общем случае система не имеет решения, поэтому будем
искать такие ux и uy , которые минимизируют ошибку:
E(ux , uy ) =
X
i,j
∂I(xi , yj , t)
∂I(xi , yj , t)
∂I(xi , yj , t)
g(xi , yj )
+ ux
+ uy
∂t
∂x
∂y
9
2
Здесь g — это функция, определяющая весовые коэффициенты для
пикселей. Самые распространенный вариант — двухмерная гауссиана, которая дает наибольший вес центральному пикселю и все меньший по мере
удаления от центра.
Чтобы найти минимум E(ux , uy ) воспользуемся методом наименьших
квадратов, найдем её частные производные по ux и uy :
∂E(ux , uy )
=
∂ux
X
∂I(xi , yj , t)
∂I(xi , yj , t)
∂I(xi , yj , t) ∂I(xi , yj , t)
=2
g(xi , yj )
+ ux
+ uy
=
∂t
∂x
∂y
∂x
i,j
2
X
∂I(xi , yj , t)
∂I(xi , yj , t) ∂I(xi , yj , t)
+
=2
g(xi , yj )[ux
+ uy
∂x
∂y
∂x
i,j
+
∂I(xi , yj , t) ∂I(xi , yj , t)
]
∂t
∂x
∂E(ux , uy )
=
∂uy
∂I(xi , yj , t)
∂I(xi , yj , t) ∂I(xi , yj , t)
∂I(xi , yj , t)
=2
g(xi , yj )
+ ux
+ uy
=
∂t
∂x
∂y
∂y
i,j
2
X
∂I(xi , yj , t) ∂I(xi , yj , t)
∂I(xi , yj , t)
=2
g(xi , yj )[ux
+ uy
+
∂x
∂y
∂y
i,j
X
+
∂I(xi , yj , t) ∂I(xi , yj , t)
]
∂t
∂y
Перепишем в более компактной форме и приравняем к нулю:
"
∂E(ux , uy )
=
∂ux
X
∂E(ux , uy )
=
∂uy
X
g(xi , yj ) ux
i,j
"
i,j
g(xi , yj ) ux
∂I
∂x
2
+ uy
∂I ∂I
+ uy
∂x ∂y
10
∂I ∂I ∂I ∂I
+
∂y ∂x ∂t ∂x
∂I
∂y
2
+
∂I ∂I
∂t ∂y
#
#
Перепишем эти два уравнения в матричной форме:
Mu = b
Где
P
P
∂I 2
∂I ∂I
i,j g(xi , yj ) ∂x
i,j g(xi , yj ) ∂x ∂y
2
M =P
P
∂I ∂I
∂I
i,j g(xi , yj ) ∂x ∂y
i,j g(xi , yj ) ∂y
P
∂I ∂I
g(xi , yj ) ∂t ∂x
b = − Pi,j
∂I ∂I
i,j g(xi , yj ) ∂t ∂y
ux
u = −
uy
Если матрица обратима (имеет ранг 2), можем вычислить ux и uy ,
которые минимизируют ошибку E:
u
e = M −1 b
Мы знаем приблизительное смещение пикселей между двумя соседними кадрами.
Поскольку в нахождении смещения каждого пикселя участвуют также соседние с ним пиксели, при реализации данного метода целесообразно
предварительно посчитать производные кадра по горизонтали и вертикали.
Описанный выше метод основан на трех значительных допущениях,
которые с одной стороны дают нам принципиальную возможность определить оптический поток, но с другой стороны вносят погрешность. Хорошая
новость для перфекционистов состоит в том, что одно допущение нужно
нам только для упрощения метода, и с его последствиями мы можем бороться. Мы предполагали, что для аппроксимации смещения нам будет
достаточно первой производной. В общем случае это конечно же не так
(рисунок слева). Для достижение требуемой точности смещение для каждой пары кадров (назовём их Fi и Fi+1 ) можно вычислять итеративно. В
11
литературе это называется искажением (warping). На практике это означает, что, вычислив смещения на первой итерации, мы перемещаем каждый
пиксель кадра Fi+1 в противоположную сторону так, чтобы это смещение
компенсировать. На следующей итерации вместо исходного кадра Fi+1 мы
0
будем использовать его искаженный вариант Fi+1
. И так далее, пока на
очередной итерации все полученные смещения не окажутся меньше заданного порогового значения. Итоговое смещение для каждого конкретного
пикселя мы получаем как сумму его смещений на всех итерациях.
Рис. 3: Пример линейных изображений.
По своей природе данный метод является локальным, то есть при
определении смещения конкретного пикселя принимается во внимание только область вокруг этого пикселя — локальная окрестность. Как следствие,
невозможно определить смещения внутри достаточно больших (больше
размера локальной окрестности) равномерно окрашенных участков кадра. К счастью на реальных кадрах такие участки встречаются не часто,
но эта особенность все же вносит дополнительное отклонение от истинного
смещения.
Ещё одна проблема связана с тем, что некоторые текстуры в изображении дают вырожденную матрицу , для которой не может быть найдена обратная матрица. Соответственно, для таких текстур мы не сможем
12
определить смещение. То есть движение вроде есть, но непонятно в какую
сторону. В общем-то, от этой проблемы страдает не только рассмотренный
метод. Даже глаз человека воспринимает такое движение не однозначно.
13
§2 Алгоритм Horn-Schunck
Следующий дифференциальный метод оценки поля скоростей, который будет нами рассмотрен — алгоритм Horn-Schunck
Будем использовать введенные ранее обозначения области в пространстве плоскости изображений и времени Ω ⊂ R3 и подобласти, отвечающей
моменту времени t0 Ωt0 . Плоскость изображений - модель экрана, на который проецируется движение. Экран по форме с течением времени не
меняется, следовательно Ωt1 = Ωt2 , ∀t1 , t2 ∈ R. Проецирование будем полагать ортогональным. Если наблюдатель перемещается в пространстве, то
эта плоскость будет перемещаться вместе с ним, следовательно, при отсутствии движения объектов в пространстве на плоскости изображений будет
просматриваться движение объектов, обратное движению наблюдателя.
Внешние условия, отвечающие окружающей среде, безусловно, будут
влиять на изображение движения. Мы будем определять проекцию движения, на плоскость изображений, так как у нас нет данных о движении,
ортогональном этой плоскости.
Основной характеристикой, которой мы будем пользоваться, будет
являться яркость точек изображений. На нее кроме текстуры передвигающихся объектов будет влиять освещенность.
В самом простом случае, можно предположить постоянство освещения, чему будет отвечать постоянство яркости, это можно выразить следующим образом:
u
I(x + u) − I(x) → 0, x, u ∈ Ωt0 .
Если освещение изменяется равномерно по всему пространству, т.е.
по всей плоскости изображений, то мы можем перейти к градиенту яркости
чтобы рассмотреть рельефность изображения и записать это в виде
14
u
5I(x + u) − 5I(x) → 0, x, u ∈ Ωt0 .
Если учитывать возможность локального изменения освещенности,
то изображение рельефности изображения может искажаться, при точечной модели освещения и условии равномерного распространения цвета, это
искажение будет только в форме круга, чем дальше от центра которого тем
будет меньше искажение. Переход к матрицам Гессе для яркости позволяет
учитывать такую модель, это можно записать в виде:
u
Ixi xj (x + u) − Ixi xj (x) → 0, i, j = 1, 2, x, u ∈ Ωt0 .
Для большей эффективности, эти условия можно комбинировать вместе
Для ввода зависимости между соседними точками плоскости изображений введем в качестве ограничения гладкость потока по переменным
x1 , x2 , для этого будем минимизировать разницу скоростей соседних точек на плоскости изображений. С помощью градиента потока это можно
выразить в следующем виде:
u
k 5 u1 k2 + k 5 u2 k2 → 0, u ∈ Ωt0 .
Здесь и далее в этой главе будет рассматриваться Евклидова норма
матрицы
v
u (n,m)
uX
ai,j , где матрицаA = {ai,j }i=1...n,j=1...m .
kAk = t
i,j=1,1
Далее рассмотрим формализацию для пар условий постоянства яркости и гладкости потока, условий постоянства Гессиана яркости и гладкости
потока.
Запишем условия и в интегральной форме для всех точек плоскости
изображений в момент времени t0 :
15
Z
Ωt0
([I(x + u) − I(x)]2 )dx1 dx2 + α
Z
u
(k 5 u1 k2 + k 5 u2 k2 )dx1 dx2 → 0.
Ωt 0
Здесь α > 0 - параметр регуляризации. В общем случае для произвольной последовательности изображений нам ничего о них неизвестно,
поэтому будем полагать граничные условия нулевыми.
Переходя к уравнениям Эйлера-Лагранжа, мы получим
I (x + u)[I(x + u) − I(x)] + α 52 u = 0
x1
1
.
Ix (x + u)[I(x + u) − I(x)] + α 52 u2 = 0
2
Далее, производя аппроксимацию по формулам Тейлора, придем к
следующей системе
I (x)[I (x) + I (x)u + I (x)u ] − α 52 u = 0
x1
t
x1
1
x2
2
1
.
Ix (x)[It (x) + Ix (x)u1 + Ix (x)u2 ] − α 52 u2 = 0
2
1
2
К дискретной форме перейдем следующим образом. Положим u1ij и
u2ij как проекции скоростей по осям в точке с координатами i, j, а x1ij и
x2ij как координаты соответствующей точки xij . i ∈ 1..N − 1, j ∈ 1..M − 1.
Для границ i = 0, i = N, j = 0, j = M положим u1ij и u2ij равным граничным условиям, задаваемым изначально, обычно они берутся равными
нулю. Раскроем лапласиан в точке с координатами i, j по следующей формуле:
52 u1ij = u1i+1j + u1i−1j + u1ij−1 + u1ij−1 − 4u1ij
52 u2ij = u2i+1j + u2i−1j + u2ij−1 + u2ij−1 − 4u2ij
Тогда система уравнений принимает вид линейной системы уравне-
16
ний:
Ix1 (xij )[It (xij ) + Ix1 (xij )u1ij + Ix2 (xij )u2ij ]−
− 1 α(u
+u
+u
+u
− 4u ) = 0
4
1i+1j
1i−1j
1ij−1
1ij−1
1ij
.
Ix2 (xij )[It (xij ) + Ix1 (xij )u1ij + Ix2 (xij )u2ij ]−
− 1 α(u2i+1j + u2i−1j + u2ij−1 + u2ij−1 − 4u2ij ) = 0
4
Введем обозначения:
A11ij = Ix1 (xij )Ix1 (xij ), A12ij = Ix1 (xij )Ix2 (xij ), A22ij = Ix2 (xij )Ix2 (xij ),
B1ij = Ix1 (xij )It (xij ), B2ij = Ix1 (xij )It (xij ).
Сгруппируем u1ij и u2ij :
(A
11ij
+ α)u1ij + A12ij u1ij = α(u1i+1j + u1i−1j + u1ij−1 + u1ij−1 ) − B1ij
.
A12ij u2ij + (A22ij + α)u2ij = α(u2i+1j + u2i−1j + u2ij−1 + u2ij−1 ) − B2ij
Слева получилась невырожденная матрица, так как определитель:
(A11ij + α)(A22ij + α) − A12ij A12ij =
= Ix1 (xij )Ix1 (xij )Ix2 (xij )Ix2 (xij ) + (Ix1 (xij )Ix1 (xij ) + Ix2 (xij )Ix2 (xij )) +
+α2 − Ix1 (xij )Ix1 (xij )Ix2 (xij )Ix2 (xij ) =
= Ix1 (xij )2 + Ix2 (xij )2 + α2 > 0
Будем решать систему итерационным методом:
(A
11ij
k−1
k−1
k−1
+ α)uk1ij + A12ij uk1ij = α(uk−1
1i+1j + u1i−1j + u1ij−1 + u1ij−1 ) − B1ij
A12ij uk + (A22ij + α)uk = α(uk−1 + uk−1 + uk−1 + uk−1 ) − B2ij
2ij
2ij
2i+1j
2i−1j
2ij−1
2ij−1
.
Пусть:
k−1
k−1
k−1
u
ek−1
= u1i+1j
+ uk−1
1
1i−1j + u1ij−1 + u1ij−1
k−1
k−1
k−1
u
ek−1
= uk−1
2
2i+1j + u2i−1j + u2ij−1 + u2ij−1
Тогда окончательное решение системы линейных уравнений запишется в
виде:
uk = u
ek−1
− (A11ij + A22ij + α2 )−1 (A11ij u
e1k−1 + A12ij u
ek−1
+ B1ij )
1ij
1
2
.
k−1
k−1
k−1
2 −1
uk = u
e2 − (A11ij + A22ij + α ) (A12ij u
e1 + A22ij u
e2 + B2ij )
2ij
17
При k = 0 u1ij и u2ij полагаются равным начальному приближению,
рационально брать это приближение равным нулю или равным результату,
полученным для предыдущего момента t.
18
§3 Обобщения для предположений постоянства производных
Приведенный в предыдущем параграфе алгоритм можно обобщить
на случаи рассмотрения производных от функции I(x, y, t). На практике
это имеет смысл, если цвет объекта меняется с течением времени и предположение постоянства яркости теряет смысл. Например, это возможно
при изменении освещенности сцены с течением времени. При одинаковом
изменении освещенности с течением времени имеет смысл рассматривать
вместо постоянства яркости I(x, y, t) = I(x + ux , y + uy , t + 1) постоянство
ее производных:
∂I(x, y, t) ∂I(x + ux , y + uy , t + 1)
=
∂x
∂x
∂I(x, y, t) ∂I(x + ux , y + uy , t + 1)
=
∂y
∂y
Если освещенность изменяется в плоскости изображения с течением
времени не равномерно, т.е. при наличии динамических источников света,
то имеет смысл рассматривать предположения постоянства производных
более высокого порядка: постоянство гессиана, лапласиана, определителя
матрицы Гессе:
H(I(x, y, t)) = H(I(x + ux , y + uy , t + 1))
∂ 2 I(x, y, t) ∂ 2 I(x, y, t)
+
=
∂x2
∂y 2
=
∂ 2 I(x + ux , y + uy , t + 1) ∂ 2 I(x + ux , y + uy , t + 1)
+
∂x2
∂y 2
det(H(I(x, y, t)) = det(H(I(x + ux , y + uy , t + 1))
где
H(I(x, y, t) =
∂ 2 I(x+ux ,y+uy ,t+1)
∂x2
2
∂ I(x+ux ,y+uy ,t+1)
∂x∂y
19
∂ 2 I(x+ux ,y+uy ,t+1)
∂x∂y
∂ 2 I(x+ux ,y+uy ,t+1)
∂y 2
В случае постоянства гессиана обобщение можно сделать следующим
образом:
X
Ωt0
(2,2)
Z
(Ixn xm (x + u) − Ixn xm (x))2 dx1 dx2 +
n,m=1,1
Z
+α
u
(k 5 u1 k2 + k 5 u2 k2 )dx1 dx2 → 0.
Ωt 0
Здесь α > 0 - параметр регуляризации. В общем случае для произвольной последовательности изображений нам ничего о них неизвестно,
поэтому будем полагать граничные условия нулевыми.
Переходя к уравнениям Эйлера-Лагранжа, мы получим
P(2,2)
n,m=1,1 (Ixn xm x1 (x + u)[(Ixn xm (x + u) − Ixn xm (x)]) +
+α 52 u1 = 0
.
P(2,2)
n,m=1,1 (Ixn xm x2 (x + u)[(Ixn xm (x + u) − Ixn xm (x)]) +
+α 52 u2 = 0
Далее, производя аппроксимацию по формулам Тейлора, придем к
следующей системе
P(2,2)
n,m=1,1 (Ixn xm x1 (x)[Ixn xm t (x) + Ixn xm x1 (x)u1 + Ixn xm x2 (x)u2 ]) +
+α 52 u1 = 0
.
P(2,2)
n,m=1,1 (Ixn xm x2 (x)[Ixn xm t (x) + Ixn xm x1 (x)u1 + Ixn xm x2 (x)u2 ]) +
+α 52 u2 = 0
Эту систему можно решать тем же образом, который рассматривался
для решения системы в предыдущем параграфе, только матрицы M теперь
20
будет записываться следующим образом:
A11ij =
A11ij =
P(2,2)
n,m=1,1 Ixn xm x1 (xij )Ixn xm x1 (xij ),
P(2,2)
n,m=1,1 Ixn xm x1 (xij )Ixn xm x2 (xij ),
A11ij =
P(2,2)
B1ij =
P(2,2)
B2ij =
n,m=1,1 Ixn xm x2 (xij )Ixn xm x2 (xij ),
n,m=1,1 Ixn xm x1 (xij )Ixn xm t (xij ),
P(2,2)
n,m=1,1 Ixn xm x2 (xij )Ixn xm t (xij ).
21
§4 Детали дискретизации
Рассмотрим случай предположения постоянства яркости точки и более подробно опишем, как можно осуществлять переход от системы дифференциальных уравнений к системе алгебраических уравнений.
Есть система дифференциальных уравнений:
I (x)[I (x) + I (x)u + I (x)u ] − α 52 u = 0
x1
t
x1
1
x2
2
1
.
Ix (x)[It (x) + Ix (x)u1 + Ix (x)u2 ] − α 52 u2 = 0
2
1
2
Как и раньше, положим u1ij и u2ij как проекции скоростей по осям
в точке с координатами i, j, а x1ij и x2ij как координаты соответствующей
точки xij . i ∈ 1..N − 1, j ∈ 1..M − 1. Для границ i = 0, i = N, j = 0, j = M
положим u1ij и u2ij равным граничным условиям, задаваемым изначально,
обычно они берутся равными нулю.
Возможно несколько способов раскрытия лапласиана. Рассматриваемый ранее:
52 u1ij ≈ u1i+1j + u1i−1j + u1ij−1 + u1ij−1 − 4u1ij
52 u2ij ≈ u2i+1j + u2i−1j + u2ij−1 + u2ij−1 − 4u2ij
Оценка производных функции I(x, y, t) тоже может быть осуществлена несколькими способами:
∂I( xij ) 1
≈ (I(xi+1j ) − I(xi−1j ))
∂x
2
∂I(xij )
≈ I(xij ) − I(xi−1j )
∂x
∂I(xij )
≈ I(xi+1j ) − I(xij )
∂x
Анализ точности этих формул можно найти в [14]. Возможно использование более сложных формул, которые требуют большего количества
точек, но при оценке производных по времени это тесно связано с количеством доступных изображений.
22
Использование более сложных конструкций приводит к размытию
оригинальной производной.
Например, рассмотрим движение черного прямоугольника на белом
фоне вдоль оси y. Пусть момент времени t + 1 отличается от момента времени t сдвигом прямоугольника на один пиксель изображения. Таким образом, ожидается наличие производной по оси y не равной нулю только
для двух фиксированных значений x, соответствующих верхней и нижней
грани прямоугольника. При использовании первой из трех указанных разностных схем происходит размытие производной вдоль оси y, движению
верхней грани будет соответствовать не одна ненулевая производная, а две
в два раза меньшие оригинальной.
В оригинальной статье для метода Horn-Schunk используется оценка
производных с усреднением:
∂I(xijk ) 1
≈ (I(xi+1jk ) − I(xijk ) + I(xi+1j+1k ) − I(xij+1k )+
∂x
4
+I(xi+1jk+1 ) − I(xijk+1 ) + I(xi+1j+1k+1 ) − I(xij+1k+1 ))
∂I(xijk ) 1
≈ (I(xij+1k ) − I(xijk ) + I(xi+1j+1k ) − I(xi+1jk )+
∂y
4
+I(xij+1k+1 ) − I(xijk+1 ) + I(xi+1j+1k+1 ) − I(xi+1jk+1 ))
∂I(xijk ) 1
≈ (I(xijk+1 ) − I(xijk ) + I(xij+1k+1 ) − I(xij+1k )+
∂t
4
+I(xi+1jk+1 ) − I(xi+1jk ) + I(xi+1j+1k+1 ) − I(xi+1j+1k ))
где k — дискретный момент времени, при разбиении для имеющихся изображений. Такой подход дополнительно размывает изображение.
23
Глава 2 Алгоритм построения поля скоростей
§1 Смешаный алгоритм Лукаса-Канаде и Horn Schunck
В рамках данной работы предлагается модификация алгоритма HornSchunck, которая заключается в объединении оригинального алгоритма
Horn-Schunck и алгоритма Лукаса-Канаде, и модификация выделения границ. Далее в этом параграфе приводится математическая формализация
этого смешанного алгоритма.
Будем использовать введенные ранее в описании алгоритма Horn Schuncd
обозначения.
Ранее использовалось условие постоянства яркости в точке:
u
I(x + u) − I(x) → 0, x, u ∈ Ωt0 .
Теперь заменим его на следующие r условий:
u
I(x[d] + u) − I(x[d] ) → 0, x[d] , u ∈ Ωt0 .
Здесь множество точек {x[d] , d = 1...r} - точки из окрестности точки
x. Запишем это условие вместе с условием гладкости потока для всех точек
плоскости изображений в момент времени t0 :
Z
Z
r h
i2
X
( I(x[d] + u) − I(x[d] ) )dx1 dx2 +α
Ωt0 d=1
u
(k5u1 k2 +k5u2 k2 )dx1 dx2 → 0.
Ωt 0
Здесь α > 0 - параметр регуляризации. В общем случае для произвольной последовательности изображений нам ничего о них неизвестно,
поэтому будем полагать граничные условия нулевыми.
Переходя к уравнениям Эйлера-Лагранжа, мы получим
24
Pr I (x[d] + u)[I(x[d] + u) − I(x[d] )] + α 52 u = 0
1
d=1 x1
.
P
r Ix (x[d] + u)[I(x[d] + u) − I(x[d] )] + α 52 u2 = 0
d=1 2
Далее, производя аппроксимацию по формулам Тейлора, придем к
следующей системе
Pr I (x[d] )[I (x[d] ) + I (x[d] )u + I (x[d] )u ] − α 52 u = 0
t
x1
1
x2
2
1
d=1 x1
.
P
r Ix (x[d] )[It (x[d] ) + Ix (x[d] )u1 + Ix (x[d] )u2 ] − α 52 u2 = 0
1
2
d=1 2
К дискретной форме перейдем следующим образом. Положим u1ij и
[d]
u2ij как проекции скоростей по осям в точке с координатами i, j, а x1ij и
[d]
[d]
x2ij как координаты соответствующей точки xij . i ∈ 1..N − 1, j ∈ 1..M − 1.
Для границ i = 0, i = N, j = 0, j = M положим u1ij и u2ij равным граничным условиям, задаваемым изначально, обычно они берутся равными
нулю. Раскроем лапласиан в точке с координатами i, j по следующей формуле:
52 u1ij = u1i+1j + u1i−1j + u1ij−1 + u1ij−1 − 4u1ij
52 u2ij = u2i+1j + u2i−1j + u2ij−1 + u2ij−1 − 4u2ij
Тогда система уравнений принимает вид линейной системы уравнений:
Pr
[d]
[d]
[d]
[d]
Ix1 (xij )[It (xij ) + Ix1 (xij )u1ij + Ix2 (xij )u2ij ]−
d=1
− 1 α(u
1i+1j + u1i−1j + u1ij−1 + u1ij−1 − 4u1ij ) = 0
4
.
P
[d]
[d]
[d]
[d]
r
d=1 Ix2 (xij )[It (xij ) + Ix1 (xij )u1ij + Ix2 (xij )u2ij ]−
− 1 α(u2i+1j + u2i−1j + u2ij−1 + u2ij−1 − 4u2ij ) = 0
4
25
Введем обозначения:
A11ij =
A12ij =
A22ij =
B1ij =
B2ij =
[d]
[d]
d=1 Ix1 (xij )Ix1 (xij ),
Pr
[d]
[d]
d=1 Ix1 (xij )Ix2 (xij ),
Pr
[d]
[d]
d=1 Ix2 (xij )Ix2 (xij ),
Pr
[d]
[d]
d=1 Ix1 (xij )It (xij ),
Pr
[d]
[d]
d=1 Ix1 (xij )It (xij ).
Pr
Сгруппируем u1ij и u2ij :
(A
11ij + α)u1ij + A12ij u1ij = α(u1i+1j + u1i−1j + u1ij−1 + u1ij−1 ) − B1ij
.
A12ij u2ij + (A22ij + α)u2ij = α(u2i+1j + u2i−1j + u2ij−1 + u2ij−1 ) − B2ij
Будем решать систему итерационным методом:
(A
11ij
k−1
k−1
k−1
+ α)uk1ij + A12ij uk1ij = α(uk−1
1i+1j + u1i−1j + u1ij−1 + u1ij−1 ) − B1ij
A12ij uk + (A22ij + α)uk = α(uk−1 + uk−1 + uk−1 + uk−1 ) − B2ij
2ij
2ij
2i+1j
2i−1j
2ij−1
2ij−1
.
Пусть:
k−1
k−1
k−1
u
ek−1
= u1i+1j
+ uk−1
1
1i−1j + u1ij−1 + u1ij−1
k−1
k−1
k−1
u
ek−1
= uk−1
2
2i+1j + u2i−1j + u2ij−1 + u2ij−1
Тогда окончательное решение системы линейных уравнений запишется в
виде:
uk = u
ek−1
− (A11ij + A22ij + α2 )−1 (A11ij u
e1k−1 + A12ij u
ek−1
+ B1ij )
1ij
1
2
.
k−1
k−1
k−1
2 −1
uk = u
e2 − (A11ij + A22ij + α ) (A12ij u
e1 + A22ij u
e2 + B2ij )
2ij
При k = 0 u1ij и u2ij полагаются равным начальному приближению,
рационально брать это приближение равным нулю или равным результату,
полученным для предыдущего момента t.
Можно сделать следующее обобщение: считать, что r условий
u
I(x[d] + u) − I(x[d] ) → 0, x[d] , u ∈ Ωt0 .
имеют разные веса ω [d] , тогда задача минимизации запишется следующим
образом:
26
R
[d]
[d]
[d] 2
ω
(
I(x
+
u)
−
I(x
)
)dx1 dx2 +
d=1
Ωt 0
R
u
+α Ωt (k 5 u1 k2 + k 5 u2 k2 )dx1 dx2 → 0.
Pr
0
При дискретизации изменится только способ расчета коэффициентов
системы:
A11ij =
A12ij =
A22ij =
B1ij =
B2ij =
[d]
[d]
[d]
d=1 ω Ix1 (xij )Ix1 (xij ),
Pr
[d]
[d]
[d]
d=1 ω Ix1 (xij )Ix2 (xij ),
Pr
[d]
[d]
[d]
d=1 ω Ix2 (xij )Ix2 (xij ),
Pr
[d]
[d]
[d]
d=1 ω Ix1 (xij )It (xij ),
Pr
[d]
[d]
[d]
d=1 ω Ix1 (xij )It (xij ).
Pr
Для определенности будем полагать, что
r
X
ω [d] = 1
d=1
Эти коэффициенты можно выбирать различными способами. Самый простой ω [d] = 1/r. Если точки x[d] лежат в окрестности точки x, то более
оправдано брать для точек, ближних к x большие веса, а для дальних
меньшие.
На сетке, соответствующей пикселям изображения множество точек
n
o
[d]
xij , d = 1...r для точки xij можно выбрать следующим образом:
{xnm , n = −1...1, m = −1...1}
тогда формулы расчета коэффициентов системы линейных уравнений примут вид:
A11ij =
A12ij =
P1,1
n=−1,m=−1 Ix1 (xi+n,j+m )Ix1 (xi+n,j+m ),
P1,1
n=−1,m=−1 Ix1 (xi+n,j+m )Ix2 (xi+n,j+m ),
A22ij =
P1,1
B1ij =
P1,1
B2ij =
n=−1,m=−1 Ix2 (xi+n,j+m )Ix2 (xi+n,j+m ),
n=−1,m=−1 Ix1 (xi+n,j+m )It (xi+n,j+m ),
P1,1
n=−1,m=−1 Ix1 (xi+n,j+m )It (xi+n,j+m ).
27
§2 Модификация для выделения границ
Приведем математическое описание предлагаемой модификации для
выделения границ. Условие гладкости потока
u
k 5 u1 k2 + k 5 u2 k2 → 0, u ∈ Ωt0 .
позволяет в последствии получить в итерационном подходе невырожденную систему линейных уравнений. Но для этого требуется гладкость потока, что приводит к сильному размытию границ движущихся объектов на
изображении.
Z
Рассмотрим задачу минимизации:
Z
u
([I(x + u) − I(x)]2 )dx1 dx2 + α
(k 5 u1 k2 + k 5 u2 k2 )dx1 dx2 → 0.
Ωt0
Ωt 0
Ей соответствует система дифференциальных уравнений
I (x)[I (x) + I (x)u + I (x)u ] − α 52 u = 0
x1
t
x1
1
x2
2
1
.
Ix (x)[It (x) + Ix (x)u1 + Ix (x)u2 ] − α 52 u2 = 0
2
1
2
Для перехода к системе линейных уравнений в описанном ранее подходе использовались следующие формулы раскрытия лапласианов:
52 u1ij ≈ u1i+1j + u1i−1j + u1ij−1 + u1ij−1 − 4u1ij
52 u2ij ≈ u2i+1j + u2i−1j + u2ij−1 + u2ij−1 − 4u2ij
Модифицируем эти формулы. Будем считать, что для каждой точки известно являются ли соседние к ней точки граничные, т.е. известны
отображения для определения наличия границы у конкретной точки:
Gu = {Ωt0 → 0, 1} − для границы сверху,
Gd = {Ωt0 → 0, 1} − для границы снизу,
Gl = {Ωt0 → 0, 1} − для границы слева,
28
Gr = {Ωt0 → 0, 1} − для границы справа.
Значение 1 соответствует тому, что граница есть, 0 тому, что ее нет. Теперь
можно модифицировать формулы раскрытия лапласианов следующим образом:
52 u1ij ≈ Gd(xij )u1i+1j + Gu(xij )u1i−1j + Gl(xij )u1ij−1 + Gr(xij )u1ij−1 −
−(Gd(xij ) + Gu(xij ) + Gl(xij ) + Gr(xij ))u1ij
52 u2ij ≈ Gd(xij )u2i+1j + Gu(xij )u2i−1j + Gl(xij )u2ij−1 + Gr(xij )u2ij−1 −
−(Gd(xij ) + Gu(xij ) + Gl(xij ) + Gr(xij ))u2ij
если для краткости обозначить Gd(xij ) = Gdij , Gd(xij ) = Guij ,
Gd(xij ) = Glij и Gr(xij ) = Grij , а их сумму
Gsij = Gdij + Guij + Glij + Grij
можно нормировать Gsij = 1 то формулы раскрытия лапласианов примут
вид
52 u1ij ≈ Gdij u1i+1j + Guij u1i−1j + Glij u1ij−1 + Grij u1ij−1 − u1ij
52 u2ij ≈ Gdij u2i+1j + Guij u2i−1j + Glij u2ij−1 + Grij u2ij−1 − u2ij
В результате значения 52 u1ij и 52 u2ij будут не равны нулю, если
не равны нулю одновременно Gdij , Guij , Glij , Grij . Эти значения можно
вычислить для определенного момента времени t заранее и использовать в
итерационном решении системы
k−1
k−1
k−1
k−1
u
ek−1
1ij = Gdij u1i+1j + Guij u1i−1j + Glij u1ij−1 + Grij u1ij−1
k−1
k−1
k−1
k−1
u
ek−1
2ij = Gdij u2i+1j + Guij u2i−1j + Glij u2ij−1 + Grij u2ij−1
uk = u
ek−1
− (A11ij + A22ij + α2 )−1 (A11ij u
ek−1
+ A12ij u
ek−1
+ B1ij )
1ij
1
1
2
.
uk = u
ek−1 − (A11ij + A22ij + α2 )−1 (A12ij u
ek−1 + A22ij u
ek−1 + B2ij )
2ij
2
1
29
2
§3 Система линейных уравнений
Выпишем систему уравнений для u
B + A u + A u − α 52 u = 0
1ij
11ij 1ij
12ij 2ij
1ij
.
B2ij + A12ij u1ij + A22ij u2ij − α 52 u2ij = 0
рассмотрим способ раскрытия лапласианов, введенный в предыдущем пераграфе
B1ij + A11ij u1ij + A12ij u2ij −
−α(Gd u
ij 1i+1j + Guij u1i−1j + Glij u1ij−1 + Grij u1ij−1 − Gsij u1ij ) = 0
.
B
+
A
u
+
A
u
−
2ij
12ij 1ij
22ij 2ij
−α(Gdij u2i+1j + Guij u2i−1j + Glij u2ij−1 + Grij u2ij−1 − Gsij u2ij ) = 0
Таким образом мы пришли системе линейных уравнений порядка 2 ∗
w ∗h. Обозначим эту систему уравнений как H0 z0 = q0 . В случае w = h = 4
схема матрицы системы
..
A . B
H0 = . . .
. . . , A = A11ij , B = A12ij , C = A22ij , i, j = 1, wh
..
B . C
при первом способе раскрытия Лапласиана имеет вид, показанный на рисунке 4.
Первые w ∗ h уравнений отвечают u1ij , вторые w ∗ h уравнений отвечают u2ij . На схеме нулям соответствуют незакрашенные цвета, сплошной
штриховке отвечают ai,i = A11ij + Gsij α и ci,i = A22ij + Gsij α в верхней левой и правой нижней четверти соответственно. Диагональной штриховке
отвечает элемент −Gij α, Gij - один из Gdij ,Guij ,Glij ,Grij , а кругам элемент
bi,i = A12ij Вектор свободных членов будет иметь вид
q0 = (−B111 , . . . , −B1wh , −B211 , . . . , −B2wh )T .
30
Рис. 4: Схема матрицы.
Искомый вектор:
z0 = (u111 , . . . , u1wh , u211 , . . . , u2wh )T .
Уравнения и переменные можно перегруппировать следующим образом Hz = q:
H = {Hi,j }, i ∈ 1, wh, j ∈ 1, wh
aij bij
,
Hi,j =
bij cij
положим n = wh, тогда
z = (u111 , u211 )T , . . . , (u1wh , u2wh )T
T
= (z1 , . . . , zn )T ,
q = (−B111 , −B211 )T , . . . , (−B1wh , −B2wh )T
31
T
= (q1 , . . . , qn )T .
§4 Численные методы решения систем линейных уравнений
Для решения системы Hz = q можно использовать блочные численные методы [11-14], построенные на основе известных численных методов
Якоби, Гаусса-Зейделя и метода последовательной верхней релаксации.
Метод Якоби для данной задачи примет вид:
Hii zik+1
=−
n
X
Hij zjk + qi , i = 1, n, k = 0, ∞, z 0 = q,
j6=i,j=1
погрешность z k относительно точного решения z ∗ определяется формулой
||z k − z ∗ || < Qk ||z 0 − z 8 ||, где Q = ||H||
или с помощью нормы невязки
||rk || = ||q − Hz k || ≤ .
Метод Гаусса-Зейделя определяется формулой
Hii zik+1
=−
n
X
Hij zjk+1
j<i,j=1
−
n
X
Hij zjk + qi , i = 1, n, k = 0, ∞, z 0 = q,
j>i,j=1
погрешность:
||z k − z ∗ || ≤
Qr
||z k − z k−1 ||, где Q = ||H||, Qr = ||H∗ ||, H∗ = H −1 Hr ,
1−Q
при разложении матрицы на диагональную, верхнюю и нижнюю треугольную H = HL + HD + Hr , или с помощью нормы невязки
||rk || = ||q − Hz k || ≤ .
Метод последовательной верхней релаксации:
!
n
n
X
X
Hii zik+1 = ω −
Hij zjk+1 −
Hij zjk + qi +(1−ω)Hii zik , i = 1, n, k = 0, ∞,
j<i,j=1
j>i,j=1
Погрешность
||rk || = ||q − Hz k || ≤ .
32
§5 Сходимость блочных методов
Будем рассматривать блочную систему, введенную в третьем параграфе этой главы Hz = q.
Для доказательства сходимости блочных итерационных методов Якоби, Гаусса-Зейделя и метода последовательной верхней релаксации применительно к этой системе будем пользоваться следующими условиями сходимости [15]:
1. aii cii − b2ii > 0, aii > 0, cii > 0, i = 1, n
q
Pn
aii +cii
aii −cii 2
2. 2 ≥
+ b2ii и для некоторого i это
r=1,r6=i ||air || +
2
неравенство выполняется в строгом виде
3. матрица системы является блочно неприводимой
В [15] доказывается справедливость условий 1. и 2. для стандартной системы, которая используется в стандартном случае алгоритма HornSchunck. В обоих случаях матрицы системы линейных уравнений будут
блочно неприводимыми.
Приведем доказательство справедливости условий 1. и 2. для нашего
смешенного алгоритма. Для упрощения выкладок положим
Ix (xij ) = Uij , Iy (xij ) = Vij , i, j ∈ 1, 2.
Условие 1.
aii cii − b2ii =
2
X
!
Uij2 + Gsij α
i,j=1
=
2
X
i,j=1
Uij2 Gsij α+
2
X
Vij2 + Gsij α
Vij2 Gsij α+Gs2ij α2 +
i,j=1
Uij2
2
X
Vij2 −
!2
Uij Vij
=
i,j=1
2
X
Uij2
i,j=1
2
X
2
X
−
i,j=1
i,j=1
>
!
2
X
2
X
2
X
Vij2 −−
i,j=1
2
X
!2
Uij Vij
>
i,j=1
!2
Uij Vij
=
i,j=1
i,j=1
2
2
2
2
= (U11
+U21
+U12
+U22
)(V112 +V212 +V122 +V222 )−(U11 V11 +U21 V21 +U12 V12 +U22 V22 ) =
33
= (U11 V12 − U12 V11 )2 + (U11 V21 − U21 V11 )2 + (U11 V22 − U22 V11 )2 +
+(U12 V21 − U21 V12 )2 + (U12 V22 − U22 V12 )2 + (U21 V22 − U22 V21 )2 ≥ 0,
aii =
2
X
Uij2
+ Gsij α > 0, cii =
i,j=1
2
X
Vij2 + Gsij α > 0.
i,j=1
Условие 2.
aii + cii
≥
2
n
X
s
||air || +
r=1,r6=i
aii − cii
2
2
+ b2ii
Для точек, расположенных на границе изображения
Pn
r=1,r6=i ||air ||
<
Gsij α, что будет обеспечивать случай строгого выполнения неравенства в
P
условии, и nr=1,r6=i ||air || < Gsij α для точек, расположенных внутри изображения.
P2
2
2
+
U
i,j=1 Vij
i,j=1 ij
≥
Gsij α +
2
v
!2
!
u P2
P2
2
2 −
2 2
u
X
U
V
i,j=1 ij
i,j=1 ij
≥ Gsij α + t
+
Uij Vij
2
i,j=1
P2
+
2
P2
2
i,j=1 Uij +
2
P2
P2
2
i,j=1 Uij
P2
2
X
i,j=1
Uij2
2
X
i,j=1
2
i,j=1 Vij
2
i,j=1 Vij
Vij2 ≥
2
X
!2
−
2
P2
2
i,j=1 Uij −
2
P2
P2
2
i,j=1 Uij
≥
!2
P2
−
2
i,j=1 Vij
!2
+
2
X
!2
Uij Vij
i,j=1
2
i,j=1 Vij
!2
≥
2
X
!2
Uij Vij
i,j=1
!2
Uij Vij
, что уже рассматривалось в условии 1.
i,j=1
Конец доказательства.
34
Глава 3 Программная реализация вычисления оптического
потока
§1 Общая схема реализации вычислительной системы
В данной главе приводится описание построенного прототипа распределенной системы для распознавания движения на видео или последовательности изображений [16].
Прототип системы реализован в виде веб-сервиса. Общая схема вычислительного процесса представлена на рис. 6. Реализованный сервис предоставляет возможность вычисления поля скоростей для набора изображений
или кадров видиофайла.
Рис. 5: Общая схема вычислительного процесса.
Пользователь на вход подает несколько изображений или видеофайл,
во втором случае по указаным моментам времени из видеофайла извлекаются изображения. Для каждой пары из набора изображений необходимо
35
вычислить поле скоростей, матрицу проекций скоростей по осям x и y.
Запрос пользователя представляет собой отдельную задачу. Процесс вычисления состоит из следующих этапов:
1. В случае обработки видеофайла, извлечение по кадров из видео по
указанным пользователем временным меткам.
2. Выбор системой балансировки вычислительных узлов для решения
задачи.
3. Загрузка данных на вычислительные узелы.
4. Вычисление поля скоростей для задачи.
5. Предоставление результата пользователю в виде набора файлов.
Для каждой пары изображений предоставляется отдельный файл.
6. Пользователь может визуализировать полученный результат средствами, предоставляемыми на сайте.
36
§2 Средства реализации
Для реализации веб-сервиса были выбраны следующие средства.
Heroku
Облачный сервис Heroku [17] предоставляет вычислительные ресурсы
для размещения проектов, созданных на базе различных языков программирования: PHP, Java, Ruby, Python и некоторые другие. Каждый проект
запускается в отдельном контейнере, что обеспечивает его защищенность
и легкую масштабируемость.
В описываемой реализаци используется вариант однопоточных проектов для организации вычислительных сервисов, каждый из которых может вычислять поле скоростей для набора изображений, и сервиса-контроллера,
который будет предоставлять информацию о состоянии вычислительных
сервисов. Отдельный сервис служит для выделения кадров из видеофайла.
Все они реализованы с использованием языка программирования JAVA.
Доступ к вычислительным сервисам (вычислительные узлы) организуется через POST запрос.
Heroku DB
Кроме вычислительных ресурсов облачный сервис Heroku предоставляет возможность размещения базы данных под СУБД PostgreSQL, эта
возможность используется для хранения состояний вычислительных узлов. Каждый, из них, в начале вычисления задачи отмечет в базе то, что он
занят. При завершении вычисления задачи он отмечает то, что он освободился. Схема базы данных состояний вычислительных узлов представлена
на рис. 6.
Столбцы таблицы nodes имеют следующее назначение:
id - номер вычислительного узла 1 - 8
isused - занятость вычислительного узла true/false
instack - дополнительное поле, не используется
37
Рис. 6: Схема базы данных состояний вычислительных узлов.
Файловый сервер и веб-сервер
Основной сайт, через который предоставляется доступ к вычислительной системе представляет собой отдельный сервис, реализованный на
PHP, также развернутый на Heroku. Временное хранилище изображений,
полученных от пользователей, и результаты вычислений размещаются на
ftp сервисе byethost.
Визуализация на стороне клиента
Пользователь системы может визуализировать результат с помощью
скрипта javascript в браузере.
38
§3 Описание жизненного цикла задачи
Схема процесса обработки отдельной задачи представлена на рис. 7.
Пользователь выбирает на основном сайте набор изображений и запускает
процесс обработки задачи.
Рис. 7: Схема цикла обработки последовательности изображений.
В случае обработки видео предварительно происходит извлечение
необходимых кадров, которые представляют собой последовательность изображений.
Далее приведено описание шагов обработки задачи:
1. Коммуникационный сервер загружает видео или изображения на
ftp сервер.
2. В случае обработки видео, оно загружается на узел для раскадровки, дополнительно
3. Результат раскадровки сохраняется на ftp-сервер, коммуникацион39
ному серверу возвращаются ссылки на имена файлов изображений подлежащих обработке.
4. Коммуникационный сервер обращается к узлу-контроллеру (heroku)
для получения информации о состоянии вычислительных узлов.
5. Узел-контроллер соединяется с базой данных состояний вычислительных узлов (Heroku postgresql) и загружает информацию о состоянии
системы.
6. Узел-контроллер возвращает коммуникационному серверу информацию о состоянии вычислительных узлов в виде строки
1 , f a l s e ;2 , true ;3 , f a l s e ;4 , f a l s e ;
5 , true ;6 , f a l s e ;7 , f a l s e ;8 , f a l s e ;
7. Коммуникационный сервер выбирает узелs, на которыt подается
запрос на обработку с параметрами с параметрами count - количество изображений в последовательности, array - массив имен файлов.
8. Вычислительный узел указывает в базе данных состояний узлов
Heroku postgresql информацию о том, что он занят.
9. По этим ссылкам вычислительный узел загружает изображения из
ftp-сервера.
10. Вычислительный узел производит решение задачи: для каждой
пары изображений вычисляет поле скоростей в виде пары матриц в одном
файле данных.
11. Вычислительный узел сохраняет файлы результата на ftp-сервер.
12. Вычислительный узел указывает в базе данных состояний узлов
Heroku postgresql информацию о том, что он свободен.
13. Коммуникационный сервер получает ответ на POST запрос от
вычислительного узла, содержащий информацию о завершении обработки
задачи. Пользователю предоставляется результат в виде набора ссылок на
файлы результата, хранящиеся на ftp сервере.
40
14. Пользователь может скачать файлы результата с ftp-сервера или
просмотреть визуализацию потока на сайте.
Протокол HTTP используется при взаимодействии коммуникационного сервера с узлами Heroku и для общения с пользователем. Протокол
FTP используется при загрузке и сохранении выйлов на ftp-сервер вычислительными узлами и узлом раскадровки.
41
Глава 4 Примеры и анализ работы алгоритмов
§1 Пример: движение прямоугольника
Продемонстрируем работу алгоритма на нескольких примерах.
Рассмотрим движение черного прямоугольника на белом фоне рис.
8. Движение происходит на единичный вектор вниз.
Рис. 8: Прямоугольник.
Результат работы стандартного алгоритма Horn Schunck показан на
рисунке 9.
42
Рис. 9: Поток для прямоугольника расчитаный стандартным алгоритмом Horn Schunck.
Цифрами отмечена координатная сетка пикселей на изображении.
Результат работы стандартного алгоритма Лукаса-Канаде показан на
рисунке 10.
43
Рис. 10: Поток для прямоугольника расчитаный стандартным алгоритмом ЛукасаКанаде. Цифрами отмечена координатная сетка пикселей на изображении.
Результат работы модифицированного смешанного алгоритма Horn
Schunck Лукаса-Канаде показан на рисунке 11.
44
Рис. 11: Поток для прямоугольника расчитаный модифицированным смешанным алгоритмом Horn Schunck Лукаса-Канаде. Цифрами отмечена координатная сетка пикселей
на изображении.
Результат работы модифицированного смешанного алгоритма Horn
Schunck Лукаса-Канаде с модификацией выделения границ показан на рисунке 12.
45
Рис. 12: Поток для прямоугольника расчитаный модифицированным смешанным алгоритмом Horn Schunck Лукаса-Канаде с модификацией выделения границ. Цифрами
отмечена координатная сетка пикселей на изображении.
.
46
§2 Пример: сдвиг изображения
Рассмотрим сдвиг изображения на вектор 1, 1 вниз-влево.
Рис. 13: Тестовое изображение. Красной рамкой выделена область, для которой будет
показано поле скоростей.
47
Результат работы стандартного алгоритма Horn Schunck показан на
рисунке 14.
Рис. 14:
Поток для прямоугольника расчитаный стандартным алгоритмом Horn
Schunck. Цифрами отмечена координатная сетка пикселей на изображении.
48
Результат работы модифицированного смешанного алгоритма Horn
Schunck Лукаса-Канаде показан на рисунке 15.
Рис. 15:
Поток для прямоугольника расчитаный стандартным алгоритмом Horn
Schunck. Цифрами отмечена координатная сетка пикселей на изображении.
49
§3 Пример: сдвиг карты
Метод оптического потока также можно применять для анализа карт
местности, получаемых в результате аэрофотосьемки.
Рассмотрим сдвиг карты на вектор 1, 1 вниз-влево.
Рис. 16: Тестовое изображение. Красной рамкой выделена область, для которой будет
показано поле скоростей.
50
Результат работы стандартного алгоритма Horn Schunck показан на
рисунке 17.
Рис. 17:
Поток для прямоугольника расчитаный стандартным алгоритмом Horn
Schunck. Цифрами отмечена координатная сетка пикселей на изображении.
51
Результат работы модифицированного смешанного алгоритма Horn
Schunck Лукаса-Канаде показан на рисунке 18.
Рис. 18:
Поток для прямоугольника расчитаный стандартным алгоритмом Horn
Schunck. Цифрами отмечена координатная сетка пикселей на изображении.
52
§4 Анализ работы алгоритмов
Для проведения анализа работы алгоритмов необходимо выбрать величины, по которым производить сравнение. Ранее были описаны три примера, для каждого из них мы знаем заданные изначально векторы сдвига
и определенные алгоритмами их оценки. В качестве таких величин выберем следующие величины, которые будут отражать средние отклонения
вычисленных величин от реальных
iw ,jh
1 X
Eφ =
(φij − φ0 )
wh i,j=i ,j
0
0
iw ,jh
1 X
Er =
(uxij − ux )2 + (uyij − uy )2
wh i,j=i ,j
0
0
здесь область i = i0 ...iw , j = j0 ...jh фрагмент изображения размером
wh, 0, ux , uy - соответственно фазовый угол, и компоненты известного вектора сдвига изображения. φij , uxij , uyij - соответственно фазовый угол, и
компоненты вектора, вычисленного для точки i, j.
Далее для каждого примера приведем сравнение результатов работы
следующих алгоритмов:
1. Орригинального алгоритма Horn Schunck,
2. Модифицированного смешанного алгоритма Horn Schunck и ЛукасаКанаде, где точке с координатами i0 , j0 соответствует множество точек
{(i0 + i, j0 + j)|i ∈ {−1, 0, 1} , j ∈ {−1, 0, 1}}
3. Модифицированного смешанного алгоритма Horn Schunck и ЛукасаКанаде с подчеркиванием границ,
4. Орригинальный алгоритм Horn Schunck с подчеркиванием границ.
Результаты сравнения методов для примера сдвиг прямоугольника
вниз представлены на таблице 1, для прмера сдвиг изображения внизвправо на таблице 2 и сдвига карты на таблице 3.
53
Eφ - средняя ошибка фазового угла вектора поля скоростей
число итераций
5
оригинальный Horn-Shunck
20
40
60
1.568 1.548 1.545 1.544
смешанный Horn-Shunck и Лукаса-Канаде 1.562 1.540 1.539 1.537
смешанный Horn-Shunck и Лукаса-Канаде 1.565 1.550 1.547 1.543
с подчеркиванием границ
оригинальный Horn-Shunck
1.525 1.468 1.461 1.458
с подчеркиванием границ
Er - средняя ошибка длины вектора поля скоростей
число итераций
5
оригинальный Horn-Shunck
20
40
60
0.614 0.560 0.477 0.434
смешанный Horn-Shunck и Лукаса-Канаде 0.474 0.419 0.355 0.328
смешанный Horn-Shunck и Лукаса-Канаде 0.469 0.415 0.356 0.335
с подчеркиванием границ
оригинальный Horn-Shunck
0.610 0.539 0.462 0.417
с подчеркиванием границ
Таблица 1: Сравнение методов для сдвига прямоугольника
54
Eφ - средняя ошибка фазового угла вектора поля скоростей
число итераций
5
оригинальный Horn-Shunck
20
40
60
0.875 0.595 0.577 0.566
смешанный Horn-Shunck и Лукаса-Канаде 0.658 0.607 0.568 0.538
смешанный Horn-Shunck и Лукаса-Канаде 0.658 0.607 0.568 0.538
с подчеркиванием границ
оригинальный Horn-Shunck
0.875 0.566 0.595 0.577
с подчеркиванием границ
Er - средняя ошибка длины вектора поля скоростей
число итераций
5
оригинальный Horn-Shunck
20
40
60
2.728 2.641 2.500 2.376
смешанный Horn-Shunck и Лукаса-Канаде 2.457 1.865 1.515 1.394
смешанный Horn-Shunck и Лукаса-Канаде 2.457 1.865 1.515 1.394
с подчеркиванием границ
оригинальный Horn-Shunck
2.728 2.376 2.641 2.500
с подчеркиванием границ
Таблица 2: Сравнение методов для сдвига изображения
55
Eφ - средняя ошибка фазового угла вектора поля скоростей
число итераций
5
оригинальный Horn-Shunck
20
40
60
0.516 0.295 0.236 0.207
смешанный Horn-Shunck и Лукаса-Канаде 0.619 0.627 0.668 0.739
смешанный Horn-Shunck и Лукаса-Канаде 0.619 0.627 0.668 0.739
с подчеркиванием границ
оригинальный Horn-Shunck
0.516 0.295 0.236 0.207
с подчеркиванием границ
Er - средняя ошибка длины вектора поля скоростей
число итераций
5
оригинальный Horn-Shunck
20
40
60
2.637 2.498 2.347 2.246
смешанный Horn-Shunck и Лукаса-Канаде 2.268 1.877 1.795 1.578
смешанный Horn-Shunck и Лукаса-Канаде 2.268 1.877 1.795 1.578
с подчеркиванием границ
оригинальный Horn-Shunck
2.637 2.498 2.347 2.246
с подчеркиванием границ
Таблица 3: Сравнение методов для сдвига катры
56
Анализируя результаты можно сказать что смешанный алгоритм дает более точные оцеки поля скоростей, а подчеркивание границ эффективно не всегда.
57
Заключение
В работе было сделано следующее:
1. Был построен прототип системы распределенного вычисления поля
скоростей по видео и последовательности изображений.
2. Было предложено две модификации метода вычисления поля скоростей Horn-Schunck: смешанный алгоритм Horn-Schunck и Лукаса-Канаде
и модификация подчеркивания границ.
3. Для смешанного алгоритма Horn-Schunck и Лукаса-Канаде была
показана сходимость итерационных методов Якоби, Гаусса-Зейделя и метода последовательной верхней релаксации.
4. На примерах было продемонстрировано, что смешанный алгоритм
дает более эффективный результат, чем стандартный метода Horn-Schunck.
Таким образом, цель работы, заключавшаяся в построении системы
распределенного вычисления поля скоростей, достигнута.
Построенный прототип системы можно улучшить несколькими способами: ускорить вычислительный процесс, добавив распределение одной задачи на несколько узлов, усовершенствовать применяемый алгоритм, усовершенствовать отказоустойчивость системы при выходе из строя вычислительных узлов.
58
Список литературы
[1] Гонсалес Дж.Ту.Р. принципы распознавания образов. М.: Мир, 1978,
414 с.
[2] Uijlings J. R. R. Selective search for object recognition// International
Journal of Computer Vision. 2013. No. 104. P. 154–171.
[3] Bruce D. Lucas, and Takeo Kanade. An iterative image registration
technique with an application to stereo vision// IJCAI, 1981, P. 674–679.
[4] Horn B. K. P., Schunck B. G. Determining Optical Flow // Artificial
Intelligence. 1981. No 17. P. 185–203.
[5] Papenberg N. Highly Accurate Optic Flow Computation with Theoretically
JustifiedWarping // International Journal of Computer Vision. 2006. Vol. 2,
No 67. P. 141–158.
[6] Barron J.L., Fleet D.J. Beauchemin S.S. Performance of Optical Flow
Techniques // International Journal of Computer Vision. 1994. No 12.
P. 43–77.
[7] Fleet D.J., Weiss Y. Optical Flow Estimation // Mathematical models for
Computer Vision. 2005. P. 25.
[8] Black M.J. Robust Incremential Optical Flow // Yale University. 1992.
P.280.
[9] Anandan P. A computational framework and an algorithm for the
measurement of visual motion // International Journal of Computer Vision.
1989. P. 283—310.
[10] Котина Е. Д., Пасечная Г. А., Определение поля скоростей в задачах
59
обработки изображений, Изв. Иркутского гос. ун-та. Сер. Математика,
2013, 48–59 c.
[11] Фаддеев Д. К., Фаддеева В. Н. Вычислительные методы линейной алгебры, Численные методы линейной алгебры. Параллельные вычисления, Л.: Наука 1975, 228 с.
[12] Форсайт Дж. Молер К. Численное решение систем линейных алгебраических уравнений М.: Мир, 1969, 167 с.
[13] Воеводин В. В. Параллельные вычисления, СПб: БХВ, 2002, 602 с.
[14] Вержбицкий В. М. Основы численных методов: Учебник для вузов/
В. М. Вержбицкий. — М.: Высш. шк., 2002. — 840 с.
[15] Котина Е. Д. О сходимости блочных итерационных методов / Е. Д.
Котина// Изв. Иркут. гос. ун-та. – 2012. – Т. 5, вып. 3. – С. 41–55
[16] Электронный ресурс: http://velocityflow.herokuapp.com/ [дата обращения: 10.05.16]
[17] Электронный ресурс: https://www.heroku.com/ [дата обращения:
10.05.16]
60
Отзывы:
Авторизуйтесь, чтобы оставить отзыв