视频流体物理驱动模型恢复及重新仿真的方法

摘要

本发明公开了一种视频流体物理驱动模型恢复及重新仿真的方法，该方法旨在利用视频流的样例帧，得到反向工程的流体仿真结果。该方法首先采用HS(Horn-Schunck)方法，利用多尺度框架，在每层对第一帧流体粒子的二维运动速度进行计算，然后在每层中利用浅水波方程做约束，得到了较为准确的流体运动速度及非归一化几何。在流体仿真过程中，利用样例帧表面的几何作为约束，调整流体粒子运动的密度属性，从而得到了流体仿真的物理驱动模型。最后利用增强的粒子密度进行重仿真，得到了具有真实感细节的重仿真结果。本发明可有效地应用于流体的反向工程的研究及应用中，可以定制任意场景，克服了基于物理的传统流体仿真中计算量大以及实时性差的问题。

说明书

技术领域

本发明涉及视频流体反向工程中物理属性及驱动模型的恢复。利用运动矢量计算的多尺度框架及浅水波方程优化得到流体表面的非归一化几何和初始帧的二维运动速度，恢复流体运动属性，得到反向工程中的驱动模型，获得具有真实感细节的重仿真结果。

背景技术

真实感流体仿真目前已经广泛地被应用于计算机游戏、动画、电影特效、军事仿真等领域，关于流体仿真计算的研究已经取得了很多的成果。限于流体物理仿真中的大量物理属性计算的问题，近些年人们开始从事于基于视频的流体仿真技术的研究。怎样充分利用流体视频中的丰富真实感信息及流体的运动信息，进一步生成实时的、具有真实感特性的流体仿真场景，这就是视频流体反向工程中亟待解决的问题。

近年来，一些研究工作已经开始从事基于视频的建模及绘制的研究，一些研究中，采用能量优化的策略，并权衡计算量及外观的特性，从而使得生成的流体仿真场景不具有流体真实的物理特性；在另外的相关研究中，人们仅仅恢复流体表面的归一化几何，由于这些几何信息仅能反映流体表面的外观特点，还不能用于对于流体表面的真实物理属性的恢复。

由于流体运动的特殊性，对于流体反向工程的研究存在着一定难度。目前人们在研究中利用流体运动速度跟踪的方法来实现流体重仿真的过程，存在的主要问题是流体运动的物理驱动模型并没有真正得到恢复；在另外的反向工程相关研究中，人们试图采用在低维空间中进行流体仿真，为了保持流体仿真的细节，如何从低维空间建立简化的仿真模型，这已经成为研究中瓶颈，尽管人们已经做了一些努力探索这个问题，目前如何在流体方向工程的低维空间中，快速生成与原流体场景相似结果以及长序列的流体重仿真结果，仍然是一个亟待解决的实际问题。

发明内容

本发明的目的是针对现有技术的不足而提供的一种视频流体反向工程中物理属性及驱动模型恢复的方法，利用该方法能够得到视频流体反向工程中流体物理属性及驱动模型，利用该模型可以得到流体重仿真的结果。

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

(1)利用多尺度的框架，计算视频流体第一帧的表面几何和流体二维运动速度

为了计算视频流体的表面几何和流体粒子二维运动速度，采用以下两个条件约束：流体表面的外观的约束和流体运动的动态行为约束。

首先，利用多尺度光流计算的框架，采用5层尺度规模，初始化流体表面粒子的二维运动速度。利用视频流体初始两帧的灰度图像I₀和I₁产生5层不同尺度的图像其中s＝{0,1,2,3,4},及对应着较粗糙的模型，而及对应着较精细的模型。将每层流体粒子运动速度记为t表示帧的序号，s表示尺度规模，s＝{0,1,2,3,4}。从第0层开始，每一层得到的二维运动速度传给下一层。每层的计算按照以下三个步骤进行：

(a)二维运动速度的初始化

利用现有多尺度光流计算方法,计算流体第一帧的二维速度

(b)非归一化几何的计算

在每层中，利用现有方法计算流体表面的归一化几何,表示为z；非归一化几何表示为Z，比例系数(我们也称之为非归一化因子)假设为S，那么，

Z＝Sz        (1)

利用浅水波方程，根据现有的研究结果中速度与高度之间的计算关系

$>V^{k+1}=\frac{1}{a_{u1}{a}_{v2}-a_{u2}{a}_{v1}}(-a_{v1}{C}_u^k+a_{u1}{C}_v^k)---\left(24\right)$>

$>u^{k+1}=\frac{1}{a_{u1}{a}_{v2}-a_{u2}{a}_{v1}}(a_{v2}{C}_u^k-a_{u2}{C}_v^k)---\left(25\right)$>

将二维运动速度的初始化的结果看作是较准确的，因此可以得到：

$>\frac{u}{v}=\frac{(a_{v2}{C}_u-a_{u2}{C}_v)}{(-a_{v1}{C}_u+a_{u1}{C}_v)}---\left(2\right)$>

其中，

$>C_u[\frac{\partial }{\partial x}{Z}^2{u}_x+Z{Z}_y{v}_x+(-{\mathrm{ZZ}}_x+\frac{\partial }{\partial x}{Z}^2)v_y+(Z^2+a^2)\tilde{u}+{\alpha }^2\bar{u}+Z^2{v}_{\mathrm{xy}}-Z_x{Z}_t+\frac{\partial }{\partial x}\left({\mathrm{ZZ}}_t\right)]---\left(3\right)$>

$>C_v=[(-{\mathrm{ZZ}}_y+\frac{\partial }{\partial y}{Z}^2)u_x+{\mathrm{ZZ}}_x{u}_y+\frac{\partial }{\partial y}{Z}^2{v}_y+Z^2{u}_{\mathrm{xy}}+{\alpha }^2\tilde{v}(Z^2+{\alpha }^2)\bar{v}-Z_y{Z}_t+\frac{\partial }{\partial y}\left({\mathrm{ZZ}}_t\right)]---\left(4\right)$>

$>a_{u1}=Z_x^2-\frac{\partial }{\partial x}\left({\mathrm{ZZ}}_x\right)+Z^2+{2\alpha }^2---\left(5\right)$>

$>a_{u2}=Z_x{Z}_y-\frac{\partial }{\partial x}\left({\mathrm{ZZ}}_y\right)---\left(6\right)$>

$>a_{v1}=Z_x{Z}_y-\frac{\partial }{\partial y}\left({\mathrm{ZZ}}_x\right)---\left(7\right)$>

$>a_{v2}=Z_y^2-\frac{\partial }{\partial y}\left({\mathrm{ZZ}}_y\right)+Z^2+2{\alpha }^2---\left(8\right)$>

$>\tilde{u}=\frac{u(x+h,y)+u(x-h,y)}{2h}---\left(9\right)$>

$>\bar{u}=\frac{u(x,y+h)+u(x,y-h)}{2h}---\left(10\right)$>

$>\tilde{v}=\frac{v(x+h,y)+v(x-h,y)}{2h}---\left(11\right)$>

$>\bar{v}=\frac{v(x,y+h)+v(x,y-h)}{2h}---\left(12\right)$>

其中，u_x、u_y表示运动速度第一个分量u在x轴和y轴方向的导数，v_x、v_y表示运动速度第二个分量v在x轴和y轴方向的导数，u_xy表示运动速度第一个分量u在x轴和y轴方向的混合偏导数，v_xy表示运动速度第二个分量v在x轴和y轴方向的混合偏导数。Z_x表示Z在x轴方向的导数，Z_y表示Z在y轴方向的导数,Z_t表示Z在时序上的导数。h表示相邻像素之间的坐标差，取1值，α是平滑系数，取7.0。

将式(1)带入到方程式(2)中，可以得到有关S的四次方程：

P₁S⁴+P₂S²+P₃＝0            (13)

利用t＝S²对上式进行换元，得到有关t的一元二次方程：

P₁t²+P₂t+P₃＝0          (14)

其中，

P₁＝CDu+DFv-ABu-BEv       (15)

$>P_2=-{\alpha }^2(\tilde{u}+\bar{u})(\mathrm{Au}+\mathrm{Ev})+{\alpha }^2(\tilde{v}+\bar{v})(\mathrm{Cu}+\mathrm{Fv})+2{\alpha }^2\mathrm{Du}-{\alpha }^2\mathrm{Bv}---\left(16\right)$>

$>P_3=2{\alpha }^{4}u(\tilde{v}+\bar{v})-2{\alpha }^{4}v(\tilde{u}+\bar{u})---\left(17\right)$>

$>A=z_x{z}_y-\frac{\partial \left({\mathrm{zz}}_x\right)}{\partial y}---\left(18\right)$>

$>B=\frac{\partial }{\partial x}{z}^2{u}_x+{\mathrm{zz}}_y{v}_x-{\mathrm{zz}}_x{v}_y+\frac{\partial }{\partial x}{z}^2{v}_y+z^2\tilde{u}+z^2{v}_{\mathrm{xy}}-z_x{z}_t+\frac{\partial }{\partial x}\left({\mathrm{zz}}_t\right)---\left(19\right)$>

$>C=z_x^2-\frac{\partial }{\partial y}\left({\mathrm{zz}}_x\right)+z^2---\left(20\right)$>

$>D=-z{z}_y{u}_x+\frac{\partial }{\partial y}\left(z^2\right)u_x+{\mathrm{zz}}_x{u}_y+\frac{\partial z^2}{\partial y}{v}_y+z^2{u}_{\mathrm{xy}}+z^2\bar{v}-z_y{z}_t+\frac{\partial }{\partial y}\left({\mathrm{zz}}_t\right)---\left(21\right)$>

$>E=z_y^2-\frac{\partial }{\partial y}\left({\mathrm{zz}}_y\right)+z^2---\left(22\right)$>

$>F=z_x{z}_y-\frac{\partial \left({\mathrm{zz}}_y\right)}{\partial x}---\left(23\right)$>

根据式(14)求得t，进一步能够得到S；在求取t和S时，一元二次方程有两个解，舍去不合理的一个解；

(c)流体表面几何和速度的优化

按照以下步骤循环两次进行优化：步骤1，将步骤(b)得到的第k次非归一化几何看作是较准确的(第一次迭代假设k＝0)，利用以下公式(24)和(25)可以得到第k+1次更准确的速度(u^k+1,v^k+1)；步骤2，将步骤1得到的速度看作是较准确的结果，再利用方程(13)计算得到非归一化因子S，从而得到更准确的表面几何；

$>v^{k+1}=\frac{1}{a_{u1}{a}_{v2}-a_{u2}{a}_{v1}}(-a_{v1}{C}_u^k+a_{u1}{C}_v^k)---\left(24\right)$>

$>u^{k+1}=\frac{1}{a_{u1}{a}_{v2}*a_{u2}{a}_{v1}}(a_{v2}{C}_u^k-a_{u2}{C}_v^k)---\left(25\right)$>

每个粒子利用上述第L-1层的步骤(c)计算得到流体表面几何的非归一化因子，以下步骤是利用粒子的非归一化因子来计算整个场景的S：如果任意两个粒子i和j计算得到的非归一化因子分别为S_i和S_j，聚类时，如果S_i和S_j满足以下条件，那么S_i和S_j就标记为同一类：

|S_i-S_j|<S_T i,j∈{0,1,2,3…N-1}       (26)

其中N为粒子总数目，S_T为S的阈值，设其值为1.0.

聚类时类的数量是自动得到的，如果聚得最大类的样本数为n，最大类中所有粒子求得的非归一化因子的和为S_s，那么，

S＝S_s/n      (27)

(2)利用步骤(1)得到流体表面几何的非归一化因子S，再计算视频样例帧非归一化几何取视频中前100帧作为样例帧，按照步骤(1)得到流体表面几何的非归一化因子S，再利用式(1)计算这些样例帧表面非归一化几何；

(3)采用LBM(Lattice Boltzmann Method)，恢复流体运动的属性及驱动模型

对于所有样例帧，利用LBM进行仿真时，利用式(28)和(29)计算并调整粒子的粒子密度和分布函数；

$>f_i(t+1)=\frac{Z}{h_L}{f}_i\left(t\right),i=\{\mathrm{0,1}...8\}---\left(28\right)$>

$>\rho (t+1)=\frac{Z}{h_L}\rho \left(t\right)---\left(29\right)$>

其中，h_L是粒子演化的高度，f_i(t+1)和f_i(t)分别是粒子在i方向的t+1时刻和t时刻的分布函数，i表示LBM的D2Q9模型中第i个方向，i＝0,1,2…8,ρ(t+1)和ρ(t)分别是粒子在t+1时刻和t时刻的密度；

(4)对流体粒子仿真密度属性进行增强

将步骤(3)中得到的粒子密度我们称其为标准密度，记为ρ_s；LBM自动推演产生的密度，称之为自动演化密度，记为ρ_a；分析ρ_s的直方图的分布，并对ρ_a进行增强，增强后的密度称之为增强密度，记为ρ_e；

对于ρ_s，求其最大值和最小值将区间等分成256个子区间：

[is_s，(i+1)s_s] i＝{0，1...255}      (30)

其中，s_s为区间尺度：

$>s_s=\frac{({\rho }_m^s-{\rho }_n^s)}{256}---\left(31\right)$>

标准密度ρ_s表示为ρ_s＝S_s×s_s+O_s，将ρ_s记为[S_s,O_s],其中S_s为区间号，O_s为区间内密度位移量，其满足：

对ρ_a，求其最大值和最小值将区间等分成256个子区间：

[is_a，(i+1)s_a] i＝{0，1...255}      (33)

其中，s_a为区间尺度：

$>s_a=\frac{({\rho }_m^a-{\rho }_n^a)}{256}---\left(34\right)$>

对于任何一个自动演化密度ρ_a，都可以表示为ρ_a＝S_a×s_a+O_a，将ρ_a记为[S_a,O_a]其中S_a为区间号，O_a为区间内密度位移量，其满足：

根据ρ_s的直方图的统计分布，利用图像增强理论中直方图规定化的方法，实现S_s和S_a的对应，如果对于任意一个自动演化密度它的区间号为如果和它匹配的区间号为它被增强后的密度ρ_e可以计算为：

$>{\rho }_e={\rho }_n^s+\frac{{\rho }_m^s-{\rho }_n^s}{256}{S}_s^m+\frac{{\rho }_m^s-{\rho }_n^s}{{\rho }_m^a-{\rho }_n^a}\times O_a^m---\left(36\right)$>

(5)重仿真过程

视频流体的重仿真过程包括以下步骤：

(a)利用步骤(1)初始化流体初始帧的运动速度及非归一化几何

(b)利用步骤(2)计算样例帧表面非归一化几何

(c)利用初始帧运动速度，推演LBM，同时利用步骤(3)和(4)增强流体粒子的密度

利用增强的密度演化仿真，得到具有细节特征的真实感流体重仿真结果。

本发明具有简单、实用的特点，能够实现流体物理属性及驱动模型的恢复，有效地实现具有细节特征的真实感流体的重仿真过程。克服了现有流体仿真的反向工程中，缺少物理模型驱动的问题。进一步实例证明了该发明方法的有效性。可以实现流体场景的放大、缩小、定制功能，并且所恢复的物理属性及驱动模型能够满足边界条件。

附图说明

图1为本发明实施例流体粒子前两帧之间运动速度初始化结果图；

图2为本发明实施例流体粒子表面非归一化几何的计算结果图；

图3为本发明实施例流体粒子的密度计算可视化的结果图；

图4为本发明实施例流体重仿真的结果图；

图5为本发明实施例中恢复流体物理属性的边界条件图；

图6为本发明实施例中烟的重仿真结果，及与现有方法的比较图；

图7为本发明实施例反向工程中场景定制尺度的应用实例图；

图8为本发明实施例反向工程中定制场景的应用实例图。

具体实施方式

实施例

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

本实施例采用DynTex动态纹理库中的流体视频进行物理属性及驱动模型的恢复。在PC机上Windows7操作系统下进行实施，其硬件配置是2.66GHz Intel Core(TM)2 Duo CPU、4GB RAM。(1)利用多尺度的框架，计算视频流体前两帧的表面几何和流体粒子二维运动速度

(a)二维运动速度的初始化

利用现有多尺度光流计算方法,采用5层尺度规模,计算流体第一帧的二维速度(u,v).

(b)非归一化几何的计算

在每层中，利用现有方法计算流体表面的归一化几何,对归一化几何我们表示为z；非归一化几何我们表示为Z,比例系数假设为S，那么，

Z＝Sz(1)

利用浅水波方程，根据现有的研究结果中速度与高度之间的计算关系，并将二维运动速度的初始化的结果看作是较准确的，因此可以得到速度的计算关系：

$>\frac{u}{v}=\frac{(a_{v2}{C}_u-a_{u2}{C}_v)}{(-a_{v1}{C}_u+a_{u1}{C}_v)}---\left(2\right)$>

其中，

$>C_u[\frac{\partial }{\partial x}{Z}^2{u}_x+Z{Z}_y{v}_x+(-{\mathrm{ZZ}}_x+\frac{\partial }{\partial x}{Z}^2)v_y+(Z^2+a^2)\tilde{u}+{\alpha }^2\bar{u}+Z^2{v}_{\mathrm{xy}}-Z_x{Z}_t+\frac{\partial }{\partial x}\left({\mathrm{ZZ}}_t\right)]---\left(3\right)$>

$>C_v=[(-{\mathrm{ZZ}}_y+\frac{\partial }{\partial y}{Z}^2)u_x+{\mathrm{ZZ}}_x{u}_y+\frac{\partial }{\partial y}{Z}^2{v}_y+Z^2{u}_{\mathrm{xy}}+{\alpha }^2\tilde{v}(Z^2+{\alpha }^2)\bar{v}-Z_y{Z}_t+\frac{\partial }{\partial y}\left({\mathrm{ZZ}}_t\right)]---\left(4\right)$>

$>a_{u1}=Z_x^2-\frac{\partial }{\partial x}\left({\mathrm{ZZ}}_x\right)+Z^2+{2\alpha }^2---\left(5\right)$>

$>a_{u2}=Z_x{Z}_y-\frac{\partial }{\partial x}\left({\mathrm{ZZ}}_y\right)---\left(6\right)$>

$>a_{v1}=Z_x{Z}_y-\frac{\partial }{\partial y}\left({\mathrm{ZZ}}_x\right)---\left(7\right)$>

$>a_{v2}=Z_y^2-\frac{\partial }{\partial y}\left({\mathrm{ZZ}}_y\right)+Z^2+2{\alpha }^2---\left(8\right)$>

$>\tilde{u}=\frac{u(x+h,y)+u(x-h,y)}{2h}---\left(9\right)$>

$>\bar{u}=\frac{u(x,y+h)+u(x,y-h)}{2h}---\left(10\right)$>

$>\tilde{v}=\frac{v(x+h,y)+v(x-h,y)}{2h}---\left(11\right)$>

$>\bar{v}=\frac{v(x,y+h)+v(x,y-h)}{2h}---\left(12\right)$>

其中，u_x、u_y表示运动速度第一个分量u在x轴和y轴方向的导数，v_x、v_y表示运动速度第二个分量v在x轴和y轴方向的导数，u_xy表示运动速度第一个分量u在x轴和y轴方向的混合偏导数，v_xy表示运动速度第二个分量v在x轴和y轴方向的混合偏导数；Z_x表示Z在x轴方向的导数，Z_y表示Z在y轴方向的导数,Z_t在时序上的导数；h表示相邻像素之间的坐标差，取1值，α是平滑系数，取7.0；

将式(1)带入到浅水波方程式(2)中，可以得到有关S的四次方程：

P₁S⁴+P₂S²+P₃＝0         (13)

我们利用t＝S²对上式进行换元，得到有关t的一元二次方程：

P_lt²+P₂t+P₃＝0    (14)

其中，

P₁＝CDu+DFv-ABu-BEv     (15)

$>P_2=-{\alpha }^2(\tilde{u}+\bar{u})(\mathrm{Au}+\mathrm{Ev})+{\alpha }^2(\tilde{v}+\bar{v})(\mathrm{Cu}+\mathrm{Fv})+2{\alpha }^2\mathrm{Du}-{\alpha }^2\mathrm{Bv}---\left(16\right)$>

$>P_3=2{\alpha }^{4}u(\tilde{v}+\bar{v})-2{\alpha }^{4}v(\tilde{u}+\bar{u})---\left(17\right)$>

$>A=z_x{z}_y-\frac{\partial \left({\mathrm{zz}}_x\right)}{\partial y}---\left(18\right)$>

$>B=\frac{\partial }{\partial x}{z}^2{u}_x+{\mathrm{zz}}_y{v}_x-{\mathrm{zz}}_x{v}_y+\frac{\partial }{\partial x}{z}^2{v}_y+z^2\tilde{u}+z^2{v}_{\mathrm{xy}}-z_x{z}_t+\frac{\partial }{\partial x}\left({\mathrm{zz}}_t\right)---\left(19\right)$>

$>C=z_x^2-\frac{\partial }{\partial y}\left({\mathrm{zz}}_x\right)+z^2---\left(20\right)$>

$>D=-z{z}_y{u}_x+\frac{\partial }{\partial y}\left(z^2\right)u_x+{\mathrm{zz}}_x{u}_y+\frac{\partial z^2}{\partial y}{v}_y+z^2{u}_{\mathrm{xy}}+z^2\bar{v}-z_y{z}_t+\frac{\partial }{\partial y}\left({\mathrm{zz}}_t\right)---\left(21\right)$>

$>E=z_y^2-\frac{\partial }{\partial y}\left({\mathrm{zz}}_y\right)+z^2---\left(22\right)$>

$>F=z_x{z}_y-\frac{\partial \left({\mathrm{zz}}_y\right)}{\partial x}---\left(23\right)$>

根据式(14)可以求得t,进一步可以得到S。

(c)流体表面几何和速度的优化

按照以下步骤循环两次进行优化：步骤1，将步骤(b)得到的第k次非归一化几何看作是较准确的(第一次迭代假设k＝0)，利用以下公式(24)和(25)可以得到第k+1次更准确的速度(u^k+1,v^k+1)；步骤2，将步骤1得到的速度看作是较准确的结果，再利用方程(13)计算得到非归一化因子S，从而得到更准确的表面几何。

$>v^{k+1}=\frac{1}{a_{u1}{a}_{v2}-a_{u2}{a}_{v1}}(-a_{v1}{C}_u^k+a_{u1}{C}_v^k)---\left(24\right)$>

$>u^{k+1}=\frac{1}{a_{u1}{a}_{v2}*a_{u2}{a}_{v1}}(a_{v2}{C}_u^k-a_{u2}{C}_v^k)---\left(25\right)$>

每个粒子利用上述第L-1层的步骤(c)计算得到流体表面几何的非归一化因子S，利用以下步骤是利用粒子的非归一化因子来计算整个场景的S；如果任意两个粒子i和j计算得到的S表示为S_i和S_j，聚类时，如果S_i和S_j满足以下条件，那么S_i和S_j就标记为同一类：

|S_i-S_j|<S_T i,j∈{0,1,2,3…N-1}      (26)

其中N为粒子总数目,S_T为S的阈值，我们设其值为1.0.

聚类时类的数量是自动得到的，如果聚得最大类的样本数为n，最大类中所有粒子求得的非归一化因子的和为S_s，那么，

S＝S_s/n     (27)

表1是不同视频计算得到的S的结果。

表1 不同视频计算得到S的结果

(2)利用步骤(1)得到流体表面几何的非归一化因子S，再计算视频样例帧非归一化几何

取视频中前100帧作为样例帧，按照步骤(1)得到流体表面几何的非归一化因子S，再利用式(1)计算这些样例帧表面非归一化几何；

图2是计算得到的流体表面非归一化几何信息。

(3)采用LBM(Lattice Boltzmann Method)，恢复流体运动的属性及驱动模型

对于所有样例帧，利用LBM进行仿真时，利用式(28)和(29)计算并调整粒子的粒子密度和分布函数；

$>f_i(t+1)=\frac{Z}{h_L}{f}_i\left(t\right),i=\{\mathrm{0,1}...8\}---\left(28\right)$>

$>\rho (t+1)=\frac{Z}{h_L}\rho \left(t\right)---\left(29\right)$>

其中，h_L是粒子演化的高度，f_i(t+1)和f_i(t)分别是粒子在i方向的t+1时刻和t时刻的分布函数，i表示LBM的D2Q9模型中第i个方向，i＝0,1,2…8,ρ(t+1)和ρ(t)分别是粒子在t+1时刻和t时刻的密度；

(4)对流体粒子仿真密度属性进行增强

将步骤(3)中得到的粒子密度我们称其为标准密度，记为ρ_s；LBM自动推演产生的密度，称之为自动演化密度，记为ρ_a；分析ρ_s的直方图的分布，并对ρ_a进行增强，增强后的密度称之为增强密度，记为ρ_e；

对于ρ_s，求其最大值和最小值将区间等分成256个子区间：

[is_s，(i+1)s_s] i＝{0，1...255}        (30)

其中，s_s为区间尺度：

$>s_s=\frac{({\rho }_m^s-{\rho }_n^s)}{256}---\left(31\right)$>

标准密度ρ_s表示为ρ_s＝S_s×s_s+O_s，将ρ_s记为[S_s,O_s],其中S_s为区间号，O_s为区间内密度位移量，其满足：

对ρ_a，求其最大值和最小值将区间等分成256个子区间：

[is_a，(i+1)s_a] i＝{0，1...255}       (33)

其中，s_a为区间尺度：

$>s_a=\frac{({\rho }_m^a-{\rho }_n^a)}{256}---\left(34\right)$>

对于任何一个自动演化密度ρ_a，都可以表示为ρ_a＝S_a×s_a+O_a，将ρ_a记为[S_a,O_a]其中S_a为区间号，O_a为区间内密度位移量，其满足：

根据ρ_s的直方图的统计分布，利用图像增强理论中直方图规定化的方法，实现S_s和S_a的对应，如果对于任意一个自动演化密度它的区间号为如果和它匹配的区间号它被增强后的密度ρ_e可以计算为：

$>{\rho }_e={\rho }_n^s+\frac{{\rho }_m^s-{\rho }_n^s}{256}{S}_s^m+\frac{{\rho }_m^s-{\rho }_n^s}{{\rho }_m^a-{\rho }_n^a}\times O_a^m---\left(36\right)$>

图3是视频流体标准密度ρ_s计算的可视化结果，从该结果可以看出该发明专利能够有效性地恢复物理驱动模型。

(5)重仿真过程

视频流体的重仿真过程包括以下四个步骤：

(a)利用步骤(1)初始化流体初始帧的运动速度及非归一化几何

(b)利用步骤(2)计算样例帧表面非归一化几何

(c)利用初始帧运动速度，推演LBM，同时利用步骤(3)和(4)增强流体粒子的密度

(d)利用增强的密度演化仿真，得到具有细节特征的真实感流体重仿真结果。

图4是重仿真的结果，每组中左上结果是原流体视频，右下结果是重仿真的结果。

利用上述的流体重仿真实施方法，进一步验证了边界条件。我们仍利用采用DynTex动态纹理库中的流体视频进行物理属性及驱动模型的恢复以后，再进行物理属性的边界条件测试，采用现有LBM的边界条件的控制方法，改变边界处的分布函数，得到了实施结果，证明边界条件的有效性。图5是边界条件的实施结果。

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

为了验证本发明中，流体运动属性及物理驱动模型恢复的有效性，采用James Gregson2014研究中烟的视频进行重仿真验证，并与现有的研究结果进行比较，图6是对烟视频进行重仿真的结果，左侧是原烟的视频，中间的是现有方法的结果，右面的结果是本发明的结果，从其可以证明流体物理属性及物理驱动模型的恢复也适用于烟流体的反向工程。与现有的基于运动跟踪的方法相比，本发明的仿真结果具有真实感强的特点，并且反向工程的结果与原视频更加相似。

为了验证该发明中物理属性恢复的准确性，与现有的方法进行比较，分别采用现有方法和本发明对流体运动速度进行初始化，然后进行驱动LBM仿真，计算演化高度和计算得到的表面几何之间的差异：

$>E_v=\sqrt{\frac{1}{n_p}\underset{i=0}{\overset{n_p-1}{\Sigma }}{(h_i-Z_i)}^2}---\left(37\right)$>

其中，n是粒子的数目，h_i是粒子i的演化高度，Z_i是粒子i的非归一化的表面几何。

表2列出了与现有方法进行比较的误差结果。

表2 与现有方法进行比较的误差结果

从表中的误差结果可以看出，本发明中恢复的运动速度更加准确一些。

在反向工程中应用的实施。

应用1：利用本发明定制场景规模。可以任意尺度的场景，实施方法：利用图像放大及缩小的理论，生成大场景(或小场景)的密度场，从而可以驱动流体物理模型进行驱动。具体地，如果场景尺度放大的倍数为2×3，采用图像放大理论中的后向映射法生成密度场，如果目标位置的坐标为(x,y),那么可以计算得到其对应的原场景的对应位置(x₀,y₀)：

x₀＝x/2 y₀＝y/3        (38)

那么可以利用对应位置(x₀,y₀)处的密度生成新的坐标处(x,y)的密度。如果计算得到的结果非整型数值，可以采用最近邻插值法或双线性插值方法计算得到(x,y)点处的密度值。图7是利用该方法得到的不同放大倍数的场景结果，每组中，左上是原来尺度重仿真结果，右上是0.7×0.7的缩放比例，即水平方向尺度缩小到0.7倍，垂直方向也缩小到0.7倍。下面的放大比例为1.5×1.5。

应用2：利用本发明的方法，定制密度分布直方图，从而可以得到定制的流体场景。在实施中，我们改变密度直方图的分布，得到了新的流体定制场景。为了可视地表现出新的流通场景高度几何与外观的对应关系，我们生成了颜色与灰度之间对应的映射关系，如果灰度值为g,那么利用原来场景，生成对应关系：C_R[g]、C_G[g]和C_B[g]。在流体仿真中可以得到定制场景的演化高度，再根据现有的基于SFS(Shape from Shading)的重建方法，可以计算得到粒子的灰度值，再利用所建立的颜色与灰度之间的映射关系，可以得到表面的颜色映射结果。图8是定制流体重仿真的结果。每组中左侧是流体重仿真的结果，右侧是定制的流体仿真结果。