KLT光流跟踪

本文主要会讲解以下的几个部分:

0. 什么是光流

1. 光流法的基本约束方程

2. KLT算法的基本思想

3. 公式推导过程

4. opencv代码实现

5. 参考部分

希望能够通过这个算法一窥光流跟踪法的基本思想与原理。
feature_points_tracking

0. 什么是光流

光流(optical flow)是空间运动物体在观察成像平面上的像素运动的瞬时速度[4]。
光流( Optical Flow )的概念是Gibson在1950年首先提出来的,它是空间运动物体在观察成像平面上的像素运动的瞬时速度,是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系,从而计算出相邻帧之间物体的运动信息的一种方法。一般而言,光流是由于场景中前景目标本身的移动、相机的运动,或者两者的共同运动所产生的。[0]

1.光流法的基本约束方程

光流法的基本思想是跟踪图像中像素点的移动来进行对物体移动的跟踪,而这种跟踪是发生在图像的像素级别的。我们用$I(x,y,t)$ 来表示图像中处于t时刻点(x,y)处的灰度值,经过$d_t$时间的移动该点的位置发生了改变,那么此时该像素点的表达方式变为$I(x +d_x,y+d_y,t+d_t)$ (未发生较大的运动,x方向上移动了$d_x$,y方向上移动了$d_y$).
我们假设图像的亮度没有发生改变,则有:
$$I(x,y,t) = I(x +d_x,y+d_y,t+d_t) \tag{1}$$
将上式泰勒展开有:
$$
I(x,y,t)=I(x,y,t)+
\frac {\partial I} {\partial x} d_x +
\frac {\partial I} {\partial y} d_y +
\frac {\partial I} {\partial t} d_t +
\xi \tag{2}
$$

$\xi$表示高阶无穷小项,可以忽略不计。再将(2)代人(1)后同除dt,有
$$
\frac {\partial I} {\partial x} {\frac {d_x} {d_t}} +
\frac {\partial I} {\partial y} {\frac {d_y} {d_t}} +
\frac {\partial I} {\partial t} {\frac {d_t} {d_t}}
$$
进行符号替换有:
$$I_x u + I_yv +I_t = 0 \tag{3}$$

(u,v)即为所求的光流矢量
式3被称为基本约束方程,如前文所述,这个方程成立有以下的两个条件:

  1. 亮度恒定
  2. 时间连续或者是运动是“小运动”,相邻帧之间位移要比较小.

2. KLT算法的基本思想

上面我们已经介绍了基本约束方程,等式(3)有两个未知量,需要引入其他的约束条件才能求解。
KLT算法就是光流法中基于梯度的一种方法,它通过增加假设空间一致,临近点有相似运动,保持相邻这一条件,来求解方程(3)进而求解出光流。
下面我们将通过学习KLT算法来探究这种像素级别的跟踪方式。

feature_points_tracking

算法的基本思想是:通过跟踪连续的图像帧中特定的点达到跟踪的效果(这种同方向上移动的点总是有很多,且有相同的移动轨迹)。这种方法最大的问题是如何求解出第二帧中的对应的点,且要求相邻帧之间的移动不能太剧烈。

将上述思想抽象成具体的数学问题就是:

“如何找到两个帧之间相似的点”

我们用$I和J$ 来表示两个不同的帧,$I(x y)$和$J(x y)$来表示图片中$[x y]^T$坐标处的灰度值(gray value),使用$u = [u_x u_y]^T$来表示第一帧$[x, y]^T$坐标处中的点,使用$v = [v_x v_y]^T$来表示第二帧$[x y]^T$坐标处中的一点,此时上述问题就转变成:

在第二帧中求解出一个$v$,使得$I(u)$和$J(v)$相近.

我们用$d=[d_x d_y]^T$来记录中图片中$u$点处的速度(或者称其为$u$点处的光流),此时有:
$$v = u + d = [u_x+d_x u_y+d_y]^T$$

通过上面的式子,我们找到了一个求解第二帧对应点的方法,通过计算第一帧某点出的$d$(即,这点处的速度)即可找到第二帧中的该点的对应点。

$d$的求解问题请查看公式推导章节.

3. 公式推导过程

//TODO ???gray 与亮度的关系
我们假设亮度恒定(如上所述,我们是根据灰度值来寻找第二帧中移动后的点),我们假设t时刻某一个点(x y)在很短时间内移动后的坐标是$(x+d_x y+d_y)$。
定义如下的剩余函数(the residual function):
$$\xi(d) = \xi(d_x d_y) = \sum_{x=u_x-w_x}^{u_x+w_x}
\sum_{x=u_y-w_y}^{u_y+w_y}
(I(x,y) - J(x+d_x,y+d_y))^2$$

  • 这个式子实际上表达的是:以两个点为中心,窗口大小为$(2w_x+1)*(w_y+1)$范围内灰度值的差。
    d_residual_function

此时我们的任务就是寻找这个函数的最小值,也就是寻找最下的$d$。
//TODO

不难证明函数$\xi(d)$取最小点处的导数为零.

上式求导数有:
$$
\frac{\partial \xi(d)}{\partial d} =
-2 \sum_{x=u_x-w_x}^{u_x+w_x}
\sum_{x=u_y-w_y}^{u_y+w_y}
(I(x,y) - J(x+d_x,y+d_y))
[\frac {\partial J} {\partial x}
\frac {\partial J} {\partial y}]
$$
将$J(x+d_x,y+d_y)$泰勒展开的:
$$\frac{\partial \xi(d)}{\partial d} ≈
-2 \sum_{x=u_x-w_x}^{u_x+w_x}
\sum_{x=u_y-w_y}^{u_y+w_y}
(I(x,y) - J(x,y) - [\frac {\partial J} {\partial x}
\frac {\partial J} {\partial y}] d)
[\frac {\partial J} {\partial x} \frac {\partial J} {\partial y}]
$$

  • $I(x,y)和J(x,y)表示的同一坐标处的灰度值$,记${\delta I} = I(x,y) - J(x,y)$
  • $[\frac {\partial J} {\partial x} \frac {\partial J} {\partial y}]$则是图像的梯度。(记为这个矩阵的转置为$\nabla I$)
  • 又记$
    I_x(x,y)=
    {\frac {\partial I(x,y)}{\partial x}} =
    {\frac {I(x+1,y) - I(x-1, y)} 2}
    $,$
    I_y(x,y)=
    {\frac {\partial J(x,y)}{\partial x}} =
    {\frac {I(x,y+1) - I(x,y-1)} 2}
    $
    ($J_x(x,y)$和$J_y(x,y)$定义类似)

符号替换有:
$$
{\frac 1 2}\frac{\partial \xi(d)}{\partial d} ≈
\sum_{x=u_x-w_x}^{u_x+w_x}
\sum_{x=u_y-w_y}^{u_y+w_y}
({\nabla I^T}d - {\delta I}){\nabla I^T}
$$
又有:
$$
{\frac 1 2} [\frac{\partial \xi(d)}{\partial d}]^T ≈
\sum_{x=u_x-w_x}^{u_x+w_x}
\sum_{x=u_y-w_y}^{u_y+w_y}
(
[
\begin{matrix}
I^2_x & I_xI_y\
I_yI_y & I^2_y \
\end{matrix}
]d
-[\begin{matrix}
\delta I & I_x\
\delta I & I_y \
\end{matrix}]
)
$$
最后一次符号替换有:
$$
{\frac 1 2} [\frac{\partial \xi(d)}{\partial d}]^T

Gd - \overline b
$$
令偏导数为零有:
$$d = G^{-1} \overline b$$

d一定可逆

4. opencv代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

Mat mask = Mat::zeros(currentFrame.size(), currentFrame.type());

//当前与下一个灰度图像
cv::Mat currentgrayImg, nextgrayImg;
cvtColor(currentFrame, currentgrayImg, COLOR_BGRA2GRAY);
cvtColor(nextframe, nextgrayImg, COLOR_BGRA2GRAY);

// 得到“好的跟踪点”
vector<Point2f> p0, p1;
cv::goodFeaturesToTrack(currentgrayImg, p0, 10, 0.3,
7,Mat(), 7,false,0.04);

vector<uchar> status;
vector<float> err;

TermCriteria criteria = TermCriteria((TermCriteria::COUNT) + (TermCriteria::EPS), 10, 0.03);
calcOpticalFlowPyrLK(currentgrayImg, nextgrayImg, p0, p1, status, err, Size(15,15), 2, criteria);

// 可视化
Mat img;
vector<Scalar> colors;
RNG rng;
for(int i = 0; i < 100; i++)
{
int r = rng.uniform(0, 256);
int g = rng.uniform(0, 256);
int b = rng.uniform(0, 256);
colors.push_back(Scalar(r,g,b));
}

for(uint i = 0; i < p0.size(); i++) {
line(mask,p1[i], p0[i], colors[i], 2);
//,circle(currentFrame, p1[i], 5, colors[i], -1);
}

add(currentFrame, mask, img);
imshow("Frame", img);
cv::waitKey(100);

5. 参考