基于改进随机游走算法的CT影像的肺裂检测

解德芳，李 彬+，刘霜纯，田联房，王立非

(1.华南理工大学 自动化科学与工程学院，广东 广州 510640；2.广东省人民医院 广东省医学科学院 广东省老年医学研究所康复学科，广东 广州 510080；3.深圳市第三人民医院放射科，广东 深圳 518112)

0 引 言

由于肺裂的边界非常模糊，灰度值的变化范围也较大，此外血管和气道等周围的肺结构、肿瘤等肺病症以及CT影像噪声干扰等，使得肺裂检测通常不能获得令人满意的效果[1]。Klinder等[2]运用线性增强滤波器原理进行检测肺裂，提出多方向假设以确定体素的线性方向，但该算法仅对低密度的肺裂具有较高的准确性；Ukil等[3]运用肺的解剖学知识定位肺裂位置，通过图搜索方法找到肺裂的最优三维曲面，但该算法涉及二次导数的计算，故对重建的假缺陷和噪声有较高的敏感性。

随机游走算法是使用最为广泛的基于图论的交互式图像分割方法之一，算法不仅对弱边界或遗漏的边界具有良好的响应，且计算速度快，分割结果好，对噪声具有较强的鲁棒性[4]。因此，对CT影像的肺裂分割，随机游走具有一定的优势，但传统的随机游走算法以图像像素之间的灰度值差来描述相邻节点的相似度，目标轮廓的分割易受到复杂纹理图像的干扰。针对上述问题，提出一种基于改进随机游走算法CT影像的肺裂检测算法，首先采用基于双边滤波的多尺度Retinex算法增强肺裂，通过增强处理使算法可对各种密度的肺裂进行检测，接着采用改进的随机游走算法对肺裂增强后的图像进行肺裂的分割，提高了随机游走算法对肺裂分割的准确率。

1 算法框架

本文提出一种肺裂检测的新算法，算法整体框架如图1所示，共分为3个主要部分。①基于自适应阈值法分割肺实质。②基于双边滤波的多尺度Retinex算法的肺裂增强。利用傅里叶变换将数据从空域转换到频域，把卷积转换为乘积的运算降低了运算量，采用双边滤波估计图像的亮度情况，并对亮度图像增加权值β，使图像更加自然协调，同时有效补偿了图像在估计亮度时产生的损失。③基于改进的随机游走算法进行肺裂的分割。对肺裂增强图像先进行Candy边缘检测，将肺裂所在区域进行展示，有利于快速准确选取种子点，提高肺裂分割的准确性。在传统随机游走算法基础上，提出了一种鲁棒性更强的权值函数融合灰度、位置和梯度方向信息的改进算法，完善了分割效果，扩展其在具体分割问题中的应用；以像素的8邻域内的梯度信息为边的权重赋值，使随机游走者能在更大范围内寻找最大权重所对应的边并引导随机游走者正确的选择种子点的走向，提高图像分割的精确性。

图1 算法整体框架

2 基于自适应阈值法的肺实质分割

最大类间OTSU阈值分割法是一种自适应确定阈值的肺实质分割方法，但OTSU只是对肺实质的粗分割，不能有效分割肺实质和背景[5]。针对上述问题，本文提出了一种基于自适应阈值法与形态学方法相结合的肺实质分割方案。首先，采用自适应阈值算法对图像进行二值化处理，以二值化的图像所体现的特征对分割效果进行检测评估，之后判断图像是否需要重新进行阈值分割处理，最后，运用形态学方法对分割后的初始肺实质图像进行修补和提取。算法框架如图2所示。

图2 肺实质提取流程

3 基于双边滤波的多尺度Retinex算法的肺裂图像增强

Retinex理论具有动态范围压缩和颜色恒常等特性，常被用于图像增强领域[6]。但Retinex理论对图像局部细节部分的处理仍然不理想，且在高对比度边缘区域易产生“光晕”现象。为解决上述问题，本文提出基于双边滤波的多尺度Retinex算法MSR(multi-scale Retinex)增强肺裂，采用双边滤波代替传统单尺度Retinex算法SSR(single-scale Retinex)中的高斯滤波并对亮度图像增加权值β，使增强后的肺裂图像整体更加自然协调。

3.1 肺裂增强算法框架

基于双边滤波的多尺度Retinex算法主要包含图像预处理、改进SSR算法和采用MSR算法3部分，流程如图3所示。

图3 肺裂增强流程

3.2 肺裂增强算法的实现

3.2.1 原始图像的预处理

图像预处理是指对图像进行取对数和傅里叶变换。由于本文算法涉及模板矩阵的卷积运算，若直接在空域上进行计算则运算量非常大，所以采用傅里叶变换先将数据从空域转换到频域，把卷积运算转换为乘积运算降低运算量；取对数不仅可以减小图像的动态范围，还能提高程序运行速度。

3.2.2 基于改进的Retinex算法的肺裂图像增强

Retinex理论的目标是将一幅图I(x,y)分解为决定图像动态范围的照射分量L(x,y)和决定图像内在性质的反射分量R(x,y)[7]，如式(1)所示

I(x,y)=L(x,y)*R(x,y)

(1)

本文针对反射分量R(x,y)进行局部对比度增强处理；针对照射分量L(x,y)采用改进的SSR算法处理。接着对3个SSR算法进行加权平均，即采用MSR算法。改进的Retinex算法不仅可以有效地去除“光晕”现象和解决图像过亮问题，还能在恢复阴影区域的过程中保持图像高亮部分的细节信息，使增强后的肺裂图像整体更加自然协调。

3.3 单尺度Retinex算法及改进

将式(1)两边取对数，整理得如式(2)、式(3)所示

lg[R(x,y)]=lg[I(x,y)]-lg[L(x,y)]

(2)

R(x,y)=lg[I(x,y)]-lg[F(x,y)*I(x,y)]=
lg{[I(x,y)/[F(x,y)*I(x,y)]}

(3)

式中：R(x,y)表示经Retinex算法输出的结果；*表示卷积算子；F(x,y)表示中心环绕函数，可由式(4)高斯函数来实现

(4)

其中，K由式(5)归一化函数来决定

(5)

基于传统SSR算法本文进行了如下改进：

(1)用双边滤波代替传统SSR算法中的高斯滤波实现对图像亮度的估计。

肺裂的边界非常模糊，灰度值的变化范围也较大，用各向同性的高斯滤波器估计光照会产生“光晕”现象，而双边滤波由于考虑了当前像素点和邻域像素点的灰度相似性和空间相近性故能在一定程度上改善“光晕”现象，双边滤波进行光照估计的公式如式(6)所示

(6)

其中，I(x,y)为输入图像经双边滤波后的照射分量，k(x)为归一化常数，可由式(7)得到，W表示窗口大小，f(x,y,m,n)、λ(x,y,m,n)和d(x,y,m,n)分别表示以当前点(x,y)为中心的邻域内点的像素值、亮度相似度和距离相似度，其中，d(x,y,m,n)可由式(8)得到

(7)

(8)

式(8)中，σd表示距离差尺度。用双边滤波对亮度估计的关键是在图像W×W窗口计算邻域内的点与当前点之间的亮度差，根据亮度差值进行判断比较，从而得到以当前点(x,y)为中心的邻域内点的亮度相似度[8]，其中λ(x,y,m,n) 可由式(9)得到

(9)

式中：f(x,y)为输出亮度值，M为[0,1]范围内的阈值，通常情况下，M值越小，图像亮度越大，σr表示亮度差尺度。

(2)对光照估计的图像增加权值β

传统SSR算法将光照估计的结果从输入图像中完全去除，如式(2)所示。针对该问题，本文在去除光照估计时加入权值β，如式(10)所示，β参数确定和增强效果将在5.1节阐述

lg[R(x,y)]=lg[I(x,y)]-βlg[L(x,y)]

(10)

4 基于改进随机游走算法分割肺裂

传统随机游走算法对弱边界或遗漏的边界具有良好的响应，计算速度快，分割结果好，但在构造权值函数时由于只考虑了相邻像素的相似度，使得分割的目标轮廓易受复杂纹理的干扰而出现误分割现象[11]。为解决上述问题，本文提出改进的随机游走算法，通过Candy边缘检测将肺裂所在区域进行展示，降低种子点的误选率，并提出将像素间的位置和梯度差信息融合进权值函数，选择像素的8领域描述像素之间的差别，增强分割准确性。

4.1 肺裂分割算法框架

基于改进的随机游走肺裂分割算法主要包含Candy边缘检测和改进随机游走算法两部分，算法流程如图4所示。

图4 肺裂分割流程

4.2 肺裂分割算法的实现

针对肺裂增强后的图像，本文算法先进行Candy边缘检测，将肺裂所在区域更直接形象的进行展示，为选取种子点提供参考，提高肺裂分割的准确性和快速性。在传统随机游走算法基础上，提出了一种鲁棒性更强的权值函数融合灰度、位置和梯度方向信息的改进算法，完善了分割效果。其次，在梯度函数融合梯度信息基础上，以像素的8邻域内的梯度信息为边的权重赋值，使随机游走者能在更大范围内寻找最大权重所对应的边并引导随机游走者正确的选择种子点的走向，提高了肺裂分割的精确度。

4.3 随机游走算法的改进

改进随机游走算法主要包含4个部分：将图像转换成有权无向图、选取种子点、建立权值函数和Dirichlet问题的转换。

4.3.1 将图像转换成有权无向图

图像转换后，像素点对应于图的顶点集，像素之间的连接关系对应于图的边集，顶点的属性对应于图像像素的特征信息，如灰度值、颜色和纹理等，边的大小由加权函数确定，权值代表像素之间的差异性或相似性[9]。

4.3.2 选取种子点

以Candy边缘检测为参照，在肺裂区域交互式选取一定数目的种子点，本文选取12个种子点。图5为肺裂分割过程示意图，其中图5(a)为肺裂增强后的图像，图5(b)为 Candy边缘检测结果，图5(c)为选取的种子点(用白色实心圆点表示)分布图，图5(d)肺裂分割结果(肺裂区域用黑线标注)。

图5 肺裂分割过程

4.3.3 建立权值函数

传统的随机游走算法在构造权值函数时只考虑了相邻像素的相似度，如式(11)所示

wij=exp(-ρ(gi-gj)2)

(11)

式中：顶点Vi和Vj处的灰度值分别为gi和gj，ρ表示权值信息的比例系数，是算法中唯一的可改变的参数[10]。

传统的随机游走算法抑制了随机游走者沿着与种子点灰度相近的边向种子点前进，从而易导致错分与漏分[11]。针对上述问题，本文提出一种鲁棒性更强的权值函数融合灰度、梯度方向和位置信息的改进算法。

(1)权值函数融合位置信息，如式(12)所示

(12)

式中：hi、hj分别为像素i、j的坐标位置；ε是为了保证wij>0而加入的，通常取ε=10-6。

(2)权值函数融合梯度信息

给定一个连续图像的函数f(i,j)，将坐标(i,j)处的梯度定义为一个矢量，如式(13)所示

(13)

其中，(i,j)方向上的梯度分别用Gi和Gj表示，利用模板函数在每个像素上进行卷积运算来近似计算上式的偏导数，把Gi和Gj各自所对应的模板结合起来便构成了一个梯度算子。

4.3.4 Dirichlet问题的转换

5 实验结果与分析

本文实验数据来自LIDC(lung image database consortium)数据库和广东省人民医院提供的临床医学肺部影像，选取4组CT序列作为实验对象，每组序列含50张512×512图像，其中前两组来自LIDC数据库，后两组来自广东省人民医院。实验环境为Windows 10 pro;MatlabR2015;CPU:Intel(R)Pentium(R);RAM:4GB。

5.1 光照估计时权值β取值确定

在基于双边滤波的多尺度Retinex肺裂增强算法中，为解决“光晕”现象和图像过亮问题，本文提出对光照估计的图像增加权值β，为确定最佳β值，对50幅CT图像选取不同β值进行肺裂增强处理，图6为其中4幅CT图像的肺裂增强效果。由图6可以看出，β值越小甚至为零时，图像越模糊，会出现严重的“光晕”现象和图像过亮问题，当β取值范围在[0.5,1]时，增强后的肺裂图像自然协调，“光晕”现象和图像过亮问题都得到了较好的解决。

图6 不同β值的肺裂增强效果对比

对应于图6(a)～图6(d)，选用均值、信息熵和清晰度进行定量分析，表1为对50幅图像统计结果的平均值。从表1可以看出，随着β值的增大，图像的平均信息熵增大，说明图像携带的信息量越多，即图像的细节部分如肺裂保持的最好；β取0.7和1.0时的清晰度和均值都能达到理想的效果，但β取0.7时的均值更接近最优均值(灰度值128附近)，同时，β取0.7的清晰度高于β取1.0时的清晰度。综上，本文β取值0.7使增强后的图像不仅视觉效果最好，同时在客观评价指标取得最优。

表1 不同β值的客观评价指标对比

5.2 肺裂增强效果分析

为验证本文算法增强肺裂性能，选用直方图均衡化和SSR算法作为对比，图7为选自各个序列的4幅CT影像的肺裂增强效果图，其中图7(a)为原始图像，图7(b)为采用自适应阈值法分割的肺实质结果，图7(c)为直方图均衡化增强肺裂结果，图7(d)为SSR算法增强肺裂结果，图7(e)为本文算法增强结果。

图7 肺裂增强结果

由图7可以看出，直方图均衡化算法只是对图像的灰度级进行了均匀化调整，效果是使图像的亮度整体提高，且在去除光照时产生了噪声。SSR算法虽然能够将受光照影响的肺部图像整体增强，但是对图像细节的处理仍然不理想，局部区域会产生“光晕”现象。本文算法不仅有效去除了“光晕”现象和图像过亮问题，还能在恢复阴影区域的过程中保持图像高亮部分的细节信息，使增强后的肺裂图像整体更加自然协调。实验结果表明，本文算法不仅能增强肺裂，而且具有良好的视觉效果。

为验证本文算法的客观评价指标取得最优，分别从均值、信息熵、清晰度以及算法平均执行时间来定量评估算法的优劣，表2为对第1组序列共50幅图像客观评价指标结果的平均值。

表2 算法的客观评价指标对比

从表2可以看出，本文算法的信息熵和清晰度都明显高于其它算法，均值更接近最优值(灰度值128附近)，且增强处理使得本文算法对中高密度的肺裂同样具有较高的准确性。由于本文算法只涉及模板矩阵的卷积运算，不涉及二次导数和肺血管分割，使得本文算法不仅降低了运算量同时减弱了对噪声的敏感度。综上，本文算法除了执行效率略低于SSR算法，均值、信息熵和清晰度都能达到满意的效果，增强的肺裂最清晰、效果最好，细节部分最突出。

5.3 肺裂检测结果分析

5.3.1 效果对比

为验证本文肺裂检测算法的有效性，对每组序列中的50幅增强处理后的图像进行本文肺裂分割算法处理，采用随机游走算法作为对比，以专家分割结果为依据，图8为第1组序列中的21幅CT图像的肺裂检测过程及效果对比图。其中，图8(a)为原始图像，图8(b)为肺实质分割结果，图8(c)为本文算法的肺裂增强结果，图8(d)为本文算法的肺裂分割结果，图8(e)为本文算法的肺裂检测结果，图8(f)为随机游走算法的肺裂分割结果，图8(g)为专家分割肺裂结果。由图8可以看出，传统随机游走算法不能完整的分割肺裂，肺裂出现断裂状态，且分割结果易受肺血管等组织的影响；本文算法分割的肺裂轮廓清晰准确，能反映肺裂的真实形态，更接近专家分割结果。

5.3.2 客观评价结果对比

本文选用Dice系数作为定量评估肺裂分割准确度的指标，分割性能对比见表3。令A表示本文算法或随机游走算法分割的结果，B表示专家手动分割的结果，即金标准，则Dice系数定义为式(14)

(14)

图8 肺裂检测过程与结果

通常，Dice系数越接近1，表明其分割准确率越高[13]。由表3数据可知，本文算法的Dice系数明显高于传统随机游走算法，分割结果更加接近金标准，肺裂分割精度较传统随机游走算法有了较大的提升。由于本文肺裂分割算法首先进行Candy边缘检测，将肺裂所在区域更直接形象进行展示，有利于我们快速准确选取种子点，提高了肺裂分割的准确性，同时提出将像素间的位置和梯度差信息融合进权值函数，选择像素的8邻域描述像素之间的差别，提高了随机游走者选择种子点走向的准确率。

表3 传统随机游走算法和本文算法的Dice系数对比(每组序列含50张CT图像)

另外，对传统随机游走算法和本文算法进行算法执行效率对比，对比结果见表4。从表4可以看出，尽管本文算法执行时间约是传统随机游走算法的两倍，但传统的随机游走算法不能实现同时分割多个肺裂，否则会出现严重的过分割现象，即会将多余的肺实质连同分割，且传统随机游走算法需要的种子点数目远大于本文算法，肺裂区域和背景区域均需选取种子点，使之操作较为复杂，另外由于传统的随机游走算法没有肺裂增强处理，故对选取种子点的技术要求较高。综上，本文肺裂检测算法可实现对多组临床医学肺部CT影像的肺裂检测，结果精确度较高。

表4 传统随机游走算法和本文算法的平均执行时间对比(每组序列含50张CT图像)

6 结束语

对于肺裂检测，首先基于自适应阈值法分割肺实质，接着基于双边滤波的多尺度Retinex算法进行肺裂的增强，之后将像素间的位置和梯度差信息融合进权值函数，采用改进的随机游走算法对肺裂增强后的图像进行肺裂的分割。实验结果表明，本文算法不仅可以抑制噪声，较好地保留边缘信息，且能准确快速的检测肺裂的弱边缘，分割的肺裂轮廓线比较清晰。但随机游走算法需要人工指定种子点，选取种子点的个数和位置对检测结果有一定的影响，因此下一步的研究工作的重点是如何提高算法的自动化程度。

[1]LIU Yang.Computer-aided method of detecting lung fissure based on CT image[D].Shenyang:Shenyang University,2012:1-8(in Chinese).[刘洋.基于CT影像的肺裂计算机辅助检测方法研究[D].沈阳:沈阳大学,2012:1-8]

[2]Klinder T,Wendland H,Wiemker R.Lobar fissure detection using line enhancing filters[J].Proceedings of SPIE-The International Society for Optical Engineering,2013,8669(6):598-608.

[3]Van R E M,Goldin J G,Maya G A,et al.A method for the automatic quantification of the completeness of pulmonary fissures:Evaluation in a database of subjects with severe emphysema[J].European Radiology,2012,22(2):302-309.

[4]Shen J,Du Y,Wang W,et al.Lazy random walks for superpixel segmentation[J].IEEE Transactions on Image Proces-sing,2014,23(4):1451.

[5]XU Changxin.PENG Guohua.Fast algorithm for 2D Otsu thresholding algorithm[J].Journal of Computer Applications,2012,32(5):1258-1260(in Chinese).[徐长新,彭国华.二维Otsu阈值法的快速算法[J].计算机应用,2012,32(5):1258-1260.]

[6]QIN Xujia,WANG Huiling,DU Yicheng,et al.Structured light image enhancement algorithm based on Retinex in HSV color space[J].Journal of Computer-Aided Design & Compu-ter Graphics,2013,25(4):488-493(in Chinese).[秦绪佳,王慧玲,杜轶诚,等.HSV色彩空间的Retinex结构光图像增强算法[J].计算机辅助设计与图形学学报,2013,25(4):488-493.]

[7]He K,Sun J,Tang X.Guided image filtering[J].IEEE Transactions on Pattern Analysis & Machine Intelligence,2013,35(6):1397-1409.

[8]ZHOU Yuwei,CHEN Qiang,SUN Quansen,et al.Remote sensing image enhancement based on dark channel prior and bilateral filtering[J].Journal of Image and Graphics,2014,19(2):313-321(in Chinese).[周雨薇,陈强,孙权森,等.结合暗通道原理和双边滤波的遥感图像增强[J].中国图象图形学报,2014,19(2):313-321.]

[9]Benjamini I,Kozma G,Wormald N.The mixing time of the giant component of a random graph[J].Random Structures & Algorithms,2016,45(3):383-407.

[10]CHEN Shengguo,SUN Zhengxing,ZHOU Jie.A segmentation method for stratum image based on FCM and random walks[J].ACTA Electronica Sinica,2013,41(3):526-531(in Chinese).[陈圣国,孙正兴,周杰.基于FCM和随机游走的地层图像分割方法[J].电子学报,2013,41(3):526-531.]

[11]GUO Pengfei,LIU Wanjun,LIN Lin,et al.Brain image segmentation method based on FCM and random walk[J].Computer Science,2014,41(7):322-325(in Chinese).[郭鹏飞,刘万军,林琳,等.结合随机游走与FCM的脑图像分割方法[J].计算机科学,2014,41(7):322-325.]

[12]He K,Sun J,Tang X.Guided image filtering[J].IEEE Transactions on Pattern Analysis & Machine Intelligence,2013,35(6):1397-1409.

[13]ZHU Wenbo,LI Bin,TIAN Lianfang,et al.A new coronary artery skeleton extraction method based on layered multiple hypothesis tracking[J].ACTA Automatica Sinica,2014,40(8):1783-1792(in Chinese).[朱文博,李彬,田联房,等.新型基于分层多假设跟踪的冠脉骨架提取算法[J].自动化学报,2014,40(8):1783-1792.]