技术领域
[0001] 本
发明属于
生物医学图像重建技术领域，具体涉及一种基于GTV(Graph Total Variation，图全变分)的心脏左心室运动分析方法。
背景技术
[0002] 在现代社会，心脏类
疾病已经日益成为人类健康的最主要威胁，显著影响人类的生命
质量；全球心
血管疾病(CVD)是导致死亡的主要原因之一，预计将在接下来继续威胁人类健康，心肌缺血或冠心病可以通过分析心脏的运动和
变形来识别和
定位。为了定量地识别和测量感兴趣区域中的病变的
位置和程度，近年来，如何从
图像序列中获得准确和稳健的心肌
运动估计已经成为挑战。
[0003] 左心室作为心脏最主要的组成部分，是躯体供
氧的动
力枢纽，是心脏影像学的主要研究对象；左心室心肌的运动情况为心脏疾病的诊断提供依据，因此成为当前研究的热点。对心脏左心室内外膜之间的心肌的定点
跟踪和左心室心肌的应变分析及其直观描述的研究尤为重要，对于心脏疾病的临床诊断具有重要的参考价值；从近几年的研究成果及文献来看，基于心脏电影MR图像的左心室运动分析的研究引起越来越多的研究机构的关注。
[0004] 由于心脏是一个非刚体材料，在心肌运动分析中应变通常作为要估计的运动参数，但由于应变是位移场导数的函数，对噪声敏感，因此如何获得准确并且完整的心肌位移场具有重要的生理和临床意义。一般策略是将从心脏的
动态图像序列中导出的稀疏和噪声数据(中间壁速度或边界位移)作为先验知识，例如通用的图像标记方法、MR的
相位对比速度等，但由于导出的特征点的不完整性和噪声干扰性，从稀疏地标位移场中恢复成密集位移场估计是一个病态问题，需要与相应的平滑和插值运动模型相结合，作为求解唯一性的约束。
[0005] 在现有的左心室运动重建模型中，这种约束模型可以是数学正则化法也可以是
连续介质力学最小化法。正则化方法历来被应用于
计算机视觉中的许多问题，当正则化模型应用于非刚体运动估计时，可以平稳地恢复运动场的连续性和平滑性，如卡尔曼
滤波器法、Tikhonov正则化法、全变分方法等；而有限元方法和无网格法常常应用在心肌几何结构中作为数值分析技术来解决工程中的近似问题。
[0006] 然而在对不准确、不完整数据的恢复过程中，由于相邻
节点之间的节点值的变化非常小并且对噪声敏感，所以先前的正则化方法可能会引发解的过平滑、丢失边缘或信息锐化等问题；而连续介质力学最小化法不仅仅强依赖于解剖学和图像学的参数，还需要手动调节参数来平衡噪声误差。在基于几何学的微分求解过程中，有限元单元法相比无网格法，在物理力学模型上具有更有效、稳定和精确的构造方法，但在边界处的连续性和光滑性也会受到邻点位移值约束的影响。
发明内容
[0007] 鉴于上述，本发明提供了一种基于GTV的心脏左心室运动分析方法，提出一种将图全变分方法与生物力学模型相结合的有限元
框架，利用KNN
算法对图的位移信息进行相似度分类、平滑，在求解步骤中应用前向-后向对偶算法，将原优化问题转变为联合求解可微分问题和不可微分问题的最小化问题，调整求解目标变量变化的可变步长参数，
加速了目标函数的收敛速度，提高了
分辨率。
[0008] 一种基于GTV的心脏左心室运动分析方法，包括如下步骤：
[0009] (1)对左心室MRI图像序列进行动态分割，获取心肌表面(心内膜和心外膜)轮廓信息；
[0010] (2)构造用于重建心肌运动信息的物理模型；
[0011] (3)根据所述物理模型利用有限元数值分析方法建立心肌的生物力学模型；
[0012] (4)将图全变分的正则化约束与心肌的生物力学模型结合，构造出左心室运动分析重建模型；
[0013] (5)对上述重建模型进行优化求解，得到心脏左心室的运动信息。
[0014] 进一步地，所述步骤(2)中由于所获得的轮廓信息是一系列有噪声且稀疏的边界位移点，需从中恢复出鲁棒且稠密的位移场，故采用连续介质力学方法利用生物解剖得到的材料参数作为力学的先验知识，模拟出心肌的非刚性变形，从而构建得到用于重建心肌运动信息的物理模型。
[0015] 进一步地，所述步骤(3)中考虑到心肌的非刚体变形，将心肌全局节点位移与全局节点力的关系由以下
刚度方程联立起来，即构建得到心肌的生物力学模型；
[0016] F＝KU
[0017] 其中：K为心肌的刚度矩阵且为n×n大小，U为心肌全局节点的位移矩阵且为n×1大小，F为心肌全局节点的力矩阵且为n×1大小，n为心肌全局节点(包括轮廓边界位移点以及心肌内部位移点)的数量。
[0018] 进一步地，所述步骤(4)中的左心室运动分析重建模型表达式如下：
[0019]
[0020] 其中：K为心肌的刚度矩阵且为n×n大小，U为心肌全局节点的位移矩阵且为n×1大小，F为心肌全局节点的力矩阵且为n×1大小，n为心肌全局节点的数量，μ为正则化参数，||||2表示2范数，GTV(U)为关于U的图全变分算子。
[0021] 进一步地，所述步骤(4)中图全变分的概念是将量化原始位移的图
信号和滤波后的位移图信号之间的差作为图形变化项进行平滑。
[0022] 进一步地，所述图全变分算子GTV(U)的表达式如下：
[0023]
[0024]
[0025] 其中：i为任一心肌全局节点，ηi为与节点i最邻近的节点集合，j为节点集合ηi中的任一心肌全局节点且i≠j，ui和uj分别为位移矩阵U中关于节点i和j的位移量，||||1表示1范数，ai为节点集合ηi中所有点与节点i的平均距离。
[0026] 进一步地，所述点集合ηi通过KNN算法(the K-nearest neighbor algorithm，K最邻近算法)确定。
[0027] 进一步地，所述步骤(5)中采用前向-后向算法(Forward-Backward Algorithm)结合原始-对偶方法(primal-dual method)对重建模型进行优化求解，以解决特殊凸优化问题。
[0028] 本发明通过生物力学模型和有限元几何分析法重建心脏恢复框架，利用图全变分方法引入图的信号概念，建立平滑模型以重构来自噪声样本的信号；图
信号处理中可利用KNN算法对图的位移信息进行相似度分类、平滑，以直观、灵活的方式获得数据中的相似关系，集合有限元和无网格法的优势，使信号在无约束的空间内平稳地减少噪声。在此过程中，图全变分将找出该距离内所有节点之间总变化最小的图形信号，优于其它正则化方法，噪声直接并入计算进行平滑去噪。在求解步骤中本发明应用前向-后向对偶算法，将原优化问题转变为联合求解可微分问题和不可微分问题的最小化问题，调整求解目标变量变化的可变步长参数，加速了目标函数的收敛速度，提高了分辨率。
附图说明
[0029] 图1为本发明基于GTV心脏左心室运动分析方法的步骤流程示意图。
具体实施方式
[0030] 为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。
[0031] 如图1所示，本发明基于GTV的心脏左心室运动分析方法，包括如下步骤：
[0032] (1)利用电影相位
对比度梯度回波序列采集具有心肌缺血症状的成年杂种犬的16
帧磁共振图像。
[0033] (2)对左室心内外膜的边缘轮廓进行分割，并选择合适的
采样点建立整个心肌的边界模型。本实施方式采用一种简单的区域填充方法，得到了心肌边界和内部的所有采样点，并利用这些采样点建立了有限元的力学分析模型。
[0034] 考虑到心肌是一个持续的不可压缩材料，考虑
惯性力和阻尼力，系统的动力学方程如下，其中F是全局位移，M质量矩阵，C阻尼矩阵；阻尼矩阵满足C＝0.01M+0.01K，U是全局位移， 和 是加速度和速度。
[0035]
[0036]
有限元分析是一种基于离散逼近法逼近连续复变函数的数值方法，考虑到心肌的非刚体变形，可将全局节点位移U与全局节点力F的关系由刚度方程K简单的联立起来，将上面的连续介质力学方程简化，建立了位移U和
载荷F的之间的直接目标方程KU＝F，并施加GTV约束，建立目标方程：
[0037]
[0038] 其中： 作为正则化项来处理噪声平滑，μ为GTV的正则化参数，||||2表示二范数。
[0039] (3)计算目标函数的一个重要步骤是构造GTV正则化步骤中的图模型G；图是一种自然的数据结构，其中对象由节点表示，它们之间的关系可以由图中的边来确定，而每条边通常表示它所连接的两个
顶点之间的相似性。为了从噪声样本中重建信号，需要建立平滑模型，平滑度的另一个概念是图的变化，它量化了原始图信号和滤波信号之间的差异。
[0040] GTV参见文献《Mahmood,Faisal,et al."Adaptive Graph-based Total Variation for Tomographic Reconstructions."IEEE Signal Processing Letters 1-1(2018).》，GTV算法的图概念将量化原始位移的图信号和滤波后的位移图信号之间的差作为图形变化项进行平滑，对所述的最优化求解表达式进行求解。
[0041] (4)通过KNN算法对数据相似度进行分类来计算图全变分，并以直观、灵活的方式获得数据中的相似关系，使信号在无约束的空间内平稳地减少噪声。
[0042] 图权重矩阵Wij采用和使用高斯
内核权重计算方案， 将计算i点的所有邻居ζi间的连接，α为连接样本间的平均距离，本发明图的全变差正则化数学模型为：
[0043]
[0044]
[0045] (5)采用前向-后向原对偶法代替传统的优化方法，即结合
梯度下降法和相邻算子求解模型中的不可微项，原对偶问题可用于求解特殊凸优化问题，该方法可同时用于求解微分问题和非微分问题，前者通过计算梯度直接计算出显式步长，后者通过逼近算子的隐式步长逼近。
[0046]
[0047] 其中：
[0048] 对于一个关于U的基本优化问题，上述式子可以视之为原始问题，则这一原始问题的对偶问题如下：
[0049]
[0050] 我们定义对偶式子是一个复凸函数，GTV项定义为g(D)＝μ||D||1， g*(V)表示其共轭；这种双重形式可以分裂问题，提高计算效率，尤其是当Γ值小于Λ值，近端操作符°表示阿达玛积。
[0051]
[0052] 主要的近似分裂算法通常包括前向-后向算法，它是Douglas-Rachford算法的复形，将该算法与流行的原对偶理论相结合，算法主要参数为两个步骤参数τ和松σ，弛因子λn和正则化参数μ，通过调整参数，以达到最佳的
迭代算法的收敛性。初始值U0＝ULSM通过基于有限元的最小二乘法计算， 则迭代优化步骤如下：
[0053] Scheme：
[0054]
[0055] tn＝Vn+σL(2pn-Un)
[0056]
[0057] (Un+1，Vn+1)＝(Un，Vn)+λn((pn-qn)-(Un，Vn))
[0058] Algorithm：giveUn，Vn，λncompute pn，qn
[0059] give pn，qncompute Un+1，Vn+1
[0060] update
[0061] 其中：Un，Vn互为原对偶问题，prox为近端操作算子，*为其共轭。
[0062] 我们通过以下方法对比数据来验证本实施方式的优越性和可靠性，表1比较了图全变分算法和最小二乘法LSM在5dB、10dB和20dB不同
信噪比下的误差，将xt真值数据和xr算法数据代入下式，计算相对误差和相关系数。
[0063]
[0064]
[0065] 表1
[0066]
[0067] 我们可以对比数据看出，与传统的有限元拟合最小二乘法相比，本发明GTV模型具有更高的计算速度和收敛速度。
[0068] 上述对
实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种
修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。
