From 5354e68e4c7ec92f8efe6fef83e258cb8124ea27 Mon Sep 17 00:00:00 2001 From: OErbil Date: Mon, 23 Feb 2026 13:48:00 +0100 Subject: [PATCH 1/3] 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/3] [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/3] 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)