CN109191391B - 一种衰减参数自适应非局部均值的图像降噪方法 - Google Patents

一种衰减参数自适应非局部均值的图像降噪方法 Download PDF

Info

Publication number
CN109191391B
CN109191391B CN201810898880.7A CN201810898880A CN109191391B CN 109191391 B CN109191391 B CN 109191391B CN 201810898880 A CN201810898880 A CN 201810898880A CN 109191391 B CN109191391 B CN 109191391B
Authority
CN
China
Prior art keywords
image
noise
attenuation parameter
optimal
value
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
CN201810898880.7A
Other languages
English (en)
Other versions
CN109191391A (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.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and 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 Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN201810898880.7A priority Critical patent/CN109191391B/zh
Publication of CN109191391A publication Critical patent/CN109191391A/zh
Application granted granted Critical
Publication of CN109191391B publication Critical patent/CN109191391B/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

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)

Abstract

本发明提供了一种衰减参数自适应非局部均值的图像降噪方法，涉及图像去噪领域。在对噪声图像进行非局部均值降噪处理中确定衰减参数与全局均方误差的关系，利用迭代算法获得衰减参数的最优值，并根据该值确定最优预降噪图像；然后利用最优预降噪图像估计降噪图像的像素级均方误差，再利用最速下降法获得像素级均方误差取得最小值时对应的最终降噪图像。本发明提供的图像降噪方法，可对受高斯噪声和斑点噪声污染的图像进行有效的恢复，与现有非局部均值方法相比，在图像噪声滤除、伪影消除和细节保留方面可获得更佳的效果。

Description

一种衰减参数自适应非局部均值的图像降噪方法
技术领域
本发明涉及图像分析与处理中的图像去噪领域，尤其涉及一种衰减参数自适应非局部均值的图像降噪方法。
背景技术
噪声会降低图像质量，并且使模式识别变得困难。降噪有助于后续图像处理任务，如分割、配准及可视化等处理任务，因此研究图像降噪方法具有重要理论价值和实际意义。
非局部均值(Nonlocal Means,NLM)是图像降噪的主流方法之一，在图像降噪领域引起了广泛关注。该方法最初是为消除高斯噪声而设计，其核心思想是利用图像在全局范围内的自相似性实现图像降噪，其中自相似性则利用像素邻域(即图像块)而非单个像素间的灰度差异进行评价。近年来，有学者将NLM推广至受斑点噪声污染图像(如医学超声图像、合成孔径雷达图像)的降噪处理中，提出了基于概率块(probabilistic patch-based,PPB)的滤波方法及优化的贝叶斯非局部均值(optimized Bayesian nonlocal means，OBNLM)方法，其取得了一定的效果。
NLM方法的恢复性能与衰减参数h有直接关系，衰减参数的自适应选择是一个极具挑战的难题。有学者提出了衰减参数全局自适应方法，这些方法采用Stein无偏风险估计(Stein’s unbiased risk estimate，SURE)来确定h，或简单地认为h与图像中噪声标准差σ之间呈线性关系，将h设计为h＝c·σ，其中c为常数。这些方法虽然能在一定程度上提高NLM的鲁棒性，但对整幅图像选择相同的衰减参数h很难保证对图像中的边缘区和平滑区皆取得良好的恢复效果。因此，研究衰减参数h的局部自适应方法非常重要。Duval等在NLM方法中引入了偏差-方差平衡原理，利用像素级SURE(pixel-wise SURE,PSURE)来确定每个像素的最优衰减参数，该方法易在图像边缘区域产生衰减参数的误估计，导致在一些边缘区域的恢复性能降低。Dore等通过在牛顿方法中嵌入Cp统计方法，提出了Cp-NLM方法，以确定最佳的衰减参数。该方法直接利用受噪声污染的图像来估计局部均方误差，当图像中噪声污染程度较严重时，因难以准确确定每个像素对应的衰减参数而导致图像恢复性能的不佳。
总体来说，现有非局部均值方法无法准确地确定衰减系数，因此会影响该方法对高斯噪声或斑点噪声污染图像进行降噪处理时的恢复性能，特别是在图像中噪声污染程度较高时，上述方法难以在有效抑制图像噪声的同时很好地保护图像细节信息。
发明内容
为了克服现有技术中的不足，本发明提供一种衰减参数自适应非局部均值的图像降噪方法，利用该方法可在图像处理中获取更为理想的衰减参数，从而有效提高图像处理过程中的噪声滤除、伪影规避和细节保留性能。
本发明是通过以下技术方案实现的：
一种衰减参数自适应非局部均值的图像降噪方法，所述方法步骤包括：
步骤S1、对噪声图像I进行非局部均值降噪处理，获得降噪图像μ，并确定在降噪处理中的衰减参数h与降噪图像μ对应的全局均方误差M[μ]的关系；
步骤S2、根据所述衰减参数h与所述全局均方误差M[μ]之间的关系，计算出所述衰减参数h的最优值hopt，并根据该最优值hopt确定最优预降噪图像
Figure GDA0002905287760000031
步骤S3、计算出所述最优预降噪图像
Figure GDA0002905287760000032
对应的方法噪声
Figure GDA0002905287760000033
并对该方法噪声
Figure GDA0002905287760000034
进行降噪处理，获得所述方法噪声
Figure GDA0002905287760000035
的残余图像细节图
Figure GDA0002905287760000036
步骤S4、利用所述最优预降噪图像
Figure GDA0002905287760000037
所述方法噪声
Figure GDA0002905287760000038
以及所述残余图像细节图
Figure GDA0002905287760000039
对所述噪声图像I的原始图像s和所述噪声图像I的噪声η进行估算，获得所述原始图像s和所述噪声η的估算值分别为
Figure GDA00029052877600000310
和
Figure GDA00029052877600000311
步骤S5、基于所述
Figure GDA00029052877600000312
和
Figure GDA00029052877600000313
计算出所述降噪图像μ的像素级均方误差，并求解获得当所述像素级均方误差取得最小值时对应的最终降噪图像。
在本发明的优选实施方式中，所述步骤S1中，确定所述衰减参数h与所述全局均方误差M[μ]的关系由以下公式确定：
Figure GDA00029052877600000314
Figure GDA00029052877600000315
Figure GDA00029052877600000316
在本发明的优选实施方式中，所述步骤S2中，利用迭代算法计算出所述衰减参数h的最优值hopt，并根据该最优值hopt确定最优预降噪图像
Figure GDA00029052877600000317
在本发明的优选实施方式中，所述步骤S3中，对所述方法噪声
Figure GDA00029052877600000318
进行降噪处理包括非局部均值方法降噪处理和均值滤波降噪处理。
在本发明的优选实施方式中，在求解获得所述最终降噪图像过程中，利用最速下降法求解获得当所述像素级均方误差取得最小值时对应的最终降噪图像；具体通过利用最速下降法获得所述降噪图像μ的每个像素(x,y)对应的最优衰减参数h′opt(x,y)，并根据该最优衰减参数h′opt(x,y)获得对应的最终降噪图像。
在本发明的优选实施方式中，在获得所述最优衰减参数h′opt(x,y)的过程中，首先将所述降噪图像μ的每个像素(x,y)对应的均方误差
Figure GDA0002905287760000041
估计为与该像素对应衰减系数h(x,y)相关的代价函数J[h(x,y)]，再利用最速下降法对J[h(x,y)]进行优化求解，从而获得所述噪声图像I的每个像素(x,y)对应的最优衰减参数h′opt(x,y)。
与现有技术相比，本发明提供的图像降噪方法可对受高斯噪声和斑点噪声污染的图像进行更有效的恢复，在图像噪声滤除、伪影消除和细节保留方面可获得更佳的效果。
附图说明
图1为本发明实施例中所述图像降噪方法的流程示意图。
图2是五幅标准图像，其中图(2a)是Lena图像，图(2b)是Barbara图像,图(2c)是Boat图像，图(2d)是Peppers图像，图(2e)是Airplane图像。
图3是各方法对受方差为40的高斯噪声污染的Boat图像进行去噪处理后的结果；其中，图(3a)为噪声图像，图(3b)为原始图像，图(3c)为NLM方法的去噪结果，图(3d)为SURE方法的去噪结果，图(3e)为PSURE方法的去噪结果，图(3f)为Cp-NLM方法的去噪结果，图(3g)为PNLM方法的去噪结果。
图4是对图3中选取的感兴趣区域的放大图；其中，图(4a)为噪声图像，图(4b)为感兴趣区域图像，图(4c)为NLM方法的去噪结果，图(4d)为SURE方法的去噪结果，图(4e)为PSURE方法的去噪结果，图(4f)为Cp-NLM方法的去噪结果，图(4g)为PNLM方法的去噪结果。
图5是各方法对受方差为0.4的斑点噪声污染的Peppers图像进行降噪处理后对应结果中感兴趣区的放大图；其中，图(5a)为噪声图像，图(5b)为感兴趣区图像，图(5c)为NLM方法的去噪结果，图(5d)为PPB方法的去噪结果，图(5e)为OBNLM方法的去噪结果，图(5f)为Cp-NLM方法的去噪结果，图(5g)为PNLM方法的去噪结果。
具体实施方式
为了便于本领域技术人员的理解，下面结合附图对本发明作进一步的描述。
参见图1，一种衰减参数自适应非局部均值的图像降噪方法，包括：
步骤S1、对噪声图像I进行非局部均值降噪处理，获得降噪图像μ，并确定在降噪处理中的衰减参数h与降噪图像μ对应的全局均方误差M[μ]的关系；
步骤S2、根据所述衰减参数h与所述全局均方误差M[μ]之间的关系，计算出所述衰减参数h的最优值hopt(即hopt为最优衰减参数)，并根据该最优值hopt确定最优预降噪图像
Figure GDA0002905287760000051
步骤S3、计算出所述最优预降噪图像
Figure GDA0002905287760000052
对应的方法噪声
Figure GDA0002905287760000053
并对该方法噪声
Figure GDA0002905287760000054
进行降噪处理，获得所述方法噪声
Figure GDA0002905287760000055
中的残余图像细节图
Figure GDA0002905287760000056
步骤S4、利用所述最优预降噪图像
Figure GDA0002905287760000061
所述方法噪声
Figure GDA0002905287760000062
以及所述残余图像细节图
Figure GDA0002905287760000063
对所述噪声图像I的原始图像s和所述噪声图像I的噪声η进行估算，获得所述原始图像s和所述噪声η的估算值分别为
Figure GDA0002905287760000064
和
Figure GDA0002905287760000065
步骤S5、基于所述
Figure GDA0002905287760000066
和
Figure GDA0002905287760000067
计算出所述降噪图像μ的像素级均方误差，并求解获得当所述像素级均方误差取得最小值时对应的最终降噪图像。其中，所述降噪图像μ的像素级均方误差，是指基于所述降噪图像μ的每个像素(x,y)对应的均方误差。
作为本发明的优选实施例，所述步骤S2中，利用迭代算法计算出所述衰减参数h的最优值hopt，并根据该最优值hopt确定最优预降噪图像
Figure GDA0002905287760000068
所述步骤S3中，对所述方法噪声
Figure GDA0002905287760000069
进行降噪处理包括非局部均值方法降噪处理和均值滤波降噪处理；
所述步骤S5中，在求解获得所述最终降噪图像过程中，利用最速下降法求解获得当所述像素级均方误差取得最小值时对应的最终降噪图像，在求解获得所述最终降噪图像过程中，具体通过利用最速下降法获得所述降噪图像μ的每个像素(x,y)对应的最优衰减参数h′opt(x,y)，并根据该最优衰减参数h′opt(x,y)获得对应的最终降噪图像。
本发明实施例的图像降噪方法的流程如附图1所示。为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。
本发明实施例的衰减参数自适应非局部均值的图像降噪方法，其主要包括两大处理环节：处理环节一，获得衰减参数h的最优值hopt(即最优衰减参数hopt)，并根据最优值hopt确定最优预降噪图像
Figure GDA00029052877600000610
处理环节二，获得所述降噪图像μ的像素级均方误差，并求解获得当所述像素级均方误差取得最小值时对应的最终降噪图像，即获得所述降噪图像μ的每个像素(x,y)对应的最优衰减参数h′opt(x,y)，并根据最优衰减参数h′opt(x,y)获得最终降噪图像。
在对上述的两大处理节作详细的说明前，为便于描述和理解，本发明实施例中，假定受污染的噪声图像I为I＝s+sγu＝s+η,其中s为噪声图像I的原始图像，η为噪声图像I中的噪声，γ为常数，γ＝0表示图像中的噪声为高斯噪声，γ＝1则表示图像中的噪声为斑点噪声，μ为零均值高斯随机变量。对于噪声图像I，假定该噪声图像I经非局部均值降噪后获得的降噪图像为μ，相应的残余分量为v，则有I＝s+η＝μ+v。另外，在本发明中，若无特别的说明，本发明中相同的字母或符号代表着相同的含义，或相同的字母或符号表示同一参数。
以下对所述的处理环节一作详细的说明：
对于降噪图像μ，其像素(x,y)的对应灰度值μ(x,y)为：
Figure GDA0002905287760000071
式(4)中，Ω代表搜索窗口，ω(x,y,p,q)表示(x,y)和(p,q)两像素的相似性，定义为：
Figure GDA0002905287760000072
式(5)中，G代表高斯核，*代表卷积运算，P(I(x,y))和P(I(p,q))分别表示以(x,y)和(p,q)两像素为中心的邻域(即图像块)，h代表衰减参数。
而对于降噪图像μ对应的全局均方误差M[μ]，则有：
Figure GDA0002905287760000081
式(6)中，E,var和cov分别代表期望，方差和协方差。
M[μ]对衰减参数h的偏导可表达为：
Figure GDA0002905287760000082
式(7)等于0时，M[μ]取得最小值，此时h取得最优值hopt。考虑到
Figure GDA0002905287760000083
因此h取得最优化的条件变为
Figure GDA0002905287760000084
式(1)两边对var[v]求偏导，结合条件
Figure GDA0002905287760000085
可得到如下条件：
Figure GDA0002905287760000086
考虑到在实际降噪处理过程中，上述等式
Figure GDA0002905287760000087
不一定会出现，因此将其修正为：
Figure GDA0002905287760000088
因此，根据上述公式(6)(7)(8)即可确定所述衰减参数h与所述全局均方误差M[μ]的关系。
在确定所述衰减参数h与所述全局均方误差M[μ]的关系后，即可根据相关算法计算出所述衰减参数h的最优值hopt，即最优衰减参数hopt，从而可根据该最优衰减参数hopt获得最优预降噪图像
Figure GDA0002905287760000089
在本实施例中，利用迭代算法计算出所述最优衰减参数hopt。当然，也可以用其他算法计算出所述最优衰减参数hopt。
本实施例中，利用迭代算法求解所述最优衰减参数hopt并获得最优预降噪图像
Figure GDA00029052877600000810
的具体实现步骤如下：
步骤S201、设定cov[η,v0]＝0，var[v0]＝0，E{v0}＝0，n＝1，δh＝10，
Figure GDA0002905287760000091
其中，δh代表每次迭代过程中衰减参数h的变化量，
Figure GDA0002905287760000092
代表噪声η的标准差的估计值。估计值
Figure GDA0002905287760000093
可利用基于伪残差的方法来求得:先人工从噪声图像I中选择一定大小的平滑区φ，该区域内任意像素(x,y),其对应的伪残差R(x,y)为:
Figure GDA0002905287760000094
然后利用平滑区φ中各像素的残差来获得估计值
Figure GDA0002905287760000095
其中,
Figure GDA0002905287760000096
这里|φ|代表平滑区φ的大小，本实施例中，取|φ|＝80×80；所述
Figure GDA0002905287760000097
代表平滑区φ内所有像素对应伪残差的均值。
步骤S202、基于所述估计值
Figure GDA0002905287760000098
生成与噪声图像I相同大小的同步噪声图，该同步噪声图实际上是固定值与随机变量相加的结果，而随机变量的标准差即为
Figure GDA0002905287760000099
本发明实施例中，所述固定值选为140。
步骤S203、通过公式(hn)2＝(hn-1)2+δh来计算hn，然后将hn取代环节1中所述公式(5)里的h,借助公式(4)和(5)获得第n次迭代时的非局部均值降噪图像μn和残差vn(vn＝I-μn)；
步骤S204、计算
Figure GDA00029052877600000910
与
Figure GDA00029052877600000911
这里
Figure GDA00029052877600000912
和
Figure GDA00029052877600000913
分别是vn和vn-1的平均值，该步骤中var[vn]，var[vn-1]，E{vn}和E{vn-1}皆可从噪声图像I中计算获得；
步骤S205、计算
Figure GDA00029052877600000914
这里cov[η,vn]和cov[η,vn-1]可从上述同步噪声图上计算获得；
步骤S206、如果不等式
Figure GDA00029052877600000915
成立，那么就停止迭代过程，并且输出最优衰减参数hopt＝hn和对应的最优预降噪图像
Figure GDA0002905287760000101
否则，令n＝n+1，然后返回步骤3。
所述的处理环节二的具体实现步骤如下：
在获得所述降噪图像μ的像素级均方误差之前，首先计算出所述最优预降噪图像
Figure GDA0002905287760000102
对应的方法噪声
Figure GDA0002905287760000103
然后再通过对该方法噪声进行降噪处理来获得所述方法噪声
Figure GDA0002905287760000104
中的残余图像细节图
Figure GDA0002905287760000105
之后再估算出所述原始图像s的估算值
Figure GDA0002905287760000106
和所述噪声η的估算值
Figure GDA0002905287760000107
最后再通过所述噪声η的估算值和
Figure GDA0002905287760000108
来计算出所述降噪图像μ的像素级均方误差，通过梯度下降法对像素级均方误差进行优化求解，获得所述降噪图像μ的每个像素(x,y)对应的最优衰减参数h′opt(x,y)。
其中，获得所述方法噪声
Figure GDA0002905287760000109
中的残余图像细节图
Figure GDA00029052877600001010
的具体步骤包括：
步骤S301、利用等式
Figure GDA00029052877600001011
计算出最优预降噪图像
Figure GDA00029052877600001012
对应的方法噪声
Figure GDA00029052877600001013
步骤S302、利用非局部均值方法对所述
Figure GDA00029052877600001014
进行降噪处理，得到
Figure GDA00029052877600001015
的降噪结果r(x,y):
Figure GDA00029052877600001016
其中，
Figure GDA00029052877600001017
为当所述ω(x,y,p,q)中的衰减参数h取最优值hopt时的权值；
步骤S303、考虑到r(x,y)仍然包括噪声，因此进一步利用均值滤波方法对所述r(x,y)进行降噪处理，得到所述方法噪声
Figure GDA00029052877600001018
中的残余图像细节图
Figure GDA00029052877600001019
本实施例中，在采用均值滤波方法进行降噪处理时，滤波窗口大小为3×3；当然，实际中可根据具体需要选用其他合适大小的滤波窗口。
本实施例中，在估算所述原始图像s的估算值
Figure GDA00029052877600001020
和所述噪声η的估算值和
Figure GDA0002905287760000111
的过程中，利用公式
Figure GDA0002905287760000112
和
Figure GDA0002905287760000113
来对所述噪声图像的原始图像s和所述噪声图像I的噪声η进行估算，从而获得所述原始图像s的估算值
Figure GDA0002905287760000114
和所述噪声η的估算值
Figure GDA0002905287760000115
在获得所述原始图像s的估算值
Figure GDA0002905287760000116
和所述噪声η的估算值
Figure GDA0002905287760000117
后，即可通过相关算法来获得所述降噪图像μ的每个像素(x,y)对应的最优衰减参数h′opt(x,y)。在本实施例中，在获得所述最优衰减参数h′opt(x,y)的过程中，首先将所述降噪图像μ的每个像素(x,y)对应的均方误差
Figure GDA0002905287760000118
估计为与该像素对应衰减系数h(x,y)相关的代价函数J[h(x,y)]，再利用梯度下降法对J[h(x,y)]进行优化求解，从而获得所述噪声图像I的每个像素(x,y)对应的最优衰减参数h′opt(x,y)。当然，在对J[h(x,y)]进行优化求解时，也可采用其他方法进行优化求解，只要可获得所述最优衰减参数h′opt(x,y)即可。
其中，所述代价函数J[h(x,y)]存在着以下关系：
Figure GDA0002905287760000119
式(10)中，Bias[μ(x,y)])2和Var[μ(x,y)]分别代表偏差和方差项，分别定义为：
Figure GDA00029052877600001110
Figure GDA00029052877600001111
其中，权值
Figure GDA00029052877600001112
定义为：
Figure GDA00029052877600001113
其中，G代表高斯核，*代表卷积运算，P(I(x,y))和P(I(p,q))分别表示以(x,y)和(p,q)两像素为中心的邻域，h(x,y)代表像素(x,y)对应的衰减参数。
本实施例中，利用最速下降法对J[h(x,y)]进行优化求解，以获得所述噪声图像I的每个像素(x,y)对应的最优衰减参数h′opt(x,y)和及其对应的最终降噪图像
Figure GDA0002905287760000126
的具体实现步骤如下：
步骤S501、对噪声图像I中的待考虑的像素(x,y)，初始化x＝1,y＝1,n＝0，
Figure GDA0002905287760000121
步骤S502、在第n次迭代计算中，根据所述公式(10)-公式(13)计算像素(x,y)对应的代价函数J[hn(x,y)]，同时计算J[hn(x,y)]对hn(x,y)的偏导
Figure GDA0002905287760000122
Figure GDA0002905287760000123
步骤S503、使用下面的公式获得hn+1(x,y)：
Figure GDA0002905287760000124
其中α是一个常数，代表步长，对于高斯噪声和斑点噪声，α分别选为8和1。
步骤S504、如果迭代次数达到最大迭代次数N(本专利中，N取为30)，那么输出h′opt(x,y)＝hn+1(x,y)，将h′opt(x,y)取代所述公式(13)中的h(x,y),并用求得的权值
Figure GDA0002905287760000125
取代所述处理环节一中的所述公式(4)里的ω(x,y,p,q)，获得最终降噪图像
Figure GDA0002905287760000131
转入步骤S505，否则n＝n+1，转步骤S502。
步骤S505、判断y是否位于噪声图像I的右边界及x是否位于噪声图像I的下边界，若y位于噪声图像I的右边界且x处于噪声图像I的下边界，则表明噪声图像I中所有像素已处理完毕，因此最速下降法求解过程终止；否则按从左向右、再从上向下的扫描顺序找到噪声图像I的下一个像素，即按以下方式对(x,y)的位置进行更新：若y未处于噪声图像I的右边界，则y＝y+1,转步骤S501，否则y＝1,x＝x+1，转步骤S501。
利用本发明实施例提供的图像降噪方法，可在图像处理中获取更为理想的衰减参数，从而有效提高图像处理过程中的噪声滤除、伪影规避和细节保留等方面的性能。
为了更好地说明本发明实施例提供的图像降噪方法的优越性，下面将传统NLM方法(非局部均值方法)、SURE方法(Stein无偏风险估计方法)、PSURE方法(像素级SURE方法)、Cp-NLM方法(嵌入Cp统计方法的非局部均值方法)、PPB滤波方法(基于概率块的滤波方法)及OBNLM方法作为比较对象。为便于描述，本发明实施例中，将所述的衰减参数自适应非局部均值的图像降噪方法简称为PNLM方法。
为实现公平比较，上述比较方法搜索窗大小统一为17×17，除PNLM方法外的其余方法其图像块大小统一选为9×9，而PNLM方法在计算最优预降噪图像时图像块大小选为9×9，获取最终降噪图像时图像块大小选为17×17。为客观比较各种方法的恢复性能，采用峰值信噪比(peak signal-to-noise ratio,PSNR)和结构相似性(structuralsimilarity,SSIM)作为评价指标。其中，PSNR的单位为dB，PSNR值越大，就代表失真越少；SSIM为衡量两幅图像相似度的指标，其值越大越好，最大为1。
请参阅附图3-5，从图3可看出，对于受方差为40的高斯噪声污染的Boat图像，NLM方法和SURE方法采用了全局固定的衰减系数，因而不能有效地抑制噪声或引起了图像细节的过平滑，PSURE和Cp-NLM方法与上面两种方法相比能更好地保护图像细节信息，但在降噪图像中仍保留了较多的噪声。相比之下，本发明实施例中的PNLM方法能在有效抑制噪声的同时很好地保护图像细节。另外，从图4给出的感兴趣区降噪结果可进一步看到，NLM，SURE，PSURE和Cp-NLM等各方法在感兴趣区产生伪影或引起图像细节的模糊，而本发明实施例中的PNLM方法则能有效消除高斯噪声，能在不引入伪影的情况下保持图像细节的锐度。上述主观比较表明：与NLM，SURE，PSURE和Cp-NLM方法相比，本发明实施例中的PNLM方法在高斯噪声污染图像恢复方面具有更优的性能。
从图5可看出，对于受方差为0.4的斑点噪声污染的Peppers图像，NLM方法和SURE方法在降噪图像中残留了一定的斑点噪声或引起了图像细节的过平滑，PPB方法和OBNLM方法能有效抑制斑点噪声，但在降噪图像局部区域引入了伪影，且一定程度上破坏了图像细节。相比之下，本发明实施例中的PNLM方法不仅能有效斑点噪声，而且能很好地保护图像边缘和纹理等信息。上述主观比较表明：与NLM，SURE，PPB和OBNLM方法相比，PNLM方法在斑点噪声污染图像恢复方面具有更优的性能。
表1、表2分别是各方法对受不同高斯噪声污染的标准图像去噪结果的PSNR/SSIM的对比表格，以及各方法对受不同斑点噪声污染的标准图像去噪结果的PSNR/SSIM的对比表格；其中表1中的ση代表加入的高斯噪声的标准差，而表2中的σu则代表斑点噪声模型
Figure GDA0002905287760000151
中
Figure GDA0002905287760000152
的标准差。从表1和表2可看出：在高斯噪声污染图像降噪方面，本发明实施例提供的PNLM方法可提供较其它方法更高的PSNR和SSIM；在斑点污染噪声图像降噪方面，本发明实施例提供的PNLM方法可提供较其它方法更高的PSNR；就SSIM这一指标而言，仅仅在Barbara图像上，当σu＝0.2和0.3时，本发明实施例中的PNLM方法提供的SSIM略低于OBNLM方法，其余情况下皆比其它方法的SSIM高。因此，上述客观比较表明：整体而言，在高斯噪声和斑点噪声降噪方面，本发明实施例中的PNLM方法较现有衰减参数自适应非局部均值方法具有更优的恢复性能。
表1各方法受不同高斯噪声污染的标准图像去噪结果的PSNR(dB)/SSIM
Figure GDA0002905287760000153
表2各方法受不同斑点噪声污染的标准图像去噪结果的PSNR(dB)/SSIM
Figure GDA0002905287760000161
需要说明的是，上述实施例中提到的内容为本发明较佳的实施方式，具体的实施例仅用以解释本发明，而并非是用于对本发明的限定，在不脱离本发明构思的前提下，任何显而易见的替换均在本发明的保护范围之内。此外，上述各实施例中所涉及到的技术特征，只要彼此之间未构成冲突就可以相互组合。

Claims (5)

1.一种衰减参数自适应非局部均值的图像降噪方法，所述方法步骤包括：
步骤S1、对噪声图像I进行非局部均值降噪处理，获得降噪图像μ，并确定在降噪处理中的衰减参数h与降噪图像μ对应的全局均方误差M[μ]的关系；
所述衰减参数h与所述全局均方误差M[μ]的关系由以下公式确定：
Figure FDA0002726827960000011
Figure FDA0002726827960000012
Figure FDA0002726827960000013
其中：E、var和cov分别代表期望、方差和协方差；s为噪声图像I的原始图像，η为噪声图像I中的噪声，μ为在对噪声图像I进行非局部均值降噪处理后获得的降噪图像，v为相应的残余分量；
步骤S2、根据所述衰减参数h与所述全局均方误差M[μ]之间的关系，计算出所述衰减参数h的最优值hopt，进而根据该最优值hopt确定最优预降噪图像
Figure FDA0002726827960000014
步骤S3、计算出所述最优预降噪图像
Figure FDA0002726827960000015
对应的方法噪声
Figure FDA0002726827960000016
并对该方法噪声
Figure FDA0002726827960000017
进行降噪处理，获得所述方法噪声
Figure FDA0002726827960000018
中的残余图像细节图
Figure FDA0002726827960000019
步骤S4、利用所述最优预降噪图像
Figure FDA0002726827960000021
所述方法噪声
Figure FDA0002726827960000022
以及所述残余图像细节图
Figure FDA0002726827960000023
对所述噪声图像I中的噪声及其对应的原始图像进行估算，获得两者对应的估算值
Figure FDA0002726827960000024
和
Figure FDA0002726827960000025
步骤S5、基于所述
Figure FDA0002726827960000026
和
Figure FDA0002726827960000027
计算出所述降噪图像μ的像素级均方误差，并求解所述像素级均方误差取得最小值时对应的最终降噪图像；
在求解获得所述最终降噪图像过程中，利用最速下降法获得所述降噪图像μ的每个像素(x,y)对应的最优衰减参数h′opt(x,y)，并根据该最优衰减参数h′opt(x,y)获得对应的最终降噪图像；
在获得所述最优衰减参数h′opt(x,y)的过程中，首先将所述降噪图像μ的每个像素(x,y)对应的均方误差
Figure FDA0002726827960000028
估计为与该像素对应衰减系数h(x,y)相关的代价函数J[h(x,y)]，再利用最速下降法对J[h(x,y)]进行优化求解，从而获得所述噪声图像I的每个像素(x,y)对应的最优衰减参数h′opt(x,y)；
其中：
Figure FDA0002726827960000029
μ(x,y)为降噪图像μ中的像素(x,y)的对应灰度值，定义为
Figure FDA00027268279600000210
式(5)中，Ω代表搜索窗口，ω(x,y,p,q)表示(x,y)和(p,q)两像素的相似性，定义为：
Figure FDA00027268279600000211
式(6)中，G代表高斯核，*代表卷积运算，P(I(x,y))和P(I(p,q))分别表示以(x,y)和(p,q)两像素为中心的邻域，h代表衰减参数；
Bias[μ(x,y)])2和Var[μ(x,y)]分别代表偏差和方差项，分别定义为：
Figure FDA0002726827960000031
Figure FDA0002726827960000032
其中，权值
Figure FDA0002726827960000033
定义为：
Figure FDA0002726827960000034
2.根据权利要求1所述的方法，其特征在于：所述步骤S2中，利用迭代算法计算出所述衰减参数h的最优值hopt，并根据该最优值hopt确定最优预降噪图像
Figure FDA0002726827960000035
3.根据权利要求1所述的方法，其特征在于：所述步骤S3中，对所述方法噪声
Figure FDA0002726827960000036
进行降噪处理采用非局部均值方法降噪处理或均值滤波降噪处理。
4.根据权利要求1～3中任意一项所述的方法，其特征在于，所述步骤S4中，利用公式
Figure FDA0002726827960000037
和
Figure FDA0002726827960000038
对所述噪声图像I的原始图像s和所述噪声图像I的噪声η进行估算，从而获得所述原始图像s的估算值
Figure FDA0002726827960000039
和所述噪声η的估算值
Figure FDA00027268279600000310
5.根据权利要求1所述的方法，其特征在于，所述步骤S3中，获得所述方法噪声
Figure FDA00027268279600000311
中的残余图像细节图
Figure FDA00027268279600000312
的具体步骤包括：
步骤S301、利用等式
Figure FDA0002726827960000041
计算出最优预降噪图像
Figure FDA0002726827960000042
对应的方法噪声
Figure FDA0002726827960000043
步骤S302、利用非局部均值方法对所述
Figure FDA0002726827960000044
进行降噪处理，得到
Figure FDA0002726827960000045
的降噪结果r(x,y)：
Figure FDA0002726827960000046
其中，
Figure FDA0002726827960000047
为当所述ω(x,y,p,q)中的衰减参数h取最优值hopt时的权值；
步骤S303、利用均值滤波方法对所述r(x,y)进行降噪处理，得到所述方法噪声
Figure FDA0002726827960000048
中的残余图像细节图
Figure FDA0002726827960000049
CN201810898880.7A 2018-08-08 2018-08-08 一种衰减参数自适应非局部均值的图像降噪方法 Active CN109191391B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810898880.7A CN109191391B (zh) 2018-08-08 2018-08-08 一种衰减参数自适应非局部均值的图像降噪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810898880.7A CN109191391B (zh) 2018-08-08 2018-08-08 一种衰减参数自适应非局部均值的图像降噪方法

Publications (2)

Publication Number Publication Date
CN109191391A CN109191391A (zh) 2019-01-11
CN109191391B true CN109191391B (zh) 2021-03-26

Family

ID=64920751

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810898880.7A Active CN109191391B (zh) 2018-08-08 2018-08-08 一种衰减参数自适应非局部均值的图像降噪方法

Country Status (1)

Country Link
CN (1) CN109191391B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110490869B (zh) * 2019-08-23 2022-03-08 北京机械设备研究所 一种超声图像对比度和横向分辨率优化方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101847257A (zh) * 2010-06-10 2010-09-29 上海电力学院 基于非局部均值与多级定向图像的图像降噪方法
CN102930508A (zh) * 2012-08-30 2013-02-13 西安电子科技大学 基于图像残余信号的非局部均值图像去噪方法
CN103544682A (zh) * 2013-09-17 2014-01-29 华中科技大学 一种三维超声图像非局部均值滤波方法
CN106991661A (zh) * 2017-03-31 2017-07-28 重庆大学 融合kl变换与灰色关联度的非局部均值去噪方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8928813B2 (en) * 2010-10-28 2015-01-06 Microsoft Corporation Methods and apparatus for reducing structured noise in video
US20170109867A1 (en) * 2015-10-16 2017-04-20 Motorola Mobility Llc Camera array for performing non-local means image processing over multiple sequential images

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101847257A (zh) * 2010-06-10 2010-09-29 上海电力学院 基于非局部均值与多级定向图像的图像降噪方法
CN102930508A (zh) * 2012-08-30 2013-02-13 西安电子科技大学 基于图像残余信号的非局部均值图像去噪方法
CN103544682A (zh) * 2013-09-17 2014-01-29 华中科技大学 一种三维超声图像非局部均值滤波方法
CN106991661A (zh) * 2017-03-31 2017-07-28 重庆大学 融合kl变换与灰色关联度的非局部均值去噪方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
An Improved Image Denoising Model Based on Nonlocal Means Filter;Yan Jin等;《Mathematical Problems in Engineering》;20180704;1-13 *
Noise Based Computation of Decay Control Parameter in Nonlocal Means Filter for MRI Restoration;Justin Joseph 等;《Journal of Medical Imaging and Health Informatics》;20160826;第6卷(第4期);1-11 *
Pixel-wise decay parameter adaption for nonlocal means image denoising;Yi Zhan 等;《ournal of Electronic Imaging》;20131219;第22卷(第4期);043034-1至043034-7 *
利用无参考图像内容量度优化非局部均值去噪方法的参数;侯鑫 等;《中国图象图形学报》;20140831;第19卷(第8期);1142-1148 *
混合鲁棒权重和改进方法噪声的两级非局部均值去噪;陆海青、葛洪伟;《计算机工程与科学》;20180731;第40卷(第7期);1227-1236 *

Also Published As

Publication number Publication date
CN109191391A (zh) 2019-01-11

Similar Documents

Publication Publication Date Title
CN108921800B (zh) 基于形状自适应搜索窗口的非局部均值去噪方法
CN109872288B (zh) 用于图像去噪的网络训练方法、装置、终端及存储介质
Thanh et al. An adaptive method for image restoration based on high-order total variation and inverse gradient
US7783125B2 (en) Multi-resolution processing of digital signals
US8509559B2 (en) Hierarchical motion deblurring method for single image
CN113222853B (zh) 一种基于噪声估计的渐进式红外图像降噪方法
JP2007018379A (ja) 画像処理方法及び画像処理装置
JP2004362578A (ja) 重み付きオーバーコンプリートな雑音除去
US10229479B2 (en) Image signal processing apparatus, image signal processing method and image signal processing program
Saadia et al. Fractional order integration and fuzzy logic based filter for denoising of echocardiographic image
CN103971345A (zh) 一种基于改进双边滤波的图像去噪方法
CN113793278A (zh) 一种改进的加权核范数最小化和Laplace算子选择性增强的遥感影像去噪方法
EP3072104B1 (en) Image de-noising method
CN109191391B (zh) 一种衰减参数自适应非局部均值的图像降噪方法
Zhang et al. Diffusion scheme using mean filter and wavelet coefficient magnitude for image denoising
CN111652810A (zh) 一种基于小波域奇异值差分模型的图像去噪方法
CN107085839B (zh) 基于纹理增强与稀疏编码的sar图像降斑方法
CN115829967A (zh) 一种工业金属表面缺陷图像去噪和增强方法
Zhang et al. Hybrid gradient-domain image denoising
Choi et al. False contour reduction using directional dilation and edge-preserving filtering
Kadri et al. Colour Image Denoising using Curvelets and Scale Dependent Shrinkage
Tsai et al. An improved adaptive deconvolution algorithm for single image deblurring
Xu et al. Non-iterative, robust monte carlo noise reduction
Hsia et al. Adaptive Wavelet Shrinkage Based On Intelligent FIS Learned Thresholding
Neole et al. Denoising of Digital Images Using Spatial Domain Edge Detection Approach

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