% Reconstruct a lotus slice from real X-ray data using Landweber iteration.
%
% Samuli Siltanen September 2016

% Load the measurement matrix and the sinogram from
load LotusData128
N = 128;
A = A/normA; % Normalize measurement matrix so that its norm is one
m = m/normA; % Normalize measurement as well
disp('Data loaded')

% Number of iterations
Niter = 500;

% Landweber iteration relaxation factor
omega = .8;

% Initialize iteration
recn = zeros(N);

% Perform iteration
for nnn = 1:Niter
    
    % Forward projection of the current iterate
    Arecn = A*recn(:);
    
    % Back-projection of the residual
    tmp = A.'*(m(:) - Arecn);
    tmp = reshape(tmp,N,N);
    
    % Landweber iteration
    recn = recn + omega*tmp;
    
    % Projection to non-negative "quadrant"
    recn = max(0,recn);
    
    if (mod(nnn,37)==0)|(nnn==Niter)
        figure(1)
        clf
        plotim = scale01(recn);
        plotim = plotim.^(.5); % gamma correction for better image shades
        imagesc(plotim)
        axis equal
        axis off
        colormap gray
        title(['Iteration number ', num2str(nnn)],'fontsize',16)
        
        %pause
        %keyboard
    end
end

