diff --git a/examples/_toc.yml b/examples/_toc.yml index 14a0273..c6dc87c 100644 --- a/examples/_toc.yml +++ b/examples/_toc.yml @@ -14,4 +14,7 @@ chapters: - file: t2star_multi_echo_flash - file: t1rho_se_single_line - file: t1_t2_spiral_cmrf +- file: system.md + sections: + - file: b0_b1_wasabi diff --git a/examples/b0_b1_wasabi.ipynb b/examples/b0_b1_wasabi.ipynb new file mode 100644 index 0000000..bfc253e --- /dev/null +++ b/examples/b0_b1_wasabi.ipynb @@ -0,0 +1,275 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# B0 and B1 Mapping using WASABI\n", + "\n", + "Simultaneous mapping of water shift (B0) and B1 (WASABI) using using an off-resonant preparation pulse. This induces Rabi oscillations which can be seen as sinc-like oscillations in the frequency-offset dimension. B0 can then be estimated by its symmetry axis and B1 by its oscillation frequency." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import tempfile\n", + "from pathlib import Path\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import MRzeroCore as mr0\n", + "import numpy as np\n", + "import torch\n", + "from einops import rearrange\n", + "from mrpro.algorithms.reconstruction import DirectReconstruction\n", + "from mrpro.data import KData\n", + "from mrpro.data import SpatialDimension\n", + "from mrpro.data.traj_calculators import KTrajectoryCartesian\n", + "from mrpro.operators import DictionaryMatchOp\n", + "from mrpro.operators.models import WASABI\n", + "\n", + "from mrseq.scripts.b0_b1_wasabi import main as create_seq\n", + "from mrseq.utils import sys_defaults" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Settings\n", + "We are going to use a numerical phantom with a matrix size of 128 x 128. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "image_matrix_size = [128, 128]\n", + "\n", + "frequency_offsets = np.linspace(-240, 240, 31)\n", + "norm_offset = -35e3\n", + "\n", + "tmp = tempfile.TemporaryDirectory()\n", + "fname_mrd = Path(tmp.name) / 'b0_b1_wasabi.mrd'" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### Create the digital phantom\n", + "\n", + "We use the standard Brainweb phantom from [MRzero](https://github.com/MRsources/MRzero-Core)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "phantom = mr0.util.load_phantom(image_matrix_size)" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "### Create the WASABI sequence\n", + "\n", + "To create the WASABI sequence, we use the previously imported [b0_b1_wasabi script](../src/mrseq/scripts/b0_b1_wasabi.py).\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "sequence, fname_seq = create_seq(\n", + " system=sys_defaults,\n", + " test_report=False,\n", + " timing_check=False,\n", + " fov_xy=float(phantom.size.numpy()[0]),\n", + " n_readout=image_matrix_size[0],\n", + " n_phase_encoding=image_matrix_size[0],\n", + " frequency_offsets=frequency_offsets,\n", + " norm_offset=norm_offset,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### Simulate the sequence\n", + "Now, we pass the sequence and the phantom to the MRzero simulation and save the simulated signal as an (ISMR)MRD file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))\n", + "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e-1)\n", + "mr0.sig_to_mrd(fname_mrd, signal, sequence)" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "### Reconstruct the images at different inversion times\n", + "\n", + "We use [MRpro](https://github.com/PTB-MR/MRpro) for the image reconstruction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "kdata = KData.from_file(fname_mrd, trajectory=KTrajectoryCartesian())\n", + "kdata.header.encoding_matrix = SpatialDimension(z=1, y=image_matrix_size[1], x=image_matrix_size[0] * 2)\n", + "kdata.header.recon_matrix = SpatialDimension(z=1, y=image_matrix_size[1], x=image_matrix_size[0])\n", + "recon = DirectReconstruction(kdata, csm=None)\n", + "idata = recon(kdata)" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "We can now plot the images at different offset frequencies." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "idat = idata.rss().abs().numpy().squeeze()\n", + "offsets = np.concatenate(([norm_offset], frequency_offsets))\n", + "fig, ax = plt.subplots(4, idat.shape[0] // 4, figsize=(3 * idat.shape[0] // 4, 3 * 4))\n", + "ax = ax.flatten()\n", + "for i in range(idat.shape[0]):\n", + " ax[i].imshow(idat[i, :, :], cmap='gray')\n", + " ax[i].set_title(f'Offset = {int(offsets[i])} Hz')\n", + " ax[i].set_xticks([])\n", + " ax[i].set_yticks([])" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "### Estimate B0 and B1\n", + "We use a dictionary matching approach to estimate the T1 maps. Afterward, we compare them to the input and ensure they match." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "dictionary = DictionaryMatchOp(WASABI(offsets=torch.as_tensor(offsets, dtype=torch.float32)))\n", + "dictionary.append(\n", + " torch.linspace(-100.0, 100.0, 100),\n", + " torch.linspace(0.1, 1.5, 100)[None, :],\n", + " torch.linspace(-100.0, 100.0, 100)[None, None, :],\n", + " torch.linspace(-100.0, 100.0, 100)[None, None, None, :],\n", + ")\n", + "b0_match, rb1_match, c_match, a_match = dictionary(idata.data)\n", + "\n", + "b0_input = np.roll(rearrange(phantom.B0.numpy().squeeze()[::-1, ::-1], 'x y -> y x'), shift=(1, 1), axis=(0, 1))\n", + "obj_mask = np.zeros_like(b0_input)\n", + "obj_mask[b0_input > 0] = 1\n", + "b0_measured = b0_match.numpy().squeeze() * obj_mask\n", + "\n", + "rb1_input = np.abs(\n", + " np.roll(rearrange(phantom.B1.numpy().squeeze()[::-1, ::-1], 'x y -> y x'), shift=(1, 1), axis=(0, 1))\n", + ")\n", + "rb1_measured = rb1_match.abs().numpy().squeeze() * obj_mask\n", + "\n", + "fig, ax = plt.subplots(1, 4)\n", + "ax[0].imshow(b0_match.squeeze().numpy() * obj_mask)\n", + "ax[1].imshow(rb1_match.abs().squeeze().numpy() * obj_mask)\n", + "ax[2].imshow(c_match.squeeze().numpy() * obj_mask)\n", + "ax[3].imshow(a_match.squeeze().numpy() * obj_mask)\n", + "\n", + "fig, ax = plt.subplots(2, 3, figsize=(15, 6))\n", + "for cax in ax.flatten():\n", + " cax.set_xticks([])\n", + " cax.set_yticks([])\n", + "\n", + "im = ax[0, 0].imshow(b0_input, vmin=-40, vmax=40, cmap='bwr')\n", + "fig.colorbar(im, ax=ax[0, 0], label='Input B0 (Hz)')\n", + "\n", + "im = ax[0, 1].imshow(b0_measured, vmin=-1, vmax=1, cmap='bwr')\n", + "fig.colorbar(im, ax=ax[0, 1], label='Measured B0 (Hz)')\n", + "\n", + "im = ax[0, 2].imshow(b0_measured - b0_input, vmin=-1.8, vmax=1.8, cmap='bwr')\n", + "fig.colorbar(im, ax=ax[0, 2], label='Difference B0 (Hz)')\n", + "\n", + "relative_error = np.sum(np.abs(b0_input - b0_measured)) / np.sum(np.abs(b0_input))\n", + "print(f'Relative error {relative_error}')\n", + "# assert relative_error < 0.08\n", + "\n", + "im = ax[1, 0].imshow(rb1_input, vmin=0, vmax=1.2, cmap='magma')\n", + "fig.colorbar(im, ax=ax[1, 0], label='Input relative B1 ')\n", + "\n", + "im = ax[1, 1].imshow(rb1_measured, vmin=-1, vmax=1, cmap='magma')\n", + "fig.colorbar(im, ax=ax[1, 1], label='Measured relative B1')\n", + "\n", + "im = ax[1, 2].imshow(rb1_measured - rb1_input, vmin=-1.8, vmax=1.8, cmap='bwr')\n", + "fig.colorbar(im, ax=ax[1, 2], label='Difference relative B1')\n", + "\n", + "relative_error = np.sum(np.abs(rb1_input - rb1_measured)) / np.sum(np.abs(rb1_input))\n", + "print(f'Relative error {relative_error}')\n", + "# assert relative_error < 0.08" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mrseq_mrzero", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/system.md b/examples/system.md new file mode 100644 index 0000000..3cf58df --- /dev/null +++ b/examples/system.md @@ -0,0 +1,5 @@ +# System Characterization + +Sequences to characterize the MR system. Available examples are: + +- [B0 and B1 mapping using WASABI](b0_b1_wasabi.ipynb) diff --git a/examples/t1_molli_bssfp.ipynb b/examples/t1_molli_bssfp.ipynb index b4631fd..6156a9c 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." ] }, { @@ -64,6 +64,14 @@ "id": "4", "metadata": {}, "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], "source": [ "image_matrix_size = [64, 64]\n", "\n", @@ -73,7 +81,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "6", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -84,7 +92,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "7", "metadata": {}, "outputs": [], "source": [ @@ -98,7 +106,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "8", "metadata": {}, "source": [ "### Create the MOLLI bSSFP sequence\n", @@ -109,7 +117,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "9", "metadata": {}, "outputs": [], "source": [ @@ -125,7 +133,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "10", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -135,7 +143,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "11", "metadata": {}, "outputs": [], "source": [ @@ -146,7 +154,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "12", "metadata": {}, "source": [ "### Reconstruct the images at different inversion times\n", @@ -157,7 +165,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "13", "metadata": {}, "outputs": [], "source": [ @@ -178,7 +186,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "14", "metadata": {}, "source": [ "We visualize the coil sensitivity maps which we used for the iterative SENSE reconstruction." @@ -187,7 +195,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "15", "metadata": {}, "outputs": [], "source": [ @@ -198,7 +206,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "16", "metadata": {}, "source": [ "We can now plot the images at different inversion times." @@ -207,7 +215,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "17", "metadata": {}, "outputs": [], "source": [ @@ -229,7 +237,7 @@ }, { "cell_type": "markdown", - "id": "17", + "id": "18", "metadata": {}, "source": [ "### Estimate the T1 maps\n", @@ -239,7 +247,7 @@ { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "19", "metadata": {}, "outputs": [], "source": [ diff --git a/src/mrseq/scripts/b0_b1_wasabi.py b/src/mrseq/scripts/b0_b1_wasabi.py new file mode 100644 index 0000000..9185dca --- /dev/null +++ b/src/mrseq/scripts/b0_b1_wasabi.py @@ -0,0 +1,349 @@ +"""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 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, + 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. + 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) + + # 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 + gx = pp.make_trapezoid( + channel='x', flat_area=n_readout * delta_k, flat_time=adc_dwell_time * n_readout, system=system + ) + adc = pp.make_adc(num_samples=n_readout, duration=gx.flat_time, delay=gx.rise_time, system=system) + + # create frequency encoding pre-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) + + # calculate gradient areas for (linear) phase encoding direction + k0_center_id = np.where((np.arange(n_readout) - n_readout / 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()) + + 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, + rf_prep_duration: float = 5e-3, + rf_prep_amplitude: float = 3.75, + rf_duration: float = 1.28e-3, + rf_flip_angle: float = 10.0, + rf_bwt: float = 4.0, + rf_apodization: float = 0.5, + rf_spoiling_inc: float = 117.0, + adc_dwell_time: float = 1e-5, + 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). + 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 + ADC dwell time (in seconds). + 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) + + 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, + 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..c85bbc7 --- /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-05 + + +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)