基于一致性法的肺部CT图像血管提取

摘要：在肺部CT图像中，血管与背景的对比度较低，很难分辨肿瘤和血管结节。为了解决这个问题，本文提出具有旋转不变性的一致性法。根据肺部CT影像细小血管具有局部亮度和结构光滑的纹理特征检测细小血管。首先采用OTSU等算法快速定位CT图像的肺部区域并对定位出的结果进行去噪，应用图像二值化方法分割出粗大血管，然后对没有粗大血管肺部区域的各个子区域采用一致性法进行分类计算，最后根据细小血管的纹理特征值使用支持向量机对有细小血管的子区域及有肺结节的子区域进行样本训练，判断是否是细小血管然后将其提取。实验说明该方法是有效的。
论文关键词：肺部CT图像,一致性法,OTSU,纹理特征,支持向量机

　　随着CT扫描技术的发展，针对肺部疾病的诊断，越来越广泛使用肺部CT计算机辅助诊断系统(CAD)[1,2,3,4], 在肺部CT疾病众多，如肺栓塞，肺结节等，应用计算机辅助诊断系统对这些疾病进行早期诊断达到早期治疗的目的是最好的解决办法[5,6,7]。对于早期肿瘤以及与肺血管粘连的肿瘤, 往往很难分辨肿瘤和血管结节, 因此计算机辅助诊断系统对这种肿瘤的识别率较低。

　　现有的血管提取方法包括阈值法和增强血管象素平衡法及跟踪算法。阈值方法简单易行, 但是鉴于肺部CT 图像中血管与背景的低对比度, 因此提取效果不理想。象素平衡法将力学平衡的观点引入象素灰度中, 给出了象素平衡计算公式, 当图像中所有象素均平衡时, 整幅图像平衡, 此时图像中目标与背景对比度最大。象素平衡法大量采用公式迭代运算, 计算量较大, 因此不适用于数据量较大的CT 图像。跟踪算法是在阈值细化方法处理后, 以检测出的血管为基础, 利用血管与背景梯度信息跟踪检测出在阈值细化方法中漏检的血管。采用OTSU 算法二值化肺部CT 图像, 形态学细化算法提取血管骨架, 依据血管光滑性确定跟踪方向并根据灰度梯度信息跟踪出更多的细小血管，跟踪算法是从单一方向进行跟踪, 而血管有很多分支, 并不是单一方向的, 因此还可以考虑从多个方向进行跟踪，所以该算法还有待进一步改善和提高。

　　本文的研究目的是快速提取肺部CT图像的大血管，然后针对细小血管与肺结节及肿瘤的纹理特征进行诊断与提取。根据研究动机，本文提出了一个新的特征提取系统即一致性法用于检测肺部大血管及细小血管。运用于肺部子区域的纹理特征信息精确诊断提取肺部的血管，大大提高肺部细小血管的提取的精确性，提高了肺部疾病的早期诊断的准确性。

　　1相关工作

　　在肺部CT图像中，血管与背景的对比度较低，很难分辨肿瘤和血管结节。为了解决这个问题，根据肺部CT影像细小血管具有局部亮度和结构光滑的纹理特征检测细小血管。本文的算法是能够提取具有旋转不变性的纹理特征来实现血管提取。局部二进制模式是一种有效的纹理描述算子，由于其对图像局部纹理特征的卓越描绘能力而获得了广泛的运用。LBP[10,11]是子区域中3 × 3中心像素阈值相对于相邻像素值生成二进制代码的算法。如果相邻的像素值比中心像素的小，它产生一个二进制代码0，否则，它会生成一个二进制代码1。这些二进制代码乘以相应的权便可得出LBP生成的代码，其计算方法如下：

　　

　　其中(xc, yc)是中心点位置，gc是中心点像素值，gp中心点相邻像素点的像素值，P是中心点相邻像素点个数，R为半径

　　

　　图1显示了一个LBP代码产生的过程。

　　

　　图1. LBP代码生成过程

　　Figure.1 The generation process of LBP code

　　为了突出旋转不变特性，将LBP扩展为一个圆形“统一”模型[10]，八个相邻像素点组成一个半径为R的圆形区域。该方法方便计算，但相邻像素灰度值不一定与该像素位置完全符合。这种旋转不变的LBP可以计算如下：

　　

　　

　　参数U是用来估计对应于空间过渡的一致性，即二进制数按位0、1之间变化的次数。因此，U的值越大，局部区域亮度发生转变越多。图2是一个圆环“统一”模型。

　　

　　图2. 一个统一模型中不同U值的情况

　　Figure.2 A uniformity model of different U-values

　　考虑到肺部CT影像的复杂结构，干扰噪声，本文提出重新定义方程（3）～（5）如下：

　　

　　

　　

　　从式子（7）可以看出，在本文中心点像素与周边像素之间亮度的关系有三种情况。不同的结果（-1，0和1）代表了不同的情况。本文提出这种新的模式可以分清中心点和更细节的相邻像素亮度的关系。与LBP相比，本文定义的NLBP在估计亮度均匀性具有更好的效果，同时具有旋转不变性。

　　2 特征系统的提取

　　要提取纹理特征，并识别肺部的正常血管区域和病变区域。必须先定位CT图像中的肺部区域，并对它进行相应的处理。

　　2.1 肺部区域的定位

　　根据肺部CT图像的影像学和解剖学特点，首先利用最大类间方差法（OTSU 法）对图像进行预分割，然后利用区域生长及小面积消除方法剔除干扰信息，同时生成掩模图像，最后运用数学形态学方法对模板进行细化，将原始图像与掩模图像进行数学运算即可得到肺部区域。具体方法如图 3 所示。

　　

　　图3 肺部区域提取示意图

　　Figure.3 The steps of lung region location phase

　　2.1.1 快速分割算法

　　最大类间方差法由日本学者Nobuyuki Otsu[8]首先提出，是一种自适应的阈值确定方法，又叫大津法，简称 OTSU 法。该方法应用类判别法寻找最佳阈值，以获得最好的分离特性。利用 OTSU 法依据类间距离极大准则来确定区域分割阈值就意味着错分概率最小。

　　本文利用大量图像进行了实验，在此以一幅正常的肺部CT图像为例。图4(a)为原始图像，(b) 为利用 OTSU 法预分割得到的结果。

　　

　　(a) 肺部CT图像 (b) OTSU法分割结果

　　（a）Lung CT image (b) enhanced image

　　图4 肺部CT图像OTSU法分割结果

　　Figure.4 Lung CT image after OTSU enhancement

　　经过 OTSU 法分割后，肺部CT图像根据灰度值分布情况将肺部区域和背景大致分开，但是图像中还存在检查床、心脏和血管等高密度区域都会对提取完整的肺部区域形成干扰，为了去除这些无关信息，我们利用基于区域生长的方法和小面积计算的方法继续进行分割。可以用图 5表示。

　　

　　图5. 区域生长流程图

　　Figure.5 Flow diagram of adaptive threshold

　　经过区域生长再分割和反色变换后得到的图像如图 6(a) 所示。除了大部分肺部组织外，气管、支气管因内部充满空气，也显示为低密度影区，而原本属于肺部组织的部分血管、结节、纤维化等则显示为高密度影区。我们分别对图像的两个密度区域进行连通域标记，测得各连通区域的面积，同时选取合适的面积阈值，并对面积小于相应阈值的区域内的像素值取反，从而弥补二值化过程带来的分类误差。处理结果如图 6(b) 所示。

　　

　　(a) 区域生长结果 (b) 小面积消除结果

　　图6. 去除无关信息的结果

　　Figure.6 Images after adaptive threshold

　　通过上述一系列的处理， 肺部模板已经基本成形， 但是，对图 6(b) 进行分析可以看到，由于肺实质边缘密度和周围组织非常相近，在肺部区域预分割时常常将其误分为背景，因此，本文利用形态学的闭运算对模板进行细化。图7(a) 为利用闭运算细化得到的最终模板，将原图与模板做减运算即得到了肺部区域的完整图像，如图 7(b) 所示。

　　

　　(a) 模板细化结果 (b) 肺部区域像

　　(a) Image after the morphology process （b）Image of Lung region

　　图7. 分割最终结果

　　Figure7. Image result

　　2.1.2 初步定位并剔除粗血管

　　先通过直方图均衡化提高上面得到的肺部区域图像的对比度，然后将图片二值化，就得到了粗血管，如图8（a）。通过图8（a）找出的粗血管图，映射到原图，将粗血管区域剔除掉。图8（b）为剔除粗血管后的肺部区域。

　　

　　（a）找出的粗血管图 （b）剔除粗血管后的图像

　　（a）the vessels （b）Lung region without vessels

　　图8.剔除肺部区域粗血管

　　Figure 9. Lung region image without big vessels

　　把最后对得到的图片进行分割，如图9所示，分成若干个25*25（单位：像素）的子区域。接下来的工作只对肺部区域占70%以上的小区域进行处理。在接下来的步骤中，没有粗血管的肺部区域被分成子区域。对于每个子区域，我们计算了基于灰度的NLBP值后计算梯度方向差,以得出条件概率密度函数，进行6个纹理特征量的计算。最后，根据6个纹理特征量使用支持向量机（SVM）对子区域进行分类 。

　　

　　图9. 对没有粗大血管的肺部区域进行分割

　　Figure 9. Sub-regions of the lung without big vessels

　　2.2 肺部纹理特征提取

　　NLBP能够描述中心像素与其邻边区域灰度值的关系，但是它不能说明这些血管像素光滑的局部结构信息。图10（a）和（b）是样品图像中两个3 × 3的局部区域模型。我们可以看到这两种模式的NLBP值是相同的,然而这两种模式的局部构造是不同的。

　　描述血管像素光滑的局部结构，提取纹理特征的有效性方面,统计家族GLCM方法是公认的有效方法,具有较强的适应能力和鲁棒性。该方法是建立在估计图像的二阶组合条件概率密度基础上的。GLCM是描述在θ方向上, 相隔d像元距离的一对像元分别具有灰度层i和j的出现概率。显然GLCM是一个对称矩阵,是距离和方向的函数, 其阶数由图像中的灰度级Ng 决定。尽管由GLCM提取的纹理特征具有较好的鉴别能力,但是这个方法在正确率依赖角度θ的选择。

　　本文对GLCM方法进行改进，提出新的不依赖角度具有旋转不变性的一致性法。

　　

　　（a）

　　

　　（b）

　　图10. 中心点的梯度方向

　　Figure center gradient orientation

　　设图像强度函数为S（x，y），计算点P（x，y）梯度方向如下：

　　（10）

　　如公式（10）所示，梯度方向表示点P（x，y）的局部结构。图10（a）和（b）表示两个中心点梯度方向。可以看出梯度方向取决于中心点周围的局部结构。在这种情况下，尽管这两中心点的NLBP值相同，但其梯度方向是不同的。因此，本文提出计算两点之间的梯度方向的差异，以评估其局部结构的相似性。图11中箭头表示模型各点的梯度方向的差异。

　　

　　图11. 梯度方向的差异

　　Figure 11. Gradient orientation differences

　　因此，点P（x，y） 和它的各邻点 P (xn, yn)之间的梯度差可通过以下计算得到：（11）

　　可知点P（x，y）的灰度值表示它的亮度，NLBP值表示点P（x，y）周围的亮度均匀性，以及梯度方向差异表示点P（x，y）局部结构的均匀性。假设邻点的个数为N，半径为R，邻点坐标分别为P（x1，y1），P（x2，y2），…，P（xn，yn），定义一个基于给定点P（x，y）与它的邻点P（x1，y1）之间的一致性如下：

　　（12）

　　其中gc是点P（x，y）的灰度值。同理， 以P（x1，y1），P（x2，y2），…，P（xn，yn）的值为基础通过公式（12）求得一致性。GNLBP是灰度值和NLBP的组合。因此，如果P（x，y）周围点灰度值均匀，那么GNLBP值较低。通过公式（12）定义条件概率密度函数f（GNLBP，GOD|N，R），其中0≤GNLBP