CN112617791B

CN112617791B - 一种模拟血管搏动的方法、装置及存储介质

Info

Publication number: CN112617791B
Application number: CN202011531414.9A
Authority: CN
Inventors: 高琪; 林志洪; 刘星利; 魏润杰
Original assignee: Hangzhou Shengshi Technology Co ltd
Current assignee: Hangzhou Shengshi Technology Co ltd
Priority date: 2020-12-22
Filing date: 2020-12-22
Publication date: 2022-09-20
Anticipated expiration: 2040-12-22
Also published as: CN112617791A

Abstract

本发明公开了一种模拟血管搏动的方法、装置及存储介质。据医学图像提取血管的点云提取中心线；对中心线进行包含切向量获取的前处理；使用从医学图像中提取血管的几何模型，网格划分处理获取网格节点的坐标，对任意网格节点利用前处理后的中心线数据获得控制网格节点运动的运动矢量和运动幅度函数；结合计算流体力学方法，输入冠脉几何模型的网格文件，输入运动矢量和运动幅度函数，设置边界条件，参数和物理模型，模拟血管搏动及血管内部流场。本发明根据血管中心线获取运动矢量控制网格节点运动，模拟在血管搏动下的血液流动状态，避免了采用流固耦合方法和测量血管壁材料属性，减少由刚性壁假设引入的误差，节省计算资源，提高计算速度。

Description

一种模拟血管搏动的方法、装置及存储介质

技术领域

本发明属于生物医学工程领域的一种处理模拟血管脉动方向和幅度的方法，具体是涉及了一种模拟血管搏动方法。

背景技术

血管是具有弹性的管状腔体，在一个心动周期内，血管在血压作用下进行舒张或收缩的运动来适应血压的变化。当前，由于计算流体力学具有较高的时空分辨率，且其准确性和可用性得到了相关研究的验证，因此，借助计算流体力学模拟血管内部血液流动成为认知血管内部流动状况和获取血流动力学参数常用方法。然而在大部分的计算流体力学模拟中，血管壁都被假设为刚性壁，忽略其运动对血液流动产生的影响，也无法模拟血管出现大幅度运动的病症，如心肌桥。

有的模拟为了还原血管的运动，采取了流固耦合的方法，然而流固耦合有两个明显的缺点：首先，流固耦合需要耗费较多的计算资源，计算速度慢；然后，流固耦合的方法需要定义血管壁的力学参数，如杨氏模量或泊松比等，当前已知的参数测量方式是对取样组织进行拉伸实验，但这种方法是有创的且效率比较低。

发明内容

为了能够在计算流体力学模拟中还原血管搏动，且在计算中避免使用流固耦合这种计算资源需求大的方法，本发明提出一种模拟血管脉动的方法。

本发明方法针对用于计算流体力学模拟的网格，根据血管的中心线获取血管壁网格节点的运动矢量，控制每个网格节点在不同时间点的位置，最终模拟出在血管脉动影响下的内部血液流动状态。

本发明方法的使用既能减少由刚性壁假设引入的误差，也能节省计算资源，提高计算速度。血管搏动是通过网格节点的运动实现的，需要规定好网格节点的运动方向来避免网格发生交错，发生网格交错的话计算就会出错。

本发明采用的技术方案是：

本发明方法主要涉及两个主要步骤，中心线数据的前处理和网格节点运动方向的计算，中心线数据的前处理流程图如图1所示。

如图1所示，方法包括：

S1、根据医学图像提取血管的点云，进一步提取中心线；具体实施的血管区域的中心线为冠脉。

S2、对中心线进行包含切向量获取的前处理；具体包含分段、光滑、切向量获取等步骤。

S3、从医学图像中提取血管的几何模型，进行网格划分处理，获取网格节点的坐标，针对其中任意网格节点，利用前处理后的中心线数据获得用于控制网格节点运动的运动矢量和运动幅度函数，作为网格节点运动方向；

S4、输入冠脉几何模型的网格文件，输入上述获得网格节点的运动矢量和运动幅度函数，设置边界条件，参数和物理模型，结合计算流体力学方法模拟搏动的血管的内部流场。

所述S2具体如下：

S2.1、中心线数据的结构：

中心线是由按照血液流动路径排序的若干个离散在三维空间的空间离散点构成的，对于沿血液流动方向的每两个相邻节点构成父子节点，其中以沿血液流动方向处于上游的一个节点作为父节点，以沿血液流动方向处于下游的另一节点作为子节点，中心线的每个节点包含自身节点的坐标和自身节点的父节点索引；

即对于每一个节点，作为子节点时，其父节点是指沿着与血液流动相反的方向，距离子节点最近的节点。

所述的血液流动方向具体是指血液从血管近心端流向远心端的方向。

中心线的每一个节点都有一个编号，父节点索引就是父节点的编号，每个子节点都有唯一的父节点索引。

S2.2、计算切向量：

针对每个中心线的节点，计算切向量，切向量处理如图2所示；

S2.3、对中心线进行分段：

以入口端节点、分岔点、出口端节点为基准对中心线进行分段，在获得的各个段内对节点进行重新编号，形成分段编号；分段效果如图3所示，以光滑中心线。

每一段的起点是中心线的出口端节点或分岔点，从起点开始，根据父节点索引对中心线进行遍历，遍历过程中把节点划分到对应的分段，直到当前节点为下一个分岔点或中心线入口端节点。最终每个分段的起点或终点分别可能是中心线的入口端节点，分岔点或出口端节点。

S2.4、计算分段后每一段的切向量振荡值；

S2.5、针对每一条中心线分段后的每一段，进行光滑拟合；

S2.6、再次计算切向量：

针对光滑后的中心线，再采用和S2.2相同方式处理计算更新获得新的切向量。

至此，中心线数据的前处理结束。

所述S2.2具体如下：

首先，节点的三维坐标随着节点编号变化，节点的三维坐标写成关于节点编号的函数，具体且对每个中心线的各个节点的三维坐标进行泰勒展开：

其中，f(p_i+1)表示节点p_i+1的三维坐标在某一方向的分量(如x,y或z)，p_i表示第i个节点，p_i+1为p_i的子节点，h_pi代表当前节点p_i与其子节点之间的距离；n 表示第n阶导数，f”(p_i)表示函数f的二阶导数，f⁽ⁿ⁾(p_i)表示f的n阶导数；

函数f是节点的x,y或z坐标，x,y或z坐标的一阶导数分别是切向量的三个分量。

然后，计算各个节点的切向量：

对于中心线上没有父节点的节点，如中心线入口端节点，则采用以下公式处理获得节点的切向量：

其中，f'(p_i)表示节点p_i+1的三维坐标在方向分量的一阶导数，即切向量在该方向的分量；当前节点为p_i，节点p_i+1的父节点为节点p_i，节点p_i+2的父节点为节点p_i+1；

对于没有子节点的节点(中心线的末端)或有不少于两个子节点的节点(如中心线的分岔点)，则采用以下公式处理获得节点的切向量：

其中，当前节点为p_i，节点p_i的父节点为节点p_i-1，节点p_i-1的父节点为节点 p_i-2。

对于除了中心线上没有父节点的节点、没有子节点的节点和有不少于两个子节点的节点之外的其他节点，均采用以下公式处理获得节点的切向量：

其中，当前节点编号为p_i，节点p_i+1的父节点为节点p_i，节点p_i的父节点为节点p_i-1。

所述S2.4具体为：

按照以下公式处理获得每一段的切向量振荡值：

其中，OSC表示切向量振荡值，

为分段中第i个节点的切向量，

为分段中所有节点的切向量矢量和的模，

为分段中所有节点的切向量模的和；

计算的切向量振荡值OSC用于表征当前分段的切向量变化幅度，最小值为 0，最大值为0.5，值越大则说明分段的中心线走向越不稳定，中心线经常出现拐弯；反之则分段的中心线走向趋于一致，中心线没有明显变向。

所述S2.5具体为：

针对每一条中心线分段后的每一段，以每个节点的分段编号为自变量，以每个节点的坐标分量(x,y或z值)为因变量，进行拟合光滑处理(如傅里叶拟合，多项式拟合等)，从而对中心线进行光滑，光滑效果如图4所示；

具体是建立坐标分量关于分段编号的表达式，重新代入分段编号，获得光滑后的节点坐标值：

(x,y,z)_smooth＝(f_x(n),f_y(n),f_z(n))

其中，f_x(n)，f_y(n)，f_z(n)分别是第n段中节点的x,y,z坐标值，n为分段节点编号，(x,y,z)_smooth则光滑后的节点的坐标值。

对于某些情况，拟合前可能需要对分段进行再次划分，获取子分段，是否进行这一步操作由分段的OSC值和分段中的节点个数决定。然后进行判断：

若当前段中的节点个数小于等于7或当前段的切向量振荡值小于等于0.05，则当前段不进行再分段；

若当前段的切向量振荡值大于0.05且当前段中的节点个数大于7，则当前段进行再分段，子分段的个数视具体情况而定，使得再分段获得的子分段中的节点个数不少于7；

对于再分段获得的子分段，再分别进行拟合光滑处理，获取光滑后的各个节点坐标。

所述S3具体为：

从医学图像中提取血管的几何模型作为初始网格，如图5所示。选取任意一个网格节点FN，寻找与网格节点FN距离最近的邻近中心线节点N2，并且寻找中心线节点N2的邻近父节点N1；根据网格节点FN、邻近中心线节点N2、邻近父节点N1的坐标，按照步骤S2中的相同方式处理获得邻近中心线节点N2 的切向量N2TAN、邻近父节点N1的切向量N1TAN，进而处理获得网格节点FN 到中心线的法向量FNNOR，作为网格节点FN的运动矢量，使得网格节点FN 沿着该法向量的进行运动，其方向指向中心线，其长度等于初始网格节点到中心线的垂直距离。

所述处理获得网格节点FN到中心线的法向量FNNOR，具体采用以下四种方式的其中一种：

下面介绍四种计算FNNOR的方式，图6-图9分别为这四种方法的示意图。

如图6所示，第一方法：

1)计算网格节点FN分别到邻近父节点N1和邻近中心线节点N2的距离，分别记为h1和h2；

2)根据邻近中心线节点N2的切向量N2TAN、邻近父节点N1的切向量 N1TAN计算网格节点FN的切向量FNTAN，表达式为：

其中，h1表示网格节点FN到邻近父节点N1的距离，h2表示网格节点FN 到邻近中心线节点N2的距离；

3)记邻近父节点N1在网格节点FNTAN的切向量FNTAN上的投影为点P，根据网格节点FN的切向量FNTAN计算点P的矢量FNP：

其中，P,FN,N1分别是点P、网格节点FN和邻近父节点N1的坐标，以矢量的形式来表示，||FNTAN||为FNTAN的模；

4)最后根据矢量计算FNNOR，表达式为：

FNNOR＝(N1-FN)-FNP。

如图7所示，第二方法：

3)由网格节点FN、邻近父节点N1、邻近中心线节点N2三点确定一个平面，求取平面的法向量NOR：

NOR＝(N1-FN)×(N2-FN)

其中，FN,N1,N2分别是网格节点FN、邻近父节点N1和邻近中心线节点N2 的坐标，以矢量的形式来表示；

4)将平面的法向量NOR与网格节点FN的切向量FNTAN叉乘，获得既在平面内又垂直于切向量FNTAN的矢量NOR1：NOR1＝FNTAN×NOR；

5)根据矢量NOR1的长度获得FNNOR：

其中，||NOR1||为NOR1的模；邻近父节点N1在NOR1上的投影为点P， P＝FNNOR+FN。

如图8所示，第三方法：

1)计算邻近父节点N1和邻近中心线节点N2之间的连接向量N1N2，表达式为：

N1N2＝N2-N1

其中，N1,N2分别是邻近父节点N1和邻近中心线节点N2的坐标；

2)平移矢量N1N2，使得矢量N1N2的起点为FN，记邻近父节点N1在平移后的N1N2上投影为点P，根据连接向量N1N2计算点P的矢量FNP：

其中，P,FN,N1分别是点P、网格节点FN和邻近父节点N1的坐标，以矢量的形式来表示；

3)最后根据矢量计算FNNOR，表达式为：

FNNOR＝(N1-FN)-FNP。

在中心线节点密度足够高的情况下，第三方法的效果与第一方法是一致的。

如图9所示，第四方法：

1)由网格节点FN、邻近父节点N1、邻近中心线节点N2三点确定一个平面，求取平面的法向量NOR：

NOR＝(N1-FN)×(N2-FN)

2)将平面的法向量NOR与连接向量N1N2叉乘，获得既在平面内又垂直于连接向量N1N2的矢量NOR1：NOR1＝N1N2×NOR；

3)根据矢量NOR1的长度获得FNNOR：

其中，||NOR1||为NOR1的模；邻近父节点N1在NOR1上的投影为点P， P＝FNNOR+FN。在中心线节点密度足够高的情况下，第四方法的效果与第二方法是一致的。

所述S3中，处理获得以下运动幅度函数g(x_P,t_k)：

g(x_P,t_k)＝(1+B(x_P))·A(t_k)

其中，B(x_P)表示和刻度x_P有关的空间相关函数，A(t_k)表示和时刻t_k有关的时间相关函数；函数A(t_k)用控制与时间相关的运动幅度，函数B(x_P)用于控制与空间位置相关的运动幅度。人体血压随时间的变化是一个周期函数。

所述S4中，按照以下公式利用运动幅度函数g(x_P,t_k)，在第k个时刻对网格节点FN的坐标FN_k进行调整运动：

FN_k＝FN+g(x_P,t_k)·FNNOR

其中，x_P与求解运动矢量过程中得到的投影点P相关，中心线的节点坐标编号也为中心线上的刻度，x_P是网格节点FN的中心线投影在中心线上的刻度，以此来表征网格节点FN的中心线投影在中心线上的相对位置；t_k表示第k个时刻，g(x_P,t_k)表示在第k个时刻中心线刻度为x_P的网格节点的运动幅度，FNNOR 表示网格节点FN到中心线的法向量。

所述的时间相关函数A(t_k)经过傅里叶分解为若干个余弦函数的线性组合：

其中，N代表自然数集，j代表线性组合中的第j项，a_j,ω_j,b_j都是分解过程中得到的第一、第二、第三常数。

通常情况下认为血管的运动幅度是均匀的，所以空间相关函数B(x_P)取常数 1，即g(x_P,t_k)＝A(t_k)；在特殊的情况，血管的运动幅度才会随空间位置变化，例如心肌桥区域但不限于心肌桥区域。以心肌桥为例，心肌桥血管受心肌包裹，在心动周期的收缩期被挤压并形成狭窄，根据心肌桥位置构建获得空间相关函数B(x_P)。

例如在心肌桥区域但不限于心肌桥区域，所述的空间相关函数B(x_P)具体为：

其中，

和

分别是心肌桥区域的网格节点中刻度x_P的最小值和最大值；c为心肌桥收缩最严重时对应的狭窄率。

对于不同的法向量FNNOR处理，x_P的计算方法不同，分别采用不同方式计算获得网格节点FN的中心线投影在中心线上的刻度x_P为：

当使用第一方法计算FNNOR时，对于网格节点FN，其x_P表达式为：

其中，num(N1)为邻近父节点N1在心肌桥中心线的节点编号，num(NC)为最严重位置中心线节点NC在心肌桥中心线的节点编号；N1,N2分别是邻近父节点 N1和邻近中心线节点N2的坐标；FNP表示从节点FN指向点P的矢量，|| ||代表矢量的模，abs为绝对值；

当使用第二方法和第四方法计算FNNOR时，对于节点FN，其x_P表达式为：

其中，num(N1)为邻近父节点N1在心肌桥中心线的节点编号，num(NC)为最严重位置中心线节点NC在心肌桥中心线的节点编号；N1,N2,P分别是邻近父节点N1，邻近中心线节点N2和投影点P的坐标；N1P表示从节点N1指向点P 的矢量，|| ||代表矢量的模，abs为绝对值；

当使用第三方法计算FNNOR时，对于网格节点FN，其x_P表达式为：

其中，num(N1)为邻近父节点N1在心肌桥中心线的节点编号，num(NC)为最严重位置中心线节点NC在心肌桥中心线的节点编号；N1,N2分别是邻近父节点 N1和邻近中心线节点N2的坐标；FNP表示从节点N1指向点P的矢量，|| ||代表矢量的模，abs为绝对值。

所述的最严重位置中心线节点NC按照以下方式确定。

根据医学图像，确定心肌桥区域对应的网格和中心线，并在心肌桥区域独立划分为一个新的分段，心肌桥区域对应的中心线作为心肌桥中心线；以心肌桥中心线的近心端节点为分段起点，远心端终点为分段终点，重新对内部的中心线节点进行编号；已知x_P表征网格节点FN的中心线投影在中心线上的相对位置，中心线节点编号是中心线的刻度，x_P的值借助已知的重新编号的中心线节点来确认。然后确认心肌桥收缩最严重的位置，记录最严重的位置对应的最严重位置中心线节点NC，实验发现心肌桥区域血管的运动是距离中心点NC越近，则运动幅度越大。

二、一种模拟血管搏动处理系统：

医学图像采集设备，采集血管的医学图像并提取获得中心线；

前处理模块，对中心线进行包含切向量获取的前处理；

运动矢量模块，从医学图像中提取血管的几何模型进行网格划分获取网格节点，利用前处理后的中心线对网格节点处理获得运动矢量和运动幅度函数；

计算流体力学模块，输入冠脉几何模型的网格文件，输入上述获得网格节点的运动矢量和运动幅度函数模拟血管搏动的运动，并获得搏动的血管内部流场。

所述的前处理模块中具体以下处理：

S2.1、建立中心线数据的父子节点关系和索引结构；

S2.2、计算切向量：针对每个中心线的节点，计算切向量；

S2.3、对中心线进行分段：

以入口端节点、分岔点、出口端节点为基准对中心线进行分段，在获得的各个段内对节点进行重新编号，形成分段编号；

S2.4、计算分段后每一段的切向量振荡值；

S2.5、针对每一条中心线分段后的每一段，进行光滑拟合；

S2.6、针对光滑后的中心线，再采用和S2.2相同方式处理计算更新获得新的切向量。

所述运动矢量模块具体处理为：选取任意一个网格节点FN，寻找与网格节点FN距离最近的邻近中心线节点N2，并且寻找中心线节点N2的邻近父节点 N1；根据网格节点FN、邻近中心线节点N2、邻近父节点N1的坐标，按照步骤S2中的相同方式处理获得邻近中心线节点N2的切向量N2TAN、邻近父节点 N1的切向量N1TAN，进而处理获得网格节点FN到中心线的法向量FNNOR，作为网格节点FN的运动矢量，且记录下不同方法计算得到的投影点P坐标。

三、一种存储介质，存储有计算机程序，所述计算机程序被处理器执行时实现上述方法。

所述计算机程序为对应实现模拟血管搏动处理方法的指令。

本发明的有益效果是：

本发明方法的计算结果(网格节点法向量)可与计算流体力学方法结合，将计算结果用于控制网格的运动，从而模拟出搏动的血管的内部流场。

将本发明方法处理获得的运动矢量和运动幅度函数能代替流固耦合方法，只需要较少的输入数据(无需考虑血管材料属性)，就能模拟受血管脉动影响下的血管内部流场；而且与流固耦合方法相比，因为省去用于模拟血管壁运动的有限元计算，避免了流固耦合方法需要提供血管壁的材料属性才能处理的问题(实际中血管壁的材料属性是难以获得的)，能大幅度缩短数据前处理时间和流场计算时间。

本发明针对每个独立的案例只需执行一次，计算结果不需要在计算流体力学模拟的迭代过程重复更新，可直接调用，所以能节省计算时间而且算法稳定性也比较高。

本发明方案适用于网格文件中的所有节点，包括壁面节点和内部节点，使用者可以选择控制全部网格节点的运动或只控制壁面网格的运动。如果只控制壁面网格运动，则内部节点坐标的更新需要调用额外的动网格算法，计算时间较长；如果同时控制所有网格节点运动，则不需要调用额外的动网格算法，直接使用本发明方案的结果更新坐标，能进一步地节省时间。具体情况视使用者的选择而定。

附图说明

图1：中心线数据的前处理流程图。

图2：中心线节点的切向量示意图。

图3：中心线分段示意图，不同分段的节点使用不同的符号表示。

图4：中心线的光滑效果，光滑前的中心线节点用“.”表示，光滑后的中心线节点用“+”表示。

图5：提取出血管的几何模型后进行网格划分处理情况图。

图6：使用切向量计算网格节点法向量的第一方法示意图

图7：使用向量叉乘和切向量计算网格节点法向量的第二方法示意图。

图8：直接使用N1N2计算网格节点切向量的第三方法示意图。

图9：使用向量叉乘和N1N2计算网格节点切向量的第四方法示意图。

图10：从CT医学图像中提取点云，最后获取中心线数据和几何模型。

图11：模拟冠脉搏动的结果图。外层为冠脉收缩前的边界，内部为收缩后的边界，冠脉血管壁的运动引起了压力分布的变化图。

图12：模拟心肌桥情况下的血管运动，外层为心肌收缩前的冠脉形态；内层为心肌收缩后冠脉被压缩的形态；形态的变化引起压力分布的变化图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细说明。

本发明的实施例如下：

实施例一：模拟冠脉整体运动

首先，需要准备中心线和网格数据。

根据医学图像提取冠脉点云，在点云基础上提取冠脉中心线数据(包括中心线节点和父节点索引)和冠脉几何模型，图10中以CT图像为例，展示从医学图像到最后得到中心线和几何模型的过程。根据技术方案的主要步骤一，处理中心线数据，获取光滑后的中心线节点坐标及对应的切向量。

接着，根据技术方案中的主要步骤二，创建计算流体力学模拟的计算案例。

本发明提供的方法，可结合计算流体力学方法模拟血管的搏动；本发明方案的计算结果能为计算流体力学模拟过程提供每个时刻的网格节点坐标，通过更新坐标的方式使网格动起来，而每个时刻的坐标计算，都需要使用上述技术方案中主要步骤二中得到的FNNOR。

本实施例中，阐述如何使用第一方法计算运动矢量。对于任意一个初始网格节点FN，其坐标为FN＝(x,y,z)_FN，遍历中心线节点，寻找距离FN最近的节点N2，其坐标为N2＝(x,y,z)_N2；N2的父节点为N1，其坐标为N1＝(x,y,z)_N1。FN 与N1之间的距离为h1，FN与N2之间的距离为h2，根据N1的切向量N1TAN和N2的切向量N2TAN，插值出FN的切向量FNTAN，

节点FN到切向量FNTAN的投影为P，P的坐标为P＝(x,y,z)_P,计算矢量FNP:

进一步地，计算运动矢量FNNOR, FNNOR＝(N1-FN)-FNP。

遍历每个初始网格节点，计算每个网格节点对应的FNNOR。

在获得所有初始网格节点的运动矢量后，设置计算流体力学模拟案例，输入冠脉几何模型的网格文件；对冠脉出入口设置边界条件，边界条件可以设置为压力，流量，速度分布或者耦合降阶模型等；设置血液的物理特性，包括密度和粘度等；设置湍流模型。

输入上述计算的初始网格节点运动矢量，模拟血管壁的运动。本发明方法可直接用于全体网格的坐标更新或只用于壁面网格的坐标更新，如果只控制壁面网格的坐标更新，则需要额外调用动网格算法。

本实施例中，使用更新网格节点的方式控制血管模拟运动，因为本实施例模拟正常的血管搏动，空间位置对运动幅度没有影响，因此使用的幅度控制函数g(x_P,t_k)可简化为只与时间相关的A(t_k)。已知原始坐标为FN＝(x,y,z)_FN，网格节点FN的运动矢量为FNNOR，设置运动控制函数在第k个时刻的网格节点 FN的坐标FN_k调整为FN_k＝FN+A(t_k)·FNNOR，其中A(t_k)为第ik时刻t_k的常态运动控制函数，运动控制函数大小的变化控制了网格节点位移的幅度。

运动矢量FNNOR根据原始网格计算，运动矢量FNNOR的模等于节点到中心线的距离，当常态运动控制函数A(t_k)＝0.8时，则在第k个时刻t_k，网格节点到中心线的距离为原始距离的0.8。常态运动控制函数A(t_k)的具体形式要考虑到人体血管内部的压力时周期性变化。

具体是的常态运动控制函数A(t_k)为周期函数，经过傅里叶分解，转化为若干个余弦函数的线性组合，表示为：

其中，N代表自然数集，j代表线性组合中的第j项，a_j,ω_j,b_j都是分解过程中得到的第一、第二、第三常数；

本实施例中，

为例，在图11中展示两个时刻下网格的形状变化和计算结果，其中T为一个心动周期的时间，取0.8s。

实施例二：心肌桥

冠脉内部随时间变化的压力分布引起血管正常的收缩和舒张运动。然而，在特殊情况下，血管会出现异常的形变，比如心肌桥。有时冠脉被埋在心肌中，当心室收缩的时候，冠脉被挤压导致出现狭窄。

实施例二针对心肌桥的模拟，在实施例一的基础上进行。

根据医学图像，确定心肌桥区域，如图12所示。

计算获得每个节点的运动矢量后，在心肌桥区域，针对心肌桥区域使用与时间和空间相关的函数g(x_P,t_k)控制心肌桥区域的网格节点运动，而剩余的区域保留使用只与时间相关的常态运动控制函数A(t_k)来控制冠脉中的血管模拟运动。

在心肌桥区域的网格节点更新方式改为：FN_k＝FN+g(x_P,t_k)·FNNOR，在非心肌桥区域保留FN_k＝FN+A(t_k)·FNNOR的网格节点更新，g(x_P,t_k)如下表达式：

g(x_P,t_k)＝(1+B(x_P))·A(t_k)

其中，A(t_k)为控制血管运动的第k个时刻t_k的常态运动控制函数；c为常数，代表心肌桥收缩最严重时的狭窄率；x_P表示在心肌桥区域，网格节点到心肌桥中心线的投影在心肌桥中心线的相对位置，与网格节点的位置相关，x代表反映了投影点P到中心点的距离；根据不同的FNNOR计算方法，x_P的计算方式也有区别，以第一方法为例，x_P的表达式为：

其中，NC为心肌桥的中心点，是血管被挤压最严重的位置；num(N1)为N1 在心肌桥中心线的节点编号，num(NC)为NC在心肌桥中心线的节点编号；N1,N2 分别是邻近父节点N1和邻近中心线节点N2的坐标；FNP表示从节点FN指向点P的矢量，|| ||代表矢量的模，abs为绝对值。

因此，相对坐标x越大，到中心点NC的距离就越小，反之则距离越大。

通过心肌桥运动控制函数g(x_P,t_i)控制心肌桥区域的网格节点运动，因为函数B的基本形式为余弦函数，所以最终的模拟结果体现为：网格节点的x_P最大或最小(距离中心点NC最远)时B(x_P)值为0，网格节点的x_P等于0(网格节点的投影与中心点NC重合)时位移幅度最大即B(x_P)值为1。图12展示了两个时刻下心肌桥区域网格的形状变化及计算结果。

Claims

1.一种模拟血管搏动的方法，其特征在于：方法包括：

S1、根据医学图像提取血管的点云，进一步提取中心线；

S2、对中心线进行包含切向量获取的前处理；

S3、从医学图像中提取血管的几何模型，进行网格划分处理，获取网格节点的坐标，针对其中任意网格节点，利用前处理后的中心线数据获得用于控制网格节点运动的运动矢量和运动幅度函数；

S4、输入冠脉几何模型的网格文件，输入上述获得网格节点的运动矢量和运动幅度函数，设置边界条件，参数和物理模型，结合计算流体力学方法模拟搏动的血管的内部流场；

所述S3具体为：选取任意一个网格节点FN，寻找与网格节点FN距离最近的邻近中心线节点N2，并且寻找中心线节点N2的邻近父节点N1；根据网格节点FN、邻近中心线节点N2、邻近父节点N1的坐标，按照步骤S2中的相同方式处理获得邻近中心线节点N2的切向量N2TAN、邻近父节点N1的切向量N1TAN，进而处理获得网格节点FN到中心线的法向量FNNOR，作为网格节点FN的运动矢量；

所述S3中，处理获得以下运动幅度函数g(x_P,t_k)：

g(x_P,t_k)＝(1+B(x_P))×A(t_k)

其中，B(x_P)表示和刻度x_P有关的空间相关函数，A(t_k)表示和时刻t_k有关的时间相关函数；

FN_k＝FN+g(x_P,t_k)×FNNOR

其中，x_P是网格节点FN的中心线投影在中心线上的刻度，以此来表征网格节点FN的中心线投影在中心线上的相对位置；t_k表示第k个时刻，g(x_P,t_k)表示在第k个时刻中心线刻度为x_P的网格节点的运动幅度，FNNOR表示网格节点FN到中心线的法向量。

2.根据权利要求1所述的一种模拟血管搏动的方法，其特征在于：所述S2具体如下：

S2.1、中心线数据的结构：

中心线是由按照血液流动路径排序的若干个空间离散点构成的，对于沿血液流动方向的每两个相邻节点构成父子节点，其中以沿血液流动方向处于上游的一个节点作为父节点，以沿血液流动方向处于下游的另一节点作为子节点，中心线的每个节点包含自身节点的坐标和自身节点的父节点索引；

S2.2、计算切向量：针对每个中心线的节点，计算切向量；

S2.3、对中心线进行分段：以入口端节点、分岔点、出口端节点为基准对中心线进行分段，在获得的各个段内对节点进行重新编号，形成分段编号；

S2.4、计算分段后每一段的切向量振荡值；

S2.5、针对每一条中心线分段后的每一段，进行光滑拟合；

S2.6、再次计算切向量：针对光滑后的中心线，再采用和S2.2相同方式处理计算更新获得新的切向量。

3.根据权利要求2所述的一种模拟血管搏动的方法，其特征在于：所述S2.2具体如下：

首先，且对每个中心线的各个节点的三维坐标进行泰勒展开：

其中，f(p_i+1)表示节点p_i+1的三维坐标在x,y或z方向的分量，p_i表示第i个节点，p_i+1为p_i的子节点，h_pi代表当前节点p_i与其子节点之间的距离；n表示第n阶导数，f”(p_i)表示函数f的二阶导数，f⁽ⁿ⁾(p_i)表示f的n阶导数；

然后，计算各个节点的切向量：

对于中心线上没有父节点的节点，则采用以下公式处理获得节点的切向量：

其中，f'(p_i)表示节点p_i+1的三维坐标在方向分量的一阶导数；

对于没有子节点的节点或有不少于两个子节点的节点，则采用以下公式处理获得节点的切向量：

4.根据权利要求2所述的一种模拟血管搏动的方法，其特征在于：所述S2.4具体为：按照以下公式处理获得每一段的切向量振荡值：

其中，OSC表示切向量振荡值，

为分段中第i个节点的切向量，

为分段中所有节点的切向量矢量和的模，

为分段中所有节点的切向量模的和。

5.根据权利要求2所述的一种模拟血管搏动的方法，其特征在于：所述S2.5具体为：

针对每一条中心线分段后的每一段，以每个节点的分段编号为自变量，以每个节点的坐标分量为因变量，进行拟合光滑处理，从而对中心线进行光滑；

然后进行判断：

若当前段的切向量振荡值大于0.05且当前段中的节点个数大于7，则当前段进行再分段，使得再分段获得的子分段中的节点个数不少于7；

6.根据权利要求1所述的一种模拟血管搏动的方法，其特征在于：在受压迫的血管区域，所述的空间相关函数B(x_P)具体为：

其中，

和

分别是受到压迫的血管区域的网格节点中刻度x_P的最小值和最大值；c为受压迫的流道血管收缩最严重时对应的狭窄率。

7.采用权利要求1～6任一所述方法的一种模拟血管搏动处理系统，其特征在于包括：

前处理模块，对中心线进行包含切向量获取的前处理；

8.根据权利要求7所述的一种模拟血管搏动处理系统，其特征在于：

所述的前处理模块中具体以下处理：

S2.1、建立中心线数据的父子节点关系和索引结构；

S2.2、计算切向量：针对每个中心线的节点，计算切向量；

S2.3、对中心线进行分段：

S2.4、计算分段后每一段的切向量振荡值；

S2.5、针对每一条中心线分段后的每一段，进行光滑拟合；

9.根据权利要求8所述的一种模拟血管搏动处理系统，其特征在于：所述运动矢量模块具体处理为：选取任意一个网格节点FN，寻找与网格节点FN距离最近的邻近中心线节点N2，并且寻找中心线节点N2的邻近父节点N1；根据网格节点FN、邻近中心线节点N2、邻近父节点N1的坐标，按照步骤S2中的相同方式处理获得邻近中心线节点N2的切向量N2TAN、邻近父节点N1的切向量N1TAN，进而处理获得网格节点FN到中心线的法向量FNNOR，作为网格节点FN的运动矢量，且记录下不同方法计算得到的投影点P坐标。

10.一种存储介质，存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1～6任一所述的方法。

11.根据权利要求10所述的一种存储介质，其特征在于，所述计算机程序为对应实现权利要求1～6任一所述的一种模拟血管搏动的方法的指令。