Aliasing artifacts

[estab1]

Hi, so I was trying to use grappa which worked as intented at first. But as soon as I started to create error images abs(reconstructed_image - reference_image) , I realized that grappa is producing aliasing artifacts. The error is also reproducible with the example provided in examples/basic_grappa.py. Therefore you only have to add the following code at the end of the example:

reference_img = np.zeros((2*N, 2*N))
kk = 0
for idx in np.ndindex((2, 2)):
    ii, jj = idx[:]
    reference_img[ii*N:(ii+1)*N, jj*N:(jj+1)*N] = abs(imspace[..., kk])
    kk += 1
plt.imshow(abs(reference_img-res0), cmap='gray')
plt.show()

Before using pygrappa I was using PULSAR which is a grappa implementation for matlab and aliasing artifacts weren't a problem. So my question is, am I missing something or is this a bug.

(I also tried using cgrappa and mdgrappa which led to the same results.)

[mckib2]

Hi @estab1 , thanks for reporting! I'll have more time to look at this a little later today, but I think I've run into relative scaling issues in the past, so normalizing both the reference and reconstruction before comparing may help. I've not used PULSAR before, but I can take a look and see if there are any obvious differences in reconstruction technique -- is the demo recon.m script a good place to start?

[estab1]

Thanks for your quick reply. recon.m contains all the parameters needed for the reconstruction. And grappa.m contains the algorithm itself. So you probably need both to understand the reconstruction technique.

[mckib2]

So I took a look at that recon.m and grappa.m and I think these are similar parameters to what's happening there:

• R=2
• 32 ACS lines
• 90% of FE lines used

The following modified basic_grappa.py script attempts to replicate these:

PULSAR_recon.py

```python '''Replicate PULSAR recon.''' import numpy as np import matplotlib.pyplot as plt from phantominator import shepp_logan from pygrappa import grappa if __name__ == '__main__': # Generate fake sensitivity maps: mps N = 128 ncoils = 4 xx = np.linspace(0, 1, N) x, y = np.meshgrid(xx, xx) mps = np.zeros((N, N, ncoils)) mps[..., 0] = x**2 mps[..., 1] = 1 - x**2 mps[..., 2] = y**2 mps[..., 3] = 1 - y**2 # generate 4 coil phantom ph = shepp_logan(N) imspace = ph[..., None]*mps imspace = imspace.astype('complex') ax = (0, 1) kspace = np.fft.fftshift(np.fft.fft2( np.fft.ifftshift(imspace, axes=ax), axes=ax), axes=ax) # crop 32x90% window from the center of k-space for calibration pd = 16 ctr = int(N/2) calib = kspace[ctr-pd:ctr+pd, int(.05*N):-int(.05*N), :].copy() # calibrate a kernel kernel_size = (7, 7) # undersample by a factor of 2 in ky kspace[:, 1::2, ...] = 0 # reconstruct: res = grappa( kspace, calib, kernel_size, coil_axis=-1, lamda=0.01) # replace calib #res[ctr-pd:ctr+pd, int(.05*N):-int(.05*N), :] = calib # SOS recon res = np.sqrt(np.sum(np.abs(np.fft.fftshift(np.fft.ifft2( np.fft.ifftshift(res, axes=ax), axes=ax), axes=ax))**2, axis=-1)) # SOS truth truth = np.sqrt(np.sum(np.abs(imspace)**2, axis=-1)) residual = np.abs(res - truth) plt.title('Residual | recon - truth |') plt.imshow(residual, vmin=0, vmax=residual.flatten().max(), cmap='gray') plt.colorbar() plt.show() ```

The results I'm getting seem alright -- I was not able to run PULSAR because I don't have easy access to MATLAB. Do you know if it runs in Octave? A couple notes:

• The unmodified basic_grappa.py is reconstructing a different scenario: R=4 (Rx=Ry=2) with a smaller calibration region (20x20)
• you are comparing the coil reconstructions to truth -- this in general will not look good because the SNR of any single coil is worse than the combined reconstruction. I recommend doing SOS on the reconstructed coil images for a good comparison
• try tuning the kernel size: larger kernels run slower but produce less residual, I found (7, 7) to be a good compromise
• replacing the calibration region into your reconstructed k-space will help (if the calibration region came from your collected data, that is)
• The quality of the coil sensitivities might also be in question -- I generated 4 very simple non-complex coil sensitivities. You might try experimenting with more/better coils

EDIT: screen shots for posterity

Cranking kernel size and replacing calibration region:

[estab1]

Thank you very much for your reply. I replicated your results with PULSAR_recon.py using my own phantom data and realized that PULSAR for matlab is producing the same results. I got confused about the results using basic_grappa.py because I forgot to add Gaussian noise to the k-space like I was doing in my matlab script. And thank you for the tip how to improve the image quality.