SimpleITK的部分介绍及使用代码（python） 简单介绍了学习中SimplelTK的部分使用方法，包括读取数据，重采样，图像分割，图像形态学操作，连通域分析等
SimpleITK的部分介绍及使用代码
1.简单介绍
SimplelTK是在Insight Toolkit(ITK)的基础上构建的简化工具，常用于处理医学图像。它支持8种编程语言，包括C++、Python、R、Java、C#、Lua、 Ruby和TCL。SimplelTK可以很便捷地读取医学图像的格式，如二维的.dicom格式以及三维图像的.nii格式，并对图像中的诸如体素、方向等等基本信息进行处理。SimplelTK中的图像是一种物理实体，图像中的每一个像素都是物理空间中的一个点，不仅有像素值，还有坐标、间距、方向等概念。
2.SimpleITK基本概念
参考SimpleITK文档：Fundamental Concepts — SimpleITK 2.0rc2 documentation
如果看得懂的话，请以官方文档为准，以下皆为个人理解和学习时的笔记，仅供参考。
2.1 SimpleITK图像基础概念
SimpleITK默认设定世界坐标系下的长度单位为mm。
世界坐标系：用于医学成像仪器的坐标系。坐标用元素是连续值的向量表示。
图像坐标系：用于表示像素点的坐标系。坐标用元素是离散值的向量表示。
基础属性：
Origin (vector like type) - 世界坐标系下，图像相对于坐标系原点的位置。
Spacing - 两个像素之间在世界坐标系下的间隔。因为图像数据（即像素）以离散值存储，而像素的间隔会因成像仪器的不同而有所不同。需要注意的是各个方向的Spacing通常是不相等的，例如(0.67777,0.67777,1.33333)，其实也就是体素的大小。
Size - 每个维度的大小。用元组表示。
Direction cosine matrix - 方向余弦矩阵。世界坐标系坐标轴对应于图像像素矩阵的方向。用一组矩阵空间内的基表示。通常每个向量形似(1,0,0).T。
官方演示图（2D）
3.SimpleITK的基本使用
在使用SimpleITK库之前，需要将SimpleITK库导入进来，如下：
import SimpleITK as sitk3.1 读取和保存图像
SimpleITK可以使用ReadImage()读取如.mhd , .nii, .nrrd等图像数据格式。
相关参数：
fileName : str or list of strs ： 读取图像的路径。
outputPixelType : 枚举类型 ： 将被作为返回值的图像自动转换成规定的像素类型。默认为sitkUnknown，用于未定义或错误的像素ID，值为-1。
imageIO : str ： 读取器的名字，类型是字符串。可用的读取器通过GetRegisteredImageIOs()查看。默认为""。不指定时，将自动使用对应的读取器读取图像，所以一般可以不指定。
关于GetRegisteredImageIOs()可查看官方说明：读取和写入图像和转换 — SimpleITK 2.0rc2 文档
返回值:
SimpleITK.SimpleITK.Image ： ITK图像类型。
代码实例
#以读取和保存nii图像为例
imagepath = "xxx.nii" #图像路径
writepath = "xxx.nii" #保存图像路径
image = sitk.ReadImage(imagepath) # 根据路径读取图像
sitk.WriteImage(image, writepath) # 根据路径保存图像
image1 = sitk.ReadImage(imagepath, sitk.sitkFloat32） #你可以指定读取的像素类型读取的像素类型表示为枚举类型，下面是部分枚举类型的表。
根据SimpleITK库的版本不同，所能支持的像素类型可能有所不同，大部分的都支持sitkUInt8、16、32、64，sitkInt8、16、32、64，sitkFloat32、64。
3.2 获取图像属性
前面我们提到过，ITK图像具有Origin，Spacing，Size，Direction cosine matrix，共四个基本属性。我们可以分别使用Get的方式获取，代码如下所示：
# 获取图像的大小，size为图像的每一个维度的长度，即每个维度像素点的个数
print(image.GetSize())
# 获取图像的原点坐标
print(image.GetOrigin())
# 获取每个维度上像素或体素之间的间距，单位mm，其中对于二维图像是像素，三维图像是体素
print(image.GetSpacing())
# 获取图像的方向，即图像坐标系相对世界坐标系的角度，角度采用的是方向余弦矩阵
print(image.GetDirection())其返回值皆为一个元组，在SimpleITK中，各术语对应如下：
Width: 宽度，X轴，矢状面（Sagittal）Height: 高度，Y轴，冠状面（Coronal）Depth: 深度， Z轴，横断面（Axial）
SimpleITK图像顺序是x，y，z三个方向。
例子：若假设提取的图像image通过Getsize()函数获得的大小为(224, 256, 176)，则可知该图像矢状面（x轴方向）切片数为224，冠状面（y轴方向）切片数为256，横断面（z轴方向）切片数为176。
3.3 SimpleITK图像数据与Numpy矩阵数据之间的转换
3.3.1 SimpleITK to Numpy:
GetArrayFromImage()：返回图像数据的副本。然后可以自由地修改数据，因为它对原始SimpleITK图像没有影响。
注：将SimpleITK转换为Numpy时，其x轴会与z轴发生调换。例如：
print("ITK image size:", image.GetSize())
#ITK image size:(224, 256, 176)
Numpy_image = sitk.GetArrayFromImage(image) # 将SimpleITK对象转换为ndarray
print("Numpy image size:", Numpy_image.shape)
# Numpy_image size:(176, 256, 224)3.3.2 Numpy to SimpleITK:
GetImageFromArray()：返回一个SimpleITK图像，在大多数情况下需要重新设置适当的元数据值。未重新设置时，其原点默认为0.0，所有维度的间距都默认为1.0，方向余弦矩阵默认为
[1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0]，而强度数据会从numpy数组中复制。
代码实例：
imagesitk_image = sitk.GetImageFromArray(np_array)
print(imagesitk_image.GetSize()) # 图像的大小不变
print(imagesitk_image.GetOrigin()) # 3D图像的原点坐标默认为(0.0, 0.0, 0.0)
print(imagesitk_image.GetSpacing()) # 每个维度上体素之间的间距默认为(1.0, 1.0, 1.0)
print(imagesitk_image.GetDirection())
# 3D图像的方向默认为
[1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0]
# 从array转成image一般需要重新设置原点,维度间距以及方向。
origin = (1.0, 1.0, 1.0)
spacing = (2.0, 2.0, 2.0)
direction = [2.0, 0.0, 0.0,
0.0, 2.0, 0.0,
0.0, 0.0, 2.0]
sitk_image.SetOrigin(origin) # 设置原点坐标
sitk_image.SetSpacing(spacing) # 设置每个维度上体素之间的间距
sitk_image.SetDirection(direction) # 设置图像的方向3.4访问像素和切片
两种方式：一是使用GetPixel()和SetPixel()函数，二是使用python的切片操作符。例子如下：
#1.使用GetPixel()和SetPixel()函数
image_3D = sitk.Image([256, 128, 64], sitk.sitkInt16)
print(image_3D.GetPixel(0, 0, 0))
image_3D.SetPixel(0, 0, 0, 1)
print(image_3D.GetPixel(0, 0, 0))
#2.使用python的切片操作
print(image_3D[0,0,1])
image_3D[0,0,1] = 2
print(image_3D[0,0,1])运行过后可以知道，这两种操作的效果是一样的。
3.5图像重采样
重采样目的是将ITK图像中不同的体素通过坐标变换后归一化到相同的大小，可分为两种方案：
方案一：按目标Spacing缩放。
方案二：按目标Size缩放。
这两种方案具体操作都可分为三个步骤：
1.使用SimpleITK读取数据，获得原始图像的Spacing以及Size；
2.如果是方案一，则图像原始Size乘以原始Spacing除以新Spacing得到新Size，如果是方案二，则图像原始Size乘以原始Spacing除以新Size得到新Spacing；
3.最后将新Spacing和新Size赋值到读取的数据即可。
在开始重采样时，我们需要了解一下重采样过滤器Resampler，其通过某些坐标变换对现有图像进行重新采样，并通过某些图像函数进行插值处理。在使用时应设置输出图像的输出信息（间距、原点坐标、输出大小和方向）。或者可以使用SetReferenceImage()直接从参考图像中获取输出信息进行变换。
然后我们就可以尝试用一个函数进行处理，以下代码按方案一实行：
def resampleImage(image, targetSpacing, resamplemethod=sitk.sitkNearestNeighbor):
"""
将体数据重采样的指定的spacing大小
paras：
image：sitk读取的image信息,这里是体数据
targetSpacing：指定的spacing,例如[1.0,1.0,1.0]
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] / targetSpacing[0])
targetsize[1] = round(ori_size[1] * ori_spacing[1] / targetSpacing[1])
targetsize[2] = round(ori_size[2] * ori_spacing[2] / targetSpacing[2])
# 设定重采样的一些参数
resampler = sitk.ResampleImageFilter() # 初始化图像重采样滤波器Resampler
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用最近邻插值,保存为Int16
resampler.SetOutputPixelType(sitk.sitkInt16)
else：
# 体数据用线性插值,保存为float32
resampler.SetOutputPixelType(sitk.sitkFloat32)
newimage = resampler.Execute(image) # 开始启动过滤器
return newimage3.6图像分割
1.二值化分割：
图像的二值化分割是分割方法中最基础的，可通过使用stik.BinaryThreshold函数设置lowerThreshold和upperThreshold这两个参数，只要图片的像素值在lowerThreshold和upperThreshold两者之间，则该像素值改为insideValue，否则改为outsideValue。这种方法只是简单的基于灰度范围标记图像像素。
image = sitk.ReadImage("./***.nii.gz")
outimage = sitk.BinaryThreshold(
image,
lowerThreshold=0,
upperThreshold=200, # [0,200),注意范围
insideValue=1, # 将0-255之间的像素设置为1
outsideValue=0 # 否则设置为0
)2.区域生长分割：
区域生长分割需要确定种子点、生长准则和终止条件。具体来说，对每一个需要分割的区域找一个种子像素作为生长的起点，根据生长准则将种子像素邻域中与种子像素具有相同或相似的像素合并到种子像素所在的区域，直到没有满足条件的像素可以被包进来就终止。
在SimpleITK中，首先会计算当前区域中包含的所有像素点灰度值的标准差和期望，通过定义multiplier因子（乘以标准差）来计算以期望为中心的灰度值范围，如果initialNeighborhoodRadius邻域半径内的灰度值位于这个范围就被包进来，灰度值改为replaceValue，当遍历了所有邻域像素，即认为完成了一次迭代，下一次迭代时，像素点的灰度值期望和标准差是以新的像素区域为基础进行计算的，一共要迭代numberOfIterations次。
image = sitk.ReadImage("./***.nii.gz")
sitk.ConfidenceConnected(
image,
seedList=[seed], # 种子点列表(坐标),该点通常是手动选择或基于某种预处理自动选择的
# 并且应该位于期望分割的目标区域内。
numberOfIterations=1, # 迭代次数,算法将在指定次数的迭代中尝试扩展或收缩分割区域
multiplier=2.5,
initialNeighborhoodRadius=1, # 初始邻域半径,定义了从种子点开始搜索的初始范围
replaceValue=1) # 将灰度值改为13.7图像的形态学操作
图像形态学操作一般常用的是腐蚀，膨胀，开，闭操作： 在进行这些操作时，都需要手动给一个kernel结构元素，2D中可以理解为一个正方形，3D中可以理解为一个正方体。事实上，结构元素可以根据需求具有不同的形状和大小，比如圆形，球等。
注：此时进行形态学操作的图像必须要为一个二值化的图像。
腐蚀：可以理解为结构元素会在图像上滑动，只有当结构元素范围内全部像素点的值都为1时，才为1，否则将该范围内的全部像素点的值都改为0。
膨胀：与腐蚀相反，当结构元素在图像上滑动时，当结构元素范围内有一个像素点的值1时，该范围内的全部像素值都改为1，否则为0。
开操作：先腐蚀，后膨胀。
闭操作：先膨胀，后腐蚀。
image = sitk.ReadImage("./***.nii.gz")
kernel = [3, 3, 3] # 定义腐蚀操作的核大小(3D)(可以为列表或元组)
sitk.BinaryErode(image, kernel) # 腐蚀操作
sitk.BinaryDilate(image, kernel) # 膨胀操作
sitk.BinaryMorphologicalOpening(image, kernel) # 开操作
sitk.BinaryMorphologicalClosing(image, kernel) # 闭操作3.8连通域分析
连通域分析一般是针对二值图像，将具有相同像素值且相邻的像素找出来并标记，其中二维图像连通域一般包括4连通和8连通，三维图像连通域包括6连通、18连通和26连通。
# ConnectedComponentImageFilter类,用于连通域划分
ccfilter = sitk.ConnectedComponentImageFilter()
ccfilter.SetFullyConnected(True)
# image为要进行连通域划分的二值图像,output为输出被标记的图像
output = ccfilter.Execute(image)
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 # 获取最大连通域4.参考
SimpleITK Sphinx Documentation — SimpleITK 2.0rc2 documentation