feat(get_trt): Allow estimated and fallback TotalReadoutTime

Author: effigies

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"] = (

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,6 +237,12 @@ 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 +253,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 +282,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")