采用改进Snake模型的胸片肺野自动分割方法

胡俊，李平

(华侨大学 信息科学与工程学院，福建 厦门 361021)

X光胸片是临床上广泛使用的医学影像，对于肺部疾病的分析、诊断具有非常重要的临床意义.随着计算机辅助诊断(CAD)技术的普及，其高效性和实用性受到人们的关注，将计算机辅助诊断技术应用到X光胸片上，可有效地帮助医生检查关于肺部的各种疾病，如气胸、肺气肿、肺炎、肺癌等[1].因此，X光胸片肺野的精确分割成了计算机辅助医生进行高效、准确地诊断肺部疾病的重要基础[2].

一直以来，X光胸片肺野分割就是医学图像处理的研究热点，许多分割方法被应用于这一任务中.曹新华等[3]通过阈值法进行分割，并利用直方图特征自动确定阈值.张继武等[4]采用了水平集方法，分割效果要优于阈值法.刘炎等[5]采用了多种分割方法联合使用的策略，达到了提高分割精度的目的.近年来，深度学习日益受到人们的重视，文献[6-7]均采用神经网络对X光胸片肺野进行分割，虽然能得到较好的分割效果，但由于计算量巨大，对数据集、算力的要求颇高，难以推广到实际应用中去.就目前而言，传统方法也能实现较好的分割效果，对数据集、算力的要求又相对较低，所以将传统方法应用于医学图像分割时更具有优势.Kass等[8]提出的Snake模型作为一种经典的传统图像分割方法，被广泛应用于图像分割、运动跟踪、目标识别等领域[9].

传统的Snake模型对轮廓初始位置选择十分敏感，轮廓初始化是采用人机交互的方式完成的.然而，采用人工初始化轮廓具有极强的主观性，初始轮廓的选择完全依赖于人的判断，而初始轮廓又决定了分割结果.所以，经常会发生因初始轮廓选择的不合理而引起分割结果不够准确的情况，甚至会导致分割失败，出现欠分割、过分割的问题.针对上述问题，本文提出一种基于自动初始化Snake模型的X光胸片肺野自动分割方法.

1 X光胸片肺野分割中的问题

1.1 X光胸片图像的特点

X光成像的原理是基于人体不同组织密度、厚度的差别，当X线透过人体不同组织结构时，被吸收的程度有所不同，到达荧屏或胶片上的X线能量则存在差异，因此形成明暗或黑白对比不同的影像[10].相较于计算机断层扫描(CT)、磁共振成像(MRI)图像，X光胸片图像最大的缺陷是分辨率较低，容易引入噪声，对一些细微组织、病灶的对比度较差.由于X光成像具有投射、吸收的特点，导致成像时各组织投影呈现重叠的状态，往往使得各组织结构模糊不清.

1.2 X光胸片肺野分割存在的问题

图1是一张典型的X光正位胸片图像.因为背景、人体肺部内存在大量空气，密度较低，吸收的X线能量较少，在图像上呈现为黑影；而人体内其他组织、器官密度相对较大，吸收的X线能量多，在图像上则呈现为白影.从图1可知：虽然肺野与其他组织也表现出一定的灰度差，但由于组织投影彼此重叠，使得肺野与其他组织之间并没有非常明确的边界；此外，在肺野内部还包含大量骨骼、血管等组织信息.

图1 X光胸片图像

正是因为X光胸片图像存在上述特点，这使得肺野的精确分割成了一个难题.想要完成对肺野的精确分割，需要解决3个方面的问题：1)如何克服弱边界、伪边界的影响，避免分割陷入局部极小值；2)如何排除肺野内其他组织信息的干扰，防止收敛到错误的边界；3)如何增加对凹陷区域、角点的分割能力，提高分割的精度.只有解决了这几个方面的问题，才能获得较精确的分割结果.

2 X光胸片肺野的自动分割

所提出的X光胸片肺野自动分割方法，主要分为3个步骤:1)对图像进行预处理；2)采用Otsu法对图像进行二值化，再结合边缘检测和形态学方法完成Snake模型轮廓的自动初始化；3)利用自动初始化的轮廓，进行Snake模型分割，得到最终的分割结果.X光胸片肺野的自动分割流程，如图2所示.

图2 分割流程图

2.1 图像预处理

在X光胸片成像、存储及传输的过程中，通常会混入一些噪声.对于计算机辅助诊断系统来说，为了获得准确的分割输出，在分割之前必须消除噪声，而且在消除噪声的同时，还需要保留图像的基本特征.高斯滤波是比较常用的滤波技术，它在滤除噪声的同时不会丢失边缘信息等基本特征.

如图2所示，采用二维零均值、标准差σ=1的高斯滤波器对尺寸为1 024 px×1 024 px的原始X光胸片图像I1(x，y)进行去噪，得到输出图像I2(x，y).

二维零均值、标准差为σ的高斯函数表达式为

(1)

设二维零均值、标准差σ=1的高斯滤波器响应为H(r，s)，去噪过程可用离散卷积表示，即

(2)

式(2)中：k，l的值根据所选滤波器窗口的大小来确定，通常选择k，l=3或k，l=5的窗口，文中选择的滤波器窗口大小为k，l=3.

为了减小计算量，提高分割效率，需要对图像尺寸进行归一化.尺寸归一化是通过图像缩放来实现，尺寸归一化后得到尺寸缩小的图像.在尺寸归一化的过程中，会损失图像中部分局部信息，但保留了全局信息，这有助于提高算法对大目标的检测、分割能力，也有助于增强算法的鲁棒性.此外，为了提升算法对凹陷区域的分割能力，需要对肺野进行边界增强.边界增强能增加肺野边缘像素点的灰度值，以此来提高肺野的边缘梯度，从而增加肺野边缘对分割轮廓线的吸引力，这有助于轮廓线收敛到肺野的凹陷区域，可以提高分割精度.

如图2所示，对高斯去噪后的图像I2(x，y)进行缩放，得到尺寸为512 px×512 px的图像I3(x，y)；然后，利用图像金字塔[11]得到肺野边缘图像I4(x，y)；最后，通过将图像I3(x，y)与图像I4(x，y)进行融合，就可得到边界增强的图像I5(x，y).图像增强前后，分割结果对比，如图3所示.

从图3可知：对图像进行边界增强后，文中算法对凹陷区域的分割更为准确.其中，红色(虚线)为初始轮廓线，蓝色(实线)为分割结束时的轮廓线.

(a)原图 (b)对原图的分割结果 (c)边缘增强图 (d)对增强图的分割结果

2.2 基于改进Snake模型的X光胸片肺野分割

2.2.1 Snake模型算法简介 考虑到Snake模型在应用于图像分割时，具有只关注分割轮廓线邻域内图像信息的特点，而初始分割轮廓通常又位于待分割目标的真实边界附近.所以，采用Snake模型算法对X光胸片肺野进行分割，这样可以有效降低X光胸片肺野内部及背景区域中其他组织信息对分割过程的干扰.传统的Snake模型依赖人工来初始化轮廓，分割效果在很大程度上取决于人对初始轮廓的选择.这种初始化方式既不能保证分割结果的准确性和稳定性，也不利于批量化的处理.所以，文中在应用Snake模型进行肺野分割之前，对Snake模型轮廓进行了自动初始化.

作为一种经典的活动轮廓模型，Snake模型是Kass等[12]于1987年提出的，最开始用于追踪人嘴部的运动.由于其具有能将图像底层特征(灰度、纹理等)和上层信息有机结合起来的特点，受到国内外大量学者的关注与研究.时至今日，它已经被广泛的应用于目标分割、运动追踪、边缘检测等领域.

Snake模型通过设定一条可变形的闭合参数曲线，并利用曲线变形的具体规律来定义相应的能量函数[13]；然后，通过最小化能量函数，来驱动参数曲线发生形变.在构造能量函数时，把待分割目标轮廓处设置为能量最小处；而在最小化能量函数的过程中，参数曲线会逐步收敛到待分割目标轮廓处，以此完成对目标的分割.

2.2.2 Snake模型轮廓的自动初始化 文中Snake模型轮廓的自动初始化分为背景区域的消除和肺野区域外轮廓的提取两个部分.背景消除旨在于消除肺野以外的区域，只保留肺野区域.这个过程是为了确保在进行肺野轮廓提取的时候不受背景信息的干扰，以得到较准确的肺野区域外轮廓.在进行背景消除之前，需要对图像进行二值化处理.图像二值化就是通过一个合适的阈值，把图像分割成只包含前景、背景这两部分的黑白图像.这个过程会损失图像的纹理、颜色等特征，但也突出了形状特征.文中Snake模型轮廓的自动初始化，只需要利用图像中的形状特征.因此，对X光胸片图像进行二值化处理，能起到简化图像的作用，有利于Snake模型轮廓的自动初始化.

假设图像中当前像素点的灰度值为G(x，y)，通过遍历图像，对图像中每一像素点的灰度值进行操作，即

(3)

式(3)中：T为所选择的灰度阈值.

选择阈值的方法有很多，最常用的有双峰直方图法[14]、模糊熵法[15]、Otsu法[16].文中分别用这3种方法对边界增强后的图像I5(x，y)进行二值化，得到效果图，如图4所示.

从图4可知：双峰直方图法和模糊熵法都不能避免骨骼带来的影响，二值化后的肺野区域内部包含大量的骨骼边缘，若此时进行边缘提取，难以得到精确的肺野区域外轮廓；而Otsu法则不受骨骼纹理的干扰，相较于前两种方法，Otsu法的二值化效果最好.因此，选择Otsu法对图像进行二值化处理，得到二值图像S1(x，y).

(a)原图 (b)双峰直方图法 (c)模糊熵法 (d)Otsu法

一般来说，对于二值图像，通常把白色区域称为前景区域，黑色区域称为背景区域，目标区域一般归属于前景区域.但对于X光胸片来说，肺野区域像素点的灰度值通常比其他区域的要低，在经过二值化以后，肺野区域内像素点灰度值被置0，所以在图像中显示为黑色.为了方便对肺野区域进一步的处理，需要对二值化后的X光胸片图像进行取反，即把原本显示为白色的区域转换成黑色，原本显示为黑色的区域转换成白色.

假设二值图像中当前像素点的灰度值为G(x，y)，因为二值图像中所有像素的灰度值只有0和255这两种，故只需要经过如下操作就可以得到取反后的灰度值，其表达式为

G(x，y)=255-G(x，y).

(4)

通过遍历图像S1(x，y)中所有像素点，并对每个像素点执行式(4)的操作，就可得到取反后的二值图像S2(x，y)，图像S2(x，y)中肺野区域此时显示为白色.然而，在整张图像中，除了肺野区域以外，还有大片背景区域也显示为白色，这些背景区域会给肺野轮廓的提取带来干扰，所以需要消除这些背景区域.对于二值图像S2(x，y)，通过连通域检测，标记出图像中所有连通域，然后通过消除背景所在连通域，就得到只含有肺野区域的二值图像S3(x，y).

接下来，需要对图像进行边缘检测，其目的在于获取肺野区域的外轮廓，用以Snake模型轮廓的自动初始化.

一般来说，Snake模型的初始轮廓需要满足两个条件：一是初始轮廓需要包裹目标区域；二是初始轮廓需要尽量接近真实的目标边界.满足这两个条件的初始轮廓，更有利于肺部区域的分割.然而，由于图像S3(x，y)的肺野区域内存在孔洞，若直接进行边缘检测，孔洞的边缘也会被检测出来，这将干扰初始轮廓的确定.所以，在边缘检测之前，还需要消除孔洞.通过形态学膨胀，不仅可以消除这些孔洞，还能略微扩大肺野区域，此时再通过边缘检测所得的轮廓更满足Snake模型对初始轮廓的要求.由图2可知：对图像S3(x，y)进行形态学膨胀，得到肺野区域内不含孔洞的图像S4(x，y).

边缘检测的方法多种多样，最常用的方法是使用边缘检测算子来检测边缘.边缘检测算子是利用图像边缘的突变性质来检测边缘，主要分成两种类型[17]：一类是以一阶导数为基础，通过计算图像的梯度值来检测边缘，如Robert算子、Sobel算子、Prewitt算子等；另一类是以二阶导数为基础，通过寻找二阶导数的过零点来检测边缘，如Laplacian算子、LOG算子、Canny算子等.Canny算子作为最常用的边缘检测算子，具有结构简单，检测出的边缘清晰、连续的优势，因此选择使用Canny算子对二值图像进行边缘检测.

用Canny算子对形态学膨胀前后的图像进行边缘检测，结果如图5所示.从图5可知：对于未经形态学膨胀的图像，其肺野区域内部的孔洞边缘也被检测出来；而在经过形态学膨胀后，肺野区域内部的孔洞被消除，通过边缘检测，得到唯一的肺野外轮廓.

(a)未经形态学膨胀 (b)膨胀前的检测结果 (c)经过形态学膨胀 (d)膨胀后的检测结果

通过对图像S4(x，y)进行边缘检测，得到了肺野区域的外轮廓，然后在所得外轮廓上等距离地选择控制点，用以构造Snake模型的初始轮廓.经此步骤，即完成了Snake模型轮廓的自动初始化.初始轮廓可视化，如图6所示.图6中：红色(虚线)为Snake模型自动初始化轮廓的所在处.通过实验可知，在外轮廓上等距离地选择一半的点所构造出的初始轮廓，既能满足分割精度的需要，也能保证分割效率.

图6 初始轮廓可视化

2.2.1 X光胸片肺野分割过程 在完成Snake模型轮廓的自动初始化后，利用Snake模型算法对X光胸片进行肺野分割.X光胸片肺野的分割过程，实际上是通过最小化能量函数驱动Snake模型参数曲线移动的过程.

在实验中，文中将自动初始化的轮廓线作为Snake模型算法的参数曲线，将初始轮廓上各像素点作为控制点，表示为

v(s)=(x(s)，y(s))s∈(0，1).

(5)

式(5)中：x(s)，y(s)是每个控制点的坐标位置；s是以傅里叶变换形式描述曲线弧长的自变量.

在参数曲线上定义Snake模型的能量函数，表示为

(6)

由式(6)可知：Snake模型的能量函数由两部分组成，一部分称之为内部能量，另一部分称为外部能量.内部能量仅与参数曲线的形状有关，而外部能量则由图像特征构成，例如灰度值、梯度等.文中使用弹性能量和弯曲能量构造内部能量，用梯度来构造外部能量.

内部能量表示为

(7)

式(7)中：一阶导数项称为弹性能量，反映曲线的连续性，弹性能量越大，曲线越不易被拉伸；α是弹性系数，它的值越大，曲线收缩越快，若其值为零，允许曲线产生不连续的点；二阶导数项称为弯曲能量，反映曲线的平滑性，弯曲能量越大，曲线越不易变形；β是弯曲系数，它的值越大，曲线越平滑，若其值为零，允许曲线产生拐角.

外部能量表示为

(8)

式(8)中：∇为梯度算子.由于外部能量和梯度值有关，且此项是非正项，所以梯度值越大，对应的外部能量越小.当参数曲线靠近目标边缘时，梯度值会增大，外部能量则会减小；而当参数曲线收敛到待分割目标轮廓处时，梯度值达到最大，外部能量则达到最小.

在定义了能量函数以后，就需要通过最小化能量函数来驱动参数曲线移动，以达到分割的目的.因为图像中像素点都是离散的，所以最小化能量函数的是一个典型的变分问题.通过变分法求解时，参数曲线需满足欧拉方程，即有

-αv″+βv″″+∇Eext=0.

(9)

对于图像而言，可以用差分近似代替微分，即有

vs=(xs，ys)，

(10)

v″≈vs+1+vs-1-2vs，

(11)

v″″≈(vs+2+vs-2vs+1)+(vs+vs-2-2vs-1)-2(vs+1+vs-1-2vs).

(12)

把式(8)，(10)，(11)，(12)代入式(9)中，整理可得

(13)

把Snake模型的参数曲线看成关于时间t的函数，当参数曲线移动到目标边缘处后，不再随时间发生变化.此时v(t)=(xt，yt)=0，代入式(13)中，求解可得

(14)

式(14)中：γ为时间步长.

通过对Snake模型参数曲线上的每一控制点执行上述计算，得到每个控制点新位置的坐标，以此达到驱使参数曲线移动的目的.通过算法的迭代运算，驱使参数曲线向肺野轮廓处移动，当达到最大迭代次数或者所有控制点都满足停止条件时，迭代停止.

停止条件为

|Eext(v(s))+Eint(v(s))|<δ.

(15)

式(15)中：δ为停止阈值，通常为趋近于0的正数.停止阈值需要根据算法需求来确定.一般来说，选择较小的阈值有利于提升分割精度，选择较大的阈值有利于提升分割效率.经过大量实验，文中将此阈值选为0.1，即可满足本文对分割算法要同时拥有高精度、高效率的需求.

由式(15)可知：当迭代停止时，每个控制点的总能量近似为0，保证了剩余的能量无法再使任何一个控制点移动，这代表了分割过程的结束.

X光胸片肺野的分割结果，如图7所示.图7(a)中的蓝色曲线为分割结束时的Snake模型轮廓所在位置；图7(b)为采用文中所提方法分割X光胸片肺野的结果.

(a)最终轮廓可视化 (b)肺野分割结果的二值图

3 实验结果与分析

3.1 不同方法的分割效果对比

实验所使用的X光胸片数据集来自AI研习社所发布的肺炎X光病灶识别竞赛，图像尺寸大小均为1 024 px×1 024px.文中进行了两组实验.第一组是采用不同方法对同一张X光胸片进行肺野分割，然后对分割结果进行比较、分析，用以评判文中所提方法的分割效果；第二组是采用文中所提方法对不同X光胸片进行肺野分割，通过对分割结果进行比较、分析，用以评判方法的适用性.

分别使用文献[5]中所提的柔性形态学与聚类分割方法、文献[18]中所提的基于ASM和GrabCut分割方法、文献[19]中所提的密度特征匹配分割方法、人工初始化轮廓的传统Snake模型分割方法和文中所提的自动初始化轮廓Snake模型分割方法对同一张X光胸片原图进行肺野分割，得到的效果图，如图8所示.

从图8可知：文献[5]方法对于X光胸片肺野的分割不够精确，所得到的肺野轮廓也不够平滑.主要原因是X光胸片中的弱、伪边界及骨骼纹理对分割造成了很大干扰，使得分割轮廓线难以收敛到真实的边界.文献[18]方法虽然能排除弱、伪边界以及骨骼纹理的干扰，但对凹陷区域的分割却不够准确.文献[19]方法具有良好的凹陷区域分割能力，但容易受到血管以及骨骼纹理的影响，导致分割精度的下降.传统Snake模型方法能得到平滑的肺野轮廓，但不能准确分割出凹陷区域，也无法排除弱、伪边界造成的干扰，容易陷入局部极小值.相较于以上方法，文中所提方法既能有效地克服弱、伪边界及骨骼纹理的干扰，又能完成对凹陷区域的精准分割.通过对比可知，文中所提方法的分割效果更好.

为了定量的分析不同方法对X光胸片肺野的分割性能，引入了灵敏度和特异性作为评判标准.灵敏度是指在X光胸片肺野区域内，正确分类的肺野像素占整个肺野区域像素的比例；特异性是指在X光胸片背景区域内，正确分类的背景像素占整个背景区域像素的比例.灵敏度越高，代表对肺野的分割越准确；特异性越高，代表把背景误分割成肺野的比率越小.灵敏度和特异性的表达式分别为

IS=NTP/(NTP+NFN)，

(16)

IP=NTN/(NTN+NFP).

(17)

式(16)中：IS表示灵敏度；IP表示特异性；NTP表示正确分类的肺野像素数量，NFN表示将肺野像素误判成背景像素的数量；NTN表示正确分类的背景像素数量，NFP表示将背景像素误判成肺野像素的数量.

对上述5种方法所得的X光胸片肺野分割结果数据进行整理，结果如表1所示.由表1可知，相较于其他方法，文中所提方法的灵敏度和特异性均更高，表明其具备更好的分割性能，对X光胸片肺野的分割更为准确.

表1 5种分割方法的性能对比

3.2 文中所提方法的分割效果对比

为了测试文中所提方法的适用性，从上述X光胸片数据集中随机选择4张肺野形态大小各异，亮度、对比度各不相同的图像，按照顺序进行编号；然后，通过采用文中所提方法进行测试，得到分割结果，如图9所示.

(a)图像1原图 (b)图像1分割结果 (c)图像2原图 (d)图像2分割结果

由图9可知：文中所提方法对肺野形态、大小各不相同的X光胸片图像均能完成有效的分割，所分割出的肺野边缘连续、平滑.对于亮度较大、对比度较低的X光胸片图像，文中所提方法也能完成较精确的分割.当图像中肺野边缘存在凹陷区域时，分割曲线能有效收敛到凹陷处；此外，对于图像中弱边界、伪边界的干扰，文中所提方法能有效排除，不易陷入局部极小值，大大提高了对肺野的分割精度.

对上述所得的X光胸片肺野分割结果数据进行整理，结果表2所示.由表2可知，对于不同的X光胸片图像，文中所提方法都获得了较高的灵敏度和特异性，表明文中所提方法对于不同的X光胸片都能完成较精确的分割，具有较好的适用性.

表2 文中所提方法对不同X光胸片的分割结果

4 结论

针对传统Snake模型依赖人工初始化轮廓而导致X光胸片肺野分割不准确的问题，提出了一种自动初始化Snake模型轮廓的X光胸片肺野分割方法.首先，利用Otsu法二值化X光胸片图像，分割出肺野区域和背景区域，再通过背景消除和边缘提取，完成对Snake模型轮廓的自动初始化，最后利用Snake模型算法对初始轮廓进行迭代，得到分割结果.

该方法摆脱了Snake模型对人工初始化轮廓的依赖，通过自动初始化轮廓，有效解决了传统Snake模型算法存在的一些问题，提高了分割效率和分割精度.此外，与人工初始化轮廓相比，自动初始化轮廓还具备一致性和稳定性的优点，提高了分割的鲁棒性.实验结果表明：相较于传统的Snake模型算法以及一些其他的分割方法，该方法简单实用，不仅提高了分割结果的准确性，而且还具有良好的适用性和稳定性.