图像噪声抑制

多媒体技术教程
实验报告
实验名称：图像噪声抑制
系（科）：信息系电子信息工程
班级：0932
学号：
姓名：
完成时间：2012年5月4日
南京师范大学中北学院
实验四图像噪声抑制
一、实验目的：
1、熟悉图像噪声的类型。

2、掌握空间域当中图像噪声抑制的基本方法，如均值滤波、中值滤波、边界保持类平滑滤波的算法。

二、实验原理：
1、图像的噪声：
所谓的图像噪声，是图像在摄取时或是传输时所受到的随机干扰信号。

这些干扰信号的抑制称为图像的噪声抑制。

对噪声进行分类，从统计数学的观点来看可将噪声分为平稳噪声和非平稳噪声。

平稳噪声：统计特性不随时间变化的噪声；非平稳噪声：统计特性随时间变化的噪声。

两种典型的图像噪声：其一为噪声的幅值基本相同，但出现的位置是随机分布的椒盐噪声，其二幅值分布是随机的，但每一点都存在噪声如高斯噪声、瑞利噪声、白噪声（频谱均匀分布）等。

对噪声的描述用统计意义上的均值和方差，噪声均值（图像中噪声的总体强度），噪声方差（图像中噪声分布的强弱）
∑∑==⨯==M x N y y x n N M y x n E n 11)
,(1)],([∑∑==-⨯=-=M x N Y n n y x n N M n y x n E 112
22),([1}]),({[σ按照噪声的模型按照对信号的影响：加性噪声、乘性噪声。

在MATLAB 中提供了给图像加入噪声的函数imnoise ，imnoise 的语法格式为J =imnoise(I,type)
J =imnoise(I,type,parameters)
其中J =imnoise(I,type)返回对原始图像I 添加典型噪声的有噪图像J 。

参数type 和parameters 用于确定噪声的类型和相应的参数。

下面的命令是对图像eight.tif 分别加入高斯噪声、椒盐噪声和乘性噪声，其结果如图所示：
例：
I=imread('eight.tif');
J1=imnoise(I,'gaussian',0,0.02);
J2=imnoise(I,'salt&pepper',0.02);
J3=imnoise(I,'speckle',0.02);
subplot(2,2,1),imshow(I),title('原图像');
subplot(2,2,2),imshow(J1),title('加高斯噪声');
subplot(2,2,3),imshow(J2),title('加椒盐噪声');
subplot(2,2,4),imshow(J3),title('加乘性噪声');
在上面的例子中使用了一个函数subplot。

其作用就是将多幅图像显示再同一幅图像显示对话框中。

其语法格式为：subplot(m,n,p)
其作用就是将一个图像显示对话框分成m行n列，并显示第p幅图像。

在MATLAB程序语言中，分号的用处为不显示程序运算中的中间结果，这在一定程度上使系统运算的效率增高，因此在不需知道中间结果的情况下，可以用分号作为一个句子的结尾，而不显示该句运算的中间结果。

2、滤波的算法
a)均值滤波：
均值滤波的原理是：在图像上，对待处理的像素给定一个模板，该模板包括了其周围的邻近像素。

将模板中的全体像素的均值来替代原来的像素值的方法。

均值滤波可以用来对椒盐噪声和高斯噪声进行滤波。

以模块运算系数表示即：⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡=111111111910H 将以上的均值滤波器加以修正，可以得到加权平均滤波器。

⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡=1111211111011H ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=111101111813H ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=1212421211612H ，⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡=0010041414141214H 。

b)中值滤波：均值滤波往往使得边界模糊，中值滤波对模板中的像素值由小到大排列，最终待处理像素的灰度取这个模板中排在中间位置上的像素的灰度值。

如做3*3的模板，对9个数排序，取第5个数替代原来的像素值。

实现方法见附录示例2：
c)边界保持类平滑的算法：经过平滑（特别是均值）滤波处理之后，图像就会变得模糊。

分析原因，在图像上的景物之所以可以辨认清楚是因为目标物之间存在边界。

而边界点与噪声点有一个共同的特点是，都具有灰度的跃变特性。

所以平滑处理会同时将边界也处理了。

为了解决图像模糊问题，一个自然的想法就是，在进行平滑处理时，首先判别当前像素是否为边界上的点，如果是，则不进行平滑处理；如果不是，则进行平滑处理。

对应算法有灰度最小方差的均值滤波器、K 近邻(KNN)平滑滤波器、对称近邻均值滤波器、西戈玛平滑滤波器等。

示例采用灰度最小方差的均值滤波器算法。

Matlab 中给出了相应函数，如conv2,imfilter 等，大家可以参考help 文件。

三、实验要求：
1、掌握以上噪声抑制的基本算法。

2、读出图像rice.tif ，分别加椒盐噪声和高斯噪声，得到两幅带有噪声的图像。

3、分别两幅图像进行均值滤波（使用⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡=1212421211612H ）和中值滤波，比较实验结果。

对加高斯噪声的图像进行均值滤波得到滤波前后的图像比较：
滤波前
滤波后
对加椒盐噪声的图像进行均值滤波得到滤波前后的图像比较：滤波前滤波后
对加高斯噪声的图像进行中值滤波的代码如下：
对加高斯噪声的图像进行中值滤波得到滤波前后的图像比较：
对加椒盐噪声的图像进行中值滤波的代码如下：
对加椒盐噪声的图像进行中值滤波得到滤波前后的图像比较：
4、如果进行K近邻(KNN)平滑滤波，请编码实现。

附录：
示例1
b=imread('pout.tif');
a1=imnoise(b,'gaussian',0.05);%加噪声
a=double(a1);
[m n]=size(a);
newimage=a;
for x=2:m-1
for y=2:n-1
s=0;
xin=[a(x-1,y-1),a(x,y-1),a(x+1,y-1),a(x-1,y),a(x,y),a(x+1,y),a(x-1,y+1),a(x,y+1),a(x+1,y+1)];
for z=1:9
s=s+xin(z);
end
newimage(x,y)=s/9;
end
end
imshow(a1)
figure,imshow(uint8(newimage))
示例2%用3*3模板
clear
f1=imread('pout.tif');
c=imnoise(f1,'salt&pepper',0.05);%加噪声
f=double(c);
[m n]=size(f);
f2=f;
f3=f;
subplot(2,2,1),imshow('pout.tif')
subplot(2,2,2),imshow(c)
for x=2:m-1
for y=2:n-1
a=0;%待排序矩阵a=[f(x-1,y-1),f(x-1,y),f(x-1,y+1),f(x,y-1),f(x,y),f(x,y+1),f(x+1,y-1),f(x+1,y),f(x+1,y+1)];
%利用中值滤波的算法，下面要排序找出第5个值，赋给当前像素b=0;
c=0;
for z=1:8%以下是排序算法
b=a(z);
c=z;
for zz=z+1:9
if a(zz)<b
b=a(zz);
c=zz;
end
end
a(c)=a(z);
a(z)=b;
end
f3(x,y)=a(5);%取出中值
end
end
subplot(2,2,4),imshow(uint8(f3))
示例3：灰度最小方差的均值滤波器，涉及模板、选最小值等操作
clear
b=imread('pout.tif');
a1=imnoise(b,'gaussian',0.05);
a=double(a1);
bbb=a;
[m n]=size(a);
for x=3:m-2
for y=3:n-2
s1=0;
s2=0;
s3=0;
s4=0;
s5=0;
s6=0;
s7=0;
s8=0;
s9=0;
%以下是四个屋形模板
xing1=[a(x-2,y-1),a(x-1,y-1),a(x-2,y),a(x-1,y),a(x-2,y+1),a(x-1,y+1),a(x,y)];
pmean(1)=mean(xing1);
for z=1:7
s1=s1+(xing1(z)-pmean(1))^2;%计算方差
end
fang(1)=s1/7;
xing2=[a(x-1,y+2),a(x-1,y+1),a(x,y+1),a(x,y+2),a(x,y+1),a(x,y+2),a(x,y)];
pmean(2)=mean(xing2);
for z=1:7
s2=s2+(xing2(z)-pmean(2))^2;
end
fang(2)=s2/7;
xing3=[a(x+2,y-1),a(x+2,y),a(x+2,y+1),a(x+1,y-1),a(x+1,y),a(x+1,y+1),a(x,y)];
pmean(3)=mean(xing3);
for z=1:7
s3=s1+(xing3(z)-pmean(3))^2;
end
fang(3)=s3/7;
xing4=[a(x-1,y-2),a(x-1,y-1),a(x,y-2),a(x,y-1),a(x,y),a(x+1,y-2),a(x,y-1)];
pmean(4)=mean(xing4);
for z=1:7
s4=s4+(xing4(z)-pmean(4))^2;
end
fang(4)=s4/7;
%以下是四个六边形模板
xing5=[a(x-2,y-2),a(x-1,y-2),a(x-2,y-1),a(x-1,y-1),a(x,y-1),a(x-1,y),a(x,y)];
pmean(5)=mean(xing5);
for z=1:7
s5=s4+(xing5(z)-pmean(5))^2;
end
fang(5)=s5/7;
xing6=[a(x+2,y+2),a(x+2,y+1),a(x+1,y+2),a(x+1,y+1),a(x+1,y),a(x,y+1),a(x,y)];
pmean(6)=mean(xing6);
for z=1:7
s6=s6+(xing6(z)-pmean(6))^2;
end
fang(6)=s6/7;
xing7=[a(x+2,y-2),a(x+2,y-1),a(x+1,y-2),a(x+1,y-1),a(x+1,y),a(x,y-1),a(x,y)];
pmean(7)=mean(xing7);
for z=1:7
s7=s7+(xing7(z)-pmean(7))^2;
end
fang(7)=s7/7;
xing8=[a(x-2,y-2),a(x-2,y-1),a(x-1,y-2),a(x-1,y-1),a(x-1,y),a(x,y-1),a(x,y)];
pmean(8)=mean(xing8);
for z=1:7
s8=s8+(xing8(z)-pmean(8))^2;
end
fang(8)=s8/7;
%以下是一个矩形模板
xing9=[a(x-1,y-1),a(x-1,y),a(x-1,y+1),a(x,y-1),a(x,y),a(x,y+1),a(x+1,y-1),a(x+1,y),a(x+1,y+1)];
pmean(9)=mean(xing9);
for z=1:9
s9=s9+(xing9(z)-pmean(9))^2;
end
fang(9)=s9/9;
minfang=min(fang);
for z=1:9
if fang(z)==minfang%得到最小的方差对应的模板
bbb(x,y)=round(pmean(z));%将最小方差对应的模板的均值赋给处理的像素end
end
end
end
subplot(1,3,1),imshow(b),title('原始图像')
subplot(1,3,2),imshow(a1),title('噪声图像')
subplot(1,3,3),imshow(uint8(bbb)),title('最小方差均值滤波图像')。