CN102360495B - 基于平均密度投影和平移高斯模型的肺结节分割方法 - Google Patents

Abstract

本发明公开了一种基于平均密度投影和平移高斯模型的肺结节分割方法，主要解决血管粘连型肺结节难以分割、分割精度低、假阳率高和处理数据量大的问题。其实现步骤是：利用平均密度投影AIP技术，对原始CT数据进行稀疏采样；其次，对AIP图像进行肺实质提取、阈值分割、形态学处理和几何特征粗提取，确定出结节大致轮廓；最后，建立平移高斯模型实现对肺结节的精确分割。本发明不但能够提取出血管粘连型肺结节，而且具有较强的鲁棒性，可有效抑制假阳，与传统方法相比分割精度更高，同时有效降低数据处理量，能够实现批处理。

Description

基于平均密度投影和平移高斯模型的肺结节分割方法

技术领域

本发明属于图像处理技术领域，特别涉及医学肺结节图像处理，具体地说是一种基于平均密度投影和平移高斯模型的肺结节分割方法，可用于血管粘连型肺结节以及一般类型的肺结节自动精确分割。

背景技术

肺结节是肺癌前期的主要表现之一，其边缘毛刺是判断结节良性和恶性的重要指标。此外，血管粘连型结节恶性概率最大，由于血管与结节在图像上的像素灰度值非常接近，给结节提取造成极大的干扰。因此，对血管粘连型肺结节的精确分割是当前研究的重点和难点。

目前，分割血管粘连型的肺结节的方法主要有：利用形态学方法分割肺结节，其中由于腐蚀膨胀的参数不易控制，导致欠分割或过分割，且结节的边缘毛刺也可能被消除；三维重建的方法是在二维序列图像上进行三维重建，能够还原结节的空间特征，但对图像的清晰度要求较高，计算量大，且不能有效解决血管与结节大面积相连的分割问题；基于EM和Mean-shift的肺结节分割的方法，运用梯度方向判别血管和结节的区别，能够提取出与血管粘连较为简单的结节，但不适用于血管中粘连大于等于2个结节的情况。

以上方法难以精确分割结节的主要原因在于：(1)没有充分利用二维图像信息，即未根据结节的灰度分布和几何特点进行建模，造成原始图像信息在一定程度上的缺失；(2)没有利用肺结节与血管在三维成像特点的不同，对肺结节的检测带来的很大干扰，导致检测效果有所欠缺。

发明内容

针对上述问题，本发明提出一种基于平均密度投影和平移高斯模型的肺结节分割方法，以实现对肺结节图像的精确分割。

本发明的具体实现步骤包括如下：

(1)输入待检测图像CT序列，读取CT序列的层数、层厚度和分辨率信息；

(2)根据读取的检测图像CT图像序列信息，通过平均密度投影AIP的方法，确定重建层厚度S_r和重建层间距I_r，得到待检测图像CT序列的平均密度投影AIP序列，记为B；

(3)对平均密度投影AIP序列B进行肺实质提取、阈值分割和形态学处理，得到B中第i幅图像的分割结果SA_i；利用平均密度投影AIP方法的对应关系，将分割结果SA_i对应的待检测图像CT序列，记为

其中i表示B中第i幅图像，j表示重建层序号；对

进行肺结节几何特征粗提取，得出B中第i幅图像的粗分割结果

(4)通过骨架的方法提取出粗分割结果的骨架，将

图像横坐标方向上骨架的灰度曲线记为g_ij(x)；根据该灰度曲线g_ij(x)的灰度分布，建立一维平移高斯模型f_ij(x)，计算该一维平移高斯模型f_ij(x)的参数；

(5)根据一维平移高斯模型f_ij(x)中计算所得模型参数，建立二维平移高斯模型F_ij(x，y)，并利用二维平移高斯模型F_ij(x，y)对粗分割结果

进行精确分割；

(6)对精确分割所得结果，利用圆形度的判别方法去除假阳，即去除不是肺结节的区域，获得最终的检测结果。

本方法具有以下优点：

(1)本发明利用平均密度投影AIP方法，充分利用局部三维信息，还原了血管、结节的结构特点，在弱化血管干扰的同时有效减少了处理数据量；

(2)本发明通过建立平移高斯模型实现对肺结节的精确分割，此模型不但适用于血管粘连型结节，而且同样适用于一般类型的结节，具有广泛的适用性，并可以在保证肺结节精确分割的前提下，完成对整个肺部CT图像的批处理；

(3)本发明利用结节的灰度、几何特征，首先进行粗分割，再通过建立平移高斯模型进行精确分割，这样多层次的处理方式可以保证最大限度的还原结节在待检测CT图像上的特点，充分利用了待检测CT图像信息。

附图说明

图1是本发明的整体流程图；

图2是本发明中的平均密度投影AIP成像原理示意图；

图3是本发明中的AIP处理效果图；

图4是本发明中生成的AIP序列示意图；

图5是本发明中的形态学处理子流程图；

图6是本发明中的肺结节粗提取区域对应的灰度曲面图；

图7是本发明中的一维平移高斯模型拟合效果图；

图8是本发明的疑似区域数目曲线图；

图9是本发明与传统分割结果对比图。

具体实施方式

一、基础理论介绍

1.肺结节与血管的几何特征

在待检测CT图像上，结节灰度类似服从高斯分布，形状类似圆形；血管灰度沿中心线基本不变，并以中心线上灰度值为基准，向血管壁呈逐次递减的分布，形状多为长条形。在三维空间上，结节在空间上类似球形，血管则为圆柱形；轮廓特征方面，结节边缘可能出现许多毛刺，而血管边缘则相对平滑；与血管粘连型结节的半径大于血的半径，所以结节不会被血管完全遮掩。

2.平均密度投影

平均密度投影AIP是将连续多帧图像沿射线方向上所有点灰度的平均值作为投影图像上对应点的灰度值，如图2所示，定义如下：

$\mathrm{IP}(x,y)=\frac{1}{{\mathrm{SN}}_r}\underset{k=1}{\overset{{\mathrm{SN}}_r}{\mathrm{\Σ}}}{I}_k(x,y),1\≤x\≤H,1,\≤y\≤W---1)$

其中，IP(x，y)为平均密度投影AIP图像中点(x，y)处的灰度值，SN_r为投影的层数，I_k(x，y)为原始CT序列图像中第k层图像上点(x，y)处的灰度值，H、W分别原始序列图像在x，y横、纵坐标分辨率。利用平均密度投影，可以还原肺结节与血管的几何特征，处理效果如图3所示，其中图组3(a)为平均密度投影AIP方法对血管的弱化效果图，图组3(b)为平均密度投影AIP方法对血管的几何特征恢复效果图。

二、基于平均密度投影和平移高斯模型的肺结节分割方法

参照图1，本发明的具体实施过程包括以下步骤：

步骤1.输入待检测图像CT序列，读取CT图像信息U、N_o、H、W，其中，U为原始CT序列图像层厚度，N_o为总层数，H、W分别为图像横、纵坐标分辨率。

步骤2.初始化平均密度投影AIP方法的参数，获得平均密度投影AIP序列。

对于原始CT序列图像采用AIP进行稀疏采样，其对应关系如图4所示，通过平均密度投影AIP的方法确定重建层厚度S_r和重建层间距I_r，得到待检测图像CT序列的平均密度投影AIP序列如下式：

N_r＝(N_o-SN_r)/IN_r    2)

其中，N_r为待检测图像CT序列通过AIP方法所得重建序列图像的总数，

      N_o为待检测图像CT序列图像的总数，

      SN_r为待检测图像CT序列通过AIP方法所得重建区域层数，

      SN_r＝S_r/U，U为待检测图像CT序列图像的层厚度，

      IN_r为待检测图像CT序列通过AIP方法所得的重复区域层数，

      IN_r＝I_r/U。

步骤3.根据平均密度投影AIP序列，进行待检测图像CT序列的粗分割。

3a)对平均密度投影AIP序列进行肺实质提取，利用梯度判别的方法，将肺部信息保存，排除检查床等无用信息的干扰，梯度判别的方法参考于“王伟.肺部CT图像中结节的检测算法研究[D].西安：西安电子科技大学，2007”；

3b)采用双高斯混合模型来拟合CT图像的灰度直方图，从而获得图像的整体阈值，并进行区域二值化，得到阈值分割图像，双高斯混合模型参考于“Li Huai，WangYue，Liu K J R，et al.Computerized radiographic mass detection part I：lesion site selectionby morphological enhancement and contextual segmentation[J].IEEE Transactions onMedical Imaging，2001，20(4)：289-301”；

3c)根据阈值分割图像连通区域像素点数的不同做形态学操作：

参照图5，本步骤操作如下：

3c1)当像素点数大于等于指定像素点数门限N时，做参数为A的腐蚀运算，如果运算后的小区域个数大于等于2，对每个小区域做参数为A+1的膨胀运算，其中A＝2，得到第i幅图像的分割结果SA_i，如果运算后的小区域个数小于2，将这小区域舍弃；

3c2)当像素点数小于指定像素点数门限N时，做参数为1的膨胀运算，得到第i幅图像的分割结果SA_i；

3d)利用平均密度投影AIP方法N_r＝(N_o-SN_r)/IN_r的对应关系，将分割结果SA_i对应的待检测图像CT序列，记为

其中i表示B中第i幅图像，j表示重建层序号；

3e)对

做二值化分割处理，所得序列图像

包含血管和结节的几何信息，最后将SA_i上的疑似区域分别与作逻辑‘与’操作，得到肺结节的粗分割结果

步骤4.建立一维高斯模型并计算参数。

4a)根据肺结节和血管的灰度分布，如图6所示，利用骨架的方法提取出粗分割结果

的骨架，并将在x方向骨架的灰度曲线记为g_ij(x)，如图7所示，骨架的方法参考于“刘俊涛，刘文予，吴彩华，等.一种提取物体线形骨架的新方法[J].自动化学报，2008，(06)：617-622”；

4b)建立一维平移高斯模型f_ij(x)，按如下公式确定：

其中，μ和σ分别为高斯模型的均值和方差，反映结节的灰度特征；C为平移量，反映杂波的灰度平均值，该杂波为血管或其他组织；

4c)根据所构建的一维平移高斯模型f_ij(x)，计算其中的参数：

4c1)取骨架灰度曲线g_ij(x)的最小值记为sg₀，作为骨架分割初始值；

4c2)令第k次所得的骨架分割值为sg_k，粗分割图像G_ij的灰度最大值为G_ij，max，利用sg_k将骨架曲线g_ij(x)分割成多个部分，计算G_ij，max所在区域与sg_k轴所围成的封闭区域面积，记为S_k；

4c3)将sg_k以δ步长进行递增，即：

sg_k+1＝sg_k+δ，sg_k+1≤G_ij，max    4)

计算G_ij，max所在区域与sg_k+1所围成的封闭区域面积，记为S_k+1；

4c4)将前后两次所得区域面积S_k与S_k+1的相对误差，作为判定收敛的依据：

$\left|\frac{S_k-S_{k+1}}{S_k}\right|\≤\mathrm{\ξ}---5)$

其中，ξ为收敛系数，如果满足上式则令C＝sg_k，如果不满足，则转到步骤4c2)进行迭代；

4c5)由骨架灰度曲线g_ij(x)的灰度值减去平移量C，将此所得结果利用非线性最小二乘法求解，得出一维平移高斯模型f_ij(x)中高斯模型的均值μ和高斯模型的方差σ。

步骤5.构建二维平移高斯模型，并进行精确分割。

5a)根据一维平移高斯模型f_ij(x)的参数μ、σ和C，构建二维平移高斯模型F_ij(x，y)：

$F_{\mathrm{ij}}(x,y)=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp \{-\frac{{(x-\mathrm{\μ})}^2+{(y-y_0)}^2}{2{\mathrm{\σ}}^2}\}+C---6)$

其中y₀为骨架曲线g_ij(x)上横坐标为μ时，在粗分割结果对应的纵坐标值，μ为高斯模型的均值，σ为高斯模型的方差，C为平移量；

5b)根据二维平移高斯模型F_ij(x，y)，对粗分割结果

进行精确分割，精确分割的方法按如下公式确定：

|G_ij(x，y)-F_ij(x，y)|＜ε    7)

其中ε为判别门限，G_ij(x，y)为粗分割结果对应坐标的灰度值，满足上式保留此点，否则去除，将处理后的图像记为精确分割图像D_ij(x，y)。

步骤6.利用圆形度对精确分割图像D_ij(x，y)进行计算并判断，去除假阳，即去除不是肺结节的区域，将最终检测结果在CT图像上输出显示。

6a)圆形度定义为：

$\mathrm{CD}=\frac{4\mathrm{\πS}}{P^2}---8)$

其中，S为所在区域面积，P为所在区域周长。圆形度定义参考于“郭浩.低剂量CT孤立肺结节检测方法研究[D].西安：西安电子科技大学，2010”；

对精确分割图像计算圆形度CD，如果满足判别式CD≥η，则保留精确分割图像区域，反之去除此区域，其中η为特定门限。

6b)将判别后所得精确分割图像D_ij(x，y)边界提取出来，在待检测序列CT图像上输出显示。

6c)返回步骤3，进行重复操作，直到对待检测序列CT的所有图像处理完为止。

本发明的效果可通过以下实验进一步说明：

1.实验条件

本实验采用层厚为1～2mm，间距0.5～1mm，分辨率为512×512像素的薄层待检测CT图像，此外，医生对实验所用数据的30个血管粘连型肺结节进行了标记，作为分割结果的参考。在本发明中，平均密度投影AIP方法重建层厚度取S_r＝2×I_r。

判断分割性能，采用面积交迭度AOM作为分割效果的评价指标，定义为：

$\mathrm{AOM}(A,B)=\frac{S\left(\mathrm{AIB}\right)}{S\left(\mathrm{AUB}\right)}\×100\%---9)$

其中，AOM为面积交迭度，A为医生标注的图像，B为分割结果图，S(·)表示对应区域的像素点数，AOM取值越大表明分割效果越好。AOM定义参考于“艾海明，吴水才，高宏建，等.基于图论的肝肿瘤CT图像自动分割方法[J].北京工业大学学报，2010，36(4)：572-576”。

2.仿真内容与结果

实验一：比较不同的平均密度投影AIP方法重构厚度S_r对肺结节提取的影响，结果如表1所示。在医生标记的疑似区域全部检测到的情况下，采用不同S_r重构，通过圆形度的检测，比较各层疑似区域的检测结果，如图8所示。

实验二：选取各种粘连形态、尺寸大小以及良/恶性的肺结节的待检测原始序列CT图像，对本发明、区域生长法和形态学法进行比较，其中本发明的平均密度投影AIP重建层厚度取S_r＝7mm，结果如图9所示，其中图组9(a)为待检测CT图像，图组9(b)为用区域生长法对各种结节的提取结果，图组9(c)为用形态学方法对各种结节的提取结果，图组9(d)为用本发明对各种结节的分割结果。

表1不同重建层厚层间距检测结果

由表1所得数据可以看出，在S_r＝5mm时，S_r过薄，血管信息没得到很好的弱化；在S_r＝9mm时，有三个点未检测到，根据观察，是由于这三个结节的面积较小且灰度值较低，当S_r过厚时，结节的信息在很大程度上被弱化了，不利于提取结节；在S_r＝7mm时，检测效果较好。

从图8中可以看出，采用平均密度投影AIP处理后的序列图像的疑似区域个数明显小于待检测序列图像的检测结果；在S_r＝5mm时，血管的空间几何特征，在平均密度投影AIP图像没有很好的体现出来，血管仍然可能以类似结节的圆形出现；在S_r＝7mm时，假阳大幅度下降，效果明显。

从图9中可以看出，区域生长法不能有效的将结节与血管分割开；形态学方法对于面积较小的良性结节，出现了漏检现象，对于恶性结节，丢失了边缘毛刺信息，效果不理想；本发明对良/恶性结节、血管粘连多个结节，都具有很好的分割效果，表现出较强的鲁棒性。

对以上三种算法均能检测到结节的实验结果统计其AOM，如表2所示。

表2不同方法提取精度

由表2所得数据可以看出，本发明的提取结果与专业医师标记的面积交迭度达到93％，高出区域生长法和形态学法10～20％，具有更好的分割效果和分割精度。

Claims (3)

1.一种基于平均密度投影和平移高斯模型的肺结节分割方法，包括如下步骤：

（1）输入待检测图像CT序列，读取CT序列的层数、层厚度和分辨率信息；

（2）根据读取的检测图像CT图像序列信息，通过平均密度投影AIP的方法，确定重建层厚度S_r和重建层间距I_r，得到待检测图像CT序列的平均密度投影AIP序列，记为B；

（3）对平均密度投影AIP序列B进行肺实质提取、阈值分割和形态学处理，得到B中第i幅图像的分割结果SA_i；利用平均密度投影AIP方法的对应关系，将分割结果SA_i对应的待检测图像CT序列，记为

其中i表示B中第i幅图像，j表示重建层序号；对

进行肺结节几何特征粗提取，得出B中第i幅图像的粗分割结果

其中SN_r为投影的层数；

（4）通过骨架的方法提取出粗分割结果

的骨架，将

图像横坐标方向上骨架的灰度曲线记为g_ij(x)；

(5)根据该灰度曲线g_ij(x)的灰度分布，建立一维平移高斯模型f_ij(x):

$f_{\mathrm{ij}}\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}{e}^{-\frac{{(x-\mathrm{\μ})}^2}{2{\mathrm{\σ}}^2}}+C,$

其中，μ和σ分别为高斯模型的均值和方差，反映结节的灰度特征；C为反映杂波的平移量，即血管或其他组织的灰度平均值;

(6)计算该一维平移高斯模型f_ij(x)的参数:

（6a）取骨架灰度曲线g_ij(x)的最小值记为sg₀,作为骨架分割初始值；

（6b）令第k次所得的骨架分割值为sg_k，粗分割图像G_ij的灰度最大值为G_ij,max，利用sg_k将骨架曲线g_ij(x)分割成多个部分，计算G_ij,max所在区域与sg_k轴所围成的封闭区域面积，记为S_k；

（6c）将sg_k以δ步长进行递增，即：

sg_k+1＝sg_k+δ,sg_k+1≤G_ij,max

计算G_ij,max所在区域与sg_k+1所围成的封闭区域面积，记为S_k+1；

（6d）将前后两次所得区域面积S_k与S_k+1的相对误差，作为判定收敛的依据：

$\left|\frac{S_k-S_{k+1}}{S_k}\right|\≤\mathrm{\ξ}$

其中，ξ为收敛系数，如果满足上式则令C＝sg_k，如果不满足，则转到步骤（6b）进行迭代；

（6e）由骨架灰度曲线g_ij(x)的灰度值减去平移量C，再利用非线性最小二乘法求得所述参数μ和σ；

（7）根据一维平移高斯模型f_ij(x)中计算所得模型参数，建立二维平移高斯模型F_ij(x,y)，并利用二维平移高斯模型F_ij(x,y)对粗分割结果

进行精确分割,

该二维平移高斯模型F_ij(x,y)表示如下：

$F_{\mathrm{ij}}(x,y)=\frac{1}{2{\mathrm{\π\σ}}^2}\exp \{-\frac{{(x-\mathrm{\μ})}^2+{(y-y_0)}^2}{{2\mathrm{\σ}}^2}\}+C$

其中y₀为骨架曲线g_ij(x)上横坐标为μ时，在粗分割结果

对应的纵坐标值；

（8）对精确分割所得结果，利用圆形度的判别方法去除假阳，即去除不是肺结节的区域，获得最终的检测结果。

2.根据权利要求1所述的分割方法，其中步骤（2）所述的通过平均密度投影AIP的方法确定重建层厚度S_r和重建层间距I_r，得到待检测图像CT序列的平均密度投影AIP序列，按如下公式确定：

N_r=(N_o-SN_r)/IN_r

其中，N_r为待检测图像CT序列通过AIP方法所得重建序列图像的总数，

N_o为待检测图像CT序列图像的总数，

SN_r为待检测图像CT序列通过AIP方法所得重建区域层数，

SN_r＝S_r/U，U为待检测图像CT序列图像的层厚度，

IN_r为待检测图像CT序列通过AIP方法所得的重复区域层数，

IN_r=I_r/U。

3.根据权利要求1所述的方法，其中步骤（7）所述的利用二维平移高斯模型F_ij(x,y)对粗分割结果

进行精确分割，按如下判别式确定：

|G_ij(x,y)-F_ij(x,y)|＜ε

其中ε为判别门限，G_ij(x,y)为粗分割结果对应坐标的灰度值，满足上式保留此点，否则去除。

Images