[ENH] add preproc workflow for anat with several sessions

demos/openneuro/Makefile

@@ -49,7 +49,7 @@ data_ds000224:
 data_ds001168:
 	mkdir -p inputs
 	cd inputs && datalad install ///openneuro/ds001168
-	cd inputs/ds001168 && datalad get ds001168/sub-0[12] -J 2
+	cd inputs/ds001168 && datalad get ds001168/sub-0[12] -J 12
 
 data_ds001734:
 	mkdir -p inputs
@@ -66,5 +66,13 @@ data_ds002799:
 	cd inputs && datalad install ///openneuro/ds002799
 	datalad get -d inputs/ds002799 inputs/ds002799/derivatives/fmriprep/sub-292/*/func/*MNI152NLin2009cAsym*
 	datalad get -d inputs/ds002799 inputs/ds002799/derivatives/fmriprep/sub-292/*/func/*tsv
-	datalad get -d inputs/ds002799 inputs/ds002799/derivatives/fmriprep/sub-30[27]/*/func/*MNI152NLin2009cAsym* -J 2
-	datalad get -d inputs/ds002799 inputs/ds002799/derivatives/fmriprep/sub-30[27]/*/func/*tsv -J 2
+	datalad get -d inputs/ds002799 inputs/ds002799/derivatives/fmriprep/sub-30[27]/*/func/*MNI152NLin2009cAsym* -J 12
+	datalad get -d inputs/ds002799 inputs/ds002799/derivatives/fmriprep/sub-30[27]/*/func/*tsv -J 12
+
+
+data_ds000201:
+	mkdir -p inputs
+	cd inputs && datalad install ///openneuro/ds000201
+	datalad get -d inputs/ds000201 inputs/ds000201/sub-900[1-5]/*/anat/*T1w* -J 12
+# datalad get -d inputs/ds000201 inputs/ds000201/sub-900[1-5]/*/func -J 12
+# datalad get -d inputs/ds000201 inputs/ds000201/sub-900[1-5]/*/fmap -J 12

demos/openneuro/ds000201_run.m

@@ -0,0 +1,26 @@
+% (C) Copyright 2024 bidspm developers
+
+% Example of anatomical preprocessing with one anat per session
+
+clear;
+clc;
+
+addpath(fullfile(pwd, '..', '..'));
+bidspm();
+
+% The directory where the data are located
+root_dir = fileparts(mfilename('fullpath'));
+bids_dir = fullfile(root_dir, 'inputs', 'ds000201');
+output_dir = fullfile(root_dir, 'outputs', 'ds000201', 'derivatives');
+
+participant_label = {'9001', '9002', '9003', '9004', '9005'};
+
+%% Preprocessing
+bidspm(bids_dir, output_dir, 'subject', ...
+       'participant_label', participant_label, ...
+       'action', 'preprocess', ...
+       'anat_only', true, ...
+       'space', {'individual'}, ...
+       'skip_validation', true, ...
+       'dry_run', true, ...
+       'verbosity', 3);

src/batches/preproc/setBatchSegmentation.m

@@ -19,11 +19,22 @@
 
   % (C) Copyright 2020 bidspm developers
 
+  % TODO
+  % implement the opt.segment.do and opt.segment.force
+  % with segmentationAlreadyDone ?
+
+  if nargin < 3
+    imageToSegment = '';
+  end
+  if iscell(imageToSegment)
+    imageToSegment = char(imageToSegment);
+  end
+
   if ~opt.segment.do
-    opt.orderBatches.segment = 0;
+    opt.orderBatches.segment = [0];
     return
   else
-    opt.orderBatches.segment = numel(matlabbatch) + 1;
+    opt.orderBatches.segment(end + 1) = numel(matlabbatch) + 1;
   end
 
   printBatchName('Segmentation anatomical image', opt);
@@ -35,29 +46,23 @@
   % save bias corrected image = true
   preproc.channel.write = [false true];
 
-  % first part assumes we are in the bidsSpatialPreproc workflow
-  if isfield(opt, 'orderBatches') && isfield(opt.orderBatches, 'selectAnat')
+  % assumes vanilla bidsSpatialPreproc workflow
+  % - not anat only
+  if ~strcmp(imageToSegment, '')
+    for iImg = 1:size(imageToSegment, 1)
+      file = validationInputFile([], deblank(imageToSegment(iImg, :)));
+      preproc.channel.vols{iImg, 1} = file;
+    end
+
+  elseif isfield(opt, 'orderBatches') && isfield(opt.orderBatches, 'selectAnat')
 
-    % SAVE BIAS CORRECTED IMAGE
     preproc.channel.vols(1) = cfg_dep('Named File Selector: Anatomical(1) - Files', ...
                                       returnDependency(opt, 'selectAnat'), ...
                                       substruct('.', 'files', '{}', {1}));
-  else
 
-    % TODO
-    % implement the opt.segment.do and opt.segment.force
-    % with segmentationAlreadyDone ?
-
-    % in case a cell was given as input
-    if iscell(imageToSegment)
-      imageToSegment = char(imageToSegment);
-    end
-
-    % add all the images to segment
-    for iImg = 1:size(imageToSegment, 1)
-      file = validationInputFile([], deblank(imageToSegment(iImg, :)));
-      preproc.channel.vols{iImg, 1} = file;
-    end
+  else
+    % ???
+    return
 
   end
 

src/batches/preproc/setBatchSkullStripping.m

@@ -25,9 +25,10 @@
   % :rtype: structure
   %
   % This function will get its inputs from the segmentation batch by reading
-  % the dependency from ``opt.orderBatches.segment``. If this field is not specified it will
-  % try to get the results from the segmentation by relying on the ``anat``
-  % image returned by ``getAnatFilename``.
+  % the dependency from ``opt.orderBatches.segment``.
+  % If this field is not specified it will try
+  % to get the results from the segmentation by relying
+  % on the ``anat`` image returned by ``getAnatFilename``.
   %
   % The threshold for inclusion in the mask can be set by::
   %
@@ -67,8 +68,8 @@
     return
   end
 
-  % if this is part of a pipeline we get the segmentation dependency to get
-  % the input from.
+  % if this is part of a pipeline
+  % we get the segmentation dependency to get the input from.
   % Otherwise the files to process are stored in a cell
   if isfield(opt, 'orderBatches') && ...
       isfield(opt.orderBatches, 'segment') && ...

src/batches/setBatchSelectAnat.m

@@ -38,9 +38,17 @@
 
   printBatchName('selecting anatomical image', opt);
 
-  [anatImage, anatDataDir] = getAnatFilename(BIDS, opt, subLabel);
+  nbImgToReturn = 1;
+  if opt.anatOnly
+    nbImgToReturn = Inf;
+  end
+  [anatImage, anatDataDir] = getAnatFilename(BIDS, opt, subLabel, nbImgToReturn);
 
   matlabbatch{end + 1}.cfg_basicio.cfg_named_file.name = 'Anatomical';
-  matlabbatch{end}.cfg_basicio.cfg_named_file.files = { {fullfile(anatDataDir, anatImage)} };
+  if ischar(anatImage)
+    matlabbatch{end}.cfg_basicio.cfg_named_file.files = { {fullfile(anatDataDir, anatImage)} };
+  else
+    matlabbatch{end}.cfg_basicio.cfg_named_file.files = { fullfile(anatDataDir, anatImage) };
+  end
 
 end

src/bids/getAnatFilename.m

@@ -107,7 +107,7 @@
 
   anatImage = unzipAndReturnsFullpathName(anat);
 
-  msg = sprintf('  selecting anat file: %s', ...
+  msg = sprintf('  selecting anat file(s): %s', ...
                 bids.internal.create_unordered_list(bids.internal.format_path(anat)));
   logger('DEBUG', msg, 'options', opt, 'filename', mfilename());
 

src/cli/cliPreprocess.m

@@ -39,6 +39,7 @@ function cliPreprocess(varargin)
   end
 
   bidsReport(opt);
+  % TODO adapt boiler plate for anat only
   boilerplate(opt, ...
               'outputPath', fullfile(opt.dir.output, 'reports'), ...
               'pipelineType', 'preprocess', ...
@@ -49,22 +50,26 @@ function cliPreprocess(varargin)
 
   bidsCopyInputFolder(opt);
 
-  if opt.dummyScans > 0
-    tmpOpt = opt;
-    tmpOpt.dir.input = tmpOpt.dir.preproc;
-    bidsRemoveDummies(tmpOpt, ...
-                      'dummyScans', tmpOpt.dummyScans, ...
-                      'force', false);
-  end
+  if ~opt.anatOnly
 
-  bidsCheckVoxelSize(opt);
+    bidsCheckVoxelSize(opt);
 
-  if opt.useFieldmaps && ~opt.anatOnly
-    bidsCreateVDM(opt);
-  end
+    if opt.dummyScans > 0
+      tmpOpt = opt;
+      tmpOpt.dir.input = tmpOpt.dir.preproc;
+      bidsRemoveDummies(tmpOpt, ...
+                        'dummyScans', tmpOpt.dummyScans, ...
+                        'force', false);
+    end
+
+    if opt.useFieldmaps
+      bidsCreateVDM(opt);
+    end
+
+    if ~opt.stc.skip
+      bidsSTC(opt);
+    end
 
-  if ~opt.stc.skip && ~opt.anatOnly
-    bidsSTC(opt);
   end
 
   bidsSpatialPrepro(opt);

src/workflows/preproc/bidsSpatialPrepro.m

@@ -79,8 +79,22 @@
 
     printProcessingSubject(iSub, subLabel, opt);
 
+    %% Anat
     matlabbatch = setBatchSelectAnat(matlabbatch, BIDS, opt, subLabel);
+    anatFiles = matlabbatch{1}.cfg_basicio.cfg_named_file.files{1};
+
+    if opt.anatOnly && numel(anatFiles) > 1
+      printBatchName('coregister all anatomical data to first anatomical', opt);
+      for iAnat = 2:numel(anatFiles)
+        matlabbatch = setBatchCoregistration(matlabbatch, ...
+                                             opt, ...
+                                             anatFiles{1}, ...
+                                             anatFiles{iAnat});
+      end
+    end
 
+    %% Func
+    % dependency from SelectAnat
     if ~opt.realign.useUnwarp
       action = 'realign';
     else
@@ -92,36 +106,21 @@
                                                                subLabel, ...
                                                                action);
 
-    % dependency from file selector ('Anatomical')
     matlabbatch = setBatchCoregistrationFuncToAnat(matlabbatch, BIDS, opt, subLabel);
 
     matlabbatch = setBatchSaveCoregistrationMatrix(matlabbatch, BIDS, opt, subLabel);
 
-    anatFile = matlabbatch{1}.cfg_basicio.cfg_named_file.files{1}{1};
-
-    % TODO refactor with bidsSegmentSkullstrip
-    %% Skip segmentation / skullstripping if done previously
-    if skullstripDo && ...
-            ~opt.skullstrip.force && ...
-            skullstrippingAlreadyDone(anatFile, BIDS)
-      opt.skullstrip.do = false;
-    end
-    if segmentDo && ~opt.segment.force && ...
-            segmentationAlreadyDone(anatFile, BIDS)
-      opt.segment.do = false;
-      % but if we must force the skullstripping
-      % then we will need some segmentation input
-    elseif opt.skullstrip.do && ...
-            ~segmentationAlreadyDone(anatFile, BIDS)
-      opt.segment.do = true;
+    %% Segmentation / Skullstripping
+    for iAnat = 1:numel(anatFiles)
+      [matlabbatch, opt] = setBatchesSegmentationAndSkullstrip(matlabbatch, ...
+                                                               BIDS, ...
+                                                               opt, ...
+                                                               subLabel, ...
+                                                               segmentDo, ...
+                                                               skullstripDo, ...
+                                                               anatFiles{iAnat});
     end
 
-    [matlabbatch, opt] = setBatchSegmentation(matlabbatch, opt);
-
-    matlabbatch = setBatchSkullStripping(matlabbatch, BIDS, opt, subLabel);
-    opt.orderBatches.skullStripping = numel(matlabbatch) - 1;
-    opt.orderBatches.skullStrippingMask = numel(matlabbatch);
-
     %% Normalization
     if ismember('IXI549Space', opt.space)
       % dependency from segmentation
@@ -210,3 +209,49 @@
   createdFiles = bidsRename(opt);
 
 end
+
+function [matlabbatch, opt] = setBatchesSegmentationAndSkullstrip(varargin)
+
+  args = inputParser;
+
+  addRequired(args, 'matlabbatch', @iscell);
+  addRequired(args, 'BIDS', @isstruct);
+  addRequired(args, 'opt', @isstruct);
+  addRequired(args, 'subLabel', @ischar);
+  addRequired(args, 'segmentDo', @islogical);
+  addRequired(args, 'skullstripDo', @islogical);
+  addRequired(args, 'anatFile', @ischar);
+
+  parse(args, varargin{:});
+
+  matlabbatch = args.Results.matlabbatch;
+  BIDS = args.Results.BIDS;
+  opt = args.Results.opt;
+  subLabel = args.Results.subLabel;
+  segmentDo = args.Results.segmentDo;
+  skullstripDo = args.Results.skullstripDo;
+  anatFile = args.Results.anatFile;
+
+  % TODO refactor with bidsSegmentSkullstrip
+  %% Skip segmentation / skullstripping if done previously
+  if skullstripDo && ...
+          ~opt.skullstrip.force && ...
+          skullstrippingAlreadyDone(anatFile, BIDS)
+    opt.skullstrip.do = false;
+  end
+  if segmentDo && ~opt.segment.force && ...
+          segmentationAlreadyDone(anatFile, BIDS)
+    opt.segment.do = false;
+    % but if we must force the skullstripping
+    % then we will need some segmentation input
+  elseif opt.skullstrip.do && ...
+          ~segmentationAlreadyDone(anatFile, BIDS)
+    opt.segment.do = true;
+  end
+
+  [matlabbatch, opt] = setBatchSegmentation(matlabbatch, opt, anatFile);
+
+  matlabbatch = setBatchSkullStripping(matlabbatch, BIDS, opt, subLabel);
+  opt.orderBatches.skullStripping = numel(matlabbatch) - 1;
+  opt.orderBatches.skullStrippingMask = numel(matlabbatch);
+end

tests/create_dummy_dataset.py

@@ -42,9 +42,9 @@ def create_raw_dataset(target_dir, subject_list, session_list):
             create_raw_func_vismotion(target_dir, sub, ses)
             create_raw_func_vislocalizer(target_dir, sub, ses)
             create_raw_fmap(target_dir, sub, ses)
+            create_raw_anat(target_dir, sub, ses)
 
         create_raw_func_rest(target_dir, sub, "02")
-        create_raw_anat(target_dir, sub)
 
 
 def create_raw_func_vismotion(target_dir, sub, ses):
@@ -150,8 +150,7 @@ def create_raw_fmap(target_dir, sub, ses):
         json.dump(content, f, indent=4)
 
 
-def create_raw_anat(target_dir, sub):
-    ses = "01"
+def create_raw_anat(target_dir, sub, ses="01"):
     suffix = "T1w"
     this_dir = target_dir / f"sub-{sub}" / f"ses-{ses}" / "anat"