CN106600629B - 一种面向大尺度运动的光流估计方法 - Google Patents

Abstract

本发明公开了一种面向大尺度运动的光流估计方法，其目的是计算包含大位移和大尺度变化场景对象的图像之间的稠密像素对应关系。首先，对输入的图像对进行预处理；然后，构建光流控制点集合，并在此基础上计算控制光流向量和控制点状态矩阵，以及构建初始光流向量和光流偏移向量；最后，在离散优化框架下实现光流向量的迭代优化。采用本发明方法可以有效处理光流估计中场景对象大位移和大尺度变化问题，显著提高大尺度运动光流估计的精度，且与已有方法相比本发明方法计算效率有显著提升。

Description

一种面向大尺度运动的光流估计方法

技术领域

本发明属于图像信息处理技术领域，特指一种面向大尺度运动的光流估计方法。

背景技术

光流估计指的是计算两幅或多幅图像之间的稠密像素对应关系的过程，光流估计是图像信息处理研究领域的基础性问题。光流估计在运动估计、视频压缩、图像配准、视图插值、三维重建等许多计算机视觉任务中都有重要应用。在过去一、二十年里光流估计方法得以广泛研究，光流估计的性能也取得了较大的进展。

目前，比较流行的光流估计方法是经典变形光流估计方法(T.Brox,A.Bruhn,N.Papenberg,and J.Weickert,“High accuracy optical flow estimation based onatheory for warping,”in ECCV,2004)或者基于该方法的改进版本。这类方法根据估计的光流计算第一幅输入图像的变形图像，然后计算该变形图像与第二幅输入图像之间的亮度、梯度差异，并通过数值方法对估计的光流进行迭代优化。对于图像中场景对象移动范围较小的情况，上述方法可以获得较高精度的光流估计结果，其缺点是不能处理包含大尺度运动的情形。

场景对象的大尺度运动可能导致场景对象在最终成像平面上的位置、大小发生较大变化，这种情况下的光流估计与传统的变化幅度较小的光流估计相比具有更大的挑战。针对大尺度运动的光流估计方法近年来开始引起学者关注。其中，大位移光流估计方法(T.Brox and J.Malik,“Large displacement optical flow:Descriptor matching invariational motionestimation,”IEEE Transactions on Pattern Analysis andMachine Intelligence,33(3):500-513,2011)利用较大的场景对象的运动来初始化整个图像的运动，然后通过多步迭代逐步求取更加细小的场景结构光流变化。上述方法的缺点是在小尺度场景对象的运动幅度超过较大的场景对象的运动幅度的情况下，该方法不能获得正确的光流估计结果。SIFT流估计方法(C.Liu,J.Yuen,and A.Torralba,“SIFT ow:dense correspondence across scenes and its applications,”IEEE Transactions onPattern Analysis and Machine Intelligence,33(5):978-994,2011)计算图像中每个像素的SIFT特征，利用特征匹配和迭代优化方法计算图像间像素的对应关系。该方法可一定程度上提高大位移光流估计性能，但其精度局限于像素级，不能实现亚像素级光流估计。运动细节保持光流估计方法(L.Xu,J.Jia,and Y.Matsushita,“Motion detail preservingoptical flow estimation,”IEEE Transactionson Pattern Analysis and MachineIntelligence,34(9):1744-1757,2012)可以处理包含复杂大尺度运动的光流估计，但该方法的运行效率远低于经典变形光流估计方法。

发明内容

为了解决现有技术的不足，本发明的目的在于提出一种面向大尺度运动的光流估计方法。本发明利用图像稀疏特征点匹配关系，构建图像初始光流和控制光流约束关系，在此基础上对光流进行迭代优化，实现包含大尺度运动的光流的精确、高效计算。

一种面向大尺度运动的光流估计方法，包括以下步骤：

第一步，图像预处理

若输入的图像为彩色图像，则先转化为灰度图。记输入的两幅图像(均为灰度图)分别为和图像宽度为W个像素、高度为H个像素，任意像素位置p＝(x,y)(x＝{1,2,…,W}，y＝{1,2,…,H})上图像灰度值记为I₁(p)＝I₁(x,y)、I₂(p)＝I₂(x,y)。其中，x为像素沿水平方向的坐标，y为像素沿垂直方向的坐标。

分别提取图像和中的ASIFT特征点(Jean-Michel Morel and GuoshenYu.ASIFT:A new framework for fully affine invariant image comparison)，并利用SIFT特征描述子(D.G Lowe.Distinctive image features from scale-invariant keypoints.International Journal of Computer Vision,60(2):91–110,2004.)对所提取的特征点进行描述。

然后，利用最近邻搜索方法计算中各ASIFT特征点在中的对应特征点，并利用特征点描述子距离比测试剔除误匹配(参见D.G Lowe.Distinctive image features fromscale-invariant key points.International Journal of Computer Vision,60(2):91–110,2004.)，其中距离比阈值取值范围是0.6～0.8。记图像中最终匹配特征点集合为

最后，计算在图像上的特征覆盖率c₁，即中所有特征点的最小外接矩形面积与图像的面积之比。

第二步，构建光流控制点集合步骤如下：

1)参数初始化：令k＝1,

2)利用k-均值方法，根据特征点的图像坐标对中的特征点进行聚类，聚类数量设定为k；

3)将各聚类中距离相应聚类中心最近的匹配特征点加入光流控制点集合并将该匹配特征点从中删除；

4)计算在图像上的特征覆盖率c₂。如果c₂＜0.8c₁，则令k←k+1，转至步骤2)；否则，结束并输出控制点集合

第三步，计算控制光流向量和控制点状态矩阵

3.1初始化控制点约束光流和状态标识。

将所有像素水平、垂直方向约束光流置为零，即令u_gt(p)＝u_gt(x,y)＝0，v_gt(p)＝v_gt(x,y)＝0(以及)，并将所有像素位置上的控制点状态设置为没有控制约束光流状态，即令ρ(p)＝ρ(x,y)＝0(以及)。其中，u_gt(p)和v_gt(p)分别表示像素位置p处控制光流向量的水平分量和垂直分量，ρ(p)表示像素位置p处的控制点标识(0表示无控制点，1表示有控制点)。

3.2根据控制点集合更新控制点约束光流和状态标识。

对于控制点集合中的每一个控制点，按照如下方式确定相应的约束光流和状态标识：

首先，根据控制点在图像中的坐标分别在中确定距离匹配特征点最近的像素位置，记为p_s＝(x_s,y_s)、p_t＝(x_t,y_t)。其中，x_s和x_t均为不大于W的正整数，y_s和y_t均为不大于H的正整数。

然后，按照如下公式更新u_gt、v_gt、ρ：

3.3构建控制光流向量U_gt、V_gt以及控制点状态矩阵G

其中，U_gt、V_gt为WH维列向量，G为WH×WH维对角方阵。U_gt为由所有像素位置上水平方向的控制光流u_gt按照从左至右、从上至下的顺序排列而成的光流向量，V_gt为由所有像素位置上垂直方向的控制光流v_gt按照从左至右、从上至下的顺序排列而成的光流向量，对角矩阵G的对角线元素由所有像素位置上的控制光流状态ρ按照从左至右、从上至下的顺序排列而成。

第四步，构建初始光流向量和光流偏移向量

4.1光流和偏移初始化。

根据控制点约束光流u_gt、v_gt(ρ≠0)，利用三线性插值方法计算图像各像素位置上的初始光流w(p)＝(u(p),v(p))。其中，w(p)表示像素p处的光流向量，u(p)和v(p)分别表示像素位置p处光流向量的水平分量和垂直分量。并将各像素位置上的初始光流偏移值置为零，即du(p)＝0、du(p)＝0。

4.2构建初始光流向量U和V

其中，U为由所有像素水平方向光流值按照从左至右、从上至下的顺序排列构建而成的光流向量，V为由所有像素垂直方向光流值按照从左至右、从上至下的顺序排列构建而成的光流向量。

4.3构建初始光流偏移向量dU和dV

dU为由所有像素水平方向光流偏移值按照从左至右、从上至下的顺序排列构建而成的光流向量，dV为由所有像素垂直方向光流偏移值按照从左至右、从上至下的顺序排列构建而成的光流向量。

第五步，光流向量迭代优化

光流向量迭代优化计算过程中，首先初始化迭代次数：N_out＝0，然后按照下述过程迭代优化光流向量。

5.1按照如下公式更新光流向量、光流偏移向量：

5.2计算最佳光流偏移向量

计算最佳光流偏移向量的过程中，首先初始化迭代：N_in＝0，然后按照下述过程迭代更新最佳光流偏移向量。

5.2.1计算权重矩阵

1)构建梯度图像以及偏移亮度差异图像图像各像素位置p的取值按照如下公式计算：

2)构建WH×WH维对角方阵I_x、I_y、I_z。

I_x、I_y、I_z对角线上各元素分别由图像中各像素按照从左至右、从上至下的顺序填充而成。

3)构建权重矩阵T_f、T_g、T_h。上述三个矩阵均为WH×WH维对角方阵。对角线上第i个元素的取值由如下公式确定：

其中，符号[·]_i,i表示矩阵第i行、第i列元素，符号[·]_i表示向量第i个元素，D_x和D_y分别表示计算图像水平、垂直方向导数的WH×WH维矩阵算子，函数ψ的定义为：常数ε的取值范围是0.001～0.002。

5.2.2通过求解如下线性方程组，计算dU和dV

其中，

参数α的取值范围是0.01～0.03，参数β的取值范围是10～15。

5.2.3更新迭代次数：N_in←N_in+1。如果N_in＜2，则更新dU和dV，转到步骤5.2.1；否则，转到步骤5.3。

5.3更新迭代次数：N_out←N_out+1。如果N_out＜8，则转到步骤5.1；否则，输出最终光流：

本发明相比以往技术方法，具有以下明显优点：

与经典的光流估计方法相比，本发明方法可以有效处理光流估计中场景对象大位移和大尺度变化问题，获得高精度的光流估计结果。与已有面向大位移场景对象的光流估计方法相比，本发明方法在计算效率上有显著提升。

附图说明

图1为本发明一种面向大尺度运动的光流估计方法的流程图；

图2为本发明光流控制点集合生成过程的流程图；

图3为本发明光流向量迭代优化计算过程的流程图；

图4为本发明计算最佳光流偏移向量过程的流程图；

图5(a)至图5(d)为四组包含大尺度运动的测试图像；

图6为本发明方法与已有方法光流估计精度的比较；

图7为本发明方法与已有方法光流估计所需时间的比较。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

参照图1，为本发明一种面向大尺度运动的光流估计方法的流程图，其目的是计算包含大位移和大尺度变化场景对象的图像之间的稠密像素对应关系。首先，对输入的图像对进行预处理；然后，构建光流控制点集合，并在此基础上计算控制光流向量和控制点状态矩阵，以及构建初始光流向量和光流偏移向量；最后，在离散优化框架下实现光流向量的迭代优化。采用本发明方法可以有效处理光流估计中场景对象大位移和大尺度变化问题，显著提高大尺度运动光流估计的精度，且与已有方法相比本发明方法计算效率有显著提升。具体地，其通过以下步骤实现：

第一步，图像预处理

若输入的图像为彩色图像，则先转化为灰度图。记输入的两幅图像(均为灰度图)分别为和图像宽度为W个像素、高度为H个像素，任意像素位置p＝(x,y)(x＝{1,2,…,W}，y＝{1,2,…,H})上图像灰度值记为I₁(p)＝I₁(x,y)、I₂(p)＝I₂(x,y)。其中，x为像素沿水平方向的坐标，y为像素沿垂直方向的坐标。

分别提取图像和中的ASIFT特征点(Jean-Michel Morel and GuoshenYu.ASIFT:A new framework for fully affine invariant image comparison)，并利用SIFT特征描述子(D.G Lowe.Distinctive image features from scale-invariant keypoints.International Journal of Computer Vision,60(2):91–110,2004.)对所提取的特征点进行描述。

然后，利用最近邻搜索方法计算中各ASIFT特征点在中的对应特征点，并利用特征点描述子距离比测试剔除误匹配，距离比阈值取值范围是0.6～0.8。记图像中最终匹配特征点集合为

最后，计算在图像上的特征覆盖率c₁，即中所有特征点的最小外接矩形面积与图像的面积之比。

第二步，构建光流控制点集合步骤如下(如图2所示)：

1)参数初始化：令k＝1,

2)利用k-均值方法，根据特征点的图像坐标对中的特征点进行聚类，聚类数量设定为k；

3)将各聚类中距离相应聚类中心最近的匹配特征点加入光流控制点集合并将该匹配特征点从中删除；

4)计算在图像上的特征覆盖率c₂。如果c₂＜0.8c₁，则令k←k+1，转至步骤2)；否则，结束并输出控制点集合

第三步，计算控制光流向量和控制点状态矩阵

3.1初始化控制点约束光流和状态标识。

将所有像素水平、垂直方向约束光流置为零，即令u_gt(p)＝u_gt(x,y)＝0，v_gt(p)＝v_gt(x,y)＝0(以及)，并将所有像素位置上的控制点状态设置为没有控制约束光流状态，即令ρ(p)＝ρ(x,y)＝0(以及)。其中，u_gt(p)和v_gt(p)分别表示像素位置p处控制光流向量的水平分量和垂直分量，ρ(p)表示像素位置p处的控制点标识(0表示无控制点，1表示有控制点)。

3.2根据控制点集合更新控制点约束光流和状态标识。

对于控制点集合中的每一个控制点，按照如下方式确定相应的约束光流和状态标识：

首先，根据控制点在图像中的坐标分别在中确定距离匹配特征点最近的像素位置，记为p_s＝(x_s,y_s)、p_t＝(x_t,y_t)。其中，x_s和x_t均为不大于W的正整数，y_s和y_t均为不大于H的正整数。

然后，按照如下公式更新u_gt、v_gt、ρ：

3.3构建控制光流向量U_gt、V_gt以及控制点状态矩阵G

其中，U_gt、V_gt为WH维列向量，G为WH×WH维对角方阵。U_gt为由所有像素位置上水平方向的控制光流u_gt按照从左至右、从上至下的顺序排列而成的光流向量，V_gt为由所有像素位置上垂直方向的控制光流v_gt按照从左至右、从上至下的顺序排列而成的光流向量，对角矩阵G的对角线元素由所有像素位置上的控制光流状态ρ按照从左至右、从上至下的顺序排列而成。

第四步，构建初始光流向量和光流偏移向量

4.1光流和偏移初始化。

根据控制点约束光流u_gt、v_gt(ρ≠0)，利用三线性插值方法计算图像各像素位置上的初始光流w(p)＝(u(p),v(p))。其中，w(p)表示像素p处的光流向量，u(p)和v(p)分别表示像素位置p处光流向量的水平分量和垂直分量。并将各像素位置上的初始光流偏移值置为零，即du(p)＝0、du(p)＝0。

4.2构建初始光流向量U和V

其中，U为由所有像素水平方向光流值按照从左至右、从上至下的顺序排列构建而成的光流向量，V为由所有像素垂直方向光流值按照从左至右、从上至下的顺序排列构建而成的光流向量。

4.3构建初始光流偏移向量dU和dV

dU为由所有像素水平方向光流偏移值按照从左至右、从上至下的顺序排列构建而成的光流向量，dV为由所有像素垂直方向光流偏移值按照从左至右、从上至下的顺序排列构建而成的光流向量。

第五步，光流向量迭代优化

光流向量迭代优化计算过程如图3所示，首先初始化迭代次数：N_out＝0，然后按照下述过程迭代优化光流向量。

5.1按照如下公式更新光流向量、光流偏移向量：

5.2计算最佳光流偏移向量

该过程如图4所示。计算最佳光流偏移向量的过程中，首先初始化迭代：N_in＝0，然后按照下述过程迭代更新最佳光流偏移向量。

5.2.1计算权重矩阵

1)构建梯度图像以及偏移亮度差异图像图像各像素位置p的取值按照如下公式计算：

2)构建WH×WH维对角方阵I_x、I_y、I_z。I_x、I_y、I_z对角线上各元素分别由图像中各像素按照从左至右、从上至下的顺序填充而成，即对角线上第i个元素的取值由如下公式确定：

其中，符号[·]_i,j表示矩阵第i行、第j列元素。

3)构建权重矩阵T_f、T_g、T_h。上述三个矩阵均为WH×WH维对角方阵。对角线上第i个元素的取值由如下公式确定：

其中，符号[·]_i,i表示矩阵第i行、第i列元素，符号[·]_i表示向量第i个元素，D_x和D_y分别表示计算图像水平、垂直方向导数的WH×WH维矩阵算子，函数ψ的定义为：在此常数ε取0.001。

5.2.2通过求解如下线性方程组，计算dU和dV

其中，

参数α＝0.02，β＝12。

5.2.3更新迭代次数：N_in←N_in+1。如果N_in＜2，则更新dU和dV，转到步骤5.2.1；否则，转到步骤5.3。

5.3更新迭代次数：N_out←N_out+1。如果N_out＜8，则转到步骤5.1；否则，输出最终光流：

图5展示四组有代表性的、包含大尺度运动的测试图像，图像分辨率均为800×600像素。其中，图5(a)和图5(b)是由软件合成的图像，因此实际光流是已知的，图5(c)和图5(d)是真实图像。图5(a)和图5(c)是典型的大尺度远近运动情形，图5(b)和图5(d)是典型的大尺度平移运动情形。

图6是本发明方法与已有方法在图5(a)和图5(b)所示图像上的光流估计精度的比较。精度评估考察平均角度误差和平均端点误差两个指标。其中，平均角度误差指的是各像素到其估计对应像素形成的向量与该像素到其真实对应像素形成的向量之间的夹角的平均值，平均端点误差指的是各像素的估计对应像素与真实对应像素之间的距离的平均值。由于本发明方法引入了稀疏特征控制光流，因此可以有效处理光流估计中场景对象大位移和大尺度变化问题，显著提高大尺度运动光流估计精度。

本发明方法与已有方法所需计算时间的比较如图7所示。注意，经典变形光流估计方法所需时间与本发明方法相差较小，但从图6可以看出其估计精度远低于本发明方法，因此图7中没有展示经典变形光流估计方法所需时间。由于本发明采用高效的离散迭代优化方法，与具有相当精度的运动细节保持光流估计方法相比，本发明方法的计算效率得以显著提升。

以上包含了本发明优选实施例的说明，这是为了详细说明本发明的技术特征，并不是想要将发明内容限制在实施例所描述的具体形式中，依据本发明内容主旨进行的其他修改和变型也受本专利保护。本发明内容的主旨是由权利要求书所界定，而非由实施例的具体描述所界定。

Claims (1)

1.一种面向大尺度运动的光流估计方法，其特征在于，包括以下步骤：

第一步，图像预处理

若输入的图像为彩色图像，则先转化为灰度图；记输入的两幅图像分别为和图像宽度为W个像素、高度为H个像素，中任意像素位置p＝(x,y)上图像灰度值记为I₁(p)＝I₁(x,y)、I₂(p)＝I₂(x,y)；其中，x＝{1,2,…,W}，y＝{1,2,…,H}，x为像素沿水平方向的坐标，y为像素沿垂直方向的坐标；

分别提取图像和中的ASIFT特征点对所提取的特征点进行描述；

然后，利用最近邻搜索方法计算中各ASIFT特征点在中的对应特征点，并利用特征点描述子距离比测试剔除误匹配，其中距离比阈值取值范围是0.6～0.8，记图像中最终匹配特征点集合为

最后，计算在图像上的特征覆盖率c₁，即中所有特征点的最小外接矩形面积与图像的面积之比；

第二步，构建光流控制点集合C，步骤如下：

1)参数初始化：令k＝1,

2)利用k-均值方法，根据特征点的图像坐标对中的特征点进行聚类，聚类数量设定为k；

3)将各聚类中距离相应聚类中心最近的匹配特征点加入光流控制点集合C，并将该匹配特征点从中删除；

4)计算C在图像上的特征覆盖率c₂；如果c₂＜0.8c₁，则令k←k+1，转至步骤2)；否则，结束并输出控制点集合C；

第三步，计算控制光流向量和控制点状态矩阵

3.1初始化控制点约束光流和状态标识；

将所有像素水平、垂直方向约束光流置为零，即令u_gt(p)＝u_gt(x,y)＝0，v_gt(p)＝v_gt(x,y)＝0，其中以及并将所有像素位置上的控制点状态设置为没有控制约束光流状态，即令ρ(p)＝ρ(x,y)＝0，其中以及u_gt(p)和v_gt(p)分别表示像素位置p处控制光流向量的水平分量和垂直分量，ρ(p)表示像素位置p处的控制点标识，0表示无控制点，1表示有控制点；

3.2根据控制点集合C更新控制点约束光流和状态标识；

首先，根据控制点在图像中的坐标分别在中确定距离匹配特征点最近的像素位置，记为p_s＝(x_s,y_s)、p_t＝(x_t,y_t)；其中，x_s和x_t均为不大于W的正整数，y_s和y_t均为不大于H的正整数；

然后，按照如下公式更新u_gt、v_gt、ρ：

3.3构建控制光流向量U_gt、V_gt以及控制点状态矩阵G；

其中，U_gt、V_gt为WH维列向量，控制点状态矩阵G为WH×WH维对角方阵；U_gt为由所有像素位置上水平方向的控制光流u_gt按照从左至右、从上至下的顺序排列而成的光流向量，V_gt为由所有像素位置上垂直方向的控制光流v_gt按照从左至右、从上至下的顺序排列而成的光流向量，控制点状态矩阵G的对角线元素由所有像素位置上的控制光流状态ρ按照从左至右、从上至下的顺序排列而成；

第四步，构建初始光流向量和光流偏移向量

4.1光流和偏移初始化

根据控制点约束光流u_gt、v_gt，利用三线性插值方法计算图像各像素位置上的初始光流w(p)＝(u(p),v(p))，其中，w(p)表示像素p处的光流向量，u(p)和v(p)分别表示像素位置p处光流向量的水平分量和垂直分量；并将各像素位置上的初始光流偏移值置为零，即du(p)＝0、du(p)＝0；

4.2构建初始光流向量U和V

其中，U为由所有像素水平方向光流值按照从左至右、从上至下的顺序排列构建而成的光流向量，V为由所有像素垂直方向光流值按照从左至右、从上至下的顺序排列构建而成的光流向量；

4.3构建初始光流偏移向量dU和dV

dU为由所有像素水平方向光流偏移值按照从左至右、从上至下的顺序排列构建而成的光流向量，dV为由所有像素垂直方向光流偏移值按照从左至右、从上至下的顺序排列构建而成的光流向量；

第五步，光流向量迭代优化

光流向量迭代优化过程中，首先初始化迭代次数：N_out＝0，然后按照以下步骤迭代优化光流向量；

5.1按照如下公式更新光流向量、光流偏移向量：

5.2计算最佳光流偏移向量

计算最佳光流偏移向量的过程中，首先初始化迭代：N_in＝0，然后按照以下步骤迭代更新最佳光流偏移向量；

5.2.1计算权重矩阵

1)构建梯度图像以及偏移亮度差异图像图像各像素位置p的取值按照如下公式计算：

2)构建WH×WH维对角方阵I_x、I_y、I_z

I_x、I_y、I_z对角线上各元素分别由图像中各像素按照从左至右、从上至下的顺序填充而成；

3)构建权重矩阵T_f、T_g、T_h，上述三个矩阵均为WH×WH维对角方阵，对角线上第i个元素的取值由如下公式确定：

其中，符号[·]_i,i表示矩阵第i行、第i列元素，符号[·]_i表示向量第i个元素，D_x和D_y分别表示计算图像水平、垂直方向导数的WH×WH维矩阵算子，函数ψ的定义为：常数ε的取值范围是0.001～0.002；

5.2.2通过求解如下线性方程组，计算dU和dV

其中，

参数α的取值范围是0.01～0.03，参数β的取值范围是10～15；

5.2.3更新迭代次数：N_in←N_in+1；如果N_in＜2，则更新dU和dV，转到步骤5.2.1；否则，转到步骤5.3；

5.3更新迭代次数：N_out←N_out+1；如果N_out＜8，则转到步骤5.1；否则，输出最终光流：