SimpleITK学习笔记

前言

SimpleITK是专门处理医学影像的软件，是ITK的简化接口，使用起来非常便捷，SimpleITK支持8种编程语言，包括c++、Python、R、Java、c#、Lua、Ruby和TCL。在SimpleITK中，图像的概念与我们在计算机视觉中常用的RGB图像差异很大，后者只是一个多维矩阵，是一个数学上的概念，而在SimpleITK中图像是一个物理实体，图像中的每一个像素都是物理空间中的一个点，不光有着像素值，还有坐标、间距、方向等概念（这些概念我们后续会介绍）。

下面我们用python对SimpleITK库中常用的函数进行举例说明。

在使用SimpleITK库之前，需要将SimpleITK库导入进来，如下

import SimpleITK as sitk

1 sitk中的常见属性值

sitk中有以下四种常见的属性，分别可以使用get的方式获取，代码如下所示：

print(image.GetSize()) #获取图像的大小，size为图像的每一个维度的长度，即每个维度像素点的个数
print(image.GetOrigin()) #获取图像的原点坐标
print(image.GetSpacing()) #获取每个维度上像素或体素之间的间距，单位mm，其中对于二维图像是像素，三维图像是体素
print(image.GetDirection()) #获取图像的方向，即图像坐标系相对世界坐标系的角度，角度采用的是方向余弦矩阵

以二维图像为例，如下图所示，左下图是世界坐标系，右下图是图像坐标系。属性值如下所示：

print(image.GetSize())
（7， 6） #在SimpleITK中，读取的顺序为[X, Y, Z], 即先宽度，后高度，再深度，对应的三维医学图像中为：先矢状位，后冠状位，再横断位。

print(image.GetOrigin())
（60.0， 70.0） #读取顺序同GetSize()方法

print(image.GetSpacing())
（20.0， 30.0） #读取顺序同GetSize()方法， 每个维度可以具有不同的间距，且每个维度不一定是正交的。

print(image.GetDirection())
（1.0， 0.0， 0.0， 1.0） 

2 读取和保存图像

SimpleITK可以读取如.mhd , .nii, .nrrd等图像数据格式。

#以读取和保存mhd图像为例
imagepath = "xxx.mhd"   #图像路径
writepath = "xxx.mhd"  #保存图像路径
image = sitk.ReadImage(imagepath）
sitk.Write(image, writepath)
image1 = sitk.ReadImage(imagepath, sitk.sitkFloat32） #你可以指定读取的数据类型

Note1：图像访问是以x， y， z顺序GetPixel(x,y,z) 或 image[x,y,z]， 从0开始索引。

Note2：默认的mask图像类型和默认值为sitkUInt8或标量图像uint8_t， 默认值为0和1，其中1代表的是mask。

3 像素类型

像素类型表示为枚举类型，下面是部分枚举类型的表。

数据类型	含义
sitkUInt8	无符号8位整数
sitkInt8	有符号的8位整数
sitkUInt16	无符号16位整数
sitkInt16	有符号的16位整数
sitkFloat32	32位浮点
sitkFloat64	64位浮点

还有sitkUnknown类型，用于未定义或错误的像素ID，值为-1。

Note：64位整数类型并非在所有的发行版上都可用，如果不可用，则值为sitkUnknown，将图像保存为nii文件，用ITKsnap读取时就会出现的错误如下：

4 SimpleITK图像数据和Numpy矩阵数据之间的转换

般我们会用SimpleITK读取图像，再转换为numpy矩阵格式，这样方便数据的处理。

Note

在SimpleITK中，各术语对应如下

Width: 宽度，X轴，矢状面（Sagittal）
Height: 高度，Y轴，冠状面（Coronal）
Depth: 深度， Z轴，横断面（Axial）

SimpleITK图像顺序是x，y，z三个方向的大小（在第一节中也讲过），而numpy矩阵的顺序是z，y，x三个方向的大小, 要注意索引位置。

举个例子：假设实验用的图片大小为320*250*80，即矢状面（x轴方向）切片数为320，冠状面（y轴方向）切片数为250，横断面（z轴方向）片数为80。

SimpleITK2Numpy：

GetArrayFromImage()：返回图像数据的副本。然后可以自由地修改数据，因为它对原始SimpleITK图像没有影响。

# sitk image to numpy 
shape_img = image.GetSize() #输出形状为：(Width, Height, Depth)，即原始SimpleITK数据的存储形式
print("image size:", shape_img)
#image size:(320, 250, 80)

datanp_array = sitk.GetArrayFromImage(image) # 将SimpleITK对象转换为ndarray，X轴与Z轴发生了对调，输出形状为：(Depth, Height, Width)
print("np_array size:", np_array.shape)
# np_array size:(80, 250, 320)

Numpy2SimpleITK：

GetImageFromArray()：返回一个SimpleITK图像，原点设置为0，所有维度的间距设置为1，方向余弦矩阵设置为identity，强度数据从numpy数组中复制。在大多数情况下需要设置适当的元数据值。

# numpy data to sitk 
imagesitk_image = sitk.GetImageFromArray(np_array) 
sitk_image.SetOrigin(origin)
sitk_image.SetSpacing(spacing)
sitk_image.SetDirection(direction)

5 访问像素和切片

两种方式：一是使用GetPixel和SetPixel函数，二是使用python的切片操作符。例子如下：

#法一：使用GetPixel和SetPixel函数
image_3D = sitk.Image(256, 128, 64, sitk.sitkInt16)
print(image_3D.GetPixel(0, 0, 0))  #注意：是从0开始索引
image_3D.SetPixel(0, 0, 0, 1)
print(image_3D.GetPixel(0, 0, 0))

#法二：使用python的切片操作符
print(image_3D[0,0,1])
image_3D[0,0,1] = 2
print(image_3D[0,0,1])

6 图像重采样

重采样目的是将医疗图像中不同的体素归一化到相同的大小，可分为两种方案，方案一：按目标Spacing缩放，方案二：按目标Size缩放。

这两种方案具体操作分为三个步骤：

1.使用SimpleITK读取数据，获得原始图像的Spacing以及Size；

2.如果是方案一，则图像原始Size乘以原始Spacing除以新Spacing得到新Size，如果是方案二，则图像原始Size乘以原始Spacing除以新Size得到新Spacing；

3.最后将新Spacing和新Size赋值到读取的数据即可。

核心代码如下：

resampler = sitk.ResampleImageFilter() #初始化一个图像重采样滤波器resampler
resampler.SetReferenceImage(itkimage) #设置需要重新采样的目标图像
resampler.SetSize(newSize.tolist()) #设置目标尺寸
resampler.SetTransform(sitk.Transform(3, sitk.sitkIdentity))
resampler.SetInterpolator(sitk.sitkNearestNeighbor) #设置插值方式,体数据采用线性插值，mask采用最近邻插值
itkimgResampled = resampler.Execute(itkimage) #得到重新采样后的图像

下面以指定Spacing大小，对原始数据进行重采样，例子如下：

def resampleImage(image, targetSpacing, resamplemethod=sitk.sitkNearestNeighbor):
    """
    将体数据重采样的指定的spacing大小
    paras：
    image：sitk读取的image信息，这里是体数据
    targetSpacing：指定的spacing，例如[1,1,1]
    resamplemethod:插值类型
    return：重采样后的数据
    """
    targetsize = [0, 0, 0]
    #读取原始数据的size和spacing信息
    ori_size = image.GetSize()
    ori_spacing = image.GetSpacing()
    transform = sitk.Transform()
    transform.SetIdentity()
    #计算改变spacing后的size，用物理尺寸/体素的大小
    targetsize[0] = round(ori_size[0] * ori_spacing[0] / targetsize[0])
    targetsize[1] = round(ori_size[1] * ori_spacing[1] / targetsize[1])
    targetsize[2] = round(ori_size[2] * ori_spacing[2] / targetsize[2])
    #设定重采样的一些参数
    resampler = sitk.ResampleImageFilter()
    resampler.SetTransform(transform)
    resampler.SetSize(targetsize)
    resampler.SetOutputOrigin(image.GetOrigin())
    resampler.SetOutputSpacing(targetSpacing)
    resampler.SetOutputDirection(image.GetDirection())
    resampler.SetInterpolator(resamplemethod)
    if resamplemethod=sitk.sitkNearestNeighbor:
        #mask用最近邻插值，保存为uint8
        resampler.SetOutputPixelType(sitk.sitkUInt8)
    else：
    	#体数据用线性插值，保存为float32
        resampler.SetOutputPixelType(sitk.sitkFloat32)
    newImage = resampler.Execute(image)
    return newImage

7 图像分割

图像二值化分割： 是分割方法中最基础的，通过定义lowerThreshold和upperThreshold两个像素临界点，只要像素值在两者之间，则该像素值改为insideValue，否则改为outsideValue。这种方法只是简单的基于灰度范围标记图像像素，不考虑几何或连通性。

#sitk API
sitk.BinaryThreshold(
    sitk_src, 
    lowerThreshold=lowervalue, 
    upperThreshold=uppervalue, 
    insideValue=255, 
    outsideValue=0) 

#例子：将0-500之间的像素置为1，其他置为0
image = sitk.ReadImage("./***.nii.gz")
outimage = sitk.BinaryThreshold(
        image,
        lowerThreshold=0,
        upperThreshold=500,
        insideValue=1,
        outsideValue=0
    )

图像区域生长分割： 需要确定种子点、生长准则和终止条件。具体来说，对每一个需要分割的区域找一个种子像素作为生长的起点，根据生长准则将种子像素邻域中与种子像素具有相同或相似的像素合并到种子像素所在的区域，直到没有满足条件的像素可以被包进来就终止。在SimpleITK中，首先会计算当前区域中包含的所有像素点灰度值的标准差和期望，通过定义multiplier因子（乘以标准差）来计算以期望为中心的灰度值范围，如果initialNeighborhoodRadius邻域半径内的灰度值位于这个范围就被包进来，灰度值改为replaceValue，当遍历了所有邻域像素，即认为完成了一次迭代，下一次迭代时，像素点的灰度值期望和标准差是以新的像素区域为基础进行计算的，一共要迭代numberOfIterations次。

sitk.ConfidenceConnected(
    image,
    seedList=[seed],
    numberOfIterations=1,
    multiplier=2.5,
    initialNeighborhoodRadius=1,
    replaceValue=1)

8 图像的形态学操作

# 图像的形态学操作：开、闭、膨胀、腐蚀
sitk.BinaryMorphologicalOpening(sitk.ReadImage() != 0, kernelsize) 
sitk.BinaryMorphologicalClosing(sitk.ReadImage() != 0, kernelsize)
sitk.BinaryDilate(sitk.ReadImage() != 0, kernelsize)
sitk.BinaryErode(sitk.ReadImage() != 0, kernelsize)

9 连通域分析

连通域分析一般是针对二值图像，将具有相同像素值且相邻的像素找出来并标记，其中二维图像连通域一般包括4连通和8连通，三维图像连通域包括6连通、18连通和26连通。

# ConnectedComponentImageFilter类，用于连通域划分
ccfilter = sitk.ConnectedComponentImageFilter()
ccfilter.SetFullyConnected(True)
output = ccfilter.Execute(image)  #image为要进行连通域划分的二值图像，output为输出被标记的图像
count = ccfilter.GetObjectCount()

# LabelShapeStatisticsImageFilter类，用于获取各label的形状属性
lssfilter = sitk.LabelShapeStatisticsImageFilter()
lssfilter.Execute(output) 

#根据连通域的体积大小保留最大连通域
area_max_label = 0  #记录最大的连通域的label
area_max = 0 #记录连通域中最大的体积
for i in range(1, count + 1):
    area = lssfilter.GetNumberOfPixels(i)  # 根据label获取连通域体积
    if area > area_max:
        area_max_label = i
        area_max = area
output_array = sitk.GetArrayFromImage(output)
res = np.zeros_like(output_array)
res[output_array == area_max_label] = 1 #获取最大连通域

OpenIssue

1 图片读取的宽高顺序

先上结论：

① 图片一般表示为(width,height)，其中原点在图片的左上角，width为column，height为row。举个例子，一张图片为1280x960的分辨率，即width为1280个像素，height为960个像素。

② 图片的存储形式为矩阵，在数学概念中，一般矩阵的表达中第一个值为row，第二个值为column，对应为(height,width)。因此图片为(width,height)转化为numpy数组后变为(height,width)。

下面分别用opencv、pillow和skimage读取图片。

import cv2
from PIL import Image
import numpy as np
from skimage.io import imread

filepath = './xxx.jpg'   #读取图片的路径

#opencv读取：opencv读进来的图片已经是一个numpy矩阵了
cv2_im = cv2.imread(filepath)
print('cv2_im shape ',cv2_im.shape) # (height, width, ch)

#pillow读取：pillow都进来的图片是一个对象，因此还需要进一步转化为numpy矩阵
im = Image.open(filepath)
print('PIL image size', im.size) # (width, height, ch)
pil_im = np.asarray(im)
print('PIL_im shape ',pil_im.shape) # (height, width, ch)

#skimage读取：同opencv
sk_im = imread(filepath)
print('sk_im shape', sk_im.shape) # (height, width, ch)

summary：

①opencv读取和存储的格式是BGR，而不是主流的RGB。

②在深度学习进行推理阶段的前处理时，一般要将RGB图像或者BGR图像转化为CHW格式，一般来说，opencv、pillow和skimage读取图片并转化为矩阵都是HWC格式，而CHW更适合CNN，因为网络是逐个通道的对图像进行卷积运算，希望在访问同一个channel的像素是连续的，一般存储为CHW格式，这样访问内存的时候是连续的，比较方便。

#yolov5读取图片
img = cv2.imread(path)
img = img[:, :, ::-1].transpose(2, 0, 1)  # img[:, :, ::-1]是将BGR to RGB，transpose(2, 0, 1)是将HWC转换成CHW格式

2 图像重采样可用的工具包

下面分别介绍PIL、opencv工具包对图像进行重采样。

#PIL重采样
from PIL import Image
image = Image.open("./xxx.jpg")
new_image = image.resize((width, height), Image.BILINEAR)
new_image.save("./xxx.jpg")

#opencv重采样
import cv2
image = cv2.imread("../xxx.jpg")
new_image= cv2.resize(image, (width,height), interpolation=cv2.INTER_CUBIC)

Note：PIL，CV2一般针对2D图像做处理，比如普通的自然图像，或者CT三维图像中的某一层。所以，并不适用于3D的医学图像处理。