diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 038c78ae..c100f40f 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -15,6 +15,16 @@ repos: exclude: ^.vscode/ - id: mixed-line-ending + - repo: https://github.com/kynan/nbstripout + rev: 0.9.0 + hooks: + # cleans the .ipynbs (removes outputs, resets all cell-ids to 0..N, cleans steps) + # also clean any kernel information left after execution + - id: nbstripout + name: clean .ipynb output + args: [--extra-keys, "metadata.language_info"] + files: examples/.*\.ipynb$ + - repo: https://github.com/astral-sh/ruff-pre-commit rev: v0.14.14 hooks: diff --git a/examples/notebooks/write_t2_se.ipynb b/examples/notebooks/write_t2_se.ipynb index d58768ed..8e667655 100644 --- a/examples/notebooks/write_t2_se.ipynb +++ b/examples/notebooks/write_t2_se.ipynb @@ -25,7 +25,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -48,7 +48,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -83,7 +83,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -116,7 +116,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -148,7 +148,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -175,7 +175,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -223,7 +223,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -249,7 +249,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -275,7 +275,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -299,7 +299,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -334,7 +334,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -382,7 +382,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -406,7 +406,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -421,7 +421,7 @@ }, { "cell_type": "code", - "execution_count": 0, + "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", @@ -442,18 +442,6 @@ "display_name": "Python 3", "language": "python", "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.6.3" } }, "nbformat": 4, diff --git a/examples/tutorials/01_introduction.ipynb b/examples/tutorials/01_introduction.ipynb new file mode 100644 index 00000000..2754eef7 --- /dev/null +++ b/examples/tutorials/01_introduction.ipynb @@ -0,0 +1,620 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Tutorial 1: Introduction to PyPulseq\n", + "\n", + "Welcome to the PyPulseq tutorial series! This first notebook introduces the core concepts\n", + "and building blocks you need to start creating MRI pulse sequences with PyPulseq.\n", + "\n", + "## What is Pulseq?\n", + "\n", + "[Pulseq](https://pulseq.github.io/) is an open-source, vendor-independent file format for\n", + "MRI pulse sequence programming. Instead of writing sequences in vendor-specific languages\n", + "(e.g., Siemens IDEA, GE EPIC, Bruker ParaVision), you describe your sequence in a `.seq`\n", + "file that can be executed on any scanner with a Pulseq interpreter.\n", + "\n", + "## What is PyPulseq?\n", + "\n", + "**PyPulseq** is a Python package for creating Pulseq `.seq` files. It provides functions to\n", + "define RF pulses, gradient waveforms, ADC readouts, and delays — and to assemble them into\n", + "complete pulse sequences.\n", + "\n", + "## What you will learn\n", + "\n", + "In this notebook, you will learn:\n", + "\n", + "1. How to install and import PyPulseq\n", + "2. How to define scanner hardware limits with `pp.Opts`\n", + "3. How to create a `Sequence` object\n", + "4. How to create basic sequence events (RF pulse, ADC, delay)\n", + "5. How to assemble events into blocks and inspect timing\n", + "6. How to check, plot, and export the sequence\n", + "7. How to investigate the `.seq` file format" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 1. Installation\n", + "\n", + "There are several ways to install PyPulseq. The recommended and most convenient way is to install the latest stable release version from [PyPI](https://pypi.org/project/pypulseq/) using `pip`:\n", + "\n", + "```bash\n", + "pip install pypulseq\n", + "```\n", + "\n", + "If you prefer to use conda, you can install the latest stable release version from [conda-forge](https://anaconda.org/conda-forge/pypulseq) using `conda`:\n", + "```bash\n", + "conda install -c conda-forge pypulseq\n", + "```\n", + "\n", + "If you want to install the latest development version, you can install it directly from the GitHub repository:\n", + "```bash\n", + "pip install git+https://github.com/pypulseq/pypulseq.git@\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('pypulseq'):\n", + " %pip install pypulseq" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 2. Importing PyPulseq\n", + "\n", + "By convention, we import PyPulseq under the alias `pp`. We usually also import NumPy as `np`, which is used throughout for numerical operations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "import pypulseq as pp\n", + "\n", + "print(f'PyPulseq version: {pp.__version__}')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 3. Scanner Hardware Limits: `pp.Opts`\n", + "\n", + "Every MRI scanner has physical limits on how fast gradients can change (slew rate), how\n", + "strong they can be (maximum gradient amplitude), how long the gaps between individual events of different types should be\n", + "(dead times), and on which time grid events should be scheduled (raster times). Before creating any sequence events, we\n", + "define these limits in an [Opts](../../src/pypulseq/opts.py) object.\n", + "\n", + "This ensures that all events we create as well as the concatenation of events (the sequence) will respect the scanner's capabilities.\n", + "\n", + "Typical values we set for individual scanners are:\n", + "\n", + "| Parameter | Description | Typical value |\n", + "|---|---|---|\n", + "| `max_grad` | Maximum gradient amplitude | 24–80 mT/m |\n", + "| `max_slew` | Maximum gradient slew rate | 100–200 T/m/s |\n", + "| `rf_dead_time` | Dead time before RF pulse | 100 µs |\n", + "| `rf_ringdown_time` | Ringdown time after RF pulse | 10–30 µs |\n", + "| `adc_dead_time` | Dead time before ADC readout | 10 µs |\n", + "\n", + "You can find a complete documentation of all possible arguments in the docstring of the [Opts](../../src/pypulseq/opts.py) class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "system = pp.Opts(\n", + " max_grad=28,\n", + " grad_unit='mT/m',\n", + " max_slew=150,\n", + " slew_unit='T/m/s',\n", + " rf_ringdown_time=20e-6,\n", + " rf_dead_time=100e-6,\n", + " adc_dead_time=10e-6,\n", + ")\n", + "\n", + "print(f'Max gradient: {system.max_grad / system.gamma * 1e3:.1f} mT/m')\n", + "print(f'Max slew rate: {system.max_slew / system.gamma:.1f} T/m/s')\n", + "print(f'Gradient raster time: {system.grad_raster_time * 1e6:.1f} µs')\n", + "print(f'RF raster time: {system.rf_raster_time * 1e6:.1f} µs')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 4. The Sequence Object\n", + "\n", + "The `Sequence` object is the central container for every (Py)Pulseq pulse sequence. Usually, we name it `seq` and create it by\n", + "passing the system limits from the previous cell. \n", + "\n", + "In (Py)Pulseq, a sequence is a collection / concatenation of non-overlapping **blocks**, which we will cover later in this tutorial." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "seq = pp.Sequence(system)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 5. Creating Basic Events\n", + "\n", + "PyPulseq provides `make_*` functions to create different types of events. In this\n", + "introductory tutorial, we'll start with the simplest ones:\n", + "\n", + "- **`make_block_pulse`** — a rectangular (block) RF pulse\n", + "- **`make_adc`** — a readout / ADC event\n", + "- **`make_delay`** — a simple delay\n", + "\n", + "These events already allow us to create the most basic sequence, a Free Induction Decay (FID) pulse sequence.\n", + "\n", + "Later tutorials will also cover different types of gradients, shaped RF pulses, and more." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 5.1 Creating an RF pulse\n", + "\n", + "The first important event for our FID pulse sequence is the RF pulse used to excite the spin system. The simplest RF pulse is a rectangular (block) pulse. \n", + "\n", + "We specify:\n", + "- **`flip_angle`** — the desired flip angle in **radians** (use `np.deg2rad()` to convert from degrees)\n", + "- **`duration`** — the pulse duration in seconds\n", + "- **`system`** — the system limits object\n", + "- **`delay`** — the delay before the pulse in seconds (this has to be >= `system.rf_dead_time`)\n", + "- **`use`** — the purpose of the pulse (e.g., `'excitation'` or `'refocusing'`)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a non-selective 90° block pulse with 1 ms duration\n", + "rf_90 = pp.make_block_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " duration=1e-3,\n", + " system=system,\n", + " delay=system.rf_dead_time,\n", + " use='excitation',\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 5.2 Creating a readout / ADC event\n", + "\n", + "The second event we need is an ADC (Analog-to-Digital Converter) event. An ADC event is used to read out the signal from the receiver coil(s).\n", + "\n", + "We specify:\n", + "\n", + "- **`num_samples`** — Number of readout samples\n", + "- **`delay`** — Delay before the ADC event in seconds (this has to be >= `system.adc_dead_time`)\n", + "- **`dwell`** — The time between two samples in seconds\n", + "- **`system`** — System limits\n", + "\n", + "Please note that there are different ways to specify the same ADC event. For example, you can either specify the number of samples and the dwell time, or the number of samples and the total duration. Additionally, parameters such as `freq_offset`, `phase_offset`, `freq_ppm`, and `phase_ppm` can be specified to control the frequency and phase of the ADC event. We don't need these for the simple FID pulse sequence, but we will cover them in later tutorials." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a simple ADC event with 1024 samples and a dwell time of 20 µs\n", + "adc = pp.make_adc(\n", + " num_samples=1024,\n", + " delay=system.adc_dead_time,\n", + " dwell=20e-6,\n", + " system=system,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 5.3 Creating a delay\n", + "\n", + "A delay is simply a period of time where nothing happens. Delays are used to control the timing of the sequence (e.g., to achieve a specific TE or TR). For our simple example, we will create a delay of 100 ms using the `make_delay` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "delay = pp.make_delay(100e-3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 6. Building a Sequence: Blocks and `add_block`\n", + "\n", + "Now let's assemble these events into a sequence. As mentioned above, a Pulseq sequence is a concatenation of non-overlapping **blocks**.\n", + "\n", + "A **block** is a group of different events (e.g. RF pulse, gradient, ADC) that overlap in time within a block. \n", + "\n", + "Each block may contain:\n", + "\n", + "- one optional gradient per axis\n", + "- one optional RF pulse\n", + "- one optional ADC event\n", + "- one optional delay or soft delay\n", + "- one optional trigger\n", + "- various optional labels\n", + "\n", + "Individual events within a block may define their own start delays. We already used this feature above when defining the delays of the RF pulse and the ADC event.\n", + "\n", + "The key method to add events to a sequence is `seq.add_block()`. It takes one or more events as arguments and adds them as a block to the sequence." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 6.1 Building our simple FID sequence\n", + "\n", + "We can now build our simple FID sequence by adding our RF pulse (`rf_90`), the ADC event (`adc`), and the 100 ms delay (`delay`) to the Sequence object (`seq`)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a new sequence using the system limits\n", + "seq = pp.Sequence(system=system)\n", + "\n", + "# Add the RF pulse to the sequence\n", + "seq.add_block(rf_90)\n", + "\n", + "# Add the ADC event and the 100 ms delay to the sequence\n", + "seq.add_block(adc, delay)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 6.2 Inspecting event and block durations\n", + "\n", + "The `add_block` method can take a single event or multiple events as arguments. If you pass multiple events, they will be added to the same block, where they overlap in time. Therefore, the duration of a block is given by the longest event in the block. \n", + "\n", + "Let's investigate the duration of all the events and blocks we created so far using the `calc_duration` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "duration_rf_90 = pp.calc_duration(rf_90)\n", + "print(f'The duration of the RF pulse is {duration_rf_90 * 1e3:.2f} ms.')\n", + "\n", + "duration_adc = pp.calc_duration(adc)\n", + "print(f'The duration of the ADC is {duration_adc * 1e3:.2f} ms.')\n", + "\n", + "duration_delay = pp.calc_duration(delay)\n", + "print(f'The duration of the delay is {duration_delay * 1e3:.2f} ms.')\n", + "\n", + "duration_adc_and_delay = pp.calc_duration(adc, delay)\n", + "print(f'The duration of the ADC and delay is {duration_adc_and_delay * 1e3:.2f} ms.')\n", + "\n", + "print('')\n", + "\n", + "duration_block1 = seq.block_durations[1]\n", + "print(f'The duration of the first block is {duration_block1 * 1e3:.2f} ms.')\n", + "\n", + "duration_block2 = seq.block_durations[2]\n", + "print(f'The duration of the second block is {duration_block2 * 1e3:.2f} ms.')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Durations of the individual events\n", + "As you can see, the duration of our 1 ms RF pulse is NOT exactly 1 ms, but 1.12 ms. This is because the total duration of the RF pulse event includes the RF dead time before the RF pulse (`rf_dead_time=100 µs`) and the RF ringdown time after the RF pulse (`rf_ringdown_time=20 µs`).\n", + "\n", + "Something similar can be seen for the total duration of the ADC event. The readout duration is expected to be `num_samples * dwell = 20.48 ms`, but the total duration returned by the `calc_duration` function is 20.50 ms. The additional 20 µs are due to the ADC delay (10 µs from `adc_dead_time`) and rounding to the block duration raster.\n", + "\n", + "For the delay event, the total duration is simply the delay time we specified.\n", + "\n", + "From the last `calc_duration` call, you can see that the function also accepts several events at once. In this case, the returned duration is the maximum duration across all provided events, which is the duration of the delay event (100 ms) in this case.\n", + "\n", + "### Durations of the different blocks\n", + "The durations of all blocks of a sequence are accessible via the `seq.block_durations` dictionary, which uses the block numbers as keys and the durations as values.\n", + "\n", + "As we expected, the duration of the first block is given by the total duration (1.12 ms) of the RF pulse event it contains.\n", + "\n", + "The duration of the second block is given by the longest event it contains, which is the delay event (100 ms)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 7. Checking, Plotting, and Exporting\n", + "\n", + "### 7.1 Checking timing\n", + "\n", + "Before using a sequence, it's good practice to check for timing errors. The `check_timing()`\n", + "method verifies that all events are properly aligned to the raster times and that there are\n", + "no overlapping or conflicting events." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ok, error_report = seq.check_timing()\n", + "\n", + "if ok:\n", + " print('Timing check passed!')\n", + "else:\n", + " print('Timing check failed:')\n", + " print(error_report)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 7.2 Advanced test report\n", + "\n", + "(Py)Pulseq also provides a more detailed test report that will give you information about the number of different events in your sequence, the total duration, the estimated echo time (TE) and the repetition time (TR), the flip angle(s), the gradient amplitudes and slew rates in the different directions and so on. This is especially useful for more complex sequences, but for completeness, we also show it here." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(seq.test_report())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 7.3 Plotting the sequence\n", + "\n", + "The `seq.plot()` method provides a visual representation of the sequence, showing RF pulses, gradients, and ADC events over time. This can be helpful for debugging and verification, as well as for educational purposes. If we don't specify a time range, the plot will show the entire sequence." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "seq.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "However, sometimes it is more helpful to focus on a specific part of the sequence, e.g. the RF excitation pulse and the beginning of the ADC readout in our simple FID sequence. Let's zoom in on this part of the sequence by setting the `time_range` parameter of `seq.plot()` to the first 1.5 ms of the sequence, which allows us to see the RF excitation pulse and its preceding delay (RF dead time) and the first few of the 1024 ADC samples." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "seq.plot(time_range=(0.0, 1.5e-3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 7.4 Sequence definitions\n", + "\n", + "Before writing a sequence, you can attach metadata using `seq.set_definition()`. This information is stored in the `.seq` file header and can be used by the scanner or\n", + "reconstruction software.\n", + "\n", + "Common definitions include the field of view (FoV), a sequence name, the slice thickness, the echo time (TE), the repetition time (TR), etc. For our FID sequence, we will only define the name to be `simple_fid`, the echo time to be `0.54 ms` (as given in the test report) and a `random_parameter` with a value of `42`. We will go into more detail about which definitions are used by the scanner or reconstruction software in later tutorials." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "seq.set_definition(key='Name', value='simple_fid')\n", + "seq.set_definition(key='TE', value=0.54e-3)\n", + "seq.set_definition(key='random_parameter', value=42)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 7.5 Writing a `.seq` file\n", + "\n", + "To export the sequence as a `.seq` file for use on a scanner or simulation software, call `seq.write()`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "seq.write('simple_fid.seq')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 8. Investigating the `.seq` File\n", + "\n", + "The `.seq` file is a human-readable text file that contains all the information needed to play out the sequence on a scanner. Let's take a look at the content of the `simple_fid.seq` file we just created:\n", + "\n", + "> **`simple_fid.seq`**\n", + "> ```\n", + "> # Pulseq sequence file\n", + "> # Created by PyPulseq\n", + "> \n", + "> [VERSION]\n", + "> major 1\n", + "> minor 5\n", + "> revision 0\n", + "> \n", + "> [DEFINITIONS]\n", + "> AdcRasterTime 1e-07 \n", + "> BlockDurationRaster 1e-05 \n", + "> GradientRasterTime 1e-05 \n", + "> Name simple_fid \n", + "> RadiofrequencyRasterTime 1e-06 \n", + "> TE 0.00054 \n", + "> TotalDuration 0.10112 \n", + "> random_parameter 42 \n", + "> \n", + "> # Format of blocks:\n", + "> # NUM DUR RF GX GY GZ ADC EXT\n", + "> [BLOCKS]\n", + "> 1 112 1 0 0 0 0 0\n", + "> 2 10000 0 0 0 0 1 0\n", + "> \n", + "> # Format of RF events:\n", + "> # id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use\n", + "> # .. Hz .. .. .. us us ppm rad/MHz Hz rad ..\n", + "> # Field \"use\" is the initial of: excitation refocusing inversion saturation preparation other undefined\n", + "> [RF]\n", + "> 1 250 1 2 3 500 100 0 0 0 0 e\n", + "> \n", + "> # Format of ADC events:\n", + "> # id num dwell delay freqPPM phasePPM freq phase phase_id\n", + "> # .. .. ns us ppm rad/MHz Hz rad ..\n", + "> [ADC]\n", + "> 1 1024 20000 10 0 0 0 0 0\n", + "> \n", + "> # Sequence Shapes\n", + "> [SHAPES]\n", + "> \n", + "> shape_id 1\n", + "> num_samples 2\n", + "> 1\n", + "> 1\n", + "> \n", + "> shape_id 2\n", + "> num_samples 2\n", + "> 0\n", + "> 0\n", + "> \n", + "> shape_id 3\n", + "> num_samples 2\n", + "> 0\n", + "> 1000\n", + "> \n", + "> \n", + "> [SIGNATURE]\n", + "> # This is the hash of the Pulseq file, calculated right before the [SIGNATURE] section was added\n", + "> # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE]\n", + "> # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification)\n", + "> Type md5\n", + "> Hash f955d3c1205157caf089defa15528e7d\n", + "> ```\n", + "\n", + "A detailed description of the `.seq` file format can be found in the [specification.pdf](https://pulseq.github.io/specification.pdf) on the Pulseq website." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## Summary\n", + "\n", + "In this tutorial, you learned the fundamental building blocks of PyPulseq:\n", + "\n", + "| Concept | Function / Class |\n", + "|---|---|\n", + "| Hardware limits | `pp.Opts(max_grad=..., max_slew=..., ...)` |\n", + "| Sequence container | `pp.Sequence(system)` |\n", + "| Block RF pulse | `pp.make_block_pulse(flip_angle=..., duration=..., system=...)` |\n", + "| ADC readout | `pp.make_adc(num_samples=..., dwell=..., system=...)` |\n", + "| Delay | `pp.make_delay(duration)` |\n", + "| Add events to sequence | `seq.add_block(event1, event2, ...)` |\n", + "| Compute event duration | `pp.calc_duration(event1, event2, ...)` |\n", + "| Check timing | `seq.check_timing()` |\n", + "| Test report | `seq.test_report()` |\n", + "| Plot | `seq.plot()` |\n", + "| Metadata | `seq.set_definition(key=..., value=...)` |\n", + "| Export | `seq.write('filename.seq')` |\n", + "\n", + "## Next steps\n", + "\n", + "In the next tutorial, we'll explore **RF pulses** in more detail — including sinc pulses\n", + "for slice-selective excitation, and how RF pulses interact with slice-selection gradients." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "pypulseq", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/tutorials/02_rf_pulses.ipynb b/examples/tutorials/02_rf_pulses.ipynb new file mode 100644 index 00000000..88227508 --- /dev/null +++ b/examples/tutorials/02_rf_pulses.ipynb @@ -0,0 +1,1040 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Tutorial 2: RF Pulses\n", + "\n", + "In this tutorial, we explore the different types of radio-frequency (RF) pulses available in\n", + "PyPulseq. RF pulses are fundamental building blocks of every MRI pulse sequence — they are\n", + "used to excite, refocus, invert, or saturate the magnetization.\n", + "\n", + "## What you will learn\n", + "\n", + "1. How to create different RF pulse shapes (block, sinc, Gaussian)\n", + "2. How pulse parameters affect the waveform (duration, time-bandwidth product, apodization)\n", + "3. How to make RF pulses slice-selective using gradients\n", + "4. How to control flip angle, phase offset, and frequency offset\n", + "5. How to specify the purpose of an RF pulse (`use` parameter)\n", + "6. How to assemble RF pulses with gradients into a sequence" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('pypulseq'):\n", + " %pip install pypulseq" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "import pypulseq as pp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "system = pp.Opts(\n", + " max_grad=28,\n", + " grad_unit='mT/m',\n", + " max_slew=150,\n", + " slew_unit='T/m/s',\n", + " rf_ringdown_time=20e-6,\n", + " rf_dead_time=100e-6,\n", + " adc_dead_time=10e-6,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Before we start, let's define a small helper function to plot the RF pulse waveform. This\n", + "will allow us to quickly visualize and compare different RF pulses throughout this tutorial." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_rf(rf, title='', ax=None):\n", + " \"\"\"Plot the magnitude and phase of an RF pulse waveform.\n", + "\n", + " The plot includes the full event timing: dead time before the pulse,\n", + " the pulse waveform itself, and the ringdown time after the pulse.\n", + " \"\"\"\n", + " if ax is None:\n", + " _, ax = plt.subplots(2, 1, figsize=(8, 4), sharex=True)\n", + "\n", + " # Build the full time axis including dead time and ringdown\n", + " t_full = np.concatenate(\n", + " ([0], [rf.delay], rf.delay + rf.t, [rf.delay + rf.shape_dur], [rf.delay + rf.shape_dur + rf.ringdown_time])\n", + " )\n", + " signal_full = np.concatenate(([0], [0], rf.signal, [0], [0]))\n", + "\n", + " t_ms = t_full * 1e3 # Convert to ms\n", + "\n", + " ax[0].plot(t_ms, np.abs(signal_full), 'b-', linewidth=1.5)\n", + " ax[0].set_ylabel('Amplitude (Hz)')\n", + " ax[0].set_title(title)\n", + " ax[0].grid(True, alpha=0.3)\n", + "\n", + " ax[1].plot(t_ms, np.angle(signal_full, deg=True), 'r-', linewidth=1.5)\n", + " ax[1].set_ylabel('Phase (°)')\n", + " ax[1].set_xlabel('Time (ms)')\n", + " ax[1].grid(True, alpha=0.3)\n", + "\n", + " return ax" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 1. Block (Rectangular) Pulses\n", + "\n", + "The simplest RF pulse is a **block pulse** (also called a rectangular or hard pulse). It has\n", + "a constant amplitude for its entire duration. Block pulses are non-selective, meaning they\n", + "excite all spins equally regardless of their spatial position.\n", + "\n", + "We already used `pp.make_block_pulse` in Tutorial 1. Let's now look at it in more detail." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_block = pp.make_block_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=1e-3,\n", + " system=system,\n", + " use='excitation',\n", + ")\n", + "\n", + "print(f'RF amplitude: {np.max(np.abs(rf_block.signal)):.1f} Hz')\n", + "print(f'Shape duration: {rf_block.shape_dur * 1e3:.2f} ms')\n", + "print(f'Total duration (incl. dead time + ringdown): {pp.calc_duration(rf_block) * 1e3:.2f} ms')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "plot_rf(rf_block, title='90° Block Pulse (1 ms)')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As expected, the block pulse has a constant amplitude and zero phase. The amplitude is\n", + "determined by the flip angle and the duration: a shorter pulse requires a higher amplitude\n", + "to achieve the same flip angle.\n", + "\n", + "### 1.1 Effect of flip angle on amplitude\n", + "\n", + "Let's compare block pulses with different flip angles but the same duration." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "flip_angles_deg = [10, 30, 60, 90, 180]\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(8, 3))\n", + "for fa_deg in flip_angles_deg:\n", + " rf = pp.make_block_pulse(\n", + " flip_angle=np.deg2rad(fa_deg),\n", + " delay=system.rf_dead_time,\n", + " duration=1e-3,\n", + " system=system,\n", + " )\n", + " rf_signal = np.concatenate(([0], rf.signal, [0]))\n", + " rf_time = np.concatenate(([0], rf.t, [rf.t[-1]]))\n", + " ax.plot(rf_time * 1e3, rf_signal, label=f'{fa_deg}°', linewidth=1.5)\n", + "\n", + "ax.set_xlabel('Time (ms)')\n", + "ax.set_ylabel('Amplitude (Hz)')\n", + "ax.set_title('Block Pulses: Effect of Flip Angle')\n", + "ax.legend()\n", + "ax.grid(True, alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The amplitude scales linearly with the flip angle. A 180° pulse has exactly twice the\n", + "amplitude of a 90° pulse of the same duration." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 2. Sinc Pulses\n", + "\n", + "In practice, block pulses are rarely used for slice-selective excitation because their\n", + "frequency profile (the Fourier transform of the pulse shape) is a sinc function, which has\n", + "poor slice selectivity.\n", + "\n", + "Instead, we use **sinc-shaped** RF pulses, whose Fourier transform approximates a\n", + "rectangular frequency profile — exactly what we need for sharp slice selection.\n", + "\n", + "The key parameters of a sinc pulse are:\n", + "\n", + "- **`duration`** — the total pulse duration\n", + "- **`time_bw_product`** (TBW) — the time-bandwidth product, which controls the number of\n", + " zero-crossings and thus the sharpness of the slice profile. Higher TBW = sharper profile\n", + " but more sidelobes.\n", + "- **`apodization`** — a windowing factor (0 to 1) that suppresses sidelobes at the cost of\n", + " a slightly wider transition band. 0 = no apodization (pure sinc), 1 = full Hanning window." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_sinc = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " apodization=0.5,\n", + " system=system,\n", + " use='excitation',\n", + ")\n", + "\n", + "print(f'Peak amplitude: {np.max(np.abs(rf_sinc.signal)):.1f} Hz')\n", + "print(f'Shape duration: {rf_sinc.shape_dur * 1e3:.2f} ms')\n", + "print(f'Total duration: {pp.calc_duration(rf_sinc) * 1e3:.2f} ms')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_rf(rf_sinc, title='90° Sinc Pulse (3 ms, TBW=4, apod=0.5)')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.1 Effect of time-bandwidth product\n", + "\n", + "The time-bandwidth product (TBW) determines the number of lobes in the sinc pulse. A higher\n", + "TBW results in more zero-crossings, which leads to a sharper slice profile but also requires\n", + "a longer pulse duration (or higher gradient amplitude) to maintain the same slice thickness." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tbw_values = [2, 4, 8]\n", + "\n", + "fig, axes = plt.subplots(1, len(tbw_values), figsize=(14, 3), sharey=True)\n", + "for ax, tbw in zip(axes, tbw_values):\n", + " rf = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=tbw,\n", + " apodization=0.0,\n", + " system=system,\n", + " )\n", + " ax.plot(rf.t * 1e3, np.abs(rf.signal), 'b-', linewidth=1.5)\n", + " ax.set_title(f'TBW = {tbw}')\n", + " ax.set_xlabel('Time (ms)')\n", + " ax.grid(True, alpha=0.3)\n", + "\n", + "axes[0].set_ylabel('Amplitude (Hz)')\n", + "fig.suptitle('Sinc Pulses: Effect of Time-Bandwidth Product (no apodization)', y=1.02)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.2 Effect of apodization\n", + "\n", + "Apodization applies a Hanning window to the sinc pulse. This suppresses the sidelobes,\n", + "reducing ringing artifacts in the slice profile at the cost of a slightly wider transition\n", + "band. The `apodization` parameter ranges from 0 (no windowing) to 1 (full Hanning window).\n", + "A value of 0.5 is commonly used as a good compromise." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "apodization_values = [0.0, 0.25, 0.5, 0.75]\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(8, 3))\n", + "for apod in apodization_values:\n", + " rf = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " apodization=apod,\n", + " system=system,\n", + " )\n", + " ax.plot(rf.t * 1e3, np.abs(rf.signal), label=f'apodization = {apod}', linewidth=1.5)\n", + "\n", + "ax.set_xlabel('Time (ms)')\n", + "ax.set_ylabel('Amplitude (Hz)')\n", + "ax.set_title('Sinc Pulses: Effect of Apodization (TBW=4)')\n", + "ax.legend()\n", + "ax.grid(True, alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice how increasing the apodization smoothly tapers the sidelobes towards zero, resulting\n", + "in a smoother pulse envelope." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 3. Gaussian Pulses\n", + "\n", + "**Gaussian pulses** have a bell-shaped envelope. They are smooth and have no sidelobes,\n", + "which makes them useful in situations where a smooth frequency response is more important\n", + "than a sharp slice profile (e.g., fat saturation, magnetization preparation).\n", + "\n", + "The `pp.make_gauss_pulse` function has a very similar interface to `pp.make_sinc_pulse`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_gauss = pp.make_gauss_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " system=system,\n", + " use='excitation',\n", + ")\n", + "\n", + "print(f'Peak amplitude: {np.max(np.abs(rf_gauss.signal)):.1f} Hz')\n", + "print(f'Shape duration: {rf_gauss.shape_dur * 1e3:.2f} ms')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_rf(rf_gauss, title='90° Gaussian Pulse (3 ms, TBW=4)')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3.1 Comparing pulse shapes\n", + "\n", + "Let's compare all three pulse shapes side by side with the same duration and flip angle." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "duration = 3e-3\n", + "flip_angle = np.deg2rad(90)\n", + "\n", + "rf_block_cmp = pp.make_block_pulse(flip_angle=flip_angle, delay=system.rf_dead_time, duration=duration, system=system)\n", + "rf_sinc_cmp = pp.make_sinc_pulse(\n", + " flip_angle=flip_angle,\n", + " delay=system.rf_dead_time,\n", + " duration=duration,\n", + " time_bw_product=4,\n", + " apodization=0.5,\n", + " system=system,\n", + ")\n", + "rf_gauss_cmp = pp.make_gauss_pulse(\n", + " flip_angle=flip_angle, delay=system.rf_dead_time, duration=duration, time_bw_product=4, system=system\n", + ")\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(8, 3))\n", + "ax.plot(\n", + " np.concatenate(([0], rf_block_cmp.t, [rf_block_cmp.t[-1]])) * 1e3,\n", + " np.concatenate(([0], np.abs(rf_block_cmp.signal), [0])),\n", + " label='Block',\n", + " linewidth=1.5,\n", + ")\n", + "ax.plot(rf_sinc_cmp.t * 1e3, np.abs(rf_sinc_cmp.signal), label='Sinc (TBW=4, apod=0.5)', linewidth=1.5)\n", + "ax.plot(rf_gauss_cmp.t * 1e3, np.abs(rf_gauss_cmp.signal), label='Gaussian (TBW=4)', linewidth=1.5)\n", + "ax.set_xlabel('Time (ms)')\n", + "ax.set_ylabel('Amplitude (Hz)')\n", + "ax.set_title('Comparison of RF Pulse Shapes (90°, 3 ms)')\n", + "ax.legend()\n", + "ax.grid(True, alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Each shape has different trade-offs:\n", + "\n", + "| Shape | Slice profile | Sidelobes | Typical use |\n", + "|---|---|---|---|\n", + "| Block | sinc-shaped (poor) | Many | Non-selective excitation |\n", + "| Sinc | Approximately rectangular | Controllable via TBW and apodization | Slice-selective excitation/refocusing |\n", + "| Gaussian | Gaussian (smooth) | None | Fat saturation, preparation |" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 4. Slice-Selective RF Pulses\n", + "\n", + "To excite only a specific slice of the object, we play a **slice-selection gradient** during\n", + "the RF pulse. The gradient creates a linear mapping between spatial position and resonance\n", + "frequency. Combined with the frequency-selective RF pulse, only spins within the desired\n", + "slice are affected.\n", + "\n", + "After the RF pulse, the spins within the slice have accumulated different phases due to the\n", + "slice-selection gradient. A **rephasing gradient** (also called rewinder) is applied to\n", + "refocus this phase dispersion.\n", + "\n", + "In PyPulseq, you can request the slice-selection gradient and its rephaser by setting\n", + "`return_gz=True` and providing a `slice_thickness`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_ss, gz_ss, gz_reph = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " apodization=0.5,\n", + " slice_thickness=5e-3,\n", + " return_gz=True,\n", + " system=system,\n", + " use='excitation',\n", + ")\n", + "\n", + "print(f'RF shape duration: {rf_ss.shape_dur * 1e3:.2f} ms')\n", + "print(f'RF delay (adjusted due to rf_dead_time): {rf_ss.delay * 1e6:.0f} µs')\n", + "print(f'RF total duration (including dead/ringdown time): {pp.calc_duration(rf_ss) * 1e3:.2f} ms')\n", + "print(f'Slice-select gradient delay: {gz_ss.delay * 1e6:.0f} µs')\n", + "print(f'Slice-select gradient rise time: {gz_ss.rise_time * 1e6:.0f} µs')\n", + "print(f'Slice-select gradient flat time: {gz_ss.flat_time * 1e3:.2f} ms')\n", + "print(f'Slice-select gradient total duration: {pp.calc_duration(gz_ss) * 1e3:.2f} ms')\n", + "print(f'Rephasing gradient duration: {pp.calc_duration(gz_reph) * 1e3:.2f} ms')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4.1 Visualizing the slice-selective excitation block\n", + "\n", + "Let's build a minimal sequence with just the slice-selective RF pulse and its gradients to\n", + "visualize how they are arranged in time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "seq = pp.Sequence(system)\n", + "seq.add_block(rf_ss, gz_ss)\n", + "seq.add_block(gz_reph)\n", + "\n", + "seq.plot(grad_disp='mT/m')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the plot above, you can see:\n", + "\n", + "- **Block 1**: The RF pulse (top) plays simultaneously with the slice-selection gradient on\n", + " the z-axis (bottom). The gradient has a trapezoidal shape with a flat top during the RF\n", + " pulse and ramps on either side.\n", + "- **Block 2**: The rephasing gradient on the z-axis, which has the opposite polarity to\n", + " compensate for the phase accumulated during the second half of the slice-selection gradient.\n", + "\n", + "Note how the RF pulse's and gradient's `delay` were automatically adjusted to account for both, the RF dead time and the gradient's\n", + "rise time, ensuring that the RF pulse starts when the gradient has reached its flat top." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4.2 Effect of slice thickness and gradient units\n", + "\n", + "The slice thickness is controlled by the amplitude of the slice-selection gradient: a\n", + "thinner slice requires a stronger gradient. Let's compare different slice thicknesses.\n", + "\n", + "**A note on gradient units in Pulseq:** Internally, (Py)Pulseq stores all gradient\n", + "amplitudes in **Hz/m** rather than the more familiar **mT/m**. This is actually very\n", + "convenient because the Larmor equation directly links frequency to position:\n", + "\n", + "$$f = \\gamma \\cdot G \\cdot x$$\n", + "\n", + "where $G$ is the gradient amplitude in T/m and $\\gamma$ is the gyromagnetic ratio (Hz/T).\n", + "By storing gradients in Hz/m, we can directly compute frequencies without extra conversion\n", + "factors. The relationship between the two units is:\n", + "\n", + "$$G_\\text{[Hz/m]} = \\gamma \\cdot G_\\text{[T/m]} = \\gamma \\cdot G_\\text{[mT/m]} \\cdot 10^{-3}$$\n", + "\n", + "To convert back to mT/m (e.g., for display), we divide by `system.gamma` and multiply by\n", + "1000:\n", + "\n", + "$$G_\\text{[mT/m]} = \\frac{G_\\text{[Hz/m]}}{\\gamma} \\cdot 10^3$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "slice_thicknesses_mm = [2, 5, 10]\n", + "\n", + "for st_mm in slice_thicknesses_mm:\n", + " _, gz, gzr = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " apodization=0.5,\n", + " slice_thickness=st_mm * 1e-3,\n", + " return_gz=True,\n", + " system=system,\n", + " )\n", + " grad_mTm = gz.amplitude / system.gamma * 1e3\n", + " print(\n", + " f'Slice thickness: {st_mm:2d} mm -> '\n", + " f'Gradient: {gz.amplitude:10.1f} Hz/m = {grad_mTm:6.2f} mT/m, '\n", + " f'Rephaser duration: {pp.calc_duration(gzr) * 1e3:.2f} ms'\n", + " )\n", + "\n", + "print(\n", + " f'\\nMax gradient amplitude of our system: {system.max_grad:.1f} Hz/m = {system.max_grad / system.gamma * 1e3:.1f} mT/m'\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As expected, thinner slices require stronger gradients. The rephasing gradient duration also\n", + "changes because its area must match the area of the second half of the slice-selection\n", + "gradient plus the ramp-down area.\n", + "\n", + "\n", + "### 4.3 What happens when the gradient limit is exceeded?\n", + "\n", + "If we request a very thin slice, the required gradient amplitude may exceed the hardware\n", + "limit (`max_grad`). In that case, PyPulseq raises an error. Let's try a 0.5 mm slice with\n", + "our current system limits (28 mT/m):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " _, gz_thin, _ = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " apodization=0.5,\n", + " slice_thickness=0.5e-3, # 0.5 mm — very thin!\n", + " return_gz=True,\n", + " system=system,\n", + " )\n", + "except ValueError as e:\n", + " # Calculate what the required gradient would be\n", + " bandwidth = 4 / 3e-3 # TBW / duration\n", + " required_grad = bandwidth / 0.5e-3 # Hz/m\n", + " required_grad_mTm = required_grad / system.gamma * 1e3\n", + " print(f'Error: {e}')\n", + " print(f'\\nRequired gradient amplitude: {required_grad:.1f} Hz/m = {required_grad_mTm:.1f} mT/m')\n", + " print(\n", + " f'System max gradient amplitude: {system.max_grad:.1f} Hz/m = {system.max_grad / system.gamma * 1e3:.1f} mT/m'\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To achieve thinner slices without exceeding the gradient limit, you can either increase the\n", + "pulse duration (which reduces the required bandwidth) or increase the `max_grad` in your\n", + "system limits (if the scanner hardware supports it)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4.4 Gaussian slice-selective pulse\n", + "\n", + "Slice selection also works with Gaussian pulses. The interface is identical." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_gauss_ss, gz_gauss_ss, gz_gauss_reph = pp.make_gauss_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " slice_thickness=5e-3,\n", + " return_gz=True,\n", + " system=system,\n", + " use='excitation',\n", + ")\n", + "\n", + "seq = pp.Sequence(system)\n", + "seq.add_block(rf_gauss_ss, gz_gauss_ss)\n", + "seq.add_block(gz_gauss_reph)\n", + "\n", + "seq.plot()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 5. Flip Angle, Phase Offset, and Frequency Offset\n", + "\n", + "Every RF pulse in PyPulseq can be configured with three important parameters that control\n", + "how it interacts with the magnetization:\n", + "\n", + "- **`flip_angle`** — the angle (in radians) by which the magnetization is rotated\n", + "- **`phase_offset`** — the phase of the RF pulse (in radians), which determines the axis of\n", + " rotation in the transverse plane\n", + "- **`freq_offset`** — a frequency offset (in Hz) applied to the RF pulse, which can be used\n", + " to shift the excitation to a different frequency (e.g., for multi-slice imaging or\n", + " fat/water-selective excitation)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 5.1 Phase offset\n", + "\n", + "The phase offset rotates the axis of excitation in the transverse plane. For example, a\n", + "phase offset of 0 excites along the x-axis, while π/2 excites along the y-axis. This is\n", + "commonly used in RF spoiling schemes and for refocusing pulses.\n", + "\n", + "Let's create sinc pulses with different phase offsets and visualize the complex waveform." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "phase_offsets = [0, np.pi / 4, np.pi / 2, np.pi]\n", + "phase_labels = ['0', 'π/4', 'π/2', 'π']\n", + "\n", + "fig, axes = plt.subplots(2, len(phase_offsets), figsize=(16, 4), sharex=True, sharey=True)\n", + "\n", + "for i, (phase, label) in enumerate(zip(phase_offsets, phase_labels)):\n", + " rf = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " apodization=0.5,\n", + " phase_offset=phase,\n", + " system=system,\n", + " )\n", + " t_ms = rf.t * 1e3\n", + " # Apply the phase offset to the signal for visualization.\n", + " # PyPulseq stores phase_offset separately (applied at playback), so the\n", + " # signal array itself is always real-valued. We rotate it here to show\n", + " # the effect of the phase offset on the complex waveform.\n", + " signal_rotated = rf.signal * np.exp(1j * rf.phase_offset)\n", + " axes[0, i].plot(t_ms, np.real(signal_rotated), 'b-', linewidth=1.5)\n", + " axes[0, i].set_title(f'Phase = {label}')\n", + " axes[0, i].grid(True, alpha=0.3)\n", + " axes[1, i].plot(t_ms, np.imag(signal_rotated), 'r-', linewidth=1.5)\n", + " axes[1, i].set_xlabel('Time (ms)')\n", + " axes[1, i].grid(True, alpha=0.3)\n", + "\n", + "axes[0, 0].set_ylabel('Real part (Hz)')\n", + "axes[1, 0].set_ylabel('Imaginary part (Hz)')\n", + "fig.suptitle('Sinc Pulses with Different Phase Offsets', y=1.02)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The phase offset rotates the complex RF waveform. At phase = 0, the signal is purely real.\n", + "At phase = π/2, it becomes purely imaginary. At intermediate values, both real and imaginary\n", + "components are present." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 5.2 Frequency offset\n", + "\n", + "The frequency offset shifts the center frequency of the RF pulse. This is useful for:\n", + "\n", + "- **Multi-slice imaging**: exciting different slices by offsetting the RF frequency while\n", + " keeping the same slice-selection gradient\n", + "- **Fat/water-selective excitation**: targeting specific spectral components\n", + "\n", + "The frequency offset is stored as a property of the RF event and is applied during playback\n", + "on the scanner. It does **not** change the waveform shape itself.\n", + "\n", + "Let's demonstrate multi-slice excitation by creating RF pulses with different frequency\n", + "offsets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "slice_positions_mm = [-10, -5, 0, 5, 10] # 5 slices, 5 mm apart\n", + "\n", + "# Create the slice-selective pulse (we need the gradient amplitude to calculate freq offsets)\n", + "rf_base, gz_base, gz_reph_base = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " apodization=0.5,\n", + " slice_thickness=5e-3,\n", + " return_gz=True,\n", + " system=system,\n", + " use='excitation',\n", + ")\n", + "\n", + "# The frequency offset for a given slice position is:\n", + "# freq_offset = gradient_amplitude * slice_position\n", + "for pos_mm in slice_positions_mm:\n", + " freq_offset = gz_base.amplitude * pos_mm * 1e-3 # Convert mm to m\n", + " print(f'Slice at {pos_mm:+3d} mm -> freq_offset = {freq_offset:+8.1f} Hz')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 6. The `use` Parameter\n", + "\n", + "Every RF pulse in (Py)Pulseq has a `use` parameter that describes its purpose in the sequence.\n", + "Since Pulseq file format version 1.5, this metadata is stored in the `.seq` file and is thus still available \n", + "when reading the sequence. \n", + "\n", + "The supported values are:\n", + "\n", + "| Value | Description |\n", + "|---|---|\n", + "| `'excitation'` | Excitation pulse (tips magnetization into the transverse plane) |\n", + "| `'refocusing'` | Refocusing pulse (e.g., 180° pulse in spin echo sequences) |\n", + "| `'inversion'` | Inversion pulse (e.g., 180° pulse for inversion recovery) |\n", + "| `'saturation'` | Saturation pulse (e.g., fat saturation, CEST) |\n", + "| `'preparation'` | Preparation pulse (e.g., T2 preparation) |\n", + "| `'other'` | Other purpose |\n", + "| `'undefined'` | Default if not specified |\n", + "\n", + "Let's create examples of the most common use cases." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 6.1 Excitation pulse\n", + "\n", + "An excitation pulse tips the longitudinal magnetization into the transverse plane. It\n", + "typically has a flip angle between 1° and 90°." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_excitation, gz_excitation, gz_excitation_reph = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " apodization=0.5,\n", + " slice_thickness=5e-3,\n", + " return_gz=True,\n", + " system=system,\n", + " use='excitation',\n", + ")\n", + "\n", + "print(f'Excitation pulse: flip angle = 90°, use = {rf_excitation.use}')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 6.2 Refocusing pulse\n", + "\n", + "A refocusing pulse is a 180° pulse used in spin echo sequences to refocus the transverse\n", + "magnetization. It is typically played with a phase offset of π/2 relative to the excitation\n", + "pulse (i.e., along the y-axis if the excitation is along the x-axis)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_refocusing, gz_refocusing, gz_refocusing_reph = pp.make_sinc_pulse(\n", + " flip_angle=np.deg2rad(180),\n", + " delay=system.rf_dead_time,\n", + " duration=3e-3,\n", + " time_bw_product=4,\n", + " apodization=0.5,\n", + " phase_offset=np.pi / 2,\n", + " slice_thickness=5e-3,\n", + " return_gz=True,\n", + " system=system,\n", + " use='refocusing',\n", + ")\n", + "\n", + "print(f'Refocusing pulse: flip angle = 180°, phase = π/2, use = {rf_refocusing.use}')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 6.3 Inversion pulse\n", + "\n", + "An inversion pulse flips the longitudinal magnetization by 180°. Unlike a refocusing pulse,\n", + "it acts on the longitudinal component and is used in inversion recovery sequences (e.g.,\n", + "MPRAGE, STIR, FLAIR).\n", + "\n", + "For inversion, adiabatic pulses are often preferred because they achieve a uniform 180° flip\n", + "even in the presence of B1 inhomogeneities. PyPulseq provides `pp.make_adiabatic_pulse` for\n", + "this purpose, supporting hyperbolic secant (`hypsec`) and WURST pulse types." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_inversion = pp.make_adiabatic_pulse(\n", + " pulse_type='hypsec',\n", + " delay=system.rf_dead_time,\n", + " duration=10e-3,\n", + " system=system,\n", + " use='inversion',\n", + ")\n", + "\n", + "print(f'Inversion pulse type: hypsec, use = {rf_inversion.use}')\n", + "print(f'Duration: {rf_inversion.shape_dur * 1e3:.1f} ms')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_rf(rf_inversion, title='Adiabatic Inversion Pulse (Hyperbolic Secant, 10 ms)')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Unlike the previous pulses, the adiabatic pulse has a **non-zero phase** that varies across\n", + "the pulse duration. This frequency sweep is what makes adiabatic pulses robust to B1\n", + "inhomogeneities. For more details about adiabatic pulses in general, see for example [https://mriquestions.com/adiabatic-excitation.html](https://mriquestions.com/adiabatic-excitation.html)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 6.4 Saturation pulse\n", + "\n", + "A saturation pulse destroys the magnetization of a specific spectral component (e.g., fat).\n", + "It typically has a flip angle around 90° followed by a spoiler gradient. For frequency-selective\n", + "saturation (e.g., fat sat), a Gaussian pulse with a specific frequency offset is commonly\n", + "used. For (Py)Pulseq versions >= 1.5.0, the frequency offset can be specified in ppm (`freq_ppm`)\n", + "to make the RF pulse field strength independent." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Fat saturation pulse: 90° Gaussian pulse at the fat frequency offset (~-440 Hz at 3T)\n", + "rf_fatsat = pp.make_gauss_pulse(\n", + " flip_angle=np.deg2rad(90),\n", + " delay=system.rf_dead_time,\n", + " duration=8e-3,\n", + " time_bw_product=2,\n", + " freq_ppm=-3.45,\n", + " system=system,\n", + " use='saturation',\n", + ")\n", + "\n", + "print(f'Fat saturation pulse: use = {rf_fatsat.use}, freq_offset = {rf_fatsat.freq_ppm} ppm')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## Summary\n", + "\n", + "In this tutorial, you learned about the different types of RF pulses available in PyPulseq:\n", + "\n", + "| Concept | Function |\n", + "|---|---|\n", + "| Block (rectangular) pulse | `pp.make_block_pulse(flip_angle, duration, ...)` |\n", + "| Sinc pulse | `pp.make_sinc_pulse(flip_angle, duration, time_bw_product, apodization, ...)` |\n", + "| Gaussian pulse | `pp.make_gauss_pulse(flip_angle, duration, time_bw_product, ...)` |\n", + "| Adiabatic pulse | `pp.make_adiabatic_pulse(pulse_type, duration, ...)` |\n", + "| Slice selection | Set `return_gz=True` and `slice_thickness=...` |\n", + "| Phase offset | Set `phase_offset=...` (in radians) |\n", + "| Frequency offset | Set `freq_offset=...` (in Hz) |\n", + "| Pulse purpose | Set `use='excitation'`, `'refocusing'`, `'inversion'`, `'saturation'`, ... |\n", + "\n", + "### Key takeaways\n", + "\n", + "- **Block pulses** are simple but non-selective. Use them for hard-pulse excitation.\n", + "- **Sinc pulses** provide good slice selectivity. Tune `time_bw_product` and `apodization`\n", + " for the desired profile sharpness.\n", + "- **Gaussian pulses** are smooth and sidelobe-free. Use them for spectral saturation or\n", + " preparation.\n", + "- **Adiabatic pulses** are robust to B1 inhomogeneities. Use them for inversion.\n", + "- Set `return_gz=True` with a `slice_thickness` to automatically generate the slice-selection\n", + " gradient and its rephaser.\n", + "- Use `phase_offset` and `freq_offset` to control the excitation axis and frequency.\n", + "\n", + "## Next steps\n", + "\n", + "In the next tutorial, we'll explore **gradient events** in detail — including trapezoids,\n", + "arbitrary gradients, and how to use them for spatial encoding." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "pypulseq", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}