传统光流算法和代码

光流(optical flow)表示连续两帧的图片中目标像素的移动, 光流是二维矢量场, 表示目标像素从第一帧到第二帧的位移, 即对相邻两帧图片的目标像素的速度估计.

光流算法分为“稠密光流”和“稀疏光流”:

• 稠密光流: 对相邻两帧图像中的每个像素的速度进行估计, 也就是计算图像中每个像素的位移矢量.
• 稀疏光流: 只对特定像素点进行跟踪, 关注特定点的位移. 这种算法更加快速高效, 成本低于稠密跟踪.

光流算法分为传统方法和深度学习相关的方法, 下面会逐一介绍.

稀疏光流的传统算法 —— Lucas-Kanade算法

简称为LK光流法, 是两帧差分的光流估计算法, 1981年发表于“An Iterative Image Registration Technique with an Application to Stereo Vision”.
LK光流法基于以下三个假设:

1. 亮度恒定: 同一个点随着时间的变化, 亮度不会发生变化. 意味着像素的亮度不会随着帧的追踪而变化. (亮度恒定是所有光流法变种都需要满足的条件)
2. 时间持续性: 随着时间的变化, 位置不会发生剧烈变化. 这样像素才可以对位置求偏导. (所有光流法都需要满足该条件)
3. 空间一致性: 图片中相邻点具有相似的运动. (只对LK光流法的特定假设).

LK光流法的算法原理:
对一段相机视频, 图像是随时间变化的. $I(x,y,t)$表示t时刻位于$(x,y)$的像素的灰度值. 根据亮度恒定的条件, 有约束方程:

$$I(x,y,t) = I(x+\delta x, y+\delta y, t+\delta t)$$

根据泰勒公式推导得到(具体看参考[1]):

$$\frac{\partial I}{\partial x} \frac{\delta x}{\delta t} + \frac{\partial I}{\partial y} \frac{\delta y}{\delta t} + \frac{\partial I}{ \partial t} = 0$$

其中$\frac{\delta x}{\delta t}$和$\frac{\delta y}{\delta t}$分别是像素沿着x和y轴的导数. 沿x和y方向的速度分量记为u和v, 则上述式子简化为(其中$I_x=\frac{\partial I}{\partial x}$表示图像I的亮度在x轴方向的变化):

$$I_x u + I_y v + I_t = 0$$

根据假设3, 在大小为$m\times m =n$的小窗口内, 所有像素具有相同的光流u和v, 即:

$$I_x^{(i)} u + I_y^{(i)} v = - I_t^{(i)}, \ for \ i=1,2...n$$

现在有n个方程, 两个未知数u和v, 可以用最小二乘法求解, 最终得到光流(具体推倒看参考[1], 或者问GPT):

$$\begin{bmatrix} u \\ v \\ \end{bmatrix} = (A^TA)^{-1}A^T(b)$$

以上是LK算法的推导, LK的局限性如下:

1. 是稀疏光流算法, 只能检测关键点, 不能得到密度很高的光流向量.
2. 存在‘孔径问题(Aperture Problem)’, 即当我们从小窗口观察纹理时, 无法判断其真正的运动方向, 只能看到垂直于边缘方向的运动分量.
3. 假设过于理想, LK算法假设像素的亮度随时间不变, 相邻帧的像素运动很小. 当物体运动太快或者有遮挡, 或者光照变化, 假设会失效.
4. 对光照理想. 若图像在亮度上随时间变化, LK算法会出问题. 若物体运动快速, 需要结合金字塔(Pyramid)方法分层计算.

金字塔光流(Pyramidal Optical Flow)模型

LK算法在处理大位移时性能下降, 为了对此改进, 我们介绍金字塔光流模型. 金字塔模型会对原始图像进行L次上采样得到金字塔. 其中最底层是原图, 假设最底层是第L层, 最顶层是第0层.

基于LK算法的金字塔光流:

1. 首先构建金字塔. 原始图像在最底层, 上一层是下一层的下采样, 上一层的图像大小是下一层的一半.
2. 对最顶层的图像进行LK光流计算. 得到初始光流$(u_0,v_0)$, 因为最顶层的图像最小, 运动也被压缩, 原始图像的大位移被压缩为小位移.
3. 逐层向下传播, 假设现在计算了第$l-1$层的光流是$(u_{l-1},v_{l-1})$, 如何得到下一层也就是第$l$层的光流$(u_l,v_l)$. 首先对$l-1$层的光流上采样放大, 然后用这个光流对第$l$层的后一帧图像$I^2_{l-1}$进行反向转换(warping),warp后的图片和前一帧$I^1_{l-1}$计算LK光流,也就是对误差进行光流修正. 准确值=估计值+残差, 其中对warp图片计算LK就是计算的残差. 具体公式如下: $$(u_{l-1}, v_{l-1}) = LK(I^1_{l-1}, warp(I^2_{l-1}, u_l,v_l)) + 2*upsample(u_l,v_l)$$

注意, 下采样(downsample)会降低图片尺寸, 上采样(upsample)增加图片尺寸.

Farneback光流(Farneback Dense Optical Flow)

Farneback方法是一种稠密光流估计算法, 而LK算法是稀疏光流. 金字塔模型既可以用于LK算法, 也可以用于Farneback算法.

Farneback算法的核心思想:
用多项式（尤其是二次多项式）近似图像局部区域的像素强度分布，然后通过观察这些多项式在两帧图像间的变化来估计运动.

1. 局部建模: 将图像每个小区域(例如5*5 patch)用一个二次多项式建模, 近似为一个‘模糊的山丘’.
2. 多项式变化跟踪: 比如前一帧是一个 “山丘” 的形状，后一帧它的位置稍微偏移了。Farneback通过比较前后两帧这个“山丘”的位置和形状变化，来估算这个区域的运动。
3. 平滑约束: 使用高斯核等平滑操作让估计更稳健(防止噪声影响).
4. 也可以添加金字塔模型, 提高对大位移的处理.

Farneback算法所需要的假设条件:

• 灰度一致性brightness constancy, 对应点的亮度值在不同帧之间需要相同.(这是对光流的基本假设)
• 小位移假设(small motion), 相邻帧之间的物体移动不能太大, 否则会拟合失效, 可通过金字塔缓解.
• 光滑性假设(spatial coherence), 相邻像素的运动趋势需要相似.

franeback算法不能处理遮挡和断裂, 不适用于大角度旋转或者强烈形变. 对亮度变化和噪声敏感, 不如深度学习光流方法的精度高.

openCV中常用的光流代码

金字塔LK方法 — cv2.calcOpticalFlowPyrLK

金字塔LK方法, 用于追踪关键点, 是稀疏光流算法.
下面是案例代码, 首先在第一帧img1寻找特征点p0, 然后使用LK算法找这些特征点在img2中对应位置, 最终得到运动轨迹, 即$p0\rightarrow p1$的位移向量.

p0 = cv2.goodFeaturesToTrack(img1, maxCorners=100, qualityLevel=0.3, minDistance=7) 
# 寻找第一帧img1的特征点(角点)p0
p1, st, err = cv2.calcOpticalFlowPyrLK(img1, img2, p0, None,
                                       winSize=(15, 15), maxLevel=2,
                                       criteria=(cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03)) 
# 找到第一帧特征点p0 在第二帧的对应位置p1, p0->p1的位移向量就是光流

其中img1是第一帧图像(灰度图), img2是第二帧图像;
maxCorners=100: 最多检测 100 个特征点;
qualityLevel=0.3: 阈值因子，0.3 表示角点响应大于最大值的 30% 才保留;
minDistance=7: 每两个角点之间的最小距离(像素);
p0: img1上找到的要跟踪的点, (N×1×2)数组，形如[x,y];

winSize=(15, 15): 在每一层图像金字塔中使用的搜索窗口大小;
maxLevel=2: 使用的金字塔层数(原图 + 2层 = 共3层);
criteria: 终止条件：
    cv2.TERM_CRITERIA_EPS：当迭代变化小于 epsilon=0.03 时停止;
    cv2.TERM_CRITERIA_COUNT：最多迭代 10 次;

p1：在第二帧上的预测点(same shape as p0);
st：状态标志(1 表示成功跟踪，0 表示失败);
err：每个点的误差估计(越小越好);

Farneback方法 — cv2.calcOpticalFlowFarneback

flow = cv2.calcOpticalFlowFarneback(
    img1, img2, None,
    pyr_scale=0.5,     # 图像金字塔每层的缩放比例
    levels=3,          # 金字塔层数（越多越能捕捉大位移）
    winsize=15,        # 光流窗口大小（越大越平滑，越小越敏感）
    iterations=3,      # 每层迭代次数
    poly_n=5,          # 多项式展开邻域像素块大小（通常为5或7）
    poly_sigma=1.2,    # 多项式展开高斯平滑参数（通常1.1~1.5）
    flags=0            # 算法选项（默认0，可加OPTFLOW_FARNEBACK_GAUSSIAN）
)

输出的flow就是一个两通道的光流向量，实际上是每个点的像素位移值.

参考:

[1] 知乎: 光流法(optical flow methods)