对图像提取ROI区域以及将输出结果扩展为原图像大小的代码实现

这篇文章总共分为准备工作、ROI提取、图像扩展三部分。

1. 准备工作

这里我使用的是SimpleITK模块对图像进行操作。
需要导入的模块如下：

import SimpleITK as sitk
import os
import matplotlib.pyplot as plt
import numpy as np

然后，需要设置根目录root_dir ，待截取文件所在文件夹image_file ，我是根据我自己的分割结果（保存在seg_result_file 内）结合金标准（保存在label_reference_file）作为参考共同确定的ROI区域的大小，因为我自己的分割结果在某些层可能没有分割，所以需要结合金标准来确保自己的ROI包含整个待分割区域。

root_dir = r'根目录'
image_file = os.path.join(root_dir, "imagesTr")
label_reference_file = os.path.join(root_dir, "labelsTr")
seg_result_file = os.path.join(root_dir, "prediction")

然后开始读取图片获取bounding box ,在labelFilter.GetBoundingBox(1)这一指令中，参数1代表标签值，意思是获取标签值为1的区域所在的bounding box。因为我要结合自己的分割结果和金标准来获取ROI所以我分别获取了两个bounding box。获取的bounding box是一个元组(x, y, z, w, h, d)，其中前三个元素表示的是bounding box的角点，后三个元素分别表示bounding box的宽度、高度和深度。
我将两个bounding box的z轴上的起点（z）和其深度做了改变。选择两个bounding box中z的最小值作为最终bounding box的起点，通过z+d-1可以算出终点值，将终点最大值作为最终bounding box的终点，使用这个终点减去起点再加上1即为最终bounding box的深度。
因为我们希望图像包含多一点背景信息，所以在bounding box的基础上添加一个pad_size.

'''读取图片'''
prostate_img = sitk.Cast(sitk.ReadImage(os.path.join(prostate_file, prostate_name)), sitk.sitkUInt8)
label_reference_image = sitk.Cast(sitk.ReadImage(os.path.join(label_reference_file, label_reference_name)), sitk.sitkUInt8)
seg_result_img = sitk.Cast(sitk.ReadImage(os.path.join(seg_result_file, seg_result_name)), sitk.sitkUInt8)
seg_result_img.SetOrigin(label_img.GetOrigin())
'''获取两个bounding box'''    
labelFilter = sitk.LabelShapeStatisticsImageFilter()
labelFilter.Execute(label_reference_image)
label_reference_bbox = labelFilter.GetBoundingBox(1)
labelFilter.Execute(seg_result_img)
seg_result_bbox = labelFilter.GetBoundingBox(1)
'''计算最终bounding box'''   
z_start = min(label_reference_bbox[2], seg_result_bbox [2])
z_end = max(label_reference_bbox[2] + label_reference_bbox[-1] - 1, seg_result_bbox [2] + seg_result_bbox [-1] - 1)
depth = z_end - z_start + 1 - seg_result_bbox [-1]
bbox _list = list(seg_result_bbox ) # 元组不能改变内部元素，所以先将其转变为列表，改变元素值之后再转换为元组。
bbox _list[2] = z_start
bbox = tuple(bbox _list)

pad_size = 20

2. ROI提取

使用sitk.RegionOfInterest（img， roi_size, roi_index）进行ROI提取。

def crop_image(image, pad_size, depth, bbox, save_path, save_name):
    roi_index = [x - y for x, y in zip(bbox[:3], [pad_size, pad_size, 0])] # 角点坐标
    roi_size = [x + y for x, y in zip(bbox[3:], [pad_size * 2, pad_size * 2, depth])] # ROI的宽度、高度和深度
    # print("roi_size: ", roi_size)
    image_crop = sitk.RegionOfInterest(image, roi_size, roi_index)
    cropped_image_save_path = save_path
    if not os.path.exists(cropped_image_save_path):
        os.mkdir(cropped_image_save_path)
    sitk.WriteImage(image_crop, os.path.join(cropped_image_save_path, save_name))

3. 图像扩展

本来考虑到使用MONAI中的pading模块对图像进行填充，但是对称填充并不能适应我裁剪过得图像，因此我的思路是：先创建一个和原图像大小相同的数组，然后根据我之前计算ROI区域的步骤进行反推，将我的ROI图像中的像素值一个一个对应填进我创建的数组内。

def expand_image(image, label_reference, pad_size, bbox, save_path, save_name):
    white_arr = np.zeros(label_reference.GetSize()).T
    image_arr = sitk.GetArrayFromImage(image)
    image_shape = image_arr.shape
    origin = np.array([bbox[2], bbox[1]-pad_size, bbox[0]-pad_size]) # 将之前准备工作中获取的最终bbox的前三个元素作为参考原点。
    # 这里要注意，simpleITK和ndarray相互转换时其尺寸会相互颠倒，
    # simpleITK读取的图片尺寸格式为（w，h，d），转换成ndarray数组的尺寸会变为（d，h，w）
    print(origin)
    for z in range(image_shape[0]):
        for y in range(image_shape[1]):
            for x in range(image_shape[-1]):
                coordinate = np.array([z, y, x])
                mixed_coordinate = origin.__add__(coordinate)
                mixed_coordinate = tuple(mixed_coordinate)
                white_arr[mixed_coordinate] = image_arr[z, y, x]
    expanded_image = sitk.GetImageFromArray(white_arr)
    expanded_image.CopyInformation(label_reference)
    expanded_image_save_path = save_path
    if not os.path.exists(expanded_image_save_path):
        os.mkdir(expanded_image_save_path)
    sitk.WriteImage(expanded_image, os.path.join(expanded_image_save_path, save_name))

以上是我对自己代码的记录，为了方便阅读和理解我对自己原来的代码做了一些改动，主要是提供思路。
如果对你有所帮助我不胜荣幸，如果大家有宝贵意见或发现我存在的错误，请向我提出，互相学习，共同进步。