这一章也是按照网萨雷斯的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函数的功能不是很了解，好像是直接将图像值压缩到指定的范围，相信以后会慢慢了解的!

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