function [f,g,Hfit,Hreg,cft] = KLTKb_obj(rec,fun_params,data)
%KLTKB_OBJ   objective function for deblurring of Poisson-corrupted images 
%with Tikhonov regularization.
%   KLTKB_OBJ evaluates the objective function
%           F(X)=KL(X,A,Y)+ETA*JR(X)
%   where KL(X,A,Y) is the Kullback-Leibler divergence of X from Y:
%           KL(X,A,Y)= sum(sum(Y.*log(Y./(A*X+BG))+A*X+BG-Y))
%   (Y is the noisy and blurred image, A is the blurring operator, 
%   BG is the background value) and the regularization term
%   JR(X)=0.5*||X||^2 is the Tikhonov regularization function.
%
%   [F]=KLTKB_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]=KLTKB_OBJ(REC,FUN_PARAMS,DATA) also returns the value G of the
%   gradient of the objective function.
%
%   [F,G,HFIT]=KLTKB_OBJ(REC,FUN_PARAMS,DATA) also returns the NX-by-NY
%   matrix HFIT such that the Hessian of KL is A'*diag(HFIT(:))*A.
%
%   [F,G,HFIT,HREG]=KLTKB_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]=KLTKB_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
bg = data.bg; % background
eta = fun_params.eta; % regularization parameter
nx = fun_params.nx; ny = fun_params.ny;
cft = 0;

% Evaluate KL functional
den = real(ifft2(TF.*fft2(rec)))+bg; cft = cft+2;
fkl = sum(sum(gn.*log(gn./den)+den-gn));

% Evaluate cost function
f = fkl+eta*0.5*rec(:)'*rec(:);

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




