function [Z,tv,res]=tv_PETquant( beta,TOL,MAXIT,search,alpha, actM, fitdata)
% tv_tik.m
% TV minimization with quadratic constraint for PET parametric imaging
% min |x|_{TV,beta}+alpha ||Ax-b||^2, where linear system Ax=b are equations from loagn or Patlak method for all voxels.
%|x|_{TV,beta}= \sum sqrt(gi'*gi+beta), gi is the gradient of voxel i.
%
% beta      % for smoothing TV calculation
% TOL       %Tolerance
% MAXIT     %maximum number of iteration
% search   %1-line search after Newton method. 0-no search.
% alpha   % the regularization parameter(constant). =0.05 for our FDG data.
% actM,     %the Mask, a binary image, zeros for background and ones for active pixels.
% fitdata   % the coefficient matrix of Patlak or Logan method. Let C_P(t)
% and C_T(t) be the input function and time activity curve respectively. For Patlak method, the left columns are integral(C_P)/C_P(t) and the right
% columns are C_T/C_P(t) for active pixels. Each column is a vector over
% time points in the assumed "equilibrium time window".
%
% Z=[estK, intercept, u]
%tv,  this vector hold toltal variations for each iteration.
%
%
% Copyright 2010, Hongbin Guo, Arizona State University
% Email: hb_guo@asu.edu
% Created: September 2007
%revised on Feb 2010
%
% The users should acknowledge this software and cite the following paper in all manuscripts:
% H. Guo,   R. Renaut, K. Chen, E. Reiman. FDG-PET Parametric Imaging by Total Variation Minimization,
% Computerized Medical Imaging and Graphics,  volume 33, issue 4, page 295-303, 2009.
%
% Note: test data and main code are also available upon request by email.

[N,pt]=size(fitdata); pt=pt/2; %N=the # of active pixels, pt=# of time points in the "equilibrium time window".
sz=size(actM);% the image size.
actind=find(actM>0);
inact=setdiff(1:sz(1)*sz(2),actind);
[bdind, neighb]=bound(actM,actind);
[locr,locl,subr,subl,idr,idl]=locrl(N,sz,actind);
xdata=fitdata(:,1:pt);
ydata=fitdata(:,pt+1:end);
%prepare data for normal equation of each pixel: ATA=[xTx, sumx; sumx, pt], ATb=[xTy;sumy];
eqdata=[sum(xdata.^2,2), sum(xdata,2), sum(xdata.*ydata,2), sum(ydata,2)]; %[xTx,sumx,xTy,sumy]
[x0,y0]=linfit(eqdata,pt);
M=bdupdate(sz,x0,neighb,actind); %boundary update
sigmasq=zeros(N,1);
u=zeros(N,1);
k=0;
Z=[x0;y0;u];
while 1
    conf=con_fun(Z,M,sigmasq,actind,fitdata);
    phi=sqrt(conf(1:N)+beta);
    tv(k+1)=sum(phi);
    res=sum(conf(N+1:2*N));
    dphi=1./(2*phi);
    ddphi=-0.25*phi.^(-3);
    lambda=[dphi;alpha*ones(N,1);zeros(N,1)];
    obj(k+1)=sum(phi)+alpha*sum(conf(N+1:2*N));
   
    J=diff_mat(Z,M,actind,eqdata,pt,locr,locl,subr,subl,idr,idl);
    grd=J(1:2*N,1:2*N)*lambda(1:2*N);
    H=hess_tik(Z,M,N,lambda,actind,eqdata,pt,locr,locl,subr,subl,J,ddphi);
    X = qmr(H,-grd,TOL,MAXIT);
    Z = Z + [X;u]; xt=Z(1:N);
    M=bdupdate(sz,xt,neighb,actind);


    %one step line search
    if search==1
        conf=con_fun(Z,M,sigmasq,actind,fitdata);
        phi=sqrt(conf(1:N)+beta);
        ddphi=-0.25*phi.^(-3);
        lambda=[1./(2*phi);alpha*ones(N,1);zeros(N,1)];
        J=diff_mat(Z,M,actind,eqdata,pt,locr,locl,subr,subl,idr,idl);
        grd=J(1:2*N,1:2*N)*lambda(1:2*N);
        H=hess_tik(Z,M,N,lambda,actind,eqdata,pt,locr,locl,subr,subl,J,ddphi);
        s=-(grd'*X )/(X'*H*X);
        sl(k+1)=s;
        Z = Z + [s*X;u]; xt=Z(1:N);
        M=bdupdate(sz,xt,neighb,actind);
    end
    %end of line search
    k=k+1;
    if k>10
        break
    end
end

    function M=bdupdate(sz,x,neighb,actind)
        M=zeros(sz);
        M(actind)=x;
        [I,J]=ind2sub(sz,neighb);
        idr_b=sub2ind(sz,I,J+1);
        idd_b=sub2ind(sz,I+1,J);
        idl_b=sub2ind(sz,I,J-1);
        idu_b=sub2ind(sz,I-1,J);
        count= ismember(idr_b,actind)+ismember(idl_b,actind)+ismember(idu_b,actind)+ismember(idd_b,actind);
        M(neighb)=(M(idr_b)+M(idl_b)+M(idu_b)+M(idd_b))./count;
    end


    function [bdind, neighb]=bound(M,actind);
        % generate boundary indices.
        sz=size(M);
        M(actind)=1;
        [I,J]=ind2sub(sz,actind);
        idr=sub2ind(sz,I,J+1);
        idl=sub2ind(sz,I+1,J);
        M0=M(idr);M1=M(idl);
        M2=M0+M1;
        ind1=find(M2<1.5);
        bdind=actind(ind1);
        id2=find(M0(ind1)<0.5);
        id3=find(M1(ind1)<0.5);
        neighb=[idr(ind1(id2));idl(ind1(id3))];
        neighb=unique(neighb);
    end

    function [v,r]=con_fun(Z,M,epsilonsq,actind,fitdata)
        %x--the solpe
        %y-interept
        [N,pt]=size(fitdata); pt=pt/2;
        N=length(Z)/3;
        x=Z(1:N); y=Z(N+1:2*N); u=Z(2*N+1:3*N);
        szrow=size(M,1);
        T1=diff(M,1,2); T1(szrow,szrow)=0;
        T2=diff(M,1,1); T2(szrow,szrow)=0;
        T=T1.^2+T2.^2;
        v=T(actind)-u.^2;
        r=(x*ones(1,pt)).*fitdata(:,1:pt)+y*ones(1,pt)-fitdata(:,pt+1:2*pt);
        v(N+1:2*N)=sum(r.^2,2)-epsilonsq;
        v(2*N+1:3*N)=-u;
    end


    function J=diff_mat(Z,M,actind,eqdata,pt,locr,locl,subr,subl,idr,idl)
        %creat Jacobian matrix
        N=length(Z)/3;
        x=Z(1:N); y=Z(N+1:2*N); u=Z(2*N+1:3*N);
        % build A11
        v0=4*M(actind)-2*M(idr)-2*M(idl);
        v1=2*M(idr)-2*M(actind); v1=v1(subr);
        v2=2*M(idl)-2*M(actind); v2=v2(subl);
        % derivative of subr-th function with respect to loc1-th variables = v1;
        % derivative of subl-th function with respect to loc2-th variables = v2;
        A11=sparse(locr,subr,v1,N,N)+sparse(locl,subl,v2,N,N)+sparse(1:N,1:N,v0);

        %build A12, A22
        %4 collumns matrix eqdata=[xTx,sumx,xTy,sumy], may from Patlak or Logan formulation.
        %solve normaleqs: ATA=[xTx, sumx; sumx, pt], ATb=[xTy;sumy]
        Pm1=2*(eqdata(:,1).*x+eqdata(:,2).*y-eqdata(:,3));
        Pm2=2*(eqdata(:,2).*x+pt*y-eqdata(:,4));
        A12=sparse(1:N,1:N,Pm1);
        A22=sparse(1:N,1:N,Pm2);
        A31=sparse(1:N,1:N,-2*u);
        A33=sparse(1:N,1:N,-ones(N,1));
        J=[A11,A12,sparse(N,N); sparse(N,N),A22, sparse(N,N); A31, sparse(N,N),A33];
    end

    function [x,y]=linfit(d,pt);
        % 4 collumns matrix d=[xTx,sumx,xTy,sumy], may from Patlak or Logan formulation.
        %solve normaleqs: ATA=[xTx, sumx; sumx, pt], ATb=[xTy;sumy] corresponding
        %to each row of d by Gaussian elemination.

        right=(d(:,4)-d(:,3).*d(:,2)./d(:,1));
        left=(pt-d(:,2).^2./d(:,1));
        y=right./left;
        x=(d(:,3)-d(:,2).*y)./d(:,1);
    end


    function [locr,locl,subr,subl,idr,idl]=locrl(N,sz,actind)
        %locations of the right and lower neighbours.
        %N=length(x);
        %sz=size(M);
        [I,J]=ind2sub(sz,actind);
        idr=sub2ind(sz,I,J+1);
        idl=sub2ind(sz,I+1,J);
        %inact=setdiff(1:sz(1)*sz(2),actind);
        TF1=ismember(idr,actind); subr=find(TF1>0);
        TF2=ismember(idl,actind); subl=find(TF2>0);
        [TF1,locr]=ismember(idr(subr),actind);
        [TF2,locl]=ismember(idl(subl),actind);
    end

    function H=hess_tik(Z,M,N,lambda,actind,eqdata,pt,locr,locl,subr,subl,J,ddphi)
        %creat sum Hessian matrices with weight lambda
        H=hess_pd(Z,M,N,lambda,actind,eqdata,pt,locr,locl,subr,subl);
        H=H(1:2*N,1:2*N);
        [i1,j1,s1] = find(J(1:N,1:N)*sparse(1:N,1:N,ddphi)*J(1:N,1:N)');
        [i2,j2,s2] = find(H);
        H=sparse([i1;i2],[j1;j2],[s1;s2]);
        H=H+sparse(1:2*N,1:2*N,1e-6*norm(H,'fro')*ones(2*N,1));

        function H=hess_pd(Z,M,N,lambda,actind,eqdata,pt,locr,locl,subr,subl)
            %creat sum Hessian matrices with weight lambda

            N=length(Z)/3;
            x=Z(1:N); y=Z(N+1:2*N); u=Z(2*N+1:3*N);
            % constraint functions 1:N
            A11=sparse(locr,locr,lambda(subr),N,N)-sparse(locr,subr,lambda(subr),N,N)-sparse(subr,locr,lambda(subr),N,N);
            A11=A11+sparse(locl,locl,lambda(subl),N,N)-sparse(locl,subl,lambda(subl),N,N)-sparse(subl,locl,lambda(subl),N,N);
            A11=2*A11+4*sparse(1:N,1:N,lambda(1:N));
            A33=-2*sparse(1:N,1:N,lambda(1:N));
            H1=[A11,sparse(N,N),sparse(N,N); sparse(N,N),sparse(N,N), sparse(N,N); sparse(N,N), sparse(N,N),A33];
            % constraint functions N+1:2N
            %4 collumns matrix eqdata=[xTx,sumx,xTy,sumy], may from Patlak or Logan formulation.
            %solve normal eqs: ATA=[xTx, sumx; sumx, pt], ATb=[xTy;sumy]
            coe=lambda(N+1:2*N);
            B11=2*sparse(1:N,1:N,coe.*eqdata(:,1));
            B12=2*sparse(1:N,1:N,coe.*eqdata(:,2));
            B22=2*pt*sparse(1:N,1:N,coe);
            H2=[B11,B12,sparse(N,N); B12,B22, sparse(N,N); sparse(N,N), sparse(N,N),sparse(N,N)];
            H=H1+H2;
        end
    end

end