CN102819662B - 一种视频流体高度的计算方法 - Google Patents

Abstract

本发明公开了一种视频流体高度的计算方法，该方法首先根据流体运动特征，对流体的运动矢量进行初始化；根据运动矢量的结果，利用LBM计算流体粒子的高度；再对高度进行去噪和平滑处理；然后，根据流体运动的连续性，利用LBM递推连续多帧的高度结果，得到具有流体物理运动特性的高度计算结果；再利用流体粒子分布函数线性插值的方法，得到连续运动流体的均匀变化的高度结果；最后，采用中间帧对流体高度进行校正。本发明具有简单、方便、快速的特点，利用本发明可得到连续变化的流体高度场，并且产生实时的高度结果，能够应用于自然景观虚拟场景的构建，并且可有效地实现虚拟场景中的流体运动实时交互，克服基于视觉方法中重建的实时性差的问题，具有一定的应用价值。

Description

一种视频流体高度的计算方法

技术领域

本发明涉及视频流体高度的计算方法。

背景技术

目前，在增强现实研究中，对于刚体虚实结合的研究已经取得了初步的成果，而对于含有水、烟雾、云、火等流体的场景，由于流体属于强纹理，在运动过程中存在着遮挡及再现的问题，所以有关流体的增强现实技术的研究存在着一定挑战性，在这一方面研究中取得的成果还不够显著。流体重建技术是研究中的关键问题，如何快速、准确地计算流体高度，建立具有真实感和交互性的流体自然景观的增强现实场景，仍是该领域中极具挑战性的课题，其研究具有重要的现实意义和应用价值。

随着图像处理及计算机视觉技术的发展，对于自然景观场景高度计算的研究出现了一些技术和方法。人们利用图像中的颜色、纹理、形状、运动等特征进行研究，取得了一些成果。现有方法从图像纹理计算场景高度，通过纹理基元的尺寸、形状、和密度等变化特性来得到场景的高度值。

由于流体属于强纹理，在运动过程中不能保持纹理不变的属性，因此利用纹理信息准确地计算出流体的高度较为困难；有些研究是利用颜色不变性规则，通过局部区域的强度相关性进行研究，计算出三维场景的高度信息，从而实现对场景的重建；现有研究中一些研究从不同视图计算时空连续的运动信息，得到两帧图像的对应点，再结合摄像机标定、三维坐标恢复等技术，获取了场景的高度信息；另一些研究中通过提取运动的物理特性，恢复并校准摄像机内外参数，从而计算出物体的高度。基于运动特征的高度计算方法中存在的主要问题是：物体运动特征的确定比较困难；摄像机参数的确定通常采用离线的批处理方法，这些方法存在复杂度高、计算量大的问题，很难利用在线方法完成，因此不能满足实时性需求。

对于流体的重建技术人们曾经进行过系统的探讨，这方面研究的典型成果就是基于明暗度恢复形状（SFS)的方法。它是利用不同照明条件下的物体表面强度变化来估算物体表面的法线方向，从而达到重建的目的。

近年来，为了提高重建的准确性和真实感，人们进行了一些研究。现有研究包括利用三次样条与SFS结合的方法计算场景的高度信息、对求解SFS问题的算法进行研究；人们也利用SFS方法对高度计算的结果进行平滑处理，得到了较好的重建结果。

一些人员利用SFS和流体的运动矢量进行研究，基于质量守恒的原则实现了流体的高度计算，这种方法对于非透明水域能够取得较满意的高度计算结果，但强光照射的情况下，例如含有倒影或强亮度的水面，高度计算结果的准确性就会受到影响。

发明内容

为解决现有技术研究中存在的上述问题，本发明的目的是提出一种有效的视频流体高度的计算方法，利用该方法能够得到反应流体连续运动特征的流体运动的高度场，快速、准确地计算出流体高度。

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

一种视频流体高度的计算方法，该算法分为五个步骤：

1）流体运动矢量的计算；

2）运动矢量和LBM相结合的高度计算；

3）对高度计算结果的去噪和平滑；

4）用流体运动的连续性方程进行约束。

5）高度校正。

本发明中具体的技术细节如下：

（1）流体运动矢量的计算

由于流体粒子所在的邻域在运动前后具有保持强度成分不变的特性。根据这一特性，采用区域强度相关性对流体运动矢量进行初始化。再进行聚类。

定义用于聚类的特征向量，其形式为vector=[x,y,u,v,sigu,sigv]，其中(x,y)为流体粒子所在的位置，u和v分别表示该粒子的运动矢量沿x方向和沿y方向的分量，sigu和sigv分别表示u和v的符号，如果为正取1，如果为负取-1，否则取0。所定义的特征向量进行聚类，并根据聚类结果的规模进行优化，最后为了得到密集的运动矢量场，采用最短距离线性插值的方法，对于每个粒子，找到距离它最近的粒子，并采用线性插值法进一步得到密集的运动矢量结果。

流体运动矢量的计算具体包括：

1）利用连续两帧图像，采用区域相关运算方式对流体运动矢量进行初始化。

2）利用运动矢量的特征向量vector=[x,y,u,v,sigu，sigv]进行聚类，具体地，如果粒子P_i和P_j的运动矢量的特征向量分别为：vector_i=[x_pi,y_pi,u_pi，v_pi,sigu_pi,sigv_pi]及vector_j=[x_pj,y_pj,u_pj,v_pj,sigu_pj,sigv_pj]，并且，它们同时满足式（1）至（4），则认为粒子P_i和P_j的运动具有相似性，并将它们聚到相同的类别之中。

$|\sqrt{{(x_{\mathrm{pi}}-x_{\mathrm{pj}})}^2+{(y_{\mathrm{pi}}-y_{\mathrm{pj}})}^2}<{\mathrm{threshold}}_{\mathrm{dis}}---\left(1\right)$

$\left|\sqrt{\mathrm{angle}({\mathrm{vec}}_{\mathrm{pi}},{\mathrm{vec}}_{\mathrm{pj}})}\right|<{\mathrm{threshold}}_{\mathrm{angle}}---\left(2\right)$

sigu_pi*sigu_pj≥0      （3）

sigv_pi*sigv_pj≥0      （4）

其中，threshold_dis表示粒子P_i和P_j之间距离的阈值；vec_pi和vec_pj分别表示粒子P_i和P_j的运动矢量(u_pi,v_pi)及(u_pj,v_pj)；threshold_angle表示它们运动矢量方向的阈值。

3）统计各类别的规模，对于聚类结果中规模较大的类别，保留其结果。设任意一个类D，它的类内粒子数为D_num，即该类的规模为D_num，如果其满足式（5）则认为类内粒子运动矢量是较为准确的，并将这些粒子的flag值置为1；否则，对于粒子的flag值置为0；

D_num>threshold_D     （5）

其中，threshold_D是粒子数目的阈值。

4）对于任意一个flag值为1的粒子P_m，在其周围的区域中选取一个与之最近且flag值为1的粒子P_n,对这两个粒子连线上的所有flag值为0的粒子采用线性插值的方法，得到它们的运动矢量结果。

（2）运动矢量和LBM相结合的高度计算

LBM是一种简化的流体运动的微观模型，它将大量流体粒子的微观热运动的统计结果表示为粒子分布函数的形式，通过分布函数的演化，进一步反映流体的宏观运动。

定义标准方向。所谓标准方向是中心粒子受到来自周围8个粒子力的作用方向，这些力对中心粒子在三维空间的高度演化起着作用。

如周围粒子的运动矢量方向与某个标准方向相同，则它在该方向上对中心粒子的作用力较大。根据这一特点，将周围粒子的运动矢量投影到标准方向上，确定周围粒子对中心粒子的作用。当周围粒子的运动方向的分量与对应标准方向的夹角大于等于90度时，则在该标准方向上对中心粒子产生的作用力为零。

流体粒子在三维空间中的高度为f_i(i＝1,...,8)之和，和越大则粒子在三维空间中高度值越大，反之越小。

运动矢量和LBM相结合的高度计算具体包括：

1）假设任意一个尚未计算高度的粒子S(x,y)，其周围粒子为T_i(x,y)(i＝1,...,8)，它们对应的运动矢量为(u_i,v_i)。粒子T_i(x,y)对粒子S(x,y)作用力的标准方向为F_i。设标准方向F_i与向量Hor(u_i，0)的夹角为θ₁，F_i与向量Ver(0，v_i)的夹角为θ₂。用式(6)计算粒子S(x，y)的分布函数f_i(x，y，t)。

2)利用式(7)计算粒子S(x，y)在三维空间中的高度h。

$h(r,t)=\underset{i=0}{\overset{d}{\mathrm{\&Sigma;}}}{f}_i(r,t)---\left(7\right)$

3)判断是否所有粒子的高度已经计算，如果是，转步骤4)，否则，转步骤1)。

4)算法结束。

(3)对高度计算结果的去噪和平滑

先利用高度变化的连续性，对高度计算的结果进行去噪，然后在局部区域内进行平滑处理。其具体步骤如下：

1)对任意一个尚未处理的粒子W(x，y)，利用式(8)和(9)计算其高度h(x，y)的二阶差分，并记为及

$\frac{{\&PartialD;}^{2}h}{\&PartialD;x}=h(x+1,y)-2*h(x,y)+h(x-1,y)---\left(8\right)$

$\frac{{\&PartialD;}^{2}h}{\&PartialD;y}=h(x,y+1)-2*h(x,y)+h(x,y-1)---\left(9\right)$

2)采用式(10)和式(11)对和进行判断。如果同时满足(10)和式(11)，置粒子W(x，y)的flagh值为1；否则，将粒子W(x，y)的高度h(x，y)和flagh值置0。

$\left|\frac{{\&PartialD;}^{2}h}{\&PartialD;x}\right|<\mathrm{threshold}---\left(10\right)$

$\left|\frac{{\&PartialD;}^{2}h}{\&PartialD;y}\right|<\mathrm{threshold}---\left(11\right)$

3)判断是否所有的粒子都进行了处理，如果是，转步骤4)；否则，转步骤1)。

4)对任意一个尚未平滑处理的粒子W(x，y)，取以它为中心s×s区域，用公式(12)对高度进行平滑。

$h(x,y)=\underset{i=1}{\overset{s\&times;s}{\mathrm{\&Sigma;}}}(\mathrm{flagh}\left(i\right)\&times;h\left(i\right))/\underset{i=1}{\overset{s\&times;s}{\mathrm{\&Sigma;}}}\mathrm{flagh}\left(i\right)---\left(12\right)$

其中，flagh(i)和h(i)分别为s×s邻域中第i个粒子的flagh值和高度。

5）判断是否所有的粒子都进行了平滑处理，如果是，转步骤6）；否则，转步骤4）。

6）算法结束。

（4）利用LBM递推流体的高度

为了满足流体实时重建的实际需要，采用递推高度的方法产生中间连续帧的高度值。采用LBM方法，利用粒子分布函数f_i来实现碰撞和平流的过程。

利用LBM方法的高度计算包括平流、碰撞和边界处理的过程。具体为：

1）在平流的过程中，在运动矢量初始化以后，可以得到当前帧的高度，然后利用公式(13)计算下一帧的各个方向上的粒子分布函数。

f_i(r+e_iδ_t,t+δ_t)=f_i(r,t)      （13）

其中，f_i(r,t)是t时刻位于r处的粒子分布函数。

2）在碰撞的过程中，利用公式（14）计算粒子分布函数，并利用公式(15)计算平衡态分布函数。

$f_i(r+e_i{\mathrm{\&delta;}}_t,t+{\mathrm{\&delta;}}_t)-f_i(r,t)=-\frac{1}{\mathrm{\&tau;}}[f_i(r,t)-f_i^{\mathrm{eq}}(r,t)]---\left(14\right)$

$f_i^{\mathrm{eq}}={\mathrm{\&omega;}}_i\mathrm{\&rho;}[1+\frac{e_i\&CenterDot;c}{e_s^2}+\frac{{(e_i\&CenterDot;c)}^2}{2e_s^4}-\frac{c^2}{2e_s^2}]---\left(15\right)$

其中，为沿i方向上的平衡态分布函数；ω_i表示沿i方向的权系数；ρ表示流体的宏观密度；e_s为常数；e_i和c的意义与前述相同。

利用公式（7）即可计算出流体高度。

通过周期边界方法对边界进行处理。通常将式(15)中的参数取为：

$e_i=\left[\begin{array}{ccccccccc}0& 1& 0& -1& 0& 1& -1& -1& 1\\ 0& 0& 1& 0& -1& 1& 1& -1& -1\end{array}\right]$

${\mathrm{\&omega;}}_i=\left\{\begin{array}{c}4/9,e_i^2=0\\ 1/9,e_i^2=1\\ 1/36,e_i^2=2\end{array}\right.$ $e_s=\frac{1}{\sqrt{3}}$

（5）高度校正

假设第T帧及其之前的高度已经由LBM递推产生，为了减少递推过程中的积累误差，在第K帧处对流体高度进行校正。第K帧的高度值可以由其运动矢量计算得到，较为准确。对于(T，K)之间的各帧的高度，采用分布函数的插值方法计算各周围粒子的分布函数，进一步用分布函数求和的方法来计算中心粒子的高度信息。流体高度校正的具体步骤如下：

1）对第T至第K帧中任意尚未校正高度的帧X，设在(x，y)处粒子i方向的粒子分布函数分别为f_Ti(x,y,t)和f_Ki(x,y,t)，如果它们之间相差过大，即不满足公式（16）的条件，则可能产生流体高度突变的情况，为此按照公式（17）调整f_Ki的大小。转步骤2）。

$\left|\frac{f_{\mathrm{Ki}}-f_{\mathrm{Ti}}}{f_{\mathrm{Ti}}}\right|\&le;\mathrm{factor}---\left(16\right)$

f_Ki=f_Ti×(1+sig_factor×factor)     （17）

公式（16）和公式（17）中，factor用来衡量粒子分布函数的增长速度。sig_factor为f_Ki-f_Ti的符号，如果为正取1，如果为负取-1，否则取0。

2）第X(X∈(T,K))帧的粒子分布函数利用公式（18）线性插值得到。

f_Xi=k×(f_Ki-f_Ti)/(K-T-1)       （18）

其中，k（k=1，2……，K-T-1），f_Xi表示第X帧的i方向的粒子分布函数。

3）判断是否存在尚未进行高度校正的帧，如果存在，转步骤1）；否则转步骤4）。

4）算法结束。

利用式（19）计算factor。

$\mathrm{factor}=\frac{K-T-1}{f_{\mathrm{Ti}}+f_{\mathrm{Ki}}}---\left(19\right)$

经过高度校正后，得到了与视频流体真实感相一致的高度计算结果。

运动矢量和LBM相结合的高度计算方法的步骤：

为了实时计算具有保持流体物理运动特性和流体连续性的流体高度，获得与视频真实感一致的流体高度结果，将视频每m+1帧作为一个周期。假设对于视频序列中的m+1帧Seg_i(i＝k0,k1,...km)，首先计算Seg_k0和Seg_km帧的运动矢量，利用运动矢量的结果并结合LBM进行高度的计算，这m+1帧高度计算的步骤如下：

1）流体运动矢量的初始化。对Seg_k0帧和Seg_km帧分别计算密集的运动矢量；

2）运动矢量和LBM相结合的高度计算。利用运动矢量和LBM结合的高度计算方法，先计算出帧Seg_k0和帧Seg_km的粒子分布函数和进一步计算这两帧的流体高度h^(k0)(x，y，t)和h^(km)(x，y，t)；

3)对高度计算结果的去噪和平滑。利用流体运动连续性，对计算得到的流体高度h^(k0)(x，y，t)和h^(km)(x，y，t)进行去噪和平滑；

4）利用LBM递推流体的高度。以作为初始分布函数，计算出连续num(num∈[k0，km])帧的高度结果。

5)高度校正。对高度结果进行校准，从而得到了较为准确的高度计算结果。

本发明在对剧烈运动的海洋区域计算时，高度计算结果能明显区分出较高的浪和平静的海面的不同运动特征；对运动平缓的水域进行高度计算时，流体整体运动平缓，计算出的流体高度变化不大，计算结果仍保留了波浪的细节，能够区分波峰和波谷。运动矢量的高度计算结果具有保持流体运动的特征，计算出的高度与实际流体运动时的高度相一致。利用本发明的算法可以得到连续变化的流体高度场，并且产生实时的高度结果，可以应用于自然景观虚拟场景的构建，并且利用该方法可以有效地实现虚拟场景中的流体运动实时交互，可以克服基于视觉方法中重建的实时性较差的问题，因此本发明具有一定的应用价值。

附图说明

图1是粒子对中心粒子作用力的标准方向示意图；

图2是本发明实施例“6482810”的第61帧运动矢量的初始化结果图。

具体实施方式

实施例

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

本实施例采用DynTex动态纹理库中的“6482810”的第61帧进行流体的高度进行计算。在PC机上Windows XP操作系统下进行，其硬件配置是2.0GHz Intel Core(TM)2 Duo CPU、2GB内存。对视频每40帧作为一个周期，利用LBM递推20帧的流体高度，并进行了中间帧的矫正。具体步骤如下：

(1)流体运动矢量的计算

由于流体粒子所在的邻域在运动前后具有保持强度成分不变的特性。根据这一特性，采用区域强度相关性对流体运动矢量进行初始化。再进行聚类，定义用于聚类的特征向量，其形式为vector＝[x，y，u，v，sigu，sigv]，其中(x，y)为流体粒子所在的位置，u和v分别表示该粒子的运动矢量沿x方向和沿y方向的分量，sigu和sigv分别表示u和v的符号，如果为正取1，如果为负取-1，否则取0。所定义的特征向量进行聚类分析，并根据聚类结果的规模进行优化，最后为了得到密集的运动矢量场，采用最短距离线性插值的方法，对于每个粒子，找到距离它最近的粒子，并采用线性插值法进一步得到密集的运动矢量结果，具体为：

1）利用“6482810”的第61帧和62帧图像，采用区域相关运算对流体运动矢量进行初始化。

2）利用运动矢量的特征向量vector＝[x,y,u,v,sigu，sigv]进行聚类，具体地，如果粒子P_i和P_j的运动矢量的特征向量分别为：vector_i=[x_pi,y_pi,u_pi,v_pi,sigu_pi,sigv_pi]及vector_j=[x_pj,y_pj,u_pj,v_pj,sigu_pj,sigv_pj]，并且，它们同时满足式（1）至（4），则我们认为粒子P_i和P_j的运动具有相似性，并将它们聚到相同的类别之中。

$|\sqrt{{(x_{\mathrm{pi}}-x_{\mathrm{pj}})}^2+{(y_{\mathrm{pi}}-y_{\mathrm{pj}})}^2}<{\mathrm{threshold}}_{\mathrm{dis}}---\left(1\right)$

$\left|\sqrt{\mathrm{angle}({\mathrm{vec}}_{\mathrm{pi}},{\mathrm{vec}}_{\mathrm{pj}})}\right|<{\mathrm{threshold}}_{\mathrm{angle}}---\left(2\right)$

sigu_pi*sigu_pj≥0       （3）

sigv_pi*sigv_pj≥0       （4）

其中，threshold_dis表示粒子P_i和P_j之间距离的阈值；vec_pi和vec_pj分别表示粒子P_i和P_j的运动矢量(u_pi，v_pi)及(u_pj,v_pj)；threshold_angle表示它们运动矢量方向的阈值。

3）统计各类别的规模，对于聚类结果中规模较大的类别，保留其结果。设任意一个类D，它的类内粒子数为D_num，即该类的规模为D_num，如果其满足式（5）则认为类内粒子运动矢量是较为准确的，并将这些粒子的flag值置为1；否则，对于粒子的flag值置为0；

D_num>threshold_D        （5）

threshold_D是粒子数目的阈值。

4）对于任意一个flag值为1的粒子P_m，在其周围的区域中选取一个与之最近且flag值为1的粒子P_n,对这两个粒子连线上的所有flag值为0的粒子采用线性插值的方法，得到它们的运动矢量结果。

在本实施例中，将threshold_dis取为42、threshold_angle取为并且threshold_D取为所有类的样本数目平均值的即可。经过以上四个步骤，得到密集的流体运动矢量的初始化结果。

（2）运动矢量和LBM相结合的高度计算

LBM是一种简化的流体运动的微观模型，它将大量流体粒子的微观热运动的统计结果表示为粒子分布函数的形式，通过分布函数的演化，进一步反映流体的宏观运动。

定义标准方向。所谓标准方向就是在图1中所示的F_i(i=1,...,8)的方向。中心粒子受到来自周围8个粒子的标准方向力的作用，这些力对中心粒子在三维空间的高度演化起着作用。

如周围粒子的运动矢量方向与某个标准方向相同，则它在该方向上对中心粒子的作用力较大。根据这一特点，将周围粒子的运动矢量投影到标准方向上，从而能够确定周围粒子对中心粒子的作用。当周围粒子的运动方向的分量与对应标准方向的夹角大于等于90度时，则认为在该标准方向上对中心粒子产生的作用力为零。

流体粒子在三维空间中的高度为f_i(i＝1,...,8)之和，和越大则粒子在三维空间中高度值越大，反之越小。

具体计算步骤如下：

1）假设任意一个尚未计算高度的粒子S(x,y)，其周围粒子为T_i(x,y)(i＝1,...,8)，它们对应的运动矢量为(u_i,v_i)。粒子T_i(x,y)对粒子S(x,y)作用力的标准方向为F_i。设标准方向F_i与向量Hor(u_i,0)的夹角为θ₁，F_i与向量Ver(0,v_i)的夹角为θ₂。用式(6)计算粒子S(x,y)的分布函数f_i（x,y,t)。

2）利用式(7)计算粒子S(x,y)在三维空间中的高度h。

$h(r,t)=\underset{i=0}{\overset{d}{\mathrm{\&Sigma;}}}{f}_i(r,t)---\left(7\right)$

3）判断是否所有粒子的高度已经计算，如果是，转步骤4，否则，转步骤1。

4）算法结束。

（3）对高度计算结果的去噪和平滑

先利用高度变化的连续性，对高度计算的结果进行去噪，然后在局部区域内进行平滑处理。其具体步骤如下：

1）对任意一个尚未处理的粒子W(x，y)，利用式(8)和(9)计算其高度h(x,y)的二阶差分，并记为及

$\frac{{\&PartialD;}^{2}h}{\&PartialD;x}=h(x+1,y)-2*h(x,y)+h(x-1,y)---\left(8\right)$

$\frac{{\&PartialD;}^{2}h}{\&PartialD;y}=h(x,y+1)-2*h(x,y)+h(x,y-1)---\left(9\right)$

2)采用式(10)和式(11)对和进行判断。如果同时满足(10)和式(11)，置粒子W(x，y)的flagh值为1；否则，将粒子W(x，y)的高度h(x，y)和flagh值置0。

$\left|\frac{{\&PartialD;}^{2}h}{\&PartialD;x}\right|<\mathrm{threshold}---\left(10\right)$

$\left|\frac{{\&PartialD;}^{2}h}{\&PartialD;y}\right|<\mathrm{threshold}---\left(11\right)$

3)判断是否所有的粒子都进行了处理，如果是，转步骤4；否则，转步骤1。

4)对任意一个尚未平滑处理的粒子W(x，y)，取以它为中心s×s区域，用公式(12)对高度进行平滑。

$h(x,y)=\underset{i=1}{\overset{s\&times;s}{\mathrm{\&Sigma;}}}(\mathrm{flagh}\left(i\right)\&times;h\left(i\right))/\underset{i=1}{\overset{s\&times;s}{\mathrm{\&Sigma;}}}\mathrm{flagh}\left(i\right)---\left(12\right)$

其中，flagh(i)和h(i)分别为s×s邻域中第i个粒子的flagh值和高度。

5)判断是否所有的粒子都进行了平滑处理，如果是，转步骤6；否则，转步骤4。

6)算法结束。

将式(10)和式(11)中的threshold取为0.5，对于式(12)中区域s取为3即可，经去噪、平滑处理后得到了连续的流体高度结果。

(4)利用LBM递推流体的高度

为了满足流体实时重建的实际需要，用递推高度的方法产生中间连续帧的高度值。采用LBM方法，利用粒子分布函数f_i来实现碰撞和平流的过程。

利用LBM方法的高度计算包括平流、碰撞和边界处理的过程。具体为：

1)在平流的过程中，在运动矢量初始化以后，可以得到当前帧的高度，然后利用公式(13)计算下一帧的各个方向上的粒子分布函数。

f_i(r+e_iδ_t，t+δ_t)＝f_i(r，t)    (13)

其中，f_i(r,t)是t时刻位于r处的粒子分布函数。

2）在碰撞的过程中，利用公式（14）计算粒子分布函数，并利用公式(15)计算公式平衡态分布函数。

$f_i(r+e_i{\mathrm{\&delta;}}_t,t+{\mathrm{\&delta;}}_t)-f_i(r,t)=-\frac{1}{\mathrm{\&tau;}}[f_i(r,t)-f_i^{\mathrm{eq}}(r,t)]---\left(14\right)$

$f_i^{\mathrm{eq}}={\mathrm{\&omega;}}_i\mathrm{\&rho;}[1+\frac{e_i\&CenterDot;c}{e_s^2}+\frac{{(e_i\&CenterDot;c)}^2}{2e_s^4}-\frac{c^2}{2e_s^2}]---\left(15\right)$

其中，为沿i方向上的平衡态分布函数；ω_i表示沿i方向的权系数；ρ表示流体的宏观密度；e_s为常数；e_i和c的意义与前述相同。递推时松弛时间τ取为0.8。

3）利用公式（7）即可计算出流体高度。

通过周期边界方法对边界进行处理。通常将式(15)中的参数取为：

$e_i=\left[\begin{array}{ccccccccc}0& 1& 0& -1& 0& 1& -1& -1& 1\\ 0& 0& 1& 0& -1& 1& 1& -1& -1\end{array}\right]$

${\mathrm{\&omega;}}_i=\left\{\begin{array}{c}4/9,e_i^2=0\\ 1/9,e_i^2=1\\ 1/36,e_i^2=2\end{array}\right.$ $e_s=\frac{1}{\sqrt{3}}$

（5）高度校正

为了减少递推过程中的积累误差，在第K帧处对流体高度进行校正。第K帧的高度值可以由其运动矢量计算得到，较为准确。对于(T,K)之间的各帧的高度，采用分布函数的插值方法计算各周围粒子的分布函数，进一步用分布函数求和的方法来计算中心粒子的高度信息。在该实施例中，取T为20、K为40。流体高度校正的具体步骤如下：

1）对第T至第K帧中任意尚未校正高度的帧X，设在(x,y)处粒子i方向的粒子分布函数分别为f_Ti(x,y,t)和f_Ki(x,y,t)，如果它们之间相差过大，即不满足公式（16）的条件，则可能产生流体高度突变的情况，为此按照公式（17）调整f_Ki的大小。转步骤2。

$\left|\frac{f_{\mathrm{Ki}}-f_{\mathrm{Ti}}}{f_{\mathrm{Ti}}}\right|\&le;\mathrm{factor}---\left(16\right)$

f_Ki=f_Ti×(1+sig_factor×factor)     （17）

公式（16）和公式（17）中，factor用来衡量粒子分布函数的增长速度。sig_factor为f_Ki-f_Ti的符号，如果为正取1，如果为负取-1，否则取0。

2)第X(X∈(T，K))帧的粒子分布函数利用公式(18)线性插值得到。

f_Xi＝k×(f_Ki-f_Ti)/(K-T-1)     (18)

其中，k(k＝1，2……，K-T-1)，f_Xi表示第X帧的i方向的粒子分布函数。

3)判断是否存在尚未进行高度校正的帧，如果存在，转步骤1；否则转步骤4。

4)算法结束。

利用式(19)计算factor。

$\mathrm{factor}=\frac{K-T-1}{f_{\mathrm{Ti}}+f_{\mathrm{Ki}}}---\left(19\right)$

经过高度校正后，得到了与视频流体真实感相一致的高度计算结果。

运动矢量和LBM相结合的高度计算方法的步骤：

为了实时计算具有保持流体物理运动特性和流体连续性的流体高度，获得与视频真实感一致的流体高度结果，将视频每m+1帧作为一个周期，m取40。假设对于视频序列中的m+1帧Seg_i(i＝k0，k1，...km)，首先计算Seg_k0和Seg_km帧的运动矢量，利用运动矢量的结果并结合LBM进行高度的计算，这m+1帧高度计算的步骤如下：

1)流体运动矢量的初始化。对Seg_k0帧和Seg_km帧分别计算密集的运动矢量；

2)运动矢量和LBM相结合的高度计算。利用运动矢量和LBM结合的高度计算方法，先计算出帧Seg_k0和帧Seg_km的粒子分布函数和进一步计算这两帧的流体高度h^(k0)(x，y，t)和h^(km)(x，y，t)；

3)对高度计算结果的去噪和平滑。利用流体运动连续性，对计算得到的流体高度h^(k0)(x，y，t)和h^(km)(x，y，t)进行去噪和平滑；

4)利用LBM递推流体的高度。以作为初始分布函数，计算出连续num(num∈[k0，km])帧的高度结果。

5)高度校正。对高度结果进行校准，从而得到了较为准确的高度计算结果。

高度计算的可视化结果可以证明利用本发明对剧烈运动的海洋区域计算时，高度计算结果能明显区分出较高的浪和平静的海面的不同运动特征；对运动平缓的水域进行高度计算时，流体整体运动平缓，计算出的流体高度变化不大，计算结果仍保留了波浪的细节，能够区分波峰和波谷。运动矢量的高度计算结果具有保持流体运动的特征，计算出的高度与实际流体运动时的高度相一致。附图2为“6482810”的第61帧运动矢量的初始化结果。

本发明与同类算法的比较。

为了说明算法的有效性，将本发明与现有方法进行比较。利用三种不同类型的场景进行实验对比：一种是含有倒影的流体场景；二是背景区域不动的流体场景；三是具有局部细节的运动场景。

1）对含有倒影的流体场景，采用现有方法计算时，由于依赖场景的亮度，不能有效的进行高度计算。由于被倒影覆盖的部分区域亮度较暗，所以计算出的高度较低，而对于较亮区域，计算出的高度值较大，但是该区域在实际场景中高度比较平缓，现有方法在处理这种场景时，精度还有待于改进。利用本发明进行计算时，计算结果比较符合实际。

2）背景区域不动的流体场景。实验结果证明，对背景区域不动的流体场景，本发明能区分出背景区域，并对流体区域进行高度计算，计算结果较为准确，而利用现有的方法进行计算时，高度计算结果误差较大。

3）具有局部细节的运动场景。实验结果证明，在对于具有局部细节的运动场景进行对比时，采用本发明中能计算出比较符合真实的高度结果；对于现有方法进行实验，存在的主要问题是，当绘制高度镜面反射的水面时，例如含有阴影或强亮度的水面，对局部细节的高度计算中，精度受到了影响，得到的高度结果与视频真实情况差别较大。而利用本发明计算时，高度计算结果比较符合实际。

与同类算法的对比实验说明了本发明在计算流体高度时，特别是在计算含有倒影的流体场景、背景区域不动的流体场景以及具有局部细节的运动场景时，本发明较为符合实际，进一步可以体现出本发明的有效性。

为了验证本发明的准确性，采用相同的视图来显示原始视频的帧以及三维重建的结果，为了验证其准确性，实验中利用这两个视图结果的平均颜色进行比较，误差error的计算方法见公式（20）。

$\mathrm{error}=\frac{1}{n}\&times;\left(\right|\underset{w=1}{\overset{n}{\mathrm{\&Sigma;}}}{r}_{2d\_w}-\underset{w=1}{\overset{n}{\mathrm{\&Sigma;}}}{r}_{3d|\_w}|+|\underset{w=1}{\overset{m}{\mathrm{\&Sigma;}}}{g}_{2d\_w}-\underset{w=1}{\overset{n}{\mathrm{\&Sigma;}}}{g}_{3d\_w}|+|\underset{w=1}{\overset{n}{\mathrm{\&Sigma;}}}{b}_{2d\_w}-\underset{w=1}{\overset{n}{\mathrm{\&Sigma;}}}{b}_{3d\_w}\left|\right)---\left(20\right)$

其中，n表示图像中像素总数。r_{2d_w}、g_{2d_w}和b_{2d_w}分别是原始视频的帧视图中w像素颜色的三个分量。r_{3d_w}、g_{3d_w}和b_{3d_w}是三维视图中w像素颜色的三个分量。分别利用本文的新方法和现有方法，对动态纹理库DynTex中一些视频进行高度计算，然后利用公式（20）分别计算误差，误差对比结果见表1。

表1

从误差的计算结果中可以明显地看出，在实验中，利用本发明得到的误差较小，高度计算的结果比较符合实际，这进一步说明本发明计算结果的准确性和有效性。

本发明的时间性能分析。为了说明本发明高度计算的时间性能，采用连续100帧的平均运行时间来测试本发明的时间性能，计算结果见表2。从表2的结果可以看出，本发明具有较低的运行时间。对所有视频的运行时间进行平均统计，得到每一帧计算高度需要的时间约0.096787秒，即帧率可以达到为10.33帧/秒，从这些运行时间的结果可以看出，本发明运行需要较少的时间，能够满足流体三维重建的实际需要。

表2 新算法每100帧的平均时间(单位：秒)

Claims (1)

1.一种视频流体高度的计算方法，其特征在于该方法利用流体运动矢量的计算结果，根据粒子运动时的相互作用，计算得到流体粒子的分布函数，再根据分布函数，计算得到流体运动的高度，具体包括以下步骤：

(1)流体运动矢量的计算；

(2)运动矢量和LBM(Lattice Boltzmann Method)相结合的高度计算；

(3)对高度计算结果的去噪和平滑；

(4)利用LBM递推流体的高度；

(5)高度校正；其中：

所述流体运动矢量的计算：利用连续两帧图像，采用区域相关运算方式对流体运动矢量进行初始化；再运用所定义的运动矢量的特征向量，对初始化的结果进行聚类，并统计各类别的规模，保留主要类别，最后用主要类别中的流体运动矢量的结果进行线性插值，得到运动矢量计算结果；其中，用于聚类的特征向量定义为：

vector＝[x,y,u,v,sigu,sigv]     (1)

(x,y)为流体粒子所在的位置，u和v分别表示该粒子的运动矢量沿x方向和沿y方向的分量，sigu和sigv分别表示u和v的符号，如果为正取1，如果为负取-1，否则取0；具体计算包括：

1)利用连续两帧图像，采用区域相关运算对流体运动矢量进行初始化；

2)利用运动矢量的特征向量vector＝[x,y,u,v,sigu,sigv]进行聚类，具体地，如果粒子P_i和P_j的运动矢量的特征向量分别为：vector_i＝[x_pi,y_pi,u_pi,v_pi,sigu_pi,sigv_pi]及

vector_j＝[x_pj,y_pj,u_pj,v_pj,sigu_pj,sigv_pj]，并且，它们同时满足式(2)至(5)，则认为粒子P_i和P_j的运动具有相似性，并将它们聚到相同的类别之中；

$\left|\sqrt{{(x_{\mathrm{pi}}-x_{\mathrm{pj}})}^2+{(y_{\mathrm{pi}}-y_{\mathrm{pj}})}^2}\right|<{\mathrm{threshold}}_{\mathrm{dis}}---\left(2\right)$

$\left|\sqrt{\mathrm{angle}({\mathrm{vec}}_{\mathrm{pi}},{\mathrm{vec}}_{\mathrm{pj}})}\right|<{\mathrm{threshold}}_{\mathrm{angle}}---\left(3\right)$

sigu_pi*sigu_pj≥0      (4)

sigv_pi*sigv_pj≥0     (5)

其中，threshold_dis表示粒子P_i和P_j之间距离的阈值；vec_pi和vec_pj分别表示粒子P_i和P_j的运动矢量(u_pi,v_pi)及(u_pj,v_pj)；threshold_angle表示它们运动矢量方向的阈值；

3)统计各类别的规模，对于聚类结果中规模较大的类别，保留其结果；设任意一个类D，它的类内粒子数为D_num，即该类的规模为D_num，如果其满足式(6)则认为类内粒子运动矢量是较为准确的，并将这些粒子的flag值置为1；否则，对于粒子的flag值置为0；

D_num>threshold_D            (6)

其中，threshold_D是粒子数目的阈值；

4)对于任意一个flag值为1的粒子P_m，在其周围的区域中选取一个与之最近且flag值为1的粒子P_n,对这两个粒子连线上的所有flag值为0的粒子采用线性插值的方法，得到它们的运动矢量结果；

所述运动矢量和LBM相结合的高度计算包括：

a)假设任意一个尚未计算高度的粒子S(x,y)，其周围粒子为T_i(x,y)(i＝1,...,8)，它们对应的运动矢量为(u_i,v_i)，粒子T_i(x,y)对粒子S(x,y)作用力的标准方向为F_i，设标准方向F_i与水平向量Hor(u_i,0)的夹角为θ₁，F_i与垂直向量Ver(0,v_i)的夹角为θ₂，用下式计算粒子S(x,y)的分布函数f_i(x,y,t)：

b)利用下式计算粒子S(x,y)在三维空间中的高度h(r,t)

$h(r,t)=\underset{i=0}{\overset{d}{\mathrm{\&Sigma;}}}{f}_i(r,t)$

其中，h(r,t)为t时刻位于r处的粒子高度；d表示周围粒子个数；f_i(r,t)表示t时刻位于r处的粒子分布函数；

c)判断是否所有粒子的高度已经计算，如果是，转步骤d)，否则，转步骤a)；

d)计算结束。