SimpleITK学习笔记

SimpleITK学习笔记

前言

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

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

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

import SimpleITK as sitk

1 sitk中的常见属性值

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

print(image.GetSize())
print(image.GetOrigin())
print(image.GetSpacing())
print(image.GetDirection())

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

print(image.GetSize())
(7, 6)

print(image.GetOrigin())
(60.0, 70.0)

print(image.GetSpacing())
(20.0, 30.0)

print(image.GetDirection())
(1.0, 0.0, 0.0, 1.0)

SimpleITK学习笔记

2 读取和保存图像

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


imagepath = "xxx.mhd"
writepath = "xxx.mhd"
image = sitk.ReadImage(imagepath)
sitk.Write(image, writepath)
image1 = sitk.ReadImage(imagepath, sitk.sitkFloat32)

Note1:图像访问是以 xyz顺序 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位整数sitkFloat3232位浮点sitkFloat6464位浮点

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

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

SimpleITK学习笔记

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

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

Note

SimpleITK中,各术语对应如下

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

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

举个例子:假设实验用的图片大小为32025080,即矢状面(x轴方向)切片数为320,冠状面(y轴方向)切片数为250,横断面(z轴方向)片数为80。

SimpleITK2Numpy:

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


shape_img = image.GetSize()
print("image size:", shape_img)

datanp_array = sitk.GetArrayFromImage(image)
print("np_array size:", np_array.shape)

Numpy2SimpleITK:

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


imagesitk_image = sitk.GetImageFromArray(np_array)
sitk_image.SetOrigin(origin)
sitk_image.SetSpacing(spacing)
sitk_image.SetDirection(direction)

5 访问像素和切片

两种方式:一是使用 GetPixelSetPixel函数,二是使用python的切片操作符。例子如下:


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))

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.SetReferenceImage(itkimage)
resampler.SetSize(newSize.tolist())
resampler.SetTransform(sitk.Transform(3, sitk.sitkIdentity))
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
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]

    ori_size = image.GetSize()
    ori_spacing = image.GetSpacing()
    transform = sitk.Transform()
    transform.SetIdentity()

    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:

        resampler.SetOutputPixelType(sitk.sitkUInt8)
    else:

        resampler.SetOutputPixelType(sitk.sitkFloat32)
    newImage = resampler.Execute(image)
    return newImage

7 图像分割

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


sitk.BinaryThreshold(
    sitk_src,
    lowerThreshold=lowervalue,
    upperThreshold=uppervalue,
    insideValue=255,
    outsideValue=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连通。


ccfilter = sitk.ConnectedComponentImageFilter()
ccfilter.SetFullyConnected(True)
output = ccfilter.Execute(image)
count = ccfilter.GetObjectCount()

lssfilter = sitk.LabelShapeStatisticsImageFilter()
lssfilter.Execute(output)

area_max_label = 0
area_max = 0
for i in range(1, count + 1):
    area = lssfilter.GetNumberOfPixels(i)
    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),其中原点在图片的左上角, widthcolumnheightrow。举个例子,一张图片为1280×960的分辨率,即 width为1280个像素, height为960个像素。

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

下面分别用 opencvpillowskimage读取图片。

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

filepath = './xxx.jpg'

cv2_im = cv2.imread(filepath)
print('cv2_im shape ',cv2_im.shape)

im = Image.open(filepath)
print('PIL image size', im.size)
pil_im = np.asarray(im)
print('PIL_im shape ',pil_im.shape)

sk_im = imread(filepath)
print('sk_im shape', sk_im.shape)

summary:

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

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


img = cv2.imread(path)
img = img[:, :, ::-1].transpose(2, 0, 1)

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

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


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

import cv2
image = cv2.imread("../xxx.jpg")
new_image= cv2.resize(image, (width,height), interpolation=cv2.INTER_CUBIC)

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

Original: https://blog.csdn.net/YYY_77/article/details/125547278
Author: Jiaxxxxxx
Title: SimpleITK学习笔记

原创文章受到原创版权保护。转载请注明出处:https://www.johngo689.com/518748/

转载文章受原作者版权保护。转载请注明原作者出处!

(0)

大家都在看

亲爱的 Coder【最近整理,可免费获取】👉 最新必读书单  | 👏 面试题下载  | 🌎 免费的AI知识星球