function [output,h]=gabor(input,CtrFreq,Angle,SigmaX,SigmaY,nstd)
%GABOR	Apply an even-symmetric 2D Gabor filter
%
%	If speed is critical, rather use SEPGABOR!
%
%	output=gabor(input,CtrFreq,Angle,SigmaX,SigmaY,nstd) applies an
%	even-symmetric 2D Gabor filter with center frequency of radius
%	CtrFreq and angle Angle, and with SigmaX and SigmaY parameters to
%	the image input. CtrFreq is consequently the radius of the center
%	frequency, and the range of the center frequency is 0 to 0.5 along
%	each direction. CtrFreq has therefore a range of 0 to sqrt(2)/2. The
%	impulse response will be of length nstd (default 2) "standard
%	deviations", i.e.
%	nstd*(SigmaX*abs(cos(Angle))+SigmaY*abs(sin(Angle))) times
%	nstd*(SigmaY*abs(cos(Angle))+SigmaX*abs(sin(Angle))).
%	The filtering is performed in the spatial domain. 
%	
%	If not SigmaX and SigmaY are specified (not given, NaN or empty),
%	they are computed as SigmaX = (3*sqrt(2*log(2)))/(2*pi*CtrFreq) and
%	SigmaY=sqrt(2*log(2))/(2*pi*tan(pi/8)*CtrFreq) . SigmaX specifies the
%	width of the impulse response along the specified angle and SigmaY
%	specifies the width along the ortogonal angle.
%
%	[output,h]=gabor(...) also returns the unit pulse response of the
%	Gabor filter.
%
%	The Farrokhnia approach applies Gabor filters with CtrFreq's of
%	2.^(-6:-2)*sqrt(2), Angle (0:3)*pi/4, and the default SigmaX and
%	SigmaY values. 
%
%	See also: SEPGABOR, GABORF, GAUSSLP, GAUSSLPF

if nargin<6,
  nstd=2;
  if nargin<5,
    SigmaY=[];
    if nargin<4,
      SigmaX=[];
    end
  end
end

if (isnan(SigmaX)==1) | isempty(SigmaX),
  SigmaX = (3*sqrt(2*log(2)))/(2*pi*CtrFreq);
end
if (isnan(SigmaY)==1) | isempty(SigmaY),
  SigmaY=sqrt(2*log(2))/(2*pi*tan(pi/8)*CtrFreq);
end

[M,N]=size(input);

% Generate the filter impulse response
xlim=round(nstd*(SigmaX*abs(cos(Angle))+SigmaY*abs(sin(Angle))));
ylim=round(nstd*(SigmaY*abs(cos(Angle))+SigmaX*abs(sin(Angle))));
[j,i]=meshgrid(-ylim:ylim,-xlim:xlim);
x=i.*cos(Angle) + j.*sin(Angle);
y=j.*cos(Angle) - i.*sin(Angle);
h=exp(-1/2*(x.^2/SigmaX^2 + y.^2/SigmaY^2)).*cos(2*pi*CtrFreq*x);
h=h./sum(sum(abs(h)));

% Filter the image
output=filter2(h,input);








