[发明专利]一种面向大尺度运动的光流估计方法

权利要求书

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₁，即中所有特征点的最小外接矩形面积与图像的面积之比；

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

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、ρ：

$\left\{\begin{array}{c}{u}_{gt}\left(p_s\right)=x_t-x_s\\ v_{gt}\left(p_s\right)=y_t-y_s\\ \mathrm{\ρ}\left(p_s\right)=1\end{array}\right.$

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按照如下公式更新光流向量、光流偏移向量：

$\left\{\begin{array}{c}U\←U+dU\\ V\←V+dV\\ dU\←0\\ dV\←0\end{array}\right.$

5.2计算最佳光流偏移向量

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

5.2.1计算权重矩阵

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

$\left\{\begin{array}{c}{I}_x\left(p\right)=\frac{\∂I_2(p+w(p\left)\right)}{\∂x}\\ I_y\left(p\right)=\frac{\∂I_2(p+w(p\left)\right)}{\∂y}\\ I_z\left(p\right)=I_2(p+w\left(p\right))-I_1\left(p\right)\end{array}\right.$

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个元素的取值由如下公式确定：

$\left\{\begin{array}{c}{\[T_f\]}_{i,i}={\mathrm{\ψ}}^{\′}\left({\[I_z+I_{x}dU+I_{y}dV\]}_i^2\right)\\ {\[T_g\]}_{i,i}={\mathrm{\ψ}}^{\′}({\[D_x(U+dU)\]}_i^2+{\[D_x(V+dV)\]}_i^2+{\[D_y(U+dU)\]}_i^2+{\[D_y(V+dV)\]}_i^2)\\ {\[T_h\]}_{i,i}={\mathrm{\ψ}}^{\′}({\[G(U+dU-U_{gt})\]}_i^2+{\[G(V+dV-V_{gt})\]}_i^2)\end{array}\right.$

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

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

$\left[\begin{array}{cc}{T}_f{I}_x^2+\mathrm{\α}M+\mathrm{\β}N& T_f{I}_x{I}_y\\ T_f{I}_x{I}_y& T_f{I}_y^2+\mathrm{\α}M+\mathrm{\β}N\end{array}\right]\left[\begin{array}{c}dU\\ dV\end{array}\right]=-\left[\begin{array}{c}{T}_f{I}_x{I}_z+\mathrm{\α}MU+\mathrm{\β}N(U-U_{gt})\\ T_f{I}_y{I}_z+\mathrm{\α}MV+\mathrm{\β}N(V-V_{gt})\end{array}\right]$

其中，

$\left\{\begin{array}{c}M=D_x^T{T}_g{D}_x+D_y^T{T}_g{D}_y\\ N=G^T{T}_{h}G\end{array}\right.$

参数α的取值范围是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；否则，输出最终光流：

$\left\{\begin{array}{c}U\←U+dU\\ V\←V+dV\end{array}\right..$