CN112070899B - 基于医学图像还原血管零压模型的方法、装置及存储介质 - Google Patents

Abstract

本发明公开了一种基于医学图像还原血管零压模型的方法、装置及存储介质。具体包括以下步骤：S1：根据医学图像提取原始血管模型，对原始血管模型进行网格划分获得网格节点编号集合及网格节点坐标集合；S2：根据原始血管模型获取初始模型的点坐标集合；S3：以步骤S2得到的初始模型的点坐标集合作为迭代的初始模型进行循环迭代处理，直至满足预设的迭代停止条件时停止进行迭代处理；并输出当前血管模型的点坐标集合作为血管零压模型的点坐标集合；S4：从血管零压模型的点坐标集合中获取血管零压状态下的网格模型，以还原零压状态下血管的形状。本发明相比传统分析方法能更真实还原血管零压模型，同时不会引入过量的形变造成误差。

Description

基于医学图像还原血管零压模型的方法、装置及存储介质

技术领域

本发明属于生物医学工程领域，具体涉及一种基于医学图像还原血管零压模型的方法。

背景技术

人体内部的血管在血压的作用下一直处于膨胀的状态，医学图像中提取出来的血管模型也已经处于膨胀的形态，血管壁内部存在一定的应力分布。对血管的医学图像进行结构力学分析时，理想的分析方法应该是首先获取血管内表面没有载荷且血管壁内部应力为0的状态下的模型，称之为零压模型或零应力模型，再用零压模型进行结构力学分析。然而常见的分析方法是忽略已经发生的膨胀，认为从医学图像中提取的原始血管模型处于零应力状态；然后直接在原始医学图像模型内部施加压力边界，如人体舒张压或收缩压，通过有限元算法或软件模拟血管在血压作用下的受力状态。这种简化的分析方法最终会导致产生过大的形变，分析结果的准确性也会下降。

发明内容

为了从血管的结构力学分析中获得更准确的计算结果，本发明提出了一种基于医学图像还原血管零压模型的方法，该方法以动脉血压为边界条件，结合正向的有限元模拟和逆向的坐标还原，经过多次迭代后逐步将医学图像中提取的原始血管模型还原为零压模型。

本发明采用的技术方案如下：

一、一种基于医学图像还原血管零压模型的方法，包括以下步骤：

步骤S1：根据医学图像重构得到原始血管模型，对原始血管模型进行网格划分获得网格节点集合N_ref及网格节点坐标集合X_ref；

步骤S2：根据原始血管模型获取初始模型的点坐标集合X₀；

步骤S3：以步骤S2得到的初始模型的点坐标集合X₀作为迭代的初始模型进行循环迭代处理，直至满足预设的迭代停止条件时停止进行迭代处理；并输出当前血管模型的点坐标集合X_j作为血管零压模型的点坐标集合X_zero；其中，j为循环次数，j的初始值为0；

步骤S4：由于网格单元各节点之间的连接关系不变，能从血管零压模型的点坐标集合X_zero中获取血管零压状态下的网格模型，以还原零压状态下血管的形状。

所述步骤S3具体为：

3.1)对当前血管模型的点坐标集合X_j进行有限元分析，得到形变后血管模型的点坐标集合迭代记数j表示当前为第j次迭代；

3.2)根据和X_ref计算残差矢量集R_j，根据本次迭代的残差矢量集R_j和上一次迭代的残差矢量集R_j-1，对每个节点i计算第一参数β_i：

判断第一参数β_i的最大值是否小于0.0001，如果是，则在完成步骤3.2)后停止计算，直接输出当前血管模型的点坐标集合X_j作为血管零压模型的点坐标集合X_zero；

否则，执行步骤3.3)；

3.3)根据第一参数β_i计算第二参数α_i：

当j＝0时，α_i等于1；

当j＞0时，若β_i的最大值小于1，则α_i＝β_i；若β_i的最大值大于1，则用β_i的最大值对β_i作归一化并得到第二参数α_i：

3.4)根据第二参数α_i和残差矢量集R_j计算用于下一次迭代的点坐标集合：

其中，X_j为当前迭代步的血管模型点坐标集合，X_j,i为第i个网格节点在当前迭代步的血管模型中的坐标；

3.5)更新迭代记数以及点坐标集合：

X_j＝X_j+1,j＝j+1

更新后返回步骤3.1)进行下一次循环。

所述步骤3.2)具体为：

3.2.1)根据和X_ref计算残差矢量集R_j：

为形变后血管模型中第i个网格节点的点坐标，X_ref,i为原始血管模型中第i个网格节点的点坐标，i为网格节点集合N_ref中的一个网格节点；

3.2.2)根据本次迭代的残差矢量集R_j和上一次迭代的残差矢量集R_j-1，对每个节点i计算第一参数β_i：

当j＝0时，不计算第一参数β_i；

当j＞0时，对每个节点i计算第一参数β_i：若残差矢量为0，则β_i＝0；反之则用以下表达式计算β_i：

其中，R_j-1,i为第i个网格节点上一次迭代的残差矢量，R_j,i为第i个网格节点本次迭代的残差矢量，θ_i是两个残差矢量之间的夹角。

施加载荷p(t)，对原始血管模型的点坐标集合X_ref进行有限元分析得到变形后血管模型的点坐标集合变形后的点坐标集合/>随时间t变化；根据变形后血管模型的点坐标集合/>和原始血管模型的点坐标集合X_ref得到每个节点的位移矢量集D(t)；位移矢量集D(t)进行x,y,z三个方向的分解后得到位移分量集合D_x(t)，D_y(t)和D_z(t)；对每个网格节点，进行位移分量与载荷之间的二阶拟合，得到以下关系式：

其中，a_k,i，b_k,i和c_k,i是二阶拟合得到的三个系数；D_k,i(t)表示第i个网格节点在x,y,z三个方向上关于载荷的二阶回归；

从x,y,z三个方向的回归结果中取各自的截距c_k,j组成截距矢量集C：将C的元素与X_ref的元素相加，得到以下初始模型的点坐标集合X₀：

其中，C_i表示第i个网格节点的截距矢量，X_ref,i表示原始血管模型中第i个网格节点的坐标。

所述的有限元分析包括以下步骤：

首先对血管模型进行网格划分，将目标几何结构离散为若干个网格单元；然后求出各网格单元的刚度矩阵；确定边界条件和载荷，构建物理问题；最后通过方程求解得出各网格单元的受力和位移。

所述网格划分的具体内容为：对血管壁表面，以三角形，四边形等简单平面几何体为网格单元，使用网格单元组合出血管壁的形状；对血管壁，以四面体，六面体等简单立体几何体为网格单元，使用网格单元填充血管壁内表面(与血液接触的一侧)和外表面(与器官组织接触的一侧)包裹的空间域。最终所有网格单元组成了血管壁模型的网格，用于有限元分析。

所述的有限元分析采用但不限于Abaqus软件或Ansys软件。

所述的网格划分采用但不限于开源的Gmsh或Hypermesh软件。

所述的有限元分析求解的物理问题为：

记血管模型所处的三维空间域为Ω，在有限元计算中，首先需要对空间域Ω进行网格划分，同时血管的内表面(与血液接触的表面)记为Γ_N边界，血管出口和入口的横截面记为Γ_D边界并限制位移为0：

其中，P和S分别为任意网格节点的第一和第二皮奥拉-基尔霍夫应力；

u为Γ_D边界上任意网格节点的位移，u_D为Γ_D边界上的位移值；

N_Γ为原始坐标系下Γ_N边界上任意网格节点的法向量，法向量方向朝内指向血液；

t₀(t)为在Γ_N边界上任意网格节点随时间变化的受力，方向与Γ_N的法向量相反，Γ_N上所有网格节点的受力大小的和等于血管内表面载荷与血管内表面的积，即：

其中，A为内表面的面积，p(t)为载荷，p(t)随时间t逐渐增加，最大值为获取医学图像的时刻下的血管压强p_artery；

F是变形梯度张量：x和x_def分别为原始血管模型中的任意一个网格节点的坐标和该网格节点在形变后血管模型中的坐标。

所述步骤4中，获得血管零压模型的点坐标集合X_zero的同时输出血管零压模型的应力结果，应力结果包括Von Mises应力。

二、采用上述方法的一种基于医学图像还原血管零压模型的装置

包括三维重构模块、有限元分析模块、零压模型迭代模块和后处理模块；

三维重构模块，用于对医学图像进行重构后获取三维的血管模型；

有限元分析模块，用于对血管模型进行网格划分，确定边界条件和载荷后，通过有限元分析获取各网格单元的受力和位移；

零压模型迭代模块，以初始模型开始循环迭代获取血管零压模型，每次循环迭代时通过有限元分析模块获得新的迭代模型；

后处理模块，用于获取血管零压模型的应力结果。

三、一种存储介质

包括处理器和存储器，存储器用于存储程序，程序为对应上述基于医学图像还原血管零压模型方法的指令。

处理器通过运行程序以执行上述基于医学图像还原血管零压模型的方法。

本发明的有益效果：

1)本发明可以用于还原主动脉(包括升主动脉、主动脉弓、弓上动脉和降主动脉)的零压模型，而现有的其他迭代方法在用于重构主动脉(包括升主动脉、主动脉弓、弓上动脉和降主动脉)零压模型的时候，无法得到合理的形状。零压模型相对原始模型有一定收缩，而现有方法会使主动脉模型过度收缩，主动脉表面发生交错如图5；又或者，如图6所示，因为网格节点不合理的位移，导致产生网格质量差的网格(圈内区域因为过度扭曲出现负网格)，而且网格表面的连续性不好，出现凹凸不平的情况。

2)本发明改进了传统分析方法，能合理地评估血管在人体内受力状态，同时不会引入过量的形变造成误差；且整体易操作。

3)本发明采用向前式迭代，在每次迭代过程中保持网格单元各节点的连接关系，不需要重新进行前处理，节省时间。

4)循环迭代流程中，α_i的最大值被限制为1，这样的处理能防止还原过程中网格节点过度位移，产生交错的网格或者过度扭曲的网格，保证计算过程的稳定性。

附图说明

图1所示为本发明的整体流程图。

图2所示为本发明循环迭代过程的流程图。

图3所示为实施例中从医学图像提取的主动脉模型的各组成部分以及进行前处理后获得的网格。

图4所示为实施例中在血管模型上设置的Γ_D边界和Γ_N边界。

图5所示为采用现有方法重构零压模型使主动脉模型过度收缩的结果图。

图6所示为采用现有方法重构零压模型导致网格过度扭曲的结果图。

图7所示为实施例的初始化步骤中，网格中每个节点进行三个方向位移d_x，d_y，d_z关于载荷的二阶回归。

图8所示为β_i绝对值的最大值在迭代过程中的变化曲线(因为第0次循环不计算β_i，所以曲线从j＝1开始)。

图9所示为本发明与原始模型之间的效果显示图。

图10所示为实施例中对血管零压模型进行应力分析的结果图。

图11所示为实施例对应的还原血管零压模型的装置。

具体实施方式

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

本实施例以主动脉为例，阐述如图1所示的还原零压模型的具体过程。

1)从医学图像中获取原始血管模型Ω并进行前处理。

医学图像可以选择CT图像或4D FLOW图像等，为了方便确认血管模型对应的心动周期时刻，本实施例中采用4D FLOW图像，重构出血管模型并记录下该心动周期时刻的主动脉压p_artery；本实施例中采用的血管模型各组成部分分别是：升主动脉，主动脉弓，无名动脉，左颈动脉，左锁骨下动脉和降主动脉；对主动脉模型进行网格划分得到节点集合X_ref；本实施例中使用的二维网格为三角形，三维网格为四面体。如图3所示展示了提取得到的主动脉模型以及进行网格划分后的升主动脉部分血管。

2)进行坐标还原的第一步，初始化。

首先对X_ref构建有限元分析的物理问题：

如图4所示，把Γ_D边界设为血管模型出口和入口的横截面，把Γ_N边界设为血管模型的内表面；在Γ_D边界施加Dirichlet边界条件，约束其位移为0；在Γ_N边界施加Neumann边界条件，即随时间线性增长的载荷p(t)，载荷在经过一个心动周期后达到最大值p_artery。

求解物理问题并得到网格节点随载荷增加而变化的位移矢量集合D(t)以及在x,y,z三个方向位移分量集合D_x(t)，D_y(t)和D_z(t)。对每个节点，都进行位移分量关于载荷变化的二阶回归，回归结果如图7所示。本实施例中，对载荷进行了归一化处理，因此横坐标范围为[0,1]。从回归结果中得到截距矢量的集合C，将C的元素与X_ref的元素相加，得到用作迭代初始模型的节点集合X₀。

3)循环迭代。

如图2所示，输入前两个步骤的数据X₀，X_ref和p(t)，把循环次数标记j初始化为0，由于网格单元各节点之间的连接关系不变，直接从X₀获取对应的初始网格。循环从j＝0开始，每次循环都根据第(2)步的方法对X_j建立物理问题，求解后得到血管形变后的坐标集合对/>和X_ref作差获取残差矢量集R_j。j＝0时，每个节点i的参数α_i等于1，根据α_i和X₀直接更新用于下一次循环的网格和循环次数标记j，j＝1且对应的点坐标集为X₁。从j＝1开始，在每次循环内求解有限元问题后得到R_j，结合上一次循环的残差矢量集R_j-1计算每个节点i的参数β_i和α_i，根据α_i和R_j更新下一次循环的点坐标集以及对应的网格。从j＝1开始，监视β_i绝对值的最大值，直到其小于0.0001。本实施例中，在j＝12时符合该条件，因此零压模型的点坐标集合X_zero＝X₁₂，β_i绝对值的最大值的变化情况如图8所示。

4)效果显示和应力分布显示。

本实施例中，为了显示本方法的合理性，在图9中对比了三种模型的轮廓线，分别为：血管零压模型膨胀后的轮廓，原始血管模型的轮廓和原始血管模型膨胀后的轮廓；血管零压模型膨胀后的轮廓与原始血管模型的轮廓接近，而对原始血管模型施加载荷后则发生过度膨胀。

确认血管零压模型后，可以同时输出应力结果，例如图10中的Von Mises应力分布，完成对医学图像的结构力学分析。

如图11所示，本实施例提供了一种基于医学图像还原血管零压模型的处理装置，包括以下模块：1、基于医学图像的主动脉三维几何模型重构模块；2、有限元分析模块；3、零压模型迭代模块；4、计算结果后处理模块。

其中模块1用于获取主动脉的医学图像，如4Dflow图像和CTA图像，根据图像信息构建主动脉模型，包括升主动脉，主动脉弓，降主动脉，无名动脉，左颈动脉，左锁骨下动脉等。模块2用于计算主动脉模型的受力和形变，根据模块1提供的几何模型进行网格划分，基于网格求解有限元问题，得到主动脉的受力和形变。模块3用于还原零压模型，根据模块2的形变结果，计算出新的主动脉模型以及对应的网格文件；然后，模块3基于计算新模型过程中的中间数据，选择进行以下操作：1)中间数据提示继续迭代计算，新模型输入到模块2，计算新的形变结果；2)中间数据提示中止计算，认为当前的主动脉模型作为零压模型，输出对应网格文件。模块4用于显示主动脉模型和受力状态，根据模块3输出的网格文件和模块2计算的受力结果，展现主动脉在血流载荷下的受力情况。

此外，本实施例提供了一种存储介质，包括处理器和存储器，存储器用于存储程序，程序对应本实施例中基于医学图像还原血管零压模型方法的指令，通过运行处理器执行程序实现血管零压模型的还原。

Claims (6)

1.一种基于医学图像还原血管零压模型的方法，其特征在于，包括以下步骤：

步骤S1：根据医学图像提取原始血管模型，对原始血管模型进行网格划分获得网格节点集合N_ref及网格节点坐标集合X_ref；

步骤S2：根据原始血管模型获取初始模型的点坐标集合X₀；

步骤S3：以步骤S2得到的初始模型的点坐标集合X₀作为迭代的初始模型进行循环迭代处理，直至满足预设的迭代停止条件时停止进行迭代处理；并输出当前血管模型的点坐标集合X_j作为血管零压模型的点坐标集合X_zero；其中，j为循环次数，j的初始值为0；

步骤S4：从血管零压模型的点坐标集合X_zero中获取血管零压状态下的网格模型，以还原零压状态下血管的形状；

所述步骤S2具体为：

施加载荷p(t)，对原始血管模型的点坐标集合X_ref进行有限元分析得到变形后血管模型的点坐标集合变形后的点坐标集合/>随时间t变化；根据变形后血管模型的点坐标集合/>和原始血管模型的点坐标集合X_ref得到每个节点的位移矢量集D(t)；位移矢量集D(t)进行x,y,z三个方向的分解后得到位移分量集合D_x(t)，D_y(t)和D_z(t)；对每个网格节点，进行位移分量与载荷之间的二阶拟合，得到以下关系式：

D_k,i(t)＝a_k,ip(t)²+b_k,ip(t)+c_k,i,k＝x,y,z

其中，a_k,i，b_k,i和c_k,i是二阶拟合得到的三个系数；D_k,i(t)表示第i个网格节点在x,y,z三个方向上关于载荷的二阶回归；

从x,y,z三个方向的回归结果中取各自的截距c_k,j组成截距矢量集C：将C的元素与X_ref的元素相加，得到以下初始模型的点坐标集合X₀：

其中，C_i表示第i个网格节点的截距矢量，X_ref,i表示原始血管模型中第i个网格节点的坐标；

所述步骤S3具体为：

3.1)对当前血管模型的点坐标集合X_j进行有限元分析，得到形变后血管模型的点坐标集合迭代记数j表示当前为第j次迭代；

3.2)根据和X_Ω计算残差矢量集R_j，根据本次迭代的残差矢量集R_j和上一次迭代的残差矢量集R_j-1，对每个节点i计算第一参数β_i：

判断第一参数β_i的最大值是否小于0.0001，如果是，则在完成步骤3.2)后停止计算，直接输出当前血管模型的点坐标集合X_j作为血管零压模型的点坐标集合X_zero；

否则，执行步骤3.3)；

3.3)根据第一参数β_i计算第二参数α_i：

当j＝0时，α_i等于1；

当j>0时，若β_i的最大值小于1，则α_i＝β_i；若β_i的最大值大于1，则用β_i的最大值对β_i作归一化并得到第二参数α_i：

3.4)根据第二参数α_i和残差矢量集R_j计算用于下一次迭代的点坐标集合：

其中，X_j为当前迭代步的血管模型点坐标集合，X_j,i为第i个网格节点在当前迭代步的血管模型中的坐标；

3.5)更新迭代记数以及点坐标集合：

X_j＝X_j+1,j＝j+1

更新后返回步骤3.1)进行下一次循环；

所述的有限元分析求解的物理问题为：

记血管模型所处的三维空间域为Ω，在有限元计算中，首先需要对空间域Ω进行网格划分，同时血管的内表面记为Γ_N边界，血管出口和入口的横截面记为Γ_D边界并限制位移为0：

其中，P和S分别为任意网格节点的第一和第二皮奥拉-基尔霍夫应力；

u为Γ_D边界上任意网格节点的位移，u_D为Γ_D边界上网格节点的位移值；

N_Γ为Γ_N边界上任意网格节点的法向量；

t₀(t)为在Γ_N边界的任意网格节点上随时间变化的受力，方向与Γ_N的法向量相反，Γ_N上所有网格节点的受力大小的和等于血管内表面载荷与血管内表面的积，即：

其中，A为内表面的面积，p(t)为载荷，p(t)随时间t逐渐增加，最大值为获取医学图像的时刻下的血管压强p_artery；

F是变形梯度张量：x和x_def分别为原始血管模型中的任意一个网格节点的坐标和该网格节点在形变后血管模型中的坐标。

2.根据权利要求1所述的一种基于医学图像还原血管零压模型的方法，其特征在于，所述步骤3.2)具体为：

3.2.1)根据和X_ref计算残差矢量集R_j：

为形变后血管模型中第i个网格节点的点坐标，X_ref,i为原始血管模型中第i个网格节点的点坐标，i为网格节点集合N_ref中的一个网格节点编号；

为形变后血管模型的点坐标集合；

3.2.2)根据本次迭代的残差矢量集R_j和上一次迭代的残差矢量集R_j-1，对每个节点i计算第一参数β_i：

当j＝0时，不计算第一参数β_i；

当j>0时，对每个节点i计算第一参数β_i：若残差矢量为0，则β_i＝0；反之则用以下表达式计算β_i：

θ_i＝∠(R_j-1,i,R_j,i)

其中，R_j-1,i为第i个网格节点上一次迭代的残差矢量，R_j,i为第i个网格节点本次迭代的残差矢量，θ_i是两个残差矢量之间的夹角。

3.根据权利要求1所述的一种基于医学图像还原血管零压模型的方法，其特征在于，所述步骤4中，获得血管零压模型的点坐标集合X_zero的同时输出血管零压模型的应力结果，应力结果包括Von Mises应力。

4.采用权利要求1～3任一所述方法的一种基于医学图像还原血管零压模型的装置，其特征在于，包括三维重构模块、有限元分析模块、零压模型迭代模块和后处理模块；

三维重构模块，用于对医学图像进行重构后获取三维的血管模型；

有限元分析模块，用于对血管模型进行网格划分，确定边界条件和载荷后，通过有限元分析获取各网格单元的受力和位移；

零压模型迭代模块，以初始模型开始循环迭代获取血管零压模型，每次循环迭代时通过有限元分析模块获得新的迭代模型；

后处理模块，用于获取血管零压模型的应力结果。

5.一种存储介质，其特征在于，包括处理器和存储器，存储器用于存储程序，程序为对应权利要求1～4中任一所述的基于医学图像还原血管零压模型方法的指令。

6.根据权利要求5所述的一种存储介质，其特征在于，所述处理器通过运行存储器的程序以执行权利要求1～4中任一所述的基于医学图像还原血管零压模型的方法。