function [f,g,Hfit,Hreg,cft]=KLTVn_obj(rec,fun_params,data)
%KLTVN_OBJ   objective function for denoising of Poisson-corrupted images 
%with Total Variation regularization.
%   KLTVN_OBJ evaluates the objective function
%           F=KL(X,Y)+eta*JR(X)
%   where KL(X,Y) is the Kullback-Leibler divergence of X from Y:
%           KL(X,Y)= sum(sum(Y.*log(Y./(X+BG))+X+BG-Y))
%   (Y is the noisy image and BG is the background value) 
%   and JR(x) is the Total Variation regularization function.
%
%   [F]=KLTVN_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]=KLTVN_OBJ(REC,FUN_PARAMS,DATA) also returns the value G of the
%   gradient of the objective function.
%
%   [F,G,HFIT]=KLTVN_OBJ(REC,FUN_PARAMS,DATA) also returns the NX-by-NY
%   matrix HFIT such that the Hessian of KL is diag(HFIT(:)).
%
%   [F,G,HFIT,HREG]=KLTVN_OBJ(REC,FUN_PARAMS,DATA) also returns the 
%   (NX*NY)-by-(NX*NY) first order approximation HREG of the hessian 
%   of the Total Variation function.
%
%   [F,G,HFIT,HREG,CFT]=KLTVN_OBJ(REC,FUN_PARAMS,DATA) also returns the
%   number CFT=0 of performed ffts.
%
%   See also SetTV.

%   The TV discretization is obtained as described in:
%   C. Vogel, "Computational Methods for Inverse Problems", SIAM, 2002.
%   The first order approximation of the TV hessian is obtained as 
%   described in:
%   C. Vogel and M.E. Oman, "Iterative methods for total variation
%   denoising", SIAM J. Sci. Comput, 17, pp. 227-238, 1996.

% Get data 
gn = data.gn; % noisy image
bg = data.bg; % background
eta = fun_params.eta; % regularization parameter
beta = fun_params.betatv; % smoothing parameter for tv
cft=0;

% TV discretization
Dx = fun_params.Dx; Dy = fun_params.Dy;
nx = fun_params.nx; ny = fun_params.ny;
hx = fun_params.hx; hy = fun_params.hy; 
hxy = fun_params.hxy;

% Evaluate TV function
Dxu = reshape(Dx*rec(:),nx-1,ny-1);
Dyu = reshape(Dy*rec(:),nx-1,ny-1); 
phi = sqrt(Dxu.^2+Dyu.^2+beta^2); 
ftv = sum(phi(:))*hxy;

% Evaluate cost function
den = rec+bg;
arg = gn./den; arg(gn==0) = 1;
fkl = sum(sum(gn.*log(arg)+den-gn));
f = fkl+eta*ftv;

if nargout == 1 
    return
else
    % Evaluate gradient
    dphi = 1./phi;
    Q = spdiags(dphi(:),0,(nx-1)*(ny-1),(nx-1)*(ny-1));
    L = (Dx'*Q*Dx+Dy'*Q*Dy)*hxy;
    gtv = reshape(L*rec(:),nx,ny); 
    temp = gn./den;
    gkl = -temp+1 ;  
    g = gkl+eta*gtv;
    % Evaluate hessian
    Hreg = L;
    Hfit = temp./(den);
end




