Transcript 2011

Оптический поток с CUDA
Михаил Смирнов (ВМК МГУ)
План
• Что такое оптический поток
• Классический метод Horn-Schunck
• Современный подход к реализации
модели Horn-Schunck
2011
2011
2011
2011
2011
Оптический поток
• Набор векторов смещения
– Плотный (для каждого пикселя)
– Разреженный (для некоторого
подмножества пикселей)
• Описывает движение в кадре
2011
Как найти оптический поток?
• Яркость объекта в кадре не меняется
– Если никто не меняет освещения
2011
Как найти оптический поток?
• Яркость объекта в кадре не меняется
– Если никто не меняет освещения
• Будем сопоставлять пиксели
одинаковой яркости
2011
2011
Не самый удачный вариант
• Нужна дополнительная информация
2011
Aperture Problem
2011
Дополнительная информация
• Посмотрим на движение в целом
2011
Дополнительная информация
• Посмотрим на движение в целом
• Соседние точки скорее всего
двигаются одинаково
2011
Гладкость поля смещений
2011
Обозначения
• Смещение по вертикали
–𝒗
• Смещение по горизонтали
–𝒖
2011
Гладкость поля смещений
𝝏𝒗
𝝏𝒚
>𝟎
<𝟎
2011
Гладкость поля смещений
𝝏𝒗
𝝏𝒚
>𝟎
<𝟎
2011
Гладкость поля смещений
• «Хорошее» поле
– Небольшие абсолютные значения
производных
2011
«Хорошее» поле смещений
• Сохраняет яркость
• Гладкое
2011
Обозначения
• Прямоугольная область изображения
–Ω
• Яркость изображения с номером 𝑡 в
точке (𝑥, 𝑦)
– 𝐼(𝑥, 𝑦, 𝑡)
2011
Классическая модель
Horn-Schunck
𝐸 𝑢, 𝑣 = 𝐸𝐷 𝑢, 𝑣 + 𝛼𝐸𝑆 𝑢, 𝑣 → 𝑖𝑛𝑓
• Сохранение яркости
– data term
𝐸𝐷 𝑢, 𝑣 =
𝐼 𝑥 + 𝑢 𝑥, 𝑦 , 𝑦 + 𝑣 𝑥, 𝑦 , 𝑡 + 1 − 𝐼 𝑥, 𝑦, 𝑡
2
𝑑𝑥𝑑𝑦
Ω
• Гладкость поля
– smoothness term
𝐸𝑆 𝑢, 𝑣 =
Ω
𝑢𝑥
2
+ 𝑢𝑦
2
+ 𝑣𝑥
2
2
+ 𝑣𝑦 𝑑𝑥𝑑𝑦
2011
Точка минимума 𝑬(𝒖, 𝒗)
• Линеаризуем подинтегральное
выражение в 𝐸𝐷 𝑢, 𝑣
• Уравнения Эйлера-Лагранжа для
линеаризованной модели
𝐼𝑥 𝑢 + 𝐼𝑦 𝑣 + 𝐼𝑡 𝐼𝑥 − 𝛼 𝑢𝑥𝑥 + 𝑢𝑦𝑦 = 0
𝐼𝑥 𝑢 + 𝐼𝑦 𝑣 + 𝐼𝑡 𝐼𝑦 − 𝛼 𝑣𝑥𝑥 + 𝑣𝑦𝑦 = 0
Граничные условия: отражение
2011
Метод Якоби для дискретных
уравнений Эйлера-Лагранжа
𝑛+1
𝑛
𝑢𝑖𝑗
= 𝑢𝑖𝑗
−
𝑛+1
𝑣𝑖𝑗
=
𝑛
𝑣𝑖𝑗
𝑛
𝑛
𝐼𝑥,𝑖𝑗 𝐼𝑥,𝑖𝑗 𝑢𝑖𝑗
+ 𝐼𝑦,𝑖𝑗 𝑣𝑖𝑗
+ 𝐼𝑡,𝑖𝑗
2
2
𝛼 + 𝐼𝑥,𝑖𝑗
+ 𝐼𝑦,𝑖𝑗
−
𝑛
𝑛
𝐼𝑦,𝑖𝑗 𝐼𝑥,𝑖𝑗 𝑢𝑖𝑗
+𝐼𝑦,𝑖𝑗 𝑣𝑖𝑗
+𝐼𝑡,𝑖𝑗
2 +𝐼 2
𝛼+𝐼𝑥,𝑖𝑗
𝑦,𝑖𝑗
где
𝑢𝑖𝑗 = 𝑢𝑖−1,𝑗 + 𝑢𝑖+1,𝑗 + 𝑢𝑖,𝑗−1 + 𝑢𝑖,𝑗+1
𝑣𝑖𝑗 = 𝑣𝑖−1,𝑗 + 𝑣𝑖+1,𝑗 + 𝑣𝑖,𝑗−1 + 𝑣𝑖,𝑗+1
2011
Реализация
Общая схема
• Копирование изображений в память
GPU
• Вычисление производных
изображения
• Решение СЛАУ
• Копирование результат в основную
память
2011
Производные изображения
• Пятиточечный шаблон:
– 1, −8,0,8, −1 ⋅
1
12
• Усреднение по времени
– 𝐼𝑥,𝑖𝑗 =
𝐼𝑥 𝑡+1 +𝐼𝑥 𝑡
2
𝑖𝑗
2011
Производные изображения
Реализация
• Чтение через текстуры
• Нормализованные координаты
• cudaAddressModeMirror
2011
Производные изображения
float t0, t1;
// x derivative
t0 = tex2D(texSource,
t0 -= tex2D(texSource,
t0 += tex2D(texSource,
t0 -= tex2D(texSource,
t0 /= 12.0f;
t1
t1
t1
t1
t1
=
-=
+=
-=
/=
tex2D(texTarget,
tex2D(texTarget,
tex2D(texTarget,
tex2D(texTarget,
12.0f;
x
x
x
x
+
+
2.0f
1.0f
1.0f
2.0f
*
*
*
*
dx,
dx,
dx,
dx,
y);
y) * 8.0f;
y) * 8.0f;
y);
x
x
x
x
+
+
2.0f
1.0f
1.0f
2.0f
*
*
*
*
dx,
dx,
dx,
dx,
y);
y) * 8.0f;
y) * 8.0f;
y);
Ix[pos] = (t0 + t1) * 0.5f;
// t derivative
Iz[pos] = tex2D(texTarget, x, y) - tex2D(texSource, x, y);
2011
Решение уравнений
Эйлера-Лагранжа
• Ядро выполняют одну итерацию
– Обновляются и 𝑢, и 𝑣
• Каждый поток обрабатывает вектор
для одного пикселя
• Исходные данные общие для всех
потоков
– Возможен выигрыш от использования
кеша
2011
Одна итерация
// handle borders
if (ix != 0) left = pos - 1;
else left = pos;
…
float sumU =
(u0[left] + u0[right] + u0[up] + u0[down]) * 0.25f;
float sumV =
(v0[left] + v0[right] + v0[up] + v0[down]) * 0.25f;
float frac = (Ix[pos] * sumU + Iy[pos] * sumV + Iz[pos])
/ (Ix[pos] * Ix[pos] + Iy[pos] * Iy[pos] + alpha);
u1[pos] = sumU - Ix[pos] * frac;
v1[pos] = sumV - Iy[pos] * frac;
2011
Современный подход
• Классический метод хорошо работает
только для небольших смещений
– Линеаризация
2011
Современный подход
• Уменьшим изображения
– «Большие» смещения станут
«маленькими»
• Решим задачу для уменьшенных
изображений
– Классический метод Horn-Schunck
• Используем найденное решение как
стартовую точку в исходной задаче
2011
Современный подход
• Можно рассмотреть пирамиду
изображений
– Каждое следующее меньше
предыдущего
• Используем решение задачи для
текущего уровня как начальное
приближение для следующего
– Стартуем с самых маленьких
изображений
2011
Современный подход
• Линеаризация работает для
небольших смещений
• Разобьем решение на текущем уровне
на два слагаемых
𝑢𝑛+1 = 𝑢𝑛 + 𝑑𝑢𝑛
𝑣 𝑛+1 = 𝑣 𝑛 + 𝑑𝑣 𝑛
𝑢𝑛 , 𝑣 𝑛 − решение с предущего уровня
2011
Деформация
• На текущем уровне вместо 𝐼(𝑥, 𝑦, 𝑡 + 1)
будем работать с 𝐼(𝑥 + 𝑢𝑛 , 𝑦 + 𝑣 𝑛 , 𝑡 + 1)
• Задача сводится к нахождению
небольших смещений 𝑑𝑢, 𝑑𝑣
– Классический Horn-Schunck
2011
Общая схема
• Копирование изображений в память GPU
• Построение пирамиды (коэффициент 0.5)
• 𝑢, 𝑣 ← 0
• Начиная с наименьшего разрешения
– Замена 𝐼 𝑡 + 1 на его деформированную версию 𝐼
– Вычисление производных (для 𝐼 𝑡 и 𝐼(𝑡 + 1))
– Нахождение 𝑑𝑢, 𝑑𝑣
– 𝑢 ← 𝑢 + 𝑑𝑢, 𝑣 ← 𝑣 + 𝑑𝑣
– Масштабирование 𝑢, 𝑣 (переход к большему разрешению)
• Копирование результатов в основную память
2011
Деформация
const int ix = threadIdx.x + blockIdx.x * blockDim.x;
const int iy = threadIdx.y + blockIdx.y * blockDim.y;
const int pos = ix + iy * stride;
if (ix >= width || iy >= height) return;
float x = ((float)ix + u[pos] + 0.5f) / (float)width;
float y = ((float)iy + v[pos] + 0.5f) / (float)height;
out[pos] = tex2D(texToWarp, x, y);
2011
Современный подход
• На каждом уровне изображение
можно деформировать несколько раз
– warping iterations
2011
Современный подход
• Создание пирамиды
• 𝑢, 𝑣 ← 0
• Для всех уровней пирамиды начиная с
наименьшего
– Деформация
– Вычисление производных
– Нахождение 𝑑𝑢, 𝑑𝑣
– 𝑢 ← 𝑢 + 𝑑𝑢, 𝑣 ← 𝑣 + 𝑑𝑣
Выполнить несколько
итераций
(warping iterations)
– Перейти к следующему уровню
2011
Размер изображения
Warping
iterations
решение
Warping
iterations
Продолжить u, v
Warping
iterations
Продолжить u, v
2011
Решение СЛАУ
Оптимизация доступа к памяти
• Выравнивание адреса начала строк
матрицы
• Использование текстур
• Использование общей памяти
мультипроцессора (shared memory)
2011
Решение СЛАУ
Оптимизация доступа к памяти
• 640x480, 500 итераций
• G92 (GeForce 8800GTX)
Обычный доступ
2.39 с
Shared memory
0.54 с
Текстуры (tex1Dfetch)
0.58 с
2011
Решение СЛАУ
Оптимизация доступа к памяти
• 640x480, 500 итераций
• GF104 (GeForce GTX 460)
Обычный доступ
0.319 с
Shared memory
0.307 с
Текстуры (tex1Dfetch)
0.359 с
2011
2011
2011
2011
Ресурсы нашего курса
• Steps3d.Narod.Ru
• Google Site CUDA.CS.MSU.SU
• Google Group CUDA.CS.MSU.SU
• Google Mail CS.MSU.SU
• Google SVN
• Tesla.Parallel.Ru
• Twirpx.Com
• Nvidia.Ru
2011
Полезные источники
• Black, M.J., & Sun, D. (2010). Secrets of Optical Flow
and Their Principles. CVPR 2010
• Horn, B.K., & Schunck, B.G. (1981) Determining
Optical Flow. Artificial Intelligence, 17, pp. 185-203.
• NVIDIA CUDA С Programming Guide
• http://vision.middlebury.edu/flow
2011
Вопросы
2011
Оптический поток
PDE = уравнения в частных
производных
( I z2   ( I xz2  I yz2 ))  ( I x I z   ( I xx I xz  I xy I yz ))    div(( 3u  3v )3u)  0,
2
2
( I z2   ( I xz2  I yz2 ))  ( I y I z   ( I yy I yz  I xy I xz ))    div(( 3u  3v )3v)  0
2
Видео = …
2
Компьютерное зрение
Одна из основных задач:
Получить информацию о движении в кадре
Компьютерное зрение
Одна из основных задач:
Получить информацию о движении в кадре
•Идентифицировать движение
•Определить его направление
•Определить скорость
Компьютерное зрение
Одна из основных задач:
Получить информацию о движении в кадре
•Идентифицировать движение
•Определить его направление
•Определить скорость
Положение и перемещения камеры обычно неизвестны
Компьютерное зрение
Одна из основных задач:
Получить информацию о движении в кадре
•Идентифицировать движение
•Определить его направление
•Определить скорость
Положение и перемещения камеры обычно неизвестны
Ограничиваем задачу поиском относительного движения объектов в
системе отсчета, связанной с камерой
Компьютерное зрение
Относительное движение между двумя последовательными кадрами можно
представить как векторное поле – оптический поток
Оптический поток
Применение
•Визуальные эффекты
•Отслеживание объектов
•Автономная навигация роботов
•Системы видеонаблюдения
Оптический поток
Методы нахождения
• Фазовая корреляция
• преобразование Фурье
• только прямолинейное движение
• все точки перемещаются одинаково
Взаимная корреляция
сигналов

( f  g )(t ) 
f
*
( y ) g ( y  t )dy

F { f  g}  ( F { f })  F {g}
*
Оптический поток
Методы нахождения
• Фазовая корреляция
• преобразование Фурье
• только прямолинейное движение
• все точки перемещаются одинаково
• Блочные методы
• поиск похожих блоков
• все точки блока перемещаются одинаково
Оптический поток
Методы нахождения
• Фазовая корреляция
• преобразование Фурье
• только прямолинейное движение
• все точки перемещаются одинаково
• Блочные методы
• поиск похожих блоков
• все точки блока перемещаются одинаково
• Вариационные методы
• минимизация некоторого функционала
E (u, v)  inf
x a  b  x a b  0
u  решение u  a  b  inf x  a  b  0
x
 E (u )  u  a  b  inf
Оптический поток
Вариационные методы
•формулировка в виде задачи оптимизации
•плотное поле скоростей
•для каждого пикселя
•высокая точность
•смещение не обязательно кратно размеру пикселя
Оптический поток
Построение модели
•Дано:
•Непрерывная последовательность изображений
I ( x, y, t ), ( x, y)  , t  0, T 
яркость в точке
в момент времени
•Надо найти:
•Поле смещений
 u ( x, y , t ) 


w ( x, y , t )   v ( x, y , t ) 


1


Оптический поток

I (,, t )
I (,, t  1)
v
u
w (,, t )
Построение модели
Предположение 1
Постоянство яркости пикселей
I ( x  u, y  v, t  1)  I ( x, y, t )  0
Предположение 2
u , v – малы, I – гладкая функция
В этом случае предыдущее равенство можно линеаризовать в точке ( x, y, t )
I ( x  u, y  v, t  1)  I ( x, y, t )  I x ( x, y, t )u  I y ( x, y, t )v  I t ( x, y, t ) 1
 I xu  I y v  I t  0
Построение модели
I xu  I y v  I t  0
•Одно уравнение для двух неизвестных
•Некорректная задача с бесконечным числом решений
•Известно как проблема апертуры(aperture problem)
•наблюдая лишь за небольшой частью кадра
невозможно однозначно определить направление
движения
Построение модели
I xu  I y v  I t  0
•Одно уравнение для двух неизвестных
•Некорректная задача с бесконечным числом решений
•Известно как проблема апертуры(aperture problem)
•наблюдая лишь за небольшой частью кадра
невозможно однозначно определить направление
движения
Как движется
шаблон?
Построение модели
I xu  I y v  I t  0
•Одно уравнение для двух неизвестных
•Некорректная задача с бесконечным числом решений
•Известно как проблема апертуры(aperture problem)
•наблюдая лишь за небольшой частью кадра
невозможно однозначно определить направление
движения
Как движется
шаблон?
По диагонали вправо?
Построение модели
I xu  I y v  I t  0
•Одно уравнение для двух неизвестных
•Некорректная задача с бесконечным числом решений
•Известно как проблема апертуры(aperture problem)
•наблюдая лишь за небольшой частью кадра
невозможно однозначно определить направление
движения
Как движется
шаблон?
По диагонали вправо?
Только вправо?
Построение модели
I xu  I y v  I t  0
•Одно уравнение для двух неизвестных
•Некорректная задача с бесконечным числом решений
•Известно как проблема апертуры(aperture problem)
•наблюдая лишь за небольшой частью кадра
невозможно однозначно определить направление
движения
Как движется
шаблон?
По диагонали вправо?
Только вправо?
Или только вниз?
Как визуализировать
векторное поле
При помощи стрелок,
предварительно понизив
разрешение
При помощи цвета:
•Направление – цвет
•Абсолютное значение – яркость
Как оценить качество
рассчитанного потока
Настоящий поток (ground truth)
Вычисленный поток (estimated)
wt
we
M N
Размер изображения
Средняя угловая ошибка (AAE – Average Angular Error)
t T
e

w
w
1
M
N
i, j
 i, j
AAE 
arccos
 
 wt we
NM i 1 j 1
 i, j i, j




Средняя ошибка по конечной точке (AEE – Average Endpoint Error)
AEE 
1
NM
 
M
i 1
N
t
e
w

w
i, j
i, j
j 1
Как оценить качество
рассчитанного потока
Ground truth
AAE = 3.14
AEE = 1.53
AAE = 32.55
AEE = 7.22
AAE = 2.76
AEE = 0.37
Вариационные методы
Основная идея
поле смещений минимизирует некоторый энергетический функционал
E(u, v)   D(u, v)    S (u, v)dxdy

D(u , v)
(data term) – штраф за отклонения от предположений о
постоянстве какой-либо величины (например, яркости)
S (u, v)
(smooth term) – штраф за отклонения от предположений
гладкости векторного поля

параметр регуляризации – определяет гладкость получаемого
векторного поля
МЕТОД HORN-SCHUNK
B.K.P. Horn and B.G. Schunck, "Determining optical
flow." Artificial Intelligence, vol 17, pp 185-203, 1981
Метод Horn-Schunck
Предположение
Поле смещений – гладкая функция в
области

Поле смещений минимизирует функционал


E (u, v)   I xu  I y v  I t    u  v dxdy

 



2
data
2
2
smooth
data term – штраф за отклонения от предположения постоянства
яркости
smooth term – штраф за отклонение от предположения гладкости
поля
Минимизация функционала
E (u, v)   F ( x, y, u, v, u x , u y , vx , v y )dxdy

Уравнения Эйлера-Лагранжа


0  Fu  Fu x  Fu y
x
y


0  Fv  Fvx  Fv y
x
y
С граничными условиями
 Fu x 
  0,
n 

 Fu y 
T
 Fv x 
n    0
 Fv y 
T
Метод Horn-Schunck
Уравнения Эйлера-Лагранжа
Граничные условия
 I x ( I x u  I y v  I t )    u  0

 I y ( I x u  I y v  I t )    v  0
2
2
 2  2
x
y
nT u  0, nT v  0
Метод Horn-Schunck
Аппроксимация оператора Лапласа
Ш (i, j )
ui , j 
ui 1, j  2ui , j  ui 1, j
h
2

ui , j 1  2ui , j  ui , j 1
h2
Метод Horn-Schunck
Пространственные производные удобно вычислить заранее
•свертка столбцов(строк) с ядром 5x1 (1x5)
f ( xi 2 , y j )  8 f ( xi 1 , y j )  8 f ( xi 1 , y j )  f ( xx  2 , y j )
f
( xi , y j ) 
x
12h
xi  2
Ядро
1
12

xi 1
xi
xi 1
8
12
0
8
12
xi  2

1
12
Метод Horn-Schunck
Дискретизация уравнений
u i , j  ui , j

2
0  I x ui , j  I x I y vi , j  I x I t    
2
h
( i , j )Ш  ( i , j )


vi , j  vi , j
0  I I u  I 2v  I I   

x y i, j
y i, j
y t
2

h
( i , j )Ш ( i , j )

Граничные условия уже учтены в дискретизованном уравнении
•Шаг сетки h обычно принимают равным 1
•Система с 2xNxM неизвестными – метод Гаусса неприменим из-за
высокой сложности и плохой устойчивости
Метод Якоби для решения
СЛАУ
Рассмотрим СЛАУ вида
Ax  b
A  L  D U
Итерационный метод
x0  0
x
k 1
1
k 1
i
 D (b  ( L  U ) x )  x
k

1 
k
  bi   aij x j 
aii 
j i

Метод Гаусса-Зейделя для
решения СЛАУ
Рассмотрим СЛАУ вида
Ax  b
A  L  D U
Итерационный метод
x0  0
x
k 1
1
k 1
i
 ( D  L) (b  Ux )  x
k

1 
k
k 1
  bi   aij xi   aij xi 
aii 
j i
j i

Метод релаксации (SOR) для
решения СЛАУ
Ax  b
Рассмотрим СЛАУ вида
A  L  D U
Итерационный метод
x0  0
x k 1  ( D  L) 1 (b  (U  (  1) D) x k ) 
k 1
i
x
 (1   ) x 
k
i
 
bi   aij x   aij x

aii 
j i
j i
k
j
k 1
j




Метод Horn-Schunck

  I x It

k 1
u 

1 k  
k

  I x I y vi , j   
u
2 i, j  
h
( i , j )Ш  ( i , j )


1
I x2   
2
( i , j )Ш  ( i , j ) h

  I y It

k 1
v 


1
k
k
  I x I y ui , j   
v 
2 i, j  
( i , j )Ш  ( i , j ) h


1
I y2   
2
( i , j )Ш  ( i , j ) h
i, j
Метод Якоби
i, j
Просто распараллеливается
Приближение на следующей итерации полностью определяется по
приближению на текущей итерации – нет зависимости между разными
пикселями
Метод Horn-Schunck

  I x It

k 1
u 
i, j
Метод
Гаусса-Зейделя

  I y It

k 1
v 
i, j


1
1
k
k 1
k
  I x I y vi , j   
u  
u 
2 i, j
2 i, j  
( i , j )Ш  ( i , j ) h
( i , j )Ш  ( i , j ) h


1
I x2   
2
( i , j )Ш  ( i , j ) h

1 k 1
1 k  
k 1

  I x I y ui , j   
v  
v
2 i, j
2 i, j  
( i , j )Ш  ( i , j ) h
( i , j )Ш  ( i , j ) h


1
I y2   
2
h

( i , j )Ш ( i , j )
i,j+1
i-1,j
i,j
i,j-1
i+1,j
Ш  (i, j )
Ш  (i, j )
Нельзя напрямую использовать для параллельных
вычислений
Красно-черная схема ГауссаЗейделя
•Новое значение в красном узле
зависит только от текущих значений в
самом узле и в его «черных соседях»
•Новое значение в черном узле
зависит только от текущих значений в
самом узле и в его «красных соседях»
Новые значения в узлах одного
цвета можно вычислять
параллельно
Схема одной итерации
•Обновить значения в красных узлах
•Обновить значения в черных узлах
Аналогично можно сделать параллельную версию SOR
Метод Horn-Schunck
•Ядро выполняет одну итерацию итерационного метода
Схема вычислений
•Вычислить производные изображения (пространственные и
временные)
•Установить начальное приближение - нулевой поток
•Выполнить некоторое количество итераций одного из итерационных
методов
Новые значения сохраняются в массив, отличный от массива
со старыми значениями. В противном случае возникает конфликт
чтения-записи.
max xin 1  xin  
i


max Axn 1 i  bi  
i
критерий останова
МЕТОД BROX ET AL
T. Brox, A. Bruhn, N. Papenberg, J. Weickert
High accuracy optical flow estimation based on a theory for
warping,
T. Pajdla and J. Matas (Eds.), European Conference on Computer
Vision (ECCV) Prague, Czech Republic, Springer, LNCS, Vol. 3024, 2536, May 2004
Метод Brox et al
Рассмотренные методы чувствительны к изменению освещения
Метод Brox et al
Рассмотренные методы чувствительны к изменению
освещения
Идея
•Рассмотреть другой data term
Предположение
Сохраняется градиент изображения
I ( x  u, y  v, t  1)  I ( x, y, t )  0
Плюс:
•градиент не чувствителен к аддитивному изменению яркости
Минус:
•чувствительность к шуму
Метод Brox et al
Рассмотренные методы чувствительны к изменению
освещения
Идея
•Рассмотреть другой data term
Предположение
Сохраняется градиент изображения
I ( x  u, y  v, t  1)  I ( x, y, t )  0
С самого начала ищутся только небольшие векторы смещения
Метод Brox et al
Рассмотренные методы чувствительны к изменению
освещения
Идея
•Рассмотреть другой data term
Предположение
Сохраняется градиент изображения
I ( x  u, y  v, t  1)  I ( x, y, t )  0
С самого начала ищутся только небольшие векторы смещения
Идея
Использовать нелинеаризованное уравнение сохранения яркости
Метод Brox et al
Тогда штраф за отклонение от предположений о сохранении яркости и
градиента яркости примет вид


EData (u, v)   I (x  w)  I (x)   I (x  w)  I (x) dxdy

2
2
Метод Brox et al
Проблема:
•квадратичный штраф чувствителен к выбросам
Метод Brox et al
Проблема:
•квадратичный штраф чувствителен к выбросам
Идея
Заменить квадратичный штраф выпуклой возрастающей функцией
( s )  s  
2

2
малый положительный параметр
2
Метод Brox et al
Получим новый функционал


EData (u, v)    I (x  w)  I (x)   I (x  w)  I (x) dxdy

2
2
Метод Brox et al
Предположение
Поток кусочно-гладкий
Штраф за отклонения от предположения о гладкости потока


ESmooth (u, v)    3u  3v dxdy

2
2
Метод Brox et al
Задача
Найти u и v, минимизирующие функционал
E(u, v)  EData (u, v)    ESmooth (u, v)


EData (u, v)    I (x  w)  I (x)   I (x  w)  I (x) dxdy

2

2

ESmooth (u, v)    3u  3v dxdy

2
2
Метод Brox et al
Уравнения Эйлера-Лагранжа
Введем обозначения:
I x :  x I (x  w )
I y :  y I (x  w )
I z : I (x  w )  I (x)
I xx :  xx I (x  w )
I yy :  yy I (x  w )
I xz :  x I (x  w )   x I (x)
I yz :  y I (x  w )   y I (x)
Использование z вместо t говорит о том, что соответствующее
выражение НЕ производная, а разность, которую будем
минимизировать
Метод Brox et al
Уравнения Эйлера-Лагранжа
( I z2   ( I xz2  I yz2 ))  ( I x I z   ( I xx I xz  I xy I yz )) 
   div((  3u   3v ) 3u )  0,
2
2
( I z2   ( I xz2  I yz2 ))  ( I y I z   ( I yy I yz  I xy I xz )) 
   div((  3u   3v ) 3v)  0
2
Граничные условия Неймана
2
Метод Brox et al
Численный метод
•Уравнения Эйлера-Лагранжа нелинейные
•Используем метод неподвижной точки для w
•Если есть смещения на расстояние, большее размера пикселя, то
функционал может иметь множество локальных минимумов
Идея
•Использовать пирамиду изображений
разного разрешения
Метод неподвижной точки
f ( x)  x
x
n 1
 f (x )
n
 lim x
n 
n 1
x
 f (x )  x
*
*
*
Метод Brox et al
Численный метод
•Совместим метод неподвижной точки с масштабированием
•Масштабировать будем с коэффициентом   (0,1)
•Для более плавного перехода от разрешения к разрешению
коэффициент выберем близким к 1
•Используем полную пирамиду изображений, начиная с
наименьшего возможного разрешения (например, с 10x10)
Метод Brox et al
Численный метод
Итерации по разрешению назовем внешними
Пусть k – номер изображения в пирамиде (0 – самое низкое разрешение)
w k  (u k , v k ,1)T – решение на k-ой внешней итерации
' ((I zk 1 ) 2   ((I xzk 1 ) 2  ( I yzk 1 ) 2 ))  ( I xk I zk 1   ( I xxk I xzk 1  I xyk I yzk 1 )) 
   div(( 3u
k 1 2
  3v
Осталась нелинейность в
k 1
•
I*
•

k 1 2
)3u k 1 )  0
Метод Brox et al
Численный метод
Устранение нелинейности в производных изображения
I zk 1  I zk  I xk duk  I yk dvk
I xzk 1  I xzk  I xxk duk  I xyk dvk
k
I yzk 1  I yzk  I xyk duk  I yy
dvk
u k 1  u k  duk , v k 1  v k  dvk
Метод Brox et al
Численный метод
Введем обозначения
() kData : ((I zk  I xk duk  I yk dvk ) 2 
k
  ((I xzk  I xxk duk  I xyk dvk ) 2  ( I yzk  I xyk duk  I yy
dvk ) 2 ))
()
k
Smooth
2
2
: ( 3 (u  du )  3 (v  dv ) )
k
k
k
k
Метод Brox et al
Численный метод
Аппроксимация слагаемого с дивергенцией
div(k  u )i , j 
ui 1, j  ui , j
ui , j  ui 1, j
1
  ki  1 , j
 ki  1 , j
2
h 2
h
h

 

ui , j 1  ui , j
ui , j  ui , j 1 
1

  ki , j  1
 ki , j  1
2
2
h
h
h

Метод Brox et al
Численный метод
Аппроксимация smooth term
Метод Brox et al
Численный метод
Аппроксимация smooth term
Метод Brox et al
Численный метод
Аппроксимация smooth term
Метод Brox et al
Численный метод
Аппроксимация smooth term
Метод Brox et al
Численный метод
Аппроксимация smooth term
u i  1 , j  ui 1, j  ui , j 
2
2
2
u i , j  1  ui , j 1  ui , j 
2
2
2
  ui 1, j 1  ui , j 1   ui 1, j 1  ui , j 1  
  
 
  

2
2






  ui 1, j 1  ui 1, j
  
2

  ui 1, j 1  ui 1, j
  
2
 

 


2
2
Метод Brox et al
Численный метод
Оставшаяся нелинейность в Ѱ устраняется повторным применением
метода неподвижной точки
Начальные значения
w0  (0,0,1)T , duk ,0  0, dvk ,0  0
,l
0  () kData
 ( I xk ( I zk  I xk duk ,l 1  I yk dvk ,l 1 ) 
k
   I xxk ( I xzk  I xxk duk ,l 1  I xyk dvk ,l 1 )    I xyk ( I yzk  I xyk duk ,l 1  I yy
dvk ,l 1 )) 
,l
   div () kSmooth
3 (u k  duk ,l 1 ) 
Полученная система уравнений
•линейная
•может быть решена одним из рассмотренных методов
Метод Brox et al
Решение линейной системы
Используем красно-черный метод релаксации (Red-Black SOR)
duik ,l ,m1  (1   )duik .l ,m 


k ,l
k
k ,l , m
k ,l
k


(

)
(
u

du
)

(

)
u

S
i

j
j
j
S
i

j
i
jШ  ( i )
jШ  ( i )
 jШ (i ) (S )ik,l j 
(D )ik ,l


( )
k ,l
D i


((I xk,i ) 2    ((I xyk ,i ) 2  ( I xxk ,i ) 2 ))
k
k ,l , m
((I xk,i I yk,i    ( I xxk ,i I xyk ,i  I xyk ,i I yy
 I xk,i I zk,i    ( I xxk ,i I xzk ,i  I xyk ,i I yzk ,i ))
,i ))dvi

jШ  ( i )
( )
k ,l
S i j

(D )ik ,l

((I xk,i ) 2    ((I xyk ,i ) 2  ( I xxk ,i ) 2 ))
Граничные условия
du k ,l ,m

 0, dv k ,l ,m

0
Метод Brox et al
Решение линейной системы
Используем красно-черный метод релаксации (Red-Black SOR)
duik ,l ,m 1  (1   )duik .l ,m 



jШ  ( i )
k ,l

(

)
S
i
j 
jШ  ( i )
(D )ik ,l

(S )ik,l j (u kj  dukj ,l ,m )   jШ (i ) (S )ik,l j uik

(D )ik ,l


((I xk,i ) 2    ((I xyk ,i ) 2  ( I xxk ,i ) 2 ))
k
k ,l , m
((I xk,i I yk,i    ( I xxk ,i I xyk ,i  I xyk ,i I yy
 I xk,i I zk,i    ( I xxk ,i I xzk ,i  I xyk ,i I yzk ,i ))
,i ))dvi

jШ  ( i )
( )
k ,l
S i j

(D )ik ,l

((I xk,i ) 2    ((I xyk ,i ) 2  ( I xxk ,i ) 2 ))
Метод Brox et al
Билинейная интерполяция для I(x+w)
Метод Brox et al
Реализация
•Кадры храним в виде текстур
•Для продолжения решения с разрешения k на разрешение k+1
используем текстуры
Общая схема метода
Подготовить текстуры с изображениями
Задать начальное значения для u и v
•Для каждого разрешения, начиная с наименьшего
•Установить начальное значение du, dv
•Подготовить данные для итераций SOR
•Выполнить некоторое число итераций SOR
•Обновить u, v
•Продолжить решение на большее разрешение
•Сохранить результат
Метод Brox et al
Реализация
•Кадры храним в виде текстур
•Для продолжения решения с разрешения k на разрешение k+1
используем текстуры
Общая схема метода
Подготовить текстуры с изображениями
Задать начальное значения для u и v
•Для каждого разрешения, начиная с наименьшего
•Установить начальное значение du, dv
•Подготовить данные для итераций SOR
•Выполнить некоторое число итераций SOR
•Обновить u, v
•Продолжить решение на большее разрешение
•Сохранить результат
warping FP iteration
lagged nonlinearity FP
iteration
solver FP iteration
Метод Brox et al
•Сложные шаблоны доступа к памяти
•Большая часть данных не меняется во время итераций
Используем текстуры
Текстуры + shared memory работает медленнее чем просто текстуры
Для объединения запросов при записи индексы начала строк округляем
до значений, кратных 16
Дополнительные материалы:
•http://www.google.com
•http://www.mia.uni-saarland.de
•http://vision.middlebury.edu/flow/