一种两步去除具有图像细节保护能力的非局部双边滤波图像椒盐噪声的方法与流程

本发明涉及基于一种两步去除具有图像细节保护能力的非局部双边滤波图像椒盐噪声的方法，属于计算机图像处理领域。

背景技术：

图像椒盐噪声尤其是高频椒盐噪声给图像的后续处理工作带来极大的挑战。有效地去除椒盐噪声并尽可能多的回复图像的细节信息仍然是一个极具挑战性的工作。

现有的方法是关注于图像的局部信息，利用局部信息估计检测到的噪声点。但是，图像的局部极值点和非平稳点无法精确的恢复图像细节信息，甚至产生错误估计。为了解决这一问题，王等人提出了利用非局部信息去除图像椒盐噪声的方法(具体参见xiaotianwang,shanshanshen,guangmingshi,yuannanxu,peiyuzhang.iterativenon-localmeansfilterforsaltandpeppernoiseremoval.j.vis.commun.imager.38(2016)440–450,2018)，该方法利用非局部区域的特征相似度寻找相似非局部区域块，通过更新权重参数提高图像恢复精度。该方法考虑了利用非局部均值迭代消除椒盐噪声，但是没有空间距离因素对于重建效果的影响。

双边滤波(bilateralfilter)是一种非线性滤波，该方法考虑了图像的空间邻近度与像素值的相似度，被广泛的用于高斯噪声的去除。veerakumaretal.提出了使用自适应双边滤波算法去除高斯噪声(具体参见t.veerakumar,s.esakkirajan,i.vennila,anapproachtominimizeveryhighdensitysaltandpeppernoisethroughtrimmedglobalmean.int.j.compt.appl.39(12):29–332012)，该方法通过修改空间临近高斯函数和像素值相似度高斯函数实现。在像素值相似度函数计算中，像素值之间的差异性通过l2范数度量，该范数可以有效的保护图像纹理细节信息。该方法不能直接适用于椒盐噪声。

并且，现有的上述两种方法以l2范数度量区域之间的相似性，l2范数无法精确地逼近边缘区域，易造成边缘模糊过平滑的现象。

技术实现要素：

针对现有技术的不足，本发明提供了一种两步去除具有图像细节保护能力的非局部双边滤波图像椒盐噪声的方法；

本发明提出了一种可以自适应保护图像边缘和纹理细节的去噪算法，该算法利用nasmf将椒盐噪声转化为类高斯噪声，然后利用自适应非局部高斯滤波算法去除噪声，根据图像非局部细节特征，自适应的趋向于l1范数或者l2范数，从而达到保护边缘和纹理细节的结果。

术语解释：

1、灰度化，在rgb模型中，如果r＝g＝b时，则彩色表示一种灰度颜色，其中r＝g＝b的值叫灰度值，因此，灰度图像每个像素只需一个字节存放灰度值(又称强度值、亮度值)，灰度范围为0-255。一般有分量法最大值法平均值法加权平均法四种方法对彩色图像进行灰度化。

2、图像细节保护能力，由于现有方法问题，在复原图像的过程中，采用临近点或者相似区域对噪声点进行估值，但是估计结果与图像整体风格有所偏差，导致图像失真。这些失真的原因并非所复原图像的结构存在问题，而是由于图像在局部区域内的边缘被过于弱化或者强化导致的。过于弱化会导致过平滑现象的产生，常见于纹理区域；过于强化会导致人工伪迹的产生，常见于强边缘区域。因此，本发明所提出的图像细节保护能力指的是对图像局部的边缘进行修复，避免所复原的图像具有人工伪迹、过平滑等现象的产生。

3、nasmf算法，nasmf算法的全称为非局部自适应中值滤波算法。参见论文：kennykalvintoh,andnorashidimatisa.noiseadaptivefuzzyswitchingmedianfilterforsalt-and-peppernoisereduction.ieeesignalprocessingletters,17(3):281-284march2010。

本发明的技术方案为：

一种两步去除具有图像细节保护能力的非局部双边滤波图像椒盐噪声的方法，包括步骤如下：

s1、遍历图像，利用mask矩阵记录噪声点位置；

s2、利用nasmf算法和mask矩阵对噪声点的定位，对噪声点进行滤波将椒盐噪声转化为高斯噪声；

s3、利用自适应非局部双边滤波方法对图像进行二次滤波，去除噪声。

根据本发明优选的，所述步骤s1，包括：

(1)图像灰度化，得到灰度图像；

(2)对于步骤(1)得到的灰度图像，nsalt＝255,npepper＝0；nsalt表示盐噪声(白点)，npepper表示椒噪声(黑点)，因此，可以得到一个二值蒙版矩阵m，该二值蒙版矩阵m的尺寸与该灰度图像相同，其元素用来标记灰度图像中相应位置的像素点是噪声点还是非噪声点，如式(ⅰ)所示：

式(ⅰ)中，i(i,j)表示灰度图像中像素点(i,j)的灰度值；m(i,j)表示二值蒙版矩阵中(i,j)位置的值，当m(i,j)＝1表示灰度图像中像素点(i,j)是非噪声点；m(i,j)＝0表示灰度图像中像素点(i,j)是噪声点。

根据本发明优选的，所述步骤s2，包括：

(3)计算滤波窗口的中值，并计算噪声点的估计值：

当m(i,j)＝0时，相应的图像像素点(i,j)要进行滤波处理，选取(2s+1)×(2s+1)的滤波窗口，s为预先设定的窗口参数；

滤波窗口形式化的表示为：w2s+1(i,j)＝{i(i+m,j+n)|m,n∈(-s,...,0,,,s)}，w2s+1(i,j)表示以(i,j)位置为中心的、边长为2s+1的窗口；m,n分别表示距离当前窗口中心点(i,j)的行、列距离；i(i+m,j+n)表示图像i中(i+m,j+n)处的像素值；

记：

t2s+1(i,j)表示蒙版矩阵m中以(i,j)位置为中心的、边长为2s+1的窗口内部非噪声点的数量，m(i+m,j+n)表示蒙版矩阵m中(i+m,j+n)处的值；

预先设定区域内非噪声点数量的阈值t＝3，若t2s+1(i,j)>t，说明以(i,j)为中心，(2s+1)×(2s+1)区域内的非噪声点具有足够的数量，则使用式(ⅱ)计算噪声点的所在区域的中值即滤波窗口的中值并将结果保存到imed中，图像imed的尺寸与含噪图像i的尺寸一致：

imed(i,j)＝median{i(i+m,j+n)|m,n∈(-s,...,0...,s)}(ⅱ)

式(ⅱ)中，median{}指的是将图像内以(i,j)为中心，边长为2s+1的窗口内部所有的点从小到大排列，记录中间数的值，并记录到图像imed(i,j)位置上；

若t2s+1(i,j)<t，这说明当前的滤波窗口内没有足够多的非噪声点,则扩展滤波窗口，即s→s+1，重新利用公式(ⅱ)计算滤波窗口的中值，并记录到imed中；

当遍历完矩阵m中所有值为0的点后，利用式(ⅲ)计算噪声点的估计值：

ifinal＝(1-f(i,j))i(i,j)+f(i,j)imed(i,j)(ⅲ)

式(ⅲ)中，ifinal指的是去噪之后的图像；

d(i,j)为滤波窗口内局部灰度变化最大值：

d(i,j)＝max{|i(i+k,j+l)-i(i,j)|,k,l∈(-2,2),(i+k,j+l)≠(i,j)}；

max{}返回的是该集合中的最大值，并记录在d(i,j)中；t1和t2为预先设定的阈值；

(4)将若干副原始不含噪图像ioriginal分别加入椒盐噪声，加入的椒盐噪声的强度为60％～90％，得到对应的若干副含噪图像i，通过步骤(3)得到对应的噪声点的估计值ifinal，对含噪图像i和原始不含噪图像ioriginal做逐点的差值运算，利用式(ⅳ)得到误差图像ierr：

ierr＝ioriginal-ifinal(ⅳ)

对误差图像ierr进行直方图统计，统计结果见附图1，差值结果服从统计意义上均值为0的高斯分布，并对这些直方图进行拟合，该过程由matlab软件自带的fit函数完成。得到一个均值为0的高斯分布，如式(ⅴ)所示：

式(ⅴ)中，表示定义为，即将ngauss定义为gauss(0,σ²)，gauss(0,σ²)表示均值为0，方差为σ²的高斯分布，将式(ⅴ)带入式(ⅳ)中，得到式(ⅵ)：

ioriginal＝ifinal+ngauss(ⅵ)

由此，将原椒盐噪声去噪问题转化为高斯噪声去除问题，即通过图像ifinal获取原始图像ioriginal。

根据本发明优选的，s的初始值为1，即初始时，所选窗口大小为3×3。

t1＝10，t2＝30。

根据本发明优选的，所述步骤s3，包括：

(5)在本发明中ap,q,r,s的计算是根据图像的非局部信息自适应完成的。由于图像的局部特征各不相同，因此，在使用来度量像素间差异性时，使用同一ap,q,r,s值，显然不符合实际。因此，度量像素间差异性，求取不同的ap,q,r,s值，如式(ⅶ)所示：

式(ⅶ)中，y(p,q)和y(r,s)分别表示以(i,j)和(k,l)为中心的图像块中的像素点的灰度值，ε是一个极小数；

下面分析函数的性质。如果y(p,q)与y(r,s)灰度值相近，即(y(p,q)-y(r,s))值较小，则ap,q,r,s较大，此时，自适应范数效果趋近于l2范数，反之，图像区域突变区域，ap,q,r,s较小，使得自适应范数效果趋近于l1范数，

(6)计算空间临近度度量函数ωd(i,j,k,l)及块像素值相似度度量函数ωr(i,j,k,l)，分别如式(ⅷ)、式(ⅸ)所示：

式(ⅷ)、式(ⅸ)中,σd和σr为平滑系数，θ(i,j)表示以(i,j)为中心、大小为3×3的图像块，y(p,q)表示θ(i,j)中的点,θ(k,l)表示以(k,l)为中心、大小为3×3的图像块，y(r,s)表示θ(k,l)中的点，且ap,q,r,s>0，exp()表示指数函数，即exp(·)＝e(·)；

(7)为了进一步对步骤s2的结果进行滤波,提出了如下非局部双边滤波定义，通过非局部双边滤波公式对图像进行二次滤波，去除噪声，得到输出图像ioutput，用来替换噪声点的估计值ifinal，非局部双边滤波公式如式(ⅹ)所示：

式(ⅶ)中,ioutput(i,j)表示在(i,j)位置上的输出估计,ifinal(k,l)表示图像ifinal中(k,l)的像素值,b(i,j)表示ifinal图像中以(i,j)为中心点图像块，图像块的尺寸为奇数，该图像块区域内的点记为(k,l)，即(k,l)∈b(i,j)。ω(i,j,k,l)表示ifinal(k,l)对于估计ioutput(i,j)的贡献度，即权重参数，贡献度ω(i,j,k,l)的计算分为两部分，即ωd(i,j,k,l)和ωr(i,j,k,l)。

根据本发明优选的，ε为10^-3，σd设置为80乘以噪声强度，σr设置为1，图像块的尺寸为7×7。

本发明的有益效果为：

1、本发明提出了一种两步去除椒盐噪声的算法，第一步利用自适应中值滤波算法去除椒盐噪声，第二步利用自适应非局部双边滤波算法进一步抑制噪声，从而达到复原图像的目的。本发明利用自适应非局部双边滤波算法来修复中值滤波算法带来的误差，该算法通过计算图像块之间的相似度和距离来自适应的求解贡献权重。

3、本发明提出了一种在自适应范数中计算尺度算子的方法，通过使用该发放，自适应范数可以根据图像非局部特征选择逼近真实值的范数，在平滑区域，逼近l2范数，在边缘区域，逼近l1范数，从而避免人工伪迹和图像过平滑现象的产生。

附图说明

图1为本发明椒盐噪声转化后的图像与原图像差值图的直方图；

图2为基于lena图像的试验结果示意图；

图3为基于peppers图像的试验结果示意图；

图4为lena图像在加有90％强度的椒盐噪声条件下，不同算法的去噪重建结果示意图；

具体实施方式

下面结合说明书附图和实施例对本发明作进一步限定，但不限于此。

实施例1

一种两步去除具有图像细节保护能力的非局部双边滤波图像椒盐噪声的方法，包括步骤如下：

s1、遍历图像，利用mask矩阵记录噪声点位置；

s2、利用nasmf算法和mask矩阵对噪声点的定位，对噪声点进行滤波将椒盐噪声转化为高斯噪声；

s3、利用自适应非局部双边滤波方法对图像进行二次滤波，去除噪声。

实施例2

根据实施例1所述的一种两步去除具有图像细节保护能力的非局部双边滤波图像椒盐噪声的方法，其区别在于：

步骤s1，包括：

(1)图像灰度化，得到灰度图像；

(2)对于步骤(1)得到的灰度图像，nsalt＝255,npepper＝0；nsalt表示盐噪声(白点)，npepper表示椒噪声(黑点)，因此，我们可以得到一个二值蒙版矩阵m，该二值蒙版矩阵m的尺寸与该灰度图像相同，其元素用来标记灰度图像中相应位置的像素点是噪声点还是非噪声点，如式(ⅰ)所示：

式(ⅰ)中，i(i,j)表示灰度图像中像素点(i,j)的灰度值；m(i,j)表示二值蒙版矩阵中(i,j)位置的值，当m(i,j)＝1表示灰度图像中像素点(i,j)是非噪声点；m(i,j)＝0表示灰度图像中像素点(i,j)是噪声点。

步骤s2，包括：

(3)计算滤波窗口的中值，并计算噪声点的估计值：

当m(i,j)＝0时，相应的图像像素点(i,j)要进行滤波处理，选取(2s+1)×(2s+1)的滤波窗口，s为预先设定的窗口参数；

滤波窗口形式化的表示为：w2s+1(i,j)＝{i(i+m,j+n)|m,n∈(-s,...,0,,,s)}，w2s+1(i,j)表示以(i,j)位置为中心的、边长为2s+1的窗口；m,n分别表示距离当前窗口中心点(i,j)的行、列距离；i(i+m,j+n)表示图像i中(i+m,j+n)处的像素值；

记：

t2s+1(i,j)表示蒙版矩阵m中以(i,j)位置为中心的、边长为2s+1的窗口内部非噪声点的数量，m(i+m,j+n)表示蒙版矩阵m中(i+m,j+n)处的值；

预先设定区域内非噪声点数量的阈值t＝3，若t2s+1(i,j)>t，说明以(i,j)为中心，(2s+1)×(2s+1)区域内的非噪声点具有足够的数量，则使用式(ⅱ)计算噪声点的所在区域的中值即滤波窗口的中值并将结果保存到imed中，图像imed的尺寸与含噪图像i的尺寸一致：

imed(i,j)＝median{i(i+m,j+n)|m,n∈(-s,...,0...,s)}(ⅱ)

式(ⅱ)中，median{}指的是将图像内以(i,j)为中心，边长为2s+1的窗口内部所有的点从小到大排列，记录中间数的值，并记录到图像imed(i,j)位置上；

若t2s+1(i,j)<t，这说明当前的滤波窗口内没有足够多的非噪声点,则扩展滤波窗口，即s→s+1，重新利用公式(ⅱ)计算滤波窗口的中值，并记录到imed中；

当遍历完矩阵m中所有值为0的点后，利用式(ⅲ)计算噪声点的估计值：

ifinal＝(1-f(i,j))i(i,j)+f(i,j)imed(i,j)(ⅲ)

式(ⅲ)中，ifinal指的是去噪之后的图像；

d(i,j)为滤波窗口内局部灰度变化最大值：

d(i,j)＝max{|i(i+k,j+l)-i(i,j)|,k,l∈(-2,2),(i+k,j+l)≠(i,j)}；

max{}返回的是该集合中的最大值，并记录在d(i,j)中；t1和t2为预先设定的阈值；

(4)将若干副原始不含噪图像ioriginal分别加入椒盐噪声，加入的椒盐噪声的强度为60％～90％，得到对应的若干副含噪图像i，通过步骤(3)得到对应的噪声点的估计值ifinal，对含噪图像i和原始不含噪图像ioriginal做逐点的差值运算，利用式(ⅳ)得到误差图像ierr：

ierr＝ioriginal-ifinal(ⅳ)

对误差图像ierr进行直方图统计，统计结果见附图1，ierr图像统计结果，横坐标为误差统计范围，纵坐标为归一化后的统计值(百分比)，图例中“-”前部分表示图片名称，“-”后部分表示噪声强度。差值结果服从统计意义上均值为0的高斯分布，并对这些直方图进行拟合，该过程由matlab软件自带的fit函数完成。得到一个均值为0的高斯分布，如式(ⅴ)所示：

式(ⅴ)中，表示定义为，即将ngauss定义为gauss(0,σ²)，gauss(0,σ²)表示均值为0，方差为σ²的高斯分布，将式(ⅴ)带入式(ⅳ)中，得到式(ⅵ)：

ioriginal＝ifinal+ngauss(ⅵ)

由此，将原椒盐噪声去噪问题转化为高斯噪声去除问题，即通过图像ifinal获取原始图像ioriginal。

s的初始值为1，即初始时，所选窗口大小为3×3。

t1＝10，t2＝30。

步骤s3，包括：

(5)在本发明中ap,q,r,s的计算是根据图像的非局部信息自适应完成的。由于图像的局部特征各不相同，因此，在使用来度量像素间差异性时，使用同一ap,q,r,s值，显然不符合实际。因此，度量像素间差异性，求取不同的ap,q,r,s值，如式(ⅶ)所示：

式(ⅶ)中，y(p,q)和y(r,s)分别表示以(i,j)和(k,l)为中心的图像块中的像素点的灰度值，ε是一个极小数；

下面分析函数的性质。如果y(p,q)与y(r,s)灰度值相近，即(y(p,q)-y(r,s))值较小，则ap,q,r,s较大，此时，自适应范数效果趋近于l2范数，反之，图像区域突变区域，ap,q,r,s较小，使得自适应范数效果趋近于l1范数，

(6)计算空间临近度度量函数ωd(i,j,k,l)及块像素值相似度度量函数ωr(i,j,k,l)，分别如式(ⅷ)、式(ⅸ)所示：

式(ⅷ)、式(ⅸ)中,σd和σr为平滑系数，θ(i,j)表示以(i,j)为中心、大小为3×3的图像块，y(p,q)表示θ(i,j)中的点,θ(k,l)表示以(k,l)为中心、大小为3×3的图像块，y(r,s)表示θ(k,l)中的点，且ap,q,r,s>0，exp()表示指数函数，即exp(·)＝e^(·)；

(7)为了进一步对步骤s2的结果进行滤波,提出了如下非局部双边滤波定义，通过非局部双边滤波公式对图像进行二次滤波，去除噪声，得到输出图像ioutput，用来替换噪声点的估计值ifinal，非局部双边滤波公式如式(ⅹ)所示：

式(ⅶ)中,ioutput(i,j)表示在(i,j)位置上的输出估计,ifinal(k,l)表示图像ifinal中(k,l)的像素值,b(i,j)表示ifinal图像中以(i,j)为中心点图像块，图像块的尺寸为奇数，该图像块区域内的点记为(k,l)，即(k,l)∈b(i,j)。ω(i,j,k,l)表示ifinal(k,l)对于估计ioutput(i,j)的贡献度，即权重参数，贡献度ω(i,j,k,l)的计算分为两部分，即ωd(i,j,k,l)和ωr(i,j,k,l)。

ε为10^-3，σd设置为80乘以噪声强度，σr设置为1，图像块的尺寸为7×7。

本实施例提出的算法在两幅常用图像(lena和peppers)上进行了测试，分别对两幅图像施加10％到90％强度的椒盐噪声，数值越大，噪声强度越强，然后利用本本实施例提出的算法对含噪图像进行去噪处理。并与现有的多种常用算法(如中值滤波及其改进算法(mf，nasmf，acwmf，damf)，基于决策的算法(dba)以及基于模糊集的去噪算法(fsap)等)进行了比对，测试结果通过psnr,ssim和ief等三种评价标准进行客观评价，表1是psnr和ssim两个指标的评价结果；

表1

表1为本实施例提出算法与传统算法的数值结果在lena和peppers图像上不同噪声强度下对比表。”/”前为psnr值，后为ssim值，这两个数值越大，说明重建效果越好。

图2和图3是ief指标的评价结果，图2中，不同的线表示不同的对比算法，横坐标是不同的噪声强度，纵坐标不同强度噪声下，不同算法去噪结果的ief评价值。图3中，不同的线表示不同的对比算法，横坐标是不同的噪声强度，纵坐标不同强度噪声下，不同算法去噪结果的ief评价值。图4为lena图像在加有90％强度的椒盐噪声条件下，不同算法的去噪重建结果示意图；图4中，第一行从左往右依次是lena原图，mf方法,acwmf方法,dba方法和nafsm方法；第二行从左往右依次是nasepf方法,nasnlm方法,damf方法,fsap方法以及本实施例提出的算法。其中，psnr值、ssim值和ief值越高说明去噪效果越好。从评价结果(表1和图2，图3、图4)来看，本实施例提出的方法在客观评价指标上，均好于对比算法。