ValueError when using complex64 with cgrappa

[zaccharieramzi]

Hi, thanks very much for this implementation.

I tried to use cgrappa with np.complex64 data and got the following error:

ValueError                                Traceback (most recent call last)
<ipython-input-4-2802ffa5c73b> in <module>
     12 # coil_axis=-1 is default, so if coil dimension is last we don't
     13 # need to explicity provide it
---> 14 res = cgrappa(kspace, calib, kernel_size=(5, 5))
     15 sx, sy, ncoils = res.shape[:]

src/cgrappa.pyx in cgrappa.cgrappa()

ValueError: Buffer dtype mismatch, expected 'double complex' but got 'complex float'

The code to reproduce the error is the following:

from pygrappa import cgrappa 
import numpy as np 

kspace = np.random.normal(size=[640, 372, 15]) 
kspace = kspace.astype(np.complex64) 


sx, sy, ncoils = kspace.shape[:] # center 20 lines are ACS 
ctr, pd = int(sy/2), 10 
calib = kspace[:, ctr-pd:ctr+pd, :].copy() # call copy()! 

# coil_axis=-1 is default, so if coil dimension is last we don't 
# need to explicity provide it 
res = cgrappa(kspace, calib, kernel_size=(5, 5)) 
sx, sy, ncoils = res.shape[:]             

[mckib2]

Hi @zaccharieramzi, thanks for reporting this. What operating system are you using and what version of python?

[zaccharieramzi]

I am on Linux Ubuntu 16.04 with Python 3.6.8.

[mckib2]

Alright, I knew this looked familiar. This is a known issue (not well documented, tests only run on complex128 because complex64 fails). It's a limitation of Cython's C++ template support. I have plans on doing a better, more flexible C++ implementation.

There's not an easy fix, but here are some options:

• Take the memory hit and use complex128
• Use mdgrappa which is a pure Python implementation with n-dimensional support until cgrappa gets an overhaul

[zaccharieramzi]

Not really in the scope of this issue but since you mention it, does mdgrappa use the 3rd dimension information or does loop over the 3rd dimension (what I am doing right now)?
In particular do you need isotropic resolution?

[mckib2]

mdgrappa assumes you want to use all spatial frequency dimensions in the reconstruction. If you have a 3D+coil dataset and you want to only reconstruct inplane, you would need to write a loop and pass mdgrappa a single slice each iteration. Would it be a useful feature for you to specify which axes to either use or exclude from the reconstruction?

Keep in mind many of the algorithms found in pygrappa currently only support 2D+coil datasets as I was trying to just get them working. They are slowly getting migrated to be able to handle 3D+coil datasets (e.g., most recently cgsense).

In particular do you need isotropic resolution?

I don't see any reason nonisotropic voxels shouldn't work, let me know if you run into issues.

[zaccharieramzi]

Yes that's what I did, looping over the 3rd dimension was my solution so far.

However, it is very slow (12 minutes for a fastMRI volume) and doesn't have acceptable results (13dB in PSNR). It probably has to do with the setting of the lamda regularisation parameter, does it sound plausible?
Have you played around with fastMRI yourself?

[mckib2]

Might be regularization term, could be something else. I've also had issues with some methods using 7T data (#74).

I haven't worked with fastMRI volumes. Do you mind sending a link to the one you are trying to use? Real datasets often behave very differently from synthetic ones and I'd like to start setting up some benchmarks with real data to make sure performance is acceptable.

[zaccharieramzi]

I am not sure I understood, do you want the link to the fastMRI dataset? In this case here you go: https://fastmri.med.nyu.edu/.

If you want like an extract from this dataset to run GRAPPA on I can zip one for you just tell me.

[mckib2]

I think the confusion is because I am not familiar with fastMRI.

If you want like an extract from this dataset to run GRAPPA on I can zip one for you just tell me.

This would probably be easiest for now. Do you mind sharing the dataset as a zip file?

[zaccharieramzi]

I can't really share the whole dataset as a zip file since it's 1TB but I ll try with just one volume.