一种基于三变量双调和B样条的流体表面的时空向量化方法与流程

本发明涉及计算机视觉、计算机计算几何、流体仿真技术领域。

背景技术：

近年来，随着虚拟现实领域的快速发展，对自然环境中复杂对象和大规模场景的仿真需求日益增加。在诸多自然场景中，流体和类流体行为是普遍存在的现象之一，例如海浪、河流、瀑布、洪涝、浮云、雾霾、烟雾、火灾等大规模场景，及露珠、水滴、气泡、沸水、油污、水彩、油画、燃烧、血流等小规模场景，均有流体或类流体参与，并呈现出绚丽多变的视觉效果。因此，对流体和类流体行为的建模与仿真，得到了大量优秀研究者的关注，逐渐成为虚拟现实领域研究的重点方向之一。

此外，随着流体仿真技术在影视制作、3d游戏制作、自然现象展示、医学手术教学(模拟血液血管的交互、流血现象等)等方面的广泛应用，考虑到计算机的计算能力与场景建模和制作成本的限制，如何在尽可能保证流体物理规律正确性的前提下提高仿真效率，如何保证流体场景整体和细节上的高视觉逼真度，成为亟需解决的技术问题，也进而促进了流体仿真技术的发展。

计算机行业中，对流体仿真技术的研究主要分为两个方向，计算流体动力学(computerfluiddynamics，cfd)与图形学流体仿真(fluidanimation)。其中，计算流体动力学，以科学严谨的方式，研究流体运动的物理规律，进行精确的数值求解，其主要目的，是解决工业界流体相关场景等高精确度仿真要求环境下的流体建模和模型求解问题。其后，随着计算机图形学技术的进步和计算机成像效果的提升，以追求仿真效果视觉逼真度为目的的图形学流体仿真进入快速发展阶段。1996年，首次实现了基于navier-stokes方程的图形学流体仿真，这种不要求精确的流体动力学方程求解，但同时又能保证视觉逼真效果的流体仿真方法，迎合了影视、游戏制作等娱乐行业的发展需要，并同时符合医学领域教学仿真需求。

尽管已经拥有二十多年的发展历史，流体仿真仍然是一个具有挑战性的计算机图形学问题，这是由流体现象的复杂性造成的。基础仿真方法，虽然可以对较简单的自然现象进行较好的仿真，但对于自然界中存在的复杂自然现象和大规模流体场景的精确仿真，现存方法还不能进行准确描述。单纯的物理建模难以满足仿真需求，因此，大量学者将数据驱动的方法引入到了流体仿真技术中，希望以真实流体的运动状态驱动流体仿真的进行。但由于流体不具有固定的几何结构，流体数据的获取和分析，又成为数据驱动思路中新的难点。

针对科研和娱乐市场的大量需求，对于高精度高效率的流体建模和仿真方法的研究具有重要意义与价值，其研究成果将在实际生活中得到广泛应用。

技术实现要素：

为了克服上述现有技术的不足，本发明基于三变量双调和b样条的流体表面的时空向量化方法，该方法通过对流体视频的分析，提取流体的高度信息和速度信息，并通过聚类、数据拟合等方法，对流体信息进行归类和简化，并以此为驱动，实现新流体场景的仿真。

为完成发明目的，本发明采用的技术方案是：一种基于三变量双调和b样条的流体表面的时空向量化方法，其具体步骤如下：

一种基于三变量双调和b样条的流体表面的时空向量化方法，包含以下步骤：

(1)基于sfs的流体视频高度场还原：利用shapefromshading方法重建流体表面高度场，在迭代过程中，使用gaussianfilter和guidedfilter对高度场进行去噪，进行保留边界的平滑优化，还原流体表面几何形状；

(2)基于浅水方程的流体速度场还原：使用浅水方程和平滑项建立流体高度场和速度场关联的优化方程，求解该优化，获得流体表面在水平面上的速度场；

(3)基于三变量双调和b样条的数据拟合：将流体表面的高度场信息和速度场信息视作三维体数据，进行超体素分割，利用分割结果选取knot点，构建三变量双调和b样条基底，对流体表面信息进行向量化，实现时空离散数据到连续数据的映射；

(4)基于粒子表示的三维流体重现：使用粒子表示流体的基本单元，利用连续化结果，构建三维流体，重现流体的运动。

进一步地，在步骤(1)中，使用shapefromshading方法，基于对图像辐射方程的迭代拟合，对流体视频的每一帧的高度场进行迭代分析，使用gaussianfilter和guidedfilter结合的方法，在每次迭代步中，对分析得到的高度场进行滤波去噪，去除分析结果中的噪声点和异常值。

进一步地，在步骤(2)中，以浅水方程为基本物理模型，分析获得流体高度场和速度场关系的描述方程，以该方程为基础，添加带权重的平滑项，构建最小化目标函数，求解目标函数的欧拉-拉格朗日方程，，使用迭代求解方法，实现目标函数的最小化，使用步骤(1)中分析获得的流体表面高度信息，即可获取流体表面的二维速度场信息。

进一步地，在步骤(3)中，将流体的速度场信息可视为以视频帧的横纵坐标为x、y轴，以时间为z轴的三维体数据，对该体数据进行超体素分割，并计算超体素内速度向量的平均向量，并将超体素中心平移到坐标原点，计算获取几组三变量双调和b样条基，对每一个超体素内的速度数据进行拟合，获得拟合系数向量组，实现对流体表面信息的向量化。

进一步地，在步骤(4)中，使用粒子表示流体的基本单元。

进一步地，在步骤(1)中，以去噪和建立更好的迭代基础为目的的guidedfilter在前几步迭代过程中使用，以去噪和平滑为目的的gaussianfilter在后几步迭代过程中使用，guidedfilter和gaussianfilter的使用次数均可由用户决定。

与现有技术相比，本发明的优点和积极效果在于：

1、本发明通过shapefromshading方法，对单目流体视频中的流体高度场进行分析，并给予gaussianfilter和guidedfilter对高度场的提取过程进行优化。同时，本发明实现了基于cuda的算法优化，利用gpu实现了对视频流体高度场的实时分析。

2、本发明通过浅水方程和欧拉-拉格朗日方程，实现了基于流体高度场的流体表面速度场分析，实现了对流体表面运动的全面描述；

3、本发明通过将基于视频获取的流体运动信息视做流体演化过程的离散采样，利用三变量双调和b样条对该离散数据进行拟合，实现了流体运动演化在时空方向的连续化表示。

4、本发明基于流体向量化结果，实现了对数据的重用，并构建了流体视频对应的三维流体场景。

附图说明

图1是本发明基于三变量双调和b样条的流体表面的时空向量化方法所述方法流程图；

图2是本发明的基于sfs及不同次filter去噪的流体视频高度场还原结果图；

图3是本发明的基于sfs的流体视频高度场还原过程和基于浅水方程的流体速度场还原过程中不同迭代次数的时间统计图；

图4是本发明的基于三变量双调和b样条的数据拟合结果示意图；

图5是本发明的基于粒子表示的三维流体重现方法的粒子初始化示意图；

图6是本发明的基于三变量双调和b样条的数据拟合前后，相邻帧之间的时间导数示意图；

图7是本发明的基于三变量双调和b样条的数据拟合的时空超体素化结果；

图8是本发明的基于三变量双调和b样条的数据拟合的空间超体素化结果；

图9是本发明的基于粒子表示的三维流体重现渲染结果(微小波纹场景)；

图10是本发明的基于粒子表示的三维流体重现渲染结果(海浪场景)。

具体实施方式

本发明提供一种基于三变量双调和b样条的流体表面的时空向量化方法，包括：输入流体样例视频，重建流体表面高度场，还原流体几何形状；利用潜水方程和平滑约束，对流体的表面速度场进行拟合；将流体表面的几何和速度信息，视作时空三维的体数据，使用三变量双调和b样条作为时空基底，分别拟合流体的几何和速度数据，实现离散空间到连续空间的映射；基于时空连续化结果，采用粒子作为流体的基本单元，对流体的运动进行重现，并提取流体表面进行渲染，获得类似输入视频中流体运动模式的流体仿真结果。

下面结合具体方法步骤，介绍本发明原理。

步骤一：基于sfs的流体视频高度场还原：以dyntex数据库中的海洋视频为流体视频来源，使用shapefromshading方法，基于对图像辐射方程的迭代拟合，对流体视频的每一帧的高度场进行迭代分析。使用gaussianfilter和guidedfilter结合的方法，在每次迭代步中，对分析得到的高度场进行滤波去噪，去除分析结果中的噪声点和异常值。其中，基于shapefromshading的流体高度场分析是通过迭代实现的，以去噪和建立更好的迭代基础为目的的guidedfilter在前几步迭代过程中使用，以去噪和平滑为目的的gaussianfilter在后几步迭代过程中使用，guidedfilter和gaussianfilter的使用次数均可由用户决定。

步骤二：基于浅水方程的流体速度场还原：以浅水方程为基本物理模型，分析获得流体高度场和速度场关系的描述方程，以该方程为基础，添加带权重的平滑项，构建最小化目标函数，求解目标函数的欧拉-拉格朗日方程(euler-lagrangianequation)，使用迭代求解方法，实现目标函数的最小化，获取流体表面的二维速度场信息。其中，“浅水方程”是描述流体运动的三维n-s方程的二维近似形式，欧拉-拉格朗日方程是变分法中的经典能量极小化的求解方法。

步骤三：基于三变量双调和b样条的数据拟合：流体的速度场信息可视为以视频帧的横纵坐标为x、y轴，以时间为z轴的三维体数据，对该体数据进行超体素分割，并计算超体素内速度向量的平均向量，并将超体素中心平移到坐标原点。计算获取几组三变量双调和b样条基，对每一个超体素内的速度数据进行拟合，获得拟合系数向量组。其中，三变量双调和b样条是一组三维空间的基，三变量双调和b样条具有局部化特性，能够实现单位分解，并具有一定的鲁棒性和稳定性。

步骤四：基于粒子表示的三维流体重现：使用粒子表示流体的基本单元，利用连续化结果，构建三维流体，重现流体的运动。

本发明只需要使用单目流体视频作为输入，自动分析视频中流体的几何和速度信息，并利用三变量双调和b样条，实现了流体表面几何和速度在时间和空间维度上的连续化，利用连续化结果，可实现流体表面在时间和空间维度的超分辨还原；此外，本发明基于粒子实现了流体的重仿真，将流体视频还原为三维流体，进一步证明了实现视频流体运动特征驱动的物理仿真的可能性。

为了使本发明的目的、技术方案及优点更加清楚明白，下面结合附图，对本发明的方法作详细解释说明。应当理解，此处所描述的具体实例仅仅用以解释本发明，并不用于限定本发明。

本发明提出了一种框架，通过对流体视频的分析，提取流体的高度信息和速度信息，并通过聚类、数据拟合等方法，对流体信息进行归类和简化，并以此为驱动，实现了流体三维场景的重建与仿真。

本发明提供了一种基于三变量双调和b样条的流体表面的时空向量化方法，参考图1整体方法流程图所示，具体实施方式如下：

步骤一：基于sfs的流体视频高度场还原：

shapefromshading方法的实质，是根据图像的辐射度方程(公式(1))，建立图像亮度值和位置高度值之间的关系，进而通过图像的亮度值，计算每个像素位置的高度值。公式(1)中，r(n(x))＝ω·n(x)表示辐射方程，ω为入射光照方向，n(x)表示表面法向，计算方法如公式(3-2)所示，p(x)和q(x)分别表示高度场在x和y方向上的空间导数。

r(n(x))＝i(x)(1)

参考tsai和shah提出的直接拟合辐射度方程的方法，使用有限差分方法，近似计算高度场在x和y方向上的空间导数p(x)和q(x)，因此，若将公式(1)中的辐射方程项，写作以p(x)和q(x)为自变量的形式，则可得公式(3)。其中，下角标i,j表示相应物理量在图像中的位置坐标为(i,j)。

r(pi,j,qi,j)＝r(ui,j-ui-1,j,ui,j-ui,j-1)＝ii,j(3)

本节中使用u表示高度场，ui,j表示(i,j)位置的高度值，用表示当前迭代步开始前对高度场的分析值，用表示当前迭代步计算后对高度场的分析值(更加精确的分析值)。公式(3)是一个一维非线性方程，可使用牛顿迭代法(newton-raphsonmethod)思想，构造公式(3)的迭代函数(如公式(4)所示)，用于实现对高度场ui,j的求解。

其中，对函数f增添角标表示对f的自变量增加相应角标，即：但需要注意的是，公式(4)不是严格的牛顿迭代法给出的迭代函数，因为ui-1,j和ui,j-1在每个迭代步结束后都会改变，因此，公式(4)只是通过借鉴牛顿迭代法的求解思想及求解形式得到的，对于这样的迭代求解是否严格收敛，目前还没有普适的证明。将迭代函数(即公式(4))写成显示的迭代方程，可得如下公式(公式(5)):

由于没有普适的证明保证公式(5)收敛，tsai和shah等人提出如下改进的迭代方程形式，如公式(6)所示。式中的应满足在非零时接近原值，否则为0，因此，使用辅助量将表示成如公式(3-7)所示的形式，式中的w取值为0.01，用于保证分母不为0，辅助量的迭代计算方法如公式(3-8)所示。

根据上述分析，可得高度场的迭代计算算法，如下所示。

为了消除视频各帧高度场的竖直平移，在迭代计算流体高度场后，需要对高度值进行归一化，将视频每帧的高度场的平均值归一化到0。

步骤一实验结果如图2所示。图中使用灰度图像和3d重建场景来表示流体的高度场。如图2中的无guassianfilter和无guidedfilter子图中的方框圈出的位置所示，高度场的重建过程会产生累计误差和微小误差。为了消除显著的累计误差和类似噪声点的微小误差，采用guidedfilter和gaussianfilter对每一次迭代结果进行优化。图2展示了使用不同次数的guidedfilter和gaussianfilter进行优化的实验结果。不使用guidedfilter时，累计误差难以得到有效消除，不适用gaussianfilter时，微小误差难以得到有效消除。

步骤一的时间统计结果如图3所示。图中给出了四个场景(波纹、雨滴、瀑布、喷泉场景)在不同迭代次数下的运行时间统计。从图中可以看出，步骤一提出的方法能够实现实时运行。

步骤二：基于浅水方程的流体速度场还原；

为了获取能够驱动流体表面运动的速度场，采用模型拟合的方法，根据上述高度场，利用浅水方程，对速度场进行还原。

浅水方程的基本方程组中，能够反映高度场和速度场关系的方程式如公式(9)所示。其中，z表示高度场，u＝(u,v)表示未知速度场。

zt＝-u·▽z-z(ux+vy)(9)

将公式(9)展开为更加一般的形式，如公式(10)所示。

zt+zxu+zyv+z(ux+vy)＝0(10)

直观理解上述公式的物理含义，流体表面高度场随时间的演化zt可由两个子项进行建模，首先，高度场在水平面(与高度场方向垂直的表面)上的速度场的平流输送下进行更新，由子项-(zxu+zyv)表达；其次，高度场的变化与速度场的二维散度场的变化比例相关，由子项-z(ux+vy)表达。

除了物理控制项外，为了保证获取的速度场的平滑性，还应在优化目标函数中加入平滑项(|▽u|²+|▽v|²)，通过调节平滑项全权重，获得预期的速度场提取效果。综上所述，优化目标函数的形式如公式(11)所示。

e＝∫∫[(zt+zxu+zyv+z(ux+vy))²+α²(|▽u|²+|▽v|²)]dxdy(11)

对公式(11)的优化可通过求解该优化方程对应的欧拉-拉格朗日方程实现。欧拉-拉格朗日方程可由如下线性系统表达：aub^t＝cu，avb^t＝cv，其中b＝(u,v)，表示待求解的速度场，au和av的展开形式如公式(12)和(13)所示，cu和cv。的展开形式如公式(14)和(15)所示；设w＝{u,v}，和分别表示w在(x,y)处沿x和y方向的平均值，其计算公式如公式(15)和(16)所示。

通过公式可得据此，使用迭代方法逐渐逼近真实解，得到如下迭代公式(公式(3-17)和公式(3-18))：

利用上述迭代公式，可根据流体表面高度场，对流体表面速度场进行拟合，同时，使用高速滤波方法，去除上述拟合结果中的异常值和噪声点。

步骤二实验时间统计结果如图3所示。图中给出了四个场景(波纹、雨滴、瀑布、喷泉场景)在不同迭代次数下的运行时间统计。从图中可以看出，步骤二提出的方法能够实现实时运行。

步骤三：基于三变量双调和b样条的数据拟合；

首先，将流体的速度场信息可视为以视频帧的横纵坐标为x、y轴，以时间为z轴的三维体数据，且每个(x,y,t)位置的数据应是一个二维向量(u,v)。参考基于线性迭代聚类的图像超像素分割算法，对速度场体数据进行超体素分割。

传统的基于线性迭代聚类的图像超像素分割算法，使用lab颜色空间坐标和二维坐标构成的五维空间上的欧式距离，来度量像素点之间的距离。类似的，本文采用速度场和三维位置组成的五维空间刻画每个体素点，如公式(19)所示。使用上的欧式距离，来度量体素之间的距离，距离ds可通过公式(20)得到。式中，m用于控制超体素的紧凑性(m越大，超体素越紧凑)，表示超体素的网格间隔，其中n表示体素点总个数，k表示超体素的大致数目。

当前步骤的预期结果为分割好的超体素结果，其可视化结果大概类似二维图像超体素分割结果，不同的是超体素分割对应于三维体素数据，而超像素分割对应于二维像素数据。

超体素分割完成后，根据超体素的位置进行knot点选择。在每个超体素内，使用随机采样方法，选取m个固定的knot点，若k表示超体素的个数，则全局来看，共有k·m个knot点，若表示knot点集合，即有

计算所有knot点位置处的三变量双调和b样条基函数。双调和b样条基函数的离散形式如公式(21)所示，其中，p(x)表示格林函数，基函数是由knot点处的格林函数的加权组合构成的，系数的求解方法如公式(3-22)所示，h的每一列即是格林函数的系数。

公式(22)中，表示knot点集合，表示第j个knot点的一环邻接knot点(即直接相邻的邻居knot点)，表示格林函数集合。

使用固定的m个基函数，对速度场值(u,v)进行拟合。用2m维的拟合向量表示每个体素值，用于表征速度场。第i个点的系数向量可通过如下优化方程(公式(23))求得，其中，

步骤三的实验结果如图4所示。图4的第一列是从原输入视频中截取的示例帧。对于每一个示例帧，图4进一步给出了步骤一重建的高度场结果，步骤二重建的速度场结果的灰度示意图和3d重建结果示意图，后续列相继给出了超体素分割结果、基于步骤三方法的高度、速度场重建结果(灰度图、3d示意图，及拟合残差示意)。

基于三变量双调和b样条的数据拟合前后，相邻帧之间的时间导数如图6所示，从图中可以看出，拟合后的结果在一定程度上消除了时间上的不连续噪声；时空超体素化结果如图7所示，图7以微小波纹为例，给出了时间和空间维度的超体素空间化重建结果，图中的数字代表相邻两帧间的时间导数值，该值反映了步骤三方法能够相对连续的重建流体的高度场和速度场；空间超体素化结果如图8所示，重建图像的分辨率列在图像正上方，空间超体素化结果显示了步骤三方法能够在保持流体细节的前提下，重建连续的流体空间。

步骤四：基于粒子表示的三维流体重现；

为了重用流体运动演化过程向量化后的数据，本发明采用使用粒子表示流体的基本单元，对流体的进行三维的重仿真，其中，使用流体表面的高度场和速度场是向量化的流体演化结果的超分辨结果作为流体粒子运动的向导。基于粒子表示的三维流体重现方法的粒子初始化结果如图5所示。使用marchingcube算法对每一帧的流体运动结果进行表面提取，并使用blender软件进行流体表面渲染，获得如本发明图9和图10所示的渲染结果示意图。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例应用于其它领域，但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。