【SimpleITK】使用区域生长法/种子填充法/形态学分割肺-CT img

区域生长

定义参考 维基百科 :
对图像分割的定义：

图像分割是对图像中的每个像素加标签的一个过程，这一过程使得具有相同标签的像素具有某种共同视觉特性。

区域生长的定义：

区域生长的基本思想是将具有相似性质的像素集合起来构成区域。具体先对每个需要分割的区域找一个种子像素作为生长的起点，然后将种子像素周围邻域中与种子像素具有相同或相似性质的像素（根据某种事先确定的生长或相似准则来判定）合并到种子像素所在的区域中。将这些新像素当做新的种子像素继续进行上面的过程，直到再没有满足条件的像素可被包括进来，这样，一个区域就长成了。

定义简单明了。

由定义可知，区域生长的三个要点：

1. 种子点
2. 生长准则
3. 终止条件

根据不同的生长准则和终止条件有不同的区域生长算法。
本次使用的是Confidence Connected。

Confidence Connected

参考官方文档： http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/Python_html/30_Segmentation_Region_Growing.html

This region growing algorithm allows the user to implicitly specify the threshold bounds based on the statistics estimated from the seed points, 𝜇±𝑐𝜎 . This algorithm has some flexibility which you should familiarize yourself with:
参数 ：

• multiplier ： the constant 𝑐 from the formula above.
  You can specify a region around each seed point
• initialNeighborhoodRadius ： he statistics are estimated, see what happens when you set it to zero.
• numberOfIterations ： to rerun the algorithm.
  • In the first run the bounds are defined by the seed voxels you specified, in the following iterations 𝜇 and 𝜎 are estimated from the segmented points and the region growing is updated accordingly.

肺窗

肺窗WW1500—2000HU 、WL-450—-600HU

代码

导入simpleITK:

import SimpleITK as sitk
import matplotlib.pyplot as plt 
WINDOW_LEVEL = (1000,-500)

导入数据：

img_path = “...”
img = sitk.ReadImage(img_path)
img.GetSize()

(512, 512, 121)
查看其中一层：

img_arr = sitk.GetArrayFromImage(img)
plt.imshow(img_arr[68])

调整图像到肺窗！

img_norm = sitk.Cast(sitk.IntensityWindowing(img,
                                  windowMinimum=WINDOW_LEVEL[1] - WINDOW_LEVEL[0] / 2.0,
                                  windowMaximum=WINDOW_LEVEL[1] + WINDOW_LEVEL[0] / 2.0),
                                  sitk.sitkUInt8)

显示下肺窗下的图像：

img_norm_arr = sitk.GetArrayFromImage(img_norm)
plt.imshow(img_norm_arr[68], cmap='bone')

使用区域生长法分割：

seed_pts = [(108, 246, 68), (390, 255, 68)]
img_grow_Confidence = sitk.ConfidenceConnected(img, seedList=seed_pts,
                                                           numberOfIterations=0,
                                                           multiplier=2,
                                                           initialNeighborhoodRadius=1,
                                                           replaceValue=1)

打印结果：

img_grow_arr_Confidence = sitk.GetArrayFromImage(img_grow_Confidence)
plt.imshow(img_grow_arr_Confidence[68])

结果上有很多的空洞需要填补。
对结果进行闭包处理：

• BinaryMorphologicalClosingImageFilter

BMC = sitk.BinaryMorphologicalClosingImageFilter()
BMC.SetKernelType(sitk.sitkBall)
BMC.SetKernelRadius(2)
BMC.SetForegroundValue(1)
OUT = BMC.Execute(img_grow_Confidence)

显示下闭包之后的结果：

ov_close = sitk.LabelOverlay(img_norm, OUT)
plt.imshow(sitk.GetArrayFromImage(ov_close)[68])

此方法的唯一缺点是需要手动 定种子点。
更好的方法在如下链接：

https://www.cnblogs.com/WAoyu/p/12037576.html

改进方法

使用：

• 阈值分割
• 种子填充
• 形态学闭包
• 图像联通域

1. 阈值分割

vol = img
size = vol.GetSize()
print('size :', size)
spacing = vol.GetSpacing()
print('spacing :', spacing)
volarray = sitk.GetArrayFromImage(vol)

size : (512, 512, 121)
spacing : (0.5625, 0.5625, 2.5)

volarray[volarray>=-300]=1
volarray[volarray<=- 300]=0
#生成阈值图像
threshold = sitk.GetImageFromArray(volarray)
threshold.SetSpacing(spacing)

显示阈值分割之后的图像：

2. 种子生成法填充空气

#利用种子生成算法，填充空气
ConnectedThresholdImageFilter = sitk.ConnectedThresholdImageFilter()
ConnectedThresholdImageFilter.SetLower(0)
ConnectedThresholdImageFilter.SetUpper(0)
ConnectedThresholdImageFilter.SetSeedList([(0,0,0),(size[0]-1,size[1]-1,0)])
    
#得到body的mask，此时body部分是0，所以反转一下
bodymask = ConnectedThresholdImageFilter.Execute(threshold)

显示结果：

像素反转：

bodymask = sitk.ShiftScale(bodymask,-1,-1)

显示结果：

3. 相减的肺mask

#用bodymask减去threshold，得到初步的lung的mask
temp = sitk.GetImageFromArray(sitk.GetArrayFromImage(bodymask)-sitk.GetArrayFromImage(threshold))
temp.SetSpacing(spacing)
temp_arr = sitk.GetArrayFromImage(temp)
plt.imshow(temp_arr[68], cmap='bone')

显示结果：

4. 对结果形态学闭包

#利用形态学来去掉一定的肺部的小区域   
bm = sitk.BinaryMorphologicalClosingImageFilter()
bm.SetKernelType(sitk.sitkBall)
bm.SetKernelRadius(2)
bm.SetForegroundValue(1)
lungmask = bm.Execute(temp)

显示结果：

5. 保留最大联通域

#利用measure来计算连通域
lungmaskarray = sitk.GetArrayFromImage(lungmask)
label = measure.label(lungmaskarray,connectivity=2)
props = measure.regionprops(label)

#计算每个连通域的体素的个数
numPix = []
for ia in range(len(props)):
    numPix += [props[ia].area]

#最大连通域的体素个数，也就是肺部
maxnum = max(numPix)
#遍历每个连通区域
for i in range(len(numPix)):
    #如果当前连通区域不是最大值所在的区域，则当前区域的值全部置为0，否则为1
    if numPix[i]!=maxnum:
        label[label==i+1]=0
    else:
        label[label==i+1]=1

label = label.astype("int16")
l = sitk.GetImageFromArray(label)
l.SetSpacing(spacing)

显示结果：