Gridding a sequence containing trapezoidal gradients results in 'cut off' corners and a time-grid shift of grad_raster_time/2

[jweine]

Describe the bug When creating a sequence that contains trapezoidal gradient lobes using make_trapezoid , the result of Sequence.gradient_waveforms() shows an incorrectly gridded waveform. When gridded (by calling the function Sequence.gridded_wave_forms() which in subsequently calls the function `pypulseq.points_to_waveform'), the waveform timing starts at grad_raster_time/2. When defining trapezoids such that the interval border of the piecewise linear function lie exactly on a regular grid with spacing grad_raster_time, this results in 'capped off' corners, a non-zero start value as well as a bend waveform in the last interval.

Not sure why the function pypulseq.points_to_waveform performs the time-grid shift of grad_raster_time/2 but for this case it should not be applied.

So, if this is intentional in some cases this conditional needs to be revised to catch the correct case for the trapezoidal gradients.

To Reproduce

import pypulseq
import matplotlib.pyplot as plt
import numpy as np

system = pypulseq.Opts(max_grad=80, grad_unit='mT/m', max_slew=200,
                       slew_unit='mT/m/ms', grad_raster_time=0.01)
seq = pypulseq.Sequence(system=system)
l1 = pypulseq.make_trapezoid(channel="x", flat_time=5.6, amplitude=1., rise_time=0.9)

seq.add_block(l1)
wf = seq.gradient_waveforms()
t = np.arange(0, seq.duration()[0]+seq.grad_raster_time, seq.grad_raster_time)

f, a = plt.subplots(1, 1)
a.plot(t, wf.T)
a.grid(True)

Expected behavior
My expectation is, that the reference-points of normalized amplitude [0, 1, 1, 0] at times [0, rise_time, rise_time+flat_time, rise_time+flat_time+fall_time] should be met exactly if all time-points match the grad_raster_time as shown in the figure below.

Desktop (please complete the following information):

• OS: [Windows 11 & Ubuntu 20.04]
  • pypulseq version: 1.3.1.post1

Additional context To catch errors like this i propose to add a test similar to the following:

import unittest
import pypulseq
import numpy as np

system = pypulseq.Opts(max_grad=80, grad_unit='mT/m', max_slew=200,
                       slew_unit='mT/m/ms', grad_raster_time=0.01)

class TestGradientWaveform(unittest.TestCase):
    def test_correct_start_trap(self):
        seq = pypulseq.Sequence(system=system)
        l1 = pypulseq.make_trapezoid(channel="x", flat_time=5.6, amplitude=1., rise_time=0.9,
                                     delay=0)
        seq.add_block(l1)
        wf = seq.gradient_waveforms()
        self.assertTrue(np.allclose(wf[:, 0], 0, rtol=1e-10),
                        msg="First gridded sample is not zero")
        self.assertTrue(np.allclose(wf[:, -1], 0, rtol=1e-10),
                        msg="Last gridded sample is not zero")

[sravan953]

@jweine Issue #72 was related to this, fixed in PR #84. Please try it out and let me know!

[btasdelen]

@sravan953 Negative, unit test still fails:

======================================================================
FAIL: test_correct_start_trap (GradWaveformUT.TestGradientWaveform)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/bilal/Insync/tasdelen@usc.edu/GoogleDrive/codebase/pulseq_sequences/bilal_seqs/GradWaveformUT.py", line 16, in test_correct_start_trap
    self.assertTrue(np.allclose(wf[:, 0], 0, rtol=1e-10),
AssertionError: False is not true : First gridded sample is not zero

----------------------------------------------------------------------

I feel like the issue arises from the use of interpolation in the points_to_waveform(), as @jweine pointed out. I will take a look at the issue and see if I can come up with a better solution.

[btasdelen]

Update: Behavior is the same between 'Pulseq' and 'PyPulseq'. 'The issue' comes from the convention Siemens uses, shifting gradient sampling by half raster time, or equivalently, sampling the gradients at the center, not the edges. I believe confusion comes from the fact that, we use shifted gradients for visualization, when in fact, it makes more sense to display non shifted ones.

If we also plot the time-axis, as the first sample will not be at time=0, but time=GRT/2, it all makes sense. So I believe there is no bug here, just a confusing feature.

[jweine]

@bilal-tasdelen Thanks for following up on this. The existence of this shift is mentioned in the pulseq file format definition when using arbitrary gradients and rf-waveforms but not Trapezoidals. But I couldn't find the actual reason for using the shift / center-sampling. I feel it should be documented that this is a Siemens specific feature, as i don't think it is intuitive without more context (especially for non-siemens sites :)). Personally I think, this behaviour should not be presented to user but handled in the Siemens sequence export. But this is more of a design decision to be discussed also with the Pulseq project. Anyway, thanks for the clarification!