diff --git a/src/mrpro/data/KNoise.py b/src/mrpro/data/KNoise.py index ec60160a..32f36a72 100644 --- a/src/mrpro/data/KNoise.py +++ b/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,6 +52,8 @@ 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: @@ -56,7 +61,7 @@ def from_file(cls, filename: str | Path, dataset_idx: int = -1) -> KNoise: 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]) diff --git a/src/mrpro/data/__init__.py b/src/mrpro/data/__init__.py index d84d9e35..59a04984 100644 --- a/src/mrpro/data/__init__.py +++ b/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 diff --git a/src/mrpro/data/_kdata/KData.py b/src/mrpro/data/_kdata/KData.py index f50159ad..6d97ffc6 100644 --- a/src/mrpro/data/_kdata/KData.py +++ b/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) diff --git a/src/mrpro/data/acq_filters.py b/src/mrpro/data/acq_filters.py new file mode 100644 index 00000000..2f1b28cb --- /dev/null +++ b/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 diff --git a/tests/data/_IsmrmrdRawTestData.py b/tests/data/_IsmrmrdRawTestData.py index fdcf06a4..dbc93c3f 100644 --- a/tests/data/_IsmrmrdRawTestData.py +++ b/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) diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index fbcd19da..a5f3affb 100644 --- a/tests/data/test_kdata.py +++ b/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())