% Reconstruct a walnut slice from real X-ray data using Besov B111
% regularization implemented via Haar wavelets.
%
% Samuli Siltanen March 2017

% Parameters for controlling plot appearance
lwidth = 2;

% Number of angles used
Nang = 20;

% Load the measurement matrix and the sinogram from file.
% Discretization of the unknown is NxN.
% Number of projection directions is Nang.
disp(['Loading precomputed data with Nang=', num2str(Nang)])
res  = 2; % This can be 1 or 2 or 4
N    = res*82;
eval(['load data/Data', num2str(N), '_', num2str(Nang), ' A m normA'])
A = A/normA; % Normalize measurement matrix so that its norm is one
m = m/normA; % Normalize measurement as well
disp('Data loaded')

% Set wavelet transform depth
D = log2(2*res);

% Number of iterations
Niter = 1000;

% Threshold mu used in the transform (kind of regularization parameter)
mu = .00007;

% Convolution kernels for the Haar wavelet
h = [ 1 1]/sqrt(2);
g = [-1 1]/sqrt(2);

% Initialize iteration
recn = zeros(N);
dfitvec = repmat(NaN,1,Niter);
sparvec = repmat(NaN,1,Niter);
diffvec = repmat(NaN,1,Niter);

% Perform iteration
for nnn = 1:Niter
    % Record previous iterate
    previous = recn;
    
    % Forward projection of the current iterate
    Arecn = A*recn(:);
    
    % Compare to measured data
    datafit = m(:) - Arecn;
    dfitvec(nnn) = norm(datafit)/norm(m(:));
    
    % Back-projection of the residual
    tmp = A.'*datafit;
    tmp = reshape(tmp,N,N);
    
    % DDD iteration f{n} = S_mu (f{n-1} + A^T(m - A f{n-1}))
    [recn,nonzero,tr,tmp] = Smu_wavelet_oper(h,g,recn+tmp,D,mu);
    sparvec(nnn) = nonzero;
    
    % Projection to non-negative "quadrant"
    recn = max(0,recn);
    
    % Record change in iterate
    diffvec(nnn) = (norm(recn(:)-previous(:)))/norm(previous(:));
    
    if (mod(nnn,30)==1)|(nnn==Niter)
        figure(1)
        clf
        plotim = [scale01(recn),scale01(log(1+abs(tr)))];
        plotim = plotim.^(.5); % gamma correction for image shades
        imagesc(plotim)
        axis equal
        axis off
        colormap gray
        title(['Angles ', num2str(Nang), ', iteration number ', num2str(nnn)],'fontsize',16)
        
        % Monitor crucial parameters of the process
        figure(2)
        clf
        subplot(3,1,1)
        plot(1:Niter,sparvec,'linewidth',lwidth)
        xlim([1 Niter])
        title('Ratio of nonzero coefficients')
        subplot(3,1,2)
        semilogy(1:Niter,diffvec,'linewidth',lwidth)
        xlim([1 Niter])
        title('Relative square norm change in iterate')
        subplot(3,1,3)
        semilogy(1:Niter,dfitvec,'linewidth',lwidth)
        xlim([1 Niter])
        title('Data fit')
        drawnow
                
    end
    disp([nnn Niter])
end

% Save results to file
eval(['save data/recon_Haar_', num2str(N), '_', num2str(Nang), ' recn Niter mu D sparvec diffvec dfitvec'])

