頻率域濾波（2）

1、頻率域濾波的基本步驟：

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

　2、從空間濾波器得到頻率域濾波器

一般，當濾波器較小時，在計算上空間濾波要比頻率域濾波更有效率，可是濾波器較大時，使用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')

 

3、繪製線框圖和表面圖

三維線框圖和表面圖對於觀察二維濾波器頗有做用，繪製大小爲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