diff --git a/pypulseq/seq_examples/scripts/write_haste.py b/pypulseq/seq_examples/scripts/write_haste.py index 5365479..2cb6898 100644 --- a/pypulseq/seq_examples/scripts/write_haste.py +++ b/pypulseq/seq_examples/scripts/write_haste.py @@ -1,205 +1,212 @@ import math import warnings import matplotlib.pyplot as plt import numpy as np from pypulseq.Sequence.sequence import Sequence from pypulseq.calc_rf_center import calc_rf_center from pypulseq.make_adc import make_adc from pypulseq.make_delay import make_delay from pypulseq.make_extended_trapezoid import make_extended_trapezoid from pypulseq.make_sinc_pulse import make_sinc_pulse from pypulseq.make_trap_pulse import make_trapezoid from pypulseq.opts import Opts # ====== # SETUP # ====== dG = 250e-6 # Set system limits system = Opts(max_grad=30, grad_unit='mT/m', max_slew=170, slew_unit='T/m/s', rf_ringdown_time=100e-6, rf_dead_time=100e-6, adc_dead_time=10e-6) seq = Sequence(system=system) # Create a new sequence object # Define FOV and resolution fov = 256e-3 Ny_pre = 8 Nx, Ny = 128, 128 n_echo = int(Ny / 2 + Ny_pre) n_slices = 1 rf_flip = 180 if isinstance(rf_flip, int): rf_flip = np.zeros(n_echo) + rf_flip slice_thickness = 5e-3 # Slice thickness TE = 12e-3 TR = 2000e-3 TE_eff = 60e-3 k0 = round(TE_eff / TE) PE_type = 'linear' readout_time = 6.4e-3 + 2 * system.adc_dead_time t_ex = 2.5e-3 t_ex_wd = t_ex + system.rf_ringdown_time + system.rf_dead_time t_ref = 2e-3 tf_ref_wd = t_ref + system.rf_ringdown_time + system.rf_dead_time t_sp = 0.5 * (TE - readout_time - tf_ref_wd) t_sp_ex = 0.5 * (TE - t_ex_wd - tf_ref_wd) fspR = 1.0 fspS = 0.5 rfex_phase = math.pi / 2 rfref_phase = 0 # ====== # CREATE EVENTS # ====== # Create 90 degree slice selection pulse and gradient flipex = 90 * math.pi / 180 rfex, gz, _ = make_sinc_pulse(flip_angle=flipex, system=system, duration=t_ex, slice_thickness=slice_thickness, apodization=0.5, time_bw_product=4, phase_offset=rfex_phase, return_gz=True) GS_ex = make_trapezoid(channel='z', system=system, amplitude=gz.amplitude, flat_time=t_ex_wd, rise_time=dG) flipref = rf_flip[0] * math.pi / 180 rfref, gz, _ = make_sinc_pulse(flip_angle=flipref, system=system, duration=t_ref, slice_thickness=slice_thickness, apodization=0.5, time_bw_product=4, phase_offset=rfref_phase, use='refocusing', return_gz=True) GS_ref = make_trapezoid(channel='z', system=system, amplitude=GS_ex.amplitude, flat_time=tf_ref_wd, rise_time=dG) AGS_ex = GS_ex.area / 2 GS_spr = make_trapezoid(channel='z', system=system, area=AGS_ex * (1 + fspS), duration=t_sp, rise_time=dG) GS_spex = make_trapezoid(channel='z', system=system, area=AGS_ex * fspS, duration=t_sp_ex, rise_time=dG) delta_k = 1 / fov k_width = Nx * delta_k GR_acq = make_trapezoid(channel='x', system=system, flat_area=k_width, flat_time=readout_time, rise_time=dG) -adc = make_adc(num_samples=Nx, duration=GR_acq.flat_time - 40e-6, delay=20e-6) +adc = make_adc(num_samples=Nx, duration=GR_acq.flat_time - 2 * system.adc_dead_time, delay=20e-6) GR_spr = make_trapezoid(channel='x', system=system, area=GR_acq.area * fspR, duration=t_sp, rise_time=dG) GR_spex = make_trapezoid(channel='x', system=system, area=GR_acq.area * (1 + fspR), duration=t_sp_ex, rise_time=dG) AGR_spr = GR_spr.area AGR_preph = GR_acq.area / 2 + AGR_spr GR_preph = make_trapezoid(channel='x', system=system, area=AGR_preph, duration=t_sp_ex, rise_time=dG) n_ex = 1 PE_order = np.arange(-Ny_pre, Ny + 1).T phase_areas = PE_order * delta_k # Split gradients and recombine into blocks GS1_times = [0, GS_ex.rise_time] GS1_amp = [0, GS_ex.amplitude] GS1 = make_extended_trapezoid(channel='z', times=GS1_times, amplitudes=GS1_amp) GS2_times = [0, GS_ex.flat_time] GS2_amp = [GS_ex.amplitude, GS_ex.amplitude] GS2 = make_extended_trapezoid(channel='z', times=GS2_times, amplitudes=GS2_amp) GS3_times = [0, GS_spex.rise_time, GS_spex.rise_time + GS_spex.flat_time, GS_spex.rise_time + GS_spex.flat_time + GS_spex.fall_time] GS3_amp = [GS_ex.amplitude, GS_spex.amplitude, GS_spex.amplitude, GS_ref.amplitude] GS3 = make_extended_trapezoid(channel='z', times=GS3_times, amplitudes=GS3_amp) GS4_times = [0, GS_ref.flat_time] GS4_amp = [GS_ref.amplitude, GS_ref.amplitude] GS4 = make_extended_trapezoid(channel='z', times=GS4_times, amplitudes=GS4_amp) GS5_times = [0, GS_spr.rise_time, GS_spr.rise_time + GS_spr.flat_time, GS_spr.rise_time + GS_spr.flat_time + GS_spr.fall_time] GS5_amp = [GS_ref.amplitude, GS_spr.amplitude, GS_spr.amplitude, 0] GS5 = make_extended_trapezoid(channel='z', times=GS5_times, amplitudes=GS5_amp) GS7_times = [0, GS_spr.rise_time, GS_spr.rise_time + GS_spr.flat_time, GS_spr.rise_time + GS_spr.flat_time + GS_spr.fall_time] GS7_amp = [0, GS_spr.amplitude, GS_spr.amplitude, GS_ref.amplitude] GS7 = make_extended_trapezoid(channel='z', times=GS7_times, amplitudes=GS7_amp) # Readout gradient GR3 = GR_preph GR5_times = [0, GR_spr.rise_time, GR_spr.rise_time + GR_spr.flat_time, GR_spr.rise_time + GR_spr.flat_time + GR_spr.fall_time] GR5_amp = [0, GR_spr.amplitude, GR_spr.amplitude, GR_acq.amplitude] GR5 = make_extended_trapezoid(channel='x', times=GR5_times, amplitudes=GR5_amp) GR6_times = [0, readout_time] GR6_amp = [GR_acq.amplitude, GR_acq.amplitude] GR6 = make_extended_trapezoid(channel='x', times=GR6_times, amplitudes=GR6_amp) GR7_times = [0, GR_spr.rise_time, GR_spr.rise_time + GR_spr.flat_time, GR_spr.rise_time + GR_spr.flat_time + GR_spr.fall_time] GR7_amp = [GR_acq.amplitude, GR_spr.amplitude, GR_spr.amplitude, 0] GR7 = make_extended_trapezoid(channel='x', times=GR7_times, amplitudes=GR7_amp) # Fill-times tex = GS1.t[-1] + GS2.t[-1] + GS3.t[-1] tref = GS4.t[-1] + GS5.t[-1] + GS7.t[-1] + readout_time tend = GS4.t[-1] + GS5.t[-1] TE_train = tex + n_echo * tref + tend TR_fill = (TR - n_slices * TE_train) / n_slices # Round to gradient raster TR_fill = system.grad_raster_time * round(TR_fill / system.grad_raster_time) if TR_fill < 0: TR_fill = 1e-3 warnings.warn(f'TR too short, adapted to include all slices to: {1000 * n_slices * (TE_train + TR_fill)} ms') else: print(f'TR fill: {1000 * TR_fill} ms') delay_TR = make_delay(TR_fill) delay_end = make_delay(5) # ====== # CONSTRUCT SEQUENCE # ====== # Define sequence blocks for k_ex in range(n_ex): for s in range(n_slices): rfex.freq_offset = GS_ex.amplitude * slice_thickness * (s - (n_slices - 1) / 2) rfref.freq_offset = GS_ref.amplitude * slice_thickness * (s - (n_slices - 1) / 2) # Align the phase for off-center slices rfex.phase_offset = rfex_phase - 2 * math.pi * rfex.freq_offset * calc_rf_center(rfex)[0] rfref.phase_offset = rfref_phase - 2 * math.pi * rfref.freq_offset * calc_rf_center(rfref)[0] seq.add_block(GS1) seq.add_block(GS2, rfex) seq.add_block(GS3, GR3) for k_ech in range(n_echo): if k_ex >= 0: phase_area = phase_areas[k_ech] else: phase_area = 0 GP_pre = make_trapezoid(channel='y', system=system, area=phase_area, duration=t_sp, rise_time=dG) GP_rew = make_trapezoid(channel='y', system=system, area=-phase_area, duration=t_sp, rise_time=dG) seq.add_block(GS4, rfref) seq.add_block(GS5, GR5, GP_pre) if k_ex >= 0: seq.add_block(GR6, adc) else: seq.add_block(GR6) seq.add_block(GS7, GR7, GP_rew) seq.add_block(GS4) seq.add_block(GS5) seq.add_block(delay_TR) seq.add_block(delay_end) +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() ktraj_adc, ktraj, t_excitation, t_refocusing, _ = seq.calculate_kspace() plt.plot(ktraj.T) # Plot the entire k-space trajectory 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.') plt.show() seq.write('haste_pypulseq.seq') diff --git a/pypulseq/seq_examples/scripts/write_tse.py b/pypulseq/seq_examples/scripts/write_tse.py index f94b8f5..b674e79 100644 --- a/pypulseq/seq_examples/scripts/write_tse.py +++ b/pypulseq/seq_examples/scripts/write_tse.py @@ -1,193 +1,200 @@ import math import warnings import numpy as np from matplotlib import pyplot as plt import pypulseq as pp # ====== # SETUP # ====== dG = 250e-6 # Set system limits system = pp.Opts(max_grad=32, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s', rf_ringdown_time=100e-6, rf_dead_time=100e-6, adc_dead_time=10e-6) seq = pp.Sequence(system) # Create a new sequence object # Define FOV and resolution fov = 256e-3 Nx, Ny = 128, 128 n_echo = 16 n_slices = 1 rf_flip = 180 if isinstance(rf_flip, int): rf_flip = np.zeros(n_echo) + rf_flip slice_thickness = 5e-3 TE = 12e-3 TR = 2000e-3 TE_eff = 60e-3 k0 = round(TE_eff / TE) pe_type = 'linear' readout_time = 6.4e-3 + 2 * system.adc_dead_time t_ex = 2.5e-3 t_exwd = t_ex + system.rf_ringdown_time + system.rf_dead_time t_ref = 2e-3 t_refwd = t_ref + system.rf_ringdown_time + system.rf_dead_time t_sp = 0.5 * (TE - readout_time - t_refwd) t_spex = 0.5 * (TE - t_exwd - t_refwd) fsp_r = 1 fsp_s = 0.5 rf_ex_phase = np.pi / 2 rf_ref_phase = 0 # ====== # CREATE EVENTS # ====== flip_ex = 90 * np.pi / 180 rf_ex, gz, _ = pp.make_sinc_pulse(flip_angle=flip_ex, system=system, duration=t_ex, slice_thickness=slice_thickness, apodization=0.5, time_bw_product=4, phase_offset=rf_ex_phase, return_gz=True) gs_ex = pp.make_trapezoid(channel='z', system=system, amplitude=gz.amplitude, flat_time=t_exwd, rise_time=dG) flip_ref = rf_flip[0] * np.pi / 180 rf_ref, gz, _ = pp.make_sinc_pulse(flip_angle=flip_ref, system=system, duration=t_ref, slice_thickness=slice_thickness, apodization=0.5, time_bw_product=4, phase_offset=rf_ref_phase, use='refocusing', return_gz=True) gs_ref = pp.make_trapezoid(channel='z', system=system, amplitude=gs_ex.amplitude, flat_time=t_refwd, rise_time=dG) ags_ex = gs_ex.area / 2 gs_spr = pp.make_trapezoid(channel='z', system=system, area=ags_ex * (1 + fsp_s), duration=t_sp, rise_time=dG) gs_spex = pp.make_trapezoid(channel='z', system=system, area=ags_ex * fsp_s, duration=t_spex, rise_time=dG) delta_k = 1 / fov k_width = Nx * delta_k gr_acq = pp.make_trapezoid(channel='x', system=system, flat_area=k_width, flat_time=readout_time, rise_time=dG) -adc = pp.make_adc(num_samples=Nx, duration=gr_acq.flat_time - 40e-6, delay=20e-6) +adc = pp.make_adc(num_samples=Nx, duration=gr_acq.flat_time - 2 * system.adc_dead_time, delay=20e-6) gr_spr = pp.make_trapezoid(channel='x', system=system, area=gr_acq.area * fsp_r, duration=t_sp, rise_time=dG) gr_spex = pp.make_trapezoid(channel='x', system=system, area=gr_acq.area * (1 + fsp_r), duration=t_spex, rise_time=dG) agr_spr = gr_spr.area agr_preph = gr_acq.area / 2 + agr_spr gr_preph = pp.make_trapezoid(channel='x', system=system, area=agr_preph, duration=t_spex, rise_time=dG) # Phase-encoding n_ex = math.floor(Ny / n_echo) pe_steps = np.arange(1, n_echo * n_ex + 1) - 0.5 * n_echo * n_ex - 1 if divmod(n_echo, 2)[1] == 0: pe_steps = np.roll(pe_steps, -round(n_ex / 2)) pe_order = pe_steps.reshape((n_ex, n_echo), order='F').T phase_areas = pe_order * delta_k # Split gradients and recombine into blocks gs1_times = [0, gs_ex.rise_time] gs1_amp = [0, gs_ex.amplitude] gs1 = pp.make_extended_trapezoid(channel='z', times=gs1_times, amplitudes=gs1_amp) gs2_times = [0, gs_ex.flat_time] gs2_amp = [gs_ex.amplitude, gs_ex.amplitude] gs2 = pp.make_extended_trapezoid(channel='z', times=gs2_times, amplitudes=gs2_amp) gs3_times = [0, gs_spex.rise_time, gs_spex.rise_time + gs_spex.flat_time, gs_spex.rise_time + gs_spex.flat_time + gs_spex.fall_time] gs3_amp = [gs_ex.amplitude, gs_spex.amplitude, gs_spex.amplitude, gs_ref.amplitude] gs3 = pp.make_extended_trapezoid(channel='z', times=gs3_times, amplitudes=gs3_amp) gs4_times = [0, gs_ref.flat_time] gs4_amp = [gs_ref.amplitude, gs_ref.amplitude] gs4 = pp.make_extended_trapezoid(channel='z', times=gs4_times, amplitudes=gs4_amp) gs5_times = [0, gs_spr.rise_time, gs_spr.rise_time + gs_spr.flat_time, gs_spr.rise_time + gs_spr.flat_time + gs_spr.fall_time] gs5_amp = [gs_ref.amplitude, gs_spr.amplitude, gs_spr.amplitude, 0] gs5 = pp.make_extended_trapezoid(channel='z', times=gs5_times, amplitudes=gs5_amp) gs7_times = [0, gs_spr.rise_time, gs_spr.rise_time + gs_spr.flat_time, gs_spr.rise_time + gs_spr.flat_time + gs_spr.fall_time] gs7_amp = [0, gs_spr.amplitude, gs_spr.amplitude, gs_ref.amplitude] gs7 = pp.make_extended_trapezoid(channel='z', times=gs7_times, amplitudes=gs7_amp) # Readout gradient gr3 = gr_preph gr5_times = [0, gr_spr.rise_time, gr_spr.rise_time + gr_spr.flat_time, gr_spr.rise_time + gr_spr.flat_time + gr_spr.fall_time] gr5_amp = [0, gr_spr.amplitude, gr_spr.amplitude, gr_acq.amplitude] gr5 = pp.make_extended_trapezoid(channel='x', times=gr5_times, amplitudes=gr5_amp) gr6_times = [0, readout_time] gr6_amp = [gr_acq.amplitude, gr_acq.amplitude] gr6 = pp.make_extended_trapezoid(channel='x', times=gr6_times, amplitudes=gr6_amp) gr7_times = [0, gr_spr.rise_time, gr_spr.rise_time + gr_spr.flat_time, gr_spr.rise_time + gr_spr.flat_time + gr_spr.fall_time] gr7_amp = [gr_acq.amplitude, gr_spr.amplitude, gr_spr.amplitude, 0] gr7 = pp.make_extended_trapezoid(channel='x', times=gr7_times, amplitudes=gr7_amp) # Fill-times t_ex = gs1.t[-1] + gs2.t[-1] + gs3.t[-1] t_ref = gs4.t[-1] + gs5.t[-1] + gs7.t[-1] + readout_time t_end = gs4.t[-1] + gs5.t[-1] TE_train = t_ex + n_echo * t_ref + t_end TR_fill = (TR - n_slices * TE_train) / n_slices TR_fill = system.grad_raster_time * round(TR_fill / system.grad_raster_time) # Round to gradient raster if TR_fill < 0: TR_fill = 1e-3 warnings.warn(f'TR too short, adapted to include all slices to: {1000 * n_slices * (TE_train + TR_fill)} ms') else: print(f'TR fill: {1000 * TR_fill} ms') delay_TR = pp.make_delay(TR_fill) # ====== # CONSTRUCT SEQUENCE # ====== for k_ex in range(n_ex + 1): for s in range(n_slices): rf_ex.freq_offset = gs_ex.amplitude * slice_thickness * (s - (n_slices - 1) / 2) rf_ref.freq_offset = gs_ref.amplitude * slice_thickness * (s - (n_slices - 1) / 2) rf_ex.phase_offset = rf_ex_phase - 2 * np.pi * rf_ex.freq_offset * pp.calc_rf_center(rf_ex)[0] rf_ref.phase_offset = rf_ref_phase - 2 * np.pi * rf_ref.freq_offset * pp.calc_rf_center(rf_ref)[0] seq.add_block(gs1) seq.add_block(gs2, rf_ex) seq.add_block(gs3, gr3) for k_echo in range(n_echo): if k_ex > 0: phase_area = phase_areas[k_echo, k_ex - 1] else: phase_area = 0.0 # 0.0 and not 0 because -phase_area should successfully result in negative zero gp_pre = pp.make_trapezoid(channel='y', system=system, area=phase_area, duration=t_sp, rise_time=dG) gp_rew = pp.make_trapezoid(channel='y', system=system, area=-phase_area, duration=t_sp, rise_time=dG) seq.add_block(gs4, rf_ref) seq.add_block(gs5, gr5, gp_pre) if k_ex > 0: seq.add_block(gr6, adc) else: seq.add_block(gr6) seq.add_block(gs7, gr7, gp_rew) seq.add_block(gs4) seq.add_block(gs5) seq.add_block(delay_TR) +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() # Plot k-spaces ktraj_adc, ktraj, t_excitation, t_refocusing, _ = seq.calculate_kspace() plt.plot(ktraj.T) # Plot the entire k-space trajectory 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.') plt.show() seq.write('tse_pypulseq.seq') diff --git a/pypulseq/seq_examples/scripts/write_ute.py b/pypulseq/seq_examples/scripts/write_ute.py index f9a406b..fdbf06f 100644 --- a/pypulseq/seq_examples/scripts/write_ute.py +++ b/pypulseq/seq_examples/scripts/write_ute.py @@ -1,135 +1,141 @@ """ A very basic UTE-like sequence, without ramp-sampling, ramp-RF. Achieves TE in the range of 300-400 us """ from copy import copy import numpy as np from matplotlib import pyplot as plt import pypulseq as pp # ====== # SETUP # ====== seq = pp.Sequence() # Create a new sequence object # Define FOV and resolution fov = 250e-3 -Nx = 256 +Nx = 250 alpha = 10 # Flip angle slice_thickness = 3e-3 # Slice thickness TR = 10e-3 # TR Nr = 128 # Number of radial spokes delta = 2 * np.pi / Nr # AAngular increment; try golden angle pi*(3-5^0.5) or 0.5 of i -ro_duration = 2.4e-3 # Read-out time: controls RO bandwidth and T2-blurring +ro_duration = 2.5e-3 # Read-out time: controls RO bandwidth and T2-blurring ro_os = 2 # Oversampling ro_asymmetry = 0.97 # 0: Fully symmetric; 1: half-echo minRF_to_ADC_time = 50e-6 # Defines TE together with the RO asymmetry rf_spoiling_inc = 117 # RF spoiling increment # Set system limits system = pp.Opts(max_grad=28, grad_unit='mT/m', max_slew=100, slew_unit='T/m/s', rf_ringdown_time=20e-6, rf_dead_time=100e-6, adc_dead_time=10e-6) # ====== # CREATE EVENTS # ====== # Create alpha-degree slice selection pulse and gradient rf, gz, gz_reph = pp.make_sinc_pulse(flip_angle=alpha * np.pi / 180, duration=1e-3, slice_thickness=slice_thickness, apodization=0.5, time_bw_product=2, center_pos=1, system=system, return_gz=True) # Align RO asymmetry to ADC samples Nxo = np.round(ro_os * Nx) ro_asymmetry = np.round(ro_asymmetry * Nxo / 2) / Nxo * 2 # Define other gradients and ADC events delta_k = 1 / fov / (1 + ro_asymmetry) ro_area = Nx * delta_k gx = pp.make_trapezoid(channel='x', flat_area=ro_area, flat_time=ro_duration, system=system) adc = pp.make_adc(num_samples=Nxo, duration=gx.flat_time, delay=gx.rise_time, system=system) -adc.delay = adc.delay - 0.5 * adc.dwell # compensate for the 0.5 samples shift gx_pre = pp.make_trapezoid(channel='x', area=-(gx.area - ro_area) / 2 - ro_area / 2 * (1 - ro_asymmetry), system=system) # Gradient spoiling gx_spoil = pp.make_trapezoid(channel='x', area=0.2 * Nx * delta_k, system=system) # Calculate timing TE = gz.fall_time + pp.calc_duration(gx_pre, gz_reph) + gx.rise_time + adc.dwell * Nxo / 2 * (1 - ro_asymmetry) delay_TR = np.ceil((TR - pp.calc_duration(gx_pre, gz_reph) - pp.calc_duration(gz) - pp.calc_duration( gx)) / seq.grad_raster_time) * seq.grad_raster_time assert np.all(delay_TR >= pp.calc_duration(gx_spoil)) print(f'TE = {TE * 1e6:.0f} us') if pp.calc_duration(gz_reph) > pp.calc_duration(gx_pre): gx_pre.delay = pp.calc_duration(gz_reph) - pp.calc_duration(gx_pre) rf_phase = 0 rf_inc = 0 # ====== # CONSTRUCT SEQUENCE # ====== for i in range(Nr): for c in range(2): rf.phase_offset = rf_phase / 180 * np.pi adc.phase_offset = rf_phase / 180 * np.pi rf_inc = np.mod(rf_inc + rf_spoiling_inc, 360.0) rf_phase = np.mod(rf_phase + rf_inc, 360.0) gz.amplitude = -gz.amplitude # Alternate GZ amplitude gz_reph.amplitude = -gz_reph.amplitude seq.add_block(rf, gz) phi = delta * i gpc = copy(gx_pre) gps = copy(gx_pre) gpc.amplitude = gx_pre.amplitude * np.cos(phi) gps.amplitude = gx_pre.amplitude * np.sin(phi) gps.channel = 'y' grc = copy(gx) grs = copy(gx) grc.amplitude = gx.amplitude * np.cos(phi) grs.amplitude = gx.amplitude * np.sin(phi) grs.channel = 'y' gsc = copy(gx_spoil) gss = copy(gx_spoil) gsc.amplitude = gx_spoil.amplitude * np.cos(phi) gss.amplitude = gx_spoil.amplitude * np.sin(phi) gss.channel = 'y' seq.add_block(gpc, gps, gz_reph) seq.add_block(grc, grs, adc) seq.add_block(gsc, gss, pp.make_delay(delay_TR)) +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() # Plot gradients to check for gaps and optimality of the timing gw = seq.gradient_waveforms() plt.plot(gw.T) # Plot the entire gradient shape # Trajectory calculation ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc = seq.calculate_kspace() # Plot k-spaces time_axis = np.arange(ktraj.shape[1]) * seq.grad_raster_time plt.figure() plt.plot(time_axis, ktraj.T) # Plot the entire k-space trajectory plt.plot(t_adc, ktraj_adc[0], '.') # 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 the sampling points plt.show() seq.set_definition('FOV', [fov, fov, slice_thickness]) seq.set_definition('Name', 'UTE') seq.write('ute_pypulseq.seq')