%%%
% Fourier2   experiments with 2D Fourier transform
%%%

clear all
ARTIFICIAL = 1;
A_SQUARE   = 1;
A_CIRCLE   = 0;


if ARTIFICIAL			% then create artificial image
  %%%
  % creat distance matrices in different metrics
  % try yourself to display them using MESH function
  % images need NOT to be square
  %%%
  im.r  = 256;	% number of rows
  im.c  = 256;	% number of columns
  vecr = [-(im.r-1)/2:(im.r-1)/2];
  vecc = [-(im.c-1)/2:(im.c-1)/2];
  matX = repmat(vecr',1,im.c);			% create matrices from vectors
  matY = repmat(vecc,im.r,1);
  matD = max(abs(matX),abs(matY));		% distances in absolute metric
  matR = sqrt(matX.^2+matY.^2);		    % distances in euclidean matrix
  if A_SQUARE
	im.or = 60;							% size of the square
	im.IM = ones(im.r,im.c);			% set background to be black
	im.IM(find(matD<=im.or))=0;			% fill square with white
  elseif A_CIRCLE
	im.or = 60;							% radius of the circle
	im.IM = ones(im.r,im.c);			% black background
	im.IM(find(matR<=im.or))=0;			% fill sphere with white
  end
else
  im.IM = imread('Images/radlines','tif')		% read an example image
  [im.r,im.c] = size(im.IM);
  vecr = [-(im.r-1)/2:(im.r-1)/2];
  vecc = [-(im.c-1)/2:(im.c-1)/2];
  matX = repmat(vecr',1,im.c);			% create matrices from vectors
  matY = repmat(vecc,im.r,1);
  matD = max(abs(matX),abs(matY));		% distances in absolute metric
  matR = sqrt(matX.^2+matY.^2);		    % distances in euclidean matrix
end

figure(1)
imshow(im.IM)
title('Original image');
% axis on
% grid

%%%
% compute FFT and shift it to be low frequecies in the centre
%%%
fd.IM = fftshift(fft2(im.IM));
fd.P  = abs(fd.IM).^2;					% power of the spectrum
figure(2)
imshow(log(fd.P)/max(log(fd.P(:))))		% display power with log scale
title('Power of Fourier spectrum of the original image');


%%%
% Create mask with different cutoff value
%%%
sd(1).t = 4;
sd(2).t = 10;
sd(3).t = 50;
sd(4).t = 100;

for i=1:4,
  %%% first create logical mask
  sd(i).mask = zeros(im.r,im.c);
  sd(i).mask(find(matD<=sd(i).t/2))=1;
  %%% use only masked part of the Fourier spectrum
  sd(i).fdIM = fd.IM.*sd(i).mask;						% masking
  sd(i).P	 = abs(sd(i).fdIM).^2;						% power
  sd(i).IM	 = real(ifft2(fftshift(sd(i).fdIM)));		% inverse FFT
  sd(i).IM	 = im2uint8(sd(i).IM/max(sd(i).IM(:)));		% just conversion
  %%% plot results
  figure(i+2)
  subplot(121)
  imshow(sd(i).IM)
  title('reconstructed image');
  subplot(122)
  imshow(log(sd(i).P)/max(log(sd(i).P(:))))
  title(sprintf('limited spectrum used, cutoff=%d',sd(i).t));
end;


  %%%
  % Low pass Butherworth filter
  %%%
for j=1:4
  %%% set parameters of the Butherworth filter
  cutoff = sd(j).t/2;
  n		 = 1;
  %%% create wieght matrix
  fdbf(j).L		= 1./(1+(matR./cutoff).^(2*n));
  %%% reconstruct image with weighted fourier transform
  fdbf(j).IMlow = real(ifft2(fftshift(fd.IM.*fdbf(j).L)));
  fdbf(j).IMlow = fdbf(j).IMlow - min(fdbf(j).IMlow(:));		% avoid <0 elements
  fdbf(j).IMlow = im2uint8(fdbf(j).IMlow/(max(max(fdbf(j).IMlow))));
  figure(j+6)
  clf
  subplot(121)
  subimage(fdbf(j).IMlow)
  title('Low pass Butherworth filter')
  subplot(122)
  mesh(fdbf(j).L(1:7:end,1:7:end));
  title(sprintf('LP Buth filter n=%d, cutoff=%d',n,cutoff))
  % rotate3d on
end












