胸部CT提取分割肺部

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