-
Notifications
You must be signed in to change notification settings - Fork 73
Description
Describe the bug
Changing the sampling duration in the example sequence write_epi.py
affects the SAR calculation. Since TE and TR are fixed for the sequence, changing the sampling duration does not affect the sequence duration and hence should not affect SAR calculation.
To Reproduce
Run the write_epi.py
script and append pp.SAR.SAR_calc.calc_SAR(seq)
at the end. Now change the sampling duration and re-run the script.
sampling_time = 6.4e-3/2
https://github.com/imr-framework/pypulseq/blob/dev/examples/scripts/write_epi.py
Expected behavior
SAR calculation should be independent on other sequence events that are not RF related and that do not change the timing of the sequence.
Additional context
What happens in the write_epi.py
is that changing the sampling duration affects the spoiler duration t_sp = 0.5 * (TE - readout_time - t_refwd)
. The SAR calculation is affected by the duration of the event following the RF event. In this case that is the duration of the spoiler. I think there is a bug in the SAR calculation. In Line 101 and further in SAR_calc.py
it reads:
for block_counter in block_events:
block = seq.get_block(block_counter)
block_dur = calc_duration(block)
t[block_counter - 1] = t_prev + block_dur
t_prev = t[block_counter - 1]
if hasattr(block, "rf"): # has rf
rf = block.rf
signal = rf.signal
# This rf could be parallel transmit as well
SAR_wbg[block_counter] = _calc_SAR(Qtmf, signal)
SAR_hg[block_counter] = _calc_SAR(Qhmf, signal)
The index of the time vector t
is decremented by 1, however the decrement is missing for the SAR_wbg
vector, which leads to a misalignment of time and SAR.
Here you can see the sequence diagram. Please be aware that I deactivated the dummy scans in write_epi.py
for sake of clarity.
And here is the corresponding time and SAR vector. The SAR for the 180 refocussing pulse is associated with a time stamp of 11.99 ms which is the end of the spoiler and not the end of the RF pulse.
IDX t[ms] RF
1: 0.25 0.0
2: 2.95 0.0
3: 6.5 16.5908
4: 8.7 0.0
5: 11.99 103.6925
6: 15.21 0.0
7: 18.5 0.0
8: 20.7 0.0
In my opinion the SAR vector should be decremented according to the time vector.
for block_counter in block_events:
block = seq.get_block(block_counter)
block_dur = calc_duration(block)
t[block_counter - 1] = t_prev + block_dur
t_prev = t[block_counter - 1]
if hasattr(block, "rf"): # has rf
rf = block.rf
signal = rf.signal
# This rf could be parallel transmit as well
SAR_wbg[block_counter - 1] = _calc_SAR(Qtmf, signal)
SAR_hg[block_counter - 1] = _calc_SAR(Qhmf, signal)
Thanks a lot for your support! I hope I understood the issue right.
Kind regards
Helge