Pulseq

examples/example_pulseq.py

@@ -0,0 +1,92 @@
+# %%
+"""
+==================================
+Create a GRE Sequence using Pulseq
+==================================
+
+Example how to create sequences using PyPulseq.
+"""
+import numpy as np
+from mrinufft.io import pulseq_gre, prepare_trajectory_for_seq
+from mrinufft.trajectories.display import display_3D_trajectory
+from mrinufft.trajectories import initialize_2D_spiral, stack
+
+import matplotlib.pyplot as plt
+
+# %%
+# Defining the sequence parameters like repetition time (TR), echo time (TE), flip angle (FA), field-of-view (FOV) and image matrix size.
+
+TR = 100  # ms
+TE = 50  # ms
+FA = 10  # degrees
+FOV = np.array([0.192, 0.192, 0.128])  # Field of View in meters
+img_size = np.array([64, 64, 48])  # Image size in pixels
+
+
+# %%
+# Create a stack of spiral for our trajectory
+
+traj = stack(
+    initialize_2D_spiral(Nc=1, Ns=4096, nb_revolutions=12, in_out=True)[:, ::-1, :],
+    nb_stacks=img_size[-1],
+)
+
+display_3D_trajectory(traj)
+
+# %%
+grads, Ns, Ne = prepare_trajectory_for_seq(traj, fov=FOV, img_size=img_size)
+
+# %%
+# Each gradient shot has now a prephaser to move to the first point of the trajectory, and a rewind to the edge of k-space.
+
+plt.plot(grads[0, 0, :], label="Gx")
+plt.plot(grads[0, 1, :], label="Gy")
+plt.plot(grads[0, 2, :], label="Gz")
+
+# %%
+
+# Create a 3D GRE sequence with the trajectory:
+
+seq = pulseq_gre(traj[:3], fov=FOV, img_size=img_size, TR=TR, TE=TE, FA=FA)
+
+
+# %%
+
+# Let's show the sequence.
+plt.rcParams["figure.figsize"] = (20, 10)
+seq.plot(show_blocks=True, grad_disp="mT/m")
+
+# %%
+from mrinufft.io import read_pulseq_traj
+
+# %%
+read_kspace = read_pulseq_traj(seq)
+
+# %%
+KMAX = 0.5
+
+# %%
+kspace_adc, _, t_exc, t_refocus, t_adc = seq.calculate_kspace()
+
+# split t_adc with t_exc and t_refocus, the index are then used to split kspace_adc
+FOV = seq.get_definition("FOV")
+t_splits = np.sort(np.concatenate([t_exc, t_refocus]))
+idx_prev = 0
+kspace_shots = []
+for t in t_splits:
+    idx_next = np.searchsorted(t_adc, t, side="left")
+    if idx_next == idx_prev:
+        continue
+    kspace_shots.append(kspace_adc[:, idx_prev:idx_next].T)
+    if idx_next == kspace_adc.shape[1] and t > t_adc[-1]:  # last useful point
+        break
+    idx_prev = idx_next
+if idx_next < kspace_adc.shape[1]:
+    kspace_shots.append(kspace_adc[:, idx_next:].T)  # add remaining gradients.
+# convert to KMAX standard.
+kspace_shots = np.ascontiguousarray(kspace_shots) * KMAX * 2 * np.asarray(FOV)
+
+
+# %%
+# and we can display the k-space trajectory, loaded from the sequence file.
+display_3D_trajectory(kspace_shots)

pyproject.toml

@@ -27,7 +27,7 @@ pynufft = ["pynufft"]
 pynufft-cpu = ["pynufft"]
 pynufft-gpu = ["pynufft"]
 
-io = ["pymapvbvd"]
+io = ["pymapvbvd", "pypulseq"]
 smaps = ["scikit-image"]
 extra = ["pymapvbvd", "scikit-image", "scikit-learn", "pywavelets"]
 autodiff = ["torch"]

src/mrinufft/io/__init__.py

@@ -3,13 +3,17 @@
 from .cfl import traj2cfl, cfl2traj
 from .nsp import read_trajectory, write_trajectory, read_arbgrad_rawdat
 from .siemens import read_siemens_rawdat
-
+from .pulseq import read_pulseq_traj, pulseq_gre
+from .utils import prepare_trajectory_for_seq
 
 __all__ = [
-    "traj2cfl",
     "cfl2traj",
-    "read_trajectory",
-    "write_trajectory",
+    "prepare_trajectory_for_seq",
+    "pulseq_gre",
     "read_arbgrad_rawdat",
+    "read_pulseq_traj",
     "read_siemens_rawdat",
+    "read_trajectory",
+    "traj2cfl",
+    "write_trajectory",
 ]

src/mrinufft/io/nsp.py

@@ -8,6 +8,7 @@
 from array import array
 from datetime import datetime
 
+from mrinufft.io.utils import prepare_trajectory_for_seq
 import numpy as np
 
 from mrinufft.trajectories.utils import (
@@ -210,8 +211,8 @@ def _pop_elements(array, num_elements=1, type=np.float32):
 
 def write_trajectory(
     trajectory: np.ndarray,
-    FOV: tuple[float, ...],
-    img_size: tuple[int, ...],
+    FOV: tuple[float, float, float],
+    img_size: tuple[int, int, int],
     grad_filename: str,
     norm_factor: float = KMAX,
     gamma: float = Gammas.HYDROGEN / 1e3,
@@ -311,14 +312,15 @@ def write_trajectory(
         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,
-            acq=acq,
-        )
-        gradients = np.hstack([gradients, end_gradients])
-        Ns_to_skip_at_end = end_gradients.shape[1]
+            end_gradients = get_gradient_amplitudes_to_travel_for_set_time(
+                kspace_end_loc=edge_locations,
+                start_gradients=gradients[:, -1],
+                kspace_start_loc=final_positions,
+                acq=acq,
+            )
+
+            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, acq)