CN105654445A - 一种基于小波变换边缘检测的手机图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于小波变换边缘检测的手机图像去噪方法，能够在有效平滑噪声的同时保留图像边缘细节，去除彩色噪声，从而得到理想的去噪效果，包括如下步骤：第一步，读入彩色噪声图像，将输入的图像从RGB彩色空间转换到YUV彩色空间；第二步，用高斯滤波器平滑图像，对Y通道图像进行三次低通滤波；第三步，利用小波变换算法检测图像边缘；第四步，采用阈值分割得到边缘的二值化图像；第五步，在图像边缘区域采用各向异性扩散平滑噪声；第六步，在图像平滑区域采用均值滤波去噪；第七步，对平滑区域图像进一步去除彩色噪声；第八步，将图像从YUV彩色空间转换到RGB彩色空间，得到最终的去噪图像。

Description

一种基于小波变换边缘检测的手机图像去噪方法

技术领域

本发明属于图像处理领域，涉及一种基于小波变换边缘检测的手机图像去噪方法。

背景技术

手机拍照由于其便捷性广泛应用于生活实际中，人们对手机图像要求越来越高。在图像处理领域，去除噪声一直是研究热点。手机图像去噪过程中易损害边缘信息，为了提高手机图像的质量，在去除噪声的同时又不能使图像的边缘轮廓变模糊，对图像边缘处和平滑区域采用不同的方法去噪。现有的边缘检测算子主要检测信号的高频部分，受噪声影响严重，并不能检测到图像的最佳边缘。基于小波变换的边缘检测算法通过构造小波系数，分离信号和噪声，用检测模极大值方法可有效检测出图像的边缘。

由于图像边缘往往携带大部分的信息，对图像边缘去噪的同时也要保留图像的细节。现有的去噪算法主要适用于特定类型的噪声，如均值滤波主要用于去除高斯噪声，中值滤波主要用于去除椒盐噪声或称脉冲噪声，基于双变量收缩函数的去噪方法主要用于去除高斯噪声，手机图像所含噪声并不属于某种特定类型，且含有彩色噪声。本发明提供一种基于小波变换边缘检测的手机图像去噪方法，通过设置扩散系数，在各个方向上扩散程度根据梯度模的不同而改变，在图像边缘实行较弱的平滑，受噪声类型影响不大。通过减小图像的色度值进一步平滑彩色噪声，从而使处理后的手机图像视觉效果更好。

发明内容

技术问题：本发明提供一种可以有效去除噪声并保留图像边缘的基于小波变换边缘检测的手机图像去噪方法。

技术方案：本发明的一种基于小波变换边缘检测的手机图像去噪方法，包括如下步骤：

步骤1：读入一帧大小为L₁×L₂×3的噪声彩色图像I₀，其中L₁和L₂为正整数，分别表示噪声彩色图像I₀的总行数和总列数，然后将输入的噪声彩色图像I₀从RGB彩色空间转换到YUV彩色空间：

$\left[\begin{array}{c}Y\\ U\\ V\end{array}\right]=\left[\begin{array}{ccc}0.299& 0.587& 0.114\\ -0.147& -0.289& 0.436\\ 0.615& -0.515& -0.100\end{array}\right]\left[\begin{array}{c}R\\ G\\ B\end{array}\right]$

其中，R表示RGB彩色空间中的红色分量，G表示RGB彩色空间中的绿色分量，B表示RGB彩色空间中的蓝色分量，Y表示YUV彩色空间中的亮度值，U表示YUV彩色空间中的饱和度，V表示YUV彩色空间中的色调值，转换后的YUV噪声图像记为f，大小为L₁×L₂×3；

步骤2：用高斯滤波器平滑YUV噪声图像f，对Y通道进行三次低通滤波，得到平滑图像f₁：

f₁＝((f_Y*h₁)*h₂)*h₃

其中，f_Y表示YUV噪声图像f在Y通道的分量，*表示卷积，h₁、h₂、h₃为低通滤波器，分别为：

h₁＝[0.1250.3750.3750.125],

h₂＝[0.12500.37500.37500.125],

h₃＝[0.1250000.3750000.3750000.125]

步骤3：对平滑图像f₁，用小波变换算法检测图像边缘，具体方法如下：

由平滑图像f₁和离散化后的小波函数μ^x和μ^y进行卷积，得到平滑图像f₁沿水平和垂直方向的梯度M^xf₁和M^yf₁：

M^xf₁＝f₁*μ^x,

M^yf₁＝f₁*μ^y

其中，x和y分别表示水平和垂直方向坐标轴，μ^x和μ^y分别为离散化后水平和垂直方向的小波函数，构建方法如下：

构建两个二维小波函数μ^x(x,y)和μ^y(x,y)：

${\mathrm{\μ}}^x(x,y)=\frac{\∂\mathrm{\θ}(x,y)}{\∂x},$

${\mathrm{\μ}}^y(x,y)=\frac{\∂\mathrm{\θ}(x,y)}{\∂y}$

其中，θ(x,y)是二维高斯平滑函数，μ^x(x,y)和μ^y(x,y)分别表示二维高斯平滑函数θ(x,y)沿x和y方向的一阶偏导数，离散化后水平和垂直方向的小波函数μ^x和μ^y分别为：

${\mathrm{\μ}}^x=\left[\begin{array}{ccccc}0.0583& 0.1306& 0.0000& -0.1306& -0.0583\\ 0.2613& 0.5855& 0.0000& -0.5855& -0.2613\\ 0.4308& 0.9653& 0.0000& -0.9653& -0.4308\\ 0.2613& 0.5855& 0.0000& -0.5855& -0.2613\\ 0.0583& 0.1306& 0.0000& -0.1306& -0.0583\end{array}\right],$

${\mathrm{\μ}}^y=\left[\begin{array}{ccccc}0.0583& 0.2613& 0.4308& 0.2613& 0.0583\\ 0.1306& 0.5855& 0.9653& 0.5855& 0.1306\\ 0.0000& 0.0000& 0.0000& 0.0000& 0.0000\\ -0.1306& -0.5855& -0.9653& -0.5855& -0.1306\\ -0.0583& -0.2613& -0.4308& -0.2613& -0.0583\end{array}\right]$

根据平滑图像f₁沿水平和垂直方向的梯度M^xf₁和M^yf₁构建图像梯度的模Cf₁和幅角Af₁：

${\mathrm{Cf}}_1=\sqrt{{\left|M^x{f}_1\right|}^2+{\left|M^y{f}_1\right|}^2},$

${\mathrm{Af}}_1=arctg\left(\frac{M^y{f}_1}{M^x{f}_1}\right)$

其中，|M^xf₁|表示M^xf₁的绝对值，|M^yf₁|表示M^yf₁的绝对值，表示的反正切函数，

沿幅角Af₁方向寻找图像梯度的模Cf₁的局部极大值点，具体方法为：

$f_2\left(i,j\right)=\left\{\begin{array}{cc}{\mathrm{Cf}}_1\left(i+1,j\right),& -\frac{1}{8}\mathrm{\&pi;}\&le;{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{1}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i+1,j+1\right),& \frac{1}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{3}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i,j+1\right),& \frac{3}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{5}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i-1,j+1\right),& \frac{5}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{7}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i-1,j\right),& \frac{7}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{9}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i-1,j-1\right),& -\frac{7}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;-\frac{5}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i,j-1\right),& -\frac{5}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;-\frac{3}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i+1,j-1\right),& -\frac{3}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;-\frac{1}{8}\mathrm{\&pi;}\end{array}\right.$

其中，i和j分别表示整幅图像的行坐标和列坐标，f₂(i,j)表示边缘图像f₂在点(i,j)处的值，Cf₁(i+1,j)、Cf₁(i+1,j+1)、Cf₁(i,j+1)、Cf₁(i-1,j+1)、Cf₁(i-1,j)、Cf₁(i-1,j-1)、Cf₁(i,j-1)、Cf₁(i+1,j-1)分别表示图像梯度的模Cf₁在点(i+1,j)、(i+1,j+1)、(i,j+1)、(i-1,j+1)、(i-1,j)、(i-1,j-1)、(i,j-1)、(i+1,j-1)处的值，Af₁(i,j)表示幅角Af₁在点(i,j)处的值，并以此找出所有图像梯度的模Cf₁的局部极大值点，构成边缘图像f₂；

步骤4：用阈值分割得到边缘的二值化图像，将YUV噪声图像f分为边缘区域和平滑区域，其具体方法为：

将小波变换后的边缘图像f₂拉伸到0～255区间，获得对比度增强后的图像f₃：

$f_3(i,j)=\frac{f_2(i,j)-f_{2min}}{f_{2\max }-f_{2min}}\&times;255$

其中，f₂(i,j)表示边缘图像f₂在点(i,j)处的值，f₃(i,j)表示对比度增强后的图像f₃在点(i,j)处的值，f_2max和f_2min分别为边缘图像f₂的灰度最大值和最小值，

设置阈值T，比T大的像素点为255，比T小的像素点为0，得到二值边缘图像f₄：

$f_4(i,j)=\left\{\begin{array}{c}255,f_3(i,j)>T\\ 0,f_3(i,j)<T\end{array}\right.$

其中，f₄(i,j)表示二值边缘图像f₄在点(i,j)处的值，阈值T根据图像灰度的不同选取不同的值：

$T=\left\{\begin{array}{cc}16,& 0\&le;f_{ave}\&le;60\\ 14,& 60<f_{ave}\&le;75\\ 12,& 75<f_{ave}\&le;90\\ 10,& 90<f_{ave}\&le;255\end{array}\right.$

其中，f_ave表示YUV噪声图像f所有像素点的灰度均值，由所有二值边缘图像f₄中灰度值为255的像素点所对应的YUV噪声图像f中的区域构成边缘区域含噪图像P，由所有二值边缘图像f₄中灰度值为0的像素点所对应的YUV噪声图像f中的区域构成平滑区域含噪图像Q；

步骤5：采用各向异性扩散方法对边缘区域含噪图像P去噪，具体算法如下：

根据各向异性扩散模型，按照如下方式，构建拉普拉斯算子8方向离散化的各向异性扩散方程：

P^t+1(i1,j1)＝P^t(i1,j1)+λ[r_N·▽_NP+r_S·▽_SP+r_W·▽_WP+r_E·▽_EP+r_NE·▽_NEP+r_SE·▽_SEP+r_NW·▽_NWP+r_SW·▽_SWP]^t(i1,j1)其中，i1和j1分别表示边缘区域图像的行坐标和列坐标，t表示迭代次数，设置为20次，P^t(i1,j1)表示边缘区域含噪图像P在点(i1,j1)处第t次迭代后的各向异性扩散去噪值，迭代初始值为边缘区域含噪图像P在Y通道的分量在点(i1,j1)处的值，即P⁰(i1,j1)＝P_Y(i1,j1)，λ表示控制扩散总体强度的常数，取λ＝1/4，r表示扩散系数，其选取遵循图像梯度越大，扩散系数越小，图像梯度越小，扩散系数越大的原则，故扩散系数r的表达式为：

$r=F(\&dtri;P)=\frac{1}{1+{\left(\frac{|\&dtri;P|}{D}\right)}^2}$

其中，D表示扩散门限，取值为10，▽表示梯度算子，F(▽P)表示关于▽P的函数，▽P表示边缘区域含噪图像P的梯度，在各个方向的分量为：

▽_NP(i1,j1)＝P(i1-1,j1)-P(i1,j1),▽_SP(i1,j1)＝P(i1+1,j1)-P(i1,j1),

▽_EP(i1,j1)＝P(i1,j1+1)-P(i1,j1),▽_WP(i1,j1)＝P(i1,j1-1)-P(i1,j1),

▽_NEP(i1,j1)＝P(i1-1,j1+1)-P(i1,j1),▽_SEP(i1,j1)＝P(i1+1,j1+1)-P(i1,j1),

▽_SWP(i1,j1)＝P(i1+1,j1-1)-P(i1,j1),▽_NWP(i1,j1)＝P(i1-1,j1-1)-P(i1,j1)

其中，N、S、E、W分别表示上、下、右、左方向，NE、SE、NW、SW则分别表示右上、右下、左上、左下方向，P(i1-1,j1)、P(i1+1,j1)、P(i1,j1+1)、P(i1,j1-1)、P(i1-1,j1+1)、P(i1+1,j1+1)、P(i1+1,j1-1)、P(i1-1,j1-1)分别表示边缘区域含噪图像P在点(i1-1,j1、(i1+1,j1)、(i1,j1+1)、(i1,j1-1)、(i1-1,j1+1)、(i1+1,j1+1)、(i1+1,j1-1)、(i1-1,j1-1)处的值，▽_NP(i1,j1)、▽_SP(i1,j1)、▽_EP(i1,j1)、▽_WP(i1,j1)、▽_NEP(i1,j1)、▽_SEP(i1,j1)、▽_SWP(i1,j1)、▽_NWP(i1,j1)分别表示边缘区域含噪图像P在点(i1,j1)处上、下、右、左、右上、右下、左上、左下方向上的梯度，

扩散系数r在每次迭代中作为边缘区域含噪图像的梯度▽P的函数被更新：

$r_N^t(i1,j1)=F\&lsqb;{\&dtri;}_N{P}^t(i1,j1)\&rsqb;,r_S^t(i1,j1)=F\&lsqb;{\&dtri;}_S{P}^t(i1,j1)\&rsqb;,$

$r_E^t\left(i1,j1\right)=F\&lsqb;{\&dtri;}_E{P}^t\left(i1,j1\right)\&rsqb;,r_W^t\left(i1,j1\right)=F\&lsqb;{\&dtri;}_W{P}^t\left(i1,j1\right)\&rsqb;,$

$r_{NE}^t\left(i1,j1\right)=F\&lsqb;{\&dtri;}_{NE}{P}^t\left(i1,j1\right)\&rsqb;,r_{SE}^t\left(i1,j1\right)=F\&lsqb;{\&dtri;}_{SE}{P}^t\left(i1,j1\right)\&rsqb;,$

$r_{NW}^t(i1,j1)=F\&lsqb;{\&dtri;}_{NW}{P}^t(i1,j1)\&rsqb;,r_{SW}^t(i1,j1)=F\&lsqb;{\&dtri;}_{SW}{P}^t(i1,j1)\&rsqb;$

其中， 分别表示扩散系数r在点(i1,j1)处第t次迭代时上、下、右、左、右上、右下、左上、左下方向上的分量，分别对应边缘区域含噪图像P在点(i1,j1)处第t次迭代时上、下、右、左、右上、右下、左上、左下方向上的梯度▽_NP^t(i1,j1)、▽_SP^t(i1,j1)、▽_EP^t(i1,j1)、▽_WP^t(i1,j1)、▽_NEP^t(i1,j1)、▽_SEP^t(i1,j1)、▽_SWP^t(i1,j1)、▽_NWP^t(i1,j1)的函数，

经过一定迭代次数，自适应平滑噪声，并保留梯度强的图像边缘；

步骤6：采用均值滤波对平滑区域含噪图像Q进行去噪，得到平滑区域均值滤波后的图像H：

$H(i2,j2)=\frac{1}{m\&times;n}{\mathrm{\&Sigma;}}_{(i2,j2)\&Element;Z_{i2,j2}}Q(i2,j2)$

其中，i2和j2分别表示平滑区域图像的行坐标和列坐标，Q(i2,j2)表示平滑区域含噪图像Q在点(i2,j2)处的值，H(i2,j2)表示平滑区域均值滤波后的图像H在点(i2,j2)处的值，Z_i2,j2表示中心在点(i2,j2)处、大小为m×n的矩形子图像邻域，取m＝n＝17个像素点；

步骤7：寻找出平滑区域均值滤波后的图像H的彩色噪声点，再对彩色噪声点进行平滑处理，其中，彩色噪声点的寻找方法为：

如果满足以下条件，则点(i2,j2)为彩色噪声点：

$\left\{\begin{array}{c}|1.140\&times;H_V(i2,j2)|<\mathrm{\&delta;}\\ |2.032\&times;H_U(i2,j2)|<\mathrm{\&delta;}\\ |0.394\&times;H_U(i2,j2)+0.581\&times;H_V(i2,j2)|<\mathrm{\&delta;}\end{array}\right.$

其中，δ表示判断彩色噪声的阈值，取δ＝30，H_U(i2,j2)和H_V(i2,j2)分别表示平滑区域均值滤波后的图像H在点(i2,j2)处的U值和V值，即饱和度和色调值，

对彩色噪声进行平滑处理的具体方法为：缩小H_U(i2,j2)和H_V(i2,j2)的值，得到平滑区域均值滤波后的图像H在点(i2,j2)处的缩小后的饱和度H_U2(i2,j2)和色调值H_V2(i2,j2)：

H_U2(i2,j2)＝H_U(i2,j2)/k,

H_V2(i2,j2)＝H_V(i2,j2)/k

其中，k表示阈值，根据平滑区域均值滤波后的图像H在点(i2,j2)处亮度值H_Y(i2,j2)的不同选取不同的值：

$k=\left\{\begin{array}{cc}2.0,& 0\&le;H_Y(i2,j2)\&le;40\\ 1.7,& 40<H_Y(i2,j2)\&le;60\\ 1.4,& 60<H_Y(i2,j2)\&le;80\\ 1.1,& 80<H_Y(i2,j2)\&le;255\end{array}\right.;$

步骤8：将经步骤7处理后的彩色图像从YUV彩色空间转换到RGB彩色空间：

$\left[\begin{array}{c}R\\ G\\ B\end{array}\right]=\left[\begin{array}{ccc}1& 0& 1.140\\ 1& -0.394& 0.581\\ 1& 2.032& 0\end{array}\right]\left[\begin{array}{c}Y\\ U\\ V\end{array}\right]$

并将转换后的图像输出。

有益效果：与现有方法相比，本发明具有如下优点：

1、由于图片受噪声干扰严重，现有的边缘检测算子对噪声较为敏感，易把噪声较强的点误判为边缘，损失细节信息。为了精确定位边缘，本发明利用小波变换检测局域突变的能力，通过构造水平和垂直方向的小波函数，与低通滤波后的图像卷积，取得平滑噪声后的图像梯度；根据图像边缘梯度大，噪声点梯度小，通过检测图像梯度的模极大值检测局域突变，减小噪声点的误判；而弱边缘的模值比较小，为了防止边缘丢失，沿图像梯度的幅角多个方向，在像素点邻域通过判断幅角所在范围检测图像梯度的模极大值，对于弱边缘也能在局部区域根据边缘的方向性被检测出来，且多个方向的划分有利于细化边缘，弱化噪声的影响。针对不同灰度值的图片设置合适的阈值定位边缘，可抑制噪声的干扰，提高了边缘检测的精度，从而有效区分图像边缘处和平滑区域。

2、为了在去噪的同时避免边缘模糊，本发明在图像边缘处采用各向异性扩散去噪，设置自适应扩散系数，在各个方向上扩散程度根据梯度模的不同而改变，噪声点梯度低，扩散系数大，平滑效果越强，而在梯度高的边缘区域，扩散系数小，从而实现在图像边缘实行较弱的平滑，保留图像边缘细节。由于待去噪的图像较大，手机图像去噪实时性要求高，本发明在各向异性扩散基础上增加梯度的方向性，有利于提高单次迭代的扩散效果，从而减少扩散的迭代次数，在保证去噪效果的同时提高了运算效率，便于应用实现。

3、现有的图像去噪方法对于固定类型的噪声去噪效果好，而手机图像受拍摄条件如光线、拍摄对象等因素的影响，图像噪声类型特殊，噪声点呈大片块状，图像亮度值越低噪声点越多，且彩色噪声严重，噪声分布不规律，随机性大，不属于高斯噪声、椒盐噪声等现有的噪声类型。本发明在图像平滑区域采用均值滤波可平滑大片块状噪声，再通过减小图像的饱和度和色调值，且根据亮度值的不同缩减程度不同，进一步平滑彩色噪声，从而使处理后的手机图像去噪效果理想，视觉效果更好。

附图说明

图1为一种基于小波变换边缘检测的手机图像去噪方法流程图。

图2为手机拍摄的彩色噪声图像。

图3为用一种基于小波变换边缘检测的手机图像去噪方法去噪后的图像。

具体实施方式

下面结合说明书附图和实施例，清楚、完整地描述本发明方法的详细过程。

步骤1：读入一帧大小为L₁×L₂×3的噪声彩色图像I₀，其中L₁和L₂为正整数，分别表示噪声彩色图像I₀的总行数和总列数，然后将输入的噪声彩色图像I₀从RGB彩色空间转换到YUV彩色空间：

$\left[\begin{array}{c}Y\\ U\\ V\end{array}\right]=\left[\begin{array}{ccc}0.299& 0.587& 0.114\\ -0.147& -0.289& 0.436\\ 0.615& -0.515& -0.100\end{array}\right]\left[\begin{array}{c}R\\ G\\ B\end{array}\right]$

其中，R表示RGB彩色空间中的红色分量，G表示RGB彩色空间中的绿色分量，B表示RGB彩色空间中的蓝色分量，Y表示YUV彩色空间中的亮度值，U表示YUV彩色空间中的饱和度，V表示YUV彩色空间中的色调值，转换后的YUV噪声图像记为f，大小为L₁×L₂×3；

步骤2：用高斯滤波器平滑YUV噪声图像f，对Y通道进行三次低通滤波，得到平滑图像f₁：

f₁＝((f_Y*h₁)*h₂)*h₃

其中，f_Y表示YUV噪声图像f在Y通道的分量，*表示卷积，h₁、h₂、h₃为低通滤波器，分别为：

h₁＝[0.1250.3750.3750.125],

h₂＝[0.12500.37500.37500.125],

h₃＝[0.1250000.3750000.3750000.125]

步骤3：对平滑图像f₁，用小波变换算法检测图像边缘，具体方法如下：

由平滑图像f₁和离散化后的小波函数μ^x和μ^y进行卷积，得到平滑图像f₁沿水平和垂直方向的梯度M^xf₁和M^yf₁：

M^xf₁＝f₁*μ^x,

M^yf₁＝f₁*μ^y

其中，x和y分别表示水平和垂直方向坐标轴，μ^x和μ^y分别为离散化后水平和垂直方向的小波函数，构建方法如下：

构建两个二维小波函数μ^x(x,y)和μ^y(x,y)：

${\mathrm{\&mu;}}^x(x,y)=\frac{\&part;\mathrm{\&theta;}(x,y)}{\&part;x},$

${\mathrm{\&mu;}}^y(x,y)=\frac{\&part;\mathrm{\&theta;}(x,y)}{\&part;y}$

其中，θ(x,y)是二维高斯平滑函数，μ^x(x,y)和μ^y(x,y)分别表示二维高斯平滑函数θ(x,y)沿x和y方向的一阶偏导数，取x＝{-2,-1,0,1,2}，y＝{-2,-1,0,1,2}，离散化后水平和垂直方向的小波函数μ^x和μ^y分别为：

${\mathrm{\&mu;}}^x=\left[\begin{array}{ccccc}0.0583& 0.1306& 0.0000& -0.1306& -0.0583\\ 0.2613& 0.5855& 0.0000& -0.5855& -0.2613\\ 0.4308& 0.9653& 0.0000& -0.9653& -0.4308\\ 0.2613& 0.5855& 0.0000& -0.5855& -0.2613\\ 0.0583& 0.1306& 0.0000& -0.1306& -0.0583\end{array}\right],$

${\mathrm{\&mu;}}^y=\left[\begin{array}{ccccc}0.0583& 0.2613& 0.4308& 0.2613& 0.0583\\ 0.1306& 0.5855& 0.9653& 0.5855& 0.1306\\ 0.0000& 0.0000& 0.0000& 0.0000& 0.0000\\ -0.1306& -0.5855& -0.9653& -0.5855& -0.1306\\ -0.0583& -0.2613& -0.4308& -0.2613& -0.0583\end{array}\right]$

根据平滑图像f₁沿水平和垂直方向的梯度M^xf₁和M^yf₁构建图像梯度的模Cf₁和幅角Af₁：

${\mathrm{Cf}}_1=\sqrt{{\left|M^x{f}_1\right|}^2+{\left|M^y{f}_1\right|}^2},$

${\mathrm{Af}}_1=arctg\left(\frac{M^y{f}_1}{M^x{f}_1}\right)$

其中，|M^xf₁|表示M^xf₁的绝对值，|M^yf₁|表示M^yf₁的绝对值，表示的反正切函数，

沿幅角Af₁方向寻找图像梯度的模Cf₁的局部极大值点，具体方法为：

$f_2\left(i,j\right)=\left\{\begin{array}{cc}{\mathrm{Cf}}_1\left(i+1,j\right),& -\frac{1}{8}\mathrm{\&pi;}\&le;{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{1}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i+1,j+1\right),& \frac{1}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{3}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i,j+1\right),& \frac{3}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{5}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i-1,j+1\right),& \frac{5}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{7}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i-1,j\right),& \frac{7}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;\frac{9}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i-1,j-1\right),& -\frac{7}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;-\frac{5}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i,j-1\right),& -\frac{5}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;-\frac{3}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1\left(i+1,j-1\right),& -\frac{3}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1\left(i,j\right)\&le;-\frac{1}{8}\mathrm{\&pi;}\end{array}\right.$

其中，i和j分别表示整幅图像的行坐标和列坐标，f₂(i,j)表示边缘图像f₂在点(i,j)处的值，Cf₁(i+1,j)、Cf₁(i+1,j+1)、Cf₁(i,j+1)、Cf₁(i-1,j+1)、Cf₁(i-1,j)、Cf₁(i-1,j-1)、Cf₁(i,j-1)、Cf₁(i+1,j-1)分别表示图像梯度的模Cf₁在点(i+1,j)、(i+1,j+1)、(i,j+1)、(i-1,j+1)、(i-1,j)、(i-1,j-1)、(i,j-1)、(i+1,j-1)处的值，Af₁(i,j)表示幅角Af₁在点(i,j)处的值，并以此找出所有图像梯度的模Cf₁的局部极大值点，构成边缘图像f₂；

步骤4：用阈值分割得到边缘的二值化图像，将YUV噪声图像f分为边缘区域和平滑区域，其具体方法为：

将小波变换后的边缘图像f₂拉伸到0～255区间，获得对比度增强后的图像f₃：

$f_3(i,j)=\frac{f_2(i,j)-f_{2min}}{f_{2\max }-f_{2min}}\&times;255$

其中，f₂(i,j)表示边缘图像f₂在点(i,j)处的值，f₃(i,j)表示对比度增强后的图像f₃在点(i,j)处的值，f_2max和f_2min分别为边缘图像f₂的灰度最大值和最小值，

设置阈值T，比T大的像素点为255，比T小的像素点为0，得到二值边缘图像f₄：

$f_4(i,j)=\left\{\begin{array}{cc}255,& f_3(i,j)>T\\ 0,& f_3(i,j)<T\end{array}\right.$

其中，f₄(i,j)表示二值边缘图像f₄在点(i,j)处的值，阈值T根据图像灰度的不同选取不同的值：

$T=\left\{\begin{array}{cc}16,& 0\&le;f_{ave}\&le;60\\ 14,& 60<f_{ave}\&le;75\\ 12,& 75<f_{ave}\&le;90\\ 10,& 90<f_{ave}\&le;255\end{array}\right.$

其中，f_ave表示YUV噪声图像f所有像素点的灰度均值，由于光照比较弱的图片噪声相对较大，设置比较大的阈值能够滤除大部分噪声，将不同图像按灰度均值的不同分为四类，取不同的阈值T，由所有二值边缘图像f₄中灰度值为255的像素点所对应的YUV噪声图像f中的区域构成边缘区域含噪图像P，由所有二值边缘图像f₄中灰度值为0的像素点所对应的YUV噪声图像f中的区域构成平滑区域含噪图像Q；

步骤5：采用各向异性扩散方法对边缘区域含噪图像P去噪，具体算法如下：

根据各向异性扩散模型，按照如下方式，构建拉普拉斯算子8方向离散化的各向异性扩散方程：

P^t+1(i1,j1)＝P^t(i1,j1)+λ[r_N·▽_NP+r_S·▽_SP+r_W·▽_WP+r_E·▽_EP+r_NE·▽_NEP+r_SE·▽_SEP+r_NW·▽_NWP+r_SW·▽_SWP]^t(i1,j1)其中，i1和j1分别表示边缘区域图像的行坐标和列坐标，t表示迭代次数，设置为20次，P^t(i1,j1)表示边缘区域含噪图像P在点(i1,j1)处第t次迭代后的各向异性扩散去噪值，迭代初始值为边缘区域含噪图像P在Y通道的分量在点(i1,j1)处的值，即P⁰(i1,j1)＝P_Y(i1,j1)，λ表示控制扩散总体强度的常数，取λ＝1/4，r表示扩散系数，其选取遵循图像梯度越大，扩散系数越小，图像梯度越小，扩散系数越大的原则，故扩散系数r的表达式为：

$r=F(\&dtri;P)=\frac{1}{1+{\left(\frac{|\&dtri;P|}{D}\right)}^2}$

其中，D表示扩散门限，取值为10，▽表示梯度算子，F(▽P)表示关于▽P的函数，▽P表示边缘区域含噪图像P的梯度，在各个方向的分量为：

▽_NP(i1,j1)＝P(i1-1,j1)-P(i1,j1),▽_SP(i1,j1)＝P(i1+1,j1)-P(i1,j1),

▽_EP(i1,j1)＝P(i1,j1+1)-P(i1,j1),▽_WP(i1,j1)＝P(i1,j1-1)-P(i1,j1),

▽_NEP(i1,j1)＝P(i1-1,j1+1)-P(i1,j1),▽_SEP(i1,j1)＝P(i1+1,j1+1)-P(i1,j1),

▽_SWP(i1,j1)＝P(i1+1,j1-1)-P(i1,j1),▽_NWP(i1,j1)＝P(i1-1,j1-1)-P(i1,j1)

其中，N、S、E、W分别表示上、下、右、左方向，NE、SE、NW、SW则分别表示右上、右下、左上、左下方向，P(i1-1,j1)、P(i1+1,j1)、P(i1,j1+1)、P(i1,j1-1)、P(i1-1,j1+1)、P(i1+1,j1+1)、P(i1+1,j1-1)、P(i1-1,j1-1)分别表示边缘区域含噪图像P在点(i1-1,j1、(i1+1,j1)、(i1,j1+1)、(i1,j1-1)、(i1-1,j1+1)、(i1+1,j1+1)、(i1+1,j1-1)、(i1-1,j1-1)处的值，▽_NP(i1,j1)、▽_SP(i1,j1)、▽_EP(i1,j1)、▽_WP(i1,j1)、▽_NEP(i1,j1)、▽_SEP(i1,j1)、▽_SWP(i1,j1)、▽_NWP(i1,j1)分别表示边缘区域含噪图像P在点(i1,j1)处上、下、右、左、右上、右下、左上、左下方向上的梯度，

扩散系数r在每次迭代中作为边缘区域含噪图像的梯度▽P的函数被更新：

$r_N^t(i1,j1)=F\&lsqb;{\&dtri;}_N{P}^t(i1,j1)\&rsqb;,r_S^t(i1,j1)=F\&lsqb;{\&dtri;}_S{P}^t(i1,j1)\&rsqb;,$

$r_E^t(i1,j1)=F\&lsqb;{\&dtri;}_E{P}^t(i1,j1)\&rsqb;,r_W^t(i1,j1)=F\&lsqb;{\&dtri;}_W{P}^t(i1,j1)\&rsqb;,$

$r_{NE}^t(i1,j1)=F\&lsqb;{\&dtri;}_{NE}{P}^t(i1,j1)\&rsqb;,r_{SE}^t(i1,j1)=F\&lsqb;{\&dtri;}_{SE}{P}^t(i1,j1)\&rsqb;,$

$r_{NW}^t(i1,j1)=F\&lsqb;{\&dtri;}_{NW}{P}^t(i1,j1)\&rsqb;,r_{SW}^t(i1,j1)=F\&lsqb;{\&dtri;}_{SW}{P}^t(i1,j1)\&rsqb;$

其中， 分别表示扩散系数r在点(i1,j1)处第t次迭代时上、下、右、左、右上、右下、左上、左下方向上的分量，分别对应边缘区域含噪图像P在点(i1,j1)处第t次迭代时上、下、右、左、右上、右下、左上、左下方向上的梯度▽_NP^t(i1,j1)、▽_SP^t(i1,j1)、▽_EP^t(i1,j1)、▽_WP^t(i1,j1)、▽_NEP^t(i1,j1)、▽_SEP^t(i1,j1)、▽_SWP^t(i1,j1)、▽_NWP^t(i1,j1)的函数，

经过一定迭代次数，自适应平滑噪声，并保留梯度强的图像边缘；

步骤6：采用均值滤波对平滑区域含噪图像Q进行去噪，得到平滑区域均值滤波后的图像H：

$H(i2,j2)=\frac{1}{m\&times;n}{\mathrm{\&Sigma;}}_{(i2,j2)\&Element;Z_{i2,j2}}Q(i2,j2)$

其中，i2和j2分别表示平滑区域图像的行坐标和列坐标，Q(i2,j2)表示平滑区域含噪图像Q在点(i2,j2)处的值，H(i2,j2)表示平滑区域均值滤波后的图像H在点(i2,j2)处的值，Z_i2,j2表示中心在点(i2,j2)处、大小为m×n的矩形子图像邻域，取m＝n＝17个像素点；

步骤7：寻找出平滑区域均值滤波后的图像H的彩色噪声点，再对彩色噪声点进行平滑处理，其中，彩色噪声点的寻找方法为：

考虑到RGB空间和YUV空间存在彩色空间转换，则RGB空间下，如果满足以下条件，则点(i2,j2)为彩色噪声点：

$\left\{\begin{array}{c}\left|H_R\right(i2,j2)-H_{ave}(i2,j2\left)\right|<\mathrm{\&delta;}\\ \left|H_G\right(i2,j2)-H_{ave}(i2,j2\left)\right|<\mathrm{\&delta;}\\ \left|H_B\right(i2,j2)-H_{ave}(i2,j2\left)\right|<\mathrm{\&delta;}\end{array}\right.$

YUV空间下，如果满足以下条件，则点(i2,j2)为彩色噪声点：

$\left\{\begin{array}{c}|1.140\&times;H_V(i2,j2\left)\right|<\mathrm{\&delta;}\\ |2.032\&times;H_U(i2,j2\left)\right|<\mathrm{\&delta;}\\ |0.394\&times;H_U(i2,j2)+0.581\&times;H_V(i2,j2\left)\right|<\mathrm{\&delta;}\end{array}\right.$

其中，δ表示判断彩色噪声的阈值，取δ＝30，H_R(i2,j2)、H_G(i2,j2)、H_B(i2,j2)分别表示平滑区域均值滤波后的图像H在点(i2,j2)处的R、G、B值，H_ave(i2,j2)表示H_B(i2,j2)值的均值，即H_ave(i2,j2)＝[H_R(i2,j2)+H_G(i2,j2)+H_B(i2,j2)]/3，H_U(i2,j2)和H_V(i2,j2)分别表示平滑区域均值滤波后的图像H在点(i2,j2)处的U值和V值，即饱和度和色调值，

对彩色噪声进行平滑处理的具体方法为：缩小H_U(i2,j2)和H_V(i2,j2)的值，得到平滑区域均值滤波后的图像H在点(i2,j2)处的缩小后的饱和度H_U2(i2,j2)和色调值H_V2(i2,j2)：

H_U2(i2,j2)＝H_U(i2,j2)/k,

H_V2(i2,j2)＝H_V(i2,j2)/k

其中，k表示阈值，由于亮度暗的区域噪声较大，当平滑区域均值滤波后的图像H在点(i2,j2)处亮度值H_Y(i2,j2)较小时，对应的k稍大，而H_Y(i2,j2)较大时对应的k稍小，根据H_Y(i2,j2)的不同k选取不同的值：

$k=\{\begin{array}{cc}2.0,& 0\&le;H_Y\left(i2,j2\right)\&le;40\\ 1.7,& 40<H_Y\left(i2,j2\right)\&le;60\\ 1.4,& 60<H_Y\left(i2,j2\right)\&le;80\\ 1.1,& 80<H_Y\left(i2,j2\right)\&le;255\end{array};$

步骤8：将经步骤7处理后的彩色图像从YUV彩色空间转换到RGB彩色空间：

$\left[\begin{array}{c}R\\ G\\ B\end{array}\right]=\left[\begin{array}{ccc}1& 0& 1.140\\ 1& -0.394& 0.581\\ 1& 2.032& 0\end{array}\right]\left[\begin{array}{c}Y\\ U\\ V\end{array}\right]$

并将转换后的图像输出。

一种基于小波变换边缘检测的手机图像去噪的应用实验：

在应用实验中，对大量实际拍摄的噪声图像进行了测试。图2中上图是一张大小为4128×3096的手机拍摄的彩色噪声图像，下图是由图2上图红框标注区域放大后的图像，图3中上图是利用本发明提出的一种基于小波变换边缘检测的手机图像去噪方法处理后的图像，下图是由图3上图红框标注区域放大后的图像，可以看出，图2中的图像含有大量噪声，彩色噪声十分明显，尤其在边缘区域受噪声干扰严重，利用本发明的方法可以避免把噪声较强的点误判为边缘，并保留梯度小的弱图像边缘，有效区分图像边缘处和平滑区域；有效去除图像中的噪声，在图像边缘实行较弱的平滑，避免边缘模糊；在图像平滑区域去除较大噪声，使平滑区域的图像更加光滑，并且进一步平滑彩色噪声，从而使处理后的手机图像视觉效果更好。

应理解上述实施例仅用于说明本发明技术方案的具体实施方式，而不用于限制本发明的范围。在阅读了本发明之后，本领域技术人员对本发明的各种等同形式的修改和替换均落于本申请权利要求所限定的保护范围。

Claims (1)

1.一种基于小波变换边缘检测的手机图像去噪方法，其特征在于，按照以下步骤进行：

步骤1：读入一帧大小为L₁×L₂×3的噪声彩色图像I₀，其中L₁和L₂为正整数，分别表示噪声彩色图像I₀的总行数和总列数，然后将输入的噪声彩色图像I₀从RGB彩色空间转换到YUV彩色空间：

$\left[\begin{array}{c}Y\\ U\\ V\end{array}\right]=\left[\begin{array}{ccc}0.299& 0.587& 0.114\\ -0.147& -0.289& 0.436\\ 0.615& -0.515& -0.100\end{array}\right]\left[\begin{array}{c}R\\ G\\ B\end{array}\right]$

其中，R表示RGB彩色空间中的红色分量，G表示RGB彩色空间中的绿色分量，B表示RGB彩色空间中的蓝色分量，Y表示YUV彩色空间中的亮度值，U表示YUV彩色空间中的饱和度，V表示YUV彩色空间中的色调值，转换后的YUV噪声图像记为f，大小为L₁×L₂×3；

步骤2：用高斯滤波器平滑YUV噪声图像f，对Y通道进行三次低通滤波，得到平滑图像f₁：

f₁＝((f_Y*h₁)*h₂)*h₃

其中，f_Y表示YUV噪声图像f在Y通道的分量，*表示卷积，h₁、h₂、h₃为低通滤波器，分别为：

h₁＝[0.1250.3750.3750.125],

h₂＝[0.12500.37500.37500.125],

h₃＝[0.1250000.3750000.3750000.125]

步骤3：对平滑图像f₁，用小波变换算法检测图像边缘，具体方法如下：

由平滑图像f₁和离散化后的小波函数μ^x和μ^y进行卷积，得到平滑图像f₁沿水平和垂直方向的梯度M^xf₁和M^yf₁：

M^xf₁＝f₁*μ^x,

M^yf₁＝f₁*μ^y

其中，x和y分别表示水平和垂直方向坐标轴，μ^x和μ^y分别为离散化后水平和垂直方向的小波函数，构建方法如下：

构建两个二维小波函数μ^x(x,y)和μ^y(x,y)：

${\mathrm{\&mu;}}^x(x,y)=\frac{\&part;\mathrm{\&theta;}(x,y)}{\&part;x},$

${\mathrm{\&mu;}}^y(x,y)=\frac{\&part;\mathrm{\&theta;}(x,y)}{\&part;y}$

其中，θ(x,y)是二维高斯平滑函数，μ^x(x,y)和μ^y(x,y)分别表示二维高斯平滑函数θ(x,y)沿x和y方向的一阶偏导数，离散化后水平和垂直方向的小波函数μ^x和μ^y分别为：

${\mathrm{\&mu;}}^x=\left[\begin{array}{ccccc}0.0583& 0.1306& 0.0000& -0.1306& -0.0583\\ 0.2613& 0.5855& 0.0000& -0.5855& -0.2613\\ 0.4308& 0.9653& 0.0000& -0.9653& -0.4308\\ 0.2613& 0.5855& 0.0000& -0.5855& -0.2613\\ 0.0583& 0.1306& 0.0000& -0.1306& -0.0583\end{array}\right],$

${\mathrm{\&mu;}}^y=\left[\begin{array}{ccccc}0.0583& 0.2613& 0.4308& 0.2613& 0.0583\\ 0.1306& 0.5855& 0.9653& 0.5855& 0.1306\\ 0.0000& 0.0000& 0.0000& 0.0000& 0.0000\\ -0.1306& -0.5855& -0.9653& -0.5855& -0.1306\\ -0.0583& -0.2613& -0.4308& -0.2613& -0.0583\end{array}\right]$

根据平滑图像f₁沿水平和垂直方向的梯度M^xf₁和M^yf₁构建图像梯度的模Cf₁和幅角Af₁：

${\mathrm{Cf}}_1=\sqrt{{\left|M^x{f}_1\right|}^2+{\left|M^y{f}_1\right|}^2},$

${\mathrm{Af}}_1=arctg\left(\frac{M^y{f}_1}{M^x{f}_1}\right)$

其中，|M^xf₁|表示M^xf₁的绝对值，|M^yf₁|表示M^yf₁的绝对值，表示的反正切函数，

沿幅角Af₁方向寻找图像梯度的模Cf₁的局部极大值点，具体方法为：

$f_2(i,j)=\left\{\begin{array}{cc}{\mathrm{Cf}}_1(i+1,j),& -\frac{1}{8}\mathrm{\&pi;}\&le;{\mathrm{Af}}_1(i,j)\&le;\frac{1}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1(i+1,j+1),& \frac{1}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1(i,j)\&le;\frac{3}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1(i,j+1),& \frac{3}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1(i,j)\&le;\frac{5}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1(i-1,j+1),& \frac{5}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1(i,j)\&le;\frac{7}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1(i-1,j),& \frac{7}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1(i,j)\&le;\frac{9}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1(i-1,j-1),& -\frac{7}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1(i,j)\&le;-\frac{5}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1(i,j-1),& -\frac{5}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1(i,j)\&le;-\frac{3}{8}\mathrm{\&pi;}\\ {\mathrm{Cf}}_1(i+1,j-1),& -\frac{3}{8}\mathrm{\&pi;}<{\mathrm{Af}}_1(i,j)\&le;-\frac{1}{8}\mathrm{\&pi;}\end{array}\right.$

其中，i和j分别表示整幅图像的行坐标和列坐标，f₂(i,j)表示边缘图像f₂在点(i,j)处的值，Cf₁(i+1,j)、Cf₁(i+1,j+1)、Cf₁(i,j+1)、Cf₁(i-1,j+1)、Cf₁(i-1,j)、Cf₁(i-1,j-1)、Cf₁(i,j-1)、Cf₁(i+1,j-1)分别表示图像梯度的模Cf₁在点(i+1,j)、(i+1,j+1)、(i,j+1)、(i-1,j+1)、(i-1,j)、(i-1,j-1)、(i,j-1)、(i+1,j-1)处的值，Af₁(i,j)表示幅角Af₁在点(i,j)处的值，并以此找出所有图像梯度的模Cf₁的局部极大值点，构成边缘图像f₂；

步骤4：用阈值分割得到边缘的二值化图像，将YUV噪声图像f分为边缘区域和平滑区域，其具体方法为：

将小波变换后的边缘图像f₂拉伸到0～255区间，获得对比度增强后的图像f₃：

$f_3(i,j)=\frac{f_2(i,j)-f_{2min}}{f_{2\max }-f_{2\min }}\&times;255$

其中，f₂(i,j)表示边缘图像f₂在点(i,j)处的值，f₃(i,j)表示对比度增强后的图像f₃在点(i,j)处的值，f_2max和f_2min分别为边缘图像f₂的灰度最大值和最小值，

设置阈值T，比T大的像素点为255，比T小的像素点为0，得到二值边缘图像f₄：

$f_4(i,j)=\left\{\begin{array}{cc}255,& f_3(i,j)>T\\ 0,& f_3(i,j)<T\end{array}\right.$

其中，f₄(i,j)表示二值边缘图像f₄在点(i,j)处的值，阈值T根据图像灰度的不同选取不同的值：

$T=\left\{\begin{array}{cc}16,& 0\&le;f_{ave}\&le;60\\ 14,& 60<f_{ave}\&le;75\\ 12,& 75<f_{ave}\&le;90\\ 10,& 90<f_{ave}\&le;255\end{array}\right.$

其中，f_ave表示YUV噪声图像f所有像素点的灰度均值，由所有二值边缘图像f₄中灰度值为255的像素点所对应的YUV噪声图像f中的区域构成边缘区域含噪图像P，由所有二值边缘图像f₄中灰度值为0的像素点所对应的YUV噪声图像f中的区域构成平滑区域含噪图像Q；

步骤5：采用各向异性扩散方法对边缘区域含噪图像P去噪，具体算法如下：

根据各向异性扩散模型，按照如下方式，构建拉普拉斯算子8方向离散化的各向异性扩散方程：

P^t+1(i1,j1)＝P^t(i1,j1)+λ[r_N·▽_NP+r_S·▽_SP+r_W·▽_WP+r_E·▽_EP+r_NE·▽_NEP+r_SE·▽_SEP+r_NW·▽_NWP+r_SW·▽_SWP]^t(i1,j1)

其中，i1和j1分别表示边缘区域图像的行坐标和列坐标，t表示迭代次数，设置为20次，P^t(i1,j1)表示边缘区域含噪图像P在点(i1,j1)处第t次迭代后的各向异性扩散去噪值，迭代初始值为边缘区域含噪图像P在Y通道的分量在点(i1,j1)处的值，即P⁰(i1,j1)＝P_Y(i1,j1)，λ表示控制扩散总体强度的常数，取λ＝1/4，r表示扩散系数，其选取遵循图像梯度越大，扩散系数越小，图像梯度越小，扩散系数越大的原则，故扩散系数r的表达式为：

$r=F(\&dtri;P)=\frac{1}{1+{\left(\frac{\left|\&dtri;P\right|}{D}\right)}^2}$

其中，D表示扩散门限，取值为10，▽表示梯度算子，F(▽P)表示关于▽P的函数，▽P表示边缘区域含噪图像P的梯度，在各个方向的分量为：

▽_NP(i1,j1)＝P(i1-1,j1)-P(i1,j1),▽_SP(i1,j1)＝P(i1+1,j1)-P(i1,j1),

▽_EP(i1,j1)＝P(i1,j1+1)-P(i1,j1),▽_WP(i1,j1)＝P(i1,j1-1)-P(i1,j1),

▽_NEP(i1,j1)＝P(i1-1,j1+1)-P(i1,j1),▽_SEP(i1,j1)＝P(i1+1,j1+1)-P(i1,j1),

▽_SWP(i1,j1)＝P(i1+1,j1-1)-P(i1,j1),▽_NWP(i1,j1)＝P(i1-1,j1-1)-P(i1,j1)

其中，N、S、E、W分别表示上、下、右、左方向，NE、SE、NW、SW则分别表示右上、右下、左上、左下方向，P(i1-1,j1)、P(i1+1,j1)、P(i1,j1+1)、P(i1,j1-1)、P(i1-1,j1+1)、P(i1+1,j1+1)、P(i1+1,j1-1)、P(i1-1,j1-1)分别表示边缘区域含噪图像P在点(i1-1,j1)、(i1+1,j1)、(i1,j1+1)、(i1,j1-1)、(i1-1,j1+1)、(i1+1,j1+1)、(i1+1,j1-1)、(i1-1,j1-1)处的值，▽_NP(i1,j1)、▽_SP(i1,j1)、▽_EP(i1,j1)、▽_WP(i1,j1)、▽_NEP(i1,j1)、▽_SEP(i1,j1)、▽_SWP(i1,j1)、▽_NWP(i1,j1)分别表示边缘区域含噪图像P在点(i1,j1)处上、下、右、左、右上、右下、左上、左下方向上的梯度，

扩散系数r在每次迭代中作为边缘区域含噪图像的梯度▽P的函数被更新：

$r_N^t(i1,j1)=F\&lsqb;{\&dtri;}_N{P}^t(i1,j1)\&rsqb;,r_S^t(i1,j1)=F\&lsqb;{\&dtri;}_S{P}^t(i1,j1)\&rsqb;,$

$r_E^t(i1,j1)=F\&lsqb;{\&dtri;}_E{P}^t(i1,j1)\&rsqb;,r_W^t(i1,j1)=F\&lsqb;{\&dtri;}_W{P}^t(i1,j1)\&rsqb;,$

$r_{NE}^t(i1,j1)=F\&lsqb;{\&dtri;}_{NE}{P}^t(i1,j1)\&rsqb;,r_{SE}^t(i1,j1)=F\&lsqb;{\&dtri;}_{SE}{P}^t(i1,j1)\&rsqb;,$

$r_{NW}^t(i1,j1)=F\&lsqb;{\&dtri;}_{NW}{P}^t(i1,j1)\&rsqb;,r_{SW}^t(i1,j1)=F\&lsqb;{\&dtri;}_{SW}{P}^t(i1,j1)\&rsqb;$

其中， 分别表示扩散系数r在点(i1,j1)处第t次迭代时上、下、右、左、右上、右下、左上、左下方向上的分量，分别对应边缘区域含噪图像P在点(i1,j1)处第t次迭代时上、下、右、左、右上、右下、左上、左下方向上的梯度▽_NP^t(i1,j1)、▽_SP^t(i1,j1)、▽_EP^t(i1,j1)、▽_WP^t(i1,j1)、▽_NEP^t(i1,j1)、▽_SEP^t(i1,j1)、▽_SWP^t(i1,j1)、▽_NWP^t(i1,j1)的函数，

经过一定迭代次数，自适应平滑噪声，并保留梯度强的图像边缘；

步骤6：采用均值滤波对平滑区域含噪图像Q进行去噪，得到平滑区域均值滤波后的图像H：

$H(i2,j2)=\frac{1}{m\&times;n}{\mathrm{\&Sigma;}}_{(i2,j2)\&Element;Z_{i2,j2}}Q(i2,j2)$

其中，i2和j2分别表示平滑区域图像的行坐标和列坐标，Q(i2,j2)表示平滑区域含噪图像Q在点(i2,j2)处的值，H(i2,j2)表示平滑区域均值滤波后的图像H在点(i2,j2)处的值，Z_i2,j2表示中心在点(i2,j2)处、大小为m×n的矩形子图像邻域，取m＝n＝17个像素点；

步骤7：寻找出平滑区域均值滤波后的图像H的彩色噪声点，再对彩色噪声点进行平滑处理，其中，彩色噪声点的寻找方法为：

如果满足以下条件，则点(i2,j2)为彩色噪声点：

$\left\{\begin{array}{c}|1.140\&times;H_V(i2,j2)|<\mathrm{\&delta;}\\ |2.032\&times;H_U(i2,j2)|<\mathrm{\&delta;}\\ |0.394\&times;H_U(i2,j2)+0.581\&times;H_V(i2,j2)|<\mathrm{\&delta;}\end{array}\right.$

其中，δ表示判断彩色噪声的阈值，取δ＝30，H_U(i2,j2)和H_V(i2,j2)分别表示平滑区域均值滤波后的图像H在点(i2,j2)处的U值和V值，即饱和度和色调值，

对彩色噪声进行平滑处理的具体方法为：缩小H_U(i2,j2)和H_V(i2,j2)的值，得到平滑区域均值滤波后的图像H在点(i2,j2)处的缩小后的饱和度H_U2(i2,j2)和色调值H_V2(i2,j2)：

H_U2(i2,j2)＝H_U(i2,j2)/k,

H_V2(i2,j2)＝H_V(i2,j2)/k

其中，k表示阈值，根据平滑区域均值滤波后的图像H在点(i2,j2)处亮度值H_Y(i2,j2)的不同选取不同的值：

$k=\{\begin{array}{cc}2.0,& 0\&le;H_Y(i2,j2)\&le;40\\ 1.7,& 40<H_Y(i2,j2)\&le;60\\ 1.4,& 60<H_Y(i2,j2)\&le;80\\ 1.1,& 80<H_Y(i2,j2)\&le;255\end{array};$

步骤8：将经步骤7处理后的彩色图像从YUV彩色空间转换到RGB彩色空间：

$\left[\begin{array}{c}R\\ G\\ B\end{array}\right]=\left[\begin{array}{ccc}1& 0& 1.140\\ 1& -0.394& 0.581\\ 1& 2.032& 0\end{array}\right]\left[\begin{array}{c}Y\\ U\\ V\end{array}\right]$

并将转换后的图像输出。