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/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/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/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)