[ENH] create an anat only spatial preprocessing

.github/workflows/miss_hit_quality.yml

@@ -29,8 +29,6 @@ jobs:
       run: |
         python -m pip install --upgrade pip setuptools
         pip3 install -r requirements.txt
-        cd tests
-        make data
 
     - name: MISS_HIT Metrics
       run: |

.github/workflows/miss_hit_style.yml

@@ -29,8 +29,6 @@ jobs:
       run: |
         python -m pip install --upgrade pip setuptools
         pip3 install -r requirements.txt
-        cd tests
-        make data
 
     - name: MISS_HIT Code style
       run: |

demos/face_repetition/face_rep_anat.m

@@ -6,7 +6,7 @@
 clear;
 clc;
 
-downloadData = true;
+downloadData = false;
 
 try
   run ../../initCppSpm.m;
@@ -20,6 +20,7 @@
 opt.pipeline.type = 'preproc';
 opt.pipeline.name = 'cpp_spm-anat';
 opt.query.modality = 'anat';
+opt.anatOnly = true;
 
 opt = checkOptions(opt);
 
@@ -32,10 +33,8 @@
 
 end
 
-%% Run batches
 reportBIDS(opt);
-bidsCopyInputFolder(opt);
 
-bidsSegmentSkullStrip(opt);
+bidsCopyInputFolder(opt);
 
-anatomicalQA(opt);
+bidsSpatialPrepro(opt);

src/QA/functionalQA.m

@@ -48,6 +48,10 @@ function functionalQA(opt)
   % TODO: - the tissue probability maps are in the "native" space of each subject
   %         and are resliced to the dimension of the functional
 
+  if opt.anatOnly
+    return
+  end
+
   if isOctave()
     warning('\nfunctionalQA is not yet supported on Octave. This step will be skipped.');
     return

src/batches/preproc/setBatchCoregistrationFuncToAnat.m

@@ -21,6 +21,10 @@
   %
   % (C) Copyright 2020 CPP_SPM developers
 
+  if opt.anatOnly
+    return
+  end
+
   printBatchName('coregister functional data to anatomical', opt);
 
   matlabbatch{end + 1}.spm.spatial.coreg.estimate.ref(1) = ...

src/batches/preproc/setBatchNormalizationSpatialPrepro.m

@@ -19,54 +19,60 @@
 
   jobsToAdd = numel(matlabbatch) + 1;
 
+  if opt.anatOnly
+    maxJobsToAdd = 4;
+  else
+    maxJobsToAdd = 5;
+  end
+
   % set the deformation field for all the images we are about to normalize
   % we reuse a previously created deformation field if it is there
   % otherwise we link to a segmentation module by using a dependency
   deformationField = getDeformationField(matlabbatch, BIDS, opt);
 
-  for iJob = jobsToAdd:(jobsToAdd + 5)
-
+  for iJob = jobsToAdd:(jobsToAdd + maxJobsToAdd)
     % we set images to be resampled at the voxel size we had at acquisition
     matlabbatch = setBatchNormalize(matlabbatch, deformationField, voxDim);
-
   end
 
-  printBatchName('normalise functional images', opt);
-
-  matlabbatch{jobsToAdd}.spm.spatial.normalise.write.subj.resample(1) = ...
-      cfg_dep('Coregister: Estimate: Coregistered Images', ...
-              returnDependency(opt, 'coregister'), ...
-              substruct('.', 'cfiles'));
+  if ~opt.anatOnly
+    printBatchName('normalise functional images', opt);
+    matlabbatch{jobsToAdd}.spm.spatial.normalise.write.subj.resample(1) = ...
+        cfg_dep('Coregister: Estimate: Coregistered Images', ...
+                returnDependency(opt, 'coregister'), ...
+                substruct('.', 'cfiles'));
+    jobsToAdd = jobsToAdd + 1;
+  end
 
   % NORMALIZE STRUCTURAL
   biasCorrectedImage = getBiasCorrectedImage(matlabbatch, BIDS, opt);
-
   printBatchName('normalise anatomical images', opt);
-  matlabbatch{jobsToAdd + 1}.spm.spatial.normalise.write.subj.resample(1) = biasCorrectedImage;
-
-  % size 3 allow to run RunQA / original voxel size at acquisition
-  matlabbatch{jobsToAdd + 1}.spm.spatial.normalise.write.woptions.vox = [1 1 1];
+  matlabbatch{jobsToAdd}.spm.spatial.normalise.write.subj.resample(1) = biasCorrectedImage;
+  % TODO why do we choose this resolution for this normalization?
+  matlabbatch{jobsToAdd}.spm.spatial.normalise.write.woptions.vox = [1 1 1];
+  jobsToAdd = jobsToAdd + 1;
 
   % NORMALIZE TISSUE PROBABILITY MAPS
   [gmTpm, wmTpm, csfTpm] = getTpms(matlabbatch, BIDS, opt);
 
   printBatchName('normalise grey matter tissue probability map', opt);
-  matlabbatch{jobsToAdd + 2}.spm.spatial.normalise.write.subj.resample(1) = gmTpm;
+  matlabbatch{jobsToAdd}.spm.spatial.normalise.write.subj.resample(1) = gmTpm;
+  jobsToAdd = jobsToAdd + 1;
 
   printBatchName('normalise white matter tissue probability map', opt);
-  matlabbatch{jobsToAdd + 3}.spm.spatial.normalise.write.subj.resample(1) = wmTpm;
+  matlabbatch{jobsToAdd}.spm.spatial.normalise.write.subj.resample(1) = wmTpm;
+  jobsToAdd = jobsToAdd + 1;
 
   printBatchName('normalise csf tissue probability map', opt);
-  matlabbatch{jobsToAdd + 4}.spm.spatial.normalise.write.subj.resample(1) = csfTpm;
+  matlabbatch{jobsToAdd}.spm.spatial.normalise.write.subj.resample(1) = csfTpm;
+  jobsToAdd = jobsToAdd + 1;
 
   % NORMALIZE SKULSTRIPPED STRUCTURAL
   skullstrippedImage = getSkullstrippedImage(matlabbatch, BIDS, opt);
-
   printBatchName('normalise skullstripped anatomical images', opt);
-  matlabbatch{jobsToAdd + 5}.spm.spatial.normalise.write.subj.resample(1) = skullstrippedImage;
-
-  % size 3 allow to run RunQA / original voxel size at acquisition
-  matlabbatch{jobsToAdd + 5}.spm.spatial.normalise.write.woptions.vox = [1 1 1];
+  matlabbatch{jobsToAdd}.spm.spatial.normalise.write.subj.resample(1) = skullstrippedImage;
+  % TODO why do we choose this resolution for this normalization?
+  matlabbatch{jobsToAdd}.spm.spatial.normalise.write.woptions.vox = [1 1 1];
 
 end
 

src/batches/preproc/setBatchRealign.m

@@ -5,21 +5,21 @@
   % USAGE::
   %
   %   [matlabbatch, voxDim] = setBatchRealign(matlabbatch, ...
-  %                                           [action = 'realign'], ...
   %                                           BIDS, ...
   %                                           opt, ...
-  %                                           subLabel)
+  %                                           subLabel, ...
+  %                                           [action = 'realign'])
   %
   % :param matlabbatch: SPM batch
-  % :type matlabbatch: structure
+  % :type matlabbatch: cell
   % :param BIDS: BIDS layout returned by ``getData``.
   % :type BIDS: structure
-  % :param action: ``realign``, ``realignReslice``, ``realignUnwarp``
-  % :type action: string
   % :param opt: Options chosen for the analysis. See ``checkOptions()``.
   % :type opt: structure
-  % :type subLabel: string
   % :param subLabel: subject label
+  % :type subLabel: string
+  % :param action: ``realign``, ``realignReslice``, ``realignUnwarp``, ``'reslice'``
+  % :type action: string
   %
   % :returns: - :matlabbatch: (structure) (dimension)
   %           - :voxDim: (array) (dimension)
@@ -29,15 +29,33 @@
 
   % TODO  make which image is resliced more consistent 'which = []'
 
-  if numel(varargin) < 5
-    [matlabbatch, BIDS, opt, subLabel] = deal(varargin{:});
-    action = '';
-  else
-    [matlabbatch, BIDS, opt, subLabel, action] = deal(varargin{:});
-  end
-
-  if isempty(action)
-    action = 'realignUnwarp';
+  p = inputParser;
+
+  default_action = 'realignUnwarp';
+  allowedActions = @(x) ischar(x) && ...
+                        ismember(lower(x), ...
+                                 {'realignunwarp', ...
+                                  'realignreslice', ...
+                                  'realign', ...
+                                  'reslice'});
+
+  addRequired(p, 'matlabbatch', @iscell);
+  addRequired(p, 'BIDS', @isstruct);
+  addRequired(p, 'opt', @isstruct);
+  addRequired(p, 'subLabel', @ischar);
+  addOptional(p, 'action', default_action, allowedActions);
+
+  parse(p, varargin{:});
+
+  matlabbatch = p.Results.matlabbatch;
+  BIDS = p.Results.BIDS;
+  opt = p.Results.opt;
+  subLabel = p.Results.subLabel;
+  action = p.Results.action;
+
+  if opt.anatOnly
+    voxDim = [];
+    return
   end
 
   msg = '';

src/batches/preproc/setBatchSaveCoregistrationMatrix.m

@@ -19,6 +19,10 @@
   %
   % (C) Copyright 2020 CPP_SPM developers
 
+  if opt.anatOnly
+    return
+  end
+
   printBatchName('saving coregistration matrix', opt);
 
   % create name of the output file based on the name

src/defaults/checkOptions.m

@@ -176,6 +176,8 @@
   fieldsToSet.pipeline.type = '';
   fieldsToSet.pipeline.name = 'cpp_spm';
 
+  fieldsToSet.anatOnly = false;
+
   fieldsToSet.useBidsSchema = false;
 
   fieldsToSet.fwhm.func = 6;

src/reports/copyFigures.m

@@ -19,7 +19,7 @@ function copyFigures(BIDS, opt, subLabel)
   imgNb = copyGraphWindownOutput(opt, subLabel, 'realign');
 
   % loop through the figures outputed for unwarp: one per run
-  if  opt.realign.useUnwarp
+  if opt.realign.useUnwarp && ~opt.anatOnly
 
     runs = bids.query(BIDS, 'runs', ...
                       'sub', subLabel, ...

src/reports/reportBIDS.m

@@ -25,9 +25,11 @@ function reportBIDS(opt)
 
     bids.util.mkdir(outputDir);
 
+    filter = struct('sub', subLabel, 'modality', {opt.query.modality});
+
     try
       bids.report(BIDS, ...
-                  'filter', struct('sub', subLabel), ...
+                  'filter', filter, ...
                   'output_path', outputDir, ...
                   'read_nifti', true, ...
                   'verbose', opt.verbosity);
@@ -39,7 +41,7 @@ function reportBIDS(opt)
       errorHandling(mfilename(), 'unspecifiedError', msg, true, opt.verbosity);
 
       bids.report(BIDS, ...
-                  'filter', struct('sub', subLabel), ...
+                  'filter', filter, ...
                   'output_path', outputDir, ...
                   'read_nifti', false, ...
                   'verbose', opt.verbosity);

src/utils/renameUnwarpParameter.m

@@ -8,6 +8,10 @@ function renameUnwarpParameter(BIDS, subLabel, opt)
   %
   % (C) Copyright 2020 CPP_SPM developers
 
+  if opt.anatOnly
+    return
+  end
+
   for iTask = 1:numel(opt.taskName)
 
     unwarpParam = spm_select('FPListRec', BIDS.pth, ...

src/workflows/preproc/bidsSpatialPrepro.m

@@ -46,9 +46,13 @@
   [BIDS, opt] = setUpWorkflow(opt, 'spatial preprocessing');
 
   opt.orderBatches.selectAnat = 1;
-  opt.orderBatches.realign = 2;
-  opt.orderBatches.coregister = 3;
-  opt.orderBatches.saveCoregistrationMatrix = 4;
+  lastBatch = 1;
+  if ~opt.anatOnly
+    opt.orderBatches.realign = 2;
+    opt.orderBatches.coregister = 3;
+    opt.orderBatches.saveCoregistrationMatrix = 4;
+    lastBatch = 4;
+  end
 
   for iSub = 1:numel(opt.subjects)
 
@@ -60,12 +64,11 @@
 
     matlabbatch = setBatchSelectAnat(matlabbatch, BIDS, opt, subLabel);
 
-    % if action is emtpy then only realign will be done
-    action = [];
     if ~opt.realign.useUnwarp
-      action = 'realign';
+      [matlabbatch, voxDim] = setBatchRealign(matlabbatch, BIDS, opt, subLabel, 'realign');
+    else
+      [matlabbatch, voxDim] = setBatchRealign(matlabbatch, BIDS, opt, subLabel);
     end
-    [matlabbatch, voxDim] = setBatchRealign(matlabbatch, BIDS, opt, subLabel, action);
 
     % dependency from file selector ('Anatomical')
     matlabbatch = setBatchCoregistrationFuncToAnat(matlabbatch, BIDS, opt, subLabel);
@@ -89,9 +92,9 @@
     doSegmentAndSkullstrip = opt.segment.force || isempty(tpm) || isempty(biasCorrectedImage);
 
     if doSegmentAndSkullstrip
-      opt.orderBatches.segment = 5;
-      opt.orderBatches.skullStripping = 6;
-      opt.orderBatches.skullStrippingMask = 7;
+      opt.orderBatches.segment = lastBatch + 1;
+      opt.orderBatches.skullStripping = lastBatch + 2;
+      opt.orderBatches.skullStrippingMask = lastBatch + 3;
       matlabbatch = setBatchSegmentation(matlabbatch, opt);
       matlabbatch = setBatchSkullStripping(matlabbatch, BIDS, opt, subLabel);
     end
@@ -145,7 +148,7 @@ function renameFile(BIDS, opt)
 
     subLabel = opt.subjects{iSub};
 
-    if ~opt.dryRun
+    if ~opt.dryRun && ~opt.anatOnly
 
       % this is only necessary if the rp_ files have not already been converted
       % and deleted by functionalQA
@@ -168,6 +171,8 @@ function renameFile(BIDS, opt)
 
   end
 
+  % TODO: when anatOnly update the res label for TPMs
+
   % TODO adapt spm_2_bids map to rename eventual files
   % that only have a "r" or "ra" prefix
   opt.query =  struct('modality', {{'anat', 'func'}});

tests/Makefile

@@ -10,7 +10,7 @@ clean:
 data:
 	sh createDummyDataSet.sh
 	rm -f bids-examples/.gitkeep
-	git clone git://github.com/bids-standard/bids-examples.git --depth 1
+	git clone https://github.com/bids-standard/bids-examples.git --depth 1
 	cp bids-examples/synthetic/dataset_description.json bids-examples/synthetic/derivatives/fmriprep
 
 clean_openneuro:

tests/tests_batches/test_setBatchCoregistrationFuncToAnat.m

@@ -8,6 +8,22 @@
   initTestSuite;
 end
 
+function test_setBatchCoregistrationFuncToAnat_anat_only()
+
+  subLabel = '^01';
+
+  opt = setOptions('vismotion', subLabel);
+  opt.anatOnly = true;
+
+  BIDS = cell(1);
+
+  matlabbatch = {};
+  matlabbatch = setBatchCoregistrationFuncToAnat(matlabbatch, BIDS, opt, subLabel);
+
+  assertEqual(matlabbatch, {});
+
+end
+
 function test_setBatchCoregistrationFuncToAnat_basic()
 
   % necessarry to deal with SPM module dependencies