CN103337053B - 一种基于开关非局部全变分的椒盐噪声污染图像滤波方法 - Google Patents

Abstract

本发明公开了一种基于开关非局部全变分的椒盐噪声污染图像滤波方法，包括以下步骤：通过一个二阶段的形态学检测算子对噪声污染图像进行预处理，得到参考图像和噪声标志位。基于噪声标志位，再用改进的非局部全变分方法对参考图像进行滤波，以得到去噪后的图像。本发明可以有效检测图像的噪声分布情况，产生很低的漏检率和误检率。此外，本发明基于参考图像，能精确计算两个图像块之间的相似度，可在有效抑制椒盐噪声的同时很好地保护图像边缘和纹理等细节信息，其提供的峰值信噪比和结构相似度优于现有椒盐噪声滤波方法。

Description

一种基于开关非局部全变分的椒盐噪声污染图像滤波方法

技术领域

本发明属于图像去噪领域，更具体地，涉及一种椒盐噪声污染图像的开关非局部全变分滤波方法。

背景技术

在图像获取或传输的过程中，由于电磁或机械误差，图像不可避免地会受到脉冲噪声的污染。椒盐噪声是一种常见的脉冲噪声，它会严重影响图像的质量，给图像分割、特征提取等图像后续处理过程带来了极大的困难。研究有效的椒盐噪声消除方法，实现在图像降噪的同时有效保护图像的边缘和纹理等细节信息具有重要理论价值和实际意义。

针对图像中的椒盐噪声，目前主流的降噪技术是基于决策的滤波方法，即先采用噪声检测算法检测图像中的噪声，然后对检测出的噪声污染像素进行滤波，而未被污染的像素则保留其恢复值。典型的基于决策的滤波方法包括开关自适应加权均值(switchingadaptiveweightedmean,简称SAWM)滤波方法、基于拉普拉斯噪声检测的开关中值(Laplaciandetector-basedswitchingmedian,简称LDSM)滤波方法及噪声自适应模糊开关中值(noisyadaptivefuzzyswitchingmedian,简称NAFSM)滤波方法等。这些方法在对某一噪声像素进行处理时，通过选定其邻域内像素点的加权或非加权平均来恢复其灰度值。由于只利用了较小的局部邻域内像素点的统计信息，上述方法往往会导致细节模糊和失真，尤其是在图像中噪声水平较高的时候。和上述方法不同，非局部均值(non-localmean，简称NLM)方法利用两个像素点邻域(也称为图像块)之间的欧式距离来度量像素间的相似度，并利基于相似度的权系数来对图像的所有像素点进行加权平均得到恢复图像。受该思想的启发，研究人员通过引入加权非局部梯度函数，将传统的全变分扩展为非局部全变分(non-localtotalvariation,简称NLTV)方法。虽然NLM方法和NLTV方法可以复原受高斯噪声污染的图像，但无法有效恢复受脉冲噪声污染的图像。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种基于开关非局部全变分的椒盐噪声污染图像滤波方法，其目的在于通过结合基于形态学的噪声检测方法和改进的NLTV方法，达到在有效抑制椒盐噪声的同时有效保护图像的边缘和纹理等细节信息。

为实现上述目的，按照本发明的一个方面，提供了一种基于开关非局部全变分的椒盐噪声污染图像滤波方法，包括以下步骤：

第1步通过一个二阶段的形态学检测算子对椒盐噪声污染的图像进行预处理，以得到参考图像和噪声标志位，具体包括以下子步骤：

(1.1)获得椒盐噪声污染的图像f的形态学腐蚀算子(fΘb)(i,j)与膨胀算子图像f的大小为P×Q，其中：

(fΘb)(i,j)＝min{f(i+s,j+t)-b(s,t)|(i+s,j+t)∈D_f,(s,t)∈D_b}(1)

$(f\⊕b)(i,j)=max\left\{f\right(i-s,j-t)+b(s,t\left)\right|(i-s,j-t)\∈D_f,(s,t)\∈D_b\}---\left(2\right)$

其中(i,j)为图像f中的像素点坐标，b为方形结构元素，大小为L_b×L_b，D_f和D_b分别为图像f与方形结构元素b中元素的取值范围，min表示取集合中最小元素，max表示取集合中最大元素。

(1.2)对形态学腐蚀算子(fΘb)(i,j)与膨胀算子进行n次迭代：

${\left(f\mathrm{\Θ}b\right)}^n\left(i,j\right)=\left\{\begin{array}{cc}\left({\left(f\mathrm{\Θ}b\right)}^{n-1}\mathrm{\Θ}b\right)\left(i,j\right)& 2\≤n\≤N\\ \left(f\mathrm{\Θ}b\right)\left(i,j\right)& n=1\end{array}\right.---\left(3\right)$

${(f\⊕b)}^n(i,j)=\left\{\begin{array}{cc}\left({\left(f\⊕b\right)}^{n-1}\⊕b\right)\left(i,j\right)& 2\≤n\≤N\\ \left(f\⊕b\right)\left(i,j\right)& n=1\end{array}\right.---\left(4\right)$

其中(fΘb)ⁿ(i,j)和分别代表腐蚀算子和膨胀算子的第n次迭代结果，N为最大迭代次数；

(1.3)根据迭代结果并采用以下公式(5)和(6)计算内部和外部形态学梯度：

$g_e^n(i,j)=f(i,j)-{\left(f\mathrm{\Θ}b\right)}^n(i,j)---\left(5\right)$

$g_d^n(i,j)={(f\⊕b)}^n(i,j)-f(i,j)---\left(6\right)$

其中和分别为第n次迭代的内部和外部形态学梯度，并获取混合形态学梯度 $g_h(i,j)=\min \{g_e^N(i,j),g_d^N(i,j)\}.$

(1.4)判断是否有g_h(i,j)＝0，若有，则设置噪声标志位η₁(i,j)＝1，表示位于(i，j)处的像素点为候选噪声点；否则设置噪声标志位η₁(i,j)＝0，表示位于(i，j)处的像素点被判定为未被噪声污染点。

(1.5)判断像素点(i，j)处大小为L_b'×L_b'的邻域内的非噪声点数目S(i,j)是否大于0，若S(i,j)≥1，则采用以下公式(7)和(8)进行形态学腐蚀与膨胀运算，否则从L_b'＝3开始，将方形结构元素b′以L_b'＝L_b'+2的长度向外扩展，直到满足S(i,j)≥1为止，然后采用以下公式(7)和(8)进行形态学腐蚀与膨胀运算：

(fΘb′)^c(i,j)＝min{f(i+s,j+t)-b′(s,t)|η₁(i+s,j+t)＝0,(i+s,j+t)∈D_f,(s,t)∈D_b'}(7)

${\left(f\⊕b^{\′}\right)}^c\left(i,j\right)=\max \left\{f\left(i-s,j-t\right)+b^{\′}\left(s,t\right)|{\mathrm{\η}}_1\left(i-s,j-t\right)=0,\left(i-s,j-t\right)\∈D_f,\left(s,t\right)\∈D_{b^{\′}}\right\}---\left(8\right)$

其中(fΘb′)^c(i,j)表示形态学腐蚀结果，表示形态学膨胀结果；

(1.6)对噪声候选点重新识别，以生成新的噪声标志位η₂(i,j)，具体为：

其中，D(i,j)为m(i,j)与f(i,j)的差值的绝对值，m(i,j)为两种形态学算子(fΘb′)^c(i,j)和的均值，T为检测阈值：

(1.7)计算像素点(i,j)处的像素点修复值v(i,j)，由下式得到：

$v\left(i,j\right)=\frac{{\left({\mathrm{f\Θb}}^{\′}\right)}^c(i,j)+{(f\⊕b^{\′})}^c(i,j)}{2}---\left(10\right)$

(1.8)根据像素点修复值v(i,j)和新的噪声标志位η₂(i,j)并利用下式计算参考图像：

r(i,j)＝η₂(i,j)·v(i,j)+(1-η₂(i,j))·f(i,j)(11)

第2步基于噪声标志位η₂(i,j)，用改进的非局部全变分方法对参考图像进行滤波，以得到去噪后的图像；具体包括以下子步骤：

(2.1)利用改进的非局部全变分方法对参考图像r进行处理，以得到待求的去噪图像u，具体采用以下公式：

$\begin{array}{c}\underset{u}{\inf }F\left(u\right)=\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\left|\right|{\▿}_{NL}u\left|\right|\left(x\right)+\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\mathrm{\λ}\left|u\left(x\right)-r\left(x\right)\right|\\ =\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\sqrt{\underset{y\∈{\mathrm{\Ω}}_x}{\mathrm{\Σ}}{\left(u\left(y\right)-u\left(x\right)\right)}^{2}w\left(x,y\right)}+\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\mathrm{\λ}\left|u\left(x\right)-r\left(x\right)\right|\end{array}---\left(12\right)$

其中表示使F(u)最小时u的取值，Ω_n代表噪声像素点的集合，u(x)是待求的去噪图像，Ω_x表示坐标为x的像素点的邻域，y代表Ω_x内的像素点的坐标； $\left|\right|{\▿}_{NL}u\left|\right|\left(x\right)=\sqrt{{\mathrm{\Σ}}_{y\∈\mathrm{\Ω}}{\left(u\right(y)-u(x\left)\right)}^2\mathrm{\ω}(x,y)}$ 表示非局部梯度的范数。

(2.2)对式(12)利用分裂的布莱格曼方法求最小化问题；

(2.3)对步骤(2.2)的结果中加上L₁惩罚项，以得到无约束问题；

(2.4)用分裂的布莱格曼方法对步骤(2.3)的结果迭代，以得到u的迭代形式；

(2.5)用高斯赛德尔方法对步骤(2.4)的结果进行处理，以得到u的快速解。

步骤(2.1)中，w(x,y)表示权重函数，可由下式定义：

$w(x,y)=e^{-\frac{G*\left|\right|f_N\left(x\right)-f_N\left(y\right)|{|}_2^2}{h^2}}---\left(13\right)$

G是高斯核函数，h是控制指数函数衰减的参数，*代表卷积运算,f_N(x)和f_N(y)分别代表以x和y为中心的窗口内所有像素的灰度值，

步骤(2.2)中是采用以下公式：

$\underset{u,d}{\inf }\underset{x\∈{\mathrm{\Ω}}_n}{\Σ}\left(\right|\left|d\right|\left|\right(x)+\mathrm{\λ}|u\left(x\right)-r\left(x\right)\left|\right)---\left(14\right)$

其中 $d={\▿}_{\mathrm{NL}}u.$

步骤(2.3)中是采用以下公式：

$\underset{u,d}{inf}\underset{x\∈{\mathrm{\Ω}}_n}{\Σ}\left(\right|\left|d\right|\left|\right(x)+\mathrm{\λ}|u\left(x\right)-r\left(x\right)|+\mathrm{\β}|d\left(x\right)-\left({\▿}_{NL}u\right)\left(x\right)\left|\right)---\left(15\right)$

其中β为惩罚系数。

步骤(2.4)中是采用以下公式：

$\left\{\begin{array}{c}\left(u^{k+1},d^{k+1}\right)=\underset{u,d}{\arg \min }\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\left(\left|\right|d\left|\right|\left(x\right)+\mathrm{\λ}|u\left(x\right)-r\left(x\right)|+\mathrm{\β}|\left({\▿}_{NL}u\right)\left(x\right)-b^k\left(x\right)|\right)\\ b^{k+1}=b^k\left(x\right)+d^{k+1}\left(x\right)-\left({\▿}_{NL}{u}^{k+1}\right)\left(x\right)\end{array}\right.---\left(16\right)$

其中，k为迭代次数，u^k+1是第k+1次迭代后得到的降噪图像，d^k+1和b^k+1均为第k+1次迭代的中间结果。

步骤(2.5)中是采用以下公式：

$\begin{array}{c}{u}_x^{k+1,n+1}=\frac{\mathrm{\β}\sqrt{{\left(u_x^{k+1,n}-r_x\right)}^2+\mathrm{\τ}}\left({\mathrm{\Σ}}_{y}w\left(x,y\right)u_y^{k+1,n}-{\mathrm{\Σ}}_y\sqrt{w\left(x,y\right)}\left(d_{xy}^{k+1,n}+b_{yx}^{k+1,n}-d_{yx}^{k+1,n}-b_{xy}^{k+1,n}\right)\right)}{\mathrm{\λ}\sqrt{{\mathrm{\Σ}}_y(d_{xy}^{k+1,n}-\sqrt{w\left(x,y\right)}{\left(u_x^{k+1,n}-b_{xy}^{k+1,n}\right)}^2+\mathrm{\τ}}+\mathrm{\β}\sqrt{{\left(u_x^{k+1,n}-r_x\right)}^2+\mathrm{\τ}}}\\ +\frac{\mathrm{\λ}\sqrt{{\mathrm{\Σ}}_y{\left(d_{xy}^{k+1,n}-\sqrt{w\left(x,y\right)}\left(u_x^{k+1,n}-u_y^{k+1,n}\right)-b_{xy}^{k+1,n}\right)}^2+{\mathrm{\τr}}_x}}{\mathrm{\λ}\sqrt{{\mathrm{\Σ}}_y(d_{xy}^{k+1,n}-\sqrt{w\left(x,y\right)}{\left(u_x^{k+1,n}-b_{xy}^{k+1,n}\right)}^2+\mathrm{\τ}}+\mathrm{\β}\sqrt{{\left(u_x^{k+1,n}-r_x\right)}^2+\mathrm{\τ}}}\end{array}---\left(17\right)$

其中τ为趋近于0的正数，与由以下两式给出：

$d_{xy}^{k+1}=\frac{\sqrt{w(x,y)}(u_y^{k+1}-u_x^{k+1})+b_{xy}^k}{\sqrt{{\mathrm{\Σ}}_{y}w(x,y){(u_y^{k+1}-u_x^{k+1})}^2+{\left(b_{xy}^k\right)}^2}}max(\sqrt{\underset{y}{\Σ}w(x,y){(u_y^{k+1}-u_x^{k+1})}^2+{\left(b_{xy}^k\right)}^2}-\mathrm{\β},0)---\left(18\right)$

$b_{xy}^{k+1}=b_{xy}^k+\sqrt{w(x,y)}(u_y^{k+1}-u_x^{k+1})-d_{xy}^k---\left(19\right)$

其中n和k分别为内迭代次数和外迭代次数。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列的有益效果：

(1)本发明采用两阶段的基于形态学的噪声检测算子，对噪声图像进行预处理，得到参考图像和噪声标志位。本发明对噪声点的检测可产生很低的漏检和误检率。

(2)本发明利用参考图像块来刻画原噪声图像中两像素点的差异，利用噪声标志位，采用改进的非局部全变分方法，只对噪声点修复处理，由此得到滤波后的图像，从而克服了现有开关算法中基于单个像素灰度值的计算方法易受噪声影响的不足，即使在噪声水平比较高的情况下，可更好的修复噪声点，并且具有更好的鲁棒性，能有效的保护图像的边缘及纹理等细节信息。

附图说明

图1是本发明基于开关非局部全变分的椒盐噪声污染图像滤波方法的流程图。

图2是采用本发明的方法和现有技术的方法对大小为512×512的8位灰度Pepper图进行测试的结果比较。

图3是采用本发明的方法和现有技术的方法对大小为512×512的8位灰度Boat图进行测试的结果比较。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1所示，本发明基于开关非局部全变分的椒盐噪声污染图像滤波方法包括以下步骤：

第1步通过一个二阶段的形态学检测算子对椒盐噪声污染的图像进行预处理，以得到参考图像和噪声标志位，具体包括以下子步骤：

(1.1)获得椒盐噪声污染的图像f的形态学腐蚀算子(fΘb)(i,j)与膨胀算子图像f的大小为P×Q，P和Q均为正整数，其中：

(fΘb)(i,j)＝min{f(i+s,j+t)-b(s,t)|(i+s,j+t)∈D_f,(s,t)∈D_b}(1)

$(f\⊕b)(i,j)=max\left\{f\right(i-s,j-t)+b(s,t\left)\right|(i-s,j-t)\∈D_f,(s,t)\∈D_b\}---\left(2\right)$

其中(i,j)为图像f中的像素点坐标，且有1≤i≤P,1≤j≤Q，b为方形结构元素，大小为L_b×L_b，L_b为大于等于3的正整数，D_f和D_b分别为图像f与方形结构元素b中的元素的取值范围。

(1.2)对形态学腐蚀算子(fΘb)(i,j)与膨胀算子进行n次迭代，表示为：

${\left(f\mathrm{\Θ}b\right)}^n\left(i,j\right)=\left\{\begin{array}{cc}\left({\left(f\mathrm{\Θ}b\right)}^{n-1}\mathrm{\Θ}b\right)\left(i,j\right)& 2\≤n\≤N\\ \left(f\mathrm{\Θ}b\right)\left(i,j\right)& n=1\end{array}\right.---\left(3\right)$

${(f\⊕b)}^n(i,j)=\left\{\begin{array}{cc}\left({\left(f\⊕b\right)}^{n-1}\⊕b\right)\left(i,j\right)& 2\≤n\≤N\\ \left(f\⊕b\right)\left(i,j\right)& n=1\end{array}\right.---\left(4\right)$

其中(fΘb)ⁿ(i,j)和分别代表腐蚀算子和膨胀算子的第n次迭代结果，N为最大迭代次数，且为大于等于1的正整数，一般N的取值越大，检测结果越好，但计算时间较长，在本实施方式中，N的取值为5。(1.3)根据迭代结果并采用以下公式(5)和(6)计算内部和外部形态学梯度：

$g_e^n(i,j)=f(i,j)-{\left(f\mathrm{\Θ}b\right)}^n(i,j)---\left(5\right)$

$g_d^n(i,j)={(f\⊕b)}^n(i,j)-f(i,j)---\left(6\right)$

其中和分别为第n次迭代的内部和外部形态学梯度，并获取混合形态学梯度 $g_h(i,j)=\min \{g_e^N(i,j),g_d^N(i,j)\}.$

(1.4)判断是否有g_h(i,j)＝0，若有，则设置噪声标志位η₁(i,j)＝1，表示位于(i，j)处的像素点为候选噪声点；否则设置噪声标志位η₁(i,j)＝0，表示位于(i，j)处的像素点被判定为未被噪声污染点。

(1.5)对于候选噪声点，用方形结构元素b′代替b，该方形结构元素b′的大小为L_b'×L_b'，L_b'为大于等于3的正整数，判断像素点(i,j)处大小为L_b'×L_b'的邻域内的非噪声点数目S(i,j)是否大于0，若S(i,j)≥1，则采用以下公式(7)和(8)进行形态学腐蚀与膨胀运算，否则从L_b'＝3开始，将方形结构元素b′以L_b'＝L_b'+2的长度向外扩展，直到满足S(i,j)≥1为止，然后采用以下公式(7)和(8)进行形态学腐蚀与膨胀运算：

(fΘb′)^c(i,j)＝min{f(i+s,j+t)-b′(s,t)|η₁(i+s,j+t)＝0,(i+s,j+t)∈D_f,(s,t)∈D_b'}(7)

${\left(f\⊕b^{\′}\right)}^c\left(i,j\right)=\max \left(f\left(i-s,j-t\right)+b^{\′}\left(s,t\right)|{\mathrm{\η}}_1\left(i-s,j-t\right)=0,\left(i-s,j-t\right)\∈D_f,\left(s,t\right)\∈D_{b^{\′}}\right)---\left(8\right)$

其中(fΘb′)^c(i,j)表示形态学腐蚀结果，表示形态学膨胀结果；

(1.6)对噪声候选点(即满足η₁(i,j)＝1的图像点)重新识别，以生成新的噪声标志位η₂(i,j)，具体为：

其中，D(i,j)为m(i,j)与f(i,j)的差值的绝对值，m(i,j)为两种形态学算子(fΘb′)^c(i,j)和的均值，T为检测阈值：T＝α·(1-R)，α为常数，其取值为正数，R为噪声比，可由(1.1)中确定的候选噪声点与总像素点相除后得到。

(1.7)计算像素点(i,j)处的像素点修复值v(i,j)，由下式得到：

$v\left(i,j\right)=\frac{{\left({\mathrm{f\Θb}}^{\′}\right)}^c(i,j)+{(f\⊕b^{\′})}^c(i,j)}{2}---\left(10\right)$

其中，(fΘb')^c(i,j)与是由将η₂(i,j)代入式(7)、(8)中替换η₁(i,j)求得。

(1.8)根据像素点修复值v(i,j)和新的噪声标志位η₂(i,j)并利用下式计算参考图像：

r(i,j)＝η₂(i,j)·v(i,j)+(1-η₂(i,j))·f(i,j)(11)

第2步基于噪声标志位η₂(i,j)，用改进的非局部全变分方法对参考图像进行滤波，以得到去噪后的图像；具体包括以下子步骤：

(2.1)利用改进的非局部全变分方法对参考图像r进行处理，以得到待求的去噪图像u，具体采用以下公式：

$\begin{array}{c}\underset{u}{\inf }F\left(u\right)=\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\left|\right|{\▿}_{NL}u\left|\right|\left(x\right)+\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\mathrm{\λ}\left|u\left(x\right)-r\left(x\right)\right|\\ =\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\sqrt{\underset{y\∈{\mathrm{\Ω}}_x}{\mathrm{\Σ}}{\left(u\left(y\right)-u\left(x\right)\right)}^{2}w\left(x,y\right)}+\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\mathrm{\λ}\left|u\left(x\right)-r\left(x\right)\right|\end{array}---\left(12\right)$

其中表示使F(u)最小时u的取值，Ω_n代表噪声像素点的集合，u(x)是待求的去噪图像，为了清楚简要描述起见，x用于表示像素点的坐标(i,j)，r(x)是第1步中得到的参考图像，Ω_x表示坐标为x的像

素点的邻域，y代表Ω_x内的像素点的坐标，其中

$\left|\right|{\▿}_{NL}u\left|\right|\left(x\right)=\sqrt{{\mathrm{\Σ}}_{y\∈\mathrm{\Ω}}{\left(u\right(y)-u(x\left)\right)}^2\mathrm{\ω}(x,y)}$ 表示非局部梯度的范数，的定义为 ${\▿}_{\mathrm{NL}}u=(u\left(y\right)-u\left(x\right))\sqrt{w(x,y)},$ λ是任意正数，本实施方

式中取λ＝0.5。w(x,y)表示权重函数，可由下式定义：

$w(x,y)=e^{-\frac{G*\left|\right|f_N\left(x\right)-f_N\left(y\right)|{|}_2^2}{h^2}}---\left(13\right)$

G是高斯核函数，h是控制指数函数衰减的参数，其为正数，在本实施方式中，h的取值是25，*代表卷积运算，f_N(x)和f_N(y)分别代表以x和y为中心的窗口内所有像素的灰度值，窗口大小为正整数，本实施方式中取3×3。

(2.2)对式(12)利用分裂的布莱格曼(SplitBregman)方法求最小化问题：

$\underset{u,d}{\inf }\underset{x\∈{\mathrm{\Ω}}_n}{\Σ}\left(\right|\left|d\right|\left|\right(x)+\mathrm{\λ}|u\left(x\right)-r\left(x\right)\left|\right)---\left(14\right)$

其中 $d={\▿}_{\mathrm{NL}}u.$

(2.3)对式(14)加上L₁惩罚项，得到无约束问题：

$\underset{u,d}{inf}\underset{x\∈{\mathrm{\Ω}}_n}{\Σ}\left(\right|\left|d\right|\left|\right(x)+\mathrm{\λ}|u\left(x\right)-r\left(x\right)|+\mathrm{\β}|d\left(x\right)-\left({\▿}_{NL}u\right)\left(x\right)\left|\right)---\left(15\right)$

其中β为惩罚系数，其取值范围是80至150，在本实施方式中其为100。

(2.4)用分裂的布莱格曼方法对式(15)迭代，以得到u的迭代形式：

$\left\{\begin{array}{c}\left(u^{k+1},d^{k+1}\right)=\underset{u,d}{\arg \min }\underset{x\∈{\mathrm{\Ω}}_n}{\mathrm{\Σ}}\left(\left|\right|d\left|\right|\left(x\right)+\mathrm{\λ}|u\left(x\right)-r\left(x\right)|+\mathrm{\β}|\left({\▿}_{NL}u\right)\left(x\right)-b^k\left(x\right)|\right)\\ b^{k+1}=b^k\left(x\right)+d^{k+1}\left(x\right)-\left({\▿}_{NL}{u}^{k+1}\right)\left(x\right)\end{array}\right.---\left(16\right)$

其中，k为迭代次数，其取值范围是k≥0，本实施方式中其取值为k＝3；u^k+1是第k+1次迭代后得到的降噪图像，d^k+1和b^k+1均为第k+1次迭代的中间结果。

(2.5)用高斯赛德尔(Gauss-Seidel)方法对式(16)进行处理，以得到u的快速解：

$\begin{array}{c}{u}_x^{k+1,n+1}=\frac{\mathrm{\β}\sqrt{{\left(u_x^{k+1,n}-r_x\right)}^2+\mathrm{\τ}}\left({\mathrm{\Σ}}_{y}w\left(x,y\right)u_y^{k+1,n}-{\mathrm{\Σ}}_y\sqrt{w\left(x,y\right)}\left(d_{xy}^{k+1,n}+b_{yx}^{k+1,n}-d_{yx}^{k+1,n}-b_{xy}^{k+1,n}\right)\right)}{\mathrm{\λ}\sqrt{{\mathrm{\Σ}}_y(d_{xy}^{k+1,n}-\sqrt{w\left(x,y\right)}{\left(u_x^{k+1,n}-b_{xy}^{k+1,n}\right)}^2+\mathrm{\τ}}+\mathrm{\β}\sqrt{{\left(u_x^{k+1,n}-r_x\right)}^2+\mathrm{\τ}}}\\ +\frac{\mathrm{\λ}\sqrt{{\mathrm{\Σ}}_y{\left(d_{xy}^{k+1,n}-\sqrt{w\left(x,y\right)}\left(u_x^{k+1,n}-u_y^{k+1,n}\right)-b_{xy}^{k+1,n}\right)}^2+{\mathrm{\τr}}_x}}{\mathrm{\λ}\sqrt{{\mathrm{\Σ}}_y(d_{xy}^{k+1,n}-\sqrt{w\left(x,y\right)}{\left(u_x^{k+1,n}-b_{xy}^{k+1,n}\right)}^2+\mathrm{\τ}}+\mathrm{\β}\sqrt{{\left(u_x^{k+1,n}-r_x\right)}^2+\mathrm{\τ}}}\end{array}---\left(17\right)$

式(17)中，τ为趋近于0的正数，本实施方式中其取值为τ＝0.0001与由以下两式给出：

$d_{xy}^{k+1}=\frac{\sqrt{w(x,y)}(u_y^{k+1}-u_x^{k+1})+b_{xy}^k}{\sqrt{{\mathrm{\Σ}}_{y}w(x,y){(u_y^{k+1}-u_x^{k+1})}^2+{\left(b_{xy}^k\right)}^2}}max(\sqrt{\underset{y}{\Σ}w(x,y){(u_y^{k+1}-u_x^{k+1})}^2+{\left(b_{xy}^k\right)}^2}-\mathrm{\β},0)---\left(18\right)$

$b_{xy}^{k+1}=b_{xy}^k+\sqrt{w(x,y)}(u_y^{k+1}-u_x^{k+1})-d_{xy}^k---\left(19\right)$

n,k分别为内迭代次数和外迭代次数，此处k与式(16)中的含义相同，n的取值范围为n≥0，本实施方式中取n＝2。实验证明，当n＝2，k＝3时得到的去噪图像u已经足够好。

表1各种方法的漏检数和误检数之和

从上面表1中可以看到本发明对不同噪声图像的噪声检测所产生的漏检数和和误检数之和最少，这为图像的有效去噪提供了基础。对于NAFSM方法，它利用模糊度确定像素属于噪声的概率，因此本专利在噪声检测方面没有与之对比。

如图2和图3所示，分别采用大小为512×512的8位灰度标准Pepper图和Boat图作为测试图。该实例中，L_b＝3，N＝5，α＝5，h＝25，β＝100,椒盐噪声图像噪声比R为10％～80％。为实现公平比较，在SAWM滤波方法中，取探测窗口长度L_d＝1，常数r＝1，阈值T＝2；在LDSM滤波方法中，当R在30％以下时,次数取K＝1，当R为40％～70％时,迭代次数取K＝3，当R在80％以上时,迭代次数取K＝5；在NAFSM滤波方法中，取阈值T1＝10，T2＝30。

图2和图3中，(a)为测试图像,(b)为加噪声比为80％的椒盐噪声图像,(c)为本发明得到的结果,(d)为SAWM方法得到的结果,(e)为LDSM方法得到的结果,(f)为NAFSM方法得到的结果。从图中可以看到SAWN方法、LDSM方法和NAFSM方法在图像的边缘纹理处会产生不同程度的模糊,而本发明较以上算法能更好地保护图像的纹理、边缘等细节信息。

表2上述方法在不同噪声比下Pepper噪声图像上的去噪结果定量比较

上面表2中可以看出，在对Pepper噪声图像去噪结果的比较上，本发明提出的降噪方法对应峰值信噪比(peaksignal-to-noiseratio,简称PSNR)和结构相似度(structuresimilarity,简称SSIM)普遍优于其它方法，特别是在噪声水平较高时，也能取得较高PSNR和SSIM值。

表3上述滤波方法在Boat噪声图像上的去噪结果定量比较

从上面表3中可以看出，在对Boat噪声图像去噪结果的比较上，本发明提出的滤波方法较其它方法具有更高的PSNR和SSIM，即使在噪声水平较高时也能取得较高PSNR和SSIM值。

从图2与图3中也可看出本发明提供的方法即使在噪声水平比较高的情况下(R＝80％)，也能更好地抑制噪声点，更有效地保护图像的边缘及纹理等细节信息。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于开关非局部全变分的椒盐噪声污染图像滤波方法，其特征在于，包括以下步骤：

第1步通过一个二阶段的形态学检测算子对椒盐噪声污染的图像进行预处理，以得到参考图像和噪声标志位，具体包括以下子步骤：

(1.1)获得椒盐噪声污染的图像f的形态学腐蚀算子(fΘb)(i,j)与膨胀算子图像f的大小为P×Q，其中：

(fΘb)(i,j)＝min{f(i+s,j+t)-b(s,t)|(i+s,j+t)∈D_f,(s,t)∈D_b}(1)

其中(i,j)为图像f中的像素点坐标，b为方形结构元素，大小为L_b′L_b，D_f和D_b分别为图像f与方形结构元素b中元素的取值范围，min表示取集合中最小元素，max表示取集合中最大元素；

(1.2)对形态学腐蚀算子(fΘb)(i,j)与膨胀算子进行n次迭代：

其中(fΘb)ⁿ(i,j)和分别代表腐蚀算子和膨胀算子的第n次迭代结果，N为最大迭代次数；

(1.3)根据迭代结果并采用以下公式(5)和(6)计算内部和外部形态学梯度：

其中和分别为第n次迭代的内部和外部形态学梯度，并获取混合形态学梯度

(1.4)判断是否有g_h(i,j)＝0，若有，则设置噪声标志位h₁(i,j)＝1，表示位于(i，j)处的像素点为候选噪声点；否则设置噪声标志位h₁(i,j)＝0，表示位于(i，j)处的像素点被判定为未被噪声污染点；

(1.5)判断像素点(i，j)处大小为L_b'×L_b'的邻域内的非噪声点数目S(i,j)是否大于0，若S(i,j)≥1，则采用以下公式(7)和(8)进行形态学腐蚀与膨胀运算，否则从L_b'＝3开始，将方形结构元素以L_b'＝L_b'+2的长度向外扩展，直到满足S(i,j)≥1为止，然后采用以下公式(7)和(8)进行形态学腐蚀与膨胀运算：

其中表示形态学腐蚀结果，表示形态学膨胀结果；

(1.6)对噪声候选点重新识别，以生成新的噪声标志位h₂(i,j)，具体为：

其中，D(i,j)为m(i,j)与f(i,j)的差值的绝对值，m(i,j)为两种形态学算子和的均值，T为检测阈值：

(1.7)计算像素点(i,j)处的像素点修复值v(i,j)，由下式得到：

(1.8)根据像素点修复值v(i,j)和新的噪声标志位h₂(i,j)并利用下式计算参考图像：

r(i,j)＝η₂(i,j)·v(i,j)+(1-η₂(i,j))·f(i,j)(11)

第2步基于噪声标志位h₂(i,j)，用改进的非局部全变分方法对参考图像进行滤波，以得到去噪后的图像；具体包括以下子步骤：

(2.1)利用改进的非局部全变分方法对参考图像r进行处理，以得到待求的去噪图像u，具体采用以下公式：

其中表示使F(u)最小时u的取值，Ω_n代表噪声像素点的集合，u(x)是待求的去噪图像，Ω_x表示坐标为x的像素点的邻域，y代表Ω_x内的像素点的坐标；表示非局部梯度的范数；

(2.2)对式(12)利用分裂的布莱格曼方法求最小化问题；

(2.3)对步骤(2.2)的结果中加上L₁惩罚项，以得到无约束问题；

(2.4)用分裂的布莱格曼方法对步骤(2.3)的结果迭代，以得到u的迭代形式；

(2.5)用高斯赛德尔方法对步骤(2.4)的结果进行处理，以得到u的快速解。

2.根据权利要求1所述的椒盐噪声污染图像滤波方法，其特征在于，步骤(2.1)中，w(x,y)表示权重函数，可由下式定义：

G是高斯核函数，h是控制指数函数衰减的参数，*代表卷积运算,f_N(x)和f_N(y)分别代表以x和y为中心的窗口内所有像素的灰度值。

3.根据权利要求1所述的椒盐噪声污染图像滤波方法，其特征在于，步骤(2.2)中是采用以下公式：

其中

4.根据权利要求1所述的椒盐噪声污染图像滤波方法，其特征在于，步骤(2.3)中是采用以下公式：

其中β为惩罚系数。

5.根据权利要求1所述的椒盐噪声污染图像滤波方法，其特征在于，步骤(2.4)中是采用以下公式：

其中，k为迭代次数，u^k+1是第k+1次迭代后得到的降噪图像，d^k+1和b^k+1均为第k+1次迭代的中间结果。

6.根据权利要求1所述的椒盐噪声污染图像滤波方法，其特征在于，步骤(2.5)中是采用以下公式：

其中τ为趋近于0的正数，与由以下两式给出：

其中n和k分别为内迭代次数和外迭代次数。