2D Echo Planar Imaging for 128 x 128 image matrix. TR = 4000 ms, TE = 72 ms, FA = 90°. The number of subvoxels is 1×1×4. The total number of subvoxels is 13,572,096. The calculation time was about 1 s.

No magnetic field inhomogeneity. Amplitude (left) and phase (right) images.

 

No magnetic field inhomogeneity with 100 Hz field offset. Amplitude (left) and phase (right) images.

 

Magnetic filed inhomogeneity ∝ x. Amplitude (left) and phase (right) images.

 

Magnetic filed inhomogeneity ∝ y. Amplitude (left) and phase (right) images.

 

Magnetic filed inhomogeneity ∝ xy. Amplitude (left) and phase (right) images.

 

Magnetic filed inhomogeneity ∝ x2-y2. Amplitude (left) and phase (right) images.

 

Pulse sequence chart visualized by the SequenceViewer:

 

Python sequence code:

from psdk import *
import numpy as np

gamma = 42.57747892 # [MHz/T]
TR = 4000.0e+3 # [us]
TE = 6.0e+3 # [us]
NR = 128 # Number of readout points
NPE1 = 128 # Number of 1st phase encoding
fov = [220.0, 220.0, 256.0] # [mm]
dwell_time = 5.0 # [us]
slice_width = 6.0 # [mm]
gx_rise_time = 180.0 # [us]
gy_rise_time = 180.0 # [us]
gz_rise_time = 180.0 # [us]
gx_value = 1e+6 / (dwell_time * gamma * fov[0]) # [mT/m]
# gy_value = 2e+6 / (dwell_time * gamma * fov[1]) * NPE1 / NR # [mT/m]
gy_value = 1e+6 / (gy_rise_time * gamma * fov[1]) # [mT/m]
gz_value = 1.25 / (slice_width * 1.0e-3) / gamma # [mT/m]
excitation_pulse_width = 3200.0 # [us]
excitation_pulse_flip_angle = 90.0 # [degree]

def sinc_with_hamming(flip_angle, pulse_width, points, *, min=-2.0*np.pi, max=2.0*np.pi):
    x0 = np.arange(min, max, (max - min) / points)
    x1 = x0 + (max - min) / points
    y = (np.sinc(x0 / np.pi) + np.sinc(x1 / np.pi)) * 0.5 * np.hamming(points)
    return flip_angle * y * points / (y.sum() * pulse_width * 360.0e-6 * gamma)

with Sequence('2D Gradient echo EPI'):

    with Block('Excitation', excitation_pulse_width + 2.0*gz_rise_time):
        GZ(0.0, gz_value, gz_rise_time)
        RF(gz_rise_time, sinc_with_hamming(excitation_pulse_flip_angle, excitation_pulse_width, 160), excitation_pulse_width / 160)
        GZ(excitation_pulse_width + gz_rise_time, 0.0, gz_rise_time)

    with Block('PhaseEncoding', gy_rise_time * 16.0 + gy_rise_time * 2.0):
        GX(gy_rise_time * 18.0 - gx_rise_time * 2.5 - NR // 2 * dwell_time, -gx_value, gx_rise_time)
        GX(gy_rise_time * 18.0 - gx_rise_time * 2.0, gx_value, gx_rise_time * 2.0)
        GY(0.0, -gy_value * 4.0, gy_rise_time)
        GY(gy_rise_time * 16.0, 0.0, gy_rise_time)        
        GZ(0.0, -gz_value, gz_rise_time)
        GZ(excitation_pulse_width * 0.5, 0.0, gz_rise_time)
        
    with Block('Readout', NR * dwell_time):
        AD(0.0, NR, dwell_time)
 
    with Block('Reverse', gx_rise_time * 2.0):
        GX(0.0, ([-gx_value, gx_value], ['Echo']), gx_rise_time * 2.0)
        GY(0.0, gy_value, gy_rise_time)
        GY(gy_rise_time, 0.0, gy_rise_time)
        
    with Block('Rewinding', NR // 2 * dwell_time + gx_rise_time):
        GX(NR // 2 * dwell_time - gx_rise_time * 0.5, 0.0, gx_rise_time)
        
    with Main():
        BlockRef('Excitation')
        WaitUntil(TE + excitation_pulse_width * 0.5 + gz_rise_time - gy_rise_time * 18.0)
        BlockRef('PhaseEncoding')
        with Loop('Echo', 127):
            BlockRef('Readout')
            BlockRef('Reverse')
        BlockRef('Readout')
        BlockRef('Rewinding')
        WaitUntil(TR)