頻率域濾波（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