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

feat: Add workflow arguments for metadata estimates and fallback TRT

Author: effigies

sdcflows/conftest.py

@@ -21,7 +21,9 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """py.test configuration."""
+
 import os
+import logging
 from pathlib import Path
 import numpy
 import nibabel
@@ -60,7 +62,7 @@ def pytest_report_header(config):
 
 
 @pytest.fixture(autouse=True)
-def doctest_fixture(doctest_namespace, request):
+def doctest_fixture(doctest_namespace, request, caplog):
     doctest_plugin = request.config.pluginmanager.getplugin("doctest")
     if isinstance(request.node, doctest_plugin.DoctestItem):
         doctest_namespace.update(
@@ -73,6 +75,8 @@ def doctest_fixture(doctest_namespace, request):
             dsB_dir=data_dir / "dsB",
             dsC_dir=data_dir / "dsC",
             data_dir=data_dir,
+            caplog=caplog,
+            logging=logging,
         )
         doctest_namespace.update((key, Path(val.root)) for key, val in layouts.items())
 

sdcflows/fieldmaps.py

@@ -24,15 +24,19 @@
 from pathlib import Path
 from enum import Enum, auto
 from collections import defaultdict
+from contextlib import suppress
 import re
 import attr
 from json import loads
 from bids.layout import parse_file_entities
 from bids.utils import listify
 from niworkflows.utils.bids import relative_to_root
 from .utils.bimap import EstimatorRegistry
+from .utils.misc import create_logger
 
 
+logger = create_logger('sdcflows.fieldmaps')
+
 _estimators = EstimatorRegistry()
 _intents = defaultdict(set)
 
@@ -109,12 +113,7 @@ class FieldmapFile:
     >>> f.suffix
     'T1w'
 
-    >>> FieldmapFile(
-    ...     dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz",
-    ...     find_meta=False
-    ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
-    Traceback (most recent call last):
-    MetadataError:
+    By default, the immediate JSON sidecar file is consulted for metadata.
 
     >>> f = FieldmapFile(
     ...     dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz",
@@ -123,11 +122,26 @@ class FieldmapFile:
     0.005
 
     >>> f = FieldmapFile(
-    ...     dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz",
-    ...     metadata={"TotalReadoutTime": None},
-    ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
-    Traceback (most recent call last):
-    MetadataError:
+    ...     dsA_dir / "sub-01" / "fmap" / "sub-01_phasediff.nii.gz"
+    ... )
+    >>> f.metadata['EchoTime2']
+    0.00746
+
+    >>> f = FieldmapFile(
+    ...     dsA_dir / "sub-01" / "fmap" / "sub-01_phase2.nii.gz"
+    ... )
+    >>> f.metadata['EchoTime']
+    0.00746
+
+    >>> f = FieldmapFile(
+    ...     dsA_dir / "sub-01" / "fmap" / "sub-01_fieldmap.nii.gz"
+    ... )
+    >>> f.metadata['Units']
+    'rad/s'
+
+    However, it is possible to provide alternative metadata.
+    It is recommended to load the metadata with an external tool that
+    fully implements the BIDS inheritance rules to ensure correctness.
 
     >>> f = FieldmapFile(
     ...     dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz",
@@ -136,44 +150,56 @@ class FieldmapFile:
     >>> f.metadata['TotalReadoutTime']
     0.006
 
+    If metadata is present, warnings or errors may be presented.
+
+    >>> FieldmapFile(
+    ...     dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz",
+    ...     find_meta=False
+    ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
+    Traceback (most recent call last):
+    MetadataError: Missing 'PhaseEncodingDirection' ...
+
     >>> FieldmapFile(
     ...     dsA_dir / "sub-01" / "fmap" / "sub-01_phasediff.nii.gz",
     ...     find_meta=False
     ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
     Traceback (most recent call last):
     MetadataError:
 
-    >>> f = FieldmapFile(
-    ...     dsA_dir / "sub-01" / "fmap" / "sub-01_phasediff.nii.gz"
-    ... )
-    >>> f.metadata['EchoTime2']
-    0.00746
-
     >>> FieldmapFile(
     ...     dsA_dir / "sub-01" / "fmap" / "sub-01_phase2.nii.gz",
     ...     find_meta=False
     ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
     Traceback (most recent call last):
     MetadataError:
 
-    >>> f = FieldmapFile(
-    ...     dsA_dir / "sub-01" / "fmap" / "sub-01_phase2.nii.gz"
-    ... )
-    >>> f.metadata['EchoTime']
-    0.00746
-
     >>> FieldmapFile(
     ...     dsA_dir / "sub-01" / "fmap" / "sub-01_fieldmap.nii.gz",
     ...     find_meta=False
     ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
     Traceback (most recent call last):
     MetadataError:
 
-    >>> f = FieldmapFile(
-    ...     dsA_dir / "sub-01" / "fmap" / "sub-01_fieldmap.nii.gz"
-    ... )
-    >>> f.metadata['Units']
-    'rad/s'
+    Readout timing information is required to estimate PEPOLAR fieldmaps,
+    but it is possible to provide a fallback values.
+    Therefore, warnings are logged if the metadata are missing.
+
+    >>> with caplog.at_level(logging.WARNING, "sdcflows.fieldmaps"):
+    ...     f = FieldmapFile(
+    ...         dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz",
+    ...         metadata={"TotalReadoutTime": None},
+    ...     )  # doctest: +IGNORE_EXCEPTION_DETAIL
+    >>> print(caplog.text)
+    WARNING ... Missing readout timing information ... Explicit fallback must be provided.
+
+    >>> with caplog.at_level(logging.WARNING, "sdcflows.fieldmaps"):
+    ...     f = FieldmapFile(
+    ...         dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz",
+    ...         metadata={"PhaseEncodingDirection": "i", "EstimatedTotalReadoutTime": 0.05},
+    ...         find_meta=False,
+    ...     )  # doctest: +IGNORE_EXCEPTION_DETAIL
+    >>> print(caplog.text)
+    WARNING ... Missing readout timing information ... Estimated timing is available.
 
     """
 
@@ -251,12 +277,19 @@ def __attrs_post_init__(self):
 
             from .utils.epimanip import get_trt
 
+            msg = "Missing readout timing information for <%s>. %s"
+            extra = "Explicit fallback must be provided."
+            have_trt = False
             try:
                 get_trt(self.metadata, in_file=self.path)
-            except ValueError as exc:
-                raise MetadataError(
-                    f"Missing readout timing information for <{self.path}>."
-                ) from exc
+                have_trt = True
+            except ValueError:
+                with suppress(ValueError):
+                    get_trt(self.metadata, in_file=self.path, use_estimate=True)
+                    extra = "Estimated timing is available."
+
+            if not have_trt:
+                logger.warning(msg, self.path, extra)
 
         elif self.suffix == "fieldmap" and "Units" not in self.metadata:
             raise MetadataError(f"Missing 'Units' for <{self.path}>.")

sdcflows/utils/epimanip.py

@@ -244,7 +244,13 @@ def get_trt(
 
         return trt
 
-    if "PhaseEncodingDirection" in in_meta:
+    if in_meta.keys() & {
+        "PhaseEncodingDirection",
+        "EffectiveEchoSpacing",
+        "EstimatedEffectiveEchoSpacing",
+        "EchoSpacing",
+        "WaterFatShift",
+    } > {"PhaseEncodingDirection"}:
         # npe = N voxels PE direction
         pe_index = "ijk".index(in_meta["PhaseEncodingDirection"][0])
         npe = nb.load(in_file).shape[pe_index]

sdcflows/utils/wrangler.py

@@ -623,10 +623,8 @@ def find_anatomical_estimators(
             if not meta.get("PhaseEncodingDirection"):
                 continue
 
-            trt = 1.0
             with suppress(ValueError):
-                trt = get_trt(meta, candidate.path)
-            meta.update({"TotalReadoutTime": trt})
+                meta.update({"TotalReadoutTime": get_trt(meta, candidate.path)})
             epi_targets.append(fm.FieldmapFile(candidate, metadata=meta))
 
         def sort_key(fmap):

sdcflows/workflows/apply/correction.py

@@ -29,6 +29,8 @@
 def init_unwarp_wf(
     *,
     jacobian=True,
+    use_metadata_estimates=False,
+    fallback_total_readout_time=None,
     free_mem=None,
     omp_nthreads=1,
     debug=False,
@@ -117,8 +119,16 @@ def init_unwarp_wf(
         name="outputnode",
     )
 
-    rotime = pe.Node(GetReadoutTime(), name="rotime")
+    rotime = pe.Node(
+        GetReadoutTime(
+            use_estimate=use_metadata_estimates,
+        ),
+        name="rotime",
+        run_without_submitting=True,
+    )
     rotime.interface._always_run = debug
+    if fallback_total_readout_time is not None:
+        rotime.inputs.fallback = fallback_total_readout_time
 
     # resample is memory-hungry; choose a smaller number of threads
     # if we know how much memory we have to work with

sdcflows/workflows/base.py

@@ -21,11 +21,14 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """Estimate fieldmaps for :abbr:`SDC (susceptibility distortion correction)`."""
+
 from nipype import logging
 from nipype.pipeline import engine as pe
 from nipype.interfaces import utility as niu
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
+from ..utils.epimanip import get_trt
+
 LOGGER = logging.getLogger("nipype.workflow")
 DEFAULT_MEMORY_MIN_GB = 0.01
 
@@ -36,6 +39,8 @@ def init_fmap_preproc_wf(
     omp_nthreads,
     output_dir,
     subject,
+    use_metadata_estimates=False,
+    fallback_total_readout_time=None,
     sloppy=False,
     debug=False,
     name="fmap_preproc_wf",
@@ -107,7 +112,22 @@ def init_fmap_preproc_wf(
     )
 
     for n, estimator in enumerate(estimators, 1):
+        if fallback_total_readout_time is None:
+            for f in estimator.sources:
+                if f.suffix in ("bold", "dwi", "epi", "sbref", "asl", "m0scan"):
+                    try:
+                        get_trt(
+                            f.metadata,
+                            in_file=f.path,
+                            use_estimate=use_metadata_estimates,
+                        )
+                    except ValueError:
+                        msg = f"Missing readout timing information for <{f.path}>."
+                        raise RuntimeError(msg)
+
         est_wf = estimator.get_workflow(
+            use_metadata_estimates=use_metadata_estimates,
+            fallback_total_readout_time=fallback_total_readout_time,
             omp_nthreads=omp_nthreads,
             debug=debug,
             sloppy=sloppy,

sdcflows/workflows/fit/pepolar.py

@@ -36,6 +36,8 @@
 
 def init_topup_wf(
     grid_reference=0,
+    use_metadata_estimates=False,
+    fallback_total_readout_time=None,
     omp_nthreads=1,
     sloppy=False,
     debug=False,
@@ -123,11 +125,15 @@ def init_topup_wf(
 
     # Calculate the total readout time of each run
     readout_time = pe.MapNode(
-        GetReadoutTime(),
+        GetReadoutTime(
+            use_estimate=use_metadata_estimates,
+        ),
         name="readout_time",
         iterfield=["metadata", "in_file"],
         run_without_submitting=True,
     )
+    if fallback_total_readout_time is not None:
+        readout_time.inputs.fallback = fallback_total_readout_time
     # Average each run so that topup is not overwhelmed (see #279)
     runwise_avg = pe.MapNode(
         RobustAverage(num_threads=omp_nthreads),

sdcflows/workflows/fit/syn.py

@@ -44,6 +44,8 @@
 
 def init_syn_sdc_wf(
     *,
+    use_metadata_estimates=False,
+    fallback_total_readout_time=None,
     sloppy=False,
     debug=False,
     name="syn_sdc_wf",
@@ -157,10 +159,14 @@ def init_syn_sdc_wf(
     outputnode.inputs.method = 'FLB ("fieldmap-less", SyN-based)'
 
     readout_time = pe.Node(
-        GetReadoutTime(),
+        GetReadoutTime(
+            use_estimate=use_metadata_estimates,
+        ),
         name="readout_time",
         run_without_submitting=True,
     )
+    if fallback_total_readout_time is not None:
+        readout_time.inputs.fallback = fallback_total_readout_time
 
     warp_dir = pe.Node(
         niu.Function(function=_warp_dir),