一种基于原对偶算法的非凸加权变分金属伪影去除方法

技术领域

本发明涉及伪影去除的技术领域，尤其涉及一种基于原对偶算法的非凸加权变分金属伪影去除方法。

背景技术

计算机断层成像(Computed Tomography,CT)通过具有穿透能力和多色能量的射线(x射线、γ射线)扫描物体，物质吸收射线从而能量逐渐衰减，随之由探测器接收，最后使用重建技术得到物质内部信息，是最常见的医学诊断和治疗手段之一。然而，由于金属植入物的普遍存在，如补牙、线圈、髋关节假体等，射线扫描至金属时，能量急剧衰减，甚至探测光子数为0，利用传统的重建算法(例如：滤波反投影(FBP)和代数重建技术(SART))得到的图像出现明暗条纹伪影，严重降低图像质量，从而可能导致误诊。最简明的一个例子是一张包含铁金属的头部CT，如图4的(a)所示，真实的头部CT图像如图4的(b)所示，经过传统重建算法得到的图像严重退化，尤其在高衰减物质金属及骨头附近含有明显的金属伪影。因此，开发有效的金属伪影减少(MAR)方法是必要且有意义的。

一种有效的方法是结合成像原理校正系统模型，该方案可以减少伪影，同时保留边界和细节。然而，当处理高水平的噪声时，上述模型与真实模型之间存在差异，这不可避免地导致令人不满意的结果。另一种众所周知的方法是将受金属影响的投影(金属轨迹，正如图3的(m)为图3的(l)所示的实际测量投影数据对应的金属轨迹)识别为缺失数据，并使用插值来恢复。一些方法是使用不同的插值算法直接修复金属轨迹，算法简单快速但容易导致边界不准确同时引入新的伪影。为了提高性能，归一化金属伪影校正(NMAR)方法利用先验图像对测量数据进行归一化，然后再对金属区域的正弦图进行插值，减少了伪影引入且保留部分组织边缘信息，但是此方法高度依赖先验图像的质量。近年来，随着深度学习在医学图像处理中的显著成功，最近已经实现了深度神经网络(DNN)解决减少金属伪影的问题。现有的研究包括利用残差学习或对抗性学习技术完成直接的图像到图像学习或者正弦图域，再或者基于图像和正弦图域的双域学习。然而，DNN通常需要大型、有代表性的训练数据集和广泛的计算资源，这限制了此类方法的广泛应用。

作为MAR的一种关键数学技术，变分方法在适度的局部资源需求下提供了可解释性、稳定性和计算效率。由于校正过程涉及解决数学上不适定的问题，因此基于正则化的方法起着至关重要的作用。大多数的变分模型一部分基于CT图像的分片光滑性，引入全变分型(TV)正则直接修复金属轨迹(所有通过金属的投影)，有效抑制金属伪影及噪声，但重建结果在高密度物质周围容易产生强扩散，甚至引入新的伪影。另一部分方法利用预分割的投影归一化Radon域，极大地改善了模型的表现，但这需要相对准确的预分割。然而，分割被金属伪影污染的CT图像也是具有挑战性的。此外，具有TV或紧帧的凸正则化模型易于降低边缘对比度。

申请号为202310251399.X的发明专利申请公开了一种CT图像的金属伪影去除方法，包括：构建基于小波变换的初始自适应迭代学习模型；根据CT图像的纯净图像区域、二元非金属区域和金属伪影区域确定CT图像表达式；根据所述CT图像表达式、正则化项和图像小波变换构建第一优化目标函数；结合近端梯度下降算法和泰勒公式求解所述第一优化目标函数，得到第一计算结果；根据所述第一计算结果，结合L2输出损失和地面真实图像优化所述初始自适应迭代学习模型，获得目标自适应迭代学习模型；通过所述目标自适应迭代学习模型对所述CT图像进行图像分解和拼接以去除CT图像的金属伪影。上述发明利用金属伪影在不同域的空间分布特征，通过网络得到的图像域的金属伪影并进行伪影去除。但是，上述方法要求大型的、具有代表性的训练数据集和广泛的计算资源。

发明内容

针对现有金属伪影去除的变分方法的重建结果一方面在高密度物质周围容易产生强扩散，另一方面高度依赖先验图像的技术问题，本发明提出一种基于原对偶算法的非凸加权变分金属伪影去除方法，通过简单的阈值操作得到金属位置后，在投影域粗略表征金属伪影区域并去除伪影，对于不同水平的伪影都可以得到较好的MAR结果，且鲁棒性较高。

为了达到上述目的，本发明的技术方案是这样实现的：一种基于原对偶算法的非凸加权变分金属伪影去除方法，其步骤如下：

步骤一：根据实测投影数据求解成像模型得到包含噪声和金属伪影的初始重建图像；

步骤二：对初始重建图像进行阈值化得到金属位置，通过结合实测投影数据的幂次方和高度衰减投影区域为0计算加权矩阵，得到数据拟合项；

步骤三：使用加权各向异性和各向同性TV建立正则化项；

步骤四：联合数据拟合项和正则化项构建具有盒约束的非凸加权模型，并基于TV预对偶形式的分裂算法求解非凸加权模型，获得最终的金属伪影校正图像。

优选地，采用单色能量假设，所述成像模型为：

Pu≈Y；

其中，u表示待重建的处于特定但未知能量水平的目标图像，P代表投影矩阵，Y表示多能级测量下的实际投影数据。

优选地，所述求解成像模型的方法为传统的梯度下降方法或者采用共轭梯度法。

优选地，所述对初始重建图像进行阈值化得到金属位置的方法为：

金属位置为：M:＝{(i,j)∈X|u

优选地，所述加权矩阵为：

其中ε是一个正常数，

其中，Ω

优选地，所述高度衰减投影的区域

优选地，所述数据项为：||W⊙(Pu-Y)||

所述正则化项

具有盒约束的非凸加权模型为：

其中，⊙表示元素相乘，常数c是重建图像的上界，λ和α是两个正则化权值参数。

优选地，所述基于TV预对偶形式的分裂算法求解非凸加权模型的方法为：

基于||·||

其中，闭集

其中，p,q∈X×X分别为关于

其中，(p

添加了一个附加的二次项得到惩罚模型：

其中，η为罚参数；

引入辅助变量v＝Pu，以拉格朗日乘子

因此，转化为求解鞍点问题：

min

优选地，利用完全分裂的原始对偶混合梯度算法求解鞍点问题：min

Λ

其中，

关于变量u、v、q和p的子问题是严格凸的，且具有闭式解：

u

其中，Proj(.；.)代表投影算子，div代表散度算子。

优选地，所述完全分裂的原始对偶混合梯度算法的求解步骤为：

选取初始的正则化权值参数λ、α、罚参数η、步长参数ρ、σ

迭代：更新拉格朗日乘子Λ

计算目标图像u

更新外插变量

计算辅助变量v

并行计算对偶变量q

满足终止条件：

与现有技术相比，本发明的有益效果：

1)本发明基于金属伪影的特点，数据拟合项融合了自适应加权范数，可以实现更有效地金属伪影校正；也就是说，权重函数是利用测量投影的混合方案构建的，而不是简单地设置为二值矩阵。

2)本发明采用基于非凸全变分L

3)本发明通过在Radon域中引入辅助变量，提出了一种更快的完全分裂原对偶混合梯度算法，该算法在非凸算法理论框架下保证了收敛性。且算法可以有效地实现，因为每个子问题都有一个闭式解。

4)本发明进行了大量实验，包含不同数量组织的CT图像，对于不同水平的伪影都可以得到较好的MAR结果，展现了本发明的有效性及鲁棒性。

附图说明

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

图1为本发明的流程图。

图2为非凸正则与凸正则的结果图，其中，(a)为非凸正则的校正图像，(b)为使用L1凸正则的校正图像，(c)为L2凸正则的校正图像。

图3为本发明与现有方法的重建结果，其中，(a)NCAT的初始重建图像，(b)为真实的CT图像，(c)为光束硬化校正算子BCMAR方法，(d)为归一化金属伪影校正NMAR，(e)为全变分正则化金属伪影校正TV-MAR，(f)为重加权双域CT金属伪影校正Reweighted JSR，(g)为本发明的校正结果，(h)为(e)的矩形局部区域图，(i)为(f)的矩形局部区域图，(j)为(g)的矩形局部区域图，(k)为(b)的矩形局部区域图，(l)为实际测量投影图，(m)为(l)对应的金属轨迹图。

图4为本发明与现有方法对“Head”图像的重建结果，其中，(a)为待校正的“Head”图像，(b)为真实的CT图像，(c)为光束硬化校正算子BCMAR方法的校正结果，(d)为归一化金属伪影校正NMAR的校正结果，(e)为全变分正则化金属伪影校正TV-MAR的校正结果，(f)为重加权双域CT金属伪影校正Reweighted JSR的校正结果，(g)为本发明的校正结果。

图5为本发明与现有方法对“Skull”图像的重建结果，其中，(a)为待校正的“Skull”图像，(b)为真实的CT图像，(c)为利用光束硬化校正算子BCMAR方法的校正结果，(d)为归一化金属伪影校正NMAR的校正结果，(e)为全变分正则化金属伪影校正TV-MAR的校正结果，(f)为重加权双域CT金属伪影校正Reweighted JSR的校正结果，(g)为本发明的校正结果。

具体实施方式

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

如图1所示，本发明提出了一种基于原对偶算法的非凸加权变分金属伪影去除方法，根据金属伪影的特征，将数据拟合项设计为正向投影的估计值与实测投影数据之间的加权范数以更有效地确定CT图像，之后，设计图像梯度的非凸正则化以消除噪声，最后结合数据项和正则项构建优化目标函数并获得最终的MAR结果，其具体步骤如下：

步骤一：通过共轭梯度法求解成像模型得到包含噪声和金属伪影的初始重建图像。

根据朗伯-比尔定律，由多色x射线测量的投影数据表示为：

Y＝-ln(∑

其中，

对于低衰减材料，衰减系数相对能量是不敏感的，采用单色能量假设可以将成像模型改写为：

Pu≈Y； (1)

其中，u表示待重建的处于特定但未知能量水平的目标图像，P代表投影矩阵，由成像设备决定，Y表示多能级测量下的实际投影数据。求解式(1)的线性方程组可以利用传统的梯度下降方法或者采用共轭梯度法，梯度下降法尽管计算量小，但是收敛速度慢，计算效率低，而共轭梯度法在保证低计算量的前提下提高了计算效率。

但当存在具有高能量依赖性的高衰减组件，例如金属植入物，此时单能假设将不再准确。直接利用共轭梯度法求解式(1)会导致严重的金属伪影，从而大大降低图像质量，正如图3的(a)所示，金属之间存在明显的伪影。

步骤二：对初始重建图像进行阈值化得到金属位置，通过结合低于阈值水平t衰减的实测投影数据的幂次方和高于阈值水平t衰减投影区域为0计算加权矩阵，并得到迭代优化模型的数据拟合项。

在图像域，具有更高的CT值的金属区域很容易被分割得到。具体地，获得初始重建图像后，通过简单的阈值操作得到金属位置M，即：

M:＝{(i,j)∈X|u

其中ε是一个正常数，以避免分母中的零，

其中，Ω

且高度衰减投影的区域

此时，数据项表达式为：||W⊙(Pu-Y)||

如果只考虑投影数据的非金属信息(即使用二值矩阵

步骤三：使用加权各向异性和各向同性TV正则化以消除噪声。

数据拟合项确保和原始图像的拟合程度，正则化则描述CT图像的先验信息，保证校正结果的光滑性和对比度。非凸正则化是一个更好的选择。正如图2所示，图2的(a)为选用非凸正则的校正图像，图2的(b)和图2的(c)为凸正则的结果，尽管都可以有效地保证图像的平滑性，但非凸正则相比之下有更高的对比度。本发明使用加权各向异性和各向同性TV(AITV)：

可以描述图像u的先验信息，即利用/>

步骤四：联合数据拟合项和正则化项构建具有盒约束的非凸加权模型，并基于TV预对偶形式的分裂算法求解非凸加权模型，获得最终的金属伪影校正图像。

基于上述理论，本发明提出下面的具有盒约束的非凸加权模型：

其中，⊙表示元素相乘，常数c是重建图像的上界，λ和α是两个正则化权值参数，经实验测试正则化权值参数α＝0.75具有更好的重建效果，λ是衡量正则化的强度的调优参数。

由于非凸正则化，非凸加权模型是一个多变量非凸优化问题。这里本发明考虑了一种基于TV预对偶形式的分裂算法。

基于||·||

其中，闭集

其中，p,q∈X×X分别为关于

其中，(p

为了表征这种非凸优化模型的分裂算法的收敛性，在式(4)中添加了一个附加的二次项，从而得到如下的惩罚模型：

其中，η为罚参数。关于变量p的上述优化问题是强凹的。注意，只有当η＝0时，它才能返回到式(4)的原始模型。这种惩罚是为了增强图像的归一化梯度的平滑度。

进一步引入辅助变量v＝Pu，以拉格朗日乘子

因此，需要如下解决鞍点问题：

第k+1次迭代中的完全分裂的原始对偶混合梯度算法(FS-PDHG)的方案可以描述为：

Λ

为用来加速算法的外插变量，正参数ρ、σ

u

其中，Proj(.；.)代表投影算子，div代表散度算子。

FS-PDHG算法总结如下：

1.初值选取：选取合适的参数λ、α、η、ρ、σ

2.迭代：

利用式(5)更新拉格朗日乘子Λ

利用式(7)计算目标图像u

利用式(6)更新外插变量

利用式(8)计算辅助变量v

利用式(9)和式(10)并行计算对偶变量q

满足终止条件：

求解非凸加权模型的经典思想为将非凸加权正则表示为两个凸正则的差，再将L2正则线性化表示求解，完全分裂的原始对偶混合梯度算法(FS-PDHG)求解非凸加权模型避免引入线性化操作，而是将其转化为预对偶形式直接求解。

本发明的实验平台：因特尔酷睿处理器11300H，CPU@3.10赫兹，16G内存个人笔记本，使用Matlab2021进行程序设计。

具体实验一

将图3的(a)(NCAT)利用本发明和现有的部分MAR方法进行重建并对比，结果见图3。图3的(b)为真实的CT图像，图3的(c)为光束硬化校正算子BCMAR方法的校正结果，图3的(d)为归一化金属伪影校正NMAR的校正结果，图3的(e)为全变分正则化金属伪影校正TV-MAR的校正结果，图3的(f)为重加权双域CT金属伪影校正Reweighted JSR的校正结果，图3的(g)为本发明的校正结果，图3(h)为图3(e)的矩形局部区域图，图3的(i)为图3的(f)的矩形局部区域图，图3的(j)为图3的(g)的矩形局部区域图，图3的(k)为图3的(b)的矩形局部区域图。

图3表明，BCMAR方法在较高的噪声水平下的校正结果仍然存在明显的伪影和噪声；NMAR校正后总体图像质量得到部分提升，但金属周围要模糊于TV-MAR；TV-MAR在抑制伪影和噪声方面表现优于之前的算法，但是由红色箭头处可观察到骨头边界和金属出现融合，且图像的对比度较低。而Reweighted JSR和本发明在保留细节的同时，有效抑制了伪影和噪声。

为了更好的说明本发明的有效性，放大使用变分方法校正图3的(b)的局部区域，即图3的(h)-图3(k)。通过细节图可以清晰看到TV-MAR对其校正后，靠近金属周围的暗伪影得到了有效抑制，但靠近金属的骨骼模糊严重，金属边缘的次级伪影也很明显。

具体实验二

而Reweighted JSR高度依赖分割结果，以下的实验进一步说明本发明具有更好的有效性和鲁棒性。

将本发明和已有的部分MAR方法分别对“Head”和“Skull”的仿真图像进行校正并对比，结果如图4和图5。

图4的(a)为待校正的“Head”图像，图4的(b)为真实的CT图像，图4的(c)为利用BCMAR方法对图4的(a)的校正结果，图4的(d)为NMAR的校正结果，图4的(e)为TV-MAR的校正结果，图4的(f)为Reweighted JSR的校正结果，图的4(g)为本发明的校正结果。图5的(a)为待校正的“Skull”图像，图5的(b)为真实的CT图像，图5的(c)为利用BCMAR方法对图5的(a)的校正结果，图5的(d)为NMAR的校正结果，图5的(e)为TV-MAR的校正结果，图5的(f)为Reweighted JSR的校正结果，图5的(g)为本发明的校正结果。

由图4可以看出，由于“Head”图像存在大量骨头导致更多的伪影，极大地降低骨头附近图像质。BCMAR、NMAR和TV-MAR与NCAT结果类似，而且NMAR由于粗略的插值在金属与骨头之间引入新的伪影，TV-MAR则是直接在投影域恢复金属区域将导致边界扩散严重且引入新的伪影。相较于NCAT，更强的伪影降低了预分割的精度，进一步使得Reweighted JSR引入了不被期望的结构。但是本发明的结果类似于NCAT在去除伪影的同时保留了细节，有最好的整体效果。

进一步，通过观察含有较多细节的图5，TV-MAR表现出严重的伪影，ReweightedJSR因为预分割的局限导致不精确的伪影修正结果，而本发明持续在减少金属伪影和恢复组织细节方面展现了最佳的性能。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。