dynamicSPECT/matlab/fitCentersPixelIVF.m

function [globalPar, cPars] = ...
    fitCentersPixelIVF(path, patientID, cax, cm, centers,sigma2)
    %path is the core directory
    %cax are the mid intervals in images
    %cm are the interval durations (1,n)
    %data are the images
    %nclass is number of typical responses sought
    npar=1;
    nglob=2;
    nclass=size(centers,1);
    
    cm2=ones(nclass,1)*cm';
    
    %w=sqrt(cm2.*centers)./sum(cm2.*centers,2);
    w=ones(nclass,1)*sqrt(cm)';
    %make sure that the function that is most similar to input function has
    %a high weight on maximum value, which comes as the second reading
    %abandonded, because decay time gets shortened and fitting of the slow
    %rise is given up in favor of meeting the high point through k1 value
    %[~,im]=max(centers(:,2));
    %w(im,2)=w(im,2)+0.1*(sqrt(max(cm))-w(im,2));
    
    
    %for individual fits
    valueGeneratorSettings1.type='loggaus';
    valueGeneratorSettings1.parameters=[-3,1];
    valueGeneratorSettings1.u=[0,Inf];
    ivar=2;
        
    valueGeneratorSettings1(ivar).type='flat';
    valueGeneratorSettings1(ivar).parameters=[0,1];
    valueGeneratorSettings1(ivar).u=[0,1];
    ivar=ivar+1;
    
    valueGeneratorSettings1(ivar).type='loggaus';
    valueGeneratorSettings1(ivar).parameters=[-7,2];
    valueGeneratorSettings1(ivar).u=[0,Inf];
    ivar=ivar+1;
    
    valueGeneratorSettings1(ivar).type='gaus';
    valueGeneratorSettings1(ivar).parameters=[10,5];
    valueGeneratorSettings1(ivar).u=[0,Inf];
        
        
    options = optimoptions(@lsqnonlin);
    options.Display='off'; 
    %[options.OptimalityTolerance,options.FunctionTolerance,options.StepTolerance];
    %holds best estimate
    n1=40;
    nFitIndividual=4;
    qfit=zeros(1,nFitIndividual*nclass);
    
    %use input functions estimate from center fits
    p=[readFitParameters(path,patientID,10); ...
        readFitParameters(path,patientID,20); ...
        readFitParameters(path,patientID,30)];
    fp=median(p);%4 values, const, fall & rise time, delay
    
    fprintf('Global fit: (%f,%f,%f,%f)\n',fp(1),fp(2),fp(3),fp(4));

    chiSum=0;
    %now fit every curve individually with the known input function
    for j=1:nclass

        cay=centers(j,:)';
        w=sqrt(cm);

        fun1=@(par)evaluateSingleDiff(cax,cay,w,fp,par);

        qpar=zeros(1,4);
        qmin=1e30;
        for k=1:n1
            [par0,ul1,ub1]=generateInitialValues(valueGeneratorSettings1);
            [par1,resnorm,residual,exitflag,output,lambda,jacobian]  = ...
                lsqnonlin(fun1,par0,ul1,ub1,options);
            fchi2=sqrt(residual'*residual);
            if fchi2<qmin
                qmin=fchi2;
                qpar=par1;
            end
        end
        qfit((j-1)*nFitIndividual+1:j*nFitIndividual)=qpar;
        fprintf('Center: %d (%f %f %f %f)/(%f %f %f %f)\n',j,...
            par0(1),par0(2),par0(3),par0(4),...
            par1(1),par1(2),par1(3),par1(4));
        %residual
        chiSum=chiSum+qmin;
    end
    %writematrix(ci,fullfile(path,sprintf('CI%d.txt',i)),'Delimiter','tab');
    %output
    chi2=chiSum;

    fprintf('Done %f \n',chi2); 
        
    %this is the result that will be used by k1 calculator
    sigmaCode=sprintf('%.2f',sigma2);
    sigmaCode=strrep(sigmaCode,'.','p');
    code=sprintf('%s_%d_%s_PixelIVF',patientID,nclass,sigmaCode);
    ofile=fullfile(path,patientID,sprintf('%s_fitParFinal.txt',code));
    try
      writematrix([chi2 fp qfit],ofile,'Delimiter','tab');
    catch ME
      dlmwrite(ofile,[chi2 fp qfit],'\t');
    end
    %ofile=fullfile(path,sprintf('CIFinal_%d_%d.txt',nclass,realizationId));
    %writematrix(fci,ofile,'Delimiter','tab');
    
    globalPar=[fp(1:3) 0 0 0];
    
    cPars=reshape(qfit,[],nclass)';
    %reset k1 for those where BVF is nearly one
    %cPars(cPars(:,2)>0.4,1)=0;
    %reset k1 for those where k2 is abnormally large
    %cPars(cPars(:,3)>1,1)=0;
    opath=fullfile(path,patientID);
    plotFunctions(opath,code,globalPar,cPars,cax,centers);
    
    
end