feat: Add flag to fallback to Estimated* metadata or a passed value for TotalReadoutTime (#3423)

This flag provides two additional ways to retrieve TotalReadoutTime. If
`EstimatedTotalReadoutTime` or `EstimatedEffectiveEchoSpacing` are
defined, the user can pass `--fallback-total-readout-time estimated` to
direct fMRIPrep to use these.

This could be turned on by default, but I don't think this is prudent.
The reason that `dcm2niix` does not simply say `TotalReadoutTime` or
`EffectiveEchoSpacing` is not to misrepresent the confidence that these
are correct. Previously, users had to actively remove the `Estimated` to
accept the estimates, and I don't think we want to eliminate the
conscious decision. However, we can make it so that users do not need to
modify their datasets to accept them.

Closes #3009.

Author: effigies

.circleci/config.yml

@@ -486,7 +486,6 @@ jobs:
 
             # Inject pretend metadata
             json_sidecar=/tmp/data/${DATASET}/task-mixedgamblestask_bold.json
-            awk 'NR==1{print; print "    \"TotalReadoutTime\": 0.05,"} NR!=1' ${json_sidecar} > tmp && mv tmp ${json_sidecar}
             awk 'NR==1{print; print "    \"PhaseEncodingDirection\": \"j\","} NR!=1' ${json_sidecar} > tmp && mv tmp ${json_sidecar}
 
             fmriprep-docker -i nipreps/fmriprep:latest \
@@ -530,7 +529,8 @@ jobs:
                 /tmp/data/${DATASET} /tmp/${DATASET}/fmriprep-partial participant \
                 --fs-subjects-dir /tmp/${DATASET}/freesurfer \
                 ${FASTRACK_ARG} \
-                --sloppy --write-graph --use-syn-sdc --mem-mb 14336 \
+                --use-syn-sdc --fallback-total-readout-time 0.03125 \
+                --sloppy --write-graph --mem-mb 14336 \
                 --output-spaces MNI152NLin2009cAsym fsaverage5 fsnative MNI152NLin6Asym anat \
                 --nthreads 4 --cifti-output --project-goodvoxels -vv
       - store_artifacts:
@@ -745,7 +745,7 @@ jobs:
             # Inject pretend metadata for SDCFlows not to crash
             # TODO / open question - do all echos need the metadata?
             chmod +w /tmp/data/${DATASET}
-            echo '{"PhaseEncodingDirection": "j", "TotalReadoutTime": 0.058}' >> /tmp/data/${DATASET}/task-cuedSGT_bold.json
+            echo '{"PhaseEncodingDirection": "j"}' >> /tmp/data/${DATASET}/task-cuedSGT_bold.json
             chmod -R -w /tmp/data/${DATASET}
 
             fmriprep-docker -i nipreps/fmriprep:latest \
@@ -754,7 +754,8 @@ jobs:
                 /tmp/data/${DATASET} /tmp/${DATASET}/fmriprep participant \
                 ${FASTRACK_ARG} \
                 --me-output-echos \
-                --fs-no-reconall --use-syn-sdc --ignore slicetiming \
+                --fs-no-reconall --ignore slicetiming \
+                --use-syn-sdc --fallback-total-readout-time 0.0625 \
                 --dummy-scans 1 --sloppy --write-graph \
                 --output-spaces MNI152NLin2009cAsym \
                 --mem-mb 14336 --nthreads 4 -vv

fmriprep/cli/parser.py

@@ -149,6 +149,16 @@ def _slice_time_ref(value, parser):
             raise parser.error(f'Slice time reference must be in range 0-1. Received {value}.')
         return value
 
+    def _fallback_trt(value, parser):
+        if value == 'estimated':
+            return value
+        try:
+            return float(value)
+        except ValueError:
+            raise parser.error(
+                f'Falling back to TRT must be a number or "estimated". Received {value}.'
+            ) from None
+
     verstr = f'fMRIPrep v{config.environment.version}'
     currentv = Version(config.environment.version)
     is_release = not any((currentv.is_devrelease, currentv.is_prerelease, currentv.is_postrelease))
@@ -163,6 +173,7 @@ def _slice_time_ref(value, parser):
     PositiveInt = partial(_min_one, parser=parser)
     BIDSFilter = partial(_bids_filter, parser=parser)
     SliceTimeRef = partial(_slice_time_ref, parser=parser)
+    FallbackTRT = partial(_fallback_trt, parser=parser)
 
     # Arguments as specified by BIDS-Apps
     # required, positional arguments
@@ -417,6 +428,15 @@ def _slice_time_ref(value, parser):
         type=int,
         help='Number of nonsteady-state volumes. Overrides automatic detection.',
     )
+    g_conf.add_argument(
+        '--fallback-total-readout-time',
+        required=False,
+        action='store',
+        default=None,
+        type=FallbackTRT,
+        help='Fallback value for Total Readout Time (TRT) calculation. '
+        'May be a number or "estimated".',
+    )
     g_conf.add_argument(
         '--random-seed',
         dest='_random_seed',

fmriprep/config.py

@@ -575,6 +575,9 @@ class workflow(_Config):
     """Remove the mean from fieldmaps."""
     force_syn = None
     """Run *fieldmap-less* susceptibility-derived distortions estimation."""
+    fallback_total_readout_time = None
+    """Infer the total readout time if unavailable from authoritative metadata.
+    This may be a number or the string "estimated"."""
     hires = None
     """Run FreeSurfer ``recon-all`` with the ``-hires`` flag."""
     fs_no_resume = None

fmriprep/interfaces/resampling.py

@@ -30,7 +30,6 @@ class ResampleSeriesInputSpec(TraitedSpec):
     ref_file = File(exists=True, mandatory=True, desc='File to resample in_file to')
     transforms = InputMultiObject(
         File(exists=True),
-        mandatory=True,
         desc='Transform files, from in_file to ref_file (image mode)',
     )
     inverse = InputMultiObject(
@@ -92,7 +91,8 @@ def _run_interface(self, runtime):
 
         nvols = source.shape[3] if source.ndim > 3 else 1
 
-        transforms = load_transforms(self.inputs.transforms, self.inputs.inverse)
+        # No transforms appear Undefined, pass as empty list
+        transforms = load_transforms(self.inputs.transforms or [], self.inputs.inverse)
 
         pe_dir = self.inputs.pe_dir
         ro_time = self.inputs.ro_time
@@ -191,6 +191,13 @@ def _run_interface(self, runtime):
 class DistortionParametersInputSpec(TraitedSpec):
     in_file = File(exists=True, desc='EPI image corresponding to the metadata')
     metadata = traits.Dict(mandatory=True, desc='metadata corresponding to the inputs')
+    fallback = traits.Either(
+        None,
+        'estimated',
+        traits.Float,
+        usedefault=True,
+        desc='Fallback value for missing metadata',
+    )
 
 
 class DistortionParametersOutputSpec(TraitedSpec):
@@ -215,6 +222,8 @@ def _run_interface(self, runtime):
             self._results['readout_time'] = get_trt(
                 self.inputs.metadata,
                 self.inputs.in_file or None,
+                use_estimate=self.inputs.fallback == 'estimated',
+                fallback=self.inputs.fallback if isinstance(self.inputs.fallback, float) else None,
             )
             self._results['pe_direction'] = self.inputs.metadata['PhaseEncodingDirection']
         except (KeyError, ValueError):

fmriprep/workflows/base.py

@@ -614,13 +614,27 @@ def init_single_subject_wf(subject_id: str):
         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,
-        )
+        fallback_trt = config.workflow.fallback_total_readout_time
+        try:
+            fmap_wf = init_fmap_preproc_wf(
+                use_metadata_estimates=fallback_trt == 'estimated',
+                fallback_total_readout_time=fallback_trt
+                if isinstance(fallback_trt, float)
+                else None,
+                debug='fieldmaps' in config.execution.debug,
+                estimators=fmap_estimators,
+                omp_nthreads=omp_nthreads,
+                output_dir=fmriprep_dir,
+                subject=subject_id,
+            )
+        except RuntimeError:
+            message = (
+                'Missing readout time information. '
+                'See documentation for `--fallback-total-readout-time`.'
+            )
+            config.loggers.workflow.critical(message, exc_info=True)
+            sys.exit(os.EX_DATAERR)
+
         fmap_wf.__desc__ = f"""
 
 Preprocessing of B<sub>0</sub> inhomogeneity mappings

fmriprep/workflows/bold/apply.py

@@ -17,6 +17,7 @@ def init_bold_volumetric_resample_wf(
     metadata: dict,
     mem_gb: dict[str, float],
     jacobian: bool,
+    fallback_total_readout_time: str | float | None = None,
     fieldmap_id: str | None = None,
     omp_nthreads: int = 1,
     name: str = 'bold_volumetric_resample_wf',
@@ -161,7 +162,10 @@ def init_bold_volumetric_resample_wf(
         run_without_submitting=True,
     )
     distortion_params = pe.Node(
-        DistortionParameters(metadata=metadata),
+        DistortionParameters(
+            metadata=metadata,
+            fallback=fallback_total_readout_time,
+        ),
         name='distortion_params',
         run_without_submitting=True,
     )

fmriprep/workflows/bold/base.py

@@ -383,6 +383,7 @@ def init_bold_wf(
     # Resample to anatomical space
     bold_anat_wf = init_bold_volumetric_resample_wf(
         metadata=all_metadata[0],
+        fallback_total_readout_time=config.workflow.fallback_total_readout_time,
         fieldmap_id=fieldmap_id if not multiecho else None,
         omp_nthreads=omp_nthreads,
         mem_gb=mem_gb,

fmriprep/workflows/bold/fit.py

@@ -31,7 +31,6 @@
 from niworkflows.interfaces.header import ValidateImage
 from niworkflows.interfaces.nitransforms import ConcatenateXFMs
 from niworkflows.interfaces.utility import KeySelect
-from sdcflows.workflows.apply.correction import init_unwarp_wf
 from sdcflows.workflows.apply.registration import init_coeff2epi_wf
 
 from ... import config
@@ -189,7 +188,6 @@ def init_bold_fit_wf(
     * :py:func:`~fmriprep.workflows.bold.hmc.init_bold_hmc_wf`
     * :py:func:`~niworkflows.func.utils.init_enhance_and_skullstrip_bold_wf`
     * :py:func:`~sdcflows.workflows.apply.registration.init_coeff2epi_wf`
-    * :py:func:`~sdcflows.workflows.apply.correction.init_unwarp_wf`
     * :py:func:`~fmriprep.workflows.bold.registration.init_bold_reg_wf`
     * :py:func:`~fmriprep.workflows.bold.outputs.init_ds_boldref_wf`
     * :py:func:`~fmriprep.workflows.bold.outputs.init_ds_hmc_wf`
@@ -547,12 +545,26 @@ def init_bold_fit_wf(
                     (ds_fmapreg_wf, fmapreg_buffer, [('outputnode.xform', 'boldref2fmap_xfm')]),
                 ])  # fmt:skip
 
-            unwarp_wf = init_unwarp_wf(
-                free_mem=config.environment.free_mem,
-                debug='fieldmaps' in config.execution.debug,
-                omp_nthreads=config.nipype.omp_nthreads,
+            boldref_fmap = pe.Node(
+                ReconstructFieldmap(inverse=[True]), name='boldref_fmap', mem_gb=1
+            )
+
+            distortion_params = pe.Node(
+                DistortionParameters(
+                    metadata=metadata,
+                    in_file=bold_file,
+                    fallback=config.workflow.fallback_total_readout_time,
+                ),
+                name='distortion_params',
+                run_without_submitting=True,
+            )
+
+            unwarp_boldref = pe.Node(
+                ResampleSeries(jacobian=jacobian),
+                name='unwarp_boldref',
+                n_procs=omp_nthreads,
+                mem_gb=mem_gb['resampled'],
             )
-            unwarp_wf.inputs.inputnode.metadata = layout.get_metadata(bold_file)
 
             skullstrip_bold_wf = init_skullstrip_bold_wf()
 
@@ -564,24 +576,26 @@ def init_bold_fit_wf(
                     ('sdc_method', 'sdc_method'),
                     ('fmap_id', 'keys'),
                 ]),
-                (fmap_select, unwarp_wf, [
-                    ('fmap_coeff', 'inputnode.fmap_coeff'),
+                (fmapref_buffer, boldref_fmap, [('out', 'target_ref_file')]),
+                (fmapreg_buffer, boldref_fmap, [('boldref2fmap_xfm', 'transforms')]),
+                (fmap_select, boldref_fmap, [
+                    ('fmap_coeff', 'in_coeffs'),
+                    ('fmap_ref', 'fmap_ref_file'),
                 ]),
-                (fmapreg_buffer, unwarp_wf, [
-                    # This looks backwards, but unwarp_wf describes transforms in
-                    # terms of points while we (and init_coeff2epi_wf) describe them
-                    # in terms of images. Mapping fieldmap coordinates into boldref
-                    # coordinates maps the boldref image onto the fieldmap image.
-                    ('boldref2fmap_xfm', 'inputnode.fmap2data_xfm'),
+                (fmapref_buffer, unwarp_boldref, [('out', 'ref_file')]),
+                (enhance_boldref_wf, unwarp_boldref, [
+                    ('outputnode.bias_corrected_file', 'in_file'),
                 ]),
-                (enhance_boldref_wf, unwarp_wf, [
-                    ('outputnode.bias_corrected_file', 'inputnode.distorted'),
+                (boldref_fmap, unwarp_boldref, [('out_file', 'fieldmap')]),
+                (distortion_params, unwarp_boldref, [
+                    ('readout_time', 'ro_time'),
+                    ('pe_direction', 'pe_dir'),
                 ]),
-                (unwarp_wf, ds_coreg_boldref_wf, [
-                    ('outputnode.corrected', 'inputnode.boldref'),
+                (unwarp_boldref, ds_coreg_boldref_wf, [
+                    ('out_file', 'inputnode.boldref'),
                 ]),
-                (unwarp_wf, skullstrip_bold_wf, [
-                    ('outputnode.corrected', 'inputnode.in_file'),
+                (ds_coreg_boldref_wf, skullstrip_bold_wf, [
+                    ('outputnode.boldref', 'inputnode.in_file'),
                 ]),
                 (skullstrip_bold_wf, ds_boldmask_wf, [
                     ('outputnode.mask_file', 'inputnode.boldmask'),
@@ -591,7 +605,7 @@ def init_bold_fit_wf(
                 (fmapreg_buffer, func_fit_reports_wf, [
                     ('boldref2fmap_xfm', 'inputnode.boldref2fmap_xfm'),
                 ]),
-                (unwarp_wf, func_fit_reports_wf, [('outputnode.fieldmap', 'inputnode.fieldmap')]),
+                (boldref_fmap, func_fit_reports_wf, [('out_file', 'inputnode.fieldmap')]),
             ])  # fmt:skip
         else:
             workflow.connect([
@@ -859,7 +873,11 @@ def init_bold_native_wf(
         )
 
         distortion_params = pe.Node(
-            DistortionParameters(metadata=metadata, in_file=bold_file),
+            DistortionParameters(
+                metadata=metadata,
+                in_file=bold_file,
+                fallback=config.workflow.fallback_total_readout_time,
+            ),
             name='distortion_params',
             run_without_submitting=True,
         )

pyproject.toml

@@ -34,8 +34,8 @@ dependencies = [
     "psutil >= 5.4",
     "pybids >= 0.16",
     "requests >= 2.27",
-    "sdcflows >= 2.11.0",
-    "smriprep >= 0.17.0",
+    "sdcflows >= 2.13.0",
+    "smriprep >= 0.18.0",
     "tedana >= 23.0.2",
     "templateflow >= 24.2.2",
     "transforms3d >= 0.4",