CN105139439A - 一种光流场计算方法 - Google Patents

Abstract

本发明公开一种光流场计算方法，通过采用4邻域模板计算像素点光流均值时，该像素点的下邻域及左邻域光流的当次迭代新值已经求出，因此计算该点光流均值时利用该新值，根据计算出来的光流均值可以求出该点的当次迭代光流，本发明的方法不仅误差小，而且加快了算法的收敛速度。本发明的方法尤其适用于具有较多边界，边界灰度梯度比较大的情况，能大大减少迭代次数，缩减迭代时间，和传统微分法相比，本发明的方法具有很好的鲁棒性，能大大的减少光流场计算时间，效率高。

Description

一种光流场计算方法

技术领域

本发明属于图像处理技术领域，具体涉及一种快速光流场计算方法。

背景技术

光流场，它是指图像中所有像素点构成的一种二维(2D)瞬时速度场，其中的二维速度矢量是景物中可见点的三维速度矢量在成像表面的投影。所以光流不仅包含了被观察物体的运动信息，而且还包含有关景物三维结构的丰富信息。对光流的研究成为计算机视觉及有关研究领域中的一个重要部分。因为在计算机视觉中，光流扮演着重要角色，在目标对象分割、识别、跟踪、机器人导航以及形状信息恢复等都有着非常重要的应用。从光流中恢复物体三维结构和运动则是计算机视觉研究所面临的最富有意义和挑战性的任务之一。

1981年，Horn和Schunck最早研究光流场的计算。到目前为止，各种各样的方法和改进方法超过几十种，这些方法在概念上和性能上都有着很大的差异，对现有方法述缺乏一个系统的分类。习惯从概念上将现有方法分为4类：基于梯度方法(微分法)、基于区域方法(匹配法)、基于能量方法(能量法)、基于相位方法(相位法)和神经动力学法五种。但是传统的光流场计算方法，迭代次数多，计算时间长，效率底下。

发明内容

本发明为解决上述技术问题，提出一种光流场计算方法，通过改进的4邻域模板，计算光流场均值，从而计算光流，有效地减少迭代次数，大大的减少了迭代时间。

本发明采用的技术方案是：一种光流场计算方法，具体包括以下步骤：

S1：初始化各像素点迭代初值；

S2：采用4邻域模版计算某像素点光流均值时，根据该像素点的下邻域及左邻域当次迭代新值计算该像素点的光流均值；

S3：有步骤S2得到的该像素点的光流均值，计算该像素点当前迭代的光流值；

S4：重复步骤S2至S3得到当次迭代各像素点的光流值，并判断各像素点当次迭代得到的光流值与前一次迭代得到的光流值差值是否小于或等于给定估计容差，若是则停止迭代，否则转至步骤S2进入下一次迭代。

进一步地，所述步骤S2计算某像素点光流场均值的表达式为：

${\stackrel{\‾}{u}}^{(n-1)}(i,j)=\left(u^{\left(n\right)}\right(i-1,j)+u^{\left(n\right)}(i,j-1)+u^{(n-1)}(i+1,j)+u^{(n-1)}(i,j+1\left)\right)/4$

${\stackrel{\‾}{v}}^{(n-1)}(i,j)=\left(v^{\left(n\right)}\right(i-1,j)+v^{\left(n\right)}(i,j-1)+v^{(n-1)}(i+1,j)+v^{(n-1)}(i,j+1\left)\right)/4$

其中，u,v表示图像中像素点(i,j)点处分别沿x轴，y轴方向的光流速，是u,v在(i,j)点4邻域的平均值，n表示迭代次数，i表示像素点x轴坐标，j表示像素点y轴坐标。

进一步地，所述步骤S3计算该像素点当前迭代的光流值的表达式为：

$u^{\left(n\right)}={\stackrel{\‾}{u}}^{(n-1)}-\frac{I_x(I_x{\stackrel{\‾}{u}}^{(n-1)}+I_y{\stackrel{\‾}{v}}^{(n-1)}+I_t)}{({\mathrm{\α}}^2+{I_x}^2+{I_y}^2)}$

$v^{\left(n\right)}={\stackrel{\‾}{v}}^{(n-1)}-\frac{I_y(I_x{\stackrel{\‾}{u}}^{(n-1)}+I_y{\stackrel{\‾}{v}}^{(n-1)}+I_t)}{({\mathrm{\α}}^2+{I_x}^2+{I_y}^2)}$

其中，I_x为图像灰度对x的微分，I_y为图像灰度对y的微分，I_t为图像灰度对t的微分，α表示权重值。

进一步地，判断各像素点当次迭代得到的光流值与前一次迭代得到的光流值差值是否小于或等于给定容差，即：

|u⁽ⁿ⁾(i,j)-u^(n-1)(i,j)|≤ε₁

|v⁽ⁿ⁾(i,j)-v^(n-1)(i,j)|≤ε₂

其中，ε₁，ε₂为给定的估计容差，|·|表示取绝对值运算。

本发明的有益效果：本发明的一种光流场计算方法，通过采用4邻域模板计算像素点光流均值时，该像素点的下邻域及左邻域当次迭代新值已经求出，利用已经求出的迭代新值，计算该像素点的光流均值，本发明的方法不仅误差小，而且加快了算法的收敛速度。本发明的方法尤其适用于具有较多边界，边界灰度梯度比较大的情况，能大大减少迭代次数，缩减迭代时间，和传统微分法相比，本发明的方法具有很好的鲁棒性，能大大的减少光流场计算时间，效率高。

附图说明

图1为Horn模板示意图。

图2为4邻域模板示意图。

图3为本发明实施例提供的3*3的数字图像。

图4为本发明提供的视频截取连续两帧图像的第一帧示意图。

图5为本发明提供的视频截取连续两帧图像的第二帧示意图。

图6为传统方法得到的光流场。

图7为本发明方法得到的光流场。

具体实施方式

为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。

本发明在基本约束方程和全局平滑性约束条件下，提出了一种光流场计算方法。在传统的微分算法基础上，采用4邻域模板计算均值时，当计算某像素点光流场均值时，它的下邻域和左邻域当次迭代的新值已经计算出来，在计算某像素点光流均值时利用已经计算出的新值，能加速算法的收敛速度。

微分法原理简单，实现起来比较容易，效率高，而且能够计算结果相当精确，所以得到了广泛的研究。下面，详细介绍一下微分法的原理。

微分法，又叫基于梯度的光流场算法，是由图像像素的时空微分得到每个像素点的速度。其中经典的是Horn‐Schunck算法，Horn等人根据连续帧图像之间的间隔时间很小(几十毫秒之内)，两幅图像中相应像素点的灰度变化也很小，得到光流场的基本方程，再由全局平滑假设计算出每个像素点的速度矢量。

记时刻t时，图像上(x,y)点处的灰度值为I(x,y,t)。在经过δt时间后，该点在图像上的坐标变为(x+δx,y+δy)，灰度值记为I(x+δx,y+δy,t+δt)，假定它与I(x,y,t)相等，即：

I(x,y,t)＝I(x+δx,y+δy,t+δt)(2‐15)

将右边用泰勒公式展开，经化简和略去两次以上的项，得：

$I(x+\mathrm{\δ}x,y+\mathrm{\δ}y,t+\mathrm{\δ}t)=I(x,y,t)+\frac{\∂I}{\∂x}\mathrm{\δ}x+\frac{\∂I}{\∂y}\mathrm{\δ}y+\frac{\∂I}{\∂t}\mathrm{\δ}t---(2-16)$

则(2‐15)可化为：

$\frac{\∂I}{\∂x}\mathrm{\δ}x+\frac{\∂I}{\∂y}\mathrm{\δ}y+\frac{\∂I}{\∂t}\mathrm{\δ}t=0---(2-17)$

两边同除以δt得：

$\frac{\∂I}{\∂x}\frac{\mathrm{\δ}x}{\mathrm{\δ}t}+\frac{\∂I}{\∂y}\frac{\mathrm{\δ}y}{\mathrm{\δ}t}+\frac{\∂I}{\∂t}=0---(2-18)$

令为沿x方向的分量，为沿y方向的分量，则有：

$\frac{\∂I}{\∂x}u+\frac{\∂I}{\∂y}v+\frac{\∂I}{\∂t}=0---(2-19)$

这就是光流场计算的基本公式。写成梯度形式为：

其中，U＝(u,v)^T即为光流，u,v为(x,y)点处分别沿x,y方向的光流速。这是一个不适定方程，要想求解u,v，有全局平滑假设和局部平滑假设两种方法。

以下是以全局平滑假设为基础的光流场的解法。

Horn和Schunch提出了同一物体引起的光流场尽可能连续的、平滑的假设，即要求偏离平滑的误差最小，则有：

${E_s}^2=\∫\∫{({\left(\frac{\∂u}{\∂x}\right)}^2+{\left(\frac{\∂u}{\∂y}\right)}^2+{\left(\frac{\∂v}{\∂x}\right)}^2+{\left(\frac{\∂v}{\∂y}\right)}^2)}^{2}dxdy=min---(2-21)$

E_s为偏离平滑性的误差。同时应该满足偏离光流基本等式的误差，即：

${E_c}^2=\∫\∫{(\▿I)}^T\·U+I_t{)}^{2}dxdy---(2-22)$

于是光流场(u,v)应满足：

$\∫\∫({\left({\left(\▿I\right)}^T\·U+I_t\right)}^2+{\mathrm{\α}}^2{\left({\left(\frac{\∂u}{\∂x}\right)}^2+{\left(\frac{\∂u}{\∂y}\right)}^2+{\left(\frac{\∂v}{\∂x}\right)}^2+{\left(\frac{\∂v}{\∂y}\right)}^2\right)}^2)dxdy=min---(2-23)$

其中，α为权值，表示图像数据和约束条件之间的权重。当原图像噪声比较大时，原图像数据受到噪声干扰，可信度不高，此时α取值较大，更多的依赖于平滑约束条件，当原图像数据比较精确时，则应更重视原图像数据，此时α取较小值。

下面只考虑离散的情况，因为在图像处理中，多是处理的数字图像。在图像上的一点(i,j)及其4邻域上，偏离光流基本方程的误差为：

E_c ²(i,j)＝(I_xu(i,j)+I_yv(i,j)+I_t)²(2‐24)

又有： $\frac{\∂u}{\∂x}=u(i+1,j)-u(i,j)---(2-25)$

则： ${\left(\frac{\∂u}{\∂x}\right)}^2=\frac{1}{2}{\left(\frac{\∂u}{\∂x}\right)}^2+\frac{1}{2}{\left(\frac{\∂u}{\∂x}\right)}^2=\frac{{\left(u\left(i+1,j\right)-u\left(i,j\right)\right)}^2}{2}+\frac{{\left(u\left(i,j\right)-u\left(i-1,j\right)\right)}^2}{2}---\left(2-26\right)$

偏离平滑的误差为：

$\begin{array}{c}{E_s}^2\left(i,j\right)=\frac{1}{4}\[{\left(u\left(i+1,j\right)-u\left(i,j\right)\right)}^2+{\left(u\left(i,j\right)-u\left(i-1,j\right)\right)}^2+\\ {\left(u\left(i,j+1\right)-u\left(i,j\right)\right)}^2+{\left(u\left(i,j\right)-u\left(i,j-1\right)\right)}^2+\\ {\left(v\left(i+1,j\right)-v\left(i,j+1\right)\right)}^2+{\left(v\left(i,j\right)-v\left(i-1,j\right)\right)}^2+\\ {\left(v\left(i,j+1\right)-v\left(i,j\right)\right)}^2+{\left(v\left(i,j\right)-v\left(i,j-1\right)\right)}^2\]\end{array}---\left(2-27\right)$

则总的误差为： $E=\underset{i}{\Σ}\underset{j}{\Σ}\left({E_c}^2\right(i.j)+{\mathrm{\α}}^2{E_S}^2(i,j\left)\right)---(2-28)$ 要取其极小值，则必有： $\frac{\∂E}{\∂u}=0---(2-29)$

$\frac{\∂E}{\∂v}=0---(2-30)$

则有： $\frac{\∂E}{\∂u}=2(I_{x}u+I_{y}v+I_t)I_x+2{\mathrm{\α}}^2(u-\stackrel{\‾}{u})=0---(2-31)$

$\frac{\∂E}{\∂v}=2(I_{x}u+I_{y}v+I_t)I_y+2{\mathrm{\α}}^2(v-\stackrel{\‾}{v})=0---(2-32)$

将(2‐31)(2‐32)式展开并整理得： $({I_x}^2+{\mathrm{\α}}^2)u+I_x{I}_{y}v={\mathrm{\α}}^2\stackrel{\‾}{u}-I_x{I}_t---(2-33)$

$I_x{I}_{y}u+({I_y}^2+{\mathrm{\α}}^2)v={\mathrm{\α}}^2\stackrel{\‾}{v}-I_y{I}_t---(2-34)$

利用克莱姆法则求解上面方程组的解，有：

$u=\frac{D_1}{D}=\frac{\left|\begin{array}{cc}{\mathrm{\α}}^2\stackrel{\‾}{u}-I_x{I}_t& I_x{I}_y\\ {\mathrm{\α}}^2\stackrel{\‾}{v}-I_y{I}_t& {I_y}^2+{\mathrm{\α}}^2\end{array}\right|}{\left|\begin{array}{cc}{I_x}^2+{\mathrm{\α}}^2& I_x{I}_y\\ I_x{I}_y& {I_y}^2+{\mathrm{\α}}^2\end{array}\right|}---(2-35)$

$v=\frac{D_2}{D}=\frac{\left|\begin{array}{cc}{I_x}^2+{\mathrm{\α}}^2& {\mathrm{\α}}^2\stackrel{\‾}{u}-I_x{I}_t\\ I_x{I}_y& {\mathrm{\α}}^2\stackrel{\‾}{v}-I_y{I}_t\end{array}\right|}{\left|\begin{array}{cc}{I_x}^2+{\mathrm{\α}}^2& I_x{I}_y\\ I_x{I}_y& {I_y}^2+{\mathrm{\α}}^2\end{array}\right|}---(2-36)$

求解得： $u=\stackrel{\‾}{u}-I_x\frac{I_x\stackrel{\‾}{u}+I_y\stackrel{\‾}{v}+I_t}{{\mathrm{\α}}^2+{I_x}^2+{I_y}^2}---(2-37)$

$v=\stackrel{\‾}{v}-I_y\frac{I_x\stackrel{\‾}{u}+I_y\stackrel{\‾}{v}+I_t}{{\mathrm{\α}}^2+{I_x}^2+{I_y}^2}---(2-38)$

利用迭代法，设方程组Ax＝b有唯一解x^*，将Ax＝b变形为等价的方程组：

x＝Bx+f(2‐39)

对于第i个方程解出x_i，得到与原方程组等价的方程组：

$x_i=\frac{1}{a_{ii}}\[-\underset{j\≠i}{\Σ}{a}_{ij}{x}_j+b_i\],(i=1,2,......,n)---(2-40)$

迭代格式为： ${x_i}^{(n+1)}=\frac{1}{a_{ii}}\[-\underset{j\≠i}{\Σ}{a}_{ij}{x_j}^{\left(n\right)}+b_i\],(i=1,2,......,n)---(2-41)$

矩阵形式为：

x⁽ⁿ⁺¹⁾＝Bx⁽ⁿ⁾+f(2‐42)

给定初值x⁽⁰⁾，按此公式计算得到近似解向量序列{x^(k)}，当迭代次数无限增加时，序列{x^(k)}都有相同的极限x^*，即： $\underset{k\→\mathrm{\∞}}{\lim }{x}^{\left(k\right)}=x^{*}---(2-43)$

则由上述迭代法得到的解如下：

$u^{(n+1)}={\stackrel{\‾}{u}}^{\left(n\right)}-\frac{I_x(I_x{\stackrel{\‾}{u}}^{\left(n\right)}+I_y{\stackrel{\‾}{v}}^{\left(n\right)}+I_t)}{({\mathrm{\α}}^2+{I_x}^2+{I_y}^2)}---(2-44)$

同理： $v^{(n+1)}={\stackrel{\‾}{v}}^{\left(n\right)}-\frac{I_y(I_x{\stackrel{\‾}{u}}^{\left(n\right)}+I_y{\stackrel{\‾}{v}}^{\left(n\right)}+I_t)}{({\mathrm{\α}}^2+{I_x}^2+{I_y}^2)}---(2-45)$

其中，I_x,I_y,I_t为图像灰度对空间时间的微分，I_i,j,t表示在(i,j,t)点的灰度值I(x,y,t)。即：

$\begin{array}{c}{I}_x=\frac{1}{4}\[(I_{i+1,j,t}-I_{i,j,t})+(I_{i+1,j+1,t}-I_{i,j+1,t})+\\ \left(I_{i+1,j,t+1}-I_{i,j,t+1}\right)+\left(I_{i+1,j+1,t+1}-I_{i,j+1,t+1}\right)\]\end{array}---(2-46)$

$\begin{array}{c}{I}_y=\frac{1}{4}\[(I_{i,j+1,t}-I_{i,j,t})+(I_{i+1,j+1,t}-I_{i+1,j,t})+\\ \left(I_{i,j+1,t+1}-I_{i,j,t+1}\right)+\left(I_{i+1,j+1,t+1}-I_{i+1,j,t+1}\right)\]\end{array}---(2-47)$

$\begin{array}{c}{I}_t=\frac{1}{4}\[(I_{i,j,t+1}-I_{i,j,t})+(I_{i+1,j,t+1}-I_{i+1,j,t})+\\ \left(I_{i,j+1,t+1}-I_{i,j+1,t}\right)+\left(I_{i+1,j+1,t+1}-I_{i+1,j+1,t}\right)\]\end{array}---(2-48)$

n+1表示第n+1次迭代，是u,v在(x,y)点4邻域的平均值。u⁽⁰⁾,v⁽⁰⁾为光流的初始值，一般取0，当迭代结果满足事先给定的估计容差ε₁，ε₂时，迭代结束，即：

|u⁽ⁿ⁾-u^(n-1)|≤ε₁(2‐49)

|v⁽ⁿ⁾-v^(n-1)|≤ε₂(2‐50)

此时得到了图像中每个像素点的光流分量u,v。

在式(2‐44)和式(2‐45)中计算时，通常有两种方法。第一种方法是采用Horn模板，第二种通常用的方法是采用4邻域模板。

Horn在计算平均值时采用了8邻域模板，模板如图1所示，即：

$\begin{array}{c}\stackrel{\‾}{u}\left(i,j\right)=\left(u\left(i-1,j\right)+u\left(i,j-1\right)+u\left(i+1,j\right)+u\left(i,j+1\right)\right)/6+\\ \left(u\left(i-1,j-1\right)+u\left(i-1,j+1\right)+u\left(i+1,j+1\right)+u\left(i+1,j-1\right)\right)/12\end{array}---\left(2-51\right)$

$\begin{array}{c}\stackrel{\‾}{v}\left(i,j\right)=\left(v\left(i-1,j\right)+v\left(i,j-1\right)+v\left(i+1,j\right)+v\left(i,j+1\right)\right)/6+\\ \left(v\left(i-1,j-1\right)+v\left(i-1,j+1\right)+v\left(i+1,j+1\right)+v\left(i+1,j-1\right)\right)/12\end{array}---\left(2-52\right)$

模板中的权值越大，表示在计算时占的比重越大，所以Horn模板可以简化为4邻域模板。

现在计算均值时经常使用到的是4邻域模板，模板如图2所示，即：

${\stackrel{\‾}{u}}^{(n-1)}(i,j)=\left(u^{(n-1)}\right(i-1,j)+u^{(n-1)}(i,j-1)+u^{(n-1)}(i+1,j)+u^{(n-1)}(i,j+1\left)\right)/4---(2-53)$

${\stackrel{\‾}{v}}^{(n-1)}(i,j)=\left(v^{(n-1)}\right(i-1,j)+v^{(n-1)}(i,j-1)+v^{(n-1)}(i+1,j)+v^{(n-1)}(i,j+1\left)\right)/4---(2-54)$

本发明的一种光流场计算方法，具体方案为：

S1：初始化初始化各像素点迭代初值，例如初始值u⁽⁰⁾,v⁽⁰⁾，一般取0。

S2：本发明对4邻域模板进行改进，在采用4邻域模板计算均值时，其实每次计算时，它的下邻域和左邻域当次迭代的新值已经计算出来。

如图5所示为3*3的数字图像，一共有9个像素点，从(1,1)开始从左往右，从下往上求得9个像素点的(u,v)，当进行第二次迭代时，例如求点(2,2)的(u⁽²⁾(2,2),v⁽²⁾(2,2))时，按照计算顺序，点(2,2)的下邻域(1,2)和左邻域(2,1)的第二次迭代的光流值已经计算出来，即(u⁽²⁾(1,2),v⁽²⁾(1,2))与(u⁽²⁾(2,1),v⁽²⁾(2,1))。

因此，在本申请中计算时，由于u⁽ⁿ⁾(i-1,j)，u⁽ⁿ⁾(i,j-1)已算出，计算时，由于v⁽ⁿ⁾(i-1,j)，v⁽ⁿ⁾(i,j-1)已算出，故可以采用最新的值。

即，将(2‐53)式做如下修改，用该像素点的下邻域与左邻域的当前迭代值u⁽ⁿ⁾(i-1,j)，u⁽ⁿ⁾(i,j-1)，而不用u^(n-1)(i-1,j)，u^(n-1)(i,j-1)，将(2‐54)式做如下修改，用该像素点的下邻域与左邻域的当前迭代值v⁽ⁿ⁾(i-1,j)，v⁽ⁿ⁾(i,j-1)，而不用v^(n-1)(i-1,j)，v^(n-1)(i,j-1)，即：

${\stackrel{\‾}{u}}^{(n-1)}(i,j)=\left(u^{\left(n\right)}\right(i-1,j)+u^{\left(n\right)}(i,j-1)+u^{(n-1)}(i+1,j)+u^{(n-1)}(i,j+1\left)\right)/4---(2-55)$

${\stackrel{\‾}{v}}^{(n-1)}(i,j)=\left(v^{\left(n\right)}\right(i-1,j)+v^{\left(n\right)}(i,j-1)+v^{(n-1)}(i+1,j)+v^{(n-1)}(i,j+1\left)\right)/4---(2-56)$

计算时，采用这种改进的方法。

S3：将计算出的带入式(2‐44)和式(2‐45)迭代，得到当前迭代的光流值，表达式如下

$u^{\left(n\right)}={\stackrel{\‾}{u}}^{(n-1)}-\frac{I_x(I_x{\stackrel{\‾}{u}}^{(n-1)}+I_y{\stackrel{\‾}{v}}^{(n-1)}+I_t)}{({\mathrm{\α}}^2+{I_x}^2+{I_y}^2)}$

$v^{\left(n\right)}={\stackrel{\‾}{v}}^{(n-1)}-\frac{I_y(I_x{\stackrel{\‾}{u}}^{(n-1)}+I_y{\stackrel{\‾}{v}}^{(n-1)}+I_t)}{({\mathrm{\α}}^2+{I_x}^2+{I_y}^2)}.$

S4：重复步骤S2至S3得到当前迭代各像素点的光流值，当各像素点当前迭代的光流值均满足式(2‐49)和式(2‐50)时，迭代停止，此时便得到每个像素点的光流分量u,v。这样能够较好的提高算法的速度，减少迭代次数，特别是当计算精度要求很高时，改进的算法的优势更明显。

在计算均值时，不论是4邻域模板还是Horn模板，在计算时背景和运动物体的模板权值都是一样的。实际上，在背景和运动物体的边缘，像素灰度值的变化是非常大的，为了防止运动物体的梯度数据向背景的梯度数据扩散，夏毓鹏等提出了选取e^x作为每个像素点的光流分量的权值，即：

$\begin{array}{c}{\stackrel{\‾}{u}}^{\left(n-1\right)}\left(i,j\right)=(e^{-d_{i-1,j}}{u}^{\left(n-1\right)}\left(i-1,j\right)+e^{-d_{i,j-1}}{u}^{\left(n-1\right)}\left(i,j-1\right)+e^{-d_{i+1,j}}{u}^{\left(n-1\right)}\left(i+1,j\right)+\\ e^{-d_{i,j+1}}{u}^{\left(n-1\right)}\left(i,j+1\right))/(4\left(e^{-d_{i-1,j}}+e^{-d_{i,j-1}}+e^{-d_{i+1,j}}+e^{-d_{i,j+1}}\right))\end{array}---\left(2-57\right)$

$\begin{array}{c}{\stackrel{\‾}{v}}^{\left(n-1\right)}\left(i,j\right)=(e^{-d_{i-1,j}}{v}^{\left(n-1\right)}\left(i-1,j\right)+e^{-d_{i,j-1}}{v}^{\left(n-1\right)}\left(i,j-1\right)+e^{-d_{i+1,j}}{v}^{\left(n-1\right)}\left(i+1,j\right)+\\ e^{-d_{i,j+1}}{v}^{\left(n-1\right)}\left(i,j+1\right))/(4\left(e^{-d_{i-1,j}}+e^{-d_{i,j-1}}+e^{-d_{i+1,j}}+e^{-d_{i,j+1}}\right))\end{array}---\left(2-58\right)$

其中：

d_i-1,j＝|I_i,j-I_i-1,j|(2‐59)

d_i,j-1＝|I_i,j-I_i,j-1|(2‐60)

d_i+1,j＝|I_i+1,j-I_i,j|(2‐61)

d_i,j+1＝|I_i,j+1-I_i,j|(2‐62)

采用公式(2‐57)和(2‐58)来计算均值从而计算光流场，能够抑制梯度数据的盲目扩散。

将改进的4模板邻域算法应用到边界收敛的微分法中，能够加快迭代，减少迭代的次数，则式(2‐57)和式(2‐58)可变为：

$\begin{array}{c}{\stackrel{\‾}{u}}^{\left(n-1\right)}\left(i,j\right)=(e^{-d_{i-1,j}}{u}^{\left(n\right)}\left(i-1,j\right)+e^{-d_{i,j-1}}{u}^{\left(n\right)}\left(i,j-1\right)+e^{-d_{i+1,j}}{u}^{\left(n\right)}\left(i+1,j\right)+\\ e^{-d_{i,j+1}}{u}^{\left(n-1\right)}\left(i,j+1\right))/(4\left(e^{-d_{i-1,j}}+e^{-d_{i,j-1}}+e^{-d_{i+1,j}}+e^{-d_{i,j+1}}\right))\end{array}---\left(2-63\right)$

$\begin{array}{c}{\stackrel{\‾}{v}}^{\left(n-1\right)}\left(i,j\right)=(e^{-d_{i-1,j}}{v}^{\left(n\right)}\left(i-1,j\right)+e^{-d_{i,j-1}}{v}^{\left(n\right)}\left(i,j-1\right)+e^{-d_{i+1,j}}{v}^{\left(n-1\right)}\left(i+1,j\right)+\\ e^{-d_{i,j+1}}{v}^{\left(n-1\right)}\left(i,j+1\right))/(4\left(e^{-d_{i-1,j}}+e^{-d_{i,j-1}}+e^{-d_{i+1,j}}+e^{-d_{i,j+1}}\right))\end{array}---\left(2-64\right)$

利用改进的边界收敛的微分法计算的结果表明：该方法迭代次数大大减少，总的迭代的时间也大大减少，并且这种边界收敛的方法尤其适用于具有较多边界，边界灰度梯度比较大的情况。

下面通过实验比较，说明本发明的一种光流场计算方法与传统的采用一般的4邻域模板计算均值从而计算光流的方法的优势。

从任一段视频中截取出一些连续帧图像，图像的大小为1024×576，现取连续的两帧图像如图3和图4所示。

利用传统方法计算得到的光流场图如图5所示，利用改进方法计算得到的光流场图如图6所示。由图5、图6可知两种方法得到的光流场误差很小，计算得，对于传统方法，u的均值为0.011728像素/秒，均方差为0.365601像素/秒，v的均值为‐0.406717像素/秒，均方差为0.810914像素/秒。对于改进方法，u的均值为0.019390像素/秒，均方差为0.350750像素/秒，v的均值为‐0.418672像素/秒，均方差为0.830516像素/秒。

利用传统方法与改进方法分别计算两帧图像的光流场。当权值α＝20，误差容限ε₁＝ε₂＝0.01时，传统方法需要迭代183次，共计迭代时间为52秒，改进方法只需迭代153次，共计迭代时间38秒。

总之，本发明的一种光流场计算方法能够有效的减少迭代次数，大大的减少了迭代时间，并且与传统方法得到的光流场误差很小。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (4)

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

S1：初始化各像素点迭代初值；

S2：采用4邻域模版计算某像素点光流均值时，根据该像素点的下邻域及左邻域当次迭代新值计算该像素点的光流均值；

S3：由步骤S2得到的该像素点的光流均值，计算该像素点当前迭代的光流值；

S4：重复步骤S2至S3得到当次迭代各像素点的光流值，并判断各像素点当次迭代得到的光流值与前一次迭代得到的光流值差值是否小于或等于给定估计容差，若是则停止迭代，否则转至步骤S2进入下一次迭代。

2.根据权利要求1所述的一种光流场计算方法，其特征在于，所述步骤S2计算某像素点光流场均值的表达式为：

${\stackrel{\‾}{u}}^{(n-1)}(i,j)=\left(u^{\left(n\right)}\right(i-1,j)+u^{\left(n\right)}(i,j-1)+u^{(n-1)}(i+1,j)+u^{(n-1)}(i,j+1\left)\right)/4$

${\stackrel{\‾}{v}}^{(n-1)}(i,j)=\left(v^{\left(n\right)}\right(i-1,j)+v^{\left(n\right)}(i,j-1)+v^{(n-1)}(i+1,j)+v^{(n-1)}(i,j+1\left)\right)/4$

其中，u,v表示图像中像素点(i,j)点处分别沿x轴，y轴方向的光流速，是u,v在(i,j)点4邻域的平均值，n表示迭代次数，i表示像素点x轴坐标，j表示像素点y轴坐标。

3.根据权利要求1所述的一种光流场计算方法，其特征在于，所述步骤S3计算该像素点当前迭代的光流值的表达式为：

$u^{\left(n\right)}={\stackrel{\‾}{u}}^{(n-1)}-\frac{I_x(I_x{\stackrel{\‾}{u}}^{(n-1)}+I_y{\stackrel{\‾}{v}}^{(n-1)}+I_t)}{({\mathrm{\α}}^2+{I_x}^2+{I_y}^2)}$

$v^{\left(n\right)}={\stackrel{\‾}{v}}^{\left(n-1\right)}-\frac{I_y(I_x{\stackrel{\‾}{u}}^{(n-1)}+I_y{\stackrel{\‾}{v}}^{(n-1)}+I_t)}{({\mathrm{\α}}^2+{I_x}^2+{I_y}^2)}$

其中，I_x为图像灰度对x的微分，I_y为图像灰度对y的微分，I_t为图像灰度对t的微分，α表示权重值。

4.根据权利要求1所述的一种光流场计算方法，其特征在于，判断各像素点当次迭代得到的光流值与前一次迭代得到的光流值差值是否小于或等于给定容差，即：

|u⁽ⁿ⁾(i,j)-u^(n-1)(i,j)|≤ε₁

|v⁽ⁿ⁾(i,j)-v^(n-1)(i,j)|≤ε₂

其中，ε₁，ε₂为给定的估计容差，|·|表示取绝对值运算。