光流法的三维运动表示

Reference:

1. 高翔，张涛 《视觉SLAM十四讲》
2. GILAD ADIV Determining Three-Dimensional Motion and Structurefrom Optical Flow Generated by Several Moving Objects

1. 简介

光流是一种描述像素随时间在图像之间运动的方法。计算部分像素运动的方法称为稀疏光流；计算所有像素的称为稠密光流。稀疏光流以 Lucas-Kanade 光流为代表，稠密光流以 Horn-Schunck 光流为代表。本文中主要使用到了 LK 光流。

Lucas-Kanade 光流

单像素算法

在 LK 光流中，我们认为来自相机的图像是随时间变化的。图像可以看作时间的函数：$I(t)$。那么，一个在 $t$ 时刻，位于 $(x,y)$ 处的像素，它的灰度可以写成：
$$I(x,y,t)$$

光流法的基本假设：
灰度不变假设：同一个空间点的像素灰度值，在各个图像中是固定不变的。
对于 $t$ 时刻位于 $(x,y)$ 处的像素，我们设 $t+dt$ 时刻，它运动到 $(x+dx,y+dy)$ 处。由于灰度不变，有：
$$\mathbf{I}(x+\mathrm{d}x,y+\mathrm{d}y,t+\mathrm{d}t)=\mathbf{I}(x,y,t)$$
对左边进行泰勒展开，保留一阶项，得：
$$\mathbf{I}(x+\mathrm{d}x,y+\mathrm{d}y,t+\mathrm{d}t)\approx \mathbf{I}(x,y,t)+\frac{\partial \mathbf{I}}{\partial x}\mathrm{d}x+\frac{\partial \mathbf{I}}{\partial y}\mathrm{d}y+\frac{\partial \mathbf{I}}{\partial t}\mathrm{d}t$$

因为假设了灰度不变，于是下一个时刻的灰度等于之前的灰度，从而：
$$\frac{\partial \mathbf{I}}{\partial x}\mathrm{d}x+\frac{\partial \mathbf{I}}{\partial y}\mathrm{d}y+\frac{\partial \mathbf{I}}{\partial t}\mathrm{d}t=0$$

两边除以 $dt$，得：
$$\frac{\partial \mathbf{I}}{\partial x}\mathrm{d}x+\frac{\partial \mathbf{I}}{\partial y}\mathrm{d}y=-\frac{\partial \mathbf{I}}{\partial t}\mathrm{d}t$$

其中 $dx/dt$ 为像素在 $x$ 轴上运动速度，而 $dy/dt$ 为 $y$ 轴速度，把它们记为 $u$，$v$。同
时 $\partial I/\partial x$ 为图像在该点处 $x$ 方向的梯度，另一项则是在 $y$ 方向的梯度，记为 $I_x,I_y$。把图像灰度对时间的变化量记为 $I_t$，写成矩阵形式，有：
$$\left[\begin{array}{ll}{\mathbf{I}}_x& {\mathbf{I}}_y\end{array}\right]\left[\begin{array}{l}u\\ v\end{array}\right]=-{\mathbf{I}}_t$$

最小二乘获得运动轨迹

我们想计算的是像素的运动 u; v，但是该式是带有两个变量的一次方程，仅凭它无法
计算出 $u,v$。因此，必须引入额外的约束来计算 $u,v$。在 LK 光流中，我们假设某一个窗口内的像素具有相同的运动(注意是窗口，并不是整张图)。

考虑一个大小为 $w\times w$ 大小的窗口，它含有 $w^2$ 数量的像素。由于该窗口内像素具有
同样的运动，因此我们共有 $w^2$ 个方程：
$${\left[\begin{array}{ll}{\mathbf{I}}_x& {\mathbf{I}}_y\end{array}\right]}_k\left[\begin{array}{l}u\\ v\end{array}\right]=-{\mathbf{I}}_{tk},k=1,\dots ,w^2$$

记：
$$\mathbf{A}=\left[\begin{array}{c}{\left[{\mathbf{I}}_x,{\mathbf{I}}_y\right]}_1\\ ⋮\\ {\left[{\mathbf{I}}_x,{\mathbf{I}}_y\right]}_k\end{array}\right],\mathbf{b}=\left[\begin{array}{c}{\mathbf{I}}_{t1}\\ ⋮\\ {\mathbf{I}}_{tk}\end{array}\right]$$
于是整个方程为：
$$\mathbf{A}\left[\begin{array}{l}u\\ v\end{array}\right]=-\mathbf{b}$$

这是一个关于 u; v 的超定线性方程，传统解法是求最小二乘解。最小二乘在很多时候
都用到过：
$${\left[\begin{array}{l}u\\ v\end{array}\right]}^{\ast }=-{\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T\mathbf{b}$$

2. 三维运动表示

假设图像特征点 $p$ 的路标为 $P(X,Y,Z)$，它的图像点坐标是 $(x,y)$，它的光流向量是 $(u,v)$。相机有一个三维的移动，它的光心为 $O$，它的坐标轴为 $X_c,Y_c,Z_c$，它的图像平面为 $x_c,y_c$。路标 $P$ 位于图像平面的前面，它相对于相机的移动可以被分解成旋转移动 $\vec{\omega }=(A,B,C{)}^{\mathrm{T}}$ 和位移移动 $\vec{T}=(U,V,W{)}^{\mathrm{T}}$。

令 $P^{\prime }$ 为在时间 $t^{\prime }$ 上所对应的坐标系，对应关系可得：
$$\left(\begin{array}{l}{X}^{\prime }\\ Y^{\prime }\\ Z^{\prime }\end{array}\right)=R\left(\begin{array}{l}X\\ Y\\ Z\end{array}\right)+T$$
此时的旋转矩阵可以被估计，假设旋转参数的值较小，
$$R=\left(\begin{array}{ccc}1& -C& B\\ C& 1& -A\\ -B& A& 1\end{array}\right)$$

综合上式，可得 $X^{\prime }$, $Y^{\prime }$, 和 $Z^{\prime }$（注意，此时的图像已经矫正过了）
$$x^{\prime }=\frac{X^{\prime }}{Z^{\prime }}=\frac{x-Cy+B+U/Z}{-Bx+Ay+1+W/Z}$$

$$y^{\prime }=\frac{Y^{\prime }}{Z^{\prime }}=\frac{Cx+y-A+V/Z}{-Bx+Ay+1+W/Z}$$

已知位移向量 $(u=x^{\prime }-x,v=y^{\prime }-y)$，因此通过上面两式可以得到：
$$u=\frac{-Axy+B\left(1+x^2\right)-Cy+\left(U-Wx\right)/Z}{1+Ay-Bx+W/Z}$$

$$v=\frac{-A\left(1+y^2\right)+Bxy+Cx+\left(V-Wy\right)/Z}{1+Ay-Bx+W/Z}$$

如果 $|W/Z|<<1$ 且相机的视场不是很大，且旋转参数小，我们可以将分母近似为1。

因此图像三维运动的参考光流向量可以表示为：
$$\{\begin{array}{l}u=\frac{\dot{X}}{Z}-\frac{X\dot{Z}}{Z^2}=\left(\frac{U}{Z}+B-Cy\right)-x\left(\frac{W}{Z}+Ay-Bx\right)\\ v=\frac{\dot{Y}}{Z}-\frac{Y\dot{Z}}{Z^2}=\left(\frac{V}{Z}+Cx-A\right)-y\left(\frac{W}{Z}+Ay-Bx\right)\end{array}$$

其中，
$$\dot{X}=U+BZ-CY,\dot{Y}=V+CX-AZ,\dot{Z}=W+AY-BX$$