wip

sdcflows/conftest.py

@@ -12,6 +12,8 @@
 layouts = {p.name: BIDSLayout(str(p), validate=False, derivatives=True)
            for p in Path(test_data_env).glob('*') if p.is_dir()}
 
+data_dir = Path(__file__).parent / "tests" / "data" / "dsA"
+
 
 def pytest_report_header(config):
     msg = "Datasets found: %s" % ', '.join([v.root for v in layouts.values()])
@@ -28,6 +30,8 @@ def add_np(doctest_namespace):
     for key, val in list(layouts.items()):
         doctest_namespace[key] = Path(val.root)
 
+    doctest_namespace['testdata_dir'] = data_dir
+
 
 @pytest.fixture
 def workdir():
@@ -42,3 +46,8 @@ def output_path():
 @pytest.fixture
 def bids_layouts():
     return layouts
+
+
+@pytest.fixture
+def testdata_dir():
+    return data_dir

sdcflows/fieldmaps.py

@@ -0,0 +1,266 @@
+"""Utilities for fieldmap estimation."""
+from pathlib import Path
+from enum import Enum, auto
+import attr
+from bids.layout import BIDSFile, parse_file_entities
+from bids.utils import listify
+from niworkflows.utils.bids import relative_to_root
+
+
+class MetadataError(ValueError):
+    """A better name for a specific value error."""
+
+
+class EstimatorType(Enum):
+    """Represents different types of fieldmap estimation approach."""
+
+    UNKNOWN = auto()
+    PEPOLAR = auto()
+    PHASEDIFF = auto()
+    MAPPED = auto()
+    ANAT = auto()
+
+
+MODALITIES = {
+    "bold": EstimatorType.PEPOLAR,
+    "dwi": EstimatorType.PEPOLAR,
+    "epi": EstimatorType.PEPOLAR,
+    "fieldmap": EstimatorType.MAPPED,
+    "magnitude": None,
+    "magnitude1": None,
+    "magnitude2": None,
+    "phase1": EstimatorType.PHASEDIFF,
+    "phase2": EstimatorType.PHASEDIFF,
+    "phasediff": EstimatorType.PHASEDIFF,
+    "sbref": EstimatorType.PEPOLAR,
+    "T1w": EstimatorType.ANAT,
+    "T2w": EstimatorType.ANAT,
+}
+
+
+def _type_setter(obj, attribute, value):
+    """Make sure the type of estimation is not changed."""
+    if obj.method == value:
+        return value
+
+    if obj.method != EstimatorType.UNKNOWN and obj.method != value:
+        raise TypeError(
+            f"Cannot change determined method {obj.method} to {value}."
+        )
+
+    if value not in (
+        EstimatorType.PEPOLAR,
+        EstimatorType.PHASEDIFF,
+        EstimatorType.MAPPED,
+        EstimatorType.ANAT,
+    ):
+        raise ValueError(f"Invalid estimation method type {value}.")
+
+    return value
+
+
+@attr.s(slots=True)
+class FieldmapFile:
+    """
+    Represent a file that can be used in some fieldmap estimation method.
+
+    Examples
+    --------
+    >>> f = FieldmapFile(testdata_dir / "sub-01" / "anat" / "sub-01_T1w.nii.gz")
+    >>> f.suffix
+    'T1w'
+
+    >>> FieldmapFile(
+    ...     testdata_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz"
+    ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
+    Traceback (most recent call last):
+    MetadataError:
+
+    >>> f = FieldmapFile(
+    ...     testdata_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz",
+    ...     metadata={'PhaseEncodingDirection': 'i', 'TotalReadoutTime': 0.006}
+    ... )
+    >>> f.metadata['PhaseEncodingDirection']
+    'i'
+
+    >>> FieldmapFile(
+    ...     testdata_dir / "sub-01" / "fmap" / "sub-01_phasediff.nii.gz"
+    ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
+    Traceback (most recent call last):
+    MetadataError:
+
+    >>> f = FieldmapFile(
+    ...     testdata_dir / "sub-01" / "fmap" / "sub-01_phasediff.nii.gz",
+    ...     metadata={'EchoTime1': 0.005, 'EchoTime2': 0.00746}
+    ... )
+    >>> f.metadata['EchoTime2']
+    0.00746
+
+    >>> FieldmapFile(
+    ...     testdata_dir / "sub-01" / "fmap" / "sub-01_phase2.nii.gz"
+    ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
+    Traceback (most recent call last):
+    MetadataError:
+
+    >>> f = FieldmapFile(
+    ...     testdata_dir / "sub-01" / "fmap" / "sub-01_phase2.nii.gz",
+    ...     metadata={'EchoTime': 0.00746}
+    ... )
+    >>> f.metadata['EchoTime']
+    0.00746
+
+    >>> FieldmapFile(
+    ...     testdata_dir / "sub-01" / "fmap" / "sub-01_fieldmap.nii.gz"
+    ... )  # doctest: +IGNORE_EXCEPTION_DETAIL
+    Traceback (most recent call last):
+    MetadataError:
+
+    >>> f = FieldmapFile(
+    ...     testdata_dir / "sub-01" / "fmap" / "sub-01_fieldmap.nii.gz",
+    ...     metadata={'Units': "rad/s"}
+    ... )
+    >>> f.metadata['Units']
+    'rad/s'
+
+    """
+
+    path = attr.ib(converter=Path, repr=str, on_setattr=attr.setters.NO_OP)
+    """Path to a fieldmap file."""
+
+    entities = attr.ib(init=False, repr=False)
+    """BIDS entities extracted from filepath."""
+
+    suffix = attr.ib(init=False, repr=False)
+    """Extracted suffix from input file."""
+
+    bids_root = attr.ib(init=False, default=None, repr=False)
+    """Path of the BIDS root."""
+
+    metadata = attr.ib(kw_only=True, default=attr.Factory(dict))
+    """Metadata associated to this file."""
+
+    @path.validator
+    def check_path(self, attribute, value):
+        """Validate a fieldmap path."""
+        if isinstance(value, BIDSFile):
+            value = Path(value.path)
+        if isinstance(value, str):
+            value = Path(value)
+
+        if not value.is_file():
+            raise FileNotFoundError(
+                f"File path <{value}> does not exist, is a broken link, or it is not a file"
+            )
+
+        if not str(value).endswith((".nii", ".nii.gz")):
+            raise ValueError(f"File path <{value}> does not look like a NIfTI file.")
+
+    def __attrs_post_init__(self):
+        """Initialize the suffix property."""
+        self.entities = parse_file_entities(str(self.path))
+        suffix = self.entities.pop('suffix')
+        if suffix not in tuple(MODALITIES.keys()):
+            raise ValueError(
+                f"File path <{self.path}> with suffix <{suffix}> is not a valid "
+                "fieldmap sourcefile."
+            )
+        self.suffix = suffix
+
+        relative_path = relative_to_root(self.path)
+        if str(relative_path) != str(self.path):
+            self.bids_root = Path(str(self.path)[:-len(str(relative_path))])
+
+        if suffix in ("bold", "dwi", "epi", "sbref"):
+            if "PhaseEncodingDirection" not in self.metadata:
+                raise MetadataError(f"Missing 'PhaseEncodingDirection' for <{self.path}>.")
+            if not (
+                set(("TotalReadoutTime", "EffectiveEchoSpacing")).intersection(
+                    self.metadata.keys())
+            ):
+                raise MetadataError(f"Missing timing information for <{self.path}>.")
+
+        if suffix == "fieldmap" and "Units" not in self.metadata:
+            raise MetadataError(f"Missing 'Units' for <{self.path}>.")
+
+        if (
+            suffix == "phasediff"
+            and ("EchoTime1" not in self.metadata or "EchoTime2" not in self.metadata)
+        ):
+            raise MetadataError(f"Missing 'EchoTime1' and/or 'EchoTime2' for <{self.path}>.")
+
+        if suffix in ("phase1", "phase2") and ("EchoTime" not in self.metadata):
+            raise MetadataError(f"Missing 'EchoTime' for <{self.path}>.")
+
+
+@attr.s(slots=True)
+class FieldmapEstimation:
+    """
+    Represent fieldmap estimation strategies.
+
+    This class provides a consistent interface to all types of fieldmap estimation
+    strategies.
+    The actual type of method for estimation is inferred from the ``sources`` input,
+    and collects all the available metadata.
+    This class also checks all the data and metadata:
+
+      - Check that all files in ``sources`` exist
+      - Check that all files necessary for the inferred estimation procedure are
+        found within the ``sources`` attribute.
+      - Combine all metadata (and allow setting custom metadata), as well as provide a
+        ``validate()`` member that checks all necessary metadata for estimation is set.
+
+
+    """
+
+    sources = attr.ib(
+        default=None,
+        converter=lambda v: [
+            FieldmapFile(f) if not isinstance(f, FieldmapFile) else f
+            for f in listify(v)
+        ],
+        repr=lambda v: f"<{len(v)} files>",
+    )
+    """File path or list of paths indicating the source data to estimate a fieldmap."""
+
+    method = attr.ib(init=False, default=EstimatorType.UNKNOWN, on_setattr=_type_setter)
+
+    def __attrs_post_init__(self):
+        """Determine the inteded fieldmap estimation type and check for data completeness."""
+        suffix_list = [f.suffix for f in self.sources]
+        suffix_set = set(suffix_list)
+        fmap_types = suffix_set.intersection(("fieldmap", "phasediff", "phase1"))
+        if len(fmap_types) > 1:
+            raise ValueError(f"Incompatible suffices found: <{','.join(fmap_types)}>.")
+
+        if fmap_types:
+            self.method = MODALITIES[fmap_types.pop()]
+
+        pepolar_types = suffix_set.intersection(("bold", "dwi", "epi", "sbref"))
+        _pepolar_estimation = len([
+            f for f in suffix_list if f in ("bold", "dwi", "epi", "sbref")
+        ]) > 1
+
+        if _pepolar_estimation:
+            self.method = MODALITIES[pepolar_types.pop()]
+
+        anat_types = suffix_set.intersection(("T1w", "T2w"))
+        if anat_types:
+            self.method = MODALITIES[anat_types.pop()]
+
+        if fmap_types and (_pepolar_estimation or anat_types):
+            raise ValueError("Incompatible suffices found in combination with a fieldmap type.")
+        if fmap_types and not suffix_set.intersection(("magnitude", "magnitude1")):
+            raise ValueError(
+                "A phase-difference estimation type was found, "
+                "but an anatomical reference (magnitude file) is missing."
+            )
+        if anat_types and not pepolar_types:
+            raise ValueError("Only anatomical sources were found, cannot estimate fieldmap.")
+
+        if pepolar_types and not (anat_types or _pepolar_estimation):
+            raise ValueError(
+                "One EPI file was provided, but more data is necessary for estimation."
+            )
+
+        if self.method == EstimatorType.UNKNOWN:
+            raise ValueError("Insufficient sources to estimate a fieldmap.")

sdcflows/tests/test_fieldmaps.py

@@ -0,0 +1,68 @@
+"""test_fieldmaps."""
+import pytest
+from ..fieldmaps import FieldmapFile, FieldmapEstimation
+
+
+def test_FieldmapFile(testdata_dir):
+    """Test one existing file."""
+    FieldmapFile(testdata_dir / "sub-01" / "anat" / "sub-01_T1w.nii.gz")
+
+
+def test_FieldmapEstimation(testdata_dir):
+    """Test errors."""
+
+    sub_dir = testdata_dir / "sub-01"
+
+    with pytest.raises(ValueError):
+        FieldmapEstimation([
+            sub_dir / "fmap" / "sub-01_fieldmap.nii.gz",
+            sub_dir / "fmap" / "sub-01_phasediff.nii.gz",
+        ])
+
+    FieldmapEstimation([
+        sub_dir / "fmap" / "sub-01_fieldmap.nii.gz",
+        sub_dir / "fmap" / "sub-01_magnitude.nii.gz",
+    ])
+
+    with pytest.raises(ValueError):
+        FieldmapEstimation([
+            sub_dir / "fmap" / "sub-01_phase1.nii.gz",
+            sub_dir / "fmap" / "sub-01_phase2.nii.gz",
+        ])
+
+    FieldmapEstimation([
+        sub_dir / "fmap" / "sub-01_phase1.nii.gz",
+        sub_dir / "fmap" / "sub-01_phase2.nii.gz",
+        sub_dir / "fmap" / "sub-01_magnitude1.nii.gz",
+        sub_dir / "fmap" / "sub-01_magnitude2.nii.gz",
+    ])
+
+    with pytest.raises(ValueError):
+        FieldmapEstimation([
+            sub_dir / "anat" / "sub-01_T1w.nii.gz",
+        ])
+
+    with pytest.raises(ValueError):
+        FieldmapEstimation([
+            sub_dir / "anat" / "sub-01_T1w.nii.gz",
+            sub_dir / "fmap" / "sub-01_phase2.nii.gz",
+        ])
+
+    FieldmapEstimation([
+        sub_dir / "anat" / "sub-01_T1w.nii.gz",
+        sub_dir / "dwi" / "sub-01_dir-LR_sbref.nii.gz",
+    ])
+
+    with pytest.raises(ValueError):
+        FieldmapEstimation([
+            sub_dir / "fmap" / "sub-01_dir-RL_epi.nii.gz",
+        ])
+
+    FieldmapEstimation([
+        sub_dir / "fmap" / "sub-01_dir-LR_epi.nii.gz",
+        sub_dir / "fmap" / "sub-01_dir-RL_epi.nii.gz",
+    ])
+    FieldmapEstimation([
+        sub_dir / "fmap" / "sub-01_dir-LR_epi.nii.gz",
+        sub_dir / "dwi" / "sub-01_dir-RL_sbref.nii.gz",
+    ])