SimpleITK笔记本中ITK细分(Segmentation)介绍

SimpleITK笔记本中ITK细分(Segmentation)介绍

目标：熟悉ITK中可用的基本分割算法，并以交互方式探索其参数空间。

图像分割滤波器处理图像，将其分割成（希望）有意义的区域。输出通常是整数的图像，其中每个整数可以表示一个对象。值0通常用于背景，值1（有时为255）用于前景对象。

%matplotlib inline
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatSlider
import SimpleITK as sitk
# Download data to work on
%run update_path_to_download_script
from downloaddata import fetch_data as fdata
from myshow import myshow, myshow3d

<Figure size 640x480 with 0 Axes>

img_T1 = sitk.ReadImage(fdata("nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT1.nrrd")
img_T2 = sitk.ReadImage(fdata("nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT2.nrrd"))

# 要可视化RGB中的标签图像，需要0-255范围的图像
img_T1_255 = sitk.Cast(sitk.RescaleIntensity(img_T1), sitk.sitkUInt8) # 将图像转换为sitkUInt8类型


img_T2_255 = sitk.Cast(sitk.RescaleIntensity(img_T2), sitk.sitkUInt8)

myshow3d(img_T1)

Fetching nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT1.nrrd

Fetching nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT2.nrrd

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-KPqemDXI-1680056431591)(SimpleITK笔记本中ITK细分.fld/image001.png)]

阈值

阈值是最基本的分割形式。它只是根据强度范围对图像的像素进行标记，而不考虑几何形状或连通性。

seg = img_T1 > 200 # 阈值分割
myshow(sitk.LabelOverlay(img_T1_255, seg), "Basic Thresholding") # sitk.LabelOverlay()函数将标签图像与原始图像叠加

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ama6eTt2-1680056431595)(SimpleITK笔记本中ITK细分.fld/image002.png)]

seg = sitk.BinaryThreshold(
img_T1, lowerThreshold=100, upperThreshold=400, insideValue=1, outsideValue=0) # 二值化阈值分割
myshow(sitk.LabelOverlay(img_T1_255, seg), "Binary Thresholding")

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VgWusIqB-1680056431596)(SimpleITK笔记本中ITK细分.fld/image003.png)]

ITK有许多基于直方图的自动阈值滤波器，包括Huang，MaximumEntropy，Triangle和流行的Otsu方法。这些方法创建一个直方图，然后使用启发式方法来确定阈值。

Otsu’s method: https://en.wikipedia.org/wiki/Otsu%27s_method 大津算法的基本思想是：对于给定的图像，计算图像中所有像素的灰度值的直方图，然后从中找到一个阈值，使得在该阈值下前景和背景的类间方差最大。

otsu_filter = sitk.OtsuThresholdImageFilter() # Otsu阈值分割
otsu_filter.SetInsideValue(0) # 设置内部值
otsu_filter.SetOutsideValue(1) # 设置外部值
seg = otsu_filter.Execute(img_T1) # 执行Otsu阈值分割
myshow(sitk.LabelOverlay(img_T1_255, seg), "Otsu Thresholding") # sitk.LabelOverlay()函数将标签图像与原始图像叠加
print(otsu_filter.GetThreshold()) # 打印阈值

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Fa6Gl2Uc-1680056431597)(SimpleITK笔记本中ITK细分.fld/image004.png)]

236.40869140625

区域生长分割(Region Growing Segmentation)

首先改进朴素阈值的第一步是一类称为区域生长的算法。这包括：

· 连接阈值(ConnectedThreshold)

· 连接置信(ConfidenceConnected)

· 向量置信(VectorConfidenceConnected)

· 邻域连接(NeighborhoodConnected)

我们之前使用3D_Slicer来确定索引：(132,142,96)是左侧侧腔的一个好种子。

seed = (132, 142, 96) # 设置种子点
seg = sitk.Image(img_T1.GetSize(), sitk.sitkUInt8) # 创建一个与img_T1相同大小的图像
seg.CopyInformation(img_T1) # 复制图像信息
seg[seed] = 1 # 设置种子点
seg = sitk.BinaryDilate(seg, [3] * 3) # 二值膨胀
myshow(sitk.LabelOverlay(img_T1_255, seg), "Initial Seed") # sitk.LabelOverlay()函数将标签图像与原始图像叠加

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-hZzwg9f6-1680056431598)(SimpleITK笔记本中ITK细分.fld/image005.png)]

seg = sitk.ConnectedThreshold(img_T1, seedList=[seed], lower=100, upper=190) # 连通阈值分割
myshow(sitk.LabelOverlay(img_T1_255, seg), "Connected Threshold") # 【译】标图叠加

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-cIDf6vef-1680056431599)(SimpleITK笔记本中ITK细分.fld/image006.png)]

更好的是，置信连接滤波器使用初始种子或当前分割来估计阈值范围。

seg = sitk.ConfidenceConnected(img_T1, seedList=[seed], numberOfIterations=1, multiplier=2.5, initialNeighborhoodRadius=1, replaceValue=1) # 置信度连通分割
myshow(sitk.LabelOverlay(img_T1_255, seg), "ConfidenceConnected") # 标图叠加

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-AzMUTIwj-1680056431601)(SimpleITK笔记本中ITK细分.fld/image007.png)]

img_multi = sitk.Compose(img_T1, img_T2) # 将两幅图像合并为一个多通道图像
seg = sitk.VectorConfidenceConnected(img_multi, seedList=[seed], numberOfIterations=1, multiplier=2.5, initialNeighborhoodRadius=1) # 多通道置信度连通分割
myshow(sitk.LabelOverlay(img_T1_255, seg), "VectorConfidenceConnected") # 标图叠加

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-gUlYxOw7-1680056431602)(SimpleITK笔记本中ITK细分.fld/image008.png)]

快进分割(Fast Marching Segmentation)

快进图像滤波器实现了一个快速行进的简单水平集演化问题（欧几里得方程）的解决方案。在这个例子中，微分方程中使用的速度项以图像的形式提供。速度图像基于梯度幅度，并映射为有界倒数1/(1+𝑥)。

seed = (132, 142, 96) # 设置种子点
feature_img = sitk.GradientMagnitudeRecursiveGaussian(img_T1, sigma=0.5) # 计算梯度图像
speed_img = sitk.BoundedReciprocal(feature_img) # 计算速度图像, 与Sigmoid不同，这是无参数的
myshow(speed_img) # 显示速度图像

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-5Q0ylyaZ-1680056431603)(SimpleITK笔记本中ITK细分.fld/image009.png)]

The output of the FastMarchingImageFilter is a time-crossing map that indicates, for each pixel, how much time it would take for the front to arrive at the pixel location.

快进图像滤波器的输出是一个时间交叉图(time-crossing map)，它指示每个像素，前端到达像素位置需要多长时间。

fm_filter = sitk.FastMarchingBaseImageFilter() # 创建FastMarchingBaseImageFilter对象
fm_filter.SetTrialPoints([seed]) # 设置试点
fm_filter.SetStoppingValue(1000) # 设置停止值
fm_img = fm_filter.Execute(speed_img) # 执行FastMarchingBaseImageFilter
myshow( sitk.Threshold( fm_img, lower=0.0, upper=fm_filter.GetStoppingValue(), outsideValue=fm_filter.GetStoppingValue() + 1, ) ) # 显示

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-QUkyCTCS-1680056431604)(SimpleITK笔记本中ITK细分.fld/image010.png)]

def fm_callback(img, time, z): # 定义回调函数
    seg = img < time # 二值化
    myshow(sitk.LabelOverlay(img_T1_255[:, :, z], seg[:, :, z])) # 显示
interact( lambda **kwargs: fm_callback(fm_img, **kwargs), time=FloatSlider(min=0.05, max=1000.0, step=0.05, value=100.0), z=(0, fm_img.GetSize()[2] - 1), ) # 交互式显示

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-AVzpheIC-1680056431605)(SimpleITK笔记本中ITK细分.fld/image011.png)]

<function main.(**kwargs)>

水平集分割(Level-Set Segmentation)

ITK中有各种基于水平集的分割滤波器：

· 地理学活动轮廓(GeodesicActiveContour)

· 形状检测(ShapeDetection)

· 阈值分割(ThresholdSegmentation)

· 拉普拉斯分割(LaplacianSegmentation)

· 标量Chan和Vese(ScalarChanAndVese)

另外还有一个模块化水平集框架(modular Level-set framework)，它允许组合术语并在C ++中轻松扩展。

首先，我们从种子创建一个标签图像。

seed = (132, 142, 96) # 设种子
seg = sitk.Image(img_T1.GetSize(), sitk.sitkUInt8) # 创建一个与img_T1相同大小的图像
seg.CopyInformation(img_T1) # 复制图像信息
seg[seed] = 1 # 种
seg = sitk.BinaryDilate(seg, [3] * 3) # 膨胀

使用种子来估计一个合理的阈值范围。

stats = sitk.LabelStatisticsImageFilter() # 创建labelStatisticsImageFilter对象
stats.Execute(img_T1, seg) # 执行LabelStatisticsImageFilter
factor = 3.5 # 设置阈值因子
lower_threshold = stats.GetMean(1) - factor * stats.GetSigma(1) # 计算下阈值
upper_threshold = stats.GetMean(1) + factor * stats.GetSigma(1) # 计算上阈值
print(lower_threshold, upper_threshold) # 打印阈值

>>> 81.25184541308809 175.0084466827569

init_ls = sitk.SignedMaurerDistanceMap(seg, insideIsPositive=True, useImageSpacing=True) # 计算距离图像
lsFilter = sitk.ThresholdSegmentationLevelSetImageFilter() # 创建ThresholdSegmentationLevelSetImageFilter对象
lsFilter.SetLowerThreshold(lower_threshold) # 设置下阈值
lsFilter.SetUpperThreshold(upper_threshold) # 设置上阈值
lsFilter.SetMaximumRMSError(0.02) # 设置最大均方根误差
lsFilter.SetNumberOfIterations(1000) # 设置迭代次数
lsFilter.SetCurvatureScaling(0.5) # 设置曲率缩放
lsFilter.SetPropagationScaling(1) # 设置传播缩放
lsFilter.ReverseExpansionDirectionOn() # 设置反向扩展方向
ls = lsFilter.Execute(init_ls, sitk.Cast(img_T1, sitk.sitkFloat32)) # 执行ThresholdSegmentationLevelSetImageFilter
print(lsFilter) # 打印ThresholdSegmentationLevelSetImageFilter对象

itk::simple::ThresholdSegmentationLevelSetImageFilter

 LowerThreshold: 81.2518

 UpperThreshold: 175.008

 MaximumRMSError: 0.02

 PropagationScaling: 1

 CurvatureScaling: 0.5

 NumberOfIterations: 1000

 ReverseExpansionDirection: 1

 ElapsedIterations: 119

 RMSChange: 0.0180966

 Debug: 0

 NumberOfThreads: 8

 NumberOfWorkUnits: 0

 Commands: (none)

 ProgressMeasurement: 0.119

 ActiveProcess: (none)

myshow(sitk.LabelOverlay(img_T1_255, ls > 0)) # 显示

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-uTG63BbU-1680056431607)(SimpleITK笔记本中ITK细分.fld/image012.png)]

.ipynb链接