% 2D tomographic reconstruction via
% Smooth Multiplicative Algebraic Reconstruction Technique (SMART)
% This is a partial implementation. There is no goodness-of-fit, and
% therefore no active process to drive chi^2 to unity.
% 
% CCK 2011-Apr-25
% CCK 2016-Apr-30 -- Removed some old debugging statements to make code readable.
%

function object = smart(projections, angles, noise_model, Niter, a)

% projections = NxM-element array, consisting of M, N-element
%     projections through the NxN-element object,
%     at angles given by the angles array.
% noise_model = a function that can take an array corresponding
%     to expectation values of the projections array, and return
%     the 1-sigma uncertainties on all those values. This is
%     used to calculate a reduced chi-squared.
% Niter = number of iterations.
% object = the result of tomographic inversion, an NxN array.
% a = scalar smoothing parameter, typically >0, <1.
% Also see the forward() function.

[N,M] = size(projections)

object = ones(N,N); %Initial guess (not very enlightened!).

gamma = ones(1,M)/M; %Initialize array of correction strengths.

for i=1:Niter % Cycle through iterations
   object = (1-a)*object + a*smooth2(object,1); %Smooth the object.
   corrections = ones(N,N); %Initialize corrections to unity for this iteration.
   for j=1:M % Cycle through angles (projections).
      [proj, obj_rot] = forward(object, angles(j));
      correction = ( (projections(:,j) ./ max(0.1, proj)) * ones(1,N) );
      % The max operator thresholds at a minimum of 0.1, thus preventing 
      % divide by zero errors. I'm assuming here
      % that I'd never want a value smaller than 0.1 in the denominator.
      % Considering that the units are photons, I think that is a safe
      % assumption.
      correction = imrotate( correction , -angles(j), 'bicubic', 'crop', 0 );
         % That last zero means replace missing data with 0.0 rather than NA.
      correction = max(correction, 0.01); % eliminate negatives.
         % Here i'm assuming a minimum reasonable correction.
      corrections .*= correction .^ gamma(j);
      
   end
   object .*= corrections;
   i %print out the iteration we've just finished.
   correction_mean = mean( vec(corrections) )
   correction_sdev = std( vec(corrections) )
end

