在后续训练模型的时候如果直接使用原始的肺部CT图像，其中的非肺区域会对训练造成干扰，因此我们需要提取分割肺实质。这一部分参考了kaggle上的教程。

从这张CT图像中可以看出扫描器以外的部分为黑色，背景和肺内部颜色较深，非肺组织和骨头颜色较浅。我们的目的是将肺部分割出来。由于肺部区域和其他部分对比明显，我们可以设置一个合适的阈值对图像中的区域进行隔离。

import argparse
import os
import numpy as np

from medpy.filter.smoothing import anisotropic_diffusion
from scipy.ndimage import median_filter
from skimage import measure, morphology
import scipy.ndimage as ndimage
from sklearn.cluster import KMeans

def segment_lung(img):
	# 1.标准化、去噪声
	# 2.重置像素值
	# 3.Kmeans聚类、分割
	# 4.腐蚀膨胀

一.使用到的库

• SciPy：是建立在NumPy扩展之上的数学算法和方便函数的集合。
  • SciPy.ndimage：这个包包含用于多维图像处理的各种函数。
• MedPy：是一个医学图像处理库。
  • medpy.filter这个包包含各种图像过滤器和图像处理函数。
• sklearn：是一个基于NumPy、SciPy和matplotlib构建的预测数据分析工具
  • sklearn.cluster包含流行的无监督聚类算法。
• Scikit-image：一个图像处理算法的集合。
  • morphology包含各种形态学处理
  • measure

二、标准化、去噪声

首先对图像进行减均值、除方差的标准化处理：

mean = np.mean(img)
std = np.std(img)
img = img-mean
img = img/std

我们必须确保我们在肺像素值和密集组织像素值之间设置了阈值。由于肺部像素比其他组织暗，背景区域也比其他组织暗，无法设置合适的阈值。为此，我们需要先将最小值的像素重置为图像中心附近的平均像素值，再执行k=2的kmeans聚类：

middle = img[100:400,100:400] 
mean = np.mean(middle)  
img[img==max]=mean
img[img==min]=mean

聚类前进行去噪，median_filter（中值滤波）方法将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值，可以去除椒盐噪声，是一种非线性平滑技术：

img= median_filter(img,size=3)

anisotropic_diffusion（各向异性扩散）旨在减少图像噪声而不去除图像内容的重要部分的技术，通常是边缘、线条或其他对解释图像很重要的细节，在本次实验中就是微小的结节：

img= anisotropic_diffusion(img)

三、KMeans聚类获得阈值

KMeans算法简单来说就是：

• 随机初始化簇的质心

• 分配步骤：将每个观测值分配给距离（欧几里得距离）最近的簇的质心。

• 更新步骤：重新计算分配到每个集群的观测值的平均值（质心）。回到分配步骤。当分配不再改变时，算法已经收敛。该算法不能保证找到最优值。

KMeans的fit方法的输入数据是通道顺序的，所以应该将输入图像展开（np.prod方法实现连乘）。找到两类像素的质心后，将阈值设置为平均值，低于阈值的部分为肺部区域，像素值置为1，其他区域置为0：

kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))
centers = sorted(kmeans.cluster_centers_.flatten())
threshold = np.mean(centers)
thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image

结果如下图：

四、腐蚀和膨胀

可以看到二值化的图像中存在很多小黑点，会对分割结果造成影响。为消除这些区域，同时保证肺部分割的完整性，可以最图像进行先腐蚀后膨胀的开运算。腐蚀的原理就是用一个像素点周围像素的最小值替换原值，膨胀操作相反：

eroded = morphology.erosion(thresh_img,np.ones([4,4]))
dilation = morphology.dilation(eroded,np.ones([10,10]))

结果如图：

如果使用闭运算，会导致边缘细节、左右肺部的间隔的损失较大：

目前，我们得到了一张黑白图像，数据类型是浮点型。我们使用measure.label方法标记整数数组中连接的区域。当两个像素相邻且具有相同值时，它们被连接起来，默认将像素值为0的区域作为背景。返回一个只包含0和1的数组。

接下来使用measure.regionprops测量标记图像区域的属性，方法返回一个region数组，region对象的bbox属性是区域边界框(min_row, min_col, max_row, max_col)。我们根据肺部大小的大致范围设置阈值，选出合适的肺部区域，然后与原图像相乘，就可以得到分割出肺实质的图像了。

labels = measure.label(dilation)
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)
mask = np.ndarray([512,512],dtype=np.int8)
mask[:] = 0
        
for N in good_labels:
    mask = mask + np.where(labels==N,1,0)
mask = morphology.dilation(mask,np.ones([10,10])) # 
return mask*img