【深度学习】图像自动处理工具SimpleITK的使用（Python）

文章目录
1 Python安装第三方库SimpleITK
2 SimpleITK处理医学图像
3 图像增强、去噪、边缘检测（SimpleITK）
4 使用SimpleITK读取和保存Nii.gz文件

1 Python安装第三方库SimpleITK

命令行模式下，任意位置输入：

  pip install SimpleITK

第三方库会自动存在Python安装路径下的Lib\site-packages目录下。
测试代码

import SimpleITK as sitk

import sys

import os

file = r'C:\RedPanda.jpg'

example = sitk.ReadImage(file)

sitk.Show(example)

注意事项：

有时提示没有安装SimpleITK，可能是选择的解释器位置不对。将解释器设置为Python安装路径下的Python.exe.（安装了第三方库SimpleITK）

代码中.Show()函数首字母用大写，Python是大小写敏感的。不然可能会报错AttributeError: module 'SimpleITK' has no attribute 'show'

2 SimpleITK处理医学图像

做医学图像时，SimpleITK是一个很常用的库。实际上大家往往喜欢把不同类型的数据割裂开，nrrd用pynrrd处理，dicom用dicom处理，nii用nibabel处理……实际上根本没必要，SimpleITK完全可以统一处理，各种类型的读取和保存一步搞定。

1. 读取

首先是最常见的DICOM。这个文件格式是CT机器直接输出的结果，每个病人会是一个文件夹，里面有若干个DICOM文件，每个DICOM文件都是一个切片，可以被单独查阅（通过SimpleITK的ReadImage）。但是，处理起来一般我们会希望把这些DICOM组合起来形成一个3D的矩阵，这样就可以用到SimpltITK里的ImageSeriesReader()来实现：

reader = sitk.ImageSeriesReader()
img_names = reader.GetGDCMSeriesFileNames(img_path)
reader.SetFileNames(img_names)
image = reader.Execute()
image_array = sitk.GetArrayFromImage(image) # z, y, x
然后还有就是MHD。这种格式的文件由一个MHD文件和一个RAW文件组成，其中MHD里面是病人和CT扫描结果的相关信息，RAW里面存的是真正的数组，对于这种格式，SimpleITK可以这样读取：

itk_img = sitk.ReadImage(filename)
img_array = sitk.GetArrayFromImage(itk_img)
其次，还有NII格式的文件，想要读取也是类似的代码：

itk_img = sitk.ReadImage(filename)
img_array = sitk.GetArrayFromImage(itk_img)
不过，对于NII类型的文件，SimpleITK在保存上会有一些问题（读取没问题保存上竟然有问题我真是服了），如果slice的数量很大（之前有碰到过300多个slice的），保存起来会直接报错，所以如果碰到这种情况我会选择保存成nrrd（用pynrrd库）

2. 保存

说到保存，有几个天坑一定要说一下。

首先，保存的通用代码是：

savedImg = sitk.GetImageFromArray(outputs3D)
savedImg.SetSpacing(spacing)
savedImg.SetOrigin(origin)
sitk.WriteImage(savedImg, input_nii[:-4]+'_liverSegResult.nii')
这个里面，SetOrigin和SetSpacing我常常会忘记，所以导致保存的CT用ITKsnap看会变形。这里面的spacing和origin均可以从之前ReadImage()返回的对象上通过GetOrigin()和GetSpacing()获取，千万别忘了。当然，除了这两个之外还有个direction也是可以get和set的，不过一般情况下只要你代码没有搞一些矩阵的翻转等操作这个是不用管的。

3 图像增强、去噪、边缘检测（SimpleITK）

1、对数变换
图像对数变换首先将图像从SimpleITK图像数据转成Numpy矩阵数据，然后采用Numpy的log1p（）函数来计算数据的log（1+x）变换，由于1+x不能小于零，因此这里我们使用图像减去图像的最小值来计算对数变换结果。
2、幂次变换
图像对数变换首先将图像从SimpleITK图像数据转成Numpy矩阵数据，然后采用Numpy的power（）函数来计算数据的幂次变换，为了防止出现计算值溢出错误，因此这里我们使用图像减去图像均值再除以图像方差来计算图像幂次变换结果，，在这里我们计算图像3次幂变换。
3、指数变换
图像指数变换首先将图像从SimpleITK图像数据转成Numpy矩阵数据，然后采用Numpy的exp（）函数来计算数据的指数变换，为了防止出现计算值溢出错误，因此这里我们使用图像减去图像均值再除以图像方差来计算图像指数变换结果。
4、直方图均衡化
使用函数AdativeHistogramEqualizationImageFilter（）函数来实现图像的直方图均衡化，这个类函数需要人为设置三个参数：Alpha，Beta，Radius，其中Alpha参数是用来控制结果相对于经典直方图均衡化方法结果的相似程度，Beta参数用来控制图像锐化程度，Radius用来控制直方图统计时的区域大小。
5、均值滤波
均值滤波是典型的线性滤波算法，它是指在图像上对目标像素给一个模板，该模板包括了其周围的临近像素（以目标像素为中心的周围8个像素，构成一个滤波模板，即去掉目标像素本身），再用模板中的全体像素的平均值来代替原来像素值。均值滤波本身存在着固有的缺陷，即它不能很好地保护图像细节，在图像去噪的同时也破坏了图像的细节部分，从而使图像变得模糊，不能很好地去除噪声点。我们使用MeanImageFilter（）函数来对图像进行平滑去噪。参数设置就是统计均值计算的半径大小，这里可以xyz方向上都一样的半径大小，也可以设置成不同的半径大小，在这里我们设置成5，就意味着是10x10x10区域里计算均值滤波的结果。
6、中值滤波
中值滤波的原理很简答，对一副图像上的某一点作为中值滤波器，必须先将预求的像素及其邻点内的像素值排序，确定出中值，并将中值赋予该像素点，让周围的像素值接近的真实值，从而消除孤立的噪声点，中值滤波对消除椒盐噪声非常有效。我们使用MedianImageFilter（）函数来对图像进行中值滤波去噪。参数设置就是统计中值计算的半径大小，这里可以xyz方向上都一样的半径大小，也可以设置成不同的半径大小，在这里我们设置成5，就意味着是10x10x10区域里计算中值滤波的结果。7
7、高斯滤波
高斯滤波是一种线性平滑滤波，适用于消除高斯噪声，广泛应用于图像处理的减噪过程。我们使用SmoothingRecursiveGaussianImageFilter（）函数来对图像进行高斯滤波平滑去噪。参数设置就是高斯核参数Sigma大小，这里可以xyz方向上都一样的大小，也可以设置成不同的大小，还有一个参数设置是是否要对高斯进行尺度归一化。
8、Sobel算子
使用SobelEdgeDetectionImageFilter（）函数来计算图像梯度，梯度是有幅值和方向的，返回的Sobel梯度结果是梯度幅值。没有需要输入的参数，直接输入2d或3d图像就可以了，这边还要提醒一下，输入的图像必须是float类型的图像才可以计算。
9、canny算子
使用CannyEdgeDetectionImageFilter()函数提取图像的边缘。参数有4个需要人为设置，双边阈值的上限阈值和下限阈值，高斯平滑由两个参数来控制高斯核参数大小，高斯最大误差，必须设置在0~1之间。这边还要提醒一下，输入的图像必须是float类型的图像才可以计算。

直方图均衡化是一种简单有效的图像增强技术，通过改变图像的直方图来改变图像中各像素的灰度，主要用于增强动态范围偏小的图像的对比度。原始图像由于其灰度分布可能集中在较窄的区间，造成图像不够清晰。例如，过曝光图像的灰度级集中在高亮度范围内，而曝光不足将使图像灰度级集中在低亮度范围内。采用直方图均衡化，可以把原始图像的直方图变换为均匀分布（均衡）的形式，这样就增加了像素之间灰度值差别的动态范围，从而达到增强图像整体对比度的效果。换言之，直方图均衡化的基本原理是：对在图像中像素个数多的灰度值（即对画面起主要作用的灰度值）进行展宽，而对像素个数少的灰度值（即对画面不起主要作用的灰度值）进行归并，从而增大对比度，使图像清晰，达到增强的目的。举个例子，如图1所示，左图为原始图像，右图为直方图均衡化后的图像。

灰度图像增强对比度。

高斯滤波
先引入两个问题。
1.图像为什么要滤波？

答：a.消除图像在数字化过程中产生或者混入的噪声。
b.提取图片对象的特征作为图像识别的特征模式。

2.滤波器该如何去理解?
答：滤波器可以想象成一个包含加权系数的窗口或者说一个镜片，当使用滤波器去平滑处理图像的时候，就是把通过这个窗口或者镜片去看这个图像。

滤波器分为很多种，有方框滤波、均值滤波、高斯滤波等。
高斯滤波是一种线性平滑滤波，适用于消除高斯噪声。
所以在讲高斯滤波之前，先解释一下什么是高斯噪声？

高斯噪声
首先，噪声在图像当中常表现为一引起较强视觉效果的孤立像素点或像素块。简单来说，噪声的出现会给图像带来干扰，让图像变得不清楚。
高斯噪声就是它的概率密度函数服从高斯分布（即正态分布）的一类噪声。如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度又是均匀分布的，则称它为高斯白噪声。高斯白噪声的二阶矩不相关，一阶矩为常数，是指先后信号在时间上的相关性。

高斯滤波器是根据高斯函数的形状来选择权值的线性平滑滤波器

4 使用SimpleITK读取和保存Nii.gz文件

## using simpleITK to load and save data.
import SimpleITK as sitk
itk_img = sitk.ReadImage('./nifti.nii.gz')
img = sitk.GetArrayFromImage(itk_img)
print("img shape:",img.shape)

## save 
out = sitk.GetImageFromArray(img)
# # out.SetSpacing(itk_img.GetSpacing())
# # out.SetOrigin(itk_img.GetOrigin())
sitk.WriteImage(out,'simpleitk_save.nii.gz')

PS：可视化一层

下面的程序挑取一层进行可视化，并标注结节位置。

a = image_array.transpose(1,2,0)[:,:,0]  # transpose是将(z,x,y)的三维矩阵转为(x,y,z)的矩阵
plt.gca().add_patch( plt.Rectangle((147,297), 24,24, fill=False,edgecolor='r', linewidth=3))
plt.imshow(a[:,:,1]*255)  # 在图中画框
plt.show()