function [f,g,Hfit,Hreg,cft]=LSTVb_obj(rec,fun_params,data)
%LSTVB_OBJ   objective function for deblurring of Gaussian-corrupted images
%with Total Variation regularization.
%   LSTVB_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 JR(x) is the Total variation regularization function.
%
%   [F]=LSTVB_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]=LSTVB_OBJ(REC,FUN_PARAMS,DATA) also returns the value G of the
%   gradient of the objective function.
%
%   [F,G,HFIT]=LSTVB_OBJ(REC,FUN_PARAMS,DATA) also returns 
%   the scalar HFIT=1.
%
%   [F,G,HFIT,HREG]=LSTVB_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]=LSTVB_OBJ(REC,FUN_PARAMS,DATA) also returns the
%   number CFT 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 
TF = data.TF; % fft ofthe psf image
gn = data.gn; % blurred and noisy image
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
resid = real(ifft2(TF.*fft2(rec)))-gn; cft = cft+2;
fls = 0.5*resid(:)'*resid(:);
f = fls+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); 
    gls = real(ifft2(conj(TF).*fft2(resid))); cft = cft+2;
    g = gls+eta*gtv;
    % Evaluate hessian
    Hreg = L; 
    Hfit = 1;
end




