圖像分析:Gabor濾波器解析與編程


 

     在圖像分析里面,Gabor濾波器是應用非常廣泛的一種工具,也算是我本科第一次接觸圖像分析領域時用到的一種工具,剛開始學Gabor的時候,因為找不到門路,看到它的幾種形式特別傷腦筋,不知道這幾種之間到底是個什么關系,為它的復雜形式糾結了好久。后來隨着學習的深入,發現Gabor其實是一種很“傻瓜式的萬用工具”,可以應用在很多場合,就是參數設置麻煩了點。在應用的時候,只要把它當成一個濾波器來對待就行,當然,若是有時間那是最好去理解一下它公式的含義。

Gabor濾波器的歷史過程
十九世紀初,法國科學家傅里葉在向巴黎科學院呈交的關於熱傳導的著名論文中提出了Fourier級數,其意義是無法估量的。今天,Fourier分析方法已成為了各種信號數據處理中最基本的數學工具。學過《信號與系統》和《數字信號處理》的同學都知道,Fourier變換是在整體上將信號分解為不同的頻率分量,但是它只適用於統計量不隨時間變化的平穩信號,因為它並不能告之某種頻率發生在哪些時間內,而這對於非平穩信號是十分重要的。
為了研究信號在局部范圍內的頻率特性,Dennis Gabor於1946年在“Theory of communication”一文中提出了著名的“窗口”傅里葉變換(也叫短時Fourier變換,STFT),即Gabor變換。
Gabor濾波器解析
Gabor濾波器實質上一種加了高斯窗的傅里葉變換,對於圖像來說,窗函數決定了它在空域的局部性,所以可以通過移動窗口的中心來獲得不同位置的空間域信息。此外,由於髙斯函數在經過傅里葉變換以后仍然是高斯函數,這就使得Gabor變換在頻域上仍然是局部的。
1985年,Daugman將Gabor函數擴展為2維形式,並在此基礎上構造了2D Gabor濾波器。人們驚奇地發現,2D Gabor濾波器不僅同樣可以同時獲取時間和頻率域的最小不確定性,而且與哺乳動物視網膜神經細胞的接收模型相吻合。
Daugman認為,盡管Gabor濾波器的基函數不能構成一個完備的正交集,Gabor濾波器也可以看作是一種小波濾波器。他給出了2D Gabor濾波器基函數的一般形式:
從上式可以看出,Gabor函數是一個被復正弦函數調制的Gaussian函數。其中λ和θk分別是正弦波的波長和方向;θk的定義如下:
k決定了濾波器方向的個數;σ_x和σ_y分別為高斯包絡在x方向和y方向上的標准差,它們決定了高斯包絡的空間擴展;
Gabor濾波器的參數設置:目前主流的方法是基於一定的約束條件,通過實驗進行Gabor濾波器的參數選擇。
 
Gabor濾波器的缺點:計算量過大。因為(1)Gabor濾波器的實現上需要耗費時間的卷積操作,然而卻沒有行之有效的快速算法;(2)Gabor特征維數較高,Gabor濾波器的優點是多通道、多分辨率分析,通過Gabor濾波器抽取的特征通常能夠取得較高的識別率,但也產生了高維數的Gabor特征矢量,由此帶來了較大的計算量和存儲負擔。例如,對於96×48的樣本圖像來說,當選擇3尺度、8方向時,Gabor特征矢量會達到110,592(86×48×3×8)。
 
 
Gabor濾波器程序實現
matalb實現:
function [G,gabout] = gaborfilter1(I,Sx,Sy,f,theta)
%The Gabor filter is basically a Gaussian (with variances sx and sy along x and y-axes respectively)
%modulated by a complex sinusoid (with centre frequencies U and V along x and y-axes respectively) 
%described by the following equation
%%%%%%%%%%%%%%%%%%%%%%
%                                  -1     x' ^      y' ^             
%%% G(x,y,theta,f) =  exp ([----{(----) 2+(----) 2}])*cos(2*pi*f*x'); cos代表實部
%                                   2     sx'       sy'
%%% x' = x*cos(theta)+y*sin(theta);
%%% y' = y*cos(theta)-x*sin(theta);
%% Describtion :
%% I : Input image
%% Sx & Sy : Variances along x and y-axes respectively 方差
%% f : The frequency of the sinusoidal function
%% theta : The orientation of Gabor filter
%% G : The output filter as described above
%% gabout : The output filtered image
% %%isa判斷輸入參量是否為指定類型的對象
if isa(I,'double')~=1 
    I = double(I);
end
%%%%Sx,Sy在公式里分別表示Guass函數沿着x,y軸的標准差,相當於其他的gabor函數中的sigma. 
%%同時也用Sx,Sy指定了gabor濾波器的大小。(濾波器矩陣的大小)
%%這里沒有考慮到相位偏移.fix(n)是取小於n的整數(往零的方向靠)
for x = -fix(Sx):fix(Sx)
    for y = -fix(Sy):fix(Sy)
        xPrime = x * cos(theta) + y * sin(theta);
        yPrime = y * cos(theta) - x * sin(theta);
        G(fix(Sx)+x+1,fix(Sy)+y+1) = exp(-.5*((xPrime/Sx)^2+(yPrime/Sy)^2))*cos(2*pi*f*xPrime);
    end
end
Imgabout = conv2(I,double(imag(G)),'same');
Regabout = conv2(I,double(real(G)),'same');
gabout = sqrt(Imgabout.*Imgabout + Regabout.*Regabout);

OpenCV中的實現:

#include "precomp.hpp"

/*
 Gabor filters and such. To be greatly extended to have full texture analysis.
 For the formulas and the explanation of the parameters see:
 http://en.wikipedia.org/wiki/Gabor_filter
*/

cv::Mat cv::getGaborKernel( Size ksize, double sigma, double theta,
                            double lambd, double gamma, double psi, int ktype )
{
    double sigma_x = sigma;
    double sigma_y = sigma/gamma;
    int nstds = 3;
    int xmin, xmax, ymin, ymax;
    double c = cos(theta), s = sin(theta);

    if( ksize.width > 0 )
        xmax = ksize.width/2;
    else
        // cvRound返回和參數最接近的整數值;fabs返回浮點數的絕對值
        xmax = cvRound(std::max(fabs(nstds*sigma_x*c), fabs(nstds*sigma_y*s)));

    if( ksize.height > 0 )
        ymax = ksize.height/2;
    else
        ymax = cvRound(std::max(fabs(nstds*sigma_x*s), fabs(nstds*sigma_y*c)));

    xmin = -xmax;
    ymin = -ymax;
    
    //CV_Assert()若括號中的表達式值為false,則返回一個錯誤信息。
    CV_Assert( ktype == CV_32F || ktype == CV_64F );
    
    //初始化kernel矩陣
    Mat kernel(ymax - ymin + 1, xmax - xmin + 1, ktype);

    double scale = 1;
    double ex = -0.5/(sigma_x*sigma_x);
    double ey = -0.5/(sigma_y*sigma_y);
    double cscale = CV_PI*2/lambd;

    for( int y = ymin; y <= ymax; y++ )
        for( int x = xmin; x <= xmax; x++ )
        {
            double xr = x*c + y*s;
            double yr = -x*s + y*c;
            // real gabor
            double v = scale*exp(ex*xr*xr + ey*yr*yr)*cos(cscale*xr + psi);
            if( ktype == CV_32F )
                kernel.at<float>(ymax - y, xmax - x) = (float)v;
            else
                kernel.at<double>(ymax - y, xmax - x) = v;
        }

    return kernel;
}

 

Reference:

[1] Movellan J R. Tutorial on Gabor filters[J]. Open Source Document, 2002.

[2] Gabor filter for image processing and computer vision.http://matlabserver.cs.rug.nl/edgedetectionweb/web/edgedetection_params.html

[3] Prasad V S N, Domke J. Gabor filter visualization[R]. Technical Report, University of Maryland, 2005.

[4] Gabor濾波簡介和實現(Matlab,OpenCV)  http://blog.163.com/huai_jing@126/blog/static/171861983201172091718341/

[5] gabor濾波器的幾種實現方式 http://blog.csdn.NET/watkinsong/article/details/7872764

[6] 對於gabor變換和gabor小波變換理解與總結

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM