%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%
% Image Filtering
%
%       ima_out = Filter(ima,type,cutoff)
%
% This function high or low pass filter an image
%     
%       ima : real input image
%       type : 0 for high, 1 for low
%       cutoff: cutoff frequency in [0,1]. 0 is corresponds to the minimum frequency and 1
%					 to the maximum frequency.
%       ima_out : output real
%   
%
% Yvan Petillot, December 2000                                  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ima_out = Filter(ima,type,cutoff)

  % First convert ima to double
  ima = double(ima);
  % Create new figure
  figure(1)
  % Display ima
  imagesc(ima);
  axis image;
  title('Original image');
  pause(0.1);
  drawnow;
  % Get FFT:
  imafft = fftshift(fft2(fftshift(ima)));
  s = size(ima);
  maxr = 0.5*sqrt(s(1)^2+s(2)^2);
  cutoff = maxr*cutoff;
  switch(type)
    % High pass filter
    case 0
      for i = 1 : s(1)
	for j = 1 : s(2)
	  r = sqrt((i-1-s(1)/2)^2+(j-1-s(2)/2)^2);
	  if ( r < cutoff)
	    imafft(i,j) = 0;
	  end
	end
      end   
    case 1
      for i = 1 : s(1)
	for j = 1 : s(2)
	  r = sqrt((i-1-s(1)/2)^2+(j-1-s(2)/2)^2);
	  if ( r > cutoff)
	    imafft(i,j) = 0;
	  end
	end
      end
  end
  
  % Display filtered spectrum
  figure(2);
  imagesc(log(1+abs(imafft)));
  axis image;
  display('Press any key to proceed');
  pause;
  figure(3);
  ima_out = fftshift(ifft2(fftshift(imafft)));
  %keyboard;
  imagesc(abs(ima_out));
  axis image;
  title('Filtered image');