Matlab频域滤波练习

%% 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一样