医疗影像读取与写入

参考：
SimpleITK处理dicom数据 ；
【python】SimpleITK 和 Nibabel 读取医学图像 nii 数据（2D显示）

使用须知

• SimpleITK读取的图像数据sitk.ReadImage()的坐标顺序为C*H*W，即从多少张切片到单张切片的宽和高；而据SimpleITK Image获取的origin和spacing的坐标顺序则是H*W*C。

• nib.load()读取文件，会将图像向左旋转90°，坐标顺序为W*H*C（不太推荐）

SimpleITK

• ReadImage

import SimpleITK as sitk
import skimage.io as io

def read_img_sitk(filepath):
    image = sitk.ReadImage(filepath)  # <class 'SimpleITK.SimpleITK.Image'> # 支持dcm\nrrd\nii
    image_array = sitk.GetArrayFromImage(image)  # z,y,x
    return image_array

def show_img_sitk(data, slice_num=None):
    if slice_num:
        io.imshow(data[slice_num, :, :], cmap='gray')
        io.show()
    else:
        for i in range(data.shape[0]):
            io.imshow(data[i, :, :], cmap='gray')
            print(i)
            io.show()

path = "origin.nii"
data_arr = read_img_sitk(path)
show_img_sitk(data_arr, slice_num=45)

PET slice

• WriteImage

image = sitk.GetImageFromArray(array_data)
sitk.WriteImage(image, 'xxxx')

示例： nii 与 npy 相互转化

import SimpleITK as sitk
import numpy as np

# load input file
origin_data = read_img_sitk(filepath="origin.nii")
print(origin_data.shape)  # (90, 128, 128)

# save as npy file
np.save(file='origin.npy', arr=origin_data)

# load npy file
npy_data = np.load('origin.npy')
print(npy_data.shape)  # (90, 128, 128)

# save as nii file
image = sitk.GetImageFromArray(npy_data)
sitk.WriteImage(image, 'origin_npy.nii')

# load saved file
origin_npy_data = read_img_sitk(filepath="origin_npy.nii")
print(origin_npy_data.shape)  # (90, 128, 128)

左右图数据应该是一致的，左图由Mango自动翻转slice排序显示，故而显示有所差异 ⤴️

show_img_sitk(origin_data, 45)
show_img_sitk(origin_npy_data, 45)

Nibabel

• ReadImage

import nibabel as nib
import matplotlib.pyplot as plt
 
def read_img_nib(filepath):
    image_data = nib.load(filepath).get_data()  # <class 'numpy.memmap'>
    return image_data
 
def show_img_nib(data, slice_num):
    import matplotlib.pyplot as plt
    plt.imshow(data[:, :, slice_num], cmap='gray')  # channel_last
    plt.show()
 
path = 'origin.nii'
data = read_img_nib(path)  # (W, H, C, 1)
print(data.shape)  # (128, 128, 90, 1)
data_arr = np.asarray(data)  # <class 'numpy.memmap'> -> <class 'numpy.ndarray'>
data_arr = np.squeeze(data_arr)  # (W, H, C)  # or data.squeeze()
print(data_arr.shape)  # (128, 128, 90)
show_img_nib(data_arr, 45)

# correction nib offset
data_arr = data_arr.transpose((1, 0, 2))  # (H, W, C)
print(data_arr.shape)  # (128, 128, 90)
show_img_nib(data_arr, 45)

transpose PET slice

• WriteImage

# load input file
origin_data = read_img_nib(filepath="origin.nii")
print(origin_data.shape)  # (128, 128, 90, 1)  (W*H*C*1)

# save as npy file
np.save(file='origin.npy', arr=origin_data)

# load npy file
npy_data = np.load('origin.npy')
print(npy_data.shape)  # (128, 128, 90, 1)  (W*H*C*1)
npy_data = np.squeeze(npy_data)
print(npy_data.shape)  # (128, 128, 90)   (W*H*C)

# save as nii file
# 默认以(C' * W' * H')传入，即误认为传入128张128*90的切片
image = sitk.GetImageFromArray(npy_data)  
sitk.WriteImage(image, 'origin_npy.nii')

# load saved file
origin_npy_data = read_img_nib(filepath="origin_npy.nii")
print(origin_npy_data.shape)  # (90, 128, 128)   (W*H*C)  # 128张128*90的切片 ❎

显然，右图是未对origin_data做像素归一化所致 ⤴️，修改如下：

def normalization(data):
    data = np.asarray(data).squeeze()
    max, min = np.max(data), np.min(data)
    data = (data - min) / (max - min)
    return data

# load input file
origin_data = read_img_nib(filepath="origin.nii")  # (128, 128, 90, 1)  (W*H*C*1)
origin_data = normalization(origin_data)  # (128, 128, 90)  (W*H*C)

# save as npy file
np.save(file='origin.npy', arr=origin_data)

# load npy file
npy_data = np.load('origin.npy') 
print(npy_data.shape)  # (128, 128, 90)  # (W, H, C)

# save as nii file
npy_data = npy_data.transpose(2, 1, 0)
print(npy_data.shape)  # (90, 128, 128)  # (C, H, W)
image = sitk.GetImageFromArray(npy_data)
sitk.WriteImage(image, 'nib_origin_npy.nii')

# load saved file
origin_npy_data = read_img_nib(filepath="nib_origin_npy.nii")
print(origin_npy_data.shape)  # (128, 128, 90)  (W*H*C)  # 90张128*128的切片 ✅

show_img_nib(origin_data, 45)
show_img_nib(origin_npy_data, 45)

中心裁剪

• 数组操作

def center_crop_sitk(img_arr, width_target, height_target):
    """
    中心裁剪任意尺寸的图片（以中心为原点）
    """
    channel, height, width = img_arr.shape

    if width < width_target:
        raise Exception(
            "The width of image must be at least {} pixels!".format(width_target))
    elif height < height_target:
        raise Exception(
            "The height of image must be at least {} pixels!".format(height_target))

    width_crop = (width - width_target) // 2
    height_crop = (height - height_target) // 2

    if width_crop > 0:
        img_arr = img_arr[:, :, width_crop:-width_crop]
    if height_crop > 0:
        img_arr = img_arr[:, height_crop:-height_crop, :]

    return img_arr

# load input file
origin_data = read_img_sitk(filepath="origin.nii")
print(origin_data.shape)  # (90, 128, 128)

# crop with size
crop_data = center_crop_sitk(origin_data, width_target=88, height_target=112)
print(crop_data.shape)  # (90, 112, 88)  (C*W*H) sitk

# compare
show_img_sitk(origin_data, 45)
show_img_sitk(crop_data, 45)

• 自定义裁剪

origin_data = read_img_sitk(filepath="origin.nii")  # sitk
print(origin_data.shape)  # (90, 128, 128)

# height_target[0] + height_target[1] <  origin_data.shape[1]
height_target = (0, 20)  # 确实修改的是高度 ➜ origin_data.shape[1] 确实代表的是H ➜ origin_data.shape：C*H*W (可验证)
# width_target[0] + width_target[1] <  origin_data.shape[2]
width_target = (10, 10) 
crop_data = origin_data[:, height_target[0]:-height_target[1], width_target[0]:-width_target[1]]
print(crop_data.shape)  # (90, 108, 108)

# compare
show_img_sitk(origin_data, 45)
show_img_sitk(crop_data, 45)