实验二B 图像的空域与频域滤波

• 一、实验目的
• 二、实验原理
• 三、实验内容与要求
• 四、实验的具体实现

一、实验目的

1.掌握图像滤波的基本定义及目的。
2.理解空间域滤波的基本原理及方法。
3.掌握进行图像的空域滤波的方法。
4.掌握傅里叶变换及逆变换的基本原理方法。
5.理解频域滤波的基本原理及方法。
6.掌握进行图像的频域滤波的方法。

二、实验原理

1.频域增强
频域增强是利用图像变换方法将原来的图像空间中的图像以某种形式转换到其他空间中，然后利用该空间的特有性质方便地进行图像处理，最后再转换回原来的图像空间中，从而得到处理后的图像。
频域增强的主要步骤是：
1）选择变换方法，将输入图像变换到频域空间。
2）在频域空间中，根据处理目的设计一个传递函数，并进行处理。
3）将所得结果用逆变换得到增强的图像。
4）常用的频域增强方法有低通滤波和高通滤波。
2.低通滤波
图像的能量大部分集中在幅度谱的低频和中频部分，而图像的边缘和噪声对应于高频部分。因此能降低高频成分幅度的滤波器就能减弱噪声的影响。由卷积定理，在频域实现低通滤波的数学表达式：

1)理想低通滤波器（ILPF）

2)巴特沃斯低通滤波器（BLPF）

3)指数型低通滤波器（ELPF）

3.高通滤波
由于图像中的细节部分与其高频分量相对应，所以高通滤波可以对图像进行锐化处理。高通滤波与低通滤波相反，它是高频分量顺利通过，使低频分量受到削弱。高通滤波器和低通滤波器相似，其转移函数分别为:
1)理想高通滤波器（IHPF）

2)巴特沃斯高通滤波器（BLPF）

3)指数型高通滤波器（ELPF）

图像经过高通滤波处理后，会丢失许多低频信息，所以图像的平滑区基本上会消失。所以，可以采用高频加强滤波来弥补。高频加强滤波就是在设计滤波传递函数时，加上一个大于0小于1的常数，即：

三、实验内容与要求

1.傅里叶变换
1)读出一幅灰度图像，对其进行快速傅里叶变换，分别显示其幅度图像和相位图像。
2)仅对相位部分进行傅里叶逆变换后查看结果图像。
3)仅对幅度部分进行傅里叶逆变换后查看结果图像。
4)将图像的傅里叶变换置为其共轭后进行逆变换，比较新生成图像与原始图像的差异。
2.平滑频域滤波
1)设计理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器，截止频率自选。
2)读出一幅灰度图像，分别采用理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器对其进行滤波（截止频率自选），再做逆变换，观察不同的截止频率下采用不同低通滤波器得到的图像与原图像的区别，特别注意振铃效应。
3.锐化频域滤波
1)设计理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器，截止频率自选。
2)读出一幅灰度图像，分别采用理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器对其进行滤波（截止频率自选），再做逆变换，观察不同的截止频率下采用不同高通滤波器得到的图像与原图像的区别。

四、实验的具体实现

%% 傅里叶变换% 1)读出一幅灰度图像，对其进行快速傅里叶变换，分别显示其幅度图像和相位图像。img=imread('Fig0219(a).tif');f1=fft2(img);%快速傅里叶变换f2=log(1+abs(f1));%振幅谱f3=fftshift(f1);f4=angle(f1);%相位谱figure(1)subplot(1,3,1),imshow(img),title('Original Image');subplot(1,3,2),imshow(f2,[]),title('amplitude Image');subplot(1,3,3),imshow(f4),title('pahse Image');% 2)仅对相位部分进行傅里叶逆变换后查看结果图像。f=ifft2(abs(f1));figure(2)subplot(2,3,1),imshow(img),title('Original Image');subplot(2,3,2),imshow(f4),title('pahse Image');subplot(2,3,3),imshow(log(1+abs(f)),[]),title('abspahse Image');% 3)仅对幅度部分进行傅里叶逆变换后查看结果图像。f=ifft2(abs(f1));subplot(2,3,4),imshow(img),title('Original Image');subplot(2,3,5),imshow(f2,[]),title('amplitude Image');subplot(2,3,6),imshow(log(1+abs(f)),[]),title('absamplitude Image');% 4)将图像的傅里叶变换 置为其共轭后进行逆变换，比较新生成图像与原始图像的差异。%img=imread('Fig0427(a).tif');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(3)subplot(1,2,1),imshow(img),title('Original Image');subplot(1,2,2),imshow(real(f7),[]),title('inverse fourier transform');%% 2.平滑频域滤波%1)设计理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器，截止频率自选。%%%%%理想低通滤波器透视图%%%%%a=100;b=100;U=0:a;V=0:b;M=length(U);N=length(V);D0=10;%频带中心半径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    endend%figure(4);subplot(1,3,1),surf(U,V,H),title('理想低通滤波透视图');%%%%%2阶巴特沃斯低通滤波透视图%%%%%a=100;b=100;U=0:a;V=0:b;M=length(U);N=length(V);D0=10;%频带中心半径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;    endend%subplot(1,3,2),surf(U,V,H),title('2阶巴特沃斯低通滤波透视图');%%%%%高斯低通滤波透视图%%%%%a=100;b=100;U=0:a;V=0:b;M=length(U);N=length(V);D0=10;%频带中心半径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        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));%D(u,v)    endend%subplot(1,3,3),surf(U,V,H),title('高斯低通滤波透视图');%2)读出一幅灰度图像，分别采用理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器对其进行滤波（截止频率自选），再做逆变换，观察不同的截止频率下采用不同低通滤波器得到的图像与原图像的区别，特别注意振铃效应。%%%%%%%d0=15理想低通滤波%%%%%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);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));figure(5),subplot(2,2,1),imshow(img),title('Original Image');subplot(2,2,2),imshow(J2),title('d0=15 lowpass filter');%%%%%d0=50理想低通滤波%%%%%d0=50;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);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,3),imshow(J2),title('d0=50 lowpass 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=1;        else             h=0;        end        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,4),imshow(J2),title('d0=100 lowpass filter');%%%%%%%d0=15 2阶巴特沃斯低通滤波%%%%%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);        h=1/(1+0.414*(d/d0)^(2*nn));        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));figure(6),subplot(2,2,1),imshow(img),title('Original Image');subplot(2,2,2),imshow(J2),title('d0=15 Butterworth lowpss filter');%%%%%d0=50 2阶巴特沃斯低通滤波%%%%%d0=50;m=fix(M/2);n=fix(N/2);for i=1:M    for j=1:N        d=sqrt((i-m)^2+(j-n)^2);        h=1/(1+0.414*(d/d0)^(2*nn));        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,3),imshow(J2),title('d0=50 Butterworth lowpss filter');%%%%%d0=100 2阶巴特沃斯低通滤波%%%%%0;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);        h=1/(1+0.414*(d/d0)^(2*nn));        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,4),imshow(J2),title('d0=100 lButterworth lowpss filter');%%%%%%%d0=15 高斯低通滤波%%%%%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);        h=exp(-(d.^2)./(2*(d0^2)));        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));figure(7),subplot(2,2,1),imshow(img),title('Original Image');subplot(2,2,2),imshow(J2),title('d0=15 Gaussian filter');%%%%%%d0=50 高斯低通滤波%%%%%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);        h=exp(-(d.^2)./(2*(d0^2)));        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,3),imshow(J2),title('d0=50 Guassian 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);        h=exp(-(d.^2)./(2*(d0^2)));        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,4),imshow(J2),title('d0=100 Guassian filter');%% 3.锐化频域滤波%1)设计理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器，截止频率自选。%%%%%理想高通滤波器透视图%%%%%a=100;b=100;U=0:a;V=0:b;M=length(U);N=length(V);D0=15;%频带中心半径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    endendfigure(8);subplot(1,3,1),surf(U,V,H),title('理想高通滤波透视图');%%%%%2阶巴特沃斯高通滤波透视图%%%%%a=100;b=100;U=0:a;V=0:b;M=length(U);N=length(V);D0=15;%频带中心半径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;    endendsubplot(1,3,2),surf(U,V,H),title('2阶巴特沃斯高通滤波透视图');%%%%%高斯高通滤波透视图%%%%%a=100;b=100;U=0:a;V=0:b;M=length(U);N=length(V);D0=15;%频带中心半径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));%D(u,v)    endendsubplot(1,3,3),surf(U,V,H),title('高斯高通滤波透视图');%2)读出一幅灰度图像，分别采用理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器对其进行滤波（截止频率自选），再做逆变换，观察不同的截止频率下采用不同高通滤波器得到的图像与原图像的区别。%%%%%%%d0=15理想高通滤波%%%%%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);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));figure(9),subplot(2,2,1),imshow(img),title('Original Image');subplot(2,2,2),imshow(J2),title('d0=15 high 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);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,3),imshow(J2),title('d0=30 high filter');%%%%%d0=60理想高通滤波%%%%%d0=60;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);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,4),imshow(J2),title('d0=60 high filter');%%%%%%%d0=15 2阶巴特沃斯低通滤波%%%%%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);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));figure(10),subplot(2,2,1),imshow(img),title('Original Image');subplot(2,2,2),imshow(J2),title('d0=15 Butterworth high filter');%%%%%d0=30 2阶巴特沃斯低通滤波%%%%%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==0)            h=0;        else            h=1/(1+0.414*(d0/d)^(2*nn));        end        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,3),imshow(J2),title('d0=30 Butterworth high filter');%%%%%d0=60 2阶巴特沃斯高通滤波%%%%%0;d0=60;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);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,4),imshow(J2),title('d0=60 lButterworth high filter');%%%%%%%d0=15 高斯高通滤波%%%%%d0=15;n=2;m=fix(M/2);n=fix(N/2);for i=1:M    for j=1:N        d=sqrt((i-m)^2+(j-n)^2);        h=1-exp(-(d.^2)./(2*(d0^2)));        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));figure(11),subplot(2,2,1),imshow(img),title('Original Image');subplot(2,2,2),imshow(J2),title('d0=15 Gaussian 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);        h=1-exp(-(d.^2)./(2*(d0^2)));        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,3),imshow(J2),title('d0=30 Guassian filter');%%%%%%d0=60 高斯低通滤波%%%%%d0=60;m=fix(M/2);n=fix(N/2);for i=1:M    for j=1:N        d=sqrt((i-m)^2+(j-n)^2);        h=1-exp(-(d.^2)./(2*(d0^2)));        result(i,j)=h*g(i,j);    endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));subplot(2,2,4),imshow(J2),title('d0=60 Guassian filter');