肺组织包括肺实质和肺间质, 其中肺实质由肺内支气管的各级分支、血管及可能存在的肺内病变等组织构成.肺结节是肺实质内病发性的圆形或类圆形病灶.对肺结节的准确分割是提升肺CAD诊断准确率的重要前提, 目前已有诸多研究成果^[1-4], 但由于肺结节本身结构和特征的多样性, 仍存在一些问题难以解决:

1) 肺结节的CT影像形状多样且结构复杂, 大大提高了肺结节的分割难度.例如粘连型结节与周围的正常组织灰度相近, 极易发生误判和漏检.磨玻璃型结节具有密度小、显示性状模糊的特点, 检测难度大且恶性结节居多.

2) 为提升准确度, 现有的肺结节分割算法大多需要增加人工干预, 且分割目标多选取孤立肺结节, 适用范围小, 对粘连型、磨玻璃型肺结节的分割效果较差, 且此类结节的恶性概率要显著高于孤立型肺结节.

3) 现有算法多针对二维肺结节进行分割, 忽略了肺结节的空间信息, 而直接进行三维分割又普遍存在算法复杂度过高、效率低等问题.

1 研究方法

按照形状、位置等特征, 可将肺结节大致分为孤立肺结节、磨玻璃型肺结节、胸腔粘连型肺结节和血管粘连型肺结节.各类结节CT影像如图 1所示.

本文提出一种结合空间约束的图割方法来实现三维肺实质的自动分割, 旨在分割并提取肺实质内所有支气管、血管和肺结节等组织, 作为后续肺结节分类工作的必要前提.具体工作流程如图 2所示.

2 预处理 2.1 超像素构建

传统的图割方法在分割简单的低分辨率图像时具有良好的实时交互性, 但是对于高分辨率图像、多维图像或图像序列, 无论内存开销还是计算复杂度, 都严重限制了图割方法的应用.因此本文采用基于超像素构建结果构造加权图的策略来简化和加速图割模型.

超像素构建的基本原理是利用像素之间特征的相似性将图像分割成多个像素块.以超像素作为图像的基础单位可有效降低图像规模和后续处理的复杂度.从2003年Ren等^[5]首次提出超像素概念至今, 已有诸多关于超像素生成算法的研究.考虑到模型复杂度、运行时间等问题, 采用了Achanta等^[6]提出的简单线性迭代聚类(simple linear iterative cluster, SLIC)算法.相比于其他算法, SLIC在生成超像素紧凑度、运行速度和轮廓保持方面都有明显优势.

SLIC算法将图像彩色空间转化为CIELAB颜色空间, 图像中每个像素的颜色值和二维坐标组成相应的5维特征向量C_i=(l_i, a_i, b_i, x_i, y_i), 计算特征向量距离以获得聚类中心和像素间的相似性.迭代更新聚类中心直至误差收敛, 最终实现对图像像素的局部聚类.实验发现迭代10次后多数图片可以获得较好的分割结果, 如图 3所示.

2.2 空间图结构的定义

构建S-T网络, 即将图像转化为图结构, 是图割算法的首要步骤.为了充分利用肺部组织的空间信息, 在传统图结构的基础上加以改进, 将超像素视作构成图的节点, 通过连接相邻切片中的对应节点来构建肺实质的空间拓扑结构.

对于每一张超像素图像, 用G= < V, E>来表示图的构成.其中, V表示图的点集, E表示连接相邻点的边集, 图上节点与图像中的超像素一一对应.连接相邻超像素的边称为n边(n-links).除普通节点外, 增加两个终端节点:源点S和汇点T, 将连接普通节点和终端节点的边记作t边(t-links).

在此基础上, 本文新增连接相邻切片中超像素对的边, 记作n_z.限定空间距离范围为K以减小图像尺寸对计算量的影响, 连接图中每个节点z轴方向上K邻域内的相邻点, 从而构造肺实质的空间结构.对生成的子空间结构使用最小生成树进行约束, 并将每个最小生成树合并成一个节点以实现数据降维, 合并过程如图 4所示.

将每棵树视作一个小集合, 每个小集合∈ 的结构可以用一个星型结构表示.连接树和树之间边的约束极弱, 此处忽略不计, 因此树与树之间的相似性计算可以转化成计算成对的相似度问题^[7].

每个集合视作一个超像素节点, 对于每个节点, 计算该节点区域的归一化灰度直方图P(r_k):

(1)

式中:r_k表示图像中灰度级为k的像素数量; n为超像素中所含像素总数.

对于同一层切片中的超像素对, 计算P(r_k)的相似度来确定成对集合的距离.对于相邻切片中的超像素对, 选取灰度距离和超像素对质心的物理距离作为特征, 加权计算非同层超像素对的相似度:

(2)

式中:u, v分别属于相邻切片G_i, G_i+1上的超像素; d_B(u, v)表示u, v灰度特征的Bhattacharyya距离; d_E(u, v)表示u, v的欧氏距离; σ为权重系数, 决定了欧式距离对相似度的影响.

3 连续最大流分割目标区域 3.1 连续最大流

1956年, Ford等证明了最大流最小割定理:网络流中从源点S到汇点T的最大流量等于网络的最小容量和^[8].根据此定理, 可以使用最大流算法来寻找最小割能量函数的全局最小值.

离散最大流方法存在网格偏差, 且需要存储的信息量较大.相对于离散方法, 连续最大流度量误差更小且可并行实现.鉴于此, 本文采用连续最大流方法实现最小割的求解.

Yuan等^[9]提出一种连续的最大流方法求解非凸分割问题.该方法具有收敛速度快、参数选取广泛的优点.具体过程如下:

在连续2D或3D区域Ω中, 对于∀x∈Ω, 设经过x的一般流为p(x), 从源点S流向x的有向流为p_s(x), 从x流向汇点T的有向流为p_t(x).C(x), C_s(x), C_t(x)分别表示相应流的流量限制, 最大流模型为

且满足

(3)

在式(3)中引入对偶变量λ, 则该最大流模型可转化为如下鞍点问题:

(4)

其中div p表示通过局部区域x的流量总和.

3.2 连续最小割

最优化对偶变量λ, 鞍点问题(4)将等价于最大流模型(3).同样, 最优化式(4)中的流变量p_s, p_t和p可以得到与其等价的连续最小割模型:

(5)

求解连续最大流式(3)等价于求解对偶最小割式(5).对此文献[10]已给出证明, 不再赘述.因此只需求解式(4)的全局最小、最大值, 即可获得最小割的解.具体见算法1.

算法1 连续最大流算法

3.3 算法优化

对构造的三维图结构使用二维矩阵存储, 在此基础上进行连续最大流分割类似于普通二维图像的分割处理, 算法简单灵活.但由此产生的矩阵维度过大的问题亟待解决:每张切片大小为512×512, 扫描层数在100~200不等, 形成的巨型矩阵复杂度过高, 直接处理较困难.提出一种解决方案:使用正则化稀疏模型对矩阵进行稀疏表示.

采用Simon等^[10]提出的稀疏组Lasso方法对2.1节产生的大型矩阵进行稀疏化处理, 将L_{2, 1}范数罚作为正则项:

(6)

式中: y∈ R^N为响应向量; β∈ R^P表示回归系数向量^[11].稀疏约束后的矩阵再用于空间图结构构建能够有效降低数据维度, 减轻计算压力.正则化稀疏算法如算法2所示.

算法2 正则化稀疏算法

4 实验结果和分析

本文使用的肺CT影像由辽宁中医药大学附属医院提供, 业已通过伦理审查委员会(IRB)的审批, 审批通知编号为2017111CS(KT)-040-01.从辽宁中医药大学附属医院所提供的数据库中选取146组CT扫描数据, 共包含14928张肺CT切片, 进行实验.图 5展示了肺实质中肺结节的分割结果.

在空间连通性的约束下, 分割结果中的血管、支气管、肺结节等组织以连通状态存在, 更易于ROI区域的提取.然而受切片数量、厚度等客观因素影响, 一组CT扫描的分割结果中血管(肺动脉、肺静脉)的分割数量在2~20之间不等, 支气管也可能被分割成多段.对实验结果进行统计, 分割所得的肺实质组织中, 支气管数量为169个, 血管数量为677个, 肺结节数量为241个.

为了量化分割方法的性能, 参考文献[12]对三维医疗图像分割评价方法的综合评价, 采用DICE系数、敏感性(sensitivity)和特异性(specificity)作为衡量标准.TP(true positive), FP(false positive), TN(true negative), FN(false negative)分别是构成评价指标的4个基础基数.使用S_g表示金标准区域, S_t表示实验分割区域.评估指标定义如下:

(7)

(8)

(9)

作为肺结节分类工作的前期准备, 本文旨在提取所有候选结节.对照金标准提取出分割结果中的肺结节结构, 并对其采用指标式(7)~式(9)进行评价, 具体数据见表 1.分析表中数据可知, 本文算法针对不同类型肺结节的分割均取得较为满意的效果, 健壮性高且适用范围广.

5 结论

1) 本文提出一种结合空间约束和图割的三维肺实质自动分割算法, 采用超像素构造方法对肺CT图像序列进行超像素构建, 降低了数据计算的复杂度.

2) 基于所得超像素图像序列构造肺实质的三维空间结构, 充分利用了图像序列的灰度连续性和空间信息, 提高了分割的准确性.

3) 应用连续最大流方法分割肺实质, 针对内存占用和矩阵维度过高的问题, 设计了合适的优化算法, 采用合并子空间节点和矩阵稀疏化的方法进行优化, 在保证分割效果的前提下, 显著降低图像维度, 提升了算法效率.

4) 实验表明本文算法分割的平均准确度达到80%以上, 对不同类型肺结节均可取得较好的分割效果.与当前广受关注的深度学习方法相比, 在训练数据集有限的情况下, 本文算法仍可获得满意的分割结果, 且分割效果不受数据集影响, 在不同数据下均可快速获得准确的肺实质分割结果, 具有良好的稳定性和广泛的适用性.