参考：
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)
网友评论