频率域滤波（2）

一、频率域滤波的基本步骤：

1）使用函数tofloat把输入图像转换为浮点图像（im2double函数也可以）

[f,revertclass] = tofloat(f)

　2）使用函数paddedsize获得填充参数

PQ = paddedsize(size(f));

　3)得到有填充图像的傅里叶变换

F = fft2(f,PQ(1),PQ(2));

　4)生成一个大小为PQ(1)*PQ(2)的滤波器函数，用freqz2函数，如果它是居中的，那么在使用之前要反变换一下(h为空间滤波器)

H = freqz2(h,PQ(1),PQ(2));
H = ifftshift(H);

　5）用滤波器乘以该变换，然后进行傅里叶反变换

G = H.*F;
g = ifft2(G);

　6）将左上部的矩形修剪为原始大小，然后将图像转换为输入图像的类

g = g(1:size(f,1),1:size(f,2));
g = revertclass(g);   %或者 g = im2uint8(g);

 

可以把这个过程封装成一个M函数，然后用的时候直接调用，函数如下（输入为原图像和滤波器）：

%用于频率域滤波的M函数

function g = dftfilt(f,H,classout)
%f = im2double(f);  %转化为浮点图像
[f,revertclass] = tofloat(f);
F = fft2(f,size(H,1),size(H,2));   %得到有填充图像的傅里叶变换
g = ifft2(H.*F);          %滤波（相乘），傅里叶反变换
g = g(1:size(f,1),1:size(f,2));    %裁剪为原始大小
 if nargin ==2||strcmp(classout,'original')
    % g = im2uint8(g);
    g = revertclass(g);
 elseif strcmp(classout,'fltpoint')
     return
 else
     error('undefine class for the output image')
 end
 

　二、从空间滤波器获得频率域滤波器

通常，当滤波器较小时，在计算上空间滤波要比频率域滤波更有效率，但是滤波器较大时，使用FFT算法的滤波处理比空间实现快。

将空间滤波器变换为频率域滤波器需要用freqz2函数。比如下面这个例子可以比较空间滤波和频率域滤波（sobel滤波器和其对应的频率域滤波器）：

          

这里构建的滤波器类似于高通滤波器，它在变换的中心位置是0，其他位置为1，当建立乘积H(u,v)*F(u,v)时，该滤波器将抑制F(u,v)的直流项而通过其他所以项。直流项决定图像的平均灰度，因此将其置0会把输出图像的平均灰度减小为0。因此，图像会变得更暗。上溯例子的代码如下：

%空间滤波和频域滤波的比较
f=imread('G:\数字图像处理（冈萨雷斯）\DIP3E_CH04_Original_Images\DIP3E_Original_Images_CH04\Fig0438(a)(bld_600by600).tif');
 subplot(231);imshow(f);title('原图')
 f = im2double(f);
 F = fft2(f); 
 S = fftshift(log(1 + abs(F)));    %取傅里叶谱，然后进行对数变换增强细节，最后将频域中心转移到矩阵中心
 subplot(232);imshow(S,[]);title('傅里叶谱');

h = fspecial('sobel')'; %空间域sobel滤波器
PQ = paddedsize(size(f));
H = freqz2(h,PQ(1),PQ(2));   %转换为频域滤波器
H1 = ifftshift(H);            %用于重排数据，使得原点位于频率矩形左上角
subplot(233);imshow(abs(H),[]);title('滤波器的图像显示');
subplot(234);imshow(abs(H1),[]);title('反中心化滤波器的图像显示');

subplot(235);freqz2(h); title('频率域滤波器的绝对值的三维透视图');     %freqz2被写成没有输出参量的形式，显示为三维透视图
h1=ifftshift(h);
subplot(236);freqz2(h1);title('反中心化的同一滤波器的三维透视图');
figure; gs = imfilter(f,h);       %空间域滤波
subplot(221);imshow(gs,[]);title('空间滤波结果');
gf = dftfilt(f,H1);               %频域滤波
subplot(222);imshow(gf,[]);title('频率域滤波');

subplot(223);imshow(abs(gs),[]);title('空间结果的绝对值')
subplot(224);imshow(abs(gf),[]);title('频率滤波结果的绝对值（频谱）')

figure; subplot(121);imshow(abs(gs)>0.2*abs(max(gs(:))));title('阈值处理的结果1'); %阈值处理，只显示比gs、gf的最大值的20%的边缘
        subplot(122);imshow(abs(gf)>0.2*abs(max(gf(:)))); title('阈值处理结果2')

 

三、绘制线框图和表面图

三维线框图和表面图对于观察二维滤波器很有作用，绘制大小为M*N的二维函数H的最简单方法是使用mesh函数，可以每隔k个点进行绘制：

mesh(H(1:k:end,1:k:end))

 

matlab默认绘出彩色网格图，用colormap([0 0 0])将线框图设置为黑色。默认的还自带网格和坐标轴，使用grid off关闭网格，grid on打开网格，axis on打开坐标轴，axis off关闭坐标轴。还可以设置观察点的位置，view(az,el)，az和el分别代表方位角和仰角。用[az,el]=view可以查看当前的几何视图。

有时候需要绘制成表面图而不是线框图，函数surf有这个功能，语法为surf(f）；这中跟mash差不多，知识网格中的四边形会被彩色填充（小面描影），用colormap(gray)可以将彩色转换为灰色，用shading interp可以去掉网格线。简单的例子（高斯低通滤波器）如下：

                                       

具体代码如下：

%线框图，高斯低通滤波器   图3.15
H = fftshift(lpfilter('gaussian',500,500,50));
mesh(double(H(1:10:500,1:10:500)));
grid off
axis off
%colormap([0 0 0]); %将线框图设置为黑色

figure;
H = fftshift(lpfilter('gaussian',500,500,50));
surf(double(H(1:10:500,1:10:500)));      %表面图
grid off        %去掉网格
axis off         %去掉坐标轴
%colormap(gray)    %转换为灰色
shading interp        %平滑小面描影并去掉网格线

[x,y] = meshgrid(-2:0.1:2,-2:0.1:2); 
z = x.*exp(-x.^2-y.^2);    %一个函数的表面图。
figure; surf(z);
grid off
axis off
shading interp