Homec4science
Diffusion pypulseq (feature-spiral) Browse

pypulseq/gre_pypulseq_sigpy_demo2.py
1fa25251e5cdfeature-spiral

Last Change
Raw File

gre_pypulseq_sigpy_demo2.py
View Options

import math
import numpy as np
import pypulseq as pp
import sigpy.mri.rf as rf_ext
import sigpy2pulseq as sp
# ======
# SETUP
# ======
# Create a new sequence object and library options
seq = pp.Sequence()
# Define FOV and resolution
fov = 256e-3
Nx = 256
Ny = 256
alpha = 90 # flip angle
slice_thickness = 3e-3 # slice
TE = np.array([4.3e-3])
TR = 300e-3
rf_spoiling_inc = 117 # RF spoiling increment
ext_pulse_library = 1 # 0 - pypulseq, 1 - sigpy
ext_pulse_type = 'slr'
disp_pulse = 0
system = pp.Opts(max_grad=28, grad_unit='mT/m', max_slew=150, slew_unit='T/m/s', rf_ringdown_time=20e-6,
rf_dead_time=100e-6, adc_dead_time=10e-6, rf_raster_time=1e-6)
# Common pulse design parameters
duration = 3e-3 # seconds
tb = 4
# ======
# CREATE EVENTS
# ======
if (ext_pulse_library):
# Declare pulse configuration inputs here
pulse_freq_offset = 0.0 # if designing rf pulse without a FM
if ext_pulse_type == 'slr':
pulse = rf_ext.slr.dzrf(n=int(round(duration / system.rf_raster_time)), tb=tb, ptype='st', ftype='ls', d1=0.01,
d2=0.01,
cancel_alpha_phs=True)
if ext_pulse_type == 'sms':
pulse_in = rf_ext.slr.dzrf(n=int(round(duration / system.rf_raster_time)), tb=tb, ptype='st', ftype='ls',
d1=0.01, d2=0.01,
cancel_alpha_phs=False)
pulse = rf_ext.multiband.mb_rf(pulse_in, n_bands=3, band_sep=25, phs_0_pt='None')
if ext_pulse_type == 'hypsec':
# Let us design two pulses here
# - STA for excitation
pulse = rf_ext.slr.dzrf(n=int(round(duration / system.rf_raster_time)), tb=tb, ptype='st', ftype='ls', d1=0.01,
d2=0.01,
cancel_alpha_phs=True)
# - Hyperbolic secant pulse for inversion as part of magnetization prep
duration_hypsec = 4 * duration # adiabatic pulses need more time
pulse_hypsec, pulse_freq_offset = rf_ext.adiabatic.hypsec(n=2002,
beta=500,
mu=60, dur=duration_hypsec)
rf_inv = sp.sig_2_seq(pulse=pulse_hypsec, flip_angle=alpha * math.pi / 180, system=system,
duration=duration_hypsec,
slice_thickness=slice_thickness,
return_gz=False, time_bw_product=tb, rf_freq=pulse_freq_offset)
# Display if needed
if (disp_pulse):
if (ext_pulse_type == 'hypsec'):
sp.disp_pulse(pulse_hypsec, pulse_freq_offset, tb, duration_hypsec, system)
else:
sp.disp_pulse(pulse, 0, tb, duration, system)
# Convert and integrate with pypulseq
rf, gz, gzr, _ = sp.sig_2_seq(pulse=pulse, flip_angle=alpha * math.pi / 180, system=system, duration=duration,
slice_thickness=slice_thickness,
return_gz=True, time_bw_product=tb, rf_freq=0)
else:
rf, gz, gzr = pp.make_sinc_pulse(flip_angle=alpha * math.pi / 180, duration=duration,
slice_thickness=slice_thickness,
apodization=0.5, time_bw_product=tb, system=system, return_gz=True)
# Define other gradients and ADC events
delta_k = 1 / fov
gx = pp.make_trapezoid(channel='x', flat_area=Nx * delta_k, flat_time=3.2e-3, system=system)
adc = pp.make_adc(num_samples=Nx, duration=gx.flat_time, delay=gx.rise_time, system=system)
gx_pre = pp.make_trapezoid(channel='x', area=-gx.area / 2, duration=1e-3, system=system)
gz_reph = pp.make_trapezoid(channel='z', area=-gz.area / 2, duration=1e-3, system=system)
phase_areas = (np.arange(Ny) - Ny / 2) * delta_k
# gradient spoiling
gx_spoil = pp.make_trapezoid(channel='x', area=2 * Nx * delta_k, system=system)
gz_spoil = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system)
# Calculate timing
delay_TE = np.ceil((TE - pp.calc_duration(gx_pre) - gz.fall_time - gz.flat_time / 2 - pp.calc_duration(
gx) / 2) / seq.grad_raster_time) * seq.grad_raster_time
delay_TR = np.ceil((TR - pp.calc_duration(gz) - pp.calc_duration(gx_pre) - pp.calc_duration(
gx) - delay_TE) / seq.grad_raster_time) * seq.grad_raster_time
assert np.all(delay_TE >= 0)
assert np.all(delay_TR >= pp.calc_duration(gx_spoil, gz_spoil))
rf_phase = 0
rf_inc = 0
# demo slice profile
gz.channel = 'x'
# ======
# CONSTRUCT SEQUENCE
# ======
# Loop over phase encodes and define sequence blocks
# Insert adiabatic inversion here if chosen - 12 ms - for demo purposes only
for i in range(Ny):
for j in range(len(TE)):
rf.phase_offset = rf_phase / 180 * np.pi
adc.phase_offset = rf_phase / 180 * np.pi
rf_inc = divmod(rf_inc + rf_spoiling_inc, 360.0)[1]
rf_phase = divmod(rf_phase + rf_inc, 360.0)[1]
if ext_pulse_type == 'hypsec':
if np.mod(i, 64) == 0:
print('Mag prep now')
for n in range(len(rf_inv)):
seq.add_block(rf_inv[n])
seq.add_block(pp.make_delay(300e-3))
# seq.plot()
seq.add_block(rf, gz)
gy_pre = pp.make_trapezoid(channel='y', area=phase_areas[i], duration=pp.calc_duration(gx_pre), system=system)
seq.add_block(gx_pre, gy_pre, gz_reph)
seq.add_block(pp.make_delay(delay_TE[j]))
seq.add_block(gx, adc)
gy_pre.amplitude = -gy_pre.amplitude
seq.add_block(pp.make_delay(delay_TR[j]), gx_spoil, gy_pre, gz_spoil)
ok, error_report = seq.check_timing() # Check whether the timing of the sequence is correct
if ok:
print('Timing check passed successfully')
else:
print('Timing check failed. Error listing follows:')
[print(e) for e in error_report]
# ======
# VISUALIZATION
# ======
# seq.plot()
# Trajectory calculation and plotting
# ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc = seq.calculate_kspace()
# time_axis = np.arange(1, ktraj.shape[1] + 1) * system.grad_raster_time
# plt.figure()
# plt.plot(time_axis, ktraj.T) # Plot the entire k-space trajectory
# plt.plot(t_adc, ktraj_adc[0], '.') # Plot sampling points on the kx-axis
# plt.figure()
# plt.plot(ktraj[0], ktraj[1], 'b') # 2D plot
# plt.axis('equal') # Enforce aspect ratio for the correct trajectory display
# plt.plot(ktraj_adc[0], ktraj_adc[1], 'r.') # Plot sampling points
# plt.show()
# Prepare the sequence output for the scanner
seq.set_definition('FOV', [fov, fov, slice_thickness])
seq.set_definition('Name', 'gre')
seq.write('gre_pypulseq_ext_profile_' + str(alpha) + '_' + str(ext_pulse_type) + '.seq')
# Very optional slow step, but useful for testing during development e.g. for the real TE, TR or for staying within
# slew-rate limits
# rep = seq.test_report()
# print(rep)
c4science · Help