CN105279781B - 基于多精度融合的流体动画生成方法 - Google Patents

Abstract

本发明公开了基于多精度融合的流体动画生成方法，包括以下步骤：从流体计算及流体表面构建两个方面自底向上建立多精度流体动画的生成模型。首先，在底层利用自适应SPH模型计算流体方程并获取多精度粒子集。其次，提取表面粒子集并对其进行预处理。然后，在上层建立多精度表面模型。最后，通过建立多精度网格、计算无符号概率场、构建带权图并利用图切算法获得具有不同精度的流体表面，然后将其融合在一起形成最终的流体表面。最后通过使用真实感渲染方法获取对应的流体动画效果。该方法能够在保持流体动画真实感的前提下，可实现对流体动画的快速、高效的生成。

Description

基于多精度融合的流体动画生成方法

技术领域

本发明涉及一种基于多精度融合的流体动画生成方法，该方法考虑了流体运动视觉上的多尺度特征，分别从流体计算及流体表面构建两个方面自底向上建立多精度流体动画的生成模型，从而在保持流体动画真实感的前提下，实现对流体动画的快速、高效的生成。

背景技术

流体运动容易产生丰富的、具有不同精度的细节，例如广阔的水面和细小的水花等，这些细节的模拟和渲染对于流体动画的视觉真实感至关重要。然而，细节的生成需要高精度的计算和渲染，由于需要处理大规模的粒子，因此无论是在流体计算还是在流体表面构建等方面都会产生极大的计算消耗。现有的方法并没有利用流体运动在视觉上所产生的多精度特征，很难在有限计算资源的限制下进行有效的流体计算和表面生成，严重影响了动画的视觉真实感。

发明内容

本发明的目的就是为了解决上述问题，提供一种基于多精度融合的流体动画生成方法。通过将不同精度的表面融合在一起形成最终的流体表面，结合使用真实感渲染方法即可获取对应的流体动画效果。本发明方法能够在生成具有多种细节的、视觉逼真的流体动画的同时有效降低计算量。相关方法可广泛应用于大规模流体动画的生成、动漫制作、影视特效、游戏娱乐等领域中。

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

基于多精度融合的流体动画生成方法，包括如下步骤：

步骤(1)：初始化NS方程的初始条件和边界条件；所述初始条件为流体源S，流体源S由SPH粒子集P表示；所述边界条件B包括固体与液体边界、空气和液体边界；

步骤(2)：根据初始化的流体源S及其对应的粒子集P，基于NS方程，建立基于SPH的多精度流体计算方法，获得多精度粒子集；

步骤(3)：对步骤(2)的多精度粒子集进行预处理，将处理后的粒子集使用八叉树存储；

步骤(4)：从八叉树顶层开始不断重复迭代执行步骤(5)直到八叉树底层；

步骤(5)：对八叉树的每一层，计算并提取该层对应精度的流体表面网格；

步骤(6)：将不同精度的流体表面网格融合成一个流体表面，使用真实感渲染获得每帧的流体动画；

步骤(7)：如果当前帧数小于用户设定的动画帧数，则转步骤(2)继续下一帧流体动画的计算，最终得到体现多尺度细节的流体动画。

NS(Navier-Stokes，纳维-斯托克斯)

SPH(Smooth Particle Hydrodynamics,光滑粒子流体动力学)

所述场景模型包括二维或三维场景模型；

所述步骤(1)初始化NS方程的初始条件就是设定粒子集P中每个粒子的半径大小、位置和速度；

所述步骤(1)初始化边界条件B就是根据输入的流体所在的场景模型和流体源S，设定固体和液体边界；在剩余空间和流体源S之间设定空气和液体边界；

所述步骤(1)的步骤如下：

步骤(1-1)：初始化SPH粒子集P；

每个粒子具有以下属性：质量m，半径r，所在层数l，速度u，位置x。初始化时，粒子被随机的分到不同的层中l＝0,1,2,…。第0层的粒子质量和半径是最小的，令第0层粒子的质量为m⁰，半径为r⁰，那么位于第l层的质量和半径分别为m＝2^lm⁰和

根据用户输入的初始化速度向量v来初始化粒子的速度u＝v，粒子的位置x是在用户给定的流体源的内部随机生成。粒子之间的距离为r_i/h，取h＝2.5。如果两个粒子i和j满足||x_i-x_j||≤max(r_i,r_j)，则说明粒子i和j是邻居关系。将粒子用KD树存储，则该邻居关系可以通过距离查询快速获取。

用户给定的流体源是指，初始化的水源。例如模拟一个水球落入盆中，则初始化的水球即为源，它包括水球的位置和初始速度。在本发明中，用填充满水球的粒子集表示这个水源。

步骤(1-2)：初始化流体运动的边界条件；

将流体运动所在的环境模型利用体素化方法转化为NS方程计算时可识别的边界。该边界包含两个部分：固体与液体边界、空气和液体边界。本发明使用一个标记位标识不同的边界。

所述步骤(2)的多精度粒子集是指由半径有大小区别的粒子组成的集合；

所述步骤(2)的步骤如下：

步骤(2-1)：定义并计算多精度粒子半径的自适应分布函数；

步骤(2-1-1)：定义并计算由多精度粒子集所形成的二维或三维形状的中轴粒子集；

给定两个邻居粒子i和j以及精度控制常数c，如果dis(i,k)＜c且dis(j,q)＜c，其中粒子k和q分别是距离i和j最近的表面粒子，若满足

则粒子i和j为中轴粒子集中的粒子。其中，γ表示角度的阈值，取值为γ＝60；y_i为粒子k的位置，y_j为粒子q的位置。公式(1)中第一个条件表明两个互为邻居的粒子i和j位于中轴两侧，而第二个条件则说明k和q的距离需要大于互为邻居的两个粒子i和j之间的距离。如果一个粒子没有邻居，则该粒子也是中轴粒子集上的粒子。另外，dis(i,k)＜c和dis(j,q)＜c表明本发明定义的中轴上的粒子与流体表面粒子的距离小于给定常数c。c越大，中轴粒子越多；c越小，中轴粒子越少。

步骤(2-1-2)：计算流体表面M的粒子y对应的多精度函数值为MRlfs(y)＝min(||p_M-y||,c)。其中，p_M为中轴上距离y最近的粒子,c为定义中轴时的常数。

步骤(2-1-3)：将流体表面粒子y对应的自适应多精度函数值拓展到流体V的内部。表示为

即流体V内的点x位置处的局部特征尺寸MRelfs定义为体内的点到流形表面的距离与表面上该点MRlfs之和的最小值，其中，y是位于表面上的粒子。当计算获得表面粒子的MRlfs之后，应用快速步进算法获得流体V内部的粒子的局部特征尺寸函数MRelfs的值。显然，MRelfs值越小，说明点x位置处越细薄，正是流体细节较多的位置，需要粒子的半径较小；反之，MRelfs值越大，需要粒子的半径较大。

步骤(2-2)：根据粒子半径的自适应分布函数对粒子进行分裂和合并，获得多精度粒子集；

步骤(2-2-1)：判断SPH粒子集P是需要分裂、合并还是无需变化。如果MRelfs(x_i)＜αr_i则需要分裂；如果MRelfs(x_i)＞βr_i则需要合并；否则粒子无需变化；

步骤(2-2-2)：粒子分裂。对于满足分裂条件的粒子在分裂时，位于l层粒子i被分解为两个位于l-1层上的两个粒子j和k。为了保证两个粒子位于流体内部，并且与粒子i的距离最近，防止引入较大的压力作用，这两个粒子将在i的周围对称性分布，距离i为d＝r_j/(2h)＝r_k/(2h)，h＝2.5。如果存在很多可用位置，则随机选择一个。新粒子的半径和质量减半，且速度与粒子i相同。

步骤(2-2-3)：粒子合并。对于满足合并条件的粒子i在合并时，需要寻找在同一层上另一个待合并的邻居粒子j，产生新的粒子k，其层为l_k＝l_i+1＝l_j+1，半径和质量为合并前的2倍，位置为x_k＝(x_i+x_j)/2。如果该位置在流体内部且距离粒子k的半径为r_k/(2h)区域内没有其他粒子，则合并成功，h＝2.5。否则，粒子i需要寻找其他的待合并的邻居粒子合并。如果没有一个邻居粒子满足条件，则粒子i的合并就被取消并推迟到下一时间步中。

步骤(2-3)：对多精度粒子集，计算流体动力学方程，更新其速度和位置。

步骤(2-3-1)：建立基于拉格朗日法的流体动力学模型。为了计算n+1时刻粒子的动态流体速度，需要首先求解流体运动的NS方程

其中，ρ为密度，u为流体速度，p为压强，t为时间，f为外力，μ为运动粘性系数，μ▽²u表示的是粘性力，▽·为散度，▽为梯度，▽²＝▽·▽。在SPH模型中，由于粒子个数确定且质量都是常数，满足质量守恒，因此质量守恒方程(3)可以被忽略。其次，公式(4)中左侧表达式可以使用du/dt来表示，即拉格朗日法中对流项u·(▽u)可以方便的处理。因此，方程(4)变为

令f＝-▽p+ρg+μ▽²u，对于粒子i来说，根据牛顿第二定律

其中，a_i表示粒子的加速度，f_i和ρ_i分别是粒子i所在位置的力及密度。

步骤(2-3-2)：利用SPH方法估计密度ρ_i。SPH是粒子系统的一种插值方法，利用该方法，可从粒子携带的量中获取该变量在流体场中任何位置中的值。在SPH方法中，变量X在位置r处的值可根据其邻居粒子的加权和得到，即

其中，j表示粒子，m_j表示j的质量，x_j为j的位置，ρ_j为j的密度，X_j为粒子j所携带的量的值，W为各向同性的支撑半径为h的光滑的核函数

那么，粒子所在位置x_i ⁿ处的密度ρ_i为

步骤(2-3-3)：利用SPH方法估计外力f_i。这里，外力f＝-▽p+ρg+μ▽²u，故需要估计压力▽p和粘性力μ▽²u。其中压力▽p表示为

压强使用状态方程获取，即p＝k(ρ-ρ₀)，ρ₀为参考密度。为了生成对称的压力，使用(p_i+p_j)/2代替p_i。粘性力表示为

这里用u_j-u_i代替u_j其原因是粘性力与相对速度相关。

步骤(2-3-4)：更新粒子的速度和位置。粒子i的外力f_i和密度ρ_i均可被计算出来，那么粒子i在n+1时间步的新速度为

其中，Δt为时间步长。最后，根据dx_i/dt＝u_i，粒子i在n+1时间步的新位置为

所述步骤(3)的步骤如下：

步骤(3-1)：从多精度SPH粒子集P中提取表面粒子集；

步骤(3-1-1)：构建场函数，使该函数在粒子中心的取值为1，而在其他位置取0。为了保证光滑性，使用插值方法计算该场函数C_S

步骤(3-1-2)：计算每个粒子位置处的法向量。流体表面x处的法向量n可表示为

步骤(3-1-3)：根据法向量大小确定表面粒子集P_s。在流体内部，n的大小接近于零，而在流体表面位置，由于某个方向上存在粒子空隙，故n的模大于零。故给定阈值θ，如果粒子i满足||n(x_i)||≥θ，则将i并入到表面粒子集P_s中，即P_s＝P_s∪i。

步骤(3-2)：计算粒子对应的法向量、置信值和占用空间大小；

步骤(3-2-1)：计算P_s中每个粒子的法向量。对位于P_S中的粒子i，其法向量首先通过公式(15)计算，然后将其归一化，得单位向量n_i＝n(x_i)/||n(x_i)||作为该粒子的法向量。

步骤(3-2-2)：计算P_s中每个粒子的置信值。采样点的置信值c_i表示该采样点位于流体表面的概率。显然，距离流体表面越近的粒子，其置信值越大，反之则越小。为了计算P_S中的粒子i到表面的距离，需要获取距离i最近的表面粒子。使用在法向量方向上进行外推的方法，将表面粒子的编号I外推到其他粒子中，即

其中，n为表面粒子法向量，该方程可以使用快速步进法求解。这样，每个P_S中的粒子都会获取距离其表面最近的表面粒子I，故采样点i的置信值可表示为c_i＝1-α||x_i-x_I||，其中，α＞0为控制参数。

步骤(3-2-3)：计算P_c中每个粒子的占用空间大小s_i。s_i表示粒子对周围区域的影响的半径大小，为此定义为s_i＝cr_i，其中c≥1为常数。

步骤(3-2-4)：将P_c中的粒子存储在八叉树中。由于表面粒子集P_S表示计算的流体，因此包围盒B设置为整个流体计算域。将包围盒B分解为规则的体网格V，使用八叉树数据结构存储该体网格，其中层最粗糙，而最精细。多精度流体表面粒子集P_S中的粒子即采样点将根据位置被分布在不同的八叉树结点中。

所述步骤(5)的步骤如下：

步骤(5-1)：获取八叉树一层，计算包围该层所有粒子的壳；

步骤(5-1-1)：生成壳V_crust。给定八叉树的一层该层上V_crust的生成通过两个过程完成。

首先，生成表示流体表面的V_crust核。为此，将包含表面采样点的体素v的标记v_tag＝1，并将其置入到V_crust中。如果表面采样点较多，可以取较精细的层次，以尽快获得高精度的表面；反之，则取较粗糙的层次，否则很难生成封闭的V_crust。

其次，根据V_crust核中的体素单元计算位于流体表面的窄带区域上的体素单元，扩充V_crust核。具体的，执行一个标记及采样点法向外推的操作，即

其中，n为采样点的法向。根据公式(17)将当前V_crust体素往n方向膨胀τ次，膨胀的次数由用户定义，同时，将膨胀后的体素v的标记v_tag设置为距离v最近采样点的标记值。这样，V_crust即可表达流体表面及其周围的窄带。

步骤(5-1-2)：壳分割。对于体素网格集合V_crust中的每个体素v，需要确定其是位于表面的内侧还是外侧，从而实现对体素的分割。为此，需要比较体素表面的法向量与采样点的法向量。如果是流体表面的外侧，两个法向量夹角为锐角，而在流体表面的内侧，两个法向量夹角为钝角。对于体素v(v∈V_crust)来说，其表面上的法向量根据体素面的朝向容易计算，表示为对位于V_crust核上包含采样点的体素v的表面法向定义为v内采样点法向量的平均值来表示，即

对V_crust上其他体素表面法向定义为距离v最近的包含采样点的体素的法向量。在法向外推过程中，同时将体素法向量记录即可。那么，如果则体素v为外侧单元。否则体素v为内侧单元。这里ε∈(0,1)表示阈值。

步骤(5-2)：根据所得的壳，创建带权图；

步骤(5-2-1)：初始化带权图G(V,E,W)。G的顶点为体素v(v∈V_crust)的每个面f_v，即V＝{f_v}。如果在V_crust中两个面和共享一条边，则在对应的顶点之间增加一条边即E＝{e_ij}。由于V_crust中每个面有四条边，因此生成的图的每个顶点的度数为4，即规则的4连接的图。

步骤(5-2-2)：对G(V,E,W)增加源结点S和汇结点T，并在G中增加连接其他顶点到源结点S以及其他顶点汇结点T的边。具体的，对G中任意顶点f_v，如果该顶点所在体素v是V_crust的内侧且位于边界上，则增加一条连接该顶点到汇结点的边e(f_v,T)。如果该顶点所在体素v是V_crust的外侧且位于边界上，则增加一条连接该顶点到源结点的边e(f_v,S)。

步骤(5-2-3)：初始化权重。建立完图G之后，要对G中的边设置权重。对于连接两个非终点的边其权重设置为w。连接终点的边e(f_v,S)和e(f_v,T)的权重为μ，其中，w＝1，μ＝0.3。

步骤(5-3)：计算全局和局部置信值场，对带权图的权值进行变更；

步骤(5-3-1)：计算全局置信值场Γ(x)。全局置信值场为一个从R³到R的映射。该映射为V_crust中的每个体素v计算一个置信值。为此，首先定义每个采样点对周围置信值的影响。对采样点_i，其占用空间的大小为s_i，那么采样点对周围空间置信值的贡献为高斯函数

通过叠加所有采样点的空间置信值，可以得到全局置信值场。

为了提高计算效率，在具体实现中首先设定距离采样点大于其三倍s_i的值的体素，其置信度值γ_v＝0。对八叉树的每一层，计算V_crust体素中心的全局置信值Γ(x)，并将其存储在相应的八叉树节点中。当计算完一个八叉树节点的置信值之后，利用深度遍历的方式将该节点中的值加到其孩子结点中。对于采样点_i来说，如果其占用空间小于当前八叉树节点中体素的大小，则暂时增大其s_i至体素大小以计算当前的全局置信场。同时，标记该采样点，在下一层八叉树节点中进行高精度的计算。

步骤(5-3-2)：对全局置信值场Γ(x)进行自适应优化，得Γ_A(x)。由于置信值的大小要进行归一化，而最大值Γ_max可能很大，因此置信场再归一化后将导致全局置信场的局部变化很小，影响图切割算法的精度。为此，对Γ(x)进行局部自适应优化得

其中，B_edge是包围盒的边长。W(d)为权重函数

其中，D为体素中的滤波半径，D＝10，m＝7。

步骤(5-3-3)：基于优化后的置信值场Γ_A(x)对带权图G(V,E,W)的权重变更。在V_crust中的每个体素v共有六个面，而对于每个面，都有四条边与其他的面共享。因此，体素v中包含一个具有12条边的八面体图，每条边的权重为w。最佳流体表面应该是在最小割(面积最小)的基础上具有最大化的全局置信度。为此，需要在具有较高置信值的体素内设置较小的边权重，此时，修改图G(V,E,W)中2个顶点均位于体素v内的边的权重为

其中，Γ_Amax表示自适应全局置信场中的最大值，a为控制参数，调整后的权重值为w∈[a,1+a]。

步骤(5-4)：计算带权图的最小割，获得八叉树当前层对应的流体表面网格。

步骤(5-4-1)：计算流体三角形表面网格。图G(V,E,W)的最小割将产生一组切割边，这组切割边是在源结点和汇结点之间形成分割，同时也是获得的具有最高置信值的对流体表面的逼近。流体表面是在包含至少一条割边的体网格内产生。具体的，相邻的2个包含割边的体网格中心点相连，即可生成三角形的一条边，因此2×2×2的体网格内即可生成一个三角形。这样，连接相邻的、包含割边的体素中心点即可形成流体表面的三角形网格。

步骤(5-4-2)：表面网格光滑化。由于三角形的顶点只位于体素的中心，因此在获取到流体表面后，可执行一个拉普拉斯光滑化的过程，得到层上的光滑的流体表面M_i。

所述步骤(6)的步骤如下：

步骤(6-1)：将第i层得到的网格M_i融合到流体表面网格M中；

步骤(6-2)：将第i+1层产生的表面网格M_i+1加入到M中，并舍弃M_i+1所在位置的低精度三角形；

步骤(6-3)：在M_i+1与M_i的交界处，令在同一个体素内的低精度顶点代替高精度顶点完成表面的融合；

步骤(6-4)：如果i+1未达到八叉树的底层，则转步骤(6-1)继续融合；直到八叉树所有层的网格融合到M中；

步骤(6-5)：使用光子映射算法对融合得到的多精度网格进行渲染。

本发明的有益效果：

1本发明该方法通过输入流体源S、流体运动的边界B，生成多精度流体动画Frames＝{f₁,f₂,…,f_n},其中f_i表示每帧动画对应的图片，从流体运动的视觉效果出发，以多精度的方式实现对流体动画的建模和表面构建。该方法能够在有限计算资源下获取和表达流体不同尺度的流体细节，有效增强流体动画的视觉效果，提高动画生成的效率。

2本发明考虑流体运动在视觉上产生的多精度特征，自底向上的构建高精度流体动画生成模型，通过建立流体运动的多精度模型获取多精度粒子集、对粒子集进行预处理、多精度表面构建等过程，生成高精度的、视觉逼真的流体动画。

附图说明

图1是本发明模型简化的流程图。

具体实施方式

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

如图1所示，基于多精度融合的流体动画生成方法，包括以下几个过程：

过程1：初始化初始条件和边界条件；

第一步：初始化边界条件。将流体运动所在的环境模型利用体素化方法转化为方程计算时可识别的边界。它包含两个部分：固体与液体边界、空气和液体边界。本发明使用一个标记位标识不同的边界。

第二步：初始化SPH粒子集P。每个粒子具有以下属性：质量m，半径r，所在层数l，速度u，位置x。初始化时，粒子被随机的分到不同的层中l＝0,1,2,…。第0层的粒子质量和半径是最小的，令第0层粒子的质量为m⁰，半径为r⁰，那么位于第l层的质量和半径分别为m＝2^lm⁰和根据用户输入的初始化速度向量v初始化粒子的速度u＝v，粒子的位置x是在用户给定的流体源的内部随机生成。粒子之间的距离为r_i/h，取h＝2.5。如果两个粒子i和j满足||x_i-x_j||≤max(r_i,r_j)，则说明粒子i和j是邻居关系。将粒子用KD树存储，则该邻居关系可以通过距离查询快速获取。

过程2：定义多精度粒子半径的自适应分布函数；

第一步：定义并计算由多精度粒子集形成形状的中轴粒子集。给定两个邻居粒子i和j以及精度控制常数c，如果dis(i,k)＜c且dis(j,q)＜c，其中粒子k和q分别是距离i和j最近的表面粒子，若满足

则粒子i和j为中轴粒子集中的粒子。其中，γ表示角度的阈值，取值为γ＝60。公式(1)中第一个条件表明两个粒子位于中轴两侧，而第二个条件则删除接近表面的粒子。如果一个粒子没有邻居，则该粒子也是中轴上的粒子。另外，dis(i,k)＜c和dis(j,q)＜c表明本专利定义的中轴上的粒子与流体表面粒子的距离小于给定常数c。c越大，中轴粒子越多；c越小，中轴粒子越少。

第二步：计算流体表面M的粒子y对应的多精度函数值为MRlfs(y)＝min(||p_M-y||,c)。其中，p_M为中轴上距离y最近的粒子,c为定义中轴时的常数。

第三步：将流体表面粒子y对应的自适应多精度函数值拓展到流体V的内部。表示为

即V内的点x位置处的局部特征尺寸MRelfs定义为体内的点到流形表面的距离与表面上该点MRlfs之和的最小值，其中，y是位于表面上的粒子。当计算获得表面粒子的MRlfs之后，应用快速步进算法获得流体V内部的粒子的局部特征尺寸函数MRelfs的值。显然，MRelfs值越小，说明此处越细薄，正是流体细节较多的位置，需要粒子的半径较小；反之，MRelfs值越大，需要粒子的半径较大。

过程3：根据多精度粒子半径的自适应分布函数对粒子进行合并和分裂；

第一步：判断SPH粒子集P是分裂、合并还是无需变化。如果MRelfs(x_i)＜αr_i则需要分裂；如果MRelfs(x_i)＞βr_i则需要合并；否则粒子无需变化；

第二步：粒子分裂。对于满足分裂条件的粒子在分裂时，位于l层子i被分解为两个位于l-1层上的两个粒子j和k。为了保证两个粒子位于流体内部，并且与粒子i的距离最近，防止引入较大的压力作用，这两个粒子将在i的周围对称性分布，距离i为d＝r_j/(2h)＝r_k/(2h)，本发明h＝2.5。如果存在很多可用位置，则随机选择一个。新粒子的半径和质量减半，且速度与粒子i相同。

第三步：粒子合并。对于满足合并条件的粒子i在合并时，需要寻找在同一层上另一个待合并的邻居粒子j，产生新的粒子k，其层为l_k＝l_i+1＝l_j+1，半径和质量为合并前的2倍，位置为x_k＝(x_i+x_j)/2。如果该位置在流体内部且距离粒子k的半径为r_k/(2h)区域内没有其他粒子，则合并成功。本发明h＝2.5。否则，粒子i需要寻找其他的待合并的邻居粒子合并。如果没有一个邻居粒子满足条件，则粒子i的合并就被取消并推迟到下一时间步中。

过程4：SPH粒子集P执行动力学方程计算，获得粒子n+1时刻的位置和速度；

第一步：建立基于拉格朗日法的流体动力学模型。为了计算n+1时刻粒子的动态流体速度，需要首先求解流体运动的NS方程

其中，ρ为密度，u为流体速度，p为压强，t为时间，f为外力，μ为运动粘性系数，μ▽²u表示的是粘性力，▽·为散度，▽为梯度，▽²＝▽·▽。在SPH模型中，由于粒子个数确定且质量都是常数，满足质量守恒，因此质量守恒方程(3)可以被忽略。其次，公式(4)中左侧表达式可以使用du/dt来表示，即拉格朗日法中对流项u·(▽u)可以方便的处理。因此，方程(4)变为

令f＝-▽p+ρg+μ▽²u，对于粒子i来说，根据牛顿第二定律

其中，a_i表示粒子的加速度，f_i和ρ_i分别是粒子i所在位置的力及密度。

第二步：利用SPH方法估计密度ρ_i。SPH是粒子系统的一种插值方法，利用该方法，可从粒子携带的量中获取该变量在流体场中任何位置中的值。在SPH方法中，变量X在位置r处的值可根据其邻居粒子的加权和得到，即

其中，j表示粒子，m_j表示j的质量，x_j为j的位置，ρ_j为j的密度，X_j为粒子j所携带的量的值，W为各向同性的支撑半径为h的光滑的核函数

那么，粒子所在位置x_i ⁿ处的密度ρ_i为

第三步：利用SPH方法估计外力f_i。这里，外力f＝-▽p+ρg+μ▽²u，故需要估计压力▽p和粘性力μ▽²u。其中压力▽p表示为

压强使用状态方程获取，即p＝k(ρ-ρ₀)，ρ₀为参考密度。为了生成对称的压力，使用(p_i+p_j)/2代替p_i。粘性力表示为

这里用u_j-u_i代替u_j其原因是粘性力与相对速度相关。

第四步：更新粒子的速度和位置。粒子i的外力f_i和密度ρ_i均可被计算出来，那么粒子i在n+1时间步的新速度为

其中，Δt为时间步长。最后，根据dx_i/dt＝u_i，粒子i在n+1时间步的新位置为

过程5：从多精度SPH粒子集P中提取表面粒子集；

第一步：构建场函数，使该函数在粒子中心的取值为1，而在其他位置取0。为了保证光滑性，使用插值方法计算该场函数C_S

第二步：计算每个粒子位置处的法向量。流体表面x处的法向量n可表示为

第三步：根据法向大小确定表面粒子集P_s。在流体内部，n的大小接近与零，而在流体表面位置，由于某个方向上存在粒子空隙，故n的模大于零。故给定阈值θ，如果粒子i满足||n(x_i)||≥θ，则将i并入到表面粒子集P_s中，即P_s＝P_s∪i。

过程6：对多精度表面粒子集P_s进行预处理；

第一步：计算P_s中每个粒子的法向量。对位于P_S中的粒子i，其法向首先通过公式(13)计算，然后将其归一化，得单位向量n_i＝n(x_i)/||n(x_i)||作为该粒子的法向量。

第二步：计算P_s中每个粒子的置信值。采样点的置信值c_i表示该采样点位于流体表面的概率。显然，距离流体表面越近的粒子，其置信值越大，反之则越小。为了计算P_S中的粒子i到表面的距离，需要获取距离i最近的表面粒子。使用法向外推的方法，将表面粒子的编号I外推到其他粒子中，即

其中，n为表面粒子法向量，该方程可以使用快速步进法求解。这样，每个P_S中的粒子都会获取距离其表面最近的表面粒子I，故采样点i的置信值可表示为c_i＝1-α||x_i-x_I||，其中，α＞0为控制参数。

第三步：计算P_c中每个粒子的占用空间大小s_i。s_i表示粒子对周围区域的影响的半径大小，为此本发明定义为s_i＝cr_i，其中c≥1为常数。

第四步：将P_c中的粒子存储在八叉树中。由于表面粒子集P_S表示计算的流体，因此包围盒B设置为整个流体计算域。将包围盒B分解为规则的体网格V，使用八叉树数据结构存储该体网格，其中层最粗糙，而最精细。多精度流体表面粒子集P_S中的粒子即采样点将根据位置被分布在不同的八叉树结点中。

过程7：取八叉树一个层上的节点，计算包围流体表面的壳；

第一步：生成壳V_crust。给定八叉树的一层该层上V_crust的生成通过两个过程完成。首先，生成表示流体表面的V_crust核。为此，将包含表面采样点的体素v的标记v_tag＝1，并将其置入到V_crust中。一般来说，如果表面采样点较多，可以取较精细的层次，以尽快获得高精度的表面；反之，则取较粗糙的层次，否则很难生成封闭的V_crust。其次，根据V_crust核中的体素单元计算位于流体表面的窄带区域上的体素单元，扩充V_crust核。具体的，执行一个标记及采样点法向外推的操作，即

其中，n为采样点的法向。根据公式(17)将当前V_crust体素往n方向膨胀τ次，膨胀的次数由用户定义，同时，将膨胀后的体素v的标记v_tag设置为距离v最近采样点的标记值。这样，V_crust即可表达流体表面及其周围的窄带。

第二步：壳分割。对于体素网格集合V_crust中的每个体素v，需要确定其是位于表面的内侧还是外侧，从而实现对体素的分割。为此，需要比较体素表面的法向与采样点的法向。如果是流体表面的外侧，两个法向夹角为锐角，而在流体表面的内侧，两个法向夹角为钝角。对于体素v(v∈V_crust)来说，其表面上的法向根据体素面的朝向容易计算，表示为对位于V_crust核上包含采样点的体素v的表面法向定义为v内采样点法向量的平均值来表示，即

对V_crust上其他体素表面法向定义为距离v最近的包含采样点的体素的法向量。在法向外推过程中，同时将体素法向量记录即可。那么，如果则体素v为外侧单元。否则体素v为内侧单元。这里ε∈(0,1)表示阈值。

过程8：根据壳创建带权图；

第一步：初始化带权图G(V,E,W)。G的顶点为体素v(v∈V_crust)的每个面f_v，即V＝{f_v}。如果在V_crust中两个面和共享一条边，则在对应的顶点之间增加一条边即E＝{e_ij}。由于V_crust中每个面有四条边，因此生成的图都是规则的4连接的图。

第二步：对G(V,E,W)增加源结点S和汇结点T，并在G中增加连接其他顶点到源结点S以及其他顶点汇结点T的边。具体的，对G中任意顶点f_v，如果该顶点所在体素v是V_crust的内侧且位于边界上，则增加一条连接该顶点到汇结点的边e(f_v,T)。如果该顶点所在体素v是V_crust的外侧且位于边界上，则增加一条连接该顶点到源结点的边e(f_v,S)。

第三步：初始化权重。建立完图G之后，要对G中的边设置权重。对于连接两个非终点的边其权重设置为w。连接终点的边e(f_v,S)和e(f_v,T)的权重为μ，本发明中使用w＝1，μ＝0.3。

过程9：根据粒子的置信值对G(V,E,W)中的权重进行变更；

第一步：计算全局置信值场Γ(x)。全局置信值场为一个从R³到R的映射。该映射为V_crust中的每个体素v计算一个置信值。为此，首先定义每个采样点对周围置信值的影响。对采样点i，其占用空间的大小为s_i，那么采样点对周围空间置信值的贡献为高斯函数

通过叠加所有采样点的空间置信值，可以得到全局置信值场。

为了提高计算效率，在具体实现中首先设定距离采样点大于其三倍s_i的值的体素，其置信度值γ_v＝0。对八叉树的每一层，计算V_crust体素中心的全局置信值Γ(x)，并将其存储在相应的八叉树节点中。当计算完一个八叉树节点的置信值之后，利用深度遍历的方式将该节点中的值加到其孩子结点中。对于采样点i来说，如果其占用空间小于当前八叉树节点中体素的大小，则暂时增大其s_i至体素大小以计算当前的全局置信场。同时，标记该采样点，在下一层八叉树节点中进行高精度的计算。

第二步：对全局置信值场Γ(x)进行自适应优化，得Γ_A(x)。由于置信值的大小要进行归一化，而最大值Γ_max可能很大，因此置信场再归一化后将导致全局置信场的局部变化很小，影响图切割算法的精度。为此，对Γ(x)进行局部自适应优化得

其中，B_edge是包围盒的边长。W(d)为权重函数

其中，D为体素中的滤波半径，使用D＝10，m＝7。

第三步：基于优化后的置信值场Γ_A(x)对带权图G(V,E,W)的权重变更。在V_crust中的每个体素v共有六个面，而对于每个面，都有四条边与其他的面共享。因此，体素v中包含一个具有12条边的八面体图，每条边的权重为w。最佳流体表面应该是在最小割(面积最小)的基础上具有最大化的全局置信度。为此，需要在具有较高置信值的体素内设置较小的边权重，此时，修改图G(V,E,W)中2个顶点均位于体素v内的边的权重为

其中，Γ_Amax表示自适应全局置信场中的最大值，a为控制参数，调整后的权重值为w∈[a,1+a]。

过程10：计算当前层上的流体表面；

第一步：计算流体三角形表面网格。图G(V,E,W)的最小割将产生一组切割边，这组切割边是在源结点和汇结点之间形成分割，同时也是获得的具有最高置信值的对流体表面的逼近。流体表面是在包含至少一条割边的体网格内产生。具体的，相邻的2个包含割边的体网格中心点相连，即可生成三角形的一条边，因此2×2×2的体网格内即可生成一个三角形。这样，连接相邻的、包含割边的体素中心点即可形成流体表面的三角形网格。

第二步：表面网格光滑化。由于三角形的顶点只位于体素的中心，因此在获取到流体表面后，可执行一个拉普拉斯光滑化的过程，得到层上的光滑的流体表面M_i。

过程11：构建多精度流体表面；

第一步：迭代生成不同精度的流体表面网格。上述过程7-过程10是从一个八叉树层上生成对应精度的流体表面。在多精度流体表面的构建中，将不断重复上述过程，直到所有的高精度细节能够被重建出来。在过程9计算全局置信场时，对那些占用空间小于当前八叉树的节点内体素网格大小的采样点做标记。所有在当前层标记的采样点将在下一层中进行更高精度的计算，直到下一层中没有任何被标记的采样点为止。

第二步：将不同精度的流体表面网格融合在一起，形成多精度表面。具体过程是，先将初始层(第i层，i＝0)产生的表面网格M_i加入到M中，然后将第i+1层产生的表面网格M_i+1加入到M中，并舍弃M_i+1所在位置的低精度三角形。在M_i+1与M_i的交界处，令在同一个体素内的低精度顶点代替高精度顶点完成表面的融合。上述过程迭代计算直到八叉树所有层的网格融合到M中。

过程12：对多精度表面进行真实感渲染得到当前帧的流体动画效果图。

过程13：重复过程4-过程12，直到流体计算结束，即可得到体现多尺度细节的流体动画。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (3)

1.基于多精度融合的流体动画生成方法，其特征是，包括如下步骤：

步骤(1)：初始化NS方程的初始条件和边界条件；所述初始条件为流体源S，流体源S由SPH粒子集P表示；所述边界条件B包括固体与液体边界、空气和液体边界；

步骤(2)：根据初始化的流体源S及其对应的粒子集P，基于NS方程，建立基于SPH的多精度流体计算方法，获得多精度粒子集；

所述步骤(2)的步骤如下：

步骤(2-1)：定义并计算多精度粒子半径的自适应分布函数；

步骤(2-2)：根据粒子半径的自适应分布函数对粒子进行分裂和合并，获得多精度粒子集；

步骤(2-3)：对多精度粒子集，计算流体动力学方程，更新其速度和位置；

步骤(3)：对步骤(2)的多精度粒子集进行预处理，将处理后的粒子集使用八叉树存储；

步骤(4)：从八叉树顶层开始不断重复迭代执行步骤(5)直到八叉树底层；

步骤(5)：对八叉树的每一层，计算并提取该层对应精度的流体表面网格；

步骤(6)：将不同精度的流体表面网格融合成一个流体表面，使用真实感渲染获得每帧的流体动画；

步骤(7)：如果当前帧数小于用户设定的动画帧数，则转步骤(2)继续下一帧流体动画的计算，最终得到体现多尺度细节的流体动画。

2.如权利要求1所述的基于多精度融合的流体动画生成方法，其特征是，

所述步骤(1)初始化NS方程的初始条件就是设定粒子集P中每个粒子的半径大小、位置和速度；

所述步骤(1)初始化边界条件B就是根据输入的流体所在的场景模型和流体源S，设定固体和液体边界；在剩余空间和流体源S之间设定空气和液体边界。

3.如权利要求1所述的基于多精度融合的流体动画生成方法，其特征是，所述步骤(5)的步骤如下：

步骤(5-1)：获取八叉树一层，计算包围该层所有粒子的壳；

步骤(5-2)：根据所得的壳，创建带权图；

步骤(5-3)：计算全局和局部置信值场，对带权图的权值进行变更；

步骤(5-4)：计算带权图的最小割，获得八叉树当前层对应的流体表面网格。