一种基于Cross模型的多粒子混合修正的血栓模拟方法与流程

本发明涉及计算机图形学领域，具体涉及一种基于cross模型的多粒子混合修正的血栓模拟方法。

背景技术

血液是由心脏推动并在全身血管内循环流动的一种组织，具有运输、防御、调节体温和维持酸碱平衡及渗透压平衡的功能，对维持人类生命健康起着至关重要的作用。血栓是血液因为某种原因而凝固的产生的，由血栓引起的疾病可能会对人类生命健康造成巨大威胁。

流体，是气体和液体的总称，贯穿于人类生活的方方面面，如何在高效地在计算机上模拟出符合真实流体性质和运动规律的流体动画是近年来计算机图形领域一个十分热门的话题。计算机流体模拟的方法可以根据观察对象的不同来分成两大类，其中一类称为欧拉法，该类方法的主要思路是将模拟场景划分为均匀网格，将网格上的每个格点作为观察对象，通过某种数值方法来在这些格点上求解n-s方程，同时将流体的物理数值存放于格点上；另一类方法称为拉格朗日法，该类方法的主要思路是将所要模拟的流体划分为流体微团，并用流体粒子来编写流体微团，然后将每个流体粒子作为观察对象，通过某种数值方法在每个粒子上求解n-s方程，同时将流体的物理数值存放于粒子中。sph法就是一种应用广泛的拉格朗日流体模拟法。

sph法通过核函数加权法来求解n-s方程，即每个流体粒子上的n-s方程的物理项，包括密度、压力和粘性力等，由粒子固定半径内的邻域粒子的相对应的物理数值加权相加得到，其中的权重值由选取合适的核函数计算得到。进而计算出每个粒子的加速度，再选取合适的时间步长计算粒子的速度，从而可以计算出每个时刻粒子的位置。这样只要记录下每个时刻每个粒子的位置，就可以模拟出流体整体的运动情况。

自然界中的流体通常可以划分成两类，即牛顿流体和非牛顿流体。在温度不变的条件下，牛顿流体的粘性系数是常量数值，而非牛顿流体的粘性系数随着流体的切变率的改变而改变。血液是一种常见的剪切变稀非牛顿流体，即血液的粘度随着切应力的增大而减小。血液除了具有非牛顿流体常见的物理性质，还有许多独特的生化性质和流变性。

血栓的形成即是血液的一种独特的现象。血栓是一种血凝块，当血管损伤流血时，血小板、红细胞和纤维蛋白等相互作用而形成血块粘附于伤口处以阻止血液流失，但是当血栓形成于血管中，严重影响正常血液流动时，可能会引发各种疾病，对人类身体健康造成巨大危害。血栓的形成主要原因有血管内壁损伤、湍流等引起的血流异常和血液粘性过高等。正确地模拟血液中血栓的形成，在医学领域如辅助医学研究、模拟血管内的某些病理现象等具有重要的意义，甚至在影视、游戏也能有重要应用。

现有的非牛顿流体模拟的技术方案主要是使用基于cross模型的sph方法来模拟非牛顿流体。

sph法的主要思路如下：将流体离散化为流体粒子，流体的物理数值都由流体粒子携带；在每个粒子上基于核函数加权法求解n-s方程，即选取适当的核函数，对于每个粒子，找到其固定核半径内的邻域粒子，然后基于邻域粒子的相对应的物理数值加权求和，计算出该粒子的密度、所受压力和所受粘性力，再加上粒子所受重力，计算出粒子所受合力，从而计算出粒子加速度；进一步计算出粒子的速度和位置，最后由所有粒子每一时刻的位置来刻画出整个流体的运动。

cross模型是一种简化的剪切变稀非牛顿流体模型。在该模型下，在温度不变的情况下流体的粘性系数随着流体的形变率的增大而减小，而不是一个常量，再通过形变张量和粘性系数计算出流体的粘性力。因此，不同于牛顿流体的sph法的是，基于cross模型的非牛顿流体的sph法在计算完流体粒子的密度和压强后，需要先依次计算出流体粒子的形变张量、形变率，进而计算出粘性系数，接下来才可以计算流体粒子所受粘性力。

使用基于cross模型的非牛顿流体的sph法的模拟步骤如下：

1)创建流体粒子的运动空间，创建粒子并初始化粒子的各项物理数值。

2)在流体的运动空间上创建均匀网格，并记录粒子所在网格，以便于粒子固定半径内邻域搜索。

3)对于每个流体粒子，计算其固定半径内邻域粒子。选取合适的核函数，根据所有邻域粒子的密度加权求和计算出粒子的密度；然后根据理想气体状态方程计算粒子压强，若所求得的压强值为负数，则将该粒子的压强值设置为0。

4)对于每个流体粒子，计算其固定半径内邻域粒子。选取合适的核函数，求出粒子与其邻域粒子的速度差，根据所求得的速度差的梯度加权计算出粒子的形变张量，从而计算出形变率，再根据形变率计算出流体粒子的粘性系数，最后根据形变张量和粘性系数，计算出流体粒子的切应力。

5)对于每个流体粒子，计算其固定半径内邻域粒子。选取合适的核函数，根据所有邻域粒子的压强和自身压强计算出粒子所受压力；选取spline函数作为核函数，根据步骤3所求得的邻域粒子的切应力和自身的切引力，计算出粒子所受粘性力；再加上粒子所受重力，从而可以计算出粒子所受合力。

6)对于每个流体粒子，根据步骤5所求得的粒子所受的合力，计算出粒子的加速度，再根据预设的时间步长计算出粒子下一刻的速度，进而计算出粒子下一时刻的位置，最后再判断是否与固体边界发生碰撞检测，若是则修正粒子的速度和位置。

7)重复步骤2到6，直到模拟结束。

sph方法来模拟非牛顿流体具有两方面的缺陷：

1)模拟非牛顿流体时，通常只是模拟单一类型粒子的流体或是多种流体混合。然而血液是一种由血浆和各类血细胞等组合而成的复杂非牛顿流体，具有许多特殊的理化性质，血细胞本身的各种生物性质和行为也会使得血液产生其他非牛顿流体没有的特殊现象，血液既远远复杂于单一类型粒子的流体又不是简单的多种流体混合，因此现有的技术方案不能正确地模拟出血液。

2)模拟非牛顿流体时，不能模拟出流体的凝固现象。因此，仅仅使用sph法模拟血液，是无法模拟出血液中血栓的形成的。

技术实现要素：

本发明的目的是克服现有方法的不足，提出了一种基于cross模型的多粒子混合修正的血栓模拟方法。本发明解决两个核心技术问题：(1)使用一种基于cross模型的多粒子混合修正sph法来模拟血液，从而解决在模拟过程中出现的质量大的粒子往下沉、质量小的粒子往上浮的不符合正常血液生理现象的问题；(2)使用一种逐粒子聚合的方法以简化的血栓模型来模拟血液中血栓的形成，从而解决传统方法无法模拟出血液中血栓形成的问题。

为了解决上述问题，本发明提出了一种基于cross模型的多粒子混合修正的血栓模拟方法，所述方法包括：

初始化处理，创建血管模型均匀网格及血液粒子；

使用基于cross模型的多粒子混合修正sph法模拟血液；

在前述模拟血液的基础上，使用逐粒子聚合法模拟血液中血栓的形成。

优选地，所述初始化处理、创建血管模型均匀网格及血液粒子，具体包括：

初始化全局数据；

导入血管模型，作为血液粒子的运动空间；

在血管模型内壁的三角面片上采样足够密集的点，用来计算血液粒子与血管壁相碰撞时血管壁对粒子的压力；

在血管模型内创建一组血液粒子；

计算出血管模型的包围盒，并将其划分为均匀网格。

优选地，所述使用基于cross模型的多粒子混合修正sph法模拟血液，具体包括：

对于每个血液粒子，计算出其所在网格并排序其所携带数据；

对于每个血液粒子，搜索其固定半径内的邻域粒子，然后计算其密度和压强；

对于每个血液粒子，搜索其固定半径内的邻域粒子，计算其形变张量、形变率和粘性系数，从而计算出其切应力；

对于每个血液粒子，搜索其固定半径内的邻域粒子，计算其所受的压力和粘性力；

对于每个血液粒子，计算其加速度和下一时刻的速度，进而计算其下一时刻的位置，并更新其在场景中的位置。

优选地，所述逐粒子聚合法模拟血液中血栓的形成，具体为：

对于每个血液粒子，通过粒子速度或者粘性系数来判断是否血栓粒子；

对于血栓粒子，若其与邻域血栓粒子能组成新的血栓块，则计算新血栓块的体积、质量、密度；

计算新血栓块所受合力和加速度，并计算速度和位移，从而模拟出血栓块的运动。

本发明提出的基于cross模型的多粒子混合修正的血栓模拟方法能够较为真实地模拟血液，并且模拟血液中血栓的形成。正确模拟血液中血栓的形成，有利于对血栓的研究，可以应用于虚拟手术、医学辅助研究等医学应用；同时，本发明技术方案可以满足游戏与影视行业相关方面的需求。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1是本发明实施例的血栓模拟方法流程图；

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1是本发明实施例的血栓模拟方法流程图，如图1所示，该方法包括：

s1，初始化处理，创建血管模型均匀网格及血液粒子；

s2，使用基于cross模型的多粒子混合修正sph法模拟血液；

s3，在前述模拟血液的基础上，使用逐粒子聚合法模拟血液中血栓的形成。

步骤s1，具体如下：

s11，初始化全局数据：

硬编码或从配置文件中读取各种模拟数值，如模拟的时间步长、核半径、重力加速度、粘性系数的最大和最小值以及各类血液粒子的质量、静态密度和颜色等等。

s12，导入血管模型，作为血液粒子的运动空间：

血管模型由专业三维模型制作软件制作或ct数据导入并保存为obj格式文件，由模拟程序读入并提取需要使用的数据，包括血管模型的点、法向量、颜色和三角形面片顶点索引等，然后渲染并显示血管模型。

s13，在血管模型的三角面片上采样足够密集的点，用来计算血液粒子与血管壁相碰撞时血管壁对粒子的压力：

对于血管模型内壁的每个三角形面片，使用七点规则对其采样七个点，这七个点分别为三角形的重心以及每条中线由重心分割出来的两条线段的中点；若这七个点两两之间的最大距离大于核半径，则细分该三角形面片并对细分后的每个三角形面片采样点。重复这个过程直到三角形面片上的七个采样点两两之间的最大距离小于等于核半径为止。将这些采样后的点保存下来，接下来同样参与血液粒子各种物理量的加权运算。

s14，在血管模型内创建一组血液粒子：

在血管模型中，根据初始设定的血液体积以及真实血液中血浆、红细胞、白细胞和血小板所占比例来创建血液粒子，其中这4种粒子分别占总的粒子数目的55％、43％、1％和1％，然后渲染并显示血液粒子。

s15，计算出血管模型的包围盒，并将其划分为均匀网格：

分别找到血管模型顶点的最大三维坐标值xmax、ymax、zmax构成血管模型包围盒的上限点u(xmax，ymax，zmax)，分别找到血管模型顶点的最小三维坐标值xmin、ymin、zmin构成血管模型包围盒的下限点l(xmin，ymin，zmin)，从而构造出血管模型的包围盒；将包围盒划分成大小相同的立方体单元格网格，每个小立方体单元格的边长大小等于核半径的大小，分别将包围盒的长、宽和高除以单元格边长，从而构造网格分辨率向量res(xn,yn,zn)。

步骤s2，具体如下：

s21，对于每个血液粒子，计算出其所在网格并排序其所携带数据：

设当前的血液粒子为第i个粒子。将其位置坐标posi减去包围盒下限点l，并将得到的三个维度的插值除以单元格边长并向下取整，从而可得到该粒子所在单元格的三维索引idx(x,y,z)。

使用如下公式计算粒子i所在单元格的哈希值hash：

hash＝idx.x*res.yn*res.xn+idx.y*res.xn+idx.x

所有计算出来hash值相同的粒子都在同一个单元格中。按关键值对所有血液粒子所携带的物理数值进行排序，其中关键值即为粒子所在网格的哈希值。

s22，对于每个血液粒子，搜索其固定半径内的邻域粒子，然后计算其密度和压强：

设当前的血液粒子为第i个粒子。搜索粒子i的核半径内的邻域粒子的算法如下：对于粒子i所在单元格，遍历其与其相邻的26个单元格，若单元格的不在包围盒中或者不存在血液粒子，则不予以处理，否则遍历单元格中的所有血液粒子，计算每个血液粒子与粒子i的欧拉距离，若距离小于等于核半径，则血液粒子为粒子i的邻域粒子，否则不是。

使用如下公式计算粒子i的密度ρi：

其中j表示粒子i的邻域粒子索引，mj表示邻域粒子j的质量，wpoly6表示核函数为poly6函数，表示以粒子i为起点、以邻域粒子j为终点的向量，h表示核半径。

其中poly6函数如下式所示：

其中表示向量的模长，即粒子i到邻域粒子j的欧拉距离。

使用如下公式计算粒子i的压强pi：

pi＝k(ρi-ρrest)

其中k为一个与温度相关的气体常数，ρrest为血液的静态密度。

s23，对于每个血液粒子，搜索其固定半径内的邻域粒子，计算其形变张量、形变率和粘性系数，从而计算出其切应力：

设当前的血液粒子为第i个粒子，使用如下公式计算粒子i的形变张量di：

其中表示粒子i的速度的雅克比矩阵，表示克罗内克积，wspline表示使用spline函数作为核函数，表示spline函数的梯度。

其中spline函数的梯度的计算公式如下：

然后计算粒子i的所受切应力τi：

τi＝μi*di

其中dri表示粒子i的形变率，trace(di*di)表示矩阵di*di的迹，μ∞和μ0分别表示粒子i的最小和最大粘性系数。

s24，对于每个血液粒子，搜索其固定半径内的邻域粒子，计算其所受的压力和粘性力：

设当前的血液粒子为第i个粒子，

使用如下公式计算粒子i所受压力

其中表示核函数spiky函数的梯度。的计算公式如下：

然后计算粒子i所受粘性力

s25，对于每个血液粒子，计算其加速度和下一时刻的速度，进而计算其下一时刻的位置，并更新其在场景中的位置：

设当前的血液粒子为第i个粒子，根据修正的n-s方程，粒子i的加速度ai的计算公式如下：

然后根据蛙跳算法计算血液粒子的下一时刻的位置：

其中和分别表示粒子i在时刻t的速度和位移，δt为时间步长。然后再判断血液粒子是否与血管内壁发生碰撞修正粒子的速度和下一时刻的位置根据粒子i下一时刻的位置并重新渲染和显示血液粒子i。

步骤s3，具体如下：

设当前的血液粒子为第i个粒子，判断其速度是否小于某一阈值或者粘性系数是否大于某一阈值，若是，则标记粒子i为血栓粒子，否则不做处理。

若粒子i为血栓粒子，则按照步骤s22描述的固定邻域搜索算法找到粒子i的所有邻域粒子。若其邻域粒子中存在血栓粒子，则对于每个邻域血栓粒子，判断其所属血栓块是否可再生长，即粒子i的体积加上血栓块的体积是否依然不大于以血管直径为直径的球的体积；若可以再生长，则将粒子i添加到血栓块中，并相对于血栓块静止；此时，血栓块的体积增加了粒子i的体积，质量增加了粒子i的质量，血栓块的密度更新为增加了粒子i后质量和体积的商。

此后粒子i将随血栓块的运动而运动，即同一个血栓块上的所有粒子都有相同的速度和位移。血栓块所受合力等于组成血栓块的所有粒子所受力之和。对于血栓块上的每一个粒子，按照步骤s22描述的固定邻域搜索算法搜索其邻域粒子，对于邻域粒子中非血栓粒子的粒子，计算它们对血栓粒子的压力和粘性力，加上血栓粒子所受重力，即为血栓粒子所受合力。从而可以计算出整个血栓块所受合力和加速度，进而可以计算出血栓块的速度和位移，从而可以模拟出血栓块的运动。

本发明实施例提出的基于cross模型的多粒子混合修正的血栓模拟方法能够较为真实地模拟血液，并且模拟血液中血栓的形成。正确模拟血液中血栓的形成，有利于对血栓的研究，可以应用于虚拟手术、医学辅助研究等医学应用；同时，本发明技术方案可以满足游戏与影视行业相关方面的需求。

本领域普通技术人员可以理解上述实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器(rom，readonlymemory)、随机存取存储器(ram，randomaccessmemory)、磁盘或光盘等。

另外，以上对本发明实施例所提供的一种基于cross模型的多粒子混合修正的血栓模拟方法进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。