wip

Author: paquiteau

src/mrinufft/io/pulseq.py

@@ -5,19 +5,20 @@
 sequences. Requires the `pypulseq` package to be installed.
 """
 
+from types import SimpleNamespace
 import numpy as np
 import pypulseq as pp
+from numpy.typing import NDArray
+from ..trajectories.utils import convert_trajectory_to_gradients
 
 
-def read_pulseq_traj(filename):
+def read_pulseq_traj(filename, traj_delay=0.0, grad_offset=0.0):
     """Extract k-space trajectory from a Pulseq sequence file.
 
     The sequence should be a valid Pulseq `.seq` file, with arbitrary gradient
     waveforms, which all have the same length.
 
-    Unlike `Sequence.calculate_kspace`, this function returns the k-space
-    trajectory segmented in shots ("blocks" from Pulseq), and works directly
-    from the gradients waveforms stored in the `.seq` file.
+    Note
 
     Parameters
     ----------
@@ -37,5 +38,90 @@ def read_pulseq_traj(filename):
     seq = pp.Sequence()
     seq.read(filename)
 
+    kspace_adc, kspace, t_exc, t_refocus, t_adc = seq.calculate_kspace(
+        traj_delay, grad_offset
+    )
 
-    gradient_waveforms = seq.waveforms()
+    # split t_adc with t_exc and t_refocus, the index are then used to split kspace_adc
+
+    t_splits = sorted(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
+        if idx_next == len(kspace_adc) and t > t_adc[-1]: # last useful point
+            break
+        kspace_shots.append(kspace_adc[idx_prev:idx_next, :])
+        idx_prev= idx_next
+    kspace_shots = np.asarray(kspace_shots)
+    return kspace_shots
+
+
+def pulseq_grad_blocks(
+    trajectory: NDArray, system: pp.Opts = pp.Opts.default
+) -> list[SimpleNamespace]:
+    """Create Pulseq gradient blocks from a k-space trajectory.
+
+    Parameters
+    ----------
+    trajectory : np.ndarray
+        The k-space trajectory as a numpy array of shape (n_shots, n_samples, 3),
+        where the last dimension corresponds to the x, y, and z coordinates in k-space.
+
+    Returns
+    -------
+    list
+        A list of Pulseq gradient blocks corresponding to the k-space trajectory.
+    """
+    grads = convert_trajectory_to_gradients(trajectory)
+
+    # TODO add prewind/postwind stuff from #276
+    # TODO ensure amplitudes are correct
+
+    # TODO deduplicate gradients (reduces the size of the file, speed up the sequence loading).
+    # TODO Are rotation of trajectories supported in pulseq?
+
+    grads = convert_trajectory_to_gradients(trajectory)
+    for g in grads:
+        for i, c in enumerate(["x", "y", "z"]):
+            pp_g = pp.make_arbitrary_grad("x", g[:, i])
+
+
+    return blocks_pp
+
+
+def pulseq_gre_arb_grad(trajectory, pulse, system: pp.Opts = pp.Opts.default):
+    """Create a Pulseq GRE sequence with arbitrary gradient waveform.
+
+    Parameters
+    ----------
+    trajectory : np.ndarray
+        The k-space trajectory as a numpy array of shape (n_shots, n_samples, 3),
+        where the last dimension corresponds to the x, y, and z coordinates in k-space.
+    pulse : pypulseq.Pulse
+
+    system : pypulseq.Opts, optional
+        The system options for the Pulseq sequence. Default is `pp.Opts.default`.
+
+    Returns
+    -------
+    pp.Sequence
+        A Pulseq sequence object with the specified arbitrary gradient waveform.
+
+
+    See Also
+    --------
+
+    """
+    seq = pp.Sequence(system=system)
+    grad_shot_ids, grad_shots = pulseq_grad_blocks(trajectory, system)
+
+
+    # TODO
+    for grad_shot_id in grad_shot_ids:
+        pp.add_block() # RF pulse
+        pp.add_block() # pause for tuning echo time
+        pp.add_block() # Gradient block with ADC
+    return seq