Add prephasors, rewinders etc.

docs/trajectory_gradspec.rst

@@ -32,7 +32,11 @@ The binary file format is specified as follows:
 +----------------+-------+---------+---------+------------------------------------------------------------------------+
 | timestamp      | FLOAT | 1       | n.a.    | Time stamp when the binary is created                                  |
 +----------------+-------+---------+---------+------------------------------------------------------------------------+
-| Empty places   | FLOAT | 9       | n.a.    | Yet unused : Default initialized with 0                                |
+| ADC pre-skip   | FLOAT | 1       | n.a.    | Gradient samples to skip before starting ADC, for pre-phasors          |
++----------------+-------+---------+---------+------------------------------------------------------------------------+
+| ADC post-skip  | FLOAT | 1       | n.a.    | Gradient samples to skip at the end of trajectory by turning off ADC   |
++----------------+-------+---------+---------+------------------------------------------------------------------------+
+| Empty places   | FLOAT | 7       | n.a.    | Yet unused : Default initialized with 0                                |
 +----------------+-------+---------+---------+------------------------------------------------------------------------+
 | kStarts        | FLOAT | D*Nc    | 1/m     | K-space location start                                                 |
 +----------------+-------+---------+---------+------------------------------------------------------------------------+

pyproject.toml

@@ -4,7 +4,7 @@ description = "MRI Non-Cartesian Fourier Operators with multiple computation bac
 authors = [{name="Pierre-antoine Comby", email="[email protected]"}]
 
 readme = "README.md"
-dependencies = ["numpy", "scipy", "matplotlib", "tqdm"]
+dependencies = ["numpy", "scipy", "matplotlib", "tqdm", "joblib"]
 requires-python = ">=3.9"
 
 dynamic = ["version"]

src/mrinufft/io/nsp.py

@@ -8,6 +8,7 @@
 from datetime import datetime
 
 import numpy as np
+from typing import Optional
 
 from mrinufft.trajectories.utils import (
     DEFAULT_GMAX,
@@ -19,6 +20,7 @@
     convert_gradients_to_slew_rates,
     convert_trajectory_to_gradients,
 )
+from mrinufft.trajectories.tools import get_gradient_amplitudes_to_travel_for_set_time
 
 from .siemens import read_siemens_rawdat
 
@@ -29,14 +31,16 @@ def write_gradients(
     grad_filename: str,
     img_size: tuple[int, ...],
     FOV: tuple[float, ...],
-    in_out: bool = True,
+    TE_pos: float = 0.5,
     min_osf: int = 5,
     gamma: float = Gammas.HYDROGEN,
     version: float = 4.2,
     recon_tag: float = 1.1,
     timestamp: float | None = None,
     keep_txt_file: bool = False,
     final_positions: np.ndarray | None = None,
+    start_skip_samples: int = 0,
+    end_skip_samples: int = 0,
 ):
     """Create gradient file from gradients and initial positions.
 
@@ -52,8 +56,10 @@ def write_gradients(
         Image size.
     FOV : tuple[float, ...]
         Field of view.
-    in_out : bool, optional
-        Whether it is In-Out trajectory?, by default True
+    TE_pos : float, optional
+        The ratio of trajectory when TE occurs, with 0 as start of
+        trajectory and 1 as end. By default 0.5, which is the
+        center of the trajectory (in-out trajectory).
     min_osf : int, optional
         Minimum oversampling factor needed at ADC, by default 5
     gamma : float, optional
@@ -69,6 +75,12 @@ def write_gradients(
         binary file, by default False
     final_positions : np.ndarray, optional
         Final positions. Shape (num_shots, dimension), by default None
+    start_skip_samples : int, optional
+        Number of samples to skip in ADC at start of each shot, by default 0
+        This works only for version >= 5.1.
+    end_skip_samples : int, optional
+        Number of samples to skip in ADC at end of each shot, by default 0
+        This works only for version >= 5.1.
 
     """
     num_shots = gradients.shape[0]
@@ -100,14 +112,12 @@ def write_gradients(
     file.write(str(num_shots) + "\n")
     file.write(str(num_samples_per_shot) + "\n")
     if version >= 4.1:
-        if not in_out:
+        if TE_pos == 0:
             if np.sum(initial_positions) != 0:
                 warnings.warn(
                     "The initial positions are not all zero for center-out trajectory"
                 )
-            file.write("0\n")
-        else:
-            file.write("0.5\n")
+        file.write(str(TE_pos) + "\n")
         # Write the maximum Gradient
         file.write(str(max_grad) + "\n")
         # Write recon Pipeline version tag
@@ -119,6 +129,10 @@ def write_gradients(
                 timestamp = float(datetime.now().timestamp())
             file.write(str(timestamp) + "\n")
             left_over -= 1
+        if version >= 5.1:
+            file.write(str(start_skip_samples) + "\n")
+            file.write(str(end_skip_samples) + "\n")
+            left_over -= 2
         file.write(str("0\n" * left_over))
     # Write all the k0 values
     file.write(
@@ -165,7 +179,7 @@ def write_gradients(
         os.remove(grad_filename + ".txt")
 
 
-def _pop_elements(array, num_elements=1, type="float"):
+def _pop_elements(array, num_elements=1, type=np.float32):
     """Pop elements from an array.
 
     Parameters
@@ -200,8 +214,11 @@ def write_trajectory(
     gamma: float = Gammas.HYDROGEN,
     raster_time: float = DEFAULT_RASTER_TIME,
     check_constraints: bool = True,
+    TE_pos: float = 0.5,
     gmax: float = DEFAULT_GMAX,
     smax: float = DEFAULT_SMAX,
+    pregrad: str | None = None,
+    postgrad: str | None = None,
     version: float = 5,
     **kwargs,
 ):
@@ -226,10 +243,26 @@ def write_trajectory(
         Gradient raster time in ms, by default 0.01
     check_constraints : bool, optional
         Check scanner constraints, by default True
+    TE_pos : float, optional
+        The ratio of trajectory when TE occurs, with 0 as start of
+        trajectory and 1 as end. By default 0.5, which is the
+        center of the trajectory (in-out trajectory).
     gmax : float, optional
         Maximum gradient magnitude in T/m, by default 0.04
     smax : float, optional
         Maximum slew rate in T/m/ms, by default 0.1
+    pregrad : str, optional
+        Pregrad method, by default `prephase`
+        `prephase` will add a prephasing gradient to the start of the trajectory.
+    postgrad : str, optional
+        Postgrad method, by default 'slowdown_to_edge'
+        `slowdown_to_edge` will add a gradient to slow down to the edge of the k-space
+        along x-axis for all the shots i.e. go to (Kmax, 0, 0).
+        This is useful for sequences needing a spoiler at the end of the trajectory.
+        However, spoiler is still not added, it is expected that the sequence
+        handles the spoilers, which can be variable.
+        `slowdown_to_center` will add a gradient to slow down to the center
+        of the k-space.
     version: float, optional
         Trajectory versioning, by default 5
     kwargs : dict, optional
@@ -245,7 +278,46 @@ def write_trajectory(
         gamma=gamma,
         get_final_positions=True,
     )
-
+    Ns_to_skip_at_start = 0
+    Ns_to_skip_at_end = 0
+    scan_consts = {
+        "gamma": gamma,
+        "gmax": gmax,
+        "smax": smax,
+        "raster_time": raster_time,
+    }
+    if pregrad == "prephase":
+        if version < 5.1:
+            raise ValueError(
+                "pregrad is only supported for version >= 5.1, "
+                "please set version to 5.1 or higher."
+            )
+        start_gradients = get_gradient_amplitudes_to_travel_for_set_time(
+            kspace_end_loc=initial_positions,
+            end_gradients=gradients[:, 0],
+            **scan_consts,
+        )
+        initial_positions = np.zeros_like(initial_positions)
+        gradients = np.hstack([start_gradients, gradients])
+        Ns_to_skip_at_start = start_gradients.shape[1]
+    if postgrad:
+        if version < 5.1:
+            raise ValueError(
+                "postgrad is only supported for version >= 5.1, "
+                "please set version to 5.1 or higher."
+            )
+        edge_locations = np.zeros_like(final_positions)
+        if postgrad == "slowdown_to_edge":
+            # Always end at KMax, the spoilers can be handeled by the sequence.
+            edge_locations[..., 0] = img_size[0] / FOV[0] / 2
+        end_gradients = get_gradient_amplitudes_to_travel_for_set_time(
+            kspace_end_loc=edge_locations,
+            start_gradients=gradients[:, -1],
+            kspace_start_loc=final_positions,
+            **scan_consts,
+        )
+        gradients = np.hstack([gradients, end_gradients])
+        Ns_to_skip_at_end = end_gradients.shape[1]
     # Check constraints if requested
     if check_constraints:
         slewrates, _ = convert_gradients_to_slew_rates(gradients, raster_time)
@@ -261,6 +333,15 @@ def write_trajectory(
                 f"Maximum gradient amplitude: {maxG:.3f} > {gmax:.3f}"
                 f"Maximum slew rate: {maxS:.3f} > {smax:.3f}"
             )
+        if pregrad != "prephase":
+            border_slew_rate = gradients[:, 0] / raster_time
+            if np.any(np.abs(border_slew_rate) > smax):
+                warnings.warn(
+                    "Slew rate at start of trajectory exceeds maximum slew rate!"
+                    f"Maximum slew rate: {np.max(np.abs(border_slew_rate)):.3f}"
+                    f" > {smax:.3f}. Please use prephase gradient to avoid this "
+                    " issue."
+                )
 
     # Write gradients in file
     write_gradients(
@@ -270,16 +351,19 @@ def write_trajectory(
         grad_filename=grad_filename,
         img_size=img_size,
         FOV=FOV,
+        TE_pos=TE_pos,
         gamma=gamma,
         version=version,
+        start_skip_samples=Ns_to_skip_at_start,
+        end_skip_samples=Ns_to_skip_at_end,
         **kwargs,
     )
 
 
 def read_trajectory(
     grad_filename: str,
     dwell_time: float | str = DEFAULT_RASTER_TIME,
-    num_adc_samples: int = None,
+    num_adc_samples: int | None = None,
     gamma: Gammas | float = Gammas.HYDROGEN,
     raster_time: float = DEFAULT_RASTER_TIME,
     read_shots: bool = False,
@@ -336,25 +420,27 @@ def read_trajectory(
             if dwell_time == "min_osf":
                 dwell_time = raster_time / min_osf
         (num_shots, num_samples_per_shot), data = _pop_elements(data, 2, type="int")
-        if num_adc_samples is None:
-            if read_shots:
-                num_adc_samples = num_samples_per_shot + 1
-            else:
-                num_adc_samples = int(num_samples_per_shot * (raster_time / dwell_time))
-        if version >= 4.1:
+        if version > 4:
             TE, data = _pop_elements(data)
             grad_max, data = _pop_elements(data)
             recon_tag, data = _pop_elements(data)
             recon_tag = np.around(recon_tag, 2)
             left_over = 10
-            if version >= 4.2:
+            if version > 4.1:
                 timestamp, data = _pop_elements(data)
                 timestamp = datetime.fromtimestamp(float(timestamp))
                 left_over -= 1
+            if version > 5:
+                packed_skips, data = _pop_elements(data, num_elements=2, type="int")
+                start_skip_samples, end_skip_samples = packed_skips
+                left_over -= 2
+            else:
+                start_skip_samples = 0
+                end_skip_samples = 0
             _, 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:
+        if version > 4.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
@@ -365,12 +451,28 @@ def read_trajectory(
             data,
             dimension * num_samples_per_shot * num_shots,
         )
+        # Convert gradients to T/m
         gradients = np.reshape(
-            grad_max * gradients, (num_shots * num_samples_per_shot, dimension)
+            grad_max * gradients * 1e-3, (num_shots, num_samples_per_shot, dimension)
         )
-        # Convert gradients from mT/m to T/m
-        gradients = np.reshape(gradients * 1e-3, (-1, num_samples_per_shot, dimension))
-        kspace_loc = np.zeros((num_shots, num_adc_samples, dimension))
+        # Handle skipped samples
+        if start_skip_samples > 0:
+            start_location_updates = (
+                np.sum(gradients[:, :start_skip_samples], axis=1) * raster_time * gamma
+            )
+            initial_positions += start_location_updates
+            gradients = gradients[:, start_skip_samples:, :]
+            num_samples_per_shot -= start_skip_samples
+        if end_skip_samples > 0:
+            gradients = gradients[:, :-end_skip_samples, :]
+            num_samples_per_shot -= end_skip_samples
+        if num_adc_samples is None:
+            if read_shots:
+                # Acquire one extra sample at the end of each shot in read_shots mode
+                num_adc_samples = num_samples_per_shot + 1
+            else:
+                num_adc_samples = int(num_samples_per_shot * (raster_time / dwell_time))
+        kspace_loc = np.zeros((num_shots, num_adc_samples, dimension), dtype=np.float32)
         kspace_loc[:, 0, :] = initial_positions
         adc_times = dwell_time_ns * np.arange(1, num_adc_samples)
         Q, R = divmod(adc_times, gradient_raster_time_ns)