physics_based_augmentation/seqgen/utilities/standart_RF.py

# -*- coding: utf-8 -*-
"""
A subroutine functions to create different excitation and refocusing
pulses accompanied by combined SS and spoil gradients.
Requires the params structure as input.

@author: petrm
"""


# import
from pypulseq.make_sinc_pulse import make_sinc_pulse
from pypulseq.make_trapezoid import make_trapezoid
from pypulseq.make_extended_trapezoid import make_extended_trapezoid
from pypulseq.calc_duration import calc_duration
import numpy as np


# def tse_excitation_grad(params, scanner_parameters, area_gz_spoil, flip180):
#
#    return

def refocusing_grad(params, scanner_parameters, gz90_area, flip180, rf180_phase, t_refwd, spoil_duration, united: bool):
    # Create 180 degree SS refocusing pulse with SS and spoiled gradients
    rf180, gz_ref, _ = make_sinc_pulse(
        flip_angle=flip180,
        system=scanner_parameters,
        duration=params.t_ref,
        slice_thickness=params.sl_thkn,
        apodization=params.apodization,
        time_bw_product=round(params.t_BW_product_ref, 8),
        phase_offset=rf180_phase,
        use="refocusing",
        return_gz=True, 
        delay = scanner_parameters.rf_dead_time
    )


    gz180 = make_trapezoid(channel="z", system=scanner_parameters, amplitude=gz_ref.amplitude,
                           flat_time=t_refwd, rise_time=params.dG)

    if spoil_duration == 'min':
        gz_spoil_dur_min = params.dG + (gz90_area)/scanner_parameters.max_grad
        gz_spoil_dur_min = np.ceil(gz_spoil_dur_min / scanner_parameters.grad_raster_time) * scanner_parameters.grad_raster_time

        if gz_spoil_dur_min < 2 * params.dG:
            gz_spoil_dur_min = 2 * params.dG + scanner_parameters.grad_raster_time

        gz_spoil1 = make_trapezoid(channel='z', system=scanner_parameters, area= gz90_area,
                                   duration=gz_spoil_dur_min, rise_time=params.dG)
        gz_spoil2 = make_trapezoid(channel='z', system=scanner_parameters, area= gz90_area,
                                   duration=gz_spoil_dur_min, rise_time=params.dG)

    else:
        spoil_duration = float(spoil_duration)
        spoil_duration = np.ceil(spoil_duration / scanner_parameters.grad_raster_time) * scanner_parameters.grad_raster_time

        gz_spoil1 = make_trapezoid(channel='z', system=scanner_parameters, area= gz90_area,
                                   duration=spoil_duration, rise_time=params.dG)
        gz_spoil2 = make_trapezoid(channel='z', system=scanner_parameters, area= gz90_area,
                                   duration=spoil_duration, rise_time=params.dG)

    # SS refocusing gradient with spoilers

    gz_sp1_times = np.array(
        [
            0,
            gz_spoil1.rise_time,
            gz_spoil1.rise_time + gz_spoil1.flat_time,
            gz_spoil1.rise_time + gz_spoil1.flat_time + gz_spoil1.fall_time
        ]
    )
    gz_sp1_amp = np.array(
        [
            0,
            gz_spoil1.amplitude,
            gz_spoil1.amplitude,
            gz180.amplitude
        ]
    )
    gz_sp1 = make_extended_trapezoid(channel='z', system=scanner_parameters, times=gz_sp1_times, amplitudes=gz_sp1_amp)

    gz_sp2_times = np.array(
        [
            0,
            gz180.flat_time
        ]
    )
    gz_sp2_amp = np.array(
        [
            gz180.amplitude,
            gz180.amplitude
        ]
    )
    gz_sp2 = make_extended_trapezoid(channel='z', system=scanner_parameters, times=gz_sp2_times, amplitudes=gz_sp2_amp)

    gz_sp3_times = np.array(
        [
            0,
            gz_spoil2.rise_time,
            gz_spoil2.rise_time + gz_spoil2.flat_time,
            gz_spoil2.rise_time + gz_spoil2.flat_time + gz_spoil2.fall_time
        ]
    )

    gz_sp3_amp = np.array(
        [
            gz180.amplitude,
            gz_spoil2.amplitude,
            gz_spoil2.amplitude,
            0
        ]
    )
    gz_sp3 = make_extended_trapezoid(channel='z', system=scanner_parameters, times=gz_sp3_times, amplitudes=gz_sp3_amp)


    if united:
        return rf180, gz_sp1, gz_sp2, gz_sp3
    else:
        return rf180, gz_spoil1, gz180, gz_spoil2

'''
def readout_grad_x(params, scanner_parameters, spoil_duration, rewind_duration, united: bool, phi_radial):
    """
    merged with spoiler gradients readout block in x direction
    phi_radial - angle for radial aquisition
                 (0 - for cartesian)
    @author: petrm
    """

    # Create readout gradient
    readout_time = round(1 / params.BW_pixel, 8)
    k_read = np.double(params.Nf) / np.double(params.FoV_f)
    t_gx = np.ceil(readout_time / scanner_parameters.grad_raster_time) * scanner_parameters.grad_raster_time
    gx = make_trapezoid(channel='x', system=scanner_parameters, flat_area=k_read * np.cos(phi_radial),
                        flat_time=t_gx + 2 * scanner_parameters.adc_dead_time, rise_time=params.dG)

    # generate gx spoiler gradient - G_crr
    if spoil_duration == 'min':
        gx_spoil = make_trapezoid(channel='x', system=scanner_parameters, area=gx.area,
                                  rise_time=params.dG)
        if gx_spoil.flat_time == 0:
            gx_spoil = make_trapezoid(channel='x', system=scanner_parameters, area=gx.area,
                                      rise_time=params.dG, flat_time=params.dG)

    else:
        spoil_duration = float(spoil_duration)
        spoil_duration = (np.ceil(spoil_duration / scanner_parameters.grad_raster_time)
                          * scanner_parameters.grad_raster_time)
        gx_spoil = make_trapezoid(channel='x', system=scanner_parameters, area=gx.area,
                                  duration=spoil_duration, rise_time=params.dG)

    # readout gradient with spoilers

    gx_sp1_times = np.array(
        [
            0,
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil.rise_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil.rise_time + gx_spoil.flat_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil.rise_time + gx_spoil.flat_time + gx_spoil.fall_time) / scanner_parameters.grad_raster_time)
        ]
    )
    gx_sp1_amp = np.array(
        [
            0,
            gx_spoil.amplitude,
            gx_spoil.amplitude,
            gx.amplitude
        ]
    )
    gx_sp1 = make_extended_trapezoid(channel='x', system=scanner_parameters, times=gx_sp1_times, amplitudes=gx_sp1_amp)

    gx_sp2_times = np.array(
        [
            0,
            scanner_parameters.grad_raster_time * np.round(
                (gx.flat_time) / scanner_parameters.grad_raster_time)
        ]
    )
    gx_sp2_amp = np.array(
        [
            gx.amplitude,
            gx.amplitude
        ]
    )
    gx_sp2 = make_extended_trapezoid(channel='x', system=scanner_parameters, times=gx_sp2_times, amplitudes=gx_sp2_amp)

    gx_sp3_times = np.array(
        [
            0,
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil.rise_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil.rise_time + gx_spoil.flat_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil.rise_time + gx_spoil.flat_time + gx_spoil.fall_time) / scanner_parameters.grad_raster_time)
        ]
    )

    gx_sp3_amp = np.array(
        [
            gx.amplitude,
            gx_spoil.amplitude,
            gx_spoil.amplitude,
            0
        ]
    )
    gx_sp3 = make_extended_trapezoid(channel='x', system=scanner_parameters, times=gx_sp3_times, amplitudes=gx_sp3_amp)

    if rewind_duration == 'min':
        gx_pre = make_trapezoid(channel="x", system=scanner_parameters, area=gx.area * 1.5,
                                rise_time=params.dG)
        if gx_pre.flat_time == 0:
            gx_pre = make_trapezoid(channel='x', system=scanner_parameters, area=gx.area * 1.5,
                                      rise_time=params.dG, flat_time=params.dG)

    else:
        rewind_duration = float(rewind_duration)
        rewind_duration = (np.ceil(rewind_duration / scanner_parameters.grad_raster_time)
                          * scanner_parameters.grad_raster_time)
        gx_pre = make_trapezoid(channel="x", system=scanner_parameters, area=gx.area * 1.5,
                                rise_time=params.dG, duration=rewind_duration)

    if united:
        return gx_sp1, gx_sp2, gx_sp3, t_gx, gx_pre
    else:
        return gx_spoil, gx, gx_spoil, t_gx, gx_pre
'''

def readout_grad_y(params, scanner_parameters, spoil_duration, rewind_duration,
                   united: bool, phi_radial, phase_encoding_area):
    """
    merged with spoiler gradients readout block in y direction
    phi_radial - angle for radial acquisition
                 (0 - for cartesian)
    phase_encoding_area - additional gradient of phase encoding for spoiler
    @author: petrm
    """

    phase_encoding_area = phase_encoding_area * np.cos(phi_radial)
    # Create readout gradient
    readout_time = round(1 / params.BW_pixel, 8)
    k_read = np.double(params.Nf) / np.double(params.FoV_f)
    t_gx = np.ceil(readout_time / scanner_parameters.grad_raster_time) * scanner_parameters.grad_raster_time
    if phi_radial == 0.0:
        gx = make_trapezoid(channel='y', system=scanner_parameters, area=0,
                            flat_time=t_gx + 2 * scanner_parameters.adc_dead_time, rise_time=params.dG)
    else:
        gx = make_trapezoid(channel='y', system=scanner_parameters, flat_area=k_read * np.sin(phi_radial),
                            flat_time=t_gx + 2 * scanner_parameters.adc_dead_time, rise_time=params.dG)

    # generate gx spoiler gradient - G_crr
    if spoil_duration == 'min':
        gx_spoil_dur_min = params.dG + (gx.area + abs(phase_encoding_area))/scanner_parameters.max_grad
        gx_spoil_dur_min = np.ceil(gx_spoil_dur_min / scanner_parameters.grad_raster_time) * scanner_parameters.grad_raster_time
        if gx_spoil_dur_min < 2 * params.dG:
            gx_spoil_dur_min = 2 * params.dG + scanner_parameters.grad_raster_time

        gx_spoil1 = make_trapezoid(channel='y', system=scanner_parameters, area=gx.area + phase_encoding_area,
                                   duration=gx_spoil_dur_min, rise_time=params.dG)
        gx_spoil2 = make_trapezoid(channel='y', system=scanner_parameters, area=gx.area - phase_encoding_area,
                                   duration=gx_spoil_dur_min, rise_time=params.dG)

    else:
        spoil_duration = float(spoil_duration)
        spoil_duration = (np.ceil(spoil_duration / scanner_parameters.grad_raster_time)
                          * scanner_parameters.grad_raster_time)

        gx_spoil1 = make_trapezoid(channel='y', system=scanner_parameters, area=gx.area + phase_encoding_area,
                                   duration=spoil_duration, rise_time=params.dG)
        gx_spoil2 = make_trapezoid(channel='y', system=scanner_parameters, area=gx.area - phase_encoding_area,
                                   duration=spoil_duration, rise_time=params.dG)

    # readout gradient with spoilers

    gx_sp1_times = np.array(
        [
            0,
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil1.rise_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil1.rise_time + gx_spoil1.flat_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil1.rise_time + gx_spoil1.flat_time + gx_spoil1.fall_time) / scanner_parameters.grad_raster_time)
        ]
    )
    gx_sp1_amp = np.array(
        [
            0,
            gx_spoil1.amplitude,
            gx_spoil1.amplitude,
            gx.amplitude
        ]
    )
    gx_sp1 = make_extended_trapezoid(channel='y', system=scanner_parameters, times=gx_sp1_times, amplitudes=gx_sp1_amp)

    gx_sp2_times = np.array(
        [
            0,
            scanner_parameters.grad_raster_time * np.round(
                (gx.flat_time) / scanner_parameters.grad_raster_time)
        ]
    )
    gx_sp2_amp = np.array(
        [
            gx.amplitude,
            gx.amplitude
        ]
    )
    gx_sp2 = make_extended_trapezoid(channel='y', system=scanner_parameters, times=gx_sp2_times, amplitudes=gx_sp2_amp)

    gx_sp3_times = np.array(
        [
            0,
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil2.rise_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil2.rise_time + gx_spoil2.flat_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil2.rise_time + gx_spoil2.flat_time + gx_spoil2.fall_time) / scanner_parameters.grad_raster_time)
        ]
    )

    gx_sp3_amp = np.array(
        [
            gx.amplitude,
            gx_spoil2.amplitude,
            gx_spoil2.amplitude,
            0
        ]
    )
    gx_sp3 = make_extended_trapezoid(channel='y', system=scanner_parameters, times=gx_sp3_times, amplitudes=gx_sp3_amp)

    if rewind_duration == 'min':
        gx_pre_dur_min = params.dG + (gx.area * 1.5) / scanner_parameters.max_grad
        gx_pre_dur_min = np.ceil(gx_pre_dur_min / scanner_parameters.grad_raster_time) * scanner_parameters.grad_raster_time
        if gx_pre_dur_min < 2 * params.dG:
            gx_pre_dur_min = 2 * params.dG + scanner_parameters.grad_raster_time

        gx_pre = make_trapezoid(channel="y", system=scanner_parameters, area=gx.area * 1.50,
                                duration=gx_pre_dur_min, rise_time=params.dG)

    else:
        rewind_duration = float(rewind_duration)
        rewind_duration = (np.ceil(rewind_duration / scanner_parameters.grad_raster_time)
                          * scanner_parameters.grad_raster_time)
        gx_pre = make_trapezoid(channel="y", system=scanner_parameters, area=gx.area * 1.50,
                                rise_time=params.dG, duration=rewind_duration)

    if united:
        return gx_sp1, gx_sp2, gx_sp3, t_gx, gx_pre
    else:
        return gx_spoil1, gx, gx_spoil2, t_gx, gx_pre


def readout_grad_x(params, scanner_parameters, spoil_duration, rewind_duration,
                   united: bool, phi_radial, phase_encoding_area, **kwargs):
    """
    merged with spoiler gradients readout block in y direction
    phi_radial - angle for radial acquisition
                 (0 - for cartesian)
    phase_encoding_area - additional gradient of phase encoding for spoiler
    @author: petrm
    """
    if 'dixon_delta_te' in kwargs:
        dixon_delta_te = kwargs['dixon_delta_te'] #TODO: add exceptions??
    else:
        dixon_delta_te = 0

    phase_encoding_area = phase_encoding_area * np.sin(phi_radial) * (-1)

    # Create readout gradient
    readout_time = round(1 / params.BW_pixel, 8)
    k_read = np.double(params.Nf) / np.double(params.FoV_f)
    t_gx = np.ceil(readout_time / scanner_parameters.grad_raster_time) * scanner_parameters.grad_raster_time
    gx = make_trapezoid(channel='x', system=scanner_parameters, flat_area=k_read * np.cos(phi_radial),
                        flat_time=t_gx + 2 * scanner_parameters.adc_dead_time, rise_time=params.dG)
    if vars(params).get('Dixon') is not None:
        if params.Dixon == True:
            Dixon_prephasor_area = gx.amplitude * dixon_delta_te
        else:
            Dixon_prephasor_area = 0

    # generate gx spoiler gradient - G_crr
    if spoil_duration == 'min':
        gx_spoil_dur_min = 2 * params.dG + (gx.area + abs(phase_encoding_area) + abs(Dixon_prephasor_area))/scanner_parameters.max_grad
        gx_spoil_dur_min = np.ceil(gx_spoil_dur_min / scanner_parameters.grad_raster_time) * scanner_parameters.grad_raster_time
        if gx_spoil_dur_min < 2 * params.dG:
            gx_spoil_dur_min = 2 * params.dG + scanner_parameters.grad_raster_time



        gx_spoil2 = make_trapezoid(channel='x', system=scanner_parameters, area=gx.area + phase_encoding_area + Dixon_prephasor_area,
                                   rise_time=params.dG, duration=gx_spoil_dur_min)
        gx_spoil1 = make_trapezoid(channel='x', system=scanner_parameters, area=gx.area - phase_encoding_area - Dixon_prephasor_area,
                                   rise_time=params.dG, duration=gx_spoil_dur_min)

    else:
        spoil_duration = float(spoil_duration)
        spoil_duration = (np.ceil(spoil_duration / scanner_parameters.grad_raster_time)
                          * scanner_parameters.grad_raster_time)
        gx_spoil1 = make_trapezoid(channel='x', system=scanner_parameters, area=gx.area + phase_encoding_area + Dixon_prephasor_area,
                                   duration=spoil_duration, rise_time=params.dG)
        gx_spoil2 = make_trapezoid(channel='x', system=scanner_parameters, area=gx.area - phase_encoding_area - Dixon_prephasor_area,
                                   duration=spoil_duration, rise_time=params.dG)

    # readout gradient with spoilers

    gx_sp1_times = np.array(
        [
            0,
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil1.rise_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil1.rise_time + gx_spoil1.flat_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil1.rise_time + gx_spoil1.flat_time + gx_spoil1.fall_time) / scanner_parameters.grad_raster_time)
        ]
    )
    gx_sp1_amp = np.array(
        [
            0,
            gx_spoil1.amplitude,
            gx_spoil1.amplitude,
            gx.amplitude
        ]
    )
    gx_sp1 = make_extended_trapezoid(channel='x', system=scanner_parameters, times=gx_sp1_times, amplitudes=gx_sp1_amp)

    gx_sp2_times = np.array(
        [
            0,
            scanner_parameters.grad_raster_time * np.round(
                (gx.flat_time) / scanner_parameters.grad_raster_time)
        ]
    )
    gx_sp2_amp = np.array(
        [
            gx.amplitude,
            gx.amplitude
        ]
    )
    gx_sp2 = make_extended_trapezoid(channel='x', system=scanner_parameters, times=gx_sp2_times, amplitudes=gx_sp2_amp)

    gx_sp3_times = np.array(
        [
            0,
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil2.rise_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil2.rise_time + gx_spoil2.flat_time) / scanner_parameters.grad_raster_time),
            scanner_parameters.grad_raster_time * np.round(
                (gx_spoil2.rise_time + gx_spoil2.flat_time + gx_spoil2.fall_time) / scanner_parameters.grad_raster_time)
        ]
    )

    gx_sp3_amp = np.array(
        [
            gx.amplitude,
            gx_spoil2.amplitude,
            gx_spoil2.amplitude,
            0
        ]
    )
    gx_sp3 = make_extended_trapezoid(channel='x', system=scanner_parameters, times=gx_sp3_times, amplitudes=gx_sp3_amp)

    if rewind_duration == 'min':

        gx_pre_dur_min = params.dG + (gx.area * 1.5)/scanner_parameters.max_grad
        gx_pre_dur_min = np.ceil(gx_pre_dur_min / scanner_parameters.grad_raster_time) * scanner_parameters.grad_raster_time
        if gx_pre_dur_min < 2 * params.dG:
            gx_pre_dur_min = 2 * params.dG + scanner_parameters.grad_raster_time

        gx_pre = make_trapezoid(channel="x", system=scanner_parameters, area=gx.area * 1.50,
                                duration = gx_pre_dur_min, rise_time=params.dG)

    else:
        rewind_duration = float(rewind_duration)
        rewind_duration = (np.ceil(rewind_duration / scanner_parameters.grad_raster_time)
                          * scanner_parameters.grad_raster_time)
        gx_pre = make_trapezoid(channel="x", system=scanner_parameters, area=gx.area * 1.50,
                                rise_time=params.dG, duration=rewind_duration)

    if united:
        return gx_sp1, gx_sp2, gx_sp3, t_gx, gx_pre
    else:
        return gx_spoil1, gx, gx_spoil2, t_gx, gx_pre