diff --git a/.github/workflows/deployment.yml b/.github/workflows/deployment.yml index ed460e6..6ea41a1 100644 --- a/.github/workflows/deployment.yml +++ b/.github/workflows/deployment.yml @@ -66,7 +66,7 @@ jobs: python -m build - name: Store PyPI Distribution - uses: actions/upload-artifact@v6 + uses: actions/upload-artifact@v7 with: name: pypi-package-distribution path: dist/ @@ -88,7 +88,7 @@ jobs: steps: - name: Download PyPI Distribution - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: name: pypi-package-distribution path: dist/ diff --git a/.github/workflows/examples.yml b/.github/workflows/examples.yml index 28105cf..9efafdf 100644 --- a/.github/workflows/examples.yml +++ b/.github/workflows/examples.yml @@ -33,7 +33,7 @@ jobs: # Upload the book's HTML as an artifact (useful for PR previews) - name: Upload built HTML - uses: actions/upload-artifact@v6 + uses: actions/upload-artifact@v7 with: name: built-book path: examples/_build/html @@ -48,7 +48,7 @@ jobs: steps: - name: Download built HTML - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: name: built-book path: examples/_build/html diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 40d8e61..0e96dfe 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -35,7 +35,7 @@ jobs: - name: Upload PyTest Coverage Report if: always() - uses: actions/upload-artifact@v6 + uses: actions/upload-artifact@v7 with: name: pytest-report-${{ matrix.python-version }} path: | diff --git a/.github/workflows/report_pytest.yml b/.github/workflows/report_pytest.yml index 11f00d2..5342767 100644 --- a/.github/workflows/report_pytest.yml +++ b/.github/workflows/report_pytest.yml @@ -21,7 +21,7 @@ jobs: steps: # use the results of python 3.12, consider this as target platform - name: Download PyTest report artifact for Python 3.12 - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: name: pytest-report-3.12 run-id: ${{ github.event.workflow_run.id }} @@ -49,7 +49,7 @@ jobs: - name: Post PyTest Coverage Comment id: coverage_comment - uses: MishaKav/pytest-coverage-comment@v1.4.0 + uses: MishaKav/pytest-coverage-comment@v1.6.0 with: issue-number: ${{ steps.pr-context.outputs.number }} pytest-coverage-path: pytest-coverage.txt @@ -82,7 +82,7 @@ jobs: actions: read # Required for downloading artifacts steps: - name: Download artifacts - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: run-id: ${{ github.event.workflow_run.id }} github-token: ${{ secrets.GITHUB_TOKEN }} @@ -90,7 +90,7 @@ jobs: - name: Post PyTest Coverage Comment on push id: coverage_comment - uses: MishaKav/pytest-coverage-comment@v1.4.0 + uses: MishaKav/pytest-coverage-comment@v1.6.0 with: pytest-coverage-path: pytest-coverage.txt junitxml-path: pytest.xml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 5b67024..fec4759 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -15,7 +15,7 @@ repos: - id: mixed-line-ending - repo: https://github.com/kynan/nbstripout - rev: 0.9.0 + rev: 0.9.1 hooks: # cleans the .ipynbs (removes outputs, resets all cell-ids to 0..N, cleans steps) # also clean any kernel information left after execution @@ -25,14 +25,14 @@ repos: files: examples/.*\.ipynb$ - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.14.14 + rev: v0.15.4 hooks: - id: ruff # linter args: [--fix] - id: ruff-format # formatter - repo: https://github.com/crate-ci/typos - rev: v1.43.0 + rev: v1.44.0 hooks: - id: typos diff --git a/TriangleExample/duyn_triangle.py b/TriangleExample/duyn_triangle.py new file mode 100644 index 0000000..7b46913 --- /dev/null +++ b/TriangleExample/duyn_triangle.py @@ -0,0 +1,227 @@ +import os +import time + +import numpy as np +import pypulseq as pp +import yaml + +# This file generates a sequence for the Duyn method. It also incorporates the MFC, so that one can acquire the same +# measurement from both the MFC and the Scanner (provided the MFC is set up) + +# Main folder to keep everything in +main_dir = 'TriangleExample/' + +# Main YAML +yaml_dir = 'TriangleExample/triangle_parameters.yaml' + +# Add a custom name to the file to differentiate between measurements +add = 'Versuch_2' + +# Every measurement with a new name gets saved as in a folder with a date stamp +output_path = main_dir + f'/{time.strftime("%Y%m%d")}_' + str(add) +print('output path: ', output_path) + +# Check and make if directory doesn't exist +if os.path.exists(output_path) == False: + os.mkdir(f'{output_path}') + + +# Get Parameters from YAML +with open(yaml_dir) as stream: + try: + # print(yaml.safe_load(stream)) + parameters = yaml.safe_load(stream) + except yaml.YAMLError as exc: + print(exc) + + +# Step 1 - Define our sequence parameters according to our scanner +# Create PyPulseq Sequence Object and Set System Limits +seq = pp.Sequence() +system = pp.Opts( + max_grad=parameters['max_grad'] * 1e3, + grad_unit='mT/m', + max_slew=parameters['slew_rate'], + slew_unit='T/m/s', + rf_ringdown_time=parameters['rf_ringdown_time'], + rf_dead_time=parameters['rf_dead_time'], + adc_dead_time=parameters['adc_dead_time'], +) + +method = ['Duyn'] + +# Quick Dynamics Calculation (necessary for the Field Cam but also helpful to see how many measurement instances) +CamNrDynamics = ( + len(parameters['enumerate_coeff']) + * len(parameters['SLICE_POS']) + * len(parameters['GRAD']) + * len(parameters['RISE_TIMES']) + * parameters['N_AVG'] +) +print(f'CameraNrDynamics: {CamNrDynamics}') + + +# Create and Set Sinc Pulse Parameters +rf, g_sl, g_sl_r = pp.make_sinc_pulse( + flip_angle=parameters['RF_PHI'] / 180 * np.pi, + delay=system.rf_dead_time, + duration=parameters['RF_DUR'], + slice_thickness=parameters['SLICE_THICKNESS'], + apodization=parameters['apodization'], + time_bw_product=parameters['RF_BWT_PROD'], + system=system, + return_gz=True, +) + + +# Necessary delays such as TR, EddyCurrentComp, Trig for MFC +grad_free_time = pp.make_delay(parameters['grad_free']) # from skope_gtf.m +delay_tr = pp.make_delay(parameters['TR']) # TR delay +delay_ec = pp.make_delay(parameters['ec_delay']) # eddy current compensation delay + + +# Step 6 - Actually making the measurement method take in data, so trig for MFC and ADC for Scanner (our beloved) +trig = pp.make_digital_output_pulse( + channel='ext1', duration=parameters['trig_duration'], delay=0 +) # creating the MFC trigger pulse + +# Automatic ADC samples and duration +adc_num_samples = int(np.round(parameters['adc_duration'] / parameters['dwell_time'])) +print('adc_num: ', adc_num_samples) + +# ADC Block +adc = pp.make_adc( + num_samples=adc_num_samples, duration=parameters['adc_duration'], system=system, delay=system.adc_dead_time +) # Analog Digital Converter + + +# Step 7 - Creating the loop + +# Starting with averages +for avg in range(parameters['N_AVG']): + avg_str = 'N_AVG/' # Labels for loop order + avg_label = pp.make_label(type='SET', label='AVG', value=avg) + + # We enumerate over rise times to increase the amplitude of the encode gradient + for idx_rise_time, rise_time in enumerate(parameters['RISE_TIMES']): + rise_times_str = 'RiseTimes/' + seg_label = pp.make_label(type='SET', label='SEG', value=idx_rise_time) + rise_time = np.round( + rise_time, decimals=6 + ) # Not rounding the rise time here causes an error during amplitude calculation + + # Create ADC labels with label SEG + for idx_grad, grad in enumerate(parameters['GRAD']): + grad_str = 'Grad/' + grad_label = pp.make_label(type='SET', label='SET', value=idx_grad) + + # Make labels of each axis with the index of the gradients + for idx_slice_pos, slice_pos in enumerate(parameters['SLICE_POS']): + slice_pos_str = 'SlicePos/' + # Changing SLICE_POS to differentiate between Duyn and Rahmer (this is irrelevant for Duyn method) + slice_label = pp.make_label(type='SET', label='SLC', value=idx_slice_pos) + + # amp_fac goes between -1 and 0 for Duyn + for idx_amp_fac, amp_fac in enumerate(parameters['enumerate_coeff']): + enco_str = 'EnCo/' + amp_fac_label = pp.make_label(type='SET', label='REP', value=idx_amp_fac) + + amp = amp_fac * system.max_slew * (rise_time) + # Starting amp factor * how fast amp rises * how long it rises for = amplitude value at the end of rise + + # Add RF Pulse with Slice Selection + rf.freq_offset = g_sl.amplitude * slice_pos + # grad_slice_select amplitude * slice pos => + # Gr(t)*Dr = freq difference of the rf pulse + g_sl.channel = grad + g_sl_r.channel = grad + seq.add_block(rf, g_sl) # RF Pulse and the Slice Select Grad together + seq.add_block(g_sl_r) # Slice Select Rewinder + + # List of labels for everything + label_contents = [ + avg_label, + seg_label, + grad_label, + slice_label, + amp_fac_label, + ] # All labels, average, segment, gradient, slice and amp_factor + + # Add Eddy Current Compensation Delay + + seq.add_block(delay_ec) + + # Add ADC and Triangle Gradient + if amp_fac != 0: # If there is an amp (if gradient present) + # Create Triangle Gradient + + g_triangle = pp.make_trapezoid( + channel=grad, system=system, rise_time=rise_time, flat_time=0, amplitude=amp, delay=0 + ) + # trapezoid with no flat time = triangle + # enumerating over rise times gives differently sized triangles + # Watch out for delay in the recon + + seq.add_block(trig) # MFC trigger + # starts CameraAcqDuration, not relevant if no MFC present + + seq.add_block(grad_free_time) # Delay between MFC trig and ADC + + seq.add_block(adc, g_triangle) # ADC and gradients during the same block + + else: + seq.add_block(trig) # MFC trigger + + seq.add_block(grad_free_time) # Delay between MFC trig and ADC + seq.add_block(adc) # ADC without the gradients + + # Add TR Delay + seq.add_block(delay_tr) + + seq.add_block(avg_label) + +# Total loop order string for recon +label_str = avg_str + rise_times_str + grad_str + slice_pos_str + enco_str + +# Step 8 - Timing Check (Also good for sanity) +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) + + +# Step 9 - Setting Definitions for the seq file (to be used by the scanner and the MFC) +seq.set_definition('rise_time', parameters['RISE_TIMES']) +seq.set_definition('grad', parameters['GRAD']) +seq.set_definition('slice_pos', parameters['SLICE_POS']) +seq.set_definition('avg', parameters['N_AVG']) +seq.set_definition('tr_delay', parameters['TR']) +seq.set_definition('Acquisition Method: ', method) +# +seq.set_definition('CameraNrDynamics', CamNrDynamics) +# CameraNrDynamics: Maximum number of acquisitions performed by the Camera Acquisition System. +seq.set_definition('CameraNrSyncDynamics', parameters['cam_nr_sync_dyn']) +seq.set_definition('CameraAcqDuration', parameters['cam_acq_duration']) +seq.set_definition('CameraInterleaveTR', parameters['cam_interleave_tr']) +seq.set_definition('CameraAqDelay', parameters['cam_acq_delay']) +# I find it helpful to print out some important parameters +# Especially ones you need to differentiate between measurements + +# Step 10 - Save File (and maybe even plot) +filename = f'{output_path}/{time.strftime("%Y%m%d")}_Triangles' + + +print(f"\nSaving sequence file '{filename}.seq' in 'output' folder.") +seq.write(str(filename), create_signature=True) + +# seq.plot(label="avg", save=True) +seq.plot(save=False, grad_disp='mT/m') # Plot the sequence and display the gradients in mT/m + +# Write the loop order into the YAML file to be stored +with open(f'{output_path}/{add}.yaml', 'w+') as f: + yaml.dump(parameters, f, sort_keys=False) + +# In the end, you will have specifically named folders containing a sequence and the YAML parameters used +# to create that exact sequence, so that there is no confusion during reconstruction diff --git a/TriangleExample/triangle_parameters.yaml b/TriangleExample/triangle_parameters.yaml new file mode 100644 index 0000000..3ad27f7 --- /dev/null +++ b/TriangleExample/triangle_parameters.yaml @@ -0,0 +1,56 @@ +--- +# Scanner Parameters: +# Gradient Raster Time: +dwell_time : 0.00001 +# Maximum Slew Rate: +slew_rate : 120 # in [T/m/s] +# Maximum Gradient Amplitude in [T/m] +max_grad : 0.03 +# +# Sequence Parameters: +GAMMA : 267515315.1 #e8 +# Number of averages: +N_AVG : 3 +# Duration of the Scanner ADC in [sec] +adc_duration : 0.06 +SLICE_THICKNESS : 0.0015 # in [m] +SLICE_POS : [0.04] # in [m] +# For the "Duyn method": the sign before the gradient +# [-1,0] means the gradient is firstly played inverse, then the same slice is selected but without a gradient +enumerate_coeff : [-1, 0] +TR : 1 # Repetition time in [sec] +ec_delay : 0.001 # Eddy Current delay in [sec] +RISE_TIMES : [0.00005, 0.00006, 0.00007, 0.00008, 0.00009, 0.0001, 0.00011, 0.00012, 0.00013, 0.00014, 0.00015, 0.00016] # Rise times for the different triangles, in [sec] +# Order of axis: +# Change the order to change axis order, can even do xyy to double y's and leave out z +GRAD : ['x', 'y', 'z'] +# +# SINC Parameters +RF_PHI : 90 # Flip angle +RF_DUR : 0.0084 # Duration of the pulse T, in [sec] 8.4e-3 +RF_BWT_PROD : 4 # T*f, "quality" of slice profile, +apodization : 0.5 # Constant for apodization +# +# System Parameters +rf_ringdown_time : 0.00003 #30e-6 +rf_dead_time : 0.0001 #100e-6 +adc_dead_time : 0.00001 # By how much to pad the ADC duration with, in [sec] +# + + + + + + + + + +# No need to touch the MFC parameters, they only contribute a small delay in the seq (grad_free) +# MFC Parameters +grad_free : 0.0005 # Duration of gradient-free period after trigger signal in [sec] +trig_duration : 0.00001 # Duration of trigger RF signal in [sec] +cam_acq_duration : 0.07 # How long the MFC will acquire data for after the trigger in [sec] +cam_interleave_tr : 0.4 # "Interleave" duration, unimportant but don't touch +cam_acq_delay : 0 # Acquisition delay of the MFC in [sec] +cam_nr_sync_dyn : 0 # no idea (for now) +# \ No newline at end of file diff --git a/examples/_toc.yml b/examples/_toc.yml index 9f362f2..4b4992e 100644 --- a/examples/_toc.yml +++ b/examples/_toc.yml @@ -5,6 +5,7 @@ chapters: sections: - file: radial_flash - file: spiral_flash + - file: grpe_flash_dixon - file: relaxometry.md sections: - file: t1_inv_rec_gre_single_line diff --git a/examples/girf_dyn_triangle.ipynb b/examples/girf_dyn_triangle.ipynb new file mode 100644 index 0000000..64730b1 --- /dev/null +++ b/examples/girf_dyn_triangle.ipynb @@ -0,0 +1,499 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# GIRF - Dyn-method with triangular pulses\n", + "\n", + "Estimate the gradient impulse response function using the Dyn-method with triangular pulses" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import tempfile\n", + "from collections.abc import Sequence\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 mrpro.data import KData\n", + "from mrpro.data.traj_calculators import KTrajectoryCartesian\n", + "\n", + "from mrseq.scripts.girf_dyn_triangle 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", + "tmp = tempfile.TemporaryDirectory()\n", + "fname_mrd = Path(tmp.name) / 'girf_dyn_triangle.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)\n", + "phantom.T1[:] = 400e-3\n", + "phantom.T2[:] = 40e-3\n", + "phantom.B1[:] = 1.0\n", + "phantom.B0[:] = 0.0" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "### Create the GIRF Dyn sequence\n", + "\n", + "To create the GIRF Dyn sequence, we use the previously imported [girf_dyn_triangle script](../src/mrseq/scripts/girf_dyn_triangle.py)." + ] + }, + { + "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", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### Simulate the sequence\n", + "\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": [ + "### Estimate GIRF" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "#### Sort data and compress coils" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "kdata = KData.from_file(fname_mrd, trajectory=KTrajectoryCartesian())\n", + "idx = np.lexsort(\n", + " (\n", + " kdata.header.acq_info.idx.phase.squeeze(),\n", + " kdata.header.acq_info.idx.repetition.squeeze(),\n", + " kdata.header.acq_info.idx.average.squeeze(),\n", + " )\n", + ")\n", + "kdata_sorted = kdata[torch.as_tensor(idx)]\n", + "kdata_sorted = kdata_sorted.rearrange(\n", + " '(avg rep ph) ... -> avg rep ph ...',\n", + " rep=kdata.header.acq_info.idx.repetition.max() + 1,\n", + " avg=kdata.header.acq_info.idx.average.max() + 1,\n", + " ph=kdata.header.acq_info.idx.phase.max() + 1,\n", + ")\n", + "kdat_single_coil = kdata_sorted.compress_coils(n_compressed_coils=1).data.squeeze()" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "#### Get information about the scan from the seq file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "adc_duration = sequence.get_definition('AdcDuration')\n", + "dwell_time = sequence.get_definition('DwellTime')\n", + "slice_pos = sequence.get_definition('SlicePos')[0]\n", + "gamma = sequence.system.gamma * 2 * torch.pi\n", + "rise_times = sequence.get_definition('RiseTimes')\n", + "slew_rate = sequence.get_definition('SlewRate') / sequence.system.gamma\n", + "g_delay = sequence.get_definition('AdcDelay')\n", + "g_amplitude_coeff = sequence.get_definition('GradAmplitudeCoeff')\n", + "\n", + "output_delay = 2" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "#### Unwrapped measured gradient shapes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "def unwrap_phase_difference(data: torch.Tensor) -> torch.Tensor:\n", + " \"\"\"Return doubly-unwrapped phase of polarity[0] - polarity[1].\n", + "\n", + " Input shape: ``(n_avg, 2, n_axes, n_rise, n_samples)``\n", + " Output shape: ``(n_avg, n_axes, n_rise, n_samples)``\n", + " \"\"\"\n", + " diff = data[:, 0, ...].angle() - data[:, 1, ...].angle()\n", + " return torch.from_numpy(np.unwrap(np.unwrap(diff.numpy())))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "def phase_to_gradient(\n", + " phase_mean: torch.Tensor,\n", + " phase_std: torch.Tensor,\n", + " slice_pos: float,\n", + " gamma: float,\n", + " dwell_time: float,\n", + " output_delay: int,\n", + ") -> tuple[torch.Tensor, torch.Tensor]:\n", + " \"\"\"Convert unwrapped-phase arrays to gradient waveforms via finite difference.\"\"\"\n", + " scale = slice_pos * gamma * dwell_time\n", + " grad_mean = torch.roll(phase_mean.diff(dim=-1) / scale, output_delay, dims=-1)\n", + " grad_std = phase_std.diff(dim=-1) / scale\n", + " return grad_mean, grad_std" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "kdata_sorted.data.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.plot(kdata_sorted.data[0, 1, 0, 0, 0].T.abs())\n", + "\n", + "plt.figure()\n", + "plt.plot(kdata_sorted.data[0, 1, 0, 0, 0].T.angle())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "n_samples = round(adc_duration / dwell_time)\n", + "time = torch.linspace(0.0, adc_duration, n_samples)\n", + "\n", + "phase = unwrap_phase_difference(kdat_single_coil)\n", + "phase_mean = phase.mean(dim=0)\n", + "phase_std = phase.std(dim=0)\n", + "\n", + "grad_output_mean, grad_output_std = phase_to_gradient(phase_mean, phase_std, slice_pos, gamma, dwell_time, output_delay)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot unwrapped phase (k-space trajectories) for each gradient axis\n", + "grad_axes = ['x', 'y', 'z']\n", + "fig, axes = plt.subplots(len(grad_axes), 1, sharex=True) # (avg, ax, rise, samp)\n", + "for ax_idx, (ax, label) in enumerate(zip(axes, grad_axes, strict=True)):\n", + " for rise_idx in range(phase.shape[2]):\n", + " ax.plot(time, phase_mean.numpy()[ax_idx, rise_idx, :])\n", + " ax.set_title(f'K-space — {label}')\n", + " ax.set_xlabel('Time [s]')\n", + " ax.grid()\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "#### Create ideal gradient triangles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "def build_input_triangles(\n", + " rise_times: Sequence[float],\n", + " slew_rate: float,\n", + " g_delay: float,\n", + " enumerate_coeff: Sequence[float],\n", + " dwell_time: float,\n", + " n_samples: int,\n", + ") -> torch.Tensor:\n", + " \"\"\"Build ideal triangular input waveforms.\n", + "\n", + " Returns\n", + " -------\n", + " Tensor of shape ``(n_rise, n_samples)`` [T/m].\n", + " \"\"\"\n", + " triangles = torch.zeros(len(rise_times), n_samples)\n", + " for i, rise_time in enumerate(rise_times):\n", + " n_rise = round(rise_time / dwell_time)\n", + " n_pre = round(g_delay / dwell_time)\n", + " amplitude = enumerate_coeff[0] * slew_rate * rise_time\n", + "\n", + " slope_up = torch.linspace(0.0, amplitude, n_rise + 1)\n", + " slope_down = torch.linspace(amplitude, 0.0, n_rise + 1)\n", + "\n", + " triangles[i, n_pre : n_pre + n_rise + 1] = slope_up\n", + " triangles[i, n_pre + n_rise + 1 : n_pre + 2 * n_rise + 1] = slope_down[1:]\n", + "\n", + " return triangles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "grad_input = build_input_triangles(rise_times, slew_rate, g_delay, g_amplitude_coeff, dwell_time, n_samples)" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "#### Compare ideal and measured gradient triangles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27", + "metadata": {}, + "outputs": [], + "source": [ + "# Overlay ideal input triangles against measured X-axis output\n", + "plt.figure()\n", + "for i in range(grad_input.shape[0]):\n", + " line = plt.plot(time * 1e3, grad_input[i, :], '-')\n", + " plt.plot(time[:-1] * 1e3, grad_output_mean[0, i, :], color=line[0].get_color(), linestyle='dashed')\n", + "\n", + "plt.xlim((0.0, 0.4))\n", + "plt.xlabel('Time [ms]')\n", + "plt.ylabel('Gradient Strength [T/m]')\n", + "plt.title('Input (solid) vs Measurement (dashed)')\n", + "plt.grid()\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "28", + "metadata": {}, + "source": [ + "#### Calculate GIRF and plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "def estimate_gmtf(\n", + " grad_input: torch.Tensor,\n", + " grad_output_mean: torch.Tensor,\n", + ") -> torch.Tensor:\n", + " \"\"\"Least-squares GMTF estimate over all rise times.\n", + "\n", + " Parameters\n", + " ----------\n", + " grad_input:\n", + " Ideal waveforms, shape ``(n_rise, n_samples)``.\n", + " grad_output_mean:\n", + " Measured waveforms, shape ``(n_axes, n_rise, n_samples-1)``.\n", + "\n", + " Returns\n", + " -------\n", + " Complex tensor of shape ``(n_axes, n_freq)``.\n", + " \"\"\"\n", + " n_fft_in = 2 * (grad_input.shape[-1] - 1)\n", + " n_fft_out = 2 * grad_output_mean.shape[-1]\n", + "\n", + " in_spec = torch.fft.fftshift(torch.fft.fft(grad_input, n=n_fft_in, dim=-1), dim=-1)\n", + " out_spec = torch.fft.fftshift(torch.fft.fft(grad_output_mean, n=n_fft_out, dim=-1), dim=-1)\n", + "\n", + " # Least-squares: sum_rt conj(H_in) * H_out / sum_rt |H_in|^2\n", + " numerator = (in_spec.conj() * out_spec).sum(dim=-2) # sum over rise dim\n", + " denominator = (in_spec.abs() ** 2).sum(dim=-2)\n", + " return numerator / denominator" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30", + "metadata": {}, + "outputs": [], + "source": [ + "gmtf = estimate_gmtf(grad_input, grad_output_mean)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot magnitude and phase of the GMTF for all axes\n", + "grad_axes = ['x', 'y', 'z']\n", + "\n", + "n_freq = n_samples - 1\n", + "frequency_vector = (1.0 / dwell_time) * torch.arange(-n_freq, n_freq) / (2 * n_samples)\n", + "half = len(frequency_vector) // 2\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, figsize=(10, 5))\n", + "for k, label in enumerate(grad_axes):\n", + " ax1.plot(frequency_vector[half:] * 1e-3, abs(gmtf.numpy()[k, half:]), label=label)\n", + " ax2.plot(frequency_vector[half:] * 1e-3, np.unwrap(np.angle(gmtf.numpy()[k, half:])), label=label)\n", + "\n", + "ax1.set(xlabel='Frequency (kHz)', ylabel='|GMTF|', title='Magnitude', xlim=(0, 10), ylim=(0.8, 1.025))\n", + "ax2.set(xlabel='Frequency (kHz)', ylabel='arg(GMTF) (rad)', title='Phase', xlim=(0, 10), ylim=(-0.1, 0.1))\n", + "\n", + "for ax in (ax1, ax2):\n", + " ax.legend()\n", + " ax.grid(visible=True, which='major', color='#666666')\n", + " ax.minorticks_on()\n", + " ax.grid(visible=True, which='minor', color='#999999', linestyle='-', alpha=0.2)\n", + "\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mrseq_mrzero", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/grpe_flash_dixon.ipynb b/examples/grpe_flash_dixon.ipynb new file mode 100644 index 0000000..6a49333 --- /dev/null +++ b/examples/grpe_flash_dixon.ipynb @@ -0,0 +1,434 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# GRPE FLASH Dixon\n", + "\n", + "Golden radial phase encoding acquisition with a 3-point Dixon gradient echo readout. More detailed information on the acquisition can be found here:\n", + "\n", + "Mayer et al 2022. Cardio‐respiratory motion‐corrected 3D cardiac water‐fat MRI using model‐based image reconstruction. MRM. https://doi.org/10.1002/mrm.29284" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import tempfile\n", + "import time\n", + "from pathlib import Path\n", + "from urllib.request import urlretrieve\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import MRzeroCore as mr0\n", + "import numpy as np\n", + "import scipy as sp\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.enums import AcqFlags\n", + "from mrpro.data.traj_calculators import KTrajectoryIsmrmrd\n", + "from mrpro.operators.models import SpoiledGRE\n", + "from scipy.interpolate import interp1d\n", + "\n", + "from mrseq.scripts.grpe_flash_dixon import main as create_seq\n", + "from mrseq.utils import combine_ismrmrd_files\n", + "from mrseq.utils import sys_defaults" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Settings\n", + "\n", + "We are going to use a small numerical phantom with a matrix size of 32 x 32 x 32 to reduce run times.\n", + "To ensure we really acquire all data in the steady state, we use a large number of dummy excitations before the actual \n", + "image acquisition. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "# Acquisition parameters\n", + "flip_angle_degree = 12\n", + "n_dummy_spokes = 10\n", + "image_matrix_size = [32, 32, 32] # x,y,z\n", + "\n", + "# Output path\n", + "tmp = tempfile.TemporaryDirectory()\n", + "fname_mrd = Path(tmp.name) / 'grpe_flash_dixon.mrd'" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### Digital Phantom\n", + "\n", + "We use the Brainweb phantom from [MRzero](https://github.com/MRsources/MRzero-Core), interpolated to our desired matrix size, with uniform B1 field." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "# Load 3D brainweb phantom\n", + "urlretrieve(\n", + " 'https://github.com/MRsources/MRzero-Core/raw/main/documentation/playground_mr0/subject05.npz', 'subject05.npz'\n", + ")\n", + "\n", + "obj_p = mr0.VoxelGridPhantom.brainweb('subject05.npz')\n", + "phantom = obj_p.interpolate(*image_matrix_size)\n", + "phantom.size[2] = phantom.size[1] # GRPE assumes isotropic voxels in y-z plane\n", + "phantom.B1[:] = 1.0" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "### Create the 3D GRPE FLASH Sequence\n", + "\n", + "Generate a golden radial phase encoding FLASH sequence with partial Fourier along the RPE lines and partial echo along the readout." + ] + }, + { + "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", + " rf_flip_angle=flip_angle_degree,\n", + " fov_x=float(phantom.size.numpy()[0]),\n", + " fov_y=float(phantom.size.numpy()[1]),\n", + " fov_z=float(phantom.size.numpy()[1]),\n", + " n_readout=32,\n", + " n_rpe_points=32,\n", + " n_rpe_points_per_shot=4,\n", + " n_rpe_spokes=48,\n", + " partial_echo_factor=0.7,\n", + " partial_fourier_factor=0.7,\n", + " n_dummy_spokes=n_dummy_spokes,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### Simulate the Sequence\n", + "\n", + "Pass the sequence and phantom to the MRzero simulation engine and save the signal as an ISMRMRD 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", + "tstart = time.time()\n", + "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e0)\n", + "print(f'Simulation time: {(time.time() - tstart) / 60:.2f} min')\n", + "mr0.sig_to_mrd(fname_mrd, signal, sequence)\n", + "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "### Gradient correction\n", + "\n", + "We are using a bipolar readout gradient to achieve short echo times and efficient sampling. Due to hardware limitations, \n", + "the positive and negative gradient lobes are not exactly the same. This can lead to small shifts of the k-space \n", + "trajectory between positive and negative readout gradients. If these are not corrected for, there are linear phase\n", + "errors in image space along the readout direction, which can strongly impair the fat-water separation.\n", + "\n", + "There are several approaches that estimate this phase error from the imaging data. For this sequence, we acquire\n", + "additional data without phase encoding, where we change the polarity of the readout gradient. First, we acquire with \n", + "positive - negative - positive readout gradients and then negative - positive - negative readout gradients. We can then \n", + "use a cross-correlation between the pairs of positive and negative readout gradients at the same echo times to \n", + "estimate this shift and correct for it. " + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "As a first step, we get the additional data that is labeled as ACQ_IS_PHASECORR_DATA." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "kdata_corr = KData.from_file(\n", + " str(fname_mrd).replace('.mrd', '_with_traj.mrd'),\n", + " trajectory=KTrajectoryIsmrmrd(),\n", + " acquisition_filter_criterion=lambda acquisition: AcqFlags.ACQ_IS_PHASECORR_DATA.value & acquisition.flags,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "To distinguish between the \"positive - negative - positive\" and \"negative - positive - negative\" acquisition, the first \n", + "are labeled with the acquisition index repetition = 0 and the second with repetition = 1. We also use multiple averages \n", + "to reduce noise. So let's split it up into the different components. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "sort_indices = np.lexsort(\n", + " (\n", + " kdata_corr.header.acq_info.idx.contrast.squeeze(),\n", + " kdata_corr.header.acq_info.idx.repetition.squeeze(),\n", + " kdata_corr.header.acq_info.idx.average.squeeze(),\n", + " )\n", + ")\n", + "kdata_corr = kdata_corr[sort_indices.tolist(), ...].rearrange(\n", + " '(average repetition contrast)... -> average repetition contrast ...',\n", + " average=kdata_corr.header.acq_info.idx.average.unique().numel(),\n", + " repetition=kdata_corr.header.acq_info.idx.repetition.unique().numel(),\n", + " contrast=kdata_corr.header.acq_info.idx.contrast.unique().numel(),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "We are going to use the first echo for the estimation of the shift. The code below will also work for the other echoes\n", + "and should lead to the same shift." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "echo = 0\n", + "\n", + "# Get data for the current echo and compute mean over averages\n", + "kdata = torch.abs(torch.mean(kdata_corr.data[:, :, echo], dim=0)).squeeze().clone()\n", + "ktraj = kdata_corr.traj.kx[0, :, echo].squeeze().clone()\n", + "\n", + "if echo % 2 == 1: # switch positive and negative lobe to get the correct sign for the shift\n", + " kdata = kdata[(1, 0), ...]\n", + " ktraj = ktraj[(1, 0), ...]\n", + "\n", + "kspace_start = int(ktraj[0, :].min().abs().item())\n", + "ktraj_interpolated = np.linspace(-kspace_start, kspace_start - 1, 20 * kspace_start)\n", + "kdata_interp = [\n", + " interp1d(pos.numpy(), signal.numpy(), kind='linear', fill_value='extrapolate')(ktraj_interpolated)\n", + " for pos, signal in zip(ktraj, kdata, strict=True)\n", + "]\n", + "\n", + "plt.figure()\n", + "plt.plot(ktraj_interpolated, kdata_interp[0], label='positive gradient')\n", + "plt.plot(ktraj_interpolated, kdata_interp[1], label='negative gradient')\n", + "plt.legend()\n", + "\n", + "cross_correlation = sp.signal.correlate(kdata_interp[0], kdata_interp[1], mode='same')\n", + "kspace_shift = (np.argmax(cross_correlation) - len(kdata_interp[0]) // 2) / 10 # divide by interpolation factor\n", + "print(f'K-space shift (in k-space samples): {kspace_shift}')\n", + "assert kspace_shift == 0" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "### Reconstruction\n", + "\n", + "Reconstruct the 3D image from the k-space data using direct reconstruction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "# Load k-space data and reconstruct\n", + "kdata = KData.from_file(\n", + " str(fname_mrd).replace('.mrd', '_with_traj.mrd'),\n", + " trajectory=KTrajectoryIsmrmrd(),\n", + ")\n", + "\n", + "kdata.traj.kx[1] += kspace_shift # delta_k (set to actual value if needed)\n", + "\n", + "recon = DirectReconstruction(kdata, csm=None)\n", + "idata = recon(kdata)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(4, 3, figsize=(12, 10))\n", + "for cax in ax.flatten():\n", + " cax.set_xticks([])\n", + " cax.set_yticks([])\n", + "\n", + "idat = idata.data.numpy().squeeze()\n", + "z_mid, x_mid = idat.shape[-3] // 2, idat.shape[-1] // 2\n", + "\n", + "for echo in range(idat.shape[0]):\n", + " ax[0, echo].imshow(np.abs(idat[echo, :, :, x_mid]), cmap='gray')\n", + " ax[1, echo].imshow(np.abs(idat[echo, z_mid, :, :]), cmap='gray')\n", + " ax[2, echo].imshow(np.angle(idat[echo, :, :, x_mid]), cmap='bwr', vmin=-np.pi, vmax=np.pi)\n", + " ax[3, echo].imshow(np.angle(idat[echo, z_mid, :, :]), cmap='bwr', vmin=-np.pi, vmax=np.pi)\n", + " ax[0, echo].set_title(f'Echo {echo}')" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "### Compare to Theoretical Model\n", + "\n", + "Compare the reconstructed images to a theoretical spoiled GRE signal model. We calculate $T2^*$ using $1/T2^* = 1/T2 + 1/T2'$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate T2* and generate theoretical model\n", + "t2star = 1 / (1 / phantom.T2 + 1 / phantom.T2dash)\n", + "\n", + "model = SpoiledGRE(\n", + " flip_angle=np.deg2rad(flip_angle_degree),\n", + " echo_time=kdata.header.te,\n", + " repetition_time=kdata.header.tr,\n", + ")\n", + "idat_model = model(m0=phantom.PD, t1=phantom.T1, t2star=t2star)[0]\n", + "idat_model = idat_model.detach().abs().numpy().squeeze()\n", + "idat_model /= idat_model.max()\n", + "idat_model = np.roll(\n", + " rearrange(idat_model[:, ::-1, ::-1, ::-1], 'echo x y z -> echo z y x'),\n", + " shift=(1, 1, 1),\n", + " axis=(-1, -2, -3),\n", + ")\n", + "\n", + "# Create object mask\n", + "obj_mask = np.zeros_like(idat_model)\n", + "obj_mask[idat_model > 0.01] = 1\n", + "\n", + "# Normalize reconstructed image\n", + "idat = idata.data.abs().squeeze().numpy()\n", + "idat = idat * obj_mask\n", + "idat /= np.percentile(idat, 99.9)\n", + "\n", + "idat_diff = idat - idat_model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "# Comparison: axial view (middle z-slice)\n", + "fig, ax = plt.subplots(3, 3, figsize=(12, 10))\n", + "x_mid = idat.shape[-1] // 2\n", + "\n", + "for echo in range(3):\n", + " im = ax[0, echo].imshow(idat[echo, :, :, x_mid], cmap='grey', vmin=0, vmax=1)\n", + " ax[0, echo].set_title(f'Echo {echo} - Reconstructed')\n", + " ax[0, echo].set_xticks([])\n", + " ax[0, echo].set_yticks([])\n", + " fig.colorbar(im, ax=ax[0, echo])\n", + "\n", + " im = ax[1, echo].imshow(idat_model[echo, :, :, x_mid], cmap='grey', vmin=0, vmax=1)\n", + " ax[1, echo].set_title(f'Echo {echo} - Model')\n", + " ax[1, echo].set_xticks([])\n", + " ax[1, echo].set_yticks([])\n", + " fig.colorbar(im, ax=ax[1, echo])\n", + "\n", + " im = ax[2, echo].imshow(idat_diff[echo, :, :, x_mid], cmap='bwr', vmin=-0.2, vmax=0.2)\n", + " ax[2, echo].set_title(f'Echo {echo} - Difference')\n", + " ax[2, echo].set_xticks([])\n", + " ax[2, echo].set_yticks([])\n", + " fig.colorbar(im, ax=ax[2, echo])\n", + "\n", + "plt.tight_layout()\n", + "plt.show()\n", + "\n", + "\n", + "relative_error = np.sum(np.abs(idat_diff)) / np.sum(np.abs(idat_model))\n", + "print(f'Relative error: {relative_error:.4f}')\n", + "assert relative_error < 0.09" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mrseq_mrzero", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/imaging.md b/examples/imaging.md index 91c8959..9ade8b1 100644 --- a/examples/imaging.md +++ b/examples/imaging.md @@ -3,4 +3,5 @@ Sequences for qualitative MRI. Available examples are: - [M2D Radial Flash](radial_flash.ipynb) -- [M2D Spiral Flash](spiral_flash.ipynb) \ No newline at end of file +- [M2D Spiral Flash](spiral_flash.ipynb) +- [3D Golden Radial Phase Encoding with 3-point Dixon readout](grpe_flash_dixon.ipynb) \ No newline at end of file diff --git a/examples/radial_flash.ipynb b/examples/radial_flash.ipynb index 775364e..e94c724 100644 --- a/examples/radial_flash.ipynb +++ b/examples/radial_flash.ipynb @@ -49,7 +49,7 @@ "source": [ "### Settings\n", "We are going to use a numerical phantom with a matrix size of 128 x 128. \n", - "To ensure we really acquire all radial lines in the steaady state, we use a large number of dummy excitations before the actual image acquisition. " + "To ensure we really acquire all radial lines in the steady state, we use a large number of dummy excitations before the actual image acquisition. " ] }, { diff --git a/examples/t1_molli_bssfp.ipynb b/examples/t1_molli_bssfp.ipynb index b4631fd..199d4c3 100644 --- a/examples/t1_molli_bssfp.ipynb +++ b/examples/t1_molli_bssfp.ipynb @@ -245,7 +245,7 @@ "source": [ "dictionary = DictionaryMatchOp(MOLLI(ti=torch.as_tensor(ti, dtype=torch.float32)), index_of_scaling_parameter=0)\n", "dictionary.append(\n", - " torch.tensor(1.0), torch.linspace(0.1, 5.0, 1000)[None, :], torch.linspace(0.1, 5.0, 1000)[None, None, :]\n", + " torch.tensor(1.0), torch.linspace(0.1, 5.0, 300)[None, :, None], torch.linspace(0.1, 5.0, 300)[None, None, :]\n", ")\n", "m0_match, c_match, t1_match = dictionary(idata.data)\n", "\n", @@ -270,7 +270,7 @@ "\n", "relative_error = np.sum(np.abs(t1_input - t1_measured)) / np.sum(np.abs(t1_input))\n", "print(f'Relative error {relative_error}')\n", - "assert relative_error < 0.08" + "assert relative_error < 0.12" ] } ], diff --git a/src/mrseq/scripts/girf_dyn_triangle.py b/src/mrseq/scripts/girf_dyn_triangle.py new file mode 100644 index 0000000..e66e095 --- /dev/null +++ b/src/mrseq/scripts/girf_dyn_triangle.py @@ -0,0 +1,338 @@ +"""GIRF sequence creator with triangular gradients.""" + +from collections.abc import Sequence +from pathlib import Path + +import numpy as np +import pypulseq as pp + +from mrseq.utils import sys_defaults +from mrseq.utils import write_sequence + + +def girf_triangle_kernel( + system: pp.Opts, + rf_flip_angle: float, + rf_duration: float, + rf_bwt: float, + apodization: float, + n_avg: int, + adc_duration: float, + dwell_time: float, + slice_thickness: float, + slice_pos: Sequence[float], + g_amplitude_coeff: Sequence[float], + tr: float, + adc_delay: float, + rise_times: Sequence[float], + grad_free: float, + trig_duration: float, + cam_acq_duration: float, + cam_interleave_tr: float, + cam_acq_delay: float, + cam_nr_sync_dyn: int, +) -> pp.Sequence: + """Generate a GIRF sequence with triangular gradients. + + Parameters + ---------- + system + PyPulseq system limits object. + rf_flip_angle + Flip angle of RF excitation pulse (in degrees). + rf_duration + Duration of RF excitation pulse (in seconds). + rf_bwt + Bandwidth-time product of RF excitation pulse (Hz * seconds). + apodization + Apodization factor of RF excitation pulse. + n_avg + Number of averages. + adc_duration + Duration of ADC acquisition (in seconds). + dwell_time + ADC dwell time (in seconds). + slice_thickness + Slice thickness (in meters). + slice_pos + List of slice positions (in meters). + g_amplitude_coeff + List of amplitude coefficients for gradient encoding. + tr + Repetition time (in seconds). + adc_delay + Delay between gradient and ADC to minimize eddy current effects (in seconds). + rise_times + List of gradient rise times (in seconds). + grad_free + Gradient-free period after trigger (in seconds). + trig_duration + Duration of trigger pulse (in seconds). + cam_acq_duration + Camera acquisition duration (in seconds). + cam_interleave_tr + Camera interleave TR (in seconds). + cam_acq_delay + Camera acquisition delay (in seconds). + cam_nr_sync_dyn + Number of camera sync dynamics. + + Returns + ------- + seq + PyPulseq Sequence object. + """ + if n_avg < 0: + raise ValueError('Number of averages must be >= 0.') + + if not rise_times: + raise ValueError('Rise times list cannot be empty.') + + # Create PyPulseq Sequence object + seq = pp.Sequence(system=system) + + # Create and set sinc pulse parameters + rf, gz, gzr = pp.make_sinc_pulse( + flip_angle=rf_flip_angle / 180 * np.pi, + delay=system.rf_dead_time, + duration=rf_duration, + slice_thickness=slice_thickness, + apodization=apodization, + time_bw_product=rf_bwt, + system=system, + return_gz=True, + ) + + # Create necessary delays + grad_free_time = pp.make_delay(grad_free) + delay_tr = pp.make_delay(tr) + delay_ec = pp.make_delay(adc_delay) + + # Create trigger pulse for MFC + trig = pp.make_digital_output_pulse(channel='ext1', duration=trig_duration, delay=0) + + # Calculate ADC parameters + n_readout = int(np.round(adc_duration / dwell_time)) + adc = pp.make_adc( + num_samples=n_readout, + duration=adc_duration, + system=system, + delay=system.adc_dead_time, + ) + + # Calculate total number of dynamics for camera + grad_channels = ['x', 'y', 'z'] + cam_nr_dynamics = len(g_amplitude_coeff) * len(slice_pos) * len(grad_channels) * len(rise_times) * n_avg + print(f'Total dynamics for camera: {cam_nr_dynamics}') + + # Build sequence with nested loops + for avg in range(n_avg): + avg_label = pp.make_label(type='SET', label='AVG', value=avg) + + for idx_rise_time, rise_time_val in enumerate(rise_times): + seg_label = pp.make_label(type='SET', label='SEG', value=idx_rise_time) + rise_time_val = np.round(rise_time_val, decimals=6) + + for idx_grad, grad_channel in enumerate(['x', 'y', 'z']): + grad_label = pp.make_label(type='SET', label='PHS', value=idx_grad) + + for idx_slice_pos, slice_pos_val in enumerate(slice_pos): + slice_label = pp.make_label(type='SET', label='SLC', value=idx_slice_pos) + + for idx_amp_fac, amp_fac in enumerate(g_amplitude_coeff): + amp_fac_label = pp.make_label(type='SET', label='REP', value=idx_amp_fac) + + # Calculate amplitude from coefficient, slew rate, and rise time + amp = amp_fac * system.max_slew * rise_time_val + + # Set RF frequency offset for current slice + rf.freq_offset = gz.amplitude * slice_pos_val + gz.channel = grad_channel + gzr.channel = grad_channel + + # Add RF pulse with slice selection + seq.add_block(rf, gz, avg_label, seg_label, grad_label, slice_label, amp_fac_label) + seq.add_block(gzr) + + # Add eddy current compensation delay + seq.add_block(delay_ec) + seq.add_block(trig) + seq.add_block(grad_free_time) + + # Create triangle gradient (trapezoid with no flat time) + g_triangle = pp.make_trapezoid( + channel=grad_channel, + system=system, + rise_time=rise_time_val, + flat_time=0, + amplitude=amp, + delay=0, + ) + seq.add_block(adc, g_triangle) + + # Add TR delay + seq.add_block(delay_tr) + + seq.add_block(avg_label) + + # Set sequence definitions + seq.set_definition('FOV', [0.5, 0.5, slice_thickness]) + seq.set_definition('ReconMatrix', (n_readout, 1, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('SlicePos', slice_pos) + seq.set_definition('TE', 0) + seq.set_definition('TR', tr) + seq.set_definition('RiseTimes', rise_times) + seq.set_definition('CameraNrDynamics', cam_nr_dynamics) + seq.set_definition('CameraNrSyncDynamics', cam_nr_sync_dyn) + seq.set_definition('CameraAcqDuration', cam_acq_duration) + seq.set_definition('CameraInterleaveTR', cam_interleave_tr) + seq.set_definition('CameraAcqDelay', cam_acq_delay) + + seq.set_definition('AdcDuration', adc_duration) + seq.set_definition('DwellTime', dwell_time) + seq.set_definition('SlewRate', system.max_slew) + seq.set_definition('AdcDelay', adc_delay) + seq.set_definition('GradAmplitudeCoeff', g_amplitude_coeff) + + return seq + + +def main( + system: pp.Opts | None = None, + n_avg: int = 3, + tr: float = 2.0, + slice_thickness: float = 1.5e-3, + slice_pos: list[float] | None = None, + show_plots: bool = True, + test_report: bool = True, + timing_check: bool = True, + v141_compatibility: bool = True, +) -> tuple[pp.Sequence, Path]: + """Generate a sequence with triangular gradients to estimate the GIRF. + + This sequence allows for the calculation of a gradient impulse response function (GIRF) using the Dyn + method [DUY1998]_ or the use of a field camera [VAN2013]_ . + + + .. [DUY1998] Dyn J, Yang Y, Frank JA, and van der Veen JW (1998), Simple Correction Method fork-Space Trajectory + Deviations in MRI. JMR 132, 150-153. https://doi.org/10.1006/jmre.1998.1396 + + .. [VAN2013] Vannesjo SJ, Haeberlin M, Kasper L, Pavan M, Wilm, BJ, Barnet C, and PRuessman KP (2013), + Gradient system characterization by impulse response measurements with a dynamic field camera. MRM 69. 583-593. + https://doi.org/10.1002/mrm.24263 + + + Parameters + ---------- + system + PyPulseq system limits object. Uses default system if None. + n_avg + Number of averages. + tr + Repetition time (in seconds). + slice_thickness + Slice thickness (in meters). + slice_pos + List of slice positions (in meters). + 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 GIRF triangle sequence. + file_path + Path to the sequence file. + + """ + if system is None: + system = sys_defaults + + if slice_pos is None: + slice_pos = [0.04] + + # Define RF excitation pulse parameters + rf_flip_angle = 90 # flip angle [degrees] + rf_duration = 8.4e-3 # duration of the RF excitation pulse [s] + rf_bwt = 4 # bandwidth-time product of RF excitation pulse [Hz*s] + apodization = 0.5 # apodization factor of RF excitation pulse + + # Define ADC and gradient timing + adc_duration = 0.06 # ADC acquisition duration [s] + dwell_time = 10e-6 # ADC dwell time [s] + + # Define sequence parameters + g_amplitude_coeff = [-1.0, 0.0] # amplitude coefficients for gradient encoding + adc_delay = 10e-6 # eddy current compensation delay [s] + rise_times = [5e-5, 6e-5, 7e-5, 8e-5, 9e-5, 1e-4, 1.1e-4, 1.2e-4, 1.3e-4, 1.4e-4, 1.5e-4, 1.6e-4] # rise times [s] + + # Define camera and trigger parameters + grad_free = 0.5e-3 # gradient-free period after trigger [s] + trig_duration = 10e-6 # trigger pulse duration [s] + cam_acq_duration = 0.07 # camera acquisition duration [s] + cam_interleave_tr = 0.4 # camera interleave TR [s] + cam_acq_delay = 0.0 # camera acquisition delay [s] + cam_nr_sync_dyn = 0 # number of camera sync dynamics + + # Define sequence filename + filename = f'{Path(__file__).stem}' + + output_path = Path.cwd() / 'output' + output_path.mkdir(parents=True, exist_ok=True) + + seq = girf_triangle_kernel( + system=system, + rf_flip_angle=rf_flip_angle, + rf_duration=rf_duration, + rf_bwt=rf_bwt, + apodization=apodization, + n_avg=n_avg, + adc_duration=adc_duration, + dwell_time=dwell_time, + slice_thickness=slice_thickness, + slice_pos=slice_pos, + g_amplitude_coeff=g_amplitude_coeff, + tr=tr, + adc_delay=adc_delay, + rise_times=rise_times, + grad_free=grad_free, + trig_duration=trig_duration, + cam_acq_duration=cam_acq_duration, + cam_interleave_tr=cam_interleave_tr, + cam_acq_delay=cam_acq_delay, + cam_nr_sync_dyn=cam_nr_sync_dyn, + ) + + # Check sequence timing + 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 test report + if test_report: + print('\nCreating advanced test report...') + print(seq.test_report()) + + # Save sequence 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(grad_disp='mT/m') + + return seq, output_path / filename + + +if __name__ == '__main__': + main() diff --git a/src/mrseq/scripts/grpe_flash_dixon.py b/src/mrseq/scripts/grpe_flash_dixon.py new file mode 100644 index 0000000..2265a91 --- /dev/null +++ b/src/mrseq/scripts/grpe_flash_dixon.py @@ -0,0 +1,602 @@ +"""3D FLASH sequence with radial phase encoding.""" + +from pathlib import Path + +import ismrmrd +import numpy as np +import pypulseq as pp +from pypulseq.rotate import rotate + +from mrseq.utils import MultiEchoAcquisition +from mrseq.utils import find_gx_flat_time_on_adc_raster +from mrseq.utils import round_to_raster +from mrseq.utils import sys_defaults +from mrseq.utils import write_sequence +from mrseq.utils.constants import GOLDEN_ANGLE_HALF_CIRCLE +from mrseq.utils.ismrmrd import Fov +from mrseq.utils.ismrmrd import Limits +from mrseq.utils.ismrmrd import MatrixSize +from mrseq.utils.ismrmrd import create_header + + +def grpe_flash_dixon_kernel( + system: pp.Opts, + te: float | None, + delta_te: float | None, + n_echoes: int, + tr: float | None, + fov_x: float, + fov_y: float, + fov_z: float, + n_readout: int, + n_rpe_points: int, + n_rpe_points_per_shot: int, + n_rpe_spokes: int, + readout_oversampling: float, + partial_echo_factor: float, + partial_fourier_factor: float, + n_dummy_spokes: int, + gx_pre_duration: float, + gx_flat_time: float, + rf_duration: float, + rf_flip_angle: float, + rf_bwt: float, + rf_apodization: float, + rf_spoiling_phase_increment: float, + gx_spoil_duration: float, + gx_spoil_area: float, + mrd_header_file: str | Path | None, +) -> tuple[pp.Sequence, float, float, float]: + """Generate a 3D FLASH sequence with golden radial phase encoding. + + Parameters + ---------- + system + PyPulseq system limits object. + te + Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. + delta_te + Desired echo spacing (in seconds). Minimum echo spacing is used if set to None. + n_echoes + Number of echoes. + tr + Desired repetition time (TR) (in seconds). + fov_x + Field of view in x direction (in meters). + fov_y + Field of view in y direction (in meters). + fov_z + Field of view in z direction (in meters). + n_readout + Number of frequency encoding steps. + n_rpe_points + Number of radial phase encoding points (points along one RPE line). + n_rpe_points_per_shot + Shots are interleaved groups of points along RPE lines. Each shot obtains the k-space center at the cost of a + point of the highest k-space frequency. Fat-saturation pulses are applied prior to each shot. + n_rpe_spokes + Number of radial phase encoding spokes (number of RPE lines). + readout_oversampling + Readout oversampling factor, commonly 2. This reduces aliasing artifacts. + partial_echo_factor + Partial echo factor along readout (between 0.5 and 1). + partial_fourier_factor + Partial Fourier factor along RPE lines (between 0.5 and 1). + n_dummy_spokes + Number of dummy RPE spokes before data acquisition to ensure steady state. + gx_pre_duration + Duration of readout pre-winder gradient (in seconds) + gx_flat_time + Flat time of readout 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. + gx_spoil_duration + Duration of spoiler gradient (in seconds) + gx_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. + delta_te + Time between echoes. + + """ + if readout_oversampling < 1: + raise ValueError('Readout oversampling factor must be >= 1.') + + if n_dummy_spokes < 0: + raise ValueError('Number of dummy spokes must be >= 0.') + + spoke_angle = GOLDEN_ANGLE_HALF_CIRCLE + rpe_radial_shift = [0, 0.5, 0.25, 0.75] + + if partial_echo_factor > 1 or partial_echo_factor < 0.5: + raise ValueError('Partial echo factor has to be within 0.5 and 1') + if partial_fourier_factor > 1 or partial_fourier_factor < 0.5: + raise ValueError('Partial Fourier factor has to be within 0.5 and 1') + + # create PyPulseq Sequence object and set system limits + seq = pp.Sequence(system=system) + + # create slab 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=fov_z, + apodization=rf_apodization, + time_bw_product=rf_bwt, + delay=system.rf_dead_time, + system=system, + return_gz=True, + use='excitation', + ) + + multi_echo_gradient = MultiEchoAcquisition( + system=system, + delta_te=delta_te, + fov=fov_x, + n_readout=n_readout, + readout_oversampling=readout_oversampling, + partial_echo_factor=partial_echo_factor, + gx_flat_time=gx_flat_time, + gx_pre_duration=gx_pre_duration, + ) + + # calculate gradient areas for phase encoding along each RPE line in a low-high order + delta_ky = 1 / fov_y + n_rpe_points_with_partial_fourier = int( + ((n_rpe_points * partial_fourier_factor) // n_rpe_points_per_shot) * n_rpe_points_per_shot + ) + n_shots_per_rpe_spoke = n_rpe_points_with_partial_fourier // n_rpe_points_per_shot + print( + f'Number of phase encoding points {n_rpe_points_with_partial_fourier}', + f'with partial Fourier factor {partial_fourier_factor}', + ) + enc_steps_pe = np.arange(0, n_rpe_points_with_partial_fourier) + phase_areas = (enc_steps_pe - n_rpe_points / 2) * delta_ky + centric_idx = np.argsort(np.abs(phase_areas), kind='stable') + enc_steps_pe = enc_steps_pe[centric_idx] + + # Interleave shots + enc_steps_pe_interleaved = [] + for step in range(n_shots_per_rpe_spoke): + enc_steps_pe_interleaved.append(enc_steps_pe[step::n_shots_per_rpe_spoke]) + enc_steps_pe = np.concatenate(enc_steps_pe_interleaved) + + # create spoiler gradients + gx_spoil = pp.make_trapezoid(channel='x', system=system, area=gx_spoil_area, duration=gx_spoil_duration) + + # calculate minimum echo time + if te is None: + gzr_gx_dur = pp.calc_duration(gzr, multi_echo_gradient._gx_pre) # gzr and gx_pre are applied simultaneously + else: + gzr_gx_dur = pp.calc_duration(gzr) + pp.calc_duration( + multi_echo_gradient._gx_pre + ) # gzr and gx_pre are applied sequentially + + 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 + + gzr_gx_dur # slice selection re-phasing gradient and readout pre-winder + + multi_echo_gradient._gx.delay # potential delay of readout gradient + + multi_echo_gradient._gx.rise_time # rise time of readout gradient + + (multi_echo_gradient._n_readout_pre_echo + 0.5) * multi_echo_gradient._adc.dwell + ).item() + + # 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 + + gzr_gx_dur # slice selection re-phasing gradient and readout pre-winder + + pp.calc_duration(multi_echo_gradient._gx) * n_echoes # readout gradient + + pp.calc_duration(multi_echo_gradient._gx_between) * (n_echoes - 1) # readout gradient + + pp.calc_duration(gx_spoil, multi_echo_gradient._gx_post) # gradient spoiler or readout-re-winder + ).item() + + # 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 not 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_x * readout_oversampling, y=fov_y, z=fov_y), + recon_fov=Fov(x=fov_x, y=fov_y, z=fov_y), + encoding_matrix=MatrixSize(n_x=int(n_readout * readout_oversampling), n_y=n_rpe_points, n_z=n_rpe_points), + recon_matrix=MatrixSize(n_x=n_readout, n_y=n_rpe_points, n_z=n_rpe_points), + dwell_time=multi_echo_gradient._adc.dwell, + slice_limits=Limits(min=0, max=0, center=0), + k1_limits=Limits(min=0, max=n_rpe_points, center=0), + k2_limits=Limits(min=0, max=n_rpe_spokes, 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='PAR', type='SET', value=0)) + seq.add_block(multi_echo_gradient._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=3, number_of_samples=multi_echo_gradient._adc.num_samples) + prot.append_acquisition(acq) + + # Define sequence blocks + for se_index in np.arange(-n_dummy_spokes, n_rpe_spokes): + se_label = pp.make_label(type='SET', label='PAR', value=int(se_index)) + for shot_index in range(n_shots_per_rpe_spoke): + for pe_index in range(n_rpe_points_per_shot): + pe = int(enc_steps_pe[shot_index * n_rpe_points_per_shot + pe_index]) + pe_label = pp.make_label(type='SET', label='LIN', value=pe) + + # set rf pulse properties and add rf pulse block event + if rf_spoiling_phase_increment > 0: + rf.phase_offset = rf_phase / 180 * np.pi + multi_echo_gradient._adc.phase_offset = rf_phase / 180 * np.pi + + # add rf pulse + seq.add_block(rf, gz) + + # update rf pulse properties for the next loop + rf_inc = divmod(rf_inc + rf_spoiling_phase_increment, 360.0)[1] + rf_phase = divmod(rf_phase + rf_inc, 360.0)[1] + + # area of phase encoding gradient + current_phase_area = phase_areas[pe] + # do not shift k-space center for respiratory navigator calculation + if np.abs(current_phase_area) > 1e-5: + pe_shift = rpe_radial_shift[int(np.mod(se_index, len(rpe_radial_shift)))] + else: + pe_shift = 0.0 + + # phase encoding along pe + gy_pre = pp.make_trapezoid( + channel='y', + area=current_phase_area + pe_shift * delta_ky, + duration=pp.calc_duration(multi_echo_gradient._gx_pre), + system=system, + ) + + # calculate rotated phase encoding gradient + rotation_angle_rad = spoke_angle * se_index + gy_pre_rotated = rotate(gy_pre, angle=rotation_angle_rad, axis='x') + gy_pre = gy_pre_rotated[0] + if len(gy_pre_rotated) == 2: + gz_pre = gy_pre_rotated[1] + else: + gz_pre = pp.make_trapezoid( + channel='z', + area=0, + duration=pp.calc_duration(multi_echo_gradient._gx_pre), + system=system, + ) + assert gy_pre.channel == 'y' and gz_pre.channel == 'z' + + # combine slice rephaser and slice encoding gradient + gz_r_pre = pp.make_trapezoid( + channel='z', + area=gz_pre.area + gzr.area, + duration=pp.calc_duration(multi_echo_gradient._gx_pre), + system=system, + ) + + label_contents = [pe_label, se_label] + seq.add_block(multi_echo_gradient._gx_pre, gy_pre, gz_r_pre, *label_contents) + + # add delay due to TE + if te_delay > 0: + seq.add_block(pp.make_delay(te_delay)) + + # add readout gradients and ADCs + if se_index >= 0: + seq, _ = multi_echo_gradient.add_to_seq_without_pre_post_gradient(seq, n_echoes) + else: + # the most accurate way to get the duration of the readout block is to add it to a dummy sequence + seq_dummy = pp.Sequence(system=system) + seq_dummy, _ = multi_echo_gradient.add_to_seq_without_pre_post_gradient(seq_dummy, n_echoes) + readout_duration = sum(seq_dummy.block_durations.values()) + seq.add_block(pp.make_delay(readout_duration)) + + # add re-winder and spoiler gradients + gy_pre.amplitude = -gy_pre.amplitude + gz_pre.amplitude = -gz_pre.amplitude + seq.add_block(gx_spoil, gy_pre, gz_pre) + + # add delay in case TR > min_TR + if tr_delay > 0: + seq.add_block(pp.make_delay(tr_delay)) + + if mrd_header_file and se_index >= 0: + # add acquisitions to metadata + k0_trajectory = np.linspace( + -multi_echo_gradient._n_readout_pre_echo, + multi_echo_gradient._n_readout_post_echo, + multi_echo_gradient._n_readout_with_partial_echo, + ) + grpe_trajectory = np.zeros((multi_echo_gradient._n_readout_with_partial_echo, 3), dtype=np.float32) + + for echo_ in range(n_echoes): + gx_sign = (-1) ** echo_ + grpe_trajectory[:, 0] = k0_trajectory * gx_sign + grpe_trajectory[:, 1] = (pe - n_rpe_points / 2 + pe_shift) * np.cos(rotation_angle_rad) + grpe_trajectory[:, 2] = (pe - n_rpe_points / 2 + pe_shift) * np.sin(rotation_angle_rad) + + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=3, number_of_samples=multi_echo_gradient._adc.num_samples) + acq.traj[:] = grpe_trajectory + prot.append_acquisition(acq) + + # obtain echoes with positive and negative gradient polarity to be able to correct for any differences between + # positive and negative gradient waveforms in the bi-polar readout + seq.add_block(pp.make_label(type='SET', label='NAV', value=True)) + for average_idx in range(4): + average_label = pp.make_label(type='SET', label='AVG', value=average_idx) + for polarity_idx, polarity in enumerate(['positive', 'negative']): + polarity_label = pp.make_label(type='SET', label='REP', value=polarity_idx) + # set rf pulse properties and add rf pulse block event + if rf_spoiling_phase_increment > 0: + rf.phase_offset = rf_phase / 180 * np.pi + multi_echo_gradient._adc.phase_offset = rf_phase / 180 * np.pi + + # add rf pulse + seq.add_block(rf, gz) + + # update rf pulse properties for the next loop + rf_inc = divmod(rf_inc + rf_spoiling_phase_increment, 360.0)[1] + rf_phase = divmod(rf_phase + rf_inc, 360.0)[1] + + seq.add_block( + multi_echo_gradient._gx_pre + if polarity == 'positive' + else pp.scale_grad(multi_echo_gradient._gx_pre, -1), + average_label, + polarity_label, + ) + + # add delay due to TE + if te_delay > 0: + seq.add_block(pp.make_delay(te_delay)) + + # add readout gradients and ADCs + seq, time_to_echoes = multi_echo_gradient.add_to_seq_without_pre_post_gradient(seq, n_echoes, polarity) # type: ignore[arg-type] + + # add spoiler gradients + seq.add_block(gx_spoil) + + # add delay in case TR > min_TR + if tr_delay > 0: + seq.add_block(pp.make_delay(tr_delay)) + + if mrd_header_file: + # add acquisitions to metadata + k0_trajectory = np.linspace( + -multi_echo_gradient._n_readout_pre_echo, + multi_echo_gradient._n_readout_post_echo, + multi_echo_gradient._n_readout_with_partial_echo, + ) + grpe_trajectory = np.zeros((multi_echo_gradient._n_readout_with_partial_echo, 3), dtype=np.float32) + + for echo_ in range(n_echoes): + gx_sign = (-1) ** echo_ if polarity == 'positive' else (-1) ** (echo_ + 1) + grpe_trajectory[:, 0] = k0_trajectory * gx_sign + + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=3, number_of_samples=multi_echo_gradient._adc.num_samples) + acq.traj[:] = grpe_trajectory + prot.append_acquisition(acq) + + seq.add_block(pp.make_label(label='NAV', type='SET', value=False)) + + # close ISMRMRD file + if mrd_header_file: + prot.close() + + delta_te_array = np.diff(time_to_echoes) + return seq, float(min_te), float(min_tr), float(delta_te_array[0]) + + +def main( + system: pp.Opts | None = None, + tr: float | None = None, + rf_flip_angle: float = 12, + fov_x: float = 128e-3, + fov_y: float = 128e-3, + fov_z: float = 128e-3, + n_readout: int = 128, + n_rpe_points: int = 128, + n_rpe_points_per_shot: int = 8, + n_rpe_spokes: int = 16, + partial_echo_factor: float = 0.7, + partial_fourier_factor: float = 0.7, + receiver_bandwidth_per_pixel: float = 1200, # Hz/pixel + n_dummy_spokes: int = 2, + show_plots: bool = True, + test_report: bool = True, + timing_check: bool = True, + v141_compatibility: bool = True, +) -> pp.Sequence: + """Generate a 3D FLASH sequence with radial phase encoding. + + Parameters + ---------- + system + PyPulseq system limits object. + 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_x + Field of view in x direction (in meters). + fov_y + Field of view in y direction (in meters). + fov_z + Field of view in z direction (in meters). + n_readout + Number of frequency encoding steps. + n_rpe_points + Number of radial phase encoding points (points along one RPE line). + n_rpe_points_per_shot + Shots are interleaved groups of points along RPE lines. Each shot obtains the k-space center at the cost of a + point of the highest k-space frequency. Fat-saturation pulses are applied prior to each shot. + n_rpe_spokes + Number of radial phase encoding spokes (number of RPE lines). + partial_echo_factor + Partial echo factor along the readout (between 0.5 and 1). + partial_fourier_factor + Partial Fourier factor along RPE lines (between 0.5 and 1). + receiver_bandwidth_per_pixel + Desired receiver bandwidth per pixel (in Hz/pixel). This is used to calculate the readout duration. + n_dummy_spokes + Number of dummy RPE spokes 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. + """ + if system is None: + system = sys_defaults + + # define settings of rf excitation pulse + rf_duration = 0.6e-3 # duration of the rf excitation pulse [s] + rf_bwt = 2 # 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. + + # this is just approximately, the final calculation is done in the kernel + n_readout_with_oversampling = int(n_readout * readout_oversampling * partial_echo_factor) + # define ADC and gradient timing + adc_dwell_time = 1.0 / (receiver_bandwidth_per_pixel * n_readout_with_oversampling) + gx_pre_duration = 1.0e-3 # duration of readout pre-winder gradient [s] + gx_flat_time, adc_dwell_time = find_gx_flat_time_on_adc_raster( + n_readout_with_oversampling, adc_dwell_time, system.grad_raster_time, system.adc_raster_time + ) + + te = None # shortest possible echo time + delta_te = None # shortest possible delta echo time + n_echoes = 3 + + # define spoiling + gx_spoil_duration = 1.9e-3 # duration of spoiler gradient [s] + gx_spoil_area = readout_oversampling * n_readout * 1 / fov_x # 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}_fov{int(fov_x * 1000)}_{int(fov_y * 1000)}_{int(fov_z * 1000)}mm_' + filename += f'{n_readout}_{n_rpe_points}_{n_rpe_spokes}_3ne' + + 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, delta_te = grpe_flash_dixon_kernel( + system=system, + te=te, + delta_te=delta_te, + n_echoes=n_echoes, + tr=tr, + fov_x=fov_x, + fov_y=fov_y, + fov_z=fov_z, + n_readout=n_readout, + n_rpe_points=n_rpe_points, + n_rpe_points_per_shot=n_rpe_points_per_shot, + n_rpe_spokes=n_rpe_spokes, + readout_oversampling=readout_oversampling, + partial_echo_factor=partial_echo_factor, + partial_fourier_factor=partial_fourier_factor, + n_dummy_spokes=n_dummy_spokes, + gx_pre_duration=gx_pre_duration, + gx_flat_time=gx_flat_time, + 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, + gx_spoil_duration=gx_spoil_duration, + gx_spoil_area=gx_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()) + + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_x, fov_y, fov_z]) + seq.set_definition('ReconMatrix', (n_readout, n_rpe_points, n_rpe_points)) + seq.set_definition('SliceThickness', fov_z) + seq.set_definition('TE', [(te or min_te) + idx * delta_te for idx in range(n_echoes)]) + seq.set_definition('TR', tr or min_tr) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) + + # 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, (n_dummy_spokes * n_rpe_points + 20) * (tr or min_tr))) + + return seq, output_path / filename + + +if __name__ == '__main__': + main() diff --git a/src/mrseq/utils/__init__.py b/src/mrseq/utils/__init__.py index 00d413e..31bd0a0 100644 --- a/src/mrseq/utils/__init__.py +++ b/src/mrseq/utils/__init__.py @@ -1,5 +1,5 @@ from mrseq.utils.sequence_helper import find_gx_flat_time_on_adc_raster, round_to_raster, write_sequence from mrseq.utils.system_defaults import sys_defaults +from mrseq.utils.trajectory import cartesian_phase_encoding, spiral_acquisition, MultiEchoAcquisition from mrseq.utils.vds import variable_density_spiral_trajectory -from mrseq.utils.trajectory import spiral_acquisition from mrseq.utils.ismrmrd import create_header, combine_ismrmrd_files diff --git a/src/mrseq/utils/trajectory.py b/src/mrseq/utils/trajectory.py index 91beacd..c40b5ca 100644 --- a/src/mrseq/utils/trajectory.py +++ b/src/mrseq/utils/trajectory.py @@ -7,10 +7,10 @@ import numpy as np import pypulseq as pp -from mrseq.utils import find_gx_flat_time_on_adc_raster -from mrseq.utils import round_to_raster -from mrseq.utils import sys_defaults -from mrseq.utils import variable_density_spiral_trajectory +from mrseq.utils.sequence_helper import find_gx_flat_time_on_adc_raster +from mrseq.utils.sequence_helper import round_to_raster +from mrseq.utils.system_defaults import sys_defaults +from mrseq.utils.vds import variable_density_spiral_trajectory def cartesian_phase_encoding( @@ -63,10 +63,10 @@ def cartesian_phase_encoding( if n_phase_encoding_per_shot and sampling_order != 'random': kpe_extended = np.arange(-n_phase_encoding, n_phase_encoding) kpe_extended = kpe_extended[np.argsort(np.abs(kpe_extended), kind='stable')] - idx = 0 + kidx = 0 while np.mod(len(kpe), n_phase_encoding_per_shot) > 0: - kpe = np.unique(np.concatenate((kpe, (kpe_extended[idx],)))) - idx += 1 + kpe = np.unique(np.concatenate((kpe, (kpe_extended[kidx],)))) + kidx += 1 # Different temporal orders of phase encoding points if sampling_order == 'random': @@ -78,11 +78,11 @@ def cartesian_phase_encoding( elif sampling_order == 'linear': kpe = np.sort(kpe) elif sampling_order == 'low_high': - sort_idx = np.argsort(np.abs(kpe), kind='stable') - kpe = kpe[sort_idx] + idx = np.argsort(np.abs(kpe), kind='stable') + kpe = kpe[idx] elif sampling_order == 'high_low': - sort_idx = np.argsort(-np.abs(kpe), kind='stable') - kpe = kpe[sort_idx] + idx = np.argsort(-np.abs(kpe), kind='stable') + kpe = kpe[idx] else: raise ValueError(f'sampling order {sampling_order} not supported.') @@ -205,16 +205,15 @@ def __init__( min_delta_te = pp.calc_duration(self._gx) + pp.calc_duration(self._gx_between) if delta_te is None: - self._te_delay = 0.0 + self._delta_te_delay = 0.0 else: - self._te_delay = round_to_raster(delta_te - min_delta_te, self._system.block_duration_raster) - if not self._te_delay >= 0: + self._delta_te_delay = round_to_raster(delta_te - min_delta_te, self._system.block_duration_raster) + if self._delta_te_delay < 0: raise ValueError( - f'Delta TE must be larger than {min_delta_te * 1000:.2f} ms. ' - f'Current value is {delta_te * 1000:.2f} ms.' + f'TE must be larger than {min_delta_te * 1000:.3f} ms. Current value is {delta_te * 1000:.3f} ms.' ) - def add_to_seq(self, seq: pp.Sequence, n_echoes: int) -> tuple[pp.Sequence, list[float]]: + def add_to_seq(self, seq: pp.Sequence, n_echoes: int, polarity: Literal['positive', 'negative'] = 'positive'): """Add all gradients and adc to sequence. Parameters @@ -223,6 +222,8 @@ def add_to_seq(self, seq: pp.Sequence, n_echoes: int) -> tuple[pp.Sequence, list PyPulseq Sequence object. n_echoes Number of echoes + polarity + Polarity of first readout gradient Returns ------- @@ -232,17 +233,19 @@ def add_to_seq(self, seq: pp.Sequence, n_echoes: int) -> tuple[pp.Sequence, list Time from beginning of sequence to echoes. """ # readout pre-winder - seq.add_block(self._gx_pre) + seq.add_block(self._gx_pre if polarity == 'positive' else pp.scale_grad(self._gx_pre, -1)) # add readout gradients and ADCs - seq, time_to_echoes = self.add_to_seq_without_pre_post_gradient(seq, n_echoes) + seq, time_to_echoes = self.add_to_seq_without_pre_post_gradient(seq, n_echoes, polarity) # readout re-winder - seq.add_block(self._gx_post) + seq.add_block(self._gx_post if polarity == 'positive' else pp.scale_grad(self._gx_post, -1)) return seq, time_to_echoes - def add_to_seq_without_pre_post_gradient(self, seq: pp.Sequence, n_echoes: int) -> tuple[pp.Sequence, list[float]]: + def add_to_seq_without_pre_post_gradient( + self, seq: pp.Sequence, n_echoes: int, polarity: Literal['positive', 'negative'] = 'positive' + ): """Add readout gradients without pre- and re-winder gradients. Often the pre- and re-winder gradients are played out at the same time as phase encoding gradients or spoiler @@ -254,6 +257,8 @@ def add_to_seq_without_pre_post_gradient(self, seq: pp.Sequence, n_echoes: int) PyPulseq Sequence object. n_echoes Number of echoes + polarity + Polarity of first readout gradient Returns ------- @@ -266,9 +271,8 @@ def add_to_seq_without_pre_post_gradient(self, seq: pp.Sequence, n_echoes: int) time_to_echoes = [] for echo_ in range(n_echoes): start_of_current_gx = sum(seq.block_durations.values()) - gx_sign = (-1) ** echo_ + gx_sign = (-1) ** echo_ if polarity == 'positive' else (-1) ** (echo_ + 1) labels = [] - labels.append(pp.make_label(type='SET', label='REV', value=gx_sign == -1)) labels.append(pp.make_label(label='REV', type='SET', value=gx_sign == -1)) labels.append(pp.make_label(label='ECO', type='SET', value=echo_)) seq.add_block(pp.scale_grad(self._gx, gx_sign), self._adc, *labels) @@ -277,8 +281,8 @@ def add_to_seq_without_pre_post_gradient(self, seq: pp.Sequence, n_echoes: int) ) start_of_current_gx = sum(seq.block_durations.values()) if echo_ < n_echoes - 1: - if self._te_delay > 0: - seq.add_block(pp.make_delay(self._te_delay)) + if self._delta_te_delay > 0: + seq.add_block(pp.make_delay(self._delta_te_delay)) seq.add_block(pp.scale_grad(self._gx_between, -gx_sign)) return seq, time_to_echoes diff --git a/tests/scripts/test_girf_dyn_triangle.py b/tests/scripts/test_girf_dyn_triangle.py new file mode 100644 index 0000000..2b7d6be --- /dev/null +++ b/tests/scripts/test_girf_dyn_triangle.py @@ -0,0 +1,13 @@ +"""Tests for GIRF Dyn triangle sequence.""" + +import pytest +from mrseq.scripts.girf_dyn_triangle import main as create_seq + +EXPECTED_DUR = 447.42456 # defined 2026-03-19 + + +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) diff --git a/tests/scripts/test_grpe_flash_dixon.py b/tests/scripts/test_grpe_flash_dixon.py new file mode 100644 index 0000000..34759eb --- /dev/null +++ b/tests/scripts/test_grpe_flash_dixon.py @@ -0,0 +1,48 @@ +"""Tests for 3D FLASH sequence with golden radial phase encoding.""" + +import pytest +from mrseq.scripts.grpe_flash_dixon import main as create_seq + +EXPECTED_DUR = 17.84831 # defined 2025-11-09 + + +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) + + +def test_seq_creation_error_on_short_tr(system_defaults): + """Test if error is raised on too short repetition time.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, tr=5e-3, show_plots=False) + + +def test_seq_duration_vary_params_without_effect(system_defaults): + """Test if sequence duration is as expected.""" + seq, _ = create_seq( + system=system_defaults, + n_rpe_points_per_shot=4, # default 8 + show_plots=False, + test_report=False, + timing_check=False, + ) + duration = seq.duration()[0] + assert duration == pytest.approx(EXPECTED_DUR) + + +def test_seq_creation_error_on_wrong_partial_echo_factor(system_defaults): + """Test if error is raised on wrong partial echo factor.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, partial_echo_factor=1.1, show_plots=False) + with pytest.raises(ValueError): + create_seq(system=system_defaults, partial_echo_factor=0.4, show_plots=False) + + +def test_seq_creation_error_on_wrong_partial_fourier_factor(system_defaults): + """Test if error is raised on wrong partial Fourier factor.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, partial_fourier_factor=1.1, show_plots=False) + with pytest.raises(ValueError): + create_seq(system=system_defaults, partial_fourier_factor=0.4, show_plots=False) diff --git a/tests/utils/test_trajectory.py b/tests/utils/test_trajectory.py index 6a3cfdf..352c852 100644 --- a/tests/utils/test_trajectory.py +++ b/tests/utils/test_trajectory.py @@ -235,7 +235,10 @@ def test_multi_gradient_echo_error_on_short_delta_te(system_defaults): @pytest.mark.parametrize('readout_oversampling', [1, 1.5, 2]) @pytest.mark.parametrize('n_readout', [64, 128, 200]) @pytest.mark.parametrize('partial_echo_factor', [1.0, 0.8, 0.7]) -def test_multi_gradient_echo_timing(n_echoes, readout_oversampling, n_readout, partial_echo_factor, system_defaults): +@pytest.mark.parametrize('polarity', ['positive', 'negative']) +def test_multi_gradient_echo_timing( + n_echoes, readout_oversampling, n_readout, partial_echo_factor, polarity, system_defaults +): """Test that zero crossing of gradient moment coincides with echo time and correct adc sample.""" seq = pp.Sequence(system=system_defaults) mecho = MultiEchoAcquisition( @@ -244,10 +247,7 @@ def test_multi_gradient_echo_timing(n_echoes, readout_oversampling, n_readout, p readout_oversampling=readout_oversampling, partial_echo_factor=partial_echo_factor, ) - seq, time_to_echoes = mecho.add_to_seq(seq, n_echoes) - if n_echoes > 1: - # Ensure that the delta TE is the same between all echoes - assert len(np.unique(np.round(np.diff(time_to_echoes), decimals=6))) == 1 + seq, time_to_echoes = mecho.add_to_seq(seq, n_echoes, polarity) # Get full waveform for readout gradient w = seq.waveforms_and_times() @@ -262,7 +262,7 @@ def test_multi_gradient_echo_timing(n_echoes, readout_oversampling, n_readout, p k0_idx = argrelextrema(np.abs(m0_intp), np.less, order=100)[0] # Remove k0-crossings at the beginning and end of the block - k0_idx = np.array([ki for ki in k0_idx if (ki > 100 and ki < len(dt) - 100)]) + k0_idx = np.asarray([ki for ki in k0_idx if (ki > 100 and ki < len(dt) - 100)]) # Zero-crossing of the readout gradient should be within +/- .5 adc dwell time of the k-space center sample which is # the (_n_readout_pre_echo + 1)th sample. This should also be the same as the time_to_echo - time