CN110047085B - 一种针对肺ct图像阈值分割结果的肺膜粘连结节区域精确修复方法 - Google Patents

Abstract

本发明公开了一种针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，通过对肺部二值图像进行SIFT特征检测获得图像中的特征点，并对上述获得的非边界特征点计算其最近邻替代点，并使用计算出的替代点替换原非边界特征点后，获得对应的边界特征点集合；针对上述边界特征点集合中的每一个边界特征点，基于特定的规则提取支撑边界；并根据计算得到的支撑边界图像中各支撑边界的联通性对其分段，获得图中彼此独立的支撑边界；对每个支撑边界分段,使用基于弯曲度编码对待修复的边界进行识别；对识别出的部分使用基于GVF的方法驱动光滑轮廓曲线以修复凹陷区域及边界。本发明应用本发明方法可以精确修复阈值分割结果中遗漏的肺膜黏连结节区域，将肺膜结节区域重新纳入肺实质区域中，从而保障了肺部CT图像病变组织分析算法的准确性。

Description

一种针对肺CT图像阈值分割结果的肺膜粘连结节区域精确修 复方法

技术领域

本发明属于肺部CT图像分析处理技术领域，具体涉及一种针对肺CT图像阈值分割结果的肺膜粘连结节区域精确修复方法。

背景技术

近二十年以来，我国肺癌发病人数急剧增加。另外，因为肺癌在所有癌症中致死率较高，它成为危害我国居民健康的重要疾病。对于肺癌，使用CT设备进行影像学检查是当前重要的诊断手段。因此，针对肺部CT图像的识别技术是值得重视的。

随着计算机硬件计算能力的提高以及相关算法理论的完善，针对肺CT图像使用计算机硬件和算法辅助检测病变部位的方案已越来越成熟。其中，从肺部CT图像中精确提取肺区域成为节省计算资源、提高识别准确性的重要和首要手段。

在图像处理算法中，提取肺部区域的基本依据是CT图像中不同组织间的HU值范围存在差异，故大部分肺区域自动提取算法常设定一个经验阈值或通过计算获得阈值。

使用阈值遍历所有的像素点后，根据每一个像素点的CT值与阈值的大小关系，将像素点划分为肺区域或者其他区域两类。待所有点被处理完毕后，将提取到的像素点组成肺区域，并送入下一个处理流程识别其中的病变组织。

但使用上述基于阈值的方法提取肺区域有时存在问题，特别是在与肺膜粘结的肺结节上该问题更为严重。

该问题表现在CT图像上与肺膜粘连的肺结节和周围其他组织具有相近的CT值，因此基于阈值的方法常常将这些肺膜粘连型肺结节排除在肺区域之外。

但是，肺膜粘连肺结节具有较高的阳性肺结节可能性。因此，若将这些结节区域排除在肺区域外，通过算法识别出的病变部位中将缺少这部分区域，从而对算法的可靠性和准确性产生影响。

另外，一些肺边界修复算法将两个端点的连接直线段作为修复后的边界。这种方案获得的修复边界太过平直，与周围正常的肺边界存在差异而且不符合正常肺边界的形态学特征。同时，使用直线段作为边界后，使肺膜结节外部边界形状发生改变，从而影响了对肺膜结节各种形态学、统计学特征的计算结果，最终导致后续识别病变算法中出现结果识别错误。

所以，有必要开发新的算法将肺膜结节区域再次处理，以精确地将肺膜结节所在的区域重新纳入肺区域中。

综上所述，现有肺区域提取分割算法存在的主要问题是：

基于阈值提取后的肺实质区域中常常缺失肺膜粘连肺结节所在的区域。

现有技术采用直线修复的方式，修复的结果与周围正常的肺边界存在差异、不符合正常肺边界的形态学特征。

现有的部分肺结节区域处理算法存在依赖算法初始条件的问题。在初始条件较差时，需要二次或再次处理，使算法存在处理效率不高的问题。

部分肺结节区域处理算法对病变区域的位置处理方案不完善，存在对病变位置识别精度不高的问题。

发明内容

为了将二值图像中的肺膜结节区域准确地重新纳入肺实质区域中，本发明提供了一种基于弯曲度编码和梯度向量流的粘连肺结节区域修复算法。

为了达到上述目的，本发明采用的技术方案如下：

本算法是这样实现的，一种针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，包括以下步骤：

S1、获取肺部区域二值图像；

S2、对上述二值图像进行SIFT特征检测获得图像中的特征点，并对上述获得的非边界特征点计算其最近邻替代点，并使用计算出的替代点替换原非边界特征点后，获得对应的边界特征点集合；

S3、针对边界特征点集合中的每一个边界特征点，基于特定的规则提取支撑边界；并根据计算得到的支撑边界图像中各支撑边界的联通性对其分段，获得图中彼此独立的支撑边界；

S4、对上述每一段支撑边界计算其最小外接矩形并基于综合码、差值码、连续码、弯曲度编码对待修复的边界进行识别；

S5、对待修复的边界进行修复。

进一步地，所述步骤S5中，对识别为凹陷区域的部分使用基于GVF-Snake的方法驱动光滑轮廓曲线以修复凹陷区域以及相关边界。

进一步地，所述步骤S1中，输入经过阈值方法从CT图像中直接提取到的肺部区域二值图像。

进一步地，步骤S2具体包括：

S21、首先基于特征点检测、特征点描述、特征点匹配和消除错配点的SIFT特征点检测步骤后，得到二值图像中的特征点位置；

S22、对上述检测到的每一个SIFT特征点，若其不在肺部边界上则寻找距离该特征点最近的边界点作为其替代点，从而将特征点集合转变为边界特征点集合。

其中步骤S21中，SIFT特征点是尺度空间中满足位置、尺度和旋转不变量的极值点，

首先，建立尺寸空间以及相关极值检测：基于二维高斯函数以及高斯金字塔建立尺度空间，使用高斯差分算子检测尺度空间中的关键点。

然后，使用拟合三维二次方程精确确定关键点的位置和尺度，并去除其中低对比度和不稳定的点，以增强稳定性和抗噪能力。

接着，对每一个关键点，计算其金字塔图像中3σ邻域内像素的梯度以及方向分布，统计方向分布后为当前点分配一个或多个主方向。

最后，在关键点的4×4邻域内计算其8个方向的梯度信息，获得当前关键点的SIFT描述子。

通过特征点的位置和其像素值，判断当前特征点是否属于边界点。如果不是边界点，则使用距离该特征点最近的边界点替代该特征点。

进一步地，所述步骤S3具体包括：

S31、对于步骤S2中获得的每个边界特征点，比较一定范围内普通边界点与当前特征点的欧氏距离，并比较该距离与给定距离阈值R的关系，根据大小关系确定该边界点是否为支撑边界点，遍历该范围内所有边界点后获得当前特征点的支撑边界；若获得的欧氏距离小于设定的阈值R，则将该普通边界点视作支撑边界上的点。否则，不作任何处理。

S32、遍历所有的边界特征点，并将其各自的支撑边界像素点在空白图中标出，待所有支撑点都标出后即获得仅包含支撑边界的二值图像；

S33、对上述步骤中获得的二值图像进行连通性分析，将所有支撑边界分为彼此不连通的多个边界分段。

进一步地，所述步骤S4具体包括：

S41、根据每一条支撑边界上所有像素点的位置确定当前支撑边界的最小外接矩形，使用该外接矩形获取原肺分割图像中对应局部范围内的小块图像；

S42、使用二值图像边界检测算法检测出该小块图像对应的封闭边界；

S43、对封闭边界进行方向数计算得到其初始的Freeman边界链码，并基于公式计算封闭边界的综合码。在得到综合码后，通过对相邻的两个综合码求差得到对应的差值码；

S44、针对得到的差值码，基于一定的噪声筛选规则移除其中的噪声点；

S45、对去噪后的差值码重新编码得到R码，分别对R码前后相邻的两个码值分别相乘得到连续码，对连续码取最小值后得到弯曲度码；

S46、比较弯曲度码中最小码值和给定阈值的关系，从而判断当前外界矩形块是否需要被修复。

对获得的每一条支撑边界计算其最小外接矩形，其获取方式如下：

遍历每一条支撑边界，统计支撑边界上的所有像素点位置的水平坐标、垂直坐标的最大、最小值。获得的两对坐标值确定了四个角点，这四个点所圈定的矩形范围即为当前支撑边界最小外接矩形的范围。

然后同样使用上述四个角点，提取肺二值图像中上述相同矩形区域中的像素点，获得了二值图像中的局部小块图像。统计其中的目标物体的个数，若其中的目标物体个数大于1，则将面积最大的保留，其余目标的像素点的像素值转化为背景区域像素值。

使用基于弯曲度的编码对待处理的边界进行识别，其获取方式如下：

A.对当前的小块图像进行尺寸扩展，在其上下左右四个方向均扩展出一定像素宽度的背景区域。

B.对扩展后的小块图像使用Freeman方法获得其8方向边界链码，链码开始位置为目标的左上角点。

C.通过公式计算当前区域块的综合码cmb(i)。

D.基于计算所得的综合码cmb(i)，通过公式计算当前局部区域块的差值码{d′(i)|i＝1，2，…，N′}。并基于给定规则，在整个差值码编码队列中移除其中存在的噪声点。

E.根据去噪后差值码序列{d″(i)|i＝1，2，…，N″}中相邻两项的符号性质，将差值码序列转换为新的连续码序列R(j)。

F.基于上述获得的连续码序列，将其前后两个值分别相乘后求得其弯曲度码M(i)。

G.为当前图像设定阈值T，使用T与上述M(i)中的值作比较。如果其最小值小于阈值T，则认定当前小块图像区域中存在凹陷边界并需要后续的区域修复。否则，对当前图像块任何处理。

进一步地，所述步骤S5具体包括：

S51、对于需要修复的矩形块进行尺寸扩张，在小块图像的上下左右四个方向扩张出一定宽度的白色背景；

S51、将未扩张前小块图像的四条边框线设定为初始曲线轮廓，基于移除气球力和终端能量的梯度矢量流-蛇模型(GVF-Snake的方法)确定演化方程；

S53、进行迭代计算，其中根据演化方程更新坐标位置以驱动光滑曲线运动；

S54、待轮廓曲线演化结束，提取轮廓曲线包围的区域作为修复后的区域；

S55、使用上述修复区域替换原肺分割图像中对应的区域，待所有检测出的区域被处理后，肺二值图像即被修复。

对识别为凹陷的区域使用基于GVF-Snake的方法驱动光滑轮廓曲线以修复对该区域进行修复，其获取方式如下：

A.对上述步骤中尺寸扩展后的局部图像，计算图像内部目标(肺区域)的最小外接矩形，并顺时针提取该最小外接矩形的上下左右边框上的像素点，形成初始轮廓曲线K的坐标向量x，y。

B.基于上述尺寸扩展后的局部图像，计算图像中的外部力综合能量，该综合能量为线性能量与边缘能量的综合。

C.基于公式计算外部力综合能量的梯度。

D.基于公式计算图像中的梯度向量流(Gradient Vector Flow，GVF)。

E.迭代更新光滑轮廓曲线上像素点的位置。

F.使用上述得到的最终轮廓点{(x,y)}，设定轮廓线所圈定的内部区域为修复后的局部区域。

G.移除区域中尺寸拓展出的背景像素点，使用剩余像素点替换原肺二值图像中的像素点，从而完成了对单张图像上肺部区域的修复。

与现有技术相比，本发明的优点在于：

本发明方法通过计算经阈值处理后肺二值图像中存在的SIFT特征点，并基于特征点计算其替代点从而形成边界特征点集合。针对每一个边界特征点，使用半径阈值确定各自对应的支撑边界。获得所有支撑边界后，使用图像连通性分析方法获得图像中多条联通支撑边界。然后，针对每一条支撑边界取其最小外接矩形，并计算其顺时针方向的原始Freeman链码。根据原始链码，先后计算其均值码、综合码、差值码、连续码和弯曲度码。并根据弯曲度码的最小值和给定阈值T的关系决定当前局部图像块是否需要进行修复。对需要修复的局部区域块，先尺寸扩展。然后，基于最小外接矩形区域获得其初始轮廓上各像素点的坐标，并将它们作为初始轮廓曲线K的初始坐标向量。然后，基于图像中的信息，先后计算其图像外部力、梯度向量流。最后根据不带气球力和终端能量的迭代演化公式，更新轮廓上各像素点的坐标，使曲线以光滑的方式向理想的目标位置移动。最终，使用修复后的局部区域块替换原图像中的相关区域即可得到对单张CT图像进行凹陷区域修复后的结果。

应用本发明方法可以准确识别肺图像边界上有问题的部分，从而将肺边界修复工作集中于识别出的问题区域上，从而使本专利减少了处理目标数量并减少了肺区域修复步骤中的计算量，大大提高专利的计算效率。

应用本发明方法可以精确修复阈值分割结果中遗漏的肺膜黏连结节区域，将肺膜结节区域重新纳入肺实质区域中，从而保障了肺部CT图像病变组织分析算法和系统的准确性。同时，本发明方法克服了通过直线连接获取边界方式带来的修复后曲线太平直、不平滑的问题，更加符合肺部边界的形态学特性。最后，本专利方法计算的边界更加准确，确定的肺膜结节特征更精确，从而确保了后续步骤中对肺结节识别结果的准确性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明一种针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法的基本流程图；

图2为本发明实施例的处理对象：a是原始CT图像(左侧)的结构示意图；b是肺CT的阈值二值化图像(右侧)的结构示意图；

图3为本发明实施例使用SIFT算法检测二值图像中的特征点的结构示意图；

图4为本发明实施例计算支撑边界的算法原理图示；

图5为本发明实施例获取单条支撑边界最小外接矩形区域块的示意图；

图6为本发明实施例使用GVG-Snake算法修复局部区域块的算法原理图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，结合以下实施例，对本发明进一步详细说明。应当理解，此处所描述的具体实施案例仅用以解释本发明，并不用于限定本发明。

参阅图1所示，本发明提供一种针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，包括以下部分：

本算法是这样实现的，一种针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，所述对肺膜粘连结节区域修复方法包括以下步骤：

S1、参阅图2所示，获取肺部区域二值图像：输入经过阈值方法从CT图像中直接提取到的肺部区域二值图像；

S2、对经阈值法处理后的二值图像进行SIFT特征检测以获取图像中的特征点，并对上述特征点中的非边界点使用其最近邻边界点替代，使用替代点替换所有原非边界特征点后，即获得该图的边界特征点集合；

S3、对边界特征点集合中的每一个边界特征点，使用基于距离的规则提取肺部区域上的支撑边界；并根据计算得到的支撑边界图像中各支撑边界的联通性对其分段，获得图中彼此独立的支撑边界；

S4、对上述每一段支撑边界计算其最小外接矩形并基于综合码、差值码、连续码、弯曲度编码对待修复的边界进行识别；

S5、对待修复的边界进行修复，优选的，对识别为凹陷区域的部分使用基于GVF-Snake的方法驱动光滑轮廓曲线以修复修复凹陷区域以及相关边界。

上述步骤中每一步具体的实施方案为：

1.参阅图3所示，SIFT特征检测以及替代点计算

1.1步骤S2中对二值图像中的SIFT特征点是尺度空间中满足位置、尺度和旋转不变量的极值点，其获取方法如下：

首先，尺寸空间的建立以及相关极值检测：基于二维高斯函数以及高斯金字塔建立尺度空间，使用高斯差分算子检测尺度空间中的关键点；

然后，使用拟合三维二次方程精确确定关键点的位置和尺度，并去除其中低对比度和不稳定的点，以增强稳定性和抗噪能力；

接着，对每一个关键点，计算其金字塔图像中3σ邻域内像素的梯度以及方向分布，统计方向分布后为当前点分配一个或多个主方向；

最后，在关键点的4×4邻域内计算其8个方向的梯度信息，获得当前关键点的SIFT描述子。

1.2步骤S2对特征点中的非边界点使用其最近邻边界点替代，获得边界特征点集合，其获取方式如下：

通过特征点的位置和其在肺二值图像中的像素值，判断特征点是否属于边界点。如果不是边界点，则使用该特征点的最近邻边界点替代该特征点。判定方法如下：

使用二值图像检测算法检测肺二值图像中的所有边界点，所有边界点的集合为B。集合B＝{b_i(x，y)|i＝1，2，…，N′}，其中b_i(x，y)为肺二值图像上的边界点。如果s_i∈B，则不需执行任何操作。其中s_i(x，y)为检测出的肺部特征点。

如果

则执行下面操作。令S＝{s_i(x，y)|i＝1，2，…，N}为上述步骤检测出的特征点的集合，则

规则1：特征点集合中的s_i将被下面满足下面公式的b_i′代替。

遍历特征点集合S中的每一个特征点并按照规则1处理后，即可得到新的边界特征点集合S′。

2.支撑边界计算

2.1如图4所示，步骤S3中对每一个边界特征点，使用基于欧氏距离的规则提取肺部区域上的支撑边界，其获取方式如下：

计算一定范围内的普通边界点与当前特征点的欧氏距离。并比较该距离与给定距离阈值R的关系，根据大小关系确定该边界点是否为支撑边界点。若获得的欧氏距离小于设定的阈值，则将该普通边界点界点设定为支撑边界上的点。否则，不作任何处理。其具体实施方式可表达为：

令b_i(x，y)为集合B中的某个边界点，集合B′为由图像中的支撑边界像素点组成的集合，对于

若

其中，x_bi，y_bi，x_si，以及y_si，分别为点b_i和点s_i的水平、垂直坐标。

根据上述公式(2)，图像边界上所有像素点被分成了两类：支撑边界上的点以及非支撑边界上的点。

遍历所有的边界特征点，并将其各自的支撑边界像素点在空白图中标出后，获得仅包含支撑边界的二值图像。

2.2步骤S3中根据支撑边界图像中各支撑边界的联通性对其分段，获得图中彼此独立的支撑边界，其计算过程如下：

对上述得到的仅包含支撑边界的二值图像进行连通性分析后，根据图像中的支撑边界之间彼此的连通情况，图中所有的支撑边界被分成了几条彼此独立的分段。

3.凹陷区域识别

3.1步骤S4中对获得的每一条支撑边界计算其最小外接矩形，其获取方式如下：

如图5所示，遍历每一条支撑边界，统计每条支撑边界上的所有像素点位置的水平坐标、垂直坐标的最大、最小值。

最大值、最小值两对坐标确定了四个角点，这四个角点所圈定的矩形即为当前支撑边界最小外接矩形。

同样使用上述四个角点确定的矩形区域，提取肺二值图像中对应区域的像素点，因此获得了二值图像中的局部小块图像。统计其中的目标物体的个数，若其中的目标物体个数大于1，则将面积最大的保留，其余目标的像素点的像素值转化为背景区域像素值。其具体实施表达如下：

令x_min＝min{x|b_i(x，y)∈B_i′}，

x_max＝max{x|b_i(x，y)∈B_i′}，

y_min＝min{y|b_i(x，y)∈B_i′}，

y_max＝max{y|b_i(x，y)∈B_i′}，

其中B_i′为上述某一条支撑边界上的像素点集合。

则该支撑边界对应的最小外接矩形区域的范围为：[x_min→x_max]×[y_min→y_max]。取原肺二值图像中相同范围内的像素点，获得局部图像块。

步骤S4中使用基于弯曲度的编码对待处理的边界进行识别，其获取方式如下：

3.2对上述步骤获得的局部图像块进行尺寸扩展，在其上下左右四个方向均扩展出一定像素宽度的背景区域。于是，扩展后的新区域四周存在一定像素宽度的背景区域，而中间部分为待处理的肺局部区域块。

3.3对扩展后的局部区域块使用Freeman方法获得其8方向边界原始链码c(i)，i＝1，2，…，P，该原始链码开始位置为目标的左上角点。

3.4.计算当前局部区域块的综合码，其计算方式如下：

其中，i＝2，3，…，P-3。

由上述公式可知，该综合码反映了当前边界点及其邻域内的平均方向。

3.5基于上述结果cmb(i)，计算当前局部区域块的差值码，其计算方式如下：

d(i)＝cmb(i-1)-cmb(i)，i＝2，3，4，…，P (4)

d(1)＝cmb(P)-cmb (1)

由上述公式可知，该差值码反映了当前边界点及其邻域内方向数的平均变化趋势。

3.6通过公式(4)计算所得的差值码中存在大量的噪声点，这些噪声点是检测曲线的凸凹性过程中的干扰点。对算法而言，这些干扰点要从整个编码队列中移除，其对应的移除规则可表达如下：

规则2

d(i-1)	d(i)	d(i+1)	d(i+2)	d(i+3)	是否噪声点？
						＝0	＝0	≠0			不是
≠0	＝0	≠0			是
						＝0	＝0	＝0	≠0		不是
≠0	＝0	＝0	≠0		是
						＝0	＝0	＝0	＝0	≠0	不是
≠0	＝0	＝0	＝0	≠0	是

根据上述规则2，可判定差值码中某一个点所对应的编码d(i)是否为噪声点。

如果d(i)被判定为噪声点，则将d(i)从整个编码队列中移除。如果d(i)被判定为非噪声点，则不做任何操作。

将所有d(i)点按照规则2处理完毕后，原差值码序列{d(i)}变成了新的差值码序列{d′(i)}。

3.7根据去噪声后差值码序列{d′(i)|i＝1，2，…，N′}中相邻两项的符号性质，将差值码序列转换为新的连续码序列R(j)，该转换过程的基本公式为：

对于i＝1，2，…，K-1，p＝1，2，…，N′-2，如果

则连续码{R(j)|j＝1，2，…，N′}中的R(p)＝K，R(p+1)＝K，R(p+2)＝K，…，R(p+k-1)＝K。/>

对于i＝1，2，…，K-1，p＝1，2，…，N′-2，如果

则连续码{R(j)|j＝1，2，…，N′}中的R(p)＝-K，R(p+1)＝-K，R(p+2)＝-K，…，R(p+k-1)＝-K。

如果d′(p)＝0，那么R(p)＝0。 (5)

如果p∈{N′-2，N′-1，N′}，考虑数量较少，其干扰效果不大，故对这三个点不作处理。

3.8基于上述获得的连续码序列R(j)，将其前后连续的两个值相乘后求得其弯曲度码M(i)，其计算的基本方法为：

M₁(i)＝R(i)R(i-1) (6)

M₁(1)＝R(1)R(N′)

M₂(i)＝R(i)R(i+1) (7)

M₂(N′)＝R(N′)R(1)

M(i)＝min{M₁(i)，M₂(i)}，i＝1，2，…，N′ (8)

3.9为当前图像设定阈值T，使用T与上述M(i)中最小值作比较。如果其最小值小于T，则认定当前局部区域块存在凹陷边界并需要后续的区域修复。否则，不予任何处理。其计算公式为

其中，对由上述公式(9)识别出需要修复的局部区域块使用下面的方法进行处理。

4.如图6所示，步骤S5中对识别为凹陷的区域使用基于GVF-Snake的方法驱动光滑轮廓曲线以修复对该区域进行修复，其获取方式如下：

4.1在尺寸扩展后的局部区域块中，将当前支撑边界分段的最小外接矩形作为初始轮廓曲线K。

按照顺时针方向提取该最小外接矩形的上下左右边框上的像素点，形成初始轮廓曲线K的水平和垂直坐标向量x，y。

4.2基于尺寸扩展后的局部区域块，计算图像对应的的外部力综合能量。该综合能量为线性能量与边缘能量的线性组合，其计算方式如下：

E_line＝I(x，y)*G(x，y，ε)，其中E_line为线性能量，I(x，y)为尺寸扩展后的局部图像，G(x，y，ε)为标准差为ε的二维高斯函数。

(10)

其中E_edge为边缘能量，/>

为图像I(x，y)的梯度。

(11)

E_ext＝E_line+E_edge，其中E_ext为外部力综合能量。

4.3计算图像外部力综合能量的梯度，其计算公式如下：

/>

4.4计算图像中的梯度向量流(Gradient Vector Flow，GVF)，令u＝F_x，v＝F_y，通过如下的公式进行计算得到

其中，u_t-1，v_t-1分别为第t-1次迭代中u，v的值，u_t，v_t分别为第t次迭代中u，v的值。

为Laplacian算子，该项能够保障u_t和v_t的平滑性并能阻止其曲线扭曲或者发散。

4.5基于迭代算法更新轮廓上像素点的位置，其迭代更新公式如下：

其中，I为单位矩阵，α为轮廓能量对应的权重，γ为时间步长，x_t，y_t分别为第t次迭代中轮廓曲线上像素点的水平和垂直坐标值。

4.6上述迭代公式稳定后，可得到的最终轮廓点集合{(x，y)}。另外，这些轮廓点确定的曲线所圈定的内部区域即为修复后的局部区域L。

4.7对当前处理的局部区域块移除拓展区域，将局部区块中剩余像素点替换原肺二值图像中的像素点，即完成了对单张图像上肺部区域的修复。

本发明方法通过计算经阈值处理后肺二值图像中存在的SIFT特征点，并基于特征点计算其替代点从而形成边界特征点集合。针对每一个边界特征点，使用半径阈值确定各自对应的支撑边界。获得所有支撑边界后，使用图像连通性分析方法获得图像中多条联通支撑边界。然后，针对每一条支撑边界取其最小外接矩形，并计算其顺时针方向的原始Freeman链码。根据原始链码，先后计算其均值码、综合码、差值码、连续码和弯曲度码。并根据弯曲度码的最小值和给定阈值T的关系决定当前局部图像块是否需要进行修复。对需要修复的局部区域块，先尺寸扩展。然后，基于最小外接矩形区域获得其初始轮廓上各像素点的坐标，并将它们作为初始轮廓曲线K的初始坐标向量。然后，基于图像中的信息，先后计算其图像外部力、梯度向量流。最后根据不带气球力和终端能量的迭代演化公式，更新轮廓上各像素点的坐标，使曲线以光滑的方式向理想的目标位置移动。最终，使用修复后的局部区域块替换原图像中的相关区域即可得到对单张CT图像进行凹陷区域修复后的结果。

应用本发明方法可以准确识别肺图像边界上有问题的部分，从而将肺边界修复工作集中于识别出的问题区域上，从而使本专利减少了处理目标数量并减少了肺区域修复步骤中的计算量，大大提高专利的计算效率。

应用本发明方法可以精确修复阈值分割结果中遗漏的肺膜黏连结节区域，将肺膜结节区域重新纳入肺实质区域中，从而保障了肺部CT图像病变组织分析算法和系统的准确性。同时，本发明方法克服了通过直线连接获取边界方式带来的修复后曲线太平直、不平滑的问题，更加符合肺部边界的形态学特性。最后，本专利方法计算的边界更加准确，确定的肺膜结节特征更精确，从而确保了后续步骤中对肺结节识别结果的准确性。

本发明的实现并不局限于上述实例所公开的范围，可以采用不同于上述实例的方式实现上述技术方案。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书内容所做的等效结构或等效流程变换，或直接或间接运用在其他相关技术领域，均同理包括在本发明的专利保护范围内。

Claims (6)

1.一种针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，其特征在于，包括以下步骤：

S1、获取肺部区域二值图像；

S2、对上述二值图像进行SIFT特征检测获得图像中的特征点，并对上述获得的非边界特征点计算其最近邻替代点，并使用计算出的替代点替换原非边界特征点后，获得对应的边界特征点集合；

S3、针对边界特征点集合中的每一个边界特征点，基于特定的规则提取支撑边界并根据支撑边界；并根据计算得到的支撑边界图像中各支撑边界的联通性对其分段，获得图中彼此独立的支撑边界；

S4、对上述每一段支撑边界计算其最小外接矩形并基于综合码、差值码、连续码、弯曲度编码对待修复的边界进行识别；

S5、对待修复的边界进行修复；

所述步骤S3具体包括：

S31、对于步骤S2中获得的每个边界特征点，比较一定范围内普通边界点与当前特征点的欧氏距离，并比较该距离与给定距离阈值R的关系，根据大小关系确定该边界点是否为支撑边界点，遍历该范围内所有边界点后获得当前特征点的支撑边界；

S32、遍历所有的边界特征点，并将其各自的支撑边界像素点在空白图中标出，待所有支撑点都标出后即获得仅包含支撑边界的二值图像；

S33、对上述步骤中获得的二值图像进行连通性分析，将所有支撑边界分为彼此不连通的多个边界分段。

2.根据权利要求1所述的针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，其特征在于：所述步骤S5中，对识别为凹陷区域的部分使用基于GVF-Snake的方法驱动光滑轮廓曲线以修复凹陷区域以及相关边界。

3.根据权利要求2所述的针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，其特征在于：所述步骤S1中，输入经过阈值方法从CT图像中直接提取到的肺部区域二值图像。

4.根据权利要求3所述的针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，其特征在于：步骤S2具体包括：

S21、首先基于特征点检测、特征点描述、特征点匹配和消除错配点的SIFT特征点检测步骤后，得到二值图像中的特征点位置；

S22、对上述检测到的每一个SIFT特征点，若其不在肺部边界上则寻找距离该特征点最近的边界点作为其替代点，从而将特征点集合转变为边界特征点集合。

5.根据权利要求1所述的针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，其特征在于，所述步骤S4具体包括：

S41、根据每一条支撑边界上所有像素点的位置确定当前支撑边界的最小外接矩形，使用该外接矩形获取原肺分割图像中对应局部范围内的小块图像；

S42、使用二值图像边界检测算法检测出该小块图像对应的封闭边界；

S43、对封闭边界进行方向数计算得到其初始的Freeman边界链码，并基于公式计算封闭边界的综合码，在得到综合码后，通过对相邻的两个综合码求差得到对应的差值码；

S44、针对得到的差值码，基于一定的噪声筛选规则移除其中的噪声点；

S45、对去噪后的差值码重新编码得到R码，分别对R码前后相邻的两个码值分别相乘得到连续码，对连续码取最小值后得到弯曲度码；

S46、比较弯曲度码中最小码值和给定阈值的关系，从而判断当前外界矩形块是否需要被修复；

计算当前局部区域块的综合码，其计算方式如下：

其中，i＝2，3，…，P-3。

6.根据权利要求5所述的针对肺CT图像阈值分割结果的肺膜粘连结节区域修复方法，其特征在于，所述步骤S5具体包括：

S51、对于需要修复的矩形块进行尺寸扩张，在小块图像的上下左右四个方向扩张出一定宽度的白色背景；

S51、将未扩张前小块图像的四条边框线设定为初始曲线轮廓，基于移除气球力和终端能量的梯度矢量流-蛇模型确定演化方程；

S53、进行迭代计算，其中根据演化方程更新坐标位置以驱动光滑曲线运动；

S54、待轮廓曲线演化结束，提取轮廓曲线包围的区域作为修复后的区域；

S55、使用上述修复区域替换原肺分割图像中对应的区域后，肺二值图像即被修复。

Images