function batch_spm2_VOI(VOIname,VOIdef,VOIradius,nsessions,prefix,amr_prefix,amr_postfix,first_lev_dir,rootDIR,directories,Contrast_names,VOIx,VOIy,VOIz)
% script for multisubject VOI extraction
%=======================================================================
%             ________  ____  ____  ___
%            / ___/ _ \/ __ \/ __ \/ _ \
%           (__  )  __/ /_/ / /_/ /  __/
%          /____/\___/ .___/ .___/\___/
%                   /_/   /_/
%
% script written by Seppe Santens, 99% of the code stolen from Friston
%=======================================================================

%rootDIR         = '/media/Data/fMRI/data';
%data_dir = 'first_level_no_rp_m_n1_s8';
%prefix = 'Sub_';
% make sure the scriptdir is in the path
addpath(pwd);

%[directories]={'Stats_sub_17_nback','Stats_sub_217_nback','Stats_sub_223_nback','Stats_sub_243_nback'};

%[directories]={'Stats_sub_9_nback', 'Stats_sub_17_nback', 'Stats_sub_44_nback', 'Stats_sub_124_nback', 'Stats_sub_217_nback', 'Stats_sub_223_nback', 'Stats_sub_243_nback', 'Stats_sub_259_nback', 'Stats_sub_269_nback', 'Stats_sub_289_nback', 'Stats_sub_364_nback', 'Stats_sub_379_nback', 'Stats_sub_380_nback', 'Stats_sub_400_nback', 'Stats_sub_409_nback', 'Stats_sub_441_nback', 'Stats_sub_458_nback', 'Stats_sub_487_nback', 'Stats_sub_530_nback', 'Stats_sub_611_nback', 'Stats_sub_631_nback', 'Stats_sub_634_nback', 'Stats_sub_670_nback', 'Stats_sub_676_nback', 'Stats_sub_687_nback', 'Stats_sub_706_nback', 'Stats_sub_717_nback', 'Stats_sub_763_nback', 'Stats_sub_764_nback', 'Stats_sub_780_nback', 'Stats_sub_796_nback', 'Stats_sub_797_nback', 'Stats_sub_801_nback', 'Stats_sub_802_nback', 'Stats_sub_809_nback', 'Stats_sub_866_nback', 'Stats_sub_881_nback', 'Stats_sub_929_nback', 'Stats_sub_933_nback', 'Stats_sub_955_nback', 'Stats_sub_973_nback', 'Stats_sub_980_nback', 'Stats_sub_1000_nback', 'Stats_sub_1007_nback', 'Stats_sub_1026_nback', 'Stats_sub_1035_nback', 'Stats_sub_1092_nback', 'Stats_sub_1095_nback', 'Stats_sub_1096_nback', 'Stats_sub_1141_nback', 'Stats_sub_1188_nback', 'Stats_sub_1198_nback', 'Stats_sub_1433_nback', 'Stats_sub_5011_nback', 'Stats_sub_5023_nback', 'Stats_sub_5093_nback', 'Stats_sub_5098_nback', 'Stats_sub_5102_nback', 'Stats_sub_5107_nback', 'Stats_sub_5116_nback', 'Stats_sub_5127_nback', 'Stats_sub_5141_nback', 'Stats_sub_5149_nback', 'Stats_sub_5157_nback', 'Stats_sub_5158_nback', 'Stats_sub_5159_nback', 'Stats_sub_5166_nback', 'Stats_sub_5167_nback', 'Stats_sub_5180_nback', 'Stats_sub_5181_nback', 'Stats_sub_5191_nback', 'Stats_sub_5196_nback', 'Stats_sub_5212_nback', 'Stats_sub_5223_nback', 'Stats_sub_5230_nback', 'Stats_sub_5232_nback', 'Stats_sub_5238_nback', 'Stats_sub_5250_nback', 'Stats_sub_5256_nback', 'Stats_sub_5258_nback', 'Stats_sub_5273_nback', 'Stats_sub_5282_nback'};
%[directories]={'Sub_17'};

%'Stats_sub_9_nback', 
%'Stats_sub_400_nback', 

%directories = sort(directories);

% set parameters
%nsessions = 8; %this refers to the number of contrasts, at least here
nsubjects = length(directories);

% Contrast_names={'Positive_rACC_Group','Neutral_rACC_Group','Negative_rACC_Group','NegVNeu_rACC_Group','PosVNeu_rACC_Group','Valence_rACC_Group','Arousal_rACC_Group';};
% VOIx = 4;                                         % x-coordinate
% VOIy = 52;                                         % y-coordinate
% VOIz = 16;                                          % z-coordinate

VOIxyz = [VOIx VOIy VOIz];

%VOIname = 'test';
%VOIdef = 'sphere';                                  % type of VOI
%VOIradius = 8;                                      % radius of sphere (mm)
% directories=dir(rootDIR);

% create a subject loop
subj = 1:nsubjects;
for x = 1:nsubjects
    
    
    
    %fprintf('Working on participant %s\n',int2str(directories(x)));

    %SJ = sprintf('%s',int2str(directories(x)));
    % take subject directory
    %swd = [rootDIR filesep int2str(directories(x))];
    %swd = fullfile(rootDIR,directories{x});
    %change to the subjects directory
    
    swd = fullfile(rootDIR,first_lev_dir,strcat(amr_prefix,num2str(directories(x)),amr_postfix));
    cd(swd);
   
    %=======================================================================
    % - C O P Y   P A S T I N G   F R O M   spm_getSPM.m
    %=======================================================================

    %-Load SPM.mat
    %-----------------------------------------------------------------------

    load(fullfile(swd,'SPM.mat'));
    
%     for i = 1:480
%         
%         SPM.xY.VY(i).fname = strrep(SPM.xY.VY(i).fname,'\','/');
%         SPM.xY.VY(i).fname = strrep(SPM.xY.VY(i).fname,'D:/EWM_1/N-back_Task/fMRI_data_nback/Phase_3','/media2/Data/fmri/data');
%         SPM.xY.VY(i).fname = strrep(SPM.xY.VY(i).fname,'G:/EWM_1/N-back_Task/fMRI_data_nback/Phase_3','/media2/Data/fmri/data');
%         SPM.xY.VY(i).fname = strrep(SPM.xY.VY(i).fname,'E:/EWM_1/N-back_Task/fMRI_data_nback/Phase_3','/media2/Data/fmri/data');
                
%         str = SPM.xY.P(i,:);
%         str = strrep(str,'\','/');
%         str = strrep(str,'//','/');
%         str = strrep(str,'D:/EWM_1/N-back_Task/fMRI_data_nback/Phase_3','/media2/Data/fmri/data');
%         str = strrep(str,'G:/EWM_1/N-back_Task/fMRI_data_nback/Phase_3','/media2/Data/fmri/data');
%         str = strrep(str,'E:/EWM_1/N-back_Task/fMRI_data_nback/Phase_3','/media2/Data/fmri/data');
%         
%         len = length(str);
%         
%         SPM.xY.P(i,1:len)=str;
        
        %SPM.xY.P(i,:) = strrep(SPM.xY.P(i,:),'\','/');
        %SPM.xY.P(i,:) = strrep(SPM.xY.P(i,:),'//','/');
        %SPM.xY.P(i,:) = strrep(SPM.xY.P(i,:),'D:/EWM_1/N-back_Task/fMRI_data_nback/Phase_3','/media2/Data/fmri/data');
        
        %SPM.xY.P(i,:) = temp_str;
%     end  
       
    save('SPM.mat','SPM');
    
    SPM.swd = swd;

    %-Get volumetric data from SPM.mat
    %-----------------------------------------------------------------------
    try
        xX   = SPM.xX;                  %-Design definition structure
        XYZ  = SPM.xVol.XYZ;			%-XYZ coordinates
        S    = SPM.xVol.S;              %-search Volume {voxels}
        R    = SPM.xVol.R;              %-search Volume {resels}
        M    = SPM.xVol.M(1:3,1:3);		%-voxels to mm matrix
        VOX  = sqrt(diag(M'*M))';		%-voxel dimensions
    catch
        error ('This model has not been estimated.');
    end

    %-----------------------------------------------------------------------
    % - SKIPPING ALL (HOPEFULLY) UNNECESSARY STEPS
    %-----------------------------------------------------------------------

    %-Assemble ONLY NECESSARY output structures of unfiltered data
    %-----------------------------------------------------------------------
    xSPM   = struct('swd',swd,... % SPM working directory                                        
                    'XYZ',XYZ,... % location of voxels {voxel coords}
                    'XYZmm',SPM.xVol.M(1:3,:)*[XYZ; ones(1,size(XYZ,2))]);  % location of voxels {mm}

    for VOIsession=1:nsessions

        fprintf('\tWorking on session %s\n',char(Contrast_names(VOIsession)));
        VOIname=char(Contrast_names(VOIsession));

        
%=======================================================================
        % - C O P Y   P A S T I N G   F R O M   spm_regions.m
        
%=======================================================================

        % xY     - VOI structure
        
%-----------------------------------------------------------------------
        try
            xY;
        catch
            xY = {};
        end

        %-Find nearest voxel [Euclidean distance]
        
%-----------------------------------------------------------------------
        if isempty(xSPM.XYZmm)
            spm('alert!','No suprathreshold voxels!',mfilename,0);
            Y = []; xY = [];
            return
        end
        [xyz,i] = spm_XYZreg('NearestXYZ',VOIxyz,xSPM.XYZmm);
        xY.xyz  = xyz;


        %-Get adjustment options and VOI name
        %------------------------------------------------------------------
        xY.name = VOIname;

        xY.Ic = 1;                      % effects of interest

        %-if fMRI data get sessions and filtering options
        
%-----------------------------------------------------------------------
        if isfield(SPM,'Sess')

            s         = length(SPM.Sess);
            if s > 1
                s = VOIsession;
            end
            xY.Sess   = s;

        end

        %-Specify VOI
        
%-----------------------------------------------------------------------
        if ~isfield(xY,'def')
            xY.def    = VOIdef;
        end
        Q       = ones(1,size(xSPM.XYZmm,2));

        switch xY.def

            case 'sphere'
            %---------------------------------------------------------------
            if ~isfield(xY,'spec')
                xY.spec = VOIradius;
            end
            d     = [xSPM.XYZmm(1,:) - xyz(1);
                 xSPM.XYZmm(2,:) - xyz(2);
                 xSPM.XYZmm(3,:) - xyz(3)];
            Q     = find(sum(d.^2) <= xY.spec^2);

            case 'box'
            %---------------------------------------------------------------
            if ~isfield(xY,'spec')
                xY.spec = spm_input('box dimensions [x y z] {mm}',...
                    '!+0','r','0 0 0',3);
            end
            Q     = find(all(abs(xSPM.XYZmm - xyz*Q) <= xY.spec(:)*Q/2));

            case 'cluster'
            %---------------------------------------------------------------
            [x,i] = spm_XYZreg('NearestXYZ',xyz,xSPM.XYZmm);
            A     = spm_clusters(xSPM.XYZ);
            Q     = find(A == A(i));
        end

        %-Extract required data from results files
        
%=======================================================================

        %-Get raw data, whiten and filter
        
%-----------------------------------------------------------------------
        y        = spm_get_data(SPM.xY.VY,xSPM.XYZ(:,Q));
        y        = spm_filter(SPM.xX.K,SPM.xX.W*y);
        xY.XYZmm = xSPM.XYZmm(:,Q);

        %-Computation
        
%=======================================================================

        % remove null space of contrast
        
%-----------------------------------------------------------------------
        if xY.Ic

            %-Parameter estimates: beta = xX.pKX*xX.K*y
            %---------------------------------------------------------------
            beta  = spm_get_data(SPM.Vbeta,xSPM.XYZ(:,Q));

            %-subtract Y0 = XO*beta,  Y = Yc + Y0 + e
            %---------------------------------------------------------------
            y     = y - spm_FcUtil('Y0',SPM.xCon(xY.Ic),SPM.xX.xKXs,beta);

        end

        % confounds
        
%-----------------------------------------------------------------------
        xY.X0     = SPM.xX.xKXs.X(:,[SPM.xX.iB SPM.xX.iG]);

        % extract session-specific rows from data and confounds
        
%-----------------------------------------------------------------------
        try
            i     = SPM.Sess(xY.Sess).row;
            y     = y(i,:);
            xY.X0 = xY.X0(i,:);
        end

        % and add session-specific filter confounds
        
%-----------------------------------------------------------------------
        try
            xY.X0 = [xY.X0 SPM.xX.K(xY.Sess).X0];
        end

        
%=======================================================================
        try
            xY.X0 = [xY.X0 SPM.xX.K(xY.Sess).KH]; % Compatibility check
        end
        
%=======================================================================


        % Remove null space of X0
        
%-----------------------------------------------------------------------
        xY.X0   = xY.X0(:,~~any(xY.X0));


        % compute regional response in terms of first eigenvariate
        
%-----------------------------------------------------------------------
        [m n]   = size(y);
        if m > n
            [v s v] = svd(spm_atranspa(y));
            s       = diag(s);
            v       = v(:,1);
            u       = y*v/sqrt(s(1));
        else
            [u s u] = svd(spm_atranspa(y'));
            s       = diag(s);
            u       = u(:,1);
            v       = y'*u/sqrt(s(1));
        end
        d       = sign(sum(v));
        u       = u*d;
        v       = v*d;
        Y       = u*sqrt(s(1)/n);

        % set in structure
        
%-----------------------------------------------------------------------
        xY.y    = y;
        xY.u    = Y;
        xY.v    = v;
        xY.s    = s;

        % save
        
%-----------------------------------------------------------------------
        str     = ['VOI_' xY.name];
        if isfield(xY,'Sess')
            if length(xY.Sess) == 1
                str = sprintf('VOI_%s_%i',xY.name,xY.Sess);
            end
        end
        fprintf('\t\tsaving file %s\n',str);
        
        temp_path = fullfile(strcat(rootDIR,filesep,'PPI',filesep,prefix,num2str(directories(x)),filesep,char(Contrast_names(VOIsession))),str);
        try
            %mkdir(temp_path{1})
            mkdir(temp_path)
        catch
            1;
        end
        
        %save(temp_path{1},'Y','xY');
        save(temp_path,'Y','xY');

    end

    cd (rootDIR);
end

clear all;
end



