Add save kdata to ismrmrd file

src/mrpro/data/AcqInfo.py

@@ -12,7 +12,8 @@
 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.typing import TorchIndexerType
+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 write_to_ismrmrd_acquisition(
+        self, acquisition: ismrmrd.Acquisition, idx: TorchIndexerType
+    ) -> ismrmrd.Acquisition:
+        """Overwrite ISMRMRD acquisition information for single acquisition."""
+        acquisition.idx.kspace_encode_step_1 = self.idx.k1[idx]
+        acquisition.idx.kspace_encode_step_2 = self.idx.k2[idx]
+        acquisition.idx.average = self.idx.average[idx]
+        acquisition.idx.slice = self.idx.slice[idx]
+        acquisition.idx.contrast = self.idx.contrast[idx]
+        acquisition.idx.phase = self.idx.phase[idx]
+        acquisition.idx.repetition = self.idx.repetition[idx]
+        acquisition.idx.set = self.idx.set[idx]
+        acquisition.idx.segment = self.idx.segment[idx]
+        acquisition.idx.user = (
+            self.idx.user0[idx],
+            self.idx.user1[idx],
+            self.idx.user2[idx],
+            self.idx.user3[idx],
+            self.idx.user4[idx],
+            self.idx.user5[idx],
+            self.idx.user6[idx],
+            self.idx.user7[idx],
+        )
+
+        # active_channesl, number_of_samples and trajectory_dimensions are read-only and cannot be set
+        acquisition.acquisition_time_stamp = self.acquisition_time_stamp[idx].squeeze(-1)
+        acquisition.available_channels = self.available_channels[idx].squeeze(-1)
+        acquisition.center_sample = self.center_sample[idx].squeeze(-1)
+        acquisition.channel_mask = tuple(self.channel_mask[idx])
+        acquisition.discard_post = self.discard_post[idx].squeeze(-1)
+        acquisition.discard_pre = self.discard_pre[idx].squeeze(-1)
+        acquisition.encoding_space_ref = self.encoding_space_ref[idx].squeeze(-1)
+        acquisition.measurement_uid = self.measurement_uid[idx].squeeze(-1)
+        acquisition.patient_table_position = self.patient_table_position[idx][0].apply(m_to_mm).zyx[::-1]  # zyx -> xyz
+        directions = self.orientation[idx][0].as_directions()
+        acquisition.slice_dir = directions[0].zyx[::-1]  # zyx -> xyz
+        acquisition.phase_dir = directions[1].zyx[::-1]
+        acquisition.read_dir = directions[2].zyx[::-1]
+
+        acquisition.physiology_time_stamp = tuple(self.physiology_time_stamp[idx])
+        acquisition.position = self.position[idx][0].apply(m_to_mm).zyx[::-1]
+        acquisition.sample_time_us = self.sample_time_us[idx].squeeze(-1)
+        acquisition.scan_counter = self.scan_counter[idx].squeeze(-1)
+        acquisition.user_float = tuple(self.user_float[idx])
+        acquisition.user_int = tuple(self.user_int[idx])
+        acquisition.version = self.version[idx].squeeze(-1)
+        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/KData.py

@@ -265,6 +265,47 @@ 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(-1)
+        trajectory_shape = (*self.data[..., 0, :, :, :].shape, 3)  # no coils, 3d
+        trajectory = torch.broadcast_to(trajectory, trajectory_shape)
+
+        # 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 np.ndindex(self.data.shape[:-4]):
+            for k2 in range(self.data.shape[-3]):
+                for k1 in range(self.data.shape[-2]):
+                    acq.clear_all_flags()
+                    self.header.acq_info.write_to_ismrmrd_acquisition(acq, (*other, k2, k1))
+
+                    # Rearrange, switch from zyx to xyz and set trajectory.
+                    acq.traj[:] = trajectory[(*other, k2, k1)].numpy()[:, ::-1]  # zyx -> xyz
+
+                    # Set the data and append
+                    idx = other + (slice(None),) + (k2, k1)  # python 3.10 and mypy #noqa: RUF005
+                    acq.data[:] = self.data[idx].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/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 self.tr is not None:
+            seq.TR = s_to_ms(self.tr).tolist()
+        if self.te is not None:
+            seq.TE = s_to_ms(self.te).tolist()
+        if self.ti is not None:
+            seq.TI = s_to_ms(self.ti).tolist()
+        if self.fa is not None:
+            seq.flipAngle_deg = torch.rad2deg(self.fa).tolist()
+        if self.echo_spacing is not None:
+            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/SpatialDimension.py

@@ -118,6 +118,7 @@ def apply_(self, function: Callable[[T], T] | None = None, **_) -> Self:
         function
             function to apply
         """
+        # This function is mainly used for type hinting
         return super(SpatialDimension, self).apply_(function)
 
     # This function is mainly for type hinting and docstring

tests/data/_IsmrmrdRawTestData.py

@@ -221,9 +221,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
@@ -581,3 +582,27 @@ def test_KData_compress_coils_error_n_coils(consistently_shaped_kdata):
     existing_coils = consistently_shaped_kdata.data.shape[-4]
     with pytest.raises(ValueError, match='greater'):
         consistently_shaped_kdata.compress_coils(existing_coils + 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)