imlist = cell(4);
imlist{1} = 'im1.jpg';
imlist{2} = 'im2.jpg';
imlist{3} = 'im3.jpg';
imlist{4} = 'im4.jpg';

for i=4:4
    x = imread(imlist{i});
    %x = rgb2gray(x);
    [r,c,u]=size(x);
    rc = r*c;
    rcu = rc*u;
    X = zeros(size(x));
    for k=1:u;
        X(:,:,k) = fft2(x(:,:,k));
    end
    [Z,I] = sort(abs(X(:)),'descend');
    P = 10*log10(abs(X));
    P = P - min(min(min(P)));
    P = P / max(max(max(P)));
    figure(1);
    subplot(221);
    imagesc(P); axis image;
    subplot(212);
    plot(10*log10(Z));
    grid on;
    hold on;
    ax = axis;
    han = plot([0 0],ax(3:4),'r');
    hold off;

    subplot(222);
    for m = 1:100:rcu
       set(han,'XData',[m m]);
       Y = X;
       Y(abs(Y)<Z(m)) = 0;
       y = zeros(size(x));
       for k=1:u
           y(:,:,k) = real(ifft2(Y(:,:,k)));
       end
       imagesc(real(y)/256); axis image;
       title(int2str(m));
       drawnow;
       pause;
    end
end

return;

%    kr = [0:r-1].';
%    kc = [0:c-1];
%    
%    y = zeros(size(x));
%    for m = 1:rcu
%        p = I(m)-1;
%        uu = floor(p/(rc));
%        pp = mod(p,rc);
%        ss = mod(pp,r);   % row index
%        tt = floor(pp/r); % col index
%        W = exp(1j*2*pi*(bsxfun(@plus,kr*ss/r,kc*tt/c)))/(rc);
%        y(:,:,uu+1) = y(:,:,uu+1) + X(I(m))*W;
%        
%        if(mod(m,500)==0)
%        subplot(211);
%        imagesc(real(y)); axis image;
%        subplot(212);
%        imagesc(real(W)); axis image;
%        title(int2str(m));
%        drawnow;
%        %pause;
%        end
%        
%    end
%end
