function [f,g,Hfit,Hreg,cft]=LSTVn_obj(rec,fun_params,data)
%LSTVN_OBJ   objective function for denoising of Gaussian-corrupted images 
%with Total Variation regularization.
%   LSTVN_OBJ evaluates the objective function
%           F(X)=0.5*||X-Y||^2+ETA*JR(X)
%   where the data fitting term is the least-squares term ||X-Y||^2 
%   (Y is the noisy image) and JR(X) is the Total Variation 
%   regularization function.
%
%   [F]=LSTVN_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]=LSTVN_OBJ(REC,FUN_PARAMS,DATA) also returns the value G of the
%   gradient of the objective function.
%
%   [F,G,HFIT]=LSTVN_OBJ(REC,FUN_PARAMS,DATA) also returns 
%   the scalar HFIT=1.
%
%   [F,G,HFIT,HREG]=LSTVN_OBJ(REC,FUN_PARAMS,DATA) also returns the
%   (NX*NY)-by-(NX*NY) first order approximation HREG of the hessian 
%   of the Total Variation functional.
%
%   [F,G,HFIT,HREG,CFT]=LSTVN_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
eta = fun_params.eta; % regulariztion parameter
beta = fun_params.betatv; % smoothing paramter 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 = rec-gn;
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 = resid ;  
    g = gls+eta*gtv;
    % Evaluate hessian
    Hreg = L; 
    Hfit = 1;
end




