Some Points discarded

[LiangWang-bme]

When I perform zero echo time(ZTE) sequence reconstruction and use recon-lsq for reconstruction, it indicates that some points have been discarded. What could be the reason for this? My sampling matrix is 256, and my trajectory size range is 0-128. What could cause the points to be discarded, or under what circumstances would points be discarded?

[spinicist]

Hello,

Short answer: Provided the number of points discarded is small (<1%?) I think your setup is correct and riesling is also correctly discarding points at the positive edge of Cartesian k-space, because they lie outside the definition of the FFT.

Longer answer: Sample points are discarded if they lie outside the Cartesian grid. The Cartesian gridpoints are defined by those required by the FFT, but shifted so that the zero-frequency bin is in "the center".

For an even transform size the FFT bins are by definition 0 to N/2 - 1, then -N/2 to -1. After shifting this becomes -N/2 to N/2 - 1. In other words, the FFT and hence Cartesian grid are asymmetric by one sample. As a concrete example, in 1D with N=8 the Cartesian grid points are -4 -3 -2 -1 0 1 2 3. Cartesian MRI sequences are setup such that they only sample these points (at least they should be).

However, ZTE and other radial MRI trajectories are usually implemented symmetrically. By this we mean the trajectory starts at zero frequency and proceeds a fixed distance in a particular direction. Again in 1D, this means we acquire one spoke that samples 0 1 2 3 4 and another that samples 0 -1 -2 -3 -4. This means the trajectory samples both the positive and negative edges of k-space. But the positive edge doesn't exist on the matching Cartesian grid. Hence there are is a choice as to what to do with these points:

1. Use an odd sized symmetric Cartesian grid. Slow and wasteful, the even sized grid should produce the same results.
2. Wrap the positive edge to the negative edge. This is actually what the FFT algorithm expects (because the definition of the FFT is that signals repeat/wrap) but not representative of the MRI physics (which does not wrap).
3. Discard the points, which is representative of the MRI physics.

I opted for option 3.

If you use riesling lsq --ndft then riesling will perform a Non-uniform Discrete Fourier Transform, i.e. it does not do gridding+FFT, and these spatial frequencies will get used during the reconstruction. Sadly, the NDFT is O(N^2) whereas the NUFFT is O(NlogN + N) (I think? It's a while since I checked), and is hence infeasibly slow for any realistic data set.

Let me know if the explanation isn't clear. I should probably add this to the FAQ.

[LiangWang-bme]

Yes, it indicates that no more than 0.1% of the points are discarded. However, the size of the reconstructed matrix will be missing, for example, the size of the reconstructed image is 128 * 127 * 128. Is this correct? I tried using recon-lsq --ndft instead of fft and found that the reconstruction time is too long, which is usually unacceptable. I think it would be good to add this to the FAQ.

[spinicist]

Have you added the matrix attribute to the trajectory dataset? This is required to explicitly set the matrix size, see here: https://riesling.readthedocs.io/en/latest/data.html#trajectory. If you do not provide this attribute riesling will try to guess the matrix size from the maximum extents of the trajectory in each direction. This can often be wrong by 1 voxel if your trajectory doesn't have spokes pointing perfectly in each direction. The way to create the attribute in Python is here: https://github.com/spinicist/riesling/blob/main/python/riesling/io.py#L32

[LiangWang-bme]

I used MATLAB code to add a matrix to trajectory with a size of 200, as shown in the figure below

Is the addition correct?
After the addition is completed, I encountered the following error when using recon-lsq:

[spinicist]

It needs to be an array of integers, the first screenshot looks like an array of floats. Can you post the Matlab code you used to create the attribute please?

[LiangWang-bme]

I tried converting the matrix values to integers and reconstructing with recon-lsq, but it still resulted in the error shown in the screenshot above. Below is the MATLAB code for adding the matrix under trajectory, for example, with the matrix size being [200,200,200].

`fileID = 'zte.h5';
datasetName = '/trajectory';

file = H5F.open(fileID, 'H5F_ACC_RDWR', 'H5P_DEFAULT');

dset = H5D.open(file, datasetName);

type = H5T.copy('H5T_NATIVE_INT');

dims = [1 3];
space = H5S.create_simple(2, fliplr(dims), []);

attr = H5A.create(dset, 'matrix', type, space, 'H5P_DEFAULT');

matrixData = int64([200, 200, 200]);
H5A.write(attr, 'H5T_NATIVE_INT64', matrixData);
H5A.close(attr);
H5D.close(dset);
H5F.close(file);

`

[spinicist]

Is it possible to share an example file created with this code? It looks correct, I'd like to check what type the C++ code thinks this attribute ends up being.