EPI blip-up blip-down correction

Author: Zimu Huo

Date: 07.2022

use blip-up-blip-down image pairs to estimate the b0 field maps, then perform the distortion correction based on the field maps generated.

Here I used third party package (FSL). Just to double check my simulation is correct.

---

References

[1] 
Author:  M.W. Woolrich et al.
Title:  Bayesian analysis of neuroimaging data in FSL
Link: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki

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 nibabel as nib
import os
from tqdm.notebook import tqdm

import scipy.io
def get_tissue_images(slices):
    tissuetype = ['graymatter', 'deep_graymatter', 'whitematter', 'csf']
    T2 = [110, 100, 60, 1500]
    T2s = [40, 45, 50, 1000]
    mat = scipy.io.loadmat('../lib/resource/data/tissue_images/tissue_images.mat')
    tissues = mat.get("tissue_images")[:,:,slices,:]
    return np.squeeze(tissues), tissuetype

def brain_tissue(slices, coils):
    tissues,tissuetype = get_tissue_images(slices)
    if (coils== 1): return tissues,tissuetype, T2, T2s
    tissues = np.repeat(tissues[:, :,  np.newaxis,:], coils, axis=-2)
    coils = generate_birdcage_sensitivities(matrix_size = 222,number_of_coils = coils)
    tissues = tissues *  np.repeat(coils[...,np.newaxis], 4, axis = -1)
    return tissues,tissuetype
def sim_blip_up(ideal_image, field):
    ny, nx = ideal_image.shape
    image = np.zeros([ny, nx], dtype = complex)
    for y in range(ny-1,0,-1):
        modulation = np.exp(1j * t * field[ny-y])
        modulation = np.repeat(modulation[:,None], nx, -1)
        tmp = ideal_image * modulation 
        image[y] = fft2c(tmp)[y]
    extended = ifft2c(image)
    return extended
def sim_blip_down(ideal_image, field):
    ny, nx = ideal_image.shape
    image = np.zeros([ny, nx], dtype = complex)
    for y in range(ny):
        modulation = np.exp(1j * t * field[y])
        modulation = np.repeat(modulation[:,None], nx, -1)
        tmp = ideal_image * modulation 
        image[y] = fft2c(tmp)[y]
    contracted = ifft2c(image)
    return contracted

tissues,tissuetype = brain_tissue(70, 4)
ny, nx, nc, nt = tissues.shape
ns = 144

blip_up_data = np.zeros([ny,nx,ns])
blip_down_data = np.zeros([ny,nx,ns])
for sli in tqdm(range(ns)):
    tissues,tissuetype = brain_tissue(sli, 4)
    TE = 100
    ideal_image = np.zeros([ny, nx, nc], dtype = complex)
    for t in range(nt):
        ideal_image += tissues[...,t] * np.exp(TE/t2(tissuetype[t]))
    ideal_image = rsos(ideal_image) 
    G = 0.003
    fov = ny
    dwell_time = 1
    gamma = 2* np.pi * 42.58
    t = np.linspace(-1, 1, nx) * 1E-3 * fov
    phase_inc = gamma * G * dwell_time
    field = np.zeros([ny, nx], dtype = complex)
    for y in range(ny):
        field[y] =  phase_inc * y 
    blip_down_data[...,sli] = np.abs(sim_blip_down(ideal_image, field))
    blip_up_data[...,sli] = np.abs(sim_blip_up(ideal_image, field))

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

for sli in range(ns):
    blip_up_data[...,sli] = blip_up_data[...,sli].T

for sli in range(ns):
    blip_down_data[...,sli] = blip_down_data[...,sli].T

path = "../lib/EPI_DISTORTION/"
if os.path.isdir(path) is False:
    os.mkdir(path)
%cd  ../lib/EPI_DISTORTION/

/home/zimu/Desktop/mri_recon/lib/EPI_DISTORTION

corrected = nib.load('merged.nii.gz').get_fdata()
print(corrected.shape)

(144, 222, 145, 2)

new_image = nib.Nifti1Image(blip_up_data, affine=np.eye(4))
nib.save(new_image,"blip_up.nii")
new_image = nib.Nifti1Image(blip_down_data, affine=np.eye(4))
nib.save(new_image, "blip_down.nii")
with open('acqparams.txt', 'w') as f:
    f.write('0 -1 0 '+str(10E-6*(ny-1)*nx)[:5])
    f.write('\n')
    f.write('0 1 0 '+str(10E-6*(ny-1)*nx)[:5])

merged_file = "merged"

os.system("fslmerge -t "+merged_file+".nii.gz blip_down.nii blip_up.nii")

0

os.system("topup --imain="+merged_file+" --datain=acqparams.txt --config=b02b0.cnf --out=topup_bfield")

0

os.system("applytopup --imain=blip_down,blip_up --topup=topup_bfield --datain=acqparams.txt --inindex=1,2 --out=fsl_corrected")

0

fieldmap = nib.load('topup_bfield_fieldcoef.nii.gz').get_fdata()
show(fieldmap)

corrected = nib.load('fsl_corrected.nii.gz').get_fdata()
show(corrected[...,70])