Trajectory discrepancies between the sequence when is created and the same sequence when is being readed

[DelfinCalles]

Hello,

I have started using PyPulseq for the design of sequences and BART for the reconstruction of the signal coming from a Pulseq sequence. I've noticed that when trying to reconstruct a 2D image from a 2D radial sequence generated by the test case write_radial_gre.py BART would detect two pixels on the z-axis coming from the fact that the k_z frequencies were not zero.

When looking at the trajectories directly generated by the sequence with seq.calculate_kspace()[0] I've observed that while the k_z trajectories are originally close to 0 (~10^-12), once I look at the same sequence after doing seq.write('seq_name') and seq.read('seq_name'). Now the k_z trajectories have become way bigger (~10^-3).

I know there are some rounding errors when writing the .seq file, as well as some parameters that could cause some additional errors trying to mimic a true MRI scanner. I'm not sure if all of these could be the cause of this difference.

[rmz-oz]

I could reproduce this with pypulseq 1.5.0.post1 (and current main, f2c582b) using write_radial_gre.py. Before writing, max |kz| at the ADC samples is 7.9e-10 1/m (2.4e-12 once scaled by the 3 mm FOV in z, which is probably the ~1e-12 you saw). After seq.write() + seq.read() every ADC sample has the same kz = -1.265e-3 1/m. kx/ky move by at most 6.4e-4 1/m.

Where it comes from

It is the precision of the .seq text format. Trapezoid amplitudes are written with {:12g}, which is 6 significant digits:

pypulseq/src/pypulseq/Sequence/write_seq.py

Line 151 in f2c582b

id_format_str = '{:2.0f} {:12g} {:3.0f} {:4.0f} {:3.0f} {:3.0f}\n'

and remove_duplicates(), which write() calls by default, rounds them to 6 digits too:

pypulseq/src/pypulseq/Sequence/sequence.py

Line 1211 in f2c582b

seq_copy.grad_library, mapping = seq_copy.grad_library.remove_duplicates((6, -6, -6, -6, -6, -6))

MATLAB Pulseq writes [TRAP] the same way (%12g), so this is how the file format behaves in general and not something specific to the pypulseq reader:
https://github.com/pulseq/pulseq/blob/c7469123c2f381f065986e6cc3a7d09730ed16ef/matlab/+mr/@Sequence/write.m#L118

In this example the numbers are:

event	amplitude in memory	in the file	error
slice select gz	333333.333 Hz/m	333333	-0.333 Hz/m
slice rephaser gz_reph	-385416.667 Hz/m	-385417	-0.333 Hz/m

The part of gz after the RF center lasts an effective 2035 us, and the rephaser 1760 us. That gives -0.333 * 2035e-6 + (-0.333 * 1760e-6) = -6.78e-4 - 5.87e-4 = -1.265e-3 1/m, which matches the measured offset. Two more checks:

• If I patch those two lines locally to keep 12 digits, the before/after difference drops to 1.5e-9 1/m.
• Changing the slice thickness changes the offset exactly as this rounding predicts: +6.4e-4 (2 mm), +8.0e-5 (4 mm), +8.0e-4 (5 mm).

Physically it doesn't matter. For a 3 mm slice, one full kz cycle across the slice is 333 1/m, so 1.3e-3 1/m is about 4e-6 of a cycle. Also, the read-back trajectory describes the amplitudes that are actually stored in the .seq file you run.

Why BART gives 2 pixels in z

BART estimates the image size from the trajectory. A dimension only gets size 1 if its max |k| is exactly 0; otherwise it gets 2 * ceil(max|k|):
https://github.com/mrirecon/bart/blob/731bfd3ec4f09ffa59b497ded1274a2c0357c8af/src/misc/mri2.c#L260

So any non-zero kz gives 2 pixels in z, including the tiny value from before writing. I ran BART v1.0.00 with trajectories scaled as k * FOV (FOV = [0.26, 0.26, 0.003] m, as in the example's FOV definition):

• bart nufft -a with the trajectory from before writing: 64x64x2
• bart nufft -a with the read-back trajectory: 66x66x2. The extra 2 in x/y is the same rounding: max |kx * FOV| is 32.000004 instead of 32, so the ceil goes to 33.
• bart nufft -a -x 64:64:1 with either of these trajectories aborts with Assertion failed: (0. == pos[j]) (grid.c line 437): https://github.com/mrirecon/bart/blob/731bfd3ec4f09ffa59b497ded1274a2c0357c8af/src/noncart/grid.c#L437

What to do

For a 2D acquisition, set kz to exactly zero when you build the BART trajectory, and pass the matrix size so BART doesn't have to guess it:

k_adc = seq.calculate_kspace()[0]                        # (3, n_x * n_spokes), 1/m
traj = k_adc * np.array([fov, fov, 0.0])[:, None]        # BART units (1/FOV), kz forced to 0 for 2D
traj = traj.reshape(3, n_spokes, n_x).transpose(0, 2, 1) # (3, samples, spokes)

bart nufft -a -x 64:64:1 traj ksp img

With this, BART runs on both the original and the read-back trajectory and gives a 64x64 image. On the same k-space data, the two images differ by 1.1e-5 in relative norm, so the leftover kx/ky rounding doesn't matter either.

Tested with pypulseq 1.5.0.post1 (numpy 2.5.3, Python 3.12) and pypulseq main at f2c582b. BART v1.0.00 was built from mrirecon/bart at 731bfd3.