LK光流法是图像亮度的运动信息描述，它基于三个假设：
- 运动物体的灰度在很短时间内保持不变
- 时间连续或者运动是小运动（可以泰勒展开）
- 空间一致，临近点有相似运动（保证对于同一个窗口所有点的偏移量都相等）
###原理
假设图像上一个像素点 (x,y) 的像素值为 \(I(x,y,t)\)
$$
I(x+dx,y+dy,t+dt) = I(x,y,t)
$$
$$
I(x+dx,y+dy,t+dt) \approx I(x,y,t) + \frac{\partial I}{\partial x} dx + \frac{\partial I}{\partial y} dy + \frac{\partial I}{\partial t} dt
$$
可以得到：
$$
- \frac{\partial I}{\partial t} dt = \frac{\partial I}{\partial x} dx + \frac{\partial I}{\partial y} dy
$$
即：
$$ - I_t = I_x u+I_y v
$$
其中 \(I_x, I_y, I_t \)分别为图像在 \((x,y)\) 处的梯度和关于 t 的导数，通常可以这样计算：
$$
I_x = \frac{I(x+1,y,t) - I(x-1,y,t)}{2} \
I_y = \frac{I(x,y+1,t) - I(x,y-1,t)}{2} \\
I_t = I(x,y,t) - I(x,y,t-1)
$$
对于一个像素点 \((x,y)\) 存在两个未知量 \((u,v)\) ，然而我们只有一个方程，这时候要利用第三个假设，考虑一个 w x w 大小的窗口，它们在两帧图像之间具有相同的运动：
$$
\left[I_{xk}, I_{yk} \right] \begin{bmatrix}u \\ v \end{bmatrix} = -I_{tk}, k=1,…,w^2
$$
令
$$
A = \begin{bmatrix} I_{x1} & I_{y1} \\ … \\ I_{xw^2} & I_{ yw^2}\end{bmatrix}, b = \begin{bmatrix} I_{t1} \\ …\\ I_{tw^2} \end{bmatrix}
$$
则
$$
A \begin{bmatrix}u \\ v \end{bmatrix} = -b
$$
这是一个超正定方程，我们可以用最小二乘法求解，当 \(A^TA\) 可逆时：
$$
\begin{bmatrix}u \\ v \end{bmatrix}^* = -(A^TA)^{-1} A^Tb
$$
求解 LK
为方便后面内容展开，我们这里重新定义一些变量，相邻的两帧图像为 \(I , J\) ，对于 \(I\) 中的像素点 \(u = [u_x, u_y]^T\) ，我们要在 \(J\) 中找出对应的像素点 \(v = u+d = [u_x+d_x, u_y+d_y] \) 使得它们之间灰度差值最小。在上一小节中光流定义为速度，而这里我们用位移代替速度，把 \(d = [d_x,d_y]^T\) 称为在 u 处的光流，引入一个具有相同光流的邻域，大小用 \(w_x, w_y\) 两个参数表示，则求解 d 转化为使下述目标函数取最小值的优化问题：
$$
\varepsilon(d) = \varepsilon(d_x,d_y) = \sum_{x=u_x-w_x}^{x=u_x+w_x} \sum_{y=u_y-w_y}^{y=u_y+w_y} (I(x,y)-J(x+d_x,y+d_y))^2
$$
$$
\frac{\partial{\varepsilon(d)}}{\partial d} = -2 \sum_{x=u_x-w_x}^{x=u_x+w_x} \sum_{y=u_y-w_y}^{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}]
$$
一阶泰勒展开：
$$
\frac{\partial{\varepsilon(d)}}{\partial d} \approx -2 \sum_{x=u_x-w_x}^{x=u_x+w_x} \sum_{y=u_y-w_y}^{y=u_y+w_y} (I(x,y)-J(x,y) - [\frac{\partial J}{\partial x} \ \ \frac{\partial J}{\partial y}] \begin{bmatrix}d_x \\d_y \end{bmatrix}) [\frac{\partial J}{\partial x} \ \ \frac{\partial J}{\partial y}]
$$
由于 \(d = [d_x \ d_y]^T\) 足够小，所以可以把 \([\frac{\partial J}{\partial x} \ \ \frac{\partial J}{\partial y}]\) 替换为 \([\frac{\partial I}{\partial x} \ \ \frac{\partial I}{\partial y}]\)，同时定义
$$
\delta I(x,y) = I(x,y) - J(x,y) \
\nabla I = [\frac{\partial I}{\partial x} \ \ \frac{\partial I}{\partial y}]^T
$$
则公式转化为：
$$
\frac{1}{2} \frac{\partial{\varepsilon(d)}}{\partial d} = \sum_{x=u_x-w_x}^{x=u_x+w_x} \sum_{y=u_y-w_y}^{y=u_y+w_y} (\nabla I^Td-\delta I)\nabla I^T
$$
$$
\frac{1}{2} [\frac{\partial{\varepsilon(d)}}{\partial d}]^T = \sum_{x=u_x-w_x}^{x=u_x+w_x} \sum_{y=u_y-w_y}^{y=u_y+w_y}(\nabla I^Td-\delta I)\nabla I = \sum_{x=u_x-w_x}^{x=u_x+w_x} \sum_{y=u_y-w_y}^{y=u_y+w_y} \begin{bmatrix}I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{bmatrix}d - \sum_{x=u_x-w_x}^{x=u_x+w_x} \sum_{y=u_y-w_y}^{y=u_y+w_y} \begin{bmatrix} \delta I \cdot I_x \\ \delta I \cdot I_y \end{bmatrix}
$$
令
$$
G = \sum_{x=u_x-w_x}^{x=u_x+w_x} \sum_{y=u_y-w_y}^{y=u_y+w_y} \begin{bmatrix}I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{bmatrix} \
b = \sum_{x=u_x-w_x}^{x=u_x+w_x} \sum_{y=u_y-w_y}^{y=u_y+w_y} \begin{bmatrix} \delta I \cdot I_x \\ \delta I \cdot I_y \end{bmatrix}
$$
则
$$
d = G^{-1} b
$$
金字塔LK
在实际场景中，物体相邻帧间运动足够小是很难满足的，如果像素点运动较大时，那么一阶泰勒展开就不太精确。所以金字塔LK就被提出了，对于分辨率为800x800的图像，若两帧间某个像素点的运动为40x40，当把图像分辨率变为400x400时，运动变为20x20……以此类推总能在图像像素降低到一定程度时原本较大的像素间的运动变得足够小，从而可以用泰勒展开。
以原始图像为第一层，宽高缩小\(2^L\) 倍的图像作为第L层，多层图像就组成了一个金字塔
需要用低通滤波器来平滑图像来防止降采样后出现锯齿现象。当原始图像尺寸无法保证严格满足整除 \(2^L\)时，可以对原始图像调整尺寸满足整除。
总体流程：
假设金字塔有 \(L_0, L_1,…,L_{m-1}, L_m\) 共 m+1 层图像，首先计算最顶层\(L_m\)的光流，然后把它作为第\(L_m-1\) 层光流的初始值，再计算该层的光流，之后再将其作为第\(L_m-2\)层光流大小的初始值…….以此类推直至 第 \(L_0\) 层得到最终光流结果。
以第 \(L_{k+1}\) 层到第\(L_{k}\) 层为例说明计算流程，得到第\(L_{k}\) 层光流大小的初始值
$$
g^{L_k} = 2(g^{L_{k+1}} + d^{L_{k+1}})= [g_x^{L_k} \ g_y^{L_k}]^T
$$
我们要找使得下列函数最小的 \(d^{L_k} = [d_x^{L_k} \ d_y^{L_k}]^T\)
$$
\varepsilon(d^{L_k}) = \varepsilon(d_x^{L_k},d_y^{L_k}) = \sum_{x=u_x^{L_k}-w_x}^{x=u_x^{L_k}+w_x} \sum_{y=u_y^{L_k}-w_y}^{y=u_y^{L_k}+w_y} (I(x,y)-J(x+g_x^{L_k}+d_x^{L_k},y+g_y^{L_k}+d_y^{L_k}))^2
$$
$$
\frac{\partial{\varepsilon(d^{L_k})}}{\partial d} = -2 \sum_{x=u_x^{L_k}-w_x}^{x=u_x^{L_k}+w_x} \sum_{y=u_y^{L_k}-w_y}^{y=u_y^{L_k}+w_y} (I(x,y)-J(x+g_x^{L_k},y+g_y^{L_y})) - [\frac{\partial J}{\partial x} \ \ \frac{\partial J}{\partial y}]\begin{bmatrix}d_x^{L_k} \\d_y^{L_k} \end{bmatrix})
$$
$$
\delta I(x,y) = I(x,y) - J(x+g_x^{L_k},y+g_y^{L_k}) \
$$
由于 \(d = [d_x \ d_y]^T\) 足够小，所以可以把 \([\frac{\partial J}{\partial x} \ \ \frac{\partial J}{\partial y}]\) 替换为 \([\frac{\partial I}{\partial x} \ \ \frac{\partial I}{\partial y}]\)
$$
\nabla I = [\frac{\partial I}{\partial x} \ \ \frac{\partial I}{\partial y}]^T
$$
定义
$$
G = \sum_{x=u_x^{L_k}-w_x}^{x=u_x^{L_k}+w_x} \sum_{y=u_y^{L_k}-w_y}^{y=u_y^{L_k}+w_y} \begin{bmatrix}I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{bmatrix} \
b = \sum_{x=u_x^{L_k}-w_x}^{x=u_x^{L_k}+w_x} \sum_{y=u_y^{L_k}-w_y}^{y=u_y^{L_k}+w_y} \begin{bmatrix} \delta I \cdot I_x \\ \delta I \cdot I_y \end{bmatrix}
$$
则
$$
d^{L_k} = G^{-1} b
$$
下一层的光流初始值为：
$$
g^{L_{k-1}} = 2(g^{L_k}+d^{L_k})
$$
在每层图像上都要根据光流初始值先平移窗口，然后计算残差位移d 就非常小，就可以泰勒展开从而利用标准 LK 算法求解。最高层 \(g^{L_m} = [ 0 \ 0 ]^T\)，那么原始图像的光流为:
$$
d = g^0 + d^0
$$
可以看出最终的光流就是所有层的分段的叠加:
$$
d = \sum_{k=0}^{m} 2^k d^{L_k}
$$
我们可以看到使用金字塔计算光流的好处。对于每一层的光流都会保持很小，但是它们可以累积从而得到原始层的光流。
特征点的选择
光流法可以完成对所有像素点的光流计算，但是时间开销很大，所以我们需要选取一些特征点，我们可以选取使得矩阵 G 可逆的特征点，参见 \<\
参考文献：
