2D multislice Fast Spin echo imaging

T2 weighted FSE: TR = 4800 ms, TE = 96 ms (effective)

Image matrix = 256 x 256, slice thickness = 5 mm, TR = 4800 ms, TE = 96 ms, Echo spacing = 12 ms, subvoxel = 1 x 1 x 16, number of isochromats = 53,381,568, calculation time = 4,877,580 s (81.3 min)

PD weighted FSE: TR = 4800 ms, TE = 12 ms (effective)

Image matrix = 256 x 256, slice thickness = 5 mm, TR = 4800 ms, TE = 12 ms, Echo spacing = 12 ms, subvoxel = 1 x 1 x 16, number of isochromats = 53,381,568, calculation time = 4,882,239 s (81.4 min)

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 = [256.0, 256.0, 256.0] # [mm]
dwell_time = 5.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 * NEcho / NR # [mT/m]
gz_value = 5.00 / (slice_width * 1.0e-3) / gamma # [mT/m]
gx_rise_time = 300.0 # [us]
gy_rise_time = 300.0 # [us]
gz_rise_time = 300.0 # [us]
excitation_pulse_width = 1600.0 # [us]
excitation_pulse_flip_angle = 90.0 # [degree]

def sinc(flip_angle, pulse_width, points, *, min=-4.0*np.pi, max=4.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))
    return flip_angle * y * points / (y.sum() * pulse_width * 360.0e-6 * gamma)

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

def phase_correction(i):
    return (i//12)*(0.5) * np.pi

with Sequence('2D multislice FastSpinEcho'):
    with Block('Excitation', excitation_pulse_width + 2.0*gz_rise_time):   
        GZ(0.0, gz_value, gz_rise_time)
        RF(gz_rise_time, sinc(excitation_pulse_flip_angle, excitation_pulse_width, 320), excitation_pulse_width / 320,
            phase=([phase_correction(i) for i in range(24)], ['SL']),
            frequency=([-60.0, -50.0, -40.0, -30.0, -20.0, -10.0, 0.0, 10.0, 20.0, 30.0, 40.0, 50.0, -55.0, -45.0, -35.0, -25.0, -15.0, -5.0, 5.0, 15.0, 25.0, 35.0, 45.0, 55.0], ['SL']))
        GZ(excitation_pulse_width + gz_rise_time, 0.0, gz_rise_time)
    with Block('Slice_refocus', excitation_pulse_width * 0.5 + gz_rise_time * 2.0) :
        GZ(0.0, -gz_value, gz_rise_time)
        GZ(excitation_pulse_width * 0.5 + gz_rise_time, 0.0, gz_rise_time) 
    with Block('ReadPredephasing', NR // 2 * dwell_time * 2.0 + gx_rise_time ):
        GX(0.0, gx_value, gx_rise_time)
        GX(NR // 2 * dwell_time * 2.0, 0.0, gx_rise_time)         
    with Block('Refocus', excitation_pulse_width + 2.0*gz_rise_time):         
        GZ(0.0, gz_value, gz_rise_time)
        RF(gz_rise_time, 2.0 * sinc(excitation_pulse_flip_angle, excitation_pulse_width, 320), excitation_pulse_width / 320,
            phase=np.pi * 0.5,
            frequency=([-60.0, -50.0, -40.0, -30.0, -20.0, -10.0, 0.0, 10.0, 20.0, 30.0, 40.0, 50.0, -55.0, -45.0, -35.0, -25.0, -15.0, -5.0, 5.0, 15.0, 25.0, 35.0, 45.0, 55.0], ['SL']))
        GZ(excitation_pulse_width + gz_rise_time, 0.0, gz_rise_time)        
    with Block('PhaseEncoding', NR // 2 * dwell_time + gx_rise_time * 2.5):
        GY(0.0, ([gy_value * encode_number(i, j) / (NPE1 * NEcho) for i in range(NPE1) for j in range(NEcho)], ['PE1', 'Echo']), gy_rise_time)
        GY(NR // 2 * dwell_time, 0.0, gy_rise_time)
    with Block('Readout', NR * dwell_time * 2.0 + gx_rise_time):
        GX(0.0, gx_value, gx_rise_time)
        AD(NR // 2 * dwell_time + gx_rise_time * 0.5, NR, dwell_time)
        GX(NR * dwell_time * 2.0, 0, gx_rise_time)
    with Block('Rewinding', NR // 2 * dwell_time + gx_rise_time * 2.5):
        GY(0.0, ([-gy_value * encode_number(i, j) / (NPE1 * NEcho) for i in range(NPE1) for j in range(NEcho)], ['PE1', 'Echo']), gy_rise_time)
        GY(NR // 2 * dwell_time, 0.0, gy_rise_time)
        
    with Main():
        with Loop('PE1', NPE1):
            with Loop('SL', 24):
                BlockRef('Excitation')
                BlockRef('Slice_refocus')
                BlockRef('ReadPredephasing')
                WaitUntil(TE/2)
                with Loop('Echo', NEcho):
                    BlockRef('Refocus')
                    BlockRef('PhaseEncoding')
                    WaitUntil(TE * 0.5 + gz_rise_time + excitation_pulse_width * 0.5 - NR * dwell_time - gx_rise_time) 
                    BlockRef('Readout')
                    BlockRef('Rewinding')
                    WaitUntil(TE * 1.0)    
                WaitUntil(200.0e+3)                
            WaitUntil(TR)