% 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
%

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.
      if (max(max( abs(imag(object)) )) != 0)
         'smoothing messed with my object'
         error(0)
      end
   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));
      if (max(max( abs(imag(object)) )) != 0)
         'error in object here'
         error(0)
      end
      if (max(max( abs(imag(proj)) )) != 0)
         'error in proj here'
         error(0)
      end
      if (max(max( abs(imag(obj_rot)) )) != 0)
         'error in obj_rot here'
         error(0)
      end
      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.
      if (max(max( abs(imag(correction)) )) != 0)
         'no, error here'
         error(0)
      end
      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.
      if (max(max( abs(imag(correction)) )) != 0)
         'error right here'
         error(0)
      end
      if (max(max( abs(imag(correction)) )) != 0)
         'the correction is already screwed up here!!!'
         i
         j
         error(0)
      end
      if (max(max( abs(imag(corrections)) )) != 0)
         'WTF! corrections already screwed up here!!!'
         i
         j
         error(0)
      end
      pre_error1 = correction;
      pre_error2 = gamma;
      pre_error3 = gamma(j);
      pre_error4 = corrections;
      corrections .*= correction .^ gamma(j);
      %errorses = max(max( abs(imag(corrections)) ))
      if (max(max( abs(imag(corrections)) )) != 0)
         'my corrections are already screwed up here.'
         i
         j
         gamma
         gamma(j)
         gotcha1 = pre_error1(29,2) .^ pre_error3
         gotcha2 = pre_error1 .^ pre_error3
         pre_error1
         pre_error2
         pre_error3
         pre_error4
         error(0)
      end
      
   end
   object .*= corrections;
   i %print out the iteration we've just finished.
   correction_mean = mean( vec(corrections) )
   correction_sdev = std( vec(corrections) )
end

