matlab实现gabor滤波器的几种方式

转自：http://blog.csdn.net/watkinsong/article/details/7882443

方式一：

[csharp] view plaincopy

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: watkins.song@gmail.com

% 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 sigma

sigma = 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 gamma

if (gamma <= 1 && gamma > 0)

n = ceil(2.5*sigma/gamma);

else

n = 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 function

result = 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 positive

pneg = find(result < 0); %pointer list to indices of elements of result which are negative

pos = sum(result(ppos)); % sum of the positive elements of result

neg = abs(sum(result(pneg))); % abs value of sum of the negative elements of result

meansum = (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 result

end

result(pneg) = pos*result(pneg);

result(ppos) = neg*result(ppos);

end

方式二：

[csharp] view plaincopy

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: watkins.song@gmail.com

% Construct even and odd Gabor filters

sigmax = 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 filter

oddFilter = 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

方式三：

[csharp] view plaincopy

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: watkins.song@gmail.com

% 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

方式四：

[csharp] view plaincopy

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: watkins.song@gmail.com

% 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 sigma

sigma = slratio * lambda;

sigma_x = sigma;

sigma_y = sigma/gamma;

% Bounding box

nstds = 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 Function

gb= exp(-.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+phi);

end

方式五：

[csharp] view plaincopy

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: watkins.song@gmail.com

HarfH = 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;

end

end

end

最后调用方式都一样：

[csharp] view plaincopy

% 测试用程序

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),[]);

end

end