Split slice grappa with in plane acceleration

Author: Zimu Huo

Date: 05.2022

The slice grappa tries to reduce the total artefact without placing constraints on the intraslice and interslice artifacts. So the artefacts can be arbitrarily large(like vectors pointing in opposite directions can cancel out each other, also called artefact cancellation). The split slice trades interslice artefact with intraslice artefact (lower interslice artefact and higher intraslice artefact) by forming a "correlation matrix" that contrains the interslice leakage to be minimum.

---

References

[1] 
Author:  Stephen F. Cauley et al. 
Title: Interslice leakage artifact reduction technique for simultaneous multislice acquisitions
Link: https://pubmed.ncbi.nlm.nih.gov/23963964/

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.spsg import *

import numpy as np
from tqdm.notebook import tqdm

def spsg(dataR, calib,R = 1, kh = 5, kw = 5):
    # Split slice slice grappa 
    # Author: Zimu Huo
    [ny, nx, nc] = dataR.shape
    [cny, cnx, _, ns] = calib.shape
    ks = kh  * kw  * nc
    nt = cny * cnx
    data = np.zeros([ny, nx, nc, ns], dtype = complex)
    Ms = np.zeros([nt,ks,ns], dtype = complex)
    M = np.zeros([ks,ks], dtype = complex)
    for sli in range(ns):
        Ms[...,sli] = patches(calib[...,sli], kh, kw)
        M += Ms[...,sli].conj().T@Ms[...,sli]
    M = np.linalg.pinv(M)
    inMat = patches(dataR,kh, kw)
    for sli in (range(ns)):
        I = calib[...,sli].reshape(-1,nc)
        MI = Ms[...,sli].conj().T@I
        w = M@MI
        data[...,sli] = (inMat @ w).reshape(ny, nx, nc)
    return data

def patches(mat, kh, kw):
    [h, w, coil] = mat.shape
    kSize = kh * kw * coil
    inMatrix = np.zeros([h * w, kSize], dtype = complex)
    num = 0
    for y in range (h):
        ys = np.mod(np.linspace(y - int((kh/2-1)+ 1), y+int((kh/2)), kh, dtype = int), h)
        for x in range(w):
            xs = np.mod(np.linspace(x - int((kw/2-1) + 1), x+int((kw/2)), kw, dtype = int), w)
            inMatrix[num,:] = mat[ys][:,xs][:,:,:].reshape(1,-1)
            num = num + 1
    return inMatrix

numSlice = 4
R = 4
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)[...,:numSlice]
rawImage = ifft2c(data)
fovHeight, fovWidth, numCoil, _ = rawImage.shape

rawData = np.zeros(rawImage.shape, dtype = complex)
for sli in range (numSlice):
        rawData[:,:,:,sli] = fft2c(rawImage[:,:,:,sli])

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

shift = cycle*numAccq/(2*np.pi)
dataR = fft2c(simulate.multiSliceCAIPI(rawImage, cycle, R))

ncx = 32
ncy = 32
acsshift = cycle*int(numSlice* ncy /R)/(2*np.pi)
acsK = simulate.acs(rawData, (ncy, ncx))
acsIm = ifft2c(acsK)
calib = fft2c(simulate.singleSliceFov(acsIm,acsshift))

fig = plt.figure(figsize=(8, 6), dpi=80)
plt.subplot(121)
plt.imshow(np.abs(rsos(ifft2c(dataR))), cmap = "gray")
plt.subplot(122)
plt.imshow(np.abs(stitch(rsos(ifft2c(calib)),1)), cmap = "gray")

<matplotlib.image.AxesImage at 0x7fb005a94f70>

recon = spsg(dataR, calib)

fig = plt.figure(figsize=(8, 6), dpi=80)
plt.subplot(211)
plt.imshow(np.abs(stitch(rsos(rawImage), 1)), cmap = "gray")
plt.subplot(212)
plt.imshow(np.abs(stitch(rsos(ifft2c(recon)),1)), cmap = "gray")

<matplotlib.image.AxesImage at 0x7fb005c69940>

reconshift = cycle*int(numSlice* fovHeight/R)/(2*np.pi) 
recon1 = simulate.singleSliceFov(ifft2c(recon), - reconshift)
fig = plt.figure(figsize=(8, 6), dpi=80)
plt.subplot(211)
plt.imshow(np.abs(stitch(rsos(rawImage), 1)), cmap = "gray")
plt.subplot(212)
plt.imshow(np.abs(stitch(rsos((recon1)),1)), cmap = "gray")

<matplotlib.image.AxesImage at 0x7fb006169b80>