physics_based_augmentation/seqgen/MRI_sequences/py2jemris/pulseq_library.py

# Copyright of the Board of Trustees of Columbia University in the City of New York
"""
Library for generating pulseq sequences: GRE, SE, IRSE, EPI
"""

# TODO update for PyPulseq 1.2.1

import copy
from math import pi, sqrt, ceil, floor

import numpy as np
from pypulseq.Sequence.sequence import Sequence
from pypulseq.calc_duration import calc_duration
from pypulseq.make_adc import make_adc
from pypulseq.make_delay import make_delay
from pypulseq.make_sinc_pulse import make_sinc_pulse
from pypulseq.make_trap_pulse import make_trapezoid
from pypulseq.opts import Opts

GAMMA_BAR = 42.5775e6
GAMMA = 2*pi*GAMMA_BAR


def make_pulseq_gre(fov,n,thk,fa,tr,te,enc='xyz',slice_locs=None,write=False):
    """Makes a gradient-echo sequence

    2D orthogonal multi-slice gradient-echo pulse sequence with Cartesian encoding
    Orthogonal means that each of slice-selection, phase encoding, and frequency encoding
    aligns with the x, y, or z directions

    Parameters
    ----------
    fov : float
        Field-of-view in meters (isotropic)
    n : int
        Matrix size (isotropic)
    thk : float
        Slice thickness in meters
    fa : float
        Flip angle in degrees
    tr : float
        Repetition time in seconds
    te : float
        Echo time in seconds
    enc : str, optional
        Spatial encoding directions
        1st - readout; 2nd - phase encoding; 3rd - slice select
        Default 'xyz' means axial(z) slice with readout in x and phase encoding in y
    slice_locs : array_like, optional
        Slice locations from isocenter in meters
        Default is None which means a single slice at the center
    write : bool, optional
        Whether to write seq into file; default is False

    Returns
    -------
    seq : Sequence
        Pulse sequence as a Pulseq object

    """
    system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130,
                  slew_unit='T/m/s', rf_ringdown_time=30e-6,
                  rf_dead_time=100e-6, adc_dead_time=20e-6)
    seq = Sequence(system)


    Nf = n
    Np = n
    flip = fa * pi / 180
    rf, g_ss, __ = make_sinc_pulse(flip_angle=flip, system=system, duration=4e-3, slice_thickness=thk,
                               apodization=0.5, time_bw_product=4)

    g_ss.channel = enc[2]

    delta_k = 1 / fov
    kWidth = Nf * delta_k

    # Readout and ADC
    readoutTime = 6.4e-3
    g_ro= make_trapezoid(channel=enc[0],system=system,flat_area=kWidth,flat_time=readoutTime)
    adc = make_adc(num_samples=Nf, duration=g_ro.flat_time, delay=g_ro.rise_time)

    # Readout rewinder
    g_ro_pre = make_trapezoid(channel=enc[0],system=system,area=-g_ro.area/2,duration=2e-3)
    # Slice refocusing
    g_ss_reph = make_trapezoid(channel=enc[2],system=system,area=-g_ss.area/2,duration=2e-3)
    phase_areas = (np.arange(Np) - (Np / 2)) * delta_k

    # TE, TR = 10e-3, 1000e-3
    TE, TR = te,tr
    delayTE = TE - calc_duration(g_ro_pre) - calc_duration(g_ss) / 2 - calc_duration(g_ro) / 2
    delayTR = TR - calc_duration(g_ro_pre) - calc_duration(g_ss) - calc_duration(g_ro) - delayTE
    delay1 = make_delay(delayTE)
    delay2 = make_delay(delayTR)

    if slice_locs is None:
        locs = [0]
    else:
        locs = slice_locs

    for u in range(len(locs)):
        # add frequency offset
        rf.freq_offset = g_ss.amplitude * locs[u]
        for i in range(Np):
            seq.add_block(rf, g_ss)
            g_pe = make_trapezoid(channel=enc[1],system=system,area=phase_areas[i],duration=2e-3)
            seq.add_block(g_ro_pre, g_pe, g_ss_reph)
            seq.add_block(delay1)
            seq.add_block(g_ro,adc)
            seq.add_block(delay2)

    if write:
        seq.write("gre_fov{:.0f}mm_Nf{:d}_Np{:d}_TE{:.0f}ms_TR{:.0f}ms_FA{:.0f}deg.seq".format(fov * 1000, Nf, Np, TE * 1000,
                                                                                               TR * 1000, flip * 180 / pi))
    print('GRE sequence constructed')
    return seq


def make_pulseq_gre_oblique(fov,n,thk,fa,tr,te,enc='xyz',slice_locs=None,write=False):
    """Makes a gradient-echo sequence in any plane

        2D oblique multi-slice gradient-echo pulse sequence with Cartesian encoding
        Oblique means that each of slice-selection, phase encoding, and frequency encoding
        can point in any specified direction

        Parameters
        ----------
        fov : array_like
            Isotropic field-of-view, or length-2 list [fov_readout, fov_phase], in meters
        n : array_like
            Isotropic matrix size, or length-2 list [n_readout, n_phase]
        thk : float
            Slice thickness in meters
        fa : float
            Flip angle in degrees
        tr : float
            Repetition time in seconds
        te : float
            Echo time in seconds
        enc : str or array_like, optional
            Spatial encoding directions
            1st - readout; 2nd - phase encoding; 3rd - slice select
            - Use str with any permutation of x, y, and z to obtain orthogonal slices
            e.g. The default 'xyz' means axial(z) slice with readout in x and phase encoding in y
            - Use list to indicate vectors in the encoding directions for oblique slices
            They should be perpendicular to each other, but not necessarily unit vectors
            e.g. [(2,1,0),(-1,2,0),(0,0,1)] rotates the two in-plane encoding directions for an axial slice
        slice_locs : array_like, optional
            Slice locations from isocenter in meters
            Default is None which means a single slice at the center
        write : bool, optional
            Whether to write seq into file; default is False

        Returns
        -------
        seq : Sequence
            Pulse sequence as a Pulseq object

        """

    # System options
    system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130,
                  slew_unit='T/m/s', rf_ringdown_time=30e-6,
                  rf_dead_time=100e-6, adc_dead_time=20e-6)
    seq = Sequence(system)

    # Calculate unit gradients for ss, fe, pe
    ug_fe, ug_pe, ug_ss = parse_enc(enc)

    # Sequence parameters
    Nf, Np = (n,n) if isinstance(n,int) else (n[0], n[1])
    delta_k_ro, delta_k_pe = (1/fov,1/fov) if isinstance(fov,float) else (1/fov[0], 1/fov[1])
    kWidth_ro = Nf * delta_k_ro
    flip = fa * pi / 180

    # Slice select: RF and gradient
    rf, g_ss, __ = make_sinc_pulse(flip_angle=flip,system=system,duration=4e-3,slice_thickness=thk,
                               apodization=0.5, time_bw_product=4)
    g_ss_x, g_ss_y, g_ss_z = make_oblique_gradients(g_ss,ug_ss)

    # Readout and ADC
#    readoutTime = 6.4e-3
    dwell = 10e-6
    g_ro= make_trapezoid(channel='x',system=system,flat_area=kWidth_ro, flat_time=dwell*Nf)
    g_ro_x, g_ro_y, g_ro_z = make_oblique_gradients(g_ro,ug_fe)#
    adc = make_adc(num_samples=Nf, duration=g_ro.flat_time,delay=g_ro.rise_time)

    # Readout rewinder
    g_ro_pre = make_trapezoid(channel='x',system=system,area=-g_ro.area/2,duration=2e-3)
    g_ro_pre_x, g_ro_pre_y, g_ro_pre_z = make_oblique_gradients(g_ro_pre,ug_fe)#

    # Slice refocusing
    g_ss_reph = make_trapezoid(channel='z',system=system,area=-g_ss.area/2,duration=2e-3)
    g_ss_reph_x, g_ss_reph_y, g_ss_reph_z = make_oblique_gradients(g_ss_reph, ug_ss)

    # Prepare phase areas
    phase_areas = (np.arange(Np) - (Np / 2)) * delta_k_pe

    TE, TR = te,tr
    delayTE = TE - calc_duration(g_ro_pre) - calc_duration(g_ss) / 2 - calc_duration(g_ro) / 2
    delayTR = TR - calc_duration(g_ro_pre) - calc_duration(g_ss) - calc_duration(g_ro) - delayTE
    delay1 = make_delay(delayTE)
    delay2 = make_delay(delayTR)

    if slice_locs is None:
        locs = [0]
    else:
        locs = slice_locs

    # Construct sequence!
    for u in range(len(locs)):
        # add frequency offset
        rf.freq_offset = g_ss.amplitude * locs[u]
        for i in range(Np):
            seq.add_block(rf,g_ss_x, g_ss_y, g_ss_z)
            g_pe = make_trapezoid(channel='y',system=system,area=phase_areas[i],duration=2e-3)

            g_pe_x, g_pe_y, g_pe_z = make_oblique_gradients(g_pe,ug_pe)

            pre_grads_list = [g_ro_pre_x, g_ro_pre_y, g_ro_pre_z,
                             g_ss_reph_x, g_ss_reph_y, g_ss_reph_z,
                             g_pe_x, g_pe_y, g_pe_z]

            gtx, gty, gtz = combine_trap_grad_xyz(gradients=pre_grads_list,system=system, dur=2e-3)

            seq.add_block(gtx, gty, gtz)
            seq.add_block(delay1)
            seq.add_block(g_ro_x, g_ro_y, g_ro_z, adc)
            seq.add_block(delay2)

    if write:
        seq.write("gre_fov{:.0f}mm_Nf{:d}_Np{:d}_TE{:.0f}ms_TR{:.0f}ms_FA{:.0f}deg.seq".format(fov * 1000, Nf, Np, TE * 1000,
                                                                                               TR * 1000, flip * 180 / pi))
    print('GRE sequence constructed')
    return seq


def make_pulseq_irse(fov,n,thk,fa,tr,te,ti,enc='xyz',slice_locs=None,write=False):
    """Makes an Inversion Recovery Spin Echo (IRSE) sequence

        2D orthogonal multi-slice IRSE pulse sequence with Cartesian encoding
        Orthogonal means that each of slice-selection, phase encoding, and frequency encoding
        aligns with the x, y, or z directions

        Parameters
        ----------
        fov : float
            Field-of-view in meters (isotropic)
        n : int
            Matrix size (isotropic)
        thk : float
            Slice thickness in meters
        fa : float
            Flip angle in degrees
        tr : float
            Repetition time in seconds
        te : float
            Echo time in seconds
        ti : float
            Inversion time in seconds
        enc : str, optional
            Spatial encoding directions
            1st - readout; 2nd - phase encoding; 3rd - slice select
            Default 'xyz' means axial(z) slice with readout in x and phase encoding in y
        slice_locs : array_like, optional
            Slice locations from isocenter in meters
            Default is None which means a single slice at the center
        write : bool, optional
            Whether to write seq into file; default is False

        Returns
        -------
        seq : Sequence
            Pulse sequence as a Pulseq object

        """
    system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130,
                  slew_unit='T/m/s', rf_ringdown_time=30e-6,
                  rf_dead_time=100e-6, adc_dead_time=20e-6)
    seq = Sequence(system)

    # Parameters
    Nf = n
    Np = n
    delta_k = 1 / fov
    kWidth = Nf * delta_k

    TI,TE,TR = ti,te,tr

    if np.shape(TI) == ():
        TI = [TI]


    # Non-180 pulse
    flip1 = fa * pi / 180
    rf, g_ss, __ = make_sinc_pulse(flip_angle=flip1, system=system, duration=2e-3, slice_thickness=thk,
                               apodization=0.5, time_bw_product=4)

    g_ss.channel = enc[2]

    # 180 pulse
    flip2 = 180 * pi / 180
    rf180, g_ss180, __ = make_sinc_pulse(flip_angle=flip2, system=system,duration=2e-3, slice_thickness=thk,
                                     apodization=0.5, time_bw_product=4)
    g_ss180.channel = enc[2]

    # Readout gradient & ADC
    readoutTime = 6.4e-3

    g_ro = make_trapezoid(channel=enc[0],system=system, flat_area=kWidth, flat_time=readoutTime)
    adc = make_adc(num_samples=Nf, system=system, duration=g_ro.flat_time, delay=g_ro.rise_time)

    # RO rewinder gradient
    g_ro_pre = make_trapezoid(channel=enc[0],system=system,area=g_ro.area/2,duration=2e-3)

    # Slice refocusing gradient
    g_ss_reph = make_trapezoid(channel=enc[2],system=system,area=-g_ss.area/2,duration=2e-3)

    # Delays

    delayTE1 = TE / 2 - max(calc_duration(g_ss_reph), calc_duration(g_ro_pre)) - calc_duration(g_ss) / 2 - calc_duration(
        g_ss180) / 2
    delayTE2 = TE / 2 - calc_duration(g_ro) / 2 - calc_duration(g_ss180) / 2
    delayTE3 = TR - TE - calc_duration(g_ss) / 2 - calc_duration(g_ro) / 2

    print('dur rf', calc_duration(rf),'dur gss:' ,calc_duration(g_ss))

    delay1 = make_delay(delayTE1)
    delay2 = make_delay(delayTE2)
    delay3 = make_delay(delayTE3)

    # Construct sequence
    if slice_locs is None:
        locs = [0]
    else:
        locs = slice_locs

    for inv in range(len(TI)):
        for u in range(len(locs)):
            rf180.freq_offset = g_ss180.amplitude * locs[u]
            rf.freq_offset = g_ss.amplitude * locs[u]
            for i in range(Np):
                # Inversion Recovery part
                seq.add_block(rf180, g_ss180)# Selective; potentially extended to be non-selective or adiabatic
                seq.add_block(make_delay(TI[inv] - calc_duration(rf) / 2 - calc_duration(rf180) / 2))  # Inversion time delay
                # Spin echo part
                seq.add_block(rf, g_ss)  # 90-deg pulse
                g_pe_pre = make_trapezoid(channel=enc[1],system=system,area=-(Np/2-i)*delta_k,
                                          duration=2e-3)  # Phase encoding gradient
                seq.add_block(g_ro_pre, g_pe_pre, g_ss_reph)  # Add a combination of ro rewinder, phase encoding, and slice refocusing
                seq.add_block(delay1)  # Delay 1: until 180-deg pulse
                seq.add_block(rf180, g_ss180)  # 180 deg pulse for SE
                seq.add_block(delay2)  # Delay 2: until readout
                seq.add_block(g_ro, adc)  # Readout!
                seq.add_block(delay3)  # Delay 3: until next inversion pulse

    if write:
        if len(TI) == 1:
            seq.write("irse_fov{:.0f}mm_Nf{:d}_Np{:d}_TI{:.0f}ms_TE{:.0f}ms_TR{:.0f}ms.seq".format(fov * 1000, Nf, Np, TI[0] * 1000, TE * 1000, TR * 1000))
        else:
            seq.write("irse_fov{:.0f}mm_Nf{:d}_Np{:d}_multiTI_TE{:.0f}ms_TR{:.0f}ms.seq".format(fov * 1000, Nf, Np, TE * 1000, TR * 1000))

    print('IRSE sequence constructed')
    return seq

def make_pulseq_irse_oblique(fov,n,thk,fa,tr,te,ti,enc='xyz',slice_locs=None,write=False):
    """Makes an Inversion Recovery Spin Echo (IRSE) sequence in any plane

        2D oblique multi-slice IRSE pulse sequence with Cartesian encoding
        Oblique means that each of slice-selection, phase encoding, and frequency encoding
        can point in any specified direction

        Parameters
        ----------
        fov : array_like
            Isotropic field-of-view, or length-2 list [fov_readout, fov_phase], in meters
        n : array_like
            Isotropic matrix size, or length-2 list [n_readout, n_phase]
        thk : float
            Slice thickness in meters
        fa : float
            Flip angle in degrees
        tr : float
            Repetition time in seconds
        te : float
            Echo time in seconds
        ti : float
            Inversion time in seconds
        enc : str or array_like, optional
            Spatial encoding directions
            1st - readout; 2nd - phase encoding; 3rd - slice select
            - Use str with any permutation of x, y, and z to obtain orthogonal slices
            e.g. The default 'xyz' means axial(z) slice with readout in x and phase encoding in y
            - Use list to indicate vectors in the encoding directions for oblique slices
            They should be perpendicular to each other, but not necessarily unit vectors
            e.g. [(2,1,0),(-1,2,0),(0,0,1)] rotates the two in-plane encoding directions for an axial slice
        slice_locs : array_like, optional
            Slice locations from isocenter in meters
            Default is None which means a single slice at the center
        write : bool, optional
            Whether to write seq into file; default is False


        Returns
        -------
        seq : Sequence
            Pulse sequence as a Pulseq object

        """
    # System options
    system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130,
                  slew_unit='T/m/s', rf_ringdown_time=30e-6,
                  rf_dead_time=100e-6, adc_dead_time=20e-6)
    seq = Sequence(system)

    # Sequence parameters
    ug_fe, ug_pe, ug_ss = parse_enc(enc)
    Nf, Np = (n,n) if isinstance(n,int) else (n[0], n[1])
    delta_k_ro, delta_k_pe = (1/fov,1/fov) if isinstance(fov,float) else (1/fov[0], 1/fov[1])
    kWidth_ro = Nf * delta_k_ro
    TI,TE,TR = ti,te,tr

    if np.shape(TI) == ():
        TI = [TI]

    # Non-180 pulse
    flip1 = fa * pi / 180
    rf, g_ss, __ = make_sinc_pulse(flip_angle=flip1, system=system, duration=2e-3, slice_thickness=thk,
                                apodization=0.5, time_bw_product=4)
    g_ss_x, g_ss_y, g_ss_z = make_oblique_gradients(g_ss, ug_ss)

    # 180 pulse
    flip2 = 180 * pi / 180
    rf180, g_ss180, __ = make_sinc_pulse(flip_angle=flip2, system=system, duration=2e-3, slice_thickness=thk,
                                    apodization=0.5, time_bw_product=4)
    g_ss180_x, g_ss180_y, g_ss180_z = make_oblique_gradients(g_ss180, ug_ss)

    # Readout gradient & ADC
    readoutTime = 6.4e-3

    g_ro = make_trapezoid(channel='x', system=system, flat_area=kWidth_ro, flat_time= readoutTime)
    g_ro_x, g_ro_y, g_ro_z = make_oblique_gradients(g_ro, ug_fe)

    adc = make_adc(num_samples=Nf, system=system, duration=g_ro.flat_time, delay=g_ro.rise_time)

    # RO rewinder gradient
    g_ro_pre = make_trapezoid(channel=enc[0], system=system, area=g_ro.area/2,duration=2e-3)
    g_ro_pre_x, g_ro_pre_y, g_ro_pre_z = make_oblique_gradients(g_ro_pre,ug_fe)#

    # Slice refocusing gradient
    g_ss_reph = make_trapezoid(channel=enc[2],system=system,area=-g_ss.area/2,duration=2e-3)
    g_ss_reph_x, g_ss_reph_y, g_ss_reph_z = make_oblique_gradients(g_ss_reph, ug_ss)

    # Delays
    delayTE1 = TE / 2 - max(calc_duration(g_ss_reph), calc_duration(g_ro_pre))\
               - calc_duration(g_ss) / 2 - calc_duration(g_ss180) / 2
    delayTE2 = TE / 2 - calc_duration(g_ro) / 2 - calc_duration(g_ss180) / 2
    delayTE3 = TR - TE - calc_duration(g_ss) / 2 - calc_duration(g_ro) / 2

    delay1 = make_delay(delayTE1)
    delay2 = make_delay(delayTE2)
    delay3 = make_delay(delayTE3)

    # Construct sequence
    if slice_locs is None:
        locs = [0]
    else:
        locs = slice_locs

    for inv in range(len(TI)):
        for u in range(len(locs)):
            rf180.freq_offset = g_ss180.amplitude * locs[u]
            rf.freq_offset = g_ss.amplitude * locs[u]
            for i in range(Np):
                # Inversion Recovery part
                seq.add_block(rf180, g_ss180_x, g_ss180_y, g_ss180_z)# Non-selective at the moment; could be extended to make this selective/adiabatic
                seq.add_block(make_delay(TI[inv] - calc_duration(rf) / 2 - calc_duration(rf180) / 2))  # Inversion time delay
                # Spin echo part
                seq.add_block(rf, g_ss_x, g_ss_y, g_ss_z)  # 90-deg pulse
                g_pe = make_trapezoid(channel='y', system=system, area=-(Np /2 - i)*delta_k_pe, duration=2e-3)  # Phase encoding gradient
                g_pe_x, g_pe_y, g_pe_z = make_oblique_gradients(g_pe, ug_pe)

                pre_grads_list = [g_ro_pre_x, g_ro_pre_y, g_ro_pre_z,
                                  g_ss_reph_x, g_ss_reph_y, g_ss_reph_z,
                                  g_pe_x, g_pe_y, g_pe_z]
                gtx, gty, gtz = combine_trap_grad_xyz(pre_grads_list,system,2e-3)

                seq.add_block(gtx, gty, gtz)  # Add a combination of ro rewinder, phase encoding, and slice refocusing
                seq.add_block(delay1)  # Delay 1: until 180-deg pulse
                seq.add_block(rf180, g_ss180_x, g_ss180_y, g_ss180_z)  # 180 deg pulse for SE
                seq.add_block(delay2)  # Delay 2: until readout
                seq.add_block(g_ro_x, g_ro_y, g_ro_z, adc)  # Readout!
                seq.add_block(delay3)  # Delay 3: until next inversion pulse

    if write:
        if len(TI) == 1:
            seq.write("irse_fov{:.0f}mm_Nf{:d}_Np{:d}_TI{:.0f}ms_TE{:.0f}ms_TR{:.0f}ms_FA{:d}deg.seq".format(fov * 1000,
                                                                            Nf, Np, TI[0] * 1000, TE * 1000, TR * 1000, fa))
        else:
            seq.write("irse_fov{:.0f}mm_Nf{:d}_Np{:d}_multiTI_TE{:.0f}ms_TR{:.0f}ms_FA{:d}deg.seq".format(fov * 1000,
                                                                            Nf, Np, TE * 1000, TR * 1000, fa))

    print('IRSE (oblique) sequence constructed')
    return seq

def make_pulseq_se(fov,n,thk,fa,tr,te,enc='xyz',slice_locs=None,write=False):
    """Makes a Spin Echo (SE) sequence

    2D orthogonal multi-slice Spin-Echo pulse sequence with Cartesian encoding
    Orthogonal means that each of slice-selection, phase encoding, and frequency encoding
    aligns with the x, y, or z directions

    Parameters
    ----------
    fov : float
        Field-of-view in meters (isotropic)
    n : int
        Matrix size (isotropic)
    thk : float
        Slice thickness in meters
    fa : float
        Flip angle in degrees
    tr : float
        Repetition time in seconds
    te : float
        Echo time in seconds
    enc : str, optional
        Spatial encoding directions
        1st - readout; 2nd - phase encoding; 3rd - slice select
        Default 'xyz' means axial(z) slice with readout in x and phase encoding in y
    slice_locs : array_like, optional
        Slice locations from isocenter in meters
        Default is None which means a single slice at the center
    write : bool, optional
        Whether to write seq into file; default is False

    Returns
    -------
    seq : Sequence
        Pulse sequence as a Pulseq object

    """
    system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130,
                  slew_unit='T/m/s', rf_ringdown_time=30e-6,
                  rf_dead_time=100e-6, adc_dead_time=20e-6)
    seq = Sequence(system)

    # Parameters
    Nf = n
    Np = n
    delta_k = 1 / fov
    kWidth = Nf * delta_k

    TE,TR = te,tr


    # Non-180 pulse
    flip1 = fa * pi / 180
    rf, g_ss, __ = make_sinc_pulse(flip_angle=flip1, system=system, duration=2e-3, slice_thickness=thk,
                                apodization=0.5, time_bw_product=4)
    g_ss.channel = enc[2]


    # 180 pulse
    flip2 = 180 * pi / 180
    rf180, g_ss180, __ = make_sinc_pulse(flip_angle=flip2, system=system, duration=2e-3, slice_thickness=thk,
                                    apodization=0.5, time_bw_product=4)
    g_ss180.channel = enc[2]

    # Readout gradient & ADC
#    readoutTime = system.grad_raster_time * Nf
    readoutTime = 6.4e-3
    g_ro = make_trapezoid(channel=enc[0],system=system,flat_area=kWidth,flat_time=readoutTime)
    adc = make_adc(num_samples=Nf, system=system, duration=g_ro.flat_time, delay=g_ro.rise_time)

    # RO rewinder gradient
    g_ro_pre = make_trapezoid(channel=enc[0],system=system,area=g_ro.area/2,duration=2e-3)

    # Slice refocusing gradient
    g_ss_reph = make_trapezoid(channel=enc[2],system=system,area=-g_ss.area/2,duration=2e-3)

    # Delays
    delayTE1 = (TE - 2*max(calc_duration(g_ss_reph), calc_duration(g_ro_pre)) - calc_duration(g_ss) - calc_duration(
        g_ss180))/2
  # delayTE2 = TE / 2 - calc_duration(g_ro) / 2 - calc_duration(g_ss180) / 2
    delayTE2 = (TE - calc_duration(g_ro) - calc_duration(g_ss180))/2
    delayTE3 = TR - TE - (calc_duration(g_ss) + calc_duration(g_ro)) / 2


    delay1 = make_delay(delayTE1)
    delay2 = make_delay(delayTE2)
    delay3 = make_delay(delayTE3)


    # Construct sequence
    if slice_locs is None:
        locs = [0]
    else:
        locs = slice_locs

    for u in range(len(locs)):
        rf180.freq_offset = g_ss180.amplitude * locs[u]
        rf.freq_offset = g_ss.amplitude * locs[u]
        for i in range(Np):
            seq.add_block(rf, g_ss)  # 90-deg pulse
            g_pe_pre = make_trapezoid(channel=enc[1],system=system,area=-(Np/2-i)*delta_k,duration=2e-3)  # Phase encoding gradient
            seq.add_block(g_ro_pre, g_pe_pre, g_ss_reph)  # Add a combination of ro rewinder, phase encoding, and slice refocusing
            seq.add_block(delay1)  # Delay 1: until 180-deg pulse
            seq.add_block(rf180, g_ss180)  # 180 deg pulse for SE
            seq.add_block(delay2)  # Delay 2: until readout
            seq.add_block(g_ro, adc)  # Readout!
            seq.add_block(delay3)  # Delay 3: until next inversion pulse

    if write:
        seq.write("se_fov{:.0f}mm_Nf{:d}_Np{:d}_TE{:.0f}ms_TR{:.0f}ms.seq".format(fov * 1000, Nf, Np, TE * 1000, TR * 1000))


    print('Spin echo sequence constructed')
    return seq

def make_pulseq_se_oblique(fov,n,thk,fa,tr,te,enc='xyz',slice_locs=None,write=False):
    """Makes a Spin Echo (SE) sequence in any plane

        2D oblique multi-slice Spin-Echo pulse sequence with Cartesian encoding
        Oblique means that each of slice-selection, phase encoding, and frequency encoding
        can point in any specified direction

        Parameters
        ----------
        fov : array_like
            Isotropic field-of-view, or length-2 list [fov_readout, fov_phase], in meters
        n : array_like
            Isotropic matrix size, or length-2 list [n_readout, n_phase]
        thk : float
            Slice thickness in meters
        fa : float
            Flip angle in degrees
        tr : float
            Repetition time in seconds
        te : float
            Echo time in seconds
        enc : str or array_like, optional
            Spatial encoding directions
            1st - readout; 2nd - phase encoding; 3rd - slice select
            - Use str with any permutation of x, y, and z to obtain orthogonal slices
            e.g. The default 'xyz' means axial(z) slice with readout in x and phase encoding in y
            - Use list to indicate vectors in the encoding directions for oblique slices
            They should be perpendicular to each other, but not necessarily unit vectors
            e.g. [(2,1,0),(-1,2,0),(0,0,1)] rotates the two in-plane encoding directions for an axial slice
        slice_locs : array_like, optional
            Slice locations from isocenter in meters
            Default is None which means a single slice at the center
        write : bool, optional
            Whether to write seq into file; default is False

        Returns
        -------
        seq : Sequence
            Pulse sequence as a Pulseq object

        """

    # System options
    system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130,
                  slew_unit='T/m/s', rf_ringdown_time=30e-6,
                  rf_dead_time=100e-6, adc_dead_time=20e-6)
    seq = Sequence(system)

    # Sequence parameters
    ug_fe, ug_pe, ug_ss = parse_enc(enc)
    Nf, Np = (n,n) if isinstance(n,int) else (n[0], n[1])
    delta_k_ro, delta_k_pe = (1/fov,1/fov) if isinstance(fov,float) else (1/fov[0], 1/fov[1])
    kWidth_ro = Nf * delta_k_ro
    TE,TR = te,tr

    # Non-180 pulse
    flip1 = fa * pi / 180
    rf, g_ss, __ = make_sinc_pulse(flip_angle=flip1, system=system, duration=2e-3, slice_thickness=thk,
                                apodization=0.5, time_bw_product=4)
    g_ss_x, g_ss_y, g_ss_z = make_oblique_gradients(g_ss,ug_ss)

    # 180 pulse
    flip2 = 180 * pi / 180
    rf180, g_ss180, __ = make_sinc_pulse(flip_angle=flip2, system=system, duration=2e-3, slice_thickness=thk,
                                     apodization=0.5, time_bw_product=4)
    g_ss180_x, g_ss180_y, g_ss180_z = make_oblique_gradients(g_ss180,ug_ss)

    # Readout gradient & ADC
    readoutTime = 6.4e-3
    g_ro = make_trapezoid(channel='x',system=system,flat_area=kWidth_ro, flat_time=readoutTime)
    g_ro_x, g_ro_y, g_ro_z = make_oblique_gradients(g_ro, ug_fe)
    adc = make_adc(num_samples=Nf, system=system, duration=g_ro.flat_time, delay=g_ro.rise_time)

    # RO rewinder gradient
    g_ro_pre = make_trapezoid(channel='x',system=system,area=g_ro.area/2,duration=2e-3)
    g_ro_pre_x, g_ro_pre_y, g_ro_pre_z = make_oblique_gradients(g_ro_pre, ug_fe)

    # Slice refocusing gradient
    g_ss_reph = make_trapezoid(channel='z',system=system,area=-g_ss.area/2,duration=2e-3)
    g_ss_reph_x, g_ss_reph_y, g_ss_reph_z = make_oblique_gradients(g_ss_reph, ug_ss)

    # Delays
    delayTE1 = (TE - 2*max(calc_duration(g_ss_reph), calc_duration(g_ro_pre)) - calc_duration(g_ss) - calc_duration(
        g_ss180))/2
    delayTE2 = (TE - calc_duration(g_ro) - calc_duration(g_ss180))/2
    delayTE3 = TR - TE - (calc_duration(g_ss) + calc_duration(g_ro)) / 2


    delay1 = make_delay(delayTE1)
    delay2 = make_delay(delayTE2)
    delay3 = make_delay(delayTE3)

    # Construct sequence
    if slice_locs is None:
        locs = [0]
    else:
        locs = slice_locs

    for u in range(len(locs)):
        rf180.freq_offset = g_ss180.amplitude * locs[u]
        rf.freq_offset = g_ss.amplitude * locs[u]
        for i in range(Np):
            seq.add_block(rf, g_ss_x, g_ss_y, g_ss_z)  # 90-deg pulse
            g_pe = make_trapezoid(channel='y',system=system,area=-(Np/2 - i)*delta_k_pe, duration=2e-3)  # Phase encoding gradient
            g_pe_x, g_pe_y, g_pe_z = make_oblique_gradients(g_pe, ug_pe)

            pre_grads_list = [g_ro_pre_x, g_ro_pre_y, g_ro_pre_z,
                              g_ss_reph_x, g_ss_reph_y, g_ss_reph_z,
                              g_pe_x, g_pe_y, g_pe_z]
            gtx, gty, gtz = combine_trap_grad_xyz(pre_grads_list, system, 2e-3)


            seq.add_block(gtx,gty,gtz)  # Add a combination of ro rewinder, phase encoding, and slice refocusing
            seq.add_block(delay1)  # Delay 1: until 180-deg pulse
            seq.add_block(rf180, g_ss180_x, g_ss180_y, g_ss180_z)  # 180 deg pulse for SE
            seq.add_block(delay2)  # Delay 2: until readout
            seq.add_block(g_ro_x, g_ro_y, g_ro_z, adc)  # Readout!
            seq.add_block(delay3)  # Delay 3: until next inversion pulse

    if write:
        seq.write("se_fov{:.0f}mm_Nf{:d}_Np{:d}_TE{:.0f}ms_TR{:.0f}ms_FA{:d}deg.seq".format(fov * 1000, Nf, Np, TE * 1000, TR * 1000, fa))


    print('Spin echo sequence (oblique) constructed')
    return seq


# TODO multi-shot epi needs to be tested on scanner! : )
def make_pulseq_epi_oblique(fov,n,thk,fa,tr,te,enc='xyz',slice_locs=None,echo_type="se",n_shots=1,seg_type='blocked',write=False):
    """Makes an Echo Planar Imaging (EPI) sequence in any plane

        2D oblique multi-slice EPI pulse sequence with Cartesian encoding
        Oblique means that each of slice-selection, phase encoding, and frequency encoding
        can point in any specified direction

        Parameters
        ----------
        fov : array_like
            Isotropic field-of-view, or length-2 list [fov_readout, fov_phase], in meters
        n : array_like
            Isotropic matrix size, or length-2 list [n_readout, n_phase]
        thk : float
            Slice thickness in meters
        fa : float
            Flip angle in degrees
        tr : float
            Repetition time in seconds
        te : float
            Echo time in seconds
        enc : str or array_like, optional
            Spatial encoding directions
            1st - readout; 2nd - phase encoding; 3rd - slice select
            - Use str with any permutation of x, y, and z to obtain orthogonal slices
            e.g. The default 'xyz' means axial(z) slice with readout in x and phase encoding in y
            - Use list to indicate vectors in the encoding directions for oblique slices
            They should be perpendicular to each other, but not necessarily unit vectors
            e.g. [(2,1,0),(-1,2,0),(0,0,1)] rotates the two in-plane encoding directions for an axial slice
        slice_locs : array_like, optional
            Slice locations from isocenter in meters
            Default is None which means a single slice at the center
        echo_type : str, optional {'se','gre'}
            Type of echo generated
            se (default) - spin echo (an 180 deg pulse is used)
            gre - gradient echo
        n_shots : int, optional
            Number of shots used to encode each slicel; default is 1
        seg_type : str, optional {'blocked','interleaved'}
            Method to divide up k-space in the case of n_shots > 1; default is 'blocked'
            'blocked' - each shot covers a rectangle, with no overlap between shots
            'interleaved' - each shot samples the full k-space but with wider phase steps

        write : bool, optional
            Whether to write seq into file; default is False

        Returns
        -------
        seq : Sequence
            Pulse sequence as a Pulseq object
        ro_dirs : numpy.ndarray
            List of 0s and 1s indicating direction of readout
            0 - left to right
            1 - right to left (needs to be reversed at recon)
        ro_order : numpy.ndarray
            Order in which to re-arrange the readout lines
            It is [] for blocked acquisition (retain original order)

        """
    # Multi-slice, multi-shot (>=1)
    # TE is set to be where the trajectory crosses the center of k-space

    # System options
    system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130,
                  slew_unit='T/m/s', rf_ringdown_time=30e-6,
                  rf_dead_time=100e-6, adc_dead_time=20e-6)
    seq = Sequence(system)


    ug_fe, ug_pe, ug_ss = parse_enc(enc)

    # Sequence parameters
    Nf, Np = (n,n) if isinstance(n,int) else (n[0], n[1])
    delta_k_ro, delta_k_pe = (1/fov,1/fov) if isinstance(fov,float) else (1/fov[0], 1/fov[1])
    kWidth_ro = Nf * delta_k_ro
    TE,TR = te,tr
    flip = fa * pi / 180

    # RF Pulse (first)
    rf, g_ss, __ = make_sinc_pulse(flip_angle=flip, system=system,duration=2.5e-3, slice_thickness=thk,
                               apodization=0.5, time_bw_product=4)
    g_ss_x, g_ss_y, g_ss_z = make_oblique_gradients(g_ss, ug_ss)


    # Readout gradients
#    readoutTime = Nf * 4e-6
    dwell=1e-5
    readoutTime = Nf*dwell
    g_ro_pos = make_trapezoid(channel='x',system=system,flat_area=kWidth_ro,flat_time=readoutTime)
    g_ro_pos_x, g_ro_pos_y, g_ro_pos_z = make_oblique_gradients(g_ro_pos,ug_fe)
    g_ro_neg = copy.deepcopy(g_ro_pos)
    modify_gradient(g_ro_neg,scale=-1)
    g_ro_neg_x, g_ro_neg_y, g_ro_neg_z = make_oblique_gradients(g_ro_neg,ug_fe)

    # TODO make sure delay is a multiple of gradient raster time
#    adc = make_adc(num_samples=Nf, system=system, duration=g_ro_pos.flat_time, delay=g_ro_pos.rise_time+dwell/2)
    adc = make_adc(num_samples=Nf, system=system, duration=g_ro_pos.flat_time, delay=g_ro_pos.rise_time)
    print("ADC delay: ", adc.delay)

    pre_time = 8e-4

    # 180 deg pulse for SE
    if echo_type == "se":
        # RF Pulse (180 deg for SE)
        flip180 = 180 * pi / 180
        rf180, g_ss180, __ = make_sinc_pulse(flip_angle=flip180, system=system,duration=2.5e-3,slice_thickness=thk,
                                         apodization=0.5, time_bw_product=4)
        g_ss180_x, g_ss180_y, g_ss180_z = make_oblique_gradients(g_ss180, ug_ss)

        # Slice-select direction spoilers
        g_ss_spoil = make_trapezoid(channel='z',system=system,area=g_ss.area*2,duration=3*pre_time)
        ##
        modify_gradient(g_ss_spoil,0)
        ##
        g_ss_spoil_x, g_ss_spoil_y, g_ss_spoil_z = make_oblique_gradients(g_ss_spoil, ug_ss)

    # Readout rewinder
    ro_pre_area = g_ro_neg.area / 2 if echo_type == 'gre' else g_ro_pos.area / 2
    g_ro_pre = make_trapezoid(channel='x',system=system, area=ro_pre_area, duration=pre_time)
    g_ro_pre_x, g_ro_pre_y, g_ro_pre_z = make_oblique_gradients(g_ro_pre,ug_fe)

    # Slice-selective rephasing
    g_ss_reph = make_trapezoid(channel='z',system=system,area=-g_ss.area/2,duration=pre_time)
    g_ss_reph_x, g_ss_reph_y, g_ss_reph_z = make_oblique_gradients(g_ss_reph, ug_ss)

    # Phase encode rewinder
    if echo_type == 'gre':
        pe_max_area = (Np/2)*delta_k_pe
    elif echo_type == 'se':
        pe_max_area = -(Np/2)*delta_k_pe

    g_pe_max = make_trapezoid(channel='y',system=system,area=pe_max_area,duration=pre_time)

    # Phase encoding blips
    dur = ceil(2 * sqrt(delta_k_pe/ system.max_slew) / 10e-6) * 10e-6
    g_blip = make_trapezoid(channel='y',system=system,area=delta_k_pe,duration=dur)

    # Delays
    duration_to_center = (Np/ 2 ) * calc_duration(g_ro_pos) + (Np-1) / 2 * calc_duration(g_blip) # why?

    if echo_type == 'se':
        delayTE1 = TE / 2 - calc_duration(g_ss) / 2 - pre_time - calc_duration(g_ss_spoil) - calc_duration(rf180) / 2
        delayTE2 = TE / 2 - calc_duration(rf180) / 2 - calc_duration(g_ss_spoil) - duration_to_center
        delay1 = make_delay(delayTE1)
        delay2 = make_delay(delayTE2)
    elif echo_type == 'gre':
        delayTE = TE - calc_duration(g_ss)/2 - pre_time - duration_to_center
        delay12 = make_delay(delayTE)

    delayTR = TR - TE - calc_duration(rf) / 2 - duration_to_center
    delay3 = make_delay(delayTR) # This might be different for each rep though. Fix later

#####################################################################################################
    # Multi-shot calculations
    ro_dirs = []
    ro_order = []

    # Find number of lines in each block

    if seg_type == 'blocked':

        # Number of lines in each full readout block
        nl = ceil(Np / n_shots)

        # Number of k-space lines per readout
        if Np%nl == 0:
            nlines_list = nl*np.ones(n_shots)
        else:
            nlines_list = nl*np.ones(n_shots-1)
            nlines_list = np.append(nlines_list,Np%nl)

        pe_scales = 2*np.append([0],np.cumsum(nlines_list)[:-1])/Np - 1
        g_blip_x, g_blip_y, g_blip_z = make_oblique_gradients(g_blip, ug_pe)
        for nlines in nlines_list:
            ro_dirs = np.append(ro_dirs, ((-1)**(np.arange(0,nlines)+1)+1)/2)


    elif seg_type == 'interleaved':
        # Minimum number of lines per readout
        nb = floor(Np / n_shots)

        # Number of k-space lines per readout
        nlines_list = np.ones(n_shots)*nb
        nlines_list[:Np%n_shots] += 1

        # Phase encoding scales (starts from -1; i.e. bottom left combined with pre-readout)
        pe_scales = 2*np.arange(0,(Np-n_shots)/Np,1/Np)[0:n_shots]-1
        print(pe_scales)
        # Larger blips
        modify_gradient(g_blip, scale=n_shots)
        g_blip_x, g_blip_y, g_blip_z = make_oblique_gradients(g_blip, ug_pe)

#        ro_order = np.reshape(np.reshape(np.arange(0,Np),(),order='F'),(0,Np))

        ro_order = np.zeros((nb+1,n_shots))
        ro_inds = np.arange(Np)
        # Readout order for recon
        for k in range(n_shots):
            cs = int(nlines_list[k])
            ro_order[:cs,k] = ro_inds[:cs]
            ro_inds = np.delete(ro_inds,range(cs))
        ro_order = ro_order.flatten()[:Np].astype(int)

        np.save('readout_order_for_interleaving.npy', ro_order)

        print(ro_order)

        # Readout directions in original (interleaved) order
        for nlines in nlines_list:
            ro_dirs = np.append(ro_dirs, ((-1)**(np.arange(0,nlines)+1)+1)/2)



#####################################################################################################

    # Add blocks

    for u in range(len(slice_locs)): # For each slice
        # Offset rf
        rf.freq_offset = g_ss.amplitude * slice_locs[u]
        for v in range(n_shots):
            # Find init. phase encode
            g_pe = copy.deepcopy(g_pe_max)
            modify_gradient(g_pe, pe_scales[v])
            g_pe_x, g_pe_y, g_pe_z = make_oblique_gradients(g_pe, ug_pe)
            # First RF
            seq.add_block(rf, g_ss_x, g_ss_y, g_ss_z)
            # Pre-winder gradients
            pre_grads_list = [g_ro_pre_x, g_ro_pre_y, g_ro_pre_z,
                              g_pe_x, g_pe_y, g_pe_z,
                              g_ss_reph_x, g_ss_reph_y, g_ss_reph_z]
            gtx, gty, gtz = combine_trap_grad_xyz(pre_grads_list, system, pre_time)
            seq.add_block(gtx, gty, gtz)


            # 180 deg pulse and spoilers, only for Spin Echo
            if echo_type == 'se':
                # First delay
                seq.add_block(delay1)
                # Second RF : 180 deg with spoilers on both sides
                seq.add_block(g_ss_spoil_x, g_ss_spoil_y, g_ss_spoil_z)#why?
                seq.add_block(rf180, g_ss180_x, g_ss180_y, g_ss180_z)
                seq.add_block(g_ss_spoil_x, g_ss_spoil_y, g_ss_spoil_z)
                # Delay between rf180 and beginning of readout
                seq.add_block(delay2)
            # For gradient echo it's just a delay
            elif echo_type == 'gre':
                seq.add_block(delay12)



            # EPI readout with blips
            for i in range(int(nlines_list[v])):
                if i%2 == 0:
                    seq.add_block(g_ro_pos_x, g_ro_pos_y, g_ro_pos_z, adc) # ro line in the positive direction
                else:
                    seq.add_block(g_ro_neg_x, g_ro_neg_y, g_ro_neg_z, adc) # ro line backwards
                seq.add_block(g_blip_x, g_blip_y, g_blip_z) # blip

            seq.add_block(delay3)

    # Display 1 TR
    #seq.plot(time_range=(0, TR))

    if write:
        seq.write("epi_{}_FOV{:.0f}mm_Nf{:d}_Np{:d}_TE{:.0f}ms_TR{:.0f}ms_FA{:d}deg_{:d}shots.seq"\
                  .format(echo_type, fov*1000, Nf, Np, TE * 1000, TR * 1000, fa, n_shots))


    print('EPI sequence (oblique) constructed')
    return seq, ro_dirs, ro_order




def parse_enc(enc):
    """Helper function for decoding enc parameter

    Parameters
    ----------
    enc : str or array_like
        Inputted encoding scheme to parse
    Returns
    -------
    ug_fe : numpy.ndarray
        Length-3 vector of readout direction
    ug_pe : numpy.ndarray
        Length-3 vector of phase encoding direction
    ug_ss : numpy.ndarray
        Length-3 vector of slice selecting direction

    """
    if isinstance(enc, str):
        xyz_dict = {'x': (1, 0, 0), 'y': (0, 1, 0), 'z': (0, 0, 1)}
        ug_fe = xyz_dict[enc[0]]
        ug_pe = xyz_dict[enc[1]]
        ug_ss = xyz_dict[enc[2]]
    else:
        ug_fe = np.array(enc[0])
        ug_pe = np.array(enc[1])
        ug_ss = np.array(enc[2])

        ug_fe = ug_fe / np.linalg.norm(ug_fe)
        ug_pe = ug_pe / np.linalg.norm(ug_pe)
        ug_ss = ug_ss / np.linalg.norm(ug_ss)

    print('ug_fe: ', ug_fe)
    print('ug_pe: ', ug_pe)
    print('ug_ss: ', ug_ss)

    return ug_fe, ug_pe, ug_ss


def make_oblique_gradients(gradient,unit_grad):
    """Helper function to make oblique gradients

    (Gx, Gy, Gz) are generated from a single orthogonal gradient
    and a direction indicated by unit vector

    Parameters
    ----------
    gradient : Gradient
        Pulseq gradient object
    unit_grad: array_like
        Length-3 unit vector indicating direction of resulting oblique gradient

    Returns
    -------
    ngx, ngy, ngz : Gradient
        Oblique gradients in x, y, and z directions

    """
    ngx = copy.deepcopy(gradient)
    ngy = copy.deepcopy(gradient)
    ngz = copy.deepcopy(gradient)

    modify_gradient(ngx, unit_grad[0],'x')
    modify_gradient(ngy, unit_grad[1],'y')
    modify_gradient(ngz, unit_grad[2],'z')

    return ngx, ngy, ngz

def modify_gradient(gradient,scale,channel=None):
    """Helper function to modify the strength and channel of an existing gradient

    Parameters
    ----------
    gradient : Gradient
        Pulseq gradient object to be modified
    scale : float
        Scalar to multiply the gradient strength by
    channel : str, optional {None, 'x','y','z'}
        Channel to switch gradient into
        Default is None which keeps the original channel

    """
    gradient.amplitude *= scale
    gradient.area *= scale
    if gradient.type == 'trap':
        gradient.flat_area *= scale
    if channel != None:
        gradient.channel = channel


def combine_trap_grad_xyz(gradients,system,dur):
    """Helper function that merges multiple gradients

    A list of gradients are combined into one set of 3 oblique gradients (Gx, Gy, Gz) with equivalent areas
    Note that the waveforms are not preserved : the outputs will always be trapezoidal gradients

    Parameters
    ----------
    gradients : list
        List of gradients to be combined; there can be any number of x, y, or z gradients
    system : Opts
        Pulseq object that indicates system constraints for gradient parameters
    dur : float
        Duration of the output oblique gradients

    Returns
    -------
    gtx, gty, gtz : Gradient
        Oblique pulseq gradients with equivalent areas to all input gradients combined

    """
    gx_area, gy_area, gz_area = (0,0,0)
    for g in gradients:
        if g.channel == 'x':
            gx_area += g.area
        elif g.channel == 'y':
            gy_area += g.area
        elif g.channel == 'z':
            gz_area += g.area


    gtx = make_trapezoid(channel='x',system=system,area=gx_area,duration=dur)
    gty = make_trapezoid(channel='y',system=system,area=gy_area,duration=dur)
    gtz = make_trapezoid(channel='z',system=system,area=gz_area,duration=dur)

    return gtx, gty, gtz