Михаил Смирнов (ВМК МГУ) Оптический поток с CUDA

2011 План Что такое оптический поток Классический метод Horn-Schunck Современный подход к реализации модели Horn-Schunck

2011

Оптический поток Набор векторов смещения – Плотный (для каждого пикселя) – Разреженный (для некоторого подмножества пикселей) Описывает движение в кадре

2011 Как найти оптический поток? Яркость объекта в кадре не меняется – Если никто не меняет освещения

2011 Как найти оптический поток? Яркость объекта в кадре не меняется – Если никто не меняет освещения Будем сопоставлять пиксели одинаковой яркости

2011

Не самый удачный вариант Нужна дополнительная информация

2011 Aperture Problem

2011 Дополнительная информация Посмотрим на движение в целом

2011 Дополнительная информация Посмотрим на движение в целом Соседние точки скорее всего двигаются одинаково

2011 Гладкость поля смещений

2011 Обозначения

2011 Гладкость поля смещений

2011 Гладкость поля смещений

2011 Гладкость поля смещений «Хорошее» поле – Небольшие абсолютные значения производных

2011 «Хорошее» поле смещений Сохраняет яркость Гладкое

2011 Обозначения

2011 Классическая модель Horn-Schunck

2011

Метод Якоби для дискретных уравнений Эйлера-Лагранжа

2011 Реализация Общая схема Копирование изображений в память GPU Вычисление производных изображения Решение СЛАУ Копирование результат в основную память

2011 Производные изображения

2011 Производные изображения Реализация Чтение через текстуры Нормализованные координаты cudaAddressModeMirror

2011 Производные изображения float t0, t1; // x derivative t0 = tex2D(texSource, x - 2.0f * dx, y); t0 -= tex2D(texSource, x - 1.0f * dx, y) * 8.0f; t0 += tex2D(texSource, x + 1.0f * dx, y) * 8.0f; t0 -= tex2D(texSource, x + 2.0f * dx, y); t0 /= 12.0f; t1 = tex2D(texTarget, x - 2.0f * dx, y); t1 -= tex2D(texTarget, x - 1.0f * dx, y) * 8.0f; t1 += tex2D(texTarget, x + 1.0f * dx, y) * 8.0f; t1 -= tex2D(texTarget, x + 2.0f * dx, y); t1 /= 12.0f; 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 Современный подход

2011 Деформация

2011 Общая схема

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 Современный подход Выполнить несколько итераций (warping iterations)

2011 Warping iterations Продолжить u, v решение Размер изображения

2011 Решение СЛАУ Оптимизация доступа к памяти Выравнивание адреса начала строк матрицы Использование текстур Использование общей памяти мультипроцессора (shared memory)

x480, 500 итераций G92 (GeForce 8800GTX) Обычный доступ2.39 с Shared memory0.54 с Текстуры (tex1Dfetch)0.58 с Решение СЛАУ Оптимизация доступа к памяти

x480, 500 итераций GF104 (GeForce GTX 460) Обычный доступ0.319 с Shared memory0.307 с Текстуры (tex1Dfetch)0.359 с Решение СЛАУ Оптимизация доступа к памяти

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 NVIDIA CUDA С Programming Guide

2011 Вопросы

Оптический поток

PDE = уравнения в частных производных Видео = …

Компьютерное зрение Одна из основных задач: Получить информацию о движении в кадре

Компьютерное зрение Одна из основных задач: Получить информацию о движении в кадре Идентифицировать движение Определить его направление Определить скорость

Компьютерное зрение Одна из основных задач: Получить информацию о движении в кадре Положение и перемещения камеры обычно неизвестны Идентифицировать движение Определить его направление Определить скорость

Компьютерное зрение Одна из основных задач: Получить информацию о движении в кадре Положение и перемещения камеры обычно неизвестны Идентифицировать движение Определить его направление Определить скорость Ограничиваем задачу поиском относительного движения объектов в системе отсчета, связанной с камерой

Компьютерное зрение Относительное движение между двумя последовательными кадрами можно представить как векторное поле – оптический поток

Оптический поток Применение Визуальные эффекты Отслеживание объектов Автономная навигация роботов Системы видеонаблюдения

Оптический поток Методы нахождения Фазовая корреляция преобразование Фурье только прямолинейное движение все точки перемещаются одинаково

Взаимная корреляция сигналов

Оптический поток Методы нахождения Фазовая корреляция преобразование Фурье только прямолинейное движение все точки перемещаются одинаково Блочные методы поиск похожих блоков все точки блока перемещаются одинаково

Оптический поток Методы нахождения Фазовая корреляция преобразование Фурье только прямолинейное движение все точки перемещаются одинаково Блочные методы поиск похожих блоков все точки блока перемещаются одинаково Вариационные методы минимизация некоторого функционала

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

Оптический поток Построение модели Дано: Непрерывная последовательность изображений Надо найти: Поле смещений в точкев момент временияркость

Оптический поток

Построение модели Предположение 1 Постоянство яркости пикселей Предположение 2 – малы, – гладкая функция В этом случае предыдущее равенство можно линеаризовать в точке

Построение модели Одно уравнение для двух неизвестных Некорректная задача с бесконечным числом решений Известно как проблема апертуры(aperture problem) наблюдая лишь за небольшой частью кадра невозможно однозначно определить направление движения

Построение модели Как движется шаблон? Одно уравнение для двух неизвестных Некорректная задача с бесконечным числом решений Известно как проблема апертуры(aperture problem) наблюдая лишь за небольшой частью кадра невозможно однозначно определить направление движения

Построение модели Как движется шаблон? По диагонали вправо? Одно уравнение для двух неизвестных Некорректная задача с бесконечным числом решений Известно как проблема апертуры(aperture problem) наблюдая лишь за небольшой частью кадра невозможно однозначно определить направление движения

Построение модели Как движется шаблон? По диагонали вправо? Только вправо? Одно уравнение для двух неизвестных Некорректная задача с бесконечным числом решений Известно как проблема апертуры(aperture problem) наблюдая лишь за небольшой частью кадра невозможно однозначно определить направление движения

Построение модели Как движется шаблон? По диагонали вправо? Только вправо? Или только вниз? Одно уравнение для двух неизвестных Некорректная задача с бесконечным числом решений Известно как проблема апертуры(aperture problem) наблюдая лишь за небольшой частью кадра невозможно однозначно определить направление движения

Как визуализировать векторное поле При помощи стрелок, предварительно понизив разрешение При помощи цвета: Направление – цвет Абсолютное значение – яркость

Как оценить качество рассчитанного потока Настоящий поток (ground truth) Вычисленный поток (estimated) Размер изображения Средняя угловая ошибка (AAE – Average Angular Error) Средняя ошибка по конечной точке (AEE – Average Endpoint Error)

Как оценить качество рассчитанного потока AAE = AEE = 7.22 AAE = 2.76 AEE = 0.37 AAE = 3.14 AEE = 1.53 Ground truth

Вариационные методы Основная идея поле смещений минимизирует некоторый энергетический функционал (data term) – штраф за отклонения от предположений о постоянстве какой-либо величины (например, яркости) (smooth term) – штраф за отклонения от предположений гладкости векторного поля параметр регуляризации – определяет гладкость получаемого векторного поля

МЕТОД HORN-SCHUNK B.K.P. Horn and B.G. Schunck, "Determining optical flow." Artificial Intelligence, vol 17, pp , 1981

Метод Horn-Schunck Предположение Поле смещений – гладкая функция в области Поле смещений минимизирует функционал data term – штраф за отклонения от предположения постоянства яркости smooth term – штраф за отклонение от предположения гладкости поля

Минимизация функционала Уравнения Эйлера-Лагранжа С граничными условиями

Метод Horn-Schunck Уравнения Эйлера-Лагранжа Граничные условия

Метод Horn-Schunck Аппроксимация оператора Лапласа

Метод Horn-Schunck Пространственные производные удобно вычислить заранее свертка столбцов(строк) с ядром 5x1 (1x5) Ядро

Метод Horn-Schunck Дискретизация уравнений Граничные условия уже учтены в дискретизованном уравнении Шаг сетки h обычно принимают равным 1 Система с 2xNxM неизвестными – метод Гаусса неприменим из-за высокой сложности и плохой устойчивости

Метод Якоби для решения СЛАУ Рассмотрим СЛАУ вида Итерационный метод

Метод Гаусса-Зейделя для решения СЛАУ Рассмотрим СЛАУ вида Итерационный метод

Метод релаксации (SOR) для решения СЛАУ Рассмотрим СЛАУ вида Итерационный метод

Метод Horn-Schunck Метод Якоби Просто распараллеливается Приближение на следующей итерации полностью определяется по приближению на текущей итерации – нет зависимости между разными пикселями

Метод Horn-Schunck Метод Гаусса-Зейделя i,j+1 i-1,ji,ji+1,j i,j-1 Нельзя напрямую использовать для параллельных вычислений

Красно-черная схема Гаусса- Зейделя Новое значение в красном узле зависит только от текущих значений в самом узле и в его «черных соседях» Новое значение в черном узле зависит только от текущих значений в самом узле и в его «красных соседях» Новые значения в узлах одного цвета можно вычислять параллельно Схема одной итерации Обновить значения в красных узлах Обновить значения в черных узлах Аналогично можно сделать параллельную версию SOR

Метод Horn-Schunck Схема вычислений Вычислить производные изображения (пространственные и временные) Установить начальное приближение - нулевой поток Выполнить некоторое количество итераций одного из итерационных методов Новые значения сохраняются в массив, отличный от массива со старыми значениями. В противном случае возникает конфликт чтения-записи. Ядро выполняет одну итерацию итерационного метода критерий останова

МЕТОД 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, , May 2004

Метод Brox et al Рассмотренные методы чувствительны к изменению освещения

Метод Brox et al Рассмотренные методы чувствительны к изменению освещения Идея Рассмотреть другой data term Предположение Сохраняется градиент изображения Плюс: градиент не чувствителен к аддитивному изменению яркости Минус: чувствительность к шуму

Метод Brox et al Рассмотренные методы чувствительны к изменению освещения Идея Рассмотреть другой data term Предположение Сохраняется градиент изображения С самого начала ищутся только небольшие векторы смещения

Метод Brox et al Рассмотренные методы чувствительны к изменению освещения Идея Рассмотреть другой data term Предположение Сохраняется градиент изображения С самого начала ищутся только небольшие векторы смещения Идея Использовать нелинеаризованное уравнение сохранения яркости

Метод Brox et al Тогда штраф за отклонение от предположений о сохранении яркости и градиента яркости примет вид

Метод Brox et al Проблема: квадратичный штраф чувствителен к выбросам

Метод Brox et al Проблема: квадратичный штраф чувствителен к выбросам Идея Заменить квадратичный штраф выпуклой возрастающей функцией малый положительный параметр

Метод Brox et al Получим новый функционал

Метод Brox et al Штраф за отклонения от предположения о гладкости потока Предположение Поток кусочно-гладкий

Метод Brox et al Задача Найти u и v, минимизирующие функционал

Метод Brox et al Уравнения Эйлера-Лагранжа Введем обозначения: Использование z вместо t говорит о том, что соответствующее выражение НЕ производная, а разность, которую будем минимизировать

Метод Brox et al Уравнения Эйлера-Лагранжа Граничные условия Неймана

Метод Brox et al Численный метод Уравнения Эйлера-Лагранжа нелинейные Используем метод неподвижной точки для w Если есть смещения на расстояние, большее размера пикселя, то функционал может иметь множество локальных минимумов Идея Использовать пирамиду изображений разного разрешения

Метод неподвижной точки

Метод Brox et al Численный метод Совместим метод неподвижной точки с масштабированием Масштабировать будем с коэффициентом Для более плавного перехода от разрешения к разрешению коэффициент выберем близким к 1 Используем полную пирамиду изображений, начиная с наименьшего возможного разрешения (например, с 10x10)

Метод Brox et al Численный метод Итерации по разрешению назовем внешними Пусть k – номер изображения в пирамиде (0 – самое низкое разрешение) – решение на k-ой внешней итерации Осталась нелинейность в

Метод Brox et al Численный метод Устранение нелинейности в производных изображения

Метод Brox et al Численный метод Введем обозначения

Метод Brox et al Численный метод Аппроксимация слагаемого с дивергенцией

Метод 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

Метод Brox et al Численный метод Оставшаяся нелинейность в Ѱ устраняется повторным применением метода неподвижной точки Начальные значения Полученная система уравнений линейная может быть решена одним из рассмотренных методов

Метод Brox et al Решение линейной системы Используем красно-черный метод релаксации (Red-Black SOR) Граничные условия

Метод Brox et al Решение линейной системы Используем красно-черный метод релаксации (Red-Black SOR)

Метод 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

Дополнительные материалы: