function spatialFilteringDemo(fname, type, param1, filval)
%% spatialFilteringDemo(fname, type, param1, filval)
%% Created: tjn, CS, NUIM, 6 IV 2003
%% Last modified: tjn, CS, NUIM, 23 IX 2011, added round() to constructLine
%%
%% Spatial filtering demo. High or low pass filter the Fourier spectrum of a
%% two-dimensional image. Selectively remove Fourier components of particular
%% orientation.
%%
%% Usage:
%%   sfdemo(fname)
%%   sfdemo(fname, type)
%%   sfdemo(fname, type, param1)
%%   sfdemo(fname, type, param1, filval)
%%
%% Parameters:
%% fname  : complete path and filename of image file
%% type   : demo type. 'freq' means that Fourier components are removed or retained
%%          based on their frequency (default), 'orient' means that Fourier components
%%          are removed or retained based on their orientation.
%% param1 : depends on demo type. If type is
%%          freq: param1 is the radius of the circular spatial filter, normalised
%%                to the minimum dimension of the image (default = 0.1)
%%          orient: param1 is the orientation (in terms of an angle) of the linear
%%                spatial filter (default = 0 degrees)
%% filval : value the filter should have, either 0 (default) or 1

if (nargin < 4)
  filval = 0;
end
if (filval ~= 1)
  filval = 0;
end
if (nargin < 2)
  type = 'freq';
end
if (nargin < 1)
  error('Must supply filename argument.');
end
if (strcmp(type,'orient') == 1)
  if (nargin < 3)
    param1 = 0;
  end
else
  type = 'freq';
  if (nargin < 3)
    param1 = 0.1;
  end
end

%% Read image from file
a = double(imread(fname));
a = a ./ max(max(a));
figure;imshow(a);title(['Original image -- ''',fname,'''']);

%% Fourier transform image
Fa = fftshift(fft2(a));
% calculate imshowHIGH such that spectrum can be easily appreciated on low intensity-resolution displays
imshowHIGH = 0.001*max(max(abs(Fa)));
figure;imshow(abs(Fa),[0,imshowHIGH]);title('Fourier spectrum');


%%
%% Construct Fourier filter
%%
if (strcmp(type,'freq') == 1)
  H = constructAperture(size(a), param1, filval);
elseif (strcmp(type,'orient') == 1)
  H = constructLine(size(a), param1, filval);
end
figure;imshow(H);title(['Spatial filter, param1=',num2str(param1),', filval=''',int2str(filval),'''']);


%% Filter image
Fa = H .* Fa;
figure;imshow(abs(Fa),[0,imshowHIGH]);title(['Filtered spectrum, param1=',num2str(param1),', filval=''',int2str(filval),'''']);


%% Inverse Fourier transform and display filtered image
a = ifft2(Fa);
figure;imshow(abs(a));title(['Filtered image, param1=',num2str(param1),', filval=''',int2str(filval),'''']);



%%----------------------------------------------------------------------------
%%----------------------------------------------------------------------------
function filter = constructAperture(imsize, radius, filval)
%% function filter = constructAperture(imsize, radius, filval)
%% Created: tjn, CS, NUIM, 6 IV 2003
%% Last modified: tjn, CS, NUIM, 6 IV 2003

p = min(imsize) * radius; % p is the normalised radius expressed as a number of pixels

%% set up the values for everywhere OUTSIDE the filter
filter = ones(imsize) - filval;

y = 1;
while (y <= imsize(1))
  x = 1;
  while (x <= imsize(2))
    % calculate the distance of point (y,x) from the centre of the image
    dist = sqrt((x - (imsize(2)/2))^2 + (y - (imsize(1)/2))^2);
    if (dist <= p)
      filter(y,x) = filval;
    end
    x = x + 1;
  end
  y = y + 1;
end


%%----------------------------------------------------------------------------
%%----------------------------------------------------------------------------
function filter = constructLine(imsize, d, filval)
%% function filter = constructLine(imsize, radius, filval)
%% Created: tjn, CS, NUIM, 6 IV 2003
%% Last modified: tjn, CS, NUIM, 6 IV 2003

%% set up the values for everywhere OUTSIDE the filter
filter = ones(round(imsize .* 1.42)) - filval;

filter(round(size(filter,1)/2)-8:round(size(filter,1)/2)+7,:) = filval;

filter = imrotate(filter,d); % rotate anticlockwise by d degrees

%% crop filter
yran1 = floor(size(filter,1)/2)-(floor(imsize(1)/2)-1);
yran2 = floor(size(filter,1)/2)+(floor(imsize(1)/2));
xran1 = floor(size(filter,2)/2)-(floor(imsize(2)/2)-1);
xran2 = floor(size(filter,2)/2)+(floor(imsize(2)/2));
filter = filter(yran1:yran2,xran1:xran2);

if (filval == 0)
  filter = bitor(filter, constructAperture(imsize, 0.1, 1));
else
  filter = bitand(filter, constructAperture(imsize, 0.025, 0));
end