% preproc_vbm - preprocessing for voxel-based morphometry
%
% Need to specify:
% * Parent directory for directory of Template, average and mask images
% * Structural images for all subjects
% * A filename for the spreadsheet (CSV) of total tissue volumes
%
% PLEASE NOTE:
% There is a minor bug in this version: the file Template_6_2mni.mat will
% be left in the directory of the first selected structural image; not
% moved to the Templates_etc directory along with the Template image files.
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Ged Ridgway
% $Id: preproc_vbm.m 6177 2014-09-16 10:44:55Z guillaume $

if exist('preproc_vbm', 'file') && ~isdeployed
    help preproc_vbm
end
matlabbatch = {};

%% Make directory to hold Dartel Templates, average image, mask
iDir = length(matlabbatch) + 1;
matlabbatch{iDir}.cfg_basicio.file_dir.dir_ops.cfg_mkdir.parent = '<UNDEFINED>';
matlabbatch{iDir}.cfg_basicio.file_dir.dir_ops.cfg_mkdir.name = 'Templates_etc';

%% Segmentation
iSeg = length(matlabbatch) + 1;
matlabbatch{iSeg}.spm.spatial.preproc.channel.write = [0 1];
ngaus  = [1 1 2 3 4 2];
native = [
    1 0 0 0 0 0
    1 1 0 0 0 0];
warped = [
    0 0 0 0 0 0
    1 1 1 0 0 0];
for c = 1:6 % tissue class c
    matlabbatch{iSeg}.spm.spatial.preproc.tissue(c).tpm = {
        fullfile(spm('dir'), 'tpm', sprintf('TPM.nii,%d', c))};
    matlabbatch{iSeg}.spm.spatial.preproc.tissue(c).ngaus = ngaus(c);
    matlabbatch{iSeg}.spm.spatial.preproc.tissue(c).native = native(:, c)';
    matlabbatch{iSeg}.spm.spatial.preproc.tissue(c).warped = warped(:, c)';
end

%% Dartel - create template
iDar = length(matlabbatch) + 1;
matlabbatch{iDar}.spm.tools.dartel.warp.images{1}(1) = ...
    cfg_dep('Segment: rc1 Images', substruct('.','val', '{}',{iSeg}, ...
    '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','tiss', '()',{1}, '.','rc', '()',{':'}));
matlabbatch{iDar}.spm.tools.dartel.warp.images{2}(1) = ...
    cfg_dep('Segment: rc2 Images', substruct('.','val', '{}',{iSeg}, ...
    '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','tiss', '()',{2}, '.','rc', '()',{':'}));

%% Normalise segmentations to MNI space
iNS = length(matlabbatch) + 1;
matlabbatch{iNS}.spm.tools.dartel.mni_norm.template(1) = ...
    cfg_dep('Run Dartel (create Templates): Template (Iteration 6)', ...
    substruct('.','val', '{}',{iDar}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','template', '()',{7}));
matlabbatch{iNS}.spm.tools.dartel.mni_norm.data.subjs.flowfields(1) = ...
    cfg_dep('Run Dartel (create Templates): Flow Fields', ...
    substruct('.','val', '{}',{iDar}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','files', '()',{':'}));
matlabbatch{iNS}.spm.tools.dartel.mni_norm.data.subjs.images{1}(1) = ...
    cfg_dep('Segment: c1 Images', ...
    substruct('.','val', '{}',{iSeg}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','tiss', '()',{1}, '.','c', '()',{':'}));
matlabbatch{iNS}.spm.tools.dartel.mni_norm.preserve = 1;
matlabbatch{iNS}.spm.tools.dartel.mni_norm.fwhm = [8 8 8];

%% Average the smwc1 images
iAS = length(matlabbatch) + 1;
matlabbatch{iAS}.spm.util.imcalc.input(1) = ...
    cfg_dep('Normalise to MNI Space: MNI Smo. Warped - Amount (Image 1)', ...
    substruct('.','val', '{}',{iNS}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('()',{':', 1}));
matlabbatch{iAS}.spm.util.imcalc.output = 'average_smwc1';
matlabbatch{iAS}.spm.util.imcalc.outdir(1) = ...
    cfg_dep('Make Directory: Make Directory ''Templates_etc''', ...
    substruct('.','val', '{}',{iDir}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','dir'));
matlabbatch{iAS}.spm.util.imcalc.expression = 'median(X)';
matlabbatch{iAS}.spm.util.imcalc.options.dmtx = 1;

%% Threshold the tissue average to get an explicit mask
iEM = length(matlabbatch) + 1;
matlabbatch{iEM}.spm.util.imcalc.input(1) = ...
    cfg_dep('Image Calculator: Imcalc Computed Image: average_smwc1', ...
    substruct('.','val', '{}',{iAS}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','files'));
matlabbatch{iEM}.spm.util.imcalc.output = 'mask_smwc1';
matlabbatch{iEM}.spm.util.imcalc.outdir(1) = cfg_dep('Make Directory: Make Directory ''Templates_etc''', substruct('.','val', '{}',{iDir}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','dir'));
matlabbatch{iEM}.spm.util.imcalc.expression = 'i1>0.2';
matlabbatch{iEM}.spm.util.imcalc.options.dtype = 2;

%% Normalise the bias corrected images to MNI space
iNI = length(matlabbatch) + 1;
matlabbatch{iNI}.spm.tools.dartel.mni_norm.template(1) = ...
    cfg_dep('Run Dartel (create Templates): Template (Iteration 6)', ...
    substruct('.','val', '{}',{iDar}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','template', '()',{7}));
matlabbatch{iNI}.spm.tools.dartel.mni_norm.data.subjs.flowfields(1) = ...
    cfg_dep('Run Dartel (create Templates): Flow Fields', ...
    substruct('.','val', '{}',{iDar}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','files', '()',{':'}));
matlabbatch{iNI}.spm.tools.dartel.mni_norm.data.subjs.images{1}(1) = ...
    cfg_dep('Segment: Bias Corrected (1)', ...
    substruct('.','val', '{}',{iSeg}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','channel', '()',{1}, '.','biascorr', '()',{':'}));
matlabbatch{iNI}.spm.tools.dartel.mni_norm.preserve = 0;
matlabbatch{iNI}.spm.tools.dartel.mni_norm.fwhm = [0 0 0];

%% Average the normalised bias corrected images
iAI = length(matlabbatch) + 1;
matlabbatch{iAI}.spm.util.imcalc.input(1) = ...
    cfg_dep('Normalise to MNI Space: MNI Smo. Warped - Concentrations (Image 1)', ...
    substruct('.','val', '{}',{7}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('()',{':', 1}));
matlabbatch{iAI}.spm.util.imcalc.output = 'average_image';
matlabbatch{iAI}.spm.util.imcalc.outdir(1) = ...
    cfg_dep('Make Directory: Make Directory ''Templates_etc''', ...
    substruct('.','val', '{}',{iDir}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','dir'));
matlabbatch{iAI}.spm.util.imcalc.expression = 'median(X)';
matlabbatch{iAI}.spm.util.imcalc.options.dmtx = 1;

%% Show the average image and explicit mask in Check Reg
iCR = length(matlabbatch) + 1;
matlabbatch{iCR}.spm.util.checkreg.data(1) = ...
    cfg_dep('Image Calculator: Imcalc Computed Image: average_image', ...
    substruct('.','val', '{}',{iAI}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','files'));
matlabbatch{iCR}.spm.util.checkreg.data(2) = ...
    cfg_dep('Image Calculator: Imcalc Computed Image: mask_smwc1', ...
    substruct('.','val', '{}',{iEM}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','files'));

%% Compute tissue volumes
iTV = length(matlabbatch) + 1;
matlabbatch{iTV}.spm.util.tvol.matfiles(1) = ...
    cfg_dep('Segment: Seg Params', ...
    substruct('.','val', '{}',{iSeg}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','param', '()',{':'}));
matlabbatch{iTV}.spm.util.tvol.outf = '<UNDEFINED>';

%% Move the Dartel templates to the Templates_etc directory
% BUG: The mat file for the affine reg of Template_6 to MNI is not moved.
% This is not too serious, since another will be created if norm-to-MNI is
% repeated later on the relocated template. Probably best to add a new
% option for a "Template directory" to the Dartel create template module.
iMov = length(matlabbatch) + 1;
for it = 0:6
    matlabbatch{iMov}.cfg_basicio.file_dir.file_ops.file_move.files(it+1) = ...
        cfg_dep(sprintf('Run Dartel (create Templates): Template (Iteration %d)', it), ...
        substruct('.','val', '{}',{iDar}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
        substruct('.','template', '()',{it+1}));
end
matlabbatch{iMov}.cfg_basicio.file_dir.file_ops.file_move.action.moveto(1) = ...
    cfg_dep('Make Directory: Make Directory ''Templates_etc''', ...
    substruct('.','val', '{}',{iDir}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','dir'));

%% Delete (most) unnecessary files
iDel = length(matlabbatch) + 1;
for c = 1:3
    matlabbatch{iDel}.cfg_basicio.file_dir.file_ops.file_move.files(c) = ...
        cfg_dep(sprintf('Segment: mwc%d Images', c), ...
        substruct('.','val', '{}',{iSeg}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
        substruct('.','tiss', '()',{c}, '.','mwc', '()',{':'}));
end
matlabbatch{iDel}.cfg_basicio.file_dir.file_ops.file_move.files(4) = ...
    cfg_dep('Segment: Bias Corrected (1)', ...
    substruct('.','val', '{}',{iSeg}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('.','channel', '()',{1}, '.','biascorr', '()',{':'}));
matlabbatch{iDel}.cfg_basicio.file_dir.file_ops.file_move.files(5) = ...
    cfg_dep('Normalise to MNI Space: MNI Smo. Warped - Concentrations (Image 1)', ...
    substruct('.','val', '{}',{iNI}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), ...
    substruct('()',{':', 1}));
matlabbatch{iDel}.cfg_basicio.file_dir.file_ops.file_move.action.delete = false;

%% If run as script, open matlabbatch GUI
% NB mfilename is cfg_load_jobs if loaded from matlabbatch GUI
if strcmp(mfilename, 'preproc_vbm')
    spm_jobman('interactive', matlabbatch)
end