EPI phase correction

Author: Zimu Huo

Date: 05.2022

Due to the alternating direction of the EPI readout lines, hardware imperfections such as timing delays, eddy currents, and gradient coil heating can cause misalignment of the forward and reverse lines in k-space. This misalignment manifests itself in the images as a Nyquist (N/2) ghost in the phase encode direction and sinusoidal modulation of the object in the frequency encode direction.

The low rank method attempts to promote self consistency in k-space to minimize EPI misalignment artifacts

---

References

[1] 
Author: Peterson E,	Aksoy M, Maclaren J, and Bammer	R.	
Title: Acquisition-free	Nyquist	ghost correction for parallel imaging accelerated EPI. In Proceedings of ISMRM.	2015;p.0075.
Link: https://patents.google.com/patent/US5581184A/en

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

import scipy.io
def get_tissue_images(): # helper function to create realistic brain phantom
    tissuetype = ['graymatter', 'deep_graymatter', 'whitematter', 'csf']
    T2 = [110, 100, 60, 1500]
    T2s = [40, 45, 50, 1000]
    mat = scipy.io.loadmat('../lib/tissue_images.mat')
    tissues = mat.get("tissue_images")[:,:,:,:]
    return np.squeeze(tissues), tissuetype

tissues,tissuetype = get_tissue_images()
ny, nx, nz, nt = tissues.shape
nc = 4 # number of coils 

TE = 100
ideal_image = np.zeros([ny, nx, nz], dtype = complex)
for t in range(nt):
    ideal_image += tissues[...,t] * np.exp(TE/t2(tissuetype[t]))
images = ideal_image[:,:,nz//2] # we are only interested in one slice 
data = fft2c(images)

coils = generate_birdcage_sensitivities(matrix_size = ny, number_of_coils = nc)
images = np.repeat(images[...,np.newaxis],  nc, axis = -1) * coils
data = fft2c(images)

# plot the simulated ground truth image 
show(rsos(images))

# helper functions 
def sake_forward(data, k):
    [ny, nx, nc] = data.shape
    data = np.moveaxis(data, -1, 0)
    mat = np.zeros([(ny-k+1)*(nx-k+1), k * k * nc], dtype = complex)
    idx = 0
    for y in range(max(1, ny - k + 1)):
        for x in range(max(1, nx - k + 1)):
            mat[idx, :] = data[:,y:y+k, x:x+k].reshape(1,-1)
            idx += 1
    return mat
    
def epi_lowrank_cost(p, data):
    [ny, nx, nc] = data.shape
    x = np.linspace(-0.5, 0.5, nx)
    obj = object_projection_from_central_line(data)
    model = p[0] * x + p[1]
    tvec = polyfit(x, model, w = obj, order = 1)
    ghost_data = fft2c(epi_phase_correction_even_odd(data,tvec))
    mat = sake_forward(ghost_data, k = 3)
    U, S, VT = np.linalg.svd(mat,full_matrices=False)
    return np.sum(S[1:])
    
    
def epi_phase_correction_even_odd(data, tvec):
    [ny, nx, nc] = data.shape
    tvecm = np.exp(1j*tvec).reshape(-1,1)
    even = np.zeros(data.shape, dtype = complex)
    odd = np.zeros(data.shape, dtype = complex)
    even[::2] = data[::2]
    odd[1::2] = data[1::2]
    dataflag = np.zeros([nx, 1])
    dataflag[1::2] = 1
    odd = apply_epi_phase(odd, tvecm)
    ghost_data = ifft2c(even) + ifft2c(odd)
    return ghost_data
    
def object_projection_from_central_line(ref):
    ref = ref[ref.shape[0]//2-1]
    obj = (np.sqrt(np.sum(np.abs(ifft(ref))**2,1)))
    obj *= 1/obj.max()
    return obj
def polyfit(x, y, w, order):
    coef = np.polyfit(x,y,order, w=w)
    phi = np.poly1d(coef) 
    return phi(x.reshape(-1,1)).flatten()
    
def apply_epi_phase(data, tvecm):
    recon = np.zeros(data.shape, dtype = complex)
    tvec = np.ones(data.shape, dtype = complex)
    tvec[1::2,...] = np.tile(tvecm.reshape(-1,1), data.shape[2])
    data = ifft(data,1)
    recon = data * tvec
    recon = ifft(recon,0)
    return fft2c(recon) 

Hardware imperfections can cause a slight deviation between otherwise equivalent points along the frequency encoding axis between even and odd echoes in the EPI readout. Such misalignment of the echoes in K-space leads to an artifact known as Nyquist ghosting. This artifact manifests as a superposition of the true image and a copy shifted by half of the field of view (FOV) in the phase encoding direction. Additionally, a cosine modulation of the signal intensities occurs in the frequency encoding direction.

The misalignment between the forward and reverse kx lines is typically modeled as a linear phase difference in the hybrid space (after Fourier transforming the data along the frequency encoding axis) as follows:

\[ S(x,k_{y_{odd}}) = e^{i \left(\frac{\pi \alpha}{N_x}x + \beta\right)} \]

where:

• ( x ) is the index of the position in the readout direction of total size ( N_x ),
• ( a ) is the slope of the phase due to the constant time echo delay,
• ( b ) is a constant off-resonance phase shift between the even and odd readouts in K-space.

a = 3.00
b = -0.5

x = np.linspace(-0.5, 0.5, nx)
tvecm = np.exp(1j* (a*x+b)) 
data = apply_epi_phase(data, tvecm)
plt.figure()
plt.imshow(rsos(ifft2c(data)), cmap = "gray")
plt.axis('off')
plt.show()

def epi_entropy_cost(params, data):
    # This function applies a phase correction using the a and b on the input data returns the entropy of the correted image 
    # Author: Zimu Huo 
    a, b = params
    [ny, nx, nc] = data.shape
    x = np.linspace(-0.5, 0.5, nx)
    obj = object_projection_from_central_line(data)
    model = a * x + b
    tvec = polyfit(x, model, w=obj, order=1)  # A linear fit, weighted by the 1D projection of the object, is not relevant in reference-less methods.
    ghost_data = epi_phase_correction_even_odd(data, tvec)
    return entropy(ghost_data)

def epi_phasecorrection_lowrank(data, grid=50, output="partial"): 
    # I implemented a brute force algorithm to search through all combinations of a's and b's, and return the global minimum.
    # Author: Zimu Huo 
    ranges = ((-5, 5), (-1, 1))
    p = scipy.optimize.brute(epi_lowrank_cost, ranges, args=(data,), Ns=grid, workers=1, full_output=True)
    
    optimal_params = p[0]
    cost_values = p[3]

    [ny, nx, nc] = data.shape
    x = np.linspace(-0.5, 0.5, nx)
    obj = object_projection_from_central_line(data)
    model = optimal_params[0] * x + optimal_params[1]
    tvec = polyfit(x, model, w=obj, order=1)
    recon = epi_phase_correction_even_odd(data, tvec)

    a_range = np.linspace(ranges[0][0], ranges[0][1], grid)
    b_range = np.linspace(ranges[1][0], ranges[1][1], grid)
    A, B = np.meshgrid(a_range, b_range)
    Z = cost_values.reshape(grid, grid)


    if output == "full":
        return fft2c(recon), np.exp(1j * tvec).reshape(-1, 1), A,B, Z
    else:
        return fft2c(recon), Z

recon, phase_ramp, A, B, Z1 = epi_phasecorrection_lowrank(data, grid = 50, output ="full")

here we plot the reconstruction results

plt.figure(figsize =(16,12))
plt.subplot(121)
plt.title("raw")
plt.axis('off')
plt.imshow((np.abs(rsos(ifft2c(data)))), cmap ="gray")
plt.subplot(122)
plt.axis('off')
plt.title("recon")
plt.imshow((np.abs(rsos(ifft2c(recon)))), cmap ="gray")
plt.show()

The cost contour plot:

plt.figure()
plt.contourf(A, B, Z1, levels=50, cmap='jet')
plt.colorbar(label='Entropy Cost')
plt.xlabel('a')
plt.ylabel('b')
plt.title('Cost Contour Plot')
plt.show()

Lets add some noise and see what happens

def gaussian_noise(shape, L = None, mu = 0, sigma = 1):
    # generate noise, L is the correlation matrix (noise colors), can also be used to pre-whitten the noise
    [ny, nx, nc] = shape
    n = np.zeros([ny * nx, nc], dtype = complex)
    n.real = np.random.normal(mu, sigma, ny*nx*nc).reshape(ny * nx, nc)
    n.imag = np.random.normal(mu, sigma, ny*nx*nc).reshape(ny * nx, nc)
    if L is not None: 
        n = n@L
    return n.reshape(ny, nx, nc)

a = 3.00
b = -0.5

data = fft2c(images)
data = data + gaussian_noise(data.shape, sigma = 0.004)
x = np.linspace(-0.5, 0.5, nx)
tvecm = np.exp(1j* (a*x+b)) 
data = apply_epi_phase(data, tvecm)
plt.figure()
plt.imshow(rsos(ifft2c(data)), cmap = "gray")
plt.axis('off')
plt.show()

recon, phase_ramp, A, B, Z2 = epi_phasecorrection_lowrank(data, grid = 50, output ="full")

plt.figure(figsize =(16,12))
plt.subplot(121)
plt.title("raw")
plt.axis('off')
plt.imshow((np.abs(rsos(ifft2c(data)))), cmap ="gray")
plt.subplot(122)
plt.axis('off')
plt.title("recon")
plt.imshow((np.abs(rsos(ifft2c(recon)))), cmap ="gray")
plt.show()

plt.figure()
plt.contourf(A, B, Z2, levels=50, cmap='jet')
plt.colorbar(label='Entropy Cost')
plt.xlabel('a')
plt.ylabel('b')
plt.title('Cost Contour Plot')
plt.show()