CN117152291B - 一种基于原对偶算法的非凸加权变分金属伪影去除方法 - Google Patents
一种基于原对偶算法的非凸加权变分金属伪影去除方法 Download PDFInfo
- Publication number
- CN117152291B CN117152291B CN202311169123.3A CN202311169123A CN117152291B CN 117152291 B CN117152291 B CN 117152291B CN 202311169123 A CN202311169123 A CN 202311169123A CN 117152291 B CN117152291 B CN 117152291B
- Authority
- CN
- China
- Prior art keywords
- metal
- image
- convex
- dual
- algorithm
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 239000002184 metal Substances 0.000 title claims abstract description 92
- 229910052751 metal Inorganic materials 0.000 title claims abstract description 92
- 238000000034 method Methods 0.000 title claims abstract description 48
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 41
- 230000009977 dual effect Effects 0.000 title claims abstract description 36
- 238000012937 correction Methods 0.000 claims abstract description 55
- 239000011159 matrix material Substances 0.000 claims abstract description 20
- 230000002238 attenuated effect Effects 0.000 claims abstract description 11
- 238000003384 imaging method Methods 0.000 claims abstract description 11
- 230000006870 function Effects 0.000 claims description 11
- 238000002939 conjugate gradient method Methods 0.000 claims description 6
- 238000013213 extrapolation Methods 0.000 claims description 6
- 238000005259 measurement Methods 0.000 claims description 6
- 150000002739 metals Chemical class 0.000 claims description 5
- 238000011478 gradient descent method Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 239000007769 metal material Substances 0.000 claims description 3
- 230000002596 correlated effect Effects 0.000 claims description 2
- 210000000988 bone and bone Anatomy 0.000 description 7
- 238000004364 calculation method Methods 0.000 description 6
- 238000005457 optimization Methods 0.000 description 6
- 238000013459 approach Methods 0.000 description 4
- 230000000875 corresponding effect Effects 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 4
- 229910052704 radon Inorganic materials 0.000 description 4
- SYUHGPGVQRZVTB-UHFFFAOYSA-N radon atom Chemical compound [Rn] SYUHGPGVQRZVTB-UHFFFAOYSA-N 0.000 description 4
- 230000011218 segmentation Effects 0.000 description 4
- 210000003625 skull Anatomy 0.000 description 4
- 239000000126 substance Substances 0.000 description 4
- 230000001419 dependent effect Effects 0.000 description 3
- 238000009792 diffusion process Methods 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 210000001519 tissue Anatomy 0.000 description 3
- 238000002591 computed tomography Methods 0.000 description 2
- 239000007943 implant Substances 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 230000005855 radiation Effects 0.000 description 2
- 230000008439 repair process Effects 0.000 description 2
- 238000012549 training Methods 0.000 description 2
- CWYNVVGOOAEACU-UHFFFAOYSA-N Fe2+ Chemical compound [Fe+2] CWYNVVGOOAEACU-UHFFFAOYSA-N 0.000 description 1
- 230000003042 antagnostic effect Effects 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 230000000593 degrading effect Effects 0.000 description 1
- 239000004053 dental implant Substances 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 239000003226 mitogen Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000001225 therapeutic effect Effects 0.000 description 1
- 238000004800 variational method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明提出了一种基于原对偶算法的非凸加权变分金属伪影去除方法，步骤为：根据实测投影数据求解成像模型得到包含噪声和金属伪影的初始重建图像；对初始重建图像进行阈值化得到金属位置，通过结合实测投影的幂次方和高度衰减投影区域为0计算加权矩阵，得到数据拟合项；使用加权各向异性和各向同性TV建立正则化项；联合数据拟合项和正则化项构建具有盒约束的非凸加权模型，基于TV预对偶形式的分裂算法求解非凸加权模型，获得金属伪影校正图像。本发明基于金属伪影的特点，数据拟合项融合了自适应加权范数，可以实现更有效地金属伪影校正；采用正则化来惩罚分片光滑函数，在保持图像稀疏性的前提下提高了图像的对比度，保证了收敛性。
Description
技术领域
本发明涉及伪影去除的技术领域，尤其涉及一种基于原对偶算法的非凸加权变分金属伪影去除方法。
背景技术
计算机断层成像(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|uo,k≥C}，ui,j代表初始重建图像在位置(i,j)处的像素值，C为初始重建图像的金属材料的像素临界值；为图像空间，n,m分别为每行、列所包含的像素个数。
优选地，所述加权矩阵为：
其中ε是一个正常数，是一个全1矩阵；/>是与子集Ωt相关的二值矩阵，且：
其中，Ωt＝Om∪Ot是金属轨迹Ω的子集，Om为每两个分离金属在投影区域的交集，Ot为高度衰减投影的区域。
优选地，所述高度衰减投影的区域 具有阈值水平t和最大像素值其中，Yi,j与/>为实际投影数据Y分别在(i,j)、/>处的像素值。
优选地，所述数据拟合项为：||W⊙(Pu-Y)||2；
所述正则化项描述图像的先验信息；
具有盒约束的非凸加权模型为：
其中，⊙表示元素相乘，常数c是重建图像的上界，λ和α是两个正则化权值参数。
优选地，所述基于TV预对偶形式的分裂算法求解非凸加权模型的方法为：
基于||·||2,1和||·||1的预对偶形式，鞍点问题写成：
其中，闭集 是关于闭集U的指示符函数：
其中，p,q∈X×X分别为关于和/>的对偶形式引入的对偶变量，/>和分别是关于闭集S和Q的指示符函数，且闭集
其中，(px)i,j、(py)i,j、(qx)i,j、(qy)i,j分别表示对偶变量p,q在x或y方向上位于(i,j)处的像素值；
添加了一个附加的二次项得到惩罚模型：
其中，η为罚参数；
引入辅助变量v＝Pu，以拉格朗日乘子为加权系数将约束增加到目标函数，得：
因此，转化为求解鞍点问题：
minu,q,vmaxp,ΛL(u,v,q,p,Λ)。
优选地，利用完全分裂的原始对偶混合梯度算法求解鞍点问题：minu,q,vmaxp,ΛL(u,v,q,p,Λ)的方法为：第k+1次迭代中：
Λk+1＝Λk+ρ(vk-Puk)；
其中，为用来加速算法的外插变量，正参数ρ、σ1、σ2、τ和β为步长参数；
关于变量u、v、q和p的子问题是严格凸的，且具有闭式解：
uk+1＝Proj(σ1div(pk+αqk)+σ1PTΛk+1+uk；U)；
其中，Proj(；)代表投影算子，div代表散度算子。
优选地，所述完全分裂的原始对偶混合梯度算法的求解步骤为：
选取初始的正则化权值参数λ、α、罚参数η、步长参数ρ、σ1、σ2和β，并令初始值u0、v0、q0、p0＝0；
迭代：更新拉格朗日乘子Λk+1；
计算目标图像uk+1；
更新外插变量
计算辅助变量vk+1；
并行计算对偶变量qk+1、pk+1；
满足终止条件：tol为给定的正常数，迭代终止。
与现有技术相比，本发明的有益效果：
1)本发明基于金属伪影的特点，数据拟合项融合了自适应加权范数，可以实现更有效地金属伪影校正；也就是说，权重函数是利用测量投影的混合方案构建的，而不是简单地设置为二值矩阵。
2)本发明采用基于非凸全变分L1-αL2范数的正则化来惩罚分片光滑函数，在保持图像稀疏性的前提下提高了图像的对比度。
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(∑Eη(E)exp(-PuE))；
其中，P:是Radon变换，表示在m1个不同投影角并沿着m2个不同探测器的离散线积分，η(E)＝S0(E)/∑ES0(E)表示归一化后的能谱，S0(E)和uE∈X分别是能级E∈[Emin,Emax]下的入射光子的数量和线性衰减系数，Emin,Emax分别表示最小和最大的能级，图像空间/>n,m分别为每行、列所包含的像素个数。显然，相应的重建过程是非线性的，求解这样一个大规模系统非常具有挑战性。
对于低衰减材料，衰减系数相对能量是不敏感的，采用单色能量假设可以将成像模型改写为：
Pu≈Y；(1)
其中，u表示待重建的处于特定但未知能量水平的目标图像，P代表投影矩阵，由成像设备决定，Y表示多能级测量下的实际投影数据。求解式(1)的线性方程组可以利用传统的梯度下降方法或者采用共轭梯度法，梯度下降法尽管计算量小，但是收敛速度慢，计算效率低，而共轭梯度法在保证低计算量的前提下提高了计算效率。
但当存在具有高能量依赖性的高衰减组件，例如金属植入物，此时单能假设将不再准确。直接利用共轭梯度法求解式(1)会导致严重的金属伪影，从而大大降低图像质量，正如图3的(a)所示，金属之间存在明显的伪影。
步骤二：对初始重建图像进行阈值化得到金属位置，通过结合低于阈值水平t衰减的实测投影数据的幂次方和高于阈值水平t衰减投影区域为0计算加权矩阵，并得到迭代优化模型的数据拟合项。
在图像域，具有更高的CT值的金属区域很容易被分割得到。具体地，获得初始重建图像后，通过简单的阈值操作得到金属位置M，即：
M＝{(i,j)∈X|ui,j≥C}，ui,j代表初始重建图像在位置(i,j)处的像素值，C为初始重建图像的金属材料的像素临界值。本发明计算了加权矩阵
其中ε是一个正常数，以避免分母中的零，是一个全1矩阵。这里是与子集Ωt相关的二值矩阵，定义如下：
其中，Ωt＝Om∪Ot是金属轨迹Ω(金属位置的投影所对应的区域)的子集，Om为每两个分离金属在投影区域的交集，Ot为高度衰减投影的区域。
且高度衰减投影的区域具有阈值水平t和其中，Yi,j与/>为实际投影数据Y分别在(i,j)、/>处的像素值。
此时，数据拟合项表达式为：||W⊙(Pu-Y)||2。
如果只考虑投影数据的非金属信息(即使用二值矩阵作为加权矩阵而不是加权矩阵W)，图像域中的边界可能会扩散，并且还会出现新的伪影。这可能是由于在金属迹线中丢弃了过多的结构而导致的。因此，这里建议只丢弃污染最严重的投影，它们要么是两种金属的交叉投影Om，要么是与其他相对高衰减材料(例如骨骼)相关的区域Ot。测量的投影数据的平方/>倒数被组合为式(2)，以便自适应地平衡误差分布，即对于更高的衰减投影，允许在Radon域中有更大的误差。
步骤三：使用加权各向异性和各向同性TV正则化以消除噪声。
数据拟合项确保和原始图像的拟合程度，正则化则描述CT图像的先验信息，保证校正结果的光滑性和对比度。非凸正则化是一个更好的选择。正如图2所示，图2的(a)为选用非凸正则的校正图像，图2的(b)和图2的(c)为凸正则的结果，尽管都可以有效地保证图像的平滑性，但非凸正则相比之下有更高的对比度。本发明使用加权各向异性和各向同性TV(AITV)：
可以描述图像u的先验信息，即利用/>表达图像在梯度变化下的的稀疏性，而AITV(u)同时具有强调保边的能力。
步骤四：联合数据拟合项和正则化项构建具有盒约束的非凸加权模型，并基于TV预对偶形式的分裂算法求解非凸加权模型，获得最终的金属伪影校正图像。
基于上述理论，本发明提出下面的具有盒约束的非凸加权模型：
其中，⊙表示元素相乘，常数c是重建图像的上界，λ和α是两个正则化权值参数，经实验测试正则化权值参数α＝0.75具有更好的重建效果，λ是衡量正则化的强度的调优参数。
由于非凸正则化，非凸加权模型是一个多变量非凸优化问题。这里本发明考虑了一种基于TV预对偶形式的分裂算法。
基于||·||2,1和||·||1的预对偶形式，相应的鞍点问题可以写成：
其中，闭集 是关于闭集U的指示符函数：
其中，p,q∈X×X分别为关于和/>的对偶形式所引入的对偶变量，/>和/>分别是关于闭集S和Q的指示符函数，且闭集
其中，(px)i,j、(py)i,j、(qx)i,j、(qy)i,j分别表示对偶变量p,q在x或y方向上位于(i,j)处的像素值。
为了表征这种非凸优化模型的分裂算法的收敛性，在式(4)中添加了一个附加的二次项，从而得到如下的惩罚模型：
其中，η为罚参数。关于变量p的上述优化问题是强凹的。注意，只有当η＝0时，它才能返回到式(4)的原始模型。这种惩罚是为了增强图像的归一化梯度的平滑度。
进一步引入辅助变量v＝Pu，以拉格朗日乘子为加权系数将约束增加到目标函数，如下所示：
因此，需要如下解决鞍点问题：
第k+1次迭代中的完全分裂的原始对偶混合梯度算法(FS-PDHG)的方案可以描述为：
Λk+1＝Λk+ρ(vk-Puk)； (5)
为用来加速算法的外插变量，正参数ρ、σ1、σ2、τ和β为步长参数。这些关于变量u、v、q和p的子问题是严格凸的，并且具有闭式解：
uk+1＝Proj(σ1div(pk+αqk)+σ1PTΛk+1+uk；U)； (7)
其中，Proj(；)代表投影算子，div代表散度算子。
FS-PDHG算法总结如下：
1.初值选取：选取合适的参数λ、α、η、ρ、σ1、σ2和β，并令初始值u0、v0、q0、p0＝0。
2.迭代：
利用式(5)更新拉格朗日乘子Λk+1；
利用式(7)计算目标图像uk+1；
利用式(6)更新外插变量
利用式(8)计算辅助变量vk+1；
利用式(9)和式(10)并行计算对偶变量qk+1、pk+1；
满足终止条件：tol为经验给定的正常数，迭代终止。
求解非凸加权模型的经典思想为将非凸加权正则表示为两个凸正则的差，再将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因为预分割的局限导致不精确的伪影修正结果，而本发明持续在减少金属伪影和恢复组织细节方面展现了最佳的性能。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (7)
1.一种基于原对偶算法的非凸加权变分金属伪影去除方法，其特征在于，其步骤如下：
步骤一：根据实测投影数据求解成像模型得到包含噪声和金属伪影的初始重建图像；
步骤二：对初始重建图像进行阈值化得到金属位置，通过结合实测投影数据的幂次方和高度衰减投影区域为0计算加权矩阵，得到数据拟合项；
所述加权矩阵为：
其中ε是一个正常数，是一个全1矩阵；/>是与子集Ωt相关的二值矩阵，且：
其中，Ωt＝Om∪Ot是金属轨迹Ω的子集，Om为每两个分离金属在投影区域的交集，Ot为高度衰减投影的区域；
步骤三：使用加权各向异性和各向同性TV建立正则化项；
所述数据拟合项为：||W⊙(Pu-Y)||2
所述正则化项描述图像的先验信息；
步骤四：联合数据拟合项和正则化项构建具有盒约束的非凸加权模型，并基于TV预对偶形式的分裂算法求解非凸加权模型，获得最终的金属伪影校正图像；
具有盒约束的非凸加权模型为：
其中，⊙表示元素相乘，常数c是重建图像的上界，λ和α是两个正则化权值参数；u表示待重建的处于特定但未知能量水平的目标图像，P代表投影矩阵，Y表示多能级测量下的实际投影数据；
所述基于TV预对偶形式的分裂算法求解非凸加权模型的方法为：
基于||·||2,1和||·||1的预对偶形式，鞍点问题写成：
其中，闭集是关于闭集U的指示符函数：
其中，p,q∈X×X分别为关于和/>的对偶形式引入的对偶变量，/>和分别是关于闭集S和Q的指示符函数，且闭集
其中，(px)i,j、(py)i,j、(qx)i,j、(qy)i,j分别表示对偶变量p,q在x或y方向上位于(i,j)处的像素值；
添加了一个附加的二次项得到惩罚模型：
其中，η为罚参数；
引入辅助变量v＝Pu，以拉格朗日乘子为加权系数将约束增加到目标函数，得：
因此，转化为求解鞍点问题：
minu,q,vmaxp,ΛL(u,v,q,p,Λ)。
2.根据权利要求1所述的基于原对偶算法的非凸加权变分金属伪影去除方法，其特征在于，采用单色能量假设，所述成像模型为：
Pu≈Y。
3.根据权利要求2所述的基于原对偶算法的非凸加权变分金属伪影去除方法，其特征在于，所述求解成像模型的方法为传统的梯度下降方法或者采用共轭梯度法。
4.根据权利要求1-3中任意一项所述的基于原对偶算法的非凸加权变分金属伪影去除方法，其特征在于，所述对初始重建图像进行阈值化得到金属位置的方法为：
金属位置为：M＝{(i,j)∈X|ui,j≥C}，ui,j代表初始重建图像在位置(i,j)处的像素值，C为初始重建图像的金属材料的像素临界值；为图像空间，n,m分别为每行、列所包含的像素个数。
5.根据权利要求1所述的基于原对偶算法的非凸加权变分金属伪影去除方法，其特征在于，所述高度衰减投影的区域具有阈值水平t和最大像素值/>其中，Yi,j与/>为实际投影数据Y分别在(i,j)、/>处的像素值。
6.根据权利要求1所述的基于原对偶算法的非凸加权变分金属伪影去除方法，其特征在于，利用完全分裂的原始对偶混合梯度算法求解鞍点问题：minu,q,vmaxp,ΛL(u,v,q,p,Λ)的方法为：第k+1次迭代中：
Λk+1＝Λk+ρ(vk-Puk)；
其中，为用来加速算法的外插变量，正参数ρ、σ1、σ2、τ和β为步长参数；
关于变量u、v、q和p的子问题是严格凸的，且具有闭式解：
uk+1＝Proj(σ1div(pk+αqk)+σ1PTΛk+1+uk；U)；
其中，Proj(；)代表投影算子，div代表散度算子。
7.根据权利要求6所述的基于原对偶算法的非凸加权变分金属伪影去除方法，其特征在于，所述完全分裂的原始对偶混合梯度算法的求解步骤为：
选取初始的正则化权值参数λ、α、罚参数η、步长参数ρ、σ1、σ2和β，并令初始值u0、v0、q0、p0＝0；
迭代：更新拉格朗日乘子Λk+1；
计算目标图像uk+1；
更新外插变量
计算辅助变量vk+1；
并行计算对偶变量qk+1、pk+1；
满足终止条件：tol为给定的正常数，迭代终止。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311169123.3A CN117152291B (zh) | 2023-09-12 | 2023-09-12 | 一种基于原对偶算法的非凸加权变分金属伪影去除方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311169123.3A CN117152291B (zh) | 2023-09-12 | 2023-09-12 | 一种基于原对偶算法的非凸加权变分金属伪影去除方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117152291A CN117152291A (zh) | 2023-12-01 |
| CN117152291B true CN117152291B (zh) | 2024-03-22 |
Family
ID=88911849
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311169123.3A Active CN117152291B (zh) | 2023-09-12 | 2023-09-12 | 一种基于原对偶算法的非凸加权变分金属伪影去除方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117152291B (zh) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107590781A (zh) * | 2017-08-17 | 2018-01-16 | 天津大学 | 基于原始对偶算法的自适应加权tgv图像去模糊方法 |
| CN112017130A (zh) * | 2020-08-31 | 2020-12-01 | 郑州财经学院 | 新颖的基于自适应各向异性全变分正则化的图像复原方法 |
| CN112069919A (zh) * | 2020-08-17 | 2020-12-11 | 浙江工业大学 | 基于非凸低秩矩阵逼近和全变分正则化的高光谱图像去噪方法 |
| CN113469905A (zh) * | 2021-06-22 | 2021-10-01 | 吉林师范大学 | 一种基于复合正则化的低剂量ct投影域去噪方法 |
| WO2021252410A1 (en) * | 2020-06-08 | 2021-12-16 | University Of Utah Research Foundation | Reducing artifacts in computerized tomography scans |
| CN116452423A (zh) * | 2023-05-11 | 2023-07-18 | 燕山大学 | 一种同时稀疏角度ct重建及金属伪影高精度校正方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8447135B2 (en) * | 2011-03-02 | 2013-05-21 | Nanyang Technological University | Methods and systems for generating enhanced images using Euler's Elastica model |
-
2023
- 2023-09-12 CN CN202311169123.3A patent/CN117152291B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107590781A (zh) * | 2017-08-17 | 2018-01-16 | 天津大学 | 基于原始对偶算法的自适应加权tgv图像去模糊方法 |
| WO2021252410A1 (en) * | 2020-06-08 | 2021-12-16 | University Of Utah Research Foundation | Reducing artifacts in computerized tomography scans |
| CN112069919A (zh) * | 2020-08-17 | 2020-12-11 | 浙江工业大学 | 基于非凸低秩矩阵逼近和全变分正则化的高光谱图像去噪方法 |
| CN112017130A (zh) * | 2020-08-31 | 2020-12-01 | 郑州财经学院 | 新颖的基于自适应各向异性全变分正则化的图像复原方法 |
| CN113469905A (zh) * | 2021-06-22 | 2021-10-01 | 吉林师范大学 | 一种基于复合正则化的低剂量ct投影域去噪方法 |
| CN116452423A (zh) * | 2023-05-11 | 2023-07-18 | 燕山大学 | 一种同时稀疏角度ct重建及金属伪影高精度校正方法 |
Non-Patent Citations (4)
| Title |
|---|
| 基于变分正则化的混合泊松-高斯噪声图像去噪方法综述;常慧宾;张婕;;天津师范大学学报(自然科学版)(04);全文 * |
| 应用先验插值校正CT金属伪影;李铭;卢彦飞;袁刚;吴中毅;张涛;;液晶与显示(06);全文 * |
| 自适应加权全变分的低剂量CT统计迭代算法;何琳;张权;上官宏;张文;张鹏程;刘;桂志国;;计算机应用(10);全文 * |
| 自适应非凸稀疏正则化下自适应光学系统加性噪声的去除;张艳艳;陈苏婷;葛俊祥;万发雨;梅永;周晓彦;;物理学报(12);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117152291A (zh) | 2023-12-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yin et al. | Domain progressive 3D residual convolution network to improve low-dose CT imaging | |
| Zhang et al. | Regularization strategies in statistical image reconstruction of low‐dose x‐ray CT: A review | |
| Jin et al. | A model-based image reconstruction algorithm with simultaneous beam hardening correction for X-ray CT | |
| Chen et al. | Nonlocal prior Bayesian tomographic reconstruction | |
| EP2783344B1 (en) | Image domain de-noising | |
| Karimi et al. | Segmentation of artifacts and anatomy in CT metal artifact reduction | |
| US20070217566A1 (en) | System and Method For Image Reconstruction | |
| Ma et al. | Sinogram denoising via attention residual dense convolutional neural network for low-dose computed tomography | |
| US20230007835A1 (en) | Composition-guided post processing for x-ray images | |
| Li et al. | Low‐dose CT image denoising with improving WGAN and hybrid loss function | |
| Nakada et al. | Joint estimation of tissue types and linear attenuation coefficients for photon counting CT | |
| Ma et al. | Generalized Gibbs priors based positron emission tomography reconstruction | |
| Yu et al. | Weighted adaptive non-local dictionary for low-dose CT reconstruction | |
| He et al. | Downsampled imaging geometric modeling for accurate CT reconstruction via deep learning | |
| Shi et al. | Spectral CT reconstruction via low-rank representation and region-specific texture preserving Markov random field regularization | |
| Zhang et al. | Spectral CT image-domain material decomposition via sparsity residual prior and dictionary learning | |
| Chan et al. | An attention-based deep convolutional neural network for ultra-sparse-view CT reconstruction | |
| Guo et al. | Spectral2Spectral: Image-spectral similarity assisted deep spectral CT reconstruction without reference | |
| Ikuta et al. | A deep recurrent neural network with FISTA optimization for CT metal artifact reduction | |
| CN117152291B (zh) | 一种基于原对偶算法的非凸加权变分金属伪影去除方法 | |
| Humphries et al. | Superiorized method for metal artifact reduction | |
| Góes et al. | Poisson denoising under a Bayesian nonlocal approach using geodesic distances with low-dose CT applications | |
| Kang et al. | Higher-order regularization based image restoration with automatic regularization parameter selection | |
| Fu et al. | PWLS-PR: low-dose computed tomography image reconstruction using a patch-based regularization method based on the penalized weighted least squares total variation approach | |
| US20220198651A1 (en) | Multi-phase filter |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |