physics_based_augmentation/seqgen/MRI_sequences/py2jemris/pulseq_jemris_simulator.py

# Same role as pulseq_bloch_simulator.py except (1) It uses JEMRIS (2) It is run as a function, not a script

# INPUTS: sequence type, geometry parameters, contrast parameters,
#         phantom type (pre-set or custom), coil type (pre-set or custom) (Tx / Rx),
#         k-space trajectory (if noncartesian; flattened version (Nro_total, 3))
import time
import os
from sim_jemris import sim_jemris
from recon_jemris import read_jemris_output, recon_jemris
from coil2xml import coil2xml
#from virtualscanner.utils import constants
from seq2xml import seq2xml
from pypulseq.Sequence.sequence import Sequence
import phantom as pht
import numpy as np
import xml.etree.ElementTree as ET
from pulseq_library import make_pulseq_se_oblique,make_pulseq_gre_oblique, make_pulseq_irse_oblique
from scipy.io import savemat, loadmat

def simulate_pulseq_jemris(seq_path, phantom_info, coil_fov,
                           tx='uniform', rx='uniform', # TODO add input that includes sequence info for
                                                       # TODO      dimensioning the RO points into kspace
                           tx_maps=None, rx_maps=None, sim_name=None, env_option="local"):
    """Runs simulation using an already-made .seq file

    Inputs
    ------
    seq_path : str
        Path to seq file
    phantom_info : dict
        Information used to create phantom; input to create_and_save_phantom()
    coil_fov : float
        Field-of-view of coil in mm
    tx : str, optional
        Type of tx coil; default is 'uniform'; the only other option is 'custom'
    rx : str, optional
        Type of rx coil; default is 'uniform'; the only other option is 'custom'
    tx_maps : list, optional
        List of np.ndarray (dtype='complex') maps for all tx channels
        Required for 'custom' type tx
    rx_maps : list, optional
        List of np.ndarray (dtype='complex') maps for all rx channels
        Required for 'custom' type rx
    sim_name : str, optional
        Used as folder name inside sim folder
        Default is None, in which case sim_name will be set to the current timestamp

    Returns
    -------
    sim_output :
        Delivers output from sim_jemris

    """
    if sim_name is None:
        sim_name = time.strftime("%Y%m%d%H%M%S")

    if env_option == 'local':
        target_path =  'sim\\' +  sim_name
    elif env_option == 'colab':
        target_path = 'sim\\' + sim_name

    # Make target folder
    dir_str = f'sim\\{sim_name}'
    if not os.path.isdir(dir_str):
        os.system(f'mkdir {dir_str}')

    # Convert .seq to .xml
    seq = Sequence()
    seq.read(seq_path)
    print(seq.get_block(1))
    seq_name = seq_path[seq_path.rfind('/')+1:seq_path.rfind('.seq')]
    seq2xml(seq, seq_name=seq_name, out_folder=str(target_path))

    # Make phantom and save as .h5 file
    pht_name = create_and_save_phantom(phantom_info, out_folder=target_path)


    # Make sure we have the tx/rx files
    tx_filename = tx + '.xml'
    rx_filename = rx + '.xml'

    # Save Tx as xml
    if tx == 'uniform':
        cp_command = f'copy {os.getcwd()}\\sim\\{tx}.xml {os.getcwd()}\\{str(target_path)}\\{tx}.xml'
        print(cp_command)
        a = os.system(cp_command)
        print(a)
    elif tx == 'custom' and tx_maps is not None:
        coil2xml(b1maps=tx_maps, fov=coil_fov, name='custom_tx', out_folder=target_path)
        tx_filename = 'custom_tx.xml'
    else:
        raise ValueError('Tx coil type not found')

    # save Rx as xml
    if rx == 'uniform':
        os.system(f'copy sim\\{rx}.xml sim\\{str(target_path)}')
    elif rx == 'custom' and rx_maps is not None:
        coil2xml(b1maps=rx_maps, fov=coil_fov, name='custom_rx', out_folder=target_path)
        rx_filename = 'custom_rx.xml'
    else:
        raise ValueError('Rx coil type not found')

    # Run simuluation in target folder
    list_sim_files = {'seq_xml': seq_name+'.xml', 'pht_h5': pht_name + '.h5', 'tx_xml': tx_filename,
                       'rx_xml': rx_filename}
    sim_output = sim_jemris(list_sim_files=list_sim_files, working_folder=target_path)

    return sim_output


# TODO
def create_and_save_phantom(phantom_info, out_folder):
    """Generates a phantom and saves it into desired folder as .h5 file for JEMRIS purposes

    Inputs
    ------
    phantom_info : dict
        Info of phantom to be constructed

        REQUIRED
        'fov' : float, field-of-view [meters]
        'N' : int, phantom matrix size (isotropic)
        'type' : str, 'spherical', 'cylindrical' or 'custom'
        'dim' : int, either 3 or 2; 3D or 2D phantom options
        'dir' : str, {'x', 'y', 'z'}; orientation of 2D phantom

        OPTIONAL (only required for 'custom' phantom type)
        'T1' : np.ndarray, T1 map matrix
        'T2' : np.ndarray, T2 map matrix
        'PD' : np.ndarray, PD map matrix
        'dr' : float, voxel size [meters] (isotropic)
        'dBmap' : optional even for 'custom' type. If not provided, dB is set to 0 everywhere.


    out_folder : str or pathlib Path object
        Path to directory where phantom will be saved

    Returns
    -------
    pht_type : str
        phantom_info['pht_type'] (returned for naming purposes)


    """
    out_folder = str(out_folder)

    FOV = phantom_info['fov']
    N = phantom_info['N']
    pht_type = phantom_info['type']
    pht_dim = phantom_info['dim']
    pht_dir = phantom_info['dir']
    pht_loc = phantom_info['loc']

    sim_phantom = 0

    if pht_type == 'spherical':
        print('Making spherical phantom')
        T1s = [1000*1e-3]
        T2s = [100*1e-3]
        PDs = [1]
        R = 0.8*FOV/2
        Rs = [R]
        if pht_dim == 3:
            sim_phantom = pht.makeSphericalPhantom(n=N, fov=FOV, T1s=T1s, T2s=T2s, PDs=PDs, radii=Rs)

        elif pht_dim == 2:
            sim_phantom = pht.makePlanarPhantom(n=N, fov=FOV, T1s=T1s, T2s=T2s, PDs=PDs, radii=Rs,
                                                dir=pht_dir, loc=0)
    elif pht_type == 'cylindrical':
        print("Making cylindrical phantom")
        sim_phantom = pht.makeCylindricalPhantom(dim=pht_dim, n=N, dir=pht_dir, loc=pht_loc)

    elif pht_type == 'custom':
        # Use a custom file!
        T1 = phantom_info['T1']
        T2 = phantom_info['T2']
        PD = phantom_info['PD']
        dr = phantom_info['dr']
        if 'dBmap' in phantom_info.keys():
            dBmap = phantom_info['dBmap']
        else:
            dBmap = 0

        sim_phantom = pht.Phantom(T1map=T1, T2map=T2, PDmap=PD, vsize=dr, dBmap=dBmap, loc=(0,0,0))

    else:
        raise ValueError("Phantom type non-existent!")

    # Save as h5
    sim_phantom.output_h5(out_folder, pht_type)

    return pht_type

if __name__ == '__main__':
    # Define the same phantom
    phantom_info = {'fov': 0.256, 'N': 32, 'type': 'cylindrical', 'dim':2, 'dir':'z'}
    sim_names = ['test0413_GRE', 'test0413_SE', 'test0413_IRSE']
    sps = ['gre_fov256mm_Nf15_Np15_TE50ms_TR200ms_FA90deg.seq',
           'se_fov256mm_Nf15_Np15_TE50ms_TR200ms_FA90deg.seq',
           'irse_fov256mm_Nf15_Np15_TI20ms_TE50ms_TR200ms_FA90deg.seq']
    # make_pulseq_irse_oblique(fov=0.256,n=15, thk=0.005, tr=0.2, te=0.05, ti=0.02, fa=90,
    #                          enc='xyz', slice_locs=[0], write=True)
    # make_pulseq_gre_oblique(fov=0.256,n=15, thk=0.005, tr=0.2, te=0.05, fa=90,
    #                          enc='xyz', slice_locs=[0], write=True)
    # make_pulseq_se_oblique(fov=0.256,n=15, thk=0.005, tr=0.2, te=0.05, fa=90,
    #                          enc='xyz', slice_locs=[0], write=True)
    simulate_pulseq_jemris(seq_path=sps[0], phantom_info=phantom_info, sim_name=sim_names[0],
                               coil_fov=0.256)
    kk, im, images = recon_jemris(file='sim/' + sim_names[0] + '/signals.h5', dims=[15,15])
    savemat('sim/'+sim_names[0]+'/output.mat', {'images': images, 'kspace': kk, 'imspace': im})