CN107633123B - 一种用于光滑粒子流体动力学模拟出血及处理加速的方法 - Google Patents

Abstract

一种用于光滑粒子流体动力学模拟出血及处理加速的方法，利用基于GPU加速的模拟方法来模拟虚拟手术中的出血及利用虹吸原理吸出渗血的效果；同时，利用网格法实时划分问题区域，创建以支持域为边长的空间网格，通过临近网格搜索最近相邻粒子，并且通过并行计算架构(CUDA)多线程并行加速技术完成粒子控制方程的求解以及血液与固体交互的计算，大大提高了运算效率，从而提高了手术训练的实时性。此外；改进的移动立方体算法(marching cube)被用于流体表面的渲染，大大提高了手术训练的真实性。

Description

一种用于光滑粒子流体动力学模拟出血及处理加速的方法

技术领域

本发明涉及计算机图形学下的医学成像、流体动力学领域，尤其是一种基于 GPU加速的流体生成及在需要擦除时利用虹吸原理对渗血吸除的模拟分析方法。

背景技术

随着计算机硬件的发展，医学成像技术的进步，虚拟现实在医学领域中占有越来越大的影响力，拥有先进的虚拟现实技术是这个时代不可或缺的标志。虚拟现实技术是一种通过计算机图形学等学科高度仿真模拟现实的技术，在医学领域中，目前正被广泛应用于医学仿真训练的研究。虚拟现实技术的模拟真实度主要取决于模型的仿真效果以及运行的流畅度。虚拟训练系统中的对血液的模拟在传统上采用光滑粒子流体动力学(SPH)方法。光滑粒子流体动力学将流体离散成带有属性和体积的粒子，再结合有限元的思想为每个粒子属性求解流体力学方程得到下一时刻的状态，粒子的属性不是单独存在的，是通过对周围光滑核函数半径范围内的粒子属性值使用核函数加权求和得到。然而，传统血液模拟的质量一直以来都不是很理想，且由于算法计算量偏大使得成像的反应时间较长，导致模拟的效果难以满足实时性的要求。

发明内容

为了解决现有系统中血液的计算量大而导致的实时性和仿真效果的不理想的问题，本发明使用CUDA指令集来作为GPU加速的技术来提高计算的速度。同时采用一种改进的移动立方体方法对血液表面进行实时可视化渲染，增强了血液模拟的真实性。

为了解决上述技术问题提供了如下技术方案：

一种用于对光滑粒子流体动力学血液模拟过程加速的方法，所述方法包括以下步骤：

步骤一，使用PCISPH方法对粒子建模，过程如下：

1.1使用平坦和归一化内核得到的二阶精度插值设计核函数；

1.2由核函数计算出单个粒子所受的压力、粘性力、外部作用力，再根据牛顿第二定律计算出合力；

1.3通过各粒子的内部和外部的受力分析对出血做出模拟；

1.4根据对当前压力计算下一步速度和位置；

步骤二，基于CUDA对出血模拟的加速；

步骤三，对流体表面渲染，过程如下：

每进行一步的SPH模拟都要将点云转化为可着色的表面，将密度场均匀划分为大小相等的格子，使用Marching Cubes算法对格子的胞元的每个角点的密度场交点进行评估，根据其是否大于某个常量阈值来决定生成三角形的形式；使用 PN三角形方法来平滑三角形网格边缘并为顶点着色操作提供更多的采样点，根据三角形的三个顶点及其法线就生成平滑的细分三角形。

进一步，在将已经生成的血液进行清除操作的模拟时，需要将吸引器管口附近的血液粒子按照距离管口距离大小施加大小不同指向管口的吸引力，使得粒子以不同的速度向管口运动，且在运动到管口处与管口相碰撞后使该粒子“消失”。由于该过程中设计到判断有效范围，并在有效范围内分多个区域，每个区域会有不同的力施加给粒子等相关运算，计算量非常大，因此也采用CUDA进行加速运算；所述方法还包括以下步骤：

步骤四确定粒子系统中在该吸管吸引范围中的粒子，为吸引范围内不同的区域设置不同的吸引力；

以吸引器的管口为中心，根据设定的吸管吸引力的大小确定吸引半径，做出的吸引范围，在该范围内吸管才能做出有效的吸引，超出该范围则是无效吸引；按照吸管吸引的反方向为法线设置不同的角度设置不同的吸引力大小；

步骤五将吸到吸引器管口位置的粒子“吸掉”，对吸血过程基于CUDA加速

将运行到吸引器管口位置并与该管口相碰撞的粒子的生命周期设为0，删除该粒子，注销该粒子的内存空间；由于吸引范围内不同区域内吸引力大小不同使得粒子在不同的区域会有不同的指向吸引器管口的加速度。

再进一步，所述CUDA的计算包含以下参数：

a.哈密顿算子▽：对粒子的受力情况分析以及粒子速度加速度的计算时需对其所在的矢量场进行分析，因此在设计核函数的时候得计算出矢量场，即粒子所在位置的“梯度”，用“梯度”来表示标量场中粒子的变化快慢和方向；

在计算加速度时还需使用拉普拉斯算子，即▽·▽或▽²；

b.Helmholtz-Hodge定理：任意的矢量场都可以唯一的分解成一个零散度的矢量场和一个标量场的梯度之和，如下所示：

w＝u+▽p

其中u是零散度矢量场，p是标量场。对上述公式求散度得：

▽·w＝▽·u+▽²p

因为u的散度为零，所以得压强泊松方程：

▽²p＝▽·w

c.扩散项的CUDA映射内核函数对计算部分加速迭代：该内核函数对应的 Grid线程结构被设计为二维成结构，其在x维度上的block数量为：

(M+blockDim.x-1)/blockDim.x

在y维度上的block数量为：

(M+blockDim.y-1)/blockDim.y

d.对流项的计算：本发明采用的是Jos Stam提出的隐式方法求解对流项，即从当前时刻起，根据粒子的运动轨迹，计算出上一时刻粒子的位置，并将该位置粒子的q量赋值给现在的这个粒子，其计算公式为:

q(x,t+δt)＝q(x-u(x,t)δt,t)

其中，q(x,t+δt)是t+δt时刻粒子的位置，u(x,t)表示x位置处t时刻的速度；

e.Navie-Stokes方程组：用于描述流体的运动，当描述牛顿流体的不可压缩流体时，方程组定义的方程可被简化为不可压缩的质量守恒方程和动量守恒方程：

其中ρ是密度，p是压强，v是速度，μ是粘性系数，f^ext是所有外力之和；

f.邻近区域内粒子插值计算属性：流体的模拟是一种基于拉格朗日方法的流体建模方法，将流体离散为大量的粒子，是一种基于粒子的插值方法，在位置 ri处的粒子i的一个标量A(r_i)由该粒子邻近区域内粒子的属性通过插值得到：

其中，n是领域内粒子的数量，j是粒子的索引号，A_j是粒子j的属性值，m_j是粒子j的质量，ρ_j是粒子j的密度，r_i是粒子j的位置，W(r_i-r_j,h) 是插值权重函数在SPH方法中称为光滑核函数，h是光滑核函数W影响区域的半径，在支持域以外光滑核函数W被定义为零；依据该插值计算可依次求得标量A(r_i)的梯度▽A(r_i)、拉普拉新算子▽²A(r_i)，进而求得粒子的各属性值；

h.圆管流雷诺数：

表针粘性影响的相似准数，表征流体流动情况的无量钢数，雷诺数大小决定了粘性流体的流动特性，u为流速，v为运动粘度， d为吸管直径。

更进一步，所述CUDA是由CPU和GPU协同完成的，CPU作为主机端，主要处理逻辑性较强的事务及串行程序的计算，GPU扮演设备端通过kernel函数并行处理负责密集型大数据的并行运算；所需要的出血模拟是一个由粒子系统模拟的模型，在该模拟中将大量的重复的运算放置到GPU端运算处理，CUDA的每一个线程都有一个唯一的线程ID，所以CUDA线程可使用一维、二维、三维索引标识，从而形成一维、二维、三维线程块，线程块包含在线程格(grid)中，线程格的尺寸根据粒子系统的规模来指定。

所述步骤二中，基于CUDA对出血模拟加速的过程如下：

2.1初始化数据，获取流体的外力大小，创建密度场和速度场的存储空间及其相关的初始化；

2.2数据传输，将相关粒子属性数据拷贝到显存中，为GPU上流体的计算提供数据；

2.3在GPU上实现模拟，包括了扩散项、外力项、对流项、投影过程中各个步骤的计算，从而得到每一个时间步长后流体的速度场和密度场，计算出下一个时间步长后的速度和位置；

2.4将GPU计算的速度场和密度场渲染到屏幕的窗口中去。

2.5边界的处理

速度在边界处遵循无滑条件，且密度保持不变，在水平方向上的流体边界条件为：

水平边界1：

Dens(i,0)＝Dens(i,1)

U(i,0)＝U(i,1)

V(i,0)＝-V(i,1)

水平边界2：

Dens(i,M+1)＝Dens(i,M)

U(i,M+1)＝U(i,M)

V(i,M+1)＝-V(i,M)

在垂直方向上的流体边界条件为：

垂直边界3：

Dens(0,j)＝Dens(1,j)

U(0,j)＝U(1,j)

V(0,j)＝-V(1,j)

垂直边界4：

Dens(N+1,j)＝Dens(N,j)

U(N+1,j)＝U(N,j)

V(N+1,j)＝-V(N,j)

最后，在四个边角单元处的边界条件为:

左下边角 x(0,0)＝0.5f*(x(1,0)+x(0,1))

右下边角 x(M+1,0)＝0.5f*(x(M,0)+x(M+1,1))

左上边角 x(0,N+1)＝0.5f*(x(1,N+1)+x(0,N))

右上边角 x(M+1,N+1)＝0.5f*(x(M,N+1)+x(M+1,N))

其中，U代表水平方向的速度，V代表垂直方向的速度，Dens代表密度，x 代表流体的速度或密度，1<＝i<＝N,1<＝j<＝M，根据给出的流体边界处的条件，在设备端定义了一个set_bnd_kernerl的内核函数，通过该内核函数在GPU上的并行执行来加速边界处流体的模拟过程；

在定义set_bnd_kernerl函数时，首先计算水平边界的流体，其次计算垂直边界的流体，最后计算四个边角网格单元的流体模拟，该内核函数对应的线程组织结构被设计成一维的形式；

根据GPU内核函数的创建方法，先划分所要处理的任务，并确定执行该任务的线程数量及其Grid结构，将每个网格单元处流体的外力和源的计算交由一个线程来处理，且将这些线程构成的Grid网格定义成一个一维的形式，共包含的block块为：

((M+2)*(N+2)+blockDim.x-1)/blockDim.x

每个block块也设计成一维的形式，且每个block块包含256个线程。

所述步骤三中，对流体表面渲染包括以下步骤：

3.1对于每个粒子i计算其位置均值

将结果与上一时刻的位置均值做线性平均；

3.2根据

及其领域的粒子位置进行加权主成分分析，得到代表粒子领域分布形状的矩阵C_i，并根据C_i及其特征值计算放缩标量S_i；

3.3对于三维密度场中每个体素，搜索其领域的粒子，根据体素中心到粒子距离进行插值得到密度，插值时将距离乘以S_i进行调整；

3.4对三维密度场进行立方体匹配，得到三角形网格；

3.5对三角形网格进行表面细分，生成PN三角形，得到光滑的网格。

本发明的技术构思为：假设模拟的血液为均匀密度单相流，继承传统模拟方法再加以改进，并通过GPU加速的手段完成一套高效率血液模拟的实现。

使用基于GPU加速的血液模拟具有很重要的实际意义。光滑粒子流体动力学将流体离散成带有属性(如压强、速度、加速度等)和体积的粒子，结合有限元的思想为每个粒子求解方程得到该粒子下一时刻的状态，粒子的属性是通过对其周围光滑核函数半径内所有粒子的属性使用核函数加权求和得到。解决了流体的不可压性和均匀密度单相流的不可压问题。由于模拟过程中计算量大，使用 GPU加速后将大为减轻CPU的负担。这些研究使虚拟手术系统的仿真度有了极大的提高，使得操作人员的操作训练更为逼真有效。

本发明的有益效果为：使用CUDA指令集来作为GPU加速的技术来提高计算速度，大大提高了血液生成和清除的运算，减轻CPU的运算量；同时采用一种改进的移动立方体方法对血液表面进行实时渲染，增强了手术训练模拟的真实性。

附图说明

图1为本发明中的出血模拟系统框图；

图2为本发明中的吸血模拟系统框图。

具体实施方式

下面结合附图对本发明作进一步说明。

参照图1和图2，一种用于对光滑粒子流体动力学模拟的出血和血液清除全过程加速的方法，包括以下步骤：

对于出血的模拟：

步骤一，使用PCISPH方法对粒子建模，过程如下：

1.1出血模拟的高度仿真依赖于对出血系统中粒子的受力情况，所以需对每个粒子的受力分析，单个粒子的属性都会“扩散”到周围，并且随着距离的减小逐渐减小。我们使用平坦和归一化的内核因为他们能够得到二阶精度的插值。我们设计了以下核函数：

1.2单个粒子所受到的各力分别为，其中外力一般就是重力，流体的质量是由流体单元的密度决定的，所以一般用密度代替质量：

粘性力：

压力：

1.3出血模拟大多采用光滑粒子流体动力学，这些粒子相互影响，共同形成了复杂的流体运动，对于每个单独的流体微粒，依旧遵循最基本的牛顿第二定律。对于空间中任意给定一点

的任意物理量A而言，使用SPH方法可以通过如下离散求和计算：

其中，

代表粒子j的位置，m_j是指粒子j的质量，ρ_j表示密度，

表示半径为h的光滑核函数。通光滑滑核函数的导数的加权和得到物理量的空间导数，即

以及：

在欧拉公式中，等温流体由速度v、密度ρ和压力p描述。这些量随时间的变化由以下两个公式表示：

其中g表示重力加速度，μ表示液体的粘度。在SPH方法中，粒子的质量和数量是恒定的，所以外部力即重力可以省略且合力的计算公式可转化为：

f＝＝-▽p+ρg+μ▽²v＝f_p+f_g+f_v

根据以上的规则，可以得到各个物理量的计算公式，如对于粒子i，它的密度可以通过领域的加权和计算：

1.4根据当前压力计算下一步速度和位置。

使用临近粒子搜索法遍历所有粒子以提高PCISPH的效率。首先我们在问题区域构建网格，为了保证构建的速度，先搜索靠边界最近的粒子，即在空间坐标 x,y,z三个方向上找出最大和最小位置，以这些位置作为血液表面边界参考点构造网格。这种构建网格的方法不仅缩小了空间域中的问题区域，减小构建网格的数量，提高运算速度，并且所有粒子都被遍历但不影响模拟效果。再计算出某一个粒子的相近粒子。在计算时只需要遍历该粒子所在网格和相邻的网格共27个网格内的粒子即可。在建立的空间索引中根据粒子的加速度和时间步长得到下一刻的速度和位置。

步骤二，基于CUDA对出血模拟的加速

为了保证模拟效果的灵活性和真实性，粒子的数量会很大，这样增加了运算的负担，容易造成实时性的伤害。本发明将利用CPU和GPU的并行运算分担运算来提高计算效率，以提高出血模拟以及虚拟手术训练系统运行的实时性。使用 CUDA架构将在PCISPH算法中对每个粒子有相同的操作集合，包括对每个粒子的密度、压强及内力和外力的计算。我们将所有可以并行运算的步骤在CUDA 上实现。

基于CUDA加速的流体模拟计算过程如下：

2.1初始化数据，获取流体的外力大小，创建密度场和速度场的存储空间及其相关的初始化，在没有人为外力加持的。

2.2数据传输，将相关粒子属性数据拷贝到显存中，为GPU上流体的计算提供数据。

2.3在GPU上实现模拟，包括了扩散项、外力项、对流项、投影过程中各个步骤的计算。从而得到每一个时间步长后流体的速度场和密度场，计算出下一个时间步长后的速度和位置。

2.4将GPU计算的速度场和密度场渲染到屏幕的窗口中去。

2.5边界的处理：

因为流体在边界区域与非边界区域的运动状态存在着一定的差异，所以需要单独处理流体边界区域的模拟。本发明规定速度在边界处遵循无滑条件，且密度保持不变。在水平方向上的流体边界条件为(分别将相对面的x轴方向上的左右边界设为边界1和边界2):

水平边界1：

Dens(i,0)＝Dens(i,1)

U(i,0)＝U(i,1)

V(i,0)＝-V(i,1)

水平边界2：

Dens(i,M+1)＝Dens(i,M)

U(i,M+1)＝U(i,M)

V(i,M+1)＝-V(i,M)

在垂直方向上的流体边界条件为(分别将相对面y轴方向上下的边界设为边界3和边界4):

垂直边界3：

Dens(0,j)＝Dens(1,j)

U(0,j)＝U(1,j)

V(0,j)＝-V(1,j)

垂直边界4：

Dens(N+1,j)＝Dens(N,j)

U(N+1,j)＝U(N,j)

V(N+1,j)＝-V(N,j)

最后，在四个边角单元处的边界条件为:

左下边角 x(0,0)＝0.5f*(x(1,0)+x(0,1))

右下边角 x(M+1,0)＝0.5f*(x(M,0)+x(M+1,1))

左上边角 x(0,N+1)＝0.5f*(x(1,N+1)+x(0,N))

右上边角 x(M+1,N+1)＝0.5f*(x(M,N+1)+x(M+1,N))

其中，U代表水平方向的速度，V代表垂直方向的速度，Dens代表密度，x 代表流体的速度或密度，1<＝i<＝N,1<＝j<＝M。根据给出的流体边界处的条件，本发明在设备端定义了一个set_bnd_kernerl的内核函数，通过该内核函数在GPU 上的并行执行来加速边界处流体的模拟过程。

在定义set_bnd_kernerl函数时，首先计算水平边界的流体，其次计算垂直边界的流体，最后计算四个边角网格单元的流体模拟。该内核函数对应的线程组织结构被设计成一维的形式，

根据GPU内核函数的创建方法，本发明先划分所要处理的任务，并确定执行该任务的线程数量及其Grid结构。本发明将每个网格单元处流体的外力和源的计算交由一个线程来处理，且将这些线程构成的Grid网格定义成一个一维的形式，共包含的block块为：

((M+2)*(N+2)+blockDim.x-1)/blockDim.x

每个block块也设计成一维的形式，且每个block块包含256个线程。

步骤三，对流体表面渲染处理

为了得到理想的渲染结果，每个时间步长的PCISPH模拟都要将点云转化为可着色的表面。因为上一步对边界的处理时垂直方向上的表面渲染是存在连接不够平滑而导致模拟的效果不够逼真，故此需要对表面的渲染做出进一步的优化渲染操作，过程如下：

3.1对于每个粒子i计算其位置均值

为了增强时域连续性，将结果与上一时刻的位置均值做线性平均(此方法称为XSPH[16])。

3.2根据

及其领域的粒子位置进行加权主成分分析，得到代表粒子领域分布形状的矩阵C_i，并根据C_i及其特征值计算放缩标量S_i。

3.3对于三维密度场中每个体素，搜索其领域的粒子，根据体素中心到粒子距离进行插值得到密度，插值时将距离乘以S_i进行调整。

3.4对三维密度场进行立方体匹配，得到三角形网格。

3.5对三角形网格进行表面细分，生成PN三角形，得到光滑的网格。

对吸血的模拟：

步骤四，确定有效的吸血范围，并为该范围内的不同区域分配不同的吸引力，过程如下

4.1以吸引器的管口为中心O，吸引器的反方向上截取线段OA，以OA的大小为半径划球，球内的范围就是有效范围。

4.2对新引力大小影响较大的因素有到管口的距离大小和与OA夹角的大小。在有效范围内，OB是与OA夹角为θ的半径锥面的素线。由于现实中的吸血的力比较小，所以按照距离吸口的距离把OA分为大小相等的三段，再按照该三段的大小将球划分为三个球层，由此三层将吸引力按照距离管口的远近分为三个大小不等的力。从吸口到球面的三个球层的吸引力分别为f，2/3f,1/3f。θ的范围是[0,90°]时，角度对力的影响较小，故选择忽略，主要考略距离的因素。当θ的范围是[90°,180°]时，三个球层的力分别是f*sinθ,2/3f*sinθ,1/3f*sinθ。

步骤五，将血液粒子吸引到管口附近并吸掉，对吸血过程基于CUDA加速，过程如下：

5.1吸管的吸引力远大于血液内部的力，由牛顿第二定律得出在有效范围内血液将按照不同的区域不同的加速度向管口运动，逼近现实现实中的流体运动的效果。

5.2基于CUDA加速的过程类似于出血的模拟 过程中的加速过程。将相关数据拷贝到显存中，计算在吸引力和液体的内部力下粒子的运动速度，加以模拟。

Claims (4)

1.一种用于光滑粒子流体动力学模拟出血及处理加速的方法，其特征在于：所述方法包括以下步骤：

步骤一，使用PCISPH方法对粒子建模，过程如下：

1.1使用平坦和归一化的内核设计核函数；

1.2由核函数计算出单个粒子所受的压力、粘性力、外部作用力，再根据牛顿第二定律计算出合力；

1.3通过各粒子的内部和外部的受力分析对出血做出模拟；

1.4根据当前压力计算下一步速度和位置；

步骤二，基于CUDA对出血模拟的加速；

步骤三，对流体表面渲染，过程如下：

每进行一步的SPH模拟都要将点云转化为可着色的表面，将密度场均匀划分为大小相等的格子，使用Marching Cubes算法对格子的胞元的每个角点的密度场交点进行评估，根据其是否大于某个常量阈值来决定生成三角形的形式；使用PN三角形方法来平滑三角形网格边缘并为顶点着色操作提供更多的采样点，根据三角形的三个顶点及其法线生成平滑的细分三角形；

所述CUDA是由CPU和GPU协同完成的，CPU作为主机端，主要处理逻辑性较强的事务及串行程序的计算，GPU扮演设备端通过kernel函数并行处理负责密集型大数据的并行运算；所需要的出血模拟是一个由粒子系统模拟的模型，在该模拟中将大量的重复的运算放置到GPU端运算处理，CUDA的每一个线程都有一个唯一的线程ID，所以CUDA线程可使用一维、二维、三维索引标识，从而形成一维、二维、三维线程块，线程块包含在线程格中，线程格的尺寸根据粒子系统的规模来指定；

所述步骤二中，基于CUDA对出血模拟的加速的过程如下：

2.1初始化数据，获取流体的外力大小，创建密度场和速度场的存储空间及其相关的初始化；

2.2数据传输，将相关粒子属性数据拷贝到显存中，为GPU上流体的计算提供数据；

2.3在GPU上实现模拟，包括了扩散项、外力项、对流项、投影过程中各个步骤的计算，从而得到每一个时间步长后流体的速度场和密度场，计算出下一个时间步长后的速度和位置；

2.4将GPU计算的速度场和密度场渲染到屏幕的窗口中去；

2.5边界的处理

速度在边界处遵循无滑条件，且密度保持不变，在水平方向上的流体边界条件为：

水平边界1：

Dens(i,0)＝Dens(i,1)

U(i,0)＝U(i,1)

V(i,0)＝-V(i,1)

水平边界2：

Dens(i,M+1)＝Dens(i,M)

U(i,M+1)＝U(i,M)

V(i,M+1)＝-V(i,M)

在垂直方向上的流体边界条件为：

垂直边界3：

Dens(0,j)＝Dens(1,j)

U(0,j)＝U(1,j)

V(0,j)＝-V(1,j)

垂直边界4：

Dens(N+1,j)＝Dens(N,j)

U(N+1,j)＝U(N,j)

V(N+1,j)＝-V(N,j)

最后，在四个边角单元处的边界条件为:

左下边角x(0,0)＝0.5f*(x(1,0)+x(0,1))

右下边角x(M+1,0)＝0.5f*(x(M,0)+x(M+1,1))

左上边角x(0,N+1)＝0.5f*(x(1,N+1)+x(0,N))

右上边角x(M+1,N+1)＝0.5f*(x(M,N+1)+x(M+1,N))

其中，U代表水平方向的速度，V代表垂直方向的速度，Dens代表密度，x代表流体的速度或密度，1<＝i<＝N,1<＝j<＝M，根据给出的流体边界处的条件，在设备端定义了一个set_bnd_kernerl的内核函数，通过该内核函数在GPU上的并行执行来加速边界处流体的模拟过程；

在定义set_bnd_kernerl函数时，首先计算水平边界的流体，其次计算垂直边界的流体，最后计算四个边角网格单元的流体模拟，该内核函数对应的线程组织结构被设计成一维的形式；

根据GPU内核函数的创建方法，先划分所要处理的任务，并确定执行该任务的线程数量及其Grid结构，将每个网格单元处流体的外力和源的计算交由一个线程来处理，且将这些线程构成的Grid网格定义成一个一维的形式，共包含的block块为：

((M+2)*(N+2)+blockDim.x-1)/blockDim.x

每个block块也设计成一维的形式，且每个block块包含256个线程。

2.如权利要求1所述的一种用于光滑粒子流体动力学模拟出血及处理加速的方法，其特征在于：在对已经生成的血液进行清除操作的模拟时，将吸引器管口附近的血液粒子按照距离管口距离大小用不同的速度向管口运动，且在运动到管口使该粒子“消失”，计算过程也采用CUDA进行加速；所述方法还包括以下步骤：

步骤四，确定粒子系统中在该吸管吸引范围中的粒子，为吸引范围内不同的区域设置不同的吸引力；

以吸引器的管口为中心，根据设定的吸管吸引力的大小确定吸引半径，做出的吸引范围，在该范围内吸管才能做出有效的吸引，超出该范围则是无效吸引；按照吸管吸引的反方向为法线设置不同的角度和不同的吸引力大小；

步骤五，将吸到吸引器管口位置的粒子“吸掉”，对渗血吸除过程基于CUDA实现加速；

将运行到吸引器管口位置并与该管口相碰撞的粒子的生命周期设为0，删除该粒子，注销该粒子的内存空间；由于吸引范围内不同区域内吸引力大小不同使得粒子在不同的区域会有不同的加速度。

3.如权利要求2所述的一种用于光滑粒子流体动力学模拟出血及处理加速的方法，其特征在于：所述CUDA的计算包含以下参数：

a.哈密顿算子▽：对粒子的受力情况分析以及粒子速度加速度的计算时需对其所在的矢量场进行分析，因此在设计核函数的时候得计算出矢量场，即粒子所在位置的“梯度”，用“梯度”来表示标量场中粒子的变化快慢和方向；

在计算加速度时还需使用拉普拉斯算子，即▽·▽或▽²；

b.Helmholtz-Hodge定理：任意的矢量场都可以唯一的分解成一个零散度的矢量场和一个标量场的梯度之和，如下所示：

w＝u+▽p

其中u是零散度矢量场，p是标量场，对上述公式求散度得：

▽·w＝▽·u+▽²p

因为u的散度为零，所以得压强泊松方程：

▽²p＝▽·w

c.扩散项的CUDA映射内核函数对计算部分加速迭代：该内核函数对应的Grid线程结构被设计为二维结构，其在x维度上的block数量为：

(M+blockDim.x-1)/blockDim.x

在y维度上的block数量为：

(M+blockDim.y-1)/blockDim.y

d.对流项的计算：采用的是Jos Stam提出的隐式方法求解对流项，即从当前时刻起，根据粒子的运动轨迹，计算出上一时刻粒子的位置，并将该位置粒子的q量赋值给现在的这个粒子，其计算公式为:

q(x,t+δt)＝q(x-u(x,t)δt,t)

其中，q(x,t+δt)是t+δt时刻粒子的位置，u(x,t)表示x位置处t时刻的速度；

e.Navie-Stokes方程组：用于描述流体的运动，当描述牛顿流体的不可压缩流体时，方程组定义的方程可被简化为不可压缩的质量守恒方程和动量守恒方程：

其中ρ是密度，v是速度，μ是粘性系数，f^ext是所有外力之和；

f.邻近区域内粒子插值计算属性：流体的模拟是一种基于拉格朗日方法的流体建模方法，将流体离散为大量的粒子，是一种基于粒子的插值方法，在位置r_i处的粒子i的一个标量A(r_i)由该粒子邻近区域内粒子的属性通过插值得到：

其中，n是领域内粒子的数量，j是粒子的索引号，A_j是粒子j的属性值，m_j是粒子j的质量，ρ_j是粒子j的密度，r_j 是粒子j的位置，W(r_i-r_j,h)是插值权重函数在SPH方法中的光滑核函数，h是光滑核函数W影响区域的半径，在影响区域以外光滑核函数W被定义为零；依据该插值计算可依次求得标量A(r_i)的梯度▽A(r_i)、拉普拉斯算子▽²A(r_i)，进而求得粒子的各属性值；

h.圆管流雷诺数：

表征粘性影响的相似准数，相似准数为流体流动情况的无量纲数，雷诺数大小决定了粘性流体的流动特性，u为流速，γ为运动粘度，d为吸管直径。

4.如权利要求1所述的一种用于光滑粒子流体动力学模拟出血及处理加速的方法，其特征在于：所述步骤三中，对流体表面渲染包括以下步骤：

3.1对于每个粒子i计算其位置均值

将结果与上一时刻的位置均值做线性平均；

3.2根据

及其领域的粒子位置进行加权主成分分析，得到代表粒子领域分布形状的矩阵C_i，并根据C_i及其特征值计算放缩标量S_i；

3.3对于三维密度场中每个体素，搜索其领域的粒子，根据体素中心到粒子距离进行插值得到密度，插值时将距离乘以S_i进行调整；

3.4对三维密度场进行立方体匹配，得到三角形网格；

3.5对三角形网格进行表面细分，生成PN三角形，得到光滑的网格。

Images