% Toolbox for batch processing ASL perfusion based fMRI data.
% All rights reserved.
% Ze Wang @ TRC, CFN, Upenn 2004
%
% Batch realigning images.
% Get the global subject information

disp('Calculating SEN map for all subjects, just wait....');
% load spm defaults
global defaults;
spm_defaults;

 
clear P;
% dirnames, 
% get the subdirectories in the main directory
for sb =1:length(PAR.subjects) % for each subject 
    
    num=0;
    clear P;
    P=[];
    
    for c=1:PAR.ncond 
           if isempty(PAR.condirs{sb,c}), continue; end
           ana_dir = fullfile(PAR.root, PAR.subjects{sb}, PAR.rbimdir);
            if ~(exist(ana_dir))
              mdcom=['!mkdir ' ana_dir];
              eval(mdcom);
              %mkdir(ana_dir_pref,MYV.ana_dir);
            end
            %dirn = fullfile(PAR.root,PAR.fMRI,PAR.subjects{s},PAR.condirs{s,c});
            % get files in this directory
            P=spm_select('ExtFPList',PAR.condirs{sb, c},['^mask_bet.nii$'],1);
            v=spm_vol(P(1,:));
            dat=spm_read_vols(v);
                  
            mask=find(dat(:)>0.26*max(dat(:)));

            Ptmp=spm_select('ExtFPList',PAR.condirs{sb, c},['^sflt_r.*\.nii$'],1:400);

            v=spm_vol(Ptmp);
            dat=spm_read_vols(v);
            [dx, dy, dz, tlen]=size(dat);
            dat=reshape(dat, dx*dy*dz, tlen);
            dat=dat(mask,:);
%             gnoi=gauref(1,1:tlen);
            cd(PAR.condirs{sb,c})
            sum=0;
            std_d=std(dat,0,2);
            for scale=0:2
                
                if scale>0
                    mdat=mean(dat,2);
                    ndat=dat-repmat(mdat,1,tlen);
                    % data has been filtered using a bandwidth of 0.01
                    % to 0.1 Hz
                    [lb,la]=fltbutter(1,0.4*0.5^scale);    % if TR=2, omega=1~1/2TR =0.25 (fMax)
                    clear filter
                    ndat=filter(lb,la,ndat,[],2);
                else
                    ndat=dat;
                end
                ndat=ndat';   % need to be a column vector
                
                senmap=zeros(dx,dy,dz);
                tlens=size(ndat,1);
                
                inmaskdatfile='d_comp.dat';
                fprintf('\t\t\t saving inmask data for calc SEN.\n');
                fid=fopen(inmaskdatfile,'wb');
                fwrite(fid, length(mask), 'int');   % number of timecourses
                fwrite(fid, tlens, 'int');

                fwrite(fid, std_d, 'float');
                fwrite(fid, ndat, 'float');
                fclose(fid);
                
                str=['!/big1/zwang/workshop/sampen/SampEn -d ' num2str(PAR.dim) ' -r ' num2str(PAR.r) ' -i ' inmaskdatfile  ' -o SEN.dat '];
                eval(str);
                fid=fopen('SEN.dat','rb');
                slen=fread(fid, 1, 'int');
                dim=fread(fid, 1, 'int');
                ratio=fread(fid, 1, 'float');
                sSEN=fread(fid, slen, 'float');
                aSEN=fread(fid, slen, 'float');
                fclose(fid);

                senmap(mask) = sSEN*1e3;  
                senmap(isinf(senmap))=4e3;
                vo=v(1);
%                     vo.dt=[16 0];
                map1=sprintf('SEN_%s',PAR.subjects{sb});
                vo.fname=fullfile(ana_dir,sprintf('%s_SEN_S%d_D%d_r%f.nii',PAR.subjects{sb}, scale, PAR.dim, PAR.r)); %only 1 condition, so no condition name was marked

                vo=spm_write_vol(vo,senmap);

               
                sum=sum+senmap/mean(senmap(mask));
                
            end
            
            vo.fname=fullfile(ana_dir,sprintf('%s_SEN_Composite_D%d_r%f.nii',PAR.subjects{sb}, PAR.dim, PAR.r)); %only 1 condition, so no condition name was marked
%             vo=spm_write_vol(vo,sum);
  end
    
end

clear sb