% Here we compute the SVD of the sparse tomographic measurement matrix A
% using full matrix algorithms. This is not very efficient coding, and it
% is done only for educational purposes.
%
% The following routines must be run before this file:
% XRM5_SVD_comp.m and XRM3_NoCrimeData_comp.m.
%
% Samuli Siltanen October 2012

% Gamma correction value for making the plots brighter
gammacorr = .4;

% Choose resolution
N = 50;

% Plot parameters
fsize  = 16;
msize  = 20;
lwidth = 4;

% Load noisy measurements from disc.
eval(['load data/XRMC_NoCrime', num2str(N), ' N mnc mncn']);
mn = mncn;

% Load singular value decomposition of the measurement matrix
eval(['load data/XRME_SVD', num2str(N), ' U D V A measang target N P Nang']);

% Demonstration of inversion by truncated SVD with several singular vectors
STEP  = floor((N^2-1)/14);
Nsvec = [1:STEP:N^2];
for iii = 1:length(Nsvec)
    
    % Reconstruct from noisy data using only large singular values.
    Ns               = Nsvec(iii);
    [row,col]        = size(D.');
    Dplus            = sparse(row,col);
    svals            = diag(D);
    Dplus(1:Ns,1:Ns) = diag(1./svals(1:Ns));
    recn             = V*Dplus*U.'*mn(:);
    recn             = reshape(recn,N,N);
    
    % Calculate relative error of reconstruction
    relerr = round(norm(recn(:)-target(:))/norm(target(:))*100);

    % Construct the latest singular vector in image form
    singvec = reshape(V(:,Ns),N,N);
    
    % Normalize the values of the pixel images singvec and recn so that
    % the pixel values are between zero and one.
    singvec  = singvec-min(min(singvec));
    singvec  = singvec/max(max(singvec));
    plotrecn = recn;
    plotrecn = max(0,plotrecn);
    plotrecn = plotrecn/max(max(plotrecn));
    
    
    % Take a look at the reconstruction
    figure(1)
    clf
    subplot(2,2,1)
    imagesc(singvec)
    colormap gray
    axis square
    axis off
    title('Singular vector','fontsize',fsize)
    subplot(2,2,2)
    semilogy([1:length(svals)],svals,'r','linewidth',lwidth)
    hold on
    semilogy([1:Ns],svals(1:Ns),'b','linewidth',lwidth)
    semilogy([Ns Ns],[svals(end) svals(Ns)],'k--','linewidth',lwidth/2)
    semilogy([Ns],[svals(Ns)],'b.','markersize',msize)
    axis([1 length(svals) svals(end) svals(1)])
    set(gca,'xtick',[Ns],'fontsize',fsize)
    set(gca,'ytick',linspace(min(svals),max(svals),20),'fontsize',fsize)
    set(gca,'yticklabel',{},'fontsize',fsize)
    axis square
    box off
    title('Singular values','fontsize',fsize)
    subplot(2,2,3)
    eval(['load data/XRMD_', num2str(N), ' target']);
    imagesc(max(0,target).^gammacorr)
    colormap gray
    axis square
    axis off
    title('Original phantom','fontsize',fsize)    
    subplot(2,2,4)
    imagesc(plotrecn.^gammacorr)
    colormap gray
    axis square
    axis off
    title(['Error ', num2str(relerr),'%'],'fontsize',fsize)
    pause
end
