From 5354e68e4c7ec92f8efe6fef83e258cb8124ea27 Mon Sep 17 00:00:00 2001 From: OErbil Date: Mon, 23 Feb 2026 13:48:00 +0100 Subject: [PATCH 01/13] 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 02/13] [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 1c6eb07e41ece7bf793c888276034f11e9413ec9 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sat, 28 Feb 2026 19:49:14 +0000 Subject: [PATCH 03/13] Bump the all-actions group with 2 updates (#71) Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/deployment.yml | 4 ++-- .github/workflows/examples.yml | 4 ++-- .github/workflows/pytest.yml | 2 +- .github/workflows/report_pytest.yml | 4 ++-- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/.github/workflows/deployment.yml b/.github/workflows/deployment.yml index ed460e6..6ea41a1 100644 --- a/.github/workflows/deployment.yml +++ b/.github/workflows/deployment.yml @@ -66,7 +66,7 @@ jobs: python -m build - name: Store PyPI Distribution - uses: actions/upload-artifact@v6 + uses: actions/upload-artifact@v7 with: name: pypi-package-distribution path: dist/ @@ -88,7 +88,7 @@ jobs: steps: - name: Download PyPI Distribution - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: name: pypi-package-distribution path: dist/ diff --git a/.github/workflows/examples.yml b/.github/workflows/examples.yml index 28105cf..9efafdf 100644 --- a/.github/workflows/examples.yml +++ b/.github/workflows/examples.yml @@ -33,7 +33,7 @@ jobs: # Upload the book's HTML as an artifact (useful for PR previews) - name: Upload built HTML - uses: actions/upload-artifact@v6 + uses: actions/upload-artifact@v7 with: name: built-book path: examples/_build/html @@ -48,7 +48,7 @@ jobs: steps: - name: Download built HTML - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: name: built-book path: examples/_build/html diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 40d8e61..0e96dfe 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -35,7 +35,7 @@ jobs: - name: Upload PyTest Coverage Report if: always() - uses: actions/upload-artifact@v6 + uses: actions/upload-artifact@v7 with: name: pytest-report-${{ matrix.python-version }} path: | diff --git a/.github/workflows/report_pytest.yml b/.github/workflows/report_pytest.yml index 11f00d2..29c8a7e 100644 --- a/.github/workflows/report_pytest.yml +++ b/.github/workflows/report_pytest.yml @@ -21,7 +21,7 @@ jobs: steps: # use the results of python 3.12, consider this as target platform - name: Download PyTest report artifact for Python 3.12 - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: name: pytest-report-3.12 run-id: ${{ github.event.workflow_run.id }} @@ -82,7 +82,7 @@ jobs: actions: read # Required for downloading artifacts steps: - name: Download artifacts - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: run-id: ${{ github.event.workflow_run.id }} github-token: ${{ secrets.GITHUB_TOKEN }} From 07b0c079595f84629b68ca88322400cb52cae86d Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Thu, 5 Mar 2026 14:14:22 +0000 Subject: [PATCH 04/13] Radial Phase Encoding with multi-echo readout (#54) Co-authored-by: Patrick Schuenke <37338697+schuenke@users.noreply.github.com> Co-authored-by: Patrick Schuenke --- examples/_toc.yml | 1 + examples/grpe_flash_dixon.ipynb | 434 ++++++++++++++++++ examples/imaging.md | 3 +- examples/radial_flash.ipynb | 2 +- src/mrseq/scripts/grpe_flash_dixon.py | 602 +++++++++++++++++++++++++ src/mrseq/utils/__init__.py | 2 +- src/mrseq/utils/trajectory.py | 54 ++- tests/scripts/test_grpe_flash_dixon.py | 48 ++ tests/utils/test_trajectory.py | 12 +- 9 files changed, 1124 insertions(+), 34 deletions(-) create mode 100644 examples/grpe_flash_dixon.ipynb create mode 100644 src/mrseq/scripts/grpe_flash_dixon.py create mode 100644 tests/scripts/test_grpe_flash_dixon.py diff --git a/examples/_toc.yml b/examples/_toc.yml index 9f362f2..4b4992e 100644 --- a/examples/_toc.yml +++ b/examples/_toc.yml @@ -5,6 +5,7 @@ chapters: sections: - file: radial_flash - file: spiral_flash + - file: grpe_flash_dixon - file: relaxometry.md sections: - file: t1_inv_rec_gre_single_line diff --git a/examples/grpe_flash_dixon.ipynb b/examples/grpe_flash_dixon.ipynb new file mode 100644 index 0000000..6a49333 --- /dev/null +++ b/examples/grpe_flash_dixon.ipynb @@ -0,0 +1,434 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# GRPE FLASH Dixon\n", + "\n", + "Golden radial phase encoding acquisition with a 3-point Dixon gradient echo readout. More detailed information on the acquisition can be found here:\n", + "\n", + "Mayer et al 2022. Cardio‐respiratory motion‐corrected 3D cardiac water‐fat MRI using model‐based image reconstruction. MRM. https://doi.org/10.1002/mrm.29284" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import tempfile\n", + "import time\n", + "from pathlib import Path\n", + "from urllib.request import urlretrieve\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import MRzeroCore as mr0\n", + "import numpy as np\n", + "import scipy as sp\n", + "import torch\n", + "from einops import rearrange\n", + "from mrpro.algorithms.reconstruction import DirectReconstruction\n", + "from mrpro.data import KData\n", + "from mrpro.data.enums import AcqFlags\n", + "from mrpro.data.traj_calculators import KTrajectoryIsmrmrd\n", + "from mrpro.operators.models import SpoiledGRE\n", + "from scipy.interpolate import interp1d\n", + "\n", + "from mrseq.scripts.grpe_flash_dixon import main as create_seq\n", + "from mrseq.utils import combine_ismrmrd_files\n", + "from mrseq.utils import sys_defaults" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Settings\n", + "\n", + "We are going to use a small numerical phantom with a matrix size of 32 x 32 x 32 to reduce run times.\n", + "To ensure we really acquire all data in the steady state, we use a large number of dummy excitations before the actual \n", + "image acquisition. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "# Acquisition parameters\n", + "flip_angle_degree = 12\n", + "n_dummy_spokes = 10\n", + "image_matrix_size = [32, 32, 32] # x,y,z\n", + "\n", + "# Output path\n", + "tmp = tempfile.TemporaryDirectory()\n", + "fname_mrd = Path(tmp.name) / 'grpe_flash_dixon.mrd'" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### Digital Phantom\n", + "\n", + "We use the Brainweb phantom from [MRzero](https://github.com/MRsources/MRzero-Core), interpolated to our desired matrix size, with uniform B1 field." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "# Load 3D brainweb phantom\n", + "urlretrieve(\n", + " 'https://github.com/MRsources/MRzero-Core/raw/main/documentation/playground_mr0/subject05.npz', 'subject05.npz'\n", + ")\n", + "\n", + "obj_p = mr0.VoxelGridPhantom.brainweb('subject05.npz')\n", + "phantom = obj_p.interpolate(*image_matrix_size)\n", + "phantom.size[2] = phantom.size[1] # GRPE assumes isotropic voxels in y-z plane\n", + "phantom.B1[:] = 1.0" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "### Create the 3D GRPE FLASH Sequence\n", + "\n", + "Generate a golden radial phase encoding FLASH sequence with partial Fourier along the RPE lines and partial echo along the readout." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "sequence, fname_seq = create_seq(\n", + " system=sys_defaults,\n", + " test_report=False,\n", + " timing_check=False,\n", + " rf_flip_angle=flip_angle_degree,\n", + " fov_x=float(phantom.size.numpy()[0]),\n", + " fov_y=float(phantom.size.numpy()[1]),\n", + " fov_z=float(phantom.size.numpy()[1]),\n", + " n_readout=32,\n", + " n_rpe_points=32,\n", + " n_rpe_points_per_shot=4,\n", + " n_rpe_spokes=48,\n", + " partial_echo_factor=0.7,\n", + " partial_fourier_factor=0.7,\n", + " n_dummy_spokes=n_dummy_spokes,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### Simulate the Sequence\n", + "\n", + "Pass the sequence and phantom to the MRzero simulation engine and save the signal as an ISMRMRD file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))\n", + "tstart = time.time()\n", + "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e0)\n", + "print(f'Simulation time: {(time.time() - tstart) / 60:.2f} min')\n", + "mr0.sig_to_mrd(fname_mrd, signal, sequence)\n", + "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "### Gradient correction\n", + "\n", + "We are using a bipolar readout gradient to achieve short echo times and efficient sampling. Due to hardware limitations, \n", + "the positive and negative gradient lobes are not exactly the same. This can lead to small shifts of the k-space \n", + "trajectory between positive and negative readout gradients. If these are not corrected for, there are linear phase\n", + "errors in image space along the readout direction, which can strongly impair the fat-water separation.\n", + "\n", + "There are several approaches that estimate this phase error from the imaging data. For this sequence, we acquire\n", + "additional data without phase encoding, where we change the polarity of the readout gradient. First, we acquire with \n", + "positive - negative - positive readout gradients and then negative - positive - negative readout gradients. We can then \n", + "use a cross-correlation between the pairs of positive and negative readout gradients at the same echo times to \n", + "estimate this shift and correct for it. " + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "As a first step, we get the additional data that is labeled as ACQ_IS_PHASECORR_DATA." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "kdata_corr = KData.from_file(\n", + " str(fname_mrd).replace('.mrd', '_with_traj.mrd'),\n", + " trajectory=KTrajectoryIsmrmrd(),\n", + " acquisition_filter_criterion=lambda acquisition: AcqFlags.ACQ_IS_PHASECORR_DATA.value & acquisition.flags,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "To distinguish between the \"positive - negative - positive\" and \"negative - positive - negative\" acquisition, the first \n", + "are labeled with the acquisition index repetition = 0 and the second with repetition = 1. We also use multiple averages \n", + "to reduce noise. So let's split it up into the different components. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "sort_indices = np.lexsort(\n", + " (\n", + " kdata_corr.header.acq_info.idx.contrast.squeeze(),\n", + " kdata_corr.header.acq_info.idx.repetition.squeeze(),\n", + " kdata_corr.header.acq_info.idx.average.squeeze(),\n", + " )\n", + ")\n", + "kdata_corr = kdata_corr[sort_indices.tolist(), ...].rearrange(\n", + " '(average repetition contrast)... -> average repetition contrast ...',\n", + " average=kdata_corr.header.acq_info.idx.average.unique().numel(),\n", + " repetition=kdata_corr.header.acq_info.idx.repetition.unique().numel(),\n", + " contrast=kdata_corr.header.acq_info.idx.contrast.unique().numel(),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "We are going to use the first echo for the estimation of the shift. The code below will also work for the other echoes\n", + "and should lead to the same shift." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "echo = 0\n", + "\n", + "# Get data for the current echo and compute mean over averages\n", + "kdata = torch.abs(torch.mean(kdata_corr.data[:, :, echo], dim=0)).squeeze().clone()\n", + "ktraj = kdata_corr.traj.kx[0, :, echo].squeeze().clone()\n", + "\n", + "if echo % 2 == 1: # switch positive and negative lobe to get the correct sign for the shift\n", + " kdata = kdata[(1, 0), ...]\n", + " ktraj = ktraj[(1, 0), ...]\n", + "\n", + "kspace_start = int(ktraj[0, :].min().abs().item())\n", + "ktraj_interpolated = np.linspace(-kspace_start, kspace_start - 1, 20 * kspace_start)\n", + "kdata_interp = [\n", + " interp1d(pos.numpy(), signal.numpy(), kind='linear', fill_value='extrapolate')(ktraj_interpolated)\n", + " for pos, signal in zip(ktraj, kdata, strict=True)\n", + "]\n", + "\n", + "plt.figure()\n", + "plt.plot(ktraj_interpolated, kdata_interp[0], label='positive gradient')\n", + "plt.plot(ktraj_interpolated, kdata_interp[1], label='negative gradient')\n", + "plt.legend()\n", + "\n", + "cross_correlation = sp.signal.correlate(kdata_interp[0], kdata_interp[1], mode='same')\n", + "kspace_shift = (np.argmax(cross_correlation) - len(kdata_interp[0]) // 2) / 10 # divide by interpolation factor\n", + "print(f'K-space shift (in k-space samples): {kspace_shift}')\n", + "assert kspace_shift == 0" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "### Reconstruction\n", + "\n", + "Reconstruct the 3D image from the k-space data using direct reconstruction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "# Load k-space data and reconstruct\n", + "kdata = KData.from_file(\n", + " str(fname_mrd).replace('.mrd', '_with_traj.mrd'),\n", + " trajectory=KTrajectoryIsmrmrd(),\n", + ")\n", + "\n", + "kdata.traj.kx[1] += kspace_shift # delta_k (set to actual value if needed)\n", + "\n", + "recon = DirectReconstruction(kdata, csm=None)\n", + "idata = recon(kdata)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(4, 3, figsize=(12, 10))\n", + "for cax in ax.flatten():\n", + " cax.set_xticks([])\n", + " cax.set_yticks([])\n", + "\n", + "idat = idata.data.numpy().squeeze()\n", + "z_mid, x_mid = idat.shape[-3] // 2, idat.shape[-1] // 2\n", + "\n", + "for echo in range(idat.shape[0]):\n", + " ax[0, echo].imshow(np.abs(idat[echo, :, :, x_mid]), cmap='gray')\n", + " ax[1, echo].imshow(np.abs(idat[echo, z_mid, :, :]), cmap='gray')\n", + " ax[2, echo].imshow(np.angle(idat[echo, :, :, x_mid]), cmap='bwr', vmin=-np.pi, vmax=np.pi)\n", + " ax[3, echo].imshow(np.angle(idat[echo, z_mid, :, :]), cmap='bwr', vmin=-np.pi, vmax=np.pi)\n", + " ax[0, echo].set_title(f'Echo {echo}')" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "### Compare to Theoretical Model\n", + "\n", + "Compare the reconstructed images to a theoretical spoiled GRE signal model. We calculate $T2^*$ using $1/T2^* = 1/T2 + 1/T2'$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate T2* and generate theoretical model\n", + "t2star = 1 / (1 / phantom.T2 + 1 / phantom.T2dash)\n", + "\n", + "model = SpoiledGRE(\n", + " flip_angle=np.deg2rad(flip_angle_degree),\n", + " echo_time=kdata.header.te,\n", + " repetition_time=kdata.header.tr,\n", + ")\n", + "idat_model = model(m0=phantom.PD, t1=phantom.T1, t2star=t2star)[0]\n", + "idat_model = idat_model.detach().abs().numpy().squeeze()\n", + "idat_model /= idat_model.max()\n", + "idat_model = np.roll(\n", + " rearrange(idat_model[:, ::-1, ::-1, ::-1], 'echo x y z -> echo z y x'),\n", + " shift=(1, 1, 1),\n", + " axis=(-1, -2, -3),\n", + ")\n", + "\n", + "# Create object mask\n", + "obj_mask = np.zeros_like(idat_model)\n", + "obj_mask[idat_model > 0.01] = 1\n", + "\n", + "# Normalize reconstructed image\n", + "idat = idata.data.abs().squeeze().numpy()\n", + "idat = idat * obj_mask\n", + "idat /= np.percentile(idat, 99.9)\n", + "\n", + "idat_diff = idat - idat_model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "# Comparison: axial view (middle z-slice)\n", + "fig, ax = plt.subplots(3, 3, figsize=(12, 10))\n", + "x_mid = idat.shape[-1] // 2\n", + "\n", + "for echo in range(3):\n", + " im = ax[0, echo].imshow(idat[echo, :, :, x_mid], cmap='grey', vmin=0, vmax=1)\n", + " ax[0, echo].set_title(f'Echo {echo} - Reconstructed')\n", + " ax[0, echo].set_xticks([])\n", + " ax[0, echo].set_yticks([])\n", + " fig.colorbar(im, ax=ax[0, echo])\n", + "\n", + " im = ax[1, echo].imshow(idat_model[echo, :, :, x_mid], cmap='grey', vmin=0, vmax=1)\n", + " ax[1, echo].set_title(f'Echo {echo} - Model')\n", + " ax[1, echo].set_xticks([])\n", + " ax[1, echo].set_yticks([])\n", + " fig.colorbar(im, ax=ax[1, echo])\n", + "\n", + " im = ax[2, echo].imshow(idat_diff[echo, :, :, x_mid], cmap='bwr', vmin=-0.2, vmax=0.2)\n", + " ax[2, echo].set_title(f'Echo {echo} - Difference')\n", + " ax[2, echo].set_xticks([])\n", + " ax[2, echo].set_yticks([])\n", + " fig.colorbar(im, ax=ax[2, echo])\n", + "\n", + "plt.tight_layout()\n", + "plt.show()\n", + "\n", + "\n", + "relative_error = np.sum(np.abs(idat_diff)) / np.sum(np.abs(idat_model))\n", + "print(f'Relative error: {relative_error:.4f}')\n", + "assert relative_error < 0.09" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mrseq_mrzero", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/imaging.md b/examples/imaging.md index 91c8959..9ade8b1 100644 --- a/examples/imaging.md +++ b/examples/imaging.md @@ -3,4 +3,5 @@ Sequences for qualitative MRI. Available examples are: - [M2D Radial Flash](radial_flash.ipynb) -- [M2D Spiral Flash](spiral_flash.ipynb) \ No newline at end of file +- [M2D Spiral Flash](spiral_flash.ipynb) +- [3D Golden Radial Phase Encoding with 3-point Dixon readout](grpe_flash_dixon.ipynb) \ No newline at end of file diff --git a/examples/radial_flash.ipynb b/examples/radial_flash.ipynb index 775364e..e94c724 100644 --- a/examples/radial_flash.ipynb +++ b/examples/radial_flash.ipynb @@ -49,7 +49,7 @@ "source": [ "### Settings\n", "We are going to use a numerical phantom with a matrix size of 128 x 128. \n", - "To ensure we really acquire all radial lines in the steaady state, we use a large number of dummy excitations before the actual image acquisition. " + "To ensure we really acquire all radial lines in the steady state, we use a large number of dummy excitations before the actual image acquisition. " ] }, { diff --git a/src/mrseq/scripts/grpe_flash_dixon.py b/src/mrseq/scripts/grpe_flash_dixon.py new file mode 100644 index 0000000..2265a91 --- /dev/null +++ b/src/mrseq/scripts/grpe_flash_dixon.py @@ -0,0 +1,602 @@ +"""3D FLASH sequence with radial phase encoding.""" + +from pathlib import Path + +import ismrmrd +import numpy as np +import pypulseq as pp +from pypulseq.rotate import rotate + +from mrseq.utils import MultiEchoAcquisition +from mrseq.utils import find_gx_flat_time_on_adc_raster +from mrseq.utils import round_to_raster +from mrseq.utils import sys_defaults +from mrseq.utils import write_sequence +from mrseq.utils.constants import GOLDEN_ANGLE_HALF_CIRCLE +from mrseq.utils.ismrmrd import Fov +from mrseq.utils.ismrmrd import Limits +from mrseq.utils.ismrmrd import MatrixSize +from mrseq.utils.ismrmrd import create_header + + +def grpe_flash_dixon_kernel( + system: pp.Opts, + te: float | None, + delta_te: float | None, + n_echoes: int, + tr: float | None, + fov_x: float, + fov_y: float, + fov_z: float, + n_readout: int, + n_rpe_points: int, + n_rpe_points_per_shot: int, + n_rpe_spokes: int, + readout_oversampling: float, + partial_echo_factor: float, + partial_fourier_factor: float, + n_dummy_spokes: int, + gx_pre_duration: float, + gx_flat_time: float, + rf_duration: float, + rf_flip_angle: float, + rf_bwt: float, + rf_apodization: float, + rf_spoiling_phase_increment: float, + gx_spoil_duration: float, + gx_spoil_area: float, + mrd_header_file: str | Path | None, +) -> tuple[pp.Sequence, float, float, float]: + """Generate a 3D FLASH sequence with golden radial phase encoding. + + Parameters + ---------- + system + PyPulseq system limits object. + te + Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. + delta_te + Desired echo spacing (in seconds). Minimum echo spacing is used if set to None. + n_echoes + Number of echoes. + tr + Desired repetition time (TR) (in seconds). + fov_x + Field of view in x direction (in meters). + fov_y + Field of view in y direction (in meters). + fov_z + Field of view in z direction (in meters). + n_readout + Number of frequency encoding steps. + n_rpe_points + Number of radial phase encoding points (points along one RPE line). + n_rpe_points_per_shot + Shots are interleaved groups of points along RPE lines. Each shot obtains the k-space center at the cost of a + point of the highest k-space frequency. Fat-saturation pulses are applied prior to each shot. + n_rpe_spokes + Number of radial phase encoding spokes (number of RPE lines). + readout_oversampling + Readout oversampling factor, commonly 2. This reduces aliasing artifacts. + partial_echo_factor + Partial echo factor along readout (between 0.5 and 1). + partial_fourier_factor + Partial Fourier factor along RPE lines (between 0.5 and 1). + n_dummy_spokes + Number of dummy RPE spokes before data acquisition to ensure steady state. + gx_pre_duration + Duration of readout pre-winder gradient (in seconds) + gx_flat_time + Flat time of readout gradient (in seconds) + rf_duration + Duration of the rf excitation pulse (in seconds) + rf_flip_angle + Flip angle of rf excitation pulse (in degrees) + rf_bwt + Bandwidth-time product of rf excitation pulse (Hz * seconds) + rf_apodization + Apodization factor of rf excitation pulse + rf_spoiling_phase_increment + RF spoiling phase increment (in degrees). Set to 0 for no RF spoiling. + gx_spoil_duration + Duration of spoiler gradient (in seconds) + gx_spoil_area + Area of spoiler gradient (in mT/m * s) + mrd_header_file + Filename of the ISMRMRD header file to be created. If None, no header file is created. + + Returns + ------- + seq + PyPulseq Sequence object + min_te + Shortest possible echo time. + min_tr + Shortest possible repetition time. + delta_te + Time between echoes. + + """ + if readout_oversampling < 1: + raise ValueError('Readout oversampling factor must be >= 1.') + + if n_dummy_spokes < 0: + raise ValueError('Number of dummy spokes must be >= 0.') + + spoke_angle = GOLDEN_ANGLE_HALF_CIRCLE + rpe_radial_shift = [0, 0.5, 0.25, 0.75] + + if partial_echo_factor > 1 or partial_echo_factor < 0.5: + raise ValueError('Partial echo factor has to be within 0.5 and 1') + if partial_fourier_factor > 1 or partial_fourier_factor < 0.5: + raise ValueError('Partial Fourier factor has to be within 0.5 and 1') + + # create PyPulseq Sequence object and set system limits + seq = pp.Sequence(system=system) + + # create slab selective excitation pulse and gradients + rf, gz, gzr = pp.make_sinc_pulse( + flip_angle=rf_flip_angle / 180 * np.pi, + duration=rf_duration, + slice_thickness=fov_z, + apodization=rf_apodization, + time_bw_product=rf_bwt, + delay=system.rf_dead_time, + system=system, + return_gz=True, + use='excitation', + ) + + multi_echo_gradient = MultiEchoAcquisition( + system=system, + delta_te=delta_te, + fov=fov_x, + n_readout=n_readout, + readout_oversampling=readout_oversampling, + partial_echo_factor=partial_echo_factor, + gx_flat_time=gx_flat_time, + gx_pre_duration=gx_pre_duration, + ) + + # calculate gradient areas for phase encoding along each RPE line in a low-high order + delta_ky = 1 / fov_y + n_rpe_points_with_partial_fourier = int( + ((n_rpe_points * partial_fourier_factor) // n_rpe_points_per_shot) * n_rpe_points_per_shot + ) + n_shots_per_rpe_spoke = n_rpe_points_with_partial_fourier // n_rpe_points_per_shot + print( + f'Number of phase encoding points {n_rpe_points_with_partial_fourier}', + f'with partial Fourier factor {partial_fourier_factor}', + ) + enc_steps_pe = np.arange(0, n_rpe_points_with_partial_fourier) + phase_areas = (enc_steps_pe - n_rpe_points / 2) * delta_ky + centric_idx = np.argsort(np.abs(phase_areas), kind='stable') + enc_steps_pe = enc_steps_pe[centric_idx] + + # Interleave shots + enc_steps_pe_interleaved = [] + for step in range(n_shots_per_rpe_spoke): + enc_steps_pe_interleaved.append(enc_steps_pe[step::n_shots_per_rpe_spoke]) + enc_steps_pe = np.concatenate(enc_steps_pe_interleaved) + + # create spoiler gradients + gx_spoil = pp.make_trapezoid(channel='x', system=system, area=gx_spoil_area, duration=gx_spoil_duration) + + # calculate minimum echo time + if te is None: + gzr_gx_dur = pp.calc_duration(gzr, multi_echo_gradient._gx_pre) # gzr and gx_pre are applied simultaneously + else: + gzr_gx_dur = pp.calc_duration(gzr) + pp.calc_duration( + multi_echo_gradient._gx_pre + ) # gzr and gx_pre are applied sequentially + + min_te = ( + rf.shape_dur / 2 # time from center to end of RF pulse + + max(rf.ringdown_time, gz.fall_time) # RF ringdown time or gradient fall time + + gzr_gx_dur # slice selection re-phasing gradient and readout pre-winder + + multi_echo_gradient._gx.delay # potential delay of readout gradient + + multi_echo_gradient._gx.rise_time # rise time of readout gradient + + (multi_echo_gradient._n_readout_pre_echo + 0.5) * multi_echo_gradient._adc.dwell + ).item() + + # calculate echo time delay (te_delay) + if te is None: + te_delay = 0.0 + else: + te_delay = round_to_raster(te - min_te, system.block_duration_raster) + if te_delay < 0: + raise ValueError(f'TE must be larger than {min_te * 1000:.3f} ms. Current value is {te * 1000:.3f} ms.') + + # calculate minimum repetition time + min_tr = ( + pp.calc_duration(gz) # rf pulse + + gzr_gx_dur # slice selection re-phasing gradient and readout pre-winder + + pp.calc_duration(multi_echo_gradient._gx) * n_echoes # readout gradient + + pp.calc_duration(multi_echo_gradient._gx_between) * (n_echoes - 1) # readout gradient + + pp.calc_duration(gx_spoil, multi_echo_gradient._gx_post) # gradient spoiler or readout-re-winder + ).item() + + # calculate repetition time delay (tr_delay) + current_min_tr = min_tr + te_delay + if tr is None: + tr_delay = 0.0 + else: + tr_delay = round_to_raster(tr - current_min_tr, system.block_duration_raster) + if not tr_delay >= 0: + raise ValueError( + f'TR must be larger than {current_min_tr * 1000:.3f} ms. Current value is {tr * 1000:.3f} ms.' + ) + + print(f'\nCurrent echo time = {(min_te + te_delay) * 1000:.3f} ms') + print(f'Current repetition time = {(current_min_tr + tr_delay) * 1000:.3f} ms') + + # choose initial rf phase offset + rf_phase = 0.0 + rf_inc = 0.0 + + # create header + if mrd_header_file: + hdr = create_header( + traj_type='other', + encoding_fov=Fov(x=fov_x * readout_oversampling, y=fov_y, z=fov_y), + recon_fov=Fov(x=fov_x, y=fov_y, z=fov_y), + encoding_matrix=MatrixSize(n_x=int(n_readout * readout_oversampling), n_y=n_rpe_points, n_z=n_rpe_points), + recon_matrix=MatrixSize(n_x=n_readout, n_y=n_rpe_points, n_z=n_rpe_points), + dwell_time=multi_echo_gradient._adc.dwell, + slice_limits=Limits(min=0, max=0, center=0), + k1_limits=Limits(min=0, max=n_rpe_points, center=0), + k2_limits=Limits(min=0, max=n_rpe_spokes, center=0), + ) + + # write header to file + prot = ismrmrd.Dataset(mrd_header_file, 'w') + prot.write_xml_header(hdr.toXML('utf-8')) + + # obtain noise samples + seq.add_block(pp.make_label(label='LIN', type='SET', value=0), pp.make_label(label='PAR', type='SET', value=0)) + seq.add_block(multi_echo_gradient._adc, pp.make_label(label='NOISE', type='SET', value=True)) + seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) + seq.add_block(pp.make_delay(system.rf_dead_time)) + + if mrd_header_file: + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=3, number_of_samples=multi_echo_gradient._adc.num_samples) + prot.append_acquisition(acq) + + # Define sequence blocks + for se_index in np.arange(-n_dummy_spokes, n_rpe_spokes): + se_label = pp.make_label(type='SET', label='PAR', value=int(se_index)) + for shot_index in range(n_shots_per_rpe_spoke): + for pe_index in range(n_rpe_points_per_shot): + pe = int(enc_steps_pe[shot_index * n_rpe_points_per_shot + pe_index]) + pe_label = pp.make_label(type='SET', label='LIN', value=pe) + + # set rf pulse properties and add rf pulse block event + if rf_spoiling_phase_increment > 0: + rf.phase_offset = rf_phase / 180 * np.pi + multi_echo_gradient._adc.phase_offset = rf_phase / 180 * np.pi + + # add rf pulse + seq.add_block(rf, gz) + + # update rf pulse properties for the next loop + rf_inc = divmod(rf_inc + rf_spoiling_phase_increment, 360.0)[1] + rf_phase = divmod(rf_phase + rf_inc, 360.0)[1] + + # area of phase encoding gradient + current_phase_area = phase_areas[pe] + # do not shift k-space center for respiratory navigator calculation + if np.abs(current_phase_area) > 1e-5: + pe_shift = rpe_radial_shift[int(np.mod(se_index, len(rpe_radial_shift)))] + else: + pe_shift = 0.0 + + # phase encoding along pe + gy_pre = pp.make_trapezoid( + channel='y', + area=current_phase_area + pe_shift * delta_ky, + duration=pp.calc_duration(multi_echo_gradient._gx_pre), + system=system, + ) + + # calculate rotated phase encoding gradient + rotation_angle_rad = spoke_angle * se_index + gy_pre_rotated = rotate(gy_pre, angle=rotation_angle_rad, axis='x') + gy_pre = gy_pre_rotated[0] + if len(gy_pre_rotated) == 2: + gz_pre = gy_pre_rotated[1] + else: + gz_pre = pp.make_trapezoid( + channel='z', + area=0, + duration=pp.calc_duration(multi_echo_gradient._gx_pre), + system=system, + ) + assert gy_pre.channel == 'y' and gz_pre.channel == 'z' + + # combine slice rephaser and slice encoding gradient + gz_r_pre = pp.make_trapezoid( + channel='z', + area=gz_pre.area + gzr.area, + duration=pp.calc_duration(multi_echo_gradient._gx_pre), + system=system, + ) + + label_contents = [pe_label, se_label] + seq.add_block(multi_echo_gradient._gx_pre, gy_pre, gz_r_pre, *label_contents) + + # add delay due to TE + if te_delay > 0: + seq.add_block(pp.make_delay(te_delay)) + + # add readout gradients and ADCs + if se_index >= 0: + seq, _ = multi_echo_gradient.add_to_seq_without_pre_post_gradient(seq, n_echoes) + else: + # the most accurate way to get the duration of the readout block is to add it to a dummy sequence + seq_dummy = pp.Sequence(system=system) + seq_dummy, _ = multi_echo_gradient.add_to_seq_without_pre_post_gradient(seq_dummy, n_echoes) + readout_duration = sum(seq_dummy.block_durations.values()) + seq.add_block(pp.make_delay(readout_duration)) + + # add re-winder and spoiler gradients + gy_pre.amplitude = -gy_pre.amplitude + gz_pre.amplitude = -gz_pre.amplitude + seq.add_block(gx_spoil, gy_pre, gz_pre) + + # add delay in case TR > min_TR + if tr_delay > 0: + seq.add_block(pp.make_delay(tr_delay)) + + if mrd_header_file and se_index >= 0: + # add acquisitions to metadata + k0_trajectory = np.linspace( + -multi_echo_gradient._n_readout_pre_echo, + multi_echo_gradient._n_readout_post_echo, + multi_echo_gradient._n_readout_with_partial_echo, + ) + grpe_trajectory = np.zeros((multi_echo_gradient._n_readout_with_partial_echo, 3), dtype=np.float32) + + for echo_ in range(n_echoes): + gx_sign = (-1) ** echo_ + grpe_trajectory[:, 0] = k0_trajectory * gx_sign + grpe_trajectory[:, 1] = (pe - n_rpe_points / 2 + pe_shift) * np.cos(rotation_angle_rad) + grpe_trajectory[:, 2] = (pe - n_rpe_points / 2 + pe_shift) * np.sin(rotation_angle_rad) + + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=3, number_of_samples=multi_echo_gradient._adc.num_samples) + acq.traj[:] = grpe_trajectory + prot.append_acquisition(acq) + + # obtain echoes with positive and negative gradient polarity to be able to correct for any differences between + # positive and negative gradient waveforms in the bi-polar readout + seq.add_block(pp.make_label(type='SET', label='NAV', value=True)) + for average_idx in range(4): + average_label = pp.make_label(type='SET', label='AVG', value=average_idx) + for polarity_idx, polarity in enumerate(['positive', 'negative']): + polarity_label = pp.make_label(type='SET', label='REP', value=polarity_idx) + # set rf pulse properties and add rf pulse block event + if rf_spoiling_phase_increment > 0: + rf.phase_offset = rf_phase / 180 * np.pi + multi_echo_gradient._adc.phase_offset = rf_phase / 180 * np.pi + + # add rf pulse + seq.add_block(rf, gz) + + # update rf pulse properties for the next loop + rf_inc = divmod(rf_inc + rf_spoiling_phase_increment, 360.0)[1] + rf_phase = divmod(rf_phase + rf_inc, 360.0)[1] + + seq.add_block( + multi_echo_gradient._gx_pre + if polarity == 'positive' + else pp.scale_grad(multi_echo_gradient._gx_pre, -1), + average_label, + polarity_label, + ) + + # add delay due to TE + if te_delay > 0: + seq.add_block(pp.make_delay(te_delay)) + + # add readout gradients and ADCs + seq, time_to_echoes = multi_echo_gradient.add_to_seq_without_pre_post_gradient(seq, n_echoes, polarity) # type: ignore[arg-type] + + # add spoiler gradients + seq.add_block(gx_spoil) + + # add delay in case TR > min_TR + if tr_delay > 0: + seq.add_block(pp.make_delay(tr_delay)) + + if mrd_header_file: + # add acquisitions to metadata + k0_trajectory = np.linspace( + -multi_echo_gradient._n_readout_pre_echo, + multi_echo_gradient._n_readout_post_echo, + multi_echo_gradient._n_readout_with_partial_echo, + ) + grpe_trajectory = np.zeros((multi_echo_gradient._n_readout_with_partial_echo, 3), dtype=np.float32) + + for echo_ in range(n_echoes): + gx_sign = (-1) ** echo_ if polarity == 'positive' else (-1) ** (echo_ + 1) + grpe_trajectory[:, 0] = k0_trajectory * gx_sign + + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=3, number_of_samples=multi_echo_gradient._adc.num_samples) + acq.traj[:] = grpe_trajectory + prot.append_acquisition(acq) + + seq.add_block(pp.make_label(label='NAV', type='SET', value=False)) + + # close ISMRMRD file + if mrd_header_file: + prot.close() + + delta_te_array = np.diff(time_to_echoes) + return seq, float(min_te), float(min_tr), float(delta_te_array[0]) + + +def main( + system: pp.Opts | None = None, + tr: float | None = None, + rf_flip_angle: float = 12, + fov_x: float = 128e-3, + fov_y: float = 128e-3, + fov_z: float = 128e-3, + n_readout: int = 128, + n_rpe_points: int = 128, + n_rpe_points_per_shot: int = 8, + n_rpe_spokes: int = 16, + partial_echo_factor: float = 0.7, + partial_fourier_factor: float = 0.7, + receiver_bandwidth_per_pixel: float = 1200, # Hz/pixel + n_dummy_spokes: int = 2, + show_plots: bool = True, + test_report: bool = True, + timing_check: bool = True, + v141_compatibility: bool = True, +) -> pp.Sequence: + """Generate a 3D FLASH sequence with radial phase encoding. + + Parameters + ---------- + system + PyPulseq system limits object. + tr + Desired repetition time (TR) (in seconds). Minimum repetition time is used if set to None. + rf_flip_angle + Flip angle of rf excitation pulse (in degrees) + fov_x + Field of view in x direction (in meters). + fov_y + Field of view in y direction (in meters). + fov_z + Field of view in z direction (in meters). + n_readout + Number of frequency encoding steps. + n_rpe_points + Number of radial phase encoding points (points along one RPE line). + n_rpe_points_per_shot + Shots are interleaved groups of points along RPE lines. Each shot obtains the k-space center at the cost of a + point of the highest k-space frequency. Fat-saturation pulses are applied prior to each shot. + n_rpe_spokes + Number of radial phase encoding spokes (number of RPE lines). + partial_echo_factor + Partial echo factor along the readout (between 0.5 and 1). + partial_fourier_factor + Partial Fourier factor along RPE lines (between 0.5 and 1). + receiver_bandwidth_per_pixel + Desired receiver bandwidth per pixel (in Hz/pixel). This is used to calculate the readout duration. + n_dummy_spokes + Number of dummy RPE spokes before data acquisition to ensure steady state + show_plots + Toggles sequence plot. + test_report + Toggles advanced test report. + timing_check + Toggles timing check of the sequence. + v141_compatibility + Save the sequence in pulseq v1.4.1 for backwards compatibility. + """ + if system is None: + system = sys_defaults + + # define settings of rf excitation pulse + rf_duration = 0.6e-3 # duration of the rf excitation pulse [s] + rf_bwt = 2 # bandwidth-time product of rf excitation pulse [Hz*s] + rf_apodization = 0.5 # apodization factor of rf excitation pulse + readout_oversampling = 2 # readout oversampling factor, commonly 2. This reduces aliasing artifacts. + + # this is just approximately, the final calculation is done in the kernel + n_readout_with_oversampling = int(n_readout * readout_oversampling * partial_echo_factor) + # define ADC and gradient timing + adc_dwell_time = 1.0 / (receiver_bandwidth_per_pixel * n_readout_with_oversampling) + gx_pre_duration = 1.0e-3 # duration of readout pre-winder gradient [s] + gx_flat_time, adc_dwell_time = find_gx_flat_time_on_adc_raster( + n_readout_with_oversampling, adc_dwell_time, system.grad_raster_time, system.adc_raster_time + ) + + te = None # shortest possible echo time + delta_te = None # shortest possible delta echo time + n_echoes = 3 + + # define spoiling + gx_spoil_duration = 1.9e-3 # duration of spoiler gradient [s] + gx_spoil_area = readout_oversampling * n_readout * 1 / fov_x # area / zeroth gradient moment of spoiler gradient + rf_spoiling_phase_increment = 117 # RF spoiling phase increment [°]. Set to 0 for no RF spoiling. + + # define sequence filename + filename = f'{Path(__file__).stem}_fov{int(fov_x * 1000)}_{int(fov_y * 1000)}_{int(fov_z * 1000)}mm_' + filename += f'{n_readout}_{n_rpe_points}_{n_rpe_spokes}_3ne' + + output_path = Path.cwd() / 'output' + output_path.mkdir(parents=True, exist_ok=True) + + # delete existing header file + if (output_path / Path(filename + '_header.h5')).exists(): + (output_path / Path(filename + '_header.h5')).unlink() + + seq, min_te, min_tr, delta_te = grpe_flash_dixon_kernel( + system=system, + te=te, + delta_te=delta_te, + n_echoes=n_echoes, + tr=tr, + fov_x=fov_x, + fov_y=fov_y, + fov_z=fov_z, + n_readout=n_readout, + n_rpe_points=n_rpe_points, + n_rpe_points_per_shot=n_rpe_points_per_shot, + n_rpe_spokes=n_rpe_spokes, + readout_oversampling=readout_oversampling, + partial_echo_factor=partial_echo_factor, + partial_fourier_factor=partial_fourier_factor, + n_dummy_spokes=n_dummy_spokes, + gx_pre_duration=gx_pre_duration, + gx_flat_time=gx_flat_time, + rf_duration=rf_duration, + rf_flip_angle=rf_flip_angle, + rf_bwt=rf_bwt, + rf_apodization=rf_apodization, + rf_spoiling_phase_increment=rf_spoiling_phase_increment, + gx_spoil_duration=gx_spoil_duration, + gx_spoil_area=gx_spoil_area, + mrd_header_file=output_path / Path(filename + '_header.h5'), + ) + + # check timing of the sequence + if timing_check and not test_report: + ok, error_report = seq.check_timing() + if ok: + print('\nTiming check passed successfully') + else: + print('\nTiming check failed! Error listing follows\n') + print(error_report) + + # show advanced rest report + if test_report: + print('\nCreating advanced test report...') + print(seq.test_report()) + + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_x, fov_y, fov_z]) + seq.set_definition('ReconMatrix', (n_readout, n_rpe_points, n_rpe_points)) + seq.set_definition('SliceThickness', fov_z) + seq.set_definition('TE', [(te or min_te) + idx * delta_te for idx in range(n_echoes)]) + seq.set_definition('TR', tr or min_tr) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) + + # save seq-file to disk + print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") + write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) + + if show_plots: + seq.plot(time_range=(0, (n_dummy_spokes * n_rpe_points + 20) * (tr or min_tr))) + + return seq, output_path / filename + + +if __name__ == '__main__': + main() diff --git a/src/mrseq/utils/__init__.py b/src/mrseq/utils/__init__.py index 00d413e..31bd0a0 100644 --- a/src/mrseq/utils/__init__.py +++ b/src/mrseq/utils/__init__.py @@ -1,5 +1,5 @@ from mrseq.utils.sequence_helper import find_gx_flat_time_on_adc_raster, round_to_raster, write_sequence from mrseq.utils.system_defaults import sys_defaults +from mrseq.utils.trajectory import cartesian_phase_encoding, spiral_acquisition, MultiEchoAcquisition from mrseq.utils.vds import variable_density_spiral_trajectory -from mrseq.utils.trajectory import spiral_acquisition from mrseq.utils.ismrmrd import create_header, combine_ismrmrd_files diff --git a/src/mrseq/utils/trajectory.py b/src/mrseq/utils/trajectory.py index 91beacd..c40b5ca 100644 --- a/src/mrseq/utils/trajectory.py +++ b/src/mrseq/utils/trajectory.py @@ -7,10 +7,10 @@ import numpy as np import pypulseq as pp -from mrseq.utils import find_gx_flat_time_on_adc_raster -from mrseq.utils import round_to_raster -from mrseq.utils import sys_defaults -from mrseq.utils import variable_density_spiral_trajectory +from mrseq.utils.sequence_helper import find_gx_flat_time_on_adc_raster +from mrseq.utils.sequence_helper import round_to_raster +from mrseq.utils.system_defaults import sys_defaults +from mrseq.utils.vds import variable_density_spiral_trajectory def cartesian_phase_encoding( @@ -63,10 +63,10 @@ def cartesian_phase_encoding( if n_phase_encoding_per_shot and sampling_order != 'random': kpe_extended = np.arange(-n_phase_encoding, n_phase_encoding) kpe_extended = kpe_extended[np.argsort(np.abs(kpe_extended), kind='stable')] - idx = 0 + kidx = 0 while np.mod(len(kpe), n_phase_encoding_per_shot) > 0: - kpe = np.unique(np.concatenate((kpe, (kpe_extended[idx],)))) - idx += 1 + kpe = np.unique(np.concatenate((kpe, (kpe_extended[kidx],)))) + kidx += 1 # Different temporal orders of phase encoding points if sampling_order == 'random': @@ -78,11 +78,11 @@ def cartesian_phase_encoding( elif sampling_order == 'linear': kpe = np.sort(kpe) elif sampling_order == 'low_high': - sort_idx = np.argsort(np.abs(kpe), kind='stable') - kpe = kpe[sort_idx] + idx = np.argsort(np.abs(kpe), kind='stable') + kpe = kpe[idx] elif sampling_order == 'high_low': - sort_idx = np.argsort(-np.abs(kpe), kind='stable') - kpe = kpe[sort_idx] + idx = np.argsort(-np.abs(kpe), kind='stable') + kpe = kpe[idx] else: raise ValueError(f'sampling order {sampling_order} not supported.') @@ -205,16 +205,15 @@ def __init__( min_delta_te = pp.calc_duration(self._gx) + pp.calc_duration(self._gx_between) if delta_te is None: - self._te_delay = 0.0 + self._delta_te_delay = 0.0 else: - self._te_delay = round_to_raster(delta_te - min_delta_te, self._system.block_duration_raster) - if not self._te_delay >= 0: + self._delta_te_delay = round_to_raster(delta_te - min_delta_te, self._system.block_duration_raster) + if self._delta_te_delay < 0: raise ValueError( - f'Delta TE must be larger than {min_delta_te * 1000:.2f} ms. ' - f'Current value is {delta_te * 1000:.2f} ms.' + f'TE must be larger than {min_delta_te * 1000:.3f} ms. Current value is {delta_te * 1000:.3f} ms.' ) - def add_to_seq(self, seq: pp.Sequence, n_echoes: int) -> tuple[pp.Sequence, list[float]]: + def add_to_seq(self, seq: pp.Sequence, n_echoes: int, polarity: Literal['positive', 'negative'] = 'positive'): """Add all gradients and adc to sequence. Parameters @@ -223,6 +222,8 @@ def add_to_seq(self, seq: pp.Sequence, n_echoes: int) -> tuple[pp.Sequence, list PyPulseq Sequence object. n_echoes Number of echoes + polarity + Polarity of first readout gradient Returns ------- @@ -232,17 +233,19 @@ def add_to_seq(self, seq: pp.Sequence, n_echoes: int) -> tuple[pp.Sequence, list Time from beginning of sequence to echoes. """ # readout pre-winder - seq.add_block(self._gx_pre) + seq.add_block(self._gx_pre if polarity == 'positive' else pp.scale_grad(self._gx_pre, -1)) # add readout gradients and ADCs - seq, time_to_echoes = self.add_to_seq_without_pre_post_gradient(seq, n_echoes) + seq, time_to_echoes = self.add_to_seq_without_pre_post_gradient(seq, n_echoes, polarity) # readout re-winder - seq.add_block(self._gx_post) + seq.add_block(self._gx_post if polarity == 'positive' else pp.scale_grad(self._gx_post, -1)) return seq, time_to_echoes - def add_to_seq_without_pre_post_gradient(self, seq: pp.Sequence, n_echoes: int) -> tuple[pp.Sequence, list[float]]: + def add_to_seq_without_pre_post_gradient( + self, seq: pp.Sequence, n_echoes: int, polarity: Literal['positive', 'negative'] = 'positive' + ): """Add readout gradients without pre- and re-winder gradients. Often the pre- and re-winder gradients are played out at the same time as phase encoding gradients or spoiler @@ -254,6 +257,8 @@ def add_to_seq_without_pre_post_gradient(self, seq: pp.Sequence, n_echoes: int) PyPulseq Sequence object. n_echoes Number of echoes + polarity + Polarity of first readout gradient Returns ------- @@ -266,9 +271,8 @@ def add_to_seq_without_pre_post_gradient(self, seq: pp.Sequence, n_echoes: int) time_to_echoes = [] for echo_ in range(n_echoes): start_of_current_gx = sum(seq.block_durations.values()) - gx_sign = (-1) ** echo_ + gx_sign = (-1) ** echo_ if polarity == 'positive' else (-1) ** (echo_ + 1) labels = [] - labels.append(pp.make_label(type='SET', label='REV', value=gx_sign == -1)) labels.append(pp.make_label(label='REV', type='SET', value=gx_sign == -1)) labels.append(pp.make_label(label='ECO', type='SET', value=echo_)) seq.add_block(pp.scale_grad(self._gx, gx_sign), self._adc, *labels) @@ -277,8 +281,8 @@ def add_to_seq_without_pre_post_gradient(self, seq: pp.Sequence, n_echoes: int) ) start_of_current_gx = sum(seq.block_durations.values()) if echo_ < n_echoes - 1: - if self._te_delay > 0: - seq.add_block(pp.make_delay(self._te_delay)) + if self._delta_te_delay > 0: + seq.add_block(pp.make_delay(self._delta_te_delay)) seq.add_block(pp.scale_grad(self._gx_between, -gx_sign)) return seq, time_to_echoes diff --git a/tests/scripts/test_grpe_flash_dixon.py b/tests/scripts/test_grpe_flash_dixon.py new file mode 100644 index 0000000..34759eb --- /dev/null +++ b/tests/scripts/test_grpe_flash_dixon.py @@ -0,0 +1,48 @@ +"""Tests for 3D FLASH sequence with golden radial phase encoding.""" + +import pytest +from mrseq.scripts.grpe_flash_dixon import main as create_seq + +EXPECTED_DUR = 17.84831 # defined 2025-11-09 + + +def test_default_seq_duration(system_defaults): + """Test if default values result in expected sequence duration.""" + seq, _ = create_seq(system=system_defaults, show_plots=False) + duration = seq.duration()[0] + assert duration == pytest.approx(EXPECTED_DUR) + + +def test_seq_creation_error_on_short_tr(system_defaults): + """Test if error is raised on too short repetition time.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, tr=5e-3, show_plots=False) + + +def test_seq_duration_vary_params_without_effect(system_defaults): + """Test if sequence duration is as expected.""" + seq, _ = create_seq( + system=system_defaults, + n_rpe_points_per_shot=4, # default 8 + show_plots=False, + test_report=False, + timing_check=False, + ) + duration = seq.duration()[0] + assert duration == pytest.approx(EXPECTED_DUR) + + +def test_seq_creation_error_on_wrong_partial_echo_factor(system_defaults): + """Test if error is raised on wrong partial echo factor.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, partial_echo_factor=1.1, show_plots=False) + with pytest.raises(ValueError): + create_seq(system=system_defaults, partial_echo_factor=0.4, show_plots=False) + + +def test_seq_creation_error_on_wrong_partial_fourier_factor(system_defaults): + """Test if error is raised on wrong partial Fourier factor.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, partial_fourier_factor=1.1, show_plots=False) + with pytest.raises(ValueError): + create_seq(system=system_defaults, partial_fourier_factor=0.4, show_plots=False) diff --git a/tests/utils/test_trajectory.py b/tests/utils/test_trajectory.py index 6a3cfdf..352c852 100644 --- a/tests/utils/test_trajectory.py +++ b/tests/utils/test_trajectory.py @@ -235,7 +235,10 @@ def test_multi_gradient_echo_error_on_short_delta_te(system_defaults): @pytest.mark.parametrize('readout_oversampling', [1, 1.5, 2]) @pytest.mark.parametrize('n_readout', [64, 128, 200]) @pytest.mark.parametrize('partial_echo_factor', [1.0, 0.8, 0.7]) -def test_multi_gradient_echo_timing(n_echoes, readout_oversampling, n_readout, partial_echo_factor, system_defaults): +@pytest.mark.parametrize('polarity', ['positive', 'negative']) +def test_multi_gradient_echo_timing( + n_echoes, readout_oversampling, n_readout, partial_echo_factor, polarity, system_defaults +): """Test that zero crossing of gradient moment coincides with echo time and correct adc sample.""" seq = pp.Sequence(system=system_defaults) mecho = MultiEchoAcquisition( @@ -244,10 +247,7 @@ def test_multi_gradient_echo_timing(n_echoes, readout_oversampling, n_readout, p readout_oversampling=readout_oversampling, partial_echo_factor=partial_echo_factor, ) - seq, time_to_echoes = mecho.add_to_seq(seq, n_echoes) - if n_echoes > 1: - # Ensure that the delta TE is the same between all echoes - assert len(np.unique(np.round(np.diff(time_to_echoes), decimals=6))) == 1 + seq, time_to_echoes = mecho.add_to_seq(seq, n_echoes, polarity) # Get full waveform for readout gradient w = seq.waveforms_and_times() @@ -262,7 +262,7 @@ def test_multi_gradient_echo_timing(n_echoes, readout_oversampling, n_readout, p k0_idx = argrelextrema(np.abs(m0_intp), np.less, order=100)[0] # Remove k0-crossings at the beginning and end of the block - k0_idx = np.array([ki for ki in k0_idx if (ki > 100 and ki < len(dt) - 100)]) + k0_idx = np.asarray([ki for ki in k0_idx if (ki > 100 and ki < len(dt) - 100)]) # Zero-crossing of the readout gradient should be within +/- .5 adc dwell time of the k-space center sample which is # the (_n_readout_pre_echo + 1)th sample. This should also be the same as the time_to_echo - time From 63ccd3e5de2d86cfab46ad194b402daeb59f4051 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 5 Mar 2026 14:31:57 +0000 Subject: [PATCH 05/13] [pre-commit] pre-commit autoupdate (#73) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Christoph Kolbitsch --- .pre-commit-config.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 5b67024..fec4759 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -15,7 +15,7 @@ repos: - id: mixed-line-ending - repo: https://github.com/kynan/nbstripout - rev: 0.9.0 + rev: 0.9.1 hooks: # cleans the .ipynbs (removes outputs, resets all cell-ids to 0..N, cleans steps) # also clean any kernel information left after execution @@ -25,14 +25,14 @@ repos: files: examples/.*\.ipynb$ - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.14.14 + rev: v0.15.4 hooks: - id: ruff # linter args: [--fix] - id: ruff-format # formatter - repo: https://github.com/crate-ci/typos - rev: v1.43.0 + rev: v1.44.0 hooks: - id: typos From 5472e25a6e70b05f639b35415d1f6e379e7764f1 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 5 Mar 2026 14:56:20 +0000 Subject: [PATCH 06/13] Bump MishaKav/pytest-coverage-comment from 1.4.0 to 1.5.0 in the all-actions group (#72) Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/report_pytest.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/report_pytest.yml b/.github/workflows/report_pytest.yml index 29c8a7e..d4c6424 100644 --- a/.github/workflows/report_pytest.yml +++ b/.github/workflows/report_pytest.yml @@ -49,7 +49,7 @@ jobs: - name: Post PyTest Coverage Comment id: coverage_comment - uses: MishaKav/pytest-coverage-comment@v1.4.0 + uses: MishaKav/pytest-coverage-comment@v1.5.0 with: issue-number: ${{ steps.pr-context.outputs.number }} pytest-coverage-path: pytest-coverage.txt @@ -90,7 +90,7 @@ jobs: - name: Post PyTest Coverage Comment on push id: coverage_comment - uses: MishaKav/pytest-coverage-comment@v1.4.0 + uses: MishaKav/pytest-coverage-comment@v1.5.0 with: pytest-coverage-path: pytest-coverage.txt junitxml-path: pytest.xml From 1d3c4d0b3bb951a773c47a37f456ea8b60334aa7 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 9 Mar 2026 06:31:41 +0100 Subject: [PATCH 07/13] Bump MishaKav/pytest-coverage-comment from 1.5.0 to 1.6.0 in the all-actions group (#74) Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/report_pytest.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/report_pytest.yml b/.github/workflows/report_pytest.yml index d4c6424..5342767 100644 --- a/.github/workflows/report_pytest.yml +++ b/.github/workflows/report_pytest.yml @@ -49,7 +49,7 @@ jobs: - name: Post PyTest Coverage Comment id: coverage_comment - uses: MishaKav/pytest-coverage-comment@v1.5.0 + uses: MishaKav/pytest-coverage-comment@v1.6.0 with: issue-number: ${{ steps.pr-context.outputs.number }} pytest-coverage-path: pytest-coverage.txt @@ -90,7 +90,7 @@ jobs: - name: Post PyTest Coverage Comment on push id: coverage_comment - uses: MishaKav/pytest-coverage-comment@v1.5.0 + uses: MishaKav/pytest-coverage-comment@v1.6.0 with: pytest-coverage-path: pytest-coverage.txt junitxml-path: pytest.xml From 86861efb0213e706306a7eb990c8ff13ef7f7c72 Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Wed, 11 Mar 2026 16:13:41 +0000 Subject: [PATCH 08/13] Fix dictionary matching for MOLLI (#66) Co-authored-by: Patrick Schuenke <37338697+schuenke@users.noreply.github.com> --- examples/t1_molli_bssfp.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/t1_molli_bssfp.ipynb b/examples/t1_molli_bssfp.ipynb index b4631fd..199d4c3 100644 --- a/examples/t1_molli_bssfp.ipynb +++ b/examples/t1_molli_bssfp.ipynb @@ -245,7 +245,7 @@ "source": [ "dictionary = DictionaryMatchOp(MOLLI(ti=torch.as_tensor(ti, dtype=torch.float32)), index_of_scaling_parameter=0)\n", "dictionary.append(\n", - " torch.tensor(1.0), torch.linspace(0.1, 5.0, 1000)[None, :], torch.linspace(0.1, 5.0, 1000)[None, None, :]\n", + " torch.tensor(1.0), torch.linspace(0.1, 5.0, 300)[None, :, None], torch.linspace(0.1, 5.0, 300)[None, None, :]\n", ")\n", "m0_match, c_match, t1_match = dictionary(idata.data)\n", "\n", @@ -270,7 +270,7 @@ "\n", "relative_error = np.sum(np.abs(t1_input - t1_measured)) / np.sum(np.abs(t1_input))\n", "print(f'Relative error {relative_error}')\n", - "assert relative_error < 0.08" + "assert relative_error < 0.12" ] } ], From 74f3e64c72358459d42203e583d32b058c23a67c Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Fri, 13 Mar 2026 19:21:50 +0000 Subject: [PATCH 09/13] 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 10/13] 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 11/13] 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 12/13] 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 13/13] 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,