% Main program running the ART reconstruction algorithm
%----------------------------------------------------------------------------
% CALLS TO : Ibeam.m, mktomomat.m, ART.m
% 29/2/96 Erkki Somersalo
%----------------------------------------------------------------------------

% Computing the original density distribution matrix

Ibeam

N = length(X(1,:));  % Number of pixels in vertical/horizontal direction
M = 64;              % Number of illumination directions
K = 64;              % Number of rays in each direction

% Computing the matrix relating pixel values to data (Patience; this
% takes a while!)

A = mktomomat(N,M,K);

% Computing the tomographic data (Warning: The data is computed using the
% same matrix that is later used in reconstruction. The results may be
% misleadingly optimistic!)

xx = X'; xx = xx(:);
y = A*xx;
Y = reshape(y,M,K);

figure(2)
pcolor(Y'), colormap(1-gray),shading flat,axis('square'), axis('off'), colorbar

% Reconstruction with ART

x = zeros(N*N,1);
w = 1;
nrounds = 10; 
tol = 0.5;
count = 0; change = 2*tol;
while change > tol  
  % Iterate while the norm of the image changes more than tol; after 
  % nrounds iterations show the image
  
  x0 = x;
  x = ART(A,y,x0,w,1);
  change = norm(x - x0);
  count = count + 1
  if count == nrounds
    Xrec = reshape(x,N,N)';
    figure(3)
    pcolor(Xrec), colormap(1-gray),shading flat,axis('square'), axis('off')
    colorbar
    pause(1)
    count = 0;
  end  
end
