数字图像处理的课件

1、精品文档第五章数字图像的预处理技术5.1 图像退化模型与噪声分类5.2 图像噪声的空间域滤波方法1 .2.1针对加性噪声的空间域滤波方法(1)自适应均值滤波器2 .自适应中值滤波器3 .自适应梯度倒数加权滤波器4 .2.2空间域滤波方法的MATLAB实现5.3 图像空间域滤波新方法探索5.3.1 中值滤波算法的研究1噪声检测2基于£1数据逼近和细节保护规则函数的凸面代价函数与噪声恢复3方法性能的仿真实验评价5.3.2 均值滤波算法的研究1基于RADON变换的图像主纹理方向分析2概率统计模型3方法性能的仿真实验评价5.4 周期噪声与频域滤波5.5 小波滤波5.5.1 小波域去噪模型5.

2、5.2 阈值的估计5.5.3 小波滤噪的MATLAB实现5.6 偏微分方程图像去噪5.7 边缘检测5.7.1 边缘检测的梯度算子Roberts算子1 g(x,y)-Td(x,y)=<(5-116)0g(x,y)<t2 Prewitt算子3 Sobel算子4 LOG算子5 Canny算子6 .7.2基于梯度算子白边缘检测MATLAB编程7 .7.3Hough变换与线检测第五章数字图像的预处理技术数字图像的预处理技术是图像处理与分析的基础，其处理结果对后序图像的处理与分析结果有非常大的影响。图像的预处理技术主要包括图像复原、图像增强两种操作，图像复原是通过使用退化现象的先验知识试图重建

3、或恢复一幅退化图像，主要包括几何畸变的修正、图像滤噪等，本书在第三章介绍了几何畸变的修正的基本内容。本章主要介绍图像滤噪的相关内容。图像增强是根据图像处理或分析的主要目标，以增强图像中相关的部分信息为主要目标的操作，主要包括图像锐化、边缘检测及基于图像直方图的图像增强等。在第三章介绍了基于图像直方图的图像增强技术。本章主要介绍图像图像锐化、边缘检测的相关内容。5.1图像退化模型与噪声分类在图像生成与传输过程中，由于成像系统与通讯信道各种因素的影响，可能使图像质量降低。这种图像质量的降低，被称之为“退化”。图像模糊、失真、附加噪声等都是图像退化的典型表现。图像复原是利用退化过程的先验知识使已被退

4、化的图像恢复本来的面目，从而改善图像的质量。图像复原可以看成是退化的逆过程，因此了解图像的退化是进行图像复原的基础。一般来说，图像的退化过程可以描述为如下的数学表达式：(5-1)gx,y=Hfx,y式中，f(x,y)是原图像，H是综合退化因子，g(x,y)是退化后的图像。通常，图像f(x,y)可以表示成：f(xy=f_JJa,B)qxa,y户也#(5-2)式中，f(a,B发像素点的特征函数，d(x-a,y-P)为冲击响应。为了简化问题，假定成像系统是线性的。实际上，一般情况下，退化因子对图像的各个像素点所产生的影响是相同的，呈现出线性的规律。这样，式(5-2)可以表示为：g(xy=H_/vc(

5、f，B)q以5y户dd"bo=.fi：-,-H、xT,ydd：(5-3)令h(x,y,a,P户H6(xa,y-P),称为成像系统的冲击相应，也叫做点扩散函数(PSF)。若成像系统是线性移不变系统，则退化图像可以用f和h的卷积表示如下：g(xy=*fa,p)(h-xa,-yl|cdd(f)xyJhxy(5-4)此外，如果受到加性噪声n(x,y)的干扰，则图像的退化模型可表示为见图5.1：g(xy=f,x)y(h,x5n)xy(5-5)如果噪声类型为乘性噪声，则图像的退化模型可表示为：gx,y=fx,yhx,y1nx,y(5-6)从数学上来说，图像退化是卷积运算和加噪运算的组合，如果可以

6、施行去卷积和去噪声的操作，退化的图像就可以复原见图2.1。由此图像滤噪是图像复原的重要手段。f(x,y)* H(x, y)g(x, y)n（x,y）图5.1图像退化模型噪声可以理解为不可预测，只能用概率统计方法来认识的随机误差，因此将图像噪声看成是多维随机过程是合适的，描述噪声完全可以借用随机过程及其概论分布函数。图像噪声按其产生的原因可分为内部噪声与外部噪声；从统计特性可分为平稳噪声与非平稳噪声，平稳噪声统计特性不随时间变化，而非平稳噪声统计特性是时间的函数；按噪声与信号的关系可分为加性噪声与乘性噪声。基于前面的假设，所关心的空间噪声描述符是上面所提及模型的噪声分量灰度值的统计特性。它们可以

7、被认为是由概率密度函数（PDF）表示的随机变量，下面是图像处理常见的噪声分布类型。5.1.1 数字图像常见的噪声数字图像常见的噪声分布类型有Gaussian噪声分布、瑞利噪声分布、伽马（厄兰）噪声、指数分布噪声、对数正态分布噪声、脉冲噪声分布等。一般来说，高斯噪声可以近似描述在低照明水平图像传感器成像产生的噪声，椒盐噪声可以模拟不完备的开关设备产生的噪声，瑞利噪声产生于波段成像，而指数与伽马（厄兰）噪声可以描述激光成像产生的噪声。照像乳胶中的银粒大小产生的噪声可以用对数正态分布噪声描述。图5.2给出了Lena图片感染各种噪声的结果。1 .Gaussian噪声分布这是一种常见的噪声模型，大多数噪

8、声可近似认为满足高斯分布，而且高斯分布比较容易进行数学分析。设随机变量z满足高斯分布，则其概率密度函数为：(5-7)精品文档其中，z表示图像的灰度值，N表示期望值，仃表示z的均方差。如果z服从式（2-7）分布时，其值有70%落在（N。）1范围内，且有90%落在（N2仃）,（N+加）1范围内。2 .瑞利噪声分布若随机变量z满足瑞利分布，则其概率密度函数为:z - az : a(5-8)z的均值和方差为:N=a+Jnb/4=2_,bbb（4-n）43 .伽马（厄兰）噪声伽马噪声的PDF由下式给出:p(z )= «(b -1J0-azez :二 0(5-9)其中a>0,b为正整数。密

9、度的均值和方差为:尽管式（2-9）经常被用来表示伽马密度，严格地说，只有分母为r（b）时伽马函数才是正确的。当分母如上式所示时，该密度近似为厄兰密度。4 .指数分布噪声az指数噪声的PDF可由下式给出：z_0P（z尸(5-10)其中a>0。概率密度函数的期望值和方差为:指数分布的概率密度函数是当b=1时爱尔兰概率分布的特殊情况。5 .均匀分布噪声均匀分布噪声的概率密度可由下式给出：其他(5-11)概率密度函数的均值和方差可由下式给出:ab""2",2b-a二126 .脉冲噪声分布脉冲噪声的概率密度函数如下：Paz=apZ=pzb(5-12)o其他如果b>

10、;a,灰度值b在图像中将显示为一个亮点，相反，a将显示为一个暗点。若pa或pb为零，则脉冲噪声称为单极脉冲。如果pa和pb均不可能为零，尤其是它们近似相等时，脉冲噪声值类似于随机分布在图像上的胡椒和盐粉微粒。由于这个原因，双极脉冲噪声也称为椒盐噪声。有时也称为散粒和尖峰脉冲。脉冲噪声可以是正，也可以是负。因为脉冲干扰通常与图像信号的强度相比较大，因此，在一幅图像中，脉冲噪声总是数字化为最大值(纯白或纯黑)。这样，通常假设a和b是饱和值，从某种意义上看，在数字化图像中，它们等于所允许的最大值和最小值。由于这一结果，负脉冲以一个黑点(胡椒点)出现在图像中，正脉冲以一个白点(盐点)出现在图像中。对于

11、一个8位图像，这意味着a=0,b=255。7 .对数正态分布噪声对数正态分布噪声概率密度函数如下：p(z)=1_ln(z)-a2/(2b2)21bz(5-13)5.1.2 用MATLAB模拟各类噪声MATLAB工具中提供了imnoise()对图像添力口噪声，其调用格式为：>>g=imnoise(f,type,parameters);其中f为待感染的图像，g为感染后的图像，参数type指定感染噪声类型，取值为gaussian'感染的是高斯噪声，取值为salt&pepper',感染的是椒盐噪声，如果取值为locakar',则感染的是均值为零的高斯噪声，取

12、值为speckle',感染的为乘性噪声，取值为posson',则感染的是泊松噪声。参数parameters随type参数取值不同而不同，如果是高斯噪声，参数parameters指定其均值与方差，具体格式为:>>g=imnoise(f,'gaussian',m,var);其中，m为高斯噪声的均值，var为方差。由于imnoise()函数要求被感染图像必须是取值为0,1之间的double类型，所以对应参数应考虑这一点，例如，对于256灰度图像，如果感染均值为64、方差为100的高斯噪声，参数m=64/256,var=100/2562。如果类型为椒盐噪声，

13、则需指定噪声感染的概率密度，调用格式为：>>g=imnoise(f,'salt&pepper',d);其中，d为感染率，如果省略d,则默认值为0.05。type取值为locakar',则感染的是均值为零的高斯噪声，调用有如下两种方式：>>g=imnoise(f,'localvarv);或者为：>>g=imnoise(f,'localvar',image_intensity,v);第一种方式v为与图像f大小相同的数组，可以对每个像素规定其方差值。第二种方式v为与图像f的灰度值取值范围大小相同的数组，int

14、ensity为对应灰度值，可以对图像的每个灰度值规定其方差值。type取值为speckle则用方程g=f+nxf将均值为0,给定方差的均匀分布的乘性噪声感染到图像f中，调用格式为：>>g=imnoise(f,'speckle',var);0.04。a=10(I的瑞利噪声对数正态分布噪声噪声其中，参量var指定均匀分布噪声的方差，缺省值为厄兰噪声声差为0.04的乘性噪图5.2 Lena图片感染各种噪声的结果例5.1对Lena利用MATLAB编程实现感染均值0,方差0.01高斯噪声。程序如下：f=imread('E:matlab7lena.bmp');g

15、=rgb2gray(f);g1=im2double(g);g1=imnoise(g1,'gaussian',0,0.01);imshow(g1,);执行结果见图5.2(b)例5.2对Lena利用MATLAB编程实现感染噪声密度为0.05的椒盐噪声。程序如下：f=imread('E:matlab7lena.bmp');g=rgb2gray(f);g1=im2double(g);g1=imnoise(g1,'salt&pepper',0.05);imshow(g1);执行结果见图5.2(c)虽然利用MATLAB提供的imnoise()可以方便

16、的模拟一些图像噪声，但不能模拟所有图像噪声类型。哪么如何利用MATLAB模拟已知概率密度函数的所有图像噪声类型，根据概率论的一个著名结论：即若w是一个在区间(0,1)内均均分布的随机变量，则可以通过下式获得具有指定概率累积分布函数(CDF)的随机变量z：(5-14)z=Fz(w)其中Fz(w)为随机变量z的累积分布函数CDF。哪么如何产生一个在区间(0,1)内均均分布的随机变量，MATLAB提供的随机函数可以产生一个在区间(0,1)内近似均均分布的随机变量，其调用格式为：>>w=rand(n,m);其中，参量n,m为产生的随机变量组成的矩阵的大小，另外一个生成随机变量的函数为ran

17、dn(),其生成取值为(0,1)均值为0、单位方差的正态分布随机变量，其调用格式为：>>w=randn(n,m);参量n,m与rand()函数相同。由此，只有已知噪声模型的概率分布函数，就可以计算出其累积分布函数，然后利用式(5-14)可以产生符合概率分布函数随机噪声，如果是加性噪声，即可以利用下式完成对图像x的感染：x=x+r(5-15)其中，r为产生的与图像x尺度相同的随机噪声。如对于瑞利噪声，根据式(5-8)可以求得其累积分布函数为：1»5z.aFz(z)=4(5-16)0,z:二a根据式(5-14),得方程为：2z_a)2/b,、1-e=w(5-17)解方程得：z

18、=a+J-bln(1-w)(5-18)则生成瑞利噪声的MATLABm函数如下：Functionr=rayleigh(m,n,a,b)%rayleighgeneratesanarrayofrandomnumberswithrayleighPDF%withm-by-n,whoselementsarerandomnumbersbetween0,1with%parametersaandb.Ifm=n=1,itgeneratesasinglerandomnumber.r=a+(-b*log(1-rand(m,n).A0.5end5.2图像噪声的空间域滤波方法图像退化是卷积运算和加噪运算的组合，如果可以施

19、行去卷积和去噪声的操作，退化的图像就可以复原。当图像退化原因仅仅是感染噪声时，式(5.2)变为:(5-19)gx,y=fx,ynx,y由于噪声项n(x,y)是未知的，从g(x,y)减去它们不是一个现实的选择。不同的噪声处理方法也不同，存在加性噪声时，可以选择空间域滤波方法。如果是周期性噪声，频域处理方法较为适合。5.2.1 针对加性噪声的空间域滤波方法对于感染加性噪声的图像，常用的空间域滤波方法有均值滤波器与顺序统计滤波器，均值滤波器是一种线性滤波器，而顺序统计滤波器是一种非线性滤波器。1均值滤波器均值滤波器主要有算术均值滤波器、几何均值滤波器、谐波均值滤波器、逆谐波均值滤波器。(1)算术均值

20、滤波器算术均值滤波器是最简单的均值滤波器，如果Sy表示中心在(x,y)的点、尺寸为mn的矩形子图像窗口的坐标组。算术均值滤波算法过程就是计算由Sxy定义的区域中被干扰图像g(x,y)的平均值。在任意点(x,y)处复原图像1?的值就是用Sxy定义的区域的像素计算出来的算术平均值。即：11?x,ygs,tmns.t.Sxy这个操作也可以用系数为1/mn的卷积模板来实现。(2)几何均值滤波器用几何均值滤波器复原一幅图像由如下表达式给出:1mn?x,y尸ng(s,t)(s,t卢y(5-20)(5-21)其中每一个复原像素由子图像窗口中像素点的乘积并自乘到1mn次哥得出。几何均值滤波器所达到的平滑度可以

21、与算法均值滤波器相比，在滤波过程中可以丢失更少的图像细节。(3)谐波均值滤波器?x,y =mns,t $ g s,t(5-22)谐波滤波器对于“盐”噪声、高斯噪声效果比较好，但是不适合于“胡椒”噪声。(4)逆谐波均值滤波器逆谐波均值滤波器对一幅图像进行恢复操作依据如下公式:f?x, y =" g s,t Q 1|S,t -Sxyg g(s,t Qs,t 衿y(5-23)其中Q称为滤波器的阶数。这种滤波器适合减少或者在实际中消除椒盐噪声，当Q为正数时，滤波器用于滤除“胡椒”噪声，当Q为负数时，滤波器用于滤除“盐”噪声。当 Q = 0时，逆谐波均值滤波器退化为算术均值滤波器，当 Q=-1

22、时，它退化为谐波均值滤波器。2顺序统计滤波器顺序统计滤波器的响应基于由滤波窗口包围的图像区域中像素点的排序结果，滤波器在任意点的响应由排序结果决定。主要有中值滤波器、最大值和最小值滤波器、中点滤波器、修正后的阿尔法均值滤波器等 。(1)中值滤波器中值滤波器是最著名的顺序统计滤波器，它的输出响应是窗口中像素排序结果的中值:? x,y =med i agns,tIs,t - Sxy(5-24)谐波均值滤波器滤噪操作如下表示：像素的原始值包含在中值的计算结果中。因为它对很多随机噪声都有很好的去噪能力，且在相同窗口尺寸下比线性滤波器引起的模糊少，所以中值滤波器应用很普遍。中值滤波器对滤除单极或双极脉冲

23、噪声效果非常好。(2)最大值/最小值滤波器中值滤波器采用的是排序结果中的中值，而从基本的统计学上讲排序还有其它很多的可能性。当使用排序中的最大值时，得到最大值滤波器如下:?x,y )= max(s,t 产yg s,t(5-25)这种滤波器在发现图像中的亮点非常适用。同样，因为胡椒噪声是非常低的值,最大值滤波器是选择子图像Sy中的最大值，可以滤除胡椒噪声。如果选择排序中的最小值，即为最小值滤波器：?仅 y 肝in g(s，tbpt JSxy这种滤波器发现图像中的暗点非常有用，使用最小值滤波器，可以滤除盐噪声。(5-26)中点滤波器中点滤波器是取滤波器涉及范围内的最大值和最小值之间的中点:?x,

24、y =2- 1max Z g(s,t )+min 工 g(s,t )这种滤波器结合了顺序统计和求平均，对于高斯和均匀随机分布噪声有较好的效果。(5-27)(4)修正后的阿尔法均值滤波器假设在Sxy邻域内去掉g(s,t)最高灰度值的d/2和最低灰度值的d/2,用gr(st,)代表 剩余的mn-d个像素，根据这些像素的平均值形成的滤波器被称为修正后的阿尔法均值滤 波器。小=rnbs：sgr s,t(5-28)其中d可以取0到mn -1之间的任意整数。当 d =0时，修正后的阿尔法均值滤波器退化 为算术均值滤波器；当 d =(mn-1 H2时，修正后的阿尔法均值滤波器退化为中值滤波器。修正后的阿尔法

25、均值滤波器当 d取其它值时，结合了均值滤波器与中值滤波器的优点，对多种噪声感染形成的混合噪声非常有用。4自适应滤波器以上提到的滤波器应用于图像时，并没有充分考虑图像的局部特征的变化，自适应滤波器基于由mn矩形窗口 Sxy定义的区域内图像的统计特性，具有更好的滤波性能。(1)自适应均值滤波器随机变量最简单的统计参量是均值和方差。因为它们是与图像状态紧密相关的数据，这些参数是自适应滤波器的基础。均值给出了区域图像的灰度平均值的度量,而方差给出了这个区域的对比度的度量。滤波器作用于局部区域Sxy,在任何点(X,y口的响应基于以下四个量：g(Xy表示噪声图像在点(X,y)上的值2仃刈表不干扰f(X,y

26、)以形成g(x,y)的噪声万差mL表示在Sxy上像素点的局部均值2Dl表示在Sxy上像素点的局部方差滤波器的输出响应如下：2如果。”为零，滤波器应该简单地返回g(x,y)的值。(在零噪声情况下g(x,y)等于f(x,y)如果局部方差和仃2是高相关的，那么滤波器要返回一个g(x,y)的近似值。如果两个方差相等，滤波器返回区域&y上的算术平均值。(在局部面积与全部图像有相同特性的条件下，局部噪声简单地用求均值来降低)为了获得？(x,y),基于这些假设的自适应表达式可以写为：_2?x,y)尸gx,y-1gx,y-mL1-L(5-29)在实际中对于式(5-29)的实现应构建一个测试，以便如果条

27、件仃彳下。2发生，把比率设置成1。这使该滤波器为非线性的，但它可以防止由于缺乏图像噪声方差的知识而产生的无意义的结果(即负灰度值)。另一个方法是允许负灰度值，在最后重新标定灰度值。这个结果将损失图像的动态范围。2 .自适应中值滤波器中值滤波器在冲激噪声空间密度不大(根据经验，噪声感染概率小于0.2),性能很好。自适应中值滤波器可以处理更大概率的冲激噪声，另外在平滑非冲激噪声时可以更好地保存细节，这是传统中值滤波器做不到的。符号定义如下：Zxy表示在坐标(x,y)上的灰度值Sxy中灰度最小值Sxy中灰度最大值Sxy中灰度中值Sxy的最大值A层和B层，如下:Zmin表示窗口Zmax表小窗口Zmed

28、表示窗口Smax表布窗口自适用中值滤波器的工作在两个层次，A层：if(Zmin<Zmed<Zmax)and(Zmin<Zxy<Zmax)Zxy=Zxy;elseif(Zmin<Zmed<Zmax)and!(Zmin<Zxy<Zmax)Zxy=Zmed;elseif!(Zmin<Zmed<Zmax)increasethewindowsize;ifwindowsize>SmaxZxy=Zxy;elserepeatfromthebegin;B层：if(Zmin<Zmed<Zmax)and(Zmin<Zxy<Zm

29、ax)Zxy=Zxy;elseif(Zmin<Zmed<Zmax)and!(Zmin<Zxy<Zmax)Zxy=Zmed;elseif!(Zmin<Zmed<Zmax)increasethewindowsize;ifwindowsize>SmaxZxy=Zxy;elserepeatfromthebegin;该算法的主要目的是滤除椒盐噪声，平滑其它非冲激噪声，减小诸如边缘细化或粗化的失真。Zmin和Zmax的值进行统计后被算法认为是类冲激式的噪声成分，即使它们在图像中不是最高和最低的可能象素点。从算法中可以看到，A层的目的是决定中值滤波器的输出Zmed是

30、否是一个脉冲(黑或白)。如果条件Zm<Zmed<Zmax有效，则Zmed就不是脉冲。在该情况下就转到B层，检查一下窗口中心点Zxy本身是否是一个脉冲，如果条件Zm<ZXy<Zmax为真，则Zxy也不是脉冲。在这种情况下，算法输出Zxy本身。如果条件Zm<Zxy<Zmax为假，像素值Zxy是一个极值且算法输出中值Zmed。假设在A层找到一个脉冲，算法就会扩大窗口尺寸并重复A层。如果达到了最大的窗口尺寸还是能找到脉冲，算法就返回Zxy值。注意：如果噪声概率过小或者Smax在允许的范围内越大，过早推出条件的可能性就越小。这应该是合理的，随着脉冲密度的增大，需要更大

31、的窗口消除尖峰噪声。3 .自适应梯度倒数加权滤波器图像中往往存在相邻区域的变化大于区域内部的变化，在同一区域中间像素的变化小于边缘像素的变化。梯度值正比于邻域像素灰度级之差，即图像变化缓慢区域，梯度值小，反之则大。如果取梯度的倒数，以梯度倒数作权重因子，则区域内部的相邻点权重就大于外部像素点的权重。这种滤波器的主要贡献来自于区域内部的，对于保持图像的细节会更有效。下面以3M3的窗口为例，如下图所示。定义：g(i, j,m,n)=1D(i m, j n) -D(i, j)(5-30)D(i-1,j-1)D(i-1,j)D(i-1,j+1)D(i,j-1)D(i,j)D(i,j+1)D(i+1,j

32、-1)D(i+1,j)D(i+1,j+1)并规定D(i,+m,j+n)D(i,j)=0时，g(i,j,m,n)=0。则窗口内每一个像素点对应的权值为：w(i,j,m,n)=1g(i,j,m,n)(1-r),m,n=0(5-31)二二g(i,j,m,n)m-nJ最后的倒数加权滤波输出为：11y=D(i,j)l一二w(i,j,m,n)D(i,j,m,n)(5-32)m-n-式中r是中心像素占有的权值，可以通过调节r的值来控制周围各点对滤波作用。5.2.2空间域滤波方法的MATLAB实现从上一节可以看出：空间域的滤波方法主要有线性与非线性两种，各类均值滤波是线性滤波器，而统计滤波器是非线性滤波器。均

33、值滤波的主要思路是将像素(x,y)特定邻域中每个像素与相应的系数想相乘，然后将结果累加，取平均做为点(x,y)的滤波结果。若邻域大小为mmn,则需要m父n系数，由这些系数组成的m父n的矩阵称为滤波掩模、滤波窗口或为滤波模板。滤波过程是将滤波掩模的中心逐点移动到图像的每个像素处，然后将滤波掩模的系数与对应像素相乘，取其平均值即为中心点像素的滤波结果，这类似于信号分析中的相关运算，如果滤波掩模以中心对称，则此过程也可以理解成为信号分析的卷积过程。因为卷积是线性移不变系统的数学抽象，所以此类滤波方法是线性的。1.线性滤波器的MATLAB实现MATLAB工具箱提供了imfilter()实现线性空间滤波

34、，调用格式为：>>g=imfilter(f,w,mode,boundary_options,size_option);其中，参数f为待滤波图像，w为滤波掩模矩阵，mode规定了计算方式，如果为corr',则使用相关运算，如果是'conv',则使用卷积运算，默认值为conv'。boundary_options决定了对边界像素的处理方式，如果为p,则以值p来填充图像边界，如果为replicate',则以外边界像素值填充图像边界，如果为symmetric',则通过镜像反射来生成边界像素值，如果为'circular',则将图像

35、边界看成为一个二维周期函数的一个周期来扩展。参数size_option决定输出图像的大小，如果为full则输出图像与扩展图像大小相同，如果是same',则输出图像与输入图像大小相同，缺省值为same'。滤波过程的计算结果是双精度类型，但imfilter函数会将结果转换为输入图像原有的数据类型，如果输入图像为整数型，则计算结果中超出取值范围的数据将以最大值或最小值代替，小数部分将会四舍五入。如果要求更高的精度，则在使用该函数之前，应将输入图像通过imdouble()或double转换为double类型。例5.3对例5.1中感染高斯噪声的Lena图像利用MATLAB编程应用均值滤波

36、器进行滤波。编程如下：>>f=imread('E:matlab7lena.bmp');> >g=rgb2gray(f);> >g=im2double(f);>>g1=imnoise(g,'gaussian',0,0.01);> >imshow(g1);> >w=ones(3,3)/9;>>d=imfilter(g1,w,'conv','replicate','same');> >figure;> >imsho

37、w(d,);该程序是生成3M3的滤波掩模，然后调用imfilter实现滤波，由于imfilter只是计算卷积运算，所以要准确计算均值应除以滤波掩模的元素数，因此在生成由1组成的滤波掩模w时，除以滤波掩模的元素数，执行结果见图5.3(b)。算术均值滤波器可以直接调用函数imfilter()写成，但是几何均值滤波必须先通过数学变换才能调用imfilter(),对式(5-21)两边求对数，有：、log(gs,t)log(?(x,y=isS(5-33)mn则实现几何均值滤波的程序如下：>>w=ones(3,3)/9;> >warningoff>>d1=exp(imf

38、ilter(log(g1),w,'replicate','same')；> >warningon> >figure;> >imshow(d1,);执行结果见图5.3(c)o根据式(5-22),谐波均值滤波器相当于用滤波掩模的元素数除以原图像的倒数与滤波掩模的卷积运算结果，对应程序如下：> >w=ones(3,3);>>d2=3*3./imfilter(1./(g1+eps),w,'replicate','same');> >figure;> >im

39、show(d2,);其中，eps是MATLAB的最小常量，引入原因是以防某些像素值为0造成分母为零，程序的执行结果见5.3(d)o根据(5-23),逆谐波均值滤波器相当于原图像的Q次哥与滤波掩模的卷积运算结果除以原图像的Q-1次哥与滤波掩模的卷积运算结果。，对应程序如下:>>w=ones(3,3);>>q=3;>>d3=imfilter(g1).Aq,w,'replicate','same');>>d4=imfilter(g1).A(q-1),w,'replicate','same'

40、);>>d3=d3./d4;>>figure;>> imshow(d2,);执行结果见图5.3(e)。滤波结果滤波结果器滤波结果图5.3感染均值为0,方差为0.01高斯噪声的Lena图片用线性滤波器的操作结果为了配合函数imfilter()的使用，MATLAB还提供了产生线性空间滤波器滤波掩模的函数fspecial(),其调用格式为：>>w=fspecial(typeparameters);其中，参数type'表示滤波器的类型，parameters规定了对应类型滤波器的参数。type'可以取以下值：(1) 'average

41、'Average规定产生矩形均值滤波器的滤波掩模，其调用格式为：>>w=fspecial(average',r,c);结果是产生rxc大小的均值滤波矩阵，如果省略参数r,c,则缺省值为3黑3。(2) 'disk'disk规定产生圆形均值滤波器的滤波掩模，其调用格式为：>>w=fspecial(disk',r);结果是产生半径为r的圆形均值滤波掩模，如果省略参数r,则缺省值为5。(3) 'gaussian'gaussian规定产生gaussian低通滤波器的滤波掩模，其调用格式为：>>w=fspecial

42、(gaussian",r,c,sig);结果是产生大小为rxc、标准差为sig的gaussian低通滤波器的滤波掩模，如果省略参数r,c,则缺省值为3M3,标准差sig的缺省值0.5。(4) 'laplacian'laplacian规定产生大小为3M3的拉普拉斯滤波掩模，其调用格式为:>>w=fspecial(laplacianalpha);(5-34)结果是产生3父3大小的拉普拉斯滤波矩阵，形状由参量alpha规定，缺省值0.5。拉普拉斯滤波矩阵与参量alpha=a的关系是：*a1-aot1+ot1+a1+0(1-a-41-aw=11+a1+aa1-aa

43、1+a1+alog'log规定产生高斯-拉普拉斯滤波矩阵，其调用格式为：>>w=fspecial(lOg，,r,c,sig);结果是产生rc大小的高斯-拉普拉斯滤波矩阵，高斯的标准差由参量sig给出缺省值为0.5,如果省略参数r,c,则缺省值为5M5。(6) 'prewitt'prewitt规定产生3父3大小的prewitt滤波掩模，其调用格式为：>>w=fspecial(prewitt');(7) 'sobel'sobel规定产生3乂3大小的sobel滤波掩模，其调用格式为：>>w=fspecial(sobe

44、l');(8) 'unsharp'unsharp规定产生3M3大小的非锐化滤波掩模，其调用格式为：>>w=fspecial(unsharp",alpha);其中，参数alpha控制形状，取值为(0,1)，缺省值为0.2。(9) 'motion'motion规定产生一个图像运动滤波器，其调用格式为：>>w=fspecial(motion',len,theta);其中，参数alen是运动的像素数，缺省值为9,theta规定运动方向，即与水平方向的夹角，缺省值为0。例5.4对例5.1中感染高斯噪声的Lena图像利用MA

45、TLAB编程应用大小为5父5、方差为0.6的高斯滤波器进行滤波。编程如下：>>f=imread('E:matlab7lena.bmp');>>g=rgb2gray(f);>>g=im2double(f);>>g1=imnoise(g,'gaussian',0,0.01);>>imshow(g1);>>w=fspecial('gaussian',5,5,0.6);>>d=imfilter(g1,w,'replicate');>>figur

46、e;>>imshow(d,);Lena图片用大小为5M 5、方差为0.6的高斯程序的执行结果见图5.4。图5.4感染均值为0,方差为0.01高斯噪声的滤波器操作结果2.顺序统计滤波器的MATLAB实现MATLAB用于实现顺序统计滤波器的函数是ordM2(),其调用格式为：>>g=ordfilt2(f,order,domain);其中，参量f为输入图像，order指定了用邻域排序后的第多少位的像素做为输出图像的当前像素，排序的顺序是由小到大的升序。排序的邻域由参量domain给出。domain是取值为0或1,大小为mn的矩阵，参预排序的像素是在f中，以当前像素为中心，do

47、main矩阵中值为1的对应像素。例如：最小值排序滤波器可以用以下语句实现：>>g=ordfilt2(f,1,ones(m,n);其中，函数ones(m,n)是产生一个mn大小的、值为1的矩阵。最大值排序滤波器可以用以下语句实现：>>g=ordfilt2(f,m*n,ones(m,n);例5.5对例5.1中感染噪声密度为0.05的椒盐噪声的Lena图像利用MATLAB编程应用各类排序滤波器进行滤波。程序如下：f=imread('E:matlab7lena.bmp');> >g=rgb2gray(f);> >g2=imnoise(g,

48、'salt&pepper',0.05);> >imshow(g2);%最小值排序滤波器> >d1=ordfilt2(g2,1,ones(3,3);> >figure;> >imshow(d1);%最大值排序滤波器> >d2=ordfilt2(g2,3*3,ones(3,3);> >figure;> >imshow(d2);%中值滤波器> >d3=ordfilt2(g2,median(1:3*3),ones(3,3);> >figure;> >imsho

49、w(d3);%中点滤波器>>f1=ordfilt2(g2,1,ones(3,3),'symmetric');>>f2=ordfilt2(g2,3*3,ones(3,3),'symmetric');> >d4=imlincomb(0.5,f1,0.5,f2);> >figure;> >imshow(d4);%d=4的阿尔法均值滤波器g3=im2double(g2);> >d=4;> >f1=imfilter(g3,ones(3,3),'symmetric');>

50、; >fork=1:d/2f=imsubtract(f1,ordfilt2(g3,k,ones(3,3),'symmetric');end> >d=4;> >fork=(3*3-(d/2)+1):3*3f=imsubtract(f1,ordfilt2(g3,k,ones(3,3),'symmetric');end> >f=f./(3*3-d);> >figure;> >imshow(f);结果结果滤波器滤波结果图5.5感染噪声密度为0.05的椒盐噪声的Lena图像应用各类排序滤波器进行滤波结果程序

51、中函数d4=imlincomb(a,f1,b,是求两幅图像的线性加权和，权因子为a、b,函数imsubtract(f1,f2);是求两幅图像的差。为了使用的方便，MATLAB提供了专门的中值滤波函数medfilt2(),其调用格式为：>>g=medfilt2(f,mn,padopt);其中，参量f为待滤波图像，mn指定一个大小为mxn的邻域，缺省值为3M3,padopt指定了边界图像的填充模式，'zeros'用零填充，'symmetric'以镜像反射方式填充，'indexed'根据类型进行填充，如果f是double型，以1填充，否则以

52、零填充，缺省为zeros'。3.自适应中值滤波器的MATLABM函数中值与均值滤波器对图像的任何像素执行相同的操作，会损伤图像的细节，自适应中值滤波器是根据图像像素的具体情况，执行不同的操作，对于像素点，直接输出，对于噪声点，则以其窗口的中值代替，而且根据图像的局部情况采取不同大小的滤波窗口。自适应中值滤波器的M函数如下：functionf=amedian(g,smax)if(smax<=1)|(smax/2=round(smax/2)|(smax=round(smax)error('smaxmustbeanoddinteger>1');endm,n=siz

53、e(g);f=zeros(m,n);ap=false(m,n);fork=3:2:smaxzmin=ordfilt2(g,1,ones(k,k),'symmetric');zmax=ordfilt2(g,k*k,ones(k,k),'symmetric');zmed=medfilt2(g,kk,'symmetric');pub=(zmed>zmin)&(zmax>zmed)&ap;zb=(g>zmin)&(zmax>g);opzxy=pub&zb;opzmed=pub&zb;f(op

54、zxy尸g(opzxy);f(opzmed尸zmed(opzmed);ap=ap|pub;ifall(ap(:)break;endendf(ap)=zmed(ap);end其中，函数false(m,n);是生成大小为mn的值为0的矩阵。果果的椒盐噪声图5.6感染各种噪声密度的椒盐噪声的Lena图像应用自适中值滤波器滤波结果图5.6给出了感染不同密度椒盐噪声Lena图像应用自适中值滤波器滤波结果，对照图5.6与图5.5,可以看出：在小感染密度下，自适中值滤波器会获得比较好的滤波结果，但随着噪声密度的增加，其滤噪性能会明显下降。5.3图像空间域滤波新方法探索图像滤噪问题仍然是数字图像研究的热点问题

55、之一，空间域滤波方法的现状与研究方法有哪些特点，如何评价图像滤噪方法性能，本节结合两个研究实例进行阐述。5.3.1 中值滤波算法的研究中值滤波在低椒盐噪声感染率的情况时具有较好的滤噪能力和计算效率，然而当噪声感染率较高时，中值滤波器的滤波效果明显下降1。因此出现了许多改进型的中值滤波器，如自适应中值滤波器2、噪声自适应软开关(NASM)3、混合神经模糊滤波器4等被称为开关型的滤波器。这些滤波器的特点是首先对待处理的象素点进行判断，对于噪声点直接以其周围的象素点中值代替，这种方法并没有充分利用图像的局部特征，所以也会造成图像细节和边缘信息的损失，因此出现了一些利用噪声点周围局部特征来估计噪声点所

56、对应象素值的中值滤波算法，如基于直方图的模糊加权滤波器(HFF)5、中心加权中值滤波器(CWM)6等。此类方法在自适应中值的基础上，根据周围象素点加权求均值得到感染噪声点的输出，也会损失比掩膜更小的细节与边缘信息。为了在滤除噪声的同时，最大限度的保护图像的细节与边缘信息，基于二次噪声检测和细节保护规则函数的图像滤波算法将滤除噪声过程分为两个阶段7:噪声检测和噪声滤除，在噪声检测阶段，算法先通过自适应中值滤波器的原理检测出噪声点和象素点，然后使用模糊函数对噪声点进行二次判断，在噪声替换阶段，针对噪声点自适应地选择其周围最佳的象素点作为细节保护函数的迭代初值。当噪声集合的凸面代价函数值达到最小时得

57、到噪声集合的最佳恢复。1噪声检测基于二次噪声检测和细节保护规则函数的图像滤波算法噪声检测分为两个步骤，首先采用自适应中值滤波的原理进行第一次检测。令Sijw表示大小为ww,中心点是y(i,j)的窗口，Wmax父Wmax代表最大窗口。对于图像中的每一个像素点y(i,j)作如下操作：(1)设置初始值w=3(2)计算窗口Si/W中像素点灰度值的最大值Sjmax,w、最小值Si,jmin，W>中值Si,jmed，w(3)如果&jminw<s,med,wSiax成立，则转到步骤5,否则令w=w2(4)如果wEwmax成立，则转到步骤2,否则y(i,j)判断为噪声点如果Sijnw<

58、;y(i,j/Sijmaw,成立，则y(i,j)为像素点，否则y(i,j)判断为噪声点通过算法的原理可以看出有些特殊区域的象素点有可能被误判为噪声点(例如：偏黑区域的黑点或偏白区域上的白点)。为了提高判断的准确度，算法在第一步检测结果的基础上，对噪声点进行二次判断，二次判断根据包含噪声点(i,j)的检测窗口内未感染噪声的象素点得到模糊隶属度函数，进行模糊判断。假设象素点(i,j)第一次检测时在ww窗口大小时被判断为噪声点，那么就通过窗口中的信息点获得局部隶属度函数，进行二次模糊判断。窗口ww中信息点的均值为：II w,cJ i,j(5-35)其中Nw,c代表窗口w><w中信息点的个数，y(i+k,j+k)代表窗口中的信息点。信息点的方差：cww_仃=XZy。+k,j+k)-i,j,I(5-36)Nk=jw-根据大部分图像局部象素点灰度值分布的特点，本文采用高斯曲线型的隶属度函数，其数学表达式为:f(x户exp-(x-Nj/2仃2(5-37)当二次检测的象素点隶属度