"""M2D spiral FLASH sequence.""" from pathlib import Path from typing import Literal import ismrmrd import numpy as np import pypulseq as pp from raw2ismrmrd.utils import Fov from raw2ismrmrd.utils import Limits from raw2ismrmrd.utils import MatrixSize from raw2ismrmrd.utils import create_header from mrseq.utils import round_to_raster from mrseq.utils import spiral_acquisition from mrseq.utils import sys_defaults from mrseq.utils import write_sequence def spiral_flash_kernel( system: pp.Opts, te: float | None, tr: float | None, fov_xy: float, n_readout: int, readout_oversampling: Literal[1, 2, 4], spiral_undersampling: float, slice_thickness: float, n_slices: int, n_dummy_excitations: int, gx_pre_duration: float, rf_duration: float, rf_flip_angle: float, rf_bwt: float, rf_apodization: float, rf_spoiling_phase_increment: float, gz_spoil_duration: float, gz_spoil_area: float, mrd_header_file: str | Path | None, ) -> tuple[pp.Sequence, float, float]: """Generate a spiral FLASH sequence. Parameters ---------- system PyPulseq system limits object. te Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. tr Desired repetition time (TR) (in seconds). Minimum repetition time is used if set to None. fov_xy Field of view in x and y direction (in meters). n_readout Number of frequency encoding steps. readout_oversampling Readout oversampling. Determines the number of ADC samples along a spiral and the bandwidth. spiral_undersampling Angular undersampling of the spiral trajectoy. slice_thickness Slice thickness of the 2D slice (in meters). n_slices Number of slices. n_dummy_excitations Number of dummy excitations before data acquisition to ensure steady state. gx_pre_duration Duration of readout pre-winder gradient (in seconds) 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_phase_increment RF spoiling phase increment (in degrees). Set to 0 for no RF spoiling. gz_spoil_duration Duration of spoiler gradient (in seconds) gz_spoil_area Area of spoiler gradient (in mT/m * s) mrd_header_file Filename of the ISMRMRD header file to be created. If None, no header file is created. Returns ------- seq PyPulseq Sequence object min_te Shortest possible echo time. min_tr Shortest possible repetition time. """ if n_dummy_excitations < 0: raise ValueError('Number of dummy excitations must be >= 0.') # create PyPulseq Sequence object and set system limits seq = pp.Sequence(system=system) # create slice selective excitation pulse and gradients rf, gz, gzr = pp.make_sinc_pulse( flip_angle=rf_flip_angle / 180 * np.pi, 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', ) # create readout gradient and ADC gx, gy, adc, trajectory, time_to_echo = spiral_acquisition( system, n_readout, fov_xy, spiral_undersampling, readout_oversampling=readout_oversampling, n_spirals=None, max_pre_duration=gx_pre_duration, spiral_type='out', ) max_spiral_duration = max(pp.calc_duration(gx, gy) for gx, gy in zip(gx, gy, strict=True)) # create spoiler gradients gz_spoil = pp.make_trapezoid(channel='z', system=system, area=gz_spoil_area, duration=gz_spoil_duration) 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) # slice rewinder + time_to_echo ) # calculate echo time delay (te_delay) if te is None: te_delay = 0.0 else: te_delay = round_to_raster(te - min_te, system.block_duration_raster) if te_delay < 0: raise ValueError(f'TE must be larger than {min_te * 1000:.3f} ms. Current value is {te * 1000:.3f} ms.') # calculate minimum repetition time min_tr = ( pp.calc_duration(gz) # rf pulse + pp.calc_duration(gzr) # slice rewinder + max_spiral_duration # readout gradient + pp.calc_duration(gz_spoil) # gradient spoiler or readout-re-winder ) # calculate repetition time delay (tr_delay) current_min_tr = min_tr + te_delay if tr is None: tr_delay = 0.0 else: tr_delay = round_to_raster(tr - current_min_tr, system.block_duration_raster) if tr_delay < 0: raise ValueError( f'TR must be larger than {current_min_tr * 1000:.3f} ms. Current value is {tr * 1000:.3f} ms.' ) print(f'\nCurrent echo time = {(min_te + te_delay) * 1000:.3f} ms') print(f'Current repetition time = {(current_min_tr + tr_delay) * 1000:.3f} ms') # choose initial rf phase offset rf_phase = 0.0 rf_inc = 0.0 # create header if mrd_header_file: hdr = create_header( traj_type='other', encoding_fov=Fov(x=fov_xy, y=fov_xy, z=slice_thickness), recon_fov=Fov(x=fov_xy, y=fov_xy, z=slice_thickness), encoding_matrix=MatrixSize(n_x=int(n_readout), n_y=int(n_readout), n_z=1), recon_matrix=MatrixSize(n_x=n_readout, n_y=n_readout, n_z=1), dwell_time=adc.dwell, slice_limits=Limits(min=0, max=n_slices, center=0), k1_limits=Limits(min=0, max=len(gx), center=0), k2_limits=Limits(min=0, max=1, center=0), ) # write header to file prot = ismrmrd.Dataset(mrd_header_file, 'w') prot.write_xml_header(hdr.toXML('utf-8')) # obtain noise samples seq.add_block(pp.make_label(label='LIN', type='SET', value=0), pp.make_label(label='SLC', type='SET', value=0)) seq.add_block(adc, pp.make_label(label='NOISE', type='SET', value=True)) seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) seq.add_block(pp.make_delay(system.rf_dead_time)) if mrd_header_file: acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=adc.num_samples) prot.append_acquisition(acq) for slice_ in range(n_slices): for spiral_ in range(-n_dummy_excitations, len(gx)): # calculate current phase_offset if rf_spoiling is activated if rf_spoiling_phase_increment > 0: rf.phase_offset = rf_phase / 180 * np.pi adc.phase_offset = rf_phase / 180 * np.pi # set frequency offset for current slice rf.freq_offset = gz.amplitude * slice_thickness * (slice_ - (n_slices - 1) / 2) # add slice selective excitation pulse seq.add_block(rf, gz) seq.add_block(gzr) # update rf phase offset for the next excitation pulse rf_inc = divmod(rf_inc + rf_spoiling_phase_increment, 360.0)[1] rf_phase = divmod(rf_phase + rf_inc, 360.0)[1] if te_delay > 0: seq.add_block(pp.make_delay(te_delay)) if spiral_ >= 0: labels = [] labels.append(pp.make_label(label='LIN', type='SET', value=spiral_)) labels.append(pp.make_label(label='SLC', type='SET', value=slice_)) seq.add_block(gx[spiral_], gy[spiral_], adc, pp.make_delay(max_spiral_duration), *labels) else: seq.add_block(pp.make_delay(max_spiral_duration)) seq.add_block(gz_spoil) # add delay in case TR > min_TR if tr_delay > 0: seq.add_block(pp.make_delay(tr_delay)) if mrd_header_file and spiral_ >= 0: # add acquisitions to metadata spiral_trajectory = np.zeros((trajectory.shape[1], 2), dtype=np.float32) # the spiral trajectory is calculated in units of delta_k. for image reconstruction we use delta_k = 1 spiral_trajectory[:, 0] = trajectory[spiral_, :, 0] * fov_xy spiral_trajectory[:, 1] = trajectory[spiral_, :, 1] * fov_xy acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=adc.num_samples) acq.traj[:] = spiral_trajectory prot.append_acquisition(acq) # close ISMRMRD file if mrd_header_file: prot.close() # write all required parameters in the seq-file header/definitions seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness * n_slices]) seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) seq.set_definition('SliceThickness', slice_thickness) seq.set_definition('TE', te or min_te) seq.set_definition('TR', tr or min_tr) return seq, min_te, min_tr def main( system: pp.Opts | None = None, te: float | None = None, tr: float | None = None, rf_flip_angle: float = 12, fov_xy: float = 128e-3, n_readout: int = 128, readout_oversampling: Literal[1, 2, 4] = 2, n_spiral_arms: int = 128, slice_thickness: float = 8e-3, n_slices: int = 1, n_dummy_excitations: int = 20, show_plots: bool = True, test_report: bool = True, timing_check: bool = True, v141_compatibility: bool = True, ) -> tuple[pp.Sequence, Path]: """Generate a spiral FLASH sequence. Parameters ---------- system PyPulseq system limits object. te Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. tr Desired repetition time (TR) (in seconds). Minimum repetition time is used if set to None. rf_flip_angle Flip angle of rf excitation pulse (in degrees) fov_xy Field of view in x and y direction (in meters). n_readout Number of frequency encoding steps. readout_oversampling Readout oversampling. Determines the number of ADC samples along a spiral and the bandwidth. n_spiral_arms Number of spiral arms. slice_thickness Slice thickness of the 2D slice (in meters). n_slices Number of slices. n_dummy_excitations Number of dummy excitations before data acquisition to ensure steady state. show_plots Toggles sequence plot. test_report Toggles advanced test report. timing_check Toggles timing check of the sequence. v141_compatibility Save the sequence in pulseq v1.4.1 for backwards compatibility. Returns ------- seq Sequence object of spiral FLASH sequence. file_path Path to the sequence file. """ if system is None: system = sys_defaults # define settings of rf excitation pulse rf_duration = 1.28e-3 # duration of the rf excitation pulse [s] 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. # gradient timing gx_pre_duration = 1.0e-3 # duration of readout pre-winder gradient [s] # define spoiling gz_spoil_duration = 0.8e-3 # duration of spoiler gradient [s] gz_spoil_area = 4 / slice_thickness # area / zeroth gradient moment of spoiler gradient rf_spoiling_phase_increment = 117 # RF spoiling phase increment [°]. Set to 0 for no RF spoiling. # define sequence filename filename = f'{Path(__file__).stem}_{int(fov_xy * 1000)}fov_{n_readout}nx_{n_spiral_arms}na_{n_slices}ns' output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) # delete existing header file if (output_path / Path(filename + '_header.h5')).exists(): (output_path / Path(filename + '_header.h5')).unlink() seq, _min_te, min_tr = spiral_flash_kernel( system=system, te=te, tr=tr, fov_xy=fov_xy, n_readout=n_readout, readout_oversampling=readout_oversampling, spiral_undersampling=n_readout / n_spiral_arms, slice_thickness=slice_thickness, n_slices=n_slices, n_dummy_excitations=n_dummy_excitations, gx_pre_duration=gx_pre_duration, rf_duration=rf_duration, rf_flip_angle=rf_flip_angle, rf_bwt=rf_bwt, rf_apodization=rf_apodization, rf_spoiling_phase_increment=rf_spoiling_phase_increment, gz_spoil_duration=gz_spoil_duration, gz_spoil_area=gz_spoil_area, mrd_header_file=output_path / Path(filename + '_header.h5'), ) # check timing of the sequence if timing_check and not test_report: 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()) # save seq-file to disk print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) if show_plots: seq.plot(time_range=(0, 10 * (tr or min_tr))) return seq, output_path / filename if __name__ == '__main__': main()