% Toolbox for batch processing ASL perfusion based fMRI data.
% All rights reserved.
% Ze Wang @ TRC, CFN, Upenn 2004
%
%     please cite the following paper:
%             Wang Z, Li Y, Childress AR, Detre JA (2014) Brain Entropy Mapping Using fMRI. PLoS ONE 9(3): e89948. doi:10.1371/journal.pone.0089948	
%
%     several related abstracts can be found here:
%            1. Wang Z (2012) Characterizing Resting Brain Information using Voxel-based Brain Information Mapping (BIM). 2012 Annual Meeting of the Organization for Human Brain Mapping. Beijing, China.	
%            2. Wang Z (2012) Stable and Self-Organized Entropy in the Resting Brain. The Third Biennial Conference on Resting State Brain Connectivity. Magdeburg, Germany. pp. 208.
%            3. Ze Wang, Marcus Raichle, Anna Rose Childress, and John A Detre (2013) Mapping brain entropy using resting state fMRI. 2013 Annual Meeting of International Society of Magnetic Resonance in Medicine. Salt Lake City, USA. pp. 4861.
%            4. Ze Wang, Y Li, Z Singer, R Ehrman, A. V Hole, C. P O'Brien, Anna Rose Childress (2013) Human brain entropy mapping using thousands of subjects and its application in a drug addiction study. 2013 Annual Meeting of Society for Neuroscience. San Diego. pp. 7491.
%  The main function included in this toolkit is a C++ based compiled program for calculating sample entropy (SampEn) and approximate entropy (ApEn). The algorithm for both SampEn
%  and ApEn has been optimized by removing the redundant embedded vector matchings. The computation time for a 64x64x40 fMRI dataset is up to 10 secs in a regular computer, as compared to hours using a Matlab code
%            Users who are interested in getting the C++ source code can contact the author. Otherwise, we will only release the binary code. 
% Mac users should replace the file "SampEn" with the one in the folder "Mac" 
% Note: the multi-scale sample entropy (Costa M., Goldberger A.L., Peng C.-K. Multiscale entropy analysis of biological signals. Phys Rev E 2005;71:021906) is not implemented in the C++ code
%           but rather in this Matlab wrapper code.  scale 0 means no scale.  Increasing one scale is similar to taking the lower half of the current frequency band. So if there is an issue for using the multi-scale calculation
%           procedures given below (mainly because the number of time points will be too few when scale increases), you can simply use a low pass filter to filter the data with a cutoff of 0.5 each time when you increase the scale by 1.
%           The composite SampEn map is recommended to be used for across-subject analysis or any group level analysis.
%  Usage: 
% There are 3 steps: 1. loading the preprocessed data and save the intracranial voxels into a data file, 2. call SampEn, 3. read the calculation results and save them in nii file format.
% As can be seen in this file, the input data file has a format of: 1. an integer (2 bytes) indicating the number of voxels or the number of time series
%            2. an integer indicating the length of the time series (number of time points), 
%            3. an vector with each item indicating the standard deviation for each time series (each voxel for fMRI)
%            4. all data saved as time series by voxels, that is    voxel 1's time course    voxel 2' time course ...
%   The output of SampEn is saved as:
%            1. number of voxel (integer)
%            2. m   (int)
%            3. r (float)
%            4. SampEn of all voxels (in float)
%            5. ApEn of all voxels	
%       m can be either 2 or 3,  r can be between 0.3-0.6. 
%  This code can be easily modified to a antomatic batch_code without manually selecting files.  Please try the batch code in ASLtbx, which can be found at
%            https://cfn.upenn.edu/~zewang/    through the menu: software

disp('Calculating BEN map for all subjects, just wait....');

 
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.2*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);
            ndat=dat(mask,:)';
%             gnoi=gauref(1,1:tlen);
            cd(PAR.condirs{sb,c})
            sum=0;
            std_d=std(ndat,0,1)';
            for scale=1
                if scale>0
                    tlens=size(ndat,1);
                    tlens=floor(tlens/2)*2;
                    ndat=ndat(1:tlens,:);
                    ndat=(ndat(1:2:tlens, :)+ndat(2:2:tlens, :))/2;
                end
                std_d=std(ndat,0,1)';
                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);
                              
                 if isunix
	            str=['!SampEn -d ' num2str(PAR.dim) ' -r ' num2str(PAR.r) ' -i ' inmaskdatfile  ' -o SEN.dat '];
	            eval(str);
	          else
	            str=['SampEn.exe -d ' num2str(PAR.dim) ' -r ' num2str(PAR.r) ' -i ' inmaskdatfile  ' -o SEN.dat '];
	            system(str);
	          end      
                    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_BEN_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;
            end
            
            vo.fname=fullfile(ana_dir,sprintf('%s_BEN_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