Merge pull request #477 from effigies/enh/estimated-trt

feat(get_trt): Allow estimated and fallback TotalReadoutTime

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

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")