matlab实现gabor filter 多种方式汇总

方式一：

 

function result = gaborKernel2d( lambda, theta, phi, gamma, bandwidth)% GABORKERNEL2D % Version: 2012/8/17 by watkins.song% Version: 1.0%   Fills a (2N+1)*(2N+1) matrix with the values of a 2D Gabor function. %   N is computed from SIGMA.%%   LAMBDA - preferred wavelength (period of the cosine factor) [in pixels]%   SIGMA - standard deviation of the Gaussian factor [in pixels]%   THETA - preferred orientation [in radians]%   PHI   - phase offset [in radians] of the cosine factor%   GAMMA - spatial aspect ratio (of the x- and y-axis of the Gaussian elipse)%   BANDWIDTH - spatial frequency bandwidth at half response,%       *******************************************************************%      %       BANDWIDTH, SIGMA and LAMBDA are interdependent. To use BANDWIDTH, %       the input value of one of SIGMA or LAMBDA must be 0. Otherwise BANDWIDTH is ignored.%       The actual value of the parameter whose input value is 0 is computed inside the %       function from the input vallues of BANDWIDTH and the other parameter.% %                           pi               -1    x'^2+gamma^2*y'^2             %   G(x,y,theta,f) =  --------------- *exp ([----{-------------------}])*cos(2*pi*f*x'+phi);%                      2*sigma*sigma          2         sigma^2%%%% x' = x*cos(theta)+y*sin(theta);%%% y' = y*cos(theta)-x*sin(theta);%% Author: watkins.song% Email: [e-mail]% calculation of the ratio sigma/lambda from BANDWIDTH % according to Kruizinga and Petkov, 1999 IEEE Trans on Image Processing 8 (10) p.1396% note that in Matlab log means ln  slratio = (1/pi) * sqrt( (log(2)/2) ) * ( (2^bandwidth + 1) / (2^bandwidth - 1) );% calcuate sigmasigma = slratio * lambda;% compute the size of the 2n+1 x 2n+1 matrix to be filled with the values of a Gabor function% this size depends on sigma and gammaif (gamma <= 1 && gamma > 0)n = ceil(2.5*sigma/gamma);elsen = ceil(2.5*sigma);end% creation of two (2n+1) x (2n+1) matrices x and y that contain the x- and y-coordinates of% a square 2D-mesh; the rows of x and the columns of y are copies of the vector -n:n[x,y] = meshgrid(-n:n);% change direction of y-axis (In Matlab the vertical axis corresponds to the row index% of a matrix. If the y-coordinates run from -n to n, the lowest value (-n) comes% in the top row of the matrix ycoords and the highest value (n) in the% lowest row. This is oposite to the customary rendering of values on the y-axis: lowest value % in the bottom, highest on the top. Therefore the y-axis is inverted:y = -y;% rotate x and y% xp and yp are the coordinates of a point in a coordinate system rotated by theta.% They are the main axes of the elipse of the Gaussian factor of the Gabor function.% The wave vector of the Gabor function is along the xp axis.xp =  x * cos(theta) + y * sin(theta);yp = -x * sin(theta) + y * cos(theta);% precompute coefficients gamma2=gamma*gamma, b=1/(2*sigma*sigma) and spacial frequency% f = 2*pi/lambda to prevent multiple evaluations gamma2 = gamma*gamma;b = 1 / (2*sigma*sigma);a = b / pi;f = 2*pi/lambda;% filling (2n+1) x (2n+1) matrix result with the values of a 2D Gabor functionresult = a*exp(-b*(xp.*xp + gamma2*(yp.*yp))) .* cos(f*xp + phi);%%%%%%%%  NORMALIZATION  %%%%%%%%%%%%%%%%%%%%% NORMALIZATION of positive and negative values to ensure that the integral of the kernel is 0.% This is needed when phi is different from pi/2.ppos = find(result > 0); %pointer list to indices of elements of result which are positivepneg = find(result < 0); %pointer list to indices of elements of result which are negative pos =     sum(result(ppos));  % sum of the positive elements of resultneg = abs(sum(result(pneg))); % abs value of sum of the negative elements of resultmeansum = (pos+neg)/2;if (meansum > 0)     pos = pos / meansum; % normalization coefficient for negative values of result    neg = neg / meansum; % normalization coefficient for psoitive values of resultendresult(pneg) = pos*result(pneg);result(ppos) = neg*result(ppos);end

方式二：

 

function [Efilter, Ofilter, gb] = gaborKernel2d_evenodd( lambda, theta, kx, ky)%GABORKERNEL2D_EVENODD Summary of this function goes here % Usage: %  gb =  spatialgabor(im, wavelength, angle, kx, ky, showfilter) % Version: 2012/8/17 by watkins.song % Version: 1.0 % % Arguments: %         im         - Image to be processed. %         wavelength - Wavelength in pixels of Gabor filter to construct %         angle      - Angle of filter in degrees.  An angle of 0 gives a %                      filter that responds to vertical features. %         kx, ky     - Scale factors specifying the filter sigma relative %                      to the wavelength of the filter.  This is done so %                      that the shapes of the filters are invariant to the %                      scale.  kx controls the sigma in the x direction %                      which is along the filter, and hence controls the %                      bandwidth of the filter.  ky controls the sigma %                      across the filter and hence controls the %                      orientational selectivity of the filter. A value of %                      0.5 for both kx and ky is a good starting point. % %    lambda = 3;    %   theta = 90;    %   kx = 0.5;    %   ky = 0.5; %  % % Author: watkins.song % Email: [e-mail] % Construct even and odd Gabor filterssigmax = lambda*kx;sigmay = lambda*ky;     sze = round(3*max(sigmax,sigmay));[x,y] = meshgrid(-sze:sze);evenFilter = exp(-(x.^2/sigmax^2 + y.^2/sigmay^2)/2).*cos(2*pi*(1/lambda)*x);     % the imaginary part of the gabor filteroddFilter = exp(-(x.^2/sigmax^2 + y.^2/sigmay^2)/2).*sin(2*pi*(1/lambda)*x);     evenFilter = imrotate(evenFilter, theta, 'bilinear','crop');oddFilter = imrotate(oddFilter, theta, 'bilinear','crop');       gb = evenFilter;Efilter = evenFilter;Ofilter = oddFilter;end

 

方式三：

 

function gb = gaborKernel2d_gaborfilter( lambda, theta, phi, gamma, bw)%GABORKERNEL2D_GABORFILTER Summary of this function goes here% Version: 2012/8/17 by watkins.song% Version: 1.0%%   LAMBDA - preferred wavelength (period of the cosine factor) [in pixels]%   SIGMA - standard deviation of the Gaussian factor [in pixels]%   THETA - preferred orientation [in radians]%   PHI   - phase offset [in radians] of the cosine factor%   GAMMA - spatial aspect ratio (of the x- and y-axis of the Gaussian elipse)%   BANDWIDTH - spatial frequency bandwidth at half response,%       *******************************************************************%      %       BANDWIDTH, SIGMA and LAMBDA are interdependent. To use BANDWIDTH, %       the input value of one of SIGMA or LAMBDA must be 0. Otherwise BANDWIDTH is ignored.%       The actual value of the parameter whose input value is 0 is computed inside the %       function from the input vallues of BANDWIDTH and the other%       parameter.%                            -1     x'^2 + y'^2             %%% G(x,y,theta,f) =  exp ([----{-----------------})*cos(2*pi*f*x'+phi);%                             2     sigma*sigma%%% x' = x*cos(theta)+y*sin(theta);%%% y' = y*cos(theta)-x*sin(theta);%% Author: watkins.song% Email: [e-mail]% bw    = bandwidth, (1)% gamma = aspect ratio, (0.5)% psi   = phase shift, (0)% lambda= wave length, (>=2)% theta = angle in rad, [0 pi) sigma = lambda/pi*sqrt(log(2)/2)*(2^bw+1)/(2^bw-1);sigma_x = sigma;sigma_y = sigma/gamma;sz=fix(8*max(sigma_y,sigma_x));if mod(sz,2)==0    sz=sz+1;end% alternatively, use a fixed size% sz = 60; [x y]=meshgrid(-fix(sz/2):fix(sz/2),fix(sz/2):-1:fix(-sz/2));% x (right +)% y (up +)% Rotation x_theta = x*cos(theta)+y*sin(theta);y_theta = -x*sin(theta)+y*cos(theta); gb=exp(-0.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+phi);end

 

方式四：

 

function gb = gaborKernel2d_wiki( lambda, theta, phi, gamma, bandwidth)% GABORKERNEL2D_WIKI 改写的来自wiki的gabor函数% Version: 2012/8/17 by watkins.song% Version: 1.0%%   LAMBDA - preferred wavelength (period of the cosine factor) [in pixels]%   SIGMA - standard deviation of the Gaussian factor [in pixels]%   THETA - preferred orientation [in radians]%   PHI   - phase offset [in radians] of the cosine factor%   GAMMA - spatial aspect ratio (of the x- and y-axis of the Gaussian elipse)%   BANDWIDTH - spatial frequency bandwidth at half response,%       *******************************************************************%      %       BANDWIDTH, SIGMA and LAMBDA are interdependent. To use BANDWIDTH, %       the input value of one of SIGMA or LAMBDA must be 0. Otherwise BANDWIDTH is ignored.%       The actual value of the parameter whose input value is 0 is computed inside the %       function from the input vallues of BANDWIDTH and the other%       parameter.%                            -1     x'^2 + y'^2             %%% G(x,y,theta,f) =  exp ([----{-----------------})*cos(2*pi*f*x'+phi);%                             2     sigma*sigma%%% x' = x*cos(theta)+y*sin(theta);%%% y' = y*cos(theta)-x*sin(theta);%% Author: watkins.song% Email: [e-mail]% calculation of the ratio sigma/lambda from BANDWIDTH % according to Kruizinga and Petkov, 1999 IEEE Trans on Image Processing 8 (10) p.1396% note that in Matlab log means ln  slratio = (1/pi) * sqrt( (log(2)/2) ) * ( (2^bandwidth + 1) / (2^bandwidth - 1) );% calcuate sigmasigma = slratio * lambda;sigma_x = sigma;sigma_y = sigma/gamma;% Bounding boxnstds = 4;xmax = max(abs(nstds*sigma_x*cos(theta)),abs(nstds*sigma_y*sin(theta)));xmax = ceil(max(1,xmax));ymax = max(abs(nstds*sigma_x*sin(theta)),abs(nstds*sigma_y*cos(theta)));ymax = ceil(max(1,ymax));xmin = -xmax; ymin = -ymax;[x,y] = meshgrid(xmin:xmax,ymin:ymax);% Rotation x_theta = x*cos(theta) + y*sin(theta);y_theta = -x*sin(theta) + y*cos(theta);% Gabor Functiongb= exp(-.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+phi);end

 

方式五：

 

function [GaborReal, GaborImg] = gaborKernel_matlab( GaborH, GaborW, U, V, sigma)%GABORKERNEL_MATLAB generate very beautiful gabor filter% Version: 2012/8/17 by watkins.song% Version: 1.0% 用以生成 Gabor % GaborReal: 核实部 GaborImg: 虚部% GaborH,GaborW: Gabor窗口 高宽.% U,V: 方向 大小%            ||Ku,v||^2% G(Z) = ---------------- exp(-||Ku,v||^2 * Z^2)/(2*sigma*sigma)(exp(i*Ku,v*Z)-exp(-sigma*sigma/2))%          sigma*sigma%% 利用另外一个gabor函数来生成gabor filter, 通过u,v表示方向和尺度.% 这里的滤波器模板的大小是不变的，变化的只有滤波器的波长和方向% v: 代表波长% u: 代表方向% 缺省输入前2个参数,后面参数 Kmax=2.5*pi/2, f=sqrt(2), sigma=1.5*pi;% GaborH, GaborW, Gabor模板大小% U,方向因子{0,1,2,3,4,5,6,7}% V,大小因子{0,1,2,3,4}% Author: watkins.song% Email: [e-mail] = fix(GaborH/2);HarfW = fix(GaborW/2);Qu = pi*U/8;sqsigma = sigma*sigma;Kv = 2.5*pi*(2^(-(V+2)/2));%Kv = Kmax/(f^V);postmean = exp(-sqsigma/2);for j = -HarfH : HarfH    for i =  -HarfW : HarfW              tmp1 = exp(-(Kv*Kv*(j*j+i*i)/(2*sqsigma)));        tmp2 = cos(Kv*cos(Qu)*i+Kv*sin(Qu)*j) - postmean;        %tmp3 = sin(Kv*cos(Qu)*i+Kv*sin(Qu)*j) - exp(-sqsigma/2);        tmp3 = sin(Kv*cos(Qu)*i+Kv*sin(Qu)*j);               GaborReal(j+HarfH+1, i+HarfW+1) = Kv*Kv*tmp1*tmp2/sqsigma;        GaborImg(j+HarfH+1, i+HarfW+1) = Kv*Kv*tmp1*tmp3/sqsigma;    endendend

 

最后调用方式都一样：

 

% 测试用程序theta = [0 pi/8 2*pi/8 3*pi/8 4*pi/8 5*pi/8 6*pi/8 7*pi/8];lambda = [4 6 8 10 12];phi = 0;gamma = 1;bw = 0.5;% 计算每个滤波器figure;for i = 1:5    for j = 1:8        gaborFilter=gaborKernel2d(lambda(i), theta(j), phi, gamma, bw);        % 查看每一个滤波器        %figure;        %imshow(real(gaborFilter),[]);        % 将所有的滤波器放到一张图像中查看，查看滤波器组        subplot(5,8,(i-1)*8+j);        imshow(real(gaborFilter),[]);    endend