physics_based_augmentation/seqgen/utilities/TSE_k_space_fill.py

import numpy as np
import random
from copy import deepcopy


def TSE_k_space_fill(n_ex, ETL, k_steps, TE_eff_number, N_center, order, dummy, N_dummy):
    # function defines phase encoding steps for k space filling in liner order
    # with shifting according to the TE effective number

    if dummy:
        k_space_dummy_temp = ['dummy']*ETL
        k_space_dummy = []
        for i in range(N_dummy):
            k_space_dummy.append(k_space_dummy_temp)
    k_space_list_with_zero = []
    for i in range(ETL):
        #k_space_list_with_zero.append(int((ETL - 1) * n_ex - i * n_ex))
        k_space_list_with_zero.append(i * n_ex)

    if n_ex > 1:
        k_space_order_filing_temp = [k_space_list_with_zero]
        for i in range(n_ex - 1):
            k_space_list_temp = []
            for k in k_space_list_with_zero:
                k_space_list_temp.append(k + i + 1)
            k_space_order_filing_temp.append(k_space_list_temp)
    else:
        k_space_order_filing_temp = [k_space_list_with_zero]


    kk = 0
    for k_space_list_with_central_old in k_space_order_filing_temp:
        temporal = abs(k_space_list_with_central_old - N_center)
        center_index = np.argmin(temporal)

        central_line_position_old = center_index
        num_left_old = central_line_position_old
        num_right_old = len(k_space_list_with_central_old) - num_left_old - 1

        for i in range(len(k_space_list_with_central_old)):
            if k_space_list_with_central_old[i] > (k_steps - 1):
                k_space_list_with_central_old[i] = 'skip'

        k_space_list_with_central_new = ['skip'] * len(k_space_list_with_central_old)
        central_line_position_new = TE_eff_number - 1
        num_left_new = central_line_position_new
        num_right_new = len(k_space_list_with_central_new) - num_left_new - 1

        if num_left_new < num_right_new:
            min_new = num_left_new
            max_new = num_right_new
        else:
            min_new = num_right_new
            max_new = num_left_new

        k_space_list_with_central_new[central_line_position_new] = k_space_list_with_central_old[
            central_line_position_old]
        k_space_list_with_central_old[central_line_position_old] = 'passed'

        if min_new == 0: # effective time is first one
            i = 0
            k_space_list_with_central_old_abs = [0] * len(k_space_list_with_central_old)
            for ii in range(len(k_space_list_with_central_old)):
                if k_space_list_with_central_old[ii] != 'skip' and k_space_list_with_central_old[ii] != 'passed':
                    k_space_list_with_central_old_abs[ii] = abs(k_space_list_with_central_old[ii] - N_center)
                    #type_of_elements = type(k_space_list_with_central_old_abs[0])
                else:
                    k_space_list_with_central_old_abs[ii] = k_space_list_with_central_old[ii]
                    #type_of_elements = type(k_space_list_with_central_old_abs[0])
            type_of_string = type('passed')

            if num_left_new < num_right_new:
                flag = True
                j = central_line_position_new + i + 1
                while flag:
                    a = min(filter(lambda i: not isinstance(i, type_of_string), k_space_list_with_central_old_abs))
                    min_index = k_space_list_with_central_old_abs.index(a)
                    k_space_list_with_central_new[j] = k_space_list_with_central_old[min_index]
                    k_space_list_with_central_old_abs[min_index] = 'passed'
                    k_space_list_with_central_old[min_index] = 'passed'
                    j += 1
                    flag = any(not isinstance(y, (type_of_string)) for y in k_space_list_with_central_old)
                k_space_order_filing_temp[kk] = k_space_list_with_central_new
                kk += 1
            else:
                flag = True
                j = central_line_position_new - i - 1
                while flag:
                    a = min(filter(lambda i: not isinstance(i, type_of_string), k_space_list_with_central_old_abs))
                    min_index = k_space_list_with_central_old_abs.index(a)
                    k_space_list_with_central_new[j] = k_space_list_with_central_old[min_index]
                    k_space_list_with_central_old_abs[min_index] = 'passed'
                    k_space_list_with_central_old[min_index] = 'passed'
                    j -= 1
                    flag = any(not isinstance(y, (type_of_string)) for y in k_space_list_with_central_old)
                k_space_order_filing_temp[kk] = k_space_list_with_central_new
                kk += 1

        else:
            for i in range(1, min_new + 1):
                k_space_list_with_central_new[central_line_position_new + i] = k_space_list_with_central_old[
                    central_line_position_old + i]
                k_space_list_with_central_old[central_line_position_old + i] = 'passed'
                k_space_list_with_central_new[central_line_position_new - i] = k_space_list_with_central_old[
                    central_line_position_old - i]
                k_space_list_with_central_old[central_line_position_old - i] = 'passed'

            k_space_list_with_central_old_abs = [0] * len(k_space_list_with_central_old)
            for ii in range(len(k_space_list_with_central_old)):
                if k_space_list_with_central_old[ii] != 'skip' and k_space_list_with_central_old[ii] != 'passed':
                    k_space_list_with_central_old_abs[ii] = abs(k_space_list_with_central_old[ii] - N_center)
                    #type_of_elements = type(k_space_list_with_central_old_abs[0])
                else:
                    k_space_list_with_central_old_abs[ii] = k_space_list_with_central_old[ii]
                    #type_of_elements = type(k_space_list_with_central_old_abs[0])
            type_of_string = type('passed')

            if num_left_new < num_right_new:
                flag = True
                j = central_line_position_new + i + 1
                while flag:
                    a = min(filter(lambda i: not isinstance(i, type_of_string), k_space_list_with_central_old_abs))
                    min_index = k_space_list_with_central_old_abs.index(a)
                    k_space_list_with_central_new[j] = k_space_list_with_central_old[min_index]
                    k_space_list_with_central_old_abs[min_index] = 'passed'
                    k_space_list_with_central_old[min_index] = 'passed'
                    j += 1
                    flag = any(not isinstance(y, (type_of_string)) for y in k_space_list_with_central_old)
                k_space_order_filing_temp[kk] = k_space_list_with_central_new
                kk += 1
            else:
                flag = True
                j = central_line_position_new - i - 1
                while flag:
                    a = min(filter(lambda i: not isinstance(i, type_of_string), k_space_list_with_central_old_abs))
                    min_index = k_space_list_with_central_old_abs.index(a)
                    k_space_list_with_central_new[j] = k_space_list_with_central_old[min_index]
                    k_space_list_with_central_old_abs[min_index] = 'passed'
                    k_space_list_with_central_old[min_index] = 'passed'
                    j -= 1
                    flag = any(not isinstance(y, (type_of_string)) for y in k_space_list_with_central_old)
                k_space_order_filing_temp[kk] = k_space_list_with_central_new
                kk += 1
            #k_space_order_filing = k_space_order_filing.append(k_space_list_with_central_new)


    """
        k_space_order_filing = []
        n_ex_shift = np.int64(n_ex/2)-1
        a = range(0-n_ex_shift, n_ex-n_ex_shift)
        for i in a:
            k_space_list_temp = []
            for kk in k_space_list_with_central_new:
                value = kk + i
                if value < 0:
                    value = value + k_steps
                if value > k_steps - 1:
                    value = 'skip'
                k_space_list_temp.append(value)
            k_space_order_filing.append(k_space_list_temp)
    """

    if dummy:
        for list in k_space_order_filing_temp:
            k_space_dummy.append(list)
    else:
        k_space_dummy = k_space_order_filing_temp
    return k_space_dummy

def HASTE_k_space_fill(n_ex, ETL, k_steps, TE_eff_number, N_center, order, dummy, N_dummy, params):
    # function defines phase encoding steps for k space filling in liner order
    # with shifting according to the TE effective number

    if dummy:
        k_space_dummy_temp = ['dummy']*k_steps
        k_space_dummy = []
        for i in range(N_dummy):
            k_space_dummy.append(k_space_dummy_temp)
    k_space_list = [*range(0, k_steps)]
    eff_shift = TE_eff_number - N_center - 1

    if params['part_fourier_factor_phase'] == 1.0:
        k_space_list_reordered = np.roll(k_space_list, shift=eff_shift)
        if TE_eff_number < int(k_space_list_reordered.size / 2):
            first_descrim = (TE_eff_number * 2)
            temp_list_second = k_space_list_reordered[first_descrim:]
            temp_list_first = k_space_list_reordered[:first_descrim]

            temp_list_reordered = np.empty([int(temp_list_second.size)], dtype=int)
            for i in range(int(temp_list_reordered.size / 2)):
                print(i)
                temp_list_reordered[2 * i + 1] = temp_list_second[-1 - i]
                temp_list_reordered[2 * i] = temp_list_second[i]
            k_space_list_reordered = np.concatenate((temp_list_first, temp_list_reordered), axis=None)

    else:
        k_space_list_reordered = np.roll(k_space_list, shift=eff_shift)
        zero_index = k_space_list_reordered.tolist().index(0)
        temp_list = k_space_list_reordered[zero_index:]
        temp_list = temp_list[::-1]
        temp_max_list = k_space_list_reordered[:zero_index]
        k_space_list_reordered = np.concatenate((temp_max_list, temp_list), axis=None)


    if order == "non_linear":
        pass

    elif order == "random":
        q1 = params['Np'] // 3
        q2 = (params['Np'] - q1) // 2
        q3 = params['Np'] - q1 - q2
        # TODO: random sampling
        #list_q1 = k_space_list_reordered[:q1]
        #if q2 > len(k_space_list_reordered):
        #    list_q2 = k_space_list_reordered[q1:q1+q2]
        #    list_q3 = k_space_list_reordered[q2:]
        #else:
        #    list_q2 = k_space_list_reordered[q1:]
        #    list_q3 = []

        #list_q1, list_q3 = redistribute_randomly_same_size(list_q1, list_q3)
        #k_space_list_reordered = list_q1 + list_q2 + list_q3

    else:
        print("That's not a valid order.")  # Default case

    if dummy:
        k_space_dummy.append(k_space_list_reordered.tolist())
    else:
        k_space_dummy = [k_space_list_reordered.tolist()]
    return k_space_dummy


def redistribute_randomly_same_size(list1, list2):
    combined = list1 + list2

    random.shuffle(combined)

    size1 = len(list1)
    size2 = len(list2)

    new_list1 = combined[:size1]
    new_list2 = combined[size1:size1 + size2]

    return new_list1, new_list2


def circle_2D_filling_TSE(Nx, Ny, Nx_center, Ny_center, target_pixels, N_points, previous_pixels):
    # Начальный радиус (по формуле площади круга)
    radius = np.sqrt(target_pixels / np.pi)

    num_iteration = int(target_pixels / N_points)

    while True:
        # Список пикселей внутри круга
        selected_pixels = []

        for x in range(Nx):
            for y in range(Ny):
                if (x - Nx_center) ** 2 + (y - Ny_center) ** 2 <= radius ** 2:
                    selected_pixels.append((x, y))
        # print(len(selected_pixels))

        if (len(selected_pixels) >= target_pixels):
            break

        # Корректируем радиус: уменьшаем, если пикселей много, увеличиваем, если мало
        if len(selected_pixels) < target_pixels:
            radius += 0.1  # Увеличиваем радиус

    selected_pixels_final = [x for x in selected_pixels if x not in previous_pixels]
    selected_pixels_final = selected_pixels_final[:N_points]
    # print(len(selected_pixels_final))
    '''
    # Создаем изображение
    image = np.zeros((Nx, Ny))
    for x, y in selected_pixels_final:
        image[y, x] = 1  # Заполняем белым цветом

    # Визуализация
    plt.imshow(image, cmap="gray", origin="lower")
    plt.scatter(Nx_center, Ny_center, color="red", marker="x")  # Отмечаем центр
    plt.grid(visible=True, color="red", linestyle="--", linewidth=0.5)
    plt.title(f"Круг из {len(selected_pixels_final)} пикселей")
    plt.show()
    '''
    return selected_pixels_final


def circle_2D_filling_TSE_last_acq(Nx, Ny, previous_pixels, N_points):
    all_pixels = []
    for x in range(Nx):
        for y in range(Ny):
            all_pixels.append((x, y))

    selected_pixels_final = [x for x in all_pixels if x not in previous_pixels]

    '''
    image = np.zeros((Nx, Ny))
    for x, y in selected_pixels_final:
        image[y, x] = 1  # Заполняем белым цветом

    # Визуализация
    plt.imshow(image, cmap="gray", origin="lower")
    plt.scatter(Nx_center, Ny_center, color="red", marker="x")  # Отмечаем центр
    plt.grid(visible=True, color="red", linestyle="--", linewidth=0.5)
    plt.title(f"Круг из {len(selected_pixels_final)} пикселей")
    plt.show()
    '''

    if len(selected_pixels_final) == 0:
        return None

    elif (selected_pixels_final) != N_points:
        for i in range(N_points - len(selected_pixels_final)):
            selected_pixels_final.append(('skip', 'skip'))
            # print(len(selected_pixels_final))
        return selected_pixels_final

    elif (selected_pixels_final) == N_points:
        return selected_pixels_final


def effective_TE_reorder(num_te_eff, num_echo):
    if num_te_eff <= num_echo / 2:
        a_right = list(range(1, (num_te_eff - 1) * 2, 2))
        a_left = list(range((num_te_eff - 1) * 2, 1, -2))
        a = a_left + [0] + a_right
        a_end = list(range(len(a), num_echo))
        a_final = a + a_end
    elif num_te_eff > num_echo / 2:
        k = num_echo - num_te_eff + 1
        a_right = list(range(1, (k - 1) * 2, 2))
        a_left = list(range((k - 1) * 2, 1, -2))
        a = a_left + [0] + a_right
        a_end = list(range(num_echo - 1, len(a) - 1, -1))
        a_final = a_end + a
    return a_final

def k_space_filling_3D_TSE(Ny, Nz, Ny_center, Nz_center, num_echo, num_te_eff, dummy, N_dummy):
    num_echo #ETL


    #TODO: add
    if Ny * Nz % num_echo == 0:
        N_points = int(Ny * Nz / num_echo)
    else:
        N_points = int(Ny * Nz / num_echo ) + 1  # number TE
    previous_pixels = []
    #Ny = 32
    #Nz = 32
    #Ny_center = 25
    #Nz_center = 25

    num_sat = int(Ny * Nz / N_points)  # number os saturations
    print(num_sat)


    if dummy:
        k_space_dummy_temp =[]
        for i in range(num_echo):
            k_space_dummy_temp.append(('dummy', 'dummy'))
        k_space_dummy = []
        for i in range(N_dummy):
            k_space_dummy.append(k_space_dummy_temp)


    k_space_fill = []
    for i in range(num_sat):
        target_num_points = (i + 1) * N_points
        selected_pixels = circle_2D_filling_TSE(Ny, Nz, Ny_center, Nz_center, target_num_points, N_points, previous_pixels)
        k_space_fill.append(selected_pixels)
        previous_pixels = previous_pixels + selected_pixels

    last_k_space_volume = circle_2D_filling_TSE_last_acq(Ny, Nz, previous_pixels, N_points)
    if last_k_space_volume != None:
        k_space_fill.append(last_k_space_volume)

    if num_sat + 1 == num_echo:
        pass
    else:
        k_space_skip_temp = []
        for i in range( N_points):
            k_space_skip_temp.append(('skip', 'skip'))
        for i in range(num_echo - num_sat - 1):
            k_space_fill.append(k_space_skip_temp)

    eff_TE_reordered = effective_TE_reorder(num_te_eff, num_echo)

    k_space_fill_final = []
    for ii in range(len(k_space_fill[0])):
        one_line = []
        for jj in range(len(k_space_fill)):
            one_line.append(k_space_fill[eff_TE_reordered[jj]][ii])
        k_space_fill_final.append(one_line)
    if dummy:
        k_space_fill_final = k_space_dummy + k_space_fill_final

    return k_space_fill_final