EPI phase correction Low rank¶
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 entropy-based method attempts to minimize the image entropy of the magnitude image to minimize EPI misalignment artifacts
References
[1]
Author: Heid, O
Title: Method for the phase correction of nuclear magnetic resonance signals
Link: https://patents.google.com/patent/US5581184A/en
[2]
Author: Skare S, Clayton DB, Newbould R, Moseley M, and Bammer R.
Title: A fast and robust minimum entropy based non-interactive Nyquist ghost correction algorithm. In Proceedings of ISMRM. 2006;p.2349.
Link: https://cds.ismrm.org/protected/06MProceedings/PDFfiles/02349.pdfIn [1]:
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 *
In [2]:
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
In [3]:
tissues,tissuetype = get_tissue_images()
ny, nx, nz, nt = tissues.shape
nc = 4 # number of coils
In [4]:
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)
In [5]:
# plot the simulated ground truth image
show(rsos(images))
In [6]:
# helper functions
def entropy(data):
if len(data.shape) > 2:
image = np.sum(np.abs(data)**2, -1)
else:
image = np.abs(data)**2
b = np.sqrt(image)/np.sum(image)
return np.abs(np.sum(b/np.log(b)))
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.
In [20]:
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()
In [21]:
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_entropy(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_entropy_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
In [22]:
recon, phase_ramp, A, B, Z1 = epi_phasecorrection_entropy(data, grid = 50, output ="full")
here we plot the reconstruction results
In [23]:
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:
In [25]:
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
In [29]:
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)
In [42]:
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()
In [43]:
recon, phase_ramp, A, B, Z2 = epi_phasecorrection_entropy(data, grid = 50, output ="full")
In [44]:
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()
In [45]:
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()
In [ ]: