%% Tikhonov regularized reconstruction, example 2
% Compute a Tikhonov regularized reconstruction using
% conjugate gradient algorithm pcg.m
% To run this code, you need your data file and the auxiliary function
% blkmulti.m to be in the same folder.

% Salla Latva-Äijö and Vesa Kaarnioja 2018

% Load the data
% load DataStatic_1120x60 A sinogram or load DataStatic_560x60 A sinogram
load DataStatic_560x60 A sinogram
m = sinogram;
m = m(:);
tt = 80; % number of time steps (80 or 230)
%%
% Initial guess for u
u = zeros(size(A,2)*tt,1);
ss = sqrt(size(u,1)/tt); 
u = reshape(u,[ss,ss,tt]);

alpha = 0.5; % regularization parameter

ht = 1; % adjusts the amount of regularization in temporal direction

% First iteration
graddu2 = [];
graddu2=(8+4/ht)*u-2*([u(2:end,:,:);zeros(1,ss,tt)]+[u(:,2:end,:),...
zeros(ss,1,tt)]+[zeros(1,ss,tt);u(1:end-1,:,:)]+[zeros(ss,1,tt),...
u(:,1:end-1,:)])-(2/ht)*(cat(3,u(:,:,2:end),zeros(ss,ss,1))+...
cat(3,zeros(ss,ss,1),u(:,:,1:end-1)));
disp('Creating gradient function...');
graddu2 = graddu2(:);
u = u(:);

gradF1 = 2*blkmulti(A',blkmulti(A,u))-2*blkmulti(A',m)+alpha*graddu2;%gradF1 = 2*A'*(A*u)-2*A'*m+alpha*graddu2;
disp('Gradient function F1 ready!');
%%
lam = .0001; % first step size

oldu = u;
u = max(0,u-lam*gradF1);

% Iterate 
disp('Iterating... ' );
k = 3; % number of iterations
for l = 1:k 
graddu2 = [];
u = reshape(u,[ss,ss,tt]);
graddu2=(8+4/ht)*u-2*([u(2:end,:,:);zeros(1,ss,tt)]+[u(:,2:end,:),...
zeros(ss,1,tt)]+[zeros(1,ss,tt);u(1:end-1,:,:)]+[zeros(ss,1,tt),...
u(:,1:end-1,:)])-(2/ht)*(cat(3,u(:,:,2:end),zeros(ss,ss,1))+...
cat(3,zeros(ss,ss,1),u(:,:,1:end-1)));
graddu2 = graddu2(:);

u = u(:);
oldgradF1 = gradF1;
gradF1 = 2*blkmulti(A',blkmulti(A,u))-2*blkmulti(A',m)+alpha*graddu2;%gradF1 = 2*A'*(A*u)-2*A'*m+alpha*graddu2;
lam = (u-oldu)'*(u-oldu)/((u-oldu)'*(gradF1-oldgradF1));
oldu = u;
u = max(0,u-lam*gradF1);
disp([num2str(l),'/' num2str(k)]);
end
disp('Iteration ready!');

%% Show reconstruction
u = reshape(u,[ss,ss,tt]);
for ii = 1:tt
imagesc(u(:,:,ii));
title(['Frame ',num2str(ii)]);
axis square;
colormap jet;
%colorbar
pause(0.2);  
end
%% Normalize and watch as a black and white video
% u = u - min(u(:));
% u = u ./ max(u(:));
% implay(u);