cgSENSE reconstruction

Author: Zimu Huo

Date: 03.2022

Key points

• The main ideas are from [3] where he talked about iterative methods and built on [1] the original paper on SENSE

• As stated by the title, this algorithm uses conjugate gradient method, and it is implemented using equations from 45 to 49 in [4].

• Note that the sum of k-space integrals are weighted by complex conjugate coil sensitivity as stated in [3] equation 18

Formulation

The image is acquired with multiple receiver coils denoted by sensitivity encoding operator S. The Fourier encoding is represented as a discrete Fourier encoding matrix F, where each row captures the positive or negative frequency attributes within a signal.

\begin{equation} \begin{aligned} \min_{x} \quad & \quad & ||DFSx - y||^2 + \lambda ||x||_2 \end{aligned} \end{equation}

This can be solved using normal equation, but the inverse of system matrix DFS is computationally demanding, the problem is typically solved in an iterative fashion. A numerically fast method to solve problem subject to L2 norm is given by the conjugate gradient algorithm.

---

References

[1] 
Author: Klaas P. Pruessmann et al. 
Title: SENSE: Sensitivity Encoding for Fast MRI
Link: https://pubmed.ncbi.nlm.nih.gov/10542355/

[2] 
Author: Klaas P. Pruessmann et al. 
Title: 2D SENSE for faster 3D MRI*
Link: https://doc.rero.ch/record/317765/files/10334_2007_Article_BF02668182.pdf

[3] 
Author: Klaas P. Pruessmann et al. 
Title: Advances in Sensitivity Encoding With Arbitrary
        k-Space Trajectories
Link: https://onlinelibrary.wiley.com/doi/pdfdirect/10.1002/mrm.1241

[4] 
Author: Jonathan Richard Shewchuk. (what an absolute legend!)
Title: An Introduction to the Conjugate Gradient Method
      Without the Agonizing Pain Edition 1+1/4
Link: https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf

[5] 
Author: Oliver Maier et al.
Title: CG-SENSE revisited: Results from the first ISMRM reproducibility challenge
Link: https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.28569

import sys
sys.path.insert(1, '../')
import matplotlib.pyplot as plt
import util.coil as coil
from util.fft import *
import numpy as np
import util.mask as undersample
from util.cgSolver import *

data = np.load("../lib/slice1_grappa1.npy")
# estimate the sensivity maps
[ny, nx, nc] = data.shape
images = ifft2c(data) 
image = coil.rsos(images) 
sensMap = coil.cmap(images)
R = 4
mask = undersample.sense(data.shape,R)
dataR = np.zeros([ny, nx, nc], dtype=complex)
dataR = data*mask
imagesR = ifft2c(dataR)

# display all channels
show(ifft2c(data),5)

import numpy as np
from util.fft import *
from tqdm.notebook import tqdm
import time

def cg_sense(dataR, sensMap, numIter = 50):
    mask = np.where(dataR == 0, 0, 1)
    imagesR = ifft2c(dataR)
    [height, width, coil] = imagesR.shape
    sconj = np.conj(sensMap)
    B = np.sum(imagesR*sconj,axis = 2)
    B = B.flatten()
    x = 0*B
    r = B 
    d = r 
    for j in range(numIter):
        Ad = np.zeros([height,width],dtype = complex)
        for i in range(coil):  
            Ad += ifft2c(fft2c(d.reshape([height,width])*sensMap[:,:,i])*mask[:,:,i])*sconj[:,:,i] 
            # this is intuitive, the correct image is sensivity encoded, then undersampled in K space, but this took me quite a while :(
        # see equations from 45 to 49 in [4].
        Ad = Ad.flatten()
        a = np.dot(r,r)/(np.dot(d,Ad))
        '''
        This is the core idea behind steepest descent, where the new direction is orthogonal to the previous search. 
        This also intuitive, in some sense you only maximise the information obatined in this search when it can no longer 
        improve your next search. However, this will lead to zigzag-ish path, which is highly inefficient.     
        '''
        x = x + a*d
        '''
        This corresponds to a collection of images, like layers upon layers staring from the initial guess. 
        '''
        rn = r - a*Ad
        beta = np.dot(rn,rn)/np.dot(r,r)
        r=rn
        d = r + beta*d
        '''
        This is where conjugate gradient kicks in, it uses the previous hardwork along with the current search direction. 
        As a result, this is more efficient. Like climbing a hill but carries some weight, or having inertia. 
        This idea is also quite commonly used in machine learning -> momentum
        '''
    return x.reshape([height,width])

recon =cg_sense(dataR, sensMap, 100)

plt.figure(figsize=(16, 12), dpi=80)
plt.subplot(121)
plt.title("raw data")
tf = plt.imshow(np.abs(coil.rsos(imagesR)), cmap ="gray")
plt.axis("off")
plt.subplot(122)
plt.title("cg sense")
tf = plt.imshow(np.abs(recon), cmap ="gray")
plt.axis("off")
plt.show()