lf_mri_platform/services/reconstructor/reco.py

import numpy as np
import os
import matplotlib
matplotlib.use('Agg')  # headless backend — must be set before pyplot import
import matplotlib.pyplot as plt
import sys


# os.chdir('C:/')
# sys.path.append("C:/Users/iuliia/recoUI/serv")
def reconstruction(matrix,k_space_min, k_space_max, nCoil, nContrast, nSlice, phase_shift=0, return_image=True):
    
    for i in range(matrix.shape[1]):
        phase_factors = np.exp(-1j * phase_shift) ** i
        matrix[:, i] = matrix[:, i] * phase_factors

    fft_data = np.fft.ifft2(matrix)
    fft_shifted_data = np.fft.ifftshift(fft_data)
    amplitude = np.abs(fft_shifted_data)

    plt.imshow(amplitude, cmap='gray', vmin=amplitude.min(), vmax=amplitude.max())
    plt.colorbar()
    plt.savefig(
        os.path.join('data', f'im_contrast_{nCoil}_{nContrast}_k_{nSlice}.png'))
    plt.close()

    log_matrix = np.log(np.abs(matrix))
    print(k_space_min, " ", k_space_max)
    plt.imshow(log_matrix, cmap='gray', vmin=np.log(np.abs(k_space_min)), vmax=np.log(np.abs(k_space_max)))
    plt.colorbar()
    plt.savefig(
        os.path.join('data', f'k_contrast_{nCoil}_{nContrast}_k_{nSlice}.png'))
    plt.close()

    if return_image:
        return amplitude


def reco_3d(matrix, maximum_k_space, minimum_k_space, nCoil, nContrast, phase_shift=0, return_image=True):
    '''
        :param matrix: [x,y,z]
        :param nContrast: количество контрастов
        :param phase_shift: угол
        :return: 
    '''

    for i in range(matrix.shape[1]):
        phase_factors = np.exp(-1j * phase_shift) ** i
        matrix[:, i, :] = matrix[:, i, :] * phase_factors

    fft_data = np.fft.ifftn(matrix)
    fft_shifted_data = np.fft.ifftshift(fft_data)
    amplitude = np.abs(fft_shifted_data)
    # print("TTTT")

    for slice in range(amplitude.shape[2]):
        plt.imshow(amplitude[:, :, slice], cmap='gray', vmin=amplitude.min(),
                   vmax=amplitude.max())
        plt.colorbar()
        plt.savefig(
            os.path.join('data', f'im_contrast_{nCoil}_{nContrast}_k_{slice}.png'))
        plt.close()

        # print("slice")

        log_matrix = np.log(np.abs(matrix))
        plt.imshow(log_matrix[:, :, slice], cmap='gray', vmin=np.log(np.abs(minimum_k_space)),
                   vmax=np.log(np.abs(maximum_k_space)))
        plt.savefig(
            os.path.join('data', f'k_contrast_{nCoil}_{nContrast}_k_{slice}.png'))
        plt.close()

    if return_image:
        return amplitude


def save_ffft(data, spoil=0):
    '''
        :param data: [x, y, z, contrast, coil]
        :param spoil: градус сдвига
    '''
    print(data.shape)

    k_space_max = np.max(data)
    k_space_min = np.min(data)

    num_slice = data.shape[2]  # Размер третьего измерения slice
    num_contrast = data.shape[3]  # Размер четвертого измерения contrast
    num_coils = data.shape[4]  # Размер по оси coil

    all_img = np.zeros((num_contrast, num_slice, data.shape[0], data.shape[1]), dtype=np.float32)

    for number_matrix_k_space in range(num_slice):
        for contrast_nb in range(num_contrast):
            img_sum = None

            for coil_nb in range(num_coils):
                coil_data = data[:, :, number_matrix_k_space, contrast_nb, coil_nb]
                img_coil = reconstruction(coil_data, k_space_min, k_space_max,
                                          coil_nb, contrast_nb, number_matrix_k_space,
                                          spoil,
                                          return_image=True )

                img_abs_two = img_coil ** 2

                if img_sum is None:
                    img_sum = img_abs_two
                else:
                    img_sum += img_abs_two

            all_img[contrast_nb, number_matrix_k_space] \
                = img_sum.copy()  # тут сохраняется ИЗОБРАЖЕНИЕ, квадрат модуля для ВСЕХ катушек данного констраста и слоя

    global_max = np.max(all_img)
    global_min = np.min(all_img)

    save_img(all_img, global_max, global_min)


def save_ffft_3d(data, phase=0):
    '''
        :param data:  [x, y, z, contrast, coil]
        :param phase: spoil
    '''
    num_contrast = data.shape[3]  # Размер четвертого измерения
    num_coils = data.shape[4]
    num_slice = data.shape[2]

    k_space_max = np.max(data)
    k_space_min = np.min(data)

    all_img = np.zeros((num_contrast, num_slice, data.shape[0], data.shape[1]), dtype=np.float32)

    # print("phase: ", phase)
    # for coil_nb in range(num_coils):
    for contrast_nb in range(num_contrast):
        img_sum = None

        for coil_nb in range(num_coils):
            coil_data = data[:, :, :, contrast_nb, coil_nb]
            img_coil = reco_3d(coil_data, k_space_max, k_space_min,
                               coil_nb, contrast_nb,
                               phase,
                               return_image=True)
            
            print(img_coil.shape)

            img_abs_two = img_coil ** 2

            if img_sum is None:
                img_sum = img_abs_two
            else:
                img_sum += img_abs_two


        img_sum = np.transpose(img_sum, (2, 0, 1))

        all_img[contrast_nb] \
            = img_sum.copy()  # тут сохраняется ИЗОБРАЖЕНИЕ, квадрат модуля для ВСЕХ катушек данного констраста
        
        print(contrast_nb)
        print(all_img.shape)

    # all_img = np.transpose(all_img, (1, 2, 3, 0))

    global_max = np.max(all_img)
    global_min = np.min(all_img)
    
    print(global_max, global_min)
    print(all_img)

    save_img(all_img, global_max, global_min)


def save_img(data, maximum, minimum):
    '''
        :param data: all_img[contrast_nb, number_matrix_k_space] = img_sum.[x,y] = [contrast, slice, x, y] 
        :return: 
    '''
    
    contrasts = data.shape[0]
    slices = data.shape[1]

    for contrast in range(contrasts):
        for slice in range(slices):
            plt.imshow(data[contrast, slice], cmap='gray', vmin=minimum, vmax=maximum)
            plt.colorbar()
            plt.savefig(
                os.path.join('data', f'im_{contrast}_k_{slice}.png'))
            plt.close()