%%%%%%%%%%%%%%%% Filters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%convolution theorem
a=zeros(64,1);
a(1:10)=1;
figure;
subplot(3,1,1); plot(a); title('signal f(x)');
axis([0 70 0 1.5])


b=conv(a,a);
subplot(3,1,2); plot(b(1:64)); title('signal f(x)*f(x)'); 

A=fft(a);
B=A.*A;
b=real(ifft(B));
subplot(3,1,3); plot(b); title('ifft(F(x).F(x))'); 



% convolution theorem  2d
% Issues:    1) conv2 assumes a origin is at 1,1 and b origin is at center.
%               Convolution theorem assumes origins are aligned
%            2) conv2 performs zero padding, while fft + conv theorem
%               assumes cyclic padding.
a = [1 2 3  ;1 1 4 ; 0 0 7];
b =  [0 1 0; 1 1 1; 0 1 0];
% perform conv with cyclic padding
aa = [a a a; a a a; a a a];
c = conv2(aa,b,'same');
c = c(4:6,4:6)

% correct b so origin is at 1,1
ifft2(fft2(f).*fft2(circshift(h,[2 2])))


% on real images:
a = readImage('lena.tif');
showImage(a);
b = zeros(size(a));
b(128:130,128:130) = 1/9;
showImage(b); imagesc(b);
c = conv2(a,b,'same'); % dont need cyclic padding here since b is all 0 on boundary
showImage(c);
% use cyclic shift just to show it is same:
c = conv2([a a a; a a a; a a a],b,'same'); 
c = c(size(a,1)+1:2*size(a,1),size(a,2)+1:2*size(a,2));
showImage(c);

A = fft2(a);
B = fft2(b);
c = ifft2(A.*B);
showImage(c)

% Cyclic shift of result is due to the fact that conv2 assumes the  
% origin of b is at the center. Correct:
B = fft2(fftshift(b)); % can use fftshift (rather than circshift) because b is even size
c = ifft2(A.*B);
showImage(c)




%%%%%%%%%%%%%%
%filtering

img = readImage('lena.tif');

ffs = fftshift(fft2(img)); %the real fft shifted for symmetry

ffl = log(abs(ffs)+1); %the log of fft to be able to see it in good visible proportions

%this is a meshgrid to use with the calculation of the round filter
[x y] = meshgrid(-1:(2/size(img,1)):1-(2/size(img,1)), -1:(2/size(img,2)):1-(2/size(img,2)));
circledist = sqrt(x.^2+y.^2);

showImage(img);
fimg = gcf;
%set(fimg,'position',[10,10,400,400]);
% image(img); axis off; axis image;
% colormap(gray(256));

showImage(ffl);
imagesc([-64:63],[-64:63],ffl); 
ffilter = gcf;
%set(ffilter,'position',[10,420,400,400]);
% colormap(gray); axis image;

pause;

%Ideal Low pass
for i=0.5:-0.05:0.05,
   filter = circledist<i;
  figure(ffilter);
  imagesc(filter); colormap(gray);
  figure(fimg);
  imagesc(real(ifft2(fftshift(ffs.*filter))));  colormap(gray);
  pause;
  %imagesc([-64:63],[-64:63],ffl); colormap(gray); 
end;


%Gaussian low pass
for i=0.5:-0.05:0.05,
   filter = exp(-circledist.^2/(i^2));
  figure(ffilter);
  imagesc(filter); colormap(gray);
  figure(fimg);
  imagesc(real(ifft2(fftshift(ffs.*filter))));  colormap(gray);
  pause;
  %imagesc([-64:63],[-64:63],ffl); colormap(gray); 
end;

%Ideal high pass
for i=0.1:0.1:1,
   filter = circledist>(i) ; 
   figure(ffilter);
  imagesc(filter); colormap(gray);
  figure(fimg);
  imagesc(real(ifft2(fftshift(ffs.*filter))));  colormap(gray);
  pause;
  %imagesc([-64:63],[-64:63],ffl); colormap(gray); 
end;

%band pass;
for i=0.0:0.1:0.9,  
   filter = double(circledist<(i+0.1) & circledist >(i)) ; 
   filter=conv2(conv2(filter,1/7*[1 1 1 1 1 1 1],'same'),1/7*[1 1 1 1 1 1 1]','same');
  figure(ffilter);
  imagesc(filter); colormap(gray);
  figure(fimg);
  imagesc(real(ifft2(fftshift(ffs.*filter))));  colormap(gray);
  pause;
  %imagesc([-64:63],[-64:63],ffl); colormap(gray); 
end;


%%%