% to start the program write : learndft('saturn.tif',150);
% the program will then use the first 150 coefficient of the fourier transform
% of the saturn image. you can use any othe image.
% if you want the images to run faster or slower change pause on line 44 
function recon = learndft(image,coef) 

    image=imread(image);
    
    [m,n] = size(image);
    
    if mod(m,2)
        m=m-1;
    end
    
    if mod(n,2) 
        n=n-1;
    end

    image=image(1:m,1:n);
    p=pairs(image,coef);
    
    IMAGE = fftshift(fft2(image));
    recon = zeros(m,n);              
     
    figure(1), clf
    subplot(2,2,2),imagesc(image),colormap gray, axis image,axis off;%, title('image');
    subplot(2,2,1),fftshow(image);
     
    
    for N = 1:coef
	  
	  F = zeros(m,n);
	  F(p(N,1), p(N,2)) = IMAGE(p(N,1), p(N,2));
	  F(p(N,3), p(N,4)) = IMAGE(p(N,3), p(N,4));	  
	  
	  f = real(ifft2(fftshift(F)));
	  recon = recon+f;	  
	 
	  subplot(2,2,1),fftshow(image);
	  hold on, plot([p(N,2), p(N,4)], [p(N,1), p(N,3)],'k.');hold off;
      subplot(2,2,4),imagesc(recon),colormap gray, axis image, axis off;%, title('Reconstruction');
	  %subplot(2,2,4),surf(f(1:8:m,1:8:n)),colormap gray, axis image, axis off
	  subplot(2,2,3),imagesc(f),colormap gray, axis image, axis off
      pause(0.001);
  end
 
  for N=1:coef
      subplot(2,2,1)%,fftshow(image);
	  hold on, plot([p(N,2), p(N,4)], [p(N,1), p(N,3)],'y.');
  end
  subplot(2,2,1)%,fftshow(image);
	  hold on, plot([p(N,2), p(N,4)], [p(N,1), p(N,3)],'k.');
  
      
  
  
function M=pairs(image,coef)

[m n]=size(image);
x1=m/2 + 1;
x2=m/2 + 1;
y1=n/2;
y2=n/2+2;
count=1;

img = mod(1,2);
img = mod(2,2);
    
for loop=1:m/2-1
    i1=x1;
    i2=x2;
    j1=y1;
    j2=y2;
    
    if mod(loop,2)
        
        if(count <= coef)
            for i1=x1:(x1+2*loop - 1)
                M(count,:)=[i1 j1 i2 j2];
                i2=i2-1;
                count=count+1;
            end
        end
        i2=i2+1;
        y1=j1+1;
        j2=j2-1;
        
        if(count <= coef)
            for j1=y1:(y1+2*loop-1)
                M(count,:)=[i1 j1 i2 j2];
                j2=j2-1;
                count=count+1;
            end
        end
        x1=i1;
        y1=j1+1;
        x2=i2;
        y2=j2;
        
    else
        
        if(count <= coef)
            for i1=x1:-1:(x1-2*loop + 1)
                M(count,:)=[i1 j1 i2 j2];
                i2=i2+1;
                count=count+1;
            end
        end
        i2=i2-1;
        y1=j1-1;
        j2=j2+1;
        
        if(count <= coef)
            for j1=y1:-1:(y1-2*loop+1)
                M(count,:)=[i1 j1 i2 j2];
                j2=j2+1;
                count=count+1;
            end
        end
        x1=i1;
        y1=j1-1;
        x2=i2;
        y2=j2;
        
    end    
end

function IMAGE=fftshow(image)

IMAGE = fft2(image);
imagesc(fftshift(log(abs(IMAGE)))),colormap gray, axis image
    axis off;%,   title('fftshift(log(abs(fft2(image))))','FontSize',10);