[发明专利]一种快速且精确的大理石纹釉生成方法

权利要求书

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)式：

∂u∂t=-(u·▿)u-1ρ▿p+v▿2u+F]]>

▿·u=0...(1)]]>

式中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-vδt▿2)u(x,t+δt)=u(x,t)...(3)]]>

式中，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...(4)]]>

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

▿2p=▿·w...(5)]]>

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

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

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

∂d∂t=-(u·▿)d...(6)]]>

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

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

F(d)=0.5log(max(d,ξ)max(1-d,ξ))...(7)]]>

式中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)