2D radial scan. The step of the radial trajectory is π/256. TR = 12 ms, TE = 6 ms, FA = 60°. The number of subvoxel is 1×1×4. The total number of subvoxels is 135,72,096. The calculation time was 13 s.

 

Pulse sequence chart visualized by the SequenceViewer:

One data acquisition sequence.

 

Entire sequence.

Python sequence code:

from psdk import *
import numpy as np

gamma = 42.57747892 # [MHz/T]
TR = 12.0e+3 # [us]
TE = 6.0e+3 # [us]
NR = 256 # Number of readout points
NPE1 = 256 # Number of 1st phase encoding
fov = [256.0, 256.0, 256.0] # [mm]
dwell_time = 10.0 # [us]
slice_width = 5.0 # [mm]
gx_value = 1e+6 / (dwell_time * gamma * fov[0]) # [mT/m]
gy_value = 1e+6 / (dwell_time * gamma * fov[1]) # [mT/m]
gz_value = 2.5 / (slice_width * 1.0e-3) / gamma # [mT/m]
gx_rt = 300.0 # [us] gx rise time
gy_rt = 300.0 # [us] gy rise time
gz_rt = 300.0 # [us] gz rise time
PW = 1600.0 # [us]
ex_pulse_flip_angle = 60.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)

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

with Sequence('2D GRE projection'):

    with Block('Excitation', PW + 4.0 * gz_rt):
        GZ(0.0, gz_value, gz_rt)
        RF(gz_rt, sinc_with_hamming(ex_pulse_flip_angle, PW, 160), PW / 160, phase=([phase_shift_angle(i) for i in range(NPE1)], ['PE1']))
        GZ(PW + gz_rt, -gz_value, 2.0 * gz_rt)
        GZ(PW * 1.5 + 3.0 * gz_rt, 0.0, gz_rt)
       
    with Block('PhaseEncoding', NR // 2 * dwell_time + gx_rt * 2.5):
        GX(0.0,  ([-gx_value * np.cos(np.pi * i / NPE1) for i in range(NPE1)], ['PE1']), gx_rt)
        GY(0.0,  ([-gy_value * np.sin(np.pi * i / NPE1) for i in range(NPE1)], ['PE1']), gy_rt)
        GX(NR // 2 * dwell_time + gx_rt * 0.5, ([gx_value * np.cos(np.pi * i / NPE1) for i in range(NPE1)], ['PE1']), gx_rt * 2.0)
        GY(NR // 2 * dwell_time + gy_rt * 0.5, ([gy_value * np.sin(np.pi * i / NPE1) for i in range(NPE1)], ['PE1']), gy_rt * 2.0)        

    with Block('Readout', NR * dwell_time):
        AD(0.0, NR, dwell_time, phase=([phase_shift_angle(i) for i in range(NPE1)], ['PE1']))
        
    with Block('Rewinding', NR // 2 * dwell_time + gx_rt * 2.5):
        GX(0.0,  ([-gx_value * np.cos(np.pi * i / NPE1) for i in range(NPE1)], ['PE1']), gx_rt * 2.0)
        GY(0.0,  ([-gy_value * np.sin(np.pi * i / NPE1) for i in range(NPE1)], ['PE1']), gy_rt * 2.0)
        GX(NR // 2 * dwell_time + gx_rt * 1.5, 0.0, gx_rt)
        GY(NR // 2 * dwell_time + gy_rt * 1.5, 0.0, gy_rt)        
        
    with Main():
        with Loop('PE1', 10):
            BlockRef('Excitation')
            WaitUntil(TE + PW * 0.5 + gz_rt - NR // 2 * 2 * dwell_time - gx_rt * 2.5)
            BlockRef('PhaseEncoding')
            WaitFor(NR * dwell_time)
            BlockRef('Rewinding')
            WaitUntil(TR)
        with Loop('PE1', NPE1):
            BlockRef('Excitation')
            WaitUntil(TE + PW * 0.5 + gz_rt - NR // 2 * 2 * dwell_time - gx_rt * 2.5)
            BlockRef('PhaseEncoding')
            BlockRef('Readout')
            BlockRef('Rewinding')
            WaitUntil(TR)