Gabor 滤波器的原因及实现.docx
Gabor滤波器的原因及实现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 = -(GaborHeightT)/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个特征(包括均 值和方差)。由于整个变换过程代码比较复杂，这里仅提供测试代码供下载。该 代码仅计算了一个101X101尺寸的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.freqnumfprintf( 处理频率 %d n，f0);for o0=1:conf.orientnumfilter_,freq_,orient_=do_gabor(imsize,(f0-1),(o0-1),conf.kmax,conf.f,conf.sigma,conf.orient num);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)bankind.filter二fftshift(bankind.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:N3fprintf( 正在处理滤波 %d n，ind);temp = ifft2(imagefft .* bankind.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);
