Up to generate centers, including matlab files

matlab/HOWTO

@@ -0,0 +1,5 @@
+
+loadData
+#sets cax, data, caspline, cm
+xa=dynamicAnalysis(cax,data,caspline,cm)
+writeData

matlab/README.md

@@ -0,0 +1,8 @@
+## Usage ##
+
+Master routines are written in `analyze.m`. Full usage:
+
+```shell
+patientId='XXX'
+analyze
+```

matlab/analyze.m

@@ -0,0 +1,20 @@
+%reset random generator
+rng shuffle
+
+%add nrrd processing capability
+addpath('nrrd_read_write_rensonnet')
+
+path='/home/studen/temp/dynamicSPECT';
+[cax,cm]=loadTime(path,patientID);
+
+nclass=str2num(nclass);
+realizationId=str2num(realizationId);
+centers=loadCenters(path,patientID,nclass,realizationId);
+%this will write out fitPar
+[globalPar, cPars] = fitCenters(path, patientID, realizationId, cax, cm, centers);
+
+%[globalPar, cPars, xa, U]=fitFromClusters(path,cax, cm, data, nclass,realizationId);
+
+disp('DynamicAnalysis done');
+
+%writeData

matlab/analyzeIVF.m

@@ -0,0 +1,21 @@
+%reset random generator
+rng shuffle
+
+%add nrrd processing capability
+addpath('nrrd_read_write_rensonnet')
+
+path='/home/studen/temp/dynamicSPECT';
+[cax,cm]=loadTime(path,patientID);
+
+nclass=str2num(nclass);
+realizationId=str2num(realizationId);
+centers=loadCenters(path,patientID,nclass,realizationId);
+
+%this will write out fitPar
+[globalPar, cPars] = fitCentersIVF(path, patientID, realizationId, cax, cm, centers);
+
+%[globalPar, cPars, xa, U]=fitFromClusters(path,cax, cm, data, nclass,realizationId);
+
+disp('DynamicAnalysisIVF done');
+
+%writeData

matlab/analyzePixel.m

@@ -0,0 +1,26 @@
+%reset random generator
+rng shuffle
+
+%add nrrd processing capability
+addpath('nrrd_read_write_rensonnet')
+
+path='/home/studen/temp/dynamicSPECT';
+[cax,cm]=loadTime(path,patientID);
+
+segm=tdfread(fullfile(path,'Segmentation.tsv'));
+v0=loadPixels(segm,patientID);
+
+data=loadSPECTdata(path,patientID,cm);
+
+sigma2=str2num(sigma2);
+na=7;
+fcenters=extractCenters(data,v0,sigma2,na);
+
+%this will write out fitPar
+[globalPar, cPars] = fitCentersPixel(path, patientID, cax, cm, fcenters,sigma2);
+
+%[globalPar, cPars, xa, U]=fitFromClusters(path,cax, cm, data, nclass,realizationId);
+
+disp('DynamicPixelAnalysis done');
+
+%writeData

matlab/analyzePixelIVF.m

@@ -0,0 +1,26 @@
+%reset random generator
+rng shuffle
+
+%add nrrd processing capability
+addpath('nrrd_read_write_rensonnet')
+
+path='/home/studen/temp/dynamicSPECT';
+[cax,cm]=loadTime(path,patientID);
+
+segm=tdfread(fullfile(path,'Segmentation.tsv'));
+v0=loadPixels(segm,patientID);
+
+data=loadSPECTdata(path,patientID,cm);
+
+sigma2=str2num(sigma2);
+na=7;
+fcenters=extractCenters(data,v0,sigma2,na);
+
+%this will write out fitPar
+[globalPar, cPars] = fitCentersPixelIVF(path, patientID, cax, cm, fcenters,sigma2);
+
+%[globalPar, cPars, xa, U]=fitFromClusters(path,cax, cm, data, nclass,realizationId);
+
+disp('DynamicPixelAnalysis done');
+
+%writeData

matlab/evaluateCenterDiff.m

@@ -0,0 +1,5 @@
+function fg=evaluateCenterDiff(cax,centers,w,gParConst,x,pars)
+    globalPar=[gParConst(3) x(1:2) 0 gParConst(1:2)];
+    fv=evaluateCenters(globalPar,pars,cax);
+    fg=reshape((fv-centers).*w,[],1);
+end

matlab/evaluateCenters.m

@@ -0,0 +1,58 @@
+function fv=evaluateCenters(globalPar, par, cax)
+    %par is a nfunc x npar matrix of coefficients
+    %globalPar is a row vector of global coefficients
+    %cax is a column vector of time instances
+    
+    oneVector=ones(1,size(cax,1));
+    cx2=ones(1,size(par,1))'*cax';
+    oneMatrix=ones(size(cx2));
+    
+    A=globalPar(1);
+    tau=globalPar(2);
+    alpha=globalPar(3);
+    
+    t0=globalPar(5);
+    t1=globalPar(6);
+    
+    k1=par(:,1)*oneVector;
+    B=par(:,2)*oneVector;
+    k2=par(:,3)*oneVector;
+    dt=par(:,4)*oneVector;
+    
+    kt=k1./(k2-oneMatrix/tau);
+    ka=k1./(k2-oneMatrix*alpha);
+    cx2p=cx2-dt;
+    
+    e0=exp(-cx2p/tau);
+    e1=exp(-cx2p*alpha);
+    ek=exp(-cx2p.*k2);
+    
+    r0=zeros(size(cx2));
+    if alpha==1/tau
+        fv=e0.*cx2p;
+        rsel=oneMatrix/tau==k2;
+        r0(rsel)=k1(rsel).*cx2p(rsel).*cx2p(rsel)*ek(rsel);
+        r0(~rsel)=kt.*(fv-(e0-ek)./(k2-oneMatrix/tau));
+        
+    else
+        q=tau/(alpha*tau-1);
+        fv=q*(e0-e1);
+        rsel=oneMatrix/tau==k2;
+        r0(rsel)=r0(rsel)+q*k1(rsel).*cx2p(rsel);
+        r0(~rsel)=r0(~rsel)+q*kt(~rsel).*(e0(~rsel)-ek(~rsel));
+        rsel=alpha==k2;
+        r0(rsel)=r0(rsel)-q*k1(rsel).*cx2p(rsel);
+        r0(~rsel)=r0(~rsel)-q*ka(~rsel).*(e1(~rsel)-ek(~rsel));
+        
+    B1=oneMatrix-B;
+    
+    fv=B.*fv+B1.*r0;
+   
+    fv=fv*A;
+    sel=cx2p<0;
+    fv(sel)=0;
+    
+
+end
+
+%par=ones(2,

matlab/evaluateSingleDiff.m

@@ -0,0 +1,5 @@
+function fg=evaluateSingleDiff(cax,cay,w,x,par)
+    globalPar=[x(1:3) 0 0 0];
+    fv=fitSingle(globalPar,par,cax);
+    fg=(fv-cay).*w;
+end

matlab/evaluateSingleDiffZero.m

@@ -0,0 +1,5 @@
+function fg=evaluateSingleDiffZero(cax,cay,w,x)
+    globalPar=[x(1:4) 0 0];
+    fv=fitSingleZero(globalPar,cax);
+    fg=(fv-cay).*w;
+end

matlab/extractCenters.m

@@ -0,0 +1,32 @@
+function fcenters=extractCenters(data,v0,sigma2,na)
+%extract response by segment centered around point v0, taking weights of surrounding pixels as gaussian with sigma2
+%nd is number of time samples
+%v0 is nr by 3 matrix of centers
+%na is the neighborhood size (in pixels, should work great for odd na)
+    nr=size(v0,1);
+    [~,~,~,nd]=size(data);
+    fcenters=zeros(nr,nd);
+    na2=int8(0.5*(na+1));
+    for fi=1:nr
+        fs=0;
+        fv=v0(fi,:);%zero beased coordinates
+        idx0=int8(fv);
+        idx=idx0;
+        for i=1:na
+            idx(1)=idx0(1)+i-na2;
+            for j=1:na
+                idx(2)=idx0(2)+j-na2;
+                for k=1:na
+                    idx(3)=idx0(3)+k-na2;
+                    df=double(idx)-fv;
+                    fd=df*df';
+                    fw=exp(-0.5*fd/sigma2);
+                    fcenters(fi,:)=fcenters(fi,:)+fw*reshape(data(idx(1)+1,idx(2)+1,idx(3)+1,:),1,[]);
+                    fs=fs+fw;
+                end
+            end
+        end
+        fcenters(fi,:)=fcenters(fi,:)/fs;
+
+    end
+end

matlab/fitCenters.m

@@ -0,0 +1,173 @@
+function [globalPar, cPars] = ...
+    fitCenters(path, patientID, realizationId, cax, cm, centers)
+    %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));
+    
+    
+    %cax are the delays
+    %c
+    valueGeneratorSettings.type='poisson';
+    valueGeneratorSettings.parameters=[25];
+    valueGeneratorSettings.u=[0,Inf];
+        
+    ivar=2;
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[1.5,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+    
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[-1,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+
+    valueGeneratorSettings(ivar).type='gaus';
+    valueGeneratorSettings(ivar).parameters=[10,2];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[-Inf,Inf];
+    ivar=ivar+1;
+
+    %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
+    fmin=1e30;
+    n=20;
+    n1=20;
+    imin=-1;
+    nFitIndividual=4;
+    qfit=zeros(1,nFitIndividual*nclass);
+    qfitmin=qfit;
+    [x,ul,ub]=generateInitialValues(valueGeneratorSettings);
+    for i=1:n
+    
+        %first estimate input functions
+        %take function with maximum response as IVF surrogate
+        [m mi]=max(centers);
+        [fm mj]=max(m);
+        cay=centers(mi(mj),:)';
+        w=ones(size(cay));
+        
+        fun=@(x)evaluateSingleDiffZero(cax,cay,w,x);
+        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+        
+        [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
+        fprintf('A %f\n',x1(1));
+        
+        while x1(1)<1
+            [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+            [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
+            fprintf('A %f\n',x1(1));
+        end
+        %this gets A, tau and alpha
+        %A shouldn't be trivial
+        
+        fprintf('[%d] Global fit: (%f,%f,%f %f) (%f,%f,%f %f)\n',i,x0(1),x0(2),x0(3),x0(4),...
+            x1(1),x1(2),x1(3),x1(4));
+        
+        chiSum=0;
+        %now fit every curve individually with the known input function
+        for j=1:nclass
+            cay=centers(j,:)';
+            w=ones(size(cay));
+            w=sqrt(cm);
+            
+            fun1=@(par)evaluateSingleDiff(cax,cay,w,x1,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;
+        
+        if chi2<fmin
+            fmin=chi2;
+            x=x1;
+            imin=i;
+            qfitmin=qfit;
+            %fci=ci;
+        end
+        fprintf('chi2: %d,%d\n',size(chi2,1),size(chi2,2));    
+        %writematrix([chi2 x1],fullfile(path,sprintf('fitPar%d.txt',i)),'Delimiter','tab');
+        fprintf('Done %d %f (%.2f,%.2f,%.2f)/(%.2f,%.2f,%.2f)\n',...
+            i,chi2,x0(1),x0(2),x0(3),x1(1),x1(2),x1(3)); 
+        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+        
+        
+    end
+    
+    %this is the result that will be used by k1 calculator
+    code=sprintf('%s_%d_%d',patientID,nclass,realizationId);
+    ofile=fullfile(path,patientID,sprintf('%s_fitParFinal.txt',code));
+    writematrix([fmin x qfitmin],ofile,'Delimiter','tab');
+    %ofile=fullfile(path,sprintf('CIFinal_%d_%d.txt',nclass,realizationId));
+    %writematrix(fci,ofile,'Delimiter','tab');
+    disp(sprintf('Found min at %d',imin)); 
+    
+    globalPar=[x(1:3) 0 0 0];
+    
+    cPars=reshape(qfitmin,[],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
+

matlab/fitCentersIVF.m

@@ -0,0 +1,115 @@
+function [globalPar, cPars] = ...
+    fitCentersIVF(path, patientID, realizationId, cax, cm, centers)
+    %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=20;
+    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
+    code=sprintf('%s_%d_%d_IVF',patientID,nclass,realizationId);
+    ofile=fullfile(path,patientID,sprintf('%s_fitParFinal.txt',code));
+    writematrix([chi2 fp qfit],ofile,'Delimiter','tab');
+    %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
+

matlab/fitCentersPixel.asv

@@ -0,0 +1,178 @@
+function [globalPar, cPars] = ...
+    fitCentersPixel(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));
+    
+    
+    %cax are the delays
+    %c
+    valueGeneratorSettings.type='poisson';
+    valueGeneratorSettings.parameters=[25];
+    valueGeneratorSettings.u=[0,Inf];
+        
+    ivar=2;
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[1.5,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+    
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[-1,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+
+    valueGeneratorSettings(ivar).type='gaus';
+    valueGeneratorSettings(ivar).parameters=[10,2];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[-Inf,Inf];
+    ivar=ivar+1;
+
+    %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
+    fmin=1e30;
+    n=20;
+    n1=20;
+    imin=-1;
+    nFitIndividual=4;
+    qfit=zeros(1,nFitIndividual*nclass);
+    qfitmin=qfit;
+    [x,ul,ub]=generateInitialValues(valueGeneratorSettings);
+    for i=1:n
+    
+        %first estimate input functions
+        %take function with maximum response as IVF surrogate
+        %assume first row (ventricle) is the input function
+        cay=centers(1,:)';
+        w=ones(size(cay));
+        
+        fun=@(x)evaluateSingleDiffZero(cax,cay,w,x);
+        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+        
+        [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
+        fprintf('A %f\n',x1(1));
+        
+        while x1(1)<1
+            [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+            [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
+            fprintf('A %f\n',x1(1));
+        end
+        %this gets A, tau and alpha
+        %A shouldn't be trivial
+        
+        fprintf('[%d] Global fit: (%f,%f,%f %f) (%f,%f,%f %f)\n',i,x0(1),x0(2),x0(3),x0(4),...
+            x1(1),x1(2),x1(3),x1(4));
+        
+        chiSum=0;
+        %now fit every curve individually with the known input function
+        for j=1:nclass
+             [m mi]=max(centers);
+        [fm mj]=max(m);
+        cay=centers(mi(mj),:)';
+        w=ones(size(cay));
+            cay=centers(j,:)';
+            w=ones(size(cay));
+            w=sqrt(cm);
+            
+            fun1=@(par)evaluateSingleDiff(cax,cay,w,x1,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;
+        
+        if chi2<fmin
+            fmin=chi2;
+            x=x1;
+            imin=i;
+            qfitmin=qfit;
+            %fci=ci;
+        end
+        fprintf('chi2: %d,%d\n',size(chi2,1),size(chi2,2));    
+        %writematrix([chi2 x1],fullfile(path,sprintf('fitPar%d.txt',i)),'Delimiter','tab');
+        fprintf('Done %d %f (%.2f,%.2f,%.2f)/(%.2f,%.2f,%.2f)\n',...
+            i,chi2,x0(1),x0(2),x0(3),x1(1),x1(2),x1(3)); 
+        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+        
+        
+    end
+    
+    %this is the result that will be used by k1 calculator
+    sigmaCode=sprintf('%.2f',sigma2);
+    sigmaCode=strrep(sigmaCode,'.','p');
+    code=sprintf('%s_%d_%s_Pixel',patientID,nclass,sigmaCode);
+    ofile=fullfile(path,patientID,sprintf('%s_fitParFinal.txt',code));
+    writematrix([fmin x qfitmin],ofile,'Delimiter','tab');
+    %ofile=fullfile(path,sprintf('CIFinal_%d_%d.txt',nclass,realizationId));
+    %writematrix(fci,ofile,'Delimiter','tab');
+    disp(sprintf('Found min at %d',imin)); 
+    
+    globalPar=[x(1:3) 0 0 0];
+    
+    cPars=reshape(qfitmin,[],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
+

matlab/fitCentersPixel.m

@@ -0,0 +1,181 @@
+function [globalPar, cPars] = ...
+    fitCentersPixel(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));
+    
+    
+    %cax are the delays
+    %c
+    valueGeneratorSettings.type='poisson';
+    valueGeneratorSettings.parameters=[25];
+    valueGeneratorSettings.u=[0,Inf];
+        
+    ivar=2;
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[1.5,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+    
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[-1,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+
+    valueGeneratorSettings(ivar).type='gaus';
+    valueGeneratorSettings(ivar).parameters=[10,2];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[-Inf,Inf];
+    ivar=ivar+1;
+
+    %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
+    fmin=1e30;
+    n=20;
+    n1=20;
+    imin=-1;
+    nFitIndividual=4;
+    qfit=zeros(1,nFitIndividual*nclass);
+    qfitmin=qfit;
+    [x,ul,ub]=generateInitialValues(valueGeneratorSettings);
+    for i=1:n
+    
+        %first estimate input functions
+        %take function with maximum response as IVF surrogate
+        %assume first row (ventricle) is the input function
+        
+        %sometimes maximum crops up in unexpected places. Assume maximum
+        %response is the best IVF estimate
+        [m mi]=max(centers);
+        [fm mj]=max(m);
+        cay=centers(mi(mj),:)';
+            
+        %cay=centers(1,:)';
+        w=ones(size(cay));
+        
+        fun=@(x)evaluateSingleDiffZero(cax,cay,w,x);
+        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+        
+        [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
+        fprintf('A %f\n',x1(1));
+        
+        while x1(1)<1
+            [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+            [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
+            fprintf('A %f\n',x1(1));
+        end
+        %this gets A, tau and alpha
+        %A shouldn't be trivial
+        
+        fprintf('[%d] Global fit: (%f,%f,%f %f) (%f,%f,%f %f)\n',i,x0(1),x0(2),x0(3),x0(4),...
+            x1(1),x1(2),x1(3),x1(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,x1,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;
+        
+        if chi2<fmin
+            fmin=chi2;
+            x=x1;
+            imin=i;
+            qfitmin=qfit;
+            %fci=ci;
+        end
+        fprintf('chi2: %d,%d\n',size(chi2,1),size(chi2,2));    
+        %writematrix([chi2 x1],fullfile(path,sprintf('fitPar%d.txt',i)),'Delimiter','tab');
+        fprintf('Done %d %f (%.2f,%.2f,%.2f)/(%.2f,%.2f,%.2f)\n',...
+            i,chi2,x0(1),x0(2),x0(3),x1(1),x1(2),x1(3)); 
+        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+        
+        
+    end
+    
+    %this is the result that will be used by k1 calculator
+    sigmaCode=sprintf('%.2f',sigma2);
+    sigmaCode=strrep(sigmaCode,'.','p');
+    code=sprintf('%s_%d_%s_Pixel',patientID,nclass,sigmaCode);
+    ofile=fullfile(path,patientID,sprintf('%s_fitParFinal.txt',code));
+    writematrix([fmin x qfitmin],ofile,'Delimiter','tab');
+    %ofile=fullfile(path,sprintf('CIFinal_%d_%d.txt',nclass,realizationId));
+    %writematrix(fci,ofile,'Delimiter','tab');
+    disp(sprintf('Found min at %d',imin)); 
+    
+    globalPar=[x(1:3) 0 0 0];
+    
+    cPars=reshape(qfitmin,[],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
+

matlab/fitCentersPixelIVF.m

@@ -0,0 +1,116 @@
+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=20;
+    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));
+    writematrix([chi2 fp qfit],ofile,'Delimiter','tab');
+    %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
+

matlab/fitFromClusters.asv

@@ -0,0 +1,174 @@
+function [globalPar, cPars, xa, U] = ...
+    fitFromClusters(path, cax, cm, data, nclass, realizationId)
+    %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=3;
+    nglob=3;
+    
+    %correct plotFunctions and evaluateCenterDiff
+    
+    c1=size(cax,1);
+    A=reshape(data,[],c1);
+
+    %find most common responses
+    %nclass=10;
+    [centers,U]=fcm(A,nclass);
+
+    %data(ix,iy,iz,:)=A((iz-1)*i1*i2+(iy-1)*i1+ix,:)
+
+    cm2=ones(size(centers,1),1)*cm';
+    
+    %w=sqrt(cm2.*centers)./sum(cm2.*centers,2);
+    w=ones(size(centers,1),1)*sqrt(cm)';
+    
+    pars=cax(1:2)';
+    
+    fun=@(x)evaluateCenterDiff(cax,centers,w,pars,x(1:nglob),...
+        reshape(x(nglob+1:nglob+npar*nclass),[],nclass)');
+    %cax are the delays
+    %c
+    
+    valueGeneratorSettings.type='poisson';
+    valueGeneratorSettings.parameters=[25];%mu
+    valueGeneratorSettings.u=[0,Inf];
+    
+    ivar=2;
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[2,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+    
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[0,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+    
+    
+    for j=1:nclass
+        
+        valueGeneratorSettings(ivar).type='loggaus';
+        valueGeneratorSettings(ivar).parameters=[-3,1];
+        valueGeneratorSettings(ivar).u=[0,Inf];
+        ivar=ivar+1;
+        
+        valueGeneratorSettings(ivar).type='flat';
+        valueGeneratorSettings(ivar).parameters=[0,1];
+        valueGeneratorSettings(ivar).u=[0,1];
+        ivar=ivar+1;
+    
+        valueGeneratorSettings(ivar).type='loggaus';
+        valueGeneratorSettings(ivar).parameters=[-7,2];
+        valueGeneratorSettings(ivar).u=[0,Inf];
+        ivar=ivar+1;
+    end
+    
+    %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];
+        
+        
+    [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+    options = optimoptions(@lsqnonlin);
+    options.Display='off'; 
+    %[options.OptimalityTolerance,options.FunctionTolerance,options.StepTolerance];
+    %holds best estimate
+    x=x0;
+    fmin=1e30;
+    n=20;
+    imin=-1;
+    for i=1:n
+    
+        %find best fit for nclass responses. 
+        %total number of points that are fitted is 20*nclass, total number of
+        %parameters estimated is 4+3*nclass
+        [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
+        
+        %this gets tau and A
+        
+        %now fit every curve individually with the known input function
+        for j=1:nclass
+            cay=centers(j,:)';
+            w=ones(size(cay));
+            
+            fun1=@(par)evaluateSingleDiff(cax,cay,w,pars,x1,par);
+            
+            qpar=zeros(1,3);
+            qmin=1e30;
+            for k=1:n
+                [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
+            x1(nglob+(j-1)*3+1:nglob+j*3)=qpar;
+            disp(sprintf('Center: %d (%f %f %f)/(%f %f %f)',j,...
+                par0(1),par0(2),par0(3),...
+                par1(1),par1(2),par1(3)));
+            %residual
+        
+        end
+        %writematrix(ci,fullfile(path,sprintf('CI%d.txt',i)),'Delimiter','tab');
+        %output
+        chi2=fun(x1);
+        chi2=chi2'*chi2;
+        
+        if chi2<fmin
+            fmin=chi2;
+            x=x1;
+            imin=i;
+            %fci=ci;
+        end
+            
+        %writematrix([chi2 x1],fullfile(path,sprintf('fitPar%d.txt',i)),'Delimiter','tab');
+        disp(sprintf('Done %d %f (%.2f %.2f %.2f)/(%.2f,%.2f,%.2f)',...
+            i,chi2,x0(1),x0(2),x0(3),); 
+        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+        
+        
+    end
+    ofile=fullfile(path,sprintf('fitParFinal_%d_%d.txt',nclass,realizationId));
+    writematrix([fmin x],ofile,'Delimiter','tab');
+    %ofile=fullfile(path,sprintf('CIFinal_%d_%d.txt',nclass,realizationId));
+    %writematrix(fci,ofile,'Delimiter','tab');
+    disp(sprintf('Found min at %d',imin)); 
+    
+    globalPar=[x(1:3) 0 cax(1:2)'];
+    
+    cPars=reshape(x(nglob+1:nglob+npar*nclass),[],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;
+    
+    plotFunctions(path,globalPar,cPars,cax,centers);
+    
+    
+    dt=size(data);
+    
+    k1v=U'*cPars(:,1);
+    xa(:,:,:,1)=reshape(k1v,dt(1:3));
+
+    vbfv=U'*cPars(:,2);
+    xa(:,:,:,2)=reshape(vbfv,dt(1:3));
+    
+    k2v=U'*cPars(:,3);
+    xa(:,:,:,3)=reshape(k2v,dt(1:3));
+end
+

matlab/fitFromClusters.m

@@ -0,0 +1,175 @@
+function [globalPar, cPars, xa, U] = ...
+    fitFromClusters(path, cax, cm, data, nclass, realizationId)
+    %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=3;
+    nglob=3;
+    
+    %correct plotFunctions and evaluateCenterDiff
+    
+    c1=size(cax,1);
+    A=reshape(data,[],c1);
+
+    %find most common responses
+    %nclass=10;
+    [centers,U]=fcm(A,nclass);
+
+    %data(ix,iy,iz,:)=A((iz-1)*i1*i2+(iy-1)*i1+ix,:)
+
+    cm2=ones(size(centers,1),1)*cm';
+    
+    %w=sqrt(cm2.*centers)./sum(cm2.*centers,2);
+    w=ones(size(centers,1),1)*sqrt(cm)';
+    
+    pars=cax(1:2)';
+    
+    fun=@(x)evaluateCenterDiff(cax,centers,w,pars,x(1:nglob),...
+        reshape(x(nglob+1:nglob+npar*nclass),[],nclass)');
+    %cax are the delays
+    %c
+    
+    valueGeneratorSettings.type='poisson';
+    valueGeneratorSettings.parameters=[25];%mu
+    valueGeneratorSettings.u=[0,Inf];
+    
+    ivar=2;
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[1.5,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+    
+    valueGeneratorSettings(ivar).type='loggaus';
+    valueGeneratorSettings(ivar).parameters=[-1,1];%log10(mean),log10(std)
+    valueGeneratorSettings(ivar).u=[0,Inf];
+    ivar=ivar+1;
+    
+    
+    for j=1:nclass
+        
+        valueGeneratorSettings(ivar).type='loggaus';
+        valueGeneratorSettings(ivar).parameters=[-3,1];
+        valueGeneratorSettings(ivar).u=[0,Inf];
+        ivar=ivar+1;
+        
+        valueGeneratorSettings(ivar).type='flat';
+        valueGeneratorSettings(ivar).parameters=[0,1];
+        valueGeneratorSettings(ivar).u=[0,1];
+        ivar=ivar+1;
+    
+        valueGeneratorSettings(ivar).type='loggaus';
+        valueGeneratorSettings(ivar).parameters=[-7,2];
+        valueGeneratorSettings(ivar).u=[0,Inf];
+        ivar=ivar+1;
+    end
+    
+    %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];
+        
+        
+    [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+    options = optimoptions(@lsqnonlin);
+    options.Display='off'; 
+    %[options.OptimalityTolerance,options.FunctionTolerance,options.StepTolerance];
+    %holds best estimate
+    x=x0;
+    fmin=1e30;
+    n=20;
+    n1=20;
+    imin=-1;
+    for i=1:n
+    
+        %find best fit for nclass responses. 
+        %total number of points that are fitted is 20*nclass, total number of
+        %parameters estimated is 4+3*nclass
+        [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
+        
+        %this gets tau and A
+        
+        %now fit every curve individually with the known input function
+        for j=1:nclass
+            cay=centers(j,:)';
+            w=ones(size(cay));
+            
+            fun1=@(par)evaluateSingleDiff(cax,cay,w,pars,x1,par);
+            
+            qpar=zeros(1,3);
+            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
+            x1(nglob+(j-1)*3+1:nglob+j*3)=qpar;
+            fprintf('Center: %d (%f %f %f)/(%f %f %f)\n',j,...
+                par0(1),par0(2),par0(3),...
+                par1(1),par1(2),par1(3));
+            %residual
+        
+        end
+        %writematrix(ci,fullfile(path,sprintf('CI%d.txt',i)),'Delimiter','tab');
+        %output
+        chi2=fun(x1);
+        chi2=chi2'*chi2;
+        
+        if chi2<fmin
+            fmin=chi2;
+            x=x1;
+            imin=i;
+            %fci=ci;
+        end
+            
+        %writematrix([chi2 x1],fullfile(path,sprintf('fitPar%d.txt',i)),'Delimiter','tab');
+        fprintf('Done %d %f (%.2f,%.2f,%.2f)/(%.2f,%.2f,%.2f)\n',...
+            i,chi2,x0(1),x0(2),x0(3),x1(1),x1(2),x1(3)); 
+        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
+        
+        
+    end
+    ofile=fullfile(path,sprintf('fitParFinal_%d_%d.txt',nclass,realizationId));
+    writematrix([fmin x],ofile,'Delimiter','tab');
+    %ofile=fullfile(path,sprintf('CIFinal_%d_%d.txt',nclass,realizationId));
+    %writematrix(fci,ofile,'Delimiter','tab');
+    disp(sprintf('Found min at %d',imin)); 
+    
+    globalPar=[x(1:3) 0 cax(1:2)'];
+    
+    cPars=reshape(x(nglob+1:nglob+npar*nclass),[],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;
+    
+    plotFunctions(path,globalPar,cPars,cax,centers);
+    
+    
+    dt=size(data);
+    
+    k1v=U'*cPars(:,1);
+    xa(:,:,:,1)=reshape(k1v,dt(1:3));
+
+    vbfv=U'*cPars(:,2);
+    xa(:,:,:,2)=reshape(vbfv,dt(1:3));
+    
+    k2v=U'*cPars(:,3);
+    xa(:,:,:,3)=reshape(k2v,dt(1:3));
+end
+

matlab/fitSingle.m

@@ -0,0 +1,57 @@
+function fv=fitSingle(globalPar, par, cax)
+    %par is a nfunc x npar matrix of coefficients
+    %globalPar is a row vector of global coefficients
+    %cax is a column vector of time instances
+    %cay are the function values
+    
+    %oneVector=ones(1,size(cax,1));%column vector
+    
+    A=globalPar(1);
+    tau=globalPar(2);
+    alpha=globalPar(3);
+    
+    k1=par(1);
+    B=par(2);
+    k2=par(3);
+    dt=par(4);
+    
+    cax1=cax-dt;
+    kt=k1./(k2-1/tau);
+    ka=k1./(k2-alpha);
+    e0=exp(-cax1/tau);
+    e1=exp(-cax1*alpha);
+    ek=exp(-cax1*k2);
+    
+    
+    if alpha==1/tau
+        fv=cax1.*e0;
+        if k2==1/tau
+            r0=k1.*cax1.*cax.*ek;
+        else
+            r0=kt*(cax1.*e0-(e0-ek)/(k2-1/tau));
+        end
+    else
+        q=tau/(alpha*tau-1);
+        fv=q*(e0-e1);
+        r0=zeros(size(cax1));
+        if k2==1/tau
+            r0=r0+q*k1*cax1;
+        else
+            r0=r0+q*kt*(e0-ek);
+        end
+        if k2==alpha
+            r0=r0-q*k1*cax1;
+        else
+            r0=r0-q*ka*(e1-ek);
+        end
+    end
+    
+    fv=B*fv+(1-B)*r0;
+    
+    fv=fv*A;
+    sel=cax1<0;
+    fv(sel)=0;
+
+end
+
+%par=ones(2,

matlab/fitSingleRamp.m

@@ -0,0 +1,50 @@
+function fv=fitSingle(globalPar, par, cax)
+    %par is a nfunc x npar matrix of coefficients
+    %globalPar is a row vector of global coefficients
+    %cax is a column vector of time instances
+    %cay are the function values
+    
+    %oneVector=ones(1,size(cax,1));%column vector
+    
+    A=globalPar(1);
+    tau=globalPar(2);
+    t0=globalPar(5);
+    t1=globalPar(6);
+    
+    k1=par(1);
+    B=par(2);
+    k2=par(3);
+    
+    
+    kt=k1./(1/tau-k2);
+    e0=exp(-(cax-t1)/tau);
+    ek=exp(-(cax-t1).*k2);
+    
+    ek1=(1-exp(-k2.*(cax-t0)))./(k2*(t1-t0));
+    fk1=(1-exp(-k2*(t1-t0)))./(k2*(t1-t0));
+    fr1=k1./k2.*(1-fk1);
+    B1=1-B;
+    
+    dt=(cax-t0)/(t1-t0);
+    
+    if k2==1/tau
+        r0=k1*(cax-t1).*e0;
+    else
+        r0=kt*(ek-e0);
+    
+    fv=zeros(size(cax));
+    sel=cax>t0 & cax<=t1;
+    fv(sel)=B*dt(sel)+B1*k1/k2*(dt(sel)-ek1(sel));
+    
+   
+    sel=cax>t1;
+    
+    fv(sel)=B*e0(sel)+B1*(fr1+r0(sel));
+    
+    
+    fv=fv*A;
+    
+
+end
+
+%par=ones(2,

matlab/fitSingleZero.m

@@ -0,0 +1,28 @@
+function fv=fitSingleZero(globalPar, cax)
+    %par is a nfunc x npar matrix of coefficients
+    %globalPar is a row vector of global coefficients
+    %cax is a column vector of time instances
+    %cay are the function values
+    
+    %oneVector=ones(1,size(cax,1));%column vector
+    
+    A=globalPar(1);
+    tau=globalPar(2);
+    alpha=globalPar(3);
+    dt=globalPar(4);
+    
+    cax1=cax-dt;
+    if tau==1/alpha
+        fv=cax1.*exp(-cax/tau);
+    else
+        fv=tau*(exp(-alpha*cax1)-exp(-cax1/tau))/(1-alpha*tau);
+    end
+    fv=fv*A;
+    sel=cax1<0;
+    fv(sel)=0;
+    %fv=A*power(cax,alpha).*exp(-cax/tau)/power(tau,alpha);
+   
+
+end
+
+%par=ones(2,

matlab/generateCenters.m

@@ -0,0 +1,10 @@
+rng shuffle
+loadData
+nclass=str2num(nclass);
+nRealizations=str2num(nRealizations);
+%sets cax, data, caspline, cm, path
+%[xa,data1]=dynamicAnalysis(cax,data,caspline,cm,1);
+saveCenters(path,patientID,cax, cm, data, nclass,nRealizations);
+%with biexp fit, chat is equivalent to caspline
+%[xa,data1]=dynamicAnalysis(cax,data,chat,cm);
+disp(sprintf('[%s %d %d] generateCenters done',patientID,nclass,nRealizations));

matlab/generateInitialValues.m

@@ -0,0 +1,34 @@
+function [x0,ul,ub]=generateInitialValues(valueGeneratorSettings)
+%initial settings
+        nvar=size(valueGeneratorSettings,2);
+        
+        x0=zeros(1,nvar);
+        ul=zeros(1,nvar);
+        ub=zeros(1,nvar);
+        for ivar=1:nvar
+            t=valueGeneratorSettings(ivar).type;
+            p=valueGeneratorSettings(ivar).parameters;
+            u=valueGeneratorSettings(ivar).u;
+            if strcmp(t,'gaus')
+                x0(ivar)=normrnd(p(1),p(2));
+            end
+            if strcmp(t,'loggaus')
+                x0(ivar)=power(10,normrnd(p(1),p(2)));
+            end
+            if strcmp(t,'poisson')
+                x0(ivar)=poissrnd(p(1));
+            end
+            if strcmp(t,'flat')
+                x0(ivar)=p(1)+(p(2)-p(1))*rand;
+            end
+            if strcmp(t,'const')
+                x0(ivar)=p(1);
+            end
+            
+            
+            ul(ivar)=u(1);
+            ub(ivar)=u(2);
+        
+        end
+        
+end

matlab/getData.m

@@ -0,0 +1,26 @@
+function ty=getData(data,i,j,k) 
+%getData at i,j,k
+    xsize=size(data);
+    ty=squeeze(data(i,j,k,:));
+    %add neighbours for smoothing
+    tw=1;
+    fc=[i j k];
+    for na=1:6
+         ao=zeros(size(fc));
+         i2=idivide(int16(na-1),2);
+         ao(i2+1)=2*rem(na-1,2)-1;
+         %disp(ao)
+         fc1=fc+single(ao);
+         if any(fc1<1)
+            continue
+         end
+         if any(fc1>xsize(1:3))
+            continue
+         end
+         w=0.5/sum(ao.*ao);
+         %disp(fc1)
+         ty=ty+w*squeeze(data(fc1(1),fc1(2),fc1(3),:));
+         tw=tw+w;
+    end
+    ty=ty./tw;
+end

matlab/loadCenters.m

@@ -0,0 +1,16 @@
+function centers=loadCenters(path, patientID, nclass,realizationId)
+    %path is the core directory
+    %nclass is number of typical responses sought
+    %realizationId is the realization number
+   
+    centers=zeros(1);
+    for i=1:nclass
+        code=sprintf('%s_%d_%d_center%d',patientID,nclass,realizationId,i);
+        ofile=fullfile(path,patientID,sprintf('%s_center.txt',code));
+        fy=readmatrix(ofile,'Delimiter',',');
+        if size(centers,1)==1
+            centers=zeros(nclass,size(fy,2));
+        end  
+        centers(i,:)=fy;
+    end
+end

matlab/loadData.m

@@ -0,0 +1,50 @@
+addpath('nrrd_read_write_rensonnet')
+
+%patientID='2SMobr';
+
+path='/home/studen/temp/dynamicSPECT';
+
+%load ventricle data
+%caFile=[ patientID '_Ventricle.mcsv'];
+%caFile=fullfile(path,caFile);
+%opts = detectImportOptions(caFile,'FileType','text');
+%ca=readmatrix(caFile,opts);
+
+[cax,cm]=loadTime(path,patientID);
+%w=cay/sum(cay);
+%w(1)=1;
+%caspline=csaps(cax,cay,p,[],w);
+
+%x=lsqnonlin(biexp,[20 3 10 5.5 0.1],[0,0,0,0,0],[Inf,Inf,Inf, Inf, Inf]);
+%chat0=[20, 3, 10, 5.5, 0.1];
+%chat0=[10 0.8 10 100];
+%lb=zeros(size(chat0));
+%ub=Inf*ones(size(chat0));
+%chat0=[chat0 cax(1:2)'];
+%lb=[lb cax(1:2)'];
+%ub=[ub cax(1:2)'];
+%chat=lsqcurvefit(@biexp,chat0 ,cax,cay,lb,ub);
+
+[c1 c2]=size(cax);
+
+volumeNamePattern=[ patientID '_Volume'];
+s0=[volumeNamePattern '0.nrrd'];
+ofile=fullfile(path,patientID,sprintf('%s_Volume0.nrrd',patientID));
+headerInfo = nhdr_nrrd_read(ofile,1);
+a0=headerInfo.data;
+[i1,i2,i3]=size(a0);
+
+data=zeros(i1,i2,i3,c1);
+
+for i=1:c1
+    s=sprintf('%s_Volume%d.nrrd',patientID,i-1);
+    ofile=fullfile(path,patientID,s);
+    hInfo = nhdr_nrrd_read(ofile,1);
+    hInfo.data(hInfo.data<0)=0;
+    data(:,:,:,i)=hInfo.data/cm(i);
+end
+
+%A=reshape(data,[],c1);
+%find most common responses
+%[centers,U]=fcm(A,5);
+%data(ix,iy,iz,:)=A((iz-1)*i1*i2+(iy-1)*i1+ix,:)

matlab/loadPixels.m

@@ -0,0 +1,19 @@
+function v0=loadPixels(segm, patientID)
+    fs=string(segm.Alias_ID);
+    for i=1:size(fs,1)
+        fs(i)=erase(fs(i)," ");
+    end
+
+    regionIds=segm.Region_Id(fs==patientID);
+    nr=size(regionIds,1);
+    vtemp=zeros(nr,3);
+    %rotate to match convention below
+    vtemp(:,3)=segm.X(fs==patientID);
+    vtemp(:,2)=segm.Y(fs==patientID);
+    vtemp(:,1)=segm.Z(fs==patientID);
+
+    v0=vtemp;
+    for j=1:nr
+        v0(regionIds(j)+1,:)=vtemp(j,:);
+    end
+end

matlab/loadResults.m

@@ -0,0 +1,11 @@
+pars=zeros(i1,i2,i3,2);
+%files={'k1.nrrd','k2.nrrd','vbf.nrrd'};
+files={'k1.nrrd','vbf.nrrd'};
+
+for i=1:2
+    h1=nhdr_nrrd_read(fullfile(path,[patientID '_' files{i}]),1);
+    pars(:,:,:,i)=h1.data;
+    
+end
+
+

matlab/loadSPECTdata.m

@@ -0,0 +1,24 @@
+function data = loadSPECTdata(path,patientID,cm)
+%loadData Summary of this function goes here
+%   Detailed explanation goes here
+    
+    c1=size(cm,1);
+    
+    ofile=fullfile(path,patientID,sprintf('%sVolume0.nrrd',patientID));
+    headerInfo = nhdr_nrrd_read(ofile,1);
+    a0=headerInfo.data;
+    [i1,i2,i3]=size(a0);
+
+    data=zeros(i1,i2,i3,c1);
+
+    for i=1:c1
+        s=sprintf('%sVolume%d.nrrd',patientID,i-1);
+        ofile=fullfile(path,patientID,s);
+        hInfo = nhdr_nrrd_read(ofile,1);
+        hInfo.data(hInfo.data<0)=0;
+        data(:,:,:,i)=hInfo.data/cm(i);
+    end
+
+
+end
+

matlab/loadSegmentation.m

@@ -0,0 +1,22 @@
+
+
+fs=string(segm.Alias_ID);
+for i=1:size(fs,1);
+    fs(i)=erase(fs(i)," ");
+end
+
+regionIds=segm.Region_Id(fs==patientID);
+nr=size(regionIds,1);
+vtemp=zeros(nr,3);
+%rotate to match convention below
+vtemp(:,3)=segm.X(fs==patientID);
+vtemp(:,2)=segm.Y(fs==patientID);
+vtemp(:,1)=segm.Z(fs==patientID);
+
+v0=vtemp;
+for j=1:nr
+    v0(regionIds(j)+1,:)=vtemp(j,:);
+end
+
+
+

matlab/loadTime.asv

@@ -0,0 +1,44 @@
+[cax, cm]=function loadTime(path,patientID)
+    
+
+end
+
+%patientID='2SMobr';
+
+path=[ path '/' patientID];
+
+%load ventricle data
+%caFile=[ patientID '_Ventricle.mcsv'];
+%caFile=fullfile(path,caFile);
+%opts = detectImportOptions(caFile,'FileType','text');
+%ca=readmatrix(caFile,opts);
+
+caDTFile=[ patientID '_Dummy.mcsv'];
+caDTFile=fullfile(path,caDTFile);
+opts = detectImportOptions(caDTFile,'FileType','text');
+cma=readmatrix(caDTFile,opts);
+
+
+%convert to seconds->multiply time by 1e-3 and value (integral/t) with 1e3
+cax=cma(:,1)*1e-3;
+%cay=ca(:,2)*1e3;
+cm=cma(:,2)*1e-3;
+%this is highly tuned. Assume nframe=20
+%p=[1e-3 1e-2*ones(1,1) 1e-1*ones(1,3) ones(1,2) 10*ones(1,8) 100*ones(1,5)];
+
+%w=cay/sum(cay);
+%w(1)=1;
+%caspline=csaps(cax,cay,p,[],w);
+
+%x=lsqnonlin(biexp,[20 3 10 5.5 0.1],[0,0,0,0,0],[Inf,Inf,Inf, Inf, Inf]);
+%chat0=[20, 3, 10, 5.5, 0.1];
+%chat0=[10 0.8 10 100];
+%lb=zeros(size(chat0));
+%ub=Inf*ones(size(chat0));
+%chat0=[chat0 cax(1:2)'];
+%lb=[lb cax(1:2)'];
+%ub=[ub cax(1:2)'];
+%chat=lsqcurvefit(@biexp,chat0 ,cax,cay,lb,ub);
+
+centers=loadCenters(path,patientID,nclass,realizationId);
+

matlab/loadTime.m

@@ -0,0 +1,13 @@
+function [cax, cm]=loadTime(path,patientID)
+    ofile=fullfile(path,patientID,sprintf('%s_Dummy.mcsv',patientID));
+    opts = detectImportOptions(ofile,'FileType','text');
+    cma=readmatrix(ofile,opts);
+    %convert to seconds->multiply time by 1e-3 and value (integral/t) with 1e3
+    cax=cma(:,1)*1e-3;
+    %cay=ca(:,2)*1e3;
+    cm=cma(:,2)*1e-3;
+    
+
+
+end
+

matlab/nrrdWriter.m

@@ -0,0 +1,185 @@
+% ========================================================================
+% 
+% nrrdwriter_dan
+% 
+% filename  - 'myimage.ext' - 'veins.nrrd'
+% matrix    - data - Matlab matrix
+% pixelspacing - boxel size
+% origin    - point from the image is generated
+% encoding  - raw, ascii, gzip
+% 
+% ========================================================================
+
+function ok = nrrdWriter(filename, matrix, pixelspacing, origin, encoding)
+
+% This line gets the path, name and extension of our file:
+% pathf = /home/mario/.../myfile.myext
+% fname = myfile
+% ext = .myext
+[pathf, fname, ext] = fileparts(filename);
+
+format=ext(2:end); % We remove the . from .ext
+% so we extract the output format from the argument filename, instead of
+% put two different arguments
+
+matrix = permute(matrix, [2 1 3]); % so we undo permute of index in nrrdreader
+
+dims=(size(matrix));    % matrix dimensions (size NxMxP)
+ndims=length(dims);     % number of dimensions (dim n)
+
+
+
+% =====================================================================
+% Conditions to make sure our file is goint to be created succesfully.
+% 
+% First the code puts the argument 'encoding' in lowercase
+encoding = lower(encoding);
+
+encodingCond = isequal(encoding, 'ascii') || isequal(encoding, 'raw') || isequal(encoding, 'gzip');
+assert(encodingCond, 'Unsupported encoding')
+
+% The same with output format
+format = lower(format);
+formatCond = isequal(format,'nhdr') || isequal(format,'nrrd');
+assert(formatCond, 'Unexpected format');
+
+% ======================================================================
+
+% Now, if our conditions are satisfied:
+if (encodingCond && formatCond)
+    
+    % Header
+    
+    % Open, filename (which specifies output format) and write binary
+    fid = fopen(filename, 'wb');
+    fprintf(fid,'NRRD0004\n');      % NRRD type 4
+    
+    % Type of variable we're storing in our file
+    mtype=class(matrix);
+    outtype=setDatatype(mtype);
+    fprintf(fid,['type: ', outtype, '\n']);
+    
+    % 
+    fprintf(fid,['dimension: ', num2str(ndims), '\n']);
+    
+    if isequal(ndims, 2)
+        fprintf(fid,'space: left-posterior\n');
+    elseif isequal (ndims, 3)
+        fprintf(fid,'space: left-posterior-superior\n');
+    end
+
+    fprintf(fid,['sizes: ', num2str(dims), '\n']);
+    
+    if isequal(ndims, 2)
+        fprintf(fid,['space directions: (', num2str(pixelspacing(1)), ...
+            ',0) (0,', num2str(pixelspacing(2)), ')\n']);
+        fprintf(fid,'kinds: domain domain\n');
+    elseif isequal (ndims, 3)
+        fprintf(fid,['space directions: (', num2str(pixelspacing(1)), ...
+            ',0,0) (0,', num2str(pixelspacing(2)), ',0) (0,0,', ...
+            num2str(pixelspacing(3)), ')\n']);
+        fprintf(fid,'kinds: domain domain domain\n');
+    end
+    
+    fprintf(fid,['encoding: ', encoding, '\n']);
+    
+    [~,~,endian] = computer();
+    
+    if (isequal(endian, 'B'))
+        fprintf(fid,'endian: big\n');
+    else
+        fprintf(fid,'endian: little\n');
+    end
+    
+    if isequal(ndims, 2)
+        fprintf(fid,['space origin: (', num2str(origin(1)),',', num2str(origin(2)),')\n']);
+    elseif isequal (ndims, 3)
+        fprintf(fid,['space origin: (', num2str(origin(1)), ...
+            ',',num2str(origin(2)),',', num2str(origin(3)),')\n']);
+    end    
+    
+    if (isequal(format, 'nhdr')) % Si hay que separar
+        % Escribir el nombre del fichero con los datos
+        fprintf(fid, ['data file: ', [fname, '.', encoding], '\n']);
+        
+        fclose(fid);
+        if isequal(length(pathf),0)
+            fid = fopen([fname, '.', encoding], 'wb');
+        else
+            fid = fopen([pathf, filesep, fname, '.', encoding], 'wb');
+        end
+    else
+        fprintf(fid,'\n');
+    end
+    
+    ok = writeData(fid, matrix, outtype, encoding);
+    fclose(fid);
+end
+
+
+% ========================================================================
+% Determine the datatype --> From mtype (matlab) to outtype (NRRD) -->    
+% ========================================================================
+function datatype = setDatatype(metaType)
+
+% Determine the datatype
+switch (metaType)
+ case {'int8', 'uint8', 'int16', 'uint16', 'int32', 'uint32', 'int64',...
+       'uint64', 'double'}
+   datatype = metaType;
+  
+ case {'single'}
+  datatype = 'float';
+  
+ otherwise
+  assert(false, 'Unknown datatype')
+end
+   
+% HACER!!!!!!!!!!!!!!!!!!!!!!!!!
+% ========================================================================
+% writeData -->
+% fidIn is the open file we're overwriting
+% matrix - data that have to be written
+% datatype - type of data: int8, string, double...
+% encoding - raw, gzip, ascii
+% ========================================================================
+function ok = writeData(fidIn, matrix, datatype, encoding)
+
+switch (encoding)
+ case {'raw'}
+  
+  ok = fwrite(fidIn, matrix(:), datatype);
+  
+ case {'gzip'}
+     
+     % Store in a raw file before compressing
+     tmpBase = tempname(pwd);
+     tmpFile = [tmpBase '.gz'];
+     fidTmpRaw = fopen(tmpBase, 'wb');
+     assert(fidTmpRaw > 3, 'Could not open temporary file for GZIP compression');
+     
+     fwrite(fidTmpRaw, matrix(:), datatype);
+     fclose(fidTmpRaw);
+     
+     % Now we gzip our raw file
+     gzip(tmpBase);
+     
+     % Finally, we put this info into our nrrd file (fidIn)
+     fidTmpRaw = fopen(tmpFile, 'rb');
+     tmp = fread(fidTmpRaw, inf, [datatype '=>' datatype]);
+     cleaner = onCleanup(@() fclose(fidTmpRaw));
+     ok = fwrite (fidIn, tmp, datatype);
+     
+     delete (tmpBase);
+     delete (tmpFile);
+
+
+ case {'ascii'}
+  
+  ok = fprintf(fidIn,'%u ',matrix(:));
+  %ok = fprintf(fidIn,matrix(:), class(matrix));
+  
+ otherwise
+  assert(false, 'Unsupported encoding')
+end
+

matlab/nrrd_read_write_rensonnet/nhdr_nrrd_read.m

@@ -0,0 +1,968 @@
+% Read NHDR/NRRD files 
+% 
+% headerInfo = nhdr_nrrd_read(nrrdFileName, bReadData)
+%
+% Inputs:
+% * nrrdFileName (char): path to NHDR/NRRD file, either a detached header,
+% a detached header pointing to detached data files or a NRRD standalone
+% file with header and data included
+% * bReadData (bool): set to true to read the data and store it in
+% headerInfo.data, set to false to just import the header without the data
+%
+% Outputs:
+% * headerInfo (struct): contains all the fields specified in the file
+% header. If data was read it is contained in the 'data' field. That
+% structure can be fed to the nhdr_nrrd_write module as is and produce
+% valid NHDR/NRRD files.
+% 
+% Format definition:
+% http://teem.sourceforge.net/nrrd/format.html
+%
+% A few supported NRRD features:
+% - detached headers with all variants of 'data file:'
+% - raw, txt/text/ascii, gz/gzip encodings
+% - definition of space and orientation
+% - handling of diffusion-weighted MRI data with '<key>:=<value>' lines 
+%   'modality:=DWMRI', 'DWMRI_b-value:=' and 'DWMRI_gradient_0000:=', 
+%   'DWMRI_gradient_0001:=', 'DWMRI_gradient_0002:=', etc.
+%   (see https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:Nrrd_format)
+% 
+%
+% Other features:
+% - exits cleanly upon error (no file should be left open)
+% - informative error messages whenever possible
+% - warnings printed to console 
+% - NHDR/NRRD writer fills out missing fields by inspecting the data
+%   whenever possible
+% 
+%
+% Unsupported features:
+% 
+% In general, any header field that we were unable to parse is reported
+% in a message printed to the console. Specific examples of
+% unsupported features include:
+% 
+% - reading data along more than 1 dimension or along a dimension other
+%   than the slowest (last) axis specified by the optional <subdim>, as in
+%   'data file: <format> <min> <max> <step> [<subdim>]' or 'data file: LIST
+%   [<subdim>]'
+% - storing all the comments found in headers
+% - hex encoding
+% - bz2/bzip2 encoding
+% - byte skip; can only accept -1 with raw encoding, 0 for all other
+%   encodings
+% - line skip: can only accept 0
+% - checking that field specifications of the form '<field>: <desc>' 
+%   appear no more than once in the NRRD header (unlike '<key>:=<value>'
+%   lines)
+%
+% Date: July 2018
+% Author: Gaetan Rensonnet
+%
+%
+%
+% Notes :
+%
+% This is intended as a more comprehensive and accurate implementation of
+% the NRRD format specification than most of the Matlab scripts that have
+% been proposed so far. We try to fail graciously when an unsupported
+% feature of the NRRD format is encountered. One of our main contributions
+% is to propose a writer module which is compatible with the read module,
+% in that the output of one can be given as an argument to the other to
+% read or produce equivalent NHDR/NRRD files. This is still a version with
+% much room for improvement.
+%
+% The following contributions inspired parts of our Matlab read/write
+% modules:
+% 
+% 1. The body of the writer module was pretty much written from scratch but
+% the general structure of the reader's main body is based on the Matlab
+% functions maReadNrrdHeader.m and maReadNrrdData.m by marc@bwh.harvard.edu
+% and kquintus@bwh.harvard.edu (unpublished). Many additions were made and
+% a few bugs fixed.
+%
+% 2. Jeff Mather's NRRD reader
+% (http://nl.mathworks.com/matlabcentral/fileexchange/34653-nrrd-format-file-reader)
+% and
+% http://jeffmatherphotography.com/dispatches/2012/02/writing-a-file-reader-in-matlab/)
+% provided the auxiliary functions:
+% - adjustEndian
+% - getDatatype, which we renamed nrrd_getMatlabDataType and now throws
+% a gracious error if it encounters 'block'-type data,
+% - readData was adapted to include a cross-platform fix to delete
+% temporary files when using gzip encoding. David Feng's fix used a
+% Windows-specific command to delete temporary files
+% (https://nl.mathworks.com/matlabcentral/fileexchange/50830-nrrd-format-file-reader).
+%
+% 3. mdcacio's nrrdWriter
+% (https://nl.mathworks.com/matlabcentral/fileexchange/48621-nrrdwriter-filename--matrix--pixelspacing--origin--encoding-)
+% provided the auxiliary functions:
+% - writeData(): we got rid of the 'unexpected end of input stream when
+% attempting to gunzip the file' error when using gzip encoding, which we
+% later found had been fixed by Quan Chen independenlty and in a similar
+% manner.
+% - setDatatype(), which we renamed get_nrrd_datatype() and is the
+% reciprocal of getDatatype 
+%
+
+function headerInfo = nhdr_nrrd_read(nrrdFileName, bReadData)
+
+[mainFpath,mainFileName,mainFext] = fileparts(nrrdFileName);
+
+headerInfo = struct();
+
+% default header:
+headerInfo.content = mainFileName;        % default value, overwritten if content field is set
+headerInfo.data = [];
+
+
+fidr = fopen(nrrdFileName, 'r');
+
+if (fidr == -1)
+    error('ABORT: %s does not exist.\n', nrrdFileName);
+end
+
+try 
+    if ~(strcmpi(mainFext, '.nhdr') || strcmpi(mainFext, '.nrrd'))
+        warning('%s looks like a %s file, not a nhdr or nrrd file.\n', nrrdFileName, mainFext );
+    end
+    
+    % Magic line
+    cs  = fgetl(fidr);
+    assert(numel(cs) >= 8 && isequal(cs(1:4), 'NRRD'),...
+        'Bad signature. First line should be magic line of type "NRRD000X" with X an integer between 1 and 5.'); % FIXME should throw an error if a bad integer is provided
+    nrrd_version = sscanf(cs(5:end), '%d');
+    if nrrd_version > 5
+        error('This reader only supports versions of the NRRD file format up to 5. Detected %d.', nrrd_version)
+    end
+
+    % Always optional: defining orientation
+    define_orientation = 0;         % internal-use variable
+    
+    % Parse header
+    while ~feof(fidr)
+        
+        cs = fgetl(fidr);         % content string
+        
+        if isempty(cs)
+            % End of header
+            break;
+        end
+                
+        if foundKeyword( 'CONTENT:', cs )
+            
+            headerInfo.content = strtrim( cs( length('CONTENT:')+1:end ) );
+            
+        elseif foundKeyword('TYPE:', cs )
+            
+            headerInfo.type = strtrim( cs( length('TYPE:')+1:end ) );
+            
+        elseif foundKeyword('ENDIAN:', cs )
+            
+            headerInfo.endian = strtrim( cs( length('ENDIAN:')+1:end ) );
+            
+        elseif foundKeyword('ENCODING:', cs )
+            
+            headerInfo.encoding = strtrim( cs( length('ENCODING:')+1:end ) );
+            
+        elseif foundKeyword('DIMENSION:', cs )
+            
+            headerInfo.dimension = sscanf( cs( length('DIMENSION:')+1:end ), '%i' );
+            
+        elseif foundKeyword('SIZES:', cs )
+            
+            iSizes = sscanf( cs(length('SIZES:')+1:end), '%i' );
+            headerInfo.sizes = iSizes(:)';      % store as row vector
+            
+        elseif foundKeyword('KINDS:', cs )
+            
+            headerInfo.kinds = extractStringList( cs(length('KINDS:')+1:end) ); % bug fixed with extractStringList where 2 entries are present
+            % FIXME: check that axis sizes match each kind according to nrrd standard
+            
+        elseif foundKeyword('SPACE:', cs )
+            % Starts defining orientation (either space or space dimension,
+            % not both)
+            define_orientation = 1;
+            
+            if isfield(headerInfo, 'spacedimension')
+                fprintf(['WARNING nhdr_nrrd_read %s:\n ''space'' field specifier will ' ...
+                    'be checked for consistency but will be ignored afterwards ' ...
+                    'because ''space dimension'' was specified before.\n'], fopen(fidr)); 
+            end
+            
+            tmp_space = strtrim( cs( length('SPACE:')+1:end ) );
+            tmp_spacedimension = nrrd_getSpaceDimensions(tmp_space);
+            if tmp_spacedimension <= 0
+                error('%s: unrecognized ''space'' descriptor ''%s''.', fopen(fidr), tmp_space)
+            end
+            
+            if isfield(headerInfo, 'spacedimension')
+                % internal_spacedimension already set
+                if internal_spacedimension ~= tmp_spacedimension
+                    error(['%s: ''space'' field specifier implies a spatial ' ...
+                            '(world) dimension equal to %d, which differs from ' ...
+                            'the ''space dimension'' field specifier set to %d.'],...
+                            fopen(fidr), tmp_spacedimension, internal_spacedimension)
+                end
+                % if no inconsistencies found, just ignore space field
+            else
+                % Set space info for the first time:
+                headerInfo.space = tmp_space;
+                internal_spacedimension = tmp_spacedimension;  % for internal use only
+            end
+            
+        elseif foundKeyword('SPACE DIMENSION:', cs)
+            % Starts defining orientation (either space or space dimension,
+            % not both)
+            define_orientation = 1;
+            if isfield(headerInfo, 'space')
+                fprintf(['WARNING nhdr_nrrd_read %s:\n ''space dimension'' field specifier ' ...
+                    ' will be checked for consistency but will be ignored afterwards ' ...
+                    'because ''space'' was specified before.\n'],...
+                    fopen(fidr));
+            end
+            
+            tmp_spacedimension = sscanf( cs( length('SPACE DIMENSION:')+1:end), '%i' );
+            if numel(tmp_spacedimension) ~= 1
+                error(['%s: ''space dimension'' should be specified as one' ...
+                    ' integer number. Found %d element(s) instead.'],...
+                      fopen(fidr), numel(tmp_spacedimension))
+            end
+            if tmp_spacedimension <= 0
+                error('%s: ''space dimension'' should be specified as a strictly positive integer (found %d).',...
+                                fopen(fidr), tmp_spacedimension)
+            end
+            
+            if isfield(headerInfo, 'space')
+                if tmp_spacedimension ~= internal_spacedimension
+                    error(['%s: ''space dimension'' field specifier set to %d, ' ...
+                        'which differs from the space (world) dimension implied by the ' ...
+                        '''space'' field specifier which is equal to %d.'],...
+                            fopen(fidr), tmp_spacedimension, internal_spacedimension)
+                end
+                % if no inconsistencies found, just ignore space dimension
+                % field
+            else
+                % Set space info for the first time:
+                headerInfo.spacedimension = tmp_spacedimension;
+                internal_spacedimension = tmp_spacedimension;  % for internal use only
+            end
+            
+        elseif foundKeyword('SPACE DIRECTIONS:', cs )
+            % Required if orientation defined but must come after space or
+            % space dimension
+            % space directions: <vector[0]> <vector[1]> ... <vector[dim-1]>
+            % The format of the <vector> is as follows. The vector is
+            % delimited by "(" and ")", and the individual components are
+            % comma-separated. 
+            if ~define_orientation
+               error('%s: field specifier ''space directions'' cannot be set before ''space'' or ''space dimension''.',fopen(fidr))
+            end
+
+            space_dir_tmp = strtrim(cs(length('SPACE DIRECTIONS:')+1:end));     % remove leading and trailing spaces
+
+            spacedir_vecs = strsplit(space_dir_tmp); %  cell array of strings after split at {' ','\f','\n','\r','\t','\v'}
+            SD_data = zeros(internal_spacedimension, internal_spacedimension);
+            
+            % Check each vector: either none or (f1,...,f_spaceDim) with fi
+            % a floating-point number
+            cnt_space_vectors = 0;
+            for i = 1:length(spacedir_vecs)
+                none_start_index = strfind(lower(spacedir_vecs{i}), 'none');
+                if ~isempty(none_start_index)
+                    % Axis-specific entry contains substring none
+                    if ~strcmpi(spacedir_vecs{i}, 'none')
+                        fprintf(['WARNING nhdr_nrrd_read: detected %s instead of ' ...
+                            'expected none for axis %d of the per-axis field specifications "space directions:".\n' ...
+                            ' There should be no quotation marks, parentheses or any other characters, just plain none.\n'],...
+                            spacedir_vecs{i}, i);
+                        % Clean none vector specification
+                        spacedir_vecs{i} = 'none';
+                    end
+                else
+                    % Axis-specific entry is a numerical vector 
+                    cnt_space_vectors = cnt_space_vectors + 1;
+                    if cnt_space_vectors > internal_spacedimension
+                        error(['%s:\n ''space directions'' field specifier: ' ...
+                            'number of space vectors detected exceeds space (world)' ...
+                            ' dimension, which is equal to %d.'],...
+                            fopen(fidr), internal_spacedimension)
+                    end
+                    btw_parentheses = spacedir_vecs{i}(1) == '(' ...
+                                        && spacedir_vecs{i}(end) == ')';
+                    axis_vector = regexprep(spacedir_vecs{i}, '[()]', ''); % stripped off all parens
+                    if ~btw_parentheses
+                        fprintf(['WARNING nhdr_nrrd_read: vector should be delimited ' ...
+                                'by parentheses for axis %d of the per-axis field ' ...
+                                'specifications "space directions:".\n' ...
+                                ' At least one missing parenthesis in ''%s''.\n'],...
+                                i, spacedir_vecs{i})
+                    end
+                    % Clean up by forcing single enclosing parentheses:
+                    spacedir_vecs{i} = ['(',  axis_vector,  ')'];
+                    % Check vector and extract numerical data
+                    vector_entries = strsplit(axis_vector, ',');
+                    if length(vector_entries) ~= internal_spacedimension
+                        error(['%s:\n vector for data axis %d (space axis %d) of the ' ...
+                                'per-axis field specifications "space directions:" should ' ...
+                                'contain %d entries corresponding to the space (or world) dimension ' ...
+                                'specified in the "space" or "space dimension" field specification.' ...
+                                ' Found %d here.'],...
+                                fopen(fidr), i, cnt_space_vectors, internal_spacedimension, ...
+                                length(vector_entries))                            
+                    end
+                    for j = 1:length(vector_entries)
+                        vector_entry = sscanf(vector_entries{j}, '%f');
+                        if isempty(vector_entry)
+                            error(['%s\n in field specification "space directions:",  ' ...
+                                'vector for data axis %d (space axis %d) too short. ' ...
+                                'Detected %d entries instead of expected %d corresponding to ' ...
+                                'space (world) dimension.'],...
+                                fopen(fidr), i, cnt_space_vectors, j-1, internal_spacedimension)
+                        end
+                        SD_data(j, cnt_space_vectors) = vector_entry;
+                    end
+                end
+            end
+            % Store array of cleaned up strings
+            headerInfo.spacedirections = spacedir_vecs;               % cell array of strings, ideally of the form {'(1,0,0)' '(0,2,0.1)' 'none'  '(0.1,0.1,1.1)'} if internal_spacedimension==3
+
+            % Extract numerical data more leniently:
+            SD_data_chk = strrep(space_dir_tmp, 'none', '' );                       % ignore "none" entries
+            SD_data_chk = extractNumbersWithout(SD_data_chk, {'(',')',',', '"', ''''} );           % detects numbers up to first non numerical entry
+            if numel(SD_data_chk) ~= (internal_spacedimension)^2
+                   error(['Expected ''space directions'' to specify a %d-by-%d matrix ' ...
+                       '(%d elements in total) in agreement with world space dimension.' ...
+                       ' Found %d element(s) instead.\n'],...
+                            internal_spacedimension, internal_spacedimension,...
+                            (internal_spacedimension)^2, numel(SD_data_chk));
+            end
+
+            % Sanity check
+            if ~isequal(SD_data_chk(:), SD_data(:))
+               error(['%s:\n ''space directions'' field specifier: couldn''t' ...
+                   ' read space vectors. Please refer to the NRRD format definition.'],...
+                   fopen(fidr)) 
+            end
+            headerInfo.spacedirections_matrix = SD_data;
+            % Correctness of field specification is checked below after 
+            % whole header is parsed because we need to be sure that the 
+            % "dimension" basic field specification was set
+
+        elseif foundKeyword('SPACE UNITS:', cs )
+            % Always optional, must come after space or space dimension
+            if define_orientation ~= 1
+                error('Field specification ''space units'' cannot be specified before ''space'' or ''space dimension''.')
+            end
+            
+            space_units_tmp = strrep( cs(length('SPACE UNITS:')+1:end), 'none', '');        % ignore none entries
+            % FIXME:  ideally, check that the sum of elements including none and
+            % " " matches headerInfo.dimension, i.e. the total dimension, as
+            % specified in the standard. Standard a bit unclear: should
+            % unknown units be specified with "???", "none" or "" ? (quotes
+            % seem to be required). 
+            
+            headerInfo.spaceunits = extract_spaceunits_list( space_units_tmp );
+            
+            if length(headerInfo.spaceunits) ~= internal_spacedimension
+               error(['Expected ''space units'' to contain %d elements enclosed in double quotes' ...
+                   ' to match the ''space'' or ''space dimension'' field but found the following ' ...
+                   '%d element(s):\n%s'],...
+                   internal_spacedimension, length(headerInfo.spaceunits), sprintf('%s\t',headerInfo.spaceunits{:})) 
+            end
+            
+        elseif foundKeyword('SPACE ORIGIN:', cs )
+            % Always optional, must come after space or space dimension
+            
+            assert(define_orientation==1, ...
+                    sprintf(['%s: field ''space origin'' cannot be specified' ...
+                            ' before ''space'' or ''space dimension''.'], ...
+                            fopen(fidr)));
+            
+            iSO = extractNumbersWithout( cs(length('SPACE ORIGIN:')+1:end), {'(',')',','} );
+            
+            assert(numel(iSO) == internal_spacedimension,...
+                    sprintf(['%s: expected ''space origin'' to specify a ' ...
+                            '%d-element vector to match the ''space'' or ' ...
+                            '''space dimension'' field but found %d ' ...
+                            'element(s).'],...
+                            fopen(fidr),internal_spacedimension, numel(iSO)) );
+            
+            headerInfo.spaceorigin = iSO(:);
+            
+        elseif foundKeyword('MEASUREMENT FRAME:', cs )
+            % Always optional, must come after space or space dimension
+            
+            assert(define_orientation==1,...
+                    sprintf(['%s: field ''measurement frame'' cannot be ' ...
+                            'specified before ''space'' or ''space dimension''.'],...
+                            fopen(fidr)));
+            
+            measframe_str = strrep( cs(length('MEASUREMENT FRAME:')+1:end), 'none', '');
+            
+            iMF = extractNumbersWithout( measframe_str, {'(',')',','} ); % fails if non number entries are not previously removed
+            
+            assert(numel(iMF) == (internal_spacedimension)^2,...
+                sprintf(['%s: expected ''measurement frame'' to specify a ' ...
+                        '%d-by-%d matrix (%d total elements) but found %d ' ...
+                        'element(s).'],...
+                        fopen(fidr), internal_spacedimension, ...
+                        internal_spacedimension, (internal_spacedimension)^2,...
+                        numel(iMF)));
+            
+            headerInfo.measurementframe = reshape(iMF(:), [internal_spacedimension, internal_spacedimension]);
+            
+            % FIXME: this will gladly accept '1 0 0 0 1 0 0 0 1' instead of
+            % '(1,0,0) (0,1,0)  (0,0,1)' for instance. But it should have
+            % length spacedimension according to standard so this is not too bad.
+            
+        elseif foundKeyword('THICKNESSES:', cs )
+            
+            sThicknesses = extractStringList( cs(length('THICKNESSES:')+1:end) ); % fixed bug with extractStringList where 2 entries are present
+            iThicknesses = [];
+            lenThicknesses = length( sThicknesses );
+            for iI=1:lenThicknesses
+                iThicknesses = [iThicknesses, str2double(sThicknesses{iI}) ];
+            end
+            headerInfo.thicknesses = iThicknesses;
+            
+        elseif foundKeyword('CENTERINGS:', cs )
+            
+            headerInfo.centerings = extractStringList( cs(length('CENTERINGS:')+1:end ) ); % fixed bug with extractStringList where 2 entries are present
+            
+        elseif foundKeyword('LINE SKIP:', cs) || foundKeyword('LINESKIP', cs)
+            
+            if foundKeyword('LINE SKIP:', cs)
+                headerInfo.lineskip = sscanf( cs( length('LINE SKIP:')+1:end ), '%d' );
+            else
+                headerInfo.lineskip = sscanf( cs( length('LINESKIP:')+1:end ), '%d' );
+            end
+            assert(headerInfo.lineskip >= 0,...
+                    sprintf(['Field lineskip or line skip should be greater' ...
+                            ' than or equal to zero, detected %d.'], ...
+                            headerInfo.lineskip));
+
+        elseif foundKeyword('BYTE SKIP:', cs) || foundKeyword('BYTESKIP:', cs)
+            
+            if foundKeyword('BYTE SKIP:', cs)
+                headerInfo.byteskip = sscanf( cs( length('BYTE SKIP:')+1:end ), '%d' );
+            else
+                headerInfo.byteskip = sscanf( cs( length('BYTESKIP:')+1:end ), '%d' );
+            end
+            assert(headerInfo.byteskip >= -1, ...
+                    sprintf(['Field byteskip or byte skip can only take ' ...
+                            'non-negative integer values or -1, detected %d.'], ...
+                            headerInfo.byteskip));
+            
+        elseif foundKeyword('MODALITY', cs )
+            
+            headerInfo.modality = strtrim( extractKeyValueString( cs(length('MODALITY')+1:end ) ) );
+            
+        elseif foundKeyword('DWMRI_B-VALUE', cs )
+            
+            headerInfo.bvalue = str2double( extractKeyValueString( cs(length('DWMRI_B-VALUE')+1:end ) ) );
+            
+        elseif foundKeyword('DWMRI_GRADIENT_', cs )
+            
+            [iGNr, dwiGradient] = extractGradient(cs(length('DWMRI_GRADIENT_')+1:end ));
+            headerInfo.gradients(iGNr+1,:) = dwiGradient;
+            % FIXME: make it more general than numbering limited to 4 digits?
+            
+        elseif foundKeyword('DATA FILE:', cs ) || foundKeyword('DATAFILE:', cs)
+            % This tells us it is a detached header and ends it. 3 possibilities here
+            
+            field_value = strtrim( cs(length('DATA FILE:')+1:end) );
+            if foundKeyword('DATAFILE:', cs)
+                field_value = strtrim( cs(length('DATAFILE:')+1:end) );
+            end
+            
+            [filelist, LIST_mode, subdim] = extract_datafiles(field_value);
+            
+            
+            headerInfo.datafiles = filelist;
+            
+            % In LIST mode, filelist is empty and is filled by reading the
+            % rest of the header
+            if LIST_mode
+                datafile_cnt = 0;
+                while ~feof(fidr)
+                    cs = fgetl(fidr);
+                    if isempty(cs)
+                        break;
+                    end
+                    datafile_cnt = datafile_cnt + 1;
+                    headerInfo.datafiles{datafile_cnt} = cs;
+                end
+            end
+            
+            % We could technically break the loop here bc data file/datafile is
+            % supposed to close a header; however I think it does no harm to
+            % keep parsing (in LIST_mode, end of file reached anyways)
+            
+        else
+            
+            % see if we are dealing with a comment
+            
+            csTmp = strtrim( cs );
+            if csTmp(1)~='#' && ~strcmp(cs(1:4),'NRRD')
+                fprintf('WARNING nhdr_nrrd_read: Could not parse input line: ''%s'' \n', cs );
+                % REVIEW: is it better to blindly write it to the output
+                % structure in order to be able to write the exact same file
+                % later on ?
+            end
+        end
+        
+    end
+    
+    
+    % Check for required fields 
+    % REVIEW: should this be checked only when the data is read? People 
+    % might be interested in just parsing the header no matter what is 
+    % in it...)
+    assert(isfield(headerInfo, 'sizes'), 'Missing required ''sizes'' field in header');
+    assert(isfield(headerInfo, 'dimension'), 'Missing required ''dimension'' field in header');
+    assert(isfield(headerInfo, 'encoding'), 'Missing required ''encoding'' field in header');
+    assert(isfield(headerInfo, 'type'), 'Missing required ''type'' field in header');
+
+    % Check other cross-field dependencies, etc. 
+    % TODO: assert lengths of all
+    % fields specified, check that kinds and axis sizes match, etc.
+    if isfield(headerInfo, 'byteskip') && headerInfo.byteskip == -1
+        assert( strcmpi(headerInfo.encoding, 'raw'), ...
+            sprintf('byte skip value of -1 is only valid with raw encoding. See definition of NRRD File Format for more information.\n'));
+    end
+           
+    matlabdatatype = nrrd_getMatlabDataType(headerInfo.type);
+    
+    %  endian field required if data defined on 2 bytes or more
+    if ~(any(strcmpi(headerInfo.encoding,{'txt', 'text', 'ascii'})) || any(strcmpi(matlabdatatype, {'int8', 'uint8'})))
+        assert(isfield(headerInfo, 'endian'), 'Missing required ''endian'' field in header');
+    else
+        if ~isfield(headerInfo,'endian')
+            headerInfo.endian = 'little';    % useless but harmless, just for code compatibility because endian field may be accessed later on while reading data
+        end
+    end
+    
+    % Space and orientation information
+    if define_orientation
+        assert(isfield(headerInfo, 'spacedirections'), ...
+                ['Missing field ''space directions'', required if either' ...
+                ' ''space'' or ''space dimension'' is set.']);
+        if length(headerInfo.spacedirections) ~= headerInfo.dimension
+            fprintf(['WARNING nhdr_nrrd_read %s:\n Unexpected format found for ''space directions'' specifier.\n',...
+                    ' Expected none entries and vectors delimited by parentheses such as (1.2,0.2,0.25), ', ...
+                    'the number of entries being equal to the dimension field specification.\n',...
+                    ' See definition of NRRD File Format for more information.\n'], fopen(fidr));
+        end
+    end
+    
+    % TODO: add support for positive line skips 
+    if isfield(headerInfo, 'lineskip') && headerInfo.lineskip > 0
+        assert(isfield(headerInfo, 'byteskip') && headerInfo.byteskip == -1, ...
+            sprintf(['lineskip option is currently not supported and can ' ...
+                    'only be set to zero unless raw encoding is used and ' ...
+                    'byte skip is set to -1 (which cancels the effect of ' ...
+                    'lineskip altogether).']));
+    end  
+    
+    % TODO: add support for positive byte skips (see line skip)
+    if isfield(headerInfo, 'byteskip') && ~strcmpi(headerInfo.encoding,'raw')
+        assert( headerInfo.byteskip == 0, ...
+            sprintf('byte skip option with non raw encoding is currently not supported and can only be set to zero.\n'));
+    end
+    if isfield(headerInfo, 'byteskip') && strcmpi(headerInfo.encoding, 'raw')
+        assert( headerInfo.byteskip == -1, ...
+            sprintf('non-negative byte skip values with raw encoding are currently not supported; byte skip can only be set to -1.\n'));
+    end
+    
+catch me
+    % Clean up before raising error
+    fclose(fidr);
+    rethrow(me);
+end
+
+% Read the data. Detect if data is in the file or in a detached data file
+% and check what type of detached file it is if need be.
+if bReadData
+    N_data_tot = prod(headerInfo.sizes);
+    if isfield(headerInfo, 'datafiles')
+        % This is a detached header file
+        
+        % We no longer need to read the header (closing it before reading
+        % all the detached data files is safer)
+        fclose(fidr);           
+        
+        % TODO: we currently only read slices along the slowest axis (i.e.,
+        % last coordinate)
+        if ~isempty(subdim) && subdim~=headerInfo.dimension
+            error(['(detached header): reading data from slices along axis other than the last' ...
+                ' (i.e. slowest) one, is currently not supported.\n' ...
+                'Last argument [<subdim>] in ''data file'' or ''datafile'' field ' ...
+                'should be removed or set equal to %d, which is the detected' ...
+                ' ''dimension'' field value, for now.'],...
+                        headerInfo.dimension);
+        end
+        
+        % Read data chunk by chunk from detached data files
+        N_data_files = length(headerInfo.datafiles);
+        assert(mod(N_data_tot, N_data_files)==0, ...
+            sprintf(['Number of detected data files (%d) does not divide total' ...
+                    ' number of values contained in data %d obtained from prod(sizes=[%s]).\n'],...
+                    N_data_files, N_data_tot, sprintf('%d ',headerInfo.sizes)));
+        
+        N_data_per_file = N_data_tot/N_data_files;
+        
+        headerInfo.data = zeros(headerInfo.sizes, matlabdatatype);      % specify right type of zeros, otherwise double by default
+        
+        for i = 1:N_data_files
+            
+            % Check type of detached data file
+            [~,fname_data,ext_data] = fileparts(headerInfo.datafiles{i});        % redundant because done above
+            
+            data_ind = (i-1)*N_data_per_file+1:i*N_data_per_file;
+            
+            if strcmpi(ext_data,'.nhdr')
+                
+                error(['datafile %di/%d: nhdr file should not be used as ' ...
+                        'detached data file.'], ...
+                        i, length(headerInfo.datafiles));
+                
+            elseif strcmpi(ext_data, '.nrrd')
+                % Detached NRRD file
+                
+                bRead_Detached_Data = true;
+                tmp_struct = nhdr_nrrd_read(fullfile(mainFpath, ...
+                                            [fname_data, ext_data]), ...
+                                            bRead_Detached_Data);   % recursive call
+                % TODO: check the rest of the structure for
+                % inconsistencies with metadata from header
+                assert(N_data_files==1 || headerInfo.dimension == tmp_struct.dimension +1, ...
+                    sprintf(['Detached header %s: the number of dimensions in detached ' ...
+                            'nrrd data file %d/%d (%s) should be one fewer than that of' ...
+                            ' the header file.\nDetected %d instead of %d.\nDifferent ' ...
+                            'datafile dimensions as specified by the nrrd standard are' ...
+                            ' not supported as of now.\n'],...
+                            fopen(fidr), i, N_data_files, headerInfo.datafiles{i},...
+                            tmp_struct.dimension, headerInfo.dimension-1 ));
+                
+                headerInfo.data(data_ind) = tmp_struct.data(:);
+                
+                % Store detached data header for last detached file seen,
+                % assuming that all are the same. Do not store data though
+                % as it would be redundant.
+                tmp_struct = rmfield(tmp_struct, 'data');
+                headerInfo.detached_header = tmp_struct;
+
+            else
+                % e.g., detached .raw file
+                
+                fid_data = fopen( fullfile(mainFpath, [fname_data, ext_data]), 'r');
+                if( fid_data < 1 )
+                    error(['While reading detached header file %s:\ndetached ' ...
+                           'data file number %d/%d (%s) could not be opened.'], ...
+                           nrrdFileName, i, N_data_files, headerInfo.datafiles{i});
+                end
+                try
+                    tmp_data = readData(fid_data, N_data_per_file, headerInfo.encoding, matlabdatatype);
+                    fclose(fid_data);
+                catch me_detached
+                    fclose(fid_data);
+                    rethrow(me_detached);
+                end
+                
+                tmp_data = adjustEndian(tmp_data, headerInfo.endian);
+                
+                headerInfo.data(data_ind) = tmp_data(:);
+            end
+            
+        end
+        
+    else
+        % This is a NRRD standalone file: read data directly from it
+        try
+            headerInfo.data = readData(fidr, N_data_tot, headerInfo.encoding, matlabdatatype);
+            fclose(fidr);
+        catch me_detached
+            fclose(fidr);
+            rethrow(me_detached);
+        end        
+        headerInfo.data = adjustEndian(headerInfo.data, headerInfo.endian);
+        headerInfo.data = reshape(headerInfo.data, headerInfo.sizes(:)');   % data into expected form. Transpose required by reshape function bc size vectors need to be row vectors
+    end
+else
+    fclose(fidr);
+end
+
+
+
+
+end
+
+% -------------------------------------------------------------------%
+
+% ====================================================================
+% --- Auxiliary functions -------------------------------------------%
+% ====================================================================
+
+function [iGNr, dwiGradient] = extractGradient( st )
+
+% first get the gradient number
+
+iGNr = str2num( st(1:4) ); % FIX numbering limited to 4 digits (synchronize with writer module)
+
+% find where the assignment is
+
+assgnLoc = strfind( st, ':=' );
+
+if ( isempty(assgnLoc) )
+    dwiGradient = [];
+    return;
+else
+    
+    dwiGradient = sscanf( st(assgnLoc+2:end), '%f' );
+    
+end
+
+end
+
+% Return part of string after :=
+function kvs = extractKeyValueString( st )
+
+assgnLoc = strfind( st, ':=' );
+
+if ( isempty(assgnLoc) )
+    kvs = [];
+    return;
+else
+    
+    kvs = st(assgnLoc(1)+2:end);
+    
+end
+
+end
+
+% Turn space-separated list into cell array of strings. 
+function sl = extractStringList( strList )
+sl = strsplit(strtrim(strList)); % old Matlab file exchange version had a strange bug with lists of length 2
+end
+
+
+% Store in an array the list of numbers separated by the tokens listed in
+% the withoutTokens cell array 
+function iNrs = extractNumbersWithout( inputString, withoutTokens )
+
+auxStr = inputString;
+
+for iI=1:length( withoutTokens )
+    
+    auxStr = strrep( auxStr, withoutTokens{iI}, ' ' );
+    
+end
+
+iNrs = sscanf( auxStr, '%f' );
+
+end
+
+% Return true if the string keyWord is the beginning of the content string
+% cs (ignoring case)
+function fk = foundKeyword( keyWord, cs )
+lenKeyword = length( keyWord );
+fk = (lenKeyword <= length(cs)) && strcmpi( cs(1:lenKeyword), keyWord);
+end
+
+% Return cell array of strings with space units in the form {mm, km, m}
+% from input of the form "mm " "mm" " m" where undesired blank spaces may
+% have crept in. The double quotes are removed in the process but will be
+% added by the nhdr/nrrd writer function.
+function su_ca = extract_spaceunits_list( fieldValue )
+fv_trimmed = strtrim( fieldValue );
+su_ca = strsplit(fv_trimmed, '"');                              % units are delimited by double quotes
+su_ca = su_ca(~ ( strcmp(su_ca, '') | strcmp(su_ca, ' ') ) );   % remove empty or blank space strings
+for i = 1:length(su_ca)
+    su_ca{i} = strtrim( su_ca{i} );
+end
+end
+
+
+% Extract data file list into a cell array from the datafile or data file
+% field in a NHDR file, stripped off its leading and trailing spaces;
+% LIST_mode is set to one if a list of files closes the header;
+% subdim is empty by default and may be set through either of these field
+% specifications:
+% data file: <format> <min> <max> <step> [<subdim>]
+% data file: LIST [<subdim>]
+
+function [filelist, LIST_mode, subdim] = extract_datafiles( field_string)
+
+field_string = strtrim(field_string);
+filelist = {};
+LIST_mode = 0;
+subdim = [];
+if length(field_string)>= 4 && strcmpi(field_string(1:4), 'LIST')
+    % LIST mode
+    subdim = sscanf(field_string(5:end), '%d'); % assert 1<=dimensions 
+    LIST_mode = 1;
+    return;
+else
+    % single detached data file or multiple files written in concise form, non LIST mode
+    str_lst = strsplit(field_string); % without delimiters specified, splits at any sequence in the set  {' ','\f','\n','\r','\t','\v'}
+    if length(str_lst) == 1
+        % Single detached data file:
+        filelist{1} = str_lst{1};
+        % subdim does not make sense here since all of the data is
+        % contained in the deatached data file
+    elseif any(length(str_lst)==[0, 2, 3]) || length(str_lst)>5 
+        % Invalid cases (2, 3 or striclty more than 5 entries)
+        error(['error in ''data list'' (or ''data list'') field, in non' ...
+                ' ''LIST'' mode:\nexpected  ''<filename>'' or ''<format> ' ...
+                '<min> <max> <step> [<subdim>]'' but found instead ''%s'' ' ...
+                '(%d elements instead of 1, 4 or 5).'],...
+                sprintf('%s ',str_lst{:}), length(str_lst));
+    else 
+        % Multiple detached data files with filenames including integral
+        % values generated by min:step:max
+                
+        str_format = str_lst{1};
+        id_min = sscanf(str_lst{2}, '%d');
+        id_max = sscanf(str_lst{3}, '%d');
+        step = sscanf(str_lst{4}, '%d');
+        
+        % check format and indexing values
+        assert(step~=0,...
+                sprintf(['detached data files with names specified by ' ...
+                        'sprintf()-like format:  step should be strictly positive' ...
+                        ' or negative, not zero.']));
+        if step < 0
+            assert(id_min >= id_max, ...
+                    sprintf(['detached data files with names specified by ' ...
+                            'sprintf()-like format:  when step is <0, min ' ...
+                            'should be larger than or equal to max, here we ' ...
+                            'found min=%d < max=%d.'],id_min, id_max));
+        else
+            assert(id_min <= id_max,...
+                    sprintf(['detached data files with names specified by ' ...
+                            'sprintf()-like format: : when step is >0, min ' ...
+                            'should be smaller than or equal to max, here we ' ...
+                            'found min=%d > max=%d.'],id_min, id_max));
+        end
+         
+        % populate data file list
+        fileIDs = id_min:step:id_max;
+        filelist = cell(length(fileIDs), 1);
+        for i = 1:length(fileIDs)  % does not necessarily incude <max>, it cannot be larger than <max> 
+            expanded_fname = sprintf(str_format, fileIDs(i));
+            filelist{i} = expanded_fname;
+        end
+        
+        if length(str_lst) == 5
+            subdim = sscanf(str_lst{5}, '%d');
+        end   
+    end
+end
+        
+
+end
+
+
+% Read data in a column vector from open file with pointer fidIn at the
+% beginning of the data to read. The other arguments specify the number of
+% elements to read (determined in main body of reader function above), data
+% file encoding (specified in NHDR/NRRD header) and the Matlab type the
+% data should be converted to (determined in main body of reader function
+% above).
+function data = readData(fidIn, Nelems, meta_encoding, matlabdatatype)
+
+switch (meta_encoding)
+    case {'raw'}
+        % This implicitely assumes that byteskip was set to -1 in raw
+        % encoding
+        data_bytes = (Nelems)*sizeOf(matlabdatatype);   % bytes of useful info to return in data array
+        fseek(fidIn,0,'eof');                           % set position pointer to end of file
+        tot_bytes_file = ftell(fidIn);                  % current location of the position pointer (I believe as a byte offset from beginning of file)
+        fseek(fidIn,tot_bytes_file-data_bytes,'bof');   % make sure position pointer is right before useful data starts
+        
+        % Just read binary file (original version of the code)
+        data = fread(fidIn, inf, [matlabdatatype '=>' matlabdatatype]);
+        
+    case {'gzip', 'gz'}
+        
+        tmp = fread(fidIn, Inf, 'uint8=>uint8'); % byte per byte
+        
+        tmpBase = tempname(pwd);
+        tmpFile = [tmpBase '.gz'];
+        fidTmp = fopen(tmpFile, 'wb');  % this creates tmpFile
+        assert(fidTmp > 3, 'Could not open temporary file for GZIP decompression')
+        
+        try
+            fwrite(fidTmp, tmp, 'uint8');
+        catch me
+            fclose(fidTmp);
+            delete(tmpFile);
+            rethrow(me);
+        end
+        fclose(fidTmp);        
+        
+        try
+            gunzip(tmpFile); % this creates tmpBase (tmpFile without the .gz)
+        catch me
+            delete(tmpFile);
+            rethrow(me);
+        end
+        delete(tmpFile);
+        
+        fidTmp = fopen(tmpBase, 'rb');
+        % cleaner = onCleanup(@() fclose(fidTmp));
+        try
+            data = readData(fidTmp, Nelems, 'raw', matlabdatatype); % recursive call            
+        catch me
+            fclose(fidTmp);
+            delete(tmpBase);
+            rethrow(me);
+        end
+        fclose(fidTmp);
+        delete(tmpBase);
+
+    case {'txt', 'text', 'ascii'}
+        
+        data = fscanf(fidIn, '%f');
+        data = cast(data, matlabdatatype);
+        
+    otherwise
+        assert(false, 'Unsupported encoding')
+end
+% Check number of read elements (sometimes there is an offest of about 1)
+assert(Nelems == numel(data),...
+        sprintf(['Error reading binary content of %s: detected %d elements' ...
+                ' instead of %d announced in header.'],...
+                fopen(fidIn), numel(data), Nelems));
+end
+
+
+% Switch byte ordering according to endianness specified in nhdr/nrrd file.
+% The metadata should just contain a field 'endian' set to 'little' or
+% 'big'.
+function data = adjustEndian(data, meta_endian)
+
+[~,~,endian] = computer();
+
+needToSwap = (isequal(endian, 'B') && isequal(lower(meta_endian), 'little')) || ...
+    (isequal(endian, 'L') && isequal(lower(meta_endian), 'big'));
+
+if (needToSwap)
+    data = swapbytes(data);
+end
+end
+
+% Size of Matlab's numeric classes in bytes
+% See: https://stackoverflow.com/questions/16104423/size-of-a-numeric-class
+
+function numBytes = sizeOf(dataClass)
+
+    % create temporary variable of data type specified
+    eval(['var = ' dataClass '(0);']); 
+
+    % Use the functional form of whos, and get the number of bytes   
+    W = whos('var');
+    numBytes = W.bytes;
+
+end

matlab/nrrd_read_write_rensonnet/nhdr_nrrd_write.m

@@ -0,0 +1,646 @@
+% Write NHDR/NRRD files
+% 
+% headerInfo = nhdr_nrrd_write(nrrdFileName, headerInfo, bWriteData)
+%
+% Inputs:
+% * nrrdFilename (char): path to NHDR/NRRD file(s), either a detached header,
+% a detached header pointing to detached data files or a nrrd standalone
+% depending on what is specified in headerInfo
+% * headerInfo (struct): structure containing all the required and optional
+% NRRD fields, such as produced by nhdr_nrrd_read.m, as well as the data. 
+% List of accepted structure field names with corresponding NRRD fields in
+% parentheses:
+%   - byteskip (byte skip: <desc>): scalar integer. For now, can only
+%       accept -1 with raw encoding and 0 with all other encodings;
+%	- content (content: <desc>): string;
+%	- data (~): numerical array containing the raw data;
+%	- datafiles (data file: <desc>): cell array of strings containing the
+%       paths to detached data files relative to the header file. This can 
+%       be a simple string if there is only one detached data file;
+% 	- detached_header (~): Matlab structure containing all the fields 
+%       required to write valid detached NRRD files containing the data. 
+%       All detached NRRD data files are assumed to have the same header;
+% 	- dimension (dimension: <desc>): scalar integer;
+%   - encoding (encoding: <desc>): string;
+%   - kinds (kinds: <desc>): cell array of strings;
+%   - lineskip (line skip: <desc>): scalar integer. For now, can only be 
+%       set to zero unless raw encoding is used and byte skip is set to -1
+%       (which cancels the effect of line skip altogether);
+%   - measurementframe (measurement frame: <desc>): matrix of size 
+%       [Ns, Ns], where Ns in the number of space dimensions;
+%   - sizes (sizes: <desc>): matrix of size [1, d] or [d, 1] where d is the
+%       dimension;
+%   - space (space: <desc>): string;
+%   - spacedimension (space dimension: <desc>): scalar integer;
+%   - spacedirections (space directions: <desc>): cell array of strings 
+%       for code compatibility, e.g. 
+%       {'(1,0,0)' '(0,2,0.1)' 'none'  '(0.1,0.1,1.1)'} for 3 
+%       world space dimensions;
+%   - spacedirections_matrix (~): matrix of size [Ns, Ns] where Ns is the 
+%       number of space dimensions. Should match the string description 
+%       contained in spacedirections, if specified;
+%   - spaceorigin (space origin: <desc>): matrix of size [1, Ns] or [Ns, 1]
+%       where Ns is the number os space dimensions;
+%   - spaceunits (space units: <desc>): cell array of strings;
+%   - type (type: <desc>): string.
+%  
+%   - bvalue (bvalue:=<value>): floating-point scalar (nominal b-value);
+%   - modality (modality:=<value>): string
+%   - gradients (DWMRI_gradient_%04d:=<value>): matrix of size [Ng, 3] 
+%       where Ng is the number of gradients in a PGSE-type 
+%       diffusion-weighted MRI protocol. Note that the b-value associated 
+%       with each gradient [g_x, g_y, g_z] is computed as 
+%       b = b_nom*(g_x^2 + g_y^2 + g_z^2), where b_nom is the
+%       nominal b-value, so the gradients must be scaled accordingly;
+% * bWriteData (bool): set to true to write the data specified in
+% headerInfo.data, set to false to just write the header file without the
+% data
+%
+% Outputs:
+% * headerInfo (struct): same structure as that passed as argument with
+% additional (potentially mandatory) fields if inferrable from the data,
+% such as 'dimension', 'sizes', etc.
+%
+% The function creates all the required NHDR/NRRD file(s).
+%
+% See nhdr_nrrd_read.m for more information.
+%
+%
+% % Format definition:
+% http://teem.sourceforge.net/nrrd/format.html
+%
+% Date: August 2018
+% Author: Gaetan Rensonnet
+function headerInfo = nhdr_nrrd_write(nrrdFileName, headerInfo, bWriteData)
+
+[path_main_header,  ~, fext] = fileparts(nrrdFileName);
+if isfield(headerInfo, 'datafiles')
+    assert(strcmpi(fext, '.nhdr'), sprintf('Invalid filename extension ''%s''. A detached nrrd header (headerInfo structure contains ''datafiles'' field) should have extension ''.nhdr''.\n',fext));
+else
+    assert(strcmpi(fext, '.nrrd'), sprintf('Invalid filename extension ''%s''. A standalone nrrd file (no ''datafiles'' field in headerInfo structure) should have extension ''.nrrd''.\n',fext));
+end
+
+header_fnames = fieldnames(headerInfo);             % cell array of strings
+fnames_parsed = false(1,length(header_fnames));     % flags field names from header successfully parsed
+
+% --- Checking availability of required information before opening files ---
+% Guess as much from the data as possible
+
+% Check required 'sizes' information
+if isfield(headerInfo, 'sizes')
+    if isfield(headerInfo, 'data') 
+        assert(nrrd_size_check(size(headerInfo.data), headerInfo.sizes),...
+        sprintf('Sizes mismatch: [%s] in headerInfo structure is not compatible with the size of the data array [%s].\n',...
+                                    sprintf('%d ',headerInfo.sizes), sprintf('%d ',size(headerInfo.data)) ));
+    end
+else
+    assert( isfield(headerInfo, 'data'), sprintf('Missing required field ''sizes'' in headerInfo structure, cannot be deduced from data because no data was provided.\n'));
+    headerInfo.sizes = size(headerInfo.data);
+end
+
+% Check required 'dimension' information
+if isfield(headerInfo, 'dimension')
+    assert( isequal(headerInfo.dimension, length(headerInfo.sizes)), sprintf('Dimension mismatch: %d in headerInfo structure is not compatible with detected data size [%s].\n',...
+                                            headerInfo.dimension, sprintf('%d ', headerInfo.sizes)));
+else
+    headerInfo.dimension = length(headerInfo.sizes);
+end
+
+% Check required 'type' information
+if isfield(headerInfo, 'type')
+    matlabdatatype = nrrd_getMatlabDataType(headerInfo.type);
+    if isfield(headerInfo, 'data')
+        assert( isequal(class(headerInfo.data), matlabdatatype), ...
+            sprintf('Type mismatch: %s in headerInfo structure is not compatible with the Matlab type of the data array (%s).\nYou may want to look at nrrd_getMatlabDataType and at Matlab''s cast functions.',headerInfo.type, class(headerInfo.data)));
+    end
+else
+   assert( isfield(headerInfo, 'data'), sprintf('Missing required field ''type'' in headerInfo structure, cannot be deduced from data because no data was provided.\n'));
+   headerInfo.type = get_nrrd_datatype(class(headerInfo.data));
+end
+
+% Check required 'encoding' information
+if ~isfield(headerInfo, 'encoding')
+    headerInfo.encoding = 'raw';
+end
+
+% Set endianness to current platform's endianness
+[~,~,endian] = computer();
+if (isequal(endian, 'B'))
+    headerInfo.endian = 'big';
+else
+    headerInfo.endian = 'little';
+end
+
+    
+fidw = fopen( nrrdFileName, 'w');
+if fidw == -1
+    error('nhdr_nrrd_writer: could not open file %s.\n',nrrdFileName);
+end
+
+% if an error is thrown while parsing, close open file before exiting
+try     
+    % --- Traditional header and required fields ---
+
+    fprintf( fidw, 'NRRD0005\n' );  
+    fprintf( fidw, '# Complete NRRD file format specification at:\n' );
+    fprintf( fidw, '# http://teem.sourceforge.net/nrrd/format.html\n' );
+
+    if isfield(headerInfo,'content')
+        fprintf( fidw, 'content: %s\n', headerInfo.content);  % always optional
+        fnames_parsed(strcmpi('content', header_fnames)) = true;
+    end
+
+    fprintf( fidw, 'type: %s\n', headerInfo.type );
+    fnames_parsed(strcmpi('type', header_fnames)) = true;
+
+    fprintf( fidw, 'encoding: %s\n', headerInfo.encoding);
+    fnames_parsed(strcmpi('encoding', header_fnames)) = true;
+
+    fprintf( fidw, 'dimension: %d\n', headerInfo.dimension );
+    fnames_parsed(strcmpi('dimension', header_fnames)) = true;
+
+    fprintf( fidw, 'sizes: ' );
+    for iI=1:length( headerInfo.sizes )
+      fprintf( fidw, '%d', headerInfo.sizes(iI) );
+      if ( iI~=length( headerInfo.sizes ) )
+        fprintf( fidw, ' ' );
+      end
+    end
+    fprintf( fidw, '\n' );
+    fnames_parsed(strcmpi('sizes', header_fnames)) = true;
+
+    fprintf( fidw, 'endian: %s\n', headerInfo.endian );
+    fnames_parsed(strcmpi('endian', header_fnames)) = true;
+    
+    % ---- Defining orientation (always optional) -----
+    
+    define_orientation = 0;
+    if isfield(headerInfo,'space')
+        fprintf(fidw, 'space: %s\n', headerInfo.space );
+        fnames_parsed(strcmpi('space', header_fnames)) = true;
+        define_orientation = 1;
+        num_sp_dim = nrrd_getSpaceDimensions(headerInfo.space);
+    elseif isfield (headerInfo,'spacedimension')
+        fprintf(fidw, 'space dimension: %d\n', headerInfo.spacedimension);
+        fnames_parsed(strcmpi('spacedimension', header_fnames)) = true;
+        define_orientation = 1;
+        num_sp_dim = headerInfo.spacedimension;
+    end
+    % Space and space dimension fields cannot be specified simultaneously
+    if isfield(headerInfo,'space') && isfield(headerInfo, 'spacedimension')
+        error('The always optional fields ''space'' and ''spacedimension'' defining orientation cannot be specified simultaneously in headerInfo structure');
+    end
+    
+    if define_orientation
+        % So far the number of spatial dimensions should be 3 or 4
+        % (this is important for writing space directions, space origin,
+        % measurement frame)
+        assert(num_sp_dim>=3, sprintf('Number of space dimensions should be at least 3 in current nrrd format, detected %d instead in headerInfo structure.\n',num_sp_dim));
+        
+        % check for necessary field 'space direction' when orientation defined
+        if ~isfield(headerInfo, 'spacedirections')
+            error('Missing field in headerInfo structure: spacedirections should always be specified if orientation is defined via ''space'' or ''space directions'' ');
+        end
+        
+        assert(iscell(headerInfo.spacedirections), ...
+            ['spacedirections field in headerInfo should be a cell array of ' ...
+            'strings for code compatibility, e.g. {''(1,0,0)'' ''(0,2,0.1)'' ''none''  ''(0.1,0.1,1.1)''}' ...
+            ' for 3 world space dimensions.'])
+
+        % Check format and clean up when possible
+        for i = 1:length(headerInfo.spacedirections)
+            none_match = regexp(headerInfo.spacedirections{i}, 'none', 'match');
+            if ~isempty(none_match)
+                % none for non space axis
+                if ~strcmpi(headerInfo.spacedirections{i}, 'none')
+                   fprintf(['nhdr_nrrd_write WARNING: detected %s instead of ' ...
+                       'expected none in header.spacedirections{%d}.\n Wrote clean version to output file(s).\n'],...
+                       headerInfo.spacedirections{i}, i)
+                end
+                headerInfo.spacedirections{i} = 'none';
+            else
+                % numerical vector for space (world) axis
+                btw_parentheses = headerInfo.spacedirections{i}(1) == '(' ...
+                    && headerInfo.spacedirections{i}(end) == ')';
+                if ~btw_parentheses
+                   fprintf(['nhdr_nrrd_write WARNING: space vector in header.spacedirections{%d} ' ...
+                       'should contain comma-separated values enclosed in single parentheses.\n ' ...
+                       'At least one parenthesis missing here. Wrote clean version to output file(s).\n'],...
+                         i)
+                end
+                headerInfo.spacedirections{i} = ['(' regexprep(headerInfo.spacedirections{i}, '[()]', ''), ')'];
+            end
+        end
+        spacedir_str = strjoin(headerInfo.spacedirections, ' ');
+
+        % Write to file
+        fprintf(fidw, 'space directions: %s\n', spacedir_str);
+        % Flag as parsed
+        fnames_parsed(strcmpi('spacedirections', header_fnames)) = true;
+        
+        % Data defining space direction matrix
+        SD_data = strrep(spacedir_str, 'none', '');
+        SD_data = extractNumbersWithout(SD_data, {'(',')',',', '"', ''''});       % this fails if non numerical entries were not previously removed
+        assert(length(SD_data) == num_sp_dim^2, sprintf('expected spacedirections field to contain %d numbers (square of %d, the world space dimension). Found %d instead.\nspacedirections should be a cell array of strings containing vectors of the form (dx,dy,dz) or none entries.\n',...
+                                                        num_sp_dim^2, num_sp_dim, length(SD_data)))
+        % Check spacedirections_matrix field (internal, not a NRRD field)
+        if isfield(headerInfo, 'spacedirections_matrix')
+            assert(isequal(headerInfo.spacedirections_matrix(:), SD_data(:)), ...
+                sprintf('Numeric data in spacedirections and spacedirections_matrix fields of headerInfo do not match. They should contain %d identical numbers in the same order.\n', num_sp_dim^2))
+        else
+            % add it to output headerInfo structure 
+            headerInfo.spacedirections_matrix = reshape(SD_data(:), [num_sp_dim, num_sp_dim]);  
+        end
+        fnames_parsed(strcmpi('spacedirections_matrix', header_fnames)) = true;
+    end
+    
+    
+    % --- optional options for the optional definition of orientation ---
+    
+    % space origin
+    if isfield(headerInfo,'spaceorigin')
+        if define_orientation
+            so = headerInfo.spaceorigin;
+            assert(length(so)==num_sp_dim, ...
+                sprintf('Field ''spaceorigin'' in headerInfo structure expected vector with %d entries to match the defined orientation, detected %d entries instead.\n',num_sp_dim, length(so)));
+            fprintf( fidw, 'space origin: (%f%s)\n', so(1), sprintf(',%f',so(2:end)) );
+            fnames_parsed(strcmpi('spaceorigin', header_fnames)) = true;
+        else
+            error('Field ''spaceorigin'' in headerInfo structure cannot be specified if neither ''space'' nor ''spacedimension'' was specified.');
+        end
+    end
+    
+    % measurement frame
+    if isfield(headerInfo,'measurementframe')
+        if define_orientation
+            mf = headerInfo.measurementframe;
+            assert(size(mf,1)==num_sp_dim && size(mf,2)==num_sp_dim,...
+                sprintf('Field ''measurementframe'' in headerInfo structure expected a %d-by-%d matrix to match the defined orientation, detected size [%s] instead.\n',...
+                                                 num_sp_dim, num_sp_dim, sprintf('%d ',size(mf))));
+                                             
+            fprintf( fidw, 'measurement frame:');
+            for imf = 1:size(mf,2)
+               fprintf(fidw,' (%f%s)', mf(1,imf), sprintf(',%f',mf(2:end,imf)));                
+            end
+            fprintf(fidw, '\n');
+            fnames_parsed(strcmpi('measurementframe', header_fnames)) = true;
+        else
+            error('Field ''measurementframe'' in headerInfo structure cannot be specified if neither ''space'' nor ''spacedimension'' were specified.');
+        end
+    end
+    
+    % spaceunits
+    if isfield(headerInfo,'spaceunits')
+        if define_orientation
+            fprintf( fidw, 'space units: ' );
+            for iI=1:length( headerInfo.spaceunits )
+                fprintf( fidw, '\"%s\"', char(headerInfo.spaceunits(iI)) );
+                if ( iI~=length( headerInfo.spaceunits ) )
+                    fprintf( fidw, ' ' );
+                end
+            end
+            fprintf( fidw, '\n' );
+            fnames_parsed(strcmpi('spaceunits', header_fnames)) = true;
+        else
+            error('Field ''spaceunits'' in headerInfo structure cannot be specified if neither ''space'' nor ''spacedimension'' were specified.');
+        end
+    end
+    % --- end of orientation definition -----
+    
+    % --- Optional fields "<field>: <desc>" ---- 
+    
+    % kinds
+    if isfield(headerInfo, 'kinds')
+        assert(length(headerInfo.kinds)==headerInfo.dimension, sprintf('Detected %d kinds instead of expected nrrd dimension %d.\n',length(headerInfo.kinds), headerInfo.dimension));
+        fprintf( fidw, 'kinds: ' );
+        for iI=1:length( headerInfo.kinds )
+          fprintf( fidw, '%s', headerInfo.kinds{iI} );
+          if ( iI~=length( headerInfo.kinds ) )
+            fprintf( fidw, ' ' );
+          end
+        end
+        fprintf( fidw, '\n' );
+        fnames_parsed(strcmpi('kinds', header_fnames)) = true;
+    end
+    
+    % line skip
+    if isfield(headerInfo, 'lineskip')
+       assert(headerInfo.lineskip >= 0, sprintf('Field ''lineskip'' in headerInfo structure should be non-negative, detected %d.\n', headerInfo.lineskip));
+       fprintf(fidw, 'lineskip: %d\n', headerInfo.lineskip);
+       fnames_parsed(strcmpi('lineskip', header_fnames)) = true;
+    end
+    
+    % byte skip
+    if isfield(headerInfo, 'byteskip')
+       assert(headerInfo.byteskip >= -1, sprintf('Field ''byteskip'' should be -1 or a non-negative integer, detected %d.\n', headerInfo.byteskip));
+       if headerInfo.byteskip == -1
+           assert(strcmpi(headerInfo.encoding, 'raw'), sprintf('byte skip value of -1 is only valid with raw encoding.\n'));
+       end
+       fprintf(fidw, 'byteskip: %d\n', headerInfo.byteskip);
+       fnames_parsed(strcmpi('byteskip', header_fnames)) = true;
+    end
+    
+    % --- "<key>:=<value>"  pairs (always optional)
+    
+    % Modality 
+    if isfield(headerInfo, 'modality')
+       fprintf( fidw, 'modality:=%s\n',headerInfo.modality);
+       fnames_parsed( strcmpi('modality', header_fnames)) = true;
+    end
+    
+    % b-value 
+    if isfield(headerInfo, 'bvalue')
+        fprintf( fidw, 'DWMRI_b-value:=%9.8f\n', headerInfo.bvalue);
+        fnames_parsed( strcmpi('bvalue', header_fnames)) = true;
+    end
+    
+    % DW-MRI gradient
+    if isfield(headerInfo, 'gradients')
+       Ngrads = size(headerInfo.gradients, 1);
+       for igrad = 1:Ngrads
+           strprint = strcat('DWMRI_gradient_',sprintf('%04d',igrad-1)) ; % The 0 flag in the %04d format specifier requests leading zeros in the output and sets minimum width of the printed value to 4
+           strprint = [strprint,':='] ;
+           strprint = [strprint, sprintf('%9.8f ', headerInfo.gradients(igrad,:))] ;
+           fprintf(fidw,[strprint '\n']) ;       
+       end
+       fnames_parsed(strcmpi('gradients', header_fnames)) = true;
+       % FIXME: make it more general than numbering limited to 4 digits?        
+    end
+       
+    
+    % --- External datafiles (should be performed in LIST mode) ---
+    
+    if isfield(headerInfo, 'datafiles')
+        % This is a detached header, not a standalone nrrd file
+        
+        % 'Convert' to cell array for code compatibility. This only works
+        % with one filename and fails with a comma-separated list of
+        % filenames for instance.
+        if ischar(headerInfo.datafiles)
+            headerInfo.datafiles = { headerInfo.datafiles};            
+        end
+            
+        if length(headerInfo.datafiles)==1
+            % data file: <filename> 
+            fprintf(fidw, 'data file: %s\n',headerInfo.datafiles{1}); % path relative to detached header
+        else
+            % data file: LIST [<subdim>]
+            % FIXME: add support for subdim argument instead of ignoring it
+            fprintf(fidw, 'data file: LIST\n');
+            for i = 1:length(headerInfo.datafiles)
+               fprintf(fidw,'%s\n', headerInfo.datafiles{i}); 
+            end
+        end
+        fnames_parsed(strcmpi('datafiles', header_fnames)) = true;
+    end
+    
+catch me
+    fclose(fidw);
+    rethrow(me);
+end
+
+%% Write data
+if bWriteData
+    
+    if isfield(headerInfo, 'datafiles')
+        % This is a detached header, detached data files need to be
+        % written to
+        
+        fclose(fidw);       % close main file
+        
+        N_data_tot = prod(headerInfo.sizes);
+        
+        % Read data chunk by chunk from detached data files
+        N_data_files = length(headerInfo.datafiles);
+        
+        assert(mod(N_data_tot, N_data_files)==0, sprintf('Number of detected data files (%d) does not divide total number of values contained in data %d obtained from prod(sizes=[%s]).\n',...
+            N_data_files, N_data_tot, sprintf('%d ',headerInfo.sizes)));
+        
+        N_data_per_file = N_data_tot/N_data_files;
+        
+        for i = 1:N_data_files
+            
+            % Check type of detached data file
+            [~,fname_data,ext_data] = fileparts(headerInfo.datafiles{i}); 
+            
+            data_ind = (i-1)*N_data_per_file+1:i*N_data_per_file;   % indices of data to be written
+            
+            if strcmpi(ext_data,'.nhdr')
+                
+                error('datafile %di/%d: nhdr file should not be used as detached data file.\n',i,length(headerInfo.datafiles));
+                
+            elseif strcmpi(ext_data, '.nrrd')                
+                % Detached nrrd file
+                
+                assert(isfield(headerInfo, 'detached_header'), sprintf('Missing field ''detached_header'' in headerInfo structure.\nThis is only required for detached headers with data stored in detached data files of type nrrd.\nShould contain a Matlab structure describing the chunk of data stored in the detached nrrd file (the data field, if present, will be overwritten).\nIt is assumed that the header is identical for all detached nrrd data files.\n'));
+                % we have to check for all i's because detached nrrd files may be
+                % interspersed with detached .raw or .hex files                   
+                info_detached_data = headerInfo.detached_header;        % change headerInfo! put only part of the data, change dimensions, remove datalists
+                fnames_parsed(strcmpi('detached_header', header_fnames)) = true;
+                
+                % Add data chunk to be written to new file
+                info_detached_data.data = headerInfo.data(data_ind);
+                info_detached_data.data = reshape(info_detached_data.data, info_detached_data.sizes);   % for code compatibility
+                
+                % Write detached nrrd file using a recursive call
+                bWrite_Detached_Data = true;
+                tmp_struct = nhdr_nrrd_write(fullfile(path_main_header, [fname_data, ext_data]), info_detached_data, bWrite_Detached_Data);   % recursive call
+                
+            else
+                % e.g., detached .raw or .hex file
+                
+                fid_data = fopen( fullfile(path_main_header, [fname_data, ext_data]), 'w');
+                if( fid_data < 1 )
+                    error('Detached data file number %d/%d (%s) could not be opened.\n', i, N_data_files, headerInfo.datafiles{i});
+                end
+                
+                try
+                    outtype = nrrd_getMatlabDataType(headerInfo.type);
+                    writeData(fid_data, headerInfo.data(data_ind), outtype, headerInfo.encoding);
+                    fclose(fid_data);
+                catch me_detached
+                    fclose(fid_data);
+                    rethrow(me_detached);
+                end
+                
+                
+            end
+            
+        end
+        
+        
+    else
+        % Standalone nrrd file with data included after header        
+        
+        try
+            % After the header, there is a single blank line containing zero
+            % characters to separate it from data
+            fprintf(fidw, '\n');
+            outtype = nrrd_getMatlabDataType(headerInfo.type);
+            writeData(fidw, headerInfo.data, outtype, headerInfo.encoding);
+            fclose(fidw);
+        catch me_data
+            fclose(fidw);
+            rethrow(me_data);
+        end
+        
+    end
+    
+else
+    fclose(fidw);
+end
+
+% Issue warnings for unknown/unsupported fields (after reading data)
+for i = 1:length(header_fnames)
+    if ~fnames_parsed(i) && ~strcmpi(header_fnames{i},'data')
+        fprintf('nhdr_nrrd_write WARNING: could not parse and write field %s from headerInfo structure.\n', header_fnames{i});
+    end
+end
+
+
+end
+
+
+
+% --- Auxiliary functions ------------
+
+% ========================================================================
+% Store in an array the list of numbers separated by the tokens listed in
+% the withoutTokens cell array 
+% ========================================================================
+function iNrs = extractNumbersWithout( inputString, withoutTokens )
+
+auxStr = inputString;
+
+for iI=1:length( withoutTokens )
+    
+    auxStr = strrep( auxStr, withoutTokens{iI}, ' ' );
+    
+end
+
+iNrs = sscanf( auxStr, '%f' );
+
+end
+
+
+% ========================================================================
+% Determine the nrrd datatype : from matlab datatype to outtype (NRRD)
+% ========================================================================
+function nrrddatatype = get_nrrd_datatype(matlab_metaType)
+
+% Determine the datatype
+switch (matlab_metaType)
+ case {'int8', 'uint8', 'int16', 'uint16', 'int32', 'uint32', 'int64',...
+       'uint64', 'double'}
+   nrrddatatype = matlab_metaType;
+  
+ case {'single'}
+  nrrddatatype = 'float';
+  
+ case {'block'}
+  error('Data type ''block'' is currently not supported.\n');
+  
+ otherwise
+  error('Unknown matlab datatype %s', matlab_metaType)
+end
+end
+
+% ========================================================================
+% writeData -->
+% fidIn is the open file we're overwriting
+% matrix - data that have to be written
+% datatype - type of data: int8, string, double...
+% encoding - raw, gzip, txt/ascii
+% ========================================================================
+function ok = writeData(fidIn, matrix, datatype, encoding)
+
+switch (encoding)
+ case {'raw'}
+  
+  ok = fwrite(fidIn, matrix(:), datatype);
+  
+ case {'gzip', 'gz'}
+     
+     % Store in a raw file before compressing
+     tmpBase = tempname(pwd);
+     fidTmpRaw = fopen(tmpBase, 'w');
+     assert(fidTmpRaw > 3, 'Could not open temporary file for GZIP compression');
+     try
+         fwrite(fidTmpRaw, matrix(:), datatype);
+         fclose(fidTmpRaw);
+     catch me
+         fclose(fidTmpRaw);
+         delete(tmpBase);
+         rethrow(me);
+     end
+     
+     
+     % Now we gzip our raw file
+     tmpFile = [tmpBase '.gz'];
+     try
+         gzip(tmpBase);     % this actually creates tmpFile
+     catch me 
+         delete(tmpBase);
+         rethrow(me);
+     end
+     delete(tmpBase);
+     
+     % Finally, we put this info into our nrrd file (fidIn)
+     
+     fidTmpRaw = fopen(tmpFile, 'r'); % should this be in a try catch as well?
+     assert(fidTmpRaw > 3, 'Could not open temporary file for writing to nrrd file during gzip compression.');
+     try
+         % tmp = fread(fidTmpRaw, Inf, [datatype '=>' datatype]); % precision argument : from datatype (source) to datatype, how about just byte by byte ?
+         tmp = fread(fidTmpRaw, Inf, 'uint8=>uint8');   % this seems to be more robust
+         fclose(fidTmpRaw) ;
+     catch me
+         fclose(fidTmpRaw);
+         delete(tmpFile);
+         rethrow(me);
+     end
+     
+     % ok = fwrite (fidIn, tmp, datatype); % why not byte by byte here?
+     ok = fwrite(fidIn, tmp, 'uint8');      % this seems to be more robust
+     delete (tmpFile);
+
+ case {'text', 'txt', 'ascii'}
+  
+  ok = fprintf(fidIn,'%u ', matrix(:)); % FIX: better with %g ?
+  %ok = fprintf(fidIn,matrix(:), class(matrix));
+  
+ otherwise
+  error('Unsupported encoding %s', encoding)
+end
+
+end
+
+% ========================================================================
+% Compare a Matlab-type size vector (no trailing ones, minimum length 2) to
+% the size vector specified in a nrrd file/header which may contain
+% trailing ones and may have only one dimension.
+% ========================================================================
+
+function sizesok = nrrd_size_check(matlab_size, nrrd_size)
+% Make row vectors
+matlab_size = matlab_size(:)';
+nrrd_size = nrrd_size(:)';
+
+% nrrd sizes may contain only one dimension
+if length(nrrd_size)==1
+    sizesok = prod(matlab_size)==nrrd_size;
+elseif length(nrrd_size)>=2
+    ind_max = find(nrrd_size~=1,1,'last');  % find last non-1 entry
+    % ones are ok in the first two dimensions:
+    if isempty(ind_max)
+        ind_max = 2;                
+    else
+        ind_max = max(2, ind_max);  % does not work if ind_max is empty
+    end
+    % Compare matlab size to nrrd size without trailing ones after from
+    % third entry on.
+    sizesok = isequal(matlab_size, nrrd_size(1:ind_max));
+else
+    error('Argument nrrd_size is an empty matrix');
+end
+
+end

matlab/nrrd_read_write_rensonnet/nrrd_getMatlabDataType.m

@@ -0,0 +1,48 @@
+% Get matlab type corresponding to each nrrd type
+% Date: November 2017
+% Author: Gaetan Rensonnet
+function matlabdatatype = nrrd_getMatlabDataType(metaType)
+
+% Determine the datatype
+switch (metaType)
+    case {'signed char', 'int8', 'int8_t'}
+        matlabdatatype = 'int8';
+        
+    case {'uchar', 'unsigned char', 'uint8', 'uint8_t'}
+        matlabdatatype = 'uint8';
+        
+    case {'short', 'short int', 'signed short', 'signed short int', ...
+            'int16', 'int16_t'}
+        matlabdatatype = 'int16';
+        
+    case {'ushort', 'unsigned short', 'unsigned short int', 'uint16', ...
+            'uint16_t'}
+        matlabdatatype = 'uint16';
+        
+    case {'int', 'signed int', 'int32', 'int32_t'}
+        matlabdatatype = 'int32';
+        
+    case {'uint', 'unsigned int', 'uint32', 'uint32_t'}
+        matlabdatatype = 'uint32';
+        
+    case {'longlong', 'long long', 'long long int', 'signed long long', ...
+            'signed long long int', 'int64', 'int64_t'}
+        matlabdatatype = 'int64';
+        
+    case {'ulonglong', 'unsigned long long', 'unsigned long long int', ...
+            'uint64', 'uint64_t'}
+        matlabdatatype = 'uint64';
+        
+    case {'float'}
+        matlabdatatype = 'single';
+        
+    case {'double'}
+        matlabdatatype = 'double';
+        
+    case {'block'}
+        error('Data type ''block'' is currently not supported.\n');
+        
+    otherwise
+        assert(false, 'Unknown datatype')
+end
+end

matlab/nrrd_read_write_rensonnet/nrrd_getSpaceDimensions.m

@@ -0,0 +1,49 @@
+% Get number of space dimensions ("space dimension" field) according to the
+% space descriptor, following the format definition at
+% http://teem.sourceforge.net/nrrd/format.html#spacedirections.
+%
+% Returns 3 if "space" field is one of
+% "right-anterior-superior" or "RAS"
+% "left-anterior-superior" or "LAS"
+% "left-posterior-superior" or "LPS"
+% "scanner-xyz"
+% "3D-right-handed"
+% "3D-left-handed"
+%
+% Returns 4 if "space" field is one of
+% "right-anterior-superior-time" or "RAST"
+% "left-anterior-superior-time" or "LAST"
+% "left-posterior-superior-time" or "LPST"
+% "scanner-xyz-time"
+% "3D-right-handed-time"
+% "3D-left-handed-time"
+%
+% Date: October 25, 2017
+% Author: Gaetan Rensonnet
+function sd = nrrd_getSpaceDimensions(spacedescriptor)
+
+if any(strcmpi(spacedescriptor,...
+        {'right-anterior-superior', 'RAS',...
+        'left-anterior-superior', 'LAS',...
+        'left-posterior-superior', 'LPS',...
+        'scanner-xyz',...
+        '3D-right-handed',...
+        '3D-left-handed'}...
+        ))
+    sd = 3;
+    
+elseif any(strcmpi(spacedescriptor,...
+        {'right-anterior-superior-time', 'RAST',...
+        'left-anterior-superior-time', 'LAST',...
+        'left-posterior-superior-time', 'LPST', ...
+        'scanner-xyz-time', ...
+        '3D-right-handed-time', ...
+        '3D-left-handed-time'}...
+        ))
+    sd = 4;
+    
+else
+    sd = -1;
+    % Unrecognized nrrd space descriptor (grace under fire)
+end
+end

matlab/nrrdread.m

@@ -0,0 +1,185 @@
+function [X, meta] = nrrdread(filename)
+%NRRDREAD  Import NRRD imagery and metadata.
+%   [X, META] = NRRDREAD(FILENAME) reads the image volume and associated
+%   metadata from the NRRD-format file specified by FILENAME.
+%
+%   Current limitations/caveats:
+%   * "Block" datatype is not supported.
+%   * Only tested with "gzip" and "raw" file encodings.
+%   * Very limited testing on actual files.
+%   * I only spent a couple minutes reading the NRRD spec.
+%
+%   See the format specification online:
+%   http://teem.sourceforge.net/nrrd/format.html
+
+% Copyright 2012 The MathWorks, Inc.
+
+
+% Open file.
+fid = fopen(filename, 'rb');
+assert(fid > 0, 'Could not open file.');
+cleaner = onCleanup(@() fclose(fid));
+
+% Magic line.
+theLine = fgetl(fid);
+assert(numel(theLine) >= 4, 'Bad signature in file.')
+assert(isequal(theLine(1:4), 'NRRD'), 'Bad signature in file.')
+
+% The general format of a NRRD file (with attached header) is:
+% 
+%     NRRD000X
+%     <field>: <desc>
+%     <field>: <desc>
+%     # <comment>
+%     ...
+%     <field>: <desc>
+%     <key>:=<value>
+%     <key>:=<value>
+%     <key>:=<value>
+%     # <comment>
+% 
+%     <data><data><data><data><data><data>...
+
+meta = struct([]);
+
+% Parse the file a line at a time.
+while (true)
+
+  theLine = fgetl(fid);
+  
+  if (isempty(theLine) || feof(fid))
+    % End of the header.
+    break;
+  end
+  
+  if (isequal(theLine(1), '#'))
+      % Comment line.
+      continue;
+  end
+  
+  % "fieldname:= value" or "fieldname: value" or "fieldname:value"
+  parsedLine = regexp(theLine, ':=?\s*', 'split','once');
+  
+  assert(numel(parsedLine) == 2, 'Parsing error')
+  
+  field = lower(parsedLine{1});
+  value = parsedLine{2};
+  
+  field(isspace(field)) = '';
+  meta(1).(field) = value;
+  
+end
+
+datatype = getDatatype(meta.type);
+
+% Get the size of the data.
+assert(isfield(meta, 'sizes') && ...
+       isfield(meta, 'dimension') && ...
+       isfield(meta, 'encoding') && ...
+       isfield(meta, 'endian'), ...
+       'Missing required metadata fields.')
+
+dims = sscanf(meta.sizes, '%d');
+ndims = sscanf(meta.dimension, '%d');
+assert(numel(dims) == ndims);
+
+data = readData(fid, meta, datatype);
+data = adjustEndian(data, meta);
+
+% Reshape and get into MATLAB's order.
+X = reshape(data, dims');
+X = permute(X, [2 1 3]);
+
+
+
+function datatype = getDatatype(metaType)
+
+% Determine the datatype
+switch (metaType)
+ case {'signed char', 'int8', 'int8_t'}
+  datatype = 'int8';
+  
+ case {'uchar', 'unsigned char', 'uint8', 'uint8_t'}
+  datatype = 'uint8';
+
+ case {'short', 'short int', 'signed short', 'signed short int', ...
+       'int16', 'int16_t'}
+  datatype = 'int16';
+  
+ case {'ushort', 'unsigned short', 'unsigned short int', 'uint16', ...
+       'uint16_t'}
+  datatype = 'uint16';
+  
+ case {'int', 'signed int', 'int32', 'int32_t'}
+  datatype = 'int32';
+  
+ case {'uint', 'unsigned int', 'uint32', 'uint32_t'}
+  datatype = 'uint32';
+  
+ case {'longlong', 'long long', 'long long int', 'signed long long', ...
+       'signed long long int', 'int64', 'int64_t'}
+  datatype = 'int64';
+  
+ case {'ulonglong', 'unsigned long long', 'unsigned long long int', ...
+       'uint64', 'uint64_t'}
+  datatype = 'uint64';
+  
+ case {'float'}
+  datatype = 'single';
+  
+ case {'double'}
+  datatype = 'double';
+  
+ otherwise
+  assert(false, 'Unknown datatype')
+end
+
+
+
+function data = readData(fidIn, meta, datatype)
+
+switch (meta.encoding)
+ case {'raw'}
+  
+  data = fread(fidIn, inf, [datatype '=>' datatype]);
+  
+ case {'gzip', 'gz'}
+
+  tmpBase = tempname();
+  tmpFile = [tmpBase '.gz'];
+  fidTmp = fopen(tmpFile, 'wb');
+  assert(fidTmp > 3, 'Could not open temporary file for GZIP decompression')
+  
+  tmp = fread(fidIn, inf, 'uint8=>uint8');
+  fwrite(fidTmp, tmp, 'uint8');
+  fclose(fidTmp);
+  
+  gunzip(tmpFile)
+  
+  fidTmp = fopen(tmpBase, 'rb');
+  cleaner = onCleanup(@() fclose(fidTmp));
+  
+  meta.encoding = 'raw';
+  data = readData(fidTmp, meta, datatype);
+  
+ case {'txt', 'text', 'ascii'}
+  
+  data = fscanf(fidIn, '%f');
+  data = cast(data, datatype);
+  
+ otherwise
+  assert(false, 'Unsupported encoding')
+end
+
+
+
+function data = adjustEndian(data, meta)
+
+[~,~,endian] = computer();
+
+needToSwap = (isequal(endian, 'B') && isequal(lower(meta.endian), 'little')) || ...
+             (isequal(endian, 'L') && isequal(lower(meta.endian), 'big'));
+         
+if (needToSwap)
+    data = swapbytes(data);
+end

matlab/plotFunctions.m

@@ -0,0 +1,24 @@
+function plotFunctions(path,code,globalPar,cPars,cax,centers)
+    n=size(centers,1);
+    cax1=linspace(cax(1),cax(size(cax,1)),1000)';
+    fv=evaluateCenters(globalPar,cPars,cax1);
+    
+    
+    for i=1:n
+        m=max(centers(i,:));
+        m1=max(fv(i,:));
+        if m1>m
+            m=m1;
+        end
+        f=figure('visible','off');
+        plot(cax1,fv(i,:))
+        hold on
+        scatter(cax,centers(i,:))
+        text(500,0.8*m,sprintf('k1=%.2f ml/g*min',60*cPars(i,1)));
+        text(500,0.7*m,sprintf('BVF=%.2f',cPars(i,2)));
+        %text(500,0.6*m,sprintf('k2=%f',cPars(i,2)));
+        hold off
+        saveas(f,fullfile(path,sprintf('%s_centers%d.png',code,i)))
+    end
+    
+    

matlab/readFitParameters.m

@@ -0,0 +1,21 @@
+function p=readFitParameters(path,patientID,nclass)
+    nr=20;
+    pf=zeros(nr,4);
+    for realizationId=1:nr
+        fname=sprintf('%s_%d_%d_fitParFinal.txt',patientID,nclass,realizationId);
+        of=fullfile(path,patientID,fname);
+        fitPar=readmatrix(of,'Delimiter','\t');
+        ft=[fitPar(3) 1/fitPar(4)];%rise and fall time; due to symmetry, they get confused, and max is always fall time and min is rise
+        pf(realizationId,2)=max(ft);
+        pf(realizationId,3)=1/min(ft);
+        pf(realizationId,1)=fitPar(2);%constant
+        pf(realizationId,4)=fitPar(5);%delay
+        %size(fitPar)
+        %fprintf('%.2f %.2f %.2f %.2f\n',fitPar(1),fitPar(2),pf(realizationId,2),pf(realizationId,3));
+        
+    
+    end
+    p=pf;
+    %fprintf('%.2f %.2f %.2f\n',median(pf(1,:)), median(pf(2,:)),median(pf(3,:)));
+  
+end

matlab/rfp.m

@@ -0,0 +1,8 @@
+function [globalPar, cPars]=rfp(path,patientID,nclass,realizationId)
+    code=sprintf('%s_%d_%d_IVF',patientID,nclass,realizationId);
+    ofile=fullfile(path,patientID,sprintf('%s_fitParFinal.txt',code));
+    mat=readmatrix(ofile,'Delimiter','\t');
+    globalPar=[mat(2:5) 0 0 0];
+    qfit=mat(6:end);
+    cPars=reshape(qfit,[],nclass)';
+end

matlab/saveCenters.m

@@ -0,0 +1,41 @@
+function saveCenters(path, patientID, cax, cm, data, nclass, nRealizations)
+    %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
+    %nRealizations is number of realizations
+    
+    c1=size(cax,1);
+    A=reshape(data,[],c1);
+    
+    if nclass==10
+        fuzzyPar=1.5;
+    end
+    if nclass==20
+        fuzzyPar=1.45;
+    end
+    
+    if nclass==30
+        fuzzyPar=1.4;
+    end
+    
+    
+    
+    options=[fuzzyPar 100 1 true];
+    
+    dt=size(data);
+    for j=1:nRealizations
+        [centers,U]=fcm(A,nclass,options);
+        for i=1:nclass
+            code=sprintf('%s_%d_%d_center%d',patientID,nclass,j,i);
+            ofile=fullfile(path,patientID,sprintf('%s_center.txt',code));
+            writematrix(centers(i,:),ofile,'Delimiter',',');
+            u=reshape(U(i,:),dt(1:3));
+            headerInfoOut.data=single(u);
+            headerInfoOut.content=sprintf('%s_centerWeight',code);
+            nhdr_nrrd_write(fullfile(path,patientID, sprintf('%s_centerWeigth.nrrd',code)),headerInfoOut,1);
+        end
+    end
+
+   

matlab/writeData.m

@@ -0,0 +1,38 @@
+%path already set by loadData
+headerInfoOut=headerInfo;
+
+k1=squeeze(xa(:,:,:,1));
+vbf=squeeze(xa(:,:,:,2));
+
+if size(xa,4)>2
+    k2=squeeze(xa(:,:,:,3));
+end
+
+headerInfoOut.data=single(k1);
+headerInfoOut.content=[patientID '_k1'];
+nhdr_nrrd_write(fullfile(path,[patientID '_k1.nrrd']),headerInfoOut,1)
+    
+if size(xa,4)>2
+    headerInfoOut.data=single(k2);
+    headerInfoOut.content=[patientID '_k2'];
+    nhdr_nrrd_write(fullfile(path,[patientID '_k2.nrrd']),headerInfoOut,1)
+end
+
+headerInfoOut.data=single(vbf);
+headerInfoOut.content=[patientID '_vbf'];
+nhdr_nrrd_write(fullfile(path,[patientID '_vbf.nrrd']),headerInfoOut,1)
+
+dt=size(data);
+for i=1:nclass
+    u=reshape(U(i,:),dt(1:3));
+    headerInfoOut.data=single(u);
+    headerInfoOut.content=sprintf('%s_center%i',patientID,i);
+    nhdr_nrrd_write(fullfile(path,sprintf('%s_center%d.nrrd',patientID,i)),headerInfoOut,1);
+    
+end
+
+%for i=1:c1
+%    headerInfoOut.data=single(data1(:,:,:,i))*cm(i);
+%    s=sprintf('%sVolume%d_model.nrrd',patientID,i-1);
+%    nhdr_nrrd_write(fullfile(path,s),headerInfoOut,1);
+%end

pythonScripts/config.py

@@ -0,0 +1,28 @@
+import os
+
+def getPatientId(row,xconfig):
+   return row[xconfig['ParticipantField']]
+
+def getVisitId(row,xconfig):
+   return row['visitName']
+
+def getCode(row,xconfig):
+   return '{}_{}'.format(getPatientId(row,xconfig),getVisitId(row,xconfig))
+
+def getPathList(row,xconfig):
+   return [xconfig['baseDir'],getPatientId(row,xconfig),getVisitId(row,xconfig)]
+
+def getOutputDir(row,xconfig):
+   return '/'.join(getPathList(row,xconfig))
+
+def getLocalDir(row,xconfig):
+   return os.path.join(xconfig['tempDir'],getCode(row,xconfig))
+
+def getNodeName(row,xconfig,mode,i=0):
+   if mode=='CT':
+      return getCode(row,xconfig)+'_CT'
+   if mode=='NM':
+      return '{}_Volume{}'.format(getCode(row,xconfig),i)
+   return getCode(row,xconfig)+'_Dummy'
+
+

pythonScripts/generateCenters.ipynb

@@ -101,7 +101,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.10"
   }
  },
  "nbformat": 4,

scripts/analyzeAll.sh

@@ -0,0 +1,10 @@
+#!/bin/bash
+
+LIST="2SMobr 2SBMIR 3ZFOBR 3ZFMIR 5MIMIR 5MIOBR"
+LIST="$LIST 7TMMIR 7TMOBR 10KFMIR 10KFOBR 11ZMMIR 11ZMOBR"
+
+for id in $LIST;
+do
+	./doAnalysis.sh $id;
+
+done;

scripts/doAnalysis.sh

@@ -0,0 +1,23 @@
+#!/bin/bash
+
+if [ $# -lt 1 ] ; then
+	echo "Usage $0 patientID \(2SMobr\) \(id\)"
+	exit 0;
+fi;
+
+SDIR=`dirname $0`;
+PATIENTID=$1;
+NCLASS=10
+if [ $# -gt 1 ] ; then
+	NCLASS=$2;
+fi;
+
+ID=0;
+if [ $# -gt 2 ] ; then
+	ID=$3;
+fi;
+
+. $SDIR/env.sh
+
+$MATLAB -sd $MATLABSCRIPTDIR -batch "patientID='${PATIENTID}'; nclass='${NCLASS}'; realizationId='${ID}'; analyze" >& ~/logs/dynamicSPECT.log;
+

scripts/doAnalysisIVF.sh

@@ -0,0 +1,26 @@
+#!/bin/bash
+
+if [ $# -lt 1 ] ; then
+	echo "Usage $0 patientID \(2SMobr\) \(nclass\) \(realizationId\)"
+	exit 0;
+fi;
+
+SDIR=`dirname $0`;
+PATIENTID=$1;
+NCLASS=10
+if [ $# -gt 1 ] ; then
+	NCLASS=$2;
+fi;
+
+ID=0;
+if [ $# -gt 2 ] ; then
+	ID=$3;
+fi;
+
+. $SDIR/env.sh
+
+#$SDIR/loadLabkeyData.sh $PATIENTID;
+$MATLAB -sd $MATLABSCRIPTDIR -batch "patientID='${PATIENTID}'; nclass='${NCLASS}'; realizationId='${ID}'; analyzeIVF" >& ~/logs/dynamicSPECT.log;
+#$SDIR/writeLabkeyData.sh $PATIENTID;
+#$SDIR/removeData.sh $PATIENTID;
+

scripts/doPixelAnalysis.sh

@@ -0,0 +1,22 @@
+#!/bin/bash
+
+if [ $# -lt 1 ] ; then
+	echo "Usage $0 patientID \(2SMobr\) \(sigma2\)"
+	exit 0;
+fi;
+
+SDIR=`dirname $0`;
+PATIENTID=$1;
+SIGMA2=1.0
+if [ $# -gt 1 ] ; then
+	SIGMA2=$2;
+fi;
+
+
+. $SDIR/env.sh
+
+#$SDIR/loadLabkeyData.sh $PATIENTID;
+$MATLAB -sd $MATLABSCRIPTDIR -batch "patientID='${PATIENTID}'; sigma2='${SIGMA2}'; analyzePixel" >& ~/logs/dynamicSPECT.log;
+#$SDIR/writeLabkeyData.sh $PATIENTID;
+#$SDIR/removeData.sh $PATIENTID;
+

scripts/doPixelIVFAnalysis.sh

@@ -0,0 +1,22 @@
+#!/bin/bash
+
+if [ $# -lt 1 ] ; then
+	echo "Usage $0 patientID \(2SMobr\) \(sigma2\)"
+	exit 0;
+fi;
+
+SDIR=`dirname $0`;
+PATIENTID=$1;
+SIGMA2=1.0
+if [ $# -gt 1 ] ; then
+	SIGMA2=$2;
+fi;
+
+
+. $SDIR/env.sh
+
+#$SDIR/loadLabkeyData.sh $PATIENTID;
+$MATLAB -sd $MATLABSCRIPTDIR -batch "patientID='${PATIENTID}'; sigma2='${SIGMA2}'; analyzePixelIVF" >& ~/logs/dynamicSPECT.log;
+#$SDIR/writeLabkeyData.sh $PATIENTID;
+#$SDIR/removeData.sh $PATIENTID;
+

scripts/env.sh

@@ -0,0 +1,6 @@
+export DSPECTPROJECT='dinamic_spect/Patients'
+export WEBDAV=$HOME/scripts/webdav
+export DSPECTTMP=$HOME/temp/dynamicSPECT
+
+export MATLABSCRIPTDIR=$HOME/software/src/dynamicSPECT/matlab
+export MATLAB=$HOME/software/install/matlab/bin/matlab;

scripts/generateCenters.sh

@@ -0,0 +1,26 @@
+#!/bin/bash
+
+if [ $# -lt 1 ] ; then
+	echo "Usage $0 patientID \(2SMobr\) \(nclass\) \(nRealizations\)"
+	exit 0;
+fi;
+
+SDIR=`dirname $0`;
+PATIENTID=$1;
+NCLASS=10
+if [ $# -gt 1 ] ; then
+	NCLASS=$2;
+fi;
+
+NR=20;
+if [ $# -gt 2 ] ; then
+	NR=$3;
+fi;
+
+. $SDIR/env.sh
+
+#$SDIR/loadLabkeyData.sh $PATIENTID;
+$MATLAB -sd $MATLABSCRIPTDIR -batch "patientID='${PATIENTID}'; nclass='${NCLASS}'; nRealizations='${NR}'; generateCenters" >& ~/logs/dynamicSPECT.log;
+#$SDIR/writeLabkeyData.sh $PATIENTID;
+#$SDIR/removeData.sh $PATIENTID;
+

scripts/loadLabkeyData.sh

@@ -0,0 +1,36 @@
+#!/bin/bash
+
+. `dirname $0`/env.sh 
+
+
+PATIENTID=$1
+
+$WEBDAV/set_project.sh $DSPECTPROJECT
+OUTDIR=$DSPECTTMP/$PATIENTID;
+
+if [ ! -d $OUTDIR ] ; then
+	mkdir $OUTDIR;
+fi;
+
+
+for ((i=0;i<20;i++)) ; do
+	F=${PATIENTID}Volume$i.nrrd;
+	$WEBDAV/getFile.sh $PATIENTID/$F $OUTDIR/$F;
+done;
+
+F=${PATIENTID}_Ventricle.mcsv;
+$WEBDAV/getFile.sh $PATIENTID/$F $OUTDIR/$F;
+
+F=${PATIENTID}_Dummy.mcsv;
+$WEBDAV/getFile.sh $PATIENTID/$F $OUTDIR/$F;
+
+FILES="k1.nrrd vbf.nrrd";
+for f in $FILES ; 
+do
+	F=${PATIENTID}_$f;
+	$WEBDAV/getFile.sh $PATIENTID/$F $OUTDIR/$F;
+done;
+
+
+
+

scripts/removeData.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+
+. `dirname $0`/env.sh
+
+PATIENTID=$1
+
+
+rm -rf $DSPECTTMP/$1;
+
+
+
+
+
+

scripts/writeLabkeyData.sh

@@ -0,0 +1,39 @@
+#!/bin/bash
+
+
+. `dirname $0`/env.sh
+
+PATIENTID=$1
+
+
+$WEBDAV/set_project.sh $DSPECTPROJECT
+
+
+FILES="k1.nrrd vbf.nrrd";
+
+for F in $FILES ; do
+	f=${PATIENTID}_$F;
+	$WEBDAV/cp.sh $DSPECTTMP/$PATIENTID/$f $PATIENTID/$f;
+done;
+
+
+$WEBDAV/cp.sh $DSPECTTMP/$PATIENTID/fitPar.txt $PATIENTID/fitPar.txt;
+
+for ((i=1;i<11;i++)) ;
+do
+	f=centers${i}.png
+	$WEBDAV/cp.sh $DSPECTTMP/$PATIENTID/$f $PATIENTID/$f;
+
+	f=${PATIENTID}_center${i}.nrrd
+	$WEBDAV/cp.sh $DSPECTTMP/$PATIENTID/$f $PATIENTID/$f;
+
+
+done;
+
+exit 0;
+
+
+
+
+
+

slicerScripts/convertToNRRD.py

@@ -9,7 +9,7 @@ import shutil
 def main(configFile):
    print('Running with {}'.format(configFile))
    with open(configFile) as f:
-      config=json.load(f)
+      xconfig=json.load(f)
 
    fsetup=os.path.join(os.path.expanduser('~'),'.labkey','setup.json')
    with open(fsetup) as f:
@@ -29,9 +29,11 @@ def main(configFile):
    import orthancDatabaseBrowser
    import orthancFileBrowser
 
+   sys.path.append('../pythonScripts')
+   import config
 
    net=labkeyInterface.labkeyInterface()
-   fnet=os.path.join(os.path.expanduser('~'),'.labkey',config['network'])
+   fnet=os.path.join(os.path.expanduser('~'),'.labkey',xconfig['network'])
    net.init(fnet)
    net.getCSRF()
    fb=labkeyFileBrowser.labkeyFileBrowser(net)
@@ -44,7 +46,7 @@ def main(configFile):
    odb=orthancDatabaseBrowser.orthancDB(onet)
 
    qFilter=[]  
-   ds=db.selectRows(config['project'],config['schemaName'],config['queryName'],qFilter)
+   ds=db.selectRows(xconfig['project'],xconfig['schemaName'],xconfig['queryName'],qFilter)
 
 
    try:
@@ -53,14 +55,13 @@ def main(configFile):
       print('No rows returned')
       return
    for r in rows:
-      print("Loading {}/{}".format(patientId(r,config),visitId(r,config)))
-      nNames=nodeNames(r,config)
-      rPath=fb.formatPathURL(config['project'],outputDir(r,config))
-      rPath+='/'+nNames['CT']+'.nrrd'
+      print("Loading {}/{}".format(config.getPatientId(r,xconfig),config.getVisitId(r,xconfig)))
+      rPath=fb.formatPathURL(xconfig['project'],config.getOutputDir(r,xconfig))
+      rPath+='/'+config.getNodeName(r,xconfig,'CT')+'.nrrd'
       entryDone=fb.entryExists(rPath)
       if entryDone:
          try:
-            if not config['recalculate']:
+            if not xconfig['recalculate']:
                print('Entry done')
                continue
          except KeyError:
@@ -70,54 +71,29 @@ def main(configFile):
          print('Forced recalculation')
 
             
-      print('{} available:{}'.format(nNames['CT'],fb.entryExists(rPath)))
+      print('{} available:{}'.format(config.getNodeName(r,xconfig,'CT'),fb.entryExists(rPath)))
       #loadPatient into slicer
-      patient=loadPatient(ofb,r,config)
+      patient=loadPatient(ofb,r,xconfig)
       #convert to nodes
-      addCT(r,patient,config)
-      addFrames(r,patient,config)
-      addDummyInputFunction(r,patient,config)
+      addCT(r,patient,xconfig)
+      addFrames(r,patient,xconfig)
+      addDummyInputFunction(r,patient,xconfig)
       
       nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLScalarVolumeNode')
       print('Nodes')
       for node in nodes:
          print('\t{}'.format(node.GetName()))
-         storeNode(fb,r,config,node)
+         storeNode(fb,r,xconfig,node)
 
       nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLDoubleArrayNode')
       print('Nodes (double array)')
       for node in nodes:
          print('\t{}'.format(node.GetName()))
-         storeNode(fb,r,config,node)
-      clearNodes(r,config)
+         storeNode(fb,r,xconfig,node)
+      clearNodes(r,xconfig)
 
       
-def patientId(row,config):
-   return row[config['ParticipantField']]
-
-def visitId(row,config):
-   return row['visitName']
-
-def code(row,config):
-   return '{}_{}'.format(patientId(row,config),visitId(row,config))
-
-def pathList(row,config):
-   return [config['baseDir'],patientId(row,config),visitId(row,config)]
-
-def outputDir(row,config):
-   return '/'.join(pathList(row,config))
-
-def localDir(row,config):
-   return os.path.join(config['tempDir'],code(row,config))
-
-def nodeNames(row,config):
-   names={}
-   x=code(row,config)
-   names['NM']=['{}_Volume{}'.format(x,i) for i in range(20)]
-   names['CT']='{}_CT'.format(x)
-   return names
-
-def clearNodes(row,config):
+def clearNodes(row,xconfig):
    nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLScalarVolumeNode')
    nodes1=slicer.mrmlScene.GetNodesByClass('vtkMRMLDoubleArrayNode')
    for n in nodes1:
@@ -125,25 +101,30 @@ def clearNodes(row,config):
 
    res=[slicer.mrmlScene.RemoveNode(f) for f in nodes]
 
-def loadPatient(ofb,r,config):
+def loadPatient(ofb,r,xconfig):
    sys.path.append('../pythonScripts')
    import parseDicom
    import vtkInterface
 
    pd=parseDicom.parseDicom()
 
-   masterPath=downloadAndUnzip(ofb,r,"nmMaster",config)
+   masterPath=downloadAndUnzip(ofb,r,"nmMaster",xconfig)
    pd.readMasterDirectory(masterPath,False)
    print('Time [{} .. {}]'.format(pd.frame_start,pd.frame_stop))
-   nmPath=downloadAndUnzip(ofb,r,"nmCorrData",config)
+   clearUnzipDir(r,xconfig)
+
+   nmPath=downloadAndUnzip(ofb,r,"nmCorrData",xconfig)
    frame_data, frame_time, frame_duration, frame_origin, \
        frame_pixel_size, frame_orientation=\
        pd.readNMDirectory(nmPath,False)
    print('Frame time {}'.format(frame_time))
-   ctPath=downloadAndUnzip(ofb,r,"ct",config)
+   clearUnzipDir(r,xconfig)
+
+   ctPath=downloadAndUnzip(ofb,r,"ct",xconfig)
    ct_data,ct_origin,ct_pixel_size, \
       ct_orientation=pd.readCTDirectory(ctPath,False)
    print('CT pixel {}'.format(ct_pixel_size))
+   clearUnzipDir(r,xconfig)
 
    ct_orientation=vtkInterface.completeOrientation(ct_orientation)
    frame_orientation=vtkInterface.completeOrientation(frame_orientation)
@@ -158,16 +139,25 @@ def loadPatient(ofb,r,config):
    
    print('Done')
 
-def downloadAndUnzip(ofb,r,xcode,config):
-   orthancId=r[xcode+'OrthancId']
-   fileCode='{}_{}'.format(code(r,config),xcode)
-   zipFile=os.path.join(config['tempDir'],fileCode+'.zip')
-   ofb.getZip('series',orthancId,zipFile)
-   zipDir=localDir(r,config)
+def clearUnzipDir(r,xconfig):
+   import config
+
+   zipDir=config.getLocalDir(r,xconfig)
    try:
       os.mkdir(zipDir)
    except FileExistsError:
       shutil.rmtree(zipDir)
+   return zipDir
+
+def downloadAndUnzip(ofb,r,code,xconfig):
+   import config 
+
+   orthancId=r[code+'OrthancId']
+   fileCode='{}_{}'.format(config.getCode(r,xconfig),code)
+   zipFile=os.path.join(xconfig['tempDir'],fileCode+'.zip')
+   ofb.getZip('series',orthancId,zipFile)
+   zipDir=clearUnzipDir(r,xconfig)
+
    try:
       outTxt=subprocess.check_output(["unzip","-d",zipDir,"-xj",zipFile])
    except subprocess.CalledProcessError:
@@ -176,25 +166,26 @@ def downloadAndUnzip(ofb,r,xcode,config):
 
    return zipDir
 
-def addCT(r,patient,config):
+def addCT(r,patient,xconfig):
    sys.path.append('../pythonScripts')
+   import config
    import vtkInterface
    ct=patient['CT']
    vtkData=vtkInterface.numpyToVTK(ct['data'],ct['data'].shape)
-   nodeName=code(r,config)+'_CT'
+   nodeName=config.getNodeName(r,xconfig,'CT')
    addNode(nodeName,vtkData,ct['origin'],ct['pixel_size'],ct['orientation'],0)
 
-def addFrames(r,patient,config):
-   sys.path.append('../pythonScripts')
+
+def addFrames(r,patient,xconfig):
    import vtkInterface
+   import config
    #convert data from numpy.array to vtkImageData
    #use time point it
-   x=code(r,config)
    nm=patient['NM']
    print("NFrames: {}".format(nm['data'].shape[3]))
    for it in range(0,nm['data'].shape[3]):
       frame_data=nm['data'][:,:,:,it];
-      nodeName=x+'_Volume'+str(it)
+      nodeName=config.getNodeName(r,xconfig,'NM',it)
       vtkData=vtkInterface.numpyToVTK(frame_data,frame_data.shape)
       addNode(nodeName,vtkData,nm['origin'],nm['pixel_size'],nm['orientation'],1)
 
@@ -226,12 +217,14 @@ def addNode(nodeName,v,origin,pixel_size,orientation,dataType):
    newNode.SetAndObserveImageData(v)
    slicer.mrmlScene.AddNode(newNode)
 
-def addDummyInputFunction(r,patient,config):
+def addDummyInputFunction(r,patient,xconfig):
 
+   import config
+   
    nm=patient['NM']
    n=nm['data'].shape[3]
 
-   dnsNodeName=code(r,config)+'_Dummy'
+   dnsNodeName=config.getNodeName(r,xconfig,'Dummy')
    dns = slicer.mrmlScene.GetNodesByClassByName('vtkMRMLDoubleArrayNode',dnsNodeName)
    if dns.GetNumberOfItems() == 0:
       dn = slicer.mrmlScene.AddNode(slicer.vtkMRMLDoubleArrayNode())
@@ -251,7 +244,8 @@ def addDummyInputFunction(r,patient,config):
       dn.SetValue(i, 2, 0)
       print('{} ({},{})'.format(i,fx,fy))
 
-def storeNode(fb,row,config,node):
+def storeNode(fb,row,xconfig,node):
+   import config
 
    suffix=".nrrd"
    if node.__class__.__name__=="vtkMRMLDoubleArrayNode":
@@ -262,11 +256,11 @@ def storeNode(fb,row,config,node):
 
    fileName=node.GetName()+suffix
 
-   localPath=os.path.join(localDir(row,config),fileName)
+   localPath=os.path.join(config.getLocalDir(row,xconfig),fileName)
 
    slicer.util.saveNode(node,localPath)
    print("Stored to: {}".format(localPath))
-   labkeyPath=fb.buildPathURL(config['project'],pathList(row,config))
+   labkeyPath=fb.buildPathURL(xconfig['project'],config.getPathList(row,xconfig))
    print ("Remote: {}".format(labkeyPath))
    #checks if exists
    remoteFile=labkeyPath+'/'+fileName

template/cardiacSPECT.json

@@ -14,5 +14,6 @@
 	"tempDir":"/home/studen/temp/dynamicSPECT",
 	"PatientId":"3ZFMIR",
    "baseDir":"processedImages",
-   "ParticipantField":"PatientId"
+   "ParticipantField":"PatientId",
+   "recalculate":1
 }