function [y,H]=gaborf(x,CtrFreq,Angle,SigmaU,SigmaV,Domain)
%GABORF	Apply an even-symmetric 2D Gabor filter
%
%	If speed is critical, rather use SEPGABOR!
%
%	y=gaborf(x,CtrFreq,Angle,SigmaU,SigmaV,Domain) applies an
%	even-symmetric 2D Gabor filter with center frequency of radius
%	CtrFreq and angle Angle, and with SigmaU and SigmaV parameters to
%	the image x. 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 filtering is performed in the Fourier domain. If the variable
%	Domain is 'Fourier', the input and the output are the Fourier
%	transform of the input image and the output image respectively,
%	otherwise the input and the output are the regular images.
%
%	If not SigmaU and SigmaV are specified (not given, NaN or empty),
%	they are computed as SigmaU = CtrFreq / (3 * sqrt(2*log(2))) and
%	SigmaV=tan(pi/8) * CtrFreq / sqrt(2*log(2)) . SigmaU specifies the
%	spread of the passband along CtrFreq direction, while SigmaV
%	specifies the spread along the Angle direction.
%
%	[y,H]=gaborf(...) also returns the 2D frequency 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 SigmaU and
%	SigmaV values. 
%
%	See also: SEPGABOR, GABOR, GAUSSLP, GAUSSLPF

[VDim,HDim]=size(x);

if nargin<6,
  Domain='Spatial';
  if nargin<5,
    SigmaV=[];
    if nargin<4,
      SigmaU=[];
    end
  end
end

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

if strcmp(lower(Domain),'fourier'),
  X=x;
else
  % Determine if the input should be symmetrically extended prior to filtering
  % to have power of two dimensions.
  if (VDim == 2^(ceil(log2(VDim)))) & (HDim == 2^(ceil(log2(HDim)))),
    SymmExt=0;
  else
    SymmExt=1;
    [x,TopMarg,BotMarg,LftMarg,RgtMarg]=symmext(x);
  end

  % Transform the input image into the Fourier domain
  X=fft2(x);
end
[M,N]=size(X);

% Generate the filter (in the Fourier-domain)
[j,i]=meshgrid(1:N/2,1:M/2);
i=[i i]; i=[i;-flipud(i)]; i=i/M;
j=[j;j]; j=[j -fliplr(j)]; j=j/N;
u=i.*cos(Angle) + j.*sin(Angle);
v=j.*cos(Angle) - i.*sin(Angle);
H=(exp(-(u-CtrFreq).^2/(2*SigmaU^2) - (v*M/N).^2/(2*SigmaV^2))...
    + exp(-(u+CtrFreq).^2/(2*SigmaU^2) - (v*M/N).^2/(2*SigmaV^2)))...
    /(2*pi*SigmaU*SigmaV);

% Filter the image
Y=X.*H;

if strcmp(lower(Domain),'fourier'),
  y=Y;
else
  % Transform the filtered image back to the spatial domain, and scale it 
  % appropriately
  y=M*N*ifft2(Y);
  y=abs(y); % Remove "rounding errors" in the imaginary part
  
  % If the input has been symmetrically extended, remove the symmetric
  % extension 
  if SymmExt,
    y=y(TopMarg+1:TopMarg+VDim,LftMarg+1:LftMarg+HDim);
  end
end








