CN101577008A - 一种快速且精确的大理石纹釉生成方法 - Google Patents

Abstract

本发明公开了一种快速且精确的大理石纹釉生成方法，用户通过基于CUDA的仿真系统执行各种操作，仿真系统执行相应的算法来进行大理石纹釉的仿真制作，方法包括以下步骤：用户设置相应参数；用户选择需要的操作方式在二维运算域中进行相应的操作；通过求解二维流体的动力学方程得到二维运算域中大理石纹釉每层图层的速度场；利用USCIP高阶插值方法，通过转换函数以及逆转换函数对速度场进行转换计算，得到大理石纹釉每层图层的密度场；每层的密度场和对应的颜料色彩相乘，并将所有层的颜色值累加得到最后大理石纹釉图像的颜色值，生成用户需要的大理石纹釉。本发明算法明确，界面友好，结果鲁棒，可以更快更好的生成大理石纹釉图案。

Description

一种快速且精确的大理石纹釉生成方法

技术领域

本发明涉及计算机图形学中的流体仿真领域，特别是涉及一种快速且精确的大理石纹釉的生成方法。

背景技术

快速并且精确的大理石纹釉生成系统通过数值方法求解二维流体动力学的Navie-Stokes方程，提供了一种用计算机来模拟制作大理石纹釉这门传统艺术的工具。

大理石纹釉工艺起源于11世纪的西南亚地区，但至今仍是一门受全世界人民喜爱的手工艺。其制作流程的基本过程是：对悬浮在基底液表层的颜料进行滴液，搅动，梳理等操作，然后将得到的图象转移到其他具有吸附性材料上进行保存。由于在很大程度上受到各种物理条件的限制，因此制作大理石纹釉并不是一样轻松的事。近年来，涌现了很多用计算机来模拟大理石纹釉工艺的仿真系统。这些系统将人类从物理的局限性中解放出来，人们通过设置一些流体的参数，随意搭配色彩方案，定制个性化的操作方式，甚至还可以撤销或重复设计过程中的每一步，就可以轻松的设计出漂亮的大理石纹釉图案。

2001年，研究人员首次提出用物理方法来模拟传统大理石纹釉工艺，参见Suzuki，T.，Mao，X.，Imamiya，A.：Simulating marbling with computergraphics.Proc.IASTED Interna-tional Conf.Visualization，Imaging，andImaging Processing，pp.208-213.IASTED，Calgary(2001)。2003年茅晓阳等人将其扩展成用户交互的大理石纹釉生成系统，该系统基于传统的大理石纹釉工艺的物理模型，制作过程可以看作一个二维流体动力学问题，并用数值方法求解。参见Mao，X.，Suzuki，T.，Imamiya，A.：AtelierM：Aphysically based interactive system for creating traditional marbling textures.Proc.1st International Conf.Computer Graphics and Interactive Techniques inAustralasia and South East Asia(GRAPHITE’03)，pp.79-86.ACM，NewYork(2003)。此后Acar和Boulanger等研究人员提出了用多图层的流体模型来模拟这门工艺，参见Acar，R.，Boulanger，P.：Digital marbling：amultiscale fluid model.IEEE Trans.Visualization and Computer Graphics，12(4)，pp.600-614(2006)。然而由于这些方法都是建立在CPU上的，运行速度较慢，用户需要等待数秒甚至数分钟才能看到模拟的结果。后来金小刚等人提出了在GPU上求解Navie-Stokes方程，虽然在速度上基本能达到实时性的要求，但却遇到了边缘模糊的问题，在真实的大理石纹釉工艺中，牛胆汁常常被添加到颜料中以在其表面形成脂肪层，从而使各种色彩不至于混合到一起，参见Jin，X.，Chen，S.，Mao，X.：Computer-generatedmarbling textures：a GPU-based design system.IEEE Computer Graphics andApplications，27(2)，pp.78-84(2007)。此后，许佳奕等人采用复杂的高阶水平对流模式改进纹釉的锐化效果，参见Xu，J.，Mao，X.，Jin，X.：Nondissipative marbling.IEEE Computer Graphics and Applications，28(2)，pp.35-43(2008)。然后高阶的方法不仅没有彻底解决边缘模糊问题，而且引入了流体产生不稳定的新问题。除此之外，以上的大理石纹釉仿真系统均不能实时生成高分辨率的大理石纹釉图案。

发明内容

本发明提供了一种快速且精确的大理石纹釉生成方法。该方法通过用数值方法求解二维流体的动力学Navie-Stokes方程，有效地解决了现有的基于物理的实时大理石纹釉生成系统中存在的问题，本发明基于最新的图形硬件并行计算架构CUDA，为艺术家和非艺术家设计大理石纹釉图案提供了一个非常便利的工具。

一种快速且精确的大理石纹釉生成方法，用户通过基于CUDA的仿真系统执行各种操作，仿真系统执行相应的算法来进行大理石纹釉的仿真制作；

方法包括以下步骤：

(1)用户设置用于大理石纹釉生成的二维运算域和基底液的属性和大理石纹釉的图层的层数n、图层的颜色color[n]；用户可以根据需要设置大理石纹釉的图层的层数，第n层的颜色为color[n]。

用户输入液滴的半径R和搅动流体时的外力大小F两个参数的数值。

大理石纹釉的制作是在一个浅的方形平底碟中进行的，在计算机上制作仿真大理石纹釉时在仿真系统中定义一个二维运算域作为平底碟的内部底平面，二维运算域是一个长为M个单位长度，宽为N个单位宽度的方形栅格上，栅格内部由M×N个正方形单元格构成，每个单元格的大小由用户在初始的时候设定。

二维运算域上设有基底液，基底液的属性包括颜色，即纹釉的底色，以及基底液流体的粘稠度。颜料在不同粘稠度的基底液表面的扩散速度不同，大理石纹理工艺所用的基底液一般具有非常高的粘性。传统的大理石纹理工艺的基底液和普通的水不同，里面添加了一些溶剂，例如牛胆汁，是大理石纹理制作工艺所专用的。

(2)采用用户指定的操作方式在二维运算域中进行相应的操作，所述的操作包括滴入各种颜色的滴液，搅动基底液带动各种颜色的图层运动和梳理图层的纹理。滴液的颜色用户可以按照自己的需要进行选择，搅动基底液所用的力也可以按照需要设定。

步骤(1)以及步骤(2)的操作，用户可以利用现有技术对相关参数进行预设定，并将预设定的参数(带有颜色信息、操作信息、以及矢量力场信息的数据)输入到二维运算域中。滴液操作时，同时要选择所要滴液的图层及其颜色。搅动和梳理基底液操作时，同时要选择力度大小，另外是梳子的齿形。选择参数完毕后，用鼠标或压杆笔进行纹釉设计。

(3)基于步骤(2)的操作，通过求解二维流体的动力学Navie-Stokes方程得到二维运算域中大理石纹釉每层图层的速度场u包括以下步骤：

所述的二维流体的动力学Navie-Stokes方程如(1)式：

$\frac{\∂u}{\∂t}=-(u\·\▿)u-\frac{1}{\mathrm{\ρ}}\▿p+v{\▿}^{2}u+F$

$\▿\·u=0...\left(1\right)$

式中u表示流体的速度，ρ表示流体的密度，v为粘度，F是作用于流体的外力，式子右边的四个项都是加速度，

是水平对流项，是压力项，是粘性扩散项，F是外部作用力项。

1)求解二维流体的动力学Navie-Stokes方程的水平对流项；

将Navie-Stokes方程的水平对流项

表示成隐式表达方式，方程如(2)式：

u(x，t+δt)＝u(x-u(x，t)δt，t)......(2)

式中u表示速度，x表示位置，t代表某一时刻，δt表示时间步长；

为了计算二维运算域中的x位置在t+δt时刻的速度，采用准拉格朗日方法沿着粒子的运动轨迹回溯到在t时刻的位置x′，对最接近x′的四个单元格用USCIP方法进行插值得到的速度拷贝到x位置，最后得到整个二维运算域中图层的水平对流速度场u₁。

无论时间步长δt取多大，准拉格朗日方法都能使更新后的速度都不会超过整个速度场在上一时刻的速度上限，从而保证了算法的稳定性。而且USCIP是一种稳定准确而且快速的高阶插值方法，并且非常适合于在GPU上实现。

2)求解二维流体的动力学Navie-Stokes方程的粘性扩散项；

将Navie-Stokes方程的粘性扩散项

转换为隐式表达式，方程如(3)式：

$(I-\mathrm{v\δt}{\▿}^2)u(x,t+\mathrm{\δt})=u(x,t)...\left(3\right)$

式中，I为单位矩阵，v为粘度，u为速度，t代表某一时刻，δt表示时间步长；

方程(3)可以看作是泊松方程，利用多栅解法在GPU上求解方程得到整个二维运算域中图层的粘性扩散速度场u₂；

多栅解法先在原始栅格上求初始值，然后将每个栅格的余数投影到粗糙一级的栅格继续求初始值和余数直到达到第N级栅格，最后逐次向高分辨率级别插值，得到最终解。

3)求解二维流体的动力学Navie-Stokes方程的外力项，将外力场直接叠加到速度场上，方法如下：

记录用户在搅动液体时经过路径的切线方向，并将所述的切线方向作为外力F的方向，外力F的大小由用户通过仿真系统交互设置，然后直接将外力场叠加到步骤1)和步骤2)的总速度场(u₁+u₂)上得到新的速度场w；

所述步骤1)和步骤2)的总速度场为水平对流速度场u₁和粘性扩散速度场u₂的叠加，即(u₁+u₂)；

4)求解二维流体的动力学Navie-Stokes方程的压力项，将步骤3)得到的速度场w减去压力的梯度得到新的速度场u的方法如下：

经过水平对流，粘性扩散，施加外力三个步骤后得到的新的速度场w是发散的，要将w转化为一个无发散的速度场u；

由于任何一个矢量场都可以被分解成另外两个矢量场的和，其中一个为无发散的矢量场，另一个为标量场的梯度，w分解如(4)式所示：

$w=u+\▿p...\left(4\right)$

式中p是Navie-Stokes方程中的压力场，对公式(4)两边应用发散度算子后根据Navie-Stokes方程(1)中的 $\▿\·u=0$ 的条件得到压力泊松方程(5)：

${\▿}^{2}p=\▿\·w...\left(5\right)$

采用多栅解法来解这个压力泊松方程，最后用w减去p的梯度就得到二维运算域中大理石纹釉每层图层的速度场u。

(4)利用USCIP高阶插值方法，通过转换函数以及逆转换函数对速度场进行转换计算，得到大理石纹釉每层图层的密度场，具体包括以下步骤：

密度场和速度场之间的关系如公式(6)所示：

$\frac{\∂d}{\∂t}=-(u\·\▿)d...\left(6\right)$

式中d为密度场，u为速度场；

为了计算x位置在t+δt时刻的密度，先将t时刻密度场用如式(7)所示的转换函数进行转换，将d转换为F(d)：

$F\left(d\right)=0.5\log \left(\frac{\max (d,\mathrm{\ξ})}{\max (1-d,\mathrm{\ξ})}\right)...\left(7\right)$

式中d表示密度，ξ是一个很小的值，目的是为了防止出现分母为0的情况；将转换后的F(d)代入到(6)式，再采用准拉格朗日方法沿着粒子的运动轨迹回溯到在t时刻的位置x′，对最接近x′的四个单元格用USCIP方法进行插值得到的密度值拷贝到x位置，得到新的F(d)，最后将新的F(d)再通过式(8)所示的逆转换函数转换得到最终的密度场d。

F^-1(d)＝0.5(tanh(d)+1)......(8)

(5)将大理石纹釉每层图层的密度场和该层的颜色相乘得到该层的颜色值，并将大理石纹釉所有层的颜色值累加得到最后大理石纹釉图像的颜色值，生成用户需要的大理石纹釉；

大理石纹釉图像的颜色值按照式(9)进行计算：

$C_{\mathrm{output}}=C_{\mathrm{liguid}}*(1-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{d}_i)+\underset{i=0}{\overset{n}{\mathrm{\Σ}}}(C_i*d_i)...\left(9\right)$

式中，C_output是输出结果，C_liguid表示基底液的颜色，d_i表示第i层的密度值，C_i表示第i层的颜色值。

本发明针对现有的基于物理模型的大理石纹釉工艺的仿真系统存在的各种缺点，采用了一种快速并且精确的方法，可以更快的，更好的生成大理石纹釉图案，本发明算法明确，界面友好，结果鲁棒。

附图说明

图1为本发明的技术方案流程图；

图2为本发明中泊松方程多栅解法的过程图。

具体实施方式

下面结合附图对本发明的一种快速且精确的大理石纹釉生成方法的实施例进行详细说明。

如图1所示，一种快速且精确的大理石纹釉生成方法，用户通过基于CUDA的仿真系统执行各种操作，仿真系统执行相应的算法来进行大理石纹釉的仿真制作，方法包括以下步骤：用户设置相应参数；用户选择需要的操作方式在二维运算域中进行相应的操作；通过求解二维流体的动力学Navie-Stokes方程得到二维运算域中大理石纹釉每层图层的速度场；利用USCIP高阶插值方法，通过转换函数以及逆转换函数对速度场进行转换计算，得到大理石纹釉每层图层的密度场；每层的密度场和对应的颜料色彩相乘，并将所有层的颜色值累加得到最后大理石纹釉图像的颜色值，生成用户需要的大理石纹釉。

(1)用户设置相关参数，包括用于大理石纹釉生成的二维运算域和基底液的属性和大理石纹釉的图层的层数n、图层的颜色color[n]；用户可以根据需要设置大理石纹釉的图层的层数，第n层的颜色为color[n]。

用户输入液滴的半径R和搅动流体时的外力大小F两个参数的数值。

大理石纹釉的制作是在一个浅的方形平底碟中进行的，在计算机上制作仿真大理石纹釉时在仿真系统中定义一个二维运算域作为平底碟的内部底平面，二维运算域是一个长为M个单位长度，宽为N个单位宽度的方形栅格上，栅格内部由M×N个正方形单元格构成，每个单元格的大小由用户在初始的时候设定。

二维运算域上设有基底液，基底液的属性包括颜色，即纹釉的底色，以及基底液流体的粘稠度。颜料在不同粘稠度的基底液表面的扩散速度不同，大理石纹理工艺所用的基底液一般具有非常高的粘性。传统的大理石纹理工艺的基底液和普通的水不同，里面添加了一些溶剂，例如牛胆汁，是大理石纹理制作工艺所专用的。

(2)采用用户指定的操作方式在二维运算域中进行相应的操作，所述的操作包括滴入各种颜色的滴液，搅动基底液带动各种颜色的图层运动和梳理图层的纹理。滴液的颜色用户可以按照自己的需要进行选择，搅动基底液所用的力也可以按照需要设定。

(3)基于步骤(2)的操作，通过求解二维流体的动力学Navie-Stokes方程得到二维运算域中大理石纹釉每层图层的速度场u包括以下步骤：

所述的二维流体的动力学Navie-Stokes方程如(1)式：

$\frac{\∂u}{\∂t}=-(u\·\▿)u-\frac{1}{\mathrm{\ρ}}\▿\mathrm{\ρ}+v{\▿}^{2}u+F$

$\▿\·u=0...\left(1\right)$

式中u表示流体的速度，ρ表示流体的密度，v为粘度，F是作用于流体的外力，式子右边的四个项都是加速度，

是水平对流项，

是压力项，

是粘性扩散项，F是外部作用力项。

1)求解二维流体的动力学Navie-Stokes方程的水平对流项；

将Navie-Stokes方程的水平对流项表示成隐式表达方式，方程如(2)式：

u(x，t+δt)＝u(x-u(x，t)δt，t)......(2)

式中u表示速度，x表示位置，t代表某一时刻，δt表示时间步长；

为了计算二维运算域中的x位置在t+δt时刻的速度，采用准拉格朗日方法沿着粒子的运动轨迹回溯到在t时刻的位置x′，对最接近x′的四个单元格用USCIP方法进行插值得到的速度拷贝到x位置，最后得到整个二维运算域中图层的水平对流速度场u₁。无论时间步长δt取多大，准拉格朗日方法都能使更新后的速度都不会超过整个速度场在上一时刻的速度上限，从而保证了算法的稳定性。

USCIP是一种稳定准确而且快速的高阶插值方法。对于二维的液体，USCIP多项式公式如下.

$\mathrm{\φ}(x,y)=\underset{0\≤i+j\≤3}{\mathrm{\Σ}}{C}_{\mathrm{ij}}{x}^i{y}^j+C_{31}{x}^{3}y+C_{13}x{y}^3$

φ₀₀，φ₀₁，φ₁₀，和φ₁₁为相邻四个象素的物理量(这里为速度)。假设φ_x和φ_y分别为φ在x方向和y方向的偏导数，则φ_x00，φ_x01，φ_x10，φ_x11，φ_y00，φ_y01，φ_y10，φ_y11为相邻四个象素的φ的偏导值。该USCIP多项式的系数由以下方程确定：

C₀₀＝φ₀₀

C₁₀＝φ_x00

C₀₁＝φ_y00

C₂₀＝3(φ₁₀-φ₀₀)-φ_x10-2φ_x00

C₀₂＝3(φ₀₁-φ₀₀)-φ_y01-2φ_y00

C₃₀＝-2(φ₁₀-φ₀₀)+φ_x10+φ_x00

C₀₃＝-2(φ₀₁-φ₀₀)+φ_y01+φ_y00

C₂₁＝3φ₁₁-2φ_x01-φ_x11-3(C₀₀+C₀₁+C₀₂+C₀₃)-C₂₀

C₃₁＝-2φ₁₁+φ_x01-φ_x11-2(C₀₀+C₀₁+C₀₂+C₀₃)-C₃₀

C₁₂＝3φ₁₁-2φ_y10-φ_y11-3(C₀₀+C₁₀+C₂₀+C₃₀)-C₀₂

C₁₃＝-2φ₁₁+φ_y10-φ_y11-2(C₀₀+C₁₀+C₂₀+C₃₀)-C₀₃

C₁₁＝φ_x01-C₁₀-C₁₂-C₁₃

从公式中可以看出，USCIP方法很容易在GPU上并行实现。我们首先计算出在各象素的两个偏导数φ_x和φ_y然后带入公式即可求出新的φ值。为了保证水平对流计算的稳定性，我们对于插值后的值限定在相邻四个象素的φ值的最大和最小值之间。

2)求解二维流体的动力学Navie-Stokes方程的粘性扩散项；

将Navie-Stokes方程的粘性扩散项

转换为隐式表达式，方程如(3)式：

$(I-\mathrm{v\δ}{t\▿}^2)u(x,t+\mathrm{\δt})=u(x,t)...\left(3\right)$

式中，I为单位矩阵，v为粘度，u为速度，t代表某一时刻，δt表示时间步长；

方程(3)可以看作是泊松方程，利用多栅解法在GPU上求解方程得到整个二维运算域中图层的粘性扩散速度场u₂；

多栅解法先在原始栅格上求初始值，然后将每个栅格的余数投影到粗糙一级的栅格继续求初始值和余数直到达到第N级栅格，最后逐次向高分辨率级别插值，得到最终解，具体包括以下步骤，如图2所示：

a)先在原始栅格上进行数次预处理求初值，获得一个初期解x^h；

b)求出各个栅格单元的余数值r^h＝b^h-A^hx^h；

c)将r^h余数值投射到粗糙一级的栅格上，成为b^2h，该栅格对应的泊松方程为A^2hx＝b^2h，其分辨率只有原始栅格的1/4大小；

d)重复1，2，3的步骤，反复递归直到某个最粗糙的栅格上(假设为第n个栅格)；此时，对当前栅格运动雅各比迭代方法解方程A^nhx＝b^nh，得到近似解x^nh；

e)将x^nh按于前面相反的顺序，逐次向高分辨率级别插值，并进行后期处理，当达到最高级别时就是最终解。

3)求解二维流体的动力学Navie-Stokes方程的外力项，将外力场直接叠加到速度场上，方法如下：

记录用户在搅动液体时经过路径的切线方向，并将所述的切线方向作为外力F的方向，外力F的大小由用户通过仿真系统交互设置，然后直接将外力场叠加到步骤1)和步骤2)的总速度场(u₁+u₂)上得到新的速度场w；

所述步骤1)和步骤2)的总速度场为水平对流速度场u₁和粘性扩散速度场u₂的叠加，即(u₁+u₂)；

4)求解二维流体的动力学Navie-Stokes方程的压力项，将步骤3)得到的速度场w减去压力的梯度得到新的速度场u的方法如下：

经过水平对流，粘性扩散，施加外力三个步骤后得到的新的速度场w是发散的，要将w转化为一个无发散的速度场u；

由于任何一个矢量场都可以被分解成另外两个矢量场的和，其中一个为无发散的矢量场，另一个为标量场的梯度，w分解如(4)式所示：

$w=u+\▿p...\left(4\right)$

式中p是Navie-Stokes方程中的压力场，对公式(4)两边应用发散度算子后根据Navie-Stokes方程(1)中的 $\▿\·u=0$ 的条件得到压力泊松方程(5)：

${\▿}^{2}p=\▿\·w...\left(5\right)$

采用多栅解法来解这个压力泊松方程，最后用w减去p的梯度就得到二维运算域中大理石纹釉每层图层的速度场u。

(4)利用USCIP高阶插值方法，通过转换函数以及逆转换函数对速度场进行转换计算，得到大理石纹釉每层图层的密度场，具体包括以下步骤：

通过求解Navie-Stokes方程后，得到了速度场u，但要将仿真结果转换成为可视的大理石纹釉图案，最为直接的途径就是加入密度场以表示置于流场中的彩色颜料。对颜料运动的精确模拟必须包括色彩的水平对流，扩散，损耗等性质。由于大理石纹釉工艺所用的基底液具有非常高的粘性，相对的颜料色彩扩散率就比较低，基本可以忽略不计，所以只考虑水平对流项。密度场u和速度场d之间的关系如公式(6)所示：

$\frac{\∂d}{\∂t}=-(u\·\▿)d...\left(6\right)$

式中d为密度场，u为速度场；

为了计算x位置在t+δt时刻的密度，先将t时刻密度场用如式(7)所示的转换函数进行转换，将d转换为F(d)：

$F\left(d\right)=0.5\log \left(\frac{\max (d,\mathrm{\ξ})}{\max (1-d,\mathrm{\ξ})}\right)...\left(7\right)$

式中d表示密度，ξ是一个很小的值，目的是为了防止出现分母为0的情况；将转换后的F(d)代入到(6)式，再采用准拉格朗日方法沿着粒子的运动轨迹回溯到在t时刻的位置x′，对最接近x′的四个单元格用USCIP方法进行插值得到的密度值拷贝到x位置，得到新的F(d)，最后将新的F(d)再通过式(8)所示的逆转换函数转换得到最终的密度场d。

F^-1(d)＝0.5(tanh(d)+1)......(8)

(5)将大理石纹釉每层图层的密度场和该层的颜色相乘得到该层的颜色值，并将大理石纹釉所有层的颜色值累加得到最后大理石纹釉图像的颜色值，生成用户需要的大理石纹釉；

大理石纹釉图像的颜色值按照式(9)进行计算：

$C_{\mathrm{output}}=C_{\mathrm{liguid}}*(1-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{d}_i)+\underset{i=0}{\overset{n}{\mathrm{\Σ}}}(C_i*d_i)...\left(9\right)$

式中，C_output是输出结果，C_liguid表示基底液的颜色，d_i表示第i层的密度值，C_i表示第i层的颜色值。

Claims (4)

1.一种快速且精确的大理石纹釉生成方法，其特征在于，包括以下步骤：

(1)用户设置用于大理石纹釉生成的二维运算域和基底液的属性和大理石纹釉的图层的层数n、图层的颜色color[n]；

输入液滴的半径R和搅动流体时的外力大小F两个参数的数值；

(2)采用用户指定的操作方式在二维运算域中进行相应的操作，所述的操作包括滴入各种颜色的滴液，搅动基底液带动各种颜色的图层运动和梳理图层的纹理；

(3)基于步骤(2)的操作，通过求解二维流体的动力学Navie-Stokes方程得到二维运算域中大理石纹釉每层图层的速度场；

(4)利用USCIP高阶插值方法，通过转换函数以及逆转换函数对速度场进行转换计算，得到大理石纹釉每层图层的密度场；

(5)将大理石纹釉每层图层的密度场和该层的颜色相乘得到该层的颜色值，并将大理石纹釉所有层的颜色值累加得到最后大理石纹釉图像的颜色值，生成用户需要的大理石纹釉。

2.如权利要求1所述的大理石纹釉生成方法，其特征在于：步骤(1)中所述的二维运算域是一个长为M个单位长度，宽为N个单位宽度的方形栅格上，栅格内部由M×N个正方形单元格构成，每个单元格的大小由用户在初始的时候设定；

基底液的属性包括颜色以及基底液流体的粘稠度。

3.如权利要求1所述的大理石纹釉生成方法，其特征在于：所述的通过求解二维流体的动力学Navie-Stokes方程得到二维运算域中大理石纹釉每层图层的速度场u包括以下步骤：

所述的二维流体的动力学Navie-Stokes方程如(1)式：

$\frac{\∂u}{\∂t}=-(u\·\▿)u-\frac{1}{\mathrm{\ρ}}\▿p+v{\▿}^{2}u+F$

$\▿\·u=0...\left(1\right)$

式中u表示流体的速度，ρ表示流体的密度，v为粘度，F是作用于流体的外力，式子右边的四个项都是加速度，

是水平对流项，是压力项，

是粘性扩散项，F是外部作用力项；

1)求解二维流体的动力学Navie-Stokes方程的水平对流项；

将Navie-Stokes方程的水平对流项

表示成隐式表达方式，方程如(2)式：

u(x，t+δt)＝u(x-u(x，t)δt，t)......(2)

式中u表示速度，x表示位置，t代表某一时刻，δt表示时间步长；

为了计算二维运算域中的x位置在t+δt时刻的速度，采用准拉格朗日方法沿着粒子的运动轨迹回溯到在t时刻的位置x′，对最接近x′的四个单元格用USCIP方法进行插值得到的速度拷贝到x位置，最后得到整个二维运算域中图层的水平对流速度场u₁；

2)求解二维流体的动力学Navie-Stokes方程的粘性扩散项；

将Navie-Stokes方程的粘性扩散项

转换为隐式表达式，方程如(3)式：

$(I-{\mathrm{v\δt}\▿}^2)u(x,t+\mathrm{\δt})=u(x,t)...\left(3\right)$

式中，I为单位矩阵，v为粘度，u为速度，t代表某一时刻，δt表示时间步长；

方程(3)可以看作是泊松方程，利用多栅解法在GPU上求解方程得到整个二维运算域中图层的粘性扩散速度场u₂；

3)求解二维流体的动力学Navie-Stokes方程的外力项，将外力场直接叠加到速度场上，方法如下：

记录用户在搅动液体时经过路径的切线方向，并将所述的切线方向作为外力F的方向，外力F的大小由用户交互设置，然后将外力场直接叠加到步骤1)和步骤2)的总速度场(u₁+u₂)上得到新的速度场w；

所述步骤1)和步骤2)的总速度场为水平对流速度场u₁和粘性扩散速度场u₂的叠加，即(u₁+u₂)；

4)求解二维流体的动力学Navie-Stokes方程的压力项，将步骤3)得到的速度场w减去压力的梯度得到新的速度场u的方法如下：

经过水平对流，粘性扩散，施加外力三个步骤后得到的新的速度场w是发散的，要将w转化为一个无发散的速度场u；

由于任何一个矢量场都可以被分解成另外两个矢量场的和，其中一个为无发散的矢量场，另一个为标量场的梯度，w分解如(4)式所示：

$w=u+\▿p...\left(4\right)$

式中p是Navie-Stokes方程中的压力场，对公式(4)两边应用发散度算子后根据Navie-Stokes方程(1)中的 $\▿\·u=0$ 的条件得到压力泊松方程(5)：

${\▿}^{2}p=\▿\·w...\left(5\right)$

采用多栅解法来解这个压力泊松方程，最后用w减去p的梯度就得到二维运算域中大理石纹釉每层图层的速度场u。

4.如权利要求1所述的大理石纹釉生成方法，其特征在于：步骤(4)所述的利用USCIP高阶插值方法，通过转换函数以及逆转换函数对速度场进行转换计算，得到大理石纹釉每层图层的密度场的方法包括以下步骤：

密度场和速度场之间的关系如公式(6)所示：

$\frac{\∂d}{\∂t}=-(u\·\▿)d...\left(6\right)$

式中d为密度场，u为速度场；

为了计算x位置在t+δt时刻的密度，先将t时刻密度场用如式(7)所示的转换函数进行转换，将d转换为F(d)：

$F\left(d\right)=0.5\log \left(\frac{\max (d,\mathrm{\ξ})}{\max (1-d,\mathrm{\ξ})}\right)...\left(7\right)$

式中d表示密度，ξ是一个很小的值，目的是为了防止出现分母为0的情况；将转换后的F(d)代入到(6)式，再采用准拉格朗日方法沿着粒子的运动轨迹回溯到在t时刻的位置x′，对最接近x′的四个单元格用USCIP方法进行插值得到的密度值拷贝到x位置，得到新的F(d)，最后将新的F(d)再通过式(8)所示的逆转换函数转换得到最终的密度场d。

F^-1(d)＝0.5(tanh(d)+1)......(8)

Images