3D RF spoiled gradient echo imaging

Image matrix = 256 x 256 x 32, TR = 18 ms, TE = 10 ms, flip angle = 30 degree, subvoxel = s x 1 x 1 (s = 1 – 16), calculation time = 771.94 s for sv = 16, number of subvoxels = 54,577,152 for sv = 16.

Amplitude of artifact vs number of subvoxels. Artifact has local minima for s = power of 2.

Python sequence code:

from psdk import *
import numpy as np

gamma = 42.57747892 # [MHz/T]
TR = 18.0e+3 # [us]
TE = 10.0e+3 # [us]
NR = 256 # Number of readout points
NPE1 = 256 # Number of 1st phase encoding
NPE2 = 32 # Number of 2nd phase encoding
fov = [256.0, 256.0, 256.0] # [mm]
dwell_time = 10.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]
gz_value = 2e+6 / (dwell_time * gamma * fov[2]) * NPE2 / NR # [mT/m]
gx_rise_time = 300.0 # [us]
gy_rise_time = 300.0 # [us]
gz_rise_time = 300.0 # [us]
excitation_pulse_width = 100.0 # [us]
excitation_pulse_flip_angle = 30.0 # [degree]
excitation_pulse_value = excitation_pulse_flip_angle / (360.0e-6 * gamma * excitation_pulse_width) # [uT]

def phase_shift_angle(i):
    phi = i * (i + 1) / 2 * 117 / 360
    phi -= round(phi)
    return 2.0 * np.pi * phi

with Sequence('3D SPGR'):

    with Block('Excitation', excitation_pulse_width):
        RF(0.0, [excitation_pulse_value] * 10, excitation_pulse_width / 10, phase=([phase_shift_angle(i) for i in range(NPE1 * NPE2)], ['PE2', 'PE1']))

    with Block('PhaseEncoding', NR // 2 * dwell_time + gx_rise_time * 2.5):
        GX(0.0, -gx_value, gx_rise_time)
        GY(0.0, ([gy_value * (i - NPE1 // 2) / NPE1 for i in range(NPE1)], ['PE1']), gy_rise_time)
        GZ(0.0, ([gz_value * (i - NPE2 // 2) / NPE2 for i in range(NPE2)], ['PE2']), gy_rise_time)
        GY(NR // 2 * dwell_time, 0.0, gy_rise_time)
        GZ(NR // 2 * dwell_time, 0.0, gz_rise_time)
        GX(NR // 2 * dwell_time + gx_rise_time * 0.5, gx_value, gx_rise_time * 2.0)
        
    with Block('Readout', NR * dwell_time):
        AD(0.0, NR, dwell_time, phase=([phase_shift_angle(i) for i in range(NPE1 * NPE2)], ['PE2', 'PE1']))
        
    with Block('Rewinding', NR // 2 * dwell_time + gx_rise_time * 2.5):
        GY(gx_rise_time * 2.5 - gy_rise_time, ([gy_value * (NPE1 // 2 - i) / NPE1 for i in range(NPE1)], ['PE1']), gy_rise_time)
        GZ(gx_rise_time * 2.5 - gz_rise_time, ([gz_value * (NPE2 // 2 - i) / NPE2 for i in range(NPE2)], ['PE2']), gz_rise_time)
        GX(NR // 2 * dwell_time + gx_rise_time * 2.5 - gx_rise_time, 0.0, gx_rise_time)
        GY(NR // 2 * dwell_time + gx_rise_time * 2.5 - gy_rise_time, 0.0, gy_rise_time)
        GZ(NR // 2 * dwell_time + gx_rise_time * 2.5 - gz_rise_time, 0.0, gz_rise_time)
        
    with Main():
        with Loop('PE2', NPE2):
            with Loop('PE1', NPE1):
                BlockRef('Excitation')
                WaitUntil(TE - NR // 2 * 2 * dwell_time - gx_rise_time * 2.5 + excitation_pulse_width * 0.5)
                BlockRef('PhaseEncoding')
                BlockRef('Readout')
                BlockRef('Rewinding')
                WaitUntil(TR)