CN102831616B - 一种视频流体运动矢量计算方法 - Google Patents

Abstract

本发明公开了一种有效的流体运动矢量计算方法，该方法旨在利用视频流的两帧图像进行流体运动矢量的计算，目的是获得流体连续性运动的二维矢量。该方法首先统计每个像素所在区域的直方图，利用区域直方图主要成分匹配方法完成对运动矢量的初始化；然后，对初步的计算结果，利用多分辨率采样方法对其进行去噪处理；再采用统计分析方法对局部模块的运动矢量进行平滑处理；最后，对运动矢量进一步优化，得到具有保持外观一致性和运动连续性的流体运动矢量。本发明可有效地应用于流体的三维重建的应用中，具有实用的特点，能够很好地解决流体的重建问题，可以进一步应用于虚拟场景建立的应用中，克服了基于视觉的流体重建的复杂性问题。

Description

一种视频流体运动矢量计算方法

技术领域

本发明涉及一种视频流体运动矢量的计算方法，利用流体运动的连续性进行优化，使得流体运动矢量的计算结果更符合流体的运动规律。

背景技术

流体运动矢量计算是增强现实技术研究的基础，目前，在流体的物理仿真中还没有完全将动态纹理区域的识别与基于物理属性的仿真结合起来，并对流体进行实时重建；同时，目前缺乏与视频媒体真实感一致的流体建模和仿真，并且较少考虑自然景物与其它景物及环境的交互作用，这些问题产生的根本原因是流体运动矢量的计算问题没有得到解决。因此，如何从流体视频中计算出较为准确的运动矢量是研究重建与交互的关键所在，它的研究具有重要的现实意义和实用价值。

流体运动是一种典型的非刚体运动，流体运动图像的计算与分析属于动态纹理分析的技术。在物体运动的图像分析中，纹理分为弱纹理和强纹理。弱纹理具有自相似性，其图像变化具有亮度不变性；而视频中的流体部分属于强纹理，其特点是亮度保持性，即其内部的运动产生了形状的改变和遮挡，但是图像中强度的成分是不变的。由于流体具有透明性，因此可以通过计算粒子的运动矢量来间接获得流体本身的运动特征。

目前，对于流体的运动矢量的计算通常采用粒子图像测速方法实现，即通过跟踪每一个可见粒子来实现运动矢量的获取。虽然这种计算方法可以克服传统方法中的不足。但研究中存在的主要问题是由于流体粒子的运动引起的遮挡及再现，会影响运动矢量计算的准确性。在传统的变分的研究方法及粒子图像测速方法中，假设运动矢量在局部空间上是光滑的，但是这样带来了运动边界的处理及遮蔽的问题。

后来人们在研究中提出了一种基于时空多分辨率的流体速度及加速度场的计算方法，该方法利用时空多分辨率直方图对流体速度及加速度的方向进行统计分析，最后得到流体的速度及加速度场，但是如果将该方法应用于多相流体运动的场景中，特别是在视频的两帧间摄像机的旋转及平移矢量非常小的情况下，势必会影响计算的结果。对于多相流场景的直方图统计分析中，很难将得到的速度及加速度与多相流区域的运动对应起来。

人们在研究中分析了现有流体运动矢量计算的模型和方法，在目标函数中加入非局部项，并给出了启发式的方法，推导出新的目标函数，并采用邻域的中值滤波方法研究了运动矢量计算的方法，利用该方法对稀疏的粒子运动进行计算，可以得到较准确的结果，对于流体中大量粒子的情形，利用该方法进行研究，所得结果的精度受到影响。

在最近的研究中，研究人员利用光度和几何进行了块匹配，从而对流体的运动进行估计,采用了改进的H-S方法，研究中使用优化技术及多分辨率技术对流体运动矢量进行计算。该方法对于较为稀疏的流体粒子，能够取得很满意的实验结果，对于大量具有复杂物理规律的粒子运动，计算的精度会受到影响。

发明内容

本发明的目的是针对现有技术的不足而提供的一种有效的视频流体运动矢量计算方法，利用该方法能够得到具有保持外观一致性和运动连续性的视频流体运动矢量。

实现本发明目的的具体技术方案是：

一种视频流体运动矢量的计算方法，该方法包括：

（1）流体运动矢量的初始化

为了求得初始的运动矢量，首先对于运动前后两幅图像分别统计其直方图，每幅图像按照一定的邻域进行统计。在视频的前后两帧中，摄像机的旋转及平移的矢量较小，采用a×a的邻域进行，将统计的结果记为中心像素邻域的统计结果。把图像的强度分为不同的级别，按照不同级别统计出区域中像素的个数。若当前的图像中像素的灰度值为density，该像素所在的级别计算为：

其中，表示取整。对于灰度图像，density表示0～255之间的一个整数值。认定直方图中像素较大数目的级别为统计直方图的主要成分。

为了求取邻域强度的主要成分，对直方图的统计结果进行排序，同时将对应直方图统计的成分进行排序，利用区域的直方图主要成分进行匹配操作。

假定在一个小区域范围内所有粒子具有相近的运动矢量，在此基础上，寻找一个运动矢量使得区域范围内所有粒子运动前后所对应的区域具有主要成分相似的特性。

假设运动前后的两幅图像为I_k(k=1,2)，中心像素I_k(i,j)在邻域内排序的统计直方图为HistogramI_k(i,j)，对应排序后的成分为HistogramxI_k(i,j)。在a×a区域范围内找到一个矢量使得在运动前后，对应邻域内所有像素的直方图主要成分的欧氏距离最小。运动矢量初始化算法如下：

1)分别计算两幅图像I_k(k=1,2)的像素点pixel_k(i,j)的a×a邻域的统计直方图，记为Histogram_k(i,j)；

2)对每个像素点pixel_k(i,j)的邻域直方图的统计结果，按照从大到小排序，直方图的成分序号也随之排序，得到排序后的统计直方图为HistogramI_k(i,j)及排序的成分HistogramxI_k(i,j)；

3)在a×a邻域内搜索一个矢量使得对应邻域内所有像素的直方图主成分的欧氏距离d(i,j)最小，其中d(i,j)为：

$d(i,j)=\underset{p=-2}{\overset{2}{\mathrm{\Σ}}}\underset{q=-2}{\overset{2}{\mathrm{\Σ}}}\underset{h=1}{\overset{w}{\mathrm{\Σ}}}{({\mathrm{HistogramxI}}_1(i+p,j+q,h)-{\mathrm{HistogramxI}}_2(i+p+f,j+q+s,h))}^2---\left(2\right)$

其中HistogramxI_k(i,j,h)表示第k幅图像，像素点pixel(i,j)邻域直方图排序后的主要成分的第h个分量；

得到当前像素pixel_k(i,j)的运动矢量初始化结束。

（2）利用多分辨率法去除噪声

为了去除初始化过程中产生的噪声，进行去噪的处理，采用多分辨率方法去除噪声，即可以将初始化的结果经过几次降低分辨率采样，再在较低分辨率的采样结果中，利用阈值法对统计结果进行处理，其具体步骤如下：

1）算法初始化。计数器m=1，设定采样窗口为WinSize×WinSize；

对于初始化的运动矢量进行第m次采样；

2）统计WinSize×WinSize区域中运动矢量的数目为numpixel，并采用阈值方法进行判断，如果满足（3）式，则认为该像素的运动矢量初始化结果计算是无效的，作为噪声去除，并得到初步的去噪结果；

numpixel<threhold_smooth                                （3）

3）判断是否得到了满意的去噪结果？如果是，则转步骤5），否则，进行下一步；

4）m=m+1，重新设置采样窗口大小，并转步骤2）；

5）结束，得到去噪后运动矢量的结果。

去除噪声是可选的步骤，只有在场景中存在大片静止的背景区，需要进行去噪时，才执行去除噪声步骤。

（3）运动矢量的平滑

运动矢量平滑算法的步骤如下：

1）初始化

对于每个象限的累加器Accumulator_k(k=1...4)初始化为0值；

2）对于像素的邻域中各像素的运动矢量，按方向分别累计，记累计矢量为并使Accumulator_k（1≤k≤4)增值；

3）选择计数最大的累加器Accumulator_m（1≤m≤4)，并判断Accumulator_m的大小，如果满足（4）式，那么，当前的像素粒子的运动矢量置为零矢量；否则，其运动矢量用(5)式计算。

Accumulator_m<threhold_smooth                                （4）

$\stackrel{\&RightArrow;}{\mathrm{vec}}=\stackrel{\&RightArrow;}{{\mathrm{vector}}_k}/{\mathrm{Accumulator}}_m---\left(5\right)$

4）判断是否所有像素都处理完毕了。如果是，转步骤5）；如果否，转步骤1）。

5）结束，得到运动矢量平滑的结果。

（4）基于流体运动连续性的优化计算

采用流体运动的连续性方程进行优化计算，这样可以体现出流体的运动特征。流体的连续性运动方程为：

${\&Integral;}_{\mathrm{\&Omega;}}{\&PartialD;}_t\mathrm{IdA}+{\&Integral;}_{\&PartialD;\mathrm{\&Omega;}}\mathrm{If}\&CenterDot;\mathrm{ndL}=0---\left(6\right)$

其中，是图像I的强度对时间的导数，表示Ω的边界，n是的外部矢量。根据散度理论及数学推导，(6)式可以写成下面的方程形式：

I_t+uI_x+vI_y+Iu_x+Iv_y=0                                        （7）

式中I_t表示图像强度对时间的导数，I_x及I_y分别表示图像强度在空间域上的导数，流体的运动矢量为(u,v)，u_x和v_y分别表示运动矢量的偏导数，采用u和v的一阶向前差分进行计算。采用连续性方程（7）对运动矢量进行优化计算。

基于流体运动连续性优化计算的步骤为：

1）初始化。对于图像的每个a×a区域，对满足（7）式的像素进行标记flag(i,j)=1(1≤i≤nWidth,1≤j≤nHeight)；

2）对于没进行标记的像素X，优化计算X像素粒子的运动矢量中的u_X分量。对于任意一个像素X，如果其邻域中像素A和B已经被标记，利用式（8），可以重新计算出X的运动矢量中的u_X分量；

I_At+u_A(I_X-I_A）+v_A(I_B-I_A)+I_A(u_X-u_A)+I_A(v_B-v_A)=0            (8)

3）对于没进行标记的像素X，优化X像素粒子的运动矢量中的v_X分量。对于任意一个像素X，如果其邻域中像素A和B已经被标记，利用运动连续性，重新计算出X的运动矢量中v_X的分量；

4）更新标记。如果像素X的u_X以及v_X都已经进行了优化，则对X进行标记；

5）判断是否满足了终止条件？如果是，则进行下一步；否则，转步骤2）；

6）得到流体运动矢量的优化计算结果。

本发明具有简单、准确的特点，计算所得到的运动矢量既具有保持运动前后外观的一致性，也具有保持流体运动连续性的特点，即使在流体相邻两帧的运动很小的情况下，也能得到满意的计算结果，与现有方法相比，它具有更准确的特点；此外利用该算法对于流体的高度进行重建，能够取得满意的结果，这进一步说明本算法可以有效地应用于流体的三维重建的应用中，因此具有实用的特点，能够很好地解决流体的重建问题，可以进一步应用于虚拟场景建立的应用中，克服了基于视觉的流体重建的复杂性问题。

附图说明

图1为本发明实施例“6482910”的第198帧运动矢量的初始化结果图；

图2为本发明实施例“6482910”的第198帧运动矢量的平滑结果图；

图3为本发明实施例“6482910”的第198帧运动矢量的优化结果图。

具体实施方式

实施例

下面结合附图对于发明进一步说明。

本实例采用DynTex动态纹理库中的“6482910”的第198帧和199帧进行流体运动矢量的计算。在PC机上Windows7操作系统下进行，其硬件配置是2.66GHz Intel Core（TM）2Duo CPU、4GB RAM。

（1）运动矢量的初始化

假定在一个小区域范围内所有粒子具有相近的运动矢量，在此基础上，寻找一个运动矢量使得区域范围内所有粒子运动前后所对应的区域具有主要成分相似的特性。

运动前后的两幅图像为I_k(k=1,2)，中心像素I_k(i,j)在邻域内排序的统计直方图为HistogramI_k(i,j)，对应排序后的成分为HistogramxI_k(i,j)。由于视频中相继两帧的摄像机位置参数变化较小，流体的运动矢量很小，因此在a×a区域范围内找到一个矢量使得在运动前后，对应邻域内所有像素的直方图主要成分的欧氏距离最小。这里a取5，如果流体的运动速度较快，可以根据需要适当扩大运动矢量的搜索范围。

首先利用下面的算法进行运动矢量的初始化。

1）分别计算两幅图像I_k(k=1,2)的像素点pixel_k(i,j)的5×5邻域的统计直方图，记为Histogram_k(i,j)；

2）对每个像素点pixel_k(i,j)的邻域直方图的统计结果，按照从大到小排序，直方图的成分序号也随之排序，得到排序后的统计直方图为HistogramI_k(i,j)及排序的成分HistogramxI_k(i,j)；

3）在5×5邻域内搜索一个矢量使得对应邻域内所有像素的直方图主成分的欧氏距离d(i,j)最小，其中d(i,j)为：

$d(i,j)=\underset{p=-2}{\overset{2}{\mathrm{\&Sigma;}}}\underset{q=-2}{\overset{2}{\mathrm{\&Sigma;}}}\underset{h=1}{\overset{w}{\mathrm{\&Sigma;}}}{({\mathrm{HistogramxI}}_1(i+p,j+q,h)-{\mathrm{HistogramxI}}_2(i+p+f,j+q+s,h))}^2$

其中HistogramxI_k(i,j,h)表示第k幅图像，像素点pixel(i,j)邻域直方图排序后的主要成分的第h个分量；

得到当前像素pixel_k(i,j)的运动矢量算法结束。

本发明实例中利用该算法进行运动矢量初始化时，利用了邻域直方图排序后的5个主要成分进行实验，即算法中w=5。图1是“6482910”的第198帧运动矢量的初始化结果。

由于该实施例中，运动矢量的噪声较少，因此就不必要进行去噪处理了。

（2）运动矢量的平滑

为了提高计算的准确性，对初始化的结果进行平滑处理，采用局部平滑的方式。将图像中每个像素周围的5×5区域中的像素进行了统计和分析，将其中的25个像素的运动矢量，按照象限进行分类统计，并利用统计数目最大的象限中运动矢量的平均值，作为该区域的中心像素的运动矢量。

运动矢量平滑算法的步骤如下：

1)初始化。对于每个象限的累加器Accumulator_k(k=1...4)初始化为0值；

2）对于像素的邻域中各像素的运动矢量，按方向分别累计，记累计矢量为并使Accumulator_k（1≤k≤4)增值；

3）选择计数最大的累加器Accumulator_m（1≤m≤4)，并判断Accumulator_m的大小，如果满足（4）式，那么，当前的像素粒子的运动矢量置为零矢量；否则，其运动矢量用(5)式计算。

Accumulator_m<threhold_smooth                                （4）

$\stackrel{\&RightArrow;}{\mathrm{vec}}=\stackrel{\&RightArrow;}{{\mathrm{vector}}_k}/{\mathrm{Accumulator}}_m---\left(5\right)$

4）判断是否所有像素都处理完毕了。如果是，转步骤5）；如果否，转步骤1）。

5）算法结束，得到运动矢量平滑的结果。

在该实例中用threhold_smooth取值为5，得到的运动矢量平滑结果，图2是“6482910”的第198帧运动矢量的平滑结果，从图中可以看出，平滑效果比较理想。

流体运动矢量计算结果的准确与否，可以从两个方面进行衡量。第一，是从外观的一致性进行衡量；第二，是流体运动的估计要符合运动的连续性。为了满足这两个要求，并保证流体运动矢量计算的准确性，采用流体运动的连续性方程进行优化，这样可以体现出流体的运动特征。

（3）基于流体运动连续性的优化计算

采用流体运动的连续性方程进行优化，这样可以体现出流体的运动特征。流体的连续性运动方程为：

${\&Integral;}_{\mathrm{\&Omega;}}{\&PartialD;}_t\mathrm{IdA}+{\&Integral;}_{\&PartialD;\mathrm{\&Omega;}}\mathrm{If}\&CenterDot;\mathrm{ndL}=0---\left(6\right)$

其中，是图像I的强度对时间的导数，表示Ω的边界，n是的外部矢量。根据散度理论及数学推导，(6)式可以写成下面的方程形式：

I_t+uI_x+vI_y+Iu_x+Iv_y=0                                            （7）

式中I_t表示图像强度对时间的导数，I_x及I_y分别表示图像强度在空间域上的导数，流体的运动矢量为(u,v)，u_x和v_y分别表示运动矢量的偏导数，在研究中采用u和v的一阶向前差分进行计算。研究中采用连续性方程（7）对运动矢量进行优化。

基于流体运动连续性优化算法的步骤为：

1）初始化。对于图像的每个3×3区域，对满足（7）式的像素进行标记flag(i,j)=1(1≤i≤nWidth,1≤j≤nHeight)；

2）对于没进行标记的像素X，优化X像素粒子的运动矢量中的u_X分量。对于任意一个像素X，如果其邻域中像素A和B已经被标记，利用式（8），可以重新计算出X的运动矢量中的u_X分量；

I_At+u_A(I_X-I_A）+v_A(I_B-I_A)+I_A(u_X-u_A)+I_A(v_B-v_A)=0        (8)

3）对于没进行标记的像素X，优化X像素粒子的运动矢量中的v_X分量。对于任意一个像素X，如果其邻域中像素A和B已经被标记，利用运动连续性，重新计算出X的运动矢量中v_X的分量；

4）更新标记。如果像素X的u_X以及v_X都已经进行了优化，则对X进行标记；

5）判断是否满足了终止条件？如果是，则进行下一步；否则，转步骤2；

6）得到流体运动矢量的优化结果，算法结束。

计算中的图像尺度是nWidth×nHeight。将被标记的像素数目达到像素总数目作为算法的终止条件。

将本发明与现有方法进行比较：

为了进行流体运动矢量计算结果的准确性比较。将本发明的误差与现有方法的实验误差进行比较。实验中计算了两种误差：第一种是运算矢量的方向误差Errangle；另一种是运动矢量大小的误差Erruv。其中方向误差Errangle是衡量运动矢量的计算结果与标准矢量之间方向的差异，它们的计算采用式(9)和（10）进行计算。

$\mathrm{Errangle}=\frac{1}{n}\underset{i=1}{\overset{\mathrm{iWidth}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{\mathrm{iHeight}}{\mathrm{\&Sigma;}}}\left|\arccos ({\mathrm{vec}}_{\mathrm{ij}},{\mathrm{vecstd}}_{\mathrm{ij}})\right|---\left(9\right)$

$\mathrm{Erruv}=\frac{1}{n}\underset{i=1}{\overset{\mathrm{iWidth}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{\mathrm{iHeight}}{\mathrm{\&Sigma;}}}\sqrt{[{\left(\right|{\mathrm{vec}}_i|-|{\mathrm{vecstd}}_i\left|\right)}^2+{\left(\right|{\mathrm{vec}}_j|-|{\mathrm{vecstd}}_j\left|\right)}^2]}---\left(10\right)$

式中vec_ij是对图像粒子P(i,j)计算得到的运动矢量的结果，vecstd_ij（vecstdi_,vecstd_j）是P(i,j)的准确运动矢量，iWidth和iHeight分别为图像的宽和高，n表示该图像像素的总数目，即n=iHeight×iWidth；vec_i和vec_j分别表示运动矢量计算结果在i和j方向上的分量，即图像粒子p[i,j]的运动矢量计算结果为vec(vec_i,vec_j)。

表1Yosemite实验结果的误差比较

表2对OpenGL生成的流体实验结果的误差比较

从误差的计算结果中可以明显地看出，对于较复杂流体运动矢量的计算，三种方法计算得到的运动矢量方向误差Errangle及大小误差Erruv略有区别，本发明的误差略较小；对于运动情形较为简单的合成流体运动来说，采用本发明得到结果的方向误差及大小误差明显偏小，这些都说明新发明计算结果更准确。

Claims (1)

1.一种视频流体运动矢量计算方法，其特征在于该方法包括以下步骤：

⑴、运动矢量的初始化

分别计算运动前后两帧图像中每个像素所在区域的灰度直方图，利用直方图的主要成分匹配方法得到初始的运动矢量；

⑵、去噪处理

采用多分辨率法对初始化的结果进行去除噪声处理；

⑶、平滑处理

采用统计分析方法对局部区域的运动矢量进行平滑处理；

⑷、优化处理

用流体运动的连续性方程对运动矢量进一步计算，得到具有保持外观一致性和运动连续性的流体运动矢量；

其中：所述直方图的主要成分匹配方法是：对运动前后两帧图像分别进行区域直方图的统计，然后对两帧图像的每个邻域求取其统计后的主要成分，利用对应区域的主要成分匹配误差最小来完成对运动矢量的初始化，按照式(1)计算像素所在的直方图统计类别

其中，表示取整，对于灰度图像，density表示0～255之间的一个灰度值；

相邻两帧中像素pixel_k(i,j)运动矢量的初始化按照如下步骤进行：

a)分别计算视频相邻两帧I_k，k＝1,2中像素点pixel_k(i,j)的a×a邻域的统计直方图，记为Histogram_k(i,j)；

b)对每个像素点pixel_k(i,j)的邻域直方图的统计结果，按照从大到小进行排序，直方图的成分序号也随之排序，得到排序后的统计直方图为HistogramI_k(i,j)及排序的成分HistogramxI_k(i,j)；

c)在pixel_k(i,j)的a×a邻域内搜索一个矢量使得对应邻域内所有像素的直方图主成分的欧氏距离d(i,j)最小，其中d(i,j)为：

$d(i,j)=\underset{p=-a}{\overset{a}{\mathrm{\&Sigma;}}}\underset{q=-a}{\overset{a}{\mathrm{\&Sigma;}}}\underset{h=1}{\overset{w}{\mathrm{\&Sigma;}}}{(\mathrm{Histogramx}{I}_1(i+p,j+q,h)-\mathrm{Histogramx}{I}_2(i+p+f,j+q+s,h))}^2---\left(2\right)$

其中HistogramxI_k(i,j,h)表示第k幅图像，像素点pixel(i,j)邻域直方图排序后的主要成分的第h个分量；在视频中相邻帧的运动较小，因此区域的范围a取值较小；

d)得到当前像素pixel_k(i,j)的运动矢量

e)得到流体运动矢量初始化结果；

步骤(4)所述的优化计算的步骤包括：

(a)初始化

对于图像的每个a×a区域，对满足(3)式的像素进行标记

flag(i,j)＝1，1≤i≤nWidth,1≤j≤nHeight；

I_t+uI_X+vI_y+Iu_x+Iv_y＝0    (3)

式中I_t表示图像强度对时间的导数，I_X及I_y分别表示图像强度在空间域上的导数，流体的运动矢量为(u,v)，u_x和v_y分别表示运动矢量的偏导数，采用u和v的一阶向前差分进行计算；

(b)对于没进行标记的像素X，如果其邻像素已经被标记，利用式(4)，重新计算出X的运动矢量中的u_X分量；

I_At+u_A(I_X-I_A)+v_A(I_B-I_A)+I_A(u_X-u_A)+I_A(v_B-v_A)＝0    (4)

(c)对于没进行标记的像素X，优化X像素粒子的运动矢量中的v_X分量，并重新计算出X的运动矢量中v_X的分量；

(d)更新标记

如果像素X的u_X以及v_X都已经进行了计算，则对X进行标记；

(e)判断是否满足了终止条件，如果是，则进行下一步；否则，转至步骤(4)中的步骤(b)；

(f)得到流体运动矢量的计算结果

计算终止条件是被标记的像素数目达到像素总数目，如果满足，则终止计算过程。