% Load the measurement matrix and the sinogram from
% file DataStatic_128x60.mat
load DataStatic_128x60 A sinogram

m = sinogram;
% Compute a Tikhonov regularized reconstruction using
% conjugate gradient algorithm pcg.m

% Define the number of slices
T = 33;
N     = sqrt(size(A,2)/T);
alpha = 10; % regularization parameter
fun   = @(x) A.'*(A*x)+alpha*x;
b     = A.'*m(:);
x     = pcg(fun,b);

% Compute a Tikhonov regularized reconstruction from only
% 10 projections
[mm,nn] = size(m);
ind     = [];
for iii=1:nn/6
    ind = [ind,(1:mm)+(6*iii-6)*mm];
end
m2    = m(:,1:6:end);
A     = A(ind,:);
alpha = 10; % regularization parameter
fun   = @(x) A.'*(A*x)+alpha*x;
b     = A.'*m2(:);
x2    = pcg(fun,b);

% Take a look at the sinograms of the first slice and the reconstructions
figure
subplot(2,2,1)
imagesc(m(:,1:size(m,2)/T))
colormap gray
axis square
axis off
title('Sinogram, 60 projections')
subplot(2,2,3)
imagesc(m2(:,1:size(m2,2)/T))
colormap gray
axis square
axis off
title('Sinogram, 10 projections')
% reshape the reconstruction
x = reshape(x,N,N,T);
x2 = reshape(x2,N,N,T);
% Show the first slice reconstruction
subplot(2,2,2)
imagesc(imrotate(x(:,:,1),-98,'bilinear','crop'))
colormap gray
axis square
axis off
title({'Tikhonov reconstruction,'; '60 projections'})
subplot(2,2,4)
imagesc(imrotate(x2(:,:,1),-98,'bilinear','crop'))
colormap gray
axis square
axis off
title({'Tikhonov reconstruction,'; '10 projections'})