From c202ee4fa06e0cde4a6816fd16432529e1c21ca9 Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Sat, 28 Sep 2024 21:00:55 +0200 Subject: [PATCH 01/19] write to ismrmrd almost finished --- src/mrpro/data/AcqInfo.py | 60 ++++++++++ src/mrpro/data/EncodingLimits.py | 15 +++ src/mrpro/data/KHeader.py | 188 ++++++++++++++++++++++++++++++- src/mrpro/data/_kdata/KData.py | 41 +++++++ tests/data/test_kdata.py | 19 ++++ tests/data/test_kheader.py | 13 +++ 6 files changed, 331 insertions(+), 5 deletions(-) diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index 525c164e0..143ec6801 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -25,6 +25,16 @@ def mm_to_m(m: T) -> T: return m / 1000 +def s_to_ms(ms: T) -> T: + """Convert s to ms.""" + return ms * 1000 + + +def m_to_mm(m: T) -> T: + """Convert m to mm.""" + return m * 1000 + + @dataclass(slots=True) class AcqIdx(MoveDataMixin): """Acquisition index for each readout.""" @@ -263,3 +273,53 @@ def spatialdimension_2d( version=tensor_2d(headers['version']), ) return acq_info + + def add_ismrmrd_acquisition_info( + self, acquisition: ismrmrd.Acquisition, other: int, 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 = tuple( + [m_to_mm(getattr(self.patient_table_position, label)[other, k2, k1, 0]) for label in ['x', 'y', 'z']] + ) + acquisition.phase_dir = tuple([getattr(self.phase_dir, label)[other, k2, k1, 0] for label in ['x', 'y', 'z']]) + acquisition.physiology_time_stamp = tuple(self.physiology_time_stamp[other, k2, k1, :]) + acquisition.position = tuple( + [m_to_mm(getattr(self.position, label)[other, k2, k1, 0]) for label in ['x', 'y', 'z']] + ) + acquisition.read_dir = tuple([getattr(self.read_dir, label)[other, k2, k1, 0] for label in ['x', 'y', 'z']]) + 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 = tuple([getattr(self.slice_dir, label)[other, k2, k1, 0] for label in ['x', 'y', 'z']]) + 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 diff --git a/src/mrpro/data/EncodingLimits.py b/src/mrpro/data/EncodingLimits.py index 546cff883..249eabb1d 100644 --- a/src/mrpro/data/EncodingLimits.py +++ b/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) diff --git a/src/mrpro/data/KHeader.py b/src/mrpro/data/KHeader.py index 0cf5fcc87..3eaca323c 100644 --- a/src/mrpro/data/KHeader.py +++ b/src/mrpro/data/KHeader.py @@ -12,7 +12,7 @@ import torch from mrpro.data import enums -from mrpro.data.AcqInfo import AcqInfo, mm_to_m, ms_to_s +from mrpro.data.AcqInfo import AcqInfo, m_to_mm, mm_to_m, ms_to_s, s_to_ms from mrpro.data.EncodingLimits import EncodingLimits from mrpro.data.MoveDataMixin import MoveDataMixin from mrpro.data.SpatialDimension import SpatialDimension @@ -171,13 +171,13 @@ def from_ismrmrd( parameters['n_coils'] = header.acquisitionSystemInformation.receiverChannels if header.sequenceParameters is not None: - if header.sequenceParameters.TR: + if any(header.sequenceParameters.TR): parameters['tr'] = ms_to_s(torch.as_tensor(header.sequenceParameters.TR)) - if header.sequenceParameters.TE: + if any(header.sequenceParameters.TE): parameters['te'] = ms_to_s(torch.as_tensor(header.sequenceParameters.TE)) - if header.sequenceParameters.TI: + if any(header.sequenceParameters.TI): parameters['ti'] = ms_to_s(torch.as_tensor(header.sequenceParameters.TI)) - if header.sequenceParameters.flipAngle_deg: + if any(header.sequenceParameters.flipAngle_deg): parameters['fa'] = torch.deg2rad(torch.as_tensor(header.sequenceParameters.flipAngle_deg)) if header.sequenceParameters.echo_spacing: parameters['echo_spacing'] = ms_to_s(torch.as_tensor(header.sequenceParameters.echo_spacing)) @@ -210,6 +210,8 @@ def from_ismrmrd( if enc.trajectory is not None: parameters['traj_type'] = enums.TrajectoryType(enc.trajectory.value) + if enc.trajectoryDescription is not None: + parameters['trajectory_description'] = TrajectoryDescription.from_ismrmrd(enc.trajectoryDescription) # Either use the series or study time if available if header.measurementInformation is not None and header.measurementInformation.seriesTime is not None: @@ -251,6 +253,8 @@ def from_ismrmrd( # Dump everything into misc parameters['misc'] = dataclasses.asdict(header) + # Remember encoding number + parameters['misc']['encoding_number'] = encoding_number if overwrite is not None: parameters.update(overwrite) @@ -270,6 +274,180 @@ 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() + header.version = self.misc['version'] + + def create_from_misc_dictionary(xsd_type, misc_dictionary_entry: dict): # noqa: ANN001 + return xsd_type(**misc_dictionary_entry) if misc_dictionary_entry is not None else xsd_type() + + # Study information + if self.misc['studyInformation'] is not None: + study = ismrmrdschema.studyInformationType(**self.misc['studyInformation']) + header.studyInformation = study + + # Experimental conditions + exp = create_from_misc_dictionary(ismrmrdschema.experimentalConditionsType, self.misc['experimentalConditions']) + exp.H1resonanceFrequency_Hz = self.h1_freq + header.experimentalConditions = exp + + # Subject information + subj = create_from_misc_dictionary(ismrmrdschema.subjectInformationType, self.misc['subjectInformation']) + subj.patientName = self.patient_name + header.subjectInformation = subj + + # Measurement information + meas = create_from_misc_dictionary( + ismrmrdschema.measurementInformationType, self.misc['measurementInformation'] + ) + if self.misc['measurementInformation'] is not None: + if self.misc['measurementInformation']['measurementDependency'] is not None: + meas.measurementDependency = [ + ismrmrdschema.measurementDependencyType(**meas_dependency) + for meas_dependency in self.misc['measurementInformation']['measurementDependency'] + ] + if self.misc['measurementInformation']['referencedImageSequence'] is not None: + meas.referencedImageSequence = [ + ismrmrdschema.measurementDependencyType(**ref_image_sequence) + for ref_image_sequence in self.misc['measurementInformation']['referencedImageSequence'] + ] + + 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 = create_from_misc_dictionary( + ismrmrdschema.acquisitionSystemInformationType, self.misc['acquisitionSystemInformation'] + ) + if ( + self.misc['acquisitionSystemInformation'] is not None + and self.misc['acquisitionSystemInformation']['coilLabel'] is not None + ): + sys.coilLabel = [ + ismrmrdschema.coilLabelType(**coil_label) + for coil_label in self.misc['acquisitionSystemInformation']['coilLabel'] + ] + if self.n_coils is not None: + sys.receiverChannels = self.n_coils + sys.systemFieldStrength_T = self.b0 + sys.systemModel = self.model + sys.systemVendor = self.vendor + header.acquisitionSystemInformation = sys + + # Sequence information + seq = create_from_misc_dictionary(ismrmrdschema.sequenceParametersType, self.misc['sequenceParameters']) + if self.misc['sequenceParameters'] is not None: + if self.misc['sequenceParameters']['diffusionDimension'] is not None: + seq.diffusion = ismrmrdschema.diffusionDimensionType( + **self.misc['sequenceParameters']['diffusionDimension'] + ) + if self.misc['sequenceParameters']['diffusion'] is not None: + seq.diffusion = [ + ismrmrdschema.diffusionTypeType(**diffusion) + for diffusion in self.misc['sequenceParameters']['diffusion'] + ] + + 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.seq_type + header.sequenceParameters = seq + + # Waveform information + wave = ismrmrdschema.waveformInformationType() + if len(self.misc['waveformInformation']): + wave = [ismrmrdschema.waveformInformationType(**wave) for wave in self.misc['waveformInformation']] + header.waveformInformation = wave + + # Encoding + encoding_number = self.misc['encoding_number'] + encoding = create_from_misc_dictionary(ismrmrdschema.encodingType, self.misc['encoding'][encoding_number]) + encoding.echoTrainLength = self.echo_train_length + par_imaging = create_from_misc_dictionary( + ismrmrdschema.parallelImagingType, self.misc['encoding'][encoding_number]['parallelImaging'] + ) + if self.misc['encoding'][encoding_number]['parallelImaging']: + par_imaging.accelerationFactor = create_from_misc_dictionary( + ismrmrdschema.accelerationFactorType, + self.misc['encoding'][encoding_number]['parallelImaging']['accelerationFactor'], + ) + par_imaging.calibrationMode = self.calibration_mode + par_imaging.interleavingDimension = self.interleave_dim + encoding.parallelImaging = par_imaging + encoding.trajectory = self.traj_type + # encoding.trajectoryDescription = self.trajectory_description.to_ismrmrd() + + # 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) + + # User parameters + if self.misc['userParameters']: + user = ismrmrdschema.userParametersType() + user.userParameterLong = [ + ismrmrdschema.userParameterLongType(**par) for par in self.misc['userParameters']['userParameterLong'] + ] + user.userParameterDouble = [ + ismrmrdschema.userParameterDoubleType(**par) + for par in self.misc['userParameters']['userParameterDouble'] + ] + user.userParameterString = [ + ismrmrdschema.userParameterStringType(**par) + for par in self.misc['userParameters']['userParameterString'] + ] + user.userParameterBase64 = [ + ismrmrdschema.userParameterBase64Type(**par) + for par in self.misc['userParameters']['userParameterBase64'] + ] + header.userParameters = user + + return header + def __repr__(self): """Representation method for KHeader class.""" te = summarize_tensorvalues(self.te) diff --git a/src/mrpro/data/_kdata/KData.py b/src/mrpro/data/_kdata/KData.py index 81b4f9003..df835b125 100644 --- a/src/mrpro/data/_kdata/KData.py +++ b/src/mrpro/data/_kdata/KData.py @@ -235,6 +235,47 @@ def sort_and_reshape_tensor_fields(input_tensor: torch.Tensor): 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 + """ + if self.data.ndim != 5: + raise ValueError('Only KData with exactly 5 dimensions (other, coil, k2, k1, k0) can be saved.') + + # Open the dataset + dataset = ismrmrd.Dataset(filename, 'dataset', create_if_needed=True) + + # Create ISMRMRD header + header = self.header.to_ismrmrd() + 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]] + for other in range(self.data.shape[-5]): + for k2 in range(self.data.shape[-3]): + for k1 in range(self.data.shape[-2]): + acq = ismrmrd.Acquisition() + acq.resize(*acq_shape, trajectory_dimensions=3) + acq = self.header.acq_info.add_ismrmrd_acquisition_info(acq, other, k2, k1) + + # Rearrange, switch from zyx to xz 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) diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index 68f9dfd07..f785dfa8e 100644 --- a/tests/data/test_kdata.py +++ b/tests/data/test_kdata.py @@ -5,6 +5,7 @@ from einops import rearrange, repeat from mrpro.data import KData, KTrajectory, SpatialDimension from mrpro.data.acq_filters import is_coil_calibration_acquisition +from mrpro.data.traj_calculators import KTrajectoryIsmrmrd from mrpro.data.traj_calculators.KTrajectoryCalculator import DummyTrajectory from mrpro.operators import FastFourierOp from mrpro.utils import modify_acq_info, split_idx @@ -469,3 +470,21 @@ def test_KData_remove_readout_os(monkeypatch, random_kheader): # testing functions such as numpy.testing.assert_almost_equal fails because there are few voxels with high # differences along the edges of the elliptic objects. assert relative_image_difference(torch.abs(img_recon), img_tensor[:, 0, ...]) <= 0.05 + + +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()) + 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 diff --git a/tests/data/test_kheader.py b/tests/data/test_kheader.py index 55cbfdf67..4395eb86b 100644 --- a/tests/data/test_kheader.py +++ b/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.acquisitionSystemInformation.systemFieldStrength_T == kheader.b0 + 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) From daf40fafcd04228437164d4e1f7d44f12c4808a9 Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Mon, 30 Sep 2024 17:32:37 +0200 Subject: [PATCH 02/19] current tests moved --- src/mrpro/data/AcqInfo.py | 16 ++-------------- src/mrpro/data/KHeader.py | 3 ++- src/mrpro/utils/__init__.py | 2 ++ 3 files changed, 6 insertions(+), 15 deletions(-) diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index 525c164e0..98ec926ea 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -2,7 +2,7 @@ from collections.abc import Callable, Sequence from dataclasses import dataclass -from typing import Self, TypeVar +from typing import Self import ismrmrd import numpy as np @@ -10,19 +10,7 @@ from mrpro.data.MoveDataMixin import MoveDataMixin from mrpro.data.SpatialDimension import SpatialDimension - -# Conversion functions for units -T = TypeVar('T', float, torch.Tensor) - - -def ms_to_s(ms: T) -> T: - """Convert ms to s.""" - return ms / 1000 - - -def mm_to_m(m: T) -> T: - """Convert mm to m.""" - return m / 1000 +from mrpro.utils.unit_conversion import mm_to_m @dataclass(slots=True) diff --git a/src/mrpro/data/KHeader.py b/src/mrpro/data/KHeader.py index 0cf5fcc87..966da772c 100644 --- a/src/mrpro/data/KHeader.py +++ b/src/mrpro/data/KHeader.py @@ -12,12 +12,13 @@ import torch from mrpro.data import enums -from mrpro.data.AcqInfo import AcqInfo, mm_to_m, ms_to_s +from mrpro.data.AcqInfo import AcqInfo from mrpro.data.EncodingLimits import EncodingLimits from mrpro.data.MoveDataMixin import MoveDataMixin from mrpro.data.SpatialDimension import SpatialDimension from mrpro.data.TrajectoryDescription import TrajectoryDescription from mrpro.utils.summarize_tensorvalues import summarize_tensorvalues +from mrpro.utils.unit_conversion import mm_to_m, ms_to_s if TYPE_CHECKING: # avoid circular imports by importing only when type checking diff --git a/src/mrpro/utils/__init__.py b/src/mrpro/utils/__init__.py index 5c02c58f2..fc5e3ec2e 100644 --- a/src/mrpro/utils/__init__.py +++ b/src/mrpro/utils/__init__.py @@ -5,3 +5,5 @@ from mrpro.utils.split_idx import split_idx from mrpro.utils.Rotation import Rotation import mrpro.utils.slice_profiles +import mrpro.utils.unit_conversion + From 27eba4720397753b32fc2a682d6c9b0e477322dd Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Tue, 1 Oct 2024 11:38:42 +0200 Subject: [PATCH 03/19] unit conversion finished --- src/mrpro/utils/unit_conversion.py | 81 ++++++++++++++++++++++++++++ tests/utils/test_unit_conversion.py | 82 +++++++++++++++++++++++++++++ 2 files changed, 163 insertions(+) create mode 100644 src/mrpro/utils/unit_conversion.py create mode 100644 tests/utils/test_unit_conversion.py diff --git a/src/mrpro/utils/unit_conversion.py b/src/mrpro/utils/unit_conversion.py new file mode 100644 index 000000000..274dc71c8 --- /dev/null +++ b/src/mrpro/utils/unit_conversion.py @@ -0,0 +1,81 @@ +"""Conversion between different units.""" + +from typing import TypeVar + +import numpy as np +import torch + +GYROMAGNETIC_MOMENT_PROTON = 42.58 * 1e6 + +# Conversion functions for units +T = TypeVar('T', float, torch.Tensor) + + +def ms_to_s(ms: T) -> T: + """Convert ms to s.""" + return ms / 1000 + + +def s_to_ms(s: T) -> T: + """Convert s to ms.""" + return s * 1000 + + +def mm_to_m(mm: T) -> T: + """Convert mm to m.""" + return mm / 1000 + + +def m_to_mm(m: T) -> T: + """Convert m to mm.""" + return m * 1000 + + +def deg_to_rad(deg: T) -> T: + """Convert degree to radians.""" + if isinstance(deg, torch.Tensor): + return torch.deg2rad(deg) + return deg / 180.0 * np.pi + + +def rad_to_deg(deg: T) -> T: + """Convert radians to degree.""" + if isinstance(deg, torch.Tensor): + return torch.rad2deg(deg) + return deg * 180.0 / np.pi + + +def lamor_frequency_to_magnetic_field(lamor_frequency: T, gyromagnetic_ratio: float = GYROMAGNETIC_MOMENT_PROTON) -> T: + """Convert the Lamor frequency [Hz] to the magntic field strength [T]. + + Parameters + ---------- + lamor_frequency + Lamor frequency [Hz] + gyromagnetic_ratio + Gyromagnetic ratio [Hz/T], default: gyromagnetic ratio of 1H proton + + Returns + ------- + Magnetic field strength [T] + """ + return lamor_frequency / gyromagnetic_ratio + + +def magnetic_field_to_lamor_frequency( + magnetic_field_strength: T, gyromagnetic_ratio: float = GYROMAGNETIC_MOMENT_PROTON +) -> T: + """Convert the magntic field strength [T] to Lamor frequency [Hz]. + + Parameters + ---------- + magnetic_field_strength + Strength of the magnetic field [T] + gyromagnetic_ratio + Gyromagnetic ratio [Hz/T], default: gyromagnetic ratio of 1H proton + + Returns + ------- + Lamor frequency [Hz] + """ + return magnetic_field_strength * gyromagnetic_ratio diff --git a/tests/utils/test_unit_conversion.py b/tests/utils/test_unit_conversion.py new file mode 100644 index 000000000..a232a5366 --- /dev/null +++ b/tests/utils/test_unit_conversion.py @@ -0,0 +1,82 @@ +"""Tests of unit conversion.""" + +import numpy as np +import torch +from mrpro.utils.unit_conversion import ( + deg_to_rad, + lamor_frequency_to_magnetic_field, + m_to_mm, + magnetic_field_to_lamor_frequency, + mm_to_m, + ms_to_s, + rad_to_deg, + s_to_ms, +) + +from tests import RandomGenerator + + +def test_mm_to_m(): + """Verify mm to m conversion.""" + generator = RandomGenerator(seed=0) + mm_input = generator.float32_tensor((3, 4, 5)) + torch.testing.assert_close(mm_to_m(mm_input), mm_input / 1000.0) + + +def test_m_to_mm(): + """Verify m to mm conversion.""" + generator = RandomGenerator(seed=0) + m_input = generator.float32_tensor((3, 4, 5)) + torch.testing.assert_close(m_to_mm(m_input), m_input * 1000.0) + + +def test_ms_to_s(): + """Verify ms to s conversion.""" + generator = RandomGenerator(seed=0) + ms_input = generator.float32_tensor((3, 4, 5)) + torch.testing.assert_close(ms_to_s(ms_input), ms_input / 1000.0) + + +def test_s_to_ms(): + """Verify s to ms conversion.""" + generator = RandomGenerator(seed=0) + s_input = generator.float32_tensor((3, 4, 5)) + torch.testing.assert_close(s_to_ms(s_input), s_input * 1000.0) + + +def test_rad_to_deg_tensor(): + """Verify radians to degree conversion.""" + generator = RandomGenerator(seed=0) + s_input = generator.float32_tensor((3, 4, 5)) + torch.testing.assert_close(rad_to_deg(s_input), torch.rad2deg(s_input)) + + +def test_deg_to_rad_tensor(): + """Verify degree to radians conversion.""" + generator = RandomGenerator(seed=0) + s_input = generator.float32_tensor((3, 4, 5)) + torch.testing.assert_close(deg_to_rad(s_input), torch.deg2rad(s_input)) + + +def test_rad_to_deg_float(): + """Verify radians to degree conversion.""" + assert rad_to_deg(np.pi / 2) == 90.0 + + +def test_deg_to_rad_float(): + """Verify degree to radians conversion.""" + assert deg_to_rad(180.0) == np.pi + + +def test_lamor_frequency_to_magnetic_field(): + """Verify conversion of lamor frequency to magnetic field.""" + proton_gyromagnetic_ratio = 42.58 * 1e6 + proton_lamor_frequency_at_3tesla = 127.74 * 1e6 + assert lamor_frequency_to_magnetic_field(proton_lamor_frequency_at_3tesla, proton_gyromagnetic_ratio) == 3.0 + + +def test_magnetic_field_to_lamor_frequency(): + """Verify conversion of magnetic field to lamor frequency.""" + proton_gyromagnetic_ratio = 42.58 * 1e6 + magnetic_field_strength = 3.0 + assert magnetic_field_to_lamor_frequency(magnetic_field_strength, proton_gyromagnetic_ratio) == 127.74 * 1e6 From bf722c25b5ad5d7207dc4842ed5308570d1163cd Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Tue, 1 Oct 2024 11:56:09 +0200 Subject: [PATCH 04/19] kheader clean-up --- src/mrpro/data/KHeader.py | 38 ++++++++++++-------------------------- 1 file changed, 12 insertions(+), 26 deletions(-) diff --git a/src/mrpro/data/KHeader.py b/src/mrpro/data/KHeader.py index 966da772c..123036c3f 100644 --- a/src/mrpro/data/KHeader.py +++ b/src/mrpro/data/KHeader.py @@ -16,7 +16,6 @@ from mrpro.data.EncodingLimits import EncodingLimits from mrpro.data.MoveDataMixin import MoveDataMixin from mrpro.data.SpatialDimension import SpatialDimension -from mrpro.data.TrajectoryDescription import TrajectoryDescription from mrpro.utils.summarize_tensorvalues import summarize_tensorvalues from mrpro.utils.unit_conversion import mm_to_m, ms_to_s @@ -40,9 +39,6 @@ class KHeader(MoveDataMixin): trajectory: KTrajectoryCalculator """Function to calculate the k-space trajectory.""" - b0: float - """Magnetic field strength [T].""" - encoding_limits: EncodingLimits """K-space encoding limits.""" @@ -61,12 +57,9 @@ class KHeader(MoveDataMixin): acq_info: AcqInfo """Information of the acquisitions (i.e. readout lines).""" - h1_freq: float + lamor_frequency_proton: float """Lamor frequency of hydrogen nuclei [Hz].""" - n_coils: int | None = None - """Number of receiver coils.""" - datetime: datetime.datetime | None = None """Date and time of acquisition.""" @@ -88,7 +81,7 @@ class KHeader(MoveDataMixin): echo_train_length: int = 1 """Number of echoes in a multi-echo acquisition.""" - seq_type: str = UNKNOWN + sequence_type: str = UNKNOWN """Type of sequence.""" model: str = UNKNOWN @@ -100,16 +93,13 @@ class KHeader(MoveDataMixin): protocol_name: str = UNKNOWN """Name of the acquisition protocol.""" - misc: dict = dataclasses.field(default_factory=dict) # do not use {} here! - """Dictionary with miscellaneous parameters.""" - calibration_mode: enums.CalibrationMode = enums.CalibrationMode.OTHER """Mode of how calibration data is acquired. """ interleave_dim: enums.InterleavingDimension = enums.InterleavingDimension.OTHER """Interleaving dimension.""" - traj_type: enums.TrajectoryType = enums.TrajectoryType.OTHER + trajectory_type: enums.TrajectoryType = enums.TrajectoryType.OTHER """Type of trajectory.""" measurement_id: str = UNKNOWN @@ -118,8 +108,8 @@ class KHeader(MoveDataMixin): patient_name: str = UNKNOWN """Name of the patient.""" - trajectory_description: TrajectoryDescription = dataclasses.field(default_factory=TrajectoryDescription) - """Description of the trajectory.""" + _misc: dict = dataclasses.field(default_factory=dict) # do not use {} here! + """Dictionary with miscellaneous parameters.""" @property def fa_degree(self) -> torch.Tensor | None: @@ -160,7 +150,10 @@ def from_ismrmrd( enc: ismrmrdschema.encodingType = header.encoding[encoding_number] # These are guaranteed to exist - parameters = {'h1_freq': header.experimentalConditions.H1resonanceFrequency_Hz, 'acq_info': acq_info} + parameters = { + 'lamor_frequency_proton': header.experimentalConditions.H1resonanceFrequency_Hz, + 'acq_info': acq_info, + } if defaults is not None: parameters.update(defaults) @@ -184,7 +177,7 @@ def from_ismrmrd( parameters['echo_spacing'] = ms_to_s(torch.as_tensor(header.sequenceParameters.echo_spacing)) if header.sequenceParameters.sequence_type is not None: - parameters['seq_type'] = header.sequenceParameters.sequence_type + parameters['sequence_type'] = header.sequenceParameters.sequence_type if enc.reconSpace is not None: parameters['recon_fov'] = SpatialDimension[float].from_xyz(enc.reconSpace.fieldOfView_mm, mm_to_m) @@ -210,7 +203,7 @@ def from_ismrmrd( ) if enc.trajectory is not None: - parameters['traj_type'] = enums.TrajectoryType(enc.trajectory.value) + parameters['trajectory_type'] = enums.TrajectoryType(enc.trajectory.value) # Either use the series or study time if available if header.measurementInformation is not None and header.measurementInformation.seriesTime is not None: @@ -243,15 +236,8 @@ def from_ismrmrd( if header.acquisitionSystemInformation.systemModel is not None: parameters['model'] = header.acquisitionSystemInformation.systemModel - if header.acquisitionSystemInformation.systemFieldStrength_T is not None: - parameters['b0'] = header.acquisitionSystemInformation.systemFieldStrength_T - - # estimate b0 from h1_freq if not given - if 'b0' not in parameters: - parameters['b0'] = parameters['h1_freq'] / 4258e4 - # Dump everything into misc - parameters['misc'] = dataclasses.asdict(header) + parameters['_misc'] = dataclasses.asdict(header) if overwrite is not None: parameters.update(overwrite) From 4c4b7fa4b77bb7185b505200aba91f84073915f4 Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Tue, 1 Oct 2024 13:36:19 +0200 Subject: [PATCH 05/19] n_coils removed --- src/mrpro/data/KHeader.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/mrpro/data/KHeader.py b/src/mrpro/data/KHeader.py index 123036c3f..1be529494 100644 --- a/src/mrpro/data/KHeader.py +++ b/src/mrpro/data/KHeader.py @@ -158,12 +158,6 @@ def from_ismrmrd( if defaults is not None: parameters.update(defaults) - if ( - header.acquisitionSystemInformation is not None - and header.acquisitionSystemInformation.receiverChannels is not None - ): - parameters['n_coils'] = header.acquisitionSystemInformation.receiverChannels - if header.sequenceParameters is not None: if header.sequenceParameters.TR: parameters['tr'] = ms_to_s(torch.as_tensor(header.sequenceParameters.TR)) From ccd93f704705523d81f7f62b8c3b29e85e317e76 Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Tue, 1 Oct 2024 16:17:34 +0200 Subject: [PATCH 06/19] orientation as rotation --- src/mrpro/data/AcqInfo.py | 38 ++++++++++++++++++++++++++------------ tests/data/test_kdata.py | 2 +- 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index 98ec926ea..37c8f5605 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -9,6 +9,7 @@ import torch 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 @@ -109,30 +110,24 @@ class AcqInfo(MoveDataMixin): number_of_samples: torch.Tensor """Number of sample points per readout (readouts may have different number of sample points).""" + orientation: Rotation + """Rotation describing the orientation of the readout, phase and slice encoding direction.""" + patient_table_position: SpatialDimension[torch.Tensor] """Offset position of the patient table, in LPS coordinates [m].""" - phase_dir: SpatialDimension[torch.Tensor] - """Directional cosine of phase encoding (2D).""" - physiology_time_stamp: torch.Tensor """Time stamps relative to physiological triggering, e.g. ECG. Not in s but in vendor-specific time units""" position: SpatialDimension[torch.Tensor] """Center of the excited volume, in LPS coordinates relative to isocenter [m].""" - read_dir: SpatialDimension[torch.Tensor] - """Directional cosine of readout/frequency encoding.""" - sample_time_us: torch.Tensor """Readout bandwidth, as time between samples [us].""" scan_counter: torch.Tensor """Zero-indexed incrementing counter for readouts.""" - slice_dir: SpatialDimension[torch.Tensor] - """Directional cosine of slice normal, i.e. cross-product of read_dir and phase_dir.""" - trajectory_dimensions: torch.Tensor # =3. We only support 3D Trajectories: kz always exists. """Dimensionality of the k-space trajectory vector.""" @@ -224,6 +219,23 @@ def spatialdimension_2d( user7=tensor(idx['user'][:, 7]), ) + # Calculate orientation as rotation matrix from directional cosines + def orientation_from_directional_cosine( + slice_dir: SpatialDimension[torch.Tensor], + phase_dir: SpatialDimension[torch.Tensor], + read_dir: SpatialDimension[torch.Tensor], + ) -> Rotation: + return Rotation.from_matrix( + torch.stack( + ( + torch.stack((slice_dir.z, slice_dir.y, slice_dir.x), dim=-1), + torch.stack((phase_dir.z, phase_dir.y, phase_dir.x), dim=-1), + torch.stack((read_dir.z, read_dir.y, read_dir.x), dim=-1), + ), + dim=-2, + ) + ) + acq_info = cls( idx=acq_idx, acquisition_time_stamp=tensor_2d(headers['acquisition_time_stamp']), @@ -237,14 +249,16 @@ def spatialdimension_2d( flags=tensor_2d(headers['flags']), measurement_uid=tensor_2d(headers['measurement_uid']), number_of_samples=tensor_2d(headers['number_of_samples']), + orientation=orientation_from_directional_cosine( + spatialdimension_2d(headers['slice_dir']), + spatialdimension_2d(headers['phase_dir']), + spatialdimension_2d(headers['read_dir']), + ), patient_table_position=spatialdimension_2d(headers['patient_table_position'], mm_to_m), - phase_dir=spatialdimension_2d(headers['phase_dir']), physiology_time_stamp=tensor_2d(headers['physiology_time_stamp']), position=spatialdimension_2d(headers['position'], mm_to_m), - read_dir=spatialdimension_2d(headers['read_dir']), sample_time_us=tensor_2d(headers['sample_time_us']), scan_counter=tensor_2d(headers['scan_counter']), - slice_dir=spatialdimension_2d(headers['slice_dir']), trajectory_dimensions=tensor_2d(headers['trajectory_dimensions']).fill_(3), # see above user_float=tensor_2d(headers['user_float']), user_int=tensor_2d(headers['user_int']), diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index 68f9dfd07..9c1bfd9d8 100644 --- a/tests/data/test_kdata.py +++ b/tests/data/test_kdata.py @@ -162,7 +162,7 @@ def test_KData_kspace(ismrmrd_cart): assert relative_image_difference(reconstructed_img[0, 0, 0, ...], ismrmrd_cart.img_ref) <= 0.05 -@pytest.mark.parametrize(('field', 'value'), [('b0', 11.3), ('tr', torch.tensor([24.3]))]) +@pytest.mark.parametrize(('field', 'value'), [('lamor_frequency_proton', 42.88 * 1e6), ('tr', torch.tensor([24.3]))]) def test_KData_modify_header(ismrmrd_cart, field, value): """Overwrite some parameters in the header.""" parameter_dict = {field: value} From 9c3d61e997be25eddaf79af19b5f44029cf8c3ea Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Tue, 1 Oct 2024 16:22:22 +0200 Subject: [PATCH 07/19] trajectory description removed --- src/mrpro/data/TrajectoryDescription.py | 29 ------------------------- src/mrpro/data/__init__.py | 1 - 2 files changed, 30 deletions(-) delete mode 100644 src/mrpro/data/TrajectoryDescription.py diff --git a/src/mrpro/data/TrajectoryDescription.py b/src/mrpro/data/TrajectoryDescription.py deleted file mode 100644 index f1a9cce5a..000000000 --- a/src/mrpro/data/TrajectoryDescription.py +++ /dev/null @@ -1,29 +0,0 @@ -"""TrajectoryDescription dataclass.""" - -import dataclasses -from dataclasses import dataclass -from typing import Self - -from ismrmrd.xsd.ismrmrdschema.ismrmrd import trajectoryDescriptionType - - -@dataclass(slots=True) -class TrajectoryDescription: - """TrajectoryDescription dataclass.""" - - identifier: str = '' - user_parameter_long: dict[str, int] = dataclasses.field(default_factory=dict) - user_parameter_double: dict[str, float] = dataclasses.field(default_factory=dict) - user_parameter_string: dict[str, str] = dataclasses.field(default_factory=dict) - comment: str = '' - - @classmethod - def from_ismrmrd(cls, trajectory_description: trajectoryDescriptionType) -> Self: - """Create TrajectoryDescription from ismrmrd traj description.""" - return cls( - user_parameter_long={p.name: int(p.value) for p in trajectory_description.userParameterLong}, - user_parameter_double={p.name: float(p.value) for p in trajectory_description.userParameterDouble}, - user_parameter_string={p.name: str(p.value) for p in trajectory_description.userParameterString}, - comment=trajectory_description.comment or '', - identifier=trajectory_description.identifier or '', - ) diff --git a/src/mrpro/data/__init__.py b/src/mrpro/data/__init__.py index 06ef1d76e..71b35c8e0 100644 --- a/src/mrpro/data/__init__.py +++ b/src/mrpro/data/__init__.py @@ -16,4 +16,3 @@ from mrpro.data.QHeader import QHeader from mrpro.data.Rotation import Rotation from mrpro.data.SpatialDimension import SpatialDimension -from mrpro.data.TrajectoryDescription import TrajectoryDescription From c3d418760c12e6b480ef3e83faedc3f2ac8528cb Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Wed, 2 Oct 2024 16:51:04 +0200 Subject: [PATCH 08/19] _apply_ instead of modify_acq_info --- src/mrpro/data/AcqInfo.py | 39 +++++++++++++++++++- src/mrpro/data/_kdata/KData.py | 6 +-- src/mrpro/data/_kdata/KDataRearrangeMixin.py | 10 ++--- src/mrpro/data/_kdata/KDataSelectMixin.py | 16 +++++--- src/mrpro/data/_kdata/KDataSplitMixin.py | 28 ++++++++++---- src/mrpro/utils/__init__.py | 1 - src/mrpro/utils/modify_acq_info.py | 35 ------------------ tests/data/test_kdata.py | 32 +++++++++++++--- tests/utils/test_modify_acq_info.py | 18 --------- 9 files changed, 102 insertions(+), 83 deletions(-) delete mode 100644 src/mrpro/utils/modify_acq_info.py delete mode 100644 tests/utils/test_modify_acq_info.py diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index 37c8f5605..50d54befa 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -1,18 +1,37 @@ """Acquisition information dataclass.""" +import dataclasses from collections.abc import Callable, Sequence from dataclasses import dataclass -from typing import Self +from typing import Self, TypeVar import ismrmrd import numpy as np import torch +from einops import rearrange 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 +T = TypeVar('T', torch.Tensor, Rotation, SpatialDimension) + + +def rearrange_acq_info_fields(field: T, pattern: str, additional_info: dict[str, int] | None = None) -> T: + """Change the shape of the fields in AcqInfo.""" + axes_lengths = {} if additional_info is None else additional_info + if isinstance(field, Rotation): + return Rotation.from_matrix(rearrange(field.as_matrix(), pattern, **axes_lengths)) + elif isinstance(field, SpatialDimension): + return SpatialDimension( + z=rearrange(field.z, pattern, **axes_lengths), + y=rearrange(field.y, pattern, **axes_lengths), + x=rearrange(field.x, pattern, **axes_lengths), + ) + else: + return rearrange(field, pattern, **axes_lengths) + @dataclass(slots=True) class AcqIdx(MoveDataMixin): @@ -265,3 +284,21 @@ def orientation_from_directional_cosine( version=tensor_2d(headers['version']), ) return acq_info + + def _apply_(self, modify_acq_info_field: Callable) -> None: + """Go through all fields of AcqInfo object and apply function in-place. + + Parameters + ---------- + modify_acq_info_field + Function which takes AcqInfo fields as input and returns modified AcqInfo field + """ + for field in dataclasses.fields(self): + current = getattr(self, field.name) + if dataclasses.is_dataclass(current): + for subfield in dataclasses.fields(current): + subcurrent = getattr(current, subfield.name) + setattr(current, subfield.name, modify_acq_info_field(subcurrent)) + else: + setattr(self, field.name, modify_acq_info_field(current)) + return None diff --git a/src/mrpro/data/_kdata/KData.py b/src/mrpro/data/_kdata/KData.py index 81b4f9003..d0025afbe 100644 --- a/src/mrpro/data/_kdata/KData.py +++ b/src/mrpro/data/_kdata/KData.py @@ -23,9 +23,9 @@ from mrpro.data.KTrajectory import KTrajectory from mrpro.data.KTrajectoryRawShape import KTrajectoryRawShape from mrpro.data.MoveDataMixin import MoveDataMixin +from mrpro.data.Rotation import Rotation from mrpro.data.traj_calculators.KTrajectoryCalculator import KTrajectoryCalculator from mrpro.data.traj_calculators.KTrajectoryIsmrmrd import KTrajectoryIsmrmrd -from mrpro.utils import modify_acq_info KDIM_SORT_LABELS = ( 'k1', @@ -199,10 +199,10 @@ def from_file( sort_idx = np.lexsort(acq_indices) # torch does not have lexsort as of pytorch 2.2 (March 2024) # Finally, reshape and sort the tensors in acqinfo and acqinfo.idx, and kdata. - def sort_and_reshape_tensor_fields(input_tensor: torch.Tensor): + def sort_and_reshape_tensor_fields(input_tensor: torch.Tensor | Rotation): return rearrange(input_tensor[sort_idx], '(other k2 k1) ... -> other k2 k1 ...', k1=n_k1, k2=n_k2) - kheader.acq_info = modify_acq_info(sort_and_reshape_tensor_fields, kheader.acq_info) + kheader.acq_info._apply_(sort_and_reshape_tensor_fields) kdata = rearrange(kdata[sort_idx], '(other k2 k1) coils k0 -> other coils k2 k1 k0', k1=n_k1, k2=n_k2) # Calculate trajectory and check if it matches the kdata shape diff --git a/src/mrpro/data/_kdata/KDataRearrangeMixin.py b/src/mrpro/data/_kdata/KDataRearrangeMixin.py index 0ae1d413c..3eba54b60 100644 --- a/src/mrpro/data/_kdata/KDataRearrangeMixin.py +++ b/src/mrpro/data/_kdata/KDataRearrangeMixin.py @@ -5,8 +5,7 @@ from einops import rearrange from mrpro.data._kdata.KDataProtocol import _KDataProtocol -from mrpro.data.AcqInfo import AcqInfo -from mrpro.utils import modify_acq_info +from mrpro.data.AcqInfo import rearrange_acq_info_fields class KDataRearrangeMixin(_KDataProtocol): @@ -34,9 +33,8 @@ def rearrange_k2_k1_into_k1(self: Self) -> Self: kheader = copy.deepcopy(self.header) # Update shape of acquisition info index - def reshape_acq_info(info: AcqInfo): - return rearrange(info, 'other k2 k1 ... -> other 1 (k2 k1) ...') - - kheader.acq_info = modify_acq_info(reshape_acq_info, kheader.acq_info) + kheader.acq_info._apply_( + lambda field: rearrange_acq_info_fields(field, 'other k2 k1 ... -> other 1 (k2 k1) ...') + ) return type(self)(kheader, kdat, type(self.traj).from_tensor(ktraj)) diff --git a/src/mrpro/data/_kdata/KDataSelectMixin.py b/src/mrpro/data/_kdata/KDataSelectMixin.py index 351f41938..463104c2e 100644 --- a/src/mrpro/data/_kdata/KDataSelectMixin.py +++ b/src/mrpro/data/_kdata/KDataSelectMixin.py @@ -1,11 +1,14 @@ """Select subset along other dimensions of KData.""" import copy -from typing import Literal, Self +from typing import Literal, Self, TypeVar import torch from mrpro.data._kdata.KDataProtocol import _KDataProtocol -from mrpro.utils import modify_acq_info +from mrpro.data.Rotation import Rotation +from mrpro.data.SpatialDimension import SpatialDimension + +T = TypeVar('T', torch.Tensor, Rotation, SpatialDimension) class KDataSelectMixin(_KDataProtocol): @@ -49,10 +52,13 @@ def select_other_subset( other_idx = torch.cat([torch.where(idx == label_idx[:, 0, 0])[0] for idx in subset_idx], dim=0) # Adapt header - def select_acq_info(info: torch.Tensor): - return info[other_idx, ...] + def select_acq_info(field: T) -> T: + if isinstance(field, SpatialDimension): + return SpatialDimension(z=field.z[other_idx, ...], y=field.y[other_idx, ...], x=field.x[other_idx, ...]) + else: + return field[other_idx, ...] - kheader.acq_info = modify_acq_info(select_acq_info, kheader.acq_info) + kheader.acq_info._apply_(select_acq_info) # Select data kdat = self.data[other_idx, ...] diff --git a/src/mrpro/data/_kdata/KDataSplitMixin.py b/src/mrpro/data/_kdata/KDataSplitMixin.py index aa6a05481..8276070a2 100644 --- a/src/mrpro/data/_kdata/KDataSplitMixin.py +++ b/src/mrpro/data/_kdata/KDataSplitMixin.py @@ -1,13 +1,18 @@ """Mixin class to split KData into other subsets.""" import copy -from typing import Literal, Self +from typing import Literal, Self, TypeVar import torch from einops import rearrange, repeat from mrpro.data._kdata.KDataProtocol import _KDataProtocol +from mrpro.data.AcqInfo import rearrange_acq_info_fields from mrpro.data.EncodingLimits import Limits -from mrpro.utils import modify_acq_info +from mrpro.data.Rotation import Rotation +from mrpro.data.SpatialDimension import SpatialDimension + +T = TypeVar('T', torch.Tensor, Rotation, SpatialDimension) +R = TypeVar('R', torch.Tensor, Rotation) class KDataSplitMixin(_KDataProtocol): @@ -54,7 +59,7 @@ def _split_k2_or_k1_into_other( def split_data_traj(dat_traj: torch.Tensor) -> torch.Tensor: return dat_traj[:, :, :, split_idx, :] - def split_acq_info(acq_info: torch.Tensor) -> torch.Tensor: + def split_acq_info_tensor(acq_info: R) -> R: return acq_info[:, :, split_idx, ...] # Rearrange other_split and k1 dimension @@ -67,7 +72,7 @@ def split_acq_info(acq_info: torch.Tensor) -> torch.Tensor: def split_data_traj(dat_traj: torch.Tensor) -> torch.Tensor: return dat_traj[:, :, split_idx, :, :] - def split_acq_info(acq_info: torch.Tensor) -> torch.Tensor: + def split_acq_info_tensor(acq_info: R) -> R: return acq_info[:, split_idx, ...] # Rearrange other_split and k1 dimension @@ -94,10 +99,17 @@ def split_acq_info(acq_info: torch.Tensor) -> torch.Tensor: kheader = copy.deepcopy(self.header) # Update shape of acquisition info index - def reshape_acq_info(info: torch.Tensor): - return rearrange(split_acq_info(info), rearrange_pattern_acq_info) - - kheader.acq_info = modify_acq_info(reshape_acq_info, kheader.acq_info) + def split_acq_info(field: T) -> T: + if isinstance(field, SpatialDimension): + return SpatialDimension( + z=split_acq_info_tensor(field.z), y=split_acq_info_tensor(field.y), x=split_acq_info_tensor(field.x) + ) + else: + return split_acq_info_tensor(field) + + kheader.acq_info._apply_( + lambda field: rearrange_acq_info_fields(split_acq_info(field), rearrange_pattern_acq_info) + ) # Update other label limits and acquisition info setattr(kheader.encoding_limits, other_label, Limits(min=0, max=n_other - 1, center=0)) diff --git a/src/mrpro/utils/__init__.py b/src/mrpro/utils/__init__.py index 8385b21cc..eee55244c 100644 --- a/src/mrpro/utils/__init__.py +++ b/src/mrpro/utils/__init__.py @@ -3,7 +3,6 @@ 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 -from mrpro.utils.modify_acq_info import modify_acq_info from mrpro.utils.split_idx import split_idx import mrpro.utils.slice_profiles import mrpro.utils.unit_conversion diff --git a/src/mrpro/utils/modify_acq_info.py b/src/mrpro/utils/modify_acq_info.py deleted file mode 100644 index d535e53c5..000000000 --- a/src/mrpro/utils/modify_acq_info.py +++ /dev/null @@ -1,35 +0,0 @@ -"""Modify AcqInfo.""" - -from __future__ import annotations - -import dataclasses -from collections.abc import Callable -from typing import TYPE_CHECKING - -import torch - -if TYPE_CHECKING: - from mrpro.data.AcqInfo import AcqInfo - - -def modify_acq_info(fun_modify: Callable, acq_info: AcqInfo) -> AcqInfo: - """Go through all fields of AcqInfo object and apply changes. - - Parameters - ---------- - fun_modify - Function which takes AcqInfo fields as input and returns modified AcqInfo field - acq_info - AcqInfo object - """ - # Apply function to all fields of acq_info - for field in dataclasses.fields(acq_info): - current = getattr(acq_info, field.name) - if isinstance(current, torch.Tensor): - setattr(acq_info, field.name, fun_modify(current)) - elif dataclasses.is_dataclass(current): - for subfield in dataclasses.fields(current): - subcurrent = getattr(current, subfield.name) - setattr(current, subfield.name, fun_modify(subcurrent)) - - return acq_info diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index 9c1bfd9d8..d1edb06d5 100644 --- a/tests/data/test_kdata.py +++ b/tests/data/test_kdata.py @@ -2,12 +2,13 @@ import pytest import torch -from einops import rearrange, repeat +from einops import repeat from mrpro.data import KData, KTrajectory, SpatialDimension from mrpro.data.acq_filters import is_coil_calibration_acquisition +from mrpro.data.AcqInfo import rearrange_acq_info_fields from mrpro.data.traj_calculators.KTrajectoryCalculator import DummyTrajectory from mrpro.operators import FastFourierOp -from mrpro.utils import modify_acq_info, split_idx +from mrpro.utils import split_idx from tests.conftest import RandomGenerator, generate_random_data from tests.data import IsmrmrdRawTestData @@ -77,10 +78,11 @@ def consistently_shaped_kdata(request, random_kheader_shape): # Start with header kheader, n_other, n_coils, n_k2, n_k1, n_k0 = random_kheader_shape - def reshape_acq_data(data): - return rearrange(data, '(other k2 k1) ... -> other k2 k1 ...', other=n_other, k2=n_k2, k1=n_k1) - - kheader.acq_info = modify_acq_info(reshape_acq_data, kheader.acq_info) + kheader.acq_info._apply_( + lambda field: rearrange_acq_info_fields( + field, '(other k2 k1) ... -> other k2 k1 ...', {'other': n_other, 'k2': n_k2, 'k1': n_k1} + ) + ) # Create kdata with consistent shape kdata = generate_random_data(RandomGenerator(request.param['seed']), (n_other, n_coils, n_k2, n_k1, n_k0)) @@ -469,3 +471,21 @@ def test_KData_remove_readout_os(monkeypatch, random_kheader): # testing functions such as numpy.testing.assert_almost_equal fails because there are few voxels with high # differences along the edges of the elliptic objects. assert relative_image_difference(torch.abs(img_recon), img_tensor[:, 0, ...]) <= 0.05 + + +def test_modify_acq_info(random_kheader_shape): + """Test the modification of the acquisition info.""" + # Create random header where AcqInfo fields are of shape [n_k1*n_k2] and reshape to [n_other, n_k2, n_k1] + kheader, n_other, _, n_k2, n_k1, _ = random_kheader_shape + + kheader.acq_info._apply_( + lambda field: rearrange_acq_info_fields( + field, '(other k2 k1) ... -> other k2 k1 ...', {'other': n_other, 'k2': n_k2, 'k1': n_k1} + ) + ) + + # Verify shape + assert kheader.acq_info.center_sample.shape == (n_other, n_k2, n_k1, 1) + 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) diff --git a/tests/utils/test_modify_acq_info.py b/tests/utils/test_modify_acq_info.py deleted file mode 100644 index 451303d02..000000000 --- a/tests/utils/test_modify_acq_info.py +++ /dev/null @@ -1,18 +0,0 @@ -"""Tests for modification of acquisition infos.""" - -from einops import rearrange -from mrpro.utils import modify_acq_info - - -def test_modify_acq_info(random_kheader_shape): - """Test the modification of the acquisition info.""" - # Create random header where AcqInfo fields are of shape [n_k1*n_k2] and reshape to [n_other, n_k2, n_k1] - kheader, n_other, _, n_k2, n_k1, _ = random_kheader_shape - - def reshape_acq_data(data): - return rearrange(data, '(other k2 k1) ... -> other k2 k1 ...', other=n_other, k2=n_k2, k1=n_k1) - - kheader.acq_info = modify_acq_info(reshape_acq_data, kheader.acq_info) - - # Verify shape - assert kheader.acq_info.center_sample.shape == (n_other, n_k2, n_k1, 1) From 539acff734c213a1bb5fa8c1c4efc07f311ebfbc Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Wed, 2 Oct 2024 16:58:15 +0200 Subject: [PATCH 09/19] adapt _apply_ in kdata --- src/mrpro/data/_kdata/KData.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/mrpro/data/_kdata/KData.py b/src/mrpro/data/_kdata/KData.py index d0025afbe..6119bd8ca 100644 --- a/src/mrpro/data/_kdata/KData.py +++ b/src/mrpro/data/_kdata/KData.py @@ -17,13 +17,12 @@ 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.AcqInfo import AcqInfo, rearrange_acq_info_fields from mrpro.data.EncodingLimits import Limits from mrpro.data.KHeader import KHeader from mrpro.data.KTrajectory import KTrajectory from mrpro.data.KTrajectoryRawShape import KTrajectoryRawShape from mrpro.data.MoveDataMixin import MoveDataMixin -from mrpro.data.Rotation import Rotation from mrpro.data.traj_calculators.KTrajectoryCalculator import KTrajectoryCalculator from mrpro.data.traj_calculators.KTrajectoryIsmrmrd import KTrajectoryIsmrmrd @@ -199,10 +198,11 @@ def from_file( sort_idx = np.lexsort(acq_indices) # torch does not have lexsort as of pytorch 2.2 (March 2024) # Finally, reshape and sort the tensors in acqinfo and acqinfo.idx, and kdata. - def sort_and_reshape_tensor_fields(input_tensor: torch.Tensor | Rotation): - return rearrange(input_tensor[sort_idx], '(other k2 k1) ... -> other k2 k1 ...', k1=n_k1, k2=n_k2) - - kheader.acq_info._apply_(sort_and_reshape_tensor_fields) + kheader.acq_info._apply_( + lambda field: rearrange_acq_info_fields( + field, '(other k2 k1) ... -> other k2 k1 ...', {'k1': n_k1, 'k2': n_k2} + ) + ) kdata = rearrange(kdata[sort_idx], '(other k2 k1) coils k0 -> other coils k2 k1 k0', k1=n_k1, k2=n_k2) # Calculate trajectory and check if it matches the kdata shape From 69dc4f17fea81ded17be31a586c534ab16d264e6 Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Wed, 2 Oct 2024 20:13:20 +0200 Subject: [PATCH 10/19] merge conflicts and xyz for dir --- src/mrpro/data/AcqInfo.py | 10 ++++------ src/mrpro/data/_kdata/KData.py | 5 +++-- tests/data/test_kdata.py | 3 --- tests/data/test_kheader.py | 2 +- 4 files changed, 8 insertions(+), 12 deletions(-) diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index ec89846bd..dbb9ab4eb 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -285,7 +285,6 @@ def orientation_from_directional_cosine( ) return acq_info -<<<<<<< HEAD def add_to_ismrmrd_acquisition( self, acquisition: ismrmrd.Acquisition, other: int, k2: int, k1: int ) -> ismrmrd.Acquisition: @@ -324,22 +323,22 @@ def add_to_ismrmrd_acquisition( m_to_mm(self.patient_table_position.y[other, k2, k1, 0]), m_to_mm(self.patient_table_position.z[other, k2, k1, 0]), ) - acquisition.phase_dir = tuple(self.orientation.as_matrix()[other, k2, k1, 1, :]) + acquisition.phase_dir = tuple(self.orientation.as_matrix()[other, k2, k1, 0, 1, :])[::-1] # zyx -> xyz acquisition.physiology_time_stamp = tuple(self.physiology_time_stamp[other, k2, k1, :]) acquisition.position = ( m_to_mm(self.position.x[other, k2, k1, 0]), m_to_mm(self.position.y[other, k2, k1, 0]), m_to_mm(self.position.z[other, k2, k1, 0]), ) - acquisition.read_dir = tuple(self.orientation.as_matrix()[other, k2, k1, 2, :]) + acquisition.read_dir = tuple(self.orientation.as_matrix()[other, k2, k1, 0, 2, :])[::-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 = tuple(self.orientation.as_matrix()[other, k2, k1, 0, :]) + acquisition.slice_dir = tuple(self.orientation.as_matrix()[other, k2, k1, 0, 0, :])[::-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 -======= + def _apply_(self, modify_acq_info_field: Callable) -> None: """Go through all fields of AcqInfo object and apply function in-place. @@ -357,4 +356,3 @@ def _apply_(self, modify_acq_info_field: Callable) -> None: else: setattr(self, field.name, modify_acq_info_field(current)) return None ->>>>>>> unit_conversion_utils diff --git a/src/mrpro/data/_kdata/KData.py b/src/mrpro/data/_kdata/KData.py index e21b65172..32085dba8 100644 --- a/src/mrpro/data/_kdata/KData.py +++ b/src/mrpro/data/_kdata/KData.py @@ -261,11 +261,12 @@ def to_file(self, filename: str | Path) -> None: # 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 range(self.data.shape[-5]): for k2 in range(self.data.shape[-3]): for k1 in range(self.data.shape[-2]): - acq = ismrmrd.Acquisition() - acq.resize(*acq_shape, trajectory_dimensions=3) + acq.clear_all_flags() acq = self.header.acq_info.add_to_ismrmrd_acquisition(acq, other, k2, k1) # Rearrange, switch from zyx to xz and set trajectory. diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index fa25b2a17..0d5469b8f 100644 --- a/tests/data/test_kdata.py +++ b/tests/data/test_kdata.py @@ -5,11 +5,8 @@ from einops import repeat from mrpro.data import KData, KTrajectory, SpatialDimension from mrpro.data.acq_filters import is_coil_calibration_acquisition -<<<<<<< HEAD from mrpro.data.traj_calculators import KTrajectoryIsmrmrd -======= from mrpro.data.AcqInfo import rearrange_acq_info_fields ->>>>>>> unit_conversion_utils from mrpro.data.traj_calculators.KTrajectoryCalculator import DummyTrajectory from mrpro.operators import FastFourierOp from mrpro.utils import split_idx diff --git a/tests/data/test_kheader.py b/tests/data/test_kheader.py index 4395eb86b..6a0092be1 100644 --- a/tests/data/test_kheader.py +++ b/tests/data/test_kheader.py @@ -48,7 +48,7 @@ def test_kheader_to_ismrmrd(random_mandatory_ismrmrd_header, random_acq_info): 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.acquisitionSystemInformation.systemFieldStrength_T == kheader.b0 + 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) From 8193e6a84999feda60cde1e8b670f12d1ced5c1e Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Wed, 2 Oct 2024 20:50:47 +0200 Subject: [PATCH 11/19] allow multiple other --- src/mrpro/data/AcqInfo.py | 84 +++++++++++++++++----------------- src/mrpro/data/_kdata/KData.py | 12 ++--- tests/data/test_kdata.py | 5 +- 3 files changed, 50 insertions(+), 51 deletions(-) diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index dbb9ab4eb..f8d5394ff 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -3,7 +3,7 @@ import dataclasses from collections.abc import Callable, Sequence from dataclasses import dataclass -from typing import Self, TypeVar +from typing import Any, Self, TypeVar import ismrmrd import numpy as np @@ -286,57 +286,57 @@ def orientation_from_directional_cosine( return acq_info def add_to_ismrmrd_acquisition( - self, acquisition: ismrmrd.Acquisition, other: int, k2: int, k1: int + self, acquisition: ismrmrd.Acquisition, other: Any, 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.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], + 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.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 = ( - m_to_mm(self.patient_table_position.x[other, k2, k1, 0]), - m_to_mm(self.patient_table_position.y[other, k2, k1, 0]), - m_to_mm(self.patient_table_position.z[other, k2, k1, 0]), + m_to_mm(self.patient_table_position.x[*other, k2, k1, 0]), + m_to_mm(self.patient_table_position.y[*other, k2, k1, 0]), + m_to_mm(self.patient_table_position.z[*other, k2, k1, 0]), ) - acquisition.phase_dir = tuple(self.orientation.as_matrix()[other, k2, k1, 0, 1, :])[::-1] # zyx -> xyz - acquisition.physiology_time_stamp = tuple(self.physiology_time_stamp[other, k2, k1, :]) + acquisition.phase_dir = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 1, :])[::-1] # zyx -> xyz + acquisition.physiology_time_stamp = tuple(self.physiology_time_stamp[*other, k2, k1, :]) acquisition.position = ( - m_to_mm(self.position.x[other, k2, k1, 0]), - m_to_mm(self.position.y[other, k2, k1, 0]), - m_to_mm(self.position.z[other, k2, k1, 0]), + m_to_mm(self.position.x[*other, k2, k1, 0]), + m_to_mm(self.position.y[*other, k2, k1, 0]), + m_to_mm(self.position.z[*other, k2, k1, 0]), ) - acquisition.read_dir = tuple(self.orientation.as_matrix()[other, k2, k1, 0, 2, :])[::-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 = tuple(self.orientation.as_matrix()[other, k2, k1, 0, 0, :])[::-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] + acquisition.read_dir = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 2, :])[::-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 = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 0, :])[::-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 def _apply_(self, modify_acq_info_field: Callable) -> None: diff --git a/src/mrpro/data/_kdata/KData.py b/src/mrpro/data/_kdata/KData.py index 32085dba8..8f1c78a8e 100644 --- a/src/mrpro/data/_kdata/KData.py +++ b/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 from typing import Self @@ -243,9 +244,6 @@ def to_file(self, filename: str | Path) -> None: filename path to the ISMRMRD file """ - if self.data.ndim != 5: - raise ValueError('Only KData with exactly 5 dimensions (other, coil, k2, k1, k0) can be saved.') - # Open the dataset dataset = ismrmrd.Dataset(filename, 'dataset', create_if_needed=True) @@ -263,17 +261,17 @@ def to_file(self, filename: str | Path) -> None: acq_shape = [self.data.shape[-1], self.data.shape[-4]] acq = ismrmrd.Acquisition() acq.resize(*acq_shape, trajectory_dimensions=3) - for other in range(self.data.shape[-5]): + 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 xz and set trajectory. - acq.traj[:] = rearrange(trajectory[:, other, k2, k1, :].numpy(), 'dim k0->k0 dim')[:, ::-1] + # 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() + acq.data[:] = self.data[*other, :, k2, k1, :].numpy() dataset.append_acquisition(acq) dataset.close() diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index 0d5469b8f..5c054f279 100644 --- a/tests/data/test_kdata.py +++ b/tests/data/test_kdata.py @@ -5,8 +5,8 @@ from einops import repeat from mrpro.data import KData, KTrajectory, SpatialDimension from mrpro.data.acq_filters import is_coil_calibration_acquisition -from mrpro.data.traj_calculators import KTrajectoryIsmrmrd 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 @@ -473,6 +473,7 @@ def test_KData_remove_readout_os(monkeypatch, random_kheader): # differences along the edges of the elliptic objects. assert relative_image_difference(torch.abs(img_recon), img_tensor[:, 0, ...]) <= 0.05 + def test_modify_acq_info(random_kheader_shape): """Test the modification of the acquisition info.""" # Create random header where AcqInfo fields are of shape [n_k1*n_k2] and reshape to [n_other, n_k2, n_k1] @@ -490,6 +491,7 @@ def test_modify_acq_info(random_kheader_shape): 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()) @@ -506,4 +508,3 @@ def test_KData_to_file(ismrmrd_cart, tmp_path_factory): torch.testing.assert_close(kdata.traj.as_tensor(), kdata_reload.traj.as_tensor()) 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 - From 885d4e00f6e0d3d6611b206b8ca05aac55d17f44 Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Thu, 3 Oct 2024 16:55:36 +0200 Subject: [PATCH 12/19] quaternion for siemens is phase read slice --- src/mrpro/data/AcqInfo.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index f8d5394ff..294e565a9 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -248,8 +248,8 @@ def orientation_from_directional_cosine( torch.stack( ( torch.stack((slice_dir.z, slice_dir.y, slice_dir.x), dim=-1), - torch.stack((phase_dir.z, phase_dir.y, phase_dir.x), dim=-1), torch.stack((read_dir.z, read_dir.y, read_dir.x), dim=-1), + torch.stack((phase_dir.z, phase_dir.y, phase_dir.x), dim=-1), ), dim=-2, ) @@ -323,14 +323,14 @@ def add_to_ismrmrd_acquisition( m_to_mm(self.patient_table_position.y[*other, k2, k1, 0]), m_to_mm(self.patient_table_position.z[*other, k2, k1, 0]), ) - acquisition.phase_dir = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 1, :])[::-1] # zyx -> xyz + acquisition.phase_dir = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 2, :])[::-1] # zyx -> xyz acquisition.physiology_time_stamp = tuple(self.physiology_time_stamp[*other, k2, k1, :]) acquisition.position = ( m_to_mm(self.position.x[*other, k2, k1, 0]), m_to_mm(self.position.y[*other, k2, k1, 0]), m_to_mm(self.position.z[*other, k2, k1, 0]), ) - acquisition.read_dir = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 2, :])[::-1] # zyx -> xyz + acquisition.read_dir = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 1, :])[::-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 = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 0, :])[::-1] # zyx -> xyz From d304aa104ab67dd953bc1a186b9d33ce1ef089a4 Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Tue, 12 Nov 2024 15:42:33 +0100 Subject: [PATCH 13/19] review and fix of orientation to read,phase,slice dir --- src/mrpro/data/AcqInfo.py | 33 +++++++++---------------------- src/mrpro/data/KHeader.py | 8 ++++---- tests/data/_IsmrmrdRawTestData.py | 6 +++--- tests/data/test_kdata.py | 6 ++++++ 4 files changed, 22 insertions(+), 31 deletions(-) diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index 94e6c7df2..c2166f579 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -2,7 +2,6 @@ from collections.abc import Sequence from dataclasses import dataclass -from typing import Any, Self import ismrmrd import numpy as np @@ -10,6 +9,7 @@ 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 @@ -230,23 +230,6 @@ def spatialdimension_2d(data: np.ndarray) -> SpatialDimension[torch.Tensor]: user7=tensor(idx['user'][:, 7]), ) - # Calculate orientation as rotation matrix from directional cosines - def orientation_from_directional_cosine( - slice_dir: SpatialDimension[torch.Tensor], - phase_dir: SpatialDimension[torch.Tensor], - read_dir: SpatialDimension[torch.Tensor], - ) -> Rotation: - return Rotation.from_matrix( - torch.stack( - ( - torch.stack((slice_dir.z, slice_dir.y, slice_dir.x), dim=-1), - torch.stack((read_dir.z, read_dir.y, read_dir.x), dim=-1), - torch.stack((phase_dir.z, phase_dir.y, phase_dir.x), dim=-1), - ), - dim=-2, - ) - ) - acq_info = cls( idx=acq_idx, acquisition_time_stamp=tensor_2d(headers['acquisition_time_stamp']), @@ -278,7 +261,7 @@ def orientation_from_directional_cosine( return acq_info def add_to_ismrmrd_acquisition( - self, acquisition: ismrmrd.Acquisition, other: Any, k2: int, k1: int + 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] @@ -310,14 +293,16 @@ def add_to_ismrmrd_acquisition( 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) - acquisition.phase_dir = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 2, :])[::-1] # zyx -> xyz + 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) - acquisition.read_dir = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 1, :])[::-1] # zyx -> xyz + 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 = tuple(self.orientation.as_matrix()[*other, k2, k1, 0, 0, :])[::-1] # zyx -> xyz + 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] diff --git a/src/mrpro/data/KHeader.py b/src/mrpro/data/KHeader.py index 24275193b..747a55dbf 100644 --- a/src/mrpro/data/KHeader.py +++ b/src/mrpro/data/KHeader.py @@ -161,13 +161,13 @@ def from_ismrmrd( parameters.update(defaults) if header.sequenceParameters is not None: - if any(header.sequenceParameters.TR): + if header.sequenceParameters.TR: parameters['tr'] = ms_to_s(torch.as_tensor(header.sequenceParameters.TR)) - if any(header.sequenceParameters.TE): + if header.sequenceParameters.TE: parameters['te'] = ms_to_s(torch.as_tensor(header.sequenceParameters.TE)) - if any(header.sequenceParameters.TI): + if header.sequenceParameters.TI: parameters['ti'] = ms_to_s(torch.as_tensor(header.sequenceParameters.TI)) - if any(header.sequenceParameters.flipAngle_deg): + if header.sequenceParameters.flipAngle_deg: parameters['fa'] = torch.deg2rad(torch.as_tensor(header.sequenceParameters.flipAngle_deg)) if header.sequenceParameters.echo_spacing: parameters['echo_spacing'] = ms_to_s(torch.as_tensor(header.sequenceParameters.echo_spacing)) diff --git a/tests/data/_IsmrmrdRawTestData.py b/tests/data/_IsmrmrdRawTestData.py index 181be56d6..b96d1a025 100644 --- a/tests/data/_IsmrmrdRawTestData.py +++ b/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 diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index 904ddda13..8981627bc 100644 --- a/tests/data/test_kdata.py +++ b/tests/data/test_kdata.py @@ -548,5 +548,11 @@ def test_KData_to_file(ismrmrd_cart, tmp_path_factory): 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 From ad79acb964952c348476601e67826f55106beae7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 14 Nov 2024 16:29:40 +0000 Subject: [PATCH 14/19] [pre-commit] auto fixes from pre-commit hooks --- tests/data/test_kdata.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/data/test_kdata.py b/tests/data/test_kdata.py index f8b876e29..88ebbd79e 100644 --- a/tests/data/test_kdata.py +++ b/tests/data/test_kdata.py @@ -576,6 +576,7 @@ def test_KData_compress_coils_error_coil_dim(consistently_shaped_kdata): with pytest.raises(ValueError, match='Coil dimension must not'): consistently_shaped_kdata.compress_coils(n_compressed_coils=3, joint_dims=(-4,)) + def test_KData_to_file(ismrmrd_cart, tmp_path_factory): """Read in data to file.""" kdata = KData.from_file(ismrmrd_cart.filename, DummyTrajectory()) @@ -597,4 +598,4 @@ def test_KData_to_file(ismrmrd_cart, tmp_path_factory): 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 \ No newline at end of file + assert kdata.header.encoding_limits.k1.max == kdata_reload.header.encoding_limits.k1.max From 9daee3f0ef194f752333c7e453e208a01bc3f5df Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Sat, 16 Nov 2024 20:45:27 +0100 Subject: [PATCH 15/19] change indexing --- src/mrpro/data/AcqInfo.py | 77 ++++++++++++++++------------------ src/mrpro/data/_kdata/KData.py | 2 +- 2 files changed, 38 insertions(+), 41 deletions(-) diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index c2166f579..e480a6b40 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -9,7 +9,6 @@ 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 @@ -260,50 +259,48 @@ def spatialdimension_2d(data: np.ndarray) -> SpatialDimension[torch.Tensor]: ) return acq_info - def add_to_ismrmrd_acquisition( - self, acquisition: ismrmrd.Acquisition, other: type_utils.TorchIndexerType, k2: int, k1: int - ) -> ismrmrd.Acquisition: + def add_to_ismrmrd_acquisition(self, acquisition: ismrmrd.Acquisition, idx: Sequence[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.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[*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], + 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[*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.acquisition_time_stamp = self.acquisition_time_stamp[*idx][0] + acquisition.available_channels = self.available_channels[*idx][0] + acquisition.center_sample = self.center_sample[*idx][0] + acquisition.channel_mask = tuple(self.channel_mask[*idx]) + acquisition.discard_post = self.discard_post[*idx][0] + acquisition.discard_pre = self.discard_pre[*idx][0] + acquisition.encoding_space_ref = self.encoding_space_ref[*idx][0] + acquisition.measurement_uid = self.measurement_uid[*idx][0] acquisition.patient_table_position = ( - self.patient_table_position[*other, k2, k1, 0].apply_(m_to_mm).zyx[::-1] + self.patient_table_position[*idx][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] + acquisition.phase_dir = self.orientation.as_directions()[-2][*idx][0].zyx[::-1] # zyx -> xyz + acquisition.physiology_time_stamp = tuple(self.physiology_time_stamp[*idx]) + acquisition.position = self.position[*idx][0].apply_(m_to_mm).zyx[::-1] # zyx -> xyz + acquisition.read_dir = self.orientation.as_directions()[-1][*idx][0].zyx[::-1] # zyx -> xyz + acquisition.sample_time_us = self.sample_time_us[*idx][0] + acquisition.scan_counter = self.scan_counter[*idx][0] + acquisition.slice_dir = self.orientation.as_directions()[-3][*idx][0].zyx[::-1] # zyx -> xyz + acquisition.user_float = tuple(self.user_float[*idx]) + acquisition.user_int = tuple(self.user_int[*idx]) + acquisition.version = self.version[*idx][0] return acquisition diff --git a/src/mrpro/data/_kdata/KData.py b/src/mrpro/data/_kdata/KData.py index ff4c19696..7ccf6061e 100644 --- a/src/mrpro/data/_kdata/KData.py +++ b/src/mrpro/data/_kdata/KData.py @@ -293,7 +293,7 @@ def to_file(self, filename: str | Path) -> None: 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) + 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] From 1a1de6e58dcf5c6465c2458698d5d0cb73ccdacf Mon Sep 17 00:00:00 2001 From: Felix Zimmermann Date: Wed, 20 Nov 2024 12:40:30 +0100 Subject: [PATCH 16/19] add apply (out of place) --- src/mrpro/data/MoveDataMixin.py | 18 +++++++++++++++++ src/mrpro/data/SpatialDimension.py | 13 ++++++++++++ tests/data/test_movedatamixin.py | 23 ++++++++++++++++++++- tests/data/test_spatial_dimension.py | 30 +++++++++++++++++++++++++++- 4 files changed, 82 insertions(+), 2 deletions(-) diff --git a/src/mrpro/data/MoveDataMixin.py b/src/mrpro/data/MoveDataMixin.py index 8d977d0a6..2ac3c1a58 100644 --- a/src/mrpro/data/MoveDataMixin.py +++ b/src/mrpro/data/MoveDataMixin.py @@ -239,6 +239,24 @@ def _convert(data: T) -> T: new.apply_(_convert, memo=memo, recurse=False) return new + def apply( + self: Self, + function: Callable[[Any], Any] | None = None, + *, + recurse: bool = True, + ) -> Self: + """Apply a function to all children. Returns a new object. + + Parameters + ---------- + function + The function to apply to all fields. None is interpreted as a no-op. + recurse + If True, the function will be applied to all children that are MoveDataMixin instances. + """ + new = self.clone().apply_(function, recurse=recurse) + return new + def apply_( self: Self, function: Callable[[Any], Any] | None = None, diff --git a/src/mrpro/data/SpatialDimension.py b/src/mrpro/data/SpatialDimension.py index 46b1db89a..221c39553 100644 --- a/src/mrpro/data/SpatialDimension.py +++ b/src/mrpro/data/SpatialDimension.py @@ -18,6 +18,7 @@ VectorTypes = torch.Tensor ScalarTypes = int | float T = TypeVar('T', torch.Tensor, int, float) + # Covariant types, as SpatialDimension is a Container # and we want, for example, SpatialDimension[int] to also be a SpatialDimension[float] T_co = TypeVar('T_co', torch.Tensor, int, float, covariant=True) @@ -116,8 +117,20 @@ 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) + def apply(self, function: Callable[[T], T] | None = None, **_) -> Self: + """Apply a function to each z, y, x (returning a new object). + + Parameters + ---------- + function + function to apply + """ + # This function is mainly used for type hinting + return super(SpatialDimension, self).apply(function) + @property def zyx(self) -> tuple[T_co, T_co, T_co]: """Return a z,y,x tuple.""" diff --git a/tests/data/test_movedatamixin.py b/tests/data/test_movedatamixin.py index 3feb091de..c92e12b2a 100644 --- a/tests/data/test_movedatamixin.py +++ b/tests/data/test_movedatamixin.py @@ -207,7 +207,7 @@ def testchild(attribute, expected_dtype): assert new.module.module1.weight is new.module.module1.weight, 'shared module parameters should remain shared' -def test_movedatamixin_apply(): +def test_movedatamixin_apply_(): """Tests apply_ method of MoveDataMixin.""" data = B() # make one of the parameters shared to test memo behavior @@ -223,3 +223,24 @@ def multiply_by_2(obj): torch.testing.assert_close(data.floattensor, original.floattensor * 2) torch.testing.assert_close(data.child.floattensor2, original.child.floattensor2 * 2) assert data.child.floattensor is data.child.floattensor2, 'shared module parameters should remain shared' + + +def test_movedatamixin_apply(): + """Tests apply method of MoveDataMixin.""" + data = B() + # make one of the parameters shared to test memo behavior + data.child.floattensor2 = data.child.floattensor + original = data.clone() + + def multiply_by_2(obj): + if isinstance(obj, torch.Tensor): + return obj * 2 + return obj + + new = data.apply(multiply_by_2) + torch.testing.assert_close(data.floattensor, original.floattensor) + torch.testing.assert_close(data.child.floattensor2, original.child.floattensor2) + torch.testing.assert_close(new.floattensor, original.floattensor * 2) + torch.testing.assert_close(new.child.floattensor2, original.child.floattensor2 * 2) + assert data.child.floattensor is data.child.floattensor2, 'shared module parameters should remain shared' + assert new is not data, 'new object should be different from the original' diff --git a/tests/data/test_spatial_dimension.py b/tests/data/test_spatial_dimension.py index afafece04..61fd127df 100644 --- a/tests/data/test_spatial_dimension.py +++ b/tests/data/test_spatial_dimension.py @@ -93,7 +93,7 @@ def test_spatial_dimension_broadcasting(): def test_spatial_dimension_apply_(): - """Test apply_ (inplace)""" + """Test apply_ (in place)""" def conversion(x: torch.Tensor) -> torch.Tensor: assert isinstance(x, torch.Tensor), 'The argument to the conversion function should be a tensor' @@ -115,6 +115,34 @@ def conversion(x: torch.Tensor) -> torch.Tensor: assert torch.equal(spatial_dimension_inplace.z, z) +def test_spatial_dimension_apply(): + """Test apply (out of place)""" + + def conversion(x: torch.Tensor) -> torch.Tensor: + assert isinstance(x, torch.Tensor), 'The argument to the conversion function should be a tensor' + return x.swapaxes(0, 1).square() + + xyz = RandomGenerator(0).float32_tensor((1, 2, 3)) + spatial_dimension = SpatialDimension.from_array_xyz(xyz.numpy()) + spatial_dimension_outofplace = spatial_dimension.apply(conversion) + + assert spatial_dimension_outofplace is not spatial_dimension + + assert isinstance(spatial_dimension_outofplace.x, torch.Tensor) + assert isinstance(spatial_dimension_outofplace.y, torch.Tensor) + assert isinstance(spatial_dimension_outofplace.z, torch.Tensor) + + x, y, z = conversion(xyz).unbind(-1) + assert torch.equal(spatial_dimension_outofplace.x, x) + assert torch.equal(spatial_dimension_outofplace.y, y) + assert torch.equal(spatial_dimension_outofplace.z, z) + + x, y, z = xyz.unbind(-1) # original should be unmodified + assert torch.equal(spatial_dimension.x, x) + assert torch.equal(spatial_dimension.y, y) + assert torch.equal(spatial_dimension.z, z) + + def test_spatial_dimension_zyx(): """Test the zyx tuple property""" z, y, x = (2, 3, 4) From fbf8c635bfa8135285860e3f1b47c00f513a40bf Mon Sep 17 00:00:00 2001 From: Felix Zimmermann Date: Wed, 20 Nov 2024 13:35:44 +0100 Subject: [PATCH 17/19] fix --- src/mrpro/data/AcqInfo.py | 83 ++++++++++++++++++---------------- src/mrpro/data/_kdata/KData.py | 18 ++++---- 2 files changed, 52 insertions(+), 49 deletions(-) diff --git a/src/mrpro/data/AcqInfo.py b/src/mrpro/data/AcqInfo.py index e480a6b40..6c478b9ce 100644 --- a/src/mrpro/data/AcqInfo.py +++ b/src/mrpro/data/AcqInfo.py @@ -12,6 +12,7 @@ from mrpro.data.MoveDataMixin import MoveDataMixin from mrpro.data.Rotation import Rotation from mrpro.data.SpatialDimension import SpatialDimension +from mrpro.utils.typing import TorchIndexerType from mrpro.utils.unit_conversion import m_to_mm, mm_to_m @@ -259,48 +260,50 @@ def spatialdimension_2d(data: np.ndarray) -> SpatialDimension[torch.Tensor]: ) return acq_info - def add_to_ismrmrd_acquisition(self, acquisition: ismrmrd.Acquisition, idx: Sequence[int]) -> ismrmrd.Acquisition: - """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] + 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], + 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][0] - acquisition.available_channels = self.available_channels[*idx][0] - acquisition.center_sample = self.center_sample[*idx][0] - acquisition.channel_mask = tuple(self.channel_mask[*idx]) - acquisition.discard_post = self.discard_post[*idx][0] - acquisition.discard_pre = self.discard_pre[*idx][0] - acquisition.encoding_space_ref = self.encoding_space_ref[*idx][0] - acquisition.measurement_uid = self.measurement_uid[*idx][0] - acquisition.patient_table_position = ( - self.patient_table_position[*idx][0].apply_(m_to_mm).zyx[::-1] - ) # zyx -> xyz - acquisition.phase_dir = self.orientation.as_directions()[-2][*idx][0].zyx[::-1] # zyx -> xyz - acquisition.physiology_time_stamp = tuple(self.physiology_time_stamp[*idx]) - acquisition.position = self.position[*idx][0].apply_(m_to_mm).zyx[::-1] # zyx -> xyz - acquisition.read_dir = self.orientation.as_directions()[-1][*idx][0].zyx[::-1] # zyx -> xyz - acquisition.sample_time_us = self.sample_time_us[*idx][0] - acquisition.scan_counter = self.scan_counter[*idx][0] - acquisition.slice_dir = self.orientation.as_directions()[-3][*idx][0].zyx[::-1] # zyx -> xyz - acquisition.user_float = tuple(self.user_float[*idx]) - acquisition.user_int = tuple(self.user_int[*idx]) - acquisition.version = self.version[*idx][0] + 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 diff --git a/src/mrpro/data/_kdata/KData.py b/src/mrpro/data/_kdata/KData.py index 7ccf6061e..7d7edf281 100644 --- a/src/mrpro/data/_kdata/KData.py +++ b/src/mrpro/data/_kdata/KData.py @@ -4,7 +4,6 @@ import datetime import warnings from collections.abc import Callable, Sequence -from itertools import product from pathlib import Path from types import EllipsisType @@ -280,26 +279,27 @@ def to_file(self, filename: str | Path) -> None: 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)] - ) + 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 product(*[range(shape) for shape in self.data.shape[:-4]]): + + 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() - acq = self.header.acq_info.add_to_ismrmrd_acquisition(acq, (*other, k2, k1)) + self.header.acq_info.write_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] + acq.traj[:] = trajectory[(*other, k2, k1)].numpy()[:, ::-1] # zyx -> xyz # Set the data and append - acq.data[:] = self.data[*other, :, k2, k1, :].numpy() + 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() From 025822ad7b5d81266aef9e55b4d2546ddcf5194d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 22 Nov 2024 13:25:34 +0000 Subject: [PATCH 18/19] [pre-commit] auto fixes from pre-commit hooks --- src/mrpro/data/SpatialDimension.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/mrpro/data/SpatialDimension.py b/src/mrpro/data/SpatialDimension.py index 28b0d4a1c..0b7735823 100644 --- a/src/mrpro/data/SpatialDimension.py +++ b/src/mrpro/data/SpatialDimension.py @@ -121,7 +121,6 @@ def apply_(self, function: Callable[[T], T] | None = None, **_) -> Self: # This function is mainly used for type hinting return super(SpatialDimension, self).apply_(function) - # This function is mainly for type hinting and docstring def apply(self, function: Callable[[T], T] | None = None, **_) -> Self: """Apply a function to each z, y, x (returning a new object). From 62c3289a2f6a7b570c758dd22cc08c04cb8b5ad3 Mon Sep 17 00:00:00 2001 From: ckolbPTB Date: Tue, 26 Nov 2024 16:04:29 +0100 Subject: [PATCH 19/19] review --- src/mrpro/data/KHeader.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mrpro/data/KHeader.py b/src/mrpro/data/KHeader.py index 747a55dbf..05c26e8b6 100644 --- a/src/mrpro/data/KHeader.py +++ b/src/mrpro/data/KHeader.py @@ -295,15 +295,15 @@ def to_ismrmrd(self) -> ismrmrdschema.ismrmrdHeader: # Sequence information seq = ismrmrdschema.sequenceParametersType() - if isinstance(self.tr, torch.Tensor): + if self.tr is not None: seq.TR = s_to_ms(self.tr).tolist() - if isinstance(self.te, torch.Tensor): + if self.te is not None: seq.TE = s_to_ms(self.te).tolist() - if isinstance(self.ti, torch.Tensor): + if self.ti is not None: seq.TI = s_to_ms(self.ti).tolist() - if isinstance(self.fa, torch.Tensor): + if self.fa is not None: seq.flipAngle_deg = torch.rad2deg(self.fa).tolist() - if isinstance(self.echo_spacing, torch.Tensor): + 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