Select ISMRMRD acquisitions with a function instead of ignore flags (#…

…315)

src/mrpro/data/KNoise.py

@@ -17,13 +17,14 @@
 from __future__ import annotations
 
 import dataclasses
+from collections.abc import Callable
 from pathlib import Path
 
 import ismrmrd
 import torch
 from einops import rearrange
 
-from mrpro.data.enums import AcqFlags
+from mrpro.data.acq_filters import is_noise_acquisition
 from mrpro.data.MoveDataMixin import MoveDataMixin
 
 
@@ -40,7 +41,9 @@ class KNoise(MoveDataMixin):
     data: torch.Tensor
 
     @classmethod
-    def from_file(cls, filename: str | Path, dataset_idx: int = -1) -> KNoise:
+    def from_file(
+        cls, filename: str | Path, dataset_idx: int = -1, acquisition_filter_criterion: Callable = is_noise_acquisition
+    ) -> KNoise:
         """Load noise measurements from ISMRMRD file.
 
         Parameters
@@ -49,14 +52,16 @@ def from_file(cls, filename: str | Path, dataset_idx: int = -1) -> KNoise:
             Path to the ISMRMRD file
         dataset_idx
             Index of the dataset to load (converter creates dataset, dataset_1, ...), default is -1 (last)
+        acquisition_filter_criterion
+            function which returns True if an acquisition should be included in KNoise
         """
         # Can raise FileNotFoundError
         with ismrmrd.File(filename, 'r') as file:
             ds = file[list(file.keys())[dataset_idx]]
             acquisitions = ds.acquisitions[:]
 
         # Read out noise measurements
-        acquisitions = list(filter(lambda acq: (AcqFlags.ACQ_IS_NOISE_MEASUREMENT.value & acq.flags), acquisitions))
+        acquisitions = [acq for acq in acquisitions if acquisition_filter_criterion(acq)]
         if len(acquisitions) == 0:
             raise ValueError(f'No noise measurements found in {filename}')
         noise_data = torch.stack([torch.as_tensor(acq.data, dtype=torch.complex64) for acq in acquisitions])

src/mrpro/data/__init__.py

@@ -1,5 +1,6 @@
 from mrpro.data import enums
 from mrpro.data import traj_calculators
+from mrpro.data import acq_filters
 from mrpro.data.AcqInfo import AcqIdx
 from mrpro.data.AcqInfo import AcqInfo
 from mrpro.data.CsmData import CsmData

src/mrpro/data/_kdata/KData.py

@@ -19,6 +19,7 @@
 import dataclasses
 import datetime
 import warnings
+from collections.abc import Callable
 from pathlib import Path
 from typing import Protocol
 
@@ -31,9 +32,9 @@
 from mrpro.data._kdata.KDataRemoveOsMixin import KDataRemoveOsMixin
 from mrpro.data._kdata.KDataSelectMixin import KDataSelectMixin
 from mrpro.data._kdata.KDataSplitMixin import KDataSplitMixin
+from mrpro.data.acq_filters import is_image_acquisition
 from mrpro.data.AcqInfo import AcqInfo
 from mrpro.data.EncodingLimits import Limits
-from mrpro.data.enums import AcqFlags
 from mrpro.data.KHeader import KHeader
 from mrpro.data.KTrajectory import KTrajectory
 from mrpro.data.KTrajectoryRawShape import KTrajectoryRawShape
@@ -46,18 +47,6 @@
 # TODO: Consider adding the users labels here, but remember issue #32 and NOT add user5 and user6.
 OTHER_LABELS = ('average', 'slice', 'contrast', 'phase', 'repetition', 'set')
 
-# Same criteria as https://github.com/wtclarke/pymapvbvd/blob/master/mapvbvd/mapVBVD.py uses
-DEFAULT_IGNORE_FLAGS = (
-    AcqFlags.ACQ_IS_NOISE_MEASUREMENT
-    | AcqFlags.ACQ_IS_DUMMYSCAN_DATA
-    | AcqFlags.ACQ_IS_HPFEEDBACK_DATA
-    | AcqFlags.ACQ_IS_NAVIGATION_DATA
-    | AcqFlags.ACQ_IS_PHASECORR_DATA
-    | AcqFlags.ACQ_IS_PHASE_STABILIZATION
-    | AcqFlags.ACQ_IS_PHASE_STABILIZATION_REFERENCE
-    | AcqFlags.ACQ_IS_PARALLEL_CALIBRATION
-)
-
 
 @dataclasses.dataclass(slots=True, frozen=True)
 class KData(KDataSplitMixin, KDataRearrangeMixin, KDataSelectMixin, KDataRemoveOsMixin, MoveDataMixin):
@@ -74,7 +63,7 @@ def from_file(
         ktrajectory: KTrajectoryCalculator | KTrajectory | KTrajectoryIsmrmrd,
         header_overwrites: dict[str, object] | None = None,
         dataset_idx: int = -1,
-        ignore_flags: AcqFlags = DEFAULT_IGNORE_FLAGS,
+        acquisition_filter_criterion: Callable = is_image_acquisition,
     ) -> KData:
         """Load k-space data from an ISMRMRD file.
 
@@ -88,11 +77,8 @@ def from_file(
             dictionary of key-value pairs to overwrite the header
         dataset_idx
             index of the ISMRMRD dataset to load (converter creates dataset, dataset_1, ...), default is -1 (last)
-        ignore_flags
-            Acqisition flags to filter out. Defaults to all non-images as defined by pymapvbvd.
-            Use ACQ_NO_FLAG to disable the filter.
-            Note: If ACQ_IS_PARALLEL_CALIBRATION is set without also setting ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING
-                    we interpret it as: Ignore if IS_PARALLEL_CALIBRATION and not PARALLEL_CALIBRATION_AND_IMAGING
+        acquisition_filter_criterion
+            function which returns True if an acquisition should be included in KData
         """
         # Can raise FileNotFoundError
         with ismrmrd.File(filename, 'r') as file:
@@ -105,22 +91,7 @@ def from_file(
                 mtime = 0
             modification_time = datetime.datetime.fromtimestamp(mtime)
 
-        if (
-            AcqFlags.ACQ_IS_PARALLEL_CALIBRATION in ignore_flags
-            and AcqFlags.ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING not in ignore_flags
-        ):
-            # if only ACQ_IS_PARALLEL_CALIBRATION is set, reinterpret it as: ignore if
-            # ACQ_IS_PARALLEL_CALIBRATION is set and ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING is not set
-            ignore_flags = ignore_flags & ~AcqFlags.ACQ_IS_PARALLEL_CALIBRATION
-            acquisitions = list(
-                filter(
-                    lambda acq: (AcqFlags.ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING.value & acq.flags)
-                    or not (AcqFlags.ACQ_IS_PARALLEL_CALIBRATION.value & acq.flags),
-                    acquisitions,
-                ),
-            )
-
-        acquisitions = list(filter(lambda acq: not (ignore_flags.value & acq.flags), acquisitions))
+        acquisitions = [acq for acq in acquisitions if acquisition_filter_criterion(acq)]
         kdata = torch.stack([torch.as_tensor(acq.data, dtype=torch.complex64) for acq in acquisitions])
 
         acqinfo = AcqInfo.from_ismrmrd_acquisitions(acquisitions)

src/mrpro/data/acq_filters.py

@@ -0,0 +1,77 @@
+"""Test ISMRMRD acquisitions based on their flags."""
+
+# Copyright 2024 Physikalisch-Technische Bundesanstalt
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at:
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import ismrmrd
+
+from mrpro.data.enums import AcqFlags
+
+# Same criteria as https://github.com/wtclarke/pymapvbvd/blob/master/mapvbvd/mapVBVD.py uses
+DEFAULT_IGNORE_FLAGS = (
+    AcqFlags.ACQ_IS_NOISE_MEASUREMENT
+    | AcqFlags.ACQ_IS_DUMMYSCAN_DATA
+    | AcqFlags.ACQ_IS_HPFEEDBACK_DATA
+    | AcqFlags.ACQ_IS_NAVIGATION_DATA
+    | AcqFlags.ACQ_IS_PHASECORR_DATA
+    | AcqFlags.ACQ_IS_PHASE_STABILIZATION
+    | AcqFlags.ACQ_IS_PHASE_STABILIZATION_REFERENCE
+    | AcqFlags.ACQ_IS_PARALLEL_CALIBRATION
+)
+
+
+def is_image_acquisition(acquisition: ismrmrd.Acquisition) -> bool:
+    """Test if acquisition was obtained for imaging.
+
+    Parameters
+    ----------
+    acquisition
+        ISMRMRD acquisition
+
+    Returns
+    -------
+        True if the acquisition was obtained for imaging
+    """
+    return not DEFAULT_IGNORE_FLAGS.value & acquisition.flags
+
+
+def is_noise_acquisition(acquisition: ismrmrd.Acquisition) -> bool:
+    """Test if acquisition contains noise measurements.
+
+    Parameters
+    ----------
+    acquisition
+        ISMRMRD acquisition
+
+    Returns
+    -------
+        True if the acquisition contains a noise measurement
+    """
+    return AcqFlags.ACQ_IS_NOISE_MEASUREMENT.value & acquisition.flags
+
+
+def is_coil_calibration_acquisition(acquisition: ismrmrd.Acquisition) -> bool:
+    """Test if acquisitions was obtained for coil calibration.
+
+    Parameters
+    ----------
+    acquisition
+        ISMRMRD acquisition
+
+    Returns
+    -------
+        True if the acquisition contains coil calibration data
+    """
+    coil_calibration_flag = AcqFlags.ACQ_IS_PARALLEL_CALIBRATION | AcqFlags.ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING
+    return coil_calibration_flag.value & acquisition.flags

tests/data/_IsmrmrdRawTestData.py

@@ -45,25 +45,27 @@ class IsmrmrdRawTestData:
     filename
         full path and filename
     matrix_size
-        size of image matrix, by default 256
+        size of image matrix
     n_coils
-        number of coils, by default 8
+        number of coils
     oversampling
-        oversampling along readout (kx) direction, by default 2
+        oversampling along readout (kx) direction
     repetitions
-        number of repetitions, by default 1
+        number of repetitions,
     flag_invalid_reps
-        flag to indicate that number of phase encoding steps are different for repetitions, by default False
+        flag to indicate that number of phase encoding steps are different for repetitions
     acceleration
-        undersampling along phase encoding (ky), by default 1
+        undersampling along phase encoding (ky)
     noise_level
-        scaling factor for noise level, by default 0.00005
+        scaling factor for noise level
     trajectory_type
-        cartesian, by default cartesian
+        cartesian
     sampling_order
-        order how phase encoding points (ky) are obtained, by default linear
+        order how phase encoding points (ky) are obtained
     phantom
         phantom with different ellipses
+    n_separate_calibration_lines
+        number of additional calibration lines, linear Cartesian sampled
     """
 
     def __init__(
@@ -79,6 +81,7 @@ def __init__(
         trajectory_type: Literal['cartesian', 'radial'] = 'cartesian',
         sampling_order: Literal['linear', 'low_high', 'high_low', 'random'] = 'linear',
         phantom: EllipsePhantom | None = None,
+        n_separate_calibration_lines: int = 0,
     ):
         if not phantom:
             phantom = EllipsePhantom()
@@ -95,6 +98,7 @@ def __init__(
         self.sampling_order: Literal['linear', 'low_high', 'high_low', 'random'] = sampling_order
         self.phantom: EllipsePhantom = phantom
         self.img_ref: torch.Tensor
+        self.n_separate_calibration_lines: int = n_separate_calibration_lines
 
         # The number of points in image space (x,y) and kspace (fe,pe)
         n_x = self.matrix_size
@@ -246,6 +250,36 @@ def __init__(
             dataset.append_acquisition(acq)
             counter += 1  # increment the scan counter
 
+        # Calibration lines
+        if n_separate_calibration_lines > 0:
+            # we take calibration lines around the k-space center
+            traj_ky_calibration, traj_kx_calibration, kpe_calibration = self._cartesian_trajectory(
+                n_separate_calibration_lines,
+                n_freq_encoding,
+                1,
+                'linear',
+            )
+            kspace_calibration = self.phantom.kspace(traj_ky_calibration, traj_kx_calibration)
+            kspace_calibration = repeat(kspace_calibration, '... -> coils ... ', coils=self.n_coils)
+            kspace_calibration = kspace_calibration + self.noise_level * torch.randn(
+                self.n_coils, n_freq_encoding, len(kpe_calibration), dtype=torch.complex64
+            )
+
+            for pe_idx, pe_pos in enumerate(kpe_calibration):
+                # Set some fields in the header
+                acq.scan_counter = counter
+
+                # kpe is in the range [-npe//2, npe//2), the ismrmrd kspace_encoding_step_1 is in the range [0, npe)
+                kspace_encoding_step_1 = pe_pos + n_phase_encoding // 2
+                acq.idx.kspace_encode_step_1 = kspace_encoding_step_1
+                acq.clearAllFlags()
+                acq.setFlag(ismrmrd.ACQ_IS_PARALLEL_CALIBRATION)
+
+                # Set the data and append
+                acq.data[:] = kspace_calibration[:, :, pe_idx].numpy()
+                dataset.append_acquisition(acq)
+                counter += 1
+
         # Loop over the repetitions, add noise and write to disk
         for rep in range(self.repetitions):
             noise = self.noise_level * torch.randn(self.n_coils, n_freq_encoding, len(kpe[rep]), dtype=torch.complex64)

tests/data/test_kdata.py

@@ -19,6 +19,7 @@
 from mrpro.data import KData
 from mrpro.data import KTrajectory
 from mrpro.data import SpatialDimension
+from mrpro.data.acq_filters import is_coil_calibration_acquisition
 from mrpro.data.traj_calculators.KTrajectoryCalculator import DummyTrajectory
 from mrpro.operators import FastFourierOp
 from mrpro.utils import modify_acq_info
@@ -44,6 +45,21 @@ def ismrmrd_cart(ellipse_phantom, tmp_path_factory):
     return ismrmrd_kdata
 
 
+@pytest.fixture(scope='session')
+def ismrmrd_cart_with_calibration_lines(ellipse_phantom, tmp_path_factory):
+    """Undersampled Cartesian data set with calibration lines."""
+    ismrmrd_filename = tmp_path_factory.mktemp('mrpro') / 'ismrmrd_cart.h5'
+    ismrmrd_kdata = IsmrmrdRawTestData(
+        filename=ismrmrd_filename,
+        noise_level=0.0,
+        repetitions=1,
+        acceleration=2,
+        phantom=ellipse_phantom.phantom,
+        n_separate_calibration_lines=16,
+    )
+    return ismrmrd_kdata
+
+
 @pytest.fixture(scope='session')
 def ismrmrd_cart_invalid_reps(tmp_path_factory):
     """Fully sampled cartesian data set."""
@@ -133,6 +149,24 @@ def test_KData_from_file_diff_nky_for_rep(ismrmrd_cart_invalid_reps):
     assert kdata.data.shape[-3] == 1, 'k2 should be 1'
 
 
+def test_KData_calibration_lines(ismrmrd_cart_with_calibration_lines):
+    """Correct handling of calibration lines."""
+    # Exclude calibration lines
+    kdata = KData.from_file(ismrmrd_cart_with_calibration_lines.filename, DummyTrajectory())
+    assert (
+        kdata.data.shape[-2]
+        == ismrmrd_cart_with_calibration_lines.matrix_size // ismrmrd_cart_with_calibration_lines.acceleration
+    )
+
+    # Get only calibration lines
+    kdata = KData.from_file(
+        ismrmrd_cart_with_calibration_lines.filename,
+        DummyTrajectory(),
+        acquisition_filter_criterion=is_coil_calibration_acquisition,
+    )
+    assert kdata.data.shape[-2] == ismrmrd_cart_with_calibration_lines.n_separate_calibration_lines
+
+
 def test_KData_kspace(ismrmrd_cart):
     """Read in data and verify k-space by comparing reconstructed image."""
     kdata = KData.from_file(ismrmrd_cart.filename, DummyTrajectory())