融合Harris角点检测算法的肺实质分割方法

孙 红，李 晶

(上海理工大学，上海 200093) (上海现代光学系统重点实验室，上海 200093)

1 引 言

由于大气污染和吸烟，肺癌是世界上对人类健康和生命威胁最大的恶性肿瘤之一，发病率和死亡率目前均占所有恶性肿瘤的首位.降低肺癌死亡率、提高病人生存质量，最关键的是早期发现、早期诊断.早期肺癌往往表现为肺部小结节，肺部结节成为肿瘤是个渐进的过程，因此我们需要及时发现并鉴别诊断出肺结节.目前主要通过电子计算机断层扫描(CT扫描)来检测肺结节.

肺结节是指小于等于75px的肺部圆形病变.然而，由于肺内器官组织的复杂性，胸膜表面的肺结节呈现不同的特征，有些甚至仅通过肉眼无法轻易识别[1].CT扫描使得早期肺癌的发现率大大提高，而薄层CT技术有助于检测更小的肺结节.但是，在使用薄层CT技术时，CT成像数量明显增加，影像医师的读片量增加，工作量也增大，这样很可能出现对肺结节的漏检和错误判断的情况.

本文着重针对胸膜表面粘连的肺结节，运用凸包与凸缺陷算法与harris角点检测的方法对肺实质的轮廓进行修正，从而对这些的肺结节做出分割.本文采用的图像样本来源于公开数据库The Cancer Imaging Archive (TCIA)[2]中的SPIE-AAPM Lung CT Challenge样本集.

目前对肺结节自动检测的研究，主要分为：肺部图像的预处理方法、结节的特征提取、运用监督或非监督学习的方法对可疑结节进行检测.在对肺实质分割之前，通常运用最大类间方差法进行自动阈值分割得到易于分析的二值化的图像.对于与胸膜粘连的结节，Armato[3]提出了一种滚球算法.在对肺实质分割之前，对肺实质的轮廓进行补正，用来避免这些位于胸膜表面的结节的丢失.滚球法的主要的问题是需要对滚球模型设定半径，不同的半径会对分割结果产生很大的影响，容易产生过分割或边界补正不完全.Pei等[4]采用基于凸包操作和射线追踪的边界校正算法进行边界校正，此方法同样可能会造成过分割的情况.Shen等[5]提出了一种无参数的肺结节分割算法.采用双向链编码方法，结合支持向量机(SVM)分类器，选择性地平滑肺边界，同时尽量减少相邻区域的过度分割.HUANG[1]提出了一种基于边界逼近的肺实质分割方法.针对自然弯曲的肺实质内轮廓与外轮廓，通过判断自然凹陷最大深度与凹陷宽度的关系，区分自然凹陷与病灶点.由于肺部结构的复杂性以及边界的不确定性.以上两种算法[1，5]的缺点在于对于一些比较特殊的将连通的肺实质断开的胸膜结节，该算法将无法对边界进行补正.本文通过harris角点检测的方法来弥补这一缺陷.

2 方 法

本文提出的方法主要涵盖4个步骤：

1)运用最大类间方差法将原CT图像转化为二值图像；

2)初步提取肺实质部分的轮廓；

3)运用凸包与凸缺陷以及角点检测方法对边界进行校正，得到完整的模板；

4)遍历整个图像根据得到的校正后的模板分割出肺实质内部的所有结节候选点.

2.1 预处理

本文运用最大类间方差法(又称Otsu算法或大津法)对图像进行预处理以得到二值化的图像以便进一步分析.如图1(b)所示，在对一幅图像进行二值化处理时，希望得到一个最佳阈值将原图像分成前景和背景两个图像.在Otsu算法中，通过使得类间方差最大来衡量前后背景的差异，即当类间方差取得最大值时，得到的阈值就是最佳阈值.因为方差是灰度分布均匀性的一种度量，那么前景与背景间的类间方差越大就说明灰度差别越大，也就越能区分前景和背景.

图1 Otsu算法处理得到的二值化图像Fig.1 Processed by the otsu algorithm

为了分割出图1(a)中肺实质内部的肺结节可疑点，需要得到一个二值模板，这个模板对应于肺实质的部分的灰度值为255，而其余部分为0.用这个模板与原图像相乘即可得到需要分割的可疑结节点.如图2所示，先提取图1(b)中肺实质的轮廓，再对轮廓内部做区域填充，即可得到二值模板.

2.2 凸包与凸缺陷

在2.1节中，初步得到的模板与原图像相乘之后，图(1)b中圆圈位置所标注的凹陷，即胸膜表面的结节无法被分割.因为与胸膜相粘连的结节灰度值与胸膜接近，在基于灰度的分割方法中，会被误判为肺实质外的区域.为了修正这个问题，本文引入了轮廓的凸包与凸缺陷.

那么凸包是什么呢？给定平面上有一个(有限)点集，这个点集的凸包就是包含点集中所有点的最小面积的凸多边形.本文通过这种方法来分割出图1(b)中圆圈标注的凹陷.凸包最常用的凸包算法是Graham扫描法[6]和Jarvis步进法[7].本文运用的是Graham扫描法.将该方法运用在图2中的轮廓图像即可得到修正以后的模板图像图3(a)，将模板图像与原图像做乘运算后得到所有结节候选点的结果，如图3(b)所示.处于胸膜表面的结节以及内部的所有组织都被分割出来，以便进行下一步分析.

图2 轮廓提取与区域填充Fig.2 Contour extraction and area filling

图3 修正以后的结果.Fig.3 Result of correction

2.3 harris角点检测方法

针对文献[1，5]中提出的问题，由于肺实质内部结构复杂，同一名患者的CT图像在不同层的肺实质轮廓都不相同.

图4 方法的失效Fig.4 Method failed

对于一些特殊位置的结节，如图4所示.粘连于胸膜表面的结节使得左肺边缘出现“断裂”.

本文运用harris角点检测算法对该问题进行修正，得到期望的理想轮廓线与模板，从而使位于边缘的病灶不产生丢失.对于图4所产生的问题，根据大量的样本显示，这类情况多出现在肺实质区域狭长的层，因病灶的缘故使得肺实质之间产生“断裂”.

通过对图像进行行扫描可寻找到该断裂位置，并以该位置为中心设定感兴趣区域.在该断裂区域寻找到肺实质区域边缘的角点并将这些点用一个最小多边形将其连接起来，这样填充以后即可得到可以分割出病灶的模板.

Harris角点检测是由Harris在文献[8]中提出的，其主要思想是用一个模板窗口在图像上滑动，在垂直和水平方向上都产生较大变化的就是角点.首先定义自相关函数：

(1)

令D(x，y)=Iu+x，v+y-Iu，v，将D(x，y)在(0，0)按泰勒级数展开得：

(2)

有D(0，0)=Iu+0，v+0-Iu，v，高阶无穷小忽略不计：

(3)

所以(1)可以化为：

(4)

(5)

(5)式本质是一个椭圆函数.椭圆的扁率和尺寸是由矩阵H的特征值决定的，椭圆函数的特征值与图像中的角点，直线和平滑区域的关系可分为：

1)直线.一个特征值较大，另外一个特征值较小.

2)平滑区域.两个特征值都很小，且近似相等.自相关函数值在个方向上都很小.

3)角点.两个特征值都很大，且近似相等，自相关函数在所有方向上都很大.

为了判定图像中角点，根据hessian矩阵中行列式的值、矩阵的迹，将它们代入角点计算公式：

(6)

当R大于一个特定阈值时，该特征点(u，v)被认为是角点.由于通常情况下一个连通域分裂为两个区域，为了补正它们之间的缝隙，寻找曲率低的平滑的点，将这些点连接所围成的区域不足以使缺失区域补全，因此考虑寻找断裂区域周围的角点作为候选点进行补正.

3 实验结果与分析

3.1 harris角点检测方法

如图5所示，左肺的病灶出现在断裂的位置，运用凸缺陷的方法无法得到包含病灶点的模板.运用上文中的harris角点检测算法，首先通过行扫描的方法对每一行像素点的个数进行积分，若某一行积分值图变为0，则认为该位置存在断裂的区域，并将上一行中首个不为0的像素所在的列数作为ROI区域的中心点，并以该点为中心设置一个感兴趣40×40像素的区域.即仅对该断裂区域检测肺实质部分的角点，将满足条件的候选点按逆时针排列后将它们两两连结成一个多边形，如图5所示.

图5(b)是原图中感兴趣区域补正的结果，大图是结果全体.对该补正后的空洞区域做填充即得到整个连通的模板，用此模板与原图像相乘可得病灶候选点.本文对The Cancer Imaging Archive (TCIA)的SPIE-AAPM Lung CT Challenge样本集中的263张CT图像(来自不同的6例病例)，对于样本集中存在如文中所说凹陷的图像均可做出正确的补正.图6是实验结果，一例在左肺产生断裂，另一例是在右肺边缘产生病灶的样本图像，本结果来自于windows 8环境下、CPU主频2.5GHz、内存3.9G、Visual Studio 2013、opencv 2.4.9.

图5 原图像与补正后的结果Fig.5 Original image and result of the harris algorithm

图6 原图像与修正结果对比Fig.6 Original image and result of correction

3.2 滚球算法结果对比与分析

文献[3]中首次提出运用一个固定半径的圆绕着轮廓滚动一周.在这过程中，平滑区域应当与圆相切即只有一个交点.在凹陷处，即表面存在结节处应当存在一个以上的交点.将这两个交点用直线相连即可完成补正的方法.本结果来自于Image J图像处理软件自带的rolling ball算法.

图7为本文修正结果与文献[3]的滚球法的效果对比，图示说明中r表示滚球法中设定的圆半径参数.对照图2(b)修正前的图像，图(a)-(d)均正确地将结节处在的边缘进行了修复.图(b)的结果与本文方法相近.而图(c)和图(d)当增大参数半径时就会产生过度分割，即将非实质区域误认为是结节区域进行了不必要的修正.又如图(f)-(h)，在不同r参数下有不同的分割结果.由此可见，文献[3]中的分割方法非常依赖参数r的设定，该参数r并不具有自适应性.

图7 本文方法与滚球法对比Fig.7 Comparison between rolling ball algorithm and the method in this paper

3.3 DRLSE(Distance Regularized Level Set Evolution)模型的水平集方法结果与分析

DRLSE模型是由LI在文献[9]中提出的，运用该方法所产生的结果来自于 windows 8 环境下、CPU主频2.5GHz、内存3.9G、Matlab 2014b.

如图8所示，该方法相比于传统的水平集方法，提出了距离正则化项用来修正水平集函数与符号距离函数之间的差，从而保证最终演化结果精确地在边缘处停止并舍弃了重新初始化水平集函数的步骤.本文根据文献[9]通过水平集方法寻找到包含病灶的肺实质轮廓，从而得到分割结节候选点的二值模板.而由结果可见，该方法主要有二点问题，第一，水平集方法需要设定初始的轮廓位置，而CT影像在不同层的位置和大小都有所不同，所以曲线演化的初始位置不易选择.第二，由图8(c)可见，在曲线演化至接近目标边缘时，病灶位置的曲线的拓扑结构发生了变化，并未将该病灶包含在轮廓内部.

图8 水平集方法Fig.8 Level set method

3.4 边界逼近法结果与分析

将边界逼近法运用于肺实质分割是由HUANG在文献[1]中提出的，该方法中选取的样本大小均为512像素×512像素，实验环境为Windows 8.1、主频2.5 GHz、内存4GB，开发平台是VS2013、OPENCV2.4.9.

如图9(b)所示，该算法对于内轮廓上自然凹陷区域分割较好，同时也较好保留外轮廓存在病灶的区域.但是图9(d)所示结果，未检测出右肺的病灶点，粘连于胸膜表面的结节使得右肺边缘出现“断裂”.该算法主要的问题是：对如图9所示样本，因病灶的缘故使得肺实质之间产生“断裂”，将肺实质割裂成两块独立的连通域，从而并未将该病灶包含在轮廓内部.

图9 边界逼近法Fig. 9 Border approximation

4 讨论和结论

本文主要针对CT图像中位于胸膜表面的结节的检测提出了解决方法，运用graham扫描算法对肺部图像的轮廓进行修补，从而得到完整的二值模板来分割位于边缘的病灶.并通过rolling ball、水平集方法、边界逼近法与本文角点检测方法做了对比.由于rolling ball与水平集方法都需要手动设定固定的参数或初始化位置，而同一位病人的CT图像在不同层上的肺实质图像差异很大，因此难以设定初始值和初始轮廓，从而不易将边缘的病灶检测出来.同时边界逼近法也不能很好的检测出粘连于胸膜表面的结节，可能会出现“断裂”现象.Harris角点检测的方法对文献[1]和文献[5]提出的问题具有针对性，如实验结果可见，该区域无需手动设定参数，当寻找到该ROI区域后自动进行补正.因此，本文方法针对将连通的肺实质断开的胸膜结节检测具有有效性.本研究的下一步工作是运用机器学习的方法根据形态、大小、密度等特征对分割出的胸膜表面的结节进行综合分析，从而对其进行分类，以完善计算机辅助诊断系统.