西安电子科技大学
数字图像处理上机实验题
李 XX 2018-12-16
数字图像处理上机实验题
1. 产生右图所示图像 f1(m,n)，其中图像大小为256×256，中间亮条为 128×32，暗处=0，亮处=100。对其进行 FFT：(matlab 程序附于文末附录)
① 同屏显示原图 f1(m,n)和 FFT(f1)的幅度谱图；
f1的图像
FFT(f1) 的幅度频谱
② 若令 f2(m,n)=(-1)m+n f1(m,n)，重复以上过程，比较二者幅度谱的异同，简述理由；
结论：
相同点：图像 f1,f2 的幅度谱的实质没有改变，幅度等都没有发生变化。
不同点：f2 的频谱是对 f1 频谱的移位，它是 f1 的频谱
从原点(0,0)移动到中心点
(,)得到的频谱。
原因：f2 图像是对 f1 图像做频谱中心化变换得来的，这一过程不改变图像的幅度等特性。
③ 若将 f2(m,n)顺时针旋转 90 度得到 f3(m,n)，试显示 FFT(f3)的幅度谱，并与 FFT(f2)的幅度谱进行比较；
FFT(f3) 的幅度频谱 FFT(f2) 的幅度频谱
结论：
由图可知，FFT(f3)的幅度谱是由 FFT(f2)的幅度谱
顺时针旋转 90 度得到的。
④ 若将f1(m,n) 顺时针旋转90 度得到f4(m,n)，令 FFT(f5)的图像 f5(m,n)=f1(m,n)+f4(m,n)，试显示 FFT(f5)的幅度谱，并指出其与 FFT(f1)和 FFT(f4)的关系；
f2的图像 FFT(f2) 的幅度频谱
结论：
FFT(f5)的幅度谱是 FFT(f1)和 FFT(f4)幅度谱
1
叠加的结果。
⑤ 若令 f6(m,n)=f2(m,n)+f3(m,n)， 试显示 FFT(f6)的幅度谱， 并指出其与FFT(f2)和 FFT(f3)
的关系，比较 FFT(f6)和 FFT(f5)的幅度谱。
FFT(f6)的图像
结论：
FFT(f6)的幅度谱是 FFT(f2)和 FFT(f3)幅度谱叠加的结果。FFT(f6)是 FFT(f5)频谱中心化的结果。
2. 产生教材 104 页题图 4.18（右图）所示的二值图像（白为 1，黑为 0），编程实现习题 4.18 所要求的处理（3×3 的平均滤波和中值滤波）功能（图像四周边界不考虑，处理结果按四舍五入仍取 0 或 1），显示处理前后的图像，比较其异同。 (matlab 程序附于文末附录)
fa的图像 fb的图像
fa的平均滤波图像
fb的平均滤波图像
结论：
图像 fa 的平均滤波图像和平
fa的加权平均滤波图像 fb的加权平均滤波图像
均加权图像与原图像相同；fb 的平均滤波图像的处理点在白、黑块组成的4块交界处时，由于窗内1值点等于0 值点个数，0 值点的值变为1，形成如图所示图像， 加权平均滤波图像与原图像相
同。
3. 产生教材 104 页题图 4.16 所示的灰度图像（白为 255，黑为 0），分别加入高斯白噪声和椒盐噪声，再分别进行 3×3 的平均滤波和中值滤波，显示原图像、加噪图像和滤波结果图像，并比较四种滤波结果。(matlab 程序附于文末附录)
2
原图像f
加高斯白噪声后的图像 平均滤波后的图像 中值滤波后的图像
加椒盐噪声后的图像
平均滤波后的图像
中值滤波后的图像
结论：
由图可知，对于加高斯白噪声后的图像，平均滤波后对噪声的滤除总体上效果好，但在局部边缘处存在误将噪声当作图像保留和误将图像细节当作噪声滤除的现象；中值滤波则对图像的边缘信息保留的较好，但总体滤除效果不佳。对于加椒盐噪声后的图像，平均滤波后在图像边缘上的损失较多，相对的，中值滤波后边缘信息损失较少。
4. 对某一灰度图像，进行如下处理：(matlab 程序附于文末附录) （1） 分别利用 Roberts、Prewitt 和 Sobel 边缘检测算子进行边缘检测；
本题选取此图片进行处理。
3
Roberts算子处理结果图
Prewitt算子处理结果图
Sobel算子处理结果图
（2） 将 Roberts、Prewitt 和 Sobel 边缘检测算子修改为锐化算子，对原图像进
行锐化，同屏显示原图像、边缘检测结果和锐化后图像，说明三者之间的关系。
原图像 Roberts算子边缘检测结果 锐化处理结果
原图像 Prewitt算子边缘检测结果 锐化处理结果
原图像 Sobel算子边缘检测结果 锐化处理结果
结论：
Roberts 边缘检测算子获得的边缘效果较差，把一些边缘也剔除掉了。
Prewitt 边缘检测算子获得的边缘要好于 Roberts 算子，但仍有部分边缘被 剔除。
Sobel 边缘检测算子相较于 Prewitt 算子主要是对噪声的抑制能力增强，由于此处原图像没有噪声，所以对比不明显。
4
每种算子的锐化结果是将其检测到的边缘相应加强的结果。
5、编程实现教材 214 页所给图像门限化分割的迭代阈值算法，实现对某一灰度图像的二值化。(matlab 程序附于文末附录)
原图像
迭代阈值处理后的二值图像
所有题目的 matlab 程序皆附于文末附录。
5
附录
题目 1 matlab 程序:
close all;clear; f1=zeros(256,256); f2=zeros(256,256); for m=:1:192 for n=112:1:144 f1(m,n)=100/255; f2(m,n)=(-1)^(m+n)*100/255; end end % figure % imshow(f1) f3=imrotate(f2,-90,'nearest'); f4=imrotate(f1,-90,'nearest'); f5=f1+f4; f6=f2+f3; %FFT 变换 fft_f1=log(1+abs(fft2(f1))); fft_f2=log(1+abs(fft2(f2))); fft_f3=log(1+abs(fft2(f3))); fft_f5=log(1+abs(fft2(f5))); fft_f6=log(1+abs(fft2(f6))); %显示图像 figure subplot(1,2,1); imshow(f1); title('f1 的图像'); subplot(1,2,2); imshow(fft_f1,[]); title('FFT(f1) 的幅度频谱'); figure subplot(1,2,1); imshow(f2); title('f2 的图像'); subplot(1,2,2); imshow(fft_f2,[]); title('FFT(f2) 的幅度频谱'); figure subplot(1,2,1); imshow(fft_f2,[]); title('FFT(f2) 的幅度频谱'); subplot(1,2,2); imshow(fft_f3,[]); title('FFT(f3) 的幅度频谱'); figure imshow(fft_f5,[]); title('FFT(f5)的图像'); figure imshow(fft_f6,[]); title('FFT(f6)的图像');
题目 2 matlab 程序： clear;close all; fb=zeros(,); fa=[ones(,32),zeros(,32)]; for m=1:8: for n=1:8: fb(m:m+7,n:n+7)=1/2*(1+(-1)^(floor(m/8)+floor(n/8)))*ones(8,8);
1
end end %3×3 平均模板 W1=1/8*[1 1 1;1 0 1;1 1 1];%平均滤波模板 fa_avefilter=fa;fb_avefilter=fb; for m=2:1:63 for n=2:1:63 fa_avefilter(m,n)=round(W1(1,1)*fa(m-1,n-1)+W1(1,2)*fa(m-1,n)+W1(1,3)*fa(m-1,n+1)... +W1(2,1)*fa(m,n-1)+W1(2,2)*fa(m,n)+W1(2,3)*fa(m,n+1)... +W1(3,1)*fa(m+1,n-1)+W1(3,2)*fa(m+1,n)+W1(3,3)*fa(m+1,n+1)); fb_avefilter(m,n)=round(W1(1,1)*fb(m-1,n-1)+W1(1,2)*fb(m-1,n)+W1(1,3)*fb(m-1,n+1)... +W1(2,1)*fb(m,n-1)+W1(2,2)*fb(m,n)+W1(2,3)*fb(m,n+1)... +W1(3,1)*fb(m+1,n-1)+W1(3,2)*fb(m+1,n)+W1(3,3)*fb(m+1,n+1)); end end W2=1/9*[1 1 1;1 1 1;1 1 1];%加权平均滤波模板 fa_whtavefilter=fa;fb_whtavefilter=fb; for m=2:1:63 for n=2:1:63 fa_whtavefilter(m,n)=round(W2(1,1)*fa(m-1,n-1)+W2(1,2)*fa(m-1,n)+W2(1,3)*fa(m- 1,n+1)... +W2(2,1)*fa(m,n-1)+W2(2,2)*fa(m,n)+W2(2,3)*fa(m,n+1)... +W2(3,1)*fa(m+1,n-1)+W2(3,2)*fa(m+1,n)+W2(3,3)*fa(m+1,n+1)); fb_whtavefilter(m,n)=round(W2(1,1)*fb(m-1,n-1)+W2(1,2)*fb(m-1,n)+W2(1,3)*fb(m- 1,n+1)... +W2(2,1)*fb(m,n-1)+W2(2,2)*fb(m,n)+W2(2,3)*fb(m,n+1)... +W2(3,1)*fb(m+1,n-1)+W2(3,2)*fb(m+1,n)+W2(3,3)*fb(m+1,n+1)); end end %中值滤波 fa_midfilter=fa; fb_midfilter=fb;
2
for m=2:1:63 for n=2:1:63 fa_midfilter(m,n)=mid([fa(m-1,n-1),fa(m-1,n),fa(m-1,n+1),fa(m,n- 1),fa(m,n),fa(m,n+1),fa(m+1,n-1),fa(m+1,n),fa(m+1,n+1)]); fb_midfilter(m,n)=mid([fb(m-1,n-1),fb(m-1,n),fb(m-1,n+1),fb(m,n- 1),fb(m,n),fb(m,n+1),fb(m+1,n-1),fb(m+1,n),fb(m+1,n+1)]); end end figure(1) subplot(1,2,1); imshow(fa); title('fa 的图像'); subplot(1,2,2); imshow(fb); title('fb 的图像'); figure(2) subplot(1,2,1); imshow(fa_avefilter); title('fa 的平均滤波图像'); subplot(1,2,2); imshow(fb_avefilter); title('fb 的平均滤波图像'); figure(3) subplot(1,2,1); imshow(fa_whtavefilter); title('fa 的加权平均滤波图像'); subplot(1,2,2); imshow(fb_whtavefilter); title('fb 的加权平均滤波图像'); figure(4) subplot(1,2,1); imshow(fa_midfilter);
3
title('fa 的中值滤波图像'); subplot(1,2,2); imshow(fb_midfilter); title('fb 的中值滤波图像'); function MID=mid(F) N=length(F); %冒泡排序for i=1:1:N for j=i:1:N if(F(i)>F(j)) tmp=F(i);F(i)=F(j);F(j)=tmp; end end end MID=F(round(N/2)); end
问题 3 matlab 程序：
close all;clear; f=zeros(256,256); for m=28:24:220 f(23:233,m:m+7)=1; end %高斯白噪声guass_noise=0.2*randn(size(f)); f_gn=f+guass_noise; %椒盐噪声f_pep=f; k1=0.1; k2=0.3; a1=rand(size(f))4a2=rand(size(f))5figure(2) subplot(1,3,1) imshow(f_gn); title('加高斯白噪声后的图像'); subplot(1,3,2) imshow(f_gn_avefilter); title('平均滤波后的图像'); subplot(1,3,3) imshow(f_gn_midfilter); title('中值滤波后的图像'); figure(3) subplot(1,3,1) imshow(f_pep); title('加椒盐噪声后的图像'); subplot(1,3,2) imshow(f_pep_avefilter); title('平均滤波后的图像'); subplot(1,3,3) imshow(f_pep_midfilter); title('中值滤波后的图像');
问题 4 matlab 程序：
close all;clear; % 读取图像fid=fopen('lena.img','r'); Img=(fread(fid,[256,256],'uint8'))'; f=Img/255; clear fid Img imshow(f); % Roberts 边缘检测算子
6
Gh=zeros(256,256); Gv=zeros(256,256); for m=2:1:255 for n=2:1:255 Gh(m,n)=f(m,n)-f(m-1,n-1); Gv(m,n)=f(m,n-1)-f(m-1,n); end end G_Roberts=sqrt(Gh.*Gh+Gv.*Gv); g_Roberts=f+Gh+Gv;% 用 Roberts 算子做锐化 clear Gh Gv B_Roberts=zeros(256,256); for m=2:1:255 for n=2:1:255 if G_Roberts(m,n)>=0.1 B_Roberts(m,n)=1; end end end clear m n figure imshow(B_Roberts); title('Roberts 算子处理结果图'); % Prewitt 梯度算子 Gh=zeros(256,256); Gv=zeros(256,256); for m=2:1:255 for n=2:1:255 Gh(m,n)=1/3*(f(m-1,n+1)+f(m,n+1)+f(m+1,n+1)-f(m-1,n-1)-f(m,n-1)-f(m+1,n-1)); Gv(m,n)=1/3*(f(m-1,n-1)+f(m-1,n)+f(m-1,n+1)-f(m+1,n-1)-f(m+1,n)-f(m+1,n+1)); end end
7
G_Prewitt=sqrt(Gh.*Gh+Gv.*Gv); g_Prewitt=f+Gh+Gv;% 用 Prewitt 算子做锐化 clear Gh Gv B_Prewitt=zeros(256,256); for m=2:1:255 for n=2:1:255 if G_Prewitt(m,n)>=0.1 B_Prewitt(m,n)=1; end end end clear m n figure imshow(B_Prewitt); title('Prewitt 算子处理结果图'); % Sobel 梯度算子 Gh=zeros(256,256); Gv=zeros(256,256); for m=2:1:255 for n=2:1:255 Gh(m,n)=1/4*(f(m-1,n+1)+2*f(m,n+1)+f(m+1,n+1)-f(m-1,n-1)-2*f(m,n-1)-f(m+1,n-1)); Gv(m,n)=1/4*(f(m-1,n-1)+2*f(m-1,n)+f(m-1,n+1)-f(m+1,n-1)-2*f(m+1,n)-f(m+1,n+1)); end end G_Sobel=sqrt(Gh.*Gh+Gv.*Gv); g_Sobel=f+Gh+Gv;% 用 Sobel 算子做锐化 clear Gh Gv B_Sobel=zeros(256,256); for m=2:1:255 for n=2:1:255 if G_Sobel(m,n)>=0.1 B_Sobel(m,n)=1;
8
end end end clear m n figure imshow(B_Sobel); title('Sobel 算子处理结果图'); figure title('Roberts 算子'); subplot(1,3,1) imshow(f); title('原图像'); subplot(1,3,2) imshow(B_Roberts); title('Roberts 算子边缘检测结果'); subplot(1,3,3) imshow(g_Roberts); title('锐化处理结果'); figure title('Prewitt 算子'); subplot(1,3,1) imshow(f); title('原图像'); subplot(1,3,2) imshow(B_Prewitt); title('Prewitt 算子边缘检测结果'); subplot(1,3,3) imshow(g_Prewitt); title('锐化处理结果'); figure
9
title('Sobel 算子') subplot(1,3,1) imshow(f); title('原图像'); subplot(1,3,2) imshow(B_Sobel); title('Sobel 算子边缘检测结果'); subplot(1,3,3) imshow(g_Sobel); title('锐化处理结果');
问题 5 matlab 程序： close all;clear; % 读取图像 dim=256; fid=fopen('lena.img','r'); Img=(fread(fid,[dim,dim],'uint8'))'; f=Img/(dim-1); clear fid Img t1=max(max(f));% 最大灰度值点tk=min(min(f));% 最小灰度值点N=ones(dim,dim); T(1)=(t1+tk)/2;% 初始阈值clear t1 tk k=1; M=10;% 设定迭代次数 % 迭代阈值法 while 1 N1(k)=0;N2(k)=0; t_O=0;t_B=0; for m=1:1:dim for n=1:1:dim 10
if f(m,n)T(k) t_B=t_B+f(m,n)*N(m,n); N2(k)=N2(k)+1; end end end t_O=t_O/N1(k); t_B=t_B/N2(k); T(k+1)=(t_O+t_B)/2; if T(k)==T(k+1)||k>M break; end k=k+1; end clear N N1 N2 t_O t_B % 使用所得阈值处理图像 f1=f; for m=1:1:dim for n=1:1:dim if f1(m,n)11subplot(1,2,2) imshow(f1); title('迭代阈值处理后的二值图像');
12