Cartesian trajectory calculator added

src/mrpro/data/traj_calculators/_KTrajectoryCalculator.py

@@ -36,6 +36,39 @@ def __call__(self, header: KHeader) -> KTrajectory | KTrajectoryRawShape:
         """
         ...
 
+    def _kfreq(self, kheader: KHeader) -> torch.Tensor:
+        """Calculate the trajectory along one readout (k0 dimension).
+
+        Parameters
+        ----------
+        kheader
+            MR raw data header (KHeader) containing required meta data
+
+        Returns
+        -------
+            trajectory along ONE readout
+
+        Raises
+        ------
+        ValueError
+            Number of samples have to be the same for each readout
+        ValueError
+            Center sample has to be the same for each readout
+        """
+        num_samples = kheader.acq_info.number_of_samples
+        center_sample = kheader.acq_info.center_sample
+
+        # Verify that each readout has the same number of samples and same center sample
+        if len(torch.unique(num_samples)) > 1:
+            raise ValueError('Trajectory can only be calculated if each acquisition has the same number of samples')
+        if len(torch.unique(center_sample)) > 1:
+            raise ValueError('Trajectory can only be calculated if each acquisition has the same center sample')
+
+        # Calculate points along readout
+        nk0 = int(num_samples[0, 0, 0])
+        k0 = torch.linspace(0, nk0 - 1, nk0, dtype=torch.float32) - center_sample[0, 0, 0]
+        return k0
+
 
 class DummyTrajectory(KTrajectoryCalculator):
     """Simple Dummy trajectory that returns zeros.

src/mrpro/data/traj_calculators/_KTrajectoryCartesian.py

@@ -0,0 +1,50 @@
+"""Cartesian trajectory class."""
+
+# Copyright 2023 Physikalisch-Technische Bundesanstalt
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#       http://www.apache.org/licenses/LICENSE-2.0
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+import torch
+from einops import repeat
+
+from mrpro.data import KHeader
+from mrpro.data import KTrajectory
+from mrpro.data.traj_calculators import KTrajectoryCalculator
+
+
+class KTrajectoryCartesian(KTrajectoryCalculator):
+    """Cartesian trajectory."""
+
+    def __call__(self, kheader: KHeader) -> KTrajectory:
+        """Calculate Cartesian trajectory for given KHeader.
+
+        Parameters
+        ----------
+        kheader
+           MR raw data header (KHeader) containing required meta data
+
+        Returns
+        -------
+            Cartesian trajectory for given KHeader
+        """
+
+        # K-space locations along readout lines
+        kx = self._kfreq(kheader)
+
+        # Trajectory along phase and slice encoding
+        ky = (kheader.acq_info.idx.k1 - kheader.encoding_limits.k1.center).to(torch.float32)
+        kz = (kheader.acq_info.idx.k2 - kheader.encoding_limits.k2.center).to(torch.float32)
+
+        # Bring to correct dimensions
+        kx = repeat(kx, 'k0-> other k2 k1 k0', other=1, k2=1, k1=1)
+        ky = repeat(ky, 'other k2 k1-> other k2 k1 k0', k0=1)
+        kz = repeat(kz, 'other k2 k1-> other k2 k1 k0', k0=1)
+        return KTrajectory(kz, ky, kx)

src/mrpro/data/traj_calculators/_KTrajectoryIsmrmrd.py

@@ -1,4 +1,4 @@
-"""Returns the trajectory save in an ISMRMRD raw data file."""
+"""Returns the trajectory saved in an ISMRMRD raw data file."""
 
 # Copyright 2023 Physikalisch-Technische Bundesanstalt
 #

src/mrpro/data/traj_calculators/_KTrajectoryRadial2D.py

@@ -1,4 +1,17 @@
-# Imports
+"""2D radial trajectory class."""
+
+# Copyright 2023 Physikalisch-Technische Bundesanstalt
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#       http://www.apache.org/licenses/LICENSE-2.0
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
 import torch
 
 from mrpro.data import KHeader
@@ -22,43 +35,6 @@ def __init__(
         super().__init__()
         self.angle: float = angle
 
-    def _krad(self, kheader: KHeader) -> torch.Tensor:
-        """Calculate the k-space locations along the read-out encoding lines.
-
-        Parameters
-        ----------
-        kheader
-            MR raw data header (KHeader) containing required meta data
-
-        Returns
-        -------
-            trajectory along ONE readout
-
-        Raises
-        ------
-        ValueError
-            Number of samples have to be the same for each readout
-        ValueError
-            Center sample has to be the same for each readout
-        """
-
-        num_samples = kheader.acq_info.number_of_samples
-        center_sample = kheader.acq_info.center_sample
-
-        if len(torch.unique(num_samples)) > 1:
-            raise ValueError(
-                'Radial 2D trajectory can only be calculated if each acquisition has the same number of samples'
-            )
-        if len(torch.unique(center_sample)) > 1:
-            raise ValueError(
-                'Radial 2D trajectory can only be calculated if each acquisition has the same center sample'
-            )
-
-        # Calculate points along readout
-        nk0 = int(num_samples[0, 0, 0])
-        k0 = torch.linspace(0, nk0 - 1, nk0, dtype=torch.float32) - center_sample[0, 0, 0]
-        return k0
-
     def __call__(self, kheader: KHeader) -> KTrajectory:
         """Calculate radial 2D trajectory for given KHeader.
 
@@ -73,7 +49,7 @@ def __call__(self, kheader: KHeader) -> KTrajectory:
         """
 
         # K-space locations along readout lines
-        krad = self._krad(kheader)
+        krad = self._kfreq(kheader)
 
         # Angles of readout lines
         kang = kheader.acq_info.idx.k1 * self.angle

src/mrpro/data/traj_calculators/_KTrajectoryRpe.py

@@ -82,39 +82,6 @@ def _apply_shifts_between_rpe_lines(self, krad: torch.Tensor, kang_idx: torch.Te
             krad[curr_angle_idx] = curr_krad
         return krad
 
-    def _kfreq(self, kheader: KHeader) -> torch.Tensor:
-        """Calculate the trajectory along one readout (k0 dimension).
-
-        Parameters
-        ----------
-        kheader
-            MR raw data header (KHeader) containing required meta data
-
-        Returns
-        -------
-            trajectory along ONE readout
-
-        Raises
-        ------
-        ValueError
-            Number of samples have to be the same for each readout
-        ValueError
-            Center sample has to be the same for each readout
-        """
-        num_samples = kheader.acq_info.number_of_samples
-        center_sample = kheader.acq_info.center_sample
-
-        # Verify that each readout has the same number of samples and same center sample
-        if len(torch.unique(num_samples)) > 1:
-            raise ValueError('RPE trajectory can only be calculated if each acquisition has the same number of samples')
-        if len(torch.unique(center_sample)) > 1:
-            raise ValueError('RPE trajectory can only be calculated if each acquisition has the same center sample')
-
-        # Calculate points along readout
-        nk0 = int(num_samples[0, 0, 0])
-        k0 = torch.linspace(0, nk0 - 1, nk0, dtype=torch.float32) - center_sample[0, 0, 0]
-        return k0
-
     def _kang(self, kheader: KHeader) -> torch.Tensor:
         """Calculate the angles of the phase encoding lines.
 

src/mrpro/data/traj_calculators/__init__.py

@@ -4,3 +4,4 @@
 from mrpro.data.traj_calculators._KTrajectoryRadial2D import KTrajectoryRadial2D
 from mrpro.data.traj_calculators._KTrajectoryIsmrmrd import KTrajectoryIsmrmrd
 from mrpro.data.traj_calculators._KTrajectoryPulseq import KTrajectoryPulseq
+from mrpro.data.traj_calculators._KTrajectoryCartesian import KTrajectoryCartesian

tests/data/test_traj_calculators.py

@@ -17,23 +17,26 @@
 import torch
 
 from mrpro.data import KData
+from mrpro.data.traj_calculators import KTrajectoryCartesian
 from mrpro.data.traj_calculators import KTrajectoryIsmrmrd
+from mrpro.data.traj_calculators import KTrajectoryPulseq
 from mrpro.data.traj_calculators import KTrajectoryRadial2D
 from mrpro.data.traj_calculators import KTrajectoryRpe
 from mrpro.data.traj_calculators import KTrajectorySunflowerGoldenRpe
-from mrpro.data.traj_calculators._KTrajectoryPulseq import KTrajectoryPulseq
 from tests.conftest import random_kheader
 from tests.data import IsmrmrdRawTestData
 from tests.data._PulseqRadialTestSeq import PulseqRadialTestSeq
+from tests.data.test_kdata import ismrmrd_cart
 from tests.phantoms.test_ellipse_phantom import ph_ellipse
 
 
 @pytest.fixture(scope='function')
 def valid_rad2d_kheader(monkeypatch, random_kheader):
     """KHeader with all necessary parameters for radial 2D trajectories."""
     # K-space dimensions
-    nk0 = 200
-    nk1 = 20
+    nk0 = 256
+    nk1 = 10
+    nk2 = 1
 
     # List of k1 indices in the shape
     idx_k1 = torch.arange(nk1, dtype=torch.int32)[None, None, ...]
@@ -43,6 +46,11 @@ def valid_rad2d_kheader(monkeypatch, random_kheader):
     monkeypatch.setattr(random_kheader.acq_info, 'center_sample', torch.zeros_like(idx_k1) + nk0 // 2)
     monkeypatch.setattr(random_kheader.acq_info.idx, 'k1', idx_k1)
 
+    # This is only needed for Pulseq trajectory calculation
+    monkeypatch.setattr(random_kheader.encoding_matrix, 'x', nk0)
+    monkeypatch.setattr(random_kheader.encoding_matrix, 'y', nk0)  # square encoding in kx-ky plane
+    monkeypatch.setattr(random_kheader.encoding_matrix, 'z', nk2)
+
     return random_kheader
 
 
@@ -65,9 +73,6 @@ def test_KTrajectoryRadial2D_golden(valid_rad2d_kheader):
     assert ktraj.kz.shape == valid_shape[0]
 
 
-from tests.data.test_kdata import ismrmrd_cart
-
-
 @pytest.fixture(scope='function')
 def valid_rpe_kheader(monkeypatch, random_kheader):
     """KHeader with all necessary parameters for RPE trajectories."""
@@ -93,36 +98,6 @@ def valid_rpe_kheader(monkeypatch, random_kheader):
     return random_kheader
 
 
-@pytest.fixture(scope='function')
-def valid_radial_kheader(monkeypatch, random_kheader):
-    """KHeader with all necessary parameters for radial trajectories."""
-    # K-space dimensions
-    nk0 = 256
-    nk1 = 10
-    nk2 = 1
-
-    # List of k1 and k2 indices in the shape (other, k2, k1)
-    k1 = torch.linspace(0, nk1 - 1, nk1, dtype=torch.int32)
-    k2 = torch.linspace(0, nk2 - 1, nk2, dtype=torch.int32)
-    idx_k1, idx_k2 = torch.meshgrid(k1, k2, indexing='xy')
-    idx_k1 = torch.reshape(idx_k1, (1, nk2, nk1))
-    idx_k2 = torch.reshape(idx_k2, (1, nk2, nk1))
-
-    # Set parameters for radial trajectory
-    monkeypatch.setattr(random_kheader.acq_info, 'number_of_samples', torch.zeros_like(idx_k1) + nk0)
-    monkeypatch.setattr(random_kheader.acq_info, 'center_sample', torch.zeros_like(idx_k1) + nk0 // 2)
-    monkeypatch.setattr(random_kheader.acq_info.idx, 'k1', idx_k1)
-    monkeypatch.setattr(random_kheader.acq_info.idx, 'k2', idx_k2)
-    monkeypatch.setattr(random_kheader.encoding_limits.k1, 'center', int(nk1 // 2))
-    monkeypatch.setattr(random_kheader.encoding_limits.k1, 'max', int(nk1 - 1))
-
-    # This is only needed for Pulseq trajectory calculation
-    monkeypatch.setattr(random_kheader.encoding_matrix, 'x', nk0)
-    monkeypatch.setattr(random_kheader.encoding_matrix, 'y', nk0)  # square encoding in kx-ky plane
-    monkeypatch.setattr(random_kheader.encoding_matrix, 'z', nk2)
-    return random_kheader
-
-
 def rpe_traj_shape(valid_rpe_kheader):
     """Expected shape of trajectory based on KHeader."""
     nk0 = valid_rpe_kheader.acq_info.number_of_samples[0, 0, 0]
@@ -174,6 +149,52 @@ def test_KTrajectorySunflowerGoldenRpe(valid_rpe_kheader):
     assert ktraj.broadcasted_shape == np.broadcast_shapes(*rpe_traj_shape(valid_rpe_kheader))
 
 
+@pytest.fixture(scope='function')
+def valid_cartesian_kheader(monkeypatch, random_kheader):
+    """KHeader with all necessary parameters for Cartesian trajectories."""
+    # K-space dimensions
+    nk0 = 200
+    nk1 = 20
+    nk2 = 10
+
+    # List of k1 and k2 indices in the shape (other, k2, k1)
+    k1 = torch.linspace(0, nk1 - 1, nk1, dtype=torch.int32)
+    k2 = torch.linspace(0, nk2 - 1, nk2, dtype=torch.int32)
+    idx_k1, idx_k2 = torch.meshgrid(k1, k2, indexing='xy')
+    idx_k1 = torch.reshape(idx_k1, (1, nk2, nk1))
+    idx_k2 = torch.reshape(idx_k2, (1, nk2, nk1))
+
+    # Set parameters for Cartesian trajectory
+    monkeypatch.setattr(random_kheader.acq_info, 'number_of_samples', torch.zeros_like(idx_k1) + nk0)
+    monkeypatch.setattr(random_kheader.acq_info, 'center_sample', torch.zeros_like(idx_k1) + nk0 // 2)
+    monkeypatch.setattr(random_kheader.acq_info.idx, 'k1', idx_k1)
+    monkeypatch.setattr(random_kheader.acq_info.idx, 'k2', idx_k2)
+    monkeypatch.setattr(random_kheader.encoding_limits.k1, 'center', int(nk1 // 2))
+    monkeypatch.setattr(random_kheader.encoding_limits.k1, 'max', int(nk1 - 1))
+    monkeypatch.setattr(random_kheader.encoding_limits.k2, 'center', int(nk2 // 2))
+    monkeypatch.setattr(random_kheader.encoding_limits.k2, 'max', int(nk2 - 1))
+    return random_kheader
+
+
+def cartesian_traj_shape(valid_cartesian_kheader):
+    """Expected shape of trajectory based on KHeader."""
+    nk0 = valid_cartesian_kheader.acq_info.number_of_samples[0, 0, 0]
+    nk1 = valid_cartesian_kheader.acq_info.idx.k1.shape[2]
+    nk2 = valid_cartesian_kheader.acq_info.idx.k1.shape[1]
+    nother = 1
+    return (torch.Size([nother, nk2, 1, 1]), torch.Size([nother, 1, nk1, 1]), torch.Size([nother, 1, 1, nk0]))
+
+
+def test_KTrajectoryCartesian(valid_cartesian_kheader):
+    """Calculate Cartesian trajectory."""
+    ktrajectory = KTrajectoryCartesian()
+    ktraj = ktrajectory(valid_cartesian_kheader)
+    valid_shape = cartesian_traj_shape(valid_cartesian_kheader)
+    assert ktraj.kz.shape == valid_shape[0]
+    assert ktraj.ky.shape == valid_shape[1]
+    assert ktraj.kx.shape == valid_shape[2]
+
+
 @pytest.fixture(scope='session')
 def ismrmrd_rad(ph_ellipse, tmp_path_factory):
     """Data set with uniform radial k-space sampling."""
@@ -210,19 +231,19 @@ def pulseq_example_rad_seq(tmp_path_factory):
     return seq
 
 
-def test_KTrajectoryPulseq_validseq_random_header(pulseq_example_rad_seq, valid_radial_kheader):
+def test_KTrajectoryPulseq_validseq_random_header(pulseq_example_rad_seq, valid_rad2d_kheader):
     """Test pulseq File reader with valid seq File."""
     # TODO: Test with valid header
     # TODO: Test with invalid seq file
 
     ktrajectory = KTrajectoryPulseq(seq_path=pulseq_example_rad_seq.seq_filename)
-    traj = ktrajectory(kheader=valid_radial_kheader)
+    traj = ktrajectory(kheader=valid_rad2d_kheader)
 
     kx_test = pulseq_example_rad_seq.traj_analytical.kx.squeeze(0).squeeze(0)
-    kx_test *= valid_radial_kheader.encoding_matrix.x / (2 * torch.max(torch.abs(kx_test)))
+    kx_test *= valid_rad2d_kheader.encoding_matrix.x / (2 * torch.max(torch.abs(kx_test)))
 
     ky_test = pulseq_example_rad_seq.traj_analytical.ky.squeeze(0).squeeze(0)
-    ky_test *= valid_radial_kheader.encoding_matrix.y / (2 * torch.max(torch.abs(ky_test)))
+    ky_test *= valid_rad2d_kheader.encoding_matrix.y / (2 * torch.max(torch.abs(ky_test)))
 
     torch.testing.assert_close(traj.kx.to(torch.float32), kx_test.to(torch.float32), atol=1e-2, rtol=1e-3)
     torch.testing.assert_close(traj.ky.to(torch.float32), ky_test.to(torch.float32), atol=1e-2, rtol=1e-3)