% Driver for the NCGP method

clear all
close all

% --------------------------------------------
% Set up the image restoration problem
% --------------------------------------------
fprintf('Set up the image restoration problem: \n');
model_choice.pb = ...
    input(' enter b for Deblurring and n for Denoising ','s');
% deblurring
if strcmp(model_choice.pb,'b')   
    model_choice.noise = ...
        input('  Enter p for Poisson noise and g for Gaussian noise ','s');
    model_choice.reg = ...
        input('  Enter tv for Total Variation and tk for Tikhonov ','s');
    if strcmp(model_choice.noise,'p')
        if strcmp(model_choice.reg,'tv')
            % Set objective and constraint functions 
            % and save in structure fun_params
            fun_params.obj = 'KLTVb_obj'; % objective function
            fun_params.Amult = 'Amult'; % matrix-vector product with 
            % the Hessian of the objective function 
        elseif strcmp(model_choice.reg,'tk')
            % Set objective and constraint functions 
            % and save in structure fun_params
            fun_params.obj = 'KLTKb_obj'; % objective function
            fun_params.Amult = 'Amult'; % matrix-vector product with 
            % the Hessian of the objective function
        else
            fprintf('Error in entering the regularization term \n'); return
        end
    elseif strcmp(model_choice.noise,'g')
        bg = [];
        if strcmp(model_choice.reg,'tv')
            % Set objective and constraint functions 
            % and save in structure fun_params
            fun_params.obj  = 'LSTVb_obj'; % objective function
            fun_params.Amult = 'Amult'; % matrix-vector product with 
            % the Hessian of the objective function  
        elseif strcmp(model_choice.reg,'tk')
            % Set objective and constraint functions 
            % and save in structure fun_params
            fun_params.obj = 'LSTKb_obj'; % objective function
            fun_params.Amult = 'Amult'; % matrix-vector product with 
            % the Hessian of the objective function
        else
            fprintf('Error in entering the regularization term \n'); return
        end
    else
        fprintf('Error in entering the type of noise \n'); return
    end    
elseif strcmp(model_choice.pb,'n')    
    psf = [];
    model_choice.noise = ...
        input(' Enter p for Poisson noise and g for Gaussian noise ','s');
    if strcmp(model_choice.noise,'p')
        % Set objective and constraint functions 
        % and save in structure fun_params
        fun_params.obj = 'KLTVn_obj'; % objective function
        fun_params.Amult = 'Amult'; % matrix-vector product with 
        % the Hessian of the objective function  
    elseif strcmp(model_choice.noise,'g')
        bg = [];
        % Set objective and constraint functions 
        % and save in structure fun_params
        fun_params.obj = 'LSTVn_obj';   % objective function
        fun_params.Amult = 'Amult';  % matrix-vector product with 
        % the Hessian of the objective function  
    else
        fprintf('Error in entering the type of noise \n'); return
    end
else
    fprintf('Error in entering the image restoration model \n'); return      
end

% --------------------------------------------
% 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);

fprintf('Set up the CG parameters: \n');
fprintf('  enter the maximum number of CG iterations \n');
fprintf('  or enter 0 for the default value (30) \n');
params_choice.pcg_max = input('');
if params_choice.pcg_max == 0
    params.pcg_max = 30;          
else
    params.pcg_max = params_choice.pcg_max;
end
fprintf('maximum number of CG iterations = %d \n\n',params.pcg_max);

fprintf('  enter the CG relative tolerance \n');
fprintf('  or enter 0 for the default value (0.1) \n'); 
params_choice.pcg_tol = input('');
if params_choice.pcg_tol == 0
    params.pcg_tol = 0.1;          
else
    params.pcg_tol = params_choice.pcg_tol;
end
fprintf('CG relative tolerance = %d \n\n',params.pcg_tol);

fprintf('  enter the 1 for CG diagonal preconditioner \n');
fprintf('  or enter 0 \n');
params_choice.precond = input('');
if params_choice.precond==0 || params_coice.precond==1
    params.precond  = params_choice.precond;      
else
    fprintf('Error in the choice of CG preconditioner \n')
    return;
end

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

% --------------------------------------------
% Newton-Cg Projection method
% --------------------------------------------
[x,histout] = NCGP(data,fun_params,params);

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


        

    
        


