% Spatial Preprocessing for fMRI data
% Christian Kaul, July 2007
% presented @ Matlab for Cogntive Neuroscience 12/07/07
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify each subject and the number of runs acquired for them.
% this script assumes a folder structure like this:
%
% exp_name/
% exp_name/s1
% exp_name/s1/structural           % the structural folder contains
% structural, fieldmaps, segmentation,...
% exp_name/s1/structural/ROI       %the ROI folder contains all the area defs, ergo put your amygdala def here
%
% exp_name/s1/localizer_data
% exp_name/s1/localizer_data/run1  % contains functional files
% exp_name/s1/localizer_data/run2
% exp_name/s1/localizer_stats      % contains the SPM, betas, con-images
%
% exp_name/s1/experiment
% exp_name/s1/experiment/run1      % contains functional files
% exp_name/s1/experiment/run2
% exp_name/s1/experiment_stats     % contains the SPM, betas, con-images


% this batch script with set perform the following chain of pre-processing
% step:
% 1) Realignment of all functional runs  (realign: estimate & reslice)
%    1.1 alternatively you can realign and unwarp. see in code.
% 2) Coregister the structural to mean-epi produced during 1)
% 3) Segment the coregistered Structural producing a gray matter file, a
%    white matter file and the parameters to normalize epi in
% 4) Normalize all functional runs  (write only)
% 5) if desired smooth data with smooth section
%

clear all

% SPM results directory for each subject is
statsDirName = 'exp_stats';
% Specify root directory
origDir  = 'D:\V1load_motion2';

dataFolder = 'experiment\run';
structuralFolder = 'structural';

% subjects and the number of volumes collected for each run and the name of
% the structural file
subjects = { ...
    '1',  ones(1,48)*48-3, 'ED-MA04903.img'; ...
    '2',  ones(1,48)*48-3, 'GV-MA04904.img'; ...
    '3',  ones(1,48)*48-3, 'HR-sMA04830.img'; ...
    '4',  ones(1,48)*48-3, 'DB-meansMA04771.img'; ...
    '5',  ones(1,48)*48-3, 'AD.img'; ...
    };
% =========================================================================
cd(origDir);
spm fmri
% =========================================================================
% =========================================================================
% =========================================================================
% =========================== Realignment =================================
% =========================================================================
clc;
disp(['Realign: Estimate and Reslice for MEAN and ALL images ']);
for s0 = 1 : size(subjects,1)
    disp(['Realign job specification for Subject : ', subjects{s0,1}]);
    runs = size(subjects{s0,2},2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % first specify the f* scans for each run
    for run = 1 : runs
        %% define epi files in the run
        epiDir = [origDir '\s' subjects{s0,1} '\' dataFolder int2str(run)];
        %%% select scans and assign to job
        f   = spm_select('List', epiDir, '^fM.*\.img$');
        fs  = cellstr([repmat([epiDir '\'],size(f,1),1) f]);
        jobs{1}.spatial{1}.realign{1}.estwrite.data{run} = fs;
        %%%% clear temporary variables for next run
        f = []; fs = [];
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % specify Estimation Options
    % all default values taken from SPM5 batch
    jobs{1}.spatial{1}.realign{1}.estwrite.eoptions.quality = 0.9;
    jobs{1}.spatial{1}.realign{1}.estwrite.eoptions.sep     = 4.0;
    jobs{1}.spatial{1}.realign{1}.estwrite.eoptions.fwhm    = 5.0;
    jobs{1}.spatial{1}.realign{1}.estwrite.eoptions.rtm     = 1.0;
    jobs{1}.spatial{1}.realign{1}.estwrite.eoptions.interp  = 2.0;
    jobs{1}.spatial{1}.realign{1}.estwrite.eoptions.wrap    = [0 0 0];
    jobs{1}.spatial{1}.realign{1}.estwrite.eoptions.weight  = {};
    % specify Reslice Options
    jobs{1}.spatial{1}.realign{1}.estwrite.roptions.which   = [2 1]; % mean and all
    jobs{1}.spatial{1}.realign{1}.estwrite.roptions.interp  = 4.0;
    jobs{1}.spatial{1}.realign{1}.estwrite.roptions.wrap    = [0 0 0];
    jobs{1}.spatial{1}.realign{1}.estwrite.roptions.mask    = 1.0;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % save and run job
    save realign.mat jobs
    disp(['RUNNING Realign for subject ' subjects{s0,1}]);
    spm_jobman('run','realign.mat');
    clear jobs
end

% %  % alternative:
% %  % if you want to realign & unwrap you will have to specify the
% %  % fieldmaps in subjects. like this:

% % subjects = { ...
% %          '7', [110]-6  , 'sMA04752-0012-00001-000176.img', 'vdm5_scsMA04752-0010-00001-000001-01.img'; ...
% %          };

% % % =========================================================================
% % cd(origDir);
% % spm fmri
% % % =========================================================================
% % % =========================================================================
% % % =========================================================================
% % % =========================== Realignment UnWarp =================================
% % % =========================================================================
% % clc;
% % disp(['Realign: Realignment & UnWarp for MEAN and ALL images ']);
% % for s0 = 1 : size(subjects,1)
% %     disp(['Realign job specification for Subject : ', subjects{s0,1}]);
% %     runs = size(subjects{s0,2},2);
% %     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %     % first specify the f* scans for each run
% %     for run = 1 : runs
% %         %% define epi files in the run
% %         epiDir = [origDir '\s' subjects{s0,1} '\' dataFolder int2str(run)];
% %         %%% select scans and assign to job
% %         f   = spm_select('List', epiDir, '^fM.*\.img$');
% %         fs  = cellstr([repmat([epiDir '\'],size(f,1),1) f]);
% %         jobs{1}.spatial{1}.realignunwarp.data.scans{run} = fs;
% %         %%%% clear temporary variables for next run
% %         f = []; fs = [];
% %     end
% %     jobs{1}.spatial{1}.realignunwarp.data.pmscan{1} = [origDir '\s' subjects{s0,1} '\' structuralFolder '\' subjects{s0,4}];
% %     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %     % specify Estimation Options
% %     % all default values taken from SPM5 batch
% %
% %
% %     jobs{1}.spatial{1}.realignunwarp.eoptions.quality = 0.9;
% %     jobs{1}.spatial{1}.realignunwarp.eoptions.sep     = 4.0;
% %     jobs{1}.spatial{1}.realignunwarp.eoptions.fwhm    = 5.0;
% %     jobs{1}.spatial{1}.realignunwarp.eoptions.rtm     = 0;
% %     jobs{1}.spatial{1}.realignunwarp.eoptions.einterp = 2.0;
% %     jobs{1}.spatial{1}.realignunwarp.eoptions.wrap    = [0 0 0];
% %     jobs{1}.spatial{1}.realignunwarp.eoptions.weight  = {};
% %
% %     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.basfcn   = [12 12];
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.regorder = 1.0;
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.lambda   = 1e+005;
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.jm       = 0;
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.fot      = [4 5];
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.sot      = [];
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.uwfwhm   = 4;
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.rem      = 1;
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.noi      = 5;
% %     jobs{1}.spatial{1}.realignunwarp.uweoptions.expround = 'Average';
% %
% %     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% %     % specify Reslice Options
% %     jobs{1}.spatial{1}.realignunwarp.uwroptions.which   = [2 1]; % mean and all
% %     jobs{1}.spatial{1}.realignunwarp.uwroptions.interp  = 4.0;
% %     jobs{1}.spatial{1}.realignunwarp.uwroptions.wrap    = [0 0 0];
% %     jobs{1}.spatial{1}.realignunwarp.uwroptions.mask    = 1.0;
% %     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %     % save and run job
% %     save realignunwarp.mat jobs
% %     disp(['RUNNING Realign for subject ' subjects{s0,1}]);
% %     spm_jobman('run','realignunwarp.mat');
% %     clear jobs
% % end



% =========================================================================
% ========================== Coregistration ===============================
% =========================================================================
% in this case we only want to coregister the structural to the mean epi
clc;
disp(['Coregister: Estimate structural (change postition file (.hdr) to match postition of epi-mean file)']);
for s0 = 1 : size(subjects,1)
    disp(['Coregister job specification for Subject : ', subjects{s0,1}]);
    %runs = size(subjects{s0,2},2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % first specify the reference image: structural scan
    sourceDir = [origDir '\s' subjects{s0,1} '\' dataFolder '1'];
    f  = spm_select('List', sourceDir, '^mean.*\.img$');
    fs = [sourceDir '\', f ',1'];
    jobs{1}.spatial{1}.coreg{1}.estimate.ref = {fs};   % select mean epi as reference
    f = []; fs = [];
    stDir = [origDir '\s' subjects{s0,1} '\' structuralFolder '\'];
    stFile = [stDir subjects{s0,3} ',1'];
    jobs{1}.spatial{1}.coreg{1}.estimate.source = {stFile};   % select strutural as source
    stFile = [];
    jobs{1}.spatial{1}.coreg{1}.estimate.other = {};    % in this case no other files
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    jobs{1}.spatial{1}.coreg{1}.estimate.eoptions.cost_fun = 'nmi';
    jobs{1}.spatial{1}.coreg{1}.estimate.eoptions.sep      = [4 2];
    jobs{1}.spatial{1}.coreg{1}.estimate.eoptions.tol      = [ 0.0200 0.0200 0.0200 0.0010 0.0010 0.0010 ...
        0.0100 0.0100 0.0100 0.0010 0.0010 0.0010];
    jobs{1}.spatial{1}.coreg{1}.estimate.eoptions.fwhm     = [7 7];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %     jobs{1}.spatial{1}.coreg{1}.estwrite.roptions.interp = 1;
    %     jobs{1}.spatial{1}.coreg{1}.estwrite.roptions.wrap = [0 0 0];
    %     jobs{1}.spatial{1}.coreg{1}.estwrite.roptions.mask = 0;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % save and run job
    save coreg_structural.mat jobs
    disp(['RUNNING Coreg for subject ' subjects{s0,1}]);
    spm_jobman('run','coreg_structural.mat');
    clear jobs
end

% =========================================================================
% =========================================================================
% =========================================================================
% =========================== Segmentation =================================
% =========================================================================
clc;
disp(['Segmentation: Produce gray matter (native & modulated normalized) & withe matter (native) from structural image ']);
spm5path = 'D:\maltlabspace\spm5\';
for s0 = 1 : size(subjects,1)
    disp(['Segmentation job specification for Subject : ', subjects{s0,1}]);


    stDir = [origDir '\s' subjects{s0,1} '\' structuralFolder '\'];
    stFilePath = {[stDir subjects{s0,3} ',1']};  % notice that this is the COREGISTERED structural
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % all default values taken from SPM5 batch
    jobs{1}.spatial{1}.preproc.data = stFilePath;
    jobs{1}.spatial{1}.preproc.output.GM     = [1 0 1];   % gray   native + modulated normalized
    jobs{1}.spatial{1}.preproc.output.WM     = [0 0 1];   % white
    jobs{1}.spatial{1}.preproc.output.CSF    = [0 0 0];
    jobs{1}.spatial{1}.preproc.output.biascor  = 1;
    jobs{1}.spatial{1}.preproc.output.cleanup  = 1;     % this option should get rid of unlikely structures like handels or holes
    % specify Segmentation Options
    jobs{1}.spatial{1}.preproc.opts.tpm  = {[spm5path 'tpm\grey.nii'];[spm5path 'tpm\white.nii'];[spm5path 'tpm\csf.nii'];};
    jobs{1}.spatial{1}.preproc.opts.ngaus     = [2 2 2 4];
    jobs{1}.spatial{1}.preproc.opts.regtype   = 'mni';
    jobs{1}.spatial{1}.preproc.opts.warpreg   = 1;
    jobs{1}.spatial{1}.preproc.opts.warpco    = 25;
    jobs{1}.spatial{1}.preproc.opts.biasreg   = 0.0001;
    jobs{1}.spatial{1}.preproc.opts.biasfwhm  = 60;
    jobs{1}.spatial{1}.preproc.opts.samp      = 3;
    jobs{1}.spatial{1}.preproc.opts.msk       = {};
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % save and run job
    save segment.mat jobs
    disp(['RUNNING Segmentation for subject ' subjects{s0,1}]);
    spm_jobman('run','segment.mat');
    clear jobs
end

clc;



% % =========================================================================
% % =========================================================================
% % ======================== Smoothing ======================================
% % =========================================================================
% clc;
% disp(['Smoothing ... ']);
% sFWHM = input('The Smoothing Kernel size [x y z]  = ');
% for s0 = 1 : size(subjects,1)
%     disp(['smoothing job specification for Subject : ', subjects{s0,1}]);
%     runs = size(subjects{s0,2},2);
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     % First specify the files: rf* scans for
%     conCat_fs = [];
%     for run = 1 : runs
%         %% define epi files in the run
%         epiDir = [origDir '\s' subjects{s0,1} '\' dataFolder int2str(run)];
%         %%% select scans and assign to job
%         f   = spm_select('List', epiDir, '^rfM.*\.img$');
%         fs  = cellstr([repmat([epiDir '\'],size(f,1),1) f repmat([',1'],size(f,1),1)]);
%         conCat_fs = [conCat_fs; fs];
%         %%%% clear temporary variables for next run
%         f = []; fs = [];
%     end
%     jobs{1}.spatial{1}.smooth.data  = conCat_fs;
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     jobs{1}.spatial{1}.smooth.fwhm  = sFWHM;
%     jobs{1}.spatial{1}.smooth.dtype = 0; % same
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     % save and run job
%     save smooth.mat jobs
%     disp(['RUNNING Realign for subject ' subjects{s0,1}]);
%     spm_jobman('run','smooth.mat');
%     clear jobs
% end



% =========================================================================
% =========================================================================
% ======================  Normalization ===================================
% =========================================================================
% here we load the normalization parameter file produced during
% segmentation and normalize all epi data with it.
clc;
disp(['Normlizing ... ']);
templateFile = 'C:\MATLAB7\toolbox\spm5\templates\EPI.nii';
for s0 = 1 : size(subjects,1)
    disp(['Normalizing job specification for Subject : ', subjects{s0,1}]);
    runs = size(subjects{s0,2},2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % first specify the parameter file (..._sn.mat) as produced by sementation
    snDir = [origDir '\s' subjects{s0,1} '\' structuralFolder '\'];
    [pathstr, Sname] = fileparts(subjects{s0,3});
    snFile = [snDir  Sname '_seg_sn.mat' ];
    jobs{1}.spatial{1}.normalise{1}.write.subj.matname = {snFile};
    snFile = [];
    % the specify the other images: the (probably smoothed) vols all runs run
    conCat_fs = [];
    for run = 1 : runs
        %% define epi files in the run
        epiDir = [origDir '\s' subjects{s0,1} '\' dataFolder int2str(run)];
        %%% select scans and assign to job
        f   = spm_select('List', epiDir, '^rfM.*\.img$');  %add an 's' if smoothed
        fs  = cellstr([repmat([epiDir '\'],size(f,1),1) f repmat([',1'],size(f,1),1)]);
        conCat_fs = [conCat_fs; fs];
        %%%% clear temporary variables for next run
        f = []; fs = [];
    end
    jobs{1}.spatial{1}.normalise{1}.write.subj.resample = conCat_fs;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %no estimation in this case, since already done by segmentation
    %     jobs{1}.spatial{1}.normalise{1}.estwrite.eoptions.template = {templateFile};
    %     jobs{1}.spatial{1}.normalise{1}.estwrite.eoptions.weight   = {};
    %     jobs{1}.spatial{1}.normalise{1}.estwrite.eoptions.smosrc   = 8;
    %     jobs{1}.spatial{1}.normalise{1}.estwrite.eoptions.smoref   = 0;
    %     jobs{1}.spatial{1}.normalise{1}.estwrite.eoptions.regtype  = 'mni';
    %     jobs{1}.spatial{1}.normalise{1}.estwrite.eoptions.cutoff   = 25;
    %     jobs{1}.spatial{1}.normalise{1}.estwrite.eoptions.nits     = 16;
    %     jobs{1}.spatial{1}.normalise{1}.estwrite.eoptions.reg      = 1;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    jobs{1}.spatial{1}.normalise{1}.write.roptions.preserve = 0;
    jobs{1}.spatial{1}.normalise{1}.write.roptions.bb       = [ -78 -112 -50; 78 76 85];
    jobs{1}.spatial{1}.normalise{1}.write.roptions.vox      = [3 3 3];  % recommended to keep resolution of your data here
    jobs{1}.spatial{1}.normalise{1}.write.roptions.interp   = 1;
    jobs{1}.spatial{1}.normalise{1}.write.roptions.wrap     = [0 0 0];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % save and run job
    save normalise.mat jobs
    disp(['RUNNING Realign for subject ' subjects{s0,1}]);
    spm_jobman('run','normalise.mat');
    clear jobs
end
%=========================================================================
%=========================================================================
%=========================================================================