diff --git a/sdcflows/conftest.py b/sdcflows/conftest.py index bfcc77ae86..2825b6b587 100644 --- a/sdcflows/conftest.py +++ b/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()) diff --git a/sdcflows/fieldmaps.py b/sdcflows/fieldmaps.py index b087731a07..e6f35d206b 100644 --- a/sdcflows/fieldmaps.py +++ b/sdcflows/fieldmaps.py @@ -24,6 +24,7 @@ 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 @@ -31,8 +32,11 @@ 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,6 +150,15 @@ 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 @@ -143,12 +166,6 @@ class FieldmapFile: 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 @@ -156,12 +173,6 @@ class FieldmapFile: 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 @@ -169,11 +180,26 @@ class FieldmapFile: 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}>.") diff --git a/sdcflows/utils/epimanip.py b/sdcflows/utils/epimanip.py index c8f08e45b9..2980693cbc 100644 --- a/sdcflows/utils/epimanip.py +++ b/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] diff --git a/sdcflows/utils/wrangler.py b/sdcflows/utils/wrangler.py index 67743b496d..977222dc58 100644 --- a/sdcflows/utils/wrangler.py +++ b/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): diff --git a/sdcflows/workflows/apply/correction.py b/sdcflows/workflows/apply/correction.py index fb20e6765c..6c99727a01 100644 --- a/sdcflows/workflows/apply/correction.py +++ b/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 diff --git a/sdcflows/workflows/base.py b/sdcflows/workflows/base.py index 3c3d97ebe8..95d3b6273e 100644 --- a/sdcflows/workflows/base.py +++ b/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, diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py index b0f6bd59fc..a353da65d3 100644 --- a/sdcflows/workflows/fit/pepolar.py +++ b/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), diff --git a/sdcflows/workflows/fit/syn.py b/sdcflows/workflows/fit/syn.py index 1ee021d495..3d8e87d5f5 100644 --- a/sdcflows/workflows/fit/syn.py +++ b/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),