CN106097305A - 双行程区域生长结合形态学重建的肺部气管树分割方法 - Google Patents

Abstract

本发明是为了解决上述问题而进行的，目的在于提供一种基于双行程三维区域生长与形态学相结合的肺部气管树分割算法，以实现对肺部气管树进行有效的分割。本发明提供的双行程区域生长结合形态学重建的肺部气管树分割方法，包括如下步骤：自动从序列断层图像中获取三维区域生长的种子点P₁；以P₁作为种子点，进行第一行程区域生长，得到肺部主支气管树；设定多尺度形态学算子，并利用灰度形态学重建操作，得到三个横切面、冠状面和矢状面解剖位上的对应的气管标记图；以肺部主支气管树为基础、三个解剖位上的气管标记图为限定条件，进行第二行程区域生长，得到对应的分割气管树；将分割气管树进行融合，得到最终的肺部气管树分割结果F。

Description

双行程区域生长结合形态学重建的肺部气管树分割方法

技术领域

本发明涉及一种医学图像处理领域，具体涉及一种双行程区域生长结合形态学重建的肺部气管树分割方法。

背景技术

肺部气管树的解剖结构，是诊断支气管狭窄、肺气肿以及慢阻肺等多种肺部疾病的病理评判标准，因此能否从CT图像中精确地分割出肺部气管树，最终会影响到临床诊断的准确性。肺部气管是人体与外界进行气体交换的重要组织，从咽喉开始由主支气管依次延伸，形成形态复杂、走向多变的类树状结构。由于CT图像中的部分容积效应、CT图像在拍摄过程中因病人运动而产生的运动伪影以及现行临床CT机极限分辨率等的限制，产生外围高级气管壁模糊、断连的现象，造成肺部气管分割过程中易逾越气管范围泄漏至肺实质中，增加了肺部气管分割的难度。针对这些问题，目前国内外只有较少量文献报道。

Lo等人在2010年Medical Image Analysis期刊上发表的论文《Vessel-guidedairway tree segmentation:A voxel classification approach》提出的基于肺部血管引导的气管分割方法，该方法利用血管与气管在解剖结构上并行的特点，以血管作为分割引导，能够较为精准地分割出肺部气管树，但算法需确定血管走向，复杂度高，时间过长；KaiLAI等人在2009年International Conference on Bioinformatics and BiomedicalEngineering上发表的论文《Automatic 3D Segmentation of Lung Airway Tree:A NovelAdaptive Region Growing Approach》，提出基于感兴趣区域划分的分割方法，虽然利用先划分感兴趣区域再进行分割的方法，可以有效地防止分割泄漏至肺实质，同时提高算法效率，但实验结果表明探测出的高级外围气管较少，分割效果欠佳；Michael W等人在2010年IEEE Transactions on Medical Imaging期刊上发表的论文《Robust 3-D airway treesegmentation for image-guided peripheral bronchoscopy》，将区域生长与气管几何结构分析相结合，以图像局部特征为引导对肺部气管树进行分割，该方法虽然能探测出较多的高级外围气管且准确性较好，但易受到低剂量CT图像中噪声和伪影等因素的影响，且需手动设置的参数过多；Anna等人在2009年Computerized Medical Imaging and Graphics期刊上发表的论文《Two-pass region growing algorithm for segmenting airway treefrom MDCT chest scans》利用双行程区域生长方法，通过形态学腐蚀和膨胀操作探测细小气管，但对高级外围气管的探测效果不佳，在低剂量、层厚较厚的CT图像上表现较差，以气管各级分支为定量指标，最多分割出第十级气管，气管总分支数最多为145，分割效果有待提高。

发明内容

本发明是为了解决上述问题而进行的，目的在于提供一种基于双行程三维区域生长与形态学相结合的肺部气管树分割算法，以实现对肺部气管树进行有效的分割。

本发明提供了一种双行程区域生长结合形态学重建的肺部气管树分割方法，具有这样的特征，包括如下步骤：步骤一，输入待分割的胸部CT的序列断层图像；步骤二，根据肺部气管在序列断层图像中对应的位置及其CT值特点，自动从序列断层图像中获取三维区域生长的种子点P₁；步骤三，以P₁作为种子点，采用3D区域生长算法进行第一行程区域生长，得到肺部主支气管树A；步骤四，设定多尺度形态学算子，并利用灰度形态学重建操作，分别从横切面、冠状面和矢状面三个方向进行潜在气管区域检测，得到三个解剖位上的对应的气管标记图M₁，M₂和M₃；步骤五，以肺部主支气管树A为基础、三个解剖位上的气管标记图M₁，M₂和M₃为限定条件，利用3D区域生长算法进行第二行程区域生长，得到对应的分割气管树A₁，A₂和A₃；步骤六，依据空间位置关系，将分割气管树A₁，A₂和A₃进行融合，得到最终的肺部气管树分割结果F。

在本发明提供的双行程区域生长结合形态学重建的肺部气管树分割方法中，还可以具有这样的特征，其中步骤三中还包括以下子步骤：子步骤3-1，设定波动范围Index的值；子步骤3-2，从种子点P₁开始，以堆栈队列型的方式存储种子点，当当前像素点是种子点的26邻域，且该像素点的CT值以及6邻域的CT值分别与被纳入气管区域像素点平均CT值的差值百分比ΔV不超过波动范围Index时，被纳入气管区域：

$\mathrm{\Δ}V=\left|\frac{{\mathrm{value}}_c-{\mathrm{value}}_a}{{\mathrm{value}}_a}\right|$

其中，ΔV表示当前像素点CT值与气管区域平均CT值之间的差值百分比，value_c表示当前像素点CT值，value_a表示气管区域平均CT值。

在本发明提供的双行程区域生长结合形态学重建的肺部气管树分割方法中，还可以具有这样的特征，其中子步骤3-1为：以Index＝0.05为初始值，ΔIndex＝0.02为步长，采用迭代寻优的方式选取波动范围Index的值，在每次生长迭代过程中递增Index，当已分割出的气管区域平均CT值小于等于-850HU、Index的值不大于0.11、且差值百分比ΔN不大于20％时，判定气管未发生泄漏，且尚未分割完全，迭代继续；否则，判定生长发生了泄漏，此时将当前的波动范围Index的值递减0.02，并将递减后的Index的值设定为最优波动范围Index_opt的值：

$\mathrm{\Δ}N=\frac{N_c-N_f}{N_f}\≤0.2,$

其中，ΔN为像素数N_c和像素数N_f的差值百分比，像素数N_c为当前区域生长结果的像素数，像素数N_f为上次迭代中区域生长结果的像素数。

在本发明提供的双行程区域生长结合形态学重建的肺部气管树分割方法中，还可以具有这样的特征，其中，步骤四中，多尺度形态学算子为：依据肺部气管在序列断层图像上的形状特征，定义一个4连通域的单位结构元素D，单位结构元素D对本身进行n次膨胀，从而得到不同尺寸的多尺度结构元素D_n，计算公式如下式所示：

式中，多尺度结构元素D_n表示用单位结构元素D对本身进行n次膨胀，n的最小值为1，最大值为6。

在本发明提供的双行程区域生长结合形态学重建的肺部气管树分割方法中，还可以具有这样的特征，其中步骤四中，灰度形态学重建操作包括以下子步骤：子步骤4-1，使用多尺度结构元素D_n，对单层序列断层图像进行形态学闭操作处理得到标记图像X_i，计算公式为：

式中B_i表示胸部CT的序列断层图像的原图像，i代表胸部CT的序列断层图像中的任意CT切层，D表示单位结构元素，X_i表示B_i由单位结构元素D经一次膨胀与腐蚀后得到的标记图像；子步骤4-2，运用标记图像X_i与原图B_i进行重构操作，重复运算直至I_k+1等于I_k，从而得到对应的高级气管重构图I_k+1，之后利用原图B_i与重构图像I_k+1相减得到标记图Map_Gray_i，I_k+1的计算公式为：

$I_{k+1}=(I_k\⊕D)\∩B_i,$

其中I₁＝X_i；子步骤4-3，选取合适的阈值，对标记图Map_Gray_i进行阈值化处理，采用的阈值化处理公式为：

Threshold＝T_index×(MaxValue-MinValue)+MinValue，

式中Threshold代表最终的阈值，Maxvalue代表标记图Map_Gray_i中的最大值，Minvalue代表标记图Map_Gray_i中的最小值，T_index代表Maxvalue与Minvalue的差值系数，取值为0.2。

在本发明提供的双行程区域生长结合形态学重建的肺部气管树分割方法中，还可以具有这样的特征：其中，在步骤五中，在第二行程区域生长中当像素点与种子点P₁为26邻域且位于M_i(i＝1，2，3)中，同时像素点的CT值与被纳入气管区域像素点平均CT值的差值百分比不大于2倍的Index时，被纳入肺部气管树区域。

在本发明提供的双行程区域生长结合形态学重建的肺部气管树分割方法中，还可以具有这样的特征：其中，Index值选取第一行程区域生长过程中得到的最优值Index_opt。

发明的作用与效果

根据本发明所涉及的双行程区域生长结合形态学重建的肺部气管树分割方法，因为设定多尺度形态学算子，并利用灰度形态学重建操作，分别从横切面、冠状面和矢状面三个方向进行潜在气管区域检测，得到三个解剖位上的对应的气管标记图，并以肺部主支气管树为基础、三个解剖位上的气管标记图为限定条件，利用3D区域生长算法进行第二行程区域生长，所以，本发明的双行程区域生长结合形态学重建的肺部气管树分割方法具有简单易实现，其分割结果显示能够分割出更高级别的气管树同时气管树分支数大幅度提高，而且分割结果侵入肺实质的比例也较低。

附图说明

图1是本发明的实施例中双行程区域生长结合形态学重建的肺部气管树分割方法的方法流程图；

图2是本发明的实施例中双行程区域生长结合形态学重建的肺部气管树分割方法的分割路径图；

图3是本发明的实施例中双行程区域生长结合形态学重建的肺部气管树分割方法的针对不同案例获得的分割结果图；以及

图4是本发明的实施例中双行程区域生长结合形态学重建的肺部气管树分割方法与TWORG的分割结果统计对照表，其中TWORG代表文献《Two-pass region growingalgorithm for segmenting airway tree from MDCT chest scans》中的算法。

具体实施方式

为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，以下实施例结合附图对本发明双行程区域生长结合形态学重建的肺部气管树分割方法作具体阐述。

实施例

图1是本发明的实施例中双行程区域生长结合形态学重建的肺部气管树分割方法的方法流程图；图2是本发明的实施例中双行程区域生长结合形态学重建的肺部气管树分割方法的分割路径图。

如图1所示，本发明的实施例中双行程区域生长结合形态学重建的肺部气管树分割方法包括如下步骤：

步骤一，输入待分割DICOM格式的胸部CT的序列断层图像。

步骤二，通过自动选取，在序列断层图像中获取三维区域生长种子点P1。肺气管内部的CT值范围为-1024HU～-800HU，且主支气管位于图像中间区域，直径为1.0cm～2.3cm。因此，设阈值T＝-800HU，得到二值图像。在二值图像中间区域搜索位置与面积符合要求的连通域并求其质心，以此质心作为最终的三维区域生长的种子点P₁。

步骤三，以P₁作为种子点，采用3D区域生长算法进行第一行程区域生长，得到肺部主支气管树A：

子步骤3-1，设定波动范围Index的值：以Index＝0.05为初始值，ΔIndex＝0.02为步长，采用迭代寻优的方式选取波动范围Index的值，在每次生长迭代过程中递增Index，当已分割出的气管区域平均CT值小于等于-850HU、Index的值不大于0.11、且差值百分比ΔN不大于20％时，判定气管未发生泄漏，且尚未分割完全，迭代继续；否则，判定生长发生了泄漏，此时将当前的波动范围Index的值递减0.02，并将递减后的Index的值设定为最优波动范围Index_opt的值：

$\mathrm{\Δ}N=\frac{N_c-N_f}{N_f}\≤0.2,$

其中，ΔN为像素数N_c和像素数N_f的差值百分比，像素数N_c为当前区域生长结果的像素数，像素数N_f为上次迭代中区域生长结果的像素数；

子步骤3-2，从种子点P₁开始，以堆栈队列型的方式存储种子点，当当前像素点是种子点的26邻域，且该像素点的CT值以及6邻域的CT值分别与被纳入气管区域像素点平均CT值的差值百分比ΔV不超过波动范围Index时，被纳入气管区域：

$\mathrm{\Δ}V=\left|\frac{{\mathrm{value}}_c-{\mathrm{value}}_a}{{\mathrm{value}}_a}\right|$

其中，ΔV表示当前像素点CT值与气管区域平均CT值之间的差值百分比，value_c表示当前像素点CT值，value_a表示气管区域平均CT值。

步骤四，设定多尺度形态学算子，并利用灰度形态学重建操作，分别从横切面、冠状面和矢状面三个方向进行潜在气管区域检测，得到三个解剖位上的对应的气管标记图M₁，M₂和M₃：

$M_i=\underset{i=1}{\overset{s}{\Σ}}{\mathrm{Map}}_i,$

M_i表示任一解剖位的最终标记图集合，Map_i表示任一解剖位的单层标记图，S表示CT序列图的总层数。

其中，多尺度形态学算子为：依据肺部气管在序列断层图像上的形状特征，定义一个4连通域的单位结构元素D，单位结构元素D对本身进行n次膨胀，从而得到不同尺寸的多尺度结构元素D_n，计算公式如下式所示：

式中，多尺度结构元素D_n表示用单位结构元素D对本身进行n次膨胀，n的最小值为1，最大值为6；

灰度形态学重建操作包括以下子步骤：

子步骤4-1，使用多尺度结构元素D_n，对单层序列断层图像进行形态学闭操作处理得到标记图像X_i，计算公式为：

式中B_i表示胸部CT的序列断层图像的原图像，i代表胸部CT的序列断层图像中的任意CT切层，D表示单位结构元素，X_i表示B_i由单位结构元素D经一次膨胀与腐蚀后得到的标记图像；

子步骤4-2，运用标记图像X_i与原图B_i进行重构操作，重复运算直至I_k+1等于I_k，从而得到对应的高级气管重构图I_k+1，之后利用原图B_i与重构图像I_k+1相减得到标记图Map_Gray_i，I_k+1的计算公式为：

$I_{k+1}=(I_k\⊕D)\∩B_i,$

其中I₁＝X_i；

子步骤4-3，选取合适的阈值，对标记图Map_Gray_i进行阈值化处理，采用的阈值化处理公式为：

Threshold＝T_index×(MaxValue-MinValue)+MinValue，

式中Threshold代表最终的阈值，Maxvalue代表标记图Map_Gray_i中的最大值，Minvalue代表标记图Map_Gray_i中的最小值，T_index

代表Maxvalue与Minvalue的差值系数，取值为0.2。

步骤五，以肺部主支气管树A为基础、三个解剖位上的气管标记图M₁，M₂和M₃为限定条件，利用3D区域生长算法进行第二行程区域生长，得到对应的分割气管树A₁，A₂和A₃。

为了能够在不发生泄漏的情况下，尽可能地分割出较高级别的细小气管，在粗分割结果的基础上，进行第二行程区域生长。当前像素点满足以下条件时，被纳入肺部气管树区域：

在第二行程区域生长中当像素点与种子点P₁为26邻域且位于M_i(i＝1，2，3)中，同时像素点的CT值与被纳入气管区域像素点平均CT值的差值百分比不大于2倍的Index时，被纳入肺部气管树区域。Index值选取第一行程区域生长过程中得到的最优值Index_opt。

由于引入了标记图集合M_i作为区域生长的路径，导致肺部气管树的平均CT值范围扩大，所以适当的增加波动范围是必要的。本文选择第一行程区域生长过程中得到的最优Index值Index_opt，并在此基础上将波动范围扩大为2倍的Index_opt。经过实验验证，以第一次区域生长确定的Index值为基础进行CT值范围的扩展，相比设定固定的扩展范围，能够更好的依据数据的实际情况提取细小气管，同时，也能防止因波动范围过大，而导致的二次分割泄漏。

步骤六，依据空间位置关系，将分割气管树A₁，A₂和A₃进行融合，得到最终的肺部气管树分割结果F。

本实验数据来自上海肺科医院40排螺旋CT扫描图像与LIDC数据库，成像层厚0.75mm～2.00mm，空间分辨率为0.6836mm～0.7993mm，每幅断层图像都是512×512的16位DICOM格式图像。总共测试了25组胸部CT序列断层图像(即25个病例)，该分割实验是在Matlab R2015a环境下进行的。

图3是本发明的实施例中双行程区域生长结合形态学重建的肺部气管树分割方法的针对不同案例获得的分割结果图。

如图3所示，对上述胸部CT序列断层图像，用本发明所述方法对肺部气管树进行了分割。

为了定量地评价该分割方法的性能，以各级气管分支数、气管总分支数为定量指标，对分割结果进行定量比较分析。

图4是本发明的实施例中双行程区域生长结合形态学重建的肺部气管树分割方法与TWORG的分割结果统计对照表，其中TWORG代表文献《Two-pass region growingalgorithm for segmenting airway tree from MDCT chest scans》中的算法。

如图4所示，利用本发明方法对25个肺部气管树进行了分割，并根据上述定量评判标准对每个分割案例的各级气管分支数以及气管总分支数进行了统计，从实验结果可得出，本发明较现有国内外文献，能够探测出更多的有效气管区域，最多可分割出12级气管，并且各级气管数目在不同的CT数据中都有较好的表现。能够探测出更多的气管树目，气管总分支数最高为293，最低为83，平均气管分支数为156。因此，使用本发明方法可以较有效地分割出肺部气管树。

实施例的作用与效果

根据本实施例所涉及的双行程区域生长结合形态学重建的肺部气管树分割方法，因为设定多尺度形态学算子，并利用灰度形态学重建操作，分别从横切面、冠状面和矢状面三个方向进行潜在气管区域检测，得到三个解剖位上的对应的气管标记图，并以肺部主支气管树为基础、三个解剖位上的气管标记图为限定条件，利用3D区域生长算法进行第二行程区域生长，所以，本实施例的双行程区域生长结合形态学重建的肺部气管树分割方法具有简单易实现，其分割结果显示能够分割出更高级别的气管树同时气管树分支数大幅度提高，而且分割结果侵入肺实质的比例也较低。而且本实施例所涉及的双行程区域生长结合形态学重建的肺部气管树分割方法能够探测出更多的有效气管区域，最多可分割出12级气管，并且各级气管数目在不同的CT数据中都有较好的表现。

上述实施方式为本发明的优选案例，并不用来限制本发明的保护范围。

Claims (7)

1.一种双行程区域生长结合形态学重建的肺部气管树分割方法，其特征在于，包括如下步骤：

步骤一，输入待分割的胸部CT的序列断层图像；

步骤二，根据肺部气管在所述序列断层图像中对应的位置及其CT值特点，自动从所述序列断层图像中获取三维区域生长的种子点P₁；

步骤三，以P₁作为种子点，采用3D区域生长算法进行第一行程区域生长，得到肺部主支气管树A；

步骤四，设定多尺度形态学算子，并利用灰度形态学重建操作，分别从横切面、冠状面和矢状面三个方向进行潜在气管区域检测，得到三个解剖位上的对应的气管标记图M₁，M₂和M₃；

步骤五，以所述肺部主支气管树A为基础、所述三个解剖位上的气管标记图M₁，M₂和M₃为限定条件，利用3D区域生长算法进行第二行程区域生长，得到对应的分割气管树A₁，A₂和A₃；

步骤六，依据空间位置关系，将所述分割气管树A₁，A₂和A₃进行融合，得到最终的肺部气管树分割结果F。

2.根据权利要求1所述的双行程区域生长结合形态学重建的肺部气管树分割算法，其特征在于，其中步骤三中还包括以下子步骤：

子步骤3-1，设定波动范围Index的值；

子步骤3-2，从种子点P₁开始，以堆栈队列型的方式存储种子点，当当前像素点是种子点的26邻域，且该像素点的CT值以及6邻域的CT值分别与被纳入气管区域像素点平均CT值的差值百分比ΔV不超过波动范围Index时，被纳入气管区域：

$\mathrm{\Δ}V=\left|\frac{{\mathrm{value}}_c-{\mathrm{value}}_a}{{\mathrm{value}}_a}\right|$

其中，ΔV表示当前像素点CT值与气管区域平均CT值之间的差值百分比，value_c表示当前像素点CT值，value_a表示气管区域平均CT值。

3.根据权利要求2所述的双行程区域生长结合形态学重建的肺部气管树分割算法，其特征在于，其中子步骤3-1为：

以Index＝0.05为初始值，ΔIndex＝0.02为步长，采用迭代寻优的方式选取波动范围Index的值，在每次生长迭代过程中递增Index，当已分割出的气管区域平均CT值小于等于-850HU、Index的值不大于0.11、且差值百分比ΔN不大于20％时，判定气管未发生泄漏，且尚未分割完全，迭代继续；否则，判定生长发生了泄漏，此时将当前的波动范围Index的值递减0.02，并将递减后的Index的值设定为最优波动范围Index_opt的值：

$\mathrm{\Δ}N=\frac{N_c-N_f}{N_f}\≤0.2,$

其中，ΔN为像素数N_c和像素数N_f的差值百分比，像素数N_c为当前区域生长结果的像素数，像素数N_f为上次迭代中区域生长结果的像素数。

4.根据权利要求1所述的双行程区域生长结合形态学重建的肺部气管树分割算法，其特征在于，其中步骤四中，所述多尺度形态学算子为：

依据肺部气管在所述序列断层图像上的形状特征，定义一个4连通域的单位结构元素D，单位结构元素D对本身进行n次膨胀，从而得到不同尺寸的多尺度结构元素D_n，计算公式如下式所示：

式中，多尺度结构元素D_n表示用单位结构元素D对本身进行n次膨胀，n的最小值为1，最大值为6。

5.根据权利要求1所述的双行程区域生长结合形态学重建的肺部气管树分割算法，其特征在于，其中步骤四中，所述灰度形态学重建操作包括以下子步骤：

子步骤4-1，使用多尺度结构元素D_n，对单层所述序列断层图像进行形态学闭操作处理得到标记图像X_i，计算公式为：

式中B_i表示所述胸部CT的序列断层图像的原图像，i代表所述胸部CT的序列断层图像中的任意CT切层，D表示单位结构元素，X_i表示B_i由单位结构元素D经一次膨胀与腐蚀后得到的标记图像；

子步骤4-2，运用标记图像X_i与原图B_i进行重构操作，重复运算直至I_k+1等于I_k，从而得到对应的高级气管重构图I_k+1，之后利用原图B_i与重构图像I_k+1相减得到标记图Map_Gray_i，I_k+1的计算公式为：

$I_{k+1}=(I_k\⊕D)\∩B_i,$

其中I₁＝X_i；

子步骤4-3，选取合适的阈值，对所述标记图Map_Gray_i进行阈值化处理，采用的阈值化处理公式为：

Threshold＝T_index×(Max Value-Min Value)+Min Value，

式中Threshold代表最终的阈值，Maxvalue代表所述标记图Map_Gray_i中的最大值，Minvalue代表所述标记图Map_Gray_i中的最小值，T_index代表Maxvalue与Minvalue的差值系数，取值为0.2。

6.根据权利要求2所述的双行程区域生长结合形态学重建的肺部气管树分割算法，其特征在于：

其中，在步骤五中，在第二行程区域生长中当所述像素点与种子点P₁为26邻域且位于M_i(i＝1，2，3)中，同时所述像素点的CT值与被纳入气管区域像素点平均CT值的差值百分比不大于2倍的Index时，被纳入肺部气管树区域。

7.根据权利要求6所述的双行程区域生长结合形态学重建的肺部气管树分割算法，其特征在于：

其中，Index值选取第一行程区域生长过程中得到的最优值Index_opt。