CN103489209B - 一种基于流体关键帧编辑的可控流体动画生成方法 - Google Patents

Abstract

本发明公开了一种基于流体关键帧编辑的可控流体动画生成方法。初始流体动画序列中指定关键帧，对关键帧处的流体形状进行预处理，分为三种特征区域；用户输入三维三角形面片或者二维曲线，获取三维控制点集；获得形变速度场参数后计算得到形变速度场；将形变速度场的每个速度模版中心沿着速度场移动得到修改后的速度场，使用修改后的速度场对关键帧形变得到形变后的关键帧；将子序列形变得到形变后的目标子序列；施加外力构造PD控制器后生成最终子序列，将最终子序列与之前序列拼接生成编辑后的流体动画；可重复上述步骤得到最终流体动画。本发明能对流体动画进行任意多次的编辑，使其形状在指定的关键帧处符合用户施加的控制，操作可行且直观。

Description

一种基于流体关键帧编辑的可控流体动画生成方法

技术领域

本发明涉及一种可控流体动画生成方法，尤其是涉及一种基于流体关键帧编辑的可控流体动画生成方法。

技术背景

基于物理的流体动画在计算机图形学，尤其是电影特效领域有着广泛的应用，如流体模拟，烟雾及火焰特效。当前流体动画生成方法已相对成熟，然而要根据用户的控制需求对一段已有的流体动画序列施加控制，目前还存在巨大的挑战，没有成熟的方法。

目前的流体动画控制方法有各自的不足：

1.需要完整的关键帧信息。[可参考LinShiandYizhouYu，TamingLiquidsforRapidlyChangingTargets，ACMSIGGRAPH／EurographicsSymposiumonComputerAnimation，LosAngeles，July2005，pp.229-236].

2.缺乏直观的用户控制方法。[可参考Rasmussen，N.，Enright，D.，Nguyen，D.，Marino，S.，Sumner，N.，Geiger，W.，Hoon，S.，ANDFedkiw，R.2004.Directiblephotorealisticliquids.InProceedingsofthe2004ACMSIGGRAPH／EurographicsSymposiumonComputerAnimation，193-202].

3.形状编辑不符合流体的特征。[可参考Bridson，R.，Houriham，J.，andNordenstam,M.2007.Curl-noiseforproceduralfluidflow.ACMTransactiononGraphics26，3，46].

发明内容

针对背景技术的不足，本发明的目的在于提供一种基于流体关键帧编辑的可控流体动画生成方法。该方法能对一段已有的流体动画序列做任意多次的关键帧编辑，并根据编辑后的关键帧生成一段新的流体动画序列。使得该序列在指定关键帧符合用户的控制要求。

为实现上述的目的，本发明采用的技术方案为：

本发明包含如下步骤：

1)输入初始流体动画序列，指定一个需要编辑的关键帧K，对关键帧K处的流体形状进行预处理，首先对流体形状进行平滑化处理，将关键帧K处的流体形状采用隐式函数表示，然后抽取形状距离场φ和边界距离场φ_Ｓ，将流体形状区域分为三种特征区域R；

2)用户从流体形状的任意视角输入一个三维三角形面片，或者在屏幕上绘制一条二维曲线：对于用户输入的三维三角形面片，直接采样三维三角形面片的网格顶点，得到三维三角形面片的三维控制点集P_M；对于用户绘制的二维曲线，在屏幕上采样一组采样点，将每个采样点X_P与绘制视角的视点位置e连成射线，并建立如下公式1计算二维曲线的深度，即每个采样点X_P的屏幕深度d_Ｓ，然后将e+D_Xd_Ｓ作为二维曲线的三维控制点集P_S：

$d_S=\arg {\min }_d\underset{X_P\∈P}{\mathrm{\Σ}}\underset{p^{*}\∈R}{\min }\mathrm{dist}(e+D_{X_p}d,p^{*})---\left(1\right)$

其中，d为沿着射线方向的每一个采样点X_P的所有可能深度，argmin_d为d的参数极小化函数，min为值极小化函数，为射线的方向，dist(a,b)为两点的距离函数，p*为特征区域R中任意一点的坐标；

3)采用Levenberg-Marquardt方法对以下公式2进行优化，得到形变速度场参数X，C，A，f，f₀：

$\underset{X,C,A,f,f_0}{\arg \min }\mathrm{\Θ}{+\mathrm{\ω}}_r{E}_{\mathrm{reg}}+{\mathrm{\ω}}_l{E}_{\mathrm{lap}}---\left(2\right)$

其中，Θ为用户控制匹配能量，E_reg为施加外力极小化能量，E_lap为形变光滑能量，ω_r为外力极小化权重，ω_l为形变光滑权重，C为速度模版中心，A为受力方向，f为优化得到的衰减系数，f₀为用户指定的衰减系数；

上述公式2中，用户控制匹配能量Θ采用如下的公式3进行计算：

D^-1(X)=Adv^N(X,-V)

其中，X为流体模拟区域中任意一点位置坐标，为用户指定的某个特征区域的距离函数，Adv为显式欧拉积分器，V为形变速度场，D^-1为逆向形变函数，N为移动次数，三维控制点集P为步骤2)中得到的三维三角形面片的三维控制点集P_M或者二维曲线的三维控制点集P_S；

4)将优化公式2得到的形变速度场参数X，C，A，f，f₀代入以下公式4得到形变速度场V，将形变速度场V的每个速度模版中心C沿着关键帧K的速度场移动w个时间步长得到修改后的速度场V′，使用修改后的速度场V′对关键帧K进行形变得到形变后的关键帧K′，其具体公式为K′=D(X)，其中w为用户设定的序列长度，D形变函数;

$\left\{\begin{array}{c}V=\▿\×\underset{i\∈N}{\mathrm{\Σ}}{H}_i(M,C,A,f,f_0)\\ V^{\′}=\▿\×\underset{i\∈N}{\mathrm{\Σ}}{H}_i(M,{\mathrm{Adv}}^w(C,-V_K^{*}),A,f,f_0)\end{array}\right.---\left(4\right)$

其中，M为任意控制点，i为速度势函数的编号，H_i为速度势函数，为关键帧K的速度场，为梯度算子；

5)初始流体动画序列中根据用户设定的序列长度w确定子序列[K-w，K]，将子序列[K-w，K]中每个关键帧j的流体形状沿着公式5形变，得到形变后的目标子序列：

F(1-j／w)V′(5)

其中，F为平滑插值函数；

6)根据步骤4)得到的形变后的目标子序列，施加如下公式6中的外力f构造PD控制器，用PD控制器生成最终子序列：

f=f_shape+f_velocity(6)

其中，f_shape为形状匹配受力，f_velocity为速度匹配受力，f_shape和f_velocity均为无散梯度场；

7)将步骤5)中生成的最终子序列与第K-w帧之前的初始流体动画的序列拼接，生成拼接后的流体动画；

8)将拼接后的流体动画依次重复上述步骤1)～步骤7)进行不同关键帧的编辑，得到符合用户要求的最终流体动画。

所述的初始流体动画序列由数值模拟方法生成。

所述的对流体形状进行平滑化处理为将流体形状在规则网格上进行隐式转换，再用半径为ε的圆对隐式形状依次做形态学的腐蚀、膨胀操作。

所述步骤1)中的三种特征区域R分别为自由液面S、中轴区域M和边界区域B，同时满足形状距离场φ(x)＜ε和边界距离场φ_S(x)＞ε的区域为自由液面S，形状距离场φ局部极小值附近ε半径内的区域为中轴区域M，满足边界距离场φ_S(x)<ε的区域为边界区域B，其中x为流体模拟区域中的任意一点坐标，ε为用户设定的流体平滑程度参数。

所述步骤3)公式2中的施加外力极小化能量E_reg和形变光滑能量E_lap采用如下公式7进行计算：

$\left\{\begin{array}{c}{E}_{\mathrm{reg}}=\underset{M\&Element;P}{\mathrm{\&Sigma;}}{\left|\right|D^{-1}(M-M)\left|\right|}^2\\ E_{\mathrm{lap}}=\underset{M\&Element;P}{\mathrm{\&Sigma;}}\underset{Y\&Element;\mathrm{NX}}{\mathrm{\&Sigma;}}{\left|\right|D^{-1}\left(M\right)-D^{-1}\left(Y\right)\left|\right|}^2\end{array}\right.---\left(7\right)$

其中，Y为控制点M的相邻控制点集中的任意一点，NM为控制点M的相邻控制点集，对于三维三角形面片，控制点M的相邻控制点集NM采用一环临域；对于二维曲线，控制点M的相邻控制点集NM采用与控制点M相邻的两个控制点。

所述步骤4)公式4中的速度势函数H_i采取如下公式8中的旋度速度势函数H_C或平动速度势函数H_W进行计算：

$\left\{\begin{array}{c}{H}_C(M,C,A,f,f_0)=\mathrm{Aexp}(-(f^2+f_0)||M-C|\left|\right)\\ H_W(M,C,A,f,f_0)=A\&times;(M-C)\exp (-(f^2+f_0)||M-C|\left|\right)\end{array}\right.---\left(8\right)$

其中，H_C为旋度速度势函数，H_W为平动速度势函数。

所述步骤6)公式6中的形状匹配受力f_shape是由通过求解以如下公式9为边值条件的拉普拉斯方程得到：

其中，C_shape为用户可调节的参数，为目标子序列的第k帧的流体形状，为距离函数，为目标子序列的第k帧的距离函数，为修改后的第k帧的距离函数。

所述步骤6)公式6中的速度匹配受力f_velocity的表达式为如下公式10：

$f_{\mathrm{velocity}}\left(X\right)=C_{\mathrm{velocity}}(V_{{\stackrel{\&OverBar;}{F}}_k}-V_{F_k})---\left(10\right)$

其中，C_velocity为用户可调节的参数，为目标子序列的第k帧的速度场，而为修改后的第k帧的速度场。

所述的第k帧的目标速度场为一无散速度场，由通过求解以如下公式11为边值条件的拉普拉斯方程得到：

其中，为目标子序列的第k-1帧的距离函数。

所述的修改后的第k帧的初始速度场为另一无散速度场，由通过求解以如下公式12为边值条件的拉普拉斯方程得到：

其中，为初始流体动画序列的第k帧的速度场。

本发明的有益效果是：

本发明能对流体动画进行任意多次的编辑，使得其形状在指定的关键帧处符合用户施加的控制。

本发明使得对流体动画的编辑操作变得可行且直观，并提出了多种用户交互手段完成编辑操作。

附图说明

图1是本发明方法的流程图。

图2是初始流体动画序列的关键帧K处的流体形状示意图。

图3是平滑化处理后的流体形状示意图。

图4是流体形状抽取特征区域后的示意图。

图5是通过二维曲线编辑前的流体形状关键帧的示意图。

图6是通过二维曲线编辑后的流体形状关键帧的示意图。

图7是通过三维三角形面片编辑前的流体形状关键帧的示意图。

图8是通过三维三角形面片编辑后的流体形状关键帧的示意图。

图9是用户输入的二维曲线在不同特征区域下的三种类型的示意图。

图中：1、自由液面曲线，2、中轴区域曲线，3、边界区域曲线。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

如图1所示，本发明方法的步骤如下：

1)输入初始流体动画序列，指定一个需要编辑的关键帧K，对关键帧K处的流体形状进行预处理，首先对流体形状进行平滑化处理，将关键帧K处的流体形状采用隐式函数表示，然后抽取形状距离场φ和边界距离场φ_S，将流体形状区域分为三种特征区域R；

2)用户从流体形状的任意视角输入一个三维三角形面片，或者在屏幕上绘制一条二维曲线：对于用户输入的三维三角形面片，直接采样三维三角形面片的网格顶点，得到三维三角形面片的三维控制点集P_M；对于用户绘制的二维曲线，在屏幕上采样一组采样点，将每个采样点X_P与绘制视角的视点位置e连成射线，并建立如下公式1计算二维曲线的深度，即每个采样点X_P的屏幕深度d_S，然后将e+D_Xd_S作为二维曲线的三维控制点集P_S：

$d_S=\arg \mathrm{mi}{n}_d\underset{X_P\&Element;P}{\mathrm{\&Sigma;}}\underset{p^{*}\&Element;R}{\min }\mathrm{dist}(e+{D_X}_{p}d,p^{*})---\left(1\right)$

其中，d为沿着射线方向的每一个采样点X_P的所有可能深度，argmin_d为d的参数极小化函数，min为值极小化函数，为射线的方向，dist(a,b)为两点的距离函数，p^*为特征区域R中任意一点的坐标；

3)采用Levenberg-Marquardt方法对以下公式2进行优化，得到形变速度场参数X，C，A，f，f₀：

$\underset{X,C,A,f,f_0}{\arg \min }\mathrm{\&Theta;}+{\mathrm{\&omega;}}_r{E}_{\mathrm{reg}}+{\mathrm{\&omega;}}_l{E}_{\mathrm{lap}}---\left(2\right)$

其中，Θ为用户控制匹配能量，E_reg为施加外力极小化能量，E_lap为形变光滑能量，ω_r为外力极小化权重，ω_l为形变光滑权重，C为速度模版中心，A为受力方向，f为优化得到的衰减系数，f₀为用户指定的衰减系数；

上述公式2中，用户控制匹配能量Θ采用如下的公式3进行计算：

(3)

D^-1(X)=Adv^N(X，-V)

其中，X为流体模拟区域中任意一点位置坐标，为用户指定的某个特征区域的距离函数，Adv为显式欧拉积分器，V为形变速度场，D^-1为逆向形变函数，N为移动次数，三维控制点集P为步骤2)中得到的三维三角形面片的三维控制点集P_M或者二维曲线的三维控制点集P_S；

4)将优化公式2得到的形变速度场参数X，C，A，f，f₀代入以下公式4得到形变速度场V，将形变速度场V的每个速度模版中心C沿着关键帧K的速度场移动w个时间步长得到修改后的速度场V′，使用修改后的速度场V′对关键帧K进行形变得到形变后的关键帧K′，其具体公式为K′：D(X)，其中w为用户设定的序列长度，D形变函数：

$\left\{\begin{array}{c}V=\&dtri;\&times;\underset{i\&Element;N}{\mathrm{\&Sigma;}}{H}_i(M,C,A,f,f_0)\\ V^{\&prime;}=\&dtri;\&times;\underset{i\&Element;N}{\mathrm{\&Sigma;}}{H}_i(M,{\mathrm{Adv}}^w(C,-V_K^{*}),A,f,f_0)\end{array}\right.---\left(4\right)$

其中，M为任意控制点，i为速度势函数的编号，H_i为速度势函数，为关键帧K的速度场，为梯度算子；

5)初始流体动画序列中根据用户设定的序列长度w确定子序列[K-w，K]，将子序列[K-w，K]中每个关键帧j的流体形状沿着公式5形变，得到形变后的目标子序列：

F(1-j/w)V′(5)

其中，F为平滑插值函数；

6)根据步骤4)得到的形变后的目标子序列，施加如下公式6中的外力f构造PD控制器，用PD控制器生成最终子序列：

f=f_shape+f_velocity(6)

其中，f_shape为形状匹配受力，f_velocity为速度匹配受力，f_shape和f_velocity均为无散梯度场；

7)将步骤5)中生成的最终子序列与第K-w帧之前的初始流体动画的序列拼接，生成拼接后的流体动画；第K帧之后的初始流体动画的序列舍弃不用；

8)将拼接后的流体动画依次重复上述步骤1)～步骤7)进行不同关键帧的编辑，得到符合用户要求的最终流体动画。

所述的初始流体动画序列由数值模拟方法生成。

所述的对流体形状进行平滑化处理为将流体形状在规则网格上进行隐式转换，再用半径为ε的圆对隐式形状依次做形态学的腐蚀、膨胀操作。

所述步骤1)中的三种特征区域R分别为自由液面S、中轴区域M和边界区域B，同时满足形状距离场φ(x)＜ε和边界距离场φ_S(x)＞ε的区域为自由液面S，形状距离场φ局部极小值附近ε半径内的区域为中轴区域M，满足边界距离场φ_S(x)＜ε的区域为边界区域B，其中x为流体模拟区域中的任意一点坐标，ε为用户设定的流体平滑程度参数。

所述步骤3)公式2中的施加外力极小化能量E_reg和形变光滑能量E_lap采用如下公式7进行计算：

$\left\{\begin{array}{c}{E}_{\mathrm{reg}}=\underset{M\&Element;P}{\mathrm{\&Sigma;}}{\left|\right|D^{-1}\left(M\right)-M\left|\right|}^2\\ E_{\mathrm{lap}}=\underset{M\&Element;P}{\mathrm{\&Sigma;}}\underset{Y\&Element;\mathrm{NX}}{\mathrm{\&Sigma;}}{\left|\right|D^{-1}\left(M\right)-D^{-1}\left(Y\right)\left|\right|}^2\end{array}\right.---\left(7\right)$

其中，Y为控制点M的相邻控制点集中的任意一点，NM为控制点M的相邻控制点集，对于三维三角形面片，控制点M的相邻控制点集NM采用一环临域；对于二维曲线，控制点M的相邻控制点集NM采用与控制点M相邻的两个控制点。

所述步骤4)公式4中的速度势函数H_i采取如下公式8中的旋度速度势函数H_C或平动速度势函数H_W进行计算：

$\left\{\begin{array}{c}{H}_C(M,C,A,f,f_0)=\mathrm{Aexp}(-(f^2+f_0)||M-C|\left|\right)\\ H_W(M,C,A,f,f_0)=A\&times;(M-C)\exp (-(f^2+f_0)||M-C|\left|\right)\end{array}\right.---\left(8\right)$

其中，H_C为旋度速度势函数，H_W为平动速度势函数。

所述步骤6)公式6中的形状匹配受力f_shape是由通过求解以如下公式9为边值条件的拉普拉斯方程(Laplace方程)得到：

其中，C_shape为用户可调节的参数，为目标子序列的第k帧的流体形状，为距离函数，为目标子序列的第k帧的距离函数，为修改后的第k帧的距离函数。

所述步骤6)公式6中的速度匹配受力f_velocity的表达式为如下公式10：

$f_{\mathrm{velocity}}\left(X\right)=C_{\mathrm{velocity}}(V_{{\stackrel{\&OverBar;}{F}}_k}-{V_F}_k)---\left(10\right)$

其中，C_velocity为用户可调节的参数，为目标子序列的第k帧的速度场，而为修改后的第k帧的速度场。

所述的第k帧的目标速度场为一无散速度场，由通过求解以如下公式11为边值条件的拉普拉斯方程(Laplace方程)得到：

其中，为目标子序列的第k-1帧的距离函数。

所述的修改后的第k帧的初始速度场为另一无散速度场，由通过求解以如下公式12为边值条件的拉普拉斯方程(Laplace方程)得到：

其中，为初始流体动画序列的第k帧的速度场。

所述的步骤1)中对流体形状进行平滑化处理为将流体形状在规则网格上进行隐式转换，再用半径为ε的圆对隐式形状依次做形态学的腐蚀Erode、膨胀Dilate操作，可参考[GonzalezRC，WoodsRE，EddinsSL.DigitalimageprocessingusingMATLAB[M].Knoxville：GatesmarkPublishing，2009.]；

所述的步骤1)初始流体动画序列可由多种数值模拟方法生成。常用的数值模拟方法有：第一数值模拟方法可参考[Zhu，Y.，andBridson，R.2005.Animatingsandasafluid.ACMTransactiononGraphics(Proc.SIGGRAPH)24，3，965-972]；第二数值模拟方法可参考[Foster，N.，ANDFedkiw，R.2001.Practicalanimationofliquids.ACMTransactiononGraphics(Proc.SIGGRAPH)，23-30]。

对于上述两种数值模拟方法得到的初始流体动画序列，特征区域R被表示为不同的形式：

a)对于采用第一数值模拟方法得到的初始流体动画序列，特征区域被表示为流体粒子的集合。该粒子集合通过一个隐式转换最终被表示为距离场。

b)对于采用第二数值模拟方法得到的初始流体动画序列，特征区域被直接表示为距离场。

本发明中采用的隐式转换方法均为[Zhu，Y.，andBridson，R.2005.Animatingsandasafluid.ACMTransactiononGraphics(Proc.SIGGRAPH)24，3，965-972]中提出的方法。

所述步骤1)中将流体形状区域分为三种特征区域R，其目的在于允许用户对关键帧K处的流体形状进行三种不同的控制。如图9所示，本发明允许用户输入的二维曲线有：自由液面曲线1，中轴区域曲线2以及边界区域曲线3三种。

步骤2)的公式1中：为最邻近流体距离，最邻近流体距离通过如下公式13计算：

步骤2)之前需对于二维曲线的三维控制点集P_S进行重新采样，使得P_S中点的间距等于规则网格单元的宽度。

步骤3)中公式2)中的施加外力极小化能量E_rep，形变光滑能量E_lap根据流体的形变应当尽量平滑，且形变不应太过剧烈的要求建立。

步骤3)采用Levenberg-Marquardt计算方法对公式2进行优化，该计算方法可参考[DonaldW.Marquardt.AnAlgorithmforLeast-SquaresEstimationofNonlinearParameters，JournaloftheSocietyforIndustrialandAppliedMathematics，Vol.11，No.2June，1963，431-441.]。执行Levenberg-Marquardt计算方法的具体细节为：

a)执行过程中梯度的计算采用如下公式14的两阶段链式求导计算方法：

$\frac{\&PartialD;(\mathrm{\&Theta;}+{\mathrm{\&omega;}}_r{E}_{\mathrm{reg}}+{\mathrm{\&omega;}}_l{E}_{\mathrm{lap}})}{\&PartialD;X}=\frac{\&PartialD;(\mathrm{\&Theta;}+{\mathrm{\&omega;}}_r{E}_{\mathrm{reg}}+{\mathrm{\&omega;}}_l+E_{\mathrm{lap}})}{{\&PartialD;D}^{-1}\left(X\right)}\frac{{\&PartialD;D}^{-1}\left(X\right)}{X}---\left(14\right)$

其中，为偏导记号，对于右侧第二项的梯度计算采用Adjoint计算方法，Adjoint计算方法可参考[MBGiles，NAPierce.Anintroductiontotheadjointapproachtodesign.Flow，turbulenceandcombustion，2000.]。优化过程中不断执行Levenberg-Marquardt计算方法进行迭代，直到相对或绝对误差减小至1e^-3后终止。

b)针对右侧第二项梯度的计算采用图形硬件对其进行加速，每个Levenberg-Marquardt计算方法迭代过程中将形变速度场参数X，C，A，f，f₀传至图形硬件，并将结果传回内存。在图形硬件上，使用单独线程，计算每个控制点X对D^-1(X)的梯度；

步骤3)中公式2根据用户绘制或者输入后的三维控制点集、流体边界条件以及流体不可压缩性质建立。

将步骤5)中生成的最终子序列与第K-w帧之前的初始流体动画的序列拼接，生成编辑后的流体动画，而第K帧之后的初始流体动画的序列被标记为不可用后删除。

步骤6)中PD控制器可参考[LinShiandYizhouYu，TamingLiquidsforRapidlyChangingTargets，ACMSIGGRAPH/EurographicsSymposiumonComputerAnimation，LosAngeles，July2005，pp.229-236.]。

步骤4)中公式4根据流体的不可压缩性质，以及固定边界条件而建立。其目的为将流体关键帧的可行形变限制为不可压缩的，且满足固定边界条件的。

步骤4)中，若存在固定边界条件，则H_i采取如下公式15：

H_i(M)＝H(M)R(φ(M)/γ)+(1-R(φ(M)/γ))N(M)(H(M)·N(M))(15)

其中R(·)为如下公式16中的截断函数，γ为截断参数，N为距离场梯度：

$R\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;1\\ \frac{15}{8}x-\frac{10}{8}{x}^3+\frac{3}{8}{x}^5& 1>x>-1\\ -1& x\&le;-1\end{array}\right.---\left(16\right)$

本发明所述的可控性体现为：允许用户对一段输入的初始流体动画序列进行控制，以及通过输入二维曲线或三维三角面片的方法，修改关键帧K处的流体形状。

本发明实施例：

a)本实施例中所使用的参数为：N＝10，W＝50，ω_l＝5.0，ω_r＝0.0以及f₀＝3。

b)输入的流体关键帧形状如图2所示，平滑化处理后的流体形状如图3所示，以及抽取的特征区域如图4所示。在本实施例使用的ε＝0.03。

c)通过控制曲线编辑关键帧K，关键帧K的流体形状如图5所示。

d)形变后的关键帧K′处的流体形状如图6所示。本实施例中，步骤3)至步骤5)总共耗费的时间为5.7秒。

e)用户通过输入一个三维三角形面片，分别对一段输入的初始流体动画序列中两个需要编辑的关键帧K进行了编辑获得形变后的关键帧K′，形变前后变化分别如图7，图8所示。

上述具体实施方式用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。

Claims (10)

1.一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于包含如下步骤：

1)输入初始流体动画序列，指定一个需要编辑的关键帧K，对关键帧K处的流体形状进行预处理，首先对流体形状进行平滑化处理，将关键帧K处的流体形状采用隐式函数表示，然后抽取形状距离场φ和边界距离场φ_S，将流体形状区域分为三种特征区域R；

2)用户从流体形状的任意视角输入一个三维三角形面片，或者在屏幕上绘制一条二维曲线：对于用户输入的三维三角形面片，直接采样三维三角形面片的网格顶点，得到三维三角形面片的三维控制点集P_M；对于用户绘制的二维曲线，在屏幕上采样一组采样点，将每个采样点X_P与绘制视角的视点位置e连成射线，并建立如下公式1计算二维曲线的深度，即每个采样点X_P的屏幕深度d_S，然后将e+D_Xd_S作为二维曲线的三维控制点集P_S：

$d_S={\mathrm{argmin}}_d\underset{X_P\&Element;P}{\&Sigma;}\underset{p*\&Element;R}{min}dist(e+D_{X_p}d,p*)---\left(1\right)$

其中，d为沿着射线方向的每一个采样点X_P的所有可能深度，argmin_d为d的参数极小化函数，min为值极小化函数，为射线的方向，dist(a,b)为两点的距离函数，p^*为特征区域R中任意一点的坐标；

3)采用Levenberg-Marquardt方法对以下公式2进行优化，得到形变速度场参数X,C,A,f,f₀：

$\underset{X,C,A,f,f_0}{\arg \min }\mathrm{\&Theta;}+{\mathrm{\&omega;}}_r{E}_{reg}+{\mathrm{\&omega;}}_l{E}_{lap}---\left(2\right)$

其中，Θ为用户控制匹配能量，E_reg为施加外力极小化能量，E_lap为形变光滑能量，ω_r为外力极小化权重，ω_l为形变光滑权重，C为速度模版中心，A为受力方向，f为优化得到的衰减系数，f₀为用户指定的衰减系数；

上述公式2中，用户控制匹配能量Θ采用如下的公式3进行计算：

其中，X为流体模拟区域中任意一点位置坐标，为用户指定的某个特征区域的距离函数，Adv为显式欧拉积分器，V为形变速度场，D^-1为逆向形变函数，N为移动次数，三维控制点集P为步骤2)中得到的三维三角形面片的三维控制点集P_M或者二维曲线的三维控制点集P_S；

4)将优化公式2得到的形变速度场参数X,C,A,f,f₀代入以下公式4得到形变速度场V，将形变速度场V的每个速度模版中心C沿着关键帧K的速度场移动w个时间步长得到修改后的速度场V'，使用修改后的速度场V'对关键帧K进行形变得到形变后的关键帧K'，其具体公式为K'＝D(X)，其中w为用户设定的序列长度，D为形变函数；

$\{\begin{array}{c}V=\&dtri;\&times;\underset{i\&Element;N}{\mathrm{\&Sigma;}}{H}_i(M,C,A,f,f_0)\\ V^{\&prime;}=\&dtri;\&times;\underset{i\&Element;N}{\mathrm{\&Sigma;}}{H}_i(M,{\mathrm{Adv}}^w(C,-V_K^{*}),A,f,f_0)\end{array}---\left(4\right)$

其中，M为任意控制点，i为速度势函数的编号，H_i为速度势函数，为关键帧K的速度场，为梯度算子；

5)初始流体动画序列中根据用户设定的序列长度w确定子序列[K-w,K]，将子序列[K-w,K]中每个关键帧j的流体形状沿着公式5形变，得到形变后的目标子序列：

F(1-j/w)V'(5)

其中，F为平滑插值函数；

6)根据步骤4)得到的形变后的目标子序列，施加如下公式6中的外力f构造PD控制器，用PD控制器生成最终子序列：

f＝f_shape+f_velocity(6)

其中，f_shape为形状匹配受力，f_velocity为速度匹配受力，f_shape和f_velocity均为无散梯度场；

7)将步骤5)中生成的最终子序列与第K-w帧之前的初始流体动画的序列拼接，生成拼接后的流体动画；

8)将拼接后的流体动画依次重复上述步骤1)～步骤7)进行不同关键帧的编辑，得到符合用户要求的最终流体动画。

2.根据权利要求1所述的一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于：所述的初始流体动画序列由数值模拟方法生成。

3.根据权利要求1所述的一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于：所述的对流体形状进行平滑化处理为将流体形状在规则网格上进行隐式转换，再用半径为ε的圆对隐式形状依次做形态学的腐蚀、膨胀操作。

4.根据权利要求1所述的一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于：所述步骤1)中的三种特征区域R分别为自由液面S、中轴区域M和边界区域B，同时满足形状距离场φ(x)<ε和边界距离场φ_S(x)>ε的区域为自由液面S，形状距离场φ局部极小值附近ε半径内的区域为中轴区域M，满足边界距离场φ_S(x)<ε的区域为边界区域B，其中x为流体模拟区域中的任意一点坐标，ε为用户设定的流体平滑程度参数。

5.根据权利要求1所述的一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于：所述步骤3)公式2中的施加外力极小化能量E_reg和形变光滑能量E_lap采用如下公式7进行计算：

$\left\{\begin{array}{c}{E}_{reg}=\underset{M\&Element;P}{\&Sigma;}\left|\right|D^{-1}\left(M\right)-M|{|}^2\\ E_{lap}=\underset{M\&Element;P}{\&Sigma;}\underset{Y\&Element;NX}{\&Sigma;}\left|\right|D^{-1}\left(M\right)-D^{-1}\left(Y\right)|{|}^2\end{array}\right.---\left(7\right)$

其中，Y为控制点M的相邻控制点集中的任意一点，NM为控制点M的相邻控制点集，对于三维三角形面片，控制点M的相邻控制点集NM采用一环临域；对于二维曲线，控制点M的相邻控制点集NM采用与控制点M相邻的两个控制点。

6.根据权利要求1所述的一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于：所述步骤4)公式4中的速度势函数H_i采取如下公式8中的旋度速度势函数H_C或平动速度势函数H_W进行计算：

$\left\{\begin{array}{c}{H}_C(M,C,A,f,f_0)=A\exp (-(f^2+f_0)\left|\right|M-C\left|\right|)\\ H_W(M,C,A,f,f_0)=A\&times;(M-C)\exp (-(f^2+f_0)\left|\right|M-C\left|\right|)\end{array}\right.---\left(8\right)$

其中，H_C为旋度速度势函数，H_W为平动速度势函数。

7.根据权利要求1所述的一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于：所述步骤6)公式6中的形状匹配受力f_shape是由通过求解以如下公式9为边值条件的拉普拉斯方程得到：

其中，C_shape为用户可调节的参数，为目标子序列的第k帧的流体形状，为距离函数，为目标子序列的第k帧的距离函数，为修改后的第k帧的距离函数。

8.根据权利要求1所述的一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于：所述步骤6)公式6中的速度匹配受力f_velocity的表达式为如下公式10：

$f_{velocity}\left(X\right)=C_{velocity}(V_{{\stackrel{\&OverBar;}{F}}_k}-V_{F_k})---\left(10\right)$

其中，C_velocity为用户可调节的参数，为目标子序列的第k帧的速度场，而为修改后的第k帧的速度场。

9.根据权利要求8所述的一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于：所述的目标子序列的第k帧的速度场为一无散速度场，由通过求解以如下公式11为边值条件的拉普拉斯方程得到：

其中，为目标子序列的第k-1帧的距离函数。

10.根据权利要求8所述的一种基于流体关键帧编辑的可控流体动画生成方法，其特征在于：所述的修改后的第k帧的速度场为另一无散速度场，由通过求解以如下公式12为边值条件的拉普拉斯方程得到：

其中，为初始流体动画序列的第k帧的速度场。