LUNA16肺实质分割，理论+代码

一、LUNA16数据集介绍

为了验证模型的检测效果，本文在LUNA16数据集上进行了测试。LUNA16数据集是LIDC-IDRI的子集，LIDC-IDRI为肺结节检测的公共数据集，LUNA16 数 据 集 删 除 了LIDC-IDRI中肺结节小于３ｍｍ 和切片厚度大于３ｍｍ 的CT图像，最后保留了888个低剂量肺部CT图像，以及由４名影像学专家手工标记的1186个结节标签，结节直径大小在３～30ｍｍ之间。

二、肺部预处理步骤

肺部预处理通常等于肺实质分割，主要是为了消除背景因素的干扰从而提高整个肺结节检测系统的性能，得到精细的肺部图像。

注：每个代码块的解释均在代码块的下方。每个代码块的下方都是该代码的宏观解释。

1.读取原始图像（raw格式），并将原始图像转换为numpy格本

import os
import SimpleITK as sitk
import numpy as np
from glob import glob
import pandas as pd
try:
    from tqdm import tqdm # long waits are not fun
except:
    print('TQDM does make much nicer wait bars...')
    tqdm = lambda x: x

def normalizePlanes(npzarray):  #对原图像进行归一化
    maxHU = 400
    minHU = -1000
    npzarray = (npzarray - minHU)/(maxHU - minHU)
    npzarray[npzarray>1] = 1
    npzarray[npzarray<0] = 0
    npzarray *= 255

    return (npzarray.astype(int))

引入HU值将肺部与其他组织隔开，HU值是用于计算人体中特定的组织和器官密度的测量单位。在低剂量螺旋CT图像中，肺部的密度和周围的器官有一定的差异，因此可以通过这种差异来进行肺实质分割。常见的HU值如图1所示。将HU值限制在[-1000,400]中，HU值<-1000人体组织表现为黑色，CT>400的表现为白色。这样就将肺部组织筛选出来。接着对像素值归一化到[0,255]

图1 常见组织的HU

luna_path = "F:/lungsdatasets/LUNA16/subset1/"         #输入路径
output_path = "E:/luna16/lungnpyset1/"        #输出路径
file_list = glob("F:/lungsdatasets/LUNA16/subset1/" + "*.mhd")



#获取数据的每一行
def get_filename(case):
    global file_list
    for f in file_list:
        if case in f:
            return (f)

df_node = pd.read_csv('F:/lungsdatasets/LUNA16/annotations.csv')
df_node["file"] = df_node["seriesuid"].apply(get_filename)   #匹配mhd名字与annotation中seriesuid列是否相同，并创建一列取名file，file中是mhd所在路径
# df_node: 文件的路径  example：E:\LUNA16\src\subset0\1.3.6.1.4.1.14519.5.2.1……
df_node = df_node.dropna()   #删除nan值

这部分是为了挑选出只含结节的CT图像，我们只处理含有结节的肺部CT图像

for fcount, img_file in enumerate(tqdm(file_list)):  #fcount计数器，一次循环计数一次， img_file就是file_list
    print("Getting mask for image file %s" % img_file.replace(luna_path, ""))
    mini_df = df_node[df_node["file"] == img_file]  # 得到所有结节
    if(len(mini_df) > 0):    #跳过没有结节的文件
      
        #读取数据
        itk_img = sitk.ReadImage(img_file)
        img_array = sitk.GetArrayFromImage(itk_img)   #(z,y,x)
        num_z,height,width = img_array.shape    #heightXwidth constitute the transverse plane  
        origin = np.array(itk_img.GetOrigin())  #世界坐标系下的x,y,z(mm)
        spacing = np.array(itk_img.GetSpacing()) #世界坐标中的体素间隔(mm)

这部分通过SimpleITK库来读取raw格式的CT图像，获取图像的宽度，高度，切片数量等

 #遍历所有节点
        for node_idx,cur_row in mini_df.iterrows():
            node_x = cur_row["coordX"]
            node_y = cur_row["coordY"]
            node_z = cur_row["coordZ"]
            diam = cur_row["diameter_mm"]

            #保留三个切片
            imgs = np.ndarray([3,height,width],dtype=np.float32)#保留3张切片
       
            center = np.array([node_x,node_y,node_z])  #结点中心
            v_center = np.rint((center-origin)/spacing) #体素坐标系的结点中心(x,y,z)
            for i,i_z in enumerate(np.arange(int(v_center[2])-1,int(v_center[2])+2).clip(0,num_z-1)):   #clip防止超出z
                imgs[i] = normalizePlanes(img_array[i_z])
            
            np.save(os.path.join(output_path,"images_%04d_%04d.npy" % (fcount, node_idx)),imgs)
          

这部分代码保留每个样本的3张切片，并将读取的原始图像进行归一化操作，保存为npy格式。

2.读取归一化后的原始图形（npy格式），分割出肺部掩码，下图为生成肺部掩膜流程图

from matplotlib import pyplot as plt
from skimage import morphology
from skimage import measure
from sklearn.cluster import KMeans
from glob import glob
import numpy as np

working_path = "E:/luna16/lungnpy/"
file_list = glob(working_path + "images_*.npy")

for img_file in file_list:
    imgs_to_process = np.load(img_file).astype(np.float64) 
    print("on image",img_file)
    for i in range(len(imgs_to_process)):
        img = imgs_to_process[i]

        #标准化
        mean = np.mean(img)   #均值
        std = np.std(img)    #标准差
        img = (img-mean)/std

        #寻找肺部附近的平均像素,以重新调整过度曝光的图像
        middle = img[100:400,100:400]

这部分代码将处理好的图像进行标准化处理，将值限定在[0,1]之间

 #使用Kmeans算法将前景（放射性不透明组织）和背景（放射性透明组织，即肺部）分离。
        #仅在图像中心进行此操作，以尽可能避免图像的非组织部分。
        kmeans = KMeans(n_clusters=2,n_init=10).fit(np.reshape(middle,[np.prod(middle.shape),1]))
        #np.prod 计算给定数组中所有元素的乘积
        #np.reshape 低维变高维  .flatten() 高维变1维
        centers = sorted(kmeans.cluster_centers_.flatten())
        threshold = np.mean(centers)
        thresh_img = np.where(img<threshold,1.0,0.0)  #阈值化图像，二值化处理

        #腐蚀和膨胀
        eroded = morphology.erosion(thresh_img,np.ones([4,4]))
        dilation = morphology.dilation(eroded,np.ones([10,10]))

这部分代码使用 k-means算法对其进行聚类分析，将图像中心部分区域分成前景和背景；再分别生成前景和背景两个中心点的值，对两个中心点取平均值，设为阈值Ｔ，对标准化后的CT切片图片 进行遍历，如果像素大于Ｔ则属于前景，设 为１，小 于Ｔ则属于背景，设为０；再利用腐蚀、膨胀操作将图像上的孤立点消除，腐蚀减去毛刺，膨胀填补孔洞。

labels = measure.label(dilation)#对二值图像进行标记，标记连通区域
        label_vals = np.unique(labels)  #获取标记值的唯一值，即标记的数量
        regions = measure.regionprops(labels) #获取连通区域的属性
        good_labels = []
        for prop in regions:
            B = prop.bbox   #边界框
            if B[2]-B[0]<475 and B[3]-B[1]<475 and B[0]>40 and B[2]<472:
                good_labels.append(prop.label)   #边界框满足条件的加入到good_labels框内，用作下面的计算
            #np.ndarray 创建多维数组
        mask = np.ndarray([512,512],dtype=np.int8)
        mask[:] = 0  #肺部mask

        for N in good_labels:
            mask = mask + np.where(labels==N,1,0)
        mask = morphology.dilation(mask,np.ones([10,10]))
        imgs_to_process[i] = mask

这部分代码对二值化图像提取连通区域，像素值相连且相同为连通区域，将连通区域打上标签，接着提取连通区域的属性，接着将根据边界框的阈值对提取的连通 区域进行挑选，调出我们所要的左右肺连通区域，接着，将图片背景全部为黑，接着将左右肺也就是我们挑选好的连通区域，赋值为1到labels标签，也就是经过腐蚀与膨胀后的二值图像 现在，就获得一张背景为黑，左右肺为白色的图片，最后再次进行膨胀操作，完成对肺的分割。

3.将第2步分割好的肺部掩码与图像与第1步得到的归一化CT切片进行逐像素相乘，得到感兴趣区域图像

from glob import glob
from skimage import measure
from skimage.transform import resize
import numpy as np

working_path = "E:/luna16/lungnpy/"

lungmask_list = glob(working_path + "lungmask_*.npy")
out_images = []  #肺部掩膜列表  共336（112*3）
#out_nodemasks = []  #结节掩膜列表
for fname in lungmask_list:
    #print("working on file ",fname)
    imgs_to_process = np.load(fname.replace("lungmask","images"))  #原始图像
    masks = np.load(fname)                                         #肺部掩膜
    #node_masks = np.load(fname.replace("lungmask","masks"))        #结节掩膜
    for i in range(len(imgs_to_process)):
        mask = masks[i]  #肺部掩膜
        #node_mask = node_masks[i]   #结节掩膜
        img = imgs_to_process[i]   #原始图像
        #new_size = [512,512]
        img = mask*img   #肺部掩码x原始图像
        #ROI区域归一化
        new_mean = np.mean(img[mask>0]) #取提取部分平均值
        new_std  = np.std(img[mask>0])  #取提取部分标准差
        old_min = np.min(img)   #提取元素最小值
        img[img == old_min] = new_mean-1.2*new_std  #更新元素最小值，防止最小值受噪音影响
        img = (img-new_mean)/(new_std)  #归一化

这部分代码加载归一化后的原始图像与分割好的肺部掩膜，并进行相乘。接着将相乘后的感兴趣区域进行归一化。

        #创建图像的边界框
        labels = measure.label(mask)   #提取连通区域
        regions = measure.regionprops(labels)   #提取连通区域的属性

        #找所有区域中的最值边界
        min_row = 512
        max_row = 0
        min_col = 512
        max_col = 0
        for prop in regions:
            B = prop.bbox  #连通区域的边界框
            if min_row > B[0]:
                min_row = B[0]
            if min_col > B[1]:
                min_col = B[1]
            if max_row < B[2]:
                max_row = B[2]
            if max_col < B[3]:
                max_col = B[3]
        height = max_row-min_row
        width = max_col-min_col
        if width > height:
            max_row = min_row + width
        else:
            max_col = min_col + height
         img = img[min_row:max_row,min_col:max_col]   #原始图像与肺部掩膜相乘
        mask = mask[min_row:max_row,min_col:max_col]  #肺部掩膜
#        print("num_row: ",max_row-min_row)
#        print("num_col: ",max_col-min_col)
        if max_row-min_row < 5 or max_col-min_col < 5:
            pass
        else:
#            print(2)
            mean = np.mean(img)
            max = np.max(img)
            min = np.min(img)
            img = (img-mean)/(max-min)
            new_img = resize(img,[512,512])
            #new_node_mask = resize(node_mask[min_row:max_row,min_col:max_col],[512,512])
            out_images.append(new_img)
#            print(len(out_images))
            #out_nodemasks.append(new_node_mask)

 这部分代码是对感兴趣区域进行裁剪

num_images = len(out_images)

#将图像和掩膜变成单通道，便于训练
final_images = np.ndarray([num_images,1,512,512],dtype=np.float32)
#final_masks = np.ndarray([num_images,1,512,512],dtype=np.float32)
for i in range(num_images):
    final_images[i,0] = out_images[i]
   #final_masks[i,0] = out_nodemasks[i]

rand_i = np.random.choice(range(num_images),size=num_images,replace=False).astype(int)  #随机生成图形索引，为生成训练与测试集
test_i = int(0.2*num_images)   #测试集占总数据集20%

np.save(working_path+"trainImages.npy",final_images[rand_i[test_i:]])
#np.save(working_path+"trainMasks.npy",final_masks[rand_i[test_i:]])
np.save(working_path+"testImages.npy",final_images[rand_i[:test_i]])
#np.save(working_path+"testMasks.npy",final_masks[rand_i[:test_i]])

这部分代码将图像改为单通道，并生成训练集与测试集，方便我们进行训练。

from glob import glob
#import matplotlib as plt
import matplotlib.pyplot as plt
import numpy as np
from skimage import measure
import os
path = "E:/luna16/lungnpyset0/"

trainImage = np.load(path + "trainImages.npy")

output_folder = "E:/luna16/11/"

for i, img in enumerate(trainImage):
    for j, slice_img in enumerate(img):
        filename = f"image_{i}_{j}.png"
        output_path = os.path.join(output_folder,filename)
        plt.imsave(output_path,slice_img,cmap="gray")

这部分代码将我们在分割好的感兴趣区域保存为png格式到指定文件夹（这里只展示训练集）。

这是分割好的成果。

三、总结

肺实质分割是肺结节检测或分类的前期准备工作，这部分工作的好坏对于肺结节检测或分类的精度有着较大影响。实际上实质分割主要步骤为：归一化原始图像-生成肺部掩膜-将肺部掩膜与原始图像进行相乘。以上方法均为为形态学分割方法，分割效果本人认为不太理想，需要在此基础上进行改进。