function dt_BOLD_calc(fdir,con,outdir,outname)
%mjf 2016 5 3
%uploaded to nitrc 2017_09_28
%the intention of this is that iit calculates 
%the sigmoidal delay function of Henson
%dx.doi.org/10.1006/nimg.2001.0940
% for the SPM.mat and contrast specified in the input

% for this version, 'con' is the name of the condition specified in the
% design matrix. ie 'congr_right' 

%mjf 2016 11 11 changed spm_imcalc_ui to spm_imcalc (for spm 12)
%mjf 2016 11 15 changed so that the outputfilename '^' gets replaced by 'p'
%which is to help with parametric modulators, which are called things like
%xRT^1 and this filename is not particularly liked by unix type things.
%
% note the condition names get appended to outname, so call like:
% dt_BOLD_calc('d:/studies/catfield/CAT001/modRT/','easy_left',
%                 'd:/studies/catfield/modRT_dtsig','CAT001')

a=load(fullfile(fdir,'SPM.mat'));
swd=a.SPM.swd;



betas={};
form='(';
forma='(';
nn=0;
pl=' ';
%get the file names of the contrasts and the first derivative
%assumes that the HRF_dt is immediately after the HRF
for i=1:length(a.SPM.Vbeta)
        if(strfind(a.SPM.Vbeta(i).descrip,[con '*bf(1)']))
            if(strfind(a.SPM.Vbeta(i+1).descrip,[con '*bf(2)'])) 
                nn=nn+1;    
                if nn > 1
                    pl=' + ';
                end
                if max(a.SPM.xX.X(:,i))==0 && min(a.SPM.xX.X(:,i))==0 
                    disp('!!!nothing in condition!')
                end
                %filenames
                betas{nn*2-1}=fullfile(swd,a.SPM.Vbeta(i).fname);
                betas{nn*2}=fullfile(swd,a.SPM.Vbeta(i+1).fname);
                
                %sigmoidal calculation for each pair of images
                %as spm imcalc expression
                %=1/(1+exp(D b2/b1))
                form=[form pl sprintf('1./(1+exp(3.1.*i%d./i%d))',nn*2,nn*2-1)];
                %also to give average of main images
                forma=[forma pl sprintf('i%d',nn*2-1)];
            end
        end
end
%boldav = [boldav sprintf(')/%d)',length(bold))]
% should be something like
% ((i1 + i3 + i5 + i7 + i9  )/5)
%dtav = [dtav sprintf(')/%d)',length(bold))]

% average the above eq, and add on the constants to make it have the right
% values
% ie 2*1.78 * (average sigmoid) -1.78
spm_calc=[form sprintf(')./%d*3.56 -1.78',nn)];
spm_calcA=[forma sprintf(')./%d',nn)];



%
%spm_imcalc(betas,'output.nii', spm_calc1, {[],[],[16],[]});
conout=strrep(con,'^','p');
outp=fullfile(outdir,[outname '_' conout '_dtsig.nii']);
spm_imcalc(betas,outp, spm_calc, {[],[],[16],[]});
outp=fullfile(outdir,[outname '_' conout '_i.nii']);
spm_imcalc(betas,outp, spm_calcA, {[],[],[16],[]});
