Fix trajectory scaling in KTrajectoryPulseq

[schuenke]

I came across a bug in the reshape_pulseq_traj function in our KTrajectoryPulseq class.

In KTrajectoryPulseq, we use seq.calculate_kspace() to get our k-space trajectory from the seq-file. Due to some rounding errors etc it happens that we get different kz values even for single 2D acquisitions. These values are usually about 10 orders of magnitude smaller compared to typical kx and ky values.

However, in the reshape_pulseq_traj function we rescale the values with encoding_size / (2 * torch.max(torch.abs(k_traj))), which can be in the order of 1e8 for very small kz values (caused by mentioned rounding errors). For encoding_size = 1, the resulting max value of kz after rescaling was always 0.5 and thus 2 orders of magnitude larger than our repeat_detection_tolerance, which is meant to compensate for the small rounding errors.

I actually don't understand why this hasn't been a problem in the past, but @fzimmermann89 and I are sure that the scaling for encoding_size = 1 is a bug, which should be fixed by this PR.

[github-actions]

Coverage Report

File	Stmts	Miss	Cover	Missing
src/mrpro/algorithms/csm
inati.py	24	1	96%	44
walsh.py	16	1	94%	34
src/mrpro/algorithms/dcf
dcf_voronoi.py	53	4	92%	15, 48–49, 76
src/mrpro/algorithms/optimizers
adam.py	20	1	95%	69
src/mrpro/algorithms/reconstruction
DirectReconstruction.py	28	16	43%	51–71, 85
IterativeSENSEReconstruction.py	13	1	92%	76
Reconstruction.py	50	22	56%	42, 54–56, 80–87, 104–113
RegularizedIterativeSENSEReconstruction.py	41	17	59%	96–100, 114–139
src/mrpro/data
AcqInfo.py	128	3	98%	26, 169, 207
CsmData.py	29	3	90%	15, 82–84
DcfData.py	45	8	82%	18, 66, 78–83
IData.py	67	9	87%	119, 125, 129, 159–167
IHeader.py	75	7	91%	75, 109, 127–131
KHeader.py	153	17	89%	25, 119–123, 150, 199, 210, 217–218, 221, 228, 260–271
KNoise.py	31	15	52%	39–52, 56–61
KTrajectory.py	69	5	93%	178–182
MoveDataMixin.py	140	18	87%	15, 113, 129, 143–145, 207, 323–325, 338, 417, 437–438, 440, 455–456, 458
QData.py	39	7	82%	42, 65–73
Rotation.py	674	35	95%	100, 198, 335, 433, 477, 495, 581, 583, 592, 626, 628, 691, 768, 773, 776, 791, 808, 813, 889, 1077, 1082, 1085, 1109, 1113, 1240, 1242, 1250–1251, 1315, 1397, 1690, 1846, 1881, 1885, 1996
SpatialDimension.py	232	21	91%	34, 104, 141, 148, 154, 274–276, 289–291, 325, 343, 356, 369, 382, 395, 404–405, 420, 429
acq_filters.py	12	1	92%	47
src/mrpro/data/_kdata
KData.py	134	18	87%	109–110, 125, 132, 142, 150, 204–205, 243, 248–249, 268–279
KDataRemoveOsMixin.py	29	2	93%	44, 46
KDataSelectMixin.py	19	2	89%	48, 63
KDataSplitMixin.py	48	3	94%	53, 84, 93
src/mrpro/data/traj_calculators
KTrajectoryCalculator.py	25	2	92%	23, 45
KTrajectoryIsmrmrd.py	13	2	85%	41, 50
KTrajectoryPulseq.py	31	1	97%	54
src/mrpro/operators
CartesianSamplingOp.py	89	3	97%	118, 157, 280
ConstraintsOp.py	60	2	97%	46, 48
EndomorphOperator.py	65	2	97%	228, 234
FiniteDifferenceOp.py	27	2	93%	40, 105
FourierOp.py	158	3	98%	263, 381, 386
Functional.py	71	5	93%	20–22, 117, 119
GridSamplingOp.py	136	9	93%	72–73, 82–83, 90–91, 94, 96, 98
LinearOperator.py	168	10	94%	55, 91, 190, 220, 261, 270, 278, 287, 295, 320
LinearOperatorMatrix.py	158	16	90%	82, 119, 152, 161, 166, 175–178, 191–194, 203, 215, 304, 331, 359
MultiIdentityOp.py	13	2	85%	43, 48
Operator.py	78	2	97%	25, 74
ProximableFunctionalSeparableSum.py	39	3	92%	50, 103, 110
SliceProjectionOp.py	173	8	95%	44, 61, 63, 69, 206, 227, 260, 300
WaveletOp.py	120	5	96%	152, 170, 205, 210, 233
ZeroPadOp.py	16	1	94%	30
src/mrpro/utils
filters.py	62	2	97%	44, 49
slice_profiles.py	46	6	87%	20, 36, 113–116, 149
sliding_window.py	34	1	97%	34
split_idx.py	10	2	80%	43, 47
summarize_tensorvalues.py	11	9	18%	20–29
typing.py	18	11	39%	8–23
zero_pad_or_crop.py	31	6	81%	26, 30, 54, 57, 60, 63
TOTAL	4871	352	93%

Tests	Skipped	Failures	Errors	Time
2257	0 💤	0 ❌	0 🔥	2m 7s ⏱️

[github-actions]

📚 Documentation

📁 Download as zip
🔍 View online

[ckolbPTB]

I actually don't understand why this hasn't been a problem in the past

I guess for 2D trajectories it does not matter because the k2 dimension would then be singleton and we simply ignore it in the FourierOp. The only one using 3D pulseq-trajectories is probably @JoHa0811 - not sure if he encountered any problems.

[fzimmermann89]

Patrick noticed the bug because the reconstruction did not work.
The kz points were not singleton, but some noise, scaled to -.5 ... .5.
So they are only considered singleton if you increase the detection tolerance to 0.5.
(which could also be a sensible default, if the biggest difference between k space points is less than 1, it surely is singleton)

The reconstruction should still have worked, in theory, but, we also have #553 , the fourier op does not work with recon matrix smaller 6.

[ckolbPTB]

Patrick noticed the bug because the reconstruction did not work.

For 2D acquisitions? Strange, because we also have pulseq reconstructions in our examples

[schuenke]

I actually don't understand why this hasn't been a problem in the past

I guess for 2D trajectories it does not matter because the k2 dimension would then be singleton and we simply ignore it in the FourierOp. The only one using 3D pulseq-trajectories is probably @JoHa0811 - not sure if he encountered any problems.

It was a 2D trajectory in my case, but the kz values after rescaling differed in the range [-0.3, 0.5] and thus much more than the default repeat_detection_tolerance resulting in final kdata.traj dimensions

kz=[  1,   1, 402,   1]
ky=[  1,   1, 402, 256]
kx=[  1,   1, 402, 256]

insted of kz = [1, 1, 1, 1].

This actually lead to errors in the reco.

[fzimmermann89]

Patrick noticed the bug because the reconstruction did not work.

For 2D acquisitions? Strange, because we also have pulseq reconstructions in our examples

Indeed. I think it depends on the sequence if pypulseq returns 0 or something close to zero.
We should ask a pypulseq maintainer for input ;)

[fzimmermann89]

Patrick noticed the bug because the reconstruction did not work.

For 2D acquisitions? Strange, because we also have pulseq reconstructions in our examples

Indeed. I think it depends on the sequence if pypulseq returns 0 or something close to zero. We should ask a pypulseq maintainer for input ;)

Ok, in the pulseq example, kz is constant (but not 0). This avoids the possible division by zero bug and results in a singleton dimension as all kz after scaling are -0.5 ...

[schuenke]

We still want to merge this? Or is it superseeded by #560?
IMHO it's not a good idea to combine a major refactoring like #560 with a bugfix like this...