基于MATLAB图像降噪方法的研究

ＸＸ大学ＸＸ学院20ＸＸ届毕业设计说明书实用文档第27页共28页1绪论1.1数字图像的定义图像是用各种观测系统以不同形式和手段观测客观世界而获得的，可以直接或间接作用于人眼并进而产生视知觉的实体【1】。人的视觉系统就是一个观测系统，通过它得到的图像就是客观景物在人心目中形成的影像。我们生活在一个信息时代，科学研究和统计表明，人类从外界获得的信息约有75%来自视觉系统，也就是从图像中获得的，例如照片、绘画、动画等等。客观世界在空间上是三维（3-D）的，但一般从客观景物得到的图像是二维（2-D）的。一幅图像可以用一个2-D数组来表示，这里和表示表示2-D空间中一个坐标点的位置，而则代表图像在点的某种性质的数值。例如常用的图像一般是灰度图，这时表示灰度值，它常对应客观景物被观察到的亮度。日常所见图像多是连续的，即f(x,y)的值可以任意实数。为了能用计算机对图像进行加工，需要把连续的图像在坐标空间XY和性质空间F都离散化【2】。这种离散化了的图像就是数字图像，可以用I(r,c)来表示。这里I代表离散化后的f，(r,c)代表离散化后的，其中r代表图像的行(row)，c代表图像的列(column)。这里的I,r,c值都是整数。在该设计中用f(x,y)代表数字图像。噪声对图像的影响由于图像在形成、传输、接收和处理的过程中,会受到外界环境、系统性能和人为因素等诸多方面影响,不可避免地存在噪声干扰,它使图像变质,影响图像的质量【3】。如果不对噪声进行及时处理,就会对后续的处理过程乃至输出结果产生影响,甚至可能得到错误的结论。1.3图像降噪的研究现状为了解决图像在生成和传输过程中因受到各种噪声的干扰和影响使图像降质，而对后续图像的处理(如分割、压缩和图像理解等)所产生不利影响的问题。人们逐步开始重视图像噪声的滤除问题。关于图像噪声去除的研究由来已久,人们根据实际图像的特点、噪声的统计特征和频谱分布的规律，以及不同的噪声源，发展了各式各样的去噪方法。常见的方法有邻域平均法，时域均值、中值滤波法、频域高斯低通滤波法等等，其余的方法往往是这几种方法的改进和变形。1.4图像降噪的研究意义噪声对图像处理十分重要，由于它影响图像处理的输入、采集、处理的各个环节以及输出结果的全过程。特别是图像的输入、采集噪声的抑制是十分关键的问题，若输入伴有较大的噪声，必然影响处理全过程及输出的结果【4】。因此一个良好的图像处理系统，不论是模拟处理还是计算机进行的数字处理，无不把减少最前一级的噪声作为主攻目标。然而由于影响图像噪声的原因不同，噪声类型不同，现有的去噪方法多，所以选择一种合适的去噪方法也尤为重要。该设计主要着重研究均值、中值、维纳去噪法，并通过MATLAB仿真结果比较这些方法的优劣。2图像噪声的基本概念2.1噪声的来源根据噪声产生的来源，大致可以分为外部噪声和内部噪声两大类。外部噪声是指从处理系统外来的影响，如天线干扰或电磁波从电源线窜入系统的噪声。内部噪声则有以下四种最常见形式。由光和电的基本性质引起的噪声。例如电流可以看作电子或空穴运动，这些粒子运动产生随机散粒噪声，导体中电子流动的热噪声，光量子运动的光量子噪声等。由机械运动引起的噪声。例如，接头振动使电流不稳，磁头或磁带、磁盘抖动等。元器件噪声。如光学底片的颗粒噪声，磁带、磁盘缺陷噪声，光盘的疵点噪声等等。系统内部电路的噪声。从噪声的分类方法来看是多种多样的。但综合来说，噪声是随机产生的量，所以又可以从统计数学的观点来定义噪声。凡是统计特性不随时间变化的噪声称为平稳噪声，而统计特性随时间变化的噪声称作非平稳噪声。2.2噪声的分类以上所讨论的各种类型的噪声反映在图像画面上大致可以分为椒盐噪声和高斯噪声【5】。椒盐噪声的幅值基本相同，但是噪声出现的位置是随机的。高斯噪声存在于图像中的每一点，且它的的幅值是随机分布的，其密度函数为高斯型的噪声称为高斯噪声，用MATLAB实现为：I=imread('01.jpg');%读入原图像J1=imnoise(I,'gaussian',0,0.01);%加高斯噪声J2=imnoise(I,'salt&pepper',0.04);%加椒盐噪声图2.1原始图像图2.2加噪（高斯）图像图2.3加噪（椒盐）图像2.3噪声的描述一般地，对噪声的描述采用统计意义上的均值与方差【6】。设图像信号的二维灰度分布为，则噪声可看作是对亮度的干扰，以来表示。噪声的均值表明了图像中噪声的总体强度，计算公式如下：(式2.1)噪声的方差表明了图像中噪声分布的强弱差异，计算公式如下：(式2.2)其中，分别为图像的行数和列数。2.4噪声的模型噪声的模型按照对信号的影响可分为加性噪声模型和乘性噪声模型两大类。设为信号，为噪声，在信号影响下的输出为，则有加性噪声模型为(式2.3)乘性噪声模型为(式2.4)乘性噪声作用下的输出是两部分的叠加，第二个噪声项信号受的影响，越大，则第二项越大，即噪声项受信号的调制。如光量子噪声、底片颗粒噪声都会因信号增大而增大。乘性噪声模型和它的分析计算都比较复杂，通常信号变化很小时，第二项近似不变，此时可用加性噪声模型来处理。一般情况下，总是假定信号和噪声是相互独立的【7】。该设计中主要讨论加性噪声的抑制方法。3均值滤波3.1均值滤波的原理所谓均值滤波实际上就是用均值替代原图像中的各个像素值。平均滤波的思想如下：设被噪声污染后的信号为，原始信号为，噪声为，则。均值滤波公式【8】如下：(式2.5)其中，是均值滤波器的长度。如当时，称为三值平均滤波；当时，称为十三值平均滤波。因为(式2.6)等式两边取数学期望（即统计均值），得(式2.7)如果噪声为高斯噪声，则有所以即理想情况下，所得到的结果中不含噪声。如果噪声为椒盐噪声，则有所以(式2.8)按照如上的定义可知，如果出现椒盐噪声，其幅值基本相同，但是数据中，还存在一些点上没有噪声，这样，可以将未出现噪声的点等同地认为噪声幅值为0，这样就有所以，经过均值滤波之后，中所包含的噪声强度低于。椒盐噪声的抑制与高斯噪声不同的是，它只是将某点出现的噪声强度让周围的数据平均分担了。3.2均值滤波的方法有了上述所介绍的一维均值滤波的原理，将其宽展到二维信号的处理，就可以实现对图像噪声的抑制。均值滤波的方法是，对待处理的当前像素，选择一个模板，该模板为其近邻的若干像素组成，用模板中像素的均值来替代原像素值的方法【9】。如图3.1所示，序号为0的是当前像素，序号为1～8的像素是其模板中的近邻像素。求模板中所有像素的均值，再把该均值赋予当前像素点，作为处理后图像在该点的灰度，即(式3.1)其中，s为模板，M为该模板中包括当前像素在内的像素总个数。考虑到数据分布的平衡性，模板一般选择为3×3，5×5，待处理像素放在模板的中心。为了使输出像素值保持在原来的灰度值范围内，模板的权值总和应维持为1。因此，模板与模板像素的乘积要除以一个系数（通常是模板系数之和），这个过程也被称为模板的归一化。模板的描述还可以采用矩阵的形式，如3×3的均值滤波器可以描述(式3.2)该模板的相应计算为(式3.3)用MATLAB实现为：I=imread('01.jpg');%读入原图J1=imnoise(‘噪声类型’);%加噪图像H=ones(3,3)/9;%选择3×3的模板B1=conv2(J1,H);%均值滤波4中值滤波从均值滤波的介绍可知，虽然均值滤波器对噪声有抑制作用，并且算法简单，但导致图像变得模糊是采用该方法所无法避免的。经过分析画面中噪声出现时所表现的形态可知，噪声点上的像素通常比周围非噪声点的像素要亮或暗。因此，可以设想，如果在噪声点像素周围寻找一个合理的值对它进行替代，在一定程度上应该可以获得较理想的滤波结果。基于以上考虑，设计的中值滤波就是一种有效的方法。4.1中值滤波原理中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术。在一维形式下，一维中值滤波器含有奇数个数据的滑动模板，对模板中的数据由小到大排序，取排在中间位置上的数据作为最终的处理结果【10】。下面，以一个简单的一维数据序列的滤波为例，介绍中值滤波原理。设模板的长度为5，模板中的数据设为｛10，15，45，20，25｝，则Med｛10,15,45,20,25｝=20(Med(…)表示取中值函数)。如果该模板中的数据为某个图像中的一个局部的数据，从数据的分布规律来看，原来模板中心位置上的像素值为45，较其周围的像素值大，画面上一定会出现一个突变的噪声点。经过中值滤波处理后的值为20，与周围的像素值差异不大，由此就得到了抑制噪声的效果。中值滤波的核心运算是将模板中的数据进行排序【11】，这样，如果一个亮点（暗点）的噪声，就会在排序过程中被排在数据序列的最右侧或者是最左侧，因此，最终选择的数据序列中间位置上的值一般不是噪声点的值，由此便可以达到抑制噪声的目的。数字图像是以二维数据来描述的，所以只要将原理中介绍的一维中值滤波推广到二维即可，具体方法如下：取某种结构的二维滑动模板，将模板内像素按照像素值的大小进行排序，生成单调上升（或下降）的二维数据序列【12】。二维中值滤波输出为：(式4.1)其中，分别为原图像和处理后图像。W为二维模板，通常选为3×3、5×5区域。4.2中值滤波的方法将模板在图中漫游，并将模板中心与图中某个像素位置重合；读取模板下各对应像素的灰度值；将这些灰度值从小到大排成一列；找出这些值中排在中间的一个；将这个中间值赋给对应模板中心位置像素。它的主要功能是让与周围像素灰度值的差比较大的像素改取与周围像素值接近的值，从而消除孤立的噪声点。用MATLAB实现为：I=imread('01.jpg');%读入原图J1=imnoise(‘噪声类型’);%加噪图像hood=3;%选择3×3的模板K1=medfilt2(J1,[hoodhood])；%中值滤波迭代中值滤波器迭代中值滤波器是传统中值滤波器的一种改进【13】，它的方法是使原中值滤波的结果再反复进行中值滤波处理，只是方法操作上的重复，其基本原理不变。MATLAB编程为：I=imread('01.jpg');%读入原图J1=imnoise(‘噪声类型’);%加噪图像hood=3;%选择3×3的模板K1=medfilt2(J1,[hoodhood]);%迭代1次K2=medfilt2(K1,[hoodhood]);%迭代2次K3=medfilt2(K2,[hoodhood]);%迭代3次K4=medfilt2(K3,[hoodhood]);%迭代4次5维纳滤波前面提到的两种降噪方法实质是通过对图像变换实现，维纳滤波是要对图像进行恢复实现[14]。图像恢复技术是图像处理领域一类非常重要的处理技术，与图像变换等其他基本图像处理技术类似，也是以获取视觉质量某种程度的改善为目的，所不同的是图像恢复过程实际上是一种估计过程，需要根据指定的图像退化的模型，对退化图像进行恢复，以取得未经过退化的原始图像。由于引起图像退化的因素众多一且性质各不相同，目前没有统一的恢复方法，众多研究人员根据不同的应用物理环境，采用了不同的退化模型、处理技巧和估计准则，从而得到了不同的恢复方法。该设计只分析维纳滤波法。5.1图像的退化模型与图像的矩阵表示在实际应用中，通常都假定传输系统是线性系统，原始图像通过系统。是综合所有退化因素得到的系统函数，称为成像系统的冲激响应或者点扩展函数(PSF)。图5.1所示的框图就是一个基本的退化模型，为实际得到的退化图像，为噪声模型。图5.1图像退化模型根据图5.1所示图像退化框图，退化模型可以表示为：(式5.1)但在实际应用中，处理的都是数字图像，所以对式5.1采用离散化形式进行表示，如下式所示：(式5.2)其中x=0,1,2,…,M-1,y=0,1,2,…,N-1。函数和分别是周期为M和N的函数，如果函数周期不是M和N，必须对其补零延拓，以避免卷积周期的交叠。数字图像一般有两种常用表示法：矩阵法和链码法【15】。本文研究的数字图像是以矩阵或数组的方式存储的。如果以列向量分别表示，和，如下式所示：,(式5.3)这样式5.3就可以写为：(式5.4)式中为维的矩阵，可以写作由个子矩阵组成的形式：(式5.5)而每个子矩阵都是由中的第j行构成的：(式5.6)由式5.5，式5.6可知和都是循环矩阵，即矩阵的第一行末尾和第二行开头元素都相同，矩阵末尾的尾元素与矩阵首行的头元素相同，所以矩阵是分块循环矩阵。现在图像恢复的主要任务就是根据和已知的和来估计。如果直接按照式5.4来计算，那是个非常繁杂的过程，数据量太大，而利用是分块循环的特点可以进行简化。其中，最常用的方法就是将分块循环矩阵对角化。可以证明，分块循环矩阵可以写成以下形式：(式5.7)式中D是对称正交对角阵,具体取值如下:(k=0,1,2,…,M-1)(式5.8)由此，对式5.4两边分别左乘,可得到下式：(式5.9)根据定义可知，，和分别是的傅里叶变换按照行堆叠而成的向量。这样，通过以下处理，就可以将空域中的复杂的方程求解问题转成频域中的简单计算，即对下式进行处理：=0,1,2,…,;=0,1,2,…,(式5.10)5.2维纳滤波原理维纳(Wiener)滤波可以归于反卷积算法一类，它是由Wiener首先提出的，应用于一维信号，并取得很好的效果。被引入二维信号处理领域后，取得相当满意的效果，尤其在图像复原过程中。在一般情况下，图像信号可近似地认为是平稳随机过程，维纳滤波将原始图像和对原始图像的估计看作为随机变量。假设和为和的自相关矩阵，在大多数图像中，邻近的像素点是高度相关的，而距离较远的像素点的相关性却较弱。通常，和的元素之间的相关不会延伸到20～30个像素的距离之外。因此，一般来说，自相关矩阵在主对角线附近有—个非零元素带，而在右上角和左上角的区域内将为零值。如果像素之问的相关是像素之间距离的函数，而不是它们位置的函数，可将和近似为分块循环矩阵，建立相应的数学模型，得到图像复原估计值：(式5.11)式中，u,v=0，1，…，N-1，。如果时，称之为维纳滤波器。当时，并不是在约束条件下得到的最佳解，即并不一定满足，若为变数，此式为参变维纳滤波器。利用最小均方误差估计，把维纳滤波应用到图像处理中，使图像估计与原始图像误差最小时有：(式5.12)我们定义为信噪比，但通常这两个值，尤其是噪声的功率谱难于得到，并且在实际应用中，对这个值的精度要求并不非常苛刻，因此我们可以选用一个正常数c来近似信噪比的倒数，因而维纳滤波器的估计值：(式5.13)式中是的傅立叶变换。维纳滤波的原理框图如图5.2，图像被线性模糊，加入噪声后称为加噪图像，经维纳滤波器W滤波后得到复原图像。图5.2维纳滤波的原理框图用MATLAB实现为：I=imread('01.JPG');%读入原图像J=imnoise(I,'gaussian',0,0.01);%加噪图像K=wiener2(J,[33]);%维纳滤波6降噪效果及讨论6.1各种降噪方法的降噪效果比较为了比较以上几种不同图像去噪方法的优劣,进行了大量实验。实验中,我们先人为地对不同图像(包括彩色图像及灰度图像)分别添加高斯白噪声、椒盐噪声。然后分别采取均值滤波、中值滤波、维纳滤波办法对图像去噪复原。给出的部分实验结果如下图6.1～6.8所示，同时根据图像质量评价标准（均方误差）公式：，(式6.1)其中设一幅图像原图像对应的象素点为，降噪图像为，根据S的大小评价图像质量好坏。以下将以上讨论的图像降噪法的降噪效果图做出比较，并给出各种降噪图像的均方差，以评价各滤波法对噪声的影响。（1）均值滤波效果（图6.1～图6.3)（a）高斯噪声（噪声密度0.01）（b）椒盐噪声（噪声密度0.04）图6.1（3×3模板）（a）高斯噪声（噪声密度0.01）（b）椒盐噪声（噪声密度0.04）图6.2（5×5模板）（a）高斯噪声（噪声密度0.01）（b）椒盐噪声（噪声密度0.04）图6.3（13×13模板）（2）中值滤波效果（图6.4～图6.6）（a）高斯噪声（噪声密度0.01）（b）椒盐噪声（噪声密度0.04）图6.4（3×3模板）（a）高斯噪声（噪声密度0.01）（b）椒盐噪声（噪声密度0.04）图6.5（5×5模板）（a）高斯噪声（噪声密度0.01）（b）椒盐噪声（噪声密度0.04）图6.6（13×13模板）（3）迭代中值滤波（用3×3模板）效果(a)迭代1次（b）迭代2次(c)迭代3次（d）迭代4次图6.7椒盐噪声（噪声密度0.3）（4）维纳滤波效果图6.8高斯噪声(噪声密度0.01)表6.1降噪图像均方误差3×3模板5×5模板13×13模板均值滤波（高斯密度=0.01）0.01740.03270.0819均值滤波（椒盐噪声密度=0.04）0.01820.03300.0815中值滤波（高斯密度=0.01）0.00520.00600.0169中值滤波（椒盐密度=0.04）0.00170.00420.0162迭代中值滤波（椒盐密度=0.3）迭代1次迭代2次迭代3次迭代4次0.01280.00640.00640.0065维纳滤波（高斯密度=0.01）0.00386.2结论通过对含有不同噪声类型的同一图像，用同一种滤波方法进行处理，可知：（1）由图6.1图～6.3及表6.1可知，均值滤波有一个非常致命的缺点，就是在求均值的计算中，会同时将景物的边缘点也同时进行均值处理，这样就使得景物的清晰度降低，画面变得模糊。虽然模板越大降噪效果越好，但是却模糊了图像边缘，在降低噪声的同时也降低了图像的质量。（2）由图6.4～图6.6可见,对于中值滤波并不是窗口越大滤波效果越好，窗口太大，在去除噪声象素的同时把一些好的像素也过滤掉了，使图像变得模糊不清，从而也破坏了图像的精度。（3）尽管维纳滤波方法是按最小原则导出的，在理论上更精确，但实验中发现维纳滤波虽能使噪声得到抑制，并较好保存图像的边缘和高频细节信息，但去噪后的图像仍略显模糊。通过对含有不同噪声类型的同一图像，用不同的滤波方法进行处理，比较可知：（1）由图6.1～图6.3的（b）与图6.4～图6.6的（b）比较及表6.1可知，均值滤波器对椒盐噪声的滤波效果不是很理想。因为椒盐噪声的幅值都是基本相同的，只是出现噪声点的位置是随机的，所以在统计意义下的噪声均值也不为0，因此，即使是理想情况也无法完全去除。但是从模板的含义来理解，经过均值处理之后，噪声部分被弱化到周围像素点上，所以得到的结果是噪声幅值减小，但是噪声点的颗粒面积同时变大。因此，处理椒盐噪声时用中值滤波效果好。（2）当椒盐噪声污染严重时（噪声密度0.3），传统中值滤波效果不佳。可采用迭代中值滤波法，由图6.7（a）～（d）及其均方误差显示可知去噪效果可达到一个最佳状态。（3）由图6.1（a）、图6.4（a）与图6.8比较可知，维纳滤波可较好地滤除高斯白噪声,图像也不致太模糊。附录本文所设计到的图像降噪去除法的对应程序：1用MATLAB语言实现的加噪图像I=imread('01.jpg');I=rgb2gray(I);figure;subplot(2,2,1);imshow(I);title('原始图像');J1=imnoise(I,'gaussian',0,0.01);subplot(2,2,2);imshow(J1);title('加噪(高斯）图像');J1=imnoise(I,'salt&pepper',0.04);subplot(2,2,3);imshow(J1);title('加噪（椒盐）图像');用MATLAB语言实现均值滤波选用3×3模板时的滤波效果I=imread('01.jpg');I=rgb2gray(I);figure;J1=imnoise(I,'gaussian',0,0.01);J2=imnoise(I,'salt&pepper',0.04);H=ones(3,3)/9;B1=conv2(J1,H);subplot(1,2,1);imshow(B1,[]);title('高斯噪声的滤波效果');B2=conv2(J2,H); subplot(1,2,2);imshow(B2,[]);title('椒盐噪声的滤波效果');fori=337:1:338forj=1:1:450I(i,j)=1;endendforj=449:1:450fori=1:1:338I(i,j)=1;endendI1=double(I);s1=0;s2=0;fori=1:1:338;forj=1:1:450;a=[I1(i,j)-B1(i,j)].^2;b=I1(i,j).^2;s1=s1+a;s2=s2+b;endends=s1/s2;disp(s');3用MATLAB语言实现中值滤波选用3×3模板时的滤波效果：I=imread('01.jpg');I=rgb2gray(I);figure;J1=imnoise(I,'gaussian',0,0.01);J2=imnoise(I,'salt&pepper',0.04);hood=3;K1=medfilt2(J1,[hoodhood]);subplot(1,2,1);imshow(K1);title('高斯噪声的滤波效果');K2=medfilt2(J2,[hoodhood]);subplot(1,2,2);imshow(K2);title('椒盐噪声的滤波效果');I1=double(I);B1=double(J1);s1=0;s2=0;fori=1:1:336;forj=1:1:448;a=[I1(i,j)-B1(i,j)].^2;b=I1(i,j).^2;s1=s1+a;s2=s2+b;endends=s1/s2;disp(s');迭代中值滤波器椒盐噪声（噪声密度0.3）I=imread('01.jpg');I=rgb2gray(I);J1=imnoise(I,'salt&pepper',0.3);hood=3;K1=medfilt2(J1,[hoodhood]);subplot(2,2,1);imshow(K1);title('迭代1次');hood=3;K2=medfilt2(K1,[hoodhood]);subplot(2,2,2);imshow(K2);title('迭代2次');hood=3;K3=medfilt2(K2,[hoodhood]);subplot(2,2,3);imshow(K3);title('迭代3次');hood=3;K4=medfilt2(K3,[hoodhood]);subplot(2,2,4);imshow(K4);title('迭代4次');I1=double(I);B1=double(K4);s1=0;s2=0;fori=1:1:336;forj=1:1:448;a=[I1(i,j)-B1(i,j)].^2;b=I1(i,j).^2;s1=s1+a;s2=s2+b;endends=s1/s2;disp(s');5维纳滤波I=imread('01.JPG');I=rgb2gray(I);figure;J=imnoise(I,'gaussian',0,0.01);K=wiener2(J,[33]);subplot(1,2,1);imshow(J);title('高斯噪声');subplot(1,2,2);imshow(K);title('高斯噪声的滤波效果');I1=double(I);M1=double(K);s1=0;s2=0;fori=1:1:336;forj=1:1:448;a=[I1(i,j)-M1(i,j)].^2;b=I1(i,j).^2;s1=s