% Toolbox for calculating Sample Entropy.
% All rights reserved.
% Ze Wang @ TRC, CFN, Upenn 2011
%     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 in an order of:
%            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
	
%%***********************************************************************************************%%
files=spm_select;   % select the preprocessed images
v=spm_vol(files);
[pth,nam,ext,num] = spm_fileparts(files(1,:));
dat=spm_read_vols(v);
mask=squeeze(mean(dat,4));
mask=mask>0.25*max(mask(:));    % generate a brain mask. If you have a better mask, you can disable this line and use your own.
tlen=size(dat,4);
dat=reshape(dat,size(dat,1)*size(dat,2)*size(dat,3),size(dat,4));
inmaskdatfile=fullfile(pth, 'braindat.dat');
dat=dat(mask,:);            % only save the intracranial data
std_d=std(dat,0,2);
dat=dat';                         % reorganize data to be time series by voxels. This is important.
loc=find(mask(:));
r=0.3;    % if the BEN map is nearly zero everywhere, increase r to be 0.6 here
ndat=dat;     
tlen=size(ndat,1);
sum=0;
for dim=[3]
    for scale=0:1
        
        if scale>0
            tlen=floor(tlen/2)*2;
            ndat=(ndat(1:2:tlen,:)+ndat(2:2:tlen,:))/2.;
        end
        tlen=size(ndat,1);
        fprintf('\t\t\t saving inmask data for calc SEN.\n');
        fid=fopen(inmaskdatfile,'wb');
        fwrite(fid, length(loc), 'int');   % number of timecourses
        fwrite(fid, tlen, 'int');

        fwrite(fid, std_d, 'float');
        fwrite(fid, ndat, 'float');
        fclose(fid);
        if isunix
            str=['!SampEn -d ' num2str(dim) ' -r ' num2str(r) ' -i ' inmaskdatfile  ' -o SEN.dat '];
            eval(str);
        else
            str=['SampEn.exe -d ' num2str(dim) ' -r ' num2str(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=zeros(size(mask));
        senmap(mask) = sSEN*1e3;       
        vo=v(1);
        vo.dt=[16 0];
        vo.fname=fullfile(pth,sprintf('SEN_scale_%d_D%d_r%f.nii', scale,dim, r)); %only 1 condition, so no condition name was marked

        vo=spm_write_vol(vo,senmap);
        sum=sum+senmap;
    end
    sum=sum/3;
    vo.fname=fullfile(pth,sprintf('SEN_Composite_D%d_r%f.nii', dim, r));
    vo=spm_write_vol(vo,sum);
end
                