clear, clc, close all
startup
N = 256; 
[A,b_ex,x_ex,theta,p,d] = paralleltomo(N,0:2:179,round(sqrt(2)*N),sqrt(2)*N);
b = addRELnoise(b_ex, 5*10^-2);
% Perform the Cimmino iterations
it_Cimm = 120;
K = 1:it_Cimm;
options.nonneg = true;
display('beginning of Cimmino...')
[X_Cimm, Info_Cimm] = cimmino(A, b, K);
RelErr_Cimm = zeros(it_Cimm, 1);
nx = norm(x_ex(:));
for i = 1:it_Cimm
    RelErr_Cimm(i) = norm(x_ex(:)-X_Cimm(:,i))/nx;
end
% Perform the Restarted FCGLS iterations
it_in = 10; it_out = 20;
stopC_1.kind = 'discr';
stopC_1.noise = 5*10^-2;
stopC_2.kind = 'stab';
stopC_2.tau = 1e-2;
%
display('beginning of MFCGLS - first test..')
[Xiter_FCGRe1, RelErr_FCGRe1, RelRes_FCGRe1, Alpha_FCGRe1, tot_it1, cycle_it1, stopC_it1] = fcgls_rp_nn_ReSt (A, b(:), it_in, it_out, x_ex(:), zeros(N^2,1), 1, stopC_1);
display('beginning of MFCGLS - second test..')
[Xiter_FCGRe2, RelErr_FCGRe2, RelRes_FCGRe2, Alpha_FCGRe2, tot_it2, cycle_it2, stopC_it2] = fcgls_rp_nn_ReSt (A, b(:), it_in, it_out, x_ex(:), zeros(N^2,1), 1, stopC_2);
display('beginning of MRNSD...')
options = IRset('x_true', x_ex(:), 'MaxIter', 200, 'x0', 1e-4*ones(size(x_ex(:))));
[X_MRNSD, IterInfo_MRNSD] = IRmrnsd(A, b(:), options);
display('beginning of NN SD...')
options = IRset(options, 'StepLength', -1);
[X_SDNN, IterInfo_SDNN] = IRgdnn(A, b(:), options);
%
% plotting the results (Relative Errors History)
figure, semilogy(RelErr_FCGRe1)
hold on
semilogy(RelErr_Cimm, '-g')
semilogy(IterInfo_MRNSD.Enrm, '-r')
semilogy(IterInfo_SDNN.Enrm, '-c')
title('relative error history, underdetermined tomographic reconstruction problem')
legend('MFCGLS', 'Cimmino', 'MRNSD', 'SDNN')
%
% showing the best reconstructed phantoms
[minvC, indvC] = min(RelErr_Cimm);
[minvF, indvF] = min(RelErr_FCGRe1);
figure, imshow(reshape(x_ex,256,256),[]), title('exact phantom')
figure, imshow(reshape(X_Cimm(:,indvC),256,256),[]), title('Cimmino')
figure, imshow(reshape(Xiter_FCGRe1(:,indvF),256,256),[]), title('MFCGLS')
figure, imshow(reshape(X_MRNSD,256,256),[]), title('MRNSD')
figure, imshow(reshape(X_SDNN,256,256),[]), title('SDNN')
