Повторяй за мной или двигаем виртуальной камерой like a pro / forpes.ru

Главная
Повторяй за мной или двигаем виртуальной камерой like a pro

Повторяй за мной или двигаем виртуальной камерой like a pro +4

01.06.2023 17:59

Coppermine 4 1400 Источник

Длинное введение

DVD – как много в этой аббревиатуре! Уверен, что вы наверняка помните такое явление, как ларьки и палатки с дисками, исчезнувшие только к началу 2010-х годов (по крайней мере так было в столице). В один из таких ларьков в конце 2009-го заглянул десятилетний я, внимание которого тут же привлекла коробка с надписью «3D Studio Max 2010»... Аниматором я, увы, так и не стал, однако интерес к области визуальных эффектов сохранился надолго.

В поисках обучающих материалов с англо-русским словарём наперевес в один прекрасный день я забрёл на Videocopilot. В одном из уроков автор с помощью неведомой волшебной софтины под названием Boujou показал, как отследить движение камеры в отснятом материале для его дальнейшего совмещения с трёхмерной графикой, что в подростковом мозгу произвело эффект разорвавшейся бомбы. Много позже, курсе на втором-третьем, меня всё чаще посещала навязчивая мысль – а как это вообще работает?

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

Первые шаги

Отправной точкой стал фундаментальный труд Роджерса и Адамса – «Математические основы машинной графики» 1980 года. Восстановлению трёхмерных объектов по проекциям посвящена отдельная глава, однако описанный подход был очень быстро отброшен, т. к.:

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

Следующим шагом, определившим дальнейшее направление, стали документация к OpenCV и монументальный труд Multiple View Geometry in Computer Vision (Hartley R., Zisserman A.). Также мне очень помогла эта статья. Данный материал практически полностью основан на этих источниках, если не указано иное.

Однородные координаты

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

Однородный вектор можно получить, просто добавив 1 к n-мерному декартову:

$\left ( \begin{matrix} x\\ y\\ z \end{matrix} \right )\rightarrow \left ( \begin{matrix} x\\ y\\ z\\ 1 \end{matrix} \right )$

Обратное преобразование выполняется путём деления каждой из компонент вектора на его последнюю компоненту с дальнейшим её отбрасыванием:

$\left ( \begin{matrix} X\\ Y\\ W \end{matrix} \right )\rightarrow \left ( \begin{matrix} \frac{X}{W}\\ \frac{Y}{W} \end{matrix} \right )$

Разумеется, если W ≠ 0.

Это всё замечательно, – скажете вы, – но хотелось бы больше конкретики. Пожалуйста! Вспоминаем алгебру за первый курс. Чтобы получить матрицу линейного преобразования нам нужно:

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

Таким образом мы получаем растяжение/сжатие и сдвиги по осям, но не параллельный перенос; произвольное аффинное преобразование будет выглядеть примерно так:

$\begin{pmatrix} a & b\\ c & d \end{pmatrix} \begin{pmatrix} x\\ y \end{pmatrix} + \begin{pmatrix} e\\ f \end{pmatrix} = \begin{pmatrix} ax+by+e\\ cx+dy+f \end{pmatrix}$

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

$\begin{pmatrix} a_2 & b_2\\ c_2 & d_2 \end{pmatrix} \left [ \begin{pmatrix} a_1 & b_1\\ c_1 & d_1 \end{pmatrix} \begin{pmatrix} x\\ y \end{pmatrix} + \begin{pmatrix} e_1\\ f_1 \end{pmatrix} \right ] + \begin{pmatrix} e_2\\ f_2 \end{pmatrix}$

В однородных координатах вышеописанный пример с преобразованием плоскости уже будет выглядеть так:

$\begin{pmatrix} a & b & e\\ c & d & f\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x\\ y\\ 1 \end{pmatrix} = \begin{pmatrix} ax+by+e\\ cx+dy+f\\ 1 \end{pmatrix}$

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

Камера-обскура

Пусть центр камеры находится в начале координат, камера сонаправлена оси OZ, а фокальная плоскость задана уравнением

где f – фокусное расстояние. Тогда точка в пространстве с координатами (X, Y, Z) отображается в точку

$\begin{pmatrix} x=\frac{fX}{Z} & y=\frac{fY}{Z} \end{pmatrix}$

Вспоминаем курс геометрии за шестой класс (перспективная проекция через подобие треугольников с коэффициентом f/Z).

В матричном виде это можно записать следующим образом:

$\begin{pmatrix} fX\\ fY\\ Z \end{pmatrix} = \begin{pmatrix} f & 0 & 0 & 0\\ 0 & f & 0 & 0\\ 0 & 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} X\\ Y\\ Z\\ 1 \end{pmatrix}$

Или:

$\begin{pmatrix} x\\ y\\ 1 \end{pmatrix} = \begin{pmatrix} f & 0 & 0\\ 0 & f & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} X\\ Y\\ Z\\ 1 \end{pmatrix}$

Всё бы ничего, но координаты точки на изображении выражены не в метрах/миллиметрах/попугаях, а в пикселях с отсчётом от верхнего левого угла – требуется перевод в новую систему координат (сдвиг, масштабирование). Для этого нам нужно знать не только фокусное расстояние, но и размеры кадра:

$\begin{pmatrix} u\\ v\\ 1 \end{pmatrix} = \begin{pmatrix} f_x & 0 & c_x\\ 0 & f_y & c_y\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} X\\ Y\\ Z\\ 1 \end{pmatrix},$

где

$f_x=\frac{f}{W}\cdot w,f_y=\frac{f}{H}\cdot h,$

W, H – физические размеры кадра, w, h – размеры кадра в пикселях, cx, cy – главная точка (пересечение фокальной плоскости с оптической осью) в координатах изображения. Коэффициенты f_x и f_y называются фокусными расстояниями по осям и могут отличаться друг от друга в случае неквадратных пикселей.

Осталось рассмотреть случай, когда камера повёрнута и/или смещена – требуется перейти из мировых координат в координаты камеры. Всё и сразу в матричном виде:

$\begin{pmatrix} u\\ v\\ 1 \end{pmatrix} = \begin{pmatrix} f_x & 0 & c_x\\ 0 & f_y & c_y\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} r_{11} & r_{12} & r_{13} & t_{1}\\ r_{21} & r_{22} & r_{23} & t_{2}\\ r_{31} & r_{32} & r_{33} & t_{3} \end{pmatrix} \begin{pmatrix} X\\ Y\\ Z\\ 1 \end{pmatrix}$

Или более кратко:

${\begin{pmatrix} u & v & 1 \end{pmatrix} }^T = K(R|T) {\begin{pmatrix} X & Y & Z & 1 \end{pmatrix}}^T,$

где K – матрица параметров камеры, R – матрица поворота, T – вектор смещения.

Геометрия нескольких видов

Рис. 2 – Дракон, видимый под различными углами зрения

Пусть имеются две камеры, как изображено на рисунке ниже. C и C’ – центры камер. Точка в пространстве X, C и C’ формируют эпиполярную плоскость. Пересечения этой плоскости с плоскостями изображений – эпиполярные прямые. Линия, соединяющая центры проецирования, пересекает плоскости изображения в сопряжённых точках, которые называются эпиполюсами.

Пусть точка X проецируется в x на плоскость изображения левой камеры и в x’ на плоскость изображения правой камеры. Тогда точке x на левом изображении соответствует эпиполярная прямая l’ на правом. При этом пара для x на изображении правой камеры может лежать только на соответствующей эпиполярной линии. То же самое работает и в обратную сторону.

Теперь всё то же самое, но в координатах:

$\exists F: x^TFx'=0, F=K^{'-T}RK^T\left [ KR^Tt \right ]_x$

F – фундаментальная матрица, где для вектора e обозначение [e]_x вычисляется как

$e_x = \begin{pmatrix} 0 & -e_z & e_y\\ e_z & 0 & -e_x\\ -e_y & e_x & 0 \end{pmatrix}$

Кроме фундаментальной матрицы также существует сущностная (или существенная? – видел оба варианта) матрица, выражающая то же ограничение:

$\exists E: x^TEx'=0, E=K'^TFK$

Разница между ними заключается в том, что сущностная матрица оперирует не пиксельными (как фундаментальная), а нормализованными координатами (имеют начало в оптическом центре изображения, безразмерны, т. к. нормализованы по f_x и f_y), а также имеет меньше степеней свободы.
С помощью сущностной матрицы можно вычислить движение камеры; пусть

$Udiag(\begin{matrix} 1 & 1 & 0 \\ \end{matrix}) V^T$

есть сингулярное разложение E, P₁ = (I|0), где U, V – ортогональные матрицы, I – единичная матрица, а P₁ и P₂ – матрицы проекций соответствующих камер. Тогда возможны 4 решения для P₂:

где

$W=\begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\ \end{pmatrix}, u_3=U\begin{pmatrix}0 & 0 & 1\end{pmatrix}^T ,$

но нас интересует только одно из них:

Рис. 4 – Геометрическая интерпретация решений

Наконец, найти координаты точки по её проекциям можно, зная матрицы проекций каждой из камер:

$\left\{\begin{matrix} x_1 = P_1X \\ x_2 = P_2X \\ \end{matrix}\right.,$

где: x₁, x₂ – проекции (в однородных координатах), P₁, P₂ – матрицы проекций камер.

Практика

Я знаю, что все снимки были сделаны камерой с APS-C матрицей (приблизительно 24×16 мм), фокусное расстояние объектива – 18 мм, дисторсия исправлена средствами камеры. Грубо прикинем матрицу параметров камеры:

# Оценка параметров камеры
FOCAL_EQ = 18.0
W, H = 24.0, 16.0  # Физические размеры сенсора в мм

cy, cx = (x / 2 for x in cam_list[0].img.shape[:-1])  # Центр изображения
# Фокусное расстояние в пикселях (прямоугольный случай)
fx = FOCAL_EQ * cam_list[0].img.shape[1] / W
fy = FOCAL_EQ * cam_list[0].img.shape[0] / H

# Матрица параметров камеры
K = np.array(((fx, 0, cx),
              (0, fy, cy),
              (0, 0, 1)), dtype=np.float32)

Особые точки

Сначала я попробовал ORB, однако методом научного тыка выяснилось, что детектор AKAZE даёт больше соответствий, хотя ORB работает быстрее. Далее нужно отфильтровать ложные совпадения, для чего ключевые точки должны пройти два теста, по результатам которых они будут приниматься или же не приниматься к рассмотрению:

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

Рис. 5 – Соответствий так много, что они сливаются в кашу

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

Восстановление движения и триангуляция

Основная часть магии происходит здесь. Я использую следующие функции OpenCV:

findEssentialMat – без комментариев;
recoverPose – для восстановления поворота и смещения из сущностной матрицы. Она возвращает матрицу поворота R (3×3) и нормализованный вектор смещения t (3×1). Тогда положение n-й камеры можно посчитать как

$T_n = T_{n-1}\cdot\begin{pmatrix} R & t \\ 0 & 1 \\ \end{pmatrix}$

src, dst = np.array(src, dtype=np.float32), np.array(dst, dtype=np.float32)
E, _ = cv.findEssentialMat(dst, src, K)
ret, local_R, local_t, mask = cv.recoverPose(E, dst, src, K)

T = np.hstack((local_R, local_t))
T = np.vstack((T, np.array((0, 0, 0, 1), dtype=np.float32)))
curr.T = np.dot(prev.T, T)  # Накапливаем преобразование

triangulatePoints – если точку видно в двух проекциях, то её можно восстановить с помощью этой функции. В качестве параметров функция принимает координаты точек и матрицы проекций для каждой камеры. Напомню, что матрицу проекции P (3×4) для n-й камеры можно найти следующим образом:

$P_n = K(R_{n}^{T}|-R_{n}^{T}\cdot t_n)$

R, t = curr.T[0:-1, 0:-1], curr.T[0:-1, 3]
# 0 = RC + t, C = -R.T * t, где C - центр камеры
P = np.hstack((R.T, np.dot(-R.T, t).reshape((3, 1))))
curr.P = np.dot(K, P)  # Новая матрица проекции

points_4d = cv.triangulatePoints(prev.P, curr.P, src.T, dst.T).T

Масштабирование

Помним, что вектор, который нам вернула recoverPose нормализован, т. е. мы знаем только направление движения. Пора это исправить.

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

Связь с внешним миром

Удовлетворяем любопытство

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

ply
format ascii 1.0
element vertex 6414
property float32 x
property float32 y
property float32 z
end_header
-2.0083678 -1.3574674  3.4676085
-1.7956636 -1.3710014  3.5252995
-1.9139004 -1.3578932  3.496787
...

Особенно хорошо получились обложки книг, а вот череп из-за равномерности окраски получился куда менее эпично.

Экспорт в в 3ds Max

Нет, использовать 3ds Max SDK в данном случае мы не будем, ибо из пушки по воробьям, для наших целей вполне достаточно возможностей MAXScript встроенного скриптового языка.

MAXScript поддерживает матрицы, но довольно специфичным образом. Класс Matrix3 описывает матрицы размером 4×3, которые задают транс-формацию объекта относительно родительской системы координат и содержат в себе локальную систему координат объекта. Каждая строка – вектор, определяющий направление осей X, Y, Z и начала координат в системе координат матрицы.

Очевидно, использовать полученные матрицы поворота напрямую не выйдет, но есть альтернатива в виде описания вращения с помощью кватернионов. Если не вдаваться в подробности и не говорить о гиперкомплексных числах, то кватернион – набор из четырёх чисел; первое (скалярная часть) определяет угол поворота вокруг вектора, оставшиеся три определяют вектор оси поворота.

point cross:off centermarker:on pos:[0.0, 0.0, 0.0] name:"origin"
animate on (
cam = freecamera name:"Camera"
cam.nearclip = 0.0
cam.parent = $origin
at time 0f cam.rotation = quat 1. 0. 0. 0.
at time 0f cam.position = [0.,0.,0.]
at time 1f cam.rotation = quat 0.99105 -0.0168   0.12458  0.04491
at time 1f cam.position = [0.96846,0.10677,0.22512]
at time 2f cam.rotation = quat 0.99587 0.06046 0.06039 0.03078
at time 2f cam.position = [ 0.40639,-0.53893, 0.11343]
...)

Заключение

Весь скрипт целиком вместе с изображениями, использованными в статье, можно найти на Github.

К сожалению, как и во всех остальных моих проектах, тут есть существенная недоделка: судя по форме облака точек (а оно не расслаивается, как могло бы), положение камер вычислено корректно, но при экспорте в 3ds Max корректное положение камера имеет только в первом кадре. Я сломал всю голову, привлёк друзей и знакомых, но решения так и не нашёл... Бонус – по-хорошему следует перейти от Y-down к Z-up системе координат графического пакета.

Отдельно хочу поблагодарить Диму, Алёну и команду Хабра за помощь с рецензированием и редактурой.

Комментарии (4)

Coppermine Автор
01.06.2023 18:01
#25608054
+1
Не хочу никого обидеть, но новый редактор - боль. Особенно что-то непонятное творится с формулами...

engine9
01.06.2023 18:01
#25608320
+2
Вдруг, если кому тема 3Д трекинга камеры интересна с практической точки зрения, в редакторе blender это реализовано на удивление хорошо. Можно на вход подавать как видеофайл так и секвенцию фотографий.

Jury_78
01.06.2023 18:01
#25608868
для его дальнейшего совмещения с трёхмерной графикой

Далек от этого, поэтому не понял... Это из плоского изображения сделать объемное?
1. Coppermine Автор
  01.06.2023 18:01
  #25608982
  Нет, речь про наложение графики на отснятый материал. Посмотрите этот ролик, думаю, сразу поймёте, о чём речь.