Add support for v5.0 trajectories (#188)

Author: chaithyagr

docs/trajectory_gradspec.rst

@@ -10,7 +10,7 @@ The binary file format is specified as follows:
 +----------------+-------+---------+---------+------------------------------------------------------------------------+
 | Name           | Type  | Size    | Unit    | Description                                                            |
 +================+=======+=========+=========+========================================================================+
-| Version        | FLOAT | 1       | n.a.    | file version this new version would be “4.0”                           |
+| Version        | FLOAT | 1       | n.a.    | file version this new version would be “5.0”                           |
 +----------------+-------+---------+---------+------------------------------------------------------------------------+
 | Dimension      | FLOAT | 1       | n.a.    | 2 -> 2D , 3 -> 3D                                                      |
 +----------------+-------+---------+---------+------------------------------------------------------------------------+
@@ -36,6 +36,8 @@ The binary file format is specified as follows:
 +----------------+-------+---------+---------+------------------------------------------------------------------------+
 | kStarts        | FLOAT | D*Nc    | 1/m     | K-space location start                                                 |
 +----------------+-------+---------+---------+------------------------------------------------------------------------+
+| kEnds          | FLOAT | D*Nc    | 1/m     | K-space location end points. Only in version>5.0                       |
++----------------+-------+---------+---------+------------------------------------------------------------------------+
 | Gradient array | FLOAT | D*Nc*Ns | unitary | Gradient trajectory expressed in the range [-1; 1] relative to MaxGrad |
 +----------------+-------+---------+---------+------------------------------------------------------------------------+
 

src/mrinufft/io/nsp.py

@@ -36,6 +36,7 @@ def write_gradients(
     recon_tag: float = 1.1,
     timestamp: float | None = None,
     keep_txt_file: bool = False,
+    final_positions: np.ndarray | None = None,
 ):
     """Create gradient file from gradients and initial positions.
 
@@ -66,6 +67,8 @@ def write_gradients(
     keep_txt_file : bool, optional
         Whether to keep the text file used temporarily which holds data pushed to
         binary file, by default False
+    final_positions : np.ndarray, optional
+        Final positions. Shape (num_shots, dimension), by default None
 
     """
     num_shots = gradients.shape[0]
@@ -125,6 +128,20 @@ def write_gradients(
         )
         + "\n"
     )
+    if version >= 5:
+        if final_positions is None:
+            warnings.warn(
+                "Final positions not provided for version >= 5,"
+                "calculating final positions from gradients"
+            )
+            final_positions = initial_positions + np.sum(gradients, axis=1)
+        file.write(
+            "\n".join(
+                " ".join([f"{iter2:5.4f}" for iter2 in iter1])
+                for iter1 in final_positions
+            )
+            + "\n"
+        )
     if version < 4.1:
         # Write the maximum Gradient
         file.write(str(max_grad) + "\n")
@@ -185,6 +202,7 @@ def write_trajectory(
     check_constraints: bool = True,
     gmax: float = DEFAULT_GMAX,
     smax: float = DEFAULT_SMAX,
+    version: float = 5,
     **kwargs,
 ):
     """Calculate gradients from k-space points and write to file.
@@ -212,17 +230,20 @@ def write_trajectory(
         Maximum gradient magnitude in T/m, by default 0.04
     smax : float, optional
         Maximum slew rate in T/m/ms, by default 0.1
+    version: float, optional
+        Trajectory versioning, by default 5
     kwargs : dict, optional
         Additional arguments for writing the gradient file.
         These are arguments passed to write_gradients function above.
     """
     # Convert normalized trajectory to gradients
-    gradients, initial_positions = convert_trajectory_to_gradients(
+    gradients, initial_positions, final_positions = convert_trajectory_to_gradients(
         trajectory,
         norm_factor=norm_factor,
         resolution=np.asarray(FOV) / np.asarray(img_size),
         raster_time=raster_time,
         gamma=gamma,
+        get_final_positions=True,
     )
 
     # Check constraints if requested
@@ -245,10 +266,12 @@ def write_trajectory(
     write_gradients(
         gradients=gradients,
         initial_positions=initial_positions,
+        final_positions=final_positions,
         grad_filename=grad_filename,
         img_size=img_size,
         FOV=FOV,
         gamma=gamma,
+        version=version,
         **kwargs,
     )
 
@@ -326,6 +349,9 @@ def read_trajectory(
             _, data = _pop_elements(data, left_over)
         initial_positions, data = _pop_elements(data, dimension * num_shots)
         initial_positions = np.reshape(initial_positions, (num_shots, dimension))
+        if version >= 5:
+            final_positions, data = _pop_elements(data, dimension * num_shots)
+            final_positions = np.reshape(final_positions, (num_shots, dimension))
         dwell_time_ns = dwell_time * 1e6
         gradient_raster_time_ns = raster_time * 1e6
         if version < 4.1:

src/mrinufft/trajectories/utils.py

@@ -207,6 +207,7 @@ def convert_trajectory_to_gradients(
     resolution=DEFAULT_RESOLUTION,
     raster_time=DEFAULT_RASTER_TIME,
     gamma=Gammas.HYDROGEN,
+    get_final_positions=False,
 ):
     """Derive a normalized trajectory over time to provide gradients.
 
@@ -227,6 +228,9 @@ def convert_trajectory_to_gradients(
     gamma : float, optional
         Gyromagnetic ratio of the selected nucleus in kHz/T
         The default is Gammas.HYDROGEN.
+    get_final_positions : bool, optional
+        If `True`, return the final positions in k-space.
+        The default is `False`.
 
     Returns
     -------
@@ -239,6 +243,8 @@ def convert_trajectory_to_gradients(
     # Compute gradients and starting positions
     gradients = np.diff(trajectory, axis=1) / gamma / raster_time
     initial_positions = trajectory[:, 0, :]
+    if get_final_positions:
+        return gradients, initial_positions, trajectory[:, -1, :]
     return gradients, initial_positions
 
 

tests/test_io.py

@@ -5,7 +5,7 @@
 from mrinufft.io.utils import add_phase_to_kspace_with_shifts
 from mrinufft.trajectories.trajectory2D import initialize_2D_radial
 from mrinufft.trajectories.trajectory3D import initialize_3D_cones
-from pytest_cases import parametrize_with_cases
+from pytest_cases import parametrize, parametrize_with_cases
 from case_trajectories import CasesTrajectories
 
 
@@ -40,8 +40,18 @@ def case_trajectory_3D(self):
     "name, trajectory, FOV, img_size, in_out, min_osf, gamma, recon_tag",
     cases=CasesIO,
 )
+@parametrize("version", [4.2, 5.0])
 def test_write_n_read(
-    name, trajectory, FOV, img_size, in_out, min_osf, gamma, recon_tag, tmp_path
+    name,
+    trajectory,
+    FOV,
+    img_size,
+    in_out,
+    min_osf,
+    gamma,
+    recon_tag,
+    tmp_path,
+    version,
 ):
     """Test function which writes the trajectory and reads it back."""
     write_trajectory(
@@ -51,15 +61,15 @@ def test_write_n_read(
         check_constraints=True,
         grad_filename=str(tmp_path / name),
         in_out=in_out,
-        version=4.2,
+        version=version,
         min_osf=min_osf,
         recon_tag=recon_tag,
         gamma=gamma,
     )
     read_traj, params = read_trajectory(
         str((tmp_path / name).with_suffix(".bin")), gamma=gamma, read_shots=True
     )
-    assert params["version"] == 4.2
+    assert params["version"] == version
     assert params["num_shots"] == trajectory.shape[0]
     assert params["num_samples_per_shot"] == trajectory.shape[1] - 1
     assert params["TE"] == (0.5 if in_out else 0)