这篇文章实际上是笔者在学习冈萨雷斯 数字图像处理 matlab 版本 的第四章时,自己动手在matlab里敲入书上的程序实验得到的。
这个DFT的滤波步骤很重要,应该弄清楚。
1, 使用函数paddedsize()获得填充参数
FQ=paddedsize(size(I));%I为原始图像灰度矩阵
2, 得到使用填充的傅里叶变换
F=fft2(I,PQ(1,),PQ(2));
3, 使用任何一种方法,例如lpfilter()生成一个大小为PQ(1)*PQ(2)的滤波函数H。这个函数如果居中,就要在使用前令H=fftshift(H).
3, 将变换乘以滤波函数:
G=H.*F;
5, 获得G的傅里叶变换的实部:
g=real(ifft2(G));
6, 将左上部的矩形修剪为原始大小:
g=g(1:size(I,1),1:size(I,2));
%%
这条语句要注意,生成一个二维频率域滤波器。
H=freqz2(h,R,C) 其中h是一个二维空间滤波器。
%频域变换
%
I=imread(‘bw.tif‘);
figure;
imshow(I);
title(‘一副简单的图像‘);
result:
%fft
figure;
F=fft2(I);
S=abs(F);
imshow(S,[]);
title(‘fft频谱‘);
figure;
Fc=fftshift(F);
Sc=abs(Fc);
imshow(Sc,[]);
title(‘居中的fft频谱‘);
%对数增强并居中的fft
S2=log(1+S);
S2=fftshift(S2);
imshow(S2,[]);
%下面我们看下傅里叶反变换
I2=I(250:280,250:280);
I2=double(I2);
F2=fft2(I2);
I3=real(ifft2(F2));
imshow(I3,[]);
%使用填充和不填充的滤波效果
%不填充的滤波
[M,N]=size(I);
F=fft2(I);
sigma=10;
H=lpfilter(‘gaussian‘,M,N,sigma);
G=H.*F;
g=real(ifft2(G));
figure;
imshow(g,[]);
%下面是使用填充的滤波
PQ=paddedsize(size(I));
Fp=fft2(I,PQ(1),PQ(2));%compute the FFT with padding
Hp=lpfilter(‘gaussian‘,PQ(1),PQ(2),2*sigma);%这里使用2*sigma的原因是现在滤波器的大小是不使用填充时的滤波器大小的两倍
Gp=Hp.*Fp;
gp=real(ifft2(Gp));
gpc=gp(1:size(I,1),1:size(I,2));
figure;
imshow(gp,[]);
%执行如下的空间滤波可以得到类似的结果
h=fspecial(‘gaussian‘,15,7);
gs=imfilter(I,h);
figure;
imshow(gs);
%从空间域中获得频域滤波器,
f=imread(‘house.tif‘);
F=fft2(f);
S=fftshift(log(1+abs(F)));
S=gscale(S);
imshow(S);
h=fspecial(‘sobel‘)‘;
% h=[ 1 0 -1;2 0 -1;1 0 -1];
freqz2(h);
PQ=paddedsize(size(f));
H=freqz2(h,PQ(1),PQ(2));
H1=ifftshift(H);
imshow(abs(H),[]);
figure,imshow(abs(H1),[]);
gs=imfilter(double(f),h);
gf=dftfilt(f,H1);%频率域滤波函数,第一个参数为原始图像,后者为二维频率域滤器
imshow(gs,[]);%空间域滤波结果,sobel算子提取边缘,这里我们提取的是垂直边缘,要注意h是转置后的sobel 算子。
figure,imshow(gf,[]);%空间域滤波结果
%通过创立一副阈值二值图像,我们可以更清楚的看到边缘,这里我们上传水平滤波的结果,只需把上面的:h=fspecial(‘sobel‘)‘;改成h=fspecial(‘sobel‘);即可
figure,imshow(abs(gs)>0.2*abs(max(gs(:))));
figure,imshow(abs(gf)>0.2*abs(max(gf(:))))
%%下面的函数俺就不给结果图了。
%下面这个函数可以绘制x=1:M,y=1:N,[M,N]=size(H)的线框图
mesh(abs(H));
H=fftshift(lpfilter(‘gaussian‘,500,500,50));
figure(2);
mesh(H(1:10:500,1:10:500))
axis([0 50 0 50 0 1])
%上面的线框图默认为彩色,它从底部的蓝色渐变到顶部的红色。下面的几行命令可将线条从彩色转变为黑色
colormap([0 0 0])
axis off;
grid off
view(-25,30);%-25 代表方位角, 30代表仰角,这里view() 是使上面的语句而不是后面的语句起作用,类似title().
surf(H);