Add save kdata to ismrmrd file

src/mrpro/data/AcqInfo.py

@@ -9,10 +9,11 @@
 from einops import rearrange
 from typing_extensions import Self
 
+import mrpro.utils.typing as type_utils
 from mrpro.data.MoveDataMixin import MoveDataMixin
 from mrpro.data.Rotation import Rotation
 from mrpro.data.SpatialDimension import SpatialDimension
-from mrpro.utils.unit_conversion import mm_to_m
+from mrpro.utils.unit_conversion import m_to_mm, mm_to_m
 
 
 def rearrange_acq_info_fields(field: object, pattern: str, **axes_lengths: dict[str, int]) -> object:
@@ -258,3 +259,51 @@ def spatialdimension_2d(data: np.ndarray) -> SpatialDimension[torch.Tensor]:
             version=tensor_2d(headers['version']),
         )
         return acq_info
+
+    def add_to_ismrmrd_acquisition(
+        self, acquisition: ismrmrd.Acquisition, other: type_utils.TorchIndexerType, k2: int, k1: int
+    ) -> ismrmrd.Acquisition:
+        """ISMRMRD acquisition information for single acquisition."""
+        acquisition.idx.kspace_encode_step_1 = self.idx.k1[*other, k2, k1]
+        acquisition.idx.kspace_encode_step_2 = self.idx.k2[*other, k2, k1]
+        acquisition.idx.average = self.idx.average[*other, k2, k1]
+        acquisition.idx.slice = self.idx.slice[*other, k2, k1]
+        acquisition.idx.contrast = self.idx.contrast[*other, k2, k1]
+        acquisition.idx.phase = self.idx.phase[*other, k2, k1]
+        acquisition.idx.repetition = self.idx.repetition[*other, k2, k1]
+        acquisition.idx.set = self.idx.set[*other, k2, k1]
+        acquisition.idx.segment = self.idx.segment[*other, k2, k1]
+        acquisition.idx.user = (
+            self.idx.user0[*other, k2, k1],
+            self.idx.user1[*other, k2, k1],
+            self.idx.user2[*other, k2, k1],
+            self.idx.user3[*other, k2, k1],
+            self.idx.user4[*other, k2, k1],
+            self.idx.user5[*other, k2, k1],
+            self.idx.user6[*other, k2, k1],
+            self.idx.user7[*other, k2, k1],
+        )
+
+        # active_channesl, number_of_samples and trajectory_dimensions are read-only and cannot be set
+        acquisition.acquisition_time_stamp = self.acquisition_time_stamp[*other, k2, k1, 0]
+        acquisition.available_channels = self.available_channels[*other, k2, k1, 0]
+        acquisition.center_sample = self.center_sample[*other, k2, k1, 0]
+        acquisition.channel_mask = tuple(self.channel_mask[*other, k2, k1, :])
+        acquisition.discard_post = self.discard_post[*other, k2, k1, 0]
+        acquisition.discard_pre = self.discard_pre[*other, k2, k1, 0]
+        acquisition.encoding_space_ref = self.encoding_space_ref[*other, k2, k1, 0]
+        acquisition.measurement_uid = self.measurement_uid[*other, k2, k1, 0]
+        acquisition.patient_table_position = (
+            self.patient_table_position[*other, k2, k1, 0].apply_(m_to_mm).zyx[::-1]
+        )  # zyx -> xyz
+        acquisition.phase_dir = self.orientation.as_directions()[-2][*other, k2, k1, 0].zyx[::-1]  # zyx -> xyz
+        acquisition.physiology_time_stamp = tuple(self.physiology_time_stamp[*other, k2, k1, :])
+        acquisition.position = self.position[*other, k2, k1, 0].apply_(m_to_mm).zyx[::-1]  # zyx -> xyz
+        acquisition.read_dir = self.orientation.as_directions()[-1][*other, k2, k1, 0].zyx[::-1]  # zyx -> xyz
+        acquisition.sample_time_us = self.sample_time_us[*other, k2, k1, 0]
+        acquisition.scan_counter = self.scan_counter[*other, k2, k1, 0]
+        acquisition.slice_dir = self.orientation.as_directions()[-3][*other, k2, k1, 0].zyx[::-1]  # zyx -> xyz
+        acquisition.user_float = tuple(self.user_float[*other, k2, k1, :])
+        acquisition.user_int = tuple(self.user_int[*other, k2, k1, :])
+        acquisition.version = self.version[*other, k2, k1, 0]
+        return acquisition

src/mrpro/data/EncodingLimits.py

@@ -27,6 +27,10 @@ def from_ismrmrd(cls, limit_type: limitType) -> Self:
             return cls()
         return cls(*dataclasses.astuple(limit_type))
 
+    def to_ismrmrd(self) -> limitType:
+        """Convert Limits to ismsmrd.limitType."""
+        return limitType(self.min, self.max, self.center)
+
     @property
     def length(self) -> int:
         """Length of the limits."""
@@ -112,3 +116,14 @@ def from_ismrmrd_encoding_limits_type(cls, encoding_limits: encodingLimitsType):
         values['k2'] = values.pop('kspace_encoding_step_2')
 
         return cls(**values)
+
+    def to_ismrmrd_encoding_limits_type(self) -> encodingLimitsType:
+        """Convert EncodingLimits to encodingLimitsType."""
+        values = {field.name: Limits.to_ismrmrd(getattr(self, field.name)) for field in dataclasses.fields(self)}
+
+        # adjust from MRPro to ISMRMRD naming convention
+        values['kspace_encoding_step_0'] = values.pop('k0')
+        values['kspace_encoding_step_1'] = values.pop('k1')
+        values['kspace_encoding_step_2'] = values.pop('k2')
+
+        return encodingLimitsType(**values)

src/mrpro/data/KHeader.py

@@ -18,7 +18,7 @@
 from mrpro.data.MoveDataMixin import MoveDataMixin
 from mrpro.data.SpatialDimension import SpatialDimension
 from mrpro.utils.summarize_tensorvalues import summarize_tensorvalues
-from mrpro.utils.unit_conversion import mm_to_m, ms_to_s
+from mrpro.utils.unit_conversion import m_to_mm, mm_to_m, ms_to_s, s_to_ms
 
 if TYPE_CHECKING:
     # avoid circular imports by importing only when type checking
@@ -255,6 +255,97 @@ def from_ismrmrd(
             ) from None
         return instance
 
+    def to_ismrmrd(self) -> ismrmrdschema.ismrmrdHeader:
+        """Create ISMRMRD header from KHeader.
+
+        All the parameters from the ISMRMRD header of the raw file used to create the kheader are saved in the
+        dictionary kheader['misc']. We first fill the information from this dictionary into the ISMRMRD header and then
+        we update the values based on the attributes of the KHeader.
+
+        Returns
+        -------
+            ISMRMRD header
+        """
+        header = ismrmrdschema.ismrmrdHeader()
+
+        # Experimental conditions
+        exp = ismrmrdschema.experimentalConditionsType()
+        exp.H1resonanceFrequency_Hz = self.lamor_frequency_proton
+        header.experimentalConditions = exp
+
+        # Subject information
+        subj = ismrmrdschema.subjectInformationType()
+        subj.patientName = self.patient_name
+        header.subjectInformation = subj
+
+        # Measurement information
+        meas = ismrmrdschema.measurementInformationType()
+        meas.protocolName = self.protocol_name
+        meas.measurementID = self.measurement_id
+        if self.datetime is not None:
+            meas.seriesTime = str(self.datetime.time())
+            meas.seriesDate = str(self.datetime.date())
+        header.measurementInformation = meas
+
+        # Acquisition system information
+        sys = ismrmrdschema.acquisitionSystemInformationType()
+        sys.systemModel = self.model
+        sys.systemVendor = self.vendor
+        header.acquisitionSystemInformation = sys
+
+        # Sequence information
+        seq = ismrmrdschema.sequenceParametersType()
+        if isinstance(self.tr, torch.Tensor):
+            seq.TR = s_to_ms(self.tr).tolist()
+        if isinstance(self.te, torch.Tensor):
+            seq.TE = s_to_ms(self.te).tolist()
+        if isinstance(self.ti, torch.Tensor):
+            seq.TI = s_to_ms(self.ti).tolist()
+        if isinstance(self.fa, torch.Tensor):
+            seq.flipAngle_deg = torch.rad2deg(self.fa).tolist()
+        if isinstance(self.echo_spacing, torch.Tensor):
+            seq.echo_spacing = s_to_ms(self.echo_spacing).tolist()
+        seq.sequence_type = self.sequence_type
+        header.sequenceParameters = seq
+
+        # Encoding
+        encoding = ismrmrdschema.encodingType()
+        par_imaging = ismrmrdschema.parallelImagingType()
+        par_imaging.calibrationMode = self.calibration_mode
+        par_imaging.interleavingDimension = self.interleave_dim
+        encoding.parallelImaging = par_imaging
+        encoding.trajectory = self.trajectory_type
+
+        # Encoded space
+        encoding_space = ismrmrdschema.encodingSpaceType()
+        encoding_fov = ismrmrdschema.fieldOfViewMm(
+            m_to_mm(self.encoding_fov.x), m_to_mm(self.encoding_fov.y), m_to_mm(self.encoding_fov.z)
+        )
+        encoding_matrix = ismrmrdschema.matrixSizeType(
+            self.encoding_matrix.x, self.encoding_matrix.y, self.encoding_matrix.z
+        )
+        encoding_space.matrixSize = encoding_matrix
+        encoding_space.fieldOfView_mm = encoding_fov
+
+        # Recon space
+        recon_space = ismrmrdschema.encodingSpaceType()
+        recon_fov = ismrmrdschema.fieldOfViewMm(
+            m_to_mm(self.recon_fov.x), m_to_mm(self.recon_fov.y), m_to_mm(self.recon_fov.z)
+        )
+        recon_matrix = ismrmrdschema.matrixSizeType(self.recon_matrix.x, self.recon_matrix.y, self.recon_matrix.z)
+        recon_space.matrixSize = recon_matrix
+        recon_space.fieldOfView_mm = recon_fov
+
+        # Set encoded and recon spaces
+        encoding.encodedSpace = encoding_space
+        encoding.reconSpace = recon_space
+
+        # Encoding limits
+        encoding.encodingLimits = self.encoding_limits.to_ismrmrd_encoding_limits_type()
+        header.encoding.append(encoding)
+
+        return header
+
     def __repr__(self):
         """Representation method for KHeader class."""
         te = summarize_tensorvalues(self.te)

src/mrpro/data/_kdata/KData.py

@@ -4,6 +4,7 @@
 import datetime
 import warnings
 from collections.abc import Callable
+from itertools import product
 from pathlib import Path
 
 import h5py
@@ -262,6 +263,46 @@ def from_file(
 
         return cls(kheader, kdata, ktrajectory_final)
 
+    def to_file(self, filename: str | Path) -> None:
+        """Save KData as ISMRMRD dataset to file.
+
+        Parameters
+        ----------
+        filename
+            path to the ISMRMRD file
+        """
+        # Open the dataset
+        dataset = ismrmrd.Dataset(filename, 'dataset', create_if_needed=True)
+
+        # Create ISMRMRD header
+        header = self.header.to_ismrmrd()
+        header.acquisitionSystemInformation.receiverChannels = self.data.shape[-4]
+        dataset.write_xml_header(header.toXML('utf-8'))
+
+        trajectory = self.traj.as_tensor()
+        trajectory = torch.stack(
+            [torch.broadcast_to(trajectory[i, ...], self.data[..., 0, :, :, :].shape) for i in range(3)]
+        )
+
+        # Go through data and save acquisitions
+        acq_shape = [self.data.shape[-1], self.data.shape[-4]]
+        acq = ismrmrd.Acquisition()
+        acq.resize(*acq_shape, trajectory_dimensions=3)
+        for other in product(*[range(shape) for shape in self.data.shape[:-4]]):
+            for k2 in range(self.data.shape[-3]):
+                for k1 in range(self.data.shape[-2]):
+                    acq.clear_all_flags()
+                    acq = self.header.acq_info.add_to_ismrmrd_acquisition(acq, other, k2, k1)
+
+                    # Rearrange, switch from zyx to xyz and set trajectory.
+                    acq.traj[:] = rearrange(trajectory[:, *other, k2, k1, :].numpy(), 'dim k0->k0 dim')[:, ::-1]
+
+                    # Set the data and append
+                    acq.data[:] = self.data[*other, :, k2, k1, :].numpy()
+                    dataset.append_acquisition(acq)
+
+        dataset.close()
+
     def __repr__(self):
         """Representation method for KData class."""
         traj = KTrajectory(self.traj.kz, self.traj.ky, self.traj.kx)

src/mrpro/utils/__init__.py

@@ -1,6 +1,5 @@
 import mrpro.utils.slice_profiles
 import mrpro.utils.typing
-import mrpro.utils.unit_conversion
 from mrpro.utils.smap import smap
 from mrpro.utils.remove_repeat import remove_repeat
 from mrpro.utils.zero_pad_or_crop import zero_pad_or_crop

tests/data/_IsmrmrdRawTestData.py

@@ -219,9 +219,9 @@ def __init__(
         acq.version = 1
         acq.available_channels = self.n_coils
         acq.center_sample = round(n_freq_encoding / 2)
-        acq.read_dir[0] = 1.0
-        acq.phase_dir[1] = 1.0
-        acq.slice_dir[2] = 1.0
+        acq.read_dir = (-0.33, 0.38, -0.86)
+        acq.phase_dir = (0.75, 0.66, 0.0)
+        acq.slice_dir = (-0.57, 0.65, 0.5)
 
         scan_counter = 0
 

tests/data/test_kdata.py

@@ -6,6 +6,7 @@
 from mrpro.data import KData, KTrajectory, SpatialDimension
 from mrpro.data.acq_filters import has_n_coils, is_coil_calibration_acquisition, is_image_acquisition
 from mrpro.data.AcqInfo import rearrange_acq_info_fields
+from mrpro.data.traj_calculators import KTrajectoryIsmrmrd
 from mrpro.data.traj_calculators.KTrajectoryCalculator import DummyTrajectory
 from mrpro.operators import FastFourierOp
 from mrpro.utils import split_idx
@@ -518,3 +519,27 @@ def test_modify_acq_info(random_kheader_shape):
     assert kheader.acq_info.idx.k1.shape == (n_other, n_k2, n_k1)
     assert kheader.acq_info.orientation.shape == (n_other, n_k2, n_k1, 1)
     assert kheader.acq_info.position.z.shape == (n_other, n_k2, n_k1, 1)
+
+
+def test_KData_to_file(ismrmrd_cart, tmp_path_factory):
+    """Read in data to file."""
+    kdata = KData.from_file(ismrmrd_cart.filename, DummyTrajectory())
+    # We need to make sure that the trajectory fits to the data
+    random_generator = RandomGenerator(seed=0)
+    traj = random_generator.float32_tensor(size=(3, kdata.data.shape[0], *kdata.data.shape[2:]))
+    kdata = KData(header=kdata.header, data=kdata.data, traj=KTrajectory.from_tensor(traj))
+
+    ismrmrd_filename = tmp_path_factory.mktemp('mrpro') / 'ismrmrd_saved_from_mrpro.h5'
+    kdata.to_file(ismrmrd_filename)
+    kdata_reload = KData.from_file(ismrmrd_filename, KTrajectoryIsmrmrd())
+
+    torch.testing.assert_close(kdata.data, kdata_reload.data)
+    torch.testing.assert_close(kdata.traj.as_tensor(), kdata_reload.traj.as_tensor())
+    torch.testing.assert_close(
+        kdata.header.acq_info.orientation.as_matrix(), kdata_reload.header.acq_info.orientation.as_matrix()
+    )
+    torch.testing.assert_close(kdata.header.acq_info.position.z, kdata_reload.header.acq_info.position.z)
+    torch.testing.assert_close(kdata.header.acq_info.position.y, kdata_reload.header.acq_info.position.y)
+    torch.testing.assert_close(kdata.header.acq_info.position.x, kdata_reload.header.acq_info.position.x)
+    assert kdata.header.encoding_fov.x == kdata_reload.header.encoding_fov.x
+    assert kdata.header.encoding_limits.k1.max == kdata_reload.header.encoding_limits.k1.max

tests/data/test_kheader.py

@@ -39,3 +39,16 @@ def test_kheader_verify_None(random_mandatory_ismrmrd_header, random_acq_info):
     assert torch.allclose(kheader.fa, fa_default)
     # tr is not mandatory but overwritten with None
     assert kheader.tr is tr_default
+
+
+def test_kheader_to_ismrmrd(random_mandatory_ismrmrd_header, random_acq_info):
+    """Create ISMRMRD header from KHeader."""
+    fa = [2.0, 3.0, 4.0, 5.0]
+    overwrite = {'trajectory': DummyTrajectory(), 'fa': torch.deg2rad(torch.as_tensor(fa))}
+    kheader = KHeader.from_ismrmrd(random_mandatory_ismrmrd_header, random_acq_info, overwrite=overwrite)
+    ismrmrd_header = kheader.to_ismrmrd()
+    kheader_again = KHeader.from_ismrmrd(ismrmrd_header, random_acq_info, {'trajectory': DummyTrajectory()})
+    assert ismrmrd_header.experimentalConditions.H1resonanceFrequency_Hz == kheader.lamor_frequency_proton
+    assert ismrmrd_header.encoding[0].encodedSpace.matrixSize.z == kheader.encoding_matrix.zyx[0]
+    assert ismrmrd_header.sequenceParameters.flipAngle_deg == fa
+    torch.testing.assert_close(kheader_again.fa, kheader.fa)