From 5354e68e4c7ec92f8efe6fef83e258cb8124ea27 Mon Sep 17 00:00:00 2001 From: OErbil Date: Mon, 23 Feb 2026 13:48:00 +0100 Subject: [PATCH 1/7] Add GIRF sequence creator with triangles and YAML Files for simple GIRF sequence creation using triangular gradients. The YAML file is used to adjust parameters. The py files reads in the parameters, creates the desired sequence, creates a specifically named folder (named by the user) and saves the sequence and the YAML parameters used into that folder. This way, one never loses the parameters required to recreate the specific sequence. Reconstruction will be added later. --- TriangleExample/duyn_triangle.py | 220 +++++++++++++++++++++++ TriangleExample/triangle_parameters.yaml | 56 ++++++ 2 files changed, 276 insertions(+) create mode 100644 TriangleExample/duyn_triangle.py create mode 100644 TriangleExample/triangle_parameters.yaml diff --git a/TriangleExample/duyn_triangle.py b/TriangleExample/duyn_triangle.py new file mode 100644 index 0000000..562d83d --- /dev/null +++ b/TriangleExample/duyn_triangle.py @@ -0,0 +1,220 @@ +import time +import os +import yaml +import numpy as np + +import pypulseq as pp + +# 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_1' + +# 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("Acquisiton 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..917b22a --- /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 : 150 # 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 # Acquisiton delay of the MFC in [sec] +cam_nr_sync_dyn : 0 # no idea (for now) +# \ No newline at end of file From 3e7b9b5ddc91e4eaf8a41886df5296b9af419986 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 24 Feb 2026 12:06:16 +0000 Subject: [PATCH 2/7] [pre-commit] auto fixes from pre-commit hooks --- TriangleExample/duyn_triangle.py | 189 ++++++++++++----------- TriangleExample/triangle_parameters.yaml | 2 +- 2 files changed, 99 insertions(+), 92 deletions(-) diff --git a/TriangleExample/duyn_triangle.py b/TriangleExample/duyn_triangle.py index 562d83d..bde0ca5 100644 --- a/TriangleExample/duyn_triangle.py +++ b/TriangleExample/duyn_triangle.py @@ -1,9 +1,9 @@ -import time import os -import yaml -import numpy as np +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) @@ -15,7 +15,7 @@ yaml_dir = 'TriangleExample/triangle_parameters.yaml' # Add a custom name to the file to differentiate between measurements -add = 'Versuch_1' +add = 'Versuch_1' # 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) @@ -24,25 +24,25 @@ # 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)) + # 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 +# 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", + 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'], @@ -51,14 +51,20 @@ 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}") +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, + delay=system.rf_dead_time, duration=parameters['RF_DUR'], slice_thickness=parameters['SLICE_THICKNESS'], apodization=parameters['apodization'], @@ -69,20 +75,24 @@ # Necessary delays such as TR, EddyCurrentComp, Trig for MFC -grad_free_time = pp.make_delay(parameters['grad_free']) # from skope_gtf.m +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 +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'])) +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 +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 @@ -90,110 +100,117 @@ # 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) + 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']): + 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) - + 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']): + 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) - - + 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)) + 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 + 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 - + 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 - + 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) + 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) + + 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 + 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 - + + 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 - + 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 +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() +# Step 8 - Timing Check (Also good for sanity) +ok, error_report = seq.check_timing() if ok: - print("\nTiming check passed successfully") + print('\nTiming check passed successfully') else: - print("\nTiming check failed! Error listing follows\n") + 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("Acquisiton Method: ", method) +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); +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('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' +filename = f'{output_path}/{time.strftime("%Y%m%d")}_Triangles' print(f"\nSaving sequence file '{filename}.seq' in 'output' folder.") @@ -203,18 +220,8 @@ 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) - +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 index 917b22a..3eab9a2 100644 --- a/TriangleExample/triangle_parameters.yaml +++ b/TriangleExample/triangle_parameters.yaml @@ -51,6 +51,6 @@ grad_free : 0.0005 # Duration of gradient-free period after trigger signal in [s 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 # Acquisiton delay of the MFC in [sec] +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 From 74f3e64c72358459d42203e583d32b058c23a67c Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Fri, 13 Mar 2026 19:21:50 +0000 Subject: [PATCH 3/7] reformat to match mrseq style --- TriangleExample/duyn_triangle.py | 2 +- TriangleExample/triangle_parameters.yaml | 2 +- examples/girf_dyn_triangle.ipynb | 155 +++++++++++ src/mrseq/scripts/girf_dyn_triangle.py | 330 +++++++++++++++++++++++ tests/scripts/test_girf_dyn_triangle.py | 13 + 5 files changed, 500 insertions(+), 2 deletions(-) create mode 100644 examples/girf_dyn_triangle.ipynb create mode 100644 src/mrseq/scripts/girf_dyn_triangle.py create mode 100644 tests/scripts/test_girf_dyn_triangle.py diff --git a/TriangleExample/duyn_triangle.py b/TriangleExample/duyn_triangle.py index bde0ca5..7b46913 100644 --- a/TriangleExample/duyn_triangle.py +++ b/TriangleExample/duyn_triangle.py @@ -15,7 +15,7 @@ yaml_dir = 'TriangleExample/triangle_parameters.yaml' # Add a custom name to the file to differentiate between measurements -add = 'Versuch_1' +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) diff --git a/TriangleExample/triangle_parameters.yaml b/TriangleExample/triangle_parameters.yaml index 3eab9a2..3ad27f7 100644 --- a/TriangleExample/triangle_parameters.yaml +++ b/TriangleExample/triangle_parameters.yaml @@ -3,7 +3,7 @@ # Gradient Raster Time: dwell_time : 0.00001 # Maximum Slew Rate: -slew_rate : 150 # in [T/m/s] +slew_rate : 120 # in [T/m/s] # Maximum Gradient Amplitude in [T/m] max_grad : 0.03 # diff --git a/examples/girf_dyn_triangle.ipynb b/examples/girf_dyn_triangle.ipynb new file mode 100644 index 0000000..71e91e5 --- /dev/null +++ b/examples/girf_dyn_triangle.ipynb @@ -0,0 +1,155 @@ +{ + "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 pathlib import Path\n", + "\n", + "import MRzeroCore as mr0\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)" + ] + }, + { + "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=1e0)\n", + "mr0.sig_to_mrd(fname_mrd, signal, sequence)" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "### Estimate GIRF" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "kdata = KData.from_file(fname_mrd, trajectory=KTrajectoryCartesian())" + ] + } + ], + "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..9a7b38d --- /dev/null +++ b/src/mrseq/scripts/girf_dyn_triangle.py @@ -0,0 +1,330 @@ +"""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_phi: float, + rf_dur: float, + rf_bwt: float, + apodization: float, + n_avg: int, + adc_duration: float, + dwell_time: float, + slice_thickness: float, + slice_pos: Sequence[float], + enumerate_coeff: Sequence[float], + tr: float, + ec_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_phi + Flip angle of RF excitation pulse (in degrees). + rf_dur + 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). + enumerate_coeff + List of amplitude coefficients for gradient encoding. + tr + Repetition time (in seconds). + ec_delay + Eddy current compensation delay (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, g_sl, g_sl_r = pp.make_sinc_pulse( + flip_angle=rf_phi / 180 * np.pi, + delay=system.rf_dead_time, + duration=rf_dur, + 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(ec_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(enumerate_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='SET', 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(enumerate_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 = g_sl.amplitude * slice_pos_val + g_sl.channel = grad_channel + g_sl_r.channel = grad_channel + + # Add RF pulse with slice selection + seq.add_block(rf, g_sl, avg_label, seg_label, grad_label, slice_label, amp_fac_label) + seq.add_block(g_sl_r) + + # Add eddy current compensation delay + seq.add_block(delay_ec) + + # Add ADC and gradient triangle if amplitude is non-zero + if amp_fac != 0: + # 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(trig) + seq.add_block(grad_free_time) + seq.add_block(adc, g_triangle) + else: + seq.add_block(trig) + seq.add_block(grad_free_time) + seq.add_block(adc) + + # 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('TR', tr) + seq.set_definition('RiseTimes', rise_times) + seq.set_definition('GradChannels', grad_channels) + seq.set_definition('SlicePositions', slice_pos) + seq.set_definition('Averages', n_avg) + 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('TE', 0) + + return seq + + +def main( + system: pp.Opts | None = None, + n_avg: int = 3, + tr: float = 1.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 GIRF sequence with triangular gradients. + + Parameters + ---------- + system + PyPulseq system limits object. Uses default system if None. + n_avg + Number of averages. Default: 3 + tr + Repetition time (in seconds). Default: 1.0 + slice_thickness + Slice thickness (in meters). Default: 1.5e-3 + slice_pos + List of slice positions (in meters). Default: [0.04] + show_plots + Toggles sequence plot. Default: True + test_report + Toggles advanced test report. Default: True + timing_check + Toggles timing check of the sequence. Default: True + v141_compatibility + Save the sequence in pulseq v1.4.1 for backwards compatibility. Default: True + + 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_phi = 90 # flip angle [degrees] + rf_dur = 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 + enumerate_coeff = [-1.0, 0.0] # amplitude coefficients for gradient encoding + ec_delay = 1e-3 # 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_phi=rf_phi, + rf_dur=rf_dur, + 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, + enumerate_coeff=enumerate_coeff, + tr=tr, + ec_delay=ec_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..4672f7c --- /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 = 231.42456 # defined 2026-03-13 + + +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) From a135885ca512c0568c17eac8ef7f8bb37407b999 Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Thu, 19 Mar 2026 13:13:35 +0000 Subject: [PATCH 4/7] add references and simplify --- src/mrseq/scripts/girf_dyn_triangle.py | 83 +++++++++++++------------ tests/scripts/test_girf_dyn_triangle.py | 2 +- 2 files changed, 44 insertions(+), 41 deletions(-) diff --git a/src/mrseq/scripts/girf_dyn_triangle.py b/src/mrseq/scripts/girf_dyn_triangle.py index 9a7b38d..ec37a18 100644 --- a/src/mrseq/scripts/girf_dyn_triangle.py +++ b/src/mrseq/scripts/girf_dyn_triangle.py @@ -23,7 +23,7 @@ def girf_triangle_kernel( slice_pos: Sequence[float], enumerate_coeff: Sequence[float], tr: float, - ec_delay: float, + adc_delay: float, rise_times: Sequence[float], grad_free: float, trig_duration: float, @@ -60,8 +60,8 @@ def girf_triangle_kernel( List of amplitude coefficients for gradient encoding. tr Repetition time (in seconds). - ec_delay - Eddy current compensation delay (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 @@ -106,7 +106,7 @@ def girf_triangle_kernel( # Create necessary delays grad_free_time = pp.make_delay(grad_free) delay_tr = pp.make_delay(tr) - delay_ec = pp.make_delay(ec_delay) + 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) @@ -157,25 +157,19 @@ def girf_triangle_kernel( # Add eddy current compensation delay seq.add_block(delay_ec) - # Add ADC and gradient triangle if amplitude is non-zero - if amp_fac != 0: - # 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(trig) - seq.add_block(grad_free_time) - seq.add_block(adc, g_triangle) - else: - seq.add_block(trig) - seq.add_block(grad_free_time) - seq.add_block(adc) + # 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(trig) + seq.add_block(grad_free_time) + seq.add_block(adc, g_triangle) # Add TR delay seq.add_block(delay_tr) @@ -186,26 +180,23 @@ def girf_triangle_kernel( 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('GradChannels', grad_channels) - seq.set_definition('SlicePositions', slice_pos) - seq.set_definition('Averages', n_avg) 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('TE', 0) - return seq def main( system: pp.Opts | None = None, n_avg: int = 3, - tr: float = 1.0, + tr: float = 2.0, slice_thickness: float = 1.5e-3, slice_pos: list[float] | None = None, show_plots: bool = True, @@ -213,28 +204,40 @@ def main( timing_check: bool = True, v141_compatibility: bool = True, ) -> tuple[pp.Sequence, Path]: - """Generate a GIRF sequence with triangular gradients. + """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. Default: 3 + Number of averages. tr - Repetition time (in seconds). Default: 1.0 + Repetition time (in seconds). slice_thickness - Slice thickness (in meters). Default: 1.5e-3 + Slice thickness (in meters). slice_pos - List of slice positions (in meters). Default: [0.04] + List of slice positions (in meters). show_plots - Toggles sequence plot. Default: True + Toggles sequence plot. test_report - Toggles advanced test report. Default: True + Toggles advanced test report. timing_check - Toggles timing check of the sequence. Default: True + Toggles timing check of the sequence. v141_compatibility - Save the sequence in pulseq v1.4.1 for backwards compatibility. Default: True + Save the sequence in pulseq v1.4.1 for backwards compatibility. Returns ------- @@ -262,7 +265,7 @@ def main( # Define sequence parameters enumerate_coeff = [-1.0, 0.0] # amplitude coefficients for gradient encoding - ec_delay = 1e-3 # eddy current compensation delay [s] + adc_delay = 1e-3 # 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 @@ -292,7 +295,7 @@ def main( slice_pos=slice_pos, enumerate_coeff=enumerate_coeff, tr=tr, - ec_delay=ec_delay, + adc_delay=adc_delay, rise_times=rise_times, grad_free=grad_free, trig_duration=trig_duration, diff --git a/tests/scripts/test_girf_dyn_triangle.py b/tests/scripts/test_girf_dyn_triangle.py index 4672f7c..2b7d6be 100644 --- a/tests/scripts/test_girf_dyn_triangle.py +++ b/tests/scripts/test_girf_dyn_triangle.py @@ -3,7 +3,7 @@ import pytest from mrseq.scripts.girf_dyn_triangle import main as create_seq -EXPECTED_DUR = 231.42456 # defined 2026-03-13 +EXPECTED_DUR = 447.42456 # defined 2026-03-19 def test_default_seq_duration(system_defaults): From 1d47433e626d132ec4e7d5a6b3f1d4e34007a1ef Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Thu, 19 Mar 2026 13:17:06 +0000 Subject: [PATCH 5/7] adapt naming --- src/mrseq/scripts/girf_dyn_triangle.py | 32 +++++++++++++------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/mrseq/scripts/girf_dyn_triangle.py b/src/mrseq/scripts/girf_dyn_triangle.py index ec37a18..f342204 100644 --- a/src/mrseq/scripts/girf_dyn_triangle.py +++ b/src/mrseq/scripts/girf_dyn_triangle.py @@ -12,8 +12,8 @@ def girf_triangle_kernel( system: pp.Opts, - rf_phi: float, - rf_dur: float, + rf_flip_angle: float, + rf_duration: float, rf_bwt: float, apodization: float, n_avg: int, @@ -38,9 +38,9 @@ def girf_triangle_kernel( ---------- system PyPulseq system limits object. - rf_phi + rf_flip_angle Flip angle of RF excitation pulse (in degrees). - rf_dur + rf_duration Duration of RF excitation pulse (in seconds). rf_bwt Bandwidth-time product of RF excitation pulse (Hz * seconds). @@ -92,10 +92,10 @@ def girf_triangle_kernel( seq = pp.Sequence(system=system) # Create and set sinc pulse parameters - rf, g_sl, g_sl_r = pp.make_sinc_pulse( - flip_angle=rf_phi / 180 * np.pi, + rf, gz, gzr = pp.make_sinc_pulse( + flip_angle=rf_flip_angle / 180 * np.pi, delay=system.rf_dead_time, - duration=rf_dur, + duration=rf_duration, slice_thickness=slice_thickness, apodization=apodization, time_bw_product=rf_bwt, @@ -146,13 +146,13 @@ def girf_triangle_kernel( amp = amp_fac * system.max_slew * rise_time_val # Set RF frequency offset for current slice - rf.freq_offset = g_sl.amplitude * slice_pos_val - g_sl.channel = grad_channel - g_sl_r.channel = grad_channel + 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, g_sl, avg_label, seg_label, grad_label, slice_label, amp_fac_label) - seq.add_block(g_sl_r) + 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) @@ -254,8 +254,8 @@ def main( slice_pos = [0.04] # Define RF excitation pulse parameters - rf_phi = 90 # flip angle [degrees] - rf_dur = 8.4e-3 # duration of the RF excitation pulse [s] + 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 @@ -284,8 +284,8 @@ def main( seq = girf_triangle_kernel( system=system, - rf_phi=rf_phi, - rf_dur=rf_dur, + rf_flip_angle=rf_flip_angle, + rf_duration=rf_duration, rf_bwt=rf_bwt, apodization=apodization, n_avg=n_avg, From 7b633e0e6c22dfc056de6deafbee43c3761c4866 Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Thu, 19 Mar 2026 13:20:55 +0000 Subject: [PATCH 6/7] reorder --- src/mrseq/scripts/girf_dyn_triangle.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/mrseq/scripts/girf_dyn_triangle.py b/src/mrseq/scripts/girf_dyn_triangle.py index f342204..8aba113 100644 --- a/src/mrseq/scripts/girf_dyn_triangle.py +++ b/src/mrseq/scripts/girf_dyn_triangle.py @@ -156,6 +156,8 @@ def girf_triangle_kernel( # 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( @@ -166,9 +168,6 @@ def girf_triangle_kernel( amplitude=amp, delay=0, ) - - seq.add_block(trig) - seq.add_block(grad_free_time) seq.add_block(adc, g_triangle) # Add TR delay From 44ad32283c74aefd18d727a4e06b3ab1a544c9f8 Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Tue, 19 May 2026 18:28:16 +0100 Subject: [PATCH 7/7] example finished --- examples/girf_dyn_triangle.ipynb | 352 ++++++++++++++++++++++++- src/mrseq/scripts/girf_dyn_triangle.py | 22 +- 2 files changed, 362 insertions(+), 12 deletions(-) diff --git a/examples/girf_dyn_triangle.ipynb b/examples/girf_dyn_triangle.ipynb index 71e91e5..64730b1 100644 --- a/examples/girf_dyn_triangle.ipynb +++ b/examples/girf_dyn_triangle.ipynb @@ -26,9 +26,13 @@ "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", @@ -75,7 +79,11 @@ "metadata": {}, "outputs": [], "source": [ - "phantom = mr0.util.load_phantom(image_matrix_size)" + "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" ] }, { @@ -120,7 +128,7 @@ "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=1e0)\n", + "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e-1)\n", "mr0.sig_to_mrd(fname_mrd, signal, sequence)" ] }, @@ -132,15 +140,351 @@ "### Estimate GIRF" ] }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "#### Sort data and compress coils" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "12", + "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": [ - "kdata = KData.from_file(fname_mrd, trajectory=KTrajectoryCartesian())" + "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": { diff --git a/src/mrseq/scripts/girf_dyn_triangle.py b/src/mrseq/scripts/girf_dyn_triangle.py index 8aba113..e66e095 100644 --- a/src/mrseq/scripts/girf_dyn_triangle.py +++ b/src/mrseq/scripts/girf_dyn_triangle.py @@ -21,7 +21,7 @@ def girf_triangle_kernel( dwell_time: float, slice_thickness: float, slice_pos: Sequence[float], - enumerate_coeff: Sequence[float], + g_amplitude_coeff: Sequence[float], tr: float, adc_delay: float, rise_times: Sequence[float], @@ -56,7 +56,7 @@ def girf_triangle_kernel( Slice thickness (in meters). slice_pos List of slice positions (in meters). - enumerate_coeff + g_amplitude_coeff List of amplitude coefficients for gradient encoding. tr Repetition time (in seconds). @@ -122,7 +122,7 @@ def girf_triangle_kernel( # Calculate total number of dynamics for camera grad_channels = ['x', 'y', 'z'] - cam_nr_dynamics = len(enumerate_coeff) * len(slice_pos) * len(grad_channels) * len(rise_times) * n_avg + 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 @@ -134,12 +134,12 @@ def girf_triangle_kernel( 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='SET', value=idx_grad) + 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(enumerate_coeff): + 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 @@ -189,6 +189,12 @@ def girf_triangle_kernel( 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 @@ -263,8 +269,8 @@ def main( dwell_time = 10e-6 # ADC dwell time [s] # Define sequence parameters - enumerate_coeff = [-1.0, 0.0] # amplitude coefficients for gradient encoding - adc_delay = 1e-3 # eddy current compensation delay [s] + 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 @@ -292,7 +298,7 @@ def main( dwell_time=dwell_time, slice_thickness=slice_thickness, slice_pos=slice_pos, - enumerate_coeff=enumerate_coeff, + g_amplitude_coeff=g_amplitude_coeff, tr=tr, adc_delay=adc_delay, rise_times=rise_times,