Farneback 光流算法详解与 calcOpticalFlowFarneback 源码分析

光流基础概念

真实的三维空间中，描述物体运动状态的物理概念是运动场。三维空间中的每一个点，经过某段时间的运动之后会到达一个新的位置，而这个位移过程可以用运动场来描述。

而在计算机视觉的空间中，计算机所接收到的信号往往是二维图片信息。由于缺少了一个维度的信息，所以其不再适用以运动场描述。光流场（optical flow）就是用于描述三维空间中的运动物体表现到二维图像中，所反映出的像素点的运动向量场。

• 当描述部分像素时，称为：稀疏光流
• 当描述全部像素时，称为：稠密光流

光流法是利用图像序列中的像素在时间域上的变化、相邻帧之间的相关性来找到的上一帧跟当前帧间存在的对应关系，计算出相邻帧之间物体的运动信息的一种方法。光流法理解的关键点有：

• 核心问题：同一个空间中的点，在下一帧即将出现的位置
• 重要假设：光流的变化（向量场）几乎是光滑
• 角点处的光流能够通过角点的邻域完全确定下来，因此角点处的运动信息最为可靠；其次是边界的信息

光流法有着各种各样的分支，本文介绍的则是一种被广泛使用的经典稠密光流算法：Farneback 光流算法

---

Farneback 光流算法

• 作者：Gunnar Farneback
• 参考资料
  
  • Two-Frame Motion Estimation Based on Polynomial Expansion，提纲阅读
  • Polynomial Expansion for Orientation and Motion Estimation，作者博士毕业论文，深读（以下简称：博士论文）
• 源码：https://searchcode.com/file/30099587/opencv_source/src/cv/cvoptflowgf.cpp

OpenCV 主函数源码：

void calcOpticalFlowFarneback( const Mat& prev0, const Mat& next0, Mat& flow0, 
    double pyr_scale, int levels, 
    int winsize, int iterations, int poly_n, double poly_sigma, int flags )
{
    ……
    for( k = 0, scale = 1; k < levels; k++ )
    {
        scale *= pyr_scale;
        if( prev0.cols*scale < min_size || prev0.rows*scale < min_size )
            break;
    }

    levels = k;

    for( k = levels; k >= 0; k-- )
    {
        for( i = 0, scale = 1; i < k; i++ )
            scale *= pyr_scale;
        ……
        Mat R[2], I, M;
        for( i = 0; i < 2; i++ )
        {
            img[i]->convertTo(fimg, CV_32F);
            GaussianBlur(fimg, fimg, Size(smooth_sz, smooth_sz), sigma, sigma);
            resize( fimg, I, Size(width, height), CV_INTER_LINEAR );
            FarnebackPolyExp( I, R[i], poly_n, poly_sigma );
        }

        FarnebackUpdateMatrices( R[0], R[1], flow, M, 0, flow.rows );

        for( i = 0; i < iterations; i++ )
        {
            if( flags & OPTFLOW_FARNEBACK_GAUSSIAN )
                FarnebackUpdateFlow_GaussianBlur( R[0], R[1], flow, M, winsize, i < iterations - 1 );
            else
                FarnebackUpdateFlow_Blur( R[0], R[1], flow, M, winsize, i < iterations - 1 );
        }

        prevFlow = flow;
    }
}

源代码中为了解决孔径问题，主函数中引入了图像金字塔。

孔径问题（Aperture Problem）：
- http://blog.csdn.net/hankai1024/article/details/23433157
- 形象理解：从小孔中观察一块移动的黑色幕布观察不到任何变化。但实际情况是幕布一直在移动中
- 解决方案：从不同尺度（图像金字塔）上对图像进行观察，由高到低逐层利用上一层已求得信息来计算下一层信息

主函数 calcOpticalFlowFarneback 中需要的变量参数包括：
1. _prev0：输入前一帧图像
2. _next0：输入后一帧图像
3. _flow0：输出的光流
4. pyr_scale：金字塔上下两层之间的尺度关系
5. levels：金字塔层数
6. winsize：均值窗口大小，越大越能 denoise 并且能够检测快速移动目标，但会引起模糊运动区域
7. iterations：迭代次数
8. poly_n：像素邻域范围大小，一般为 5、7 等
9. poly_sigma：高斯标准差，一般为 1～1.5（函数处理中需要高斯分布权重）
10. flags：计算方法，包括 OPTFLOW_USE_INITIAL_FLOW 和 OPTFLOW_FARNEBACK_GAUSSIAN

实际的 Farneback 算法在每一层金字塔上的处理过程为：

Mat R[2], I, M;
for( i = 0; i < 2; i++ )
{
    img[i]->convertTo(fimg, CV_32F);
    GaussianBlur(fimg, fimg, Size(smooth_sz, smooth_sz), sigma, sigma);
    resize( fimg, I, Size(width, height), CV_INTER_LINEAR );
    FarnebackPolyExp( I, R[i], poly_n, poly_sigma );
}

FarnebackUpdateMatrices( R[0], R[1], flow, M, 0, flow.rows );

for( i = 0; i < iterations; i++ )
{
    if( flags & OPTFLOW_FARNEBACK_GAUSSIAN )
        FarnebackUpdateFlow_GaussianBlur( R[0], R[1], flow, M, winsize, i < iterations - 1 );
    else
        FarnebackUpdateFlow_Blur( R[0], R[1], flow, M, winsize, i < iterations - 1 );
}    

输入的图像默认为灰度图片，默认只有亮度值

图像输入与输出时进行的高斯模糊操作（Gaussian Blur）操作都是使得光流场结果平滑，从而满足假设：光流的变化几乎是光滑的

所以需要关注的子函数有 4 个：
1. FarnebackPolyExp
2. FarnebackUpdateMatrices
3. FarnebackUpdateFlow_GaussianBlur
4. FarnebackUpdateFlow_Blur

OpenCV 子函数 FarnebackPolyExp：

static void FarnebackPolyExp( const Mat& src, Mat& dst, int n, double sigma )

理论基础

图像建模

将图像视为二维信号的函数（输出图像是灰度图像），因变量是二维坐标位置 $\mathbf{x}=\left( \begin{matrix} x & y \end{matrix} \right)^T$，并且利用二次多项式对于图像进行近似建模的话，会得到:

$$f(\mathbf{x})\sim \mathbf{x^TAx+b^Tx+}c$$
其中， $\mathbf{A}$ 是一个2 $\times$2的矩阵， $\mathbf{b}$是一个2 $\times$1的矩阵

因此系数化之后，以上公式等号右侧可以写为：

(xy)T(r4r6/2r6/2r5)(xy)+(r2r3)T(xy)+r1=r1+r2x+r3y+r4x2+r5y2+r6xy