2D multislice T2 weighted fast spin echo. TR = 4800 ms, TE = 96 ms, slice thickness = 5 mm. The image matrix is 256×256. The number of subvoxels was 1×1×16. The total number of subvoxels was 54,381,568. The calculation time was 896.7 s.

 

Image intensity variation:

 

Sequence visualized by the SequenceViewer:

 

Python sequence code:

from psdk import *
import numpy as np

gamma = 42.57747892 # [MHz/T]
TR = 4800.0e+3 # [us]
TE = 12.0e+3 # [us]
NR = 256 # Number of readout points
NPE1 = 32 # Number of 1st phase encoding
NEcho = 8 # Number of echoes
fov = [220.0, 220.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 = 2e+6 / (dwell_time * gamma * fov[1]) * NPE1 * 8 / NR # [mT/m]
gz_value = 1.25 / (slice_width * 1.0e-3) / gamma # [mT/m]
gx_rt = 300.0 # [us]
gy_rt = 300.0 # [us]
gz_rt = 300.0 # [us]
PW = 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)

def phase_correction(i):
    return ((i - 6) * (0.75) * np.pi + (i // 12) * (0.375) * np.pi)

def encode_number(i, j):
#    return i - ((NPE1) // 2) - j * ((NPE1) // 2) * (-1) ** (i // ((NPE1) // 2))
    return i - ((NPE1) // 2) - (7 - j) * ((NPE1) // 2) * (-1) ** (i // ((NPE1) // 2))

with Sequence('2D multislice multiple SpinEcho'):

    with Block('Excitation', PW + 2.0 * gz_rt):
        GZ(0.0, gz_value, gz_rt)
        RF(gz_rt, sinc_with_hamming(excitation_pulse_flip_angle, PW, 160), PW / 160,\
            phase=([phase_correction(i) for i in range(24)], ['SL']),\
            frequency=([-15.0, -12.5, -10.0, -7.5, -5.0, -2.5, 0.0, 2.5, 5.0, 7.5, 10.0, 12.5, -13.75, -11.25, -8.75, -6.25, -3.75, -1.25, 1.25, 3.75, 6.25, 8.75, 11.25, 13.75], ['SL']))
        GZ(PW + gz_rt, 0.0, gz_rt)

    with Block('Slice_refocus+Prephasing', PW * 0.5 + gz_rt * 2.0 ):    
         GX(0.0, gx_value * 2.0, gx_rt)
         GX((NR // 2) * dwell_time, 0.0, gx_rt)   
         GZ(0.0, -gz_value, gz_rt)
         GZ(PW * 0.5 + gz_rt, 0.0, gz_rt) 

    with Block('Refocus', PW + 2.0 * gz_rt):
        GZ(0.0, gz_value, gz_rt)
        RF(gz_rt, 2.0 * sinc_with_hamming(excitation_pulse_flip_angle, PW, 160), PW / 160,\
            phase=np.pi * 0.5,\
            frequency=([-15.0, -12.5, -10.0, -7.5, -5.0, -2.5, 0.0, 2.5, 5.0, 7.5, 10.0, 12.5, -13.75, -11.25, -8.75, -6.25, -3.75, -1.25, 1.25, 3.75, 6.25, 8.75, 11.25, 13.75], ['SL']))
        GZ(PW + gz_rt, 0.0, gz_rt)        

    with Block('Phase_encoding+Acquisition', TE/2 + 1.5 * NR * dwell_time - PW * 0.5 - gz_rt - gx_rt + gy_rt):
        GY(0.0, ([gy_value * encode_number(i, j) / (NPE1 * NEcho) for i in range(NPE1) for j in range(NEcho)], ['PE1', 'Echo']), gy_rt)
        GY(NR // 2 * dwell_time, 0.0, gy_rt)
        GX(TE/2 - PW * 0.5 - gz_rt - NR * dwell_time - 0.5 * gx_rt, gx_value, gx_rt)
        AD(TE/2 - PW * 0.5 - gz_rt - NR // 2 * dwell_time, NR, dwell_time)
        GX(TE/2 - PW * 0.5 - gz_rt + NR * dwell_time - gx_rt * 0.5, 0, gx_rt)
        GY(TE/2 - PW * 0.5 - gz_rt + NR * dwell_time - gx_rt * 0.5, ([-gy_value * encode_number(i, j) / (NPE1 * NEcho) for i in range(NPE1) for j in range(NEcho)], ['PE1', 'Echo']), gy_rt)
        GY(TE/2 - PW * 0.5 - gz_rt + NR * dwell_time - gx_rt * 0.5 + NR // 2 * dwell_time, 0.0, gy_rt)

    with Main():
        with Loop('PE1', NPE1):
            with Loop('SL', 24):
                BlockRef('Excitation')
                BlockRef("Slice_refocus+Prephasing")
                WaitUntil(TE/2)
                with Loop('Echo', NEcho):
                    BlockRef('Refocus')
                    BlockRef('Phase_encoding+Acquisition')
                    WaitUntil(TE * 1.0)    
                WaitUntil(200.0e+3)                
            WaitUntil(TR)