Skip to content

Commit 8319832

Browse files
paquiteaupaquiteau
committed
setup GRE sequence from arbitrary trajectories.
1 parent 86140df commit 8319832

File tree

2 files changed

+278
-8
lines changed

2 files changed

+278
-8
lines changed

‎src/mrinufft/io/pulseq.py‎

Lines changed: 169 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -5,19 +5,22 @@
55
sequences. Requires the `pypulseq` package to be installed.
66
"""
77

8+
import warnings
9+
from types import SimpleNamespace
10+
from mrinufft.trajectories.tools import get_gradient_amplitudes_to_travel_for_set_time
811
import numpy as np
912
import pypulseq as pp
13+
from numpy.typing import NDArray
14+
from .utils import prepare_trajectory_for_seq
1015

1116

12-
def read_pulseq_traj(filename):
17+
def read_pulseq_traj(filename, traj_delay=0.0, grad_offset=0.0):
1318
"""Extract k-space trajectory from a Pulseq sequence file.
1419
1520
The sequence should be a valid Pulseq `.seq` file, with arbitrary gradient
1621
waveforms, which all have the same length.
1722
18-
Unlike `Sequence.calculate_kspace`, this function returns the k-space
19-
trajectory segmented in shots ("blocks" from Pulseq), and works directly
20-
from the gradients waveforms stored in the `.seq` file.
23+
Note
2124
2225
Parameters
2326
----------
@@ -37,5 +40,166 @@ def read_pulseq_traj(filename):
3740
seq = pp.Sequence()
3841
seq.read(filename)
3942

43+
kspace_adc, kspace, t_exc, t_refocus, t_adc = seq.calculate_kspace(
44+
traj_delay, grad_offset
45+
)
4046

41-
gradient_waveforms = seq.waveforms()
47+
# split t_adc with t_exc and t_refocus, the index are then used to split kspace_adc
48+
49+
t_splits = sorted(np.concatenate([t_exc, t_refocus]))
50+
idx_prev = 0
51+
kspace_shots = []
52+
for t in t_splits:
53+
idx_next = np.searchsorted(t_adc, t, side="left")
54+
if idx_next == idx_prev:
55+
continue
56+
if idx_next == len(kspace_adc) and t > t_adc[-1]: # last useful point
57+
break
58+
kspace_shots.append(kspace_adc[idx_prev:idx_next, :])
59+
idx_prev = idx_next
60+
kspace_shots = np.asarray(kspace_shots)
61+
return kspace_shots
62+
63+
64+
def _check_timings(seq):
65+
# Check whether the timing of the sequence is correct
66+
ok, error_report = seq.check_timing()
67+
if ok:
68+
return None
69+
else:
70+
warnings.warn("Timing check failed. Error listing follows:" + str(error_report))
71+
72+
73+
def gre_3D(
74+
trajectory: NDArray,
75+
fov: tuple[float, float, float],
76+
img_size: tuple[int, int, int],
77+
TR: float,
78+
TE: float,
79+
TE_pos: float = 0.5,
80+
FA: float | None = None,
81+
rf_pulse: SimpleNamespace | None =None,
82+
rf_spoiling_inc: float = 0.0,
83+
system: pp.Opts = pp.Opts.default,
84+
osf: int = 1,
85+
) -> pp.Sequence:
86+
"""Create a Pulseq 3D-GRE sequence with arbitrary gradient waveform designed
87+
to follow a given k-space trajectory.
88+
89+
The sequence is designed to
90+
The
91+
92+
Parameters
93+
----------
94+
trajectory : np.ndarray
95+
The k-space trajectory as a numpy array of shape (n_shots, n_samples, 3),
96+
where the last dimension corresponds to the x, y, and z coordinates in k-space.
97+
TR: float
98+
The repetition time in milliseconds (ms).
99+
TE: float
100+
The echo time in milliseconds (ms).
101+
FA: float, optional, incompatible with `rf_pulse`
102+
The flip angle in degrees (°).
103+
TE_pos: float, optional
104+
The relative (0-1) position of the echo time within each kspace_shot.
105+
rf_pulse: SimpleNamespace, optional, incompatible with `FA`
106+
A custom radio-frequency pulse object. If not provided, a block pulse with the specified flip angle and a duration of 4 ms will be created. see `pypulseq.make_block_pulse` or `pypulseq.make_arbitrary_rf` for more details.
107+
rf_spoiling_inc: float, optional
108+
The increment in the RF phase (in degree) for spoiling. Default is 0.0, which means no spoiling.
109+
110+
osf: int, optional
111+
The oversampling factor for the ADC. Default is 1, which means no oversampling.
112+
113+
system : pypulseq.Opts, optional
114+
The system options for the Pulseq sequence. Default is `pp.Opts.default`.
115+
116+
Notes
117+
-----
118+
The Sequence cycle can be summarized as follows:
119+
1. RF pulse
120+
2. Delay to sync TE
121+
3. Gradients plays: The gradients consist in a prewind to the first point of the trajectory, the trajectory itself, and a rewind to the edge of k-space.
122+
3bis. The ADC is opened on the trajectory points (ignoring the prewind and rewinds parts)
123+
4. Gradient spoilers
124+
5. Delay to sync the next TR
125+
126+
Returns
127+
-------
128+
pp.Sequence
129+
A Pulseq sequence object with the specified arbitrary gradient waveform.
130+
131+
132+
See Also
133+
--------
134+
135+
"""
136+
TR = TR / 1000.0 # convert to seconds
137+
TE = TE / 1000.0 # convert to seconds
138+
seq = pp.Sequence(system=system)
139+
140+
if rf_pulse is None and FA is not None:
141+
rf_pulse = pp.make_block_pulse(flip_angle=FA, system=system, use="excitation")
142+
elif rf_pulse is not None and FA is not None:
143+
raise ValueError(
144+
"Cannot specify both `rf_pulse` and `FA`. Please provide only one."
145+
)
146+
elif rf_pulse is None and FA is None:
147+
raise ValueError("Either `rf_pulse` or `FA` must be provided.")
148+
if not isinstance(rf_pulse, SimpleNamespace):
149+
raise TypeError(
150+
"The `rf_pulse` parameter must be a SimpleNamespace object, "
151+
"as returned by `pypulseq.make_block_pulse` or `pypulseq.make_arbitrary_rf`."
152+
)
153+
154+
full_grads, skip_start, skip_end = prepare_trajectory_for_seq(trajectory, fov=fov, img_size=img_size, raster_time=system.grad_raster_time, gamma=system.gamma, gmax=system.max_grad, smax=system.max_slew, pregrad="prephase", postgrad="slowdown_to_edge")
155+
traj_length = full_grads.shape[1] - skip_start - skip_end
156+
157+
shot_duration = system.grad_raster_time * traj_length # in seconds
158+
adc = pp.make_adc(num_samples=traj_length*osf, system=system, duration=shot_duration, delay=skip_start * system.grad_raster_time)
159+
160+
# Add a spoiler gradient that move to twice the edge of k-space in the x direction.
161+
spoiler = pp.make_trapezoid(channel="x", area=2*img_size[0]/fov[0], system=system)
162+
163+
delay_before_grad = pp.make_delay(TE - pp.calc_rf_center(rf_pulse)[0] - rf_pulse.delay - TE_pos*shot_duration - skip_start*system.grad_raster_time)
164+
delay_end_TR = pp.make_delay(TR - rf_pulse.duration - delay_before_grad - full_grads.shape[1] * system.grad_raster_time)
165+
166+
rf_phase = 0.0
167+
for grad_xyz in full_grads:
168+
rf_phase = divmod(rf_phase + rf_spoiling_inc, 360.0)[1]
169+
rf_pulse.phase_offset = rf_phase / 180 * np.pi
170+
adc.phase_offset = rf_phase / 180 * np.pi
171+
seq.add_block(rf_pulse) # RF pulse
172+
seq.add_block(delay_before_grad) # delay to sync TE
173+
seq.add_block(*[pp.make_arbitrary_grad(channel=c, waveform=grad_xyz[i], system=system) for i,c in enumerate("xyz")], adc) # pause for tuning echo time
174+
seq.add_block(spoiler, delay_end_TR)
175+
176+
_check_timings(seq)
177+
return seq
178+
179+
180+
def gre_2D(trajectory, TR, TE, FA, pulse, system=pp.Opts.default):
181+
"""Create a Pulseq 2D-GRE sequence with arbitrary gradient waveform designed
182+
to follow a given k-space trajectory.
183+
184+
Parameters
185+
----------
186+
trajectory : np.ndarray
187+
The k-space trajectory as a numpy array of shape (n_shots, n_samples, 3),
188+
where the last dimension corresponds to the x, y, and z coordinates in k-space.
189+
TR: float
190+
The repetition time in milliseconds (ms).
191+
TE: float
192+
The echo time in milliseconds (ms).
193+
FA: float
194+
The flip angle in degrees (°).
195+
pulse: SimpleNamespace
196+
A custom radio-frequency pulse object.
197+
system : pypulseq.Opts, optional
198+
The system options for the Pulseq sequence. Default is `pp.Opts.default`.
199+
200+
Returns
201+
-------
202+
pp.Sequence
203+
A Pulseq sequence object with the specified arbitrary gradient waveform.
204+
"""
205+
...

‎src/mrinufft/io/utils.py‎

Lines changed: 109 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,22 @@
11
"""Module containing utility functions for IO in MRI NUFFT."""
22

33
import numpy as np
4+
from numpy.typing import NDArray
45
from scipy.spatial.transform import Rotation
6+
from ..trajectories.utils import (
7+
convert_trajectory_to_gradients,
8+
Gammas,
9+
DEFAULT_SMAX,
10+
DEFAULT_GMAX,
11+
DEFAULT_RASTER_TIME,
12+
KMAX,
13+
)
14+
from ..trajectories.tools import get_gradient_amplitudes_to_travel_for_set_time
515

616

7-
def add_phase_to_kspace_with_shifts(kspace_data, kspace_loc, normalized_shifts):
17+
def add_phase_to_kspace_with_shifts(
18+
kspace_data: NDArray, kspace_loc: NDArray, normalized_shifts: NDArray
19+
):
820
"""
921
Add phase shifts to k-space data.
1022
@@ -38,7 +50,7 @@ def add_phase_to_kspace_with_shifts(kspace_data, kspace_loc, normalized_shifts):
3850
return kspace_data * phase
3951

4052

41-
def siemens_quat_to_rot_mat(quat, return_det=False):
53+
def siemens_quat_to_rot_mat(quat: NDArray, return_det=False):
4254
"""
4355
Calculate the rotation matrix from Siemens Twix quaternion.
4456
@@ -134,7 +146,7 @@ def nifti_affine(twixObj):
134146
return full_mat
135147

136148

137-
def remove_extra_kspace_samples(kspace_data, num_samples_per_shot):
149+
def remove_extra_kspace_samples(kspace_data: NDArray, num_samples_per_shot: int):
138150
"""Remove extra samples from k-space data.
139151
140152
This function is useful when the k-space data has extra samples
@@ -159,3 +171,97 @@ def remove_extra_kspace_samples(kspace_data, num_samples_per_shot):
159171
if n_extra_samples > 0:
160172
kspace_data = kspace_data[..., :-n_extra_samples]
161173
return kspace_data
174+
175+
176+
def prepare_trajectory_for_seq(
177+
trajectory: NDArray,
178+
fov: tuple[float, float, float],
179+
img_size: tuple[int, int, int],
180+
norm_factor: float = KMAX,
181+
pregrad: str = "prephase",
182+
postgrad: str = "slowdown_to_edge",
183+
gamma=Gammas.HYDROGEN,
184+
raster_time=DEFAULT_RASTER_TIME,
185+
gmax=DEFAULT_GMAX,
186+
smax=DEFAULT_SMAX,
187+
):
188+
"""Prepare gradients from trajectory.
189+
190+
Parameters
191+
----------
192+
trajectory : np.ndarray
193+
The k-space trajectory as a numpy array of shape (n_shots, n_samples, 3),
194+
where the last dimension corresponds to the x, y, and z coordinates in k-space.
195+
norm_factor : float
196+
The normalization factor for the trajectory. (default is 0.5)
197+
fov : tuple[float, float, float]
198+
The field of view in the x, y, and z dimensions, in meters.
199+
img_size : tuple[int, int, int]
200+
The image size in the x, y, and z dimensions, in pixels.
201+
pregrad : str, optional
202+
The type of pre-gradient to apply. ONly "prephase" is supported currently.
203+
postgrad : str, optional
204+
The type of post-gradient to apply. Defaults to "slowdown_to_edge".
205+
206+
Returns
207+
-------
208+
np.ndarray
209+
The full gradients as a numpy array of shape (n_shots, n_samples, 3),
210+
where the last dimension corresponds to the x, y, and z gradient amplitudes.
211+
int
212+
The number of samples to skip at the start of the trajectory.
213+
int
214+
The number of samples to skip at the end of the trajectory.
215+
"""
216+
# from #276 : We need to prewind the gradients to the first point of the
217+
# trajectory, and rewind them to the edge of k-space.
218+
219+
# We will move from init_pos -[prewind]-> start_pos -> trajectory -> end_pos -[postgrad]-> final_pos
220+
221+
222+
grads, start_pos, end_pos = convert_trajectory_to_gradients(
223+
trajectory,
224+
norm_factor=0.5,
225+
resolution=fov,
226+
raster_time=raster_time,
227+
gamma=gamma,
228+
get_final_positions=True,
229+
)
230+
231+
# prewind the gradients to their first point:
232+
if pregrad == "prephase":
233+
init_pos = np.zeros_like(start_pos)
234+
else:
235+
raise ValueError("Only 'prephase' is supported for pregrad.")
236+
start_grads = get_gradient_amplitudes_to_travel_for_set_time(
237+
kspace_end_loc=start_pos,
238+
kspace_start_loc=init_pos,
239+
end_gradients=grads[:, 0, :],
240+
gamma=gamma,
241+
raster_time=raster_time,
242+
gmax=gmax,
243+
smax=smax,
244+
)
245+
skip_start = start_grads.shape[1]
246+
247+
final_pos = np.zeros_like(end_pos)
248+
if postgrad == "slowdown_to_edge":
249+
# Set the edge location to [Kmax, 0,0], to prepare for gradient spoiling.
250+
final_pos[..., 0] = img_size[0] / fov[0] / 2
251+
else:
252+
raise ValueError("Only 'slowdown_to_edge' is supported for postgrad.")
253+
254+
end_grads = get_gradient_amplitudes_to_travel_for_set_time(
255+
kspace_start_loc=end_pos,
256+
kspace_end_loc=final_pos,
257+
start_gradients=grads[:, -1, :],
258+
gamma=gamma,
259+
raster_time=raster_time,
260+
gmax=gmax,
261+
smax=smax,
262+
)
263+
264+
skip_end = end_grads.shape[1]
265+
full_gradients = np.hstack([start_grads, grads, end_grads])
266+
267+
return full_gradients, skip_start, skip_end

0 commit comments

Comments
 (0)