diff --git a/sdcflows/interfaces/epi.py b/sdcflows/interfaces/epi.py index 3d3bc7570a..b1931a9327 100644 --- a/sdcflows/interfaces/epi.py +++ b/sdcflows/interfaces/epi.py @@ -37,6 +37,10 @@ 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, usedefault=True, desc='Use "Estimated*" fields to calculate TotalReadoutTime' + ) + fallback = traits.Float(desc="A fallback value, in seconds.") class _GetReadoutTimeOutputSpec(TraitedSpec): @@ -57,6 +61,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..449c94b5db 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', use_estimate=True):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,42 +237,54 @@ 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 + + if "PhaseEncodingDirection" in in_meta: + # npe = N voxels PE direction + pe_index = "ijk".index(in_meta["PhaseEncodingDirection"][0]) + npe = nb.load(in_file).shape[pe_index] - # npe = N voxels PE direction - pe_index = "ijk".index(in_meta["PhaseEncodingDirection"][0]) - npe = nb.load(in_file).shape[pe_index] - - # Use case 2: EES is defined - ees = in_meta.get("EffectiveEchoSpacing") - if ees: - # Effective echo spacing means that acceleration factors have been accounted for. - return ees * (npe - 1) - - try: - echospacing = in_meta["EchoSpacing"] - acc_factor = in_meta["ParallelReductionFactorInPlane"] - except KeyError: - pass - else: - # etl = effective train length - etl = npe // acc_factor - return echospacing * (etl - 1) - - # Use case 3 (Philips scans) - try: - wfs = in_meta["WaterFatShift"] - epifactor = in_meta["EPIFactor"] - except KeyError: - pass - else: - wfs_hz = ( - (in_meta.get("ImagingFrequency", 0) * 3.39941) - or (in_meta.get("MagneticFieldStrength", 0) * 144.7383333) - or None - ) - if wfs_hz: - ees = wfs / (wfs_hz * (epifactor + 1)) + # Use case 2: EES is defined + ees = in_meta.get("EffectiveEchoSpacing") + 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"] + acc_factor = in_meta["ParallelReductionFactorInPlane"] + except KeyError: + pass + else: + # etl = effective train length + etl = npe // acc_factor + return echospacing * (etl - 1) + + # Use case 3 (Philips scans) + try: + wfs = in_meta["WaterFatShift"] + epifactor = in_meta["EPIFactor"] + except KeyError: + pass + else: + wfs_hz = ( + (in_meta.get("ImagingFrequency", 0) * 3.39941) + or (in_meta.get("MagneticFieldStrength", 0) * 144.7383333) + or None + ) + if wfs_hz: + ees = wfs / (wfs_hz * (epifactor + 1)) + return ees * (npe - 1) + + if fallback: + return fallback raise ValueError("Unknown total-readout time specification")