diff --git a/.circleci/config.yml b/.circleci/config.yml index 177b5fd4f..562fa7b14 100644 --- a/.circleci/config.yml +++ b/.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 diff --git a/fmriprep/cli/parser.py b/fmriprep/cli/parser.py index e29ec3867..d0fe44abd 100644 --- a/fmriprep/cli/parser.py +++ b/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', diff --git a/fmriprep/config.py b/fmriprep/config.py index 29221dc56..c9f6efd49 100644 --- a/fmriprep/config.py +++ b/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 diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py index 4f10d3260..e6162b0e2 100644 --- a/fmriprep/interfaces/resampling.py +++ b/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): diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index 8d174bb96..626cbd85c 100644 --- a/fmriprep/workflows/base.py +++ b/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 B0 inhomogeneity mappings diff --git a/fmriprep/workflows/bold/apply.py b/fmriprep/workflows/bold/apply.py index f019f13a0..408c4b289 100644 --- a/fmriprep/workflows/bold/apply.py +++ b/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, ) diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index 060db8b08..bcd282d50 100644 --- a/fmriprep/workflows/bold/base.py +++ b/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, diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 737bf9fca..3c74eccc2 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/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, ) diff --git a/pyproject.toml b/pyproject.toml index 89e1f8e4e..ef50325ad 100644 --- a/pyproject.toml +++ b/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",