CN104077752A - 一种基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法 - Google Patents

Abstract

本发明公开了一种基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法。首先假设灰度值为0和255的像素是被噪声污染的像素，得到掩膜图像，计算噪声密度；然后循环确定每个像素的复原值，若是噪声点则将带权估计值赋给目标复原图像，否则直接将当前像素值赋值给目标复原图像，并计算当前像素为窗口权值矩阵，通过高斯极大似然估计计算出当前像素的估计值，每次迭代过程都重新计算图像中未被噪声污染的像素灰度值，并计算其对于噪声图像中对应位置像素灰度值的峰值信噪比，如果该峰值信噪比不再增大，则迭代停止。本发明在有效抑制冲激噪声的同时保存局部细节，使得局部结构有更好的对比度，取得更好的图像复原效果。

Description

一种基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法

技术领域

本发明涉及计算机图像处理领域，具体涉及一种基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法。

背景技术

冲激噪声由于在信号获取过程中传感器损坏、故障或者时间错误所导致。中值滤波器就是一种用来处理冲激噪声的有效方法。抑制冲激噪声通常包含两个步骤——噪声检测和噪声去除。标准的中值滤波器(Standard Median Filter,SMF)在冲激噪声密度很高的情况下并不是很有效，并且当冲激噪声密度大于50％时，经过中值滤波器处理过的图像经常会由于细节丢失和新引入的伪影导致复原图像严重的退化。

至今人们发明了很多方法来克服标准中值滤波器的缺点。自适应的中值滤波器(Adaptive Median Filter,AMF)通过只用中值替代那些被噪声污染的像素同时不改变没有被污染的像素值从而提升中值滤波器的效果。此外，自适应的中值滤波器递归的扩大滤波窗口直到确认到一个未被污染的中值。基于决策的算法(Decision Based Algorithm,DBA)和自适应中值滤波器相似除了当中值被标记为一个被污染的像素时DBA是用一个邻域像素值来替代被污染的像素。改进的基于决策的非对称的中值滤波器(Modified Decision Based Unsymmetric TrimmedMedian Filter,MDBUTMF)能提高DBA的效果，它在中值估算时去除未被污染的像素并且当窗口中所有的像素都被污染的时候将被污染的像素值设成所选窗口的平均值。为了在抑制噪声的同时达到保存边缘的效果，提出了一个二阶段算法，它包含一个基于中值检测的滤波器和一个规格最优化的保存边缘的项。还有人提出多尺度的自适应中值滤波器来抑制椒盐噪声。

在抑制冲激噪声方面，中值方法易产生细节模糊和阶梯伪影从而导致图像退化。这个问题对于高密度噪声情况和结构模糊情况更为严重。原因可能归结于中值方法只是把从一个窗口像素得到的中值信息作为中值滤波并且没有利用到结构连续性信息。

发明内容

针对现有技术的不足，本发明在假设灰度值符合高斯分布的基础上，通过递归的高斯极大似然估计方法(Recursive Gaussian Maximum Likelihood Estimation，RGMLE)，在RGMLE算法的框架下，使用合并的确定性和相似性(certainty&similarity)信息来估计方差，提出了基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法(RGMLE-CS)，是一种新的冲激噪声抑制方法。本发明专注于噪声去除方面，只考虑已知灰度的椒盐噪声，通过递归计算像素的高斯极大似然估计值来复原图像。通过分析未被噪声污染的像素的估计误差可以得到一个控制递归次数的鲁棒算法。

为了实现上述目的，本发明采用以下技术方案：

本发明基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法，包括以下步骤：

步骤1、扫描大小为(m×n)的噪声图像x，检测椒盐噪声，得到掩膜图像M^o；

步骤2、由M^o计算噪声密度ρ，定义M^u为递归中标记未被噪声污染像素的位置的矩阵，并将M^u初始化为M^o；定义PSNR_U为噪声图像中未被噪声污染的像素对应位置处估计的像素值关于未被污染像素的峰值信噪比，初始化为0；定义Γ为第k次迭代的PSNR_U(k)与第k-1次迭代的PSNR_U(k-1)的差值，并初始化为无穷大，其中k＝1,2，…,50；

步骤3、若Γ大于0且ρ<＝70％，或者ρ>＝70％，则依次对噪声图像x中每个像素x_j(j＝1,2,3,…m×n)按照步骤4和步骤5进行处理；

步骤4、若像素x_j对应M^o中相同位置的值为0，则判断该像素x_j为被噪声污染的像素，定义像素x_ij为以x_j为中心的窗口N_j范围内的任意一个像素，计算以像素x_ij为中心的像素块与以像素x_j为中心的像素块之间的距离d_ij，然后由d_ij计算x_ij对应的权值w_ij，用窗口中所有的像素x_ij和所述w_ij通过高斯极大似然估计来计算目标复原图像μ中的像素值μ_j，然后将M^u中相应位置的元素置为1；

步骤5、若像素x_j对应M^o中相同位置的值为1，则判断该像素x_j为未被噪声污染的像素，直接将x_j赋给目标复原图像μ中的像素值μ_j；同理，在以像素x_j为中心的窗口N_j范围内，计算以像素x_ij为中心的像素块与以x_j为中心的像素块之间的距离然后由计算x_ij对应的权值用窗口中所有的像素x_ij和所述的通过高斯极大似然估计来计算图像复原值然后将M^u中相应位置的元素置为1；

步骤6、将x更新为图像μ，求步骤5中所有估计的关于原噪声图像中相应位置的像素值x_j的信噪比PSNR_U；由此计算Γ，若Γ大于0，转到步骤3，若Γ小于0，转到步骤7；

步骤7、迭代停止，得到最终复原图像。

进一步地，步骤4中采用L1范数计算两个图像块之间的距离。

进一步地，步骤6中是通过分析未被噪声污染的像素的估计误差来控制循环次数。

与现有技术相比，本发明具有以下有益效果：

1)、计算图像复原值时合并了基于L1范数的块相似性信息，使得复原结果具有很好的保存边缘的特性；

2)、计算w_ij时只利用了图像中未被噪声污染的点，减少了噪声像素对复原值的影响，使得目标复原值具有更高的精度。

3)、在估计像素值时采用了高斯极大似然估计使得估计值更接近图像原始值，复原结果更好；

4)、通过分析未被噪声污染的像素的估计误差来控制循环次数，这种迭代停止策略具有很好的鲁棒性。

5)、本发明可以在有效抑制冲激噪声的同时保存局部细节，使得局部结构有更好的对比度，取得良好的图像复原效果。

附图说明

图1为本发明实施例中原lena图像；

图2为本发明加入80％椒盐噪声的lena图像；

图3为本发明方法处理的lena复原图像；

图4为现有技术SKR方法处理的lena复原图像；

图5为本发明实施例中原pepper图像；

图6为本发明加入80％椒盐噪声的pepper图像；

图7为本发明方法处理的pepper复原图像；

图8为现有技术SKR方法处理的pepper复原图像。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1是原图像，在图1加入80％噪声得到图2，图2即为待处理的噪声图像，下面采用图2来阐明本发明基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法，包括以下步骤：

步骤1、扫描大小为(m×n)的噪声图像x，检测椒盐噪声，得到掩膜图像M^o；

具体的，扫描大小为(m×n)的噪声图像x，利用灰度图像椒盐噪声特性，将待处理图像2中的灰度为0和255的像素所对应M^o中相同位置处的元素都赋为0，其余都赋为1；

步骤2、由M^o计算噪声密度ρ，定义M^u为递归中标记未被噪声污染像素的位置的矩阵，并将M^u初始化为M^o；定义PSNR_U为噪声图像中未被噪声污染的像素对应位置处估计的像素值关于未被污染像素的峰值信噪比，初始化PSNR_U(0)为0；定义Γ为第k次迭代的PSNR_U(k)与第k-1次迭代的PSNR_U(k-1)的差值，并初始化为无穷大，其中k＝1,2，…,50；

具体的，计算M^o中零元素个数占总元素的百分比，得到噪声密度ρ，并且初始化PSNR_U＝0,Γ＝+∞；

步骤3、若Γ大于0且ρ<＝70％，或者ρ>＝70％，则依次对噪声图像x中每个像素x_j(j＝1,2,3,…m×n)按照步骤4和步骤5进行处理；

步骤4、若像素x_j对应M^o中相同位置的值为0，则判断该像素x_j为被噪声污染的像素，定义像素x_ij为以x_j为中心的窗口N_j范围内的任意一个像素，计算以x_ij为中心的像素块与以像素x_j为中心的像素块之间的距离d_ij，然后由d_ij计算x_ij对应的权值w_ij，用窗口中所有的像素x_ij和所述w_ij通过高斯极大似然估计来计算目标复原图像μ中的像素值μ_j，然后将M^u中相应位置的元素置为1；

具体的，若像素x_j对应M^o中相同位置的值为0，则判断该像素x_j为未被噪声污染的像素，定义像素x_ij为以x_j为中心的窗口N_j范围内的任意一个像素，采用L1范数来计算窗口中以x_ij为中心的像素块与以x_j为中心的像素块之间的距离d_ij：

$d_{\mathrm{ij}}={\left|\right|P_{\mathrm{ij}}^u\&CenterDot;P_{\mathrm{ij}}-P_j^u\&CenterDot;P_j\left|\right|}_1---\left(1\right)$

其中，Ρ_j以x_j为中心的像素块，Ρ_ij是以x_ij为中心的像素块，同时也是Ρ_j的邻域块，如图2所示，和是M^u中与像素块Ρ_j和Ρ_ij相同位置的元素块。然后在以像素x_j为中心的窗口N_j范围内，对窗口内所有像素x_ij采用高斯极大似然估计来计算目标复原图像μ中的像素值μ_j,计算方法如下：

对于目标复原图像μ，在一个窗口N_j范围内，对目标复原图像μ中的像素值μ_j，采用窗口内的观察值x_ij写出如下极大似然估计式：

$L(x_{1j},....,x_{\mathrm{Nj}};{\mathrm{\&mu;}}_j)=\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Pi;}}f(x_{\mathrm{ij}}-{\mathrm{\&mu;}}_j)={\left(\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_{\mathrm{ij}}^2}\right)}^{\left|N_j\right|/2}\exp (-\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Sigma;}}{(x_{\mathrm{ij}}-{\mathrm{\&mu;}}_j)}^2/2{\mathrm{\&sigma;}}_{\mathrm{ij}}^2)---\left(2\right)$

其中，|N_j|表示窗口范围内的像素个数，表示x_ij对应的的方差。

以μ_j为未知变量，将上述公式(2)最小化,可以得到如下公式(3):

${\mathrm{\&mu;}}_j=\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Sigma;}}{x}_{\mathrm{ij}}/{\mathrm{\&sigma;}}_{\mathrm{ij}}^2/\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Sigma;}}1/{\mathrm{\&sigma;}}_{\mathrm{ij}}^2---\left(3\right)$

上述式子中用x_ij来估计，当x_ij在原噪声图像中是被污染的像素且还未更新时x_ij对应的方差应为+∞，当x_ij为经过迭代估计获得一些确定性的像素时应将赋为1，而当x_ij为原噪声图像中未被污染的像素，其相应的应赋为λ，并按照公式(4)来计算：

${\mathrm{\&sigma;}}_{\mathrm{ij}}^2=\left\{\begin{array}{cc}+\&infin;& M_{\mathrm{ij}}^u=0\\ 1& (M_{\mathrm{ij}}^u=1)\mathrm{and}(M_{\mathrm{ij}}^o=0)\\ \mathrm{\&lambda;}& M_{\mathrm{ij}}^o=1\end{array}\right.---\left(4\right)$

和分别表示二元值矩阵M^o和M^u中的元素，λ是一个确定参(0＜λ＜1)，在上述公式(4)引入块相似性信息d_ij有：

${\mathrm{\&sigma;}}_{\mathrm{ij}}^2=\left\{\begin{array}{cc}+\&infin;& M_{\mathrm{ij}}^u=0\\ \exp (d_{\mathrm{ij}}/h)& (M_{\mathrm{ij}}^u=1)\mathrm{and}(M_{\mathrm{ij}}^o=0)\\ \mathrm{\&lambda;}\&CenterDot;\exp (d_{\mathrm{ij}}/h)& M_{\mathrm{ij}}^o=1\end{array}\right.---\left(5\right)$

已知参数h是用来控制距离d_ij衰减的，x_ij对应位置的权值w_ij用的倒数来估计，有：

$w_{\mathrm{ij}}=\left\{\begin{array}{cc}0& M_{\mathrm{ij}}^u=0\\ \exp (-d_{\mathrm{ij}}/h)& (M_{\mathrm{ij}}^u=1)\mathrm{and}(M_{\mathrm{ij}}^o=0)\\ \exp (-d_{\mathrm{ij}}/h)/\mathrm{\&lambda;}& M_{\mathrm{ij}}^o=1\end{array}\right.---\left(6\right)$

w_ij的计算只利用了图像中未被噪声污染的点，减少了噪声像素对复原值的影响，并且采用L1范数来计算，使得复原图像很好的保存的原图像细节，由此估计出目标复原图像μ中的像素值μ_j，

${\mathrm{\&mu;}}_j=\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Sigma;}}{x}_{\mathrm{ij}}/{\mathrm{\&sigma;}}_{\mathrm{ij}}^2/\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Sigma;}}1/{\mathrm{\&sigma;}}_{\mathrm{ij}}^2=\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}{x}_{\mathrm{ij}}/\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}---\left(7\right)$

然后将M^u中相应位置的元素置为1；

步骤5、若x_j对应M^o中相同位置的值为1，则判断该像素x_j为未被噪声污染的像素，直接将x_j赋给μ_j，定义像素x_ij为以x_j为中心的窗口N_j范围内的任意一个像素，计算窗口中以每个像素x_ij为中心的像素块与以x_j为中心的像素块之间的距离然后由此计算x_ij对应的权值同理用高斯极大似然估计计算图像复原值同时将M^u中相应位置的元素置为1；

具体的，若x_j对应M^o中相同位置的值为1，则判断该像素x_j为未被污染的像素，直接将x_j赋给μ_j，定义像素x_ij为以x_j为中心的窗口N_j范围内的任意一个像素，计算窗口中以x_ij为中心的像素块与以x_j为中心的像素块之间的距离采用如下公式进行计算：

其中， 分别为像素块Ρ_ij和像素块Ρ_j对应于M^o中相同位置的像素块，用于排除原噪声图中未被噪声污染的像素；采用上述距离来估算是为了能精确反应迭代复原效果，因此将原噪声图像中未被污染的像素去除，只用那些复原的像素值。

由计算权值如下所示(9):

采用高斯极大似然估计计算图像复原值对图像复原值用窗口内的观察值x_ij写出如下极大似然估计式：

$L(x_{1j},....,x_{\mathrm{Nj}};{\widehat{\mathrm{\&mu;}}}_j)=\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Pi;}}f(x_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_j)={\left(\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_{\mathrm{ij}}^2}\right)}^{\left|N_j\right|/2}\exp (-\underset{\mathrm{ij}\&Element;N_j}{\mathrm{\&Sigma;}}{(x_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_j)}^2/2{\mathrm{\&sigma;}}_{\mathrm{ij}}^2)---\left(10\right)$

以μ_j为未知变量，将上述式子(10)最小化,并用替代方差的倒数可以得到(11):

同时将M^u中相应位置的元素置为1；

步骤6、将x更新为图像μ，求步骤5中所有估算的关于原噪声图像中相应的像素值x_j的信噪比PSNR_U；由此计算Γ，若Γ大于0，转到步骤3，若Γ小于0，转到步骤7；

具体的，将x更新为图像μ，并按照如下公式(12)计算步骤5所有估算的关于原噪声图像中相应位置的像素值x_j的信噪比PSNR_U：

$\mathrm{PSNR}\_U\left(k\right)=10{\log }_{10}\left(\right|M^o|\&CenterDot;{225}^2/\underset{j=1}{\overset{\left|M^o\right|}{\mathrm{\&Sigma;}}}{M}_j^o{({\widehat{\mathrm{\&mu;}}}_j-x_j)}^2)---\left(12\right)$

其中k表示第k次迭代,|M^o|表示M^o中非零元素个数。由此计算Γ，Γ＝PSNR_U(k)-PSNR_U(k-1)；若Γ大于零，转到步骤3，否则,即PSNR_U不再增加时，转到步骤7；

步骤7、迭代停止，得到最终复原图像，如图3所示。

效果评估准准则

实验中选择了两个灰度图像Lena，Pepper，考虑了噪声密度从10％到90％(每次增加10％)的情况。分别用RGMLE-CS,AMF,DBA,MDBUTMF和控制核回归方法(steering kernel regression,SKR)处理不同椒盐噪声密度下的图像并比较其实验结果。对于AMF,DBA,MDBUTMF等方法，滤波窗口大小设置成3×3。对于SKR方法，通过调节全局平滑参数和核大小使得处理效果达到最好。对于RGMLE-CS方法，滤波窗口大小分别设置为5×5,块大小设置成7×7以有效描述图像局部结构，平滑参数h设置为3.5，λ设置成0.2。迭代次数T设为50。实验结果表明，RGMLE-CS可以在抑制椒盐噪声的同时保存图像细节，取得良好的复原效果。

视觉评估

通过观察lena的原图和实验图(图1至图4,分别是原图像，噪声图像，SKR处理结果，RGMLE-CS以及局部选定区域的放大图像)和Pepper的原图和实验图(图5至图8),可以看出虽然SKR方法可以有效的抑制冲激噪声，但是采用了块相似性信息的RGMLE-CS方法比SKR方法要好很多，局部细节方面可以取得更高的对比度。

量化评估

为了量化的验证本发明方法的有效性，我们计算原图和复原图像的峰值信噪比(peak signal-to-noise ratio,PSNR)和结构相似度(structural similarity index,SSIM)，这里PSNR和SSIM分别定义为：

$\mathrm{PSNR}(\mathrm{\&mu;}/I)=10{\log }_{10}\left(\right|M|\&CenterDot;{255}^2/\underset{j=1}{\overset{\left|M\right|}{\mathrm{\&Sigma;}}}{({\mathrm{\&mu;}}_j-I_j)}^2)$

$\mathrm{SSIM}(\mathrm{\&mu;},I)=\frac{\left(2m_{\mathrm{\&mu;}}{m}_I\right)(2{\mathrm{\&sigma;}}_{\mathrm{\&mu;I}}+c_2)}{(m_{\mathrm{\&mu;}}^2+m_I^2+c_1)({\mathrm{\&sigma;}}_{\mathrm{\&mu;}}^2+{\mathrm{\&sigma;}}_I^2+c_2)}$

其中，μ和I表示复原图像和原未添加噪声的图像，μ_j表示图像μ中的一个像素，I_j表示图像I的一个像素，|M|表示图像中总的像素个数，m_μ和m_I分别表示图像μ和I的平均灰度，σ_μ和σ_I分别是图像μ和I的标准差，σ_μI是图像μ和I的相关系数，c₁＝(K₁L)²，c₂＝(K₂L)²，L为像素值的动态范围(对于灰度图像为255),K₁和K₂设为0.01和0.03。

表1是部分量化评估结果，从表中可以看出SKR方法和RGMLE-CS方法复原结果的PSNR值和SSIM值都比噪声图像提高很多，表明这两种方法都可以有效抑制冲激噪声。RGMLE-CS关于Lena处理结果的PSNR值比SKR方法的处理结果略低，但是SSIM作为一种更好的图像质量度量，RGMLE-CS方法的复原结果的SSIM值比SKR方法都要高出不少。由此我们可以看出RGMLE-CS方法可以在有效抑制冲激噪声的同时保存图像局部结构信息。实验结果表明，RGMLE-CS在充分考虑块相似性信息后可以取得良好的图像复原效果。

表1

原图像	噪声图像	SKR处理结果	RGMLE-CS处理结果
				Lena	psnr＝6.43dB,ssim＝0.0091	psnr＝31.59dB,ssim＝0.8793	psnr＝31.54dB,ssim＝0.8872
Pepper	psnr＝6.26dB,ssim＝0.0093	psnr＝31.16dB,ssim＝0.8094	psnr＝32.01dB,ssim＝0.8381

Claims (3)

1.一种基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法，其特征在于，包括以下步骤：

步骤1、扫描大小为(m×n)的噪声图像x，检测椒盐噪声，得到掩膜图像M^o；

步骤2、由M^o计算噪声密度ρ，定义M^u为递归中标记未被噪声污染像素的位置的矩阵，并将M^u初始化为M^o；定义PSNR_U为噪声图像中未被噪声污染的像素对应位置处估计的像素值关于未被污染像素的峰值信噪比，初始化为0；定义Γ为第k次迭代的PSNR_U(k)与第k-1次迭代的PSNR_U(k-1)的差值，并初始化为无穷大，其中k＝1,2，…,50；

步骤3、若Γ大于0且ρ<＝70％，或者ρ>＝70％，则依次对噪声图像x中每个像素x_j(j＝1,2,3,…m×n)按照步骤4和步骤5进行处理；

步骤4、若像素x_j对应M^o中相同位置的值为0，则判断该像素x_j为被噪声污染的像素，定义像素x_ij为以x_j为中心的窗口N_j范围内的任意一个像素，计算以像素x_ij为中心的像素块与以像素x_j为中心的像素块之间的距离d_ij，然后由d_ij计算x_ij对应的权值w_ij，用窗口中所有的像素x_ij和所述w_ij通过高斯极大似然估计来计算目标复原图像μ中的像素值μ_j，然后将M^u中相应位置的元素置为1；

步骤5、若像素x_j对应M^o中相同位置的值为1，则判断该像素x_j为未被噪声污染的像素，直接将x_j赋给目标复原图像μ中的像素值μ_j；同理，在以像素xj为中心的窗口N_j范围内，计算以像素x_ij为中心的像素块与以x_j为中心的像素块之间的距离然后由计算x_ij对应的权值用窗口中所有的像素x_ij和所述的通过高斯极大似然估计来计算图像复原值然后将M^u中相应位置的元素置为1；

步骤6、将x更新为图像μ，求步骤5中所有估计的关于原噪声图像中相应位置的像素值x_j的信噪比PSNR_U；由此计算Γ，若Γ大于0，转到步骤3，若Γ小于0，转到步骤7；

步骤7、迭代停止，得到最终复原图像。

2.如权利要求1所述基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法，其特征在于，步骤4中采用L1范数计算两个图像块之间的距离。

3.如权利要求1所述基于置信相似度的递归高斯极大似然估计的冲激噪声抑制方法，其特征在于，步骤6中是通过分析未被噪声污染的像素的估计误差来控制循环次数。