基于LoG和因子分解的集成ACM的亚实性肺结节分割

陈晓楠 蒋 辉 刘晓凯 王凯欣

(大连海事大学信息科学技术学院 辽宁 大连 116026)

0 引 言

肺癌是世界上发病率和死亡率最高的癌症之一[1]，据国家癌症中心统计的数据，2019年全国肺癌的死亡率达到26.99%，占据我国十大癌症死亡率的首位，而且发病率和致死率仍在持续增长。当前，计算机断层扫描(Computed Tomography, CT)技术是检测肺癌最常用的手段之一，肺结节是肺癌在CT上的前期表现形式[2]，根据密度可以将肺结节划分为实质性肺结节(Solid Pulmonary Nodule, SPN)和亚实性肺结节(Sub-Solid Nodule, SSN)[3]。同SPN相比，SSN中既有实性成分也有周围磨玻璃影(Ground Glass Opacity, GGO)，其结构复杂，恶化程度也表现得比实质性肺结节高。因此，准确分割SSN对于前期的肺癌病变诊断有重要的临床价值。

关于提高肺结节的分割精度的研究，国内外有很多的文献记载，但是大多数都是基于SPN的分割。在SSN的分割中，大多数单一方法都不能准确而高效地分割出完整的肺结节。在近30年的医学图像分割的研究中，活动轮廓模型(Active Contour Model, ACM)作为一种水平集的能量泛型，因能量泛函广泛结合的理论性被作为医学图像分割方面的一个重要方法。Manickavasagam等[4]提出了一种基于梯度的活动轮廓模型(GACM)，通过构建梯度能量作为活动轮廓模型的边界检测函数，然而采用梯度信息容易丢失边界。文献[5-7]提出了小波能量和ACM结合的方法对肺结节进行分割，小波可以增强图像亮度，由小波能量构建ACM区域项，加强了肺结节与周围毛玻璃影的区分。Heewon等[8]提出CV活动轮廓模型结合贝叶斯方法来分割胸膜旁肺结节，在相似率上高于其他算法，但是该方法复杂、处理时间较长。Ding等[9]用高斯拉普拉斯(Laplacian of Gaussian, LoG)结合RSF模型的方法解决了弱边缘、区域亮度不均的图像分割问题。

为了解决SSN中区域亮度不均或区域不平滑、模糊边界等问题，本文提出基于LoG和因子分解的集成活动轮廓模型。首先，要得到肺实质的分割图像，这是肺结节分割的关键步骤；其次，用本文提供的方法来分割肺结节，LGIF模型[10]是一种全局区域和局部区域融合的能量框架模型，在此基础上，引入优化的LoG作为LGIF模型的一个能量项，解决亮度不均、边缘模糊的问题；然后，利用因子分解模型来构造能量函数，驱动边界轮廓曲线的演变，并平滑图像；最后，通过实验数据来验证本文算法对SSN分割的高精度性。

1 集成活动轮廓模型

如图1所示，本文对SSN的分割包含4个步骤：1) 分割原始CT图像，获取完整的肺实质图像；2) 构建优化的LoG作为一种能量项，增强较暗的灰度区域；3) 构建因子分解模型作为一种能量项来驱动轮廓曲线演变；4) 将优化的LoG能量项和因子分解能量项融入到LGIF模型中，对SSN进行分割。

图1 集成ACM模型的SSN分割流程

1.1 肺实质分割

由于肺部组织结构的复杂性，同时肺部CT图像可能受仪器设备等噪声影响，因此，为了SSN分割的高效性，需要对肺部CT图像进行肺实质分割。要将肺实质区域完整提取出来，需要去除肺气管、躯干、肺床、弧影等无关组织信息。

首先，需要对肺部CT图像粗分割，利用拉普拉斯算子获取轮廓边缘信息，并设定边界像素的阈值，通过轮廓像素的方法得到初始的肺部轮廓区域，如图2(b)所示；其次，需要对左右肺进行分割，通过设定一个圆形结构模板，以该模板进行闭运算和取反操作，获取左右肺的最大连通区域，图2(c)表示为带肺气管的左右肺；接着，需要填充孔洞和去除肺气管，用二值图像连通域标记来填充孔洞，并设定连通域的面积阈值获取完整的左右肺实质掩模图像，如图2(d)所示；最后，把原CT图像和掩模图像相乘，获取到完整的肺实质,如图2(e)所示。肺实质分割结果如图2所示。

(a) 原始CT图像(b) 初始轮廓图

(c) 左右肺分割(d) 肺实质掩模

(e) 提取的肺实质图2 肺实质整体分割过程

1.2 构建高斯拉普拉斯能量的活动轮廓模型项

由于SSN边缘的灰度值不连续，且结节亮度不均匀的特点，本文提出一种优化的高斯拉普拉斯能量项，来加强肺结节实性成分与周围背景的区分。首先，用高斯滤波器平滑处理肺实质图像I(x,y)。然后，用拉普拉斯算子处理肺结节的边缘，得到LoG算子与图像的卷积数学式为：

式中：Gσ(x,y)是标准差为σ的高斯核函数。经过LoG处理后的图像如图3所示，若区域的灰度值是固定值的时候，其LoG变换的响应值为0。在I(x,y)的灰度发生变化的区域，在较暗的一侧，响应值为负值；而在较亮的一侧，响应值则为正值。因此，肺结节边缘可以通过过零点的跳跃状态来检测。

图3 LoG处理的图像

如果图像受到噪声影响或者强度不均匀，则容易产生假的过零点边缘信息，因此，需要增强肺结节边缘和对区域平滑处理，提出优化的LoG能量泛函如下：

(L-β×Δ(Gσ×I))2dxdy

(2)

(a) 肺结节图(b) 原始的LoG图

(c) 优化的LoG图

(d) LoG响应对比图图4 原始LoG和优化LoG的图像对比

通过水平集的变分法，用水平集函数φ，简化式(2)，将优化的LoG构建为ACM的一个能量项EOL(φ)，表示如下：

EOL(φ)=mLΔ(φ)=m∬ΩH(φ)L(x,y)dxdy

(3)

式中：LΔ(φ)是优化的LoG算子项；m是权重项。本文用优化的高斯拉普拉斯处理后的图像是L(x,y)，而传统的拉普拉斯处理后的是图像是ΔI(x,y)，相关的公式为：

H(φ)是Heaviside函数，其公式为：

式中：ε>0是一个很小的数。

1.3 基于因子分解的ACM

传统的区域模型都是以灰度特征来作为描述因子，对于SSN模糊的边界信息可能分割不准确，而本文提出的因子分解模型是基于像素来分割的，引入该模型可以获得更平滑、更精确的目标边界。

因子分解模型是基于图像的纹理特征[11]，将图像区域以M×N矩阵表示，该矩阵是由每个以像素为中心的局部窗口组成，因此该区域可以表示成具有代表性特征的局部窗口的线性组合，以Y表示：

Y=R1QQ-1β1

(6)

式中：R1是M×L矩阵，R1的列表示待分割SSN区域的特征；β1是L×N矩阵，β1的列表示每个像素区域的权重向量；通过区域代表性特征归类和像素的聚类，可以得到R1的线性变换矩阵Q，用Q来表示某区域的像素特征。

该因子分解模型是将图像I(x,y)分为目标和背景两个区域，它们都是以像素为中心的局部窗口组成。设Ω是二维图像域(x,y)，Ωo和Ωb分别表示目标区域和背景区域，则像素特征模型Hw(x,y)为：

Hw(x,y)=wo(x,y)Ro+wb(x,y)Rb

(7)

式中：Rb、Rb是M维向量，它们分别是目标区域、背景区域的代表性特征矩阵；wo和wb是相应的权重向量矩阵，在(0,1)范围内。

R1的代表性特征是其线性变换的解，其特征向量在R1的L维子空间上，Q是L×L矩阵，它的列表示R1子空间上的每个代表特征的笛卡尔坐标，而本文提出目标和背景区域，即L=2，假设同一区域的像素特征是同类的，对于给定的输入图像I(x,y)：

H(φ)在式(5)中已被定义。基于特征矩阵Q可定义特征矩阵R如下：

R=[Ro,Rb]=R1×Q

(9)

wo和wo权重向量矩阵的计算方式如下：

w=[wo,wb]=(RRT)-1RTY

(10)

当wo(或wb)中某一个权重接近1时，另一个wb(wo)就接近0，通过上述的公式，创建关于因子分解的ACM模型的一个能量项Ed(φ)：

Ed(φ)=-∬Ω[H(φ)wo(x,y,R)+

(1-H(φ))wb(x,y,R)]dxdy

(11)

从式(11)得出，当能量曲线被吸引到目标与背景之间的边缘时，得到最小的能量项值。

图5是磨玻璃影SSN的灰度和像素的分布图。图5(a)是SSN的原始CT图像，而图5(b)显示的是磨玻璃影SSN的放大图像。从图5(c)和图5(d)的对比发现，肺结节的灰度差异并不明显，仅仅使用灰度特征描述不能达到很好的分割效果，而图5(d)的肺结节像素值可以直观显示像素的差异对比，以像素作为描述特征可以更好地保留肺结节的模糊成分。

(a) SSN原始CT图 (b) 磨玻璃影SSN放大图

(c) 灰度分布图

(d) 像素分布图图5 SSN的灰度分布和像素分布

1.4 集成ACM的能量演化

1) 构建ACM能量函数。本文的集成活动轮廓模型是在LGIF的能量泛型下结合高斯拉普拉斯能量项和因子分解特征矩阵能量项构建的，该集成ACM模型可以定义为：

E(φ)=μEOL(φ)+λEd(φ)+γELGIF(φ)

(12)

式中：μ、λ、γ为各项的权重参数，在本文的取值都为1；ELGIF是CV模型和LBF模型的结合[12]。ELGIF是融合了全局和局部信息的模型[13]，综合了二者的优点，该LGIF模型对边界模糊、灰度不均匀的图像分割能够达到较好的效果，公式如下：

ELGIF(φ)=θECV(φ)+(1-θ)ELBF(φ)

(13)

式中：θ是常量，0≤θ≤1。ECV(φ)和ELBF(φ)的计算公式为：

(14)

(15)

2) ACM的求解。在CV模型中，c1、c2的求法：

在LBF模型中，f1、f2的求法：

根据水平集函数φ和Euler-Lagrange方程理论：

最后，通过E(φ)对φ的变分法原理可得到φ的梯度下降流方程，即优化的LoG和因子分解结合的集成ACM的方程：

其中：

2 实验结果与分析

本文采用的亚实性肺结节的数据集都是来源于大连医科大学第二附属医院，该数据集共有19个序列，每个序列平均有300幅肺部CT图像，其中，包含亚实性肺结节的序列都是由专业医生标注。实验环境采用MATLAB R2018b，64位Windows10操作系统，3.8 GHz处理器。

2.1 定性验证

为了验证本文算法对SSN分割的有效性，从19个序列中挑选了6个典型样本，其中共包括90幅含肺结节的CT影像，大小均为512×512像素，对这些影像数据进行肺结节分割，随机选择6个样本中的单例进行展示，结果如图6所示。可以看出，第一、第二行的SSN整体都是磨玻璃状；第三行的SSN内部是实性成分，而周围是形态各异的磨玻璃影；第六行的SSN形态是边界磨玻璃影与血管黏连的情况。

(a) 原始肺部CT图像(b) 医生手动分割结果(c) 基于边界的ACM方法(d) 基于LGIF模型的方法(e) LGIF+因子分解模型(f) 本文分割算法图6 各种方法对SSN的分割结果对比

从图6(c)可以发现，基于边界的ACM方法对SSN分割时，容易丢失边界信息。从图6(d)和图6(e)对比发现，加了因子分解的LGIF模型分割SSN获得的边缘更完整，但是也存在复杂SSN过多分割背景的情况。从图6(e)和图6(f)对比看，在图6(e)的方法上加了优化的LoG能量项，可以增强SSN的边缘和区域的亮度，本文算法获得了更为完整的SSN轮廓。从图6(e)和(f)的第六行对比可以发现，本文算法可以滤除黏连的血管，直接对肺结节分割。从图6第二行的磨玻璃肺结节分割结果看，因肺结节与周围背景对比度低且内部亮度不均匀，其余三种模型都会产生欠分割与过分割的情况，只有本文方法对SSN的分割结果与医生手动标注结果接近。由图6(b)和图6(f)对比来看，本文算法对SSN的分割更接近于专业医生的手动分割结果。

对基于边界的ACM方法、基于LGIF模型的方法和LGIF与因子分解融合的方法，不能精确分割SSN的主要原因如下：1) 基于边界的ACM方法和基于LGIF模型的方法在轮廓曲线演变时是基于灰度信息来分割，并没有考虑周围模糊的磨玻璃影信息，因此，导致边界信息泄露。2) 优化的LoG能量可以加强结节亮度与边缘信息，提高SSN与其周围背景的对比度，而上述3种方法不能很好地解决亮度不均匀、假边缘等SSN的分割问题。

2.2 定量验证

为了进一步验证本文算法的分割精度，分别采用相似度系数(Dice coefficient,DICE)、灵敏度(Sensitivity)和差异度(Relative Volume Difference,RVD)这3个指标进行评价。通常将医生的手工标注结果作为真值结果，即“金标准”[14]，用“金标准”分别与文中所用算法进行比较。下列指标公式为：

式中：Seg为文中各种分割算法的结果；Sgt为医生手动分割结果；ρDICE表示各种方法和医生手动方法的分割结果的相似度，当相似度接近1时，表示分割效果越好；ρSensitivity表示实际分割SSN区域的像素点正确判断为“金标准”区域像素的比例，当比例值近似等于1时，效果越好；ρRVD反映的是分割区域差异度，当实际分割轮廓与医生标注轮廓接近时，差异度越小。

本文选择对90幅SSN的分割情况进行比较，结果如表1所示。其中，基于边界的ACM方法的各项指标性能都低于其他方法，其RVD值达到了70.60%；基于LGIF模型方法和LGIF模型与因子分解的模型方法虽然DICE值都能高达70%以上，但是Sensitivity指标都偏低，说明实际区域的像素点在真实区域中占比少；虽然，本文算法的RVD指标比LGIF模型方法和LGIF模型与因子分解的模型方法要高，但是，DICE指标达到了90.23%，比其他算法提高了近20%，其像素占比值也高(83.99%)。

表1 本文算法与其他算法的SSN分割结果评价(%)

3 结 语

本文提出一种基于LoG和因子分解的集成ACM算法对SSN进行分割。为了解决肺结节亮度不均、边界模糊的问题，本文引入了LoG算子并对其优化，构建了LoG能量项用来检测模糊的边界；其后，通过因子分解以像素值作为描述特征，构建一个特征向量的矩阵，提出因子分解的活动轮廓能量项，使分割曲线更接近于完整的肺结节。实验结果表明，本文算法对SSN的分割与专业医生的手动分割结果更为接近，比其他三种模型方法分割效果要好。