基于动态DR的肺区域智能分割方法

陈卓 邓兴广 缪晓强 郑晗

1 广东省药品监督管理局审评认证中心 （广东 广州 510080）

2 深圳技术大学4D医学动态与功能成像联合实验室 （广东 深圳 518055）

3 深圳市安健科技股份有限公司 （广东 深圳 518055）

内容提要：多功能数字化动态数字X射线摄影（Digital Radiography，DR）具有高清影像、适时点片、实时回放、连续点片等多种功能，在临床诊断中得到广泛应用。动态DR可获得肺部呼吸运动变化的图像，对呼吸康复的效果测定、肺癌、慢性阻塞性肺疾病、支气管哮喘、间质性肺炎等肺部疾病的诊断起到有用的价值，还可以观察到肺部的血流分布和通气分布情况。文章基于动态DR全呼吸周期的图像，通过传统方法和深度学习的结合进行量化分析，包括肺野识别、左右肺划分、肺部面积计算、膈肌定位、实时输出肺面积变化曲线和膈肌运动曲线。为临床医生提供多项量化指标，并且为临床动态DR诊断提供辅助支持。

随着生活质量日益提高，人们更加重视自身的健康状况，对医疗资源的需求也逐渐增大。目前，大多数人体检采用的静态数字X射线摄影（Digital Radiography，DR）拍摄，动态DR拍摄尚未普及。动态DR有着巨大的潜力，通过对动态DR影像的分析可以观察到肺部血流和通气的分布情况，以及肺面积、横膈膜和肺密度的动态变化情况。本文分别采用传统方法和深度学习方法研究动态DR图像的肺分割以及左右横膈膜的获取，目标是达成对动态DR影像肺区域的实时分割以及得到左右肺面积和横膈膜位置的实时变化曲线。

1.国内外研究进展

1.1 肺分割

基于动态DR的肺分割，可通过单帧的肺分割加上帧与帧之间的平滑实现。为了精确地分割肺野，在过去的四十年里，人们已经发展出了许多方法[1]。根据Mittal等[1]的研究，他们对肺分割的方法进行了全面的调查，重点关注了它们的基本原理、所使用的数据集、分割的性能和相对的优缺点。Van Ginneken等[2]将肺野分割方法分为四大类：①基于规则的方法；②基于像素分类（PC）的方法；③基于可变形模型的方法；④混合方法。Van Ginneken等[2]调查的目的是对胸部图像计算机分割的文献进行分类和简要回顾，其中包括了过去30年发表的150多篇论文。

1.2 横膈膜运动分析

Baltruschat等[3]通过对172名志愿者的实验，使用动态DR对站立呼吸时膈肌的时间分辨率定量分析，结果表明，在站立式潮汐呼吸时，横膈肌的平均移动距离分别为11.0mm（右）和14.9mm（左）。左侧膈肌运动明显大于右侧。较高的BMI和潮气量与双侧横膈肌的短途活动增加有关。本文将基于此项研究继续探索横膈膜的运动问题。

2.动态DR图像的传统方法肺分割

本文围绕动态DR图像的肺分割，分别采用传统方法进行处理。包括基于阈值的肺分割的全阈值分割、最大类间方差法分割、自适应阈值法分割和基于非刚性配准的分割。

2.1 基于传统方法的肺分割

图像阈值分割是一种传统的最常用的图像分割方法，因其实现简单、计算量小、性能较稳定而成为图像分割中最基本和应用最广泛的分割技术。它特别适用于目标和背景占据不同灰度级范围的图像。它不仅可以大大压缩数据量，而且也大大简化了分析和处理步骤，因此，在很多情况下图像阈值分割是进行图像分析、特征提取与模式识别之前的必要的图像预处理过程。

图像阈值的目的是要按照灰度级，对像素集合进行一个划分，得到的每个子集形成一个与现实景物相对应的区域，各个区域内部具有一致的属性，而相邻区域不具有这种一致属性。这样的划分可以通过从灰度级出发选取一个或多个阈值来实现，从而实现目标与背景的分离。在医学图像上而言，组织器官的病灶的分割给临床诊断提供了一种有效的辅助手段，下文将分别介绍三种常见的基于阈值的肺分割算法。

2.1.1 全阈值分割

全阈值分割指将灰度值大于阈值的像素设为一种颜色，小于或等于阈值的像素设为另外一种颜色，在OpenCV库中可以实现全阈值分割。

用数学表达式来表示，则可设T为阈值，使得原始图像f(x,y)转换为二值图g(x,y)，其中x和y表示像素点的位置，全阈值分割图像时则满足公式（1）。

基于全阈值肺分割的总体方案如图1所示，并且每个步骤的结果如图2所示，并且制作了Matlab可操作交互式UI界面，界面中包含左右肺面积变化曲线的实时分割结果，UI界面如图3所示。

图1. 基于全阈值肺分割的总体方案

图2. 基于阈值的肺分割的分步结果

图3. 肺分割的Matlab UI界面

2.1.2 最大类间方差法分割

最大类间方差法（OTSU）也叫大津法，是一种确定图像二值化分割阈值的算法，由日本学者大津于1979年提出。从大津法的原理上来讲，该方法又称作最大类间方差法，因为按照大津法求得的阈值进行图像二值化分割后，前景与背景图像的类间方差最大。

OTSU是一种基于全局的二值化算法，它是根据图像的灰度特性，将图像分为前景和背景两个部分。当取最佳阈值时，两部分之间的差别应该是最大的，在OTSU算法中所采用的衡量差别的标准就是较为常见的最大类间方差。前景和背景之间的类间方差如果越大，就说明构成图像的两个部分之间的差别越大，当部分目标被错分为背景或部分背景被错分为目标，都会导致两部分差别变小，当所取阈值的分割使类间方差最大时就意味着错分概率最小。

记T为前景与背景的分割阈值，前景点数占图像比例为w0，平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1，图像的总平均灰度为u，前景和背景图像方差为g，具体计算见公式（2）～（4）。

由公式（2）和（3）可得公式（4）：

当方差g最大时，可以认为此时前景和背景差异最大，此时的灰度T是最佳阈值。类间方差法对噪声以及目标大小十分敏感，它仅对类间方差函数为单峰的图像产生较好的分割效果。当目标与背景的大小比例悬殊时（例如受光照不均、反光或背景复杂等因素影响），类间方差函数可能呈现双峰或多峰，此时效果不好。

在OpenCV中实现OTSU算法使用的是cv.threshold()函数，指定使用cv.THRESH_OTSU即可。

2.1.3 自适应阈值法分割

自适应阈值分割也叫做局部阈值化，它是根据像素的邻域块的像素值分布来确定该像素位置上的阈值。这样做的好处在于每个像素位置处的阈值不是固定不变的，而是由其周围邻域像素的分布来决定的。亮度较高的图像区域的阈值通常会较高，而亮度较低的图像区域的阈值则会相适应地变小。不同亮度、对比度、纹理的局部图像区域将会拥有相对应的局部阈值。

常用的局部自适应阈值有局部邻域块的均值和局部邻域块的高斯加权和。一种较为简单的自适应阈值选取方法是对每个像素确定以其自身为中心的一个窗口，寻找窗口内像素的最大值与最小值，并取二者的平均值作为阈值，或者将窗口内所有像素的平均值作为阈值，亦或者将窗口内的所有像素的高斯卷积作为阈值。

全局阈值化使用唯一的阈值，对于亮度分布差异较大的图像，难以找到一个合适的阈值。而使用自适应的阈值分割时，因为阈值是自适应的，所以可以将物体分割出来，效果要优于全阈值分割。

2.2 基于非刚性配准的肺分割

此方法是Sema Candemir等[4]提出来的一种巧妙的方法，使用到了非刚性配准技术，和预先带有人工分割标签的数据集。非刚性配准系统由以下三个阶段组成。

2.2.1 Ⅰ阶段

在训练数据库中识别5个与患者最相似的图像，并使用这个训练图像和训练图像的标签来扭曲出患者特定的肺边界。通过比较患者X射线与人工分割的肺图像，能显著改善之后图像配准的准确率。

首先使用Partial Radon Transforms，它是一个积分变换，将定义在二维平面上的一个函数f(x,y)沿着平面上的任意一条直线做线积分，来比较和排序两个患者的肺图像之间的相似性，Partial Radon Transforms如公式（5）所示。

注：I(·)是二维图像，δ(·)是二维脉冲函数，ρ=xcosθ+ysinθ，θ是法向量和x轴的夹角。

训练数据库中的所有图像预先计算水平和垂直投影轮廓。计算直方图在垂直和水平方向上的投影，然后用Bhattacharyya系数来测量训练数据库和患者胸部X射线之间的每个投影轮廓的相似性，计算相似性如公式（6）所示。

注：p1(x)、p2(x)是图像I1和I2的水平投影，q1(x)、q2(x)是图像I1和I2的垂直投影，x和y分别是投影配置文件的直方条，n和m是投影配置文件的直方条的数量，α=n/(n+m)是每个配置文件的相对比重。

2.2.2 Ⅱ阶段

选用SIFT特征提取算法进行特征点提取。该算法实现的主要步骤包括：检测尺度、计算空间极值点、确定特征点位置、提取特征点主方向、生成SIFT特征描述符。

在I阶段从数据库找到与患者最相似的数据，分别在数据库和患者的X射线图上找出SIFT特征点，通过对应匹配之间的空间距离移动所有像素，将变换映射应用于所有像素，在人工分割标签的基础上扭曲出患者特定的掩码图。图像扭曲获取掩码图的公式如（7）所示。

注：P是X射线的像素集，N是空间邻域集（Spatial Neighborhood Set），S1、S2是由SIFT描述符向量表示的SIFT图像，是关于P的光流向量（Flow Vectors），t和d是截断的阈值。

2.2.3 Ⅲ阶段

该系统利用图像特性和上一阶段计算出的肺模型概率，来检测X射线图像的肺边界。使用图割算法进行图像分割，并用目标函数建模分割过程。

综上所述，本文通过借鉴Sema Candemir等[4]使用非刚性配准的数据集在胸片中进行肺分割，将此技术运用到动态DR图像中的肺分割中，完成了现有十四套DR数据肺区域的实时分割。但是由于数据集较小，鲁棒性不是特别好，并且分割速度较慢。

部分结果的UI界面如图4所示，界面中有四个坐标图，左上表为左肺面积变化曲线、右上表为右肺面积变化曲线、左下角为左横膈膜位移变化曲线、右下角为右横膈膜位移变化曲线。

图4. 基于非刚性配准的肺分割结果

根据分割的结果来看，部分数据出现突变点或者分割效果不稳定，这与标签的数量有关，标签的多样性和丰富性对结果有影响，后期可以通过添加数据集中的样本数量和多样化来改善分割效果。

基于非刚性配准肺分割的结果，可以提取到横膈膜的位移信息。以第一帧图像横膈膜的最高点为基准，通过固定第一帧的x轴，来确定后续帧的横膈膜位置，位置结果如图5所示。

图5. 横膈膜位置的自动选取

3.基于深度学习的肺分割

3.1 环境搭建

本文采用基于U-Net改编的网络进行实验，具体实验平台的相关配置如表1所示。

表1. 实验平台的相关配置

PyTorch是一个开源的Python机器学习库，用于自然语言处理、深度学习等领域。2017年1月，由Facebook人工智能研究院基于Torch推出了PyTorch。它是一个基于Python的可续计算包，具有强大的GPU加速的张量计算（如NumPy），并且包含自动求导系统的深度神经网络。

计算统一设备体系结构（Compute Unified Device Architecture，CUDA），是显卡厂商NVIDIA推出的运算平台。CUDA是一种通用并行计算架构，该架构使GPU能够解决复杂的计算问题。

cuDNN是用于深度神经网络的GPU加速库。它强调性能、易用性和低内存开销。NVIDIA cuDNN可以集成到更高级别的机器学习框架中，如Tensor flow和Caffe。简单的插入式设计可以让开发人员专注于设计和实现神经网络模型，而不是简单调整性能，同时还可以在GPU上实现高性能现代并行计算。

3.2 数据准备

本文采用的数据集来自Kaggle网站上关于深圳医院的肺部DR图像，数据集包含已经分割好的标签，为了提高模型的鲁棒性和泛化能力，本研究在现有的样本上进行数据增强。

数据增强也叫数据扩增，即在不实质性的增加数据的情况下，让有限的数据产生等价于更多数据的价值。数据增强可以分为有监督的数据增强和无监督的数据增强方法。其中，有监督的数据增强又可以分为单样本数据增强和多样本数据增强方法，无监督的数据增强分为生成新的数据和学习增强策略两个方向。

本文采用的数据增强的方法包括随机旋转、随机裁剪和均衡直方图化。通过调用PyTorch框架里的torchvision.transforms来实现，然后对数据增强后的图像进行训练。

数据增强的本质是为了增强模型的泛化能力，它与其他的一些方法比如DropOut，权重衰减有一定的区别。①权重衰减、DropOut等方法是专门设计来限制模型的有效容量的，用于减少过拟合，这一类是显式的正则化方法。研究表明，这一类方法可以提高泛化能力，但并非必要且能力有限，同时参数高度依赖于网络结构。②数据增强，没有降低网络的容量，也不增加计算复杂度和调参工程量，是隐式的规整化方法，在实际应用中更具有意义。

3.3 BD-U-Net网络结构

本文模型基于U-Net框架在卷积块中加入了批标准化（Batch Normalization）和DropOut来提高模型的泛化能力。U-Net模型是菲兹堡大学Olaf Ronneberger等[5]在2015年参加ISBI竞赛时，提出的一种全卷积神经网络，最初用于完成细胞分割任务。U-Net模型结构优雅精巧，对于数据较少的医学图像，端到端的分割任务效果较好。

Olaf Ronneberger等[5]提出的U-Net的U形结构如图2～6所示。网络是一个经典的全卷积网络（即网络中没有全连接操作）。网络的输入是一张572×572的边缘经过镜像操作的图片。

图6. U-Net网络结构[10]

网络的左侧在论文中将这一部分叫做压缩路径。压缩路径由4个块组成，每个块使用了2个有效卷积和1个Max Pooling降采样，每次降采样之后Feature Map的个数乘2，因此有了图中所示的Feature Map尺寸变化。最终得到了尺寸为32×32的Feature Map。

网络的右侧部分在论文中叫做扩展路径。同样由4个块组成，每个块开始之前通过反卷积将Feature Map的尺寸乘2，同时将其个数减半（最后一层略有不同），然后和左侧对称的压缩路径的Feature Map合并，由于左侧压缩路径和右侧扩展路径的Feature Map的尺寸不一样，U-Net是通过将压缩路径的Feature Map裁剪到和扩展路径相同尺寸的Feature Map进行归一化的（即图2～6中左侧虚线部分）。扩展路径的卷积操作依旧使用的是有效卷积操作，最终得到的Feature Map的尺寸是388×388。由于该任务是一个二分类任务，所以网络有两个输出Feature Map。

本文在卷积过后增加Batch Normalization和dropout操作，以此提高模型的泛化能力，并且防止模型过拟合。其中为了简洁表格，定义了两个模块，分别为Conv_BN_DP和ConvT_BN_DP，如表2所示。修改过的模型网络结构及每层的参数和每层输出图像尺寸如表3所示，并且BD-U-Net网络的总计算量如表4所示。

表2. BD-U-Net模块

表3. BD-U-Net 结构

表4. BD-U-Net总参数

Batch Normalization是Sergey Ioffe等[6]在2015年研发的，普通的归一化操作（Z-score标准化）归一化到均值为0、方差为1的公式如下，即每个样本减去样本数据的平均值后再除以样本数据的标准差，如公式（8）所示。

注：ε是个很小的数，防止分母为0，γ和β是可训练的参数。

3.4 训练与验证

本文采用交叉熵损失函数，Golik等[7]提出在随机初始化权重的情况下，基于平方误差的神经网络不能收敛到一个好的局部最优解，而交叉熵更加快速地使网络收敛。

学习率是一个至关重要的参数，影响模型的收敛速度，本文使用动态学习率，通过调用torch.optim.lrscheduler.StepLR，可以实现每30轮训练学习率递减为原来的10%，本文采用的训练参数如表5所示。

表5. 训练参数

本研究制作的数据集包含700张肺部图片，将其分为训练集：验证集：测试集=7:2:1。每训练完一轮后，就进行一次模型效果的验证，全部轮数训练完后，再用测试集中的数据进行测试。训练损失与训练轮数的关系如图7所示，验证准确率与训练轮数的关系如图8所示。

图7. 训练损失与训练轮数的关系图

图8. 验证准确率与训练轮数的关系

3.5 测试结果

将测试数据输入模型后得到肺的掩码图，通过图像形态学腐蚀、膨胀的操作，可以光滑肺边缘和去除小连通域，也可以通过掩码图减去膨胀的掩码图，得到肺掩码图的轮廓线。肺部分割效果如图9所示。

图9. 肺部分割效果

通过计算每一帧肺部图像的肺面积和肺密度，就可以的到肺面积变化曲线和肺密度变化曲线，曲线的形态以及动态区间，可以辅助医生进行诊断。动态曲线如图10所示。

图10. 肺部动态曲线

综上所述，通过传统方法和深度学习的肺分割结果来看，深度学习的速度比传统方法快，其鲁棒性比传统方法好，其准确率比传统方法高。

4.小结

医学影像是辅助临床疾病诊疗重要手段，其中DR是最为普及的医疗仪器之一，动态DR在此基础上额外提供了不同时间点的动态信息。深度学习的模型对复杂图像有高效的处理能力，在肺区域分割过程中，能够在多帧DR处理时提高效率，并且准确率和鲁棒性都很好。本研究实现了基于动态DR影像的实时左右肺分割以及横膈膜提取，输出了动态肺面积曲线，以及横膈膜运动曲线。此曲线能够辅助诊断患者的肺容量，可以辅助排查出慢性阻塞性肺疾病等肺部疾病，提升临床诊疗质量以及效率。