3D RF spoiled gradient echo sequence. TR = 6 ms, TE = 2.2 ms, FA = 30°. The image matrix is 256×256×256. The number of subvoxel is 32×1×1. The total number of subvoxel is 108,793,856. The calculation time was 6145.1 s for 256×256×256 voxels.

 

Subvoxel number dependence for 256×256×32 voxel image acquisition:

32 subvoxels were required to suppress the ghosting artifact. This was because phase accumulation contained the magnetization in the voxel in the readout direction (x)  is proportional to 2π times the number of the repetition (TR). The amplitude variation caused by the phase modulation decreases with T2, which determines the required number of subvoxels (32 subvoxels corresponds to 192 ms time interval (6 ms × 32)).

 

FA dependence of the central slice for 256×256×32 voxel images (subvoxel = 32×1×1):

 

Pulse sequence chart visualized by the SequenceViewer:

One data acquisition sequence. TR = 6 ms. TE = 2.2 ms. Dwell_time = 6 μs. Maximum amplitudes for Gy and Gz are 18 mT/m. The maximum slew rate for Gy and Gz is 150 mT/m/ms (= 150 T/m/s).

 

The inner loop sequence. RF phase modulation for the RF spoiling is visualized.

 

Entire sequence for the 256×256×32 voxel acquisition.

 

Python sequence code:

from psdk import *
import numpy as np

gamma = 42.57747892 # [MHz/T]
TR = 6.0e+3 # [us]
TE = 2.2e+3 # [us]
NR = 256 # Number of readout points
NPE1 = 256 # Number of 1st phase encoding
NPE2 = 32 # Number of 2nd phase encoding
fov = [220.0, 220.0, 256.0] # [mm]
dwell_time = 6.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 = 120.0 # [us]
gy_rise_time = 120.0 # [us]
gz_rise_time = 120.0 # [us]
PW = 300.0 # [us]
excitation_pulse_flip_angle = 30.0 # [degree]
excitation_pulse_value = excitation_pulse_flip_angle / (360.0e-6 * gamma * PW) # [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', PW):
        RF(0.0, [excitation_pulse_value] * 10, PW / 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 + gy_rise_time):
        GY(0.0, ([gy_value * (NPE1 // 2 - i) / NPE1 for i in range(NPE1)], ['PE1']), gy_rise_time)
        GZ(0.0, ([gz_value * (NPE2 // 2 - i) / NPE2 for i in range(NPE2)], ['PE2']), gz_rise_time)
        GX(NR // 2 * dwell_time - gx_rise_time * 0.5, 0.0, gx_rise_time)
        GY(NR // 2 * dwell_time, 0.0, gy_rise_time)
        GZ(NR // 2 * dwell_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)
                BlockRef('PhaseEncoding')
                BlockRef('Readout')
                BlockRef('Rewinding')
                WaitUntil(TR)