Gabor濾波原理和matlab實現

1. 傅里葉變換的缺點

傅里葉變換的公式為

從公式中可以看出，傅里葉變換對信號在整個時域做了積分處理，因此其結果對時域信號在整個時間軸上進行了信息平均。這對於平穩信號來說是可行的，然而對於在時間上具有顯著變化的非平穩信號來說，這樣的做法顯然不能滿足我們對信號進行精確分析的要求。我們希望將信號分解到不同頻率成分上來研究組成該信號的各頻率成分的含量的同時，也能看到在信號的時變過程中，到底在哪一個時間段某一頻率成分含量較多。（摘自2-D Gabor 濾波器特征提取）

2. Gabor變換

Gabor變換是D.Gabor 1946年提出的。為了由信號的傅里葉變換提取局部信息，引入了時間局部化的窗函數，得到了窗口傅里葉變換。由於窗口傅里葉變換只依賴於部分時間的信號，所以，現在窗口傅里葉變換又稱為短時傅里葉變換。Gabor變換的基本思想：把信號划分成許多小的時間間隔，用傅里葉變換分析每一個時間間隔，以便確定信號在該時間間隔存在的頻率。其處理方法是對 f(t)加一個滑動窗，再作傅里葉變換。Gabor變換所用的窗口函數是高斯函數，二維Gabor變換公式為(摘自 Gabor濾波器學習)

參數含義：

λ：正弦函數波長，它的值以像素為單位指定，通常大於等於2，但不能大於輸入圖像尺寸的1/5.

θ：Gabor核函數(濾波器)的方向,這個參數指定了Gabor函數並行條紋的方向，他的取值為0到360度 

ψ：相位偏移，調諧函數的相位偏移，取值-180到180。

σ：帶寬，高斯函數的標准差，通常取2π

γ： 空間的寬高比，決定了Gabor函數形狀的橢圓率，當γ=1時，形狀是圓的，當γ<1時，形狀隨着平行條紋方向而拉長。通常該值為0.5

在特征提取方面，Gabor小波變換與其它方法相比：一方面其處理的數據量較少，能滿足系統的實時性要求；另一方面，小波變換對光照變化不敏感，且能容忍一定程度的圖像旋轉和變形，當采用基於歐氏距離進行識別時，特征模式與待測特征不需要嚴格的對應，故能提高系統的魯棒性。

3. 用matlab手動實現Gabor濾波器

代碼：

clear;
ksize = 100;    % kernel size
d = ksize/2;
lambda = 8;     % wavelength
theta = 0;      % orientation
phase = 0;
sigma = 8;      % variation
ratio = 0.5;    % spatial aspect ratio

g = zeros(ksize, ksize);
for x = 1:ksize
    xd = x - d;  % distance from the center
    for y = 1:ksize
        yd = y - d;
        xn = xd*cos(theta) + yd*sin(theta);
        yn = -xd*sin(theta) + yd*cos(theta);
        g(x,y) = exp(-(xn^2+ratio^2*yn^2)/(2*sigma^2))*exp(1i*(2*pi*xn/lambda+phase));
    end
end

subplot(121);
gReal = real(g);
gReal = gReal/sum(sum(gReal));
mesh(1:ksize,1:ksize,gReal);
title('real part');
xlabel('x');ylabel('y');

subplot(122);
gImag = imag(g);
gImag = gImag/sum(sum(gImag));
mesh(1:ksize,1:ksize,gImag);
title('imaginary part');
xlabel('x');ylabel('y');

運行結果：

 

如果稍微轉動一下圖像，我們會發現，Gabor濾波其實就是一種小波變換（從公式上也能看出來）：

4. 用gabor濾波進行圖像特征提取

gabor函數：

function g=gabor_func_peng(ksize,lambda,theta,phase,sigma,ratio)
% input
%   ksize: kernel size
%   lambda: wavelength
%   theta: orientation
%   phase: pahse angle
%   sigma: variation
%   ratio: spatial aspect ratio
% output
%   g: gabor filter

d = ksize/2;
g = zeros(ksize, ksize);
for x = 1:ksize
    xd = x - d;  % distance from the center
    for y = 1:ksize
        yd = y - d;
        xn = xd*cos(theta) + yd*sin(theta);
        yn = -xd*sin(theta) + yd*cos(theta);
        g(x,y) = exp(-(xn^2+ratio^2*yn^2)/(2*sigma^2))*exp(1i*(2*pi*xn/lambda+phase));
    end
end

end

圖像處理函數：

function Ig=gabor_imgProcess_peng(I,ksize,lambda,theta,phase,sigma,ratio)
% input
%   I: input gray image
%   ksize: kernel size
%   lambda: wavelength
%   theta: orientation
%   phase: pahse angle
%   sigma: variation
%   ratio: spatial aspect ratio
% output
%   g: gabor filter

[m,n] = size(I);
d = ksize/2;

% pad image
Ip = zeros(m+ksize, n+ksize);
Ip(d+1:d+m, d+1:d+n)=I;

g = gabor_func_peng(ksize,lambda,theta,phase,sigma,ratio);
g = real(g); % only use the real part of gabor filter
Ig = zeros(m,n);
disp('get gabor');
for x = 1:m
    for y = 1:n
        Ig(x,y) = sum(sum(Ip(x:x+ksize-1,y:y+ksize-1).*g));
    end
end

Ig = uint8(Ig);
Ig = min(255, max(0, Ig));

主程序：

clear;

ksize = 50;     % kernel size
d = ksize/2;
lambda = 6;     % wavelength
theta = [0, pi/2];      % orientation
phase = 0;
sigma = 4;      % variation
ratio = 0.5;    % spatial aspect ratio

I = imread('../lena.jpg');
I = rgb2gray(I);

Ig_cell = cell(1,length(theta));
for k = 1:length(theta)
    Ig_cell{1,k} = gabor_imgProcess_peng(I,ksize,lambda,theta(k),phase,sigma,ratio);
end

subplot(141);
imshow(I);title('original image');
subplot(142);
imshow(Ig_cell{1,1});title('gabor image, 0 degree');
subplot(143);
imshow(Ig_cell{1,2});title('gabor image, 90 degree');
subplot(144);
imshow(Ig_cell{1,1}+Ig_cell{1,2});title('gabor image, 0 and 90 degree');

運行結果：

 

 

參考資料

[1] 代碼參考了 手動實現Gabor濾波器

[2] Gabor濾波的公式講解以及Python實現可參考油管視頻 What are Gabor filters?

[3]  Gabor濾波