Simultaneous multislice SENSE/GRAPPA

Author: Zimu Huo

Date: 05.2022

In this simulation of the SMS experiment, we concatenate multiple slices into an extended field of view. Techniques like SENSE and GRAPPA treat this extended field as a single image, applying the reconstruction in their original, unmodified form. However, GRAPPA techniques encounter difficulties in this scenario due to the CAIPI shift causing a phase discontinuity between the images. This sharp change cannot be captured by a small GRAPPA kernel. Since the GRAPPA kernel is small and bandlimited, its image domain representation is smooth and slowly varying, making it inadequate for capturing the sharp phase change.

---

References

[1] 
Author: Blaimer M et al. 
Title: Accelerated volumetric MRI with a SENSE/GRAPPA combination
Link: https://onlinelibrary.wiley.com/doi/10.1002/jmri.20632

[2] 
Author: Markus Barth et al. 
Title: Simultaneous multislice (SMS) imaging techniques
Link: https://onlinelibrary.wiley.com/doi/10.1002/mrm.25897#mrm25897-bib-0064

import sys
sys.path.insert(1, '../')
import matplotlib.pyplot as plt
import numpy as np
from scipy import io
import util.simulator as simulate
from util.coil import *
from util.fft import *
import util.mask as undersample
from util.sg import *
from tqdm.notebook import tqdm
from util.grappa import * 

slice1 = np.load("../lib/slice1_grappa1.npy")
slice2 = np.load("../lib/slice2_grappa1.npy")
slice3 = np.load("../lib/slice3_grappa1.npy")
slice4 = np.load("../lib/slice4_grappa1.npy")
data = np.concatenate((slice1[...,None], slice2[...,None], slice3[...,None], slice4[...,None]), -1)

# the coil number is high, so simulation with more than 2 slices will be long 
numSlice = 2
R = 4
data = data[...,:numSlice]
rawImage = ifft2c(data)
fovHeight, fovWidth, numCoil, _ = rawImage.shape

rawData= fft2c(rawImage)
acsK = simulate.acs(rawData, (32, 32))

# input data 
show(rsos(ifft2c(rawData), -2))
# calibration signal 
show(rsos(ifft2c(acsK), -2))

# optimal caipi shift
cycle = np.arange(0,1,1/numSlice) * 2* np.pi
numAccq = int(numSlice*fovHeight/R)
shift = cycle*numAccq/(2*np.pi)

dataR = np.zeros([fovHeight * numSlice, fovWidth, numCoil], dtype = complex)
dataR[::R] = fft2c(simulate.multiSliceCAIPI(rawImage, cycle, R))
acsIm  = ifft2c(acsK)
acsshift = cycle*(32*4/8)/(2*np.pi)
acsImLarge = simulate.multiSliceFov(acsIm, acsshift)
acs = fft2c(acsImLarge)

# undersampled with slice factor 
show(rsos(ifft2c(dataR)).T)

recon = grappa(dataR, acs, combine = False)

1934.4000735431273

# We need to unapply the CAIPI shift to the reconstructed image, but I skipped this step in this case. The key point is to demonstrate that GRAPPA cannot capture sharp phase changes.
fig = plt.figure(figsize=(16, 12), dpi=80)
plt.title("reconstruction",fontsize = 20 )
plt.imshow(np.abs(rsos(recon)).T, cmap = "gray")
plt.axis("off")

(-0.5, 511.5, 255.5, -0.5)