CN104992445B

CN104992445B - 一种ct图像肺实质的自动分割方法

Info

Publication number: CN104992445B
Application number: CN201510426463.9A
Authority: CN
Inventors: 王兵; 顾潇蒙; 顾力栩; 张欣; 田学东
Original assignee: Hebei University
Current assignee: Hebei University
Priority date: 2015-07-20
Filing date: 2015-07-20
Publication date: 2017-10-20
Anticipated expiration: 2035-07-20
Also published as: CN104992445A

Abstract

本发明提供了一种CT图像肺实质的自动分割方法。本发明通过两次随机游走算法对CT图像进行分割以得到准确的肺实质，第一次采用随机游走算法分割，得到近似肺实质掩膜；第二次采用随机游走算法对肺边界处缺陷进行修补，分割后得到准确的肺实质结果。采用随机游走算法分割时所设置的种子点是通过大津阈值和数学形态学等方法快速且自动地获取，无需手工标定，大大减少了医生的工作量和操作时间。本发明提出“两次选取种子点，两次分割”的过程，是由粗到细的自动分割过程，最终分割结果对初始种子点选择的依赖度降低，确保了分割结果的准确性、完整性、实时性和鲁棒性。本发明得到国家自然科学基金资助（编号：61375075）。

Description

一种CT图像肺实质的自动分割方法

技术领域

本发明涉及一种肺实质分割方法，具体地说是一种基于随机游走算法的CT图像肺实质的自动分割方法。

背景技术

目前，肺癌是发病率和死亡率增长最快，对人类健康和生命威胁最大的恶性肿瘤之一。根据世界卫生组织2014年5月发布的最新研究资料表明，肺癌在全球前十位死亡原因中占第五位，居癌症死亡的首位。而我国以吸烟人数广大、空气污染严重等原因，肺癌死亡率增加明显，每年平均以4.4％的速度上升，自1996年以来肺癌上升为我国癌症患者的第一杀手。如果肺癌能在早期被诊断和治疗，其五年生存率可达40％～70％。

计算机断层扫描(Computed Tomography，CT)能够提供高清晰度的图像，并且为图像中各组织提供很高的对比度，因此被广泛用于胸肺部疾病检测诊断中。临床表明肺结节是肺癌在CT图像中最为常见的表现形式，而肺结节形状多样，大小不同，分布位置也不固定，易与其他组织紧密连接，密度与肺部血管相近，在CT图像中大多表现为圆形或近似圆形的致密斑点，仅凭人眼很难区别。而且，一次肺部CT扫描将产生近百张CT图像，其中有结节的图像仅有几张。对于所有影像数据进行分析非常枯燥繁琐，即使有经验的医生也难免因诊断疲劳导致漏诊或者误诊，特别是在数据量急剧增长的情况下。因而，肺部疾病的计算机辅助诊断(Computer-Aided Diagnosis，CAD)系统应运而生，在此过程中，肺实质的分割是实现自动量化诊断的前提，同时也是进行肺结节检测和肺功能评估的基础，其分割结果的好坏直接影响到后续的分析过程。因此，自动、准确、快速地分割肺实质可为临床治疗和病理学研究提供可靠依据，具有重要的现实意义和临床价值。

近年来，国内外已提出很多肺实质分割方法，主要基于阈值法、聚类法、区域生长、边缘检测、主动轮廓模型等，这些方法能够充分利用肺部CT图像影像学的特点，针对正常的肺组织取得有效的分割结果，但对于病灶粘连胸膜、肺门区域存在高密度血管、肺区域边界模糊、病灶或组织引起的肺区断开以及左右肺连接在一起等复杂情况，采用上述分割方法均难以获得理想的分割效果。

发明内容

本发明的目的就是提供一种CT图像肺实质的自动分割方法，以解决现有肺实质分割方法无法对复杂情况下的CT图像肺实质进行有效分割的问题。

本发明是这样实现的：一种CT图像肺实质的自动分割方法，包括如下步骤：

a、对CT图像进行预处理，获得目标掩膜和背景掩膜；

b、在目标掩膜对应的CT图像中肺实质区域内选取目标种子点，在背景掩膜对应的CT图像中胸骨部位以及两肺之间的区域内选取背景种子点；

c、依据步骤b中所选取的目标种子点和背景种子点，采用随机游走算法对CT图像进行分割，获取近似肺实质掩膜，在CT图像中与近似肺实质掩膜的边界对应的部位存在缺陷；

d、调整目标种子点和背景种子点位置：采用半径为r₁的结构元对近似肺实质掩膜做膨胀运算，得到膨胀后的肺实质掩膜；采用半径为r₂的结构元对膨胀后的肺实质掩膜做腐蚀运算，得到腐蚀后的肺实质掩膜；r₁＜r₂；在CT图像中选取与膨胀后的肺实质掩膜边界对应的点作为背景种子点，在CT图像中选取与腐蚀后的肺实质掩膜边界对应的点作为目标种子点；

e、依据步骤d中调整后的目标种子点和背景种子点，采用随机游走算法对CT图像进行分割，修补CT图像中与近似肺实质掩膜边界对应的缺陷，得到精确的肺实质图像。

步骤c和步骤e中采用随机游走算法对CT图像进行分割的具体过程为：

定义一个带权无向图G＝(V,E)，用图G＝(V,E)给CT图像建模，其中V是图中顶点的集合，v∈V，v对应CT图像中的每个像素，E是图中任意两个顶点间的无向边的集合，e对应CT图像中像素的连接关系；连接两个顶点v_i和v_j的边e用e_ij表示，边e_ij的权值用w_ij表示；定义顶点v_i的度为d_i＝∑w_ij，d_i是连接顶点的所有边的权值之和；通过定义边的权值，构建Laplace矩阵，求解Dirichlet问题。

边e_ij的权值w_ij定义为：

式(1)中，g_i表示顶点v_i的灰度值；h_i表示顶点v_i的坐标值；b₁和b₂是自由参数，b₁表示像素灰度差特征的权重，b₂表示像素几何距离特征的权重。

图G＝(V,E)的Laplace矩阵L定义为：

L可通过公式L＝A^TCA进行构建；其中，A是边与顶点的m×n关联矩阵，定义为：

C是图G＝(V,E)的m×m本构矩阵，定义为对角矩阵，对角矩阵中的对角元素是相对应的边的权值。

Dirichlet积分的离散形式可表示为：

求解Dirichlet问题即是求解式(2)最小化的离散调和函数x；

将CT图像的顶点V分为两个集合V_M和V_U，集合V_M中的顶点为标记点，集合V_U中的顶点为未标记点；V_M∪V_U＝V，V_M∩V_U＝φ，则式(2)可分解为：

通过对D[x_U]关于x_U求微分来寻找其极值点：

假设表示顶点v_i属于标签s的概率，标记点的标签集合定义为Q(v_j)＝s，其中s∈Z，0<s≤K，Z为整数集合，K为标签s的数量，因此有：

则求解离散Dirichlet问题即求解：L_Ux^s＝-B^Tm^s或L_UX＝-B^TM，其中

步骤a具体包括如下步骤：

a1、采用高斯模板对CT图像进行滤波处理，以降低噪声；

a2、采用大津阈值分割技术对降噪后的CT图像进行二值化处理，得到背景掩膜和目标掩膜；背景掩膜中包含有胸骨及心脏，目标掩膜中包含有肺实质、气管、支气管和胃气；

a3、去除目标掩膜中的气管、支气管和胃气，并在CT图像中相应的位置做删除标记。

步骤d中，9mm≤r₁≤10mm，11mm≤r₂≤12mm。

本发明针对复杂情况下肺实质的分割问题，提供了一种基于随机游走(RandomWalk)算法的CT图像肺实质的自动分割方法。在本发明中，首先根据胸部组织解剖学及其CT图像的影像学特征，在肺实质及其周围组织内分别确定目标种子点和背景种子点位置；然后使用随机游走算法对CT图像进行分割，提取近似肺实质掩膜；接着对近似肺实质掩膜实施数学形态学运算，来进一步调整目标种子点和背景种子点的标定位置，使其适合具体的复杂情况；最后再次使用随机游走算法分割CT图像，得到最终的肺实质分割结果。采用本发明可解决以下几种复杂情况下肺实质的分割问题：胸膜结节及肺门附近高密度血管的存在造成分割结果有遗漏病灶的、两肺距离较近或连接在一起引起分割结果错误的、病灶或组织造成的肺区断开、肺边界模糊使得分割结果不完整等复杂情况。

本发明采用了两次随机游走算法对CT图像进行分割以得到准确的肺实质，第一次采用随机游走算法分割CT图像，得到近似肺实质掩膜；第二次采用随机游走算法对CT图像中对应近似肺实质掩膜边界部位的缺陷进行修补，分割后得到准确的肺实质分割结果。随机游走算法的权值计算，包括了像素的灰度差及像素间的几何距离，使得各边的权值不仅反映图像像素的灰度变化，还反映像素间距离，提高了算法对肺部图像的适用性和有效性。随机游走算法本身无需迭代，计算速度较快，在强噪声、弱边界或没有显著边界的情况下也能获得较好的分割效果，为上述复杂情况的肺实质分割提供有效可行的方案。

本发明中种子点是通过大津阈值和数学形态学等方法快速且自动地获取，无需手工标定，大大减少了医生的工作量和操作时间，确保了分割的实时性；根据随机游走算法的特点和肺部图像的特征，提出“两次选取种子点，两次分割”的过程，是由粗到细的自动分割过程，使最终分割结果对初始种子点选择的依赖度降低，确保了分割结果的准确性、完整性、实时性和鲁棒性，采用本发明所得分割结果和分割时间均可满足临床需求。

附图说明

图1是本发明的方法流程图。

图2是去除目标掩膜中气管、支气管的示意图。

图3是去除目标掩膜中胃气的示意图。

图4是首次分割时设置的种子点的示意图。

图5是首次使用随机游走算法分割的结果示意图。

图6是调整后的种子点的位置示意图。

图7是二次使用随机游走算法分割的最终结果图。

具体实施方式

本发明实施例所用计算机的软硬件条件是：Dual-Core CPU E5800@3.20GHz，显卡为NVIDIA GeForce GT 430，内存2.0GB，操作系统为Window 2007，软件编程语言使用MATLAB。本发明得到国家自然科学基金资助(编号：61375075)。

如图1所示，本发明包括如下步骤：

步骤S1，图像预处理。

本步骤对原始CT图像进行预处理，主要包括高斯平滑去噪，采用大津阈值分割技术对图像进行二值化处理，获得目标掩膜和背景掩膜。具体如下：

(1)高斯平滑去噪：使用3×3高斯模板对图像进行滤波处理，达到降低噪声的目的。

(2)获得目标掩膜和背景掩膜：采用大津阈值分割技术(最大类间方差方法)对上述滤波后的图像进行二值化处理。

由于胸廓的密度高，在CT图像中呈现高灰度值的亮色，肺实质包含在其中，呈现低灰度值的暗色，胸廓之外主要是灰度值更低的人体之外的暗色区域。因此，可以先采用大津阈值分割方法对图像做二值处理。由于胸廓是闭合的，并与两肺之间的高灰度值区域有近似的灰度值，二值化后成为一个区域，通过搜索最大连通区域可以找到这部分区域，将这部分区域的灰度值置255(白色)，这部分区域即是胸骨及心脏等面积最大的高密度连通区域，二值化后变成一个白色连通区域，称为背景掩膜B_M。将背景掩膜B_M之外部分全部赋值为0(黑色)。背景掩膜B_M之内包裹着有肺实质、与肺实质不连通的气管、支气管以及胃气等低密度的区域，将这些区域的灰度值赋值为0(黑色)，称为目标掩膜T_M。图2(a)示出了包含有气管的目标掩膜，图2(b)示出了包含有支气管的目标掩膜，图3(a)示出了包含有胃气的目标掩膜。

(3)去除目标掩膜中的气管、支气管和胃气，并在CT图像中相应的位置做删除标记。

去除气管、支气管方法：计算出包含目标掩膜T_M的最小矩形，称为目标窗口T_W，如图2(a)和图2(b)中的黑色矩形框。在目标窗口T_W的中心位置设置一个面积小于目标窗口的矩形，称为中心窗口C_W，如图2(a)和图2(b)中的白色矩形框。由于孤立的气管和支气管掩膜位于两肺实质掩膜之间，通常面积小于等于2000像素对应的面积，本实施例中取中心窗口C_W的面积为目标窗口T_W面积的四分之一，就可将气管和支气管掩膜包含其中。删除中心窗口C_W内面积小于等于2000像素对应的气管、支气管掩膜，同时在原图中气管及支气管位置上做删除标记。将图2(a)中气管删除后的目标掩膜如图2(c)所示，将图2(b)中支气管删除后的目标掩膜如图2(d)所示。

去除胃气的方法：胃气在底层切片，肺实质呈现月牙形，其质心不在肺实质区域内。把掩膜图像均分成四个象限，如图3(a)所示。对Ⅲ象限和Ⅳ象限(下面两个象限)中的目标掩膜分别按面积大小排序。如果在Ⅲ象限或Ⅳ象限中找到一个大的目标掩膜，其质心不在该区域内部，则认为此图像为肺底层切片；那么，判断Ⅰ象限和Ⅱ象限(上面两个象限)中是否有目标掩膜，如果有，则认为是胃气掩膜，将它们删除，同时在原图中做删除标记。将图3(a)中胃气删除后的目标掩膜如图3(b)所示。

步骤S2，设置分割算法的目标种子点和背景种子点。

去除目标掩膜中的气管、支气管、胃气之后，在目标掩膜对应的CT图像中肺实质区域内选取目标种子点，在背景掩膜对应的CT图像中胸骨部位以及两肺之间的区域内选取背景种子点。

目标种子点的选取：在目标掩膜对应的CT图像中，至少有两个不连通的区域(左右两肺区域)，在每一连通区域内随机设置N个目标种子点(N通常取值≥2)，选择的目标种子点的灰度值应该小于此连通区域的平均灰度值-Δ(Δ是个正整数，通常取值5)，使得目标种子点不设置在血管或是肿瘤上。如图4所示，图4(a)为掩膜(包括目标掩膜和背景掩膜)，图4(b)为与图4(a)相对应的CT图像。在图4(a)和图4(b)的目标掩膜及对应的CT图像中，具有两个不连通的左右肺区，在每一肺区内随机选两个目标种子点，图中所选的目标种子点以白色三角形示出。

背景种子点的选取：如图4(a)所示，在目标窗口T_W内中心位置做一个并行与窗口高的直线，在直线上任意位置选择两个未做删除标记的点，作为两个背景种子点(图中以黑色圆圈示出)；并在区域内均匀设置四个点，对应原图(即CT图像)中的四个背景种子点。即背景种子点是在目标窗口之外的胸骨上、下分别取一个，左、右各取一个，两肺之间取两个，总共设置六个背景种子点。其他实施例中也可以设置八个或十个等。

随机游走算法对种子点的个数不敏感，但均匀设置的种子点会提高分割速度。

步骤S3，第一次用随机游走算法获得近似肺实质掩膜。

根据上一步设置的目标种子点和背景种子点，第一次用随机游走算法对CT图像进行分割，获得粗分割结果，获取近似肺实质掩膜。粗分割后可将图5(a)中近乎连在一起的两个肺区域正确地区分开，分开后可得到如图5(b)和图5(c)所示结果，但是，在第一次粗分割后，对于有病灶或血管等的复杂情况下，会在肺边界处形成凹陷，如图5(b)所示；图5(c)中示出了在肺边界处形成了缺口，肺区表现出了不完整的情形。肺边界处的凹陷或缺口等均构成了肺边界的缺陷。

本步骤中随机游走算法的具体运算过程为：

定义一个带权无向图G＝(V,E)，用图G＝(V,E)给CT图像建模，其中V是图中顶点(对应CT图像中的像素点)的集合，v∈V，v对应CT图像中的每个像素，E是图中任意两个顶点间的无向边的集合，e对应CT图像中像素的连接关系；连接两个顶点v_i和v_j的边e用e_ij表示，边e_ij的权值用w_ij表示；定义顶点v_i的度为d_i＝∑w_ij，d_i是连接顶点的所有边的权值之和；因为权值w_ij代表随机游走者从顶点v_i转移到顶点v_j的概率，因此必须有w_ij>0，无向图即要求w_ij＝w_ji。通过定义边的权值，构建Laplace(拉普拉斯)矩阵，随机游走者首次到达某一种子点的概率的求解过程可转化为求解Dirichlet(狄利克雷)问题。

(1)边的权值的定义：权值表示相邻像素之间的相似度，本发明使用像素间的灰度差和几何距离作为度量特征，将边的权值定义为：其中，g_i表示顶点v_i的灰度值；h_i表示顶点v_i的坐标值；b₁和b₂是自由参数，b₁表示像素灰度差特征的权重，b₂表示像素几何距离特征的权重。本实施例中令β₁＝90，β₂＝1。

(2)Laplace矩阵的构建，具体为：

图G＝(V,E)的Laplace矩阵L定义为：

(3)Dirichlet(狄利克雷)问题的求解，具体为：

Dirichlet积分的离散形式可表示为：

求解Dirichlet问题即是求解上式最小化的离散调和函数x。由于L是半正定矩阵，故D[x]存在唯一的极小值点。将CT图像的顶点V分为两个集合V_M和V_U，集合V_M中的顶点为标记点，集合V_U中的顶点为未标记点；V_M∪V_U＝V，V_M∩V_U＝φ，则上式可分解为：

通过对D[x_U]关于x_U求微分来寻找其极值点：

在本发明中，目标种子点和背景种子点可看作两个不同的标签(即本发明中标签数量K为2)，在未运算之前，目标种子点和背景种子点是标记点，其他点均是未标记点，对每一未标记点进行随机游走运算，若未标记点首次到达某一种子点的概率最大，则该点与该种子点属于同一标签。举例来说，对于其中一个未标记点，使其随机游走到所标定的种子点处，计算出该点首次到达种子点的最大概率，则对应最大概率的种子点与该点为同一标签，将该点标记为与其具有相同标签的目标种子点或背景种子点，该点也就成为了标记点。通过对所有的未标记点进行运算，即可将所有未标记点一一标记、分类，标记、分类完毕后也就实现了对图像的分割。

步骤S4，调整目标种子点和背景种子点位置。

在缺陷附近调整种子点的位置可以修补边界缺陷。对近似肺实质掩膜实施形态学的膨胀和腐蚀操作，目的是根据实际情况重新自动调整种子点的位置，并再次使用随机游走算法分割时，可以修补边界的缺陷，得到准确的分割效果。要想修补肺边界的缺陷，必须在缺陷内选择至少一个目标种子点，在缺陷外选择至少一个背景种子点。对近似肺实质掩膜实施形态学膨胀后，区域边界对应的CT图像中边界作为背景种子点；然后再对膨胀后的掩膜实施形态学腐蚀操作，腐蚀后的区域略小于近似肺实质掩膜，其边界要与肺区域凹陷部分有重叠，腐蚀后的区域边界对应的CT图像中边界作为目标种子点。

例如可以采用半径为9mm的结构元对近似肺实质掩膜做膨胀运算，得到膨胀后的肺实质掩膜，再用半径为11mm的结构元对膨胀后的肺实质掩膜做腐蚀操作，得到腐蚀后的肺实质掩膜。如图6(a)所示，在CT图像中选取与膨胀后的肺实质掩膜边界对应的点作为背景种子点(黑色边线)，在CT图像中选取与腐蚀后的肺实质掩膜边界对应的点作为目标种子点(白色边线)。图6(b)为图6(a)中黑色矩形框的局部放大示意图。这些种子点均匀地分布在肺实质内外两侧，分割后既可以修补肺实质边界凹陷，又不会造成区域失真。若两肺间最短距离小于20mm，则左右肺依次用此方法调整种子点位置，否则可以对两个肺区域同时调整种子点位置。

步骤S5，再次用随机游走算法修补边界缺陷。

第二次使用随机游走算法修补边界缺陷：根据调整后的目标种子点和背景种子点，再次使用随机游走算法对CT图像进行分割，则修补了边界缺陷，得到精确的分割结果图，如图7所示，图7(a)即是将图5(b)中肺边界具有的凹陷修补后的示意图，图7(b)即是将图5(c)中肺边界具有的缺口修补后的示意图，修补后肺实质完整了。

Claims

1.一种CT图像肺实质的自动分割方法，其特征是，包括如下步骤：

a、对CT图像进行预处理，获得目标掩膜和背景掩膜；

2.根据权利要求1所述的CT图像肺实质的自动分割方法，其特征是，步骤c和步骤e中采用随机游走算法对CT图像进行分割的具体过程为：

3.根据权利要求2所述的CT图像肺实质的自动分割方法，其特征是，边e_ij的权值w_ij定义为：

式(1)中，g_i表示顶点v_i的灰度值，g_j表示顶点v_j的灰度值；h_i表示顶点v_i的坐标值，h_j表示顶点v_j的坐标值；β₁和β₂是自由参数，β₁表示像素灰度差特征的权重，β₂表示像素几何距离特征的权重。

4.根据权利要求3所述的CT图像肺实质的自动分割方法，其特征是，图G＝(V,E)的Laplace矩阵L定义为：

L可通过公式L＝A^TCA进行构建；其中，A是边与顶点的m×n关联矩阵，m是图G中的边数，n是图G中的顶点数，定义为：

5.根据权利要求4所述的CT图像肺实质的自动分割方法，其特征是，Dirichlet积分的离散形式可表示为：

<mrow> <mi>D</mi> <mo>&lsqb;</mo> <mi>x</mi> <mo>&rsqb;</mo> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <msup> <mrow> <mo>(</mo> <mi>A</mi> <mi>x</mi> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mi>C</mi> <mrow> <mo>(</mo> <mi>A</mi> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <msup> <mi>x</mi> <mi>T</mi> </msup> <mi>L</mi> <mi>x</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <munder> <mo>&Sigma;</mo> <mrow> <msub> <mi>e</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>&Element;</mo> <mi>E</mi> </mrow> </munder> <msub> <mi>w</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

求解Dirichlet问题即是求解式(2)最小化的离散调和函数x；

通过对D[x_U]关于x_U求微分来寻找其极值点：L_Ux_U＝-B^Tx_M；其中，x_M为标记点的灰度值，x_U为未标记点的灰度值，L_M为标记点在拉普拉斯矩阵中的对应值，L_U为未标记点在拉普拉斯矩阵中的对应值，B为图G中边界点在拉普拉斯矩阵中的对应值；

6.根据权利要求1所述的CT图像肺实质的自动分割方法，其特征是，步骤a具体包括如下步骤：

a1、采用高斯模板对CT图像进行滤波处理，以降低噪声；

a2、采用大津阈值分割技术对降噪后的CT图像进行二值化处理，得到背景掩膜和目标掩膜；背景掩膜中包含有胸骨及心脏，目标掩膜中包含有肺实质、气管、支气管和胃气；所述胃气位于CT图像的底层切片中并呈现非月牙形结构；

7.根据权利要求1所述的CT图像肺实质的自动分割方法，其特征是，步骤d中，9mm≤r₁≤10mm，11mm≤r₂≤12mm。