CN103617593A - 三维流体物理动画引擎的实现方法及装置 - Google Patents

Abstract

本发明公开了一种三维流体物理动画引擎的实现方法及装置，其中，该方法包括：利用并行空间粒子搜索算法搜索流体中每一粒子的邻居粒子，并计算所述每一粒子的状态；利用复杂边界查找算法计算流体所处容器的边界，根据所述每一粒子的状态并基于碰撞机制计算与调节所述每一粒子与边界的碰撞强度；根据碰撞后的粒子分布情况并基于流体表面重构算法对流体表面重构，或者基于屏幕空间的表面重构算法对流体表面重构，获得流体表面。通过采用本发明公开的方法及装置提高了运行速度及工作效率。

Description

三维流体物理动画引擎的实现方法及装置

技术领域

本发明涉及基于物理学的计算机动画领域，尤其涉及一种三维流体物理动画引擎的实现方法及装置。 

背景技术

SPH（Smoothed Particle Hydrodynamics，光滑粒子动力学）方法是一种典型的拉格朗日方法，在模拟流体时它是将连续的流体用相互作用的质点组描述，各个质点上承载各种物理量，包括质量、速度、密度等，通过求解质点组的动力学偏微分方程组和跟踪每个质点的运动，求得整个系统的力学行为。SPH方法使用核函数对密度、压力、粘力项进行离散化，得到Navier-Stokes（纳维叶－斯托克斯）方程组的离散化计算形式，从而在每步迭代过程近似解出各个物理量，模拟流体运动。SPH方法是由Lucy、Monaghan及Gingold等人于1977年提出用来解决天体物理学中的行星运动问题，后来被引入到计算流体力学领域。 

SPH方法模拟流体具有自适应性特性。SPH方法的这种自适应性是通过在每个时间步长计算其周围的其他粒子分布近似地获得的。由于SPH方法的这种特性，SPH方法不受粒子分布的影响。因此，SPH方法能够很好的处理具有剧烈形变问题。这是SPH方法最引起人们感兴趣的特性。另外SPH方法还克服了高维拉氏方法中网格缠绕的问题和随着维数、网格数增加其复杂度急剧上升的问题。 

发明人在进行发明创造的过程中发现，现有技术中，SPH方法模拟流体的计算量较大，运行速度较慢，工作效率较低，因此，如何提高三维流体物理动画的运行速度及工作效率成为现今的研究重点。 

发明内容

本发明的目的是提供一种三维流体物理动画引擎的实现方法及装置，提高了运行速度及工作效率。 

本发明的目的是通过以下技术方案实现的： 

三维流体物理动画引擎的实现方法，该方法包括： 

利用并行空间粒子搜索算法搜索流体中每一粒子的邻居粒子，并计算所述每一粒子的状态； 

利用复杂边界查找算法计算流体所处容器的边界，根据所述每一粒子的状态并基于碰撞机制计算与调节所述每一粒子与边界的碰撞强度； 

根据碰撞后的粒子分布情况并基于流体表面重构算法对流体表面重构，或者基于屏幕空间的表面重构算法对流体表面重构，获得流体表面。 

设计多种流体光学效应的计算方法，力求在增强流体的可视化效果同时消耗较小的系统资源。 

三维流体物理动画引擎的实现装置，该装置包括： 

并行处理单元，用于计算流体中每一粒子的状态，并利用并行空间粒子搜索算法搜索所述每一粒子的邻居粒子； 

碰撞强度控制单元，用于利用复杂边界查找算法计算流体所处容器的边界，根据所述每一粒子的状态并基于碰撞机制计算与调节所述每一粒子与边界的碰撞强度； 

流体表面获取单元，用于根据碰撞后的粒子分布情况并基于流体表面重构算法对流体表面重构，或者基于屏幕空间的表面重构算法对流体表面重构，获得流体表面。 

由上述本发明提供的技术方案可以看出，通过并行化处理粒子状态更新，并对粒子的邻居粒子进行搜索，使得系统能快速、高效的进行大规模粒子的查找；通过利用复杂边界查找算法，使得系统能在大规模复杂边界的条件下流体的仿真模拟仍然能够保持较高的帧率；另外，通过针对应用场景的不同，选择不同的表面重构算法，能够准确计算流体密度场和重构流体表面，并能够快速实时重构表面。 

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。 

图1为本发明实施例一提供的一种三维流体物理动画引擎的实现方法的流程图； 

图2为本发明实施例一提供的流体厚度与观察者关系的示意图； 

图3为本发明实施例一提供的焦散效应形成原理的示意图； 

图4为本发明实施例一提供的一种焦散效应的示意图； 

图5为本发明实施例二提供的一种三维流体物理动画引擎的实现装置的示意图。 

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。 

实施例一 

图1为本发明实施例一提供的一种三维流体物理动画引擎的实现方法的流程图。如图1所示，该方法主要包括如下步骤： 

步骤11、利用并行空间粒子搜索算法搜索流体中每一粒子的邻居粒子，并计算所述每一粒子的状态。 

由于拉格朗日方法需要利用支持域内的邻居粒子，因此，本发明实施例通过对邻居粒子物理量的加权平均来计算当前粒子所承载的相应物理量。 

本发明实施例中可以使用CUDA（统一计算设备架构）并行化处理粒子状态更新计算，具体的：对流体中每一粒子进行三维空间上的哈希操作，实现当前粒子位置到存储地址的映射，并查到找到预设的核函数半径内所述当前粒子的所有邻居粒子；将所述核函邻居粒子的状数半径内所述当前粒子的所有态带入至Navier-Stokes（纳维叶－斯托克斯）方程组离散对称形式中，并使用多项式函数作为核函数，计算所述当前粒子的密度、压力及粘力等物理状态，获得所述当前粒子的加速度；使用蛙跳式更新方法计算所述当前粒子的速度和位置状态。 

示例性的，本发明实施例中利用一6次多项式核函数对密度求解，其公式为： 

求解压力项的核函数为： 

粘度项求解使用核函数为： 

上述公式中h表示核函数的支持域半径，r表示支持域内粒子到当前粒子的距离。 

通常，邻居粒子搜索的计算量较大，本实施例设计了基于CUDA平台并行空间粒子搜索算法，并将所有粒子的邻居粒子连续存储，使空间粒子搜索算法所需的存储空间得以减少，并行化加速后算法使得系统能快速、高效的进行大规模领域粒子查找。 

该搜索算法主要包括如下几个步骤： 

（1）哈希映射处理：即对粒子点云数据所在的空间进行三维网格划分，然后再针对每个单元格得到所有空间上位于该单元格中的粒子，形成一个粒子链表。因此，需预先要对流体中的粒子i的位置P_i使用空间哈希操作，得到三维坐标A_i(A_i·x，A_i·y，A_i·z)，其计算公式为： 

$A_i\·x=\left[\frac{P_i\·x}{2h}\right];$

$A_i\·y=\left[\frac{P_i\·y}{2h}\right];$

$A_i\·z=\left[\frac{P_i\·z}{2h}\right];$

其中，2h为该流体所处容器的分段长度； 

由于网格数组在内存中以一维形式组织，因此要根据上述的三维地址坐标A_i(A_i·x，A_i·y，A_i·z)来计算一维存储偏移a'_i，其计算公式为： 

a'_i＝A_i·x+A_i·y·Dx+A_i·z·Dx·Dy； 

其中，Dx，Dy分别表示x，y方向上的最大单元格的数量，即所有粒子的包围盒以2h为单位的最大分段数。 

其次，还要针对每个单元格维护一个粒子链表，该链表的建立方法是：并行遍历所有粒子，利用头插法将粒子插入所在的单元格粒子链表中。 

（2）邻居粒子数目查询：计算所述粒子i所处单元格的三维坐标 

根据所述粒子i所处单元格中所维护的粒子链表查找该单元格内的粒子，并判断该单元格内的粒子与所述粒子i之间的距离是否小于h，若是，则为所述粒子i的邻居粒子，则将当前的粒子的邻居粒子数加一，遍历所有邻居粒子后得到当前粒子的邻居粒子数。 

示例性的，所述预设的三维坐标可以包括： $(A_i^f\·x,A_i^f\·y,A_i^f\·z),$ $(A_i^f\·x+1,A_i^f\·y,A_i^f\·z),(A_i^f\·x,A_i^f\·y+1,A_i^f\·z),(A_i^f\·x+1,A_i^f\·y+1,A_i^f\·z),$ $(A_i^f\·x,A_i^f\·y,A_i^f\·z+1),(A_i^f\·x+1,A_i^f\·y,A_i^f\·z+1),(A_i^f\·x,A_i^f\·y+1,A_i^f\·z+1),$ $(A_i^f\·x+1,A_i^f\·y+1,A_i^f\·z+1).$

（3）压缩存储空间：为了提高每次迭代过程中邻居粒子查找效率，用一个连续的数组M存放所有粒子的邻居粒子，数组M的构建方法如下：并行遍历每个粒子，先取得每个粒子的邻居粒子个数（在前述步骤（2）中已经得到），存入数组A；对数组A进行并行的exclusive scan求得其前缀和，存入数组B；最后，并行遍历所有粒子，将该粒子的每个邻居粒子按B数组中该粒子位置处的元素作为访问M数组的下标，将邻居粒子的id存入M数组中。 

（4）邻居粒子搜索：B数组中存放的是每个粒子所对应的邻居粒子在M数组中的起始位置，而每个粒子的邻居粒子数目在步骤（2）中已经获得，所以在后续计算密度、压力、粘力等步骤中查找邻居粒子只需访问步骤（3）中所构建的M数组即可，这种算法的思想是以空间换取时间。 

步骤12、利用复杂边界查找算法计算流体所处容器的边界，根据所述每一粒子的状态并基于碰撞机制计算与调节所述每一粒子与边界的碰撞强度。 

流体中的粒子在接触到容器界面时与边界相互作用，一般流体模拟中使用的容器是一些简单形状的容器，边界条件也相对简单，如长方体可以使用六个矩形边界表示，进行边界作用时可以遍历查找所有的边界条件。但一些复杂的容器或者在多孔介质中，其边界数量规模较大，分布复杂，对于这些边界作用由于遍历花费时间较长，不能进行直接遍历查找。本发明实施例设计了针对复杂边界情况下（如多孔介质边界）邻居边界的查找算法，该算法使得系统能在大规模复杂边界条件下流体的仿真仍然能够获得一定的帧率。 

该复杂边界查找算法与邻居粒子查找过程类似，在使用该复杂边界查找算法时可以进行CUDA加速处理，以加快处理速度。同时，在使用该复杂边界查找算法之前还需进行预处理，具体的：将所述容器的边界模型化为点和法向表示的理想平面；对所述理想平面进行离散采样，再进行哈希映射，以确保所述理想平面经过单元格区域。 

由于在进行复杂边界查找算法时，静态的容器边界是不动的，故相对邻居粒子查找过程，省略了哈希映射过程，而将过程在预处理中进行处理。 

另一方面，本发明实施例还根据所述每一粒子的状态并基于碰撞机制计算与调节所述每一粒子与边界的碰撞强度。 

所述碰撞机制包括：u_i＝u_i-(1+c_R)(u_i·n)n，其中，u_i为粒子i的速度，n为平面法向量，c_R为动量损失系数。 

在粒子碰撞时通过调节动量损失系数来控制粒子与边界碰撞强度，从而保证粒子与边界交互时符合物理运动同时快速处理粒子与边界交互。 

步骤13、根据碰撞后的粒子分布情况并基于流体空间密度场对流体表面重构，或者基于屏幕空间的表面重算算法对流体表面重构，获得流体表面。 

本发明实施例针对应用场景的不同采用不同的表面重构算法： 

1、根据碰撞后的粒子分布情况并基于流体空间密度场对流体表面重构，可准确计算流体密度场和重构流体表面。具体的：根据所述流体中的粒子搜索局部SPH粒子，采用矩阵理论中有关SVD（奇异值）分解的研究成果，并根据所述局部SPH粒子的位置计算所述流体中粒子分布的各向异性和大小；将所述流体中粒子分布的各向异性和大小作为核函数系数计算每一粒子局部密度场分布，以减少突兀的粒子对流体密度场的影响，再叠加所有粒子密度场，获得全局密度场分布，最后使用等值面抽取算法，获得流体表面。 

所述等值面抽取算法可以采用Marching Cube(移动立方体)算法，该算法通过一个移动的匹配立方体来抽取等值面，场数据一旦确定，算法的独立性较高。但考虑到大量无效匹配立方体存在拖慢整个GPU的处理速度，所以先进行匹配立方体有效性扫描。具体过程如下： 

1）标记，扫描全局密度场中的所用体素，得到体素的标记数组T[k]，0≤k≤M-1，M为密度场体素个数。如果体素k有等值面经过则T[k]为1，否则为0。 

2）建立映射关系，对所述数组T[k]的前缀和进行计算获得数组T’[k]，则数组T[k]中被标记的总个数为B=T’[M-1]+T[M-1]，由于空体素较多，则B远小于M；建立长度为B的数组V[j]，其中，0≤j≤N-1，且V[T’[i]]=k。 

通过以上两个步骤的操作，在进行Marching Cubes算法重构表面时可以直接对数组V[j]进行，通过数组V[j]的值就可以直接确定对应的体素非空。然后由CUDA每个核心负责一个有效的匹配立方体，最后进行等值面的三角网格生成。 

2、基于屏幕空间的表面重构算法对流体表面重构，可提高显示速度。具体的：将所 述流体中的粒子替换为球体，利用球体之间遮挡关系来近似流体表面，降低重构表面的空间和时间复杂度；利用计算机图形学中的仿射变换原理，在图像空间下计算表面像素的深度与法向等；再对计算得到的所述表面像素的深度与法向进行低通滤波，得到较为光滑的流体表面。 

步骤14、光学效应计算。 

环境对流体的影响主要集中在对流体的光学效应，包括光照、折射与反射、光线穿透的厚度、焦散和阴影等。本发明实施例的仿真模拟方法是一种基于粒子的流体模拟方法，因此，在对流体进行渲染时，不仅要考虑到流体的光学特征，还要考虑其粒子特征。本发明实施例假定流体的光学效应互不影响，通过计算各种光学效应的简化效果，最后叠加在一起，以达到渲染流体的目的。 

为了便于说明，下面针对光学效应的计算方法做详细的介绍： 

1）光照计算：本发明实施例可以采用OpenGL提供的光照模型进行场景的光照计算，通常情况下其光照计算均能达到视觉要求。 

2）折射与反射系数计算：本发明实施例可以采用环境贴图来模拟折射和反射效应，针对介质根据入射角度的不同对光线的折射和反射的也不一样（菲涅尔现象），可以采用一种简化的形式来计算不同入射角下的折射和反射系数，其计算公式为： 

r'＝fb+fs·(1.0+I·G)^fp； 

r＝min(1.0，r')； 

s＝1.0-r； 

其中，r'为中间变量（通常钳位在0到1之间），I为入射光线单位矢量，G为顶点处的法线，fb、fs与fp表示与流体介质相关的多项式系数，r为反射系数，s为折射系数。 

3）光线穿透的厚度：流体中光线能量一般随着传播距离的增加不断衰弱，尤其对于具有一定浓度的流体这种现象会更加明显，因此，通过流体的厚度计算和渲染来体现光线穿透的厚度。 

如图2所示，表示了流体厚度与观察者的关系，通过计算流体的厚度，可使观察者看到A点的流体厚度比C点的更深。 

根据SPH理论，每个粒子对周围空间密度影响是对称的，为尽量保持SPH流体密度特征，本发明实施例将流体厚度计算与核函数及密度联系起来，其计算方法如下： 

$T_{\mathrm{xy}}=c_0\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_{i}W(M_{\mathrm{mvp}}{P}_i-(x,y),h);$

I_i＝c₀ρ_i； 

其中，T_xy表示屏幕空间下x,y位置的流体厚度，M_mvp表示模视投影变换矩阵，ρ_i表示粒子i的密度，c₀为常系数。 

基于以上假设每个粒子可以看成一个亮度为I_i的Billboard，并带有核函数特征的纹理贴图。在OpenGL中，启用颜色混合模式绘制所有Billboard，绘制的结果会缓冲到帧缓冲区，最终流体的厚度信息被变换成帧图像的亮度信息。GPU从缓冲区取出结果进行简单的图像光滑处理，将处理结果交给着色模块进行合成。Billboard是一个始终面向观察者的二维矩形，绘制的复杂性较小，因此绘制速度较快。 

4）焦散效应：当光线穿过一个不平整流体表面时，折射光线并没有平行穿过，而出现漫折射，在接收面上出现不同的光线聚散，形成焦散效应，如图3所示，在接收面的有些地方（如A点）光线聚集，形成较亮区域，有些地方分散（如B点），形成相对较暗区域。 

本发明实施例在计算焦散效应时，根据表面重构算法的不同，选择不同的方法进行计算： 

①基于Marching Cubes流体表面重构算法的焦散效应计算。流体经过MarchingCube算法表面重构，形成一系列的面片，输出的数据组织结构是一系列的三角形顶点和顶点处的法向量。如图4所示，为发生焦散效应时的示意图；P表示入射光线I与流体表面的交点，N为该平面的法向量，其折射率为eta，Ref表示入射光线I射入所述流体后发生折射后的光线，折射光线Ref与接收平面（由位置Pg，法向量Ng决定的平面）相交于P'。 

其主要涉及三个特征量的计算：焦散效应光子的位置P'、焦散效应的影响范围S_b及焦散效应强度C_b。 

首先，计算焦散效应光子的位置P'，计算公式为： 

t＝1-eta²·(1.0-(N·I)²)； 

$\mathrm{Ref}=\mathrm{eta}\·I-(\mathrm{eta}\·(N\·I)+\sqrt{t})\·N;$

$\mathrm{Re}{f}^{\′}=\frac{\mathrm{Ref}}{\left|\right|\mathrm{Ref}\left|\right|};$

disPtoNg＝(P-Pg)·Ng 

$P^{\′}=P+\frac{\mathrm{disPtoNg}}{-\mathrm{Ng}\·\mathrm{Ref}}\·\mathrm{Re}{f}^{\′};$

其中，t为中间变量，ref’是折射光线的方向，disPtoNg是P到由点Pg和法向Ng构成平面的距离；另外，由于光线不能与面片背面作用，所以上述计算前提是N·I＞0。 

然后，计算焦散效应影响范围S_b，其与面片面积相关，表示该面片焦散效应在位置P'的影像范围，计算公式为： 

S_b＝c₁·S_m； 

其中，c₁为比例系数，S_m为Marching Cubes三角形面片面积或像素强度。 

再计算焦散效应强度C_b，其与入射角相关，计算公式为： 

$C_b=c_2\·{\left(\right|N\·I\left|\right)}^{c_3};$

其中，c₂为比例系数，c₃为指数系数。 

同样，可以使用Billboard技术，将大小为S_b亮度为C_b的Billboard在P'处平行于接收面绘制，使用OpenGL混合模式合成。将合成结果给GPU做图像平滑处理，然后由着色程序对其进行合成处理。 

②基于屏幕空间的表面重构算法的焦散效应计算。该焦散效应与光照密切相关，光照能照射到的表面与观察者能看到的表面具有相同的原理，不同的是投影的点和投影的角度不同，因此在光照空间下可以得到表面，其原理与屏幕空间下的表面重构方法大致相同。具体的：在光照空间下绘制粒子球，获得光照表面；恢复所述光照表面在世界坐标系下的位置；在视点空间下计算焦散效应光子的位置；其中，坐标系变换过程的公式为：

v'表示投影后的归一化视点坐标，v表示执行坐标系变换前的世界坐标，M_l与P_l分别表示模式矩阵和投影矩阵。 

本发明实施例通过并行化处理粒子状态更新，并对粒子的邻居粒子进行搜索，使得系统能快速、高效的进行大规模粒子的查找；通过利用复杂边界查找算法，使得系统能在大规模复杂边界的条件下流体的仿真模拟仍然能够保持较高的帧率；另外，通过针对应用场景的不同，选择不同的表面重构算法，能够准确计算流体密度场和重构流体表面，并能够快速实时重构表面；最后设计了多种光学效应计算方法，在假设各种效应互不影响的前提下可独立计算，最后再合成以丰富流体渲染效果。 

实施例二 

图5为本发明实施例一提供的一种三维流体物理动画引擎的实现装置的示意图。如图5所示，该装置主要包括： 

并行处理单元51，用于利用并行空间粒子搜索算法搜索流体中每一粒子的邻居粒子，并计算所述每一粒子的状态； 

碰撞强度控制单元52，用于利用复杂边界查找算法计算流体所处容器的边界，根据所述每一粒子的状态并基于碰撞机制计算与调节所述每一粒子与边界的碰撞强度； 

流体表面获取单元53，用于根据碰撞后的粒子分布情况并基于流体表面重构算法对流体表面重构，或者基于屏幕空间的表面重构算法对流体表面重构，获得流体表面。 

进一步的，所述并行处理单元51包括： 

粒子状态计算单元511，用于对流体中每一粒子进行三维空间上的哈希操作，实现当前粒子位置到存储地址的映射，并查到找到预设的核函数半径内所述当前粒子的所有邻居粒子；将所述核函邻居粒子的状数半径内所述当前粒子的所有态带入至Navier-Stokes方程组离散对称形式中，并使用多项式函数作为核函数，计算所述当前粒子的密度、压力及粘力，获得所述当前粒子的加速度；使用蛙跳式更新方法计算所述当前粒子的速度和位置状态； 

搜索单元512，用于搜索所述每一粒子的邻居粒子，且该单元包括： 

哈希映射处理单元5121，用于对流体中的粒子i的位置P_i使用空间哈希操作，得到三维坐标A_i(A_i·x，A_i·y，A_i·z)，其计算公式为： 

$A_i\·x=\left[\frac{P_i\·x}{2h}\right];$

$A_i\·y=\left[\frac{P_i\·y}{2h}\right];$

$A_i\·z=\left[\frac{P_i\·z}{2h}\right];$

其中，2h为该流体所处容器的分段长度； 

根据所述三维地址坐标A_i（x，y，z），计算一维存储偏移a'_i，其计算公式为： 

a'_i＝A_i·x+A_i·y·Dx+A_i·z·Dx·Dy； 

其中，Dx，Dy分别表示x，y方向上的最大单元格的数量； 

邻居粒子数目查询单元5122，用于计算所述粒子i所处单元格的三维坐标 根据所述粒子i所处单元格中所维护的粒子链表查找该单元格内的粒子，并判断该单元格内的粒子与所述粒子i之间的距离是否小于h，若是，则为所述粒子i的邻居粒子； 

存储空间压缩单元5123，用于构建一连续的数组M存放所有粒子的邻居粒子，其构 建方法为：并行遍历每个粒子，获得每个粒子的邻居粒子个数，存入数组A；计算所述数组A的前缀和，存入数组B；并行遍历所有粒子，将每一粒子的每个邻居粒子按数组B中其位置处的元素作为访问M数组的下标，将邻居粒子的id存入M数组中； 

邻居粒子搜索单元5124，用于通过所述数组M，完成每一粒子的邻居粒子的搜索。 

进一步的，所述碰撞强度控制单元52包括： 

预处理单元521，用于在计算所述流体所处容器的边界之前，进行预处理，具体的包括：将所述容器的边界模型化为点和法向表示的理想平面；对所述理想平面进行离散采样，再进行哈希映射，以确保所述理想平面经过单元格区域； 

所述碰撞机制包括：u_i＝u_i-(1+c_R)(u_i·n)n，其中，ui为粒子i的速度，n为平面法向量，c_R为动量损失系数；通过调节动量损失系数来控制粒子与边界碰撞强度。 

进一步的，所述流体表面获取单元53包括： 

第一流体表面重构单元531，用于根据碰撞后的粒子分布情况并基于流体表面重构算法对流体表面重构包括：根据所述流体中的粒子搜索局部SPH粒子，并根据所述局部SPH粒子的位置计算所述流体中粒子分布的各向异性和大小；将所述流体中粒子分布的各向异性和大小作为核函数系数计算每一粒子局部密度场分布，并叠加所有粒子密度场，获得全局密度场分布，再使用移动立方体Marching Cube算法抽取等值面，获得流体表面； 

且所述第一流体表面重构单元531还包括：匹配立方体有效性扫描单元5311，用于在使用Marching Cube算法抽取等值面之前进行匹配立方体有效性扫描，且该单元包括：标记单元53111，用于扫描全局密度场中的所用体素，得到体素的标记数组T[k]，0≤k≤M-1，M为密度场体素个数；映射关系建立单元53112，用于对所述数组T[k]的前缀和进行计算获得数组T’[k]，则数组T[k]中被标记的总个数为B=T’[M-1]+T[M-1]；建立长度为B的数组V[j]，其中，0≤j≤N-1，且V[T’[i]]=k； 

第二流体表面重构单元532，用于基于屏幕空间的表面重构算法对流体表面重构包括：将所述流体中的粒子替换为球体；利用计算机图形学中的仿射变换原理，在图像空间下计算表面像素的深度与法向；再对计算得到的所述表面像素的深度与法向进行低通滤波。 

进一步的，该装置还包括：光学效应计算单元54；所述光学效应包括：包括光照、折射与反射、光线穿透的厚度、焦散和阴影； 

所述光学效应计算单元；还包括： 

折射与反射系数计算单元541，其计算公式为： 

r'＝fb+fs·(1.0+I·G)^fp； 

r＝min(1.0，r')； 

s＝1.0-r； 

其中，I为入射光线单位矢量，G为顶点处的法线，fb、fs与fp表示与流体介质相关的多项式系数，r为反射系数，s为折射系数； 

流体厚度计算单元542，用于计算所述光线穿透的厚度之前计算所述流体的厚度，其计算公式为： 

$T_{\mathrm{xy}}=c_0\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_{i}W(M_{\mathrm{mvp}}{P}_i-(x,y),h);$

I_i＝c₀ρ_i； 

其中，T_xy表示屏幕空间下x，y位置的流体厚度，M_mvp表示模视投影变换矩阵，ρ_i表示粒子i的密度，c₀为常系数； 

焦散效应计算单元543，用于进行焦散效应的计算，包括：基于Marching Cubes流体表面重构算法及基于屏幕空间的表面重构算法的焦散效应计算；且该单元包括： 

第一焦散效应计算单元5431，用于基于Marching Cubes流体表面重构算法计算焦散效应，且该单元包括：位置计算单元54311，用于计算焦散效应光子的位置P'，所述焦散效应光子的位置P'的计算公式为：

其中，P表示入射光线I与流体表面的交点，Ref表示入射光线I射入所述流体后发生折射后的光线，Ng表示所述焦散效应光子的位置P'所处平面的法向量，disPtoNg是P到由点Pg和法向Ng构成平面的距离；disPtoNg＝(P-Pg)·Ng，Pg为所述焦散效应光子的位置P'所处平面上的点， 

影响范围计算单元54312，用于计算焦散效应的影响范围S_b，所述焦散效应的影响范围S_b的计算公式为：S_b＝c₁·S_m；其中，c₁为比例系数，S_m为MarchingCubes三角形面片面积或像素强度；强度计算单元54313，用于计算焦散效应强度C_b，所述焦散效应强度C_b的计算公式为：C_b＝c₂·(|N·I|)^c3；其中，c₂为比例系数，c₃为指数系数，N为入射光线I与流体表面交点P所处平面的法向量，且N·I＞0； 

第二焦散效应计算单元5432，用于基于屏幕空间的表面重构算法计算焦散效应，包括：在光照空间下绘制粒子球，获得光照表面；恢复所述光照表面在世界坐标系下的位置；在视点空间下计算焦散效应光子的位置；其中，坐标系变换过程的公式为： 

v’表示投影后的归一化视点坐标，v表示执行坐标系变换前的世界坐标，M_l与P_l分别表示模式矩阵和投影矩阵。 

需要说明的是，上述装置中包含的各个功能模块所实现的功能的具体实现方式在前面的各个实施例中已经有详细描述，故在这里不再赘述。 

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，仅以上述各功能模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能模块完成，即将装置的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。 

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质（可以是CD-ROM，U盘，移动硬盘等）中，包括若干指令用以使得一台计算机设备（可以是个人计算机，服务器，或者网络设备等）执行本发明各个实施例所述的方法。 

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。 

Claims (10)

1.三维流体物理动画引擎的实现方法，其特征在于，该方法包括：

利用并行空间粒子搜索算法搜索流体中每一粒子的邻居粒子，并计算所述每一粒子的状态；

利用复杂边界查找算法计算流体所处容器的边界，根据所述每一粒子的状态并基于碰撞机制计算与调节所述每一粒子与边界的碰撞强度；

根据碰撞后的粒子分布情况并基于流体表面重构算法对流体表面重构，或者基于屏幕空间的表面重构算法对流体表面重构，获得流体表面。

2.根据权利要求1所述的方法，其特征在于，所述利用并行空间粒子搜索算法搜索流体中每一粒子的邻居粒子，并计算所述每一粒子的状态包括：

对流体中每一粒子进行三维空间上的哈希操作，实现当前粒子位置到存储地址的映射，并查到找到预设的核函数半径内所述当前粒子的所有邻居粒子；将所述核函邻居粒子的状数半径内所述当前粒子的所有态带入至Navier-Stokes方程组离散对称形式中，并使用多项式函数作为核函数，计算所述当前粒子的密度、压力及粘力，获得所述当前粒子的加速度；使用蛙跳式更新方法计算所述当前粒子的速度和位置状态；

其中，搜索所述每一粒子的邻居粒子包括：

哈希映射处理：对流体中的粒子i的位置P_i使用空间哈希操作，得到三维坐标A_i(A_i·x，A_i·y，A_i·z)，其计算公式为：

$A_i\·x=\left[\frac{P_i\·x}{2h}\right];$

$A_i\·y=\left[\frac{P_i\·y}{2h}\right];$

$A_i\·z=\left[\frac{P_i\·z}{2h}\right];$

其中，2h为该流体所处容器的分段长度；

根据所述三维地址坐标Ai（x，y，z），计算一维存储偏移a'_i，其计算公式为：

a'_i＝A_i·x+A_i·y·Dx+A_i·z·Dx·Dy；

其中，Dx，Dy分别表示x，y方向上的最大单元格的数量；

邻居粒子数目查询：计算所述粒子i所处单元格的三维坐标

根据所述粒子i所处单元格中所维护的粒子链表查找该单元格内的粒子，并判断该单元格内的粒子与所述粒子i之间的距离是否小于h，若是，则为所述粒子i的邻居粒子；

压缩存储空间：构建一连续的数组M存放所有粒子的邻居粒子，其构建方法为：并行遍历每个粒子，获得每个粒子的邻居粒子个数，存入数组A；计算所述数组A的前缀和，存入数组B；并行遍历所有粒子，将每一粒子的每个邻居粒子按数组B中其位置处的元素作为访问M数组的下标，将邻居粒子的id存入M数组中；

邻居粒子搜索：通过所述数组M，完成每一粒子的邻居粒子的搜索。

3.根据权利要求1所述的方法，其特征在于，

计算所述流体所处容器的边界之前，进行预处理，具体的包括：将所述容器的边界模型化为点和法向表示的理想平面；对所述理想平面进行离散采样，再进行哈希映射，以确保所述理想平面经过单元格区域；

所述碰撞机制包括：u_i＝u_i-(1+c_R)(u_i·n)n，其中，ui为粒子i的速度，n为平面法向量，c_R为动量损失系数；通过调节动量损失系数来控制粒子与边界碰撞强度。

4.根据权利要求1所述的方法，其特征在于，

所述根据碰撞后的粒子分布情况并基于流体表面重构算法对流体表面重构包括：根据所述流体中的粒子搜索局部SPH粒子，并根据所述局部SPH粒子的位置计算所述流体中粒子分布的各向异性和大小；将所述流体中粒子分布的各向异性和大小作为核函数系数计算每一粒子局部密度场分布，并叠加所有粒子密度场，获得全局密度场分布，再使用移动立方体Marching Cube算法抽取等值面，获得流体表面；

其中，在使用Marching Cube算法抽取等值面之前还包括：匹配立方体有效性扫描的步骤，且包括：标记，扫描全局密度场中的所用体素，得到体素的标记数组T[k]，0≤k≤M-1，M为密度场体素个数；建立映射关系，对所述数组T[k]的前缀和进行计算获得数组T’[k]，则数组T[k]中被标记的总个数为B=T’[M-1]+T[M-1]；建立长度为B的数组V[j]，其中，0≤j≤N-1，且V[T’[i]]=k；

所述基于屏幕空间的表面重构算法对流体表面重构包括：将所述流体中的粒子替换为球体；利用计算机图形学中的仿射变换原理，在图像空间下计算表面像素的深度与法向；再对计算得到的所述表面像素的深度与法向进行低通滤波。

5.根据权利要求1或4所述的方法，其特征在于，在获得所述流体表面之后还包括：光学效应计算；所述光学效应包括：包括光照、折射与反射、光线穿透的厚度、焦散和阴影；

其中，折射与反射系数的计算公式为：

r'＝fb+fs·(1.0+I·G)^fp；

r＝min(1.0，r')；

s＝1.0-r；

其中，I为入射光线单位矢量，G为顶点处的法线，fb、fs与fp表示与流体介质相关的多项式系数，r为反射系数，s为折射系数；

计算所述光线穿透的厚度之前计算所述流体的厚度，其计算公式为：

$T_{\mathrm{xy}}=c_0\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_{i}W(M_{\mathrm{mvp}}{P}_i-(x,y),h);$

I_i＝c₀ρ_i；

其中，T_xy表示屏幕空间下x,y位置的流体厚度，M_mvp表示模视投影变换矩阵，ρ_i表示粒子i的密度，c₀为常系数；

焦散效应的计算包括：基于Marching Cubes流体表面重构算法及基于屏幕空间的表面重构算法的焦散效应计算，具体的：

所述基于Marching Cubes流体表面重构算法的焦散效应计算包括：计算焦散效应光子的位置P'、焦散效应的影响范围S_b及焦散效应强度C_b；所述焦散效应光子的位置P'的计算公式为：

其中，P表示入射光线I与流体表面的交点，Ref表示入射光线I射入所述流体后发生折射后的光线，Ng表示所述焦散效应光子的位置P'所处平面的法向量，disPtoNg是P到由点Pg和法向Ng构成平面的距离，Ref’为折射光线的方向；disPtoNg＝(P-Pg)·Ng，Pg为所述焦散效应光子的位置P'所处平面上的点，

焦散效应的影响范围S_b的计算公式为：S_b＝c₁·S_m；其中，c₁为比例系数，S_m为Marching Cubes三角形面片面积或像素强度；焦散效应强度C_b的计算公式为：

其中，c₂为比例系数，c₃为指数系数，N为入射光线I与流体表面交点P所处平面的法向量，且N·I＞0；

所述基于屏幕空间的表面重构算法的焦散效应计算包括：在光照空间下绘制粒子球，获得光照表面；恢复所述光照表面在世界坐标系下的位置；在视点空间下计算焦散效应光子的位置；其中，坐标系变换过程的公式为：

v’表示投影后的归一化视点坐标，v表示执行坐标系变换前的世界坐标，M_l与P_l分别表示模式矩阵和投影矩阵。

6.三维流体物理动画引擎的实现装置，其特征在于，该装置包括：

并行处理单元，用于利用并行空间粒子搜索算法搜索流体中每一粒子的邻居粒子，并计算所述每一粒子的状态；

碰撞强度控制单元，用于利用复杂边界查找算法计算流体所处容器的边界，根据所述每一粒子的状态并基于碰撞机制计算与调节所述每一粒子与边界的碰撞强度；

流体表面获取单元，用于根据碰撞后的粒子分布情况并基于流体表面重构算法对流体表面重构，或者基于屏幕空间的表面重构算法对流体表面重构，获得流体表面。

7.根据权利要求6所述的装置，其特征在于，所述并行处理单元包括：

粒子状态计算单元，用于对流体中每一粒子进行三维空间上的哈希操作，实现当前粒子位置到存储地址的映射，并查到找到预设的核函数半径内所述当前粒子的所有邻居粒子；将所述核函邻居粒子的状数半径内所述当前粒子的所有态带入至Navier-Stokes方程组离散对称形式中，并使用多项式函数作为核函数，计算所述当前粒子的密度、压力及粘力，获得所述当前粒子的加速度；使用蛙跳式更新方法计算所述当前粒子的速度和位置状态；

搜索单元，用于搜索所述每一粒子的邻居粒子，且该单元包括：

哈希映射处理单元，用于对流体中的粒子i的位置P_i使用空间哈希操作，得到三维坐标A_i(A_i·x，A_i·y，A_i·z)，其计算公式为：

$A_i\·x=\left[\frac{P_i\·x}{2h}\right];$

$A_i\·y=\left[\frac{P_i\·y}{2h}\right];$

$A_i\·z=\left[\frac{P_i\·z}{2h}\right];$

其中，2h为该流体所处容器的分段长度；

根据所述三维地址坐标A_i（x，y，z），计算一维存储偏移a'_i，其计算公式为：

a'_i＝A_i·x+A_i·y·Dx+A_i·z·Dx·Dy；

其中，Dx，Dy分别表示x，y方向上的最大单元格的数量；

邻居粒子数目查询单元，用于计算所述粒子i所处单元格的三维坐标

根据所述粒子i所处单元格中所维护的粒子链表查找该单元格内的粒子，并判断该单元格内的粒子与所述粒子i之间的距离是否小于h，若是，则为所述粒子i的邻居粒子；

存储空间压缩单元，用于构建一连续的数组M存放所有粒子的邻居粒子，其构建方法为：并行遍历每个粒子，获得每个粒子的邻居粒子个数，存入数组A；计算所述数组A的前缀和，存入数组B；并行遍历所有粒子，将每一粒子的每个邻居粒子按数组B中其位置处的元素作为访问M数组的下标，将邻居粒子的id存入M数组中；

邻居粒子搜索单元，用于通过所述数组M，完成每一粒子的邻居粒子的搜索。

8.根据权利要求6所述的装置，其特征在于，所述碰撞强度控制单元包括：

预处理单元，用于在计算所述流体所处容器的边界之前，进行预处理，具体的包括：将所述容器的边界模型化为点和法向表示的理想平面；对所述理想平面进行离散采样，再进行哈希映射，以确保所述理想平面经过单元格区域；

所述碰撞机制包括：u_i＝u_i-(1+c_R)(u_i·n)n，其中，ui为粒子i的速度，n为平面法向量，c_R为动量损失系数；通过调节动量损失系数来控制粒子与边界碰撞强度。

9.根据权利要求6所述的装置，其特征在于，所述流体表面获取单元包括：

第一流体表面重构单元，用于根据碰撞后的粒子分布情况并基于流体表面重构算法对流体表面重构包括：根据所述流体中的粒子搜索局部SPH粒子，并根据所述局部SPH粒子的位置计算所述流体中粒子分布的各向异性和大小；将所述流体中粒子分布的各向异性和大小作为核函数系数计算每一粒子局部密度场分布，并叠加所有粒子密度场，获得全局密度场分布，再使用移动立方体Marching Cube算法抽取等值面，获得流体表面；

且所述第一流体表面重构单元还包括：匹配立方体有效性扫描单元，用于在使用Marching Cube算法抽取等值面之前进行匹配立方体有效性扫描，且该单元包括：标记单元，用于扫描全局密度场中的所用体素，得到体素的标记数组T[k]，0≤k≤M-1，M为密度场体素个数；映射关系建立单元，用于对所述数组T[k]的前缀和进行计算获得数组T’[k]，则数组T[k]中被标记的总个数为B=T’[M-1]+T[M-1]；建立长度为B的数组V[j]，其中，0≤j≤N-1，且V[T’[i]]=k；

第二流体表面重构单元，用于基于屏幕空间的表面重构算法对流体表面重构包括：将所述流体中的粒子替换为球体；利用计算机图形学中的仿射变换原理，在图像空间下计算表面像素的深度与法向；再对计算得到的所述表面像素的深度与法向进行低通滤波。

10.根据权利要求6或9所述的装置，其特征在于，该装置还包括：光学效应计算单元；所述光学效应包括：包括光照、折射与反射、光线穿透的厚度、焦散和阴影；

所述光学效应计算单元；还包括：

折射与反射系数计算单元，其计算公式为：

r'＝fb+fs·(1.0+I·G)^fp；

r＝min(1.0，r')；

s＝1.0-r；

其中，I为入射光线单位矢量，G为顶点处的法线，fb、fs与fp表示与流体介质相关的多项式系数，r为反射系数，s为折射系数；

流体厚度计算单元，用于计算所述光线穿透的厚度之前计算所述流体的厚度，其计算公式为：

$T_{\mathrm{xy}}=c_0\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_{i}W(M_{\mathrm{mvp}}{P}_i-(x,y),h);$

I_i＝c₀ρ_i；

其中，T_xy表示屏幕空间下x，y位置的流体厚度，M_mvp表示模视投影变换矩阵，ρ_i表示粒子i的密度，c₀为常系数；

焦散效应计算单元，用于进行焦散效应的计算，包括：基于Marching Cubes流体表面重构算法及基于屏幕空间的表面重构算法的焦散效应计算；且该单元包括：

第一焦散效应计算单元，用于基于Marching Cubes流体表面重构算法计算焦散效应，且该单元包括：位置计算单元，用于计算焦散效应光子的位置P'，所述焦散效应光子的位置P'的计算公式为：其中，P表示入射光线I与流体表面的交点，Ref表示入射光线I射入所述流体后发生折射后的光线，Ng表示所述焦散效应光子的位置P'所处平面的法向量，disPtoNg是P到由点Pg和法向Ng构成平面的距离，Ref’为折射光线的方向；disPtoNg＝(P-Pg)·Ng，Pg为所述焦散效应光子的位置P'所处平面上的点，影响范围计算单元，用于计算焦散效应的影响范围S_b，所述焦散效应的影响范围S_b的计算公式为：S_b＝c₁·S_m；其中，c₁为比例系数，S_m为MarchingCubes三角形面片面积或像素强度；强度计算单元，用于计算焦散效应强度C_b，所述焦散效应强度C_b的计算公式为：C_b＝c₂·(|N·I|)^c3；其中，c₂为比例系数，c₃为指数系数，N为入射光线I与流体表面交点P所处平面的法向量，且N·I＞0；

第二焦散效应计算单元，用于基于屏幕空间的表面重构算法计算焦散效应，包括：在光照空间下绘制粒子球，获得光照表面；恢复所述光照表面在世界坐标系下的位置；在视点空间下计算焦散效应光子的位置；其中，坐标系变换过程的公式为：

v’表示投影后的归一化视点坐标，v表示执行坐标系变换前的世界坐标，M_l与P_l分别表示模式矩阵和投影矩阵。

Images