function [y,h,hi,hj]=sepgabor(x,CtrFreq,Angle,SigmaX,SigmaY,nstd)
%SEPGABOR  Apply a separable even-symmetric 2D Gabor filter
%
%	y=sepgabor(x,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 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
%	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, as a separable
%	filter. If the angle is not 0, the image is rotated, in order to be
%	able to use a separable filter. If the x is rotated, it is
%	symmetrically extended first.
%	
%	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.
%	
%	[y,h]=sepgabor(...) also returns the unit pulse response of the
%	Gabor filter.
%	
%	[y,h,hi,hj]=sepgabor(...) also returns the unit pulse responses of the
%	2D Gabor filter, and the vertically and horizontally separated parts.
%	
%	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.
%	
%	Example:
%	  load clown
%	  [y,h,hi,hj]=sepgabor(X,2.^(-5)*sqrt(2),pi/4);
%	  gimage(y);
%	
%	See also: GABOR, GABORF, GAUSSLP, GAUSSLPF, SEPGABOR_TEST


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

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

[M,N]=size(x);

% Determine the impulse response length and the symmetrical extension
ilim=round(nstd*SigmaX);
iext=max(ilim, abs(round(N*sin(Angle))));
jlim=round(nstd*SigmaY);
jext=max(jlim, abs(round(M*sin(Angle))));

% Generate the filter impulse response
i=(-ilim:ilim)';
j=(-jlim:jlim);
hi=exp(-1/2*(i.^2/SigmaX^2 + 0.^2/SigmaY^2)).*cos(2*pi*CtrFreq*i);
hj=exp(-1/2*(0.^2/SigmaX^2 + j.^2/SigmaY^2)).*cos(2*pi*CtrFreq*0);
hi=hi/sum(abs(hi));
hj=hj/sum(abs(hj));

% Rotate the image
UseImgRot=1;
if M==512 & N==512
  UseImgRot=0;
end
if UseImgRot
  x=imgrot(x,-Angle,'symmext');
  x=NaN2Val(x,0);
else
  ext=ceil(sqrt(2)/2*max(M,N)/2);
  x=imrotate(symmext(x,ext,ext,ext,ext),-Angle*180/pi,'crop');
end

% Filter the image, in each dimension separately
y=filter2(hj,filter2(hi,x));
[yM,yN]=size(y);

% Rotate the image back
if UseImgRot
  if yM>=M & yN>=N,
    y=imgrot(y,Angle,'crop');
  else
    y=imgrot(y,Angle); 
  end
  y=NaN2Val(y,0);
else
  y=imrotate(y,Angle*180/pi,'crop');
  y=y(ext+(1:M),ext+(1:N));
end

% Extract the non-extended regions
[yM,yN]=size(y);
iext=(yM-M)/2; jext=(yN-N)/2;
y=y(iext+1:M+iext,jext+1:N+jext);

if nargout>=2,
  h=imgrot(hi*hj,Angle);
end
