CN110992292B - 一种增强型低秩稀疏分解模型医学ct图像去噪方法 - Google Patents

一种增强型低秩稀疏分解模型医学ct图像去噪方法 Download PDF

Info

Publication number
CN110992292B
CN110992292B CN201911253091.9A CN201911253091A CN110992292B CN 110992292 B CN110992292 B CN 110992292B CN 201911253091 A CN201911253091 A CN 201911253091A CN 110992292 B CN110992292 B CN 110992292B
Authority
CN
China
Prior art keywords
image
rank
matrix
image block
low
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
Application number
CN201911253091.9A
Other languages
English (en)
Other versions
CN110992292A (zh
Inventor
夏克文
王斯洁
王莉
张江楠
周亚同
毛评生
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hebei University of Technology
Original Assignee
Hebei University of Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hebei University of Technology filed Critical Hebei University of Technology
Priority to CN201911253091.9A priority Critical patent/CN110992292B/zh
Publication of CN110992292A publication Critical patent/CN110992292A/zh
Application granted granted Critical
Publication of CN110992292B publication Critical patent/CN110992292B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2136—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on sparsity criteria, e.g. with an overcomplete basis
    • 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]
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30061—Lung

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Image Processing (AREA)

Abstract

本发明为一种增强型低秩稀疏分解模型医学CT图像去噪方法，该方法的步骤是：以计算出的噪声强度估计值的大小为依据确定搜索窗口、图像块矩阵中相似块的数量和大小、及迭代求解时最大迭代次数，遍历原始图像，进行非局部相似块匹配，将原始图像分成多个由非局部相似块组成的图像块矩阵；对医学CT原始图像D∈Rm×n采用加权Schatten p范数进行低秩矩阵估计，并添加联合约束的L1‑2TV正则化项，构建增强型低秩稀疏分解模型；将图像块矩阵依次输入模型中，利用交替方向乘子法进行迭代求解，得到相应图像块矩阵的低秩矩阵；聚合所有图像块矩阵对应的低秩矩阵，得到去噪后的干净图像。该方法能尽可能分离出更多的混合噪声，以获得更佳的医学CT图像去噪效果。

Description

一种增强型低秩稀疏分解模型医学CT图像去噪方法
技术领域
本发明属于医学图像去噪领域，尤其涉及一种增强型低秩稀疏分解模型医学CT图像去噪方法，该方法主要利用加权Schatten p范数和L1-2TV正则化约束来对医学CT图像进行去噪。
背景技术
图像的噪声是阻碍信息理解和分析的主要因素，图像去噪长久以来备受关注。计算机断层扫描(Computed Tomography，CT)图像是计算机辅助医疗的重要影像数据之一，CT图像承载了大量临床诊疗信息，可以有效辅助医生进行病情诊断、手术计划制定、术后治疗评估。但医学CT图像在采集、压缩、传输等过程中由于环境、设备等多方面因素的影响，不可避免地导致图像信号受到噪声的污染而质量下降。这首先会影响图像的人眼主观评价，更重要地是会对医学图像后续的处理和分析带来困难，如感兴趣区域分割、特征提取、识别分类等。在医学图像的处理和分析过程中，去噪是一个不可缺少的重要环节。
目前，针对医学CT图像的去噪方法局限性较大，常规的去噪技术主包括空间域去噪和变换域去噪两大类。空间域去噪方法，如中值滤波、均值滤波等，在空间域直接对图像的灰度值进行同一种平滑处理，忽略了各像素点的自身特征，使得去噪后图像的边缘信息模糊。变换域滤波技术的基本思想是将图像从空间域转换到变换域后，对变换域中的系数进行处理再反变换回空间域，从而达到去除噪声的目的。近年来，鲁棒主成分分析(RobustPrincipal Component Analysis，RPCA)模型广泛应用于图像去噪中，Cai等人在发表的论文“A Singular Value Thresholding Algorithm for Matrix Completion”中做出了具有开创性的工作，他们针对秩函数具有非凸性和不连续性的问题，将秩最小问题凸松弛为核范数最小化问题，并提出了奇异值阈值模型进行求解。该方法首先利用图像的低秩性构建先验模型，再通过求解矩阵的最小秩函数得到去除噪声的干净图像。但是，采用核范数进行最小秩问题的近似逼近求解，得到的结果往往是秩函数最小化问题的次优解，使得去噪图像的重建效果差。另外，随着噪声强度的增大，经过RPCA方法去噪后的图像往往会出现图像边界阶梯化的“油画”现象，造成图像边缘信息丢失严重的后果。
考虑到RPCA模型奇异值求解不准确、鲁棒性差等问题，本发明将采用加权Schatten p范数替代传统模型中的核范数进行最小秩函数估计，添加对稀疏信息的L1-2TV正则化联合约束，建立增强型低秩稀疏分解模型，并基于交替方向乘子法框架设计了新模型的迭代求解步骤。增强了传统模型的低秩矩阵估计能力，尽可能分离出了更多的混合噪声，以获得更佳的医学CT图像去噪效果。
发明内容
本发明的目的是克服现有技术的不足，提出了一种基于加权Schatten p范数和L1-2TV正则化约束的增强型低秩稀疏分解模型医学CT图像去噪方法。
为实现上述目的，本发明采取的技术方案是：
一种增强型低秩稀疏分解模型医学CT图像去噪方法，该方法的步骤是：
步骤1、对医学CT原始图像进行噪声水平估计，并以计算出的噪声强度估计值的大小为依据确定搜索窗口、图像块矩阵中相似块的数量和大小、及步骤3中迭代求解时最大迭代次数，遍历原始图像，进行非局部相似块匹配，将原始图像分成多个由非局部相似块组成的图像块矩阵；
步骤2、对医学CT原始图像D∈Rm×n采用加权Schatten p范数进行低秩矩阵估计，并添加联合约束的L1-2TV正则化项，构建增强型低秩稀疏分解模型为下式；
Figure BDA0002309577590000021
式中λ1和λ2均为平衡调节参数，L∈Rm×n表示去噪后干净的低秩矩阵，S∈Rm×n表示稀疏噪声矩阵，
Figure BDA0002309577590000022
为加权Schatten p范数，0＜p≤1；||·||1表示L1范数；
Figure BDA0002309577590000023
为L1-2TV正则化项，
Figure BDA0002309577590000024
表示各项异性TV，
Figure BDA0002309577590000025
表示各项同性TV,α∈[0,1]为正则化参数；
步骤3、将步骤1得到的图像块矩阵依次输入增强型低秩稀疏分解模型中，利用交替方向乘子法进行迭代求解，得到相应图像块矩阵的低秩矩阵；
步骤4、聚合所有图像块矩阵对应的低秩矩阵，得到去噪后的干净图像。
所述步骤1具体包括以下步骤：
步骤1.1：获取原始图像D∈Rm×n，它的秩记为rD，选取t＝3rD/5，对图像D进行噪声估计，计算噪声强度估计值σn；
步骤1.1.1：对原始图像D进行奇异值分解运算，计算出后t个奇异值的平均值为PD；
步骤1.1.2：对原始图像D添加噪声强度σD1＝30的零均值加性高斯白噪声，得到图像D1，再对D1进行奇异值分解运算，计算出后t个奇异值的平均值为PD1；
步骤1.1.3：对原始图像D添加噪声强度σD2＝60的零均值加性高斯白噪声，得到图像D2，再对D2进行奇异值分解运算，计算出后t个奇异值的平均值为PD2；
步骤1.1.4：根据下述公式计算出原始图像D的噪声强度估计值σn；
Figure BDA0002309577590000026
步骤1.2：根据噪声强度估计值σn，设定搜索窗口和当前图像块的大小，按步长step遍历整幅图像D，对图像进行非局部相似块匹配操作，得到图像块矩阵Dj；
步骤1.2.1：以欧式距离d为匹配准则，计算当前图像块Y(i,j)与搜索窗口之内其他图像块Z(i,j)之间的距离，计算公式为：
Figure BDA0002309577590000031
步骤1.2.2：对计算结果d从小到大进行排序，依据噪声强度估计值σn的大小设定相应数量的图像块作为相似块进行匹配，再令匹配的相似图像块堆叠为一个图像块矩阵，记为Dj，j表示得到第j个图像块矩阵。
所述步骤3具体包括如下步骤：
步骤3.1：将相似图像块堆叠成的图像块矩阵Dj输入增强型低秩稀疏分解模型，得：
Figure BDA0002309577590000032
式中，Lj和Sj分别表示图像块矩阵Dj对应的图像块低秩矩阵和图像块稀疏矩阵；
步骤3.2：对增强型低秩稀疏分解模型引入中间变量J,X∈Rm×n，得到：
Figure BDA0002309577590000033
步骤3.3：引入拉格朗日乘子，将上述公式转换为无约束优化问题：
Figure BDA0002309577590000034
式中，μ1,μ2,μ3＞0为惩罚系数，Y1,Y2,Y3,Y4为拉格朗日乘子；||·||2 F表示Frobenius范数，Xani和Xiso分别表示X的各向异性和各项同性，
Figure BDA0002309577590000037
和
Figure BDA0002309577590000038
分别表示Lj的各向异性和各项同性，带尖角的字母指相应变量最后求得的最优解；
步骤3.4：采用交替方向乘子法对步骤3.3的式子进行迭代求解，迭代次数记为k，k从0开始取值，最大迭代次数记为kmax，kmax的取值根据噪声强度估计值σn来设定；
步骤3.4.1：固定其余变量，更新X，更新规则为：
Figure BDA0002309577590000035
Figure BDA0002309577590000036
式中，Sτ[g]表示软阈值算子，Sτ[g]＝sgn(g)max(|g|-τ,0)，g是自变量，τ是标量；
步骤3.4.2：固定其余变量，更新Sj，更新规则为：
Figure BDA0002309577590000041
步骤3.4.3：固定其余变量，更新J，更新规则为：
Figure BDA0002309577590000042
步骤3.4.4：固定其余变量，更新Lj，更新规则为：
Figure BDA0002309577590000043
步骤3.4.5：对所有拉格朗日乘子进行更新，更新规则为：
Y1＝Y1+γ1(Lj+Sj-Dj)
Y2＝Y2+γ2(Lj-J)
Y3＝Y3+γ3(Ljani-Xani)
Figure BDA0002309577590000044
式中，γ1,γ2,γ3,γ4分别为拉格朗日乘子Y1,Y2,Y3,Y4的更新步长；
步骤3.4.6：迭代次数k＝k+1；
步骤3.4.7：进行终止条件的判决，当k＞kmax时，或者满足条件||Lj (k)-Lj (k+1)||2≤10-6||Lj (k)||2时(||·||2表示L2范数)，停止迭代；否则，返回步骤3.4.1；
步骤3.4.8：输出当前图像块矩阵Dj对应的低秩矩阵Lj。
上述方法中，当噪声强度估计值σn≤30时，将搜索窗大小设置为20×20，其余情况搜索窗大小均设置为30×30，当噪声强度估计值分别在σn≤20、20＜σn≤40和σn＞40范围内时，则相似块的大小依次设置为6×6、7×7和8×8，一个图像块矩阵中相似块的数量对应的分别为70、90和120个，最大迭代次数kmax相应分别设置为8、12和14。
与现有技术相比，本发明的优点和有益效果在于：
(1)根据图像本身固有的非局部自相关性，对原始图像进行非局部相似块搜索和匹配，将得到的相似块组成多个图像块矩阵，再将一个个图像块矩阵输入增强型低秩稀疏分解模型进行去噪处理，最后将去噪后的图像块矩阵聚合，恢复为完整图像，可以更好地发挥增强型低秩稀疏分解模型的作用效果，增强去噪性能，并提升方法整体的处理效率；
(2)本发明在基本RPCA模型的基础上，采用加权Schatten p范数代替核范数进行低秩矩阵逼近，可以有效提高对矩阵最小秩的估计精度，在增强低秩性的同时，减少了核范数最小化带来的信息损失；另外，添加L1-2TV正则化项进行联合稀疏约束，可以消除经典TV正则化的阶梯效应，更加逼近矩阵非零元素的稀疏性，进一步去除包含脉冲噪声、死线噪声、条带噪声及混合噪声在内的多种噪声，同时有效保留图像的细节和边缘信息；
(3)本发明的方法利用原始图像的非局部自相关性和低秩稀疏先验信息，能够有效去除医学CT图像中存在的噪声，同时保留其中有助于诊断的重要信息，弥补了传统去噪方法带来的细节信息丢失和边缘模糊等不足。观察附图2至附图6可以发现，RPCA、K-SVD和BM3D方法并不能较好地去除噪声，可以看到去噪后图像的背景中残存有部分的噪声干扰，并且在图像恢复的过程中，噪声和肺部组织结构发生了黏连，误把噪声恢复为了图像的主要信息部分。肺部组织的边缘信息稍显模糊，存在一定程度的“油画”现象。相比之下，本发明方法去噪处理后的图像肺部有用信息与背景区分显著，噪声信息被有效去除掉，肺实质保留完整，边缘部分的恢复效果良好。通过附图7对4种去噪方法进行定量对比分析，本发明方法的峰值信噪比，相比RPCA、K-SVD和BM3D方法分别提高了11％、3.9％和1.7％。
(4)本发明方法在RPCA模型的第一部分将核范数替换为加权schatten-p范数，同时用二维L1-2TV正则化约束项，一起对低秩矩阵进行联合约束，对改进后的模型基于交替方向乘子法思想进行求解，将整个优化问题分为四个子问题，逐个进行迭代求解。该方法能够直接应用于医学CT图像去噪中，针对真实的医学CT图像自身所带的噪声进行去除，具有在真实复杂的自然捕捉图像(图像噪声情况复杂，成分多，噪声情况不可预期)上的可行性和有效性。
(5)本发明方法设计了对整个图像进行分块和聚合步骤后，可以增强新模型的约束效果，提升去噪性能。
附图说明
图1为本发明的方法流程图；
图2为原始含噪肺部CT图像及其局部放大图像；
图3为RPCA方法得到的去噪图像(a)及其局部放大图像(b)；
图4为K奇异值分解方法得到的去噪图像(a)及其局部放大图像(b)；
图5为块匹配三维滤波方法得到的去噪图像(a)及其局部放大图像(b)；
图6为本发明方法得到的去噪图像(a)及其局部放大图像(b)；
图7为RPCA方法、K奇异值分解方法、块匹配三维滤波方法和本发明方法的去噪峰值信噪比曲线图。
具体实施方式
下面结合附图及实施例进一步解释本发明，但并不以此作为对本申请保护范围的限定。
图1为本发明医学CT图像去噪方法的流程图，展示了原始图像从输入到输出去噪后干净图像的完整流程。本发明是一种基于加权Schatten p范数和L1-2TV正则化约束的增强型低秩稀疏分解模型医学CT图像去噪方法的步骤是：
步骤1、对医学CT原始图像进行噪声水平估计，并以计算出的噪声强度的大小为依据确定搜索窗口尺寸、图像块矩阵中相似块(图像块)的大小和数量及步骤3中迭代求解时最大迭代次数，遍历原始图像，进行非局部相似块匹配，将原始图像分成多个由非局部相似块组成的图像块矩阵；
步骤2、对医学CT原始图像D∈Rm×n采用加权Schatten p范数进行低秩矩阵估计，并添加联合约束的L1-2TV正则化项，构建增强型低秩稀疏分解模型为下式；
Figure BDA0002309577590000061
式中λ1和λ2均为平衡调节参数，L∈Rm×n表示去噪后干净的低秩矩阵，S∈Rm×n表示稀疏噪声矩阵，
Figure BDA0002309577590000062
为加权Schatten p范数，0＜p≤1；||·||1表示L1范数；
Figure BDA0002309577590000063
为L1-2TV正则化项，
Figure BDA0002309577590000064
表示各项异性TV，
Figure BDA0002309577590000065
表示各项同性TV；
步骤3、将步骤1得到的图像块矩阵依次输入增强型低秩稀疏分解模型中，利用交替方向乘子法进行迭代求解，得到相应图像块矩阵的低秩矩阵；
步骤4、聚合所有图像块矩阵对应的低秩矩阵，得到去噪后的干净图像。
进一步地，所述步骤1具体包括如下步骤：
步骤1.1：获取原始图像D∈Rm×n，它的秩记为rD，选取t＝3rD/5。对图像D进行噪声估计，计算噪声强度σn；
步骤1.1.1：对原始图像D进行奇异值分解运算，计算出后t个奇异值的平均值为PD；
步骤1.1.2：对原始图像D添加噪声强度σD1＝30的零均值加性高斯白噪声，得到图像D1，再对D1进行奇异值分解运算，计算出后t个奇异值的平均值为PD1；
步骤1.1.3：对原始图像D添加噪声强度σD2＝60的零均值加性高斯白噪声，得到图像D2，再对D2进行奇异值分解运算，计算出后t个奇异值的平均值为PD2；
步骤1.1.4：根据下述公式计算出原始图像的噪声强度估计值σn；
Figure BDA0002309577590000066
步骤1.2：根据噪声强度σn的不同，设定大小不同的搜索窗口c×c和当前图像块l×l(搜索窗口大小根据图像和噪声强度大小设置)，当噪声强度σn≤30时，将搜索窗大小设置为20×20，其余情况均设置为30×30，这样可以加快算法的运行速度，当前图像块的大小依次设置为6×6，7×7，8×8分别对应于噪声强度为σn≤20，20＜σn≤40和σn＞40，按统一步长step遍历整幅图像D，对图像进行非局部相似块匹配操作，得到图像块矩阵Dj；
步骤1.2.1：以欧式距离d为匹配准则，计算当前图像块Y(i,j)与搜索窗口之内其他图像块Z(i,j)之间的距离，计算公式为：
Figure BDA0002309577590000067
步骤1.2.2：对计算结果d从小到大进行排序，依据噪声强度σn的大小取相应数量的图像块作为相似块进行匹配，当噪声强度为σn≤20，20＜σn≤40和σn＞40时分别取70，90和120个图像块作为相似块进行分组，再令同组内的相似块堆叠为图像块矩阵，记为Dj(j表示得到的相似图像块矩阵的个数)。
进一步地，所述步骤2具体包括如下步骤：
步骤2.1：将原始含噪图像表示为低秩图像与稀疏图像之和，即：D＝L+S，其中L∈Rm×n表示去噪后干净的低秩矩阵，S∈Rm×n表示稀疏噪声矩阵。本发明采用的基础模型如下所示：
Figure BDA0002309577590000071
式中，λ为低秩矩阵和稀疏矩阵的平衡调节参数，||·||*表示核范数，||·||1表示L1范数；
步骤2.2：加权Schatten p范数定义为：
Figure BDA0002309577590000072
式中，0＜p≤1，r＝min{m,n}，σi表示矩阵L的第i个奇异值，ωi为非负权重向量，表示对不同大小的奇异值做不同程度的处理，定义为对应矩阵L的奇异值的倒数，即：
Figure BDA0002309577590000073
式中ε＝10-6以防出现分母为零的情况；
步骤2.3：全变差正则化方法可以有效保护图像的边界信息，而各向异性TV和各项同性TV之差形式的L1-2TV正则化能够更精确地刻画图像的稀疏先验信息，其定义为：
Figure BDA0002309577590000074
式中，
Figure BDA0002309577590000075
表示各项异性TV，
Figure BDA0002309577590000076
表示各项同性TV，α∈[0,1]为正则化参数，当α＝0时，L1-2TV退化为经典的各向异性TV；
步骤2.4：采用加权Schatten p范数代替核范数，并加入二维L1-2TV正则化项，构造增强型低秩稀疏分解模型：
Figure BDA0002309577590000077
式中的λ1和λ2分别表示稀疏矩阵和低秩矩阵的平衡调节参数。
进一步地，所述步骤3具体包括如下步骤：
步骤3.1：将相似图像块堆叠成的图像块矩阵Dj输入增强型低秩稀疏分解模型，可得：
Figure BDA0002309577590000078
式中，Lj和Sj分别表示图像块矩阵Dj对应的图像块低秩矩阵和图像块稀疏矩阵；
步骤3.2：对增强型低秩稀疏分解模型引入中间变量J,X∈Rm×n，得到：
Figure BDA0002309577590000079
步骤3.3：引入拉格朗日乘子，将上述公式转换为无约束优化问题：
Figure BDA0002309577590000081
式中，μ1,μ2,μ3＞0为惩罚系数，Y1,Y2,Y3,Y4为拉格朗日乘子，
Figure BDA0002309577590000082
表示Frobenius范数，Xani和Xiso分别表示X的各向异性和各项同性，Ljani和Ljiso分别表示Lj的各向异性和各项同性，带尖角的字母指相应变量最后求得的最优解；
步骤3.4：采用交替方向乘子法对上式进行迭代求解，迭代次数记为k，k从0开始取值，最大迭代次数记为kmax，将其值设置为8，12和14，分别对应于噪声强度σn≤20，20＜σn≤40和σn＞40；
步骤3.4.1：固定其余变量，更新X，更新规则为：
Figure BDA0002309577590000083
Figure BDA0002309577590000084
式中，Sτ[g]表示软阈值算子，Sτ[g]＝sgn(g)max(|g|-τ,0)，g是自变量，τ是标量；
步骤3.4.2：固定其余变量，更新Sj，更新规则为：
Figure BDA0002309577590000085
步骤3.4.3：固定其余变量，更新J，更新规则为：
Figure BDA0002309577590000086
为了求解J，首先对矩阵J+Lj进行奇异值分解运算，即：
SVD(J+Lj)＝UΣVΤ,Σ＝diag(σ1,…,σr)；
其次，计算此时对应矩阵J+Lj的非负权重向量：
Figure BDA0002309577590000087
最后，得到J的最佳逼近结果为：
Figure BDA0002309577590000088
其中Δ＝diag(δ1，…,δr)，δi是如下问题的优化解：
Figure BDA0002309577590000089
步骤3.4.4：固定其余变量，更新Lj，更新规则为：
Figure BDA0002309577590000091
对函数f(Lj)求解Lj的偏导数，并令结果为零，则有：
Figure BDA0002309577590000092
对上式采用共轭梯度(PCG)算法进行求解，即可得到Lj的最佳结果；
步骤3.4.5：对所有拉格朗日乘子进行更新，更新规则为：
Y1＝Y1+γ1(Lj+Sj-Dj)
Y2＝Y2+γ2(Lj-J)
Y3＝Y3+γ3(Ljani-Xani)
Figure BDA0002309577590000093
式中，γ1,γ2,γ3,γ4分别为拉格朗日乘子Y1,Y2,Y3,Y4的更新步长；
步骤3.4.6：迭代次数k＝k+1；
步骤3.4.7：进行终止条件的判决，当k＞kmax时，或者满足条件||Lj (k)-Lj (k+1)||2≤10-6||Lj (k)||2时(||·||2表示L2范数)，停止迭代；否则，返回步骤3.4.1；
步骤3.4.8：输出当前图像块矩阵Dj对应的低秩矩阵Lj。
进一步地，所述步骤4具体包括如下步骤：按步骤1的输出顺序，将所有图像块矩阵Dj执行步骤3的操作，得到对应输出的所有低秩矩阵Lj，依次聚合(把这些得到的图像块按照原来的位置放回去，得到原来的完整图像)，最终获得原始图像D∈Rm×n对应的去噪后图像L∈Rm×n。
实施例1
本实施例选用LIDC/IDRI肺部CT图像数据库作为图像数据来源。采集肺部CT图像数据库中的图像数据(图像大小为512×512)，用D＝L+S表示，其中D∈Rm×n代表未经处理的原始图像数据，L∈Rm×n表示去噪后的低秩图像数据，S∈Rm×n表示稀疏噪声数据。
步骤1：对医学CT原始图像D∈Rm×n进行噪声估计，并根据计算出的噪声强度σn的大小，遍历原始图像，进行非局部相似块匹配，将原始图像分成多个由非局部相似块组成的图像块矩阵Dj：
步骤1.1：获取原始图像D∈Rm×n，它的秩记为rD，选取t＝3rD/5。对图像D进行噪声估计，计算噪声强度σn；
步骤1.1.1：对图像D进行奇异值分解运算，计算出后t个奇异值的平均值为PD；
步骤1.1.2：对图像D添加噪声强度σD1＝30的零均值加性高斯白噪声，得到图像D1。对D1进行奇异值分解运算，计算出后t个奇异值的平均值为PD1；
步骤1.1.3：对图像D添加噪声强度σD2＝60的零均值加性高斯白噪声，得到图像D2。对D2进行奇异值分解运算，计算出后t个奇异值的平均值为PD2；
步骤1.1.4：根据下述公式计算出本实施例输入的原始图像的噪声强度估计值为σn＝24。
Figure BDA0002309577590000101
步骤1.2：根据噪声强度σn＝24，设置搜索窗口大小为20×20，当前图像块大小为7×7，步长step统一设置为3，遍历整幅图像D，对图像进行非局部相似块匹配操作，得到图像块矩阵Dj；
步骤1.2.1：以欧式距离d为匹配准则，计算当前图像块Y(i,j)与搜索窗口之内其他图像块Z(i,j)之间的距离，计算公式为：
Figure BDA0002309577590000102
步骤1.2.2：对计算结果d从小到大进行排序，根据噪声强度σn＝24，取90个图像块作为相似块进行分组，再令同组内的相似块堆叠为图像块矩阵，记为Dj(j表示得到第j个图像块矩阵)。
步骤2：构建增强型低秩稀疏分解模型：
步骤2.1：将含噪图像表示为低秩图像与稀疏图像之和，即：D＝L+S，其中L∈Rm×n表示去噪后干净的低秩矩阵，S∈Rm×n表示稀疏噪声矩阵，本发明采用的基础模型如下所示：
Figure BDA0002309577590000103
式中，λ为低秩矩阵和稀疏矩阵的平衡调节参数，||·||*表示核范数，||·||1表示L1范数；
步骤2.2：加权Schatten p范数定义为：
Figure BDA0002309577590000104
式中，0＜p≤1，r＝min{m,n}，σi表示矩阵L的第i个奇异值，ωi为非负权重向量，i为1～r范围内的整数，表示对不同大小的奇异值做不同程度的处理，其定义为对应矩阵L的奇异值的倒数，即：
Figure BDA0002309577590000105
式中ε＝10-6以防出现分母为零的情况；
步骤2.3：全变差正则化方法可以有效保护图像的边界信息，而各向异性TV和各项同性TV之差形式的L1-2TV正则化能够更精确地刻画图像的稀疏先验信息，其定义为：
Figure BDA00023095775900001110
式中，
Figure BDA0002309577590000111
表示各项异性TV，
Figure BDA0002309577590000112
表示各项同性TV，α∈[0,1]为正则化参数，当α＝0时，L1-2TV退化为经典的各向异性TV；
步骤2.4：采用加权Schatten p范数代替核范数，并加入L1-2TV正则化项，构造增强型低秩稀疏分解模型：
Figure BDA0002309577590000113
式中的λ1和λ2分别表示稀疏矩阵和低秩矩阵的平衡调节参数。
步骤3：将Dj输入模型，并采用交替方向乘子法进行求解，得到图像块矩阵对应的低秩矩阵Lj：
步骤3.1：将图像块矩阵Dj输入上述模型，可得
Figure BDA0002309577590000114
式中，Lj和Sj分别表示图像块矩阵Dj对应的图像块低秩矩阵和图像块稀疏矩阵；
步骤3.2：对新建模型引入中间变量J,X∈Rm×n，得到：
Figure BDA0002309577590000115
步骤3.3：引入拉格朗日乘子，将上述公式转换为无约束优化问题：
Figure BDA0002309577590000116
式中，μ1,μ2,μ3＞0为惩罚系数，Y1,Y2,Y3,Y4为拉格朗日乘子，
Figure BDA0002309577590000117
表示Frobenius范数，Xani和Xiso分别表示X的各向异性和各项同性，
Figure BDA00023095775900001111
和
Figure BDA00023095775900001112
分别表示Lj的各向异性和各项同性，带尖角的字母指相应变量最后求得的最优解，这里设置μ1＝μ2＝μ3＝0.5，
Figure BDA0002309577590000118
Figure BDA0002309577590000119
本实施例中，m＝n＝512，σn＝24；
步骤3.4：采用交替方向乘子法对上式进行迭代求解，迭代次数记为k，k从0开始取值，根据噪声强度σn＝24，最大迭代次数设置为kmax＝12；
步骤3.4.1：固定其余变量，更新X，更新规则为：
Figure BDA0002309577590000121
Figure BDA0002309577590000122
式中，Sτ[g]表示软阈值算子，Sτ[g]＝sgn(g)max(|g|-τ,0)(g是自变量，τ是标量)；
步骤3.4.2：固定其余变量，更新Sj，更新规则为：
Figure BDA0002309577590000123
步骤3.4.3：固定其余变量，更新J，更新规则为：
Figure BDA0002309577590000124
为了求解上式，首先对矩阵J+Lj进行奇异值分解运算，即：
SVD(J+Lj)＝UΣVΤ,Σ＝diag(σ1,…,σr)；
其次，计算此时对应矩阵J+Lj的非负权重向量：
Figure BDA0002309577590000125
最后，得到J的最佳逼近结果为：
Figure BDA0002309577590000126
其中Δ＝diag(δ1，…,δr)，δi是如下问题的优化解：
Figure BDA0002309577590000127
步骤3.4.4：固定其余变量，更新Lj，更新规则为：
Figure BDA0002309577590000128
对函数f(Lj)求解Lj的偏导数，并令结果为零，则有：
Figure BDA0002309577590000129
对上式采用共轭梯度(PCG)算法进行求解，即可得到Lj的最佳结果；
步骤3.4.5：对所有拉格朗日乘子进行更新，更新规则为：
Y1＝Y1+γ1(Lj+Sj-Dj)
Y2＝Y2+γ2(Lj-J)
Y3＝Y3+γ3(Ljani-Xani)
Figure BDA0002309577590000131
式中，γ1,γ2,γ3,γ4分别为拉格朗日乘子Y1,Y2,Y3,Y4的更新步长，这里设置γ1＝γ2＝γ3＝γ4＝1；
步骤3.4.6：迭代次数k＝k+1；
步骤3.4.7：进行终止条件的判决，当k＞kmax时，或者满足条件||Lj (k)-Lj (k+1)||2≤10-6||Lj (k)||2时(||·||2表示L2范数)，停止迭代；否则，返回步骤3.4.1；
步骤3.4.8：输出当前图像块矩阵Dj对应的低秩矩阵Lj。
步骤4：按步骤1的输出顺序，将所有图像块矩阵Dj执行步骤3的操作，得到对应输出的所有低秩矩阵Lj，依次聚合，最终获得原始图像D∈Rm×n对应的去噪后图像L∈Rm×n。
在本发明中，对带有噪声的原始肺部CT图像进行去噪，并与传统RPCA方法、块匹配三维滤波方法、K奇异值分解方法进行对比，分析不同方法的去噪性能。
图2为原始含噪图像及其局部放大图像(b是a的局部放大图)，图3、图4、图5和图6分别展示了RPCA方法、块匹配三维滤波方法、K奇异值分解方法和本发明方法去噪后的干净肺部CT图像及其局部放大图像。相比之下，经过本发明去噪处理后的图像，肺部有用信息与背景区分显著(白色的是肺组织，黑色的是背景。边缘清晰、背景中白色小细纹少的就证明效果好。)，噪声信息被有效去除掉，肺实质保留完整，边缘部分恢复效果良好。通过人眼的主观评价，可以判定本发明方法的去噪效果最佳。
图7展示了RPCA方法、K奇异值分解(K-SVD)方法、块匹配三维滤波(BM3D)方法和本发明(Proposed)方法的去噪峰值信噪比曲线图。横坐标为噪声强度，纵坐标为去噪处理后图像的峰值信噪比(PSNR)。可以发现，本发明方法具有最高的峰值信噪比数值，且随着噪声强度的增加鲁棒性较好(横坐标为噪声强度，纵坐标为去噪处理后图像的峰值信噪比(PSNR)。去噪方法均对噪声大小没有要求，只不过有的方法随着噪声的变大能力下降，称为鲁棒性差)。在客观指标评价下，本发明方法仍然具有最佳的去噪性能。
以上所述实施例仅旨在便于对本发明的理解，并不用以限制本发明，凡在本发明的精神和原则之内所做的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
本发明未述及之处适用于现有技术。

Claims (4)

1.一种增强型低秩稀疏分解模型医学CT图像去噪方法，该方法的步骤是：
步骤1、对医学CT原始图像进行噪声水平估计，并以计算出的噪声强度估计值的大小为依据确定搜索窗口、图像块矩阵中相似块的数量和大小、及步骤3中迭代求解时最大迭代次数，遍历原始图像，进行非局部相似块匹配，将原始图像分成多个由非局部相似块组成的图像块矩阵；
步骤2、对医学CT原始图像D∈Rm×n采用加权Schattenp范数进行低秩矩阵估计，并添加联合约束的L1-2TV正则化项，构建增强型低秩稀疏分解模型为下式；
Figure FDA0002309577580000011
式中λ1和λ2均为平衡调节参数，L∈Rm×n表示去噪后干净的低秩矩阵，S∈Rm×n表示稀疏噪声矩阵，
Figure FDA0002309577580000012
为加权Schattenp范数，0＜p≤1；||·||1表示L1范数；
Figure FDA0002309577580000013
为L1-2TV正则化项，
Figure FDA0002309577580000014
表示各项异性TV，
Figure FDA0002309577580000015
表示各项同性TV，α∈[0,1]为正则化参数；
步骤3、将步骤1得到的图像块矩阵依次输入增强型低秩稀疏分解模型中，利用交替方向乘子法进行迭代求解，得到相应图像块矩阵的低秩矩阵；
步骤4、聚合所有图像块矩阵对应的低秩矩阵，得到去噪后的干净图像。
2.根据权利要求1所述的方法，其特征在于，所述步骤1具体包括以下步骤：
步骤1.1：获取原始图像D∈Rm×n，它的秩记为rD，选取t＝3rD/5，对图像D进行噪声估计，计算噪声强度估计值σn；
步骤1.1.1：对原始图像D进行奇异值分解运算，计算出后t个奇异值的平均值为PD；
步骤1.1.2：对原始图像D添加噪声强度σD1＝30的零均值加性高斯白噪声，得到图像D1，再对D1进行奇异值分解运算，计算出后t个奇异值的平均值为PD1；
步骤1.1.3：对原始图像D添加噪声强度σD2＝60的零均值加性高斯白噪声，得到图像D2，再对D2进行奇异值分解运算，计算出后t个奇异值的平均值为PD2；
步骤1.1.4：根据下述公式计算出原始图像D的噪声强度估计值σn；
Figure FDA0002309577580000016
步骤1.2：根据噪声强度估计值σn，设定搜索窗口和当前图像块的大小，按步长step遍历整幅图像D，对图像进行非局部相似块匹配操作，得到图像块矩阵Dj；
步骤1.2.1：以欧式距离d为匹配准则，计算当前图像块Y(i,j)与搜索窗口之内其他图像块Z(i,j)之间的距离，计算公式为：
Figure FDA0002309577580000021
步骤1.2.2：对计算结果d从小到大进行排序，依据噪声强度估计值σn的大小设定相应数量的图像块作为相似块进行匹配，再令匹配的相似图像块堆叠为一个图像块矩阵，记为Dj，j表示得到第j个图像块矩阵。
3.根据权利要求1所述的方法，其特征在于，所述步骤3具体包括如下步骤：
步骤3.1：将相似图像块堆叠成的图像块矩阵Dj输入增强型低秩稀疏分解模型，得：
Figure FDA0002309577580000022
式中，Lj和Sj分别表示图像块矩阵Dj对应的图像块低秩矩阵和图像块稀疏矩阵；
步骤3.2：对增强型低秩稀疏分解模型引入中间变量J,X∈Rm×n，得到：
Figure FDA0002309577580000023
步骤3.3：引入拉格朗日乘子，将上述公式转换为无约束优化问题：
Figure FDA0002309577580000024
式中，μ1,μ2,μ3＞0为惩罚系数，Y1,Y2,Y3,Y4为拉格朗日乘子；
Figure FDA0002309577580000025
表示Frobenius范数，Xani和Xiso分别表示X的各向异性和各项同性，
Figure FDA0002309577580000026
和
Figure FDA0002309577580000027
分别表示Lj的各向异性和各项同性，带尖角的字母指相应变量最后求得的最优解；
步骤3.4：采用交替方向乘子法对步骤3.3的式子进行迭代求解，迭代次数记为k，k从0开始取值，最大迭代次数记为kmax，kmax的取值根据噪声强度估计值σn来设定；
步骤3.4.1：固定其余变量，更新X，更新规则为：
Figure FDA0002309577580000028
Figure FDA0002309577580000029
式中，Sτ[g]表示软阈值算子，Sτ[g]＝sgn(g)max(|g|-τ,0)，g是自变量，τ是标量；
步骤3.4.2：固定其余变量，更新Sj，更新规则为：
Figure FDA0002309577580000031
步骤3.4.3：固定其余变量，更新J，更新规则为：
Figure FDA0002309577580000032
步骤3.4.4：固定其余变量，更新Lj，更新规则为：
Figure FDA0002309577580000033
步骤3.4.5：对所有拉格朗日乘子进行更新，更新规则为：
Y1＝Y1+γ1(Lj+Sj-Dj)
Y2＝Y2+γ2(Lj-J)
Figure FDA0002309577580000034
Figure FDA0002309577580000035
式中，γ1,γ2,γ3,γ4分别为拉格朗日乘子Y1,Y2,Y3,Y4的更新步长；
步骤3.4.6：迭代次数k＝k+1；
步骤3.4.7：进行终止条件的判决，当k＞kmax时，或者满足条件||Lj (k)-Lj (k+1)||2≤10-6||Lj (k)||2时(||·||2表示L2范数)，停止迭代；否则，返回步骤3.4.1；
步骤3.4.8：输出当前图像块矩阵Dj对应的低秩矩阵Lj。
4.根据权利要求1所述的方法，其特征在于，
当噪声强度估计值σn≤30时，将搜索窗大小设置为20×20，其余情况搜索窗大小均设置为30×30，当噪声强度估计值分别在σn≤20、20＜σn≤40和σn＞40范围内时，则相似块的大小依次设置为6×6、7×7和8×8，一个图像块矩阵中相似块的数量对应的分别为70、90和120个，最大迭代次数kmax相应分别设置为8、12和14。
CN201911253091.9A 2019-12-09 2019-12-09 一种增强型低秩稀疏分解模型医学ct图像去噪方法 Active CN110992292B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911253091.9A CN110992292B (zh) 2019-12-09 2019-12-09 一种增强型低秩稀疏分解模型医学ct图像去噪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911253091.9A CN110992292B (zh) 2019-12-09 2019-12-09 一种增强型低秩稀疏分解模型医学ct图像去噪方法

Publications (2)

Publication Number Publication Date
CN110992292A CN110992292A (zh) 2020-04-10
CN110992292B true CN110992292B (zh) 2023-04-18

Family

ID=70091518

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911253091.9A Active CN110992292B (zh) 2019-12-09 2019-12-09 一种增强型低秩稀疏分解模型医学ct图像去噪方法

Country Status (1)

Country Link
CN (1) CN110992292B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112435175A (zh) * 2020-10-30 2021-03-02 西安交通大学 一种金相图像去噪方法及系统
CN112686814B (zh) * 2020-12-25 2023-03-03 国家电网有限公司 一种基于仿射低秩的图像去噪方法
CN112734875A (zh) * 2021-01-08 2021-04-30 金陵科技学院 基于非局部低秩正则化的图像重建方法
CN113009560B (zh) * 2021-03-23 2022-03-29 中国地质大学（武汉） 一种地震数据重建方法、装置、设备及存储介质
CN113331789A (zh) * 2021-05-31 2021-09-03 浙江杜比医疗科技有限公司 肿瘤细胞生长检测系统的成像方法
CN113378415B (zh) * 2021-08-12 2021-11-02 西南科技大学 基于局部和全局约束的多媒体数据自适应恢复方法和装置
CN113837958B (zh) * 2021-09-09 2023-08-04 南方医科大学 扩散加权图像去噪算法、介质及设备

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106709881A (zh) * 2016-12-14 2017-05-24 上海增容数据科技有限公司 一种基于非凸低秩矩阵分解的高光谱图像去噪方法
CN109064412A (zh) * 2018-06-20 2018-12-21 南京邮电大学 一种低秩图像的去噪方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106709881A (zh) * 2016-12-14 2017-05-24 上海增容数据科技有限公司 一种基于非凸低秩矩阵分解的高光谱图像去噪方法
CN109064412A (zh) * 2018-06-20 2018-12-21 南京邮电大学 一种低秩图像的去噪方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Yuan Xie.etc..Weighted Schatten p-Norm Minimization for Image Denoising and Background Subtraction.IEEE .2016,全文. *
张志伟 ; 马杰 ; 夏克文 ; 李昱乐 ; .一种应用于高阶数据修复的非负稀疏Tucker分解算法.光电子・激光.2017,(第07期),全文. *
蒋明峰等.基于加权 Schatten p 范数最小化的 磁共振图像重构方法研究.电子学报.2019,(第02期),全文. *

Also Published As

Publication number Publication date
CN110992292A (zh) 2020-04-10

Similar Documents

Publication Publication Date Title
CN110992292B (zh) 一种增强型低秩稀疏分解模型医学ct图像去噪方法
CN112200750B (zh) 一种超声图像去噪模型建立方法及超声图像去噪方法
Tian et al. Deep learning on image denoising: An overview
Hou et al. NLH: A blind pixel-level non-local method for real-world image denoising
Xu et al. A trilateral weighted sparse coding scheme for real-world image denoising
CN106952228B (zh) 基于图像非局部自相似性的单幅图像的超分辨率重建方法
Zoran et al. Scale invariance and noise in natural images
CN108932699B (zh) 基于变换域的三维匹配调和滤波图像去噪方法
Singh et al. A review of image fusion: Methods, applications and performance metrics
CN114820352A (zh) 一种高光谱图像去噪方法、装置及存储介质
Malladi et al. Image denoising using superpixel-based PCA
Guo et al. Agem: Solving linear inverse problems via deep priors and sampling
CN107292316A (zh) 一种基于稀疏表示的提升图像清晰度的方法
Kumar et al. Low rank poisson denoising (LRPD): A low rank approach using split bregman algorithm for poisson noise removal from images
CN112767271A (zh) 一种基于三维变分网络的高光谱图像深度降噪的方法
Xu et al. Blind image deblurring using group sparse representation
CN111915518A (zh) 基于三重低秩模型的高光谱图像去噪方法
Lin et al. Noise2Grad: Extract Image Noise to Denoise.
Antsiperov New Centre/Surround Retinex-like Method for Low-Count Image Reconstruction.
CN112801899B (zh) 基于互补结构感知的内外循环驱动图像盲去模糊方法和装置
CN107085839B (zh) 基于纹理增强与稀疏编码的sar图像降斑方法
Zhan et al. Nonlocal means image denoising with minimum MSE-based decay parameter adaptation
Jebur et al. Image denoising techniques: An overview
CN116862809A (zh) 一种低曝光条件下的图像增强方法
Ma et al. Edge-guided cnn for denoising images from portable ultrasound devices

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