%---------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 4252 $)
%---------------------------------------------------------
%
% preprocessing_job(wrapping_vec,data_path,fieldMapDir,subjects,sessions,vox_size,smooth_vec,b_box)
%
% The function goes through each subject and each session and assigns
% appropriate values and files. The spm job is run for each session.
%
% Glad Mihai, 20.02.2012

function preprocessing_job(data_path,fieldMapDir,subjects,sessions,vars,slice_timing)

% preallocate for speed
funct_data = cell(1,numel(sessions));
session_name = cell(1,numel(sessions));
matlabbatch = cell(1,7);

% compare which direction is desired and choose appropriate vector
if strcmp(vars.wrapping_dir,'x')
    wrapping_vec = [1 0 0];
elseif strcmp(vars.wrapping_dir,'y')
    wrapping_vec = [0 1 0];
elseif strcmp(vars.wrapping_dir,'z')
    wrapping_vec = [0 0 1];

elseif strcmp(vars.wrapping_dir,'xy')
    wrapping_vec = [1 1 0];
elseif strcmp(vars.wrapping_dir,'yx')
    wrapping_vec = [1 1 0];

elseif strcmp(vars.wrapping_dir,'xz')
    wrapping_vec = [1 0 1];
elseif strcmp(vars.wrapping_dir,'zx')
    wrapping_vec = [1 0 1];

elseif strcmp(vars.wrapping_dir,'yz')
    wrapping_vec = [0 1 1];
elseif strcmp(vars.wrapping_dir,'zy')
    wrapping_vec = [0 1 1];
end

% loop over subjects
for i = 1:numel(subjects)
    for j = 1:numel(sessions)
        % load functional data for each session
        funct_data =  spm_select('FPList', fullfile(data_path,subjects{i},'spm',sessions{j}), '^f.*\.img$');
        % load structural image and field map. CAREFUL! it might change over
        % sessions
        structural = spm_select('FPList', fullfile(data_path,subjects{i},'spm','t1vol'), '^s.*\.img$');
        fieldmap = spm_select('FPList', fullfile(data_path,subjects{i},'spm',fieldMapDir), '^vdm5_scs*.*\.img$');
        display(['Field map is ' fieldmap]);
        
        % display which subject and sequence is being processed
        fprintf('Preprocessing subject "%s", "%s" (%s files)\n',...
            subjects{i},sessions{j},sprintf('%d',length(funct_data)));
        
        % CHANGE WORKING DIRECTORY (useful for .ps only)
        %-----------------------------------------------
        matlabbatch{1}.cfg_basicio.cfg_cd.dir = cellstr(fullfile(data_path,subjects{i},'spm'));
        
        % slice timing
        matlabbatch{2}.spm.temporal.st.scans = {cellstr(funct_data)}';
		matlabbatch{2}.spm.temporal.st.nslices = slice_timing.nslices;
		matlabbatch{2}.spm.temporal.st.tr = slice_timing.TR;
		matlabbatch{2}.spm.temporal.st.ta = slice_timing.TR-(slice_timing.TR/slice_timing.nslices);
		matlabbatch{2}.spm.temporal.st.so = slice_timing.sliceOrder;
		matlabbatch{2}.spm.temporal.st.refslice = slice_timing.refSlice;
		matlabbatch{2}.spm.temporal.st.prefix = 'a';
        
        % Realign: Realign and Unwarp
        matlabbatch{3}.spm.spatial.realignunwarp.data.scans(1) = cfg_dep;
        matlabbatch{3}.spm.spatial.realignunwarp.data.scans(1).tname = 'Images';
        matlabbatch{3}.spm.spatial.realignunwarp.data.scans(1).tgt_spec{1}(1).name = 'filter';
        matlabbatch{3}.spm.spatial.realignunwarp.data.scans(1).tgt_spec{1}(1).value = 'image';
        matlabbatch{3}.spm.spatial.realignunwarp.data.scans(1).tgt_spec{1}(2).name = 'strtype';
        matlabbatch{3}.spm.spatial.realignunwarp.data.scans(1).tgt_spec{1}(2).value = 'e';
        matlabbatch{3}.spm.spatial.realignunwarp.data.scans(1).sname = 'Slice Timing: Slice Timing Corr. Images (Sess 1)';
        matlabbatch{3}.spm.spatial.realignunwarp.data.scans(1).src_exbranch = substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1});
        matlabbatch{3}.spm.spatial.realignunwarp.data.scans(1).src_output = substruct('()',{1}, '.','files');
        matlabbatch{3}.spm.spatial.realignunwarp.data.pmscan = cellstr(fieldmap);
        % assign wrapping vector
        matlabbatch{3}.spm.spatial.realignunwarp.eoptions.ewrap = wrapping_vec;
        matlabbatch{3}.spm.spatial.realignunwarp.uwroptions.wrap = wrapping_vec;
        matlabbatch{3}.spm.spatial.realignunwarp.eoptions.quality = 0.9;
        matlabbatch{3}.spm.spatial.realignunwarp.eoptions.sep = 4;
        matlabbatch{3}.spm.spatial.realignunwarp.eoptions.fwhm = 5;
        matlabbatch{3}.spm.spatial.realignunwarp.eoptions.rtm = 0;
        matlabbatch{3}.spm.spatial.realignunwarp.eoptions.einterp = 2;
        matlabbatch{3}.spm.spatial.realignunwarp.eoptions.weight = '';
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.basfcn = [12 12];
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.regorder = 1;
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.lambda = 100000;
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.jm = 0;
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.fot = [4 5];
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.sot = [];
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.uwfwhm = 4;
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.rem = 1;
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.noi = 5;
        matlabbatch{3}.spm.spatial.realignunwarp.uweoptions.expround = 'Average';
        matlabbatch{3}.spm.spatial.realignunwarp.uwroptions.uwwhich = [2 1];
        matlabbatch{3}.spm.spatial.realignunwarp.uwroptions.rinterp = 4;
        matlabbatch{3}.spm.spatial.realignunwarp.uwroptions.mask = 1;
        matlabbatch{3}.spm.spatial.realignunwarp.uwroptions.prefix = 'u';
        
        % COREGISTRATION
        matlabbatch{4}.spm.spatial.coreg.estimate.ref = cfg_dep;
        matlabbatch{4}.spm.spatial.coreg.estimate.ref.tname = 'Reference Image';
        matlabbatch{4}.spm.spatial.coreg.estimate.ref.tgt_spec{1}(1).name = 'filter';
        matlabbatch{4}.spm.spatial.coreg.estimate.ref.tgt_spec{1}(1).value = 'image';
        matlabbatch{4}.spm.spatial.coreg.estimate.ref.tgt_spec{1}(2).name = 'strtype';
        matlabbatch{4}.spm.spatial.coreg.estimate.ref.tgt_spec{1}(2).value = 'e';
        matlabbatch{4}.spm.spatial.coreg.estimate.ref.sname = 'Realign & Unwarp: Unwarped Mean Image';
        matlabbatch{4}.spm.spatial.coreg.estimate.ref.src_exbranch = substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1});
        matlabbatch{4}.spm.spatial.coreg.estimate.ref.src_output = substruct( '.','meanuwr');
        matlabbatch{4}.spm.spatial.coreg.estimate.source = strcat((cellstr(structural)),',1');
        matlabbatch{4}.spm.spatial.coreg.estimate.other = {''};
        matlabbatch{4}.spm.spatial.coreg.estimate.eoptions.cost_fun = 'nmi';
        matlabbatch{4}.spm.spatial.coreg.estimate.eoptions.sep = [4 2];
        matlabbatch{4}.spm.spatial.coreg.estimate.eoptions.tol = [0.02 0.02 0.02 0.001 0.001 0.001 0.01 0.01 0.01 0.001 0.001 0.001];
        matlabbatch{4}.spm.spatial.coreg.estimate.eoptions.fwhm = [7 7];
        
        % SEGMENTATION
        matlabbatch{5}.spm.spatial.preproc.data = cfg_dep;
        matlabbatch{5}.spm.spatial.preproc.data.tname = 'Data';
        matlabbatch{5}.spm.spatial.preproc.data.tgt_spec{1}(1).name = 'filter';
        matlabbatch{5}.spm.spatial.preproc.data.tgt_spec{1}(1).value = 'image';
        matlabbatch{5}.spm.spatial.preproc.data.tgt_spec{1}(2).name = 'strtype';
        matlabbatch{5}.spm.spatial.preproc.data.tgt_spec{1}(2).value = 'e';
        matlabbatch{5}.spm.spatial.preproc.data.sname = 'Coregister: Estimate: Coregistered Images';
        matlabbatch{5}.spm.spatial.preproc.data.src_exbranch = substruct('.','val', '{}',{4}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1});
        matlabbatch{5}.spm.spatial.preproc.data.src_output = substruct('.','cfiles');
        matlabbatch{5}.spm.spatial.preproc.output.GM = [0 0 1];
        matlabbatch{5}.spm.spatial.preproc.output.WM = [0 0 1];
        matlabbatch{5}.spm.spatial.preproc.output.CSF = [0 0 0];
        matlabbatch{5}.spm.spatial.preproc.output.biascor = 1;
        matlabbatch{5}.spm.spatial.preproc.output.cleanup = 0;
        matlabbatch{5}.spm.spatial.preproc.opts.tpm = {
            '/home/glad/.spm/spm8/tpm/grey.nii'
            '/home/glad/.spm/spm8/tpm/white.nii'
            '/home/glad/.spm/spm8/tpm/csf.nii'
            };
        matlabbatch{5}.spm.spatial.preproc.opts.ngaus = [2
            2
            2
            4];
        matlabbatch{5}.spm.spatial.preproc.opts.regtype = 'mni';
        matlabbatch{5}.spm.spatial.preproc.opts.warpreg = 1;
        matlabbatch{5}.spm.spatial.preproc.opts.warpco = 25;
        matlabbatch{5}.spm.spatial.preproc.opts.biasreg = 0.0001;
        matlabbatch{5}.spm.spatial.preproc.opts.biasfwhm = 60;
        matlabbatch{5}.spm.spatial.preproc.opts.samp = 3;
        matlabbatch{5}.spm.spatial.preproc.opts.msk = {''};
        
        % Normalise: Write Functional
        matlabbatch{6}.spm.spatial.normalise.write.subj.matname(1) = cfg_dep;
        matlabbatch{6}.spm.spatial.normalise.write.subj.matname(1).tname = 'Parameter File';
        matlabbatch{6}.spm.spatial.normalise.write.subj.matname(1).tgt_spec{1}(1).name = 'filter';
        matlabbatch{6}.spm.spatial.normalise.write.subj.matname(1).tgt_spec{1}(1).value = 'mat';
        matlabbatch{6}.spm.spatial.normalise.write.subj.matname(1).tgt_spec{1}(2).name = 'strtype';
        matlabbatch{6}.spm.spatial.normalise.write.subj.matname(1).tgt_spec{1}(2).value = 'e';
        matlabbatch{6}.spm.spatial.normalise.write.subj.matname(1).sname = 'Segment: Norm Params Subj->MNI';
        matlabbatch{6}.spm.spatial.normalise.write.subj.matname(1).src_exbranch = substruct('.','val', '{}',{5}, '.','val', '{}',{1}, '.','val', '{}',{1});
        matlabbatch{6}.spm.spatial.normalise.write.subj.matname(1).src_output = substruct('()',{1}, '.','snfile', '()',{':'});
        session_name{j} = strcat ('Realign & Unwarp: Unwarped Images (Sess 	1)' );
        matlabbatch{6}.spm.spatial.normalise.write.subj.resample = cfg_dep;
        matlabbatch{6}.spm.spatial.normalise.write.subj.resample.tname = 'Images to Write';
        matlabbatch{6}.spm.spatial.normalise.write.subj.resample.tgt_spec{1}(1).name = 'filter';
        matlabbatch{6}.spm.spatial.normalise.write.subj.resample.tgt_spec{1}(1).value = 'image';
        matlabbatch{6}.spm.spatial.normalise.write.subj.resample.tgt_spec{1}(2).name = 'strtype';
        matlabbatch{6}.spm.spatial.normalise.write.subj.resample.tgt_spec{1}(2).value = 'e';
        matlabbatch{6}.spm.spatial.normalise.write.subj.resample.sname = session_name{j};
        matlabbatch{6}.spm.spatial.normalise.write.subj.resample.src_exbranch = substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1});
        matlabbatch{6}.spm.spatial.normalise.write.subj.resample.src_output = substruct('.','sess', '()',{1}, '.','uwrfiles');
        matlabbatch{6}.spm.spatial.normalise.write.roptions.preserve = 0;
        matlabbatch{6}.spm.spatial.normalise.write.roptions.bb = vars.b_box;
        matlabbatch{6}.spm.spatial.normalise.write.roptions.vox = vars.vox_size;
        matlabbatch{6}.spm.spatial.normalise.write.roptions.interp = 1;
        matlabbatch{6}.spm.spatial.normalise.write.roptions.wrap = [0 0 0];
        matlabbatch{6}.spm.spatial.normalise.write.roptions.prefix = 'w';
        
        % NORMALISE:WRITE BIAS CORRECTED STRUCTURAL DATA
        matlabbatch{7}.spm.spatial.normalise.write.subj.matname = cfg_dep;
        matlabbatch{7}.spm.spatial.normalise.write.subj.matname.tname = 'Parameter File';
        matlabbatch{7}.spm.spatial.normalise.write.subj.matname.tgt_spec{1}(1).name = 'filter';
        matlabbatch{7}.spm.spatial.normalise.write.subj.matname.tgt_spec{1}(1).value = 'mat';
        matlabbatch{7}.spm.spatial.normalise.write.subj.matname.tgt_spec{1}(2).name = 'strtype';
        matlabbatch{7}.spm.spatial.normalise.write.subj.matname.tgt_spec{1}(2).value = 'e';
        matlabbatch{7}.spm.spatial.normalise.write.subj.matname.sname = 'Segment: Norm Params Subj->MNI';
        matlabbatch{7}.spm.spatial.normalise.write.subj.matname.src_exbranch = substruct('.','val', '{}',{5}, '.','val', '{}',{1}, '.','val', '{}',{1});
        matlabbatch{7}.spm.spatial.normalise.write.subj.matname.src_output = substruct('()',{1}, '.','snfile', '()',{':'});
        matlabbatch{7}.spm.spatial.normalise.write.subj.resample = cfg_dep;
        matlabbatch{7}.spm.spatial.normalise.write.subj.resample.tname = 'Images to Write';
        matlabbatch{7}.spm.spatial.normalise.write.subj.resample.tgt_spec{1}(1).name = 'filter';
        matlabbatch{7}.spm.spatial.normalise.write.subj.resample.tgt_spec{1}(1).value = 'image';
        matlabbatch{7}.spm.spatial.normalise.write.subj.resample.tgt_spec{1}(2).name = 'strtype';
        matlabbatch{7}.spm.spatial.normalise.write.subj.resample.tgt_spec{1}(2).value = 'e';
        matlabbatch{7}.spm.spatial.normalise.write.subj.resample.sname = 'Segment: Bias Corr Images';
        matlabbatch{7}.spm.spatial.normalise.write.subj.resample.src_exbranch = substruct('.','val', '{}',{5}, '.','val', '{}',{1}, '.','val', '{}',{1});
        matlabbatch{7}.spm.spatial.normalise.write.subj.resample.src_output = substruct('.','biascorr', '()',{':'});
        matlabbatch{7}.spm.spatial.normalise.write.roptions.preserve = 0;
        matlabbatch{7}.spm.spatial.normalise.write.roptions.bb = vars.b_box;
        matlabbatch{7}.spm.spatial.normalise.write.roptions.vox = [1 1 1];
        matlabbatch{7}.spm.spatial.normalise.write.roptions.interp = 1;
        % no wrapping
        matlabbatch{7}.spm.spatial.normalise.write.roptions.wrap = [0 0 0];
        matlabbatch{7}.spm.spatial.normalise.write.roptions.prefix = 'w';
        
        % SMOOTHING
        matlabbatch{8}.spm.spatial.smooth.data = cfg_dep;
        matlabbatch{8}.spm.spatial.smooth.data.tname = 'Images to Smooth';
        matlabbatch{8}.spm.spatial.smooth.data.tgt_spec{1}(1).name = 'filter';
        matlabbatch{8}.spm.spatial.smooth.data.tgt_spec{1}(1).value = 'image';
        matlabbatch{8}.spm.spatial.smooth.data.tgt_spec{1}(2).name = 'strtype';
        matlabbatch{8}.spm.spatial.smooth.data.tgt_spec{1}(2).value = 'e';
        matlabbatch{8}.spm.spatial.smooth.data.sname = 'Normalise: Write: Normalised Images (Subj 1)';
        matlabbatch{8}.spm.spatial.smooth.data.src_exbranch = substruct('.','val', '{}',{6}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1});
        matlabbatch{8}.spm.spatial.smooth.data.src_output = substruct( '.','files');
        matlabbatch{8}.spm.spatial.smooth.fwhm = vars.smooth_vec;
        matlabbatch{8}.spm.spatial.smooth.dtype = 0;
        matlabbatch{8}.spm.spatial.smooth.im = 0;
        matlabbatch{8}.spm.spatial.smooth.prefix = 's';
        
        %% SAVE AND RUN JOB
        %------------------
        save(fullfile(data_path,subjects{i},'spm',strcat(sessions{j},'_preprocessing.mat')),'matlabbatch');
        spm_jobman('serial', matlabbatch, '' );
        clear matlabbatch
    end
end

