CN106204637A - 光流计算方法 - Google Patents

Abstract

本发明提供了一种光流计算方法，首先用Fast角点检测算法确定图像中的Fast角点候选点，然后用Shi‑Tomasi算法从候选角点中提取Shi‑Tomasi强角点，最后对提取的Shi‑Tomasi强角点用L‑K稀疏光流计算方法计算光流，得到光流场。本发明方法能够提高光流精确度，减少光流计算量，大大减少运行时间，满足实时性要求。

Description

光流计算方法

技术领域

本发明涉及视觉图像处理技术领域，具体的涉及一种光流计算方法。

背景技术

快速光流计算方法是计算机视觉领域待解决的重要问题之一。光流是图像序列中像素点运动的瞬时速度，在运动跟踪，运动估算，无人机壁障等领域有着广泛应用。由于现有的大多数光流算法计算量大，计算耗时较长等原因，阻碍了光流计算在实践中的广泛应用。

现有的光流计算方法主要分为以下几类：基于梯度的方法、基于区域的匹配法、基于能量的方法、基于相位的方法和神经动力学法。现有光流计算一般采用基于梯度的光流计算方法，较经典的计算方法是Horn-Schunck(H-S)算法和Lucas-Kanade(L-K)算法两种。1981年，Horn和Schunck根据灰度一致性假设，推导出光流约束方程，在光流约束方程的基础上附加了光流场全局平滑约束，推导出了H-S光流算法，让光流算法得以发展。同年，Lucas和Kanade两位学者提出了一种迭代的图像匹配方法，该方法使用局部平滑假设，假设一个窗口内的所有像素具有相同的运动矢量，从而推导出L-K光流算法。

光流基本约束方程是在两个假设条件下推导得出的：(1)亮度恒定假设。该假设假定物体同一点的亮度在每帧之间不会发生改变；(2)像素点的运动是小运动。即时间的变化不会引起像素点位置的剧烈变化，灰度是连续变化的，这样灰度才能对位置求偏导。假设t时刻图像上某点(x,y)的灰度值为I(x,y,t)。在t+△t(△t→0)时刻，该点运动到新位置(x+△x,y+△y)，新位置的灰度记为I(x+△x,y+△y,t+△t)，根据假设(1)，同一点的灰度值应该相同，即：I(x,y,t)＝I(x+△x,y+△y,t+△t)。根据假设(2)对I(x,y,t)＝I(x+△x,y+△y,t+△t)右边进行泰勒展开式，忽略掉二阶及以上项，变形得：I_x·u+I_y·v+I_t＝0，该式即为光流基本约束方程，其中分别表示图像灰度空间对空间和时间的偏导数。分别表示光流的水平和垂直分量。光流基本约束方程将图像灰度的时空梯度与光流两个分量联系起来，但是一个方程不能求两个未知变量，需要引进其他约束，引进不同的约束就形成了不同的光流算法。如H-S算法就是在基本约束方程的基础上引进光流场全局平滑约束，该算法是将图像中的每一个像素点都与速度关联的一种稠密光流场计算方法，该方法计算所需时间较长，实时性差。L-K算法是在基本约束方程的基础上引进局部一阶约束，是一种稀疏光流计算方法。L-K算法只需要计算感兴趣点周围的小窗口局部信息就可得到感兴趣点的稀疏光流场。与稠密光流算法相比，稀疏光流算法不需要逐点计算光流，计算效率得到提高。在稀疏光流算法中，如何进一步提高计算出感兴趣点的速度，进一步降低整个光流计算过程所需时间，是我们应该探索的领域。

发明内容

本发明的目的在于提供一种光流计算方法，该发明解决了现有技术中光流计算耗时长，难以满足实时性要求的技术问题。

参见图1～2，本发明提供的光流计算方法，先通过Fast角点算法确定视频当前帧的Fast角点，然后用Shi-Tomasi算法对候选角点进行计算，快速确定强角点，最后对获取的强角点用L-k光流算法计算光流，得到光流场。

该算法与大多数需要计算每个像素点的稠密光流算法相比，大大减少了运行时间；与其他稀疏光流算法相比，在强角点的提取过程中通过先粗略选取Fast角点，再精确计算强角点，在节约了计算所需时间的同时，提高了计算所得光流结果的精确度。该算法具体流程如下：

参见图1，步骤S100：通过强角点检测算法提取所处理视频中图像I的角点，包括以下步骤：

步骤S110：初始化，从待处理视频中获取需要计算光流的图像I，调用灰度化函数将图像I灰度化；

步骤S120：对灰度化图像I用Fast角点检测算法确定灰度化图像I中的所有的Fast角点，并存储；

具体包括以下步骤：设定阈值t，t的值根据实际视频具体设置。以待求的像素点P_i为圆心，以3个像素为半径作圆，所得圆周上包含了16个像素点，设I(P_i)为P_i点的灰度值，Fast角点检测算法计算结果中：如果16个像素点中任意12个连续的像素点的像素灰度值大于I(P_i)+t或小于I(P_i)-t，则认为该像素点P_i为Fast角点，存储起来，否则，认为像素点P_i不是Fast角点。

步骤S130：用Shi-Tomasi算法逐一计算Fast角点并判断Fast角点是否为Shi-Tomasi强角点，得到最终需要的强角点，并存储：

如果Fast角点是Shi-Tomasi强角点，则将该点设置为最终需要的强角点；

如果Fast角点不是Shi-Tomasi强角点，则舍弃该点，遍历所有Fast角点，通过Shi-Tomasi角点检测算法计算得到多个最终需要的强角点；给定阈值λ。

参见图2，步骤S200：采用L-K光流算法计算所提取的图像I中的各最终需要的强角点对应的光流，遍历各最终需要的强角点，输出光流结果，求得光流场。

其中步骤S120包括以下步骤：选择灰度化后的图像I中的一个像素点P_i，设其灰度值为I(P_i)，以P_i为圆心，以半径为3个像素作圆，考虑圆周上的16个像素点，给定一个阈值t，当16个像素点上有12个连续的像素点像素灰度值大于I(p_i)+t或小于I(p_i)-t时，则确定该点为一个Fast角点。

现有的Fast角点检测算法比其他角点检测算法检测速度虽然快很多，但当所处理图像中噪声点较多时，效果并不理想。这主要是因为检测出的角点是依据阈值来确定的。而阈值的设定难以一次就获得较好的计算记过。所以，本发明提供方法中用Fast角点检测算法计算出Fast角点后，对这些Fast角点作进一步筛选，通过Shi-Tomasi角点检测算法计算出其中的强角点，得到最终需要的强角点。这样做能避免对每个像素点都进行Shi-Tomasi角点检测运算，而只需对fast算法计算出的Fast角点进行计算，这样做大大减少了Shi-Tomasi角点检测的计算量，节省了计算时间。使得本发明提供的方法中所用Fast角点检测算法比其他现有的角点检测算法对强角点的提取做到了既快又准。

该方法首先通过对当前帧图像I中的Fast角点进行粗略提取，从而实现了对计算时间的有效缩短，避免给长时间的提取角点，同时避免了现有方法中对所提取的大量角点进行大量计算的计算时间，从而提高计算效率。

优选的，步骤S130中对于Shi-Tomasi强角点的确定包括以下步骤：设图像I中像素点p(x,y)处的灰度为I(x,y)，以P点为中心建立一个n·n的窗口Ω，将窗口平移[△x,△y]，则该窗口内的灰度变化E[△x,△y]为：

$E\[\mathrm{\Δ}x,\mathrm{\Δ}y\]=\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y){\[I(x+\mathrm{\Δ}x,y+\mathrm{\Δ}y)-I(x,y)\]}^2---\left(1\right)$

其中，w(x,y)为窗口函数，在领域中心点具有较高值，边缘点值较小。对于微小的平移，将I(x+△x,y+△y)进行泰勒展开并忽略二阶及以上项，带入(1)式可得：

$E\[\mathrm{\Δ}x,\mathrm{\Δ}y\]={\mathrm{\Δx}}^2\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x^2+2\mathrm{\Δ}x\mathrm{\Δ}y\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x{I}_y+{\mathrm{\Δy}}^2\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_y^2---\left(2\right)$

其中，I_x、I_y分别表示图像灰度在x和y方向的偏导数，将(2)写成矩阵的形式：

$E(\mathrm{\Δ}x,\mathrm{\Δ}y)=\[\mathrm{\Δ}x,\mathrm{\Δ}y\]M\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\end{array}\right]---\left(3\right)$

其中M为2×2的矩阵：

$M=\left[\begin{array}{cc}\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x^2& \underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x{I}_y\\ \underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x{I}_y& \underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_y^2\end{array}\right]---\left(4\right)$

对于像素点p(x,y)，若矩阵M的两个特征值λ₁和λ₂中，相对较小的特征值大于给定的阈值λ，即λ₁≥λ₂且λ₂≥k·λ_2max(λ_2max是所有像素点较小特征值中的最大值，k为系数，k·λ_2max即为前面给定的阈值λ)，那么像素点p(x,y)为Shi-Tomasi强角点。如果像素点p(x,y)对于的矩阵M的两个特征值不满足上述条件，那么该像素点p(x,y)不是Shi-Tomasi强角点。

本发明提出的角点提取算法就是先用Fast角点算法计算得到Fast角点，然后遍历Fast角点，逐一对各Fast角点进行Shi-Tomasi角点计算，得到Shi-Tomasi强角点。强角点提取之后，利用L-K光流算法对求得的强角点计算光流。

我们知道，L-K算法是基于3个假设的，即亮度恒定、连续缓慢运动、空间一致。空间一致性即局部区域的像素运动是一致的，可以建立邻域像素的系统方程求解中心像素的运动。邻域窗口的大小可以自己选定，但是窗口太大会违背运动一致性的假设，太小会产生光流无法求解的孔径问题。优选的，L-K算法中邻域窗口为5×5。

具体的步骤S200包括以下步骤：

根据光流约束方程I_x·u+I_y·v+I_t＝0，其中u、v为像素点的光流分量，建立25个方程：

通过公式(6)来求解最小化的||Ad-b||²，即

(A^TA)d＝A^Tb (6)

由线性代数相关知识可知，当(AA^T)可逆时，方程组的解如下：

$\left[\begin{array}{c}u\\ v\end{array}\right]=\left(A^{T}A\right)A^{T}b---\left(7\right)$

当(A^TA)满秩，即秩为2时，(A^TA)有两个较大的特征向量，(A^TA)可逆。当图像中的纹理至少包含两个梯度方向时就可以满足该条件，求得方程的解，得到光流。

采用该角点提取方法的稀疏光流算法能减少计算量，提高计算效率。对提取的角点进行光流计算的具体计算方法按现有的L-K光流计算方法进行。

本发明的技术效果：

1、本发明提供的光流计算方法，采用基于Fast角点检测算法的候选角点提取，省去了逐点计算Shi-Tomasi角点的繁琐，简化了计算过程，从而使得强角点的提取既快又准，在保证计算精度的前提下，提高了计算效率。

2、本发明提供的光流计算方法，采用稀疏光流算法，对于图像特征点的稀疏光流算法比稠密光流算法计算量少，实时性好，能满足工程上实时获取计算结果的需求。

具体请参考根据本发明的光流计算方法提出的各种实施例的如下描述，将使得本发明的上述和其他方面显而易见。

附图说明

图1是本发明提供的光流计算方法中强角点检测算法的流程示意图；

图2是本发明提供的光流计算方法中对获取的强角点计算光流的流程示意图；

图3是本发明优选实施例和对比例中用于计算光流时所处理的连续两帧图像中的第一帧图像；

图4是本发明优选实施例和对比例中用于计算光流时所处理的连续两帧图像中的第二帧图像；

图5是本发明优选实施例和对比例中所得光流计算耗时结果示意图。

具体实施方式

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

下面通过实施例和对比例的比较，说明本发明提供的光流计算方法与现有光流计算法的计算速度的不同。本发明实施例中未详述部分均按现有方法进行计算和处理。

本实验的平台为Matlab7.1。从同一段视频中任取出两帧连续的图像作为本实施例和对比例中所处理图像对象(连续两帧中的第一帧为待处理图，第二帧为处理结果参考图)，如图3～4所示，图像大小为200*200像素。图5为第一帧图像计算得到的光流计算结果图，光流方向为右下方，可以知道，第二帧中的云图像较第一帧中的云图像发生了向右下方的移动。计算图像偏导数时梯度算子为：f_t为4点像素差。窗口函数w(x,y)选择高斯函数，方差为sigma＝1，窗口大小为5*5，Fast算法中的阈值t为20，Shi-Tomasi算法中的阈值k为0.04。

采用本发明提供方法处理图3～4所示图像，作为实施例。仅采用L-K光流算法进行光流计算作为对比例1，采用直接先提取Shi-Tomasi角点再用L-K光流算法计算Shi-Tomasi角点光流作为对比例2。对比例1～2中所用各方法均按背景技术中提到的该方法的常用步骤进行。对比例1：原始的L-K光流算法计算的运行时间为1.073秒。对比例2：直接先提取Shi-Tomasi角点再用L-K光流算法计算Shi-Tomasi角点光流的运行时间为0.302秒。采用本发明提供的光流计算方法，该实施例的运行时间为0.201秒。本发明提供的光流计算方法所需计算时间远低于现有光流计算方法。

实施例所得光流计算结果列于图5中，由图5可知，该方法能有效计算出光流的大小和方向。由此可见，本发明提供方法能满足工程上实时获取计算结果的要求。

本发明提出的只针对特征角点计算光流的稀疏光流计算方法大大节省了计算时间。在角点检测时，先用Fast角点检测得到Fast角点，再用Shi-Tomaasi角点检测算法对候选点进行计算，得到最终需要的角点，该方法能够快速、准确的获取特征角点，有很好的实用性。

本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。

通过对附图，说明书和权利要求书的研究，在实施本发明时本领域技术人员可以理解和实现所公开的实施例的变形。在权利要求书中，术语“包括”不排除其他步骤或元素，而不定冠词“一个”或“一种”不排除多个。在彼此不同的从属权利要求中引用的某些措施的事实不意味着这些措施的组合不能被有利地使用。权利要求书中的任何参考标记不构成对本发明的范围的限制。

Claims (4)

1.一种光流计算方法，其特征在于，包括以下步骤：

步骤S110：初始化，从待处理视频中获取需要计算光流的图像I，调用灰度化函数将图像I灰度化；

步骤S120：用Fast角点检测算法确定灰度化图像I中的所有的Fast角点，并存储；

所述步骤S120中包括以下步骤：设定阈值t，以待求的像素点P_i为圆心，以3个像素为半径作圆，所得圆周上包含了16个像素点，设I(P_i)为P_i点的灰度值，

Fast角点检测算法处理16个像素点，如果16个像素点中任意12个连续的像素点的像素灰度值均大于I(P_i)+t或小于I(P_i)-t，则设该像素点P_i为Fast角点，存储起来，否则，设该像素点P_i不是Fast角点，其中t为阈值；

步骤S130：用Shi-Tomasi算法逐一计算所述Fast角点并判断所述Fast角点是否为Shi-Tomasi强角点，得到最终需要的强角点，并存储：

如果所述Fast角点不是Shi-Tomasi强角点，则舍弃该点，遍历所有所述Fast角点，通过Shi-Tomasi角点检测算法计算得到多个所述最终需要的强角点；

步骤S200：采用L-K光流算法计算所述Shi-Tomasi强角点对应的光流，遍历各所述Shi-Tomasi强角点，输出光流结果，求得光流场。

2.根据权利要求1所述的光流计算方法，其特征在于，所述步骤S130中对于Shi-Tomasi强角点的确定包括以下步骤：

设图像I中像素点p(x,y)处的灰度为I(x,y)，以P点为中心建立一个n·n的窗口Ω，将该窗口平移[△x,△y]，则该窗口内的灰度变化E[△x,△y]为：

$E\[\mathrm{\Δ}x,\mathrm{\Δ}y\]=\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y){\[I(x+\mathrm{\Δ}x,y+\mathrm{\Δ}y)-I(x,y)\]}^2---\left(1\right)$

其中，w(x,y)为窗口函数，在领域中心点具有较高值，边缘点值较小。对于微小的平移，将I(x+△x,y+△y)进行泰勒展开并忽略二阶及以上项，带入(1)式可得：

$E\[\mathrm{\Δ}x,\mathrm{\Δ}y\]={\mathrm{\Δx}}^2\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x^2+2\mathrm{\Δ}x\mathrm{\Δ}y\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x{I}_y+{\mathrm{\Δy}}^2\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_y^2---\left(2\right)$

其中，I_x、I_y分别表示图像灰度在x和y方向的偏导数，将(2)写成矩阵的形式：

$E(\mathrm{\Δ}x,\mathrm{\Δ}y)=\[\mathrm{\Δ}x,\mathrm{\Δ}y\]M\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\end{array}\right]---\left(3\right)$

其中M为2×2的矩阵：

$M=\left[\begin{array}{cc}\underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x^2& \underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x{I}_y\\ \underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_x{I}_y& \underset{(x,y)\∈\mathrm{\Ω}}{\Σ}w(x,y)I_y^2\end{array}\right]---\left(4\right)$

对于像素点p(x，y)，若λ₁和λ₂满足λ₁≥λ₂且λ₂≥k·λ_2max，那么像素点p(x，y)为Shi-Tomasi强角点；

如果λ₁和λ₂不满足上述条件，那么该像素点p(x，y)不是Shi-Tomasi强角点；

其中，λ₁和λ₂为素点p(x，y)对于矩阵M的两个特征值，λ_2max是所有像素点较小特征值中的最大值，k为系数，k·λ_2max即为所述阈值λ。

3.根据权利要求2所述的光流计算方法，其特征在于，所述L-K算法中邻域窗口为5×5。

4.根据权利要求3所述的光流计算方法，其特征在于，所述步骤S200中包括以下步骤：

根据光流约束方程I_x·u+I_y·v+I_t＝0，其中u、v为像素点的光流分量，建立25个方程：

通过公式(6)来求解最小化的||Ad-b||²，即

(A^TA)d＝A^Tb (6)

当(AA^T)可逆时，方程组的解如下：

$\left[\begin{array}{c}u\\ v\end{array}\right]=\left(A^{T}A\right)A^{T}b---\left(7\right)$

当图像I的纹理至少包含两个梯度方向时，(A^TA)秩为2时，(A^TA)有两个较大的特征向量，(A^TA)可逆，求得方程(7)的解，得到光流。