数字图像处理实验报告.docx

数字图像处理试验报告实验二：数字图像的空间滤波和频域滤波姓名：XX学号：2XXXXXXX实验日期：2017 年 4 月 26 日1.实验目的1.掌握图像滤波的基本定义及目的。 2.理解空间域滤波的基本原理及方法。 3.掌握进行图像的空域滤波的方法。 4.掌握傅立叶变换及逆变换的基本原理方法。 5.理解频域滤波的基本原理及方法。 6.掌握进行图像的频域滤波的方法。 2.实验内容与要求1.平滑空间滤波： 1)读出一幅图像，给这幅图像分别加入椒盐噪声和高斯噪声后并与前一张图显示在同 一图像窗口中。 2)对加入噪声图像选用不同的平滑（低通）模板做运算，对比不同模板所形成的效果， 要求在同一窗口中显示。 3) 使 用 函 数 imfilter 时 ， 分 别 采 用 不 同 的 填 充 方 法 （ 或 边 界 选 项 ， 如 零 填 充、replicate、symmetric、circular）进行低通滤波，显示处理后 的图像。 4) 运用 for 循环，将加有椒盐噪声的图像进行 10 次，20 次均值滤波，查看其特点, 显示均值处理后的图像（提示:利用 fspecial 函数的average类型生成均值滤 波器）。 5)对加入椒盐噪声的图像分别采用均值滤波法，和中值滤波法对有噪声的图像做处理， 要求在同一窗口中显示结果。 6)自己设计平滑空间滤波器，并将其对噪声图像进行处理，显示处理后的图像。 2.锐化空间滤波 1)读出一幅图像，采用 33 的拉普拉斯算子 w = 1, 1, 1; 1 8 1; 1, 1, 1对其进行滤波。 2)编写函数 w = genlaplacian(n)，自动产生任一奇数尺寸 n 的拉普拉斯算子，如 55 的拉普拉斯算子 w = 1 1 1 1 1 1 1 1 1 1 1 1 -24 1 1 1 1 1 1 1 1 1 1 1 1 3)分别采用 55，99，1515 和 2525 大小的拉普拉斯算子对 blurry_moon.tif进 行锐化滤波，并利用式 g(x, y) =f (x, y) - 2 f (x, y) 完成图像的锐化增强，观察其有何不同，要求在同一窗口中显示。 4)采用不同的梯度算子对该幅图像进行锐化滤波，并比较其效果。 5)自己设计锐化空间滤波器，并将其对噪声图像进行处理，显示处理后的图像； 3.傅立叶变换 1)读出一幅图像，对其进行快速傅立叶变换，分别显示其幅度图像和相位图像。仅对 相位部分进行傅立叶反变换后查看结果图像。 2)仅对幅度部分进行傅立叶反变换后查看结果图像。 3)将图像的傅立叶变换 F 置为其共轭后进行反变换，比较新生成图像与原始图像的差 异。 4.平滑频域滤波 1)设计理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器，截至频率自选，分 别给出各种滤波器的透视图。 2) 读出一幅图像，分别采用理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器 对其进行滤波（截至频率自选），再做反变换，观察不同的截止频率下采用不同低通 滤波器得到的图像与原图像的区别，特别注意振铃效应。(提示:1)在频率域滤波同 样要注意到填充问题；2）注意到(-1)x+y;) 5.锐化频域滤波 1)设计理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器，截至频率自选，分 别给出各种滤波器的透视图。 2) 读出一幅图像，分别采用理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器 对其进行滤波（截至频率自选），再做反变换，观察不同的截止频率下采用不同高通 滤波器得到的图像与原图像的区别。 3.实验具体实现1.平滑空间滤波： （1）.读出一幅图像，给这幅图像分别加入椒盐噪声和高斯噪声后并与前一张图显示在 同一图像窗口中。 img=imread(lena.png) figure,subplot(1,3,1); imshow(img);title(原始图像); img2=imnoise(img,salt & pepper,0.02); subplot(1,3,2); imshow(img2); title(椒盐噪声图像); img3=imnoise(img,gaussian,0.02); subplot(1,3,3),imshow(img3); title(高斯噪声图像); 实验结果如下： （2）.对加入噪声图像选用不同的平滑（低通）模板做运算，对比不同模板所形成的效果， 要求在同一窗口中显示。 平滑滤波是低频增强的空间域滤波技术。它的目的有两个，一是模糊，二是消除噪声。将空 间域低通滤波按线性和非线性特点有：线性、非线性平滑滤波器，线性平滑滤波器包括均值 滤波器，非线性的平滑滤波器有最大值滤波器，中值滤波器，最小值滤波器。 代码如下： img=imread(lena.png) img=rgb2gray(img); figure,subplot(1,3,1); imshow(img);title(原始图像); img2=imnoise(img,salt & pepper,0.02); subplot(1,3,2);imshow(img2);title(椒盐噪声图像); img3=imnoise(img,gaussian,0.02); subplot(1,3,3),imshow(img3); title(高斯噪声图像); %对椒盐噪声图像进行滤波处理 h=fspecial(average,3); I1=filter2(h,img2)/255; I2=medfilt2(img2,3 3); figure,subplot(2,2,1),imshow(img),title(原图像); subplot(2,2,2),imshow(img2),title(椒盐噪声图); subplot(2,2,3),imshow(I1),title(3*3 均值滤波图); subplot(2,2,4),imshow(I2),title(3*3 中值滤波图); %对高斯噪声图像进行滤波处理 G1=filter2(h,img3)/255; G2=medfilt2(img3,3 3); figure,subplot(2,2,1),imshow(img),title(原图像); subplot(2,2,2),imshow(img3),title(高斯噪声图); subplot(2,2,3),imshow(G1),title(3*3 均值滤波图); subplot(2,2,4),imshow(G2),title(3*3 中值滤波图); 实验结果如下： (3). 使用函数 imfilter 时，分别采用不同 的填充方法（或边界选项，如零填 充、replicate、symmetric、circular）进行低通滤波，显示处理后的图像。 h=fspecial(motion,50,45); %创建一个运动模糊滤波器 filteredimg=imfilter(img,h); boundaryReplicate=imfilter(img,h,replicate); boundary0=imfilter(img,h,0); boundarysymmetric=imfilter(img,h,symmetric); boundarycircular=imfilter(img,h,circular); figure,subplot(3,2,1),imshow(img),title(Original Image); subplot(3,2,2),imshow(filteredimg),title(Motion Blurred Image); subplot(3,2,3),imshow(boundaryReplicate),title(Replicate); subplot(3,2,4),imshow(boundary0),title(0-Padding); subplot(3,2,5),imshow(boundarysymmetric),title(symmetric); subplot(3,2,6),imshow(boundarycircular),title(circular); g = imfilter(f, w, filtering_mode, boundary_options, size_options)，其中，f 为输入图像，w 为滤波掩模，g 为滤波后图像。 实验结果如下： （4）.运用 for 循环，将加有椒盐噪声的图像进行 10 次，20 次均值滤波，查看其特点,显 示均值处理后的图像（提示:利用 fspecial 函数的average类型生成均值滤波器）。 代码如下： h=fspecial(average); for i=1:10 J1=imfilter(img2,h); end for j=1:20 J2=imfilter(img2,h); end figure,subplot(1,3,1),imshow(img2),title(salt & pepper Noise); subplot(1,3,2),imshow(J1),title(10 Average Filtering); subplot(1,3,3),imshow(J2),title(20 Average Filtering); 实验结果: （5）.对加入椒盐噪声的图像分别采用均值滤波法，和中值滤波法对有噪声的图像做 处理，要求在同一窗口中显示结果。 代码如下： h1=fspecial(average); J=imfilter(img2,h1); J2=medfilt2(img2); figure,subplot(1,3,1),imshow(img2),title(salt & pepper Noise); subplot(1,3,2),imshow(J),title(Averaging Filtering); subplot(1,3,3),imshow(J2),title(Median Filtering); 实验结果为： （6）.自己设计平滑空间滤波器，并将其对噪声图像进行处理，显示处理后的图像。 代码如下： m n=size(img2); figure,subplot(1,2,1),imshow(img2); s=zeros(1,9); for i=2:1:m-1 for j=2:1:n-1 h=1; for p=i-1:1:i+1 for q=j-1:1:j+1 s(h)=img2(p,q); h=h+1; end end s=sort(s); I(i,j)=s(5); end end subplot(1,2,2),imshow(I); 实验结果： 2.锐化空间滤波 （1）读出一幅图像，采用 33 的拉普拉斯算子 w = 1, 1, 1; 1 8 1; 1, 1, 1对其 进行滤波。 代码如下： img=imread(lena.png); img=rgb2gray(img); img=im2double(img); w=1,1,1; 1,-8,1; 1,1,1 k=conv2(img,w,same); imshow(k); 实验结果为： （2）编写函数 w = genlaplacian(n)，自动产生任一奇数尺寸 n 的拉普拉斯算子，如 55 的拉普拉斯算子 w = 1 1 1 1 1 1 1 1 1 1 1 1 -24 1 1 1 1 1 1 1 1 1 1 1 1 代码如下： num=input(please enter any num:); n=num; W=ones(n,n); for i=1:n for j=1:n if(i=fix(n/2)+1 & j=fix(n/2)+1) W(i,j)=n*n-1; end end end display (W); 代码运行结果为： （3）分别采用 55，99，1515 和 2525 大小的拉普拉斯算子对 blurry_moon.tif 进代码如下： 基于上一题要求形成一函数： 行锐化滤波，并利用式 完成图像的锐化增强，观察其有何不同，要求在同一窗口中显示。 function W=lapulasi(num) n=num,W=ones(n),x=fix(n/2)+1; W(x,x)=-(n*n-1); 其他代码： f=imread(moon.tif); f=im2double(f); figure,subplot(2,3,1),imshow(f),title(Original Image); w0=lapulasi(3),w1=lapulasi(5),w2=lapulasi(9); w3=lapulasi(15),w4=lapulasi(25); f0=f-imfilter(f,w0,replicate); subplot(2,3,2),imshow(f0),title(3*3 lapulasi); f1=f-imfilter(f,w1,replicate); subplot(2,3,3),imshow(f1),title(5*5 lapulasi); f2=f-imfilter(f,w2,replicate); subplot(2,3,4),imshow(f2),title(9*9 lapulasi); f3=f-imfilter(f,w3,replicate); subplot(2,3,5),imshow(f3),title(15*15 lapulasi); f4=f-imfilter(f,w4,replicate); subplot(2,3,6),imshow(f4),title(25*25 lapulasi); 实验结果为： 代码如下： I,map=imread(moon.tif); I=double(I); figure,subplot(2,3,1),imshow(I,map),title(Original Image); Gx,Gy=gradient(I); G=sqrt(Gx.*Gx+Gy.*Gy),J1=G; subplot(2,3,2),imshow(J1,map),title(Operator1 Image); J2=I;K=find(G=7); J2(K)=G(K); subplot(2,3,3),imshow(J2,map); title(Operator2 Image); J3=I; K=find(G=7); J3(K)=255; subplot(2,3,4),imshow(J3,map),title(Operator3 Image); J4=I; K=find(G=7); J4(K)=255; subplot(2,3,5),imshow(J4,map),title(Operator4 Image); J5=I; K=find(G=7); J5(Q)=255; subplot(2,3,6),imshow(J5,map),title(Operator5 Image); （4）采用不同的梯度算子对该幅图像进行锐化滤波，并比较其效果。 实验效果如下： 代码如下： I=imread(lena.png); I=rgb2gray(I); h=fspecial(sobel); h1=h*0.5; h2=h; h3=h*1.5; z1=imfilter(I,h1); z2=imfilter(I,h2); z3=imfilter(I,h3); figure,subplot(2,2,1),imshow(I),title(Original Image); subplot(2,2,2),imshow(z1);title(Vertical filtering1); subplot(2,2,3),imshow(z2);title(Vertical filtering2); subplot(2,2,4),imshow(z3);title(Vertical filtering3) （5）自己设计锐化空间滤波器，并将其对噪声图像进行处理，显示处理后的图像； 运行结果如下： 3.傅立叶变换 （1）.读出一幅图像，对其进行快速傅立叶变换，分别显示其幅度图像和相位图像。仅对 相位部分进行傅立叶反变换后查看结果图像。 代码为： img=imread(lena.png); img=rgb2gray(img); f1=fft2(img); %快速傅里叶变换 f2=log(1+abs(f1); %振幅谱 f3=fftshift(f1); f4=angle(f1); %相位谱 figure,subplot(1,3,1),imshow(img),title(Original Image); subplot(1,3,2),imshow(log(1+abs(f3),),title(amplitude spectrum); subplot(1,3,3),imshow(f4),title(phase spectrum); 实验结果为： 代码为： f=ifft2(abs(f1); figure,subplot(1,3,1),imshow(img),title(Original Image); subplot(1,3,2),imshow(log(1+abs(f3),),title(amplitude spectrum); subplot(1,3,3),imshow(log(1+abs(f),),title(absamplitude spectrum); （2）.仅对幅度部分进行傅立叶反变换后查看结果图像。 实验结果为： （3）.将图像的傅立叶变换 F 置为其共轭后进行反变换，比较新生成图像与原始图像的差 异。 代码 f1=fft2(img); f2=log(1+abs(f1); f3=fftshift(f1); f4=angle(f1); f5=-f4; f6=double(f3*exp(f4); %傅立叶变换的复共轭 f7=ifft2(f6); %反傅立叶变换 figure,subplot(1,2,1),imshow(img),title(Original Image); subplot(1,2,2),imshow(real(f7),),title(inverse fourier transform); 实验效果： 4.平滑频域滤波 （1）.设计理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器，截至频率自 选，分别给出各种滤波器的透视图。 %理想低通滤波器的透视图% a=100; b=100; U=0:a; V=0:b; M=length(U);N=length(V); D0=10; %D0 是频带的中心半径; x1=50;y1=50; x0=-50;y0=-50; m=fix(M/2); n=fix(N/2); H=zeros(M,N); n=2; for u=1:M for v=1:N a=sqrt(U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50);%D(u,v)的值 if(a=D0)%理想滤波器 H(u,v)=1; else H(u,v)=0; end end end figure; subplot(1,3,1), surf(U,V,H),title(理想低通滤波透视图) ; %2 阶巴特沃斯低通滤波透视图%55 a=100; b=100; U=0:a; V=0:b; M=length(U);N=length(V); D0=10;%W=200;%D0 是频带的中心半径;W 是频带的宽度 x1=50;y1=50; x0=-50;y0=-50; m=fix(M/2); n=fix(N/2); H=zeros(M,N); n=2; for u=1:M for v=1:N a=sqrt(U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50);%D(u,v)的值 b=1+(a/D0)2*n; H(u,v)=1/b; end end subplot(1,3,2),surf(U,V,H); title(n=2 Butterworth lowess filter) %高斯低通滤波% a=100; b=100; U=0:a; V=0:b; M=length(U);N=length(V); D0=10; %D0 是频带的中心半径 x1=50;y1=50; x0=-50;y0=-50; m=fix(M/2); n=fix(N/2); H=zeros(M,N); for u=1:M for v=1:N D1=(u-m-x0)2+(v-n-y0).2)0.5; D2=(u-m+x0)2+(v-n+y0).2)0.5; D11=(u-m-x1)2+(v-n-y1).2)0.5; D21=(u-m+x1)2+(v-n+y1).2)0.5; H(u,v) = (U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50); end end %在绘制高斯曲面的时候，加上下述代码，显示得美观 S=50; H = -H/(2*S); H = exp(H) / (sqrt(2*pi) * sqrt(S); subplot(1,3,3),surf(U,V,H),title(Gaussian lowess filter); 实验结果： 由上图可知，当频带中心宽度相同时，理想低通滤波器为圆柱形图像，二阶巴特沃斯 低通滤波器的面线比较紧凑，高斯滤波图像最为平滑。 （2）读出一幅图像，分别采用理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器对 其进行滤波（截至频率自选），再做反变换，观察不同的截止频率下采用不同低通滤波器得 到的图像与原图像的区别，特别注意振铃效应。(提示:1)在频率域滤波同样要注意到填充问 题；2）注意到(-1)x+y;) %理想低通滤波 img=imread(lena.png); img=rgb2gray(img); f=double(img); g=fft2(f); %傅立叶变换 g=fftshift(g); M,N=size(g); d0=15; m=fix(M/2); n=fix(N/2); for i=1:M for j=1:N d=sqrt(i-m)2+(j-n)2); if(d=d0) h=1; else h=0; end result(i,j)=h*g(i,j); end end result=ifftshift(result); J1=ifft2(result); J2=uint8(real(J1); figure,subplot(2,2,1),imshow(img),title(Original Image); subplot(2,2,2),imshow(J2),title(d0=15 lowpss filter); %d0=30 的理想低通滤波% d0=30; m=fix(M/2); n=fix(N/2); for i=1:M for j=1:N d=sqrt(i-m)2+(j-n)2); if(d=d0) h=1; else h=0; end result(i,j)=h*g(i,j); end end result=ifftshift(result); J1=ifft2(result); J2=uint8(real(J1); subplot(2,2,3),imshow(J2),title(d0=30 lowpss filter); %d0=100 的理想低通滤波% d0=100; m=fix(M/2); n=fix(N/2); for i=1:M for j=1:N d=sqrt(i-m)2+(j-n)2); if(d=D0) H(u,v)=1; else H(u,v)=0; end end end figure; surf(U,V,H),title(理想高通滤波透视图) ; 高通滤波器 巴特沃斯高通滤波器 %高斯高通滤波 a=100; b=100; U=0:a; V=0:b; M=length(U);N=length(V); D0=15; %D0 是频带的中心半径; x1=50;y1=50; x0=-50;y0=-50; m=fix(M/2); n=fix(N/2); H=zeros(M,N); for u=1:M for v=1:N D1=(u-m-x0)2+(v-n-y0).2)0.5; D2=(u-m+x0)2+(v-n+y0).2)0.5; D11=(u-m-x1)2+(v-n-y1).2)0.5; D21=(u-m+x1)2+(v-n+y1).2)0.5; %高斯低通曲面 H(u,v) = (U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50); end end % 添加方差后是显示的图像更加美观 S=50; H = -H/(2*S); H =- exp(H) / (sqrt(2*pi) * sqrt(S); figure,surf(U,V,H),title(Gaussian high filter); 高斯高通滤波器 （2）读出一幅图像，分别采用理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器对 其进行滤波（截至频率自选），再做反变换，观察不同的截止频率下采用不同高通滤波器得 到的图像与原图像的区别。 %理想高通滤波 img=imread(lena.png); img=rgb2gray(img); f=double(img); g=fft2(f); %傅立叶变换 g=fftshift(g); M,N=size(g); d0=15; m=fix(M/2); n=fix(N/2); for i=1:M for j=1:N d=sqrt(i-m)2+(j-n)2); if(d=d0) h=1; else h=0; end result(i,j)=h*g(i,j); end end result=ifftshift(result); J1=ifft2(result); J2=uint8(real(J1); figure,subplot(2,2,1),imshow(img),title(Original Image); subplot(2,2,2),imshow(J2),title(d0=15 high filter); d0=30; m=fix(M/2); n=fix(N/2); for i=1:M for j=1:N d=sqrt(i-m)2+(j-n)2); if(d=d0) h=1; else h=0; end result(i,j)=h*g(i,j); end end result=ifftshift(result); J1=ifft2(result); J2=uint8(real(J1); subplot(2,2,3),imshow(J2),title(d0=30 high filter); d0=80; m=fix(M/2); n=fix(N/2); for i=1:M for j=1:N d=sqrt(i-m)2+(j-n)2); if(d=d0) h=1; else h=0; end result(i,j)=h*g(i,j); end end result=ifftshift(result); J1=ifft2(result); J2=uint8(real(J1); subplot(2,2,4),imshow(J2),title(d0=80 high filter); 理想高通滤波器得到的图片效果为： 当 d0=15 时，滤波后的图像灰度变化部分基本保留； 当 d0=30 时，图像的轮廓还比较清楚； %巴特沃斯高通滤波 img=imread(lena.png); img=rgb2gray(img); f=double(img); g=fft2(f); %傅立叶变换 g=fftshift(g); M,N=size(g); nn=2; d0=15; m=fix(M/2),n=fix(N/2); for i=1:M for j=1:N d=sqrt(i-m)2+(j-n)2); if(d=0) h=0; else h=1/(1+0.414*(d0/d)(2*nn); end result(i,j)=h*g(i,j); end 当 d0=80 时，只能看到细微的图像轮廓； end result=ifftshift(result); J1=ifft2(result); J2=uint8(real(J1); figure,subplot(2,2,1),imshow(img),title(Original Imag