Magnetization prepared rapid gradient echo (MPRAGE). TR = 4000 ms, echo spacing = 6 ms, TE = 2.2 ms, TI = 1100 ms. The number of subvoxel is 32×1×1. The total number of subvoxels is 109,297,664. The calculation time for 256×256×256 voxel image was 6141 s using RTX 2080Ti.

 

Subvoxel dependence of the central cross section for 256×256×32 voxel images:

 

Inversion time dependence of the MPRAGE image:

 

Inversion time dependence of the MPRAGE image:

 

MPRAGE sequence chart visualized by the SequenceViewer:

One data acquisition sequence. Echo spacing = 6 ms, TE = 2.2 ms, dwell time = 6 μs.

 

Hyperbolic secant adiabatic inversion pulse. Duration time is 10 ms.

 

One inversion sequence. The delay time is 1100 ms. The total length of the data acquisition sequence is 1536 ms (6 ms x 256). RF spoiling is used.

 

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 = 256 # 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]
excitation_pulse_width = 100.0 # [us]
excitation_pulse_flip_angle = 4.0 # [degree]
excitation_pulse_value = excitation_pulse_flip_angle / (360.0e-6 * gamma * excitation_pulse_width) # [uT]
PW = 10000.0

def hyperbolic_secant_pulse(FA, pulse_width, points, window, mu):
    x = np.linspace(-window, window, points)
    y = np.power((1.0 / np.cosh(x)), complex(1, mu))
    return FA * y * points / (y.sum() * pulse_width * 360.0e-6 * gamma)

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

with Sequence('MPRAGE'):

    with Block('Inversion', 100000):
        GZ(0.0, 0.0, gz_rise_time)
        RF(gz_rise_time, hyperbolic_secant_pulse(270, 10000, 500, 5, 5), PW / (500))
        GZ(PW + gz_rise_time, 0.0, gz_rise_time)

    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 + 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):
            WaitFor(2132000)
            BlockRef('Inversion')
            WaitFor(232000)
            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)