function [f,grad,Hz]=tvmz(X,eps,Z,par)
%Smoothed Total Variation (TV) of image X, the gradient of TV and Hess-vector product
%
% Michael Zibulevsky,  14.06.2006



D1X= D1(X);
D2X= D2(X);

phi=sqrt(D1X.^2+D2X.^2 +eps);

f=sum(phi(:)) - length(phi(:))*sqrt(eps);

if nargout>1,
	tmp=1./phi;
	dphi1=D1X.*tmp;
	dphi2=D2X.*tmp;
	grad = D1adj(dphi1)+D2adj(dphi2);
end

if nargout > 2,           % Compute Hessian-vector product Hz
	tmp2=tmp.^2;
	d2phi11=tmp - D1X.*dphi1.*tmp2;
	d2phi22=tmp - D2X.*dphi2.*tmp2;
	d2phi12=       - D1X.*dphi2.*tmp2;

  K=size(Z,3); % For the case of multiple images Z
  for k=1:K
    D1z=D1(Z(:,:,k));
    D2z=D2(Z(:,:,k));
    Hz(:,:,k)=D1adj(d2phi11.*D1z + d2phi12.*D2z) + D2adj(d2phi22.*D2z + d2phi12.*D1z);
  end

end

function Y=D1(X)
[m,n]=size(X);
Y=[diff(X,1,1);zeros(1,n)];
	
function Y=D2(X)
[m,n]=size(X);
Y=[diff(X,1,2) zeros(m,1)];
	
	
function  Z=D1adj(Y)
Z=[-Y(1,:); -diff(Y,1,1)];  Z(end,:)=Z(end,:)+Y(end,:);

function  Z=D2adj(Y)
Z=[-Y(:,1)  -diff(Y,1,2)];   Z(:,end)=Z(:,end)+ Y(:,end);