% Driver for the NP or QNP methods

clear all
close all

% --------------------------------------------
% Set up the image restoration problem
% --------------------------------------------
fprintf('Set up the noise type affecting the data: \n');
model_choice.noise = ...
    input('  enter p for Poisson noise and g for Gaussian noise ','s');
fprintf('\n-------------------------------------------------\n');
if strcmp(model_choice.noise,'p') 
        fprintf(' QNP method for deblurring Poisson data \n');
        fprintf(' with Tikhonov regularization \n');
elseif model_choice.noise == 'g'  
        bg = [];
        fprintf(' NP method for deblurring Gaussian data \n');
        fprintf(' with Tikhonov regularization \n');
else
        fprintf('Error in entering the type of noise \n'); return
end
fprintf('-------------------------------------------------\n\n');

% --------------------------------------------
% Set up the data problem
% --------------------------------------------
data_file = input('Enter the data .mat file ','s'); 
load(data_file); 
TF = fft2(fftshift(psf));

% --------------------------------------------
% Set up the regularization parameter
% --------------------------------------------
fun_params.eta = input('Enter the regularization parameter value ');
fprintf('regularization parameter = %e \n',fun_params.eta);

% --------------------------------------------
% Set up the initial iterate
% --------------------------------------------
fprintf('Set up the initial iterate: \n');
fprintf('  enter d for the default initial iterate \n');
fprintf('  enter p for projected aquired image \n');
fprintf('  or enter n \n');
rec0_choice = input('','s');
if strcmp(rec0_choice,'d')
    rec0 = ones(size(gn))*sum(gn(:))/numel(gn);
elseif strcmp(rec0_choice,'p')
    rec0 = max(gn,0);
elseif strcmp(rec0_choice,'n')
    rec0 = input('    Enter the initial iterate ');
else
    fprintf('Error in entering the initial iterate \n'); return
end
   
% Save data in structure array 
data.obj = obj; % exact image
data.gn = gn; % blurred and noisy image
data.bg = bg; % background
data.TF = TF; % fft of the psf image
data.rec0 = rec0; % initial guess

% --------------------------------------------
% Set up the method parameters
% --------------------------------------------
fprintf('Set up the stopping criteria parameters: \n');
fprintf('  enter the maximum number of NCGP iterations \n');
fprintf('  or enter 0 for the default value (20) \n');
params_choice.kmax = input('');
if params_choice.kmax == 0
    params.kmax = 20;          
else
    params.kmax = params_choice.kmax;
end
fprintf('maximum number of NCGP iterations = %d \n\n',params.kmax);

fprintf('  enter the projected gradient tolerance \n');
fprintf('  or enter 0 for the default value (1e-4) \n') 
params_choice.tol_gp = input('');
if params_choice.tol_gp == 0
    params.tol_gp = 1e-4;          
else
    params.tol_gp = params_choice.tol_gp;
end
fprintf('projected gradient tolerance = %e \n\n',params.tol_gp);

fprintf('  enter the objective function tolerance \n');
fprintf('  or enter 0 for the default value (1e-4) \n');                   
params_choice.tol_f = input('');
if params_choice.tol_f == 0
    params.tol_f = 1e-4;          
else
    params.tol_f = params_choice.tol_f;
end        
fprintf('objective function tolerance = %e \n\n',params.tol_f);

params.max_fft = 1e3; % maximum number of FFTs

if strcmp(model_choice.noise,'p') 
        % Quasi-Newton Projection method
        [x,histout] = QNP(data,fun_params,params);
elseif strcmp(model_choice.noise,'g')  
        % Newton Projection method
        [x,histout] = NP(data,fun_params,params);
end

% --------------------------------------------
% Print numerical results and visualize images
% --------------------------------------------
NPoutput(x,histout,data)
