diff --git a/fmriprep/workflows/fieldmap.py b/fmriprep/workflows/fieldmap.py new file mode 100644 index 000000000..50cd881d0 --- /dev/null +++ b/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 +# +# 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 '' + ) + ) + + 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 B0 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