Shift in trajectory between TR with calculate_kspace()

[felixlandmeyer]

Hello,

this is my first post here. I hope I did not overlook something in the issues section regarding this matter...

In version 1.4.0 there is a k-space shift between TR when calling the function calculate_kspace if you use arbitrary gradients.

This does not happen in version 1.3.1

I'm just wondering if this is intentional, or if this could be an issue?

I've added two examples with interleaved spirals for version 1.3.1 and 1.4.0 (with https://github.com/FrankZijlstra/pypulseq/tree/bugfixes installation).

Spirals were created with the https://github.com/mrphysics-bonn/spiraltraj package.

Code for reproduction:

V1.4.0:

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


def rotation(phi, grad):
    grad_rot = np.zeros_like(grad)
    grad_rot[:, 0] = np.cos(phi)*grad[:, 0]-np.sin(phi)*grad[:, 1]
    grad_rot[:, 1] = np.sin(phi)*grad[:, 0]+np.cos(phi)*grad[:, 1]
    
    return grad_rot

max_slew = 200
max_grad = 45
system = pp.Opts(
    max_grad=max_grad, grad_unit='mT/m', max_slew=max_slew, slew_unit='T/m/s',
    rf_ringdown_time=30e-6, rf_dead_time=100e-6)

Rf, Gz, Gzr = pp.make_sinc_pulse(
    flip_angle=20*np.pi/180, system=system,
    duration=3e-3, slice_thickness=3e-3,
    apodization=0.5, time_bw_product=4,
    use='excitation', return_gz=True)

Nlve = 5
fov = 200
res = 5

spiral = spiraltraj.calc_traj(
    nitlv=Nlve, fov=fov, res=res, max_amp=29, min_rise=6, spiraltype=1)

spiral = np.array(spiral) * 1e-3 * system.gamma

#spiral = np.concatenate((spiral, np.zeros((5,2))), axis=0)

spiral_rot = 2*np.pi/Nlve

Gx = pp.make_arbitrary_grad(channel='x', waveform=spiral[:, 0], system=system)
Gy = pp.make_arbitrary_grad(channel='y', waveform=spiral[:, 1], system=system)

dwell = 5e-6
adc_samples = int(Gx.shape_dur/dwell)
Adc = pp.make_adc(dwell=dwell, num_samples=adc_samples, system=system)

Seq = pp.Sequence()

for n_lve in range(Nlve):
    Seq.add_block(Rf, Gz)
    Seq.add_block(Gzr)
    Seq.add_block(pp.make_delay(2e-3))
    Seq.add_block(*pp.rotate(Gx, Gy, Adc, angle=spiral_rot*n_lve, axis='z'))

    
ktraj_adc, ktraj, t_ext, t_ref, t_adc = Seq.calculate_kspace()

ktraj_adc[-1] = 0
Seq.plot()

plt.figure()
plt.plot(ktraj_adc[0, :], ktraj_adc[1, :], marker='.', ls='')

V1.3.1:

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


def rotation(phi, grad):
    grad_rot = np.zeros_like(grad)
    grad_rot[:, 0] = np.cos(phi)*grad[:, 0]-np.sin(phi)*grad[:, 1]
    grad_rot[:, 1] = np.sin(phi)*grad[:, 0]+np.cos(phi)*grad[:, 1]
    
    return grad_rot

from pypulseq.make_arbitrary_grad import make_arbitrary_grad

max_slew = 200
max_grad = 45
system = pp.Opts(
    max_grad=max_grad, grad_unit='mT/m', max_slew=max_slew, slew_unit='T/m/s',
    rf_ringdown_time=30e-6, rf_dead_time=100e-6)

Rf, Gz, Gzr = pp.make_sinc_pulse(
    flip_angle=20*np.pi/180, system=system,
    duration=3e-3, slice_thickness=3e-3,
    apodization=0.5, time_bw_product=4,
    use='excitation', return_gz=True)

Nlve = 5
fov = 200
res = 5

spiral = spiraltraj.calc_traj(
    nitlv=Nlve, fov=fov, res=res, max_amp=29, min_rise=6, spiraltype=1)

spiral = np.array(spiral) * 1e-3 * system.gamma

rot = 2*np.pi/Nlve

dwell = 5e-6
adc_samples = int(spiral.shape[0]*system.grad_raster_time/dwell)
Adc = pp.make_adc(dwell=dwell, num_samples=adc_samples, system=system)

Seq = pp.Sequence()

for n_lve in range(Nlve):
    Seq.add_block(Rf, Gz)
    Seq.add_block(Gzr)
    Seq.add_block(pp.make_delay(2e-3))
    
    spiral_rot = rotation(rot*n_lve, spiral)
    Gx = make_arbitrary_grad(channel='x', waveform=spiral_rot[:, 0], system=system)
    Gy = make_arbitrary_grad(channel='y', waveform=spiral_rot[:, 1], system=system)
    
    Seq.add_block(Gx, Gy, Adc)
    
ktraj_adc, ktraj, t_ext, t_ref, t_adc = Seq.calculate_kspace()

Seq.plot()
#
plt.figure()
plt.plot(ktraj_adc[0, :], ktraj_adc[1, :], marker='.', ls='')

[btasdelen]

Hi @felixlandmeyer. I ran your test case. The result is indeed as you shared. Here is the reason:

Spirals are not starting from zero. I guess the PyPulseq 1.3.1 was adding the zeros, but v1.4.0 is not. I feel like there should be a warning for such cases, I will take a look.

If you go ahead and uncomment line 36 in your Pulseq v1.4.0 and change it to spiral = np.concatenate((np.zeros((1,2)), spiral, np.zeros((1,2))), axis=0) to add zeros to the beginning and end (you don't need 5 zeros, 1 should be enough), here is the plot I get:

[felixlandmeyer]

Thank you @btasdelen, that solved the problem!

Yes, I did not consider the beginning of the trajectory.