Matlab DIP(瓦)ch4图像频域滤波练习

     这一章也是按照网萨雷斯的matlab书练习的，主要讲的是图像的频域滤波方面的只是，这一次的代码相对于第3章稍加了些注释。有几个函数的功能还不是特别明确。练习代码如下：

%% fftshift 对数变换,所应用的图片本身很简单，就只有黑白2种颜色
clc
clear
f = imread('.\images\dipum_images_ch04\Fig0403(a)(image).tif');
imshow(f)
title('原始图像')

imfinfo('.\images\dipum_images_ch04\Fig0403(a)(image).tif');%此处如果用Imfinfo(f)就会报错fft

%没有居中的傅里叶频谱
F=fft2(f);%进行二维快速傅里叶变换，其结果和DFT的一样，只是计算机的计算速度变快了而已，因而叫fft
S=abs(F);%求傅里叶变换后的幅值
figure,subplot(121),imshow(S,[])
title('傅里叶频谱图像1');%title函数一定要放在坐标显示的下一句才有效。

subplot(122),imshow(S),title('傅里叶频谱图像2')

%居中的傅里叶频谱
Fc=fftshift(F);%将频谱图像原点移至图像矩形中间
S1=abs(Fc);
figure,subplot(121),imshow(S1,[]);%加了第二个参数后显示的图像正常
subplot(122),imshow(S1);%当没有第二个参数时，显示的图像为竖线加一些孤立的黑点，why？

%使用对数后视觉增强后的傅里叶频谱
S2=log(1+S1);
imshow(S2,[]);

%%用fftshift对数变换显示稍复杂的图片
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0409(a)(bld).tif');
imshow(f);

f=double(f);%其实这句不要试过对后面的变换结果也没有影响
F=fft2(f);
Fc=fftshift(F);
S=abs(Fc);
S2=log(1+S);%如果没有这句的话，那么根本看不到细节的图，所以一定要用对数压缩增加对比度
figure,imshow(S2,[])

%%试一下ifft功能
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0409(a)(bld).tif');
imshow(f);

f=double(f);
f(1:8,1:8)
F=fft2(f);
f=ifft2(F);%取傅里叶变换后的反傅里叶变换,直接是整数， 并不需要像某些代码一样取real部分
f(1:8,1:8)

%%先0扩充再滤波
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0405(a)(square_original).tif');
imshow(f);

[M N]=size(f);
F=fft2(f);%没有经过0扩充直接计算fft
sig=10;%高斯滤波参数
H=lpfilter('gaussian',M,N,sig);
G=H.*F;%加了.号的乘法表示对应每个元素都相乘，没加.号的乘法说明是真正的矩阵乘法
g=real(ifft2(G));
figure,imshow(g,[]);

PQ=paddedsize(size(f));%PQ为0扩展的面积，这里设置为与图像的大小相同？这里size(f)的大小为256*256
Fp=fft2(f,PQ(1),PQ(2));%如果fft2有3个参数的话，则是进行了0扩展了
Hp=lpfilter('gaussian',PQ(1),PQ(2),2*sig);%构造高斯滤波器
Gp=Hp.*Fp;%对0扩展后的图像进行高斯滤波
gp=real(ifft2(Gp));
figure,imshow(gp,[])'

gpc=gp(1:size(f,1),1:size(f,2));%取gp图像的1到f的第一维的行，1到f第二维的列部分，即取与图像大小相同的部分
figure,imshow(gpc,[]);%此时显示的结果应该与g图像一样。

%直接使用空间域滤波
h=fspecial('gaussian',15,7);
gs=imfilter(f,h);
figure,imshow(gs,[])

%%
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0405(a)(square_original).tif');
f=double(f);
imshow(f);

PQ=paddedsize(size(f));%size(f)的大小为256*256
sig=10;
H = lpfilter('gaussian',PQ(1),PQ(2),2*sig); % PQ=[512 512]
figure,mesh(abs(H(1:10:512,1:10:512)));%画出滤波器的三维空间图
g=dftfilt(f,H);
figure,imshow(g,[])

%%空间滤波和频域滤波的比较
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0409(a)(bld).tif');
imshow(f);

F=fft2(f);
S=fftshift(log(1+abs(F)));
S=gscale(S);%用gscale来进行将数据缩放到默认值0~255
figure,imshow(S);

h=fspecial('sobel');%构造sobel空间滤波器
figure,freqz2(h);%查看滤波器的图形

PQ=paddedsize(size(f));
H=freqz2(h,PQ(1),PQ(2));%将空间滤波器h转换成频域滤波器H，但是滤波器的原点在矩阵的中心处
H1=ifftshift(H);%原点迁移到左上角

figure,mesh(abs(H1(1:20:1200,1:20:1200)));
figure,subplot(121),imshow(abs(H),[]);
subplot(122),imshow(abs(H1),[]);

gs=imfilter(double(f),h);%空间滤波，其实就是边缘检测
gf=dftfilt(f,H);%频域滤波
subplot(121),imshow(gs,[]),subplot(122),imshow(gf,[]);%将2种滤波结果显示出来

figure,imshow(abs(gs)>0.2*abs(max(gs(:))));%只显示最大值的20%的边缘
figure,imshow(abs(gf)>0.2*abs(max(gf(:))));

d=abs(gs-gf);
a=max(d(:))
b=min(d(:))

%%构造低通滤波器
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0413(a)(original_test_pattern).tif');
imshow(f)

F1=fft2(f);
imshow(log(1+abs(fftshift(F1))),[]);%直接显示取对数后的傅里叶变换

PQ=paddedsize(size(f));
[U V]=dftuv(PQ(1),PQ(2));%这个函数还真不明白什么意思，help中解释的是计算网格频率矩阵U和V，这2个矩阵是用来频率滤波的
D0=0.05*PQ(2);

F=fft2(f,PQ(1),PQ(2));%用0扩充大小为PQ(1)*PQ(2)，然后fft变换
figure,imshow(log(1+abs(fftshift(F))),[]);%显示0扩充后的fft变换图

H=exp(-(U.^2+V.^2)/(2*(D0^2)));%构造频域滤波算子
figure,imshow(fftshift(H),[]);%显示平移后滤波器算子
g=dftfilt(f,H)
figure,imshow(g,[]);

%%绘制一个低通滤波器的mesh图
clc
clear
H1=lpfilter('gaussian',500,500,50)%构造一个中心点在左上角的高斯低通滤波器,size为500*500
H2=fftshift(H1);%将中心点平移至矩阵中心
mesh(H1(1:10:500,1:10:500)),figure,mesh(H2(1:10:500,1:10:500));%绘出2者的mesh图
axis([0 50 0 50 0 1]);%xy轴范围0~50，z轴范围0~1
figure,subplot(121),imshow(H1,[]),subplot(122),imshow(H2,[]);%绘出2者的灰度图

%%绘制一个低通滤波器的mesh图
clc
clear
H=fftshift(hpfilter('ideal',500,500,100));%构造理想的高通滤波器
mesh(H(1:10:500,1:10:500));
axis([1 50 1 50 0 1]);
colormap([0 0 0]);%mesh网格全部用黑色画
axis off%关掉坐标轴显示
grid off%关掉网格显示
figure,imshow(H,[])

%%使用高通滤波器
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0413(a)(original_test_pattern).tif');
imshow(f)
title('原始图像')

PQ=paddedsize(size(f));
D0=0.05*PQ(1);
H=hpfilter('gaussian',PQ(1),PQ(2),D0);
g=dftfilt(f,H);
figure,imshow(g,[])

%%将高通滤波和直方图均衡结合起来
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0419(a)(chestXray_original).tif');
imshow(f)
title('原始图像')

PQ = paddedsize(size(f));
D0 = 0.05*PQ(1);
HBW=hpfilter('btw',PQ(1),PQ(2),D0,2);%构造高通Butterworth滤波器，截断频率为D0
H=0.5+2*HBW;
gdw=dftfilt(f,HBW);
gbw=gscale(gdw);%将灰度值自动缩放到0~255之间
figure,subplot(121),imshow(gdw,[]);
title('高通滤波后图像');

gbe=histeq(gbw,256);%直方图均衡化
subplot(122),imshow(gbe,[]);
title('高通滤波且直方图均衡化后图像');

ghf=dftfilt(f,H);%H是HBW的线性变换
ghf=gscale(ghf);
figure,subplot(121),imshow(ghf,[]);
title('高频强调滤波后图像');

ghe=histeq(ghf,256);
subplot(122),imshow(ghe,[]);
title('高频强调且均衡化后图像');

%%fftshift和ifftshift的加深理解
clc
clear
A=[2 0 0 1
    0 0 0 0
    0 0 0 0
    3 0 0 4]
B=fftshift(A)%中心点平移
C=fftshift(fftshift(A))%还原成A
D=ifftshift(fftshift(A))%也同样还原成A
E=ifftshift(A)%平移结果和B一样

注：

     freqz2 生成的滤波器原点在正中央；lpfilter（低通）生成的滤波器原点在左上角；hpfilter（高通）生成的滤波器原点在左上角

     对0扩充图像的理解还不是很透，也就是函数paddedsize不是很清楚，返回值的意义也不是很了解；对gscale函数的功能不是很了解，好像是直接将图像值压缩到指定的范围，相信以后会慢慢了解的!

    估计我的注释过程有的可能有理解错误，希望一起交流！