% Example 1
I = imread('flowers.tif');
I = im2double(I(10+[1:256],222+[1:256],:)); % crop
FI = fftn(I);
mask = zeros(size(FI));
for n = 1:8:size(FI,1)/2-1,
    mask = zeros(size(FI));
    mask(end/2-n:end/2+n,end/2-n:end/2+n,:)=1;
    RFI = fftshift(mask).*FI;
    RI = abs(ifftn(RFI));
    figure(1);imshow([mask,RI./max(RI(:))],[]);drawnow;
    M = abs(ifftn(fftshift(mask)));
    figure(2);imshow(log(1+fftshift(M./max(M(:)))),[])
    P = log(1+fftshift(M(:,:,1)./max(M(:))));
    figure(3);surf(60*P(end/2-30:end/2+30,end/2-30:end/2+30));
    axis image;
    drawnow;pause; 
end

%%%%%%%%%%%%%%%
% Example 2, filtering additive harmonic `noise' 
% Generate noise:
[X,Y] = meshgrid(1:256, 1:256);
m=35; n=20;
s = (sin(m*2*pi*X./256)+sin(n*2*pi*Y./256)+2);
s = s./max(s(:));
s3 = ones(size(I)); for i=1:3, s3(:,:,i)=s; end;
figure(1);imshow(s3,[]); pause;% as an image

% Add noise to the image:
J = I+s3;
J = J./max(J(:));
imshow([I,J],[]);pause;

M = abs(fftn(J));
imshow(M./max(M(:)),[]); pause;
mask = ones(size(J));
mask(1,m+1,:)=0; mask(n+1,1,:)=0;
mask(1,end-m+1,:)=0; mask(end-n+1,1,:)=0;
imshow([mask,M/max(M(:))],[]);pause;
imshow(mask.*M/max(M(:)),[]);pause; % mask the noise
Rec = abs(ifftn(mask.*fftn(J)));
Rec=Rec-min(Rec(:)); % restrict to zero
imshow([I,J,Rec./max(Rec(:))],[]);pause; %  original, noisy, reconstructed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 3,  same as 1 but with a circular mask
I = imread('flowers.tif');
I = im2double(I(10+[1:256],222+[1:256],:)); % crop
FI = (fftn((I)));
[X,Y] = meshgrid(1:256, 1:256);
X = X-128; Y=Y-128;
%imshow(fftshift(log(1+1.6*FI./max(FI(:)))),[]);
mask = zeros(size(FI));
for n = 1:8:size(FI,1)/2-1,
    mask = zeros(size(FI));
    m=(X.^2+Y.^2)<n^2;
    for i=1:3, mask(:,:,i)=m; end %for There is a more efficient way 
                                   % for this operation, find it!
    RFI = fftshift(mask).*FI;
    RI = abs(ifftn(RFI));
    figure(1);imshow([mask,RI./max(RI(:))],[]);drawnow;
    
    M = abs(ifftn(fftshift(mask)));
    figure(2);imshow(log(1+fftshift(M./max(M(:)))),[])
    P = log(1+fftshift(M(:,:,1)./max(M(:))));
    figure(3);surf(60*P(end/2-30:end/2+30,end/2-30:end/2+30));
    axis image;
    drawnow;pause; 
end