码迷,mamicode.com
首页 > 其他好文 > 详细

胸部CT提取分割肺部

时间:2019-12-13 23:29:46      阅读:546      评论:0      收藏:0      [点我收藏+]

标签:width   生成   dsl   efi   算法   位置   range   body   measure   

1. 肺部分割提取简介

在处理胸部CT时,我们常常需要获取肺部的一个mask,也就是将肺部结构从数据中提取出来。二维图像还好说,但是三维图像就会变得复杂复杂一点。肺部的分割常常做后续操作的预处理,所以有必要提取提取一个肺部的mask,来辅助后面的操作,所以这里利用传统图像处理方法来提取了一下肺部,当时方法又很多,这里只是抛砖引玉,也许对有些数据不适用,可以对其进行改动。

2. 原理

利用阈值分割、种子填充图像形态学、图像连通域这些操作来进行肺部的分割。

2.1 阈值分割

这个好理解,一般来说CT值的范围是-1000-+1000,而基于简单的观察,肺部就是胸腔内两个大的空洞,所以可以首先对图像进行二值化处理,这里是将CT大于-300的位置置为1,小于-300置为0,这样就将数据分为了三部分,外部空气,内部空气,躯干组织。如下图所示。


技术图片

2.2 种子填充

利用种子填充算法,将外部的空气和内部的躯干分割出来,给定两个种子,一般就能分出来。分割的效果如下。然后再用阈值图像减去这个躯干就能得到初步的肺部mask。此时肺部的值是1,组织值是0。


技术图片

2.3 图像形态学

因为肺内部有许多纤维,所以看起来会有以下空洞(相对于肺部来说),要填补这些空洞,所以利用形态学里的闭操作(先膨胀,再腐蚀)。先膨胀肺部,将小的空洞填充,再腐蚀,恢复原来的大小。

2.4 连通域

最后保留最大的连通域,此时最大的连通域就是肺部。

2.5 效果

以下图片是使用3dslicer渲染出来的。

技术图片 技术图片 技术图片 技术图片 技术图片 技术图片 技术图片 技术图片

3.代码

依赖库

  • SimpleITK
  • skimage
import SimpleITK as sitk
from skimage import measure



def lungmask(vol):
    #获取体数据的尺寸
    size = sitk.Image(vol).GetSize()
    #获取体数据的空间尺寸
    spacing = sitk.Image(vol).GetSpacing()
    #将体数据转为numpy数组
    volarray = sitk.GetArrayFromImage(vol)

    #根据CT值,将数据二值化(一般来说-300以下是空气的CT值)
    volarray[volarray>=-300]=1
    volarray[volarray<=- 300]=0
    #生成阈值图像
    threshold = sitk.GetImageFromArray(volarray)
    threshold.SetSpacing(spacing)

    #利用种子生成算法,填充空气
    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)
    
    #用bodymask减去threshold,得到初步的lung的mask
    temp = sitk.GetImageFromArray(sitk.GetArrayFromImage(bodymask)-sitk.GetArrayFromImage(threshold))
    temp.SetSpacing(spacing)
    #利用形态学来去掉一定的肺部的小区域
    
    bm = sitk.BinaryMorphologicalClosingImageFilter()
    bm.SetKernelType(sitk.sitkBall)
    bm.SetKernelRadius(2)
    bm.SetForegroundValue(1)
    lungmask = bm.Execute(temp)
    
    #利用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)
    return l



def main():
    vol = sitk.ReadImage("Test.mha")
    volarray = sitk.GetArrayFromImage(vol)
    newvol = sitk.GetImageFromArray(volarray)
    newvol.SetSpacing(vol.GetSpacing())
    newvol.SetDirection(vol.GetDirection())
    newvol.SetOrigin(vol.GetOrigin())
    mask = lungmask(newvol)
    sitk.WriteImage(mask,"newlungmask.mha")

if __name__ == "__main__":
    main()

胸部CT提取分割肺部

标签:width   生成   dsl   efi   算法   位置   range   body   measure   

原文地址:https://www.cnblogs.com/WAoyu/p/12037576.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!