diff --git a/examples/t1_molli_bssfp.ipynb b/examples/t1_molli_bssfp.ipynb index b4631fd..ba79d1f 100644 --- a/examples/t1_molli_bssfp.ipynb +++ b/examples/t1_molli_bssfp.ipynb @@ -55,7 +55,7 @@ "metadata": {}, "source": [ "### Settings\n", - "We are going to use a numerical phantom with a matrix size of 128 x 128. The repetition time is set to 20 seconds to ensure we can accurately estimate long T1 times of the CSF. " + "We are going to use a numerical phantom with a matrix size of 64 x 64." ] }, { diff --git a/src/mrseq/scripts/b0_b1_wasabi.py b/src/mrseq/scripts/b0_b1_wasabi.py new file mode 100644 index 0000000..6a8fd0f --- /dev/null +++ b/src/mrseq/scripts/b0_b1_wasabi.py @@ -0,0 +1,364 @@ +"""Simultaneous B0 and B1 mapping using the WASABI method.""" + +from pathlib import Path + +import numpy as np +import pypulseq as pp + +from mrseq.utils import find_gx_flat_time_on_adc_raster +from mrseq.utils import round_to_raster +from mrseq.utils import sys_defaults +from mrseq.utils import write_sequence +from mrseq.utils.constants import GYROMAGNETIC_RATIO_PROTON +from mrseq.utils.trajectory import cartesian_phase_encoding + + +def wasabi_gre_centric_kernel( + system: pp.Opts, + frequency_offsets: np.ndarray, + norm_offset: float | None, + t_recovery: float, + t_recovery_norm: float, + fov_xy: float, + n_readout: int, + readout_oversampling: float, + n_phase_encoding: int, + slice_thickness: float, + rf_prep_duration: float, + rf_prep_amplitude: float, + rf_duration: float, + rf_flip_angle: float, + rf_bwt: float, + rf_apodization: float, + rf_spoiling_inc: float, + adc_dwell_time: float, +) -> pp.Sequence: + """Generate a WASABI sequence for simultaneous B0 and B1 mapping using a centric-out cartesian GRE readout. + + Parameters + ---------- + system + PyPulseq system limits object. + frequency_offsets + Array of frequency offsets (in Hz). + norm_offset + Frequency offset of normalization offset (in Hz). + t_recovery + Recovery time between frequency offsets (in seconds). + t_recovery_norm + Recovery time before normalization offset (in seconds). + fov_xy + Field of view in x and y direction (in meters). + n_readout + Number of frequency encoding steps. + readout_oversampling + Readout oversampling factor, commonly 2. This reduces aliasing artifacts. + n_phase_encoding + Number of phase encoding steps. + slice_thickness + Slice thickness of the 2D slice (in meters). + rf_prep_duration + Duration of WASABI block pulse (in seconds) + rf_prep_amplitude + Amplitude of WASABI block pulse (in µT) + rf_duration + Duration of the rf excitation pulse (in seconds) + rf_flip_angle + Flip angle of rf excitation pulse (in degrees) + rf_bwt + Bandwidth-time product of rf excitation pulse (Hz * seconds) + rf_apodization + Apodization factor of rf excitation pulse + rf_spoiling_inc + Phase increment used for RF spoiling. Set to 0 to disable RF spoiling. + adc_dwell_time + Dwell time of ADC. + + Returns + ------- + seq + PyPulseq Sequence object. + """ + # create PyPulseq Sequence object and set system limits + seq = pp.Sequence(system=system) + + if readout_oversampling < 1: + raise ValueError('Readout oversampling factor must be >= 1.') + + # add normalization offset to beginning of frequency offsets if specified + if norm_offset is not None: + frequency_offsets = np.concatenate(([norm_offset], frequency_offsets)) + + # create WASABI block pulse + rf_prep_flipangle_rad = rf_prep_amplitude * 1e-6 * rf_prep_duration * GYROMAGNETIC_RATIO_PROTON * 2 * np.pi + rf_prep = pp.make_block_pulse( + flip_angle=rf_prep_flipangle_rad, + duration=rf_prep_duration, + delay=system.rf_dead_time, + system=system, + use='preparation', + ) + + # create slice selective excitation pulse and gradients + rf, gz, gzr = pp.make_sinc_pulse( + flip_angle=np.deg2rad(rf_flip_angle), + duration=rf_duration, + slice_thickness=slice_thickness, + apodization=rf_apodization, + time_bw_product=rf_bwt, + delay=system.rf_dead_time, + system=system, + return_gz=True, + use='excitation', + ) + + # define centric-out phase encoding steps + kpe, _ = cartesian_phase_encoding( + n_phase_encoding=n_phase_encoding, + acceleration=1, + sampling_order='low_high', + ) + + # create readout gradient and ADC + delta_k = 1 / (fov_xy * readout_oversampling) + n_readout_with_oversampling = int(n_readout * readout_oversampling) + gx = pp.make_trapezoid( + channel='x', + flat_area=n_readout_with_oversampling * delta_k, + flat_time=n_readout_with_oversampling * adc_dwell_time, + system=system, + ) + n_readout_with_oversampling = n_readout_with_oversampling + np.mod(n_readout_with_oversampling, 2) # make even + adc = pp.make_adc(num_samples=n_readout_with_oversampling, duration=gx.flat_time, delay=gx.rise_time, system=system) + + # create frequency encoding pre- and re-winder gradient + gx_pre = pp.make_trapezoid(channel='x', area=-gx.area / 2 - delta_k / 2, system=system) + gx_post = pp.make_trapezoid(channel='x', area=-gx.area / 2 + delta_k / 2, system=system) + k0_center_id = np.where((np.arange(n_readout_with_oversampling) - n_readout_with_oversampling / 2) * delta_k == 0)[ + 0 + ][0] + + # create readout spoiler gradient + gz_spoil = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system) + + # create post preparation spoiler gradient + prep_spoil = pp.make_trapezoid(channel='z', area=5 * gz_spoil.area, system=system) + + # calculate minimum echo time + min_te = ( + rf.shape_dur / 2 # time from center to end of RF pulse + + max(rf.ringdown_time, gz.fall_time) # RF ringdown time or gradient fall time + + pp.calc_duration(gzr, gx_pre) # slice selection rewinder gradient + + gx.delay # potential delay of readout gradient + + gx.rise_time # rise time of readout gradient + + (k0_center_id + 0.5) * adc.dwell # time from beginning of ADC to time point of k-space center sample + ).item() + + # calculate minimum repetition time + min_tr = ( + pp.calc_duration(gz) # rf pulse + + pp.calc_duration(gzr, gx_pre) # slice selection re-phasing gradient and readout pre-winder + + pp.calc_duration(gx) # readout gradient + + pp.calc_duration(gz_spoil, gx_post) # gradient spoiler or readout-re-winder + ) + + # loop over frequency offsets + for rep_idx, freq_offset_hz in enumerate(frequency_offsets): + # set repetition ('REP') label for current frequency offset + rep_label = pp.make_label(type='SET', label='REP', value=int(rep_idx)) + + # add delay for normalization offset + if rep_idx == 0 and norm_offset is not None: + seq.add_block(pp.make_delay(t_recovery_norm)) + else: + seq.add_block(pp.make_delay(t_recovery)) + + # update frequency offset of WASABI block pulse and add it to sequence + rf_prep.freq_offset = freq_offset_hz + seq.add_block(rf_prep) + + # add post prep spoiler gradient + seq.add_block(prep_spoil) + + # reset rf spoiling + rf_phase = 0.0 + rf_inc = 0.0 + + # loop over phase encoding steps + for pe_idx in kpe: + # set phase encoding ('LIN') label + pe_label = pp.make_label(type='SET', label='LIN', value=int(pe_idx + np.abs(np.min(kpe)))) + + if rf_spoiling_inc > 0: + rf.phase_offset = np.deg2rad(rf_phase) + adc.phase_offset = np.deg2rad(rf_phase) + + # update rf pulse properties for the next loop + rf_inc = divmod(rf_inc + rf_spoiling_inc, 360.0)[1] + rf_phase = divmod(rf_phase + rf_inc, 360.0)[1] + + # calculate phase encoding gradient for current phase encoding step + gy_pre = pp.make_trapezoid( + channel='y', area=pe_idx * delta_k, duration=pp.calc_duration(gzr, gx_pre), system=system + ) + + # add slice-selective rf pulse + seq.add_block(rf, gz) + + # add slice-selection rewinder and readout pre-winder + seq.add_block(gzr, gx_pre, gy_pre) + + # add readout gradient, ADC and labels + seq.add_block(gx, adc, pe_label, rep_label) + + # add x and y re-winder and spoiler gradient in z-direction + gy_post = pp.make_trapezoid( + channel='y', amplitude=-gy_pre.amplitude, rise_time=gy_pre.rise_time, flat_time=gy_pre.flat_time + ) + seq.add_block(gx_post, gy_post, gz_spoil) + + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) + seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', min_te) + seq.set_definition('TR', min_tr) + seq.set_definition('frequency_offsets', frequency_offsets.tolist()) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) + + return seq + + +def main( + system: pp.Opts | None = None, + frequency_offsets: np.ndarray | None = None, + norm_offset: float | None = -35e3, + t_recovery: float = 3.0, + t_recovery_norm: float = 12.0, + fov_xy: float = 200e-3, + n_readout: int = 128, + n_phase_encoding: int = 128, + slice_thickness: float = 8e-3, + receiver_bandwidth_per_pixel: float = 1000, # Hz/pixel + show_plots: bool = True, + test_report: bool = True, + timing_check: bool = True, +) -> tuple[pp.Sequence, Path]: + """Generate a WASABI sequence for simultaneous B0 and B1 mapping. + + Parameters + ---------- + system + PyPulseq system limits object. + frequency_offsets + Array of frequency offsets (in Hz). + norm_offset + Frequency offset of normalization offset (in Hz). Defaults to -35e3 Hz. + t_recovery + Recovery time between frequency offsets (in seconds). + t_recovery_norm + Recovery time before normalization offset (in seconds). + fov_xy + Field of view in x and y direction (in meters). + n_readout + Number of frequency encoding steps. + n_phase_encoding + Number of phase encoding steps. + slice_thickness + Slice thickness of the 2D slice (in meters). + receiver_bandwidth_per_pixel + Desired receiver bandwidth per pixel (in Hz/pixel). This is used to calculate the readout duration. + show_plots + Toggles sequence plot. + test_report + Toggles advanced test report. + timing_check + Toggles timing check of the sequence. + + Returns + ------- + seq + PyPulseq Sequence object. + file_path + Path to the sequence file without suffix (append '.seq' for the actual file). + """ + if system is None: + system = sys_defaults + + if frequency_offsets is None: + frequency_offsets = np.linspace(-240, 240, 31) + + # define settings of rf excitation pulse + rf_duration = 1.28e-3 # duration of the rf excitation pulse [s] + rf_flip_angle = 10.0 # flip angle of rf excitation pulse [°] + rf_bwt = 4 # bandwidth-time product of rf excitation pulse [Hz*s] + rf_apodization = 0.5 # apodization factor of rf excitation pulse + readout_oversampling = 2 # readout oversampling factor, commonly 2. This reduces aliasing artifacts. + rf_spoiling_inc = 117 # RF spoiling phase increment [°]. Set to 0 for no RF spoiling. + + # define ADC and gradient timing + n_readout_with_oversampling = int(n_readout * readout_oversampling) + adc_dwell_time = round_to_raster( + 1.0 / (receiver_bandwidth_per_pixel * n_readout_with_oversampling), system.adc_raster_time + ) + _, adc_dwell_time = find_gx_flat_time_on_adc_raster( + n_readout_with_oversampling, adc_dwell_time, system.grad_raster_time, system.adc_raster_time + ) + + # WASABI block pulse + rf_prep_duration: float = 5e-3 + rf_prep_amplitude: float = 3.75 + + seq = wasabi_gre_centric_kernel( + system=system, + frequency_offsets=frequency_offsets, + norm_offset=norm_offset, + t_recovery=t_recovery, + t_recovery_norm=t_recovery_norm, + fov_xy=fov_xy, + n_readout=n_readout, + readout_oversampling=readout_oversampling, + n_phase_encoding=n_phase_encoding, + slice_thickness=slice_thickness, + rf_prep_duration=rf_prep_duration, + rf_prep_amplitude=rf_prep_amplitude, + rf_duration=rf_duration, + rf_flip_angle=rf_flip_angle, + rf_bwt=rf_bwt, + rf_apodization=rf_apodization, + rf_spoiling_inc=rf_spoiling_inc, + adc_dwell_time=adc_dwell_time, + ) + + # check timing of the sequence + if timing_check: + ok, error_report = seq.check_timing() + if ok: + print('\nTiming check passed successfully') + else: + print('\nTiming check failed! Error listing follows\n') + print(error_report) + + # show advanced rest report + if test_report: + print('\nCreating advanced test report...') + print(seq.test_report()) + + # define sequence filename + filename = f'{Path(__file__).stem}_{int(fov_xy * 1000)}fov_{n_readout}nx_{n_phase_encoding}ny' + filename += f'_{len(seq.definitions["frequency_offsets"])}offsets' + + # save seq-file to disk + output_path = Path.cwd() / 'output' + output_path.mkdir(parents=True, exist_ok=True) + print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") + write_sequence(seq, str(output_path / filename), create_signature=True) + + if show_plots: + seq.plot() + + return seq, output_path / filename + + +if __name__ == '__main__': + main() diff --git a/tests/scripts/test_b0_b1_wasabi.py b/tests/scripts/test_b0_b1_wasabi.py new file mode 100644 index 0000000..78efe09 --- /dev/null +++ b/tests/scripts/test_b0_b1_wasabi.py @@ -0,0 +1,13 @@ +"""Tests for WASABI sequence.""" + +import pytest +from mrseq.scripts.b0_b1_wasabi import main as create_seq + +EXPECTED_DUR = 122.1104 # defined 2026-02-10 + + +def test_default_seq_duration(system_defaults): + """Test if default values result in expected sequence duration.""" + seq, _ = create_seq(system=system_defaults, show_plots=False) + duration = seq.duration()[0] + assert duration == pytest.approx(EXPECTED_DUR)