4、加载图像cameraman.tif，图像分别加入gaussian, salt&pepper, speckle

噪声，比较其对均值和方差的影响。

5、分别用函数medfilt2, ordfilt2 和 wiener2 去除所加的噪声。将处理后

的图像与原图像比较，讨论各种方法的作用。

6、设计陷波滤波器（Notch filter）

利用以下公式编程实现陷波滤波器。

7、研究空间滤波方法和噪声特征间的关系，利用 spfilt 函数举例说明所得结论。

8、图像的空间变换

用测试图像（C = checkerboard）产生一个指定形状的变形图像，

学习使用函数：maketform，imtransform，产生类似下列各图。

 四、实验结果以及代码（分段）

a=imread('C:\PICTURE_YUWEI\Fig0704(Vase).tif.tif');
b=imread('C:\PICTURE_YUWEI\camera.tif');

%使用imadjust函数调整图像视觉效果

%使用imadjust函数调整图像视觉效果
a1=imadjust(a,[0 1],[0 1],0.5);%gamma=0.5<1时，图象变亮
a2=imadjust(a,[0 1],[0 1],1.5);%gamma=1.5>1时，图象变暗
a3=imadjust(a,[0 1],[1 0]);%负片
figure(1);
subplot(2,2,1),imshow(a);title('图1.1 原图');
subplot(2,2,2),imshow(a1);title('图1.2 gamma=0.5<1时，图象变亮');
subplot(2,2,3),imshow(a2);title('图1.3 gamma=1.5>1时，图象变暗');
subplot(2,2,4),imshow(a3);title('图1.4 负片');

%分别用 imnoise()，imnoise2()，imnoise3()函数生成含有噪声的图像

%分别用 imnoise()，imnoise2()，imnoise3()函数生成含有噪声的图像
b1=imnoise(a,'gaussian',0,0.01); %在图像中加入均值为0、方差为0.1的高斯白噪声
n1=imnoise2('gaussian',512,512,0,0.1);%产生均值为0、方差为0.1的白噪声
b2=im2double(a)+n1; %将原图转换为双精度，并且与噪声n1进行叠加
C=[6 32;-2 2];%冲击位置的坐标
A=[1 5];%1*k维向量包含振幅的冲击对
[n2 R S]=imnoise3(512,512,C,A);%产生一个大小为M*N的正弦噪声模型n2；R为傅里叶变换，S为正弦噪声模型的傅里叶的频谱
n_gray=mat2gray(n2);%将得到的数据换归化成gray数据
b3=n_gray+im2double(a);%将原图转换为双精度再与噪声进行叠加
figure(2);
subplot(2,2,1),imshow(a);title('图2.1 原图');
subplot(2,2,2),imshow(a);title('图2.2 加imnoise噪声后的图像');
subplot(2,2,3),imshow(b2);title('图2.3 加imnoise2噪声后的图像');
subplot(2,2,4),imshow(b3,[]);title('图2.4 加imnoise3周期噪声后的图像');

%噪声估计

（这一部分没有搞出实验结果，总是报错，I'm sooo blue。遂重新写了一个程序文件，把这部分代码单独放在该文件中。在本文里，把有关噪声估计的代码放在最后。 ）

%加载图像cameraman.tif，图像分别加入gaussian, salt&pepper, speckle噪声

%加载图像cameraman.tif，图像分别加入gaussian, salt&pepper, speckle噪声
nn1=imnoise(b,'gaussian',0,0.09);%加入均值为0，方差为0.09的高斯噪声
nn2=imnoise(b,'salt & pepper',0.15);%加入密度为0.1的椒盐噪声
nn3=imnoise(b,'speckle',0.09);%加入方差为0.09的乘性噪声系数的斑点噪声
figure(3);
subplot(2,2,1);imshow(b);title('图3.1 原图');
subplot(2,2,2);imshow(nn1);title('图3.2 加高斯噪声的图');
subplot(2,2,3);imshow(nn2);title('图3.3 加椒盐噪声的图');
subplot(2,2,4);imshow(nn3);title('图3.4 加斑点噪声的图');

%分别用函数medfilt2,ordfilt2和wiener2去除所加的噪声。将处理后的图像与原图像比较。

%分别用函数medfilt2,ordfilt2和wiener2去除所加的噪声。将处理后的图像与原图像比较。
%1. 中值滤波
f11=medfilt2(b,[4 4],'symmetric');%使用4*4大小的中值滤波器进行滤波
f12=medfilt2(b,[8 8],'symmetric');%8*8
figure(4);
subplot(131);imshow(g);title('图5.1 噪声图');
subplot(132);imshow(f11);title('图5.2 4*4滤波器滤过的图');
subplot(133);imshow(f12);title('图5.3 8*8滤波器滤过的图');
%2. 顺序滤波
f21=ordfilt2(g,4,ones(4,4));%即在4*4的矩阵中将值按照由大到小排序，选取第五个值
f22=ordfilt2(g,9,ones(4,4));%第九个值
figure(5);
subplot(131);imshow(g);title('图6.1 噪声原图 ');
subplot(132);imshow(f21);title('图6.2 3*3的滤波器第四大');
subplot(133);imshow(f22);title('图6.3 3*3的滤波器第九大');
%3. 维纳滤波
f31=wiener2(g,[4 4]);%4*4滤波器滤波
f33=wiener2(g,[15 15]);%15*15滤波器
figure(6);
subplot(131);imshow(g);title('图7.1 原图');
subplot(132);imshow(f31);title('图7.2 4*4维纳滤波器滤过的图');
subplot(133);imshow(f33);title('图7.3 15*15维纳滤波器滤过的图');

设计陷波滤波器（Notch filter）

%设计陷波滤波器（Notch filter）
M=60;
N=60;%大小为60*60的滤波器
n=3;%滤波器的尖峰程度/衰减速度
D0=5;%滤波器截止频率
H=ones(60,60);%生成一个60*60矩阵
for u=1:60for v=1:60
H1(u,v)=1/(1+(D0^2/(sqrt((u-M/2-10)^2+(v-N/2-10)^2).*sqrt((u-M/2+10)^2+(v-N/2+10)^2)))^n);               end
end %在给定的范围内生成两个中心为（10,10）（-10，-10）的高通滤波器
for u=1:60for v=1:60        H2(u,v)=1/(1+(D0^2/(sqrt((u-M/2+10)^2+(v-N/2-10)^2).*sqrt((u-M/2-10)^2+(v-N/2+10)^2)))^n);    end
end%生成两个中心为（-10,10）（10，-10）的高通滤波器
H=H1.*H2;%将两个滤波器合成为一个有四个中心的高通滤波器J=ones(60,60);%生成一个矩阵
for u=1:60for v=1:60
J1(u,v)=1/(1+(D0^2/(sqrt((u-M/2-25)^2+(v-N/2-25)^2).*sqrt((u-M/2-25)^2+(v-N/2-10)^2)))^n);    end%生成两个中心为（25,25）（25，10）的高通滤波器
end
for u=1:60for v=1:60
J2(u,v)=1/(1+(D0^2/(sqrt((u-M/2-0)^2+(v-N/2-10)^2).*sqrt((u-M/2-0)^2+(v-N/2-25)^2)))^n);    end%合成
end
J=J1.*J2;%将生成的两个滤波器合成为一个带有四个高通滤波器的滤波器
figure(7);
subplot(1,2,1),mesh(H);title('图8.1 3*3的滤波器，3');
subplot(1,2,2),mesh(J);title('图8.2 3*3的滤波器，3');

%研究空间滤波方法和噪声特征间的关系，利用spfilt函数举例说明所得结论

%研究空间滤波方法和噪声特征间的关系，利用spfilt函数举例说明所得结论
addpath('C:\PICTURE_YUWEI\dipum_toolbox_2.0.2');%添加.p文件
g1=imnoise(b,'gaussian',0,0.09);%加入均值为0，方差为0.09的高斯噪声
figure(8);
subplot(1,2,1);imshow(g1);title('图9.1 被高斯噪声污染的图像');
g2=spfilt(g1,'atrimmed',3,3,2);%滤波
subplot(1,2,2);imshow(g2);title('图9.2 滤波后的图像');

图像的空间变换

%图像的空间变换
C=checkerboard(64,4,4);%在生成的图像中，每个正方形的大小是64*64
T1=[1 0.8 0;0 1 0;0 0 1];%x'=x+0.8y;y'=y
tform1=maketform('affine',T1);%根据T1创建tform，表示空间变化
C1=imtransform(C,tform1);%将原图进行反向映射
T2=[1 0 0;0.8 1 0;0 0 1];%x'=x;y'=y+0.8x
tform1=maketform('affine',T2);
C2=imtransform(C,tform1);
figure(9);
subplot(131);imshow(C);title('图10.1');
subplot(132);imshow(C2);title('图10.2');
subplot(133);imshow(C1);title('图10.3');

五、conclusions

（1）在使用imadjust函数调整视觉效果时，当gamma>1时，图像变亮；gamma<1时，图像变暗；gamma=1是，图像成为负片。

（2）学会了使用imnoise,imnoise2,imnoise3进行噪声的生成。使用imnoise的时候，是直接在图像中加入一定性质的噪声，而使用imnoise2时，是直接产生一个特定性质的噪声图像，他是双精度的，需要将原图转换为双精度后与之叠加，才能得到加入噪声的图像；使用imnoise3和imnoise2相似，不同的是，imnoise还需要将噪声图像转换为灰度图像，再与转换为双精度的原图相加才可以得到加了噪声后的图像。

（3）学会了使用函数生成特定性质的高斯、椒盐、斑点噪声。并学会使用中值medfilt2滤波、顺序ordfilt2滤波和维纳wiener2滤波将图像中的噪声去除。

（4）学会了如何设计陷波滤波器。需要定义它的像素大小、尖峰程度、截止频率，需要先生成两组包含两个有特定中心的高通滤波器，再将这两组滤波器合成一个包含四个高通滤波器的陷波滤波器。

（5）学会如何使用checkerboard函数生成棋盘格图像，使用maketform函数对其进行空间变换，使用imtransform函数对其进行反向映射。

（6）学会噪声估计的基本思路。即，首先生成一个特定的噪声图像，使用roipoly函数来选取一个ROI范围。分别画出这个范围的直方图，和噪声图像的直方图，并进行比较。可惜的是，matlab提示我输入参数太多，没有成功地完成噪声估计。

六、实验代码（全部）

%第一部分：除了噪声估计之外的其他部分代码。
a=imread('C:\PICTURE_YUWEI\Fig0704(Vase).tif.tif');
b=imread('C:\PICTURE_YUWEI\camera.tif');%使用imadjust函数调整图像视觉效果
a1=imadjust(a,[0 1],[0 1],0.5);%gamma=0.5<1时，图象变亮
a2=imadjust(a,[0 1],[0 1],1.5);%gamma=1.5>1时，图象变暗
a3=imadjust(a,[0 1],[1 0]);%负片
figure(1);
subplot(2,2,1),imshow(a);title('图1.1 原图');
subplot(2,2,2),imshow(a1);title('图1.2 gamma=0.5<1时，图象变亮');
subplot(2,2,3),imshow(a2);title('图1.3 gamma=1.5>1时，图象变暗');
subplot(2,2,4),imshow(a3);title('图1.4 负片');
%分别用 imnoise()，imnoise2()，imnoise3()函数生成含有噪声的图像
b1=imnoise(a,'gaussian',0,0.01); %在图像中加入均值为0、方差为0.1的高斯白噪声
n1=imnoise2('gaussian',512,512,0,0.1);%产生均值为0、方差为0.1的白噪声
b2=im2double(a)+n1; %将原图转换为双精度，并且与噪声n1进行叠加
C=[6 32;-2 2];%冲击位置的坐标
A=[1 5];%1*k维向量包含振幅的冲击对
[n2 R S]=imnoise3(512,512,C,A);%产生一个大小为M*N的正弦噪声模型n2；R为傅里叶变换，S为正弦噪声模型的傅里叶的频谱
n_gray=mat2gray(n2);%将得到的数据换归化成gray数据
b3=n_gray+im2double(a);%将原图转换为双精度再与噪声进行叠加
figure(2);
subplot(2,2,1),imshow(a);title('图2.1 原图');
subplot(2,2,2),imshow(a);title('图2.2 加imnoise噪声后的图像');
subplot(2,2,3),imshow(b2);title('图2.3 加imnoise2噪声后的图像');
subplot(2,2,4),imshow(b3,[]);title('图2.4 加imnoise3周期噪声后的图像');%加载图像cameraman.tif，图像分别加入gaussian, salt&pepper, speckle噪声
nn1=imnoise(b,'gaussian',0,0.09);%加入均值为0，方差为0.09的高斯噪声
nn2=imnoise(b,'salt & pepper',0.15);%加入密度为0.1的椒盐噪声
nn3=imnoise(b,'speckle',0.09);%加入方差为0.09的乘性噪声系数的斑点噪声
figure(3);
subplot(2,2,1);imshow(b);title('图3.1 原图');
subplot(2,2,2);imshow(nn1);title('图3.2 加高斯噪声的图');
subplot(2,2,3);imshow(nn2);title('图3.3 加椒盐噪声的图');
subplot(2,2,4);imshow(nn3);title('图3.4 加斑点噪声的图');%分别用函数medfilt2,ordfilt2和wiener2去除所加的噪声。将处理后的图像与原图像比较。
%1. 中值滤波
f11=medfilt2(b,[4 4],'symmetric');%使用4*4大小的中值滤波器进行滤波
f12=medfilt2(b,[8 8],'symmetric');%8*8
figure(4);
subplot(131);imshow(g);title('图5.1 噪声图');
subplot(132);imshow(f11);title('图5.2 4*4滤波器滤过的图');
subplot(133);imshow(f12);title('图5.3 8*8滤波器滤过的图');
%2. 顺序滤波
f21=ordfilt2(g,4,ones(4,4));%即在4*4的矩阵中将值按照由大到小排序，选取第五个值
f22=ordfilt2(g,9,ones(4,4));%第九个值
figure(5);
subplot(131);imshow(g);title('图6.1 噪声原图 ');
subplot(132);imshow(f21);title('图6.2 3*3的滤波器第四大');
subplot(133);imshow(f22);title('图6.3 3*3的滤波器第九大');
%3. 维纳滤波
f31=wiener2(g,[4 4]);%4*4滤波器滤波
f33=wiener2(g,[15 15]);%15*15滤波器
figure(6);
subplot(131);imshow(g);title('图7.1 原图');
subplot(132);imshow(f31);title('图7.2 4*4维纳滤波器滤过的图');
subplot(133);imshow(f33);title('图7.3 15*15维纳滤波器滤过的图');%设计陷波滤波器（Notch filter）
M=60;
N=60;%大小为60*60的滤波器
n=3;%滤波器的尖峰程度/衰减速度
D0=5;%滤波器截止频率
H=ones(60,60);%生成一个60*60矩阵
for u=1:60for v=1:60
H1(u,v)=1/(1+(D0^2/(sqrt((u-M/2-10)^2+(v-N/2-10)^2).*sqrt((u-M/2+10)^2+(v-N/2+10)^2)))^n);               end
end %在给定的范围内生成两个中心为（10,10）（-10，-10）的高通滤波器
for u=1:60for v=1:60        H2(u,v)=1/(1+(D0^2/(sqrt((u-M/2+10)^2+(v-N/2-10)^2).*sqrt((u-M/2-10)^2+(v-N/2+10)^2)))^n);    end
end%生成两个中心为（-10,10）（10，-10）的高通滤波器
H=H1.*H2;%将两个滤波器合成为一个有四个中心的高通滤波器J=ones(60,60);%生成一个矩阵
for u=1:60for v=1:60
J1(u,v)=1/(1+(D0^2/(sqrt((u-M/2-25)^2+(v-N/2-25)^2).*sqrt((u-M/2-25)^2+(v-N/2-10)^2)))^n);    end%生成两个中心为（25,25）（25，10）的高通滤波器
end
for u=1:60for v=1:60
J2(u,v)=1/(1+(D0^2/(sqrt((u-M/2-0)^2+(v-N/2-10)^2).*sqrt((u-M/2-0)^2+(v-N/2-25)^2)))^n);    end%合成
end
J=J1.*J2;%将生成的两个滤波器合成为一个带有四个高通滤波器的滤波器
figure(7);
subplot(1,2,1),mesh(H);title('图8.1 3*3的滤波器，3');
subplot(1,2,2),mesh(J);title('图8.2 3*3的滤波器，3');%研究空间滤波方法和噪声特征间的关系，利用spfilt函数举例说明所得结论
addpath('C:\PICTURE_YUWEI\dipum_toolbox_2.0.2');%添加.p文件
g1=imnoise(b,'gaussian',0,0.09);%加入均值为0，方差为0.09的高斯噪声
figure(8);
subplot(1,2,1);imshow(g1);title('图9.1 被高斯噪声污染的图像');
g2=spfilt(g1,'atrimmed',3,3,2);%滤波
subplot(1,2,2);imshow(g2);title('图9.2 滤波后的图像');%图像的空间变换
C=checkerboard(64,4,4);%在生成的图像中，每个正方形的大小是64*64
T1=[1 0.8 0;0 1 0;0 0 1];%x'=x+0.8y;y'=y
tform1=maketform('affine',T1);%根据T1创建tform，表示空间变化
C1=imtransform(C,tform1);%将原图进行反向映射
T2=[1 0 0;0.8 1 0;0 0 1];%x'=x;y'=y+0.8x
tform1=maketform('affine',T2);
C2=imtransform(C,tform1);
figure(9);
subplot(131);imshow(C);title('图10.1');
subplot(132);imshow(C2);title('图10.2');
subplot(133);imshow(C1);title('图10.3');

%第二部分：噪声估计部分的代码
c=imread('C:\PICTURE_YUWEI\Fig0704(Vase).tif.tif');
g=imnoise(c,'gaussian',0,0.1);%添加均值0 方差0.1的高斯噪声图像
B=roipoly(g);%在噪声图像中交互生成的ROI npix=sum(B(:));%求ROI区域的点数
figure(1);
subplot(2,2,1);
imshow(g);%显示被高斯噪声污染的图像
subplot(2,2,2);
imshow(B);%显示ROI区域 subplot(2,2,3);
h=imhist(c(B));%画出ROI区域的直方图
bar(h,1);
axis([0,256,0,1000]);
title('ROI区域的直方图');
[v,unv]=statmoments(h,2);%计算均值与方差；
subplot(2,2,4);
X=imnoise2('gaussian',npix,1,unv(1,1),sqrt(unv(1,2)));%生成与图片中有相同均值方差的高斯噪声
hist(X,140);%输出高斯噪声的直方图
axis([0,256,0,1000]);