ENH: Add dedicated session filtering, alternative anatomical template options

docs/faq.rst

@@ -290,8 +290,8 @@ Some examples follow:
 
 * Surgery: use only pre-operation sessions for the anatomical data. This will typically be done
   by omitting post-operation sessions from the inputs to *fMRIPrep*.
-* Developing and elderly populations: there is currently no standard way of processing these.
-  However, `as suggested by U. Tooley at NeuroStars.org
+* Developing and elderly populations: ``--subject-anatomical-reference sessionwise`` will process each session independently.
+  Additionally, `as suggested by U. Tooley at NeuroStars.org
   <https://neurostars.org/t/fmriprep-how-to-reuse-longitudinal-and-pre-run-freesurfer/4585/15>`__,
   it is theoretically possible to leverage the *anatomical fast-track* along with the
   ``--bids-filters`` option to process sessions fully independently, or grouped by some study-design

docs/workflows.rst

@@ -174,15 +174,22 @@ Your ``.bidsignore`` file should include the following line::
 
 Longitudinal processing
 ~~~~~~~~~~~~~~~~~~~~~~~
-In the case of multiple T1w images (across sessions and/or runs), T1w images are
-merged into a single template image using FreeSurfer's `mri_robust_template`_.
-This template may be *unbiased*, or equidistant from all source images, or
-aligned to the first image (determined lexicographically by session label).
+In the case of multiple T1w images (across sessions and/or runs), *fMRIPrep* provides
+a few choices on how to generate the reference anatomical space.
+
+If ``--subject-anatomical-reference first-lex`` is used, all T1w images are
+merged into a single template image using FreeSurfer's `mri_robust_template`_,
+aligned to the first image (determined lexicographically by session label). This is
+the default behavior.
+
+If ``--subject-anatomical-reference unbiased`` is used, all T1w images are merged into
+an *unbiased* template, equidistant from all source images.
 For two images, the additional cost of estimating an unbiased template is trivial,
 but aligning three or more images is too expensive to justify being the default behavior.
-For consistency, in the case of multiple images, *fMRIPrep* constructs
-templates aligned to the first image, unless passed the ``--longitudinal``
-flag, which forces the estimation of an unbiased template.
+
+If ``--subject-anatomical-reference sessionwise`` is used, a reference template will be
+generated for each session independently. If multiple T1w images are found within a session,
+the images will be aligned to the first lexicographical image from that session.
 
 .. note::
 

fmriprep/cli/parser.py

@@ -23,10 +23,14 @@
 """Parser."""
 
 import sys
+import typing as ty
 from pathlib import Path
 
 from .. import config
 
+if ty.TYPE_CHECKING:
+    from bids import BIDSLayout
+
 
 def _build_parser(**kwargs):
     """Build parser object.
@@ -51,6 +55,7 @@ def _build_parser(**kwargs):
         'force_bbr': ('--force bbr', '26.0.0'),
         'force_no_bbr': ('--force no-bbr', '26.0.0'),
         'force_syn': ('--force syn-sdc', '26.0.0'),
+        'longitudinal': ('--subject-anatomical-reference unbiased', '26.1.0'),
     }
 
     class DeprecatedAction(Action):
@@ -220,9 +225,14 @@ def _fallback_trt(value, parser):
         help='A space delimited list of participant identifiers or a single '
         'identifier (the sub- prefix can be removed)',
     )
-    # Re-enable when option is actually implemented
-    # g_bids.add_argument('-s', '--session-id', action='store', default='single_session',
-    #                     help='Select a specific session to be processed')
+
+    g_bids.add_argument(
+        '--session-label',
+        nargs='+',
+        type=lambda label: label.removeprefix('ses-'),
+        help='A space delimited list of session identifiers or a single '
+        'identifier (the ses- prefix can be removed)',
+    )
     # Re-enable when option is actually implemented
     # g_bids.add_argument('-r', '--run-id', action='store', default='single_run',
     #                     help='Select a specific run to be processed')
@@ -235,6 +245,17 @@ def _fallback_trt(value, parser):
         type=int,
         help='Select a specific echo to be processed in a multiecho series',
     )
+    g_bids.add_argument(
+        '--subject-anatomical-reference',
+        choices=['first-lex', 'unbiased', 'sessionwise'],
+        default='first-lex',
+        help='Method to produce the reference anatomical space:\n'
+        '\t"first-lex" will use the first image in lexicographical order\n'
+        '\t"unbiased" will construct an unbiased template from all images '
+        '(previously "--longitudinal")\n'
+        '\t"sessionwise" will independently process each session. If multiple runs are '
+        'found, the behavior will be similar to "first-lex" for each session.',
+    )
     g_bids.add_argument(
         '--bids-filter-file',
         dest='bids_filters',
@@ -384,8 +405,8 @@ def _fallback_trt(value, parser):
     )
     g_conf.add_argument(
         '--longitudinal',
-        action='store_true',
-        help='Treat dataset as longitudinal - may increase runtime',
+        action=DeprecatedAction,
+        help='Deprecated - use `--subject-anatomical-reference unbiased` instead',
     )
     g_conf.add_argument(
         '--bold2anat-init',
@@ -769,6 +790,7 @@ def parse_args(args=None, namespace=None):
     """Parse args and run further checks on the command line."""
     import logging
 
+    from niworkflows.utils.bids import collect_participants
     from niworkflows.utils.spaces import Reference, SpatialReferences
 
     parser = _build_parser()
@@ -779,6 +801,14 @@ def parse_args(args=None, namespace=None):
         config.load(opts.config_file, skip=skip, init=False)
         config.loggers.cli.info(f'Loaded previous configuration file {opts.config_file}')
 
+    if opts.longitudinal:
+        opts.subject_anatomical_reference = 'unbiased'
+        msg = (
+            'The `--longitudinal` flag is deprecated - use '
+            '`--subject-anatomical-reference unbiased` instead.'
+        )
+        config.loggers.cli.warning(msg)
+
     config.execution.log_level = int(max(25 - 5 * opts.verbose_count, logging.DEBUG))
     config.from_dict(vars(opts), init=['nipype'])
 
@@ -927,19 +957,68 @@ def parse_args(args=None, namespace=None):
     work_dir.mkdir(exist_ok=True, parents=True)
 
     # Force initialization of the BIDSLayout
+    config.loggers.cli.debug('Initializing BIDS Layout')
     config.execution.init()
-    all_subjects = config.execution.layout.get_subjects()
+
+    # Please note this is the input folder's dataset_description.json
+    dset_desc_path = config.execution.bids_dir / 'dataset_description.json'
+    if dset_desc_path.exists():
+        from hashlib import sha256
+
+        desc_content = dset_desc_path.read_bytes()
+        config.execution.bids_description_hash = sha256(desc_content).hexdigest()
+
+    # First check that bids_dir looks like a BIDS folder
+    subject_list = collect_participants(
+        config.execution.layout, participant_label=config.execution.participant_label
+    )
     if config.execution.participant_label is None:
-        config.execution.participant_label = all_subjects
+        config.execution.participant_label = subject_list
 
-    participant_label = set(config.execution.participant_label)
-    missing_subjects = participant_label - set(all_subjects)
-    if missing_subjects:
-        parser.error(
-            'One or more participant labels were not found in the BIDS directory: {}.'.format(
-                ', '.join(missing_subjects)
+    session_list = config.execution.session_label or []
+    subject_session_list = create_processing_groups(
+        config.execution.layout,
+        subject_list,
+        session_list,
+        config.workflow.subject_anatomical_reference,
+    )
+    config.execution.processing_groups = subject_session_list
+    config.execution.participant_label = sorted(subject_list)
+    config.workflow.skull_strip_template = config.workflow.skull_strip_template[0]
+
+
+def create_processing_groups(
+    layout: 'BIDSLayout',
+    subject_list: list,
+    session_list: list | str | None,
+    subject_anatomical_reference: str,
+) -> list[tuple[str]]:
+    """Generate a list of subject-session pairs to be processed."""
+    from bids.layout import Query
+
+    subject_session_list = []
+
+    for subject in subject_list:
+        sessions = (
+            layout.get_sessions(
+                scope='raw',
+                subject=subject,
+                session=session_list or Query.OPTIONAL,
             )
+            or None
         )
 
-    config.execution.participant_label = sorted(participant_label)
-    config.workflow.skull_strip_template = config.workflow.skull_strip_template[0]
+        if subject_anatomical_reference == 'sessionwise':
+            if not sessions:
+                raise RuntimeError(
+                    '`--subject-anatomical-reference sessionwise` was requested, but no sessions '
+                    f'found for subject {subject}.'
+                )
+            for session in sessions:
+                if len(session) == 1:
+                    session = session[0]
+                subject_session_list.append((subject, session))
+        else:
+            subject_session_list.append((subject, sessions))
+
+    return subject_session_list

fmriprep/cli/workflow.py

@@ -35,15 +35,13 @@
 def build_workflow(config_file, retval):
     """Create the Nipype Workflow that supports the whole execution graph."""
 
-    from niworkflows.utils.bids import collect_participants
     from niworkflows.utils.misc import check_valid_fs_license
 
+    from fmriprep import config, data
     from fmriprep.reports.core import generate_reports
     from fmriprep.utils.bids import check_pipeline_version
-
-    from .. import config, data
-    from ..utils.misc import check_deps
-    from ..workflows.base import init_fmriprep_wf
+    from fmriprep.utils.misc import check_deps, fmt_subjects_sessions
+    from fmriprep.workflows.base import init_fmriprep_wf
 
     config.load(config_file)
     build_log = config.loggers.workflow
@@ -69,27 +67,20 @@ def build_workflow(config_file, retval):
     if msg is not None:
         build_log.warning(msg)
 
-    # Please note this is the input folder's dataset_description.json
-    dset_desc_path = config.execution.bids_dir / 'dataset_description.json'
-    if dset_desc_path.exists():
-        from hashlib import sha256
-
-        desc_content = dset_desc_path.read_bytes()
-        config.execution.bids_description_hash = sha256(desc_content).hexdigest()
-
-    # First check that bids_dir looks like a BIDS folder
-    subject_list = collect_participants(
-        config.execution.layout, participant_label=config.execution.participant_label
-    )
-
     # Called with reports only
     if config.execution.reports_only:
-        build_log.log(25, 'Running --reports-only on participants %s', ', '.join(subject_list))
-        session_list = (
-            config.execution.bids_filters.get('bold', {}).get('session')
-            if config.execution.bids_filters
-            else None
+        build_log.log(
+            25,
+            'Running --reports-only on %s',
+            fmt_subjects_sessions(config.execution.processing_groups),
         )
+        session_list = config.execution.session_label
+        if not session_list:
+            session_list = (
+                config.execution.bids_filters.get('bold', {}).get('session')
+                if config.execution.bids_filters
+                else None
+            )
 
         failed_reports = generate_reports(
             config.execution.participant_label,
@@ -110,7 +101,7 @@ def build_workflow(config_file, retval):
     init_msg = [
         "Building fMRIPrep's workflow:",
         f'BIDS dataset path: {config.execution.bids_dir}.',
-        f'Participant list: {subject_list}.',
+        f'Participants and sessions: {fmt_subjects_sessions(config.execution.processing_groups)}.',
         f'Run identifier: {config.execution.run_uuid}.',
         f'Output spaces: {config.execution.output_spaces}.',
     ]