Conjugate synthesis with phase correction

Author: Zimu Huo

Date: 05.2022

Mainly based on Prof. John Pauly's diagram.

---

References

[1] 
Author: G McGibney et al.
Title: Quantitative evaluation of several partial Fourier reconstruction algorithms used in MRI
Link: https://pubmed.ncbi.nlm.nih.gov/8371675/

[2]
Author: John Pauly 
Title: Partial k-Space Recconstruction
Link: https://ece-classes.usc.edu/ee591/library/Pauly-PartialKspace.pdf

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
import util.simulator as simulate
from util.zpad import *
from tqdm.notebook import tqdm
from util.partialFourier import *

def conjugateSynthesis(dataR, calib):
    idx = np.where(np.sum(np.abs(dataR), axis=1)==0)[0]
    cphs = np.exp(-1j*np.angle(ifft2c(calib)))
    im = np.zeros(dataR.shape, dtype = complex)
    imR =  ifft2c(dataR) * cphs
    dataR = fft2c(imR)
    dataR[idx[0]:,:] = np.flip(np.flip(np.conj(dataR[len(idx):0:-1,:]),0),1)
    return ifft2c(dataR)

data = np.load("../lib/slice1_grappa1.npy")
rawImage = ifft2c(data)
[ny,nx,nc] = data.shape

dataR = data * undersample.partialFourier(data.shape, 0.7)
calib = simulate.acs(dataR,(64, 256))
calib = zpad3(calib)

recon = np.zeros(dataR.shape, dtype = complex)
for c in tqdm(range(nc)):
    recon[...,c] = conjugateSynthesis(dataR[...,c], calib[...,c])

  0%|          | 0/34 [00:00<?, ?it/s]

plt.figure(figsize=(12, 8), dpi=80)
plt.subplot(121)
plt.title("raw data")
plt.imshow(np.log(np.abs(coil.rsos((dataR)))),cmap='gray')
plt.subplot(122)
plt.title("conjuagte synthesis")
plt.imshow(np.abs(coil.rsos(recon.real)),cmap = 'gray')
plt.show()

C:\Users\Zimu\AppData\Local\Temp\ipykernel_9148\3236454631.py:4: RuntimeWarning: divide by zero encountered in log
  plt.imshow(np.log(np.abs(coil.rsos((dataR)))),cmap='gray')