Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix sigpy tests #227

Merged
merged 3 commits into from
Feb 13, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 22 additions & 4 deletions src/pypulseq/make_sigpy_pulse.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,16 +45,14 @@ def sigpy_n_seq(
----------
flip_angle : float
Flip angle in radians.
apodization : float, optional, default=0
Apodization.
center_pos : float, optional, default=0.5
Position of peak.5 (midway).
delay : float, optional, default=0
Delay in seconds (s).
duration : float, optional, default=4e-3
Duration in seconds (s).
freq_offset : float, optional, default=0
Frequency offset in Hertz (Hz).
center_pos : float, optional, default=0.5
Position of peak.5 (midway).
max_grad : float, optional, default=0
Maximum gradient strength of accompanying slice select trapezoidal event.
max_slew : float, optional, default=0
Expand All @@ -70,6 +68,26 @@ def sigpy_n_seq(
System limits. Default is a system limits object initialized to default values.
time_bw_product : float, optional, default=4
Time-bandwidth product.
pulse_cfg: SigpyPulseOpts, optional, default=None
Pulse configuration options. Possible keys are:
- pulse_type: str, optional, default='slr'
Pulse type. Must be one of 'slr' or 'sms'.
- ptype: str, optional, default='st'
Pulse design method. Must be one of 'st', 'ex', 'inv', 'sat', 'se', 'fi', 'fs', 'se'.
- ftype: str, optional, default='ls'
Filter type. Must be one of 'ls', 'pm', 'min', 'max', 'ap'.
- d1: float, optional, default=0.01
Passband ripple.
- d2: float, optional, default=0.01
Stopband ripple.
- cancel_alpha_phs: bool, optional, default=False
Cancel alpha phase.
- n_bands: int, optional, default=3
Number of bands. SMS only.
- band_sep: float, optional, default=20
Band separation. SMS only.
- phs_0_pt: str, optional, default='None'
Phase 0 point. SMS only.
use : str, optional, default=str()
Use of radio-frequency sinc pulse. Must be one of 'excitation', 'refocusing' or 'inversion'.
plot: bool, optional, default=True
Expand Down
89 changes: 48 additions & 41 deletions tests/test_sigpy.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
# sms - check MB
# slr - check slice profile
"""Tests for SigPy pulse generation."""

import importlib.util

Expand All @@ -21,11 +20,13 @@ def test_sigpy_import():

@pytest.mark.sigpy
def test_slr():
import sigpy.mri.rf as sigpy_rf
from pypulseq.make_sigpy_pulse import sigpy_n_seq

slice_thickness = 3e-3
flip_angle = np.pi / 2
duration = 3e-3
time_bw_product = 4

system = Opts(
max_grad=32,
Expand All @@ -35,6 +36,7 @@ def test_slr():
rf_ringdown_time=30e-6,
rf_dead_time=100e-6,
)

pulse_cfg = SigpyPulseOpts(
pulse_type='slr',
ptype='st',
Expand All @@ -46,47 +48,49 @@ def test_slr():
band_sep=20,
phs_0_pt='None',
)

rfp, _, _ = sigpy_n_seq( # type: ignore
flip_angle=flip_angle,
system=system,
delay=100e-6,
duration=duration,
slice_thickness=slice_thickness,
time_bw_product=4,
return_gz=True,
slice_thickness=slice_thickness,
system=system,
time_bw_product=time_bw_product,
pulse_cfg=pulse_cfg,
plot=False,
delay=system.rf_dead_time,
)

# Check that the number of samples in the pulse is correct
assert rfp.signal.shape[0] == pytest.approx((duration + system.rf_ringdown_time) / system.rf_raster_time)

# Check that the pulse can be added to a PyPulseq Sequence
seq = pp.Sequence(system=system)
seq.add_block(rfp)

assert rfp.signal.shape[0] == pytest.approx((duration + system.rf_ringdown_time) / system.rf_raster_time)

# [a, b] = rf.sim.abrm(
# rfp.signal,
# np.arange(-20 * time_bw_product, 20 * time_bw_product, 40 * time_bw_product / 2000),
# True,
# )
# mag_xy = 2 * np.multiply(np.conj(a), b)
# # pl.LinePlot(Mxy)
# # print(np.sum(np.abs(Mxy)))
# # peaks, dict = sis.find_peaks(np.abs(Mxy),threshold=0.5, plateau_size=40)
# plateau_widths = np.sum(np.abs(mag_xy) > 0.8)
# assert plateau_widths == 29
# Check that the pulse shape is correct
pulse = (rfp.signal * system.rf_raster_time * 2 * np.pi) / flip_angle
[a, b] = sigpy_rf.sim.abrm(
pulse,
np.arange(-20 * time_bw_product, 20 * time_bw_product, 40 * time_bw_product / 2000),
True,
)
mag_xy = 2 * np.multiply(np.conj(a), b)
plateau_widths = np.sum(np.abs(mag_xy) > 0.8)
assert plateau_widths == 29


@pytest.mark.sigpy
def test_sms():
import sigpy.mri.rf as sigpy_rf
from pypulseq.make_sigpy_pulse import sigpy_n_seq

print('Testing SMS design')

slice_thickness = 3e-3 # Slice thickness
slice_thickness = 3e-3
flip_angle = np.pi / 2
duration = 3e-3
n_bands = 3
# Set system limits
time_bw_product = 4

system = Opts(
max_grad=32,
grad_unit='mT/m',
Expand All @@ -95,6 +99,7 @@ def test_sms():
rf_ringdown_time=30e-6,
rf_dead_time=100e-6,
)

pulse_cfg = SigpyPulseOpts(
pulse_type='sms',
ptype='st',
Expand All @@ -106,32 +111,34 @@ def test_sms():
band_sep=20,
phs_0_pt='None',
)

rfp, _, _ = sigpy_n_seq( # type: ignore
flip_angle=flip_angle,
system=system,
delay=100e-6,
duration=duration,
slice_thickness=slice_thickness,
time_bw_product=4,
return_gz=True,
slice_thickness=slice_thickness,
system=system,
time_bw_product=time_bw_product,
pulse_cfg=pulse_cfg,
plot=False,
delay=system.rf_dead_time,
)

# Check that the number of samples in the pulse is correct
assert rfp.signal.shape[0] == pytest.approx((duration + system.rf_ringdown_time) / system.rf_raster_time)

# Check that the pulse can be added to a PyPulseq Sequence
seq = pp.Sequence(system=system)
seq.add_block(rfp)

assert rfp.signal.shape[0] == pytest.approx((duration + system.rf_ringdown_time) / system.rf_raster_time)

# [a, b] = rf.sim.abrm(
# rfp.signal,
# np.arange(-20 * time_bw_product, 20 * time_bw_product, 40 * time_bw_product / 2000),
# True,
# )
# mag_xy = 2 * np.multiply(np.conj(a), b)
# # pl.LinePlot(Mxy)
# # print(np.sum(np.abs(Mxy)))
# # peaks, dict = sis.find_peaks(np.abs(Mxy),threshold=0.5, plateau_size=40)
# plateau_widths = np.sum(np.abs(mag_xy) > 0.8)
# # if slr has 29 > 0.8, then sms with MB = n_bands
# assert (29 * n_bands) == plateau_widths
# Check that the pulse shape is correct
pulse = (rfp.signal * system.rf_raster_time * 2 * np.pi) / flip_angle
[a, b] = sigpy_rf.sim.abrm(
pulse,
np.arange(-20 * time_bw_product, 20 * time_bw_product, 40 * time_bw_product / 2000),
True,
)
mag_xy = 2 * np.multiply(np.conj(a), b)
plateau_widths = np.sum(np.abs(mag_xy) > 0.8)
# if slr has 29 > 0.8, then sms with MB = n_bands
assert (29 * n_bands) == plateau_widths
Loading