处理CT图像

import os

import nibabel as nib
import numpy as np

"""
通过分割图像截取出来包裹分割图像的最大面积
"""


def vol(data):
    """
    获得每一层面的框并保存
    """
    seg_vol = []
    for z in range(data.shape[2]):
        layer = data[:, :, z]
        row, col = np.nonzero(layer)
        if len(row) > 0:
            min_row = np.min(row)
            max_row = np.max(row)
            min_col = np.min(col)
            max_col = np.max(col)
            seg_vol.append((z, min_row, min_col, max_row, max_col))

    mi_z = min(seg_vol, key=lambda x: x[0])[0]
    ma_z = max(seg_vol, key=lambda x: x[0])[0]
    mi_row = min(seg_vol, key=lambda x: x[1])[1]
    ma_row = max(seg_vol, key=lambda x: x[3])[3]
    mi_col = min(seg_vol, key=lambda x: x[2])[2]
    ma_col = max(seg_vol, key=lambda x: x[4])[4]

    return mi_z, ma_z, mi_row, ma_row, mi_col, ma_col


def open_file(input_file_path):
    # input_file_path = os.path.join(input_folder, filename)  # 构建完整的输入文件路径
    nii_file = nib.load(input_file_path)  # 打开NIfTI文件
    original_data = nii_file.get_fdata()  # 获取原始图像数据和头部信息
    print("图像数据形状为{}".format(original_data.shape))
    return original_data, nii_file


def resampe_data(data, nii_file, v_x, v_y, v_z, s_x, s_y, s_z):
    header = nii_file.header  # 重采样
    header['pixdim'][1:4] = [v_x, v_y, v_z]  # 更新头部信息中的像素尺寸信息
    header['qoffset_x'] = s_x
    header['qoffset_y'] = s_y
    header['qoffset_z'] = s_z
    resampled_img = nib.Nifti1Image(data, affine=nii_file.affine, header=header)  # 创建一个新的NIfTI图像对象
    return resampled_img


def save_data2nii(output_folder, outname, data):
    print("处理后的数据形状为{}".format(data.shape))
    output_file_path = os.path.join(output_folder, outname)  # 构建输出文件路径
    nib.save(data, output_file_path)


def window_transform(ct_array, windowWidth, windowCenter, normal=True):
    """
   return: trucated image according to window center and window width
   and normalized to [0,1]
   """
    minWindow = float(windowCenter) - 0.5 * float(windowWidth)
    newimg = (ct_array - minWindow) / float(windowWidth)
    newimg[newimg < 0] = 0
    newimg[newimg > 1] = 1
    min_val = newimg.min()  # 归一化像素值到0到1之间
    max_val = newimg.max()
    resampled_data = (newimg - min_val) / (max_val - min_val)
    if not normal:
        newimg = (newimg * 255).astype('uint8')
    return resampled_data


def read_nii_files(folder_path):
    nii_files = []
    for file in os.listdir(folder_path):
        if file.endswith('.nii.gz'):
            nii_files.append(os.path.join(folder_path, file))
    return nii_files


def pair_files(folder_path, img_prefix='img', seg_prefix='seg'):
    img_files = {}
    seg_files = {}

    # 遍历文件夹中的文件
    for file_name in os.listdir(folder_path):
        file_path = os.path.join(folder_path, file_name)

        # 检查文件是图像还是分割文件
        if file_name.startswith(img_prefix) and file_name.endswith('.nii.gz'):
            identifier = file_name[len(img_prefix):len(img_prefix) + 4]  # 提取编号部分，假设是4位数字
            img_files[identifier] = file_path
        elif file_name.startswith(seg_prefix) and file_name.endswith('.nii.gz'):
            identifier = file_name[len(seg_prefix):len(seg_prefix) + 4]  # 提取编号部分，假设是4位数字
            seg_files[identifier] = file_path

    # 寻找匹配的文件对
    matched_pairs = {}
    for identifier, img_path in img_files.items():
        seg_path = seg_files.get(identifier)
        if seg_path:
            matched_pairs[identifier] = {'img': img_path, 'seg': seg_path}

    return matched_pairs


if __name__ == '__main__':
    input_folder = r'F:\imaging registration\datasets\MRCT\L2R_Task1_CT\L2R_Task1_CT'
    output_folder = r'F:\imaging registration\datasets\MRCT\L2R_Task1_CT\L2R_Task1_CT_vol'
    log_file = r'F:\imaging registration\datasets\MRCT\L2R_Task1_CT\L2R_Task1_CT_log.txt'
    # img_name = 'img0001_bcv_CT.nii.gz'
    # seg_name = 'seg0001_bcv_CT.nii.gz'
    img_outname = 'imgvol_img0001_bcv_CT.nii.gz'
    seg_outname = 'segvol_seg0001_bcv_CT.nii.gz'
    voxel_x, voxel_y, voxel_z = [2, 2, 2]  # Spacing
    set_x, set_y, set_z = [0, 0, 0]  # Initial point coordinates
    window_width, window_center = [400, 40]

    matched_pairs = pair_files(input_folder)
    with open(log_file, 'a+') as f:
        f.seek(0)  # 将文件指针移动到文件开头
        if os.stat(log_file).st_size == 0:
            f.write("日志文件\n\n")
        for identifier, pair_info in matched_pairs.items():
            f.write(f"匹配的文件对 ({identifier}):\n")
            f.write("Image File: " + pair_info['img'] + '\n')
            f.write("Segmentation File: " + pair_info['seg'] + '\n')

            img_data, img_nii = open_file(pair_info['img'])
            seg_data, seg_nii = open_file(pair_info['seg'])

            min_z, max_z, min_row, max_row, min_col, max_col = vol(seg_data)
            img_data = window_transform(img_data, window_width, window_center)
            precessed_img_data = img_data[min_row - 2:max_row + 2, min_col - 2:max_col + 2, min_z - 2:max_z + 2]
            precessed_seg_data = seg_data[min_row - 2:max_row + 2, min_col - 2:max_col + 2, min_z - 2:max_z + 2]
            precessed_img_data = resampe_data(precessed_img_data, img_nii, voxel_x, voxel_y, voxel_z, set_x, set_y, set_z)
            precessed_seg_data = resampe_data(precessed_seg_data, seg_nii, voxel_x, voxel_y, voxel_z, set_x, set_y, set_z)
            f.write("处理后的数据形状{}".format(precessed_seg_data.shape)+'\n')
            save_data2nii(output_folder, f'img_vol_{identifier}_bcv_CT.nii.gz', precessed_img_data)
            save_data2nii(output_folder, f'seg_vol_{identifier}_bcv_CT.nii.gz', precessed_seg_data)
            f.write("\n")