physics_based_augmentation/seqgen/MRI_sequences_files/TSE/seqgen_TSE_NIRSII.py

# ---------------------------------------------------------------------
# imports of the libraries
# ---------------------------------------------------------------------

import numpy as np
import math

# Важно: класс Sequence конфликтует по имени с подпакетом seqgen.pypulseq.Sequence,
# что может приводить к ошибке "'module' object is not callable" при вызове Sequence(...)
# Импортируем класс напрямую из файла определения, чтобы избежать коллизии с модулем.
from seqgen.pypulseq.Sequence.sequence import Sequence
from seqgen.pypulseq.calc_rf_center import calc_rf_center
from seqgen.pypulseq.calc_duration import calc_duration
from seqgen.pypulseq.make_adc import make_adc
from seqgen.pypulseq.make_delay import make_delay
from seqgen.pypulseq.make_sinc_pulse import make_sinc_pulse
from seqgen.pypulseq.make_trapezoid import make_trapezoid
from seqgen.pypulseq.opts import Opts

from seqgen.utilities import phase_grad_utils as pgu
from seqgen.utilities.magn_prep.magn_prep import magn_prep_calc, magn_prep_add

from seqgen.utilities.TSE_k_space_fill import TSE_k_space_fill
from seqgen.utilities.standart_RF import refocusing_grad
from seqgen.utilities.standart_RF import readout_grad_x


def seqgen_TSE_NIRSII(params):


    # --------------------------
    # Set system limits
    # --------------------------

    scanner_parameters = Opts(
        max_grad=params.G_amp_max, grad_unit="Hz/m",
        max_slew=params.G_slew_max, slew_unit="Hz/m/s",
        rf_ringdown_time=params.rf_ringdown_time,
        rf_dead_time=params.rf_dead_time,
        adc_dead_time=params.adc_dead_time,
        rf_raster_time=params.rf_raster_time,
        grad_raster_time=params.grad_raster_time,
        block_duration_raster=params.grad_raster_time,
        adc_raster_time=1 / (params.BW_pixel * params.Nf)
    )
    # Защита от редких коллизий: если по какой-то причине импортировалcя модуль, а не класс,
    # явно подтянем класс Sequence и используем его.
    SequenceClass = Sequence
    try:
        # type: ignore[attr-defined]
        _is_class = hasattr(SequenceClass, '__call__') and hasattr(SequenceClass, '__dict__')
        if not _is_class:
            from seqgen.pypulseq.Sequence.sequence import Sequence as _Seq
            SequenceClass = _Seq
    except Exception:
        from seqgen.pypulseq.Sequence.sequence import Sequence as _Seq
        SequenceClass = _Seq

    seq = SequenceClass(scanner_parameters)
    # --------------------
    # additional parameters
    # ---------------------
    united = True  # splited of combined gradients

    if params.D_scans == 0:
        D_scans = False
    else:
        D_scans = True

    params.dixon_delta_te = 1.1e-3 # in ms #TODO add to GUI
    params.Dixon = True
    # --------------------
    # quantify timings
    # ---------------------

    readout_time = round(1 / params.BW_pixel, 8)

    t_exwd = params.t_ex + 2*scanner_parameters.rf_dead_time
    t_refwd = params.t_ref + 2*scanner_parameters.rf_dead_time

    t_sp = 0.5 * (params.ES - readout_time - t_refwd)
    t_sp = scanner_parameters.grad_raster_time * np.round(t_sp / scanner_parameters.grad_raster_time)
    f_sp = 0.75  # area factor for spoiler gradient

    t_spex = 0.5 * (params.ES - t_exwd - t_refwd)
    t_spex = scanner_parameters.grad_raster_time * np.round(t_spex / scanner_parameters.grad_raster_time)
    f_spex = f_sp - 0.5

    
    # =============================================================================
    # RF & SS Gradients
    # =============================================================================

    pulse_bandwidth = np.double(params.t_BW_product_ex) / np.double(params.t_ex)  # In Hz
    gradient_slice_amp = np.double(pulse_bandwidth) / np.double(params.sl_thkn)  # In Hz/m
    pulse_freq_offsets90 = (np.linspace(0.0, params.sl_nb - 1.0, np.int16(params.sl_nb)) - 0.5 * (
            np.double(params.sl_nb) - 1.0)) * (params.sl_thkn * (100.0 + params.sl_gap) / 100.0) * gradient_slice_amp
    # Common gradient parameters

    rf90_phase = np.pi / 2
    rf180_phase = 0
    flip90 = params.FA * np.pi / 180
    flip180 = params.RA * np.pi / 180  # TODO ad parameter to GUI

    rf90, gz_ex, _ = make_sinc_pulse(flip_angle=flip90, system=scanner_parameters, duration=params.t_ex,
                                     slice_thickness=params.sl_thkn, apodization=params.apodization,
                                     time_bw_product=round(params.t_BW_product_ex, 8), phase_offset=rf90_phase,
                                     return_gz=True, use="excitation", delay = scanner_parameters.rf_dead_time)
    rf90.delay = params.dG + scanner_parameters.rf_dead_time

    gz90 = make_trapezoid(channel="z", system=scanner_parameters, amplitude=gz_ex.amplitude,
                          flat_time=t_exwd, rise_time=params.dG)
    gz_reph = make_trapezoid(channel="z", system=scanner_parameters, area=-gz90.area * 0.5,
                             duration=t_spex, rise_time=params.dG)
    # spoil gradient G_sps
    gz_cr = make_trapezoid(channel='z', system=scanner_parameters, area=gz90.area * 4, rise_time=params.dG)


    #Refocusing_gradient
    rf180, gz_spoil1, gz180, gz_spoil2 = refocusing_grad(params, scanner_parameters,
                                                         f_sp * gz90.area, flip180, rf180_phase, t_refwd,
                                                         spoil_duration=t_sp, united=united)
    rf180.delay = scanner_parameters.rf_dead_time
    pulse_freq_offsets180 = (np.linspace(0.0, params.sl_nb - 1.0, np.int16(params.sl_nb)) - 0.5 * (np.double(params.sl_nb) - 1.0)) * (params.sl_thkn * (100.0 + params.sl_gap) / 100.0) * gz180.first

    gx_spoil1, gx, gx_spoil2, t_gx, gx_pre = readout_grad_x(params, scanner_parameters,
                                                            spoil_duration=t_sp,
                                                            rewind_duration=t_spex, united=united,
                                                            phi_radial=0, phase_encoding_area=0,
                                                            dixon_delta_te=0)

    adc = make_adc(num_samples=params.Nf, duration=t_gx,
                   delay=scanner_parameters.adc_dead_time, system=scanner_parameters)
    

    # =============================================================================
    # k-space filling quantification
    # =============================================================================

    k_phase = np.double(params.Np) / np.double(params.FoV_p)
    k_steps_PE = pgu.create_k_steps_PI(k_phase, np.int16(params.Np),
                                       np.int16(params.PI_factor))  # list of phase encoding gradients

    zero_indices = np.where(k_steps_PE == 0.0)[0]
    if zero_indices.size > 0:
        N_center = zero_indices[0]
    else:
        N_center = np.argmin(np.abs(k_steps_PE))


    k_steps_PE = k_steps_PE[0:np.int16(np.round(len(k_steps_PE)))]

    n_ex = math.ceil(len(k_steps_PE) / params.ETL)  # number of excitations with partial Fourier

    # k_space_order_filing = TSE_k_space_fill(n_ex, np.int32(params.ETL), np.int32(params.Np), np.int32(params.N_TE), 'non_linear')  # number of excitations w/o partial Fourier
    k_space_order_filing = TSE_k_space_fill(n_ex, np.int32(params.ETL), len(k_steps_PE),
                                            np.int32(params.N_TE), N_center, 'non_linear',
                                            dummy=D_scans, N_dummy=params.D_scans)
    k_space_save = {'k_space_order': k_space_order_filing}

    magn_prep_objects, magn_prep_dur = magn_prep_calc(vars(params), 
                                                         scanner_parameters)
    # =============================================================================
    # block duration quant
    # =============================================================================
    # TODO:quantify block duration and TE together with d
    block_duration = 0
    block_duration += magn_prep_dur
    block_duration += calc_duration(gz90)
    block_duration += calc_duration(gz_reph)
    #block_duration += params.small_delta + 2 * params.dG
    #block_duration += calc_duration(delay_diff_empt)
    block_duration += gz_spoil1.shape_dur
    block_duration += gz180.shape_dur
    block_duration += gz_spoil2.shape_dur
    #block_duration += calc_duration(delay_diff_empt)
    #block_duration += params.small_delta + 2 * params.dG
    block_duration += gz_spoil1.shape_dur
    for i in range(np.int32(params.ETL) - 1):
        block_duration += gz180.shape_dur
        block_duration += gz_spoil2.shape_dur
        block_duration += gx.shape_dur
        block_duration += gz_spoil1.shape_dur
    block_duration += gz180.shape_dur
    block_duration += gz_spoil2.shape_dur
    block_duration += gx.shape_dur
    block_duration += calc_duration(gz_cr)
    
 
    # =============================================================================
    # Construct concatinations
    # =============================================================================

    eff_time = block_duration  # equal to previous!

    max_slices_per_TR = np.floor(params.TR / eff_time)
    if max_slices_per_TR == 0: #TODO: check GUI
        max_slices_per_TR = 1
    # required_concats = np.int32(np.ceil(params.sl_nb / max_slices_per_TR)) #maximum possible
    required_concats = params.concats
    slice_list1 = list(range(0, np.int32(params.sl_nb), 2)) #TODO: linear and interleaved slices
    slice_list2 = list(range(1, np.int32(params.sl_nb), 2))
    slice_list = slice_list1 + slice_list2
    slice_list = [slice_list[x::required_concats] for x in range(required_concats)]

    # Calculate the TR fillers
    tr_pauses = [(params.TR / np.double(len(x))) - eff_time for x in slice_list]
    tr_pauses = [
        max(seq.grad_raster_time, seq.rf_raster_time) * np.floor(x / max(seq.grad_raster_time, seq.rf_raster_time)) for
        x in tr_pauses]

    # Generate the TR fillers
    tr_fillers = [make_delay(x) for x in tr_pauses]

    # -----------------------------------------------------------------------------

    # =============================================================================
    # Construct Sequence
    # =============================================================================
    timing_quant = 0
    for k in range(params.average):  # Averages
        for curr_concat in range(required_concats):
            for phase_steps in k_space_order_filing:  # instead of phase steps list of phase steps
                # for curr_slice in range(np.int32(params.sl_nb)):  # Slices
                for curr_slice in slice_list[curr_concat]:
                    seq = magn_prep_add(magn_prep_objects, seq, curr_slice,
                                        only_one='inv_prep')  # magn prep block
                    # Apply RF offsets
                    n_echo_temp = 0
                        
                    rf90.freq_offset = pulse_freq_offsets90[curr_slice]
                    rf180.freq_offset = pulse_freq_offsets180[curr_slice]
                    rf90.phase_offset = (rf90_phase - 2 * np.pi * rf90.freq_offset * calc_rf_center(rf90)[0])
                    rf180.phase_offset = (rf180_phase - 2 * np.pi * rf180.freq_offset * calc_rf_center(rf180)[0])
                    seq.add_block(rf90, gz90)

                    _ ,gz_spoil1_first, _, _ = refocusing_grad(params, scanner_parameters,
                                                               f_spex * gz90.area, flip180,
                                                               rf180_phase, t_refwd,
                                                               spoil_duration=t_spex, united=united)
                    seq.add_block(gz_spoil1_first, gx_pre)
                    
                    timing_quant += calc_duration(gx_pre)

                    for phase_step in phase_steps:
                        # print('phase step_' + str(phase_step))

                        seq.add_block(gz180, rf180)
                        timing_quant += gz180.shape_dur

                        if phase_step == 'skip' or phase_step == 'dummy':
                            gy_pre = make_trapezoid(channel='y', system=scanner_parameters,
                                                    area=0, duration=calc_duration(gz_spoil2),
                                                    rise_time=params.dG)
                        else:
                            gy_pre = make_trapezoid(channel='y', system=scanner_parameters,
                                                    area=k_steps_PE[phase_step], duration=calc_duration(gz_spoil2),
                                                    rise_time=params.dG)

                        # print(k_steps_PE[phase_step])
                        seq.add_block(gz_spoil2, gx_spoil1, gy_pre)
                        timing_quant += calc_duration(gy_pre)

                        if phase_step == 'skip' or phase_step == 'dummy':
                            # seq.add_block(gx)
                            timing_quant += gx.shape_dur
                        else:
                            seq.add_block(gx, adc)
                            timing_quant += gx.shape_dur / 2
                            timing_quant += gx.shape_dur / 2

                        n_echo_temp += 1
                        if phase_step == 'skip' or phase_step == 'dummy':
                            gy_pre = make_trapezoid(channel='y', system=scanner_parameters,

                                                    area=-0, duration=calc_duration(gz_spoil2),
                                                    rise_time=params.dG)
                        else:
                            gy_pre = make_trapezoid(channel='y', system=scanner_parameters,
                                                    area=-k_steps_PE[phase_step], duration=calc_duration(gz_spoil2),
                                                    rise_time=params.dG)
                        if n_echo_temp == np.int32(params.ETL):
                            seq.add_block(gz_cr, gx_spoil2, gy_pre)
                            timing_quant += calc_duration(gz_cr)
                        else:
                            seq.add_block(gz_spoil1, gx_spoil2, gy_pre)
                            timing_quant += gx_spoil2.shape_dur
                        
                            
                    seq.add_block(tr_fillers[curr_concat])
                    timing_quant += calc_duration(tr_fillers[curr_concat])


                    

    ok, error_report = seq.check_timing()
    if ok:
        print("Timing check passed successfully")
    else:
        print("Timing check failed. Error listing follows:")
        [print(e) for e in error_report]

        # ======
        # VISUALIZATION
        # ======
    return seq, k_space_save