基于图像的模拟血液流动的方法、装置、设备及存储介质与流程

1.本技术涉及但不限于图像处理领域，尤其涉及一种基于图像的模拟血液流动的方法、装置、设备及存储介质。


背景技术：

2.计算流体力学(computed fluid dynamic，cfd)模拟因其能够获得更丰富的流场信息，被用于模拟心血管内的血液流动状态，以期揭示血液在心血管内的流动规律。目前，虽然有采用流固耦合的方法进行cfd模拟，但流固耦合模拟需占用较大的计算资源，且耗时较长。


技术实现要素：

3.本技术实施例提供一种基于图像的模拟血液流动的方法、装置、设备及存储介质。
4.本技术实施例的技术方案是这样实现的：
5.一种基于图像的模拟血液流动的方法，所述方法包括：
6.获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、所述心血管的血液流动速度波形数据；
7.对所述多个时刻的三维模型进行网格离散，得到所述多个时刻的网格节点坐标；
8.对所述多个时刻的网格节点坐标进行点云配准，得到所述多个时刻的运动序列；
9.通过计算流体力学算法基于所述血液流动速度波形数据以及所述多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟所述心血管内的血液流动。
10.一种基于图像的模拟血液流动的装置，所述装置包括：
11.获取模块，用于获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、心血管的血液流动速度波形数据；
12.处理模块，用于对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标；
13.处理模块，用于对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列；
14.模拟模块，用于通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动。
15.本技术提供一种基于图像的模拟血液流动的设备，所述基于图像的模拟血液流动的设备包括：处理器、存储器和通信总线；
16.通信总线用于实现处理器和存储器之间的通信连接；
17.处理器用于执行存储器中存储的基于图像的模拟血液流动的程序，以实现上述的基于图像的模拟血液流动的方法的步骤。
18.本技术实施例提供一种存储介质，存储介质存储有一个或者多个程序，一个或者多个程序可被一个或者多个处理器执行，以实现上述的基于图像的模拟血液流动的方法的步骤。
19.本技术提供的基于图像的模拟血液流动的方法、装置、设备及存储介质，通过获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、心血管的血液流动速度波形数据；对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标；对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列；通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；可见，本技术在模拟血液流动的过程中，获得一个周期内多个时刻的心血管三维模型，无需考虑血管材料属性，对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标，进而对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列后，便可以通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；避免了流固耦合方法需要提供心血管如左心房及左心耳壁面的材料属性才能处理的问题，大幅度缩短数据前处理时间和计算时间。
附图说明
20.图1本技术实施例提供的基于图像的模拟血液流动的方法的流程示意图；
21.图2本技术实施例提供的某一个时刻的左心房-左心耳的三维模型的示意图；
22.图3本技术实施例提供的体网格划分后的左心房-左心耳的三维模型的示意图；
23.图4本技术实施例提供的点云配准的示意图；
24.图5本技术实施例提供的点云转换的示意图；
25.图6本技术实施例提供的点云配准前后的示意图；
26.图7本技术实施例提供的一种网格节点更新的示意图；
27.图8本技术实施例提供的另一种网格节点更新的示意图；
28.图9本技术实施例提供的基于图像的模拟血液流动的装置的结构示意图；
29.图10本技术实施例提供的基于图像的模拟血液流动的设备的结构示意图。
具体实施方式
30.为了使本技术领域的人员更好地理解本技术方案，下面将结合本技术实施例中的附图，对本技术实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本技术一部分实施例，而不是全部的实施例。基于本技术中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本技术保护的范围。
31.本技术的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或单元的过程、方法、系统、产品或设备没有限定于已列出的步骤或单元，而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、产品或设备固有的其他步骤或单元。
32.在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包
含在本技术的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。
33.目前，大部分的cfd模拟，假设心血管壁为刚性，即使给定患者特异性边界条件，模拟的结果与实际情况还是有一定差距。例如，房颤是临床上最常见的心率失衡之一，其心房呈无序颤动的状态，评估房颤下左心房-左心耳内的血液的真实流动状态，仍然是一件十分困难的事情。除了难以评估房颤下左心房-左心耳内的血液的真实流动状态，评估正常心率下左心室内的规律血液流动以及二尖瓣的运动等同样具有一定的挑战。
34.本技术实施例提供的基于图像的模拟血液流动的设备可以实施为笔记本电脑，平板电脑，台式计算机，移动设备(例如，个人数字助理，专用消息设备)，智能机器人等任意基于图像的模拟血液流动的设备，也可以实施为服务器。下面，将说明基于图像的模拟血液流动的设备实施为基于图像的模拟血液流动的设备时的示例性应用。
35.本技术实施例提供一种基于图像的模拟血液流动的方法，该方法应用于基于图像的模拟血液流动的设备，参见图1所示，该方法包括：
36.步骤101、获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、心血管的血液流动速度波形数据。
37.本技术实施例中，医学图像包括但不限于通过如下方式获得的图像：通过4d计算机断层扫描技术(4d-computed tomography，4d-ct)、核磁共振成像技术(magnetic resonance imaging，mri)采集到的图像。
38.这里，对医学图像中的心血管进行重构处理，可以通过三维重构软件执行重构处理，三维重构软件包括但不限于医学影像交互工具包(medical imaging interaction toolkit，mitk)、医学影像控制系统(materialise's interactive medical image control system，mimics)、三维图像处理软件如simpleware或3d-doctor。
39.这里，多个时刻的数量可以是任意的，同时，多个时刻在一个周期内的分布，可以是均匀的，也可以是非均匀的。
40.在一个可实现的重构处理的场景中，获得医学图像(医学图像可以由动态医学图像数据来表征)，利用三维重构软件对医学图像中的心血管进行自动分割，得到一个周期内i个时刻的心血管的三维模型，其中，i为大于等于2的正数。示例性地，心血管模型包括但不限于左心房-左心耳、左心室、主动脉及腹主动脉。
41.在一些实施例中，一个周期内多个时刻的三维模型可以是上述自动分割得到的i个时刻的心血管的三维模型；一个周期内多个时刻的三维模型也可以是得到i个时刻的心血管的三维模型后，对i个时刻的心血管的三维模型的表面进行光滑处理后，得到的处理后的i个时刻的心血管的三维模型。这里，可以将一个周期内i个时刻的心血管的三维模型分别导入自动化逆向工程软件例如geomagic studio中，将一个周期内i个时刻的心血管的三维模型的表面进行光滑处理。
42.示例性的，以得到的一个周期内多个时刻的三维模型为一个周期内i个时刻的左心房-左心耳的三维模型为例；对应的心血管的血液流动速度波形数据为一个周期内血液流动速度随时间变化的曲线，此时为二尖瓣处的边界条件。
43.在一个可实现的场景中，采用自动化逆向工程软件对重构得到的i个时刻的左心
房-左心耳三维模型进行表面光滑，肺动脉仅保留分叉前的结构，并分割入口和出口平面。入口和出口分别代表肺动脉、二尖瓣，示例性的，i＝20，当然，i还可以取其他值，本技术对此不做具体地限定。图2展示的是某一个时刻的左心房-左心耳的三维模型，本实施例中采用的模型组成部分为：肺静脉、左心房、左心耳及二尖瓣。本技术实施例中，二尖瓣采用平面来替代，二尖瓣结构形变对血液流动的影响已通过入口边界条件来体现，因此将二尖瓣简化成平面可提高计算效率。
44.步骤102、对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标。
45.本技术实施例中，对多个时刻的三维模型进行网格离散的实现方式，包括但不限于体网格离散和面网格离散中的至少一种。
46.示例性的，以得到的一个周期内多个时刻的三维模型为一个周期内i个时刻的左心房-左心耳的三维模型为例；根据预设的网格离散方式对一个周期内i个时刻的左心房-左心耳的三维模型进行网格离散，得到一个心动周期内i个时刻的左心房-左心耳的网格节点坐标；进一步地，在网格离散的过程中，仍旧以i＝20为例，示例性的，选取左心房体积最大时刻(如左心室收缩末期)的左心房-左心耳的三维模型进行体网格离散，对剩余19个时刻的左心房-左心耳的三维模型进行面网格离散，图3展示的是体网格划分后的左心房-左心耳的三维模型。
47.步骤103、对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列。
48.本技术实施例提出的基于图像的模拟血液流动的方法代替相关技术中的流固耦合方法，只需获得一个周期内i个时刻的心血管三维模型，无需考虑血管材料属性，就能模拟心血管例如左心房及左心耳内部流场；而且与流固耦合方法相比，省去用于模拟左心房及左心耳壁面运动的有限元计算，避免了流固耦合方法需要提供左心房及左心耳壁面的材料属性才能处理的问题，经本技术实验发现，实际中左心房及左心耳壁面的材料属性是难以获得的，并且个体差异较大，从而通过本技术提供的方法，能大幅度缩短数据前处理时间和计算时间。
49.这里，流固耦合作用是指流体与固体之间的相互作用。
50.步骤104、通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动。
51.本技术实施例中，计算流体力学算法包括但不限于动网格算法、浸没边界算法；本技术根据i个时刻的运动序列结合计算流体力学算法来模拟心血管内的血液流动，其相比于相关技术中的流固耦合方法，能更简单、快速地还原心血管内部血液的真实流动情况。
52.本技术实施例提供的基于图像的模拟血液流动的方法，通过获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、心血管的血液流动速度波形数据；对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标；对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列；通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；可见，本技术在模拟血液流动的过程中，获得一个周期内多个时刻的心血管三维模型，无需考虑血管材料属性，对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标，进而对多个时刻的网格节点坐标进行点
云配准，得到多个时刻的运动序列后，便可以通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；避免了流固耦合方法需要提供心血管如左心房及左心耳壁面的材料属性才能处理的问题，大幅度缩短数据前处理时间和计算时间。
53.本技术一些实施例中，多个时刻的网格节点坐标中不同时刻的壁面网格节点坐标之间互相独立，本技术通过点云配准将各个时刻的壁面网格节点坐标进行关联；对此，步骤103对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列，可以通过如下步骤实现：
54.a11、选取多个时刻的网格节点坐标中，一个时刻的壁面网格节点坐标作为配准模板点云。
55.a12、选取多个时刻的网格节点坐标中，剩余时刻的壁面网格节点坐标作为配准目标点云。
56.其中，当选取配准模板点云和配准目标点云时，配准目标点云中点的数量大于或等于配准模板中点的数量，如此，可以避免配准过程中出现配准目标点云中的点对应多个配准模板点云中的点。
57.a13、将配准模板点云注册到每一配准目标点云上，得到配准模板点云与每一配准目标点云之间的变换矩阵。
58.a14、基于每一变换矩阵将配准模板点云进行相应的变换，得到一个周期内剩余时刻对应的运动序列。
59.本技术的点云配准过程中，参照图4所示，首先，选取一个周期内i个时刻的心血管的网格节点坐标的一个时刻的心血管的壁面网格节点坐标作为配准模板点云，剩余时刻的心血管的壁面网格节点坐标作为配准目标点云，根据点云配准得到配准模板点云及配准目标点云之间的变换矩阵(tform)。其次，根据变换矩阵tform，得到心血管在一个周期内的不同时刻的运动序列(moving)。本技术实施例中，在得到不同时刻的moving的情况下，还可以将心血管在一个心动周期内的不同时刻的moving以及配准模板点云坐标转换成不同时刻网格节点坐标的二维矩阵，分别为x_c、y_c、z_c，参照图5所示。
60.本技术一些实施例中，还可以对点云配准算法做出优化，进一步提高点云配准的效率，缩短计算时间，为此，a13将配准模板点云注册到每一配准目标点云上，得到配准模板点云与每一配准目标点云之间的变换矩阵，可以通过如下步骤实现：
61.首先，将配准模板点云中的点作为高斯混合模型的中心点；
62.其次，确定配准目标点云中的点符合独立同分布时的似然函数；
63.最后，基于高斯混合模型和似然函数，得到配准模板点云与每一配准目标点云之间的变换矩阵。
64.这里，配准模板点云中的点在配准目标点云中存在对应点。
65.本技术实施例中，对上述得到变换矩阵的过程做出进一步地说明：
66.上述配准目标点云表示为上述配准模板点云表示为本技术实施例中，将配准模板点云当作高斯混合模型的中心，在高斯混合模型中加入一项均匀分布，以考虑噪声以及异常值对模型的影响，这里，可以通过权
重参数μ(0≤μ≤1)，来平衡高斯混合模型及均匀分布之间的关系，可以表示为如下公式(1)：
[0067][0068]
其中，θ表示高斯混合模型的中心点，即配准模板点云中的点；表示点在空间变换后的新坐标位置；σ2为方差；i表示单位矩阵；xn表示当下配准所用的配准目标点云；ym表示当下配准所用的配准模板点云；n表示配准目标点云中点的数量；m表示配准模板点云中点的数量，1≤n≤n，1≤m≤m，n和m为正整数。
[0069]
示例性地，假设配准目标点云中的点符合独立同分布，则似然函数可以表示为如下公式(2)：
[0070][0071]
其中，
[0072]
示例性地，可根据期望最大算法(expectation-maximization，em)来求解公式(2)。
[0073]
进一步地，求解公式(2)可以通过如下b21-b22来实现：
[0074]
b21、计算对应关系矩阵p＝{p
mn
}∈rm×n。
[0075]
其中，点f(ym)和xn之间的对应关系表示为如下公式(3)：
[0076][0077]
其中，
[0078]
b22、通过最小化q函数来恢复点云之间的非刚性变换，同时根据正则化算法，引入一个正则项来平滑插值函数，则q函数可以表示为如下公式(4)：
[0079][0080]
这里，正则化算法包括但不限于基于变分原理的tikhonov正则化算法。
[0081]
示例性地，根据移动相干理论(motion coherent theory，mct)，经过非刚性变换后模板点云的新位置可以表示为原始位置和偏移函数ν之和，即如下公式(5)：
[0082]
f(y)＝y+ν(y) 公式(5)
[0083]
示例性地，根据变分原理即变分法得出最优的偏移函数可以定义为一组高斯核函数的线性叠加，即如下公式(6)：
[0084][0085]
其中，
[0086]
示例性地，根据公式(5)、公式(6)，将公式(4)对f进行求导，并令导数为零，得到表达式如下：
[0087][0088]
其中，φ∈rm×m是一个对称的格拉姆矩阵(gram matrix)且φm∈r1×m表示矩阵φ的第m行。
[0089]
需要说明的是，本技术实施例中引入一个对应关系约束，以降低点云配准的计算时间，该对应关系约束即如下公式(8)：
[0090][0091]
本技术中的对应关系约束表示配准模板点云中的点在配准目标点云中存在对应点，即每个配准模板点云中的点必须要在目标点云中寻找对应点。
[0092]
进一步地，将公式(8)带入公式(7)，则公式(7)可以简化如下的公式(9)：
[0093][0094]
其中，是由组成的数据矩阵，可以看作是由配准目标点云x中所有的点加权得到的配准模板点云中第m个特征点的对应点。
[0095]
示例性地，对矩阵φ进行特征值分解得到φ＝uλu
t
，同时u＝u-1
及u
t
＝(u
t
)-1
，根据公式(9)得到如下公式(10)：
[0096][0097]
由上述可知，本技术通过引入公式(8)，最终在计算公式(10)时，只需要更新一个对角矩阵(λ+λσ2im×m)-1
和进行矩阵相乘，计算量为o(2dm2+dm+m)。然而，在不引入公式(8)对公式(7)进行简化的情况下，则需要o(m3)的计算量。数据维数d远远小于m，因此o(2dm2+dm+m)远小于o(m3)。
[0098]
示例性地，令则可以得到如下公式(11)：
[0099][0100]
其中，表示点云变换后的位置，d(
·
)表示向量对角化操作。
[0101]
b21和b22相互迭代，直到到达预设的停止条件为止，最终获得变换矩阵tform，可见，变换矩阵tform由偏移函数ν转换而来。不同时刻的moving由原始位置和偏移函数之和f(y)转换而来。
[0102]
在一个可实现的优化后的点云配准场景中，以得到的一个周期内多个时刻的三维模型为一个周期内i个时刻的左心房-左心耳的三维模型为例，将其中一个时刻的壁面网格节点坐标作为配准模板点云，其余时刻的壁面网格节点坐标作为配准目标点云，通过点云配准，将配准模板点云注册到配准目标点云上，得到配准模板点云与各配准目标点云间的变换矩阵tform(位移矢量场)。示例性的，本实施例选取左心房体积最大时刻的左心房-左心耳的三维模型的壁面网格节点坐标作为配准模板点云。基于各变换矩阵tform将配准模板点云进行相应的变化，得到基于配准模板点云的各配准目标点云，即一个周期的i-1个时刻的左心房-左心耳的运动序列moving。图4所示为点云配准示意图，其中n、m表示壁面网格
节点数量；i表示一个心动周期内总共重构了i个左心房-左心耳的三维模型，以i＝20为例；(x,y,z)表示壁面网格节点坐标；(q,r,s)表示位移矢量。图6所示为点云配准前后的对比示意图。将配准模板点云以及运动序列moving转换成三个坐标矩阵，分别为x_c、y_c、z_c。图5所示为点云转换示意图。本实施例中，以进行点云配准中的配准模板点云个数约4万为例，则优化前的计算量为o(6.4
×
10
13
)，优化后的计算量为大约为o(9.6
×
109)，相差了近4个数量级，显然计算时间大幅缩短。
[0103]
本技术一些实施例中，在模拟的过程中，可以获得符合模拟要求的时间分辨率的运动序列，进而来模拟血液流动；对此，步骤104通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动，可以通过如下步骤实现：
[0104]
c31、获得当前计算时间处于一个周期中的特定时间，和多个时刻的运动序列的时间间隔。
[0105]
这里，特定时间为tt，多个时刻的运动序列的时间间隔为b。
[0106]
c32、获得特定时间处于多个时刻的运动序列的时间间隔的目标时间。
[0107]
这里，目标时间为ttt。
[0108]
c33、初始化当前网格节点。
[0109]
c34、若当前网格节点未更新，且目标时间未到预设时间，将当前网格节点对应的下一计算时间的网格节点坐标赋值给当前网格节点，并将当前网格节点标记为已更新状态。
[0110]
c35、若当前网格节点不是所有壁面网格节点中的最后一个节点，则查找当前网格节点的下一个网格节点，直至所有壁面网格节点均为已更新状态，以模拟心血管内的血液流动。
[0111]
需要说明的是，c34中将当前网格节点对应的下一计算时间的网格节点坐标赋值给当前网格节点之前，该方法还包括如下步骤：
[0112]
若当前网格节点未更新，且目标时间未到预设时间，初始化网格节点坐标矩阵行索引；
[0113]
获得当前网格节点坐标与当前计算时间网格节点坐标的差值；
[0114]
若差值与面网格最大尺寸之间满足预设数值关系，将当前网格节点对应的下一计算时间的网格节点坐标赋值给当前网格节点。
[0115]
需要说明的是，c34中将当前网格节点对应的下一计算时间的网格节点坐标赋值给当前网格节点之前，该方法还包括如下步骤：通过插值算法得到下一计算时间的网格节点坐标；其中，所有壁面网格节点均为已更新状态时，多个时刻的运动序列对应的时间分辨率符合模拟要求的时间分辨率。
[0116]
也就是说，在前述的多个时刻不足以满足cfd计算所需的时间分辨率即不符合模拟要求的时间分辨率的情况下，进行插值将多个时刻的点云扩展到符合cfd模拟的时间分辨率即符合模拟要求的时间分辨率。本技术在进行点云扩展的过程中，采用内部插值，并将插值的过程集成到计算模拟中的网格节点更新步骤中，不仅可以避免外部插值带来的巨大数据量并降低计算时间，同时无需预先指定计算时间步，其灵活度更高。可见，本技术的网格节点更新是在计算过程中，每一计算时间步的壁面网格节点是基于多个时刻的运动序列
来插值更新。本技术实施例中，内部插值算法包括但不限于如下至少之一：线性插值算法、三次样条插值算法、傅里叶插值算法。
[0117]
需要说明的是，相关技术中在进行插值的过程中，先在外部计算好每个cfd计算时间步的x_c、y_c、z_c，并在计算过程中的网格节点更新步骤中进行读取，该过程称为外部插值。这种外部插值的方式二维矩阵x_c、y_c、z_c的数据量特别大，占据较大的运行内存，增加计算时间和硬件要求。然而，本技术提出的内部插值方法，大幅减少了数组的读取和搜索，减少了计算要求的运行内存，最终减少了运行的计算时间。
[0118]
在一个可实现的采用并行处理方式实现模拟心血管内的血液流动的过程中，在模拟计算初始时刻网格节点被分配到各处理器中。这里，并行处理是指同时使用多个处理器解决计算问题的过程，可以进一步提高计算机系统计算速度和处理能力。
[0119]
在cfd计算初始时刻读入上述的多个时刻网格节点坐标的二维矩阵，即网格坐标矩阵x_c、y_c、z_c，并将网格坐标矩阵x_c、y_c、z_c储存在主节点host内，后续计算中网格坐标矩阵x_c、y_c、z_c由主节点host传入子节点node_n，其中node_n中n表示预设的并行子节点数。参照图7所述，并行处理下，网格节点更新在子节点中完成，更新过程如下：
[0120]
步骤s401：子节点node_n获取由主节点host传入的坐标矩阵x_c、y_c、z_c以及当前计算时间t和计算时间步长δt。
[0121]
步骤s402：计算当前计算时间处于一个周期中的时间以及所处时间段。
[0122]
其中，总的计算时间步(number of time steps)可以根据计算时间步长δt(time step size)和计算时间来确定。总的计算时间步＝计算时间/计算时间步长δt。
[0123]
本实施例中，一个周期为0.8s为例，图4中tt表示当前计算时间处于一个周期中的时间；b表示i个时刻的运动序列的时间间隔，本实施例b为0.04s；c表示tt/b取整；ttt表示tt处于i个时刻的运动序列的时间间隔的时间。
[0124]
步骤s403：初始化当前网格节点v＝1。
[0125]
步骤s404：判断当前网格节点v是否更新，如果当前网格节点v未更新，则进行步骤s4051；如果当前网格节点v已更新，则进行步骤s4052。
[0126]
步骤s4051：判断ttt是否为0，如果ttt不为0，则进行步骤s4061；如果为ttt为0，则进行步骤s4062。
[0127]
步骤s4061/步骤s4062：初始化网格节点坐标矩阵行索引ni＝0；
[0128]
步骤s4071：判断当前网格节点处于坐标矩阵x_c、y_c、z_c的行位置。
[0129]
如果nx&&ny&&nz《＝δs，则进行步骤s40811/步骤s40821，否则进行步骤s40812/40822。其中nx、ny及nz为当前网格节点坐标(node_x，node_y，node_z)与当前计算时间网格节点坐标(x[ni][tt]，y[ni][tt]，z[ni][tt])之差；δs为面网格最大尺寸，本实施例δs为0.05mm。应理解，当前计算时间网格节点坐标(x[ni][tt]，y[ni][tt]，z[ni][tt])是通过线性插值得到的。
[0130]
步骤s40811/步骤s40821：将下一计算时间的网格节点坐标赋值给当前网格节点。其中下一计算时间的网格节点坐标(x[ni][tt+δt]，y[ni][tt+δt]，z[ni][tt+δt])是通过线性插值得到的。
[0131]
步骤s409：将当前网格节点标记为已更新。
[0132]
步骤s410：判断当前网格节点是否为最后一个，如果是，则执行步骤s4111，否则执
行步骤s4112。
[0133]
在一个可实现的采用串行处理方式实现模拟心血管内的血液流动的过程中，相比于并行处理方式，串行处理方式可根据预设的网格节点顺序直接进行网格节点更新，无需循环判断当前网格节点处于坐标矩阵x_c、y_c、z_c的行位置。参照图8所述，串行处理下网格节点更新的过程如下：
[0134]
步骤s501：子节点node_j获取由主节点host传入的坐标矩阵x_c、y_c、z_c以及当前计算时间t和δt，其中j表示第j个子节点。
[0135]
步骤s502：计算当前计算时间处于一个周期中的时间以及所处时间段。
[0136]
本实施例中，一个周期为0.8s为例，图8中tt表示当前计算时间处于一个周期中的时间；b表示i个时刻的运动序列的时间间隔，本实施例b为0.04s；c表示tt/b取整；ttt表示tt处于i个时刻的运动序列的时间间隔的时间。
[0137]
步骤s503：初始化当前网格节点v＝1。
[0138]
步骤s504：判断当前网格节点v是否更新，如果当前网格节点v未更新，则进行步骤s5051；如果当前网格节点v已更新，则进行步骤s5052。
[0139]
步骤s5051：判断ttt是否为0，如果ttt不为0，则进行步骤s5061；如果为ttt为0，则进行步骤s5062。
[0140]
步骤s5061/步骤s5062：将下一计算时间的网格节点坐标赋值给当前网格节点。其中下一计算时间的网格节点坐标(x[v][tt+δt]，y[v][tt+δt]，z[v][tt+δt])是通过线性插值得到的。
[0141]
步骤s507：将当前网格节点标记为已更新。
[0142]
步骤s508：判断当前网格节点是否为最后一个，如果是，则执行步骤s5091，否则执行步骤s5092。
[0143]
需要说明的是，网格节点更新方法可根据实际需求灵活进行选择。示例性地，当模型网格数量较大时，可使用并行处理方式下网格节点更新的方法对网格节点进行更新；当模型网格数量较小时，可使用串行处理方式下网格节点更新的方法对网格节点进行更新。
[0144]
本技术的实施例提供一种基于图像的模拟血液流动的装置，图9所示的基于图像的模拟血液流动的装置600包括：
[0145]
获取模块601，用于获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、心血管的血液流动速度波形数据；
[0146]
处理模块602，用于对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标；
[0147]
处理模块602，用于对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列；
[0148]
模拟模块603，用于通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动
[0149]
本技术其他实施例中，处理模块602，用于选取多个时刻的网格节点坐标中，一个时刻的壁面网格节点坐标作为配准模板点云；
[0150]
选取多个时刻的网格节点坐标中，剩余时刻的壁面网格节点坐标作为配准目标点
云；
[0151]
将配准模板点云注册到每一配准目标点云上，得到配准模板点云与每一配准目标点云之间的变换矩阵；
[0152]
基于每一变换矩阵将配准模板点云进行相应的变换，得到一个周期内剩余时刻对应的运动序列。
[0153]
本技术其他实施例中，处理模块602，用于将配置模板点云中的点作为高斯混合模型的中心点；
[0154]
确定配准目标点云中的点符合独立同分布时的似然函数；
[0155]
基于高斯混合模型和似然函数，得到配准模板点云与每一配准目标点云之间的变换矩阵。
[0156]
本技术其他实施例中，配准模板点云中的点在配准目标点云中存在对应点。
[0157]
本技术其他实施例中，模拟模块603，用于获得当前计算时间处于一个周期中的特定时间，和多个时刻的时间间隔；
[0158]
获得特定时间处于多个时刻的运动序列的时间间隔的目标时间；
[0159]
初始化当前网格节点；
[0160]
若当前网格节点未更新，且目标时间未到预设时间，将当前网格节点对应的下一计算时间的网格节点坐标赋值给当前网格节点，将当前网格节点标记为已更新状态；
[0161]
若当前网格节点不是所有壁面网格节点中的最后一个节点，则查找当前网格节点的下一个网格节点，直至所有壁面网格节点均为已更新状态，以模拟心血管内的血液流动。
[0162]
本技术其他实施例中，模拟模块603，用于若当前网格节点未更新，且目标时间未到预设时间，初始化网格节点坐标矩阵行索引；
[0163]
获得当前网格节点坐标与当前计算时间网格节点坐标的差值；
[0164]
若差值与面网格最大尺寸之间满足预设数值关系，将当前网格节点对应的下一计算时间的网格节点坐标赋值给当前网格节点。
[0165]
本技术其他实施例中，模拟模块603，用于通过插值算法得到下一计算时间的网格节点坐标；其中，所有壁面网格节点均为已更新状态时，多个时刻的运动序列对应的时间分辨率符合模拟要求的时间分辨率。
[0166]
本技术实施例提供的基于图像的模拟血液流动的装置，通过获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、心血管的血液流动速度波形数据；对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标；对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列；通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；可见，本技术在模拟血液流动的过程中，获得一个周期内多个时刻的心血管三维模型，无需考虑血管材料属性，对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标，进而对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列后，便可以通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；避免了流固耦合方法需要提供心血管如左心房及左心耳壁面的材料属性才能处理的问题，大幅度缩短数据前处理时间和计算时间。
[0167]
本技术的实施例提供一种基于图像的模拟血液流动的设备，图10所示的基于图像的模拟血液流动的设备700包括：处理器701、存储器702和通信总线703，其中：
[0168]
通信总线703用于实现处理器701和存储器702之间的通信连接。
[0169]
处理器701用于执行存储器702中存储的基于图像的模拟血液流动的程序，以实现如下步骤：
[0170]
获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、心血管的血液流动速度波形数据；
[0171]
对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标；
[0172]
对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列；
[0173]
通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动。
[0174]
本技术实施例中，处理器701用于执行存储器702中存储的基于图像的模拟血液流动的程序，以实现如下步骤：
[0175]
选取多个时刻的网格节点坐标中，一个时刻的壁面网格节点坐标作为配准模板点云；
[0176]
选取多个时刻的网格节点坐标中，剩余时刻的壁面网格节点坐标作为配准目标点云；
[0177]
将配准模板点云注册到每一配准目标点云上，得到配准模板点云与每一配准目标点云之间的变换矩阵；
[0178]
基于每一变换矩阵将配准模板点云进行相应的变换，得到一个周期内剩余时刻对应的运动序列。
[0179]
本技术实施例中，处理器701用于执行存储器702中存储的基于图像的模拟血液流动的程序，以实现如下步骤：
[0180]
将配置模板点云中的点作为高斯混合模型的中心点；
[0181]
确定配准目标点云中的点符合独立同分布时的似然函数；
[0182]
基于高斯混合模型和似然函数，得到配准模板点云与每一配准目标点云之间的变换矩阵。
[0183]
本技术实施例中，配准模板点云中的点在配准目标点云中存在对应点。
[0184]
本技术实施例中，处理器701用于执行存储器702中存储的基于图像的模拟血液流动的程序，以实现如下步骤：
[0185]
获得当前计算时间处于一个周期中的特定时间，和多个时刻的运动序列的时间间隔；
[0186]
获得特定时间处于多个时刻的运动序列的时间间隔的目标时间；
[0187]
初始化当前网格节点；
[0188]
若当前网格节点未更新，且目标时间未到预设时间，将当前网格节点对应的下一计算时间的网格节点坐标赋值给当前网格节点，将当前网格节点标记为已更新状态；
[0189]
若当前网格节点不是所有壁面网格节点中的最后一个节点，则查找当前网格节点的下一网格节点，直至所有壁面网格节点均为已更新状态，以模拟心血管内的血液流动。
[0190]
本技术实施例中，处理器701用于执行存储器702中存储的基于图像的模拟血液流
动的程序，以实现如下步骤：
[0191]
若当前网格节点未更新，且目标时间未到预设时间，初始化网格节点坐标矩阵行索引；
[0192]
获得当前网格节点坐标与当前计算时间网格节点坐标的差值；
[0193]
若差值与面网格最大尺寸之间满足预设数值关系，将当前网格节点对应的下一计算时间的网格节点坐标赋值给当前网格节点。
[0194]
本技术实施例中，处理器701用于执行存储器702中存储的基于图像的模拟血液流动的程序，以实现如下步骤：
[0195]
通过插值算法得到下一计算时间的网格节点坐标；其中，所有壁面网格节点均为已更新状态时，多个时刻的运动序列对应的时间分辨率符合模拟要求的时间分辨率。
[0196]
本技术实施例提供的基于图像的模拟血液流动的设备，通过获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、心血管的血液流动速度波形数据；对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标；对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列；通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；可见，本技术在模拟血液流动的过程中，获得一个周期内多个时刻的心血管三维模型，无需考虑血管材料属性，对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标，进而对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列后，便可以通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；避免了流固耦合方法需要提供心血管如左心房及左心耳壁面的材料属性才能处理的问题，大幅度缩短数据前处理时间和计算时间。
[0197]
需要说明的是，本实施例中处理器所执行的步骤的具体实现过程，可以参照图1对应的实施例提供的基于图像的模拟血液流动的方法中的实现过程，此处不再赘述。
[0198]
本技术实施例装置的描述，与上述方法实施例的描述是类似的，具有同方法实施例相似的有益效果，因此不做赘述。对于本装置实施例中未披露的技术细节，请参照本技术方法实施例的描述而理解。
[0199]
本技术实施例提供一种存储有可执行指令的存储介质，其中存储有可执行指令，当可执行指令被处理器执行时，将引起处理器执行本技术实施例提供的方法，例如，如图1示出的方法。
[0200]
本技术实施例提供的存储介质，通过获得对医学图像中的心血管进行重构处理，得到的一个周期内多个时刻的三维模型、心血管的血液流动速度波形数据；对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标；对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列；通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；可见，本技术在模拟血液流动的过程中，获得一个周期内多个时刻的心血管三维模型，无需考虑血管材料属性，对多个时刻的三维模型进行网格离散，得到多个时刻的网格节点坐标，进而对多个时刻的网格节点坐标进行点云配准，得到多个时刻的运动序列后，便可以通过计算流体力学算法基于血液流动速度波形数据以及多个时刻的运动
序列，对每一计算时间步的壁面网格节点进行网格节点更新处理，以模拟心血管内的血液流动；避免了流固耦合方法需要提供心血管如左心房及左心耳壁面的材料属性才能处理的问题，大幅度缩短数据前处理时间和计算时间。
[0201]
在一些实施例中，存储介质可以是计算机可读存储介质，例如，铁电存储器(fram，ferromagnetic random access memory)、只读存储器(rom，read only memory)、可编程只读存储器(prom，programmable read only memory)、可擦除可编程只读存储器(eprom，erasable programmable read only memory)、带电可擦可编程只读存储器(eeprom，electrically erasable programmable read only memory)、闪存、磁表面存储器、光盘、或光盘只读存储器(cd-rom，compact disk-read only memory)等存储器；也可以是包括上述存储器之一或任意组合的各种设备。
[0202]
在一些实施例中，可执行指令可以采用程序、软件、软件模块、脚本或代码的形式，按任意形式的编程语言(包括编译或解释语言，或者声明性或过程性语言)来编写，并且其可按任意形式部署，包括被部署为独立的程序或者被部署为模块、组件、子例程或者适合在计算环境中使用的其它单元。
[0203]
作为示例，可执行指令可以但不一定对应于文件系统中的文件，可以可被存储在保存其它程序或数据的文件的一部分，例如，存储在超文本标记语言(超文本标记语言，hyper text markup language)文档中的一个或多个脚本中，存储在专用于所讨论的程序的单个文件中，或者，存储在多个协同文件(例如，存储一个或多个模块、子程序或代码部分的文件)中。作为示例，可执行指令可被部署为在一个计算设备上执行，或者在位于一个地点的多个计算设备上执行，又或者，在分布在多个地点且通过通信网络互连的多个计算设备上执行。
[0204]
以上所述，仅为本技术的实施例而已，并非用于限定本技术的保护范围。凡在本技术的精神和范围之内所作的任何修改、等同替换和改进等，均包含在本技术的保护范围之内。