imr-framework / pypulseq

Pulseq in Python
https://pypulseq.readthedocs.io
GNU Affero General Public License v3.0
111 stars 58 forks source link

add_blocks interpolation issue #168

Closed fzimmermann89 closed 5 months ago

fzimmermann89 commented 5 months ago

The documentation of pypulseq.split_gradient.split_gradient reads

Splits a trapezoidal gradient into slew up, flat top and slew down. Returns the individual gradient parts (slew up, flat top and slew down) as extended trapezoid gradient objects. The delays in the individual gradient events are adapted such that addGradients(…) produces an gradient equivalent to ‘grad’.

but

g1=pypulseq.make_trapezoid('z',flat_time=1e-4,rise_time=1e-4,amplitude=1)
g2=pypulseq.add_gradients(pypulseq.split_gradient(g1))
print(g1)
print(g2)

prints

namespace(type='trap',
          channel='z',
          amplitude=1,
          rise_time=0.0001,
          flat_time=0.0001,
          fall_time=0.0001,
          area=0.0002,
          flat_area=0.0001,
          delay=0.0,
          first=0,
          last=0)

and

namespace(type='grad',
          channel='z',
          waveform=array([2., 3., 3., 2.]),
          delay=0.0,
          tt=array([0.    , 0.0001, 0.0002, 0.0003]),
          shape_dur=0.00030000000000000003,
          first=2.0,
          last=2.0)

If I am not mistaken, the "equivalent gradient" should not have a max. waveform of 3, but 1. I believe there might be a bug inadd_gradients. --> Somewhere in the interpolation logic, both rising block gets and falling block get a flat part added where they should be 0 (the flat part is a separate block..) This also happens in more realistic settings, but this is the smallest repoducing code.

For reference, the split (pypulseq.split_gradient(g1)) is:


namespace(type='grad',
           channel='z',
           waveform=array([0, 1]),
           delay=0.0,
           tt=array([0.    , 0.0001]),
           shape_dur=0.0001,
           first=0,
           last=1)

 namespace(type='grad',
           channel='z',
           waveform=array([1, 1]),
           delay=0.0001,
           tt=array([0.    , 0.0001]),
           shape_dur=0.0001,
           first=1,
           last=1)

 namespace(type='grad',
           channel='z',
           waveform=array([1, 0]),
           delay=0.00020000000000000004,
           tt=array([0.    , 0.0001]),
           shape_dur=0.0001,
           first=1,
           last=0)
fzimmermann89 commented 5 months ago

Ok, this might be solved by #167, sorry....

fzimmermann89 commented 5 months ago

Yes, it is fixed with #167,

namespace(type='grad',
          channel='z',
          waveform=array([0., 1., 1., 0.]),
          delay=0.0,
          tt=array([0.    , 0.0001, 0.0002, 0.0003]),
          shape_dur=0.00030000000000000003,
          area=0.00020000000000000004,
          first=0.0,
          last=0.0)

So thank you for the superb response time of negative 6 hours ! 🎉🎉🎉 ;-)