algebraic reconstruction technique in matlab

ObtainART.m

function ART_final = ObtainART(sinogram, theta, max_iterations)

disp('Running ART...')

[M0, NN] = size(sinogram);

error_sums = zeros(max_iterations, 1);
last_sum = 0;

reconstruction = zeros(NN);
last_iteration_reconstruction = reconstruction;

for iteration = 1:max_iterations
    for slice_index = 1:M0;
        projection_error = sinogram(slice_index,:) - Project(reconstruction, theta(slice_index));
        reconstruction = reconstruction + imrotate((ones(NN, 1) ./ NN) * projection_error, theta(slice_index), 'bilin', 'crop');
    end

    % Intermediate images:
    %{
    [mystring, errmsg] = sprintf('reconstruction_iteration_%d.png', iteration);
    imwrite(reconstruction, mystring, 'png');    
    %}
    
    % Figure out whether we should stop;
    % Stopping criteria: when the difference between the sums of the elements of 2
    % successive iterations is less than a given number, such as 100.

    error_sum = abs( sum( sum( reconstruction - last_iteration_reconstruction ) ) );
%    error_sums(iteration) = error_sum;
    error_delta = last_sum - error_sum;

   
    [mystring, errmsg] = sprintf('Iteration %d error: %d; error delta: %d', iteration, error_sum, error_delta);
    disp(mystring);    

    if (abs(error_delta) < 100000)
        disp('ART stopping criteria reached.');
        break
    end
    
    last_sum = error_sum;
    last_iteration_reconstruction = reconstruction;    
end

%{
error_sums

DisplayImg(reconstruction, x, y)
title('ART-Reconstructed Phantom');
xlabel('x coordinate');
ylabel('y coordinate');
%}


ART_final = uint8( 255*reconstruction/max(max(reconstruction)) );
[mystring, errmsg] = sprintf('algebraic_reconstruction_technique_%d_projections_%d_iterations.tif', M0, iteration);
imwrite(ART_final , mystring, 'tif');