基于细节捕获和形态校正的流体动画渲染方法

摘要

一种计算机图像处理技术领域的基于细节捕获和形态校正的流体动画渲染方法，通过在初始场景上根据Navier-Stokes方程进行流体模拟并对速度场进行优化；然后采用半拉格朗日方法根据优化高精度速度场对相应的密度场和温度场进行更新，以用于渲染和下一帧模拟；最后将更新后的密度场渲染为流体动画。本发明可以通过速度很快的离散正弦变换算子捕获很多绚丽的细节，又可以通过降采样或八叉树方法对原始模拟方法中的最耗费时间的部分降低计算规模，同时对结果进行校正，达到同时加快模拟速度又保持流体模拟细节的作用。与现有技术相比，本发明拥有更快的模拟速度，并且能够获得更加精确的流体动画细节。

说明书

技术领域

本发明涉及的是一种计算机图像处理技术领域的方法，具体是一种基于细节捕获和形态 校正的流体动画渲染方法。

背景技术

在计算机图形学领域，流体动画，包括水，烟，爆炸以及类似的现象，因为其丰富和绚 丽的细节，开始变得越来越流行，使得电影，动画，游戏等领域中都采用增加流体动画的绚丽 程度的手段来吸引眼球。而且由于流体细节的复杂，使得流体模拟变得非常困难，因而吸引了 越来越多的科学家投入流体模拟领域。

现在流体模拟中比较令人信服的方法是欧拉法，由Jos Stam在《SIGGRAPH》上发表的 《Stable Fluids》（1999年）中，对原始的模拟方法进行了改进，在当时来说，获得了更理想的 流体动画细节和更稳定的流体模拟结果，该方法的完整过程是先将流体速度场进行对流，然后 对其应用外力项，最后进行压强项操作，即根据当前速度场和不可压缩条件，构建泊松方程， 求解压强场，并利用求解出来的压强场更新速度场；这三步其实是对Navier-Stokes方程的分解， 其中压强项操作是这个方法的精髓所在，它使得流体模拟的结果更稳定并且更接近与现实中的 不可压缩流体的性质，但是同样的压强项操作也是这种方法的缺陷所在，它变成了基于这种方 法的流体动画的模拟速度最大的瓶颈。在这个方法被提出之后的近十年时间中，大部分流体动 画研究领域中的科学家都将精力集中在对这个方法的改进上，而以该流体模拟方法为基础的改 进研究主要集中在两个方面：一是如何增加流体动画的细节；一是如何加快流体模拟的速度。 到目前为止，流体模拟中并没有效果比较信服同时模拟速度更快的改进方法。

经过对现有技术的检索发现，中国专利文献号CN102339326A，公开日2012-02-01，记 载了一种“分析模拟缝洞型油藏流体流动的方法”，该方法提出了基于Navier-Stokes方程与 Darcy方程耦合的油藏数值模拟方法，主要包括:(1)建立复杂介质Navier-Stokes与Darcy耦合数 学模型;(2)Navier-Stokes与Darcy耦合模型的数值计算方法;(3)将岩体与裂缝分开造型，各自造 型不存在困难，对于裂缝的描述也清晰容易，裂缝数量及空间产状不受限制，特别是不要求裂 缝上的节点与基岩节点完全重合，使得三维工作容易很多。并编写了数值模拟程序来实现本方 法，本方法发展了缝洞型油藏数值模拟理论和方法，科学地实现了对缝洞型油藏的模拟。但该 技术因涉及现实生活的基础设施，对模拟结果的计算精度要求导致该方法的模拟速度比较慢， 不适合流体动画的要求，而且该方法局限于缝洞型油藏流体流动现象，对于一般的流体动画模 拟则适用性不足。

谭捷在《基于物理的流体动画研究》（上海交通大学，硕士论文，2009）中提出了一套 通用的基于层次化网格的多层流体动画框架，用以解决目前流体模拟中传统欧拉方法所面临的 诸多问题，如边界条件的离散、多尺度细节的捕捉等等。但该技术在不同层次之间压强场的处 理理论性不足，而且计算结果与正确解，即最高精度网格上面求解泊松方程的结果，相似性不 足，另外层次化网格的构造非常繁琐。

任威在《大规模三维云实时模拟方法》（计算机辅助设计与图形学学报，2010，22(4)） 中提出一种大规模三维云实时模拟方法。在云建模方面，利用Navier-Stokes流体力学公式模拟 云的动态生成，提出一种基于八叉树的模型化简策略，减少了云模型粒子数;在渲染阶段，提出 一种基于Cell的绘制更新策略，结合Impostor技术自动混合绘制三维云与Impostor，实现了大 规模三维云的实时模拟.实验结果表明，该方法基于物理的方法模拟云，同时在绘制阶段根据视 点的移动实时更新，效果逼真;与同类方法相比，基于Cell的绘制策略更新时计算量更小，有效 地避免了绘制更新时常见的抖动和跳变问题。但该技术中八叉树的模型化简策略只是用于加速 绘制更新步骤，并没有对流体模拟的步骤进行加速，而且该方法同样局限于云流动的现象，不 适用于一般的流体动画模拟。

发明内容

本发明针对现有技术存在的上述不足，提出一种基于细节捕获和形态校正的流体动画渲 染方法，可以通过速度很快的离散正弦变换算子捕获很多绚丽的细节，又可以通过降采样或八 叉树方法对原始模拟方法中的最耗费时间的部分降低计算规模，同时对结果进行校正，达到同 时加快模拟速度又保持流体模拟细节的作用。与现有技术相比，本发明拥有更快的模拟速度， 并且能够获得更加精确的流体动画细节。

本发明是通过以下技术方案实现的，本发明通过在初始场景上根据Navier-Stokes方程进 行流体模拟并对速度场进行优化；然后采用半拉格朗日方法根据优化高精度速度场对相应的密 度场和温度场进行更新，以用于渲染和下一帧模拟；最后将更新后的密度场渲染为流体动画。

本发明具体包括以下步骤：

步骤一，根据需要得到的流体动画的要求，对模拟空间进行网格化，并设置初始场景；

所述的初始场景包括：设置初始速度场、初始密度场和初始温度场、流体模拟过程中的 力场条件和边界条件，以及时间步长。

本步骤确定了流体模拟的过程，流体模拟根据初始速度场不断的得到下一帧符合力场和 边界条件的速度场，然后根据这个速度场去更新密度场，用于渲染，最后得到流体动画。

步骤二，在初始场景上根据Navier-Stokes方程进行流体模拟并通过细节捕获和形态校正 对速度场进行优化，得到优化高精度速度场。

所述的流体模拟是指：在初始场景的速度场上采用半拉格朗日的方法进行Navier-Stokes 方程中的对流项的计算，得到对流模拟后的速度场，并在该对流后的速度场上根据初始场景的 流体模拟过程中的力场条件对速度场进行更新。

所述的优化包括：

细节捕获：对更新后的速度场构建泊松方程，并且利用离散正弦变换算子计算结果，然 后用光滑迭代方法迭代，得到高精度速度场；

形态校正：根据迭代后的高精度速度场计算高精度散度场，然后利用降采样或者八叉树 方法对高精度散度场进行降维操作，得到低精度的散度场，在这个低精度散度场上构建偏差泊 松方程，求解偏差泊松方程，得到低精度的偏差速度场，将这个低精度的偏差速度场上采样之 后与细节捕获得到的高精度速度场进行结合，以校正细节捕获的结果，得到优化高精度速度场。

所述的细节捕获，具体包括以下步骤：

1）通过高精度速度场按照原有压强项的方式构建泊松方程；

2）用离散正弦变换算子求解泊松方程，得到高精度压强场；

3）用光滑迭代方法在高精度压强场的基础上迭代，得到离散正弦变换算子没有捕获完 全的细节，得到完全高精度压强场；

4）用完全高精度压强场对速度场进行更新，得到包含流体细节的高精度速度场。

所述的迭代次数为3-4次；

所述的形态校正，具体包括以下步骤：

i）利用高精度速度场计算每个网格的散度，得到高精度散度场；

ii）将高精度散度场降采样得到规整低精度散度场，并：

a）根据规整低精度散度场构建低精度偏差泊松方程，或

b）根据流体模拟过程中的边界条件构建八叉树，由生成的八叉树的结构分布高精度散 度场，得到以八叉树结构分布的散度场，并在这个散度场上生成偏差泊松方程，该偏差泊松方 程将会比直接在高精度散度场上构建的泊松方程规模要小许多倍；

iii）对偏差泊松方程采用预处理共轭梯度法求解，得到对应的规整低精度压强场或八叉 树结构分布压强场；

iv）利用规整低精度压强场或八叉树结构分布压强场计算偏差速度场，然后将偏差速度 场进行上采样得到偏差高精度速度场，并将其与高精度速度场结合，得到优化高精度速度场。

因为步骤iii）降低了计算规模，所以可以达到加速模拟的目的，对于这个降维之后的偏 差泊松方程，一般采用预处理共轭梯度法来进行求解。

所述的八叉树方法中八叉树的构建是由初始场景中的流体模拟过程中的边界条件确定， 在边界条件复杂的地方，分布八叉树结构的小网格，在流体中则分布大网格；

步骤三，采用半拉格朗日方法根据优化高精度速度场对相应的密度场和温度场进行更新， 以用于渲染和下一帧模拟；最后将更新后的密度场渲染为流体动画。

所述的渲染采用pbrt（基于物理的渲染引擎Physically based rendering toolkit）渲染引擎 进行。

有益效果

1、本发明采用了离散正弦变换算子和将高精度偏差泊松方程降采样或八叉树降维的操 作，这与原来的压强项操作相比，计算算法的效率变高，计算规模成倍数的减少，从而将流体 模拟的速度大幅度的增加；

2、本发明采用了细节捕获操作捕获了许多绚丽的流体动画细节，又通过形态校正来对 流体模拟动画结果进行修正，最终达到保持流体动画细节的目的。

附图说明

图1为本发明流程图。

图2为实施例效果对比图。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施， 给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例1

本实施例中采用了烟上升的流体动画场景，场景被障碍物包住，烟从场景下方的圆孔中， 不断地往上冒，在场景的半空中放置一个圆球，烟从下方圆孔往上冒的过程中遭遇到圆球，并 在圆球附近产生大量的流体细节，造成非常绚丽的流体动画效果，本实施例采用传统流体模拟 框架进行说明。

如图1所示，本实施例包括如下步骤：

步骤一，根据上述的场景设置，将整个场景空间分成128*256*128的网格，设置流体模 拟的初始速度场为0，在场景圆孔附近设置初始的密度场为1，初始的温度场为310，根据场景 要求设置边界条件，即生成一个标识数组obstacle，在圆球占有的网格设置成true，设置时间步 长为Δt=0.05，设置从密度场和温度场产生外力的系数alpha=0.1和beta=1/273.0；

步骤一的作用是设置整个流体模拟的系数和初始状态以及一些用户指定的参数，之后的 流体模拟就会根据这些参数和初始状态生成每帧的密度场，最终用于渲染。

步骤二，在步骤一所设置的速度场上，根据Navier-Stokes方程进行模拟，不断地去更 新速度场，包括如下具体步骤：

第一步，在速度场上，进行Navier-Stokes方程中的对流项的计算，采用半拉格朗日的 方法来进行对流，这一步中产生大部分的凌乱的细节，流体动画的绚丽的部分在这里生成；

半拉格朗日方法是首先将速度场u备份一份u_temp，然后对于u_temp中的每个速度， 分别做如下操作，找到该速度对应的空间位置p，然后用该速度回溯到一时间步长之前的位置 p~=p-u_temp*Δt，然后线性插值得到p～处的速度u_temp~，然后将u_temp~设置成u中对应位 置的速度；

第二步，在第一步对流之后的速度场上，根据力场来更新速度场，这里考虑外力的作用， 使得流体按照用户的需求进行改变；

这里的力场是通过密度场density和温度场temp的共同作用得到，即 force=(alpha*density-beta*（temp-ambient_temp）)*gravity，这里ambient_temp=273.0，gravity=9.8；

而从力场到速度场的更新是采用u+=force*Δt的方式来更新速度场，这种方式是解微分 方程中的显式欧拉方法，具有一阶精度；

第三步，进行细节捕获操作，从这一操作开始就是对原有的压强项操作进行修改，而所 谓的压强项操作是指从当前速度场u，在无散性要求和边界条件的限制下，算出当前 的压强场，然后利用这个压强场的梯度来计算每个网格中u的变化，更新u，这样得到的u既 满足无散性要求又满足边界条件；

所谓无散性要求对网格i，j，k来讲实际上是指u_{i-1/2，j，k}+u_{i，j-1/2，k}+u_{i，j，k-1/2}-u_{i+1/2，j，k}-u_{i，j+1/2，k}-u_{i，j，k+1/2}=0，而边界条件实际上是指速度u在边界的法向上的速度应该等于边界的速度；

压强项操作中，对于每个网格都可以列出一个方程，对于网格i，j，k，该方程的未知 数是它自己的压强p_i，j，k，和它附近网格的压强p_i+1，j，k，p_i，j+1，k，p_i，j，k+1，p_i-1，j，k，p_i，j-1，k，p_i，j，k-1，它们线性组合的结果应该是u_{i-1/2，j，k}+u_{i，j-1/2，k}+u_{i，j，k-1/2}-u_{i+1/2，j，k}-u_{i，j+1/2，k}-u_{i，j，k+1/2}，在这些网 格都为流体的情况下，这些线性系数应该分别是6，-1，-1，-1，-1，-1，-1，如果附近网格有 一个是边界，则将6减1，对应的-1变0，然后对线性组合的结果进行修改，如果该网格正好 是边界，则对应这个网格的方程不需要进行处理，这样就可以得到数目和网格数目相同的方程 组，这里采用预处理共轭梯度法来解这个方程组；

解完方程组得到的压强场p用于更新速度场u，u_{i+1/2，j，k}+=(p_i+1，j，k-p_i，j，k)*Δt，这里也是 显式欧拉方法；

细节捕获操作则是将压强项操作中的预处理共轭梯度法替换成了离散正弦变换算子，这 个算子的特点是速度非常快，而且能够捕获到流体模拟中的大部分细节，在这个的基础上，用 光滑迭代方法迭代三到四次，例如jacobi和超松弛迭代，再多捕获一些细节，但是由于算法的 缺陷，仍然需要接下来的形态校正操作；

第四步，利用第三步的结果，进行形态校正操作，该操作的具体内容是，根据压强项操 作的步骤，再次构建散度场d^f，将d^f降采样到规整的低精度的散度场d^d，低精度散度场的值为 与其网格对应的高精度网格们的值之和；

即：

d^d_i，j，k=d^f_{2*i，2*j，2*k}+d^f_{2*i+1，2*j，2*k}+d^f_{2*i，2*j+1，2*k}+d^f_{2*i，2*j，2*k+1}+d^f_{2*1+1，2*j+1，2*k}+d^f_{2*i，2*j+1，2*k+1}+d^f_{2*1+1，2*j，2*k+1}+d^f_{2*i+1，2*j+1，2*k+1}，

或者，按照构建的八叉树结构生成相应结构的低精度的散度场d^d，其值为在这个八叉树 网格内部所有的高精度网格中心的值之和；

然后，在规整的散度场或者按照八叉树网格分布的散度场上，进行压强项操作，其中散 度场的计算结果则由上述结果代替，还是一般采用预处理共轭梯度法进行求解，计算出压强场 之后，根据显示欧拉方法计算偏差速度场u^b，这个速度场是规整的或者按照八叉树结构分布的；

最后，从规整的或者按照八叉树分布的偏差速度场u^b上，线性插值得到一个规整的高 精度的偏差速度场u^bf，细节捕获的结果速度场u^dc与u^bf之和u^f，即为形态校正的结果，也是这 一帧流体模拟的结果，这一操作主要通过降采样或者八叉树的应用使计算规模下降了许多倍， 从而加快了计算速度，并且对细节捕获的缺陷进行了修正；

八叉树网格的构建决定于步骤一设置的边界条件，具体是标识边界的数组obstacle，构 建八叉树网格首先是得到一个低精度的网格，低精度的网格的大小要比高精度网格的大小要大， 对于低精度网格，如果其中一个大网格同时具有obstacle标识的高精度网格和没有被obstacle 标识的高精度网格，则将这个大网格细分，然后再对细分之后的网格判断是否同时具有obstacle 标识和没有被obstacle标识的高精度网格，如果是则对这个细分之后的网格继续细分，直到细 分到网格和高精度网格一样大小，最终得到按照边界条件分布的八叉树网格。

步骤三，根据步骤二中得到的高精度的速度场u^f，这里采用半拉格朗日的方法来更新密 度场density和温度场temp，用于渲染和下一帧模拟；

步骤四，把步骤三得到的密度场渲染成实际的流体动画，这里采用pbrt渲染引擎进行渲 染。

本实施例的效果如图2所示，图(a)是128*256*128网格上的正确解，即在此精度网格上 求解泊松方程的结果，图(b)是本方法在64*128*64网格上进行形态校正的结果，图(c)是本方法 在32*64*32网格上进行形态校正的结果，图(d)是64*128*64网格上的正确解，图(e)是Lentine 在2010年Siggraph上面发表的《A novel algorithmfor incompressible flow using only a coarsegrid  projection》论文中的方法，其中全局网格是128*256*128，而粗化网格是64*128*64，从图中 可以看出，本方法与现有技术相比更相似于高精度网格上的正确解。

本发明可以应用于流体动画生产领域，能极快地生成想要的流体动画效果且更相似于准 确解，即高精度网格上求解泊松方程的结果，它可以用来作为精确解的预览，以便提前知道最 终效果的大致形式，并调整模拟参数或者边界条件等。