function [f,g,Hfit,Hreg,cft]=LSTKb_obj(rec,fun_params,data)
%LSTKB_OBJ   objective function for deblurring of Gaussian-corrupted images
%with Tikhonov regularization.
%   LSTKB_OBJ evaluates the objective function
%           F(X)=0.5*||A*X-Y||^2+ETA*JR(X)
%   where the data fitting term is the least-squares term ||A*X-Y||^2 
%   (Y is the noisy and blurred image and A is the blurring operator) 
%   and the regularization term JR(X)=0.5*||X||^2 is the Tikhonov 
%   regularization function.
%
%   [F]=LSTKB_OBJ(REC,FUN_PARAMS,DATA) returns the objective function 
%   value F. REC is the current iterate. FUN_PARAMS is a structure array 
%   with the objective function data. DATA is a structure array with the 
%   problem data.
%
%   [F,G]=LSTKB_OBJ(REC,FUN_PARAMS,DATA) also returns the value G of the
%   gradient of the objective function.
%
%   [F,G,HFIT]=LSTKB_OBJ(REC,FUN_PARAMS,DATA) also returns 
%   the scalar HFIT=1.
%
%   [F,G,HFIT,HREG]=LSTKB_OBJ(REC,FUN_PARAMS,DATA) also returns the
%   (NX*NY)-by-(NX*NY) sparse identity matrix HREG which is the hessian of 
%   the regularization term.
%
%   [F,G,HFIT,HREG,CFT]=LSTKB_OBJ(REC,FUN_PARAMS,DATA) also returns the
%   number CFT of performed ffts.

% Get data 
TF = data.TF; % fft of the psf image
gn = data.gn; % blurred and noisy image
eta = fun_params.eta; % regularization parameter
nx = fun_params.nx; ny = fun_params.ny;
cft = 0;

% Evaluate cost functional
resid = real(ifft2(TF.*fft2(rec)))-gn; cft = cft+2;
f = 0.5*resid(:)'*resid(:)+ eta*0.5*rec(:)'*rec(:);

if nargout == 1 
    return
else
    % Evaluate gradient
    g = real(ifft2(conj(TF).*fft2(resid)))+eta*rec; cft = cft+2;
    % Evaluate hessian
    Hreg = speye(nx*ny); 
    Hfit = 1;
end




