Gabor函数和Gabor滤波器的原理和实现

原文转自：http://www.zhizhihu.com/html/y2009/381.html

Gabor函数
　　Gabor变换属于加窗傅立叶变换，Gabor函数可以在频域不同尺度、不同方向上提取相关的特征。另外Gabor函数与人眼的生物作用相仿，所以经常用作纹理识别上，并取得了较好的效果。二维Gabor函数可以表示为：
　　其中：
　　v的取值决定了Gabor滤波的波长，u的取值表示Gabor核函数的方向，K表示总的方向数。参数决定了高斯窗口的大小，这里取。程序中取4个频率（v=0, 1, ..., 3），8个方向（即K=8，u＝0, 1, ... ,7），共32个Gabor核函数。不同频率不同方向的Gabor函数可通过下图表示：
　　图片来源：GaborFilter.html
　　图片来源：http://www.bmva.ac.uk/bmvc/1997/papers/033/node2.html
　　三、代码实现
　　Gabor函数是复值函数，因此在运算过程中要分别计算其实部和虚部。代码如下：

　　private void CalculateKernel(int Orientation, int Frequency)　　{　　double real, img;　　for(int x = -(GaborWidth-1)/2; x<(GaborWidth-1)/2+1; x++)　　for(int y = -(GaborHeight-1)/2; y<(GaborHeight-1)/2+1; y++)　　{　　real = KernelRealPart(x, y, Orientation, Frequency);　　img = KernelImgPart(x, y, Orientation, Frequency);　　KernelFFT2[(x+(GaborWidth-1)/2) + 256 * (y+(GaborHeight-1)/2)].Re = real;　　KernelFFT2[(x+(GaborWidth-1)/2) + 256 * (y+(GaborHeight-1)/2)].Im = img;　　}　　}　　private double KernelRealPart(int x, int y, int Orientation, int Frequency)　　{　　double U, V;　　double Sigma, Kv, Qu;　　double tmp1, tmp2;　　U = Orientation;　　V = Frequency;　　Sigma = 2 * Math.PI * Math.PI;　　Kv = Math.PI * Math.Exp((-(V+2)/2)*Math.Log(2, Math.E));　　Qu = U * Math.PI / 8;　　tmp1 = Math.Exp(-(Kv * Kv * ( x*x + y*y)/(2 * Sigma)));　　tmp2 = Math.Cos(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(Sigma/2));　　return tmp1 * tmp2 * Kv * Kv / Sigma;　　}　　private double KernelImgPart(int x, int y, int Orientation, int Frequency)　　{　　double U, V;　　double Sigma, Kv, Qu;　　double tmp1, tmp2;　　U = Orientation;　　V = Frequency;　　Sigma = 2 * Math.PI * Math.PI;　　Kv = Math.PI * Math.Exp((-(V+2)/2)*Math.Log(2, Math.E));　　Qu = U * Math.PI / 8;　　tmp1 = Math.Exp(-(Kv * Kv * ( x*x + y*y)/(2 * Sigma)));　　tmp2 = Math.Sin(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(Sigma/2));　　return tmp1 * tmp2 * Kv * Kv / Sigma;　　}

有了Gabor核函数后就可以采用前文中提到的“离散二维叠加和卷积”或“快速傅立叶变换卷积”的方法求解Gabor变换，并对变换结果求均值和方差作为提取的特征。32个Gabor核函数对应32次变换可以提取64个特征（包括均值和方差）。由于整个变换过程代码比较复杂，这里仅提供测试代码供下载。该代码仅计算了一个101×101尺寸的Gabor函数变换，得到均值和方差。代码采用两种卷积计算方式，从结果中可以看出，快速傅立叶变换卷积的效率是离散二维叠加和卷积的近50倍。
-----------------------------------------------------------------------

最近忙着论文，需要Gabor滤波代码，可是网上总找不到合适的代码，于是就自己编了一个，不当之处请指点。参考论文为 L. Wiskott，J. M. Fellous，N. Kruger，C. v. d.Malsburg. Face Recognition by Elastic Bunch Graph Matching，IEEE Trans. On PAMI，Vol.19，No.7，pp775-779，1997

首先实现滤波器：

function [bank] = do_createfilterbank(imsize,varargin)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 函数实现：创建Gabor 滤波组%%%% 必选参数：%   imsize - 图像大小%%% 可选参数：%   freqnum — 频率数目%   orientnum — 方向数目%   f       —   频率域中的采样步长%   kmax    —   最大的采样频率%   sigma   —   高斯窗的宽度与波向量长度的比率%%%% 返回结果：%   bank%           .freq        —     滤波频率%           .orient      —     滤波方向%           .filter      —     Gabor滤波%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%conf = struct(...,    'freqnum',3,...    'orientnum',6,...    'f',sqrt(2),...    'kmax',(pi/2),...    'sigma',(sqrt(2)*pi) ...    );conf = do_getargm(conf,varargin);bank = cell(1,conf.freqnum*conf.orientnum);for f0=1:conf.freqnum    fprintf('处理频率 %d \n', f0);    for o0=1:conf.orientnum        [filter_,freq_,orient_] = do_gabor(imsize,(f0-1),(o0-1),conf.kmax,conf.f,conf.sigma,conf.orientnum);        bank{(f0-1)*conf.orientnum + o0}.freq = freq_; %以orient增序排列        bank{(f0-1)*conf.orientnum + o0}.filter = filter_;        bank{(f0-1)*conf.orientnum + o0}.orient = orient_;    endendfor ind = 1:length(bank)    bank{ind}.filter=fftshift(bank{ind}.filter);endfunction [filter,Kv,Phiu] = do_gabor (imsize,nu,mu,Kmax,f,sigma,orientnum)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 函数实现： 创建Gabor滤波%%%% 参数：%   imsize     : 滤波的大小（即图像大小）%   nu         : 频率编号 [0 ...freqnum-1];%   mu         : 方向编号 [0...orientnum-1]%   Kmax       ： 最大的采样频率%   f          ： 频率域中的采样步长%   sigma      : 高斯窗的宽度与波向量长度的比率%   orientnum : 方向总数%%%% 返回值：%   filter :    滤波%   Kv    :    频率大小%   Phiu :    方向大小%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%rows = imsize(1);cols = imsize(2);minrow = fix(-rows/2);mincol = fix(-cols/2);row = minrow + (0:rows-1);col = mincol + (0:cols-1);[X,Y] = meshgrid(col,row);Kv = Kmax/f^nu;Phiu = pi * mu /orientnum;K = Kv * exp(i * Phiu);F1 = (Kv ^ 2)/ (sigma^2) * exp(-Kv^2 * abs(X.^2 + Y.^2) / (2*sigma^2)) ;F2 = exp(i * (real(K) * X + imag(K) * Y)) - exp(-sigma^2/2);filter = F1.* F2;

Gabor 滤波实现(1)已经创建了Gabor滤波组，现在可以使用该滤波组对图像进行转换，得到振幅和相位。

function [result] = do_filterwithbank(im,bank)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 函数实现：对图像使用Gabor滤波组进行转换换%%%% 参数：%   im       — 被转换的图像%   bank        — 由函数do_createfilterbank得到的滤波组%%%% 返回：%   result                     — 图像被转换后的结果%       .amplitudes            — 不同像素点的振幅向量%       .phases                — 不同像素点的相位向量%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[N1 N2] = size(im);N3 = length(bank);phases = zeros(N1,N2,N3);amplitudes = zeros(N1,N2,N3);imagefft = fft2(im);for ind = 1:N3    fprintf('正在处理滤波 %d \n',ind);     temp = ifft2(imagefft .* bank{ind}.filter);     phases(:,:,ind) = angle(temp);     amplitudes(:,:,ind) = abs(temp);endresult.phases = phases;result.amplitudes = amplitudes;

整个程序可以如下使用。 im = imread('image.jpg'); im = rgb2gray(im); bank = do_createfilterbank(size(im)); result = do_filterwithbank(im,bank);