pypulseq
pypulseq copied to clipboard
Sub-microsecond ADC delays are dropped without notification when exporting to .seq file
When adding ADC events that are starting on a valid time according to the adc-raster time but is smaller than a micro second the delay is simply dropped when exported to a .seq file.
This consistent with the file in section 2.8.3, as the table states delays are specified in integers of microseconds.
Below I have a minimally working example that creates a readout trapezoidal and a corresponding ADC from a number of samples and a flat-duration. To obtain a valid dwell time, the actual ADC duration is recomputed to be slighlty shorter, and subsequently the adc is shifted to be centered with respect to the flat-top duration.
Expected behaviour I suspect severity of this depends on the choses raster_times, and in most cases with the default siemens raster times the differences might be neglegible. However, a Userwarning that states that this conversion is happening might be a good idea, so the user can choose e.g. a different number of samples or delibarately choose to ignore this difference.
To Reproduce
import pypulseq as pp
import numpy as np
system = pp.Opts(max_grad=32, grad_unit="mT/m", max_slew=130, slew_unit="T/m/s",
rf_ringdown_time=30e-6, rf_dead_time=100e-6)
# Showcase invalid dwell time
Nx = 121
flat_time = 1e-3
readout_time = flat_time
dwell_time = readout_time / Nx
print(f"Not a valid dwell time: {dwell_time}")
>>> Not a valid dwell time: 8.264462809917356e-06
# Define valid dwell time
Nx = 121
flat_time = 1e-3
dwell = np.floor(flat_time / Nx / system.adc_raster_time) * system.adc_raster_time
readout_time = dwell * Nx
difference = flat_time - readout_time
print(f"Difference ADC/Flat top: {difference}")
>>> Difference ADC/Flat top: 7.800000000000081e-06
# Center the adc on flat-top duration to ensure symmetric k-space sampling
adc_delay = difference / 2
print(f"Padding prior to ADC: {adc_delay}")
>>> Padding prior to ADC: 3.900000000000041e-06
fov = 220e-3
delta_k = 1 / fov
k_width = Nx * delta_k
seq = pp.Sequence()
gx = pp.make_trapezoid(channel="x", system=system,
amplitude=k_width / readout_time, flat_time=flat_time)
adc = pp.make_adc(num_samples=Nx, duration=readout_time, delay=adc_delay)
seq.add_block(gx, adc)
seq.write("dummy_adc.seq")
Which results in the following file, stating a 4us delay instead of a 3.9us delay.
Hi @jweine. As Siemens interpreter also uses microsecond for that delay, as you said this will be consistent with Siemens convention. In Siemens, as the time is an integer, that (flat_time - readout_time)/2 division will silently do the same thing.
We may be able to add a check and print a warning if the difference after integer cast is larger than eps.
The functionality to check the timing of your sequence already exists:
ok, error_report = seq.check_timing()
if ok:
print("Timing check passed successfully")
else:
print("Timing check failed. Error listing follows:")
[print(e) for e in error_report]
This works for any timing related value, e.g. gradient timing as well. You should really always use this. It might make sense for seq.write to include a timing check.
I did notice there's a little bug in check_timing in that it does not print multiple timing errors within the same block (e.g. in the example if you provide a bad flat_time, it will still only report the ADC delay time error).
I've not been happy about the way to print the report, it would make sense to either just print a report with nicer formatting indicating where the errors are, or to return a data structure that lists the specific errors (not as a string, e.g. a list of dict(block=1, event='adc', field='delay', value=..., rounded_value=...)), and provide a function to nicely print that list.
I did not realize seq.check_timing() checks for ADC timing, good to know!
I realized that it also assumes rf_raster_time and adc_raster_time are the same?
https://github.com/imr-framework/pypulseq/blob/25938e8c1c9d23dbdb0420dace48174ea103ff3b/pypulseq/check_timing.py#L58-L59
A dict is also awesome if one wants to easily search for certain events in a spam of errors.
I believe the answers and #200 addresses the issue. Closing.