How to reconstruct image from an Under Sampled MRI k-space data?

Question

Hardware

We have a homemade MRI machine (from the open source project openimaging). It has a single coil and B0 = 0.5 Tesla (enough for teaching).

Sampling Mask

We used a randomly generated sampling mask for a general single slice Spin Echo imaging:

In the sampling mask, each white line represents a phase encoding step. The other lines were zero-padded. There are 38 white lines and 90 black lines in this mask. The sampling matrice is 128 * 256 (if fully sampled).

Reconstruction Algorithm

We've referenced the open source reconstruction code, based on sigpy, to reconstruct. Here are the code and data we used:

import numpy as np
import sigpy as sp
import sigpy.mri as mr
from scipy.io import loadmat
from matplotlib.pyplot import *


k = loadmat(r"Comp121_data_psf0.3.mat")["data"]
ksp = np.expand_dims(k,axis=0) # Dim:1 * 128 *256
img_shape = k.shape
#%% FFT
F = sp.linop.FFT(ksp.shape, axes=(-1,-2))
I0 = np.squeeze(abs(F.H * ksp))

#%% CS
mask = loadmat(r"Comp121_data_psf0.3.mat")["mask"]
P = sp.linop.Multiply(ksp.shape, mask)

mps = np.ones_like(abs(ksp))
S = sp.linop.Multiply(img_shape, mps)

W = sp.linop.Wavelet(img_shape)
wav = W *S.H* F.H * ksp

A = P * F * S * W.H

lamda = 0.0005
proxg = sp.prox.L1Reg(wav.shape, lamda)
alpha = 1
wav_thresh = proxg(alpha, wav)

# optimizing
max_iter = 10
alpha = 0.001

def gradf(x):
    return A.H * (A * x - ksp)

wav_hat = np.zeros(wav.shape, complex)
alg = sp.alg.GradientMethod(gradf, wav_hat, alpha, proxg=proxg, max_iter=max_iter)

while not alg.done():
    alg.update()
    print('\rL1WaveletRecon, Iteration={}'.format(alg.iter), end='')

I_cs = abs(W.H(wav_hat))
imshow(np.hstack((I0/np.max(I0),I_cs/np.max(I_cs))),cmap='gray')
title(f'padzero [left] and L1wav [right]')
figure()
imshow(abs(I0/np.max(I0)-I_cs/np.max(I_cs)))
colorbar()
show()

Problems

In my expectation, the reconstructed image from zero-padding would be worse. But the two images (zero-padding v.s. L1 wavelet regularized reconstruction) seem to be similar, both have a visible phase ghost. L1 wavelet regularization does not look superior to zero-padding in my situation.

Changing lambda, alpha, and max_iter has no effect, the ghost still exists. Changing the fraction of sampled phase encoding lines from 0.3 to 0.7, the ghost still exists.

My questions:

1. Are there any mistakes in my reconstruction code?
2. How to set the parameters: lambda, alpha, and max iter? Are there any standards or references for this?
3. How to remove the ghost in the reconstructed image?

Answer 1

I was experiencing a similar problem (CS reco looked the same as the simple FT reco). What helped me was to normalize my data to max=1,

ksp_norm = ksp  / np.max(np.max(np.abs(ksp)))

in your case. Maybe that helps.