SimpleITK的部分介绍及使用代码（python）

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。官方演示图2D3.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、64sitkInt8、16、32、64sitkFloat32、64。3.2 获取图像属性前面我们提到过ITK图像具有OriginSpacingSizeDirection cosine matrix共四个基本属性。我们可以分别使用Get的方式获取代码如下所示# 获取图像的大小size为图像的每一个维度的长度即每个维度像素点的个数 print(image.GetSize()) # 获取图像的原点坐标 print(image.GetOrigin()) # 获取每个维度上像素或体素之间的间距单位mm其中对于二维图像是像素三维图像是体素 print(image.GetSpacing()) # 获取图像的方向即图像坐标系相对世界坐标系的角度角度采用的是方向余弦矩阵 print(image.GetDirection())其返回值皆为一个元组在SimpleITK中各术语对应如下Width: 宽度X轴矢状面SagittalHeight: 高度Y轴冠状面CoronalDepth: 深度 Z轴横断面AxialSimpleITK图像顺序是xyz三个方向。例子若假设提取的图像image通过Getsize()函数获得的大小为(224, 256, 176)则可知该图像矢状面x轴方向切片数为224冠状面y轴方向切片数为256横断面z轴方向切片数为176。3.3 SimpleITK图像数据与Numpy矩阵数据之间的转换3.3.1SimpleITK 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.2Numpyto 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以及Size2.如果是方案一则图像原始Size乘以原始Spacing除以新Spacing得到新Size如果是方案二则图像原始Size乘以原始Spacing除以新Size得到新Spacing3.最后将新Spacing和新Size赋值到读取的数据即可。在开始重采样时我们需要了解一下重采样过滤器Resampler其通过某些坐标变换对现有图像进行重新采样并通过某些图像函数进行插值处理。在使用时应设置输出图像的输出信息间距、原点坐标、输出大小和方向。或者可以使用SetReferenceImage()直接从参考图像中获取输出信息进行变换。然后我们就可以尝试用一个函数进行处理以下代码按方案一实行def resampleImage(image, targetSpacing, resamplemethodsitk.sitkNearestNeighbor): 将体数据重采样的指定的spacing大小 paras imagesitk读取的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, lowerThreshold0, upperThreshold200, # [0,200),注意范围 insideValue1, # 将0-255之间的像素设置为1 outsideValue0 # 否则设置为0 )2.区域生长分割区域生长分割需要确定种子点、生长准则和终止条件。具体来说对每一个需要分割的区域找一个种子像素作为生长的起点根据生长准则将种子像素邻域中与种子像素具有相同或相似的像素合并到种子像素所在的区域直到没有满足条件的像素可以被包进来就终止。在SimpleITK中首先会计算当前区域中包含的所有像素点灰度值的标准差和期望通过定义multiplier因子乘以标准差来计算以期望为中心的灰度值范围如果initialNeighborhoodRadius邻域半径内的灰度值位于这个范围就被包进来灰度值改为replaceValue当遍历了所有邻域像素即认为完成了一次迭代下一次迭代时像素点的灰度值期望和标准差是以新的像素区域为基础进行计算的一共要迭代numberOfIterations次。image sitk.ReadImage(./***.nii.gz) sitk.ConfidenceConnected( image, seedList[seed], # 种子点列表(坐标),该点通常是手动选择或基于某种预处理自动选择的 # 并且应该位于期望分割的目标区域内。 numberOfIterations1, # 迭代次数,算法将在指定次数的迭代中尝试扩展或收缩分割区域 multiplier2.5, initialNeighborhoodRadius1, # 初始邻域半径,定义了从种子点开始搜索的初始范围 replaceValue1) # 将灰度值改为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 documentationSimpleITK图像基础三——SimpleITK学习笔记_sitk-CSDN博客