% disp('This program uses a threshold (on absolute values)')
% disp('to discard Fourier coefficients.')
% disp('For nine different thresholds the program')
% disp('shows the image obtained after Fourier coefficients')
% disp('are discarded, and also displays the threshold,')
% disp('the percentage of discarded (i.e. set to 0)')
% disp('coefficients, and the MSE.')
% disp(['Keep in mind that 1/256 = ',num2str(1/256)])
% disp('This MSE value is reached with approximately')
% disp('50% of discarded coefficients,')
% disp('but fairly good results are obtained even with')
% disp('appr. 70% of discarded coefficients.')
% pause

dist = min(repmat([100:-1:0 1:99],[200 1]), repmat([100:-1:0 1:99]',[1 200]));

zeroims = cell(1,9);

im = imread('100_2007.jpg');
im = im2double(im);
im = 0.6*im(:,:,1)+0.3*im(:,:,2)+0.1*im(:,:,3);
im = im(51:250,101:300);

f = fft2(im);
fs1 = fftshift(f);

for k=[0 0.0001 0.0002 0.0005 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5 10]

    fs = fs1.*(1+dist*k);

    for i=1:9

    fs2 = fs;

    t = 2^(-7-i);

    zeroims{i} = (abs(fs2/(200*200))>=t);

    fs2 = fs2.*zeroims{i};
    fs2 = fs2./(1+dist*k);

    f2 = ifftshift(fs2);
    recim =ifft2(f2); 

    recerror = mse(recim(:)-im(:));

    subplot(3,3,i);

    imshow(recim);
    title(['k = ',num2str(k),' t = ',num2str(t,5),' zeros = ',num2str(100*sum(fs2(:)==0)/(200*200),3),'%   sqrt(mse)=',num2str(sqrt(recerror),5)]);
    end

    pause

    for i=1:9
      subplot(3,3,i);
      imshow(zeroims{i});
    end
    pause
    
end
