% creating the test problem
clear, clc, close all
startup
load AtmosphericBlur30
A = psfMatrix(PSF, center,'zero');
b = A*f_true;
sigma = 20;
beta = 60;
gn = normrnd(0,sigma,256,256);
pn = poissrnd(beta,256,256);
ber = b+gn+pn;
b0 = ber; b0(b0<0)=0;
level = norm(ber(:)-b(:))/norm(b(:));
xex = f_true(:);
%
% running different solvers
m_in = 20; m_out = 50;
b0 = b0(:);
display('beginning of CP-MFCGLS...')
stopC1.kind = 'stab';
stopC1.tau = 1e-4;
[Xiter_P, RelErr_P, RelRes_P, Alpha_P, tot_it_P, cycle_it_P, stopC_itP] = fcglsPoiss_rp_nn_ReSt (A, ber(:), m_in, m_out, sigma^2, beta, xex, b0, 'off', stopC1);
display('beginning of CP-MFCGLS(k)...')
[Xiter_VP, RelErr_VP, RelRes_VP, Alpha_VP, tot_it_VP, cycle_it_VP, stopC_itVP] = fcglsPoiss_rp_nn_ReSt (A, ber(:), m_in, m_out, sigma^2, beta, xex, b0, 'on', stopC1);
options = IRset('x_true', xex, 'StepLength', -1, 'MaxIter', 300, 'x0', b0);
display('beginning of WMRNSD...')
[X_WMRNSD, IterInfo_WMRNSD]  = IRwmrnsd(A, ber(:), options, sigma^2, beta);
display('beginning of KWMRNSD...')
[X_KWMRNSD, IterInfo_KWMRNSD]  = IRkwmrnsd(A, ber(:), options, sigma^2, beta);
%
% plotting the results (Relative Errors History),
figure, semilogy(RelErr_P)
hold on
semilogy(RelErr_VP, '-c')
semilogy(IterInfo_WMRNSD.Enrm, '-r')
semilogy(IterInfo_KWMRNSD.Enrm, '-m')
title('relative error history, satellite test problem, Gaussian and Poisson noise')
legend('CP-MFCGLS', 'CP-MFCGLS(k)', 'WMRNSD', 'KWMRNSD')
%
% showing the best restored images
figure, imshow(reshape(xex,256,256),[]), title('exact image')
figure, imshow(reshape(ber,256,256),[]), title('corrupted image')
figure, imshow(reshape(Xiter_P(:,end),256,256),[]), title('CP-MFCGLS')
figure, imshow(reshape(Xiter_VP(:,end),256,256),[]), title('CP-MFCGLS(k)')
figure, imshow(reshape(X_WMRNSD,256,256),[]), title('WMRNSD')
figure, imshow(reshape(X_KWMRNSD,256,256),[]), title('KWMRNSD')
