NEW: split fieldmap wf

fmriprep/workflows/fieldmap.py

@@ -0,0 +1,375 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2023 The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""
+fMRIPrep base processing workflows
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autofunction:: init_fmriprep_wf
+.. autofunction:: init_single_subject_wf
+
+"""
+
+import warnings
+
+import bids
+from nipype.interfaces import utility as niu
+from nipype.pipeline import engine as pe
+from niworkflows.utils.connections import listify
+
+from .. import config
+
+
+def init_single_subject_fieldmap_wf(subject_id: str):
+    """
+    Organize the preprocessing pipeline for a single subject.
+
+    It collects and reports information about the subject, and prepares
+    sub-workflows to perform anatomical and functional preprocessing.
+    Anatomical preprocessing is performed in a single workflow, regardless of
+    the number of sessions.
+    Functional preprocessing is performed using a separate workflow for each
+    individual BOLD series.
+
+    Parameters
+    ----------
+    subject_id : :obj:`str`
+        Subject label for this single-subject workflow.
+
+    Inputs
+    ------
+    subjects_dir : :obj:`str`
+        FreeSurfer's ``$SUBJECTS_DIR``.
+
+    """
+    from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+    from niworkflows.interfaces.nilearn import NILEARN_VERSION
+    from niworkflows.utils.bids import collect_data
+
+    from fmriprep.workflows.bold.base import init_bold_wf
+
+    workflow = Workflow(name=f'sub_{subject_id}_wf')
+    workflow.__desc__ = """
+Results included in this manuscript come from preprocessing
+performed using *fMRIPrep* {fmriprep_ver}
+(@fmriprep1; @fmriprep2; RRID:SCR_016216),
+which is based on *Nipype* {nipype_ver}
+(@nipype1; @nipype2; RRID:SCR_002502).
+
+""".format(
+        fmriprep_ver=config.environment.version, nipype_ver=config.environment.nipype_version
+    )
+    workflow.__postdesc__ = """
+
+Many internal operations of *fMRIPrep* use
+*Nilearn* {nilearn_ver} [@nilearn, RRID:SCR_001362],
+mostly within the functional processing workflow.
+For more details of the pipeline, see [the section corresponding
+to workflows in *fMRIPrep*'s documentation]\
+(https://fmriprep.readthedocs.io/en/latest/workflows.html \
+"FMRIPrep's documentation").
+
+
+### Copyright Waiver
+
+The above boilerplate text was automatically generated by fMRIPrep
+with the express intention that users should copy and paste this
+text into their manuscripts *unchanged*.
+It is released under the [CC0]\
+(https://creativecommons.org/publicdomain/zero/1.0/) license.
+
+### References
+
+""".format(
+        nilearn_ver=NILEARN_VERSION
+    )
+
+    subject_data = collect_data(
+        config.execution.layout,
+        subject_id,
+        task=config.execution.task_id,
+        echo=config.execution.echo_idx,
+        bids_filters=config.execution.bids_filters,
+    )[0]
+
+    if 'flair' in config.workflow.ignore:
+        subject_data['flair'] = []
+    if 't2w' in config.workflow.ignore:
+        subject_data['t2w'] = []
+
+    anat_only = config.workflow.anat_only
+    # Make sure we always go through these two checks
+    if not anat_only and not subject_data['bold']:
+        task_id = config.execution.task_id
+        raise RuntimeError(
+            "No BOLD images found for participant {} and task {}. "
+            "All workflows require BOLD images.".format(
+                subject_id, task_id if task_id else '<all>'
+            )
+        )
+
+    bold_runs = [
+        sorted(
+            listify(run),
+            key=lambda file: config.execution.layout.get_metadata(file).get('EchoTime', 0),
+        )
+        for run in subject_data['bold']
+    ]
+
+    if subject_data['roi']:
+        warnings.warn(
+            f"Lesion mask {subject_data['roi']} found. "
+            "Future versions of fMRIPrep will use alternative conventions. "
+            "Please refer to the documentation before upgrading.",
+            FutureWarning,
+        )
+
+    spaces = config.workflow.spaces
+
+    anatomical_cache = {}
+    if config.execution.derivatives:
+        from smriprep.utils.bids import collect_derivatives as collect_anat_derivatives
+
+        std_spaces = spaces.get_spaces(nonstandard=False, dim=(3,))
+        std_spaces.append("fsnative")
+        for deriv_dir in config.execution.derivatives:
+            anatomical_cache.update(
+                collect_anat_derivatives(
+                    derivatives_dir=deriv_dir,
+                    subject_id=subject_id,
+                    std_spaces=std_spaces,
+                )
+            )
+
+    inputnode_fields = [
+        't1w_preproc',
+        't1w_mask',
+        't1w_dseg',
+        't1w_tpms',
+        'subjects_dir',
+        'subject_id',
+        'fsnative2t1w_xfm',
+    ]
+    inputnode = pe.Node(niu.IdentityInterface(fields=inputnode_fields), name='inputnode')
+
+    # bids_root = str(config.execution.bids_dir)
+    fmriprep_dir = str(config.execution.fmriprep_dir)
+    omp_nthreads = config.nipype.omp_nthreads
+
+    fmap_estimators, estimator_map = map_fieldmap_estimation(
+        layout=config.execution.layout,
+        subject_id=subject_id,
+        bold_data=bold_runs,
+        ignore_fieldmaps="fieldmaps" in config.workflow.ignore,
+        use_syn=config.workflow.use_syn_sdc,
+        force_syn=config.workflow.force_syn,
+        filters=config.execution.get().get('bids_filters', {}).get('fmap'),
+    )
+
+    if fmap_estimators:
+        config.loggers.workflow.info(
+            "B0 field inhomogeneity map will be estimated with the following "
+            f"{len(fmap_estimators)} estimator(s): "
+            f"{[e.method for e in fmap_estimators]}."
+        )
+
+        from sdcflows import fieldmaps as fm
+        from sdcflows.workflows.base import init_fmap_preproc_wf
+
+        fmap_wf = init_fmap_preproc_wf(
+            debug="fieldmaps" in config.execution.debug,
+            estimators=fmap_estimators,
+            omp_nthreads=omp_nthreads,
+            output_dir=fmriprep_dir,
+            subject=subject_id,
+        )
+        fmap_wf.__desc__ = f"""
+
+Preprocessing of B<sub>0</sub> inhomogeneity mappings
+
+: A total of {len(fmap_estimators)} fieldmaps were found available within the input
+BIDS structure for this particular subject.
+"""
+
+        # Overwrite ``out_path_base`` of sdcflows's DataSinks
+        for node in fmap_wf.list_node_names():
+            if node.split(".")[-1].startswith("ds_"):
+                fmap_wf.get_node(node).interface.out_path_base = ""
+
+        for estimator in fmap_estimators:
+            config.loggers.workflow.info(
+                f"""\
+Setting-up fieldmap "{estimator.bids_id}" ({estimator.method}) with \
+<{', '.join(s.path.name for s in estimator.sources)}>"""
+            )
+
+            # Mapped and phasediff can be connected internally by SDCFlows
+            if estimator.method in (fm.EstimatorType.MAPPED, fm.EstimatorType.PHASEDIFF):
+                continue
+
+            suffices = [s.suffix for s in estimator.sources]
+
+            if estimator.method == fm.EstimatorType.PEPOLAR:
+                if len(suffices) == 2 and all(suf in ("epi", "bold", "sbref") for suf in suffices):
+                    wf_inputs = getattr(fmap_wf.inputs, f"in_{estimator.bids_id}")
+                    wf_inputs.in_data = [str(s.path) for s in estimator.sources]
+                    wf_inputs.metadata = [s.metadata for s in estimator.sources]
+                else:
+                    raise NotImplementedError("Sophisticated PEPOLAR schemes are unsupported.")
+
+    # Append the functional section to the existing anatomical excerpt
+    # That way we do not need to stream down the number of bold datasets
+    func_pre_desc = """
+Functional data preprocessing
+
+: For each of the {num_bold} BOLD runs found per subject (across all
+tasks and sessions), the following preprocessing was performed.
+""".format(
+        num_bold=len(bold_runs)
+    )
+
+    for bold_series in bold_runs:
+        bold_file = bold_series[0]
+        fieldmap_id = estimator_map.get(bold_file)
+
+        functional_cache = {}
+        if config.execution.derivatives:
+            from fmriprep.utils.bids import collect_derivatives, extract_entities
+
+            entities = extract_entities(bold_series)
+
+            for deriv_dir in config.execution.derivatives:
+                functional_cache.update(
+                    collect_derivatives(
+                        derivatives_dir=deriv_dir,
+                        entities=entities,
+                        fieldmap_id=fieldmap_id,
+                    )
+                )
+
+    return clean_datasinks(workflow)
+
+
+def map_fieldmap_estimation(
+    layout: bids.BIDSLayout,
+    subject_id: str,
+    bold_data: list[list[str]],
+    ignore_fieldmaps: bool,
+    use_syn: bool | str,
+    force_syn: bool,
+    filters: dict | None,
+) -> tuple[list, dict]:
+    if not any((not ignore_fieldmaps, use_syn, force_syn)):
+        return [], {}
+
+    from sdcflows import fieldmaps as fm
+    from sdcflows.utils.wrangler import find_estimators
+
+    # In the case where fieldmaps are ignored and `--use-syn-sdc` is requested,
+    # SDCFlows `find_estimators` still receives a full layout (which includes the fmap modality)
+    # and will not calculate fmapless schemes.
+    # Similarly, if fieldmaps are ignored and `--force-syn` is requested,
+    # `fmapless` should be set to True to ensure BOLD targets are found to be corrected.
+    fmap_estimators = find_estimators(
+        layout=layout,
+        subject=subject_id,
+        fmapless=bool(use_syn) or ignore_fieldmaps and force_syn,
+        force_fmapless=force_syn or ignore_fieldmaps and use_syn,
+        bids_filters=filters,
+    )
+
+    if not fmap_estimators:
+        if use_syn:
+            message = (
+                "Fieldmap-less (SyN) estimation was requested, but PhaseEncodingDirection "
+                "information appears to be absent."
+            )
+            config.loggers.workflow.error(message)
+            if use_syn == "error":
+                raise ValueError(message)
+        return [], {}
+
+    if ignore_fieldmaps and any(f.method == fm.EstimatorType.ANAT for f in fmap_estimators):
+        config.loggers.workflow.info(
+            'Option "--ignore fieldmaps" was set, but either "--use-syn-sdc" '
+            'or "--force-syn" were given, so fieldmap-less estimation will be executed.'
+        )
+        fmap_estimators = [f for f in fmap_estimators if f.method == fm.EstimatorType.ANAT]
+
+    # Pare down estimators to those that are actually used
+    # If fmap_estimators == [], all loops/comprehensions terminate immediately
+    all_ids = {fmap.bids_id for fmap in fmap_estimators}
+    bold_files = (bold_series[0] for bold_series in bold_data)
+
+    all_estimators = {
+        bold_file: [fmap_id for fmap_id in get_estimator(layout, bold_file) if fmap_id in all_ids]
+        for bold_file in bold_files
+    }
+
+    for bold_file, estimator_key in all_estimators.items():
+        if len(estimator_key) > 1:
+            config.loggers.workflow.warning(
+                f"Several fieldmaps <{', '.join(estimator_key)}> are "
+                f"'IntendedFor' <{bold_file}>, using {estimator_key[0]}"
+            )
+            estimator_key[1:] = []
+
+    # Final, 1-1 map, dropping uncorrected BOLD
+    estimator_map = {
+        bold_file: estimator_key[0]
+        for bold_file, estimator_key in all_estimators.items()
+        if estimator_key
+    }
+
+    fmap_estimators = [f for f in fmap_estimators if f.bids_id in estimator_map.values()]
+
+    return fmap_estimators, estimator_map
+
+
+def _prefix(subid):
+    return subid if subid.startswith('sub-') else f'sub-{subid}'
+
+
+def clean_datasinks(workflow: pe.Workflow) -> pe.Workflow:
+    # Overwrite ``out_path_base`` of smriprep's DataSinks
+    for node in workflow.list_node_names():
+        if node.split('.')[-1].startswith('ds_'):
+            workflow.get_node(node).interface.out_path_base = ""
+    return workflow
+
+
+def get_estimator(layout, fname):
+    field_source = layout.get_metadata(fname).get("B0FieldSource")
+    if isinstance(field_source, str):
+        field_source = (field_source,)
+
+    if field_source is None:
+        import re
+        from pathlib import Path
+
+        from sdcflows.fieldmaps import get_identifier
+
+        # Fallback to IntendedFor
+        intended_rel = re.sub(r"^sub-[a-zA-Z0-9]*/", "", str(Path(fname).relative_to(layout.root)))
+        field_source = get_identifier(intended_rel)
+
+    return field_source