diff --git a/sdcflows/interfaces/epi.py b/sdcflows/interfaces/epi.py index 3d3bc7570a..e67192b68c 100644 --- a/sdcflows/interfaces/epi.py +++ b/sdcflows/interfaces/epi.py @@ -37,6 +37,8 @@ class _GetReadoutTimeInputSpec(BaseInterfaceInputSpec): in_file = File(exists=True, desc="EPI image corresponding to the metadata") metadata = traits.Dict(mandatory=True, desc="metadata corresponding to the inputs") + use_estimate = traits.Bool(False, desc='Use "Estimated*" fields to calculate TotalReadoutTime') + fallback = traits.Float(desc="A fallback value, in seconds.") class _GetReadoutTimeOutputSpec(TraitedSpec): @@ -57,6 +59,8 @@ def _run_interface(self, runtime): self._results["readout_time"] = get_trt( self.inputs.metadata, self.inputs.in_file if isdefined(self.inputs.in_file) else None, + use_estimate=self.inputs.use_estimate, + fallback=self.inputs.fallback or None, ) self._results["pe_direction"] = self.inputs.metadata["PhaseEncodingDirection"] self._results["pe_dir_fsl"] = ( diff --git a/sdcflows/utils/epimanip.py b/sdcflows/utils/epimanip.py index a16e3a952d..875225283c 100644 --- a/sdcflows/utils/epimanip.py +++ b/sdcflows/utils/epimanip.py @@ -26,7 +26,13 @@ """ -def get_trt(in_meta, in_file=None): +def get_trt( + in_meta, + in_file=None, + *, + use_estimate=False, + fallback=None, +): r""" Obtain the *total readout time* :math:`T_\text{ro}` from available metadata. @@ -43,6 +49,26 @@ def get_trt(in_meta, in_file=None): >>> nb.Nifti1Image(np.zeros((90, 90, 60)), None, None).to_filename( ... tmpdir.join('epi.nii.gz').strpath) + Parameters + ---------- + in_meta: :class:`dict` + BIDS metadata dictionary. + in_file: :class:`str`, optional + Path to the EPI file. Used to determine the number of voxels along the + phase-encoding direction. + use_estimate: :class:`bool`, optional + Whether to use "Estimated*" fields to calculate the total readout time. + These are generated by dcm2niix when authoritative metadata is not available + but heuristic methods permit an estimation. + fallback: :class:`float`, optional + A fallback value, in seconds, to use when the total readout time cannot be + calculated. This should only be used in situations where the field is to be + determined from displacement fields, as in SyN-SDC. + A recommended "plausible" value would be 0.03125, to minimize the impact of + floating-point errors in the calculations. + + Examples + -------- >>> meta = {'TotalReadoutTime': 0.05251} >>> get_trt(meta) @@ -159,6 +185,23 @@ def get_trt(in_meta, in_file=None): Traceback (most recent call last): ValueError: + dcm2niix may provide "EstimatedTotalReadoutTime" or "EstimatedEffectiveEchoSpacing" + fields when converting Philips data. In order to use these fields, pass + ``use_estimate=True``: + + >>> get_trt({'EstimatedTotalReadoutTime': 0.05251}, use_estimate=True) + 0.05251 + >>> meta = {'EstimatedEffectiveEchoSpacing': 0.00059, + ... 'PhaseEncodingDirection': 'j-'} + >>> f"{get_trt(meta, in_file='epi.nii.gz'):g}" + '0.05251' + + Finally, if a fallback value is provided, it will be used when the total readout + time cannot be calculated by any method: + + >>> get_trt({}, fallback=0.03125) + 0.03125 + .. testcleanup:: >>> os.chdir(cwd) @@ -194,6 +237,13 @@ def get_trt(in_meta, in_file=None): raise ValueError(f"'{trt}'") return trt + elif use_estimate and "EstimatedTotalReadoutTime" in in_meta: + trt = in_meta.get("EstimatedTotalReadoutTime") + if not trt: + raise ValueError(f"'{trt}'") + + return trt + # npe = N voxels PE direction pe_index = "ijk".index(in_meta["PhaseEncodingDirection"][0]) @@ -204,6 +254,8 @@ def get_trt(in_meta, in_file=None): if ees: # Effective echo spacing means that acceleration factors have been accounted for. return ees * (npe - 1) + elif use_estimate and "EstimatedEffectiveEchoSpacing" in in_meta: + return in_meta.get("EstimatedEffectiveEchoSpacing") * (npe - 1) try: echospacing = in_meta["EchoSpacing"] @@ -231,6 +283,9 @@ def get_trt(in_meta, in_file=None): ees = wfs / (wfs_hz * (epifactor + 1)) return ees * (npe - 1) + if fallback: + return fallback + raise ValueError("Unknown total-readout time specification")