X线胸片是目前临床上广泛使用的医学成像技术，对于肺部疾病的分析与诊断具有重要的临床意义。对胸片中肺野的准确分割是对疾病的分析与诊断的重要前提。目前，随着计算机辅助诊断的普及，已经有很多肺野自动分割算法被提出来了，分割准确度与速度也越来越高，但仍存在一些难点：锁骨和肋骨窗的强边界，容易导致局部极小值；极小的胸膈角难以准确分割等。

目前胸片肺野分割算法主要有基于规则的分割^[1]、像素分类、活动形状模型等。基于规则的肺野分割，通常只得到较为粗糙的肺野轮廓，可用做活动形状模型等的初始解。基于像素分类的分割取决于分类器预测性能，一般难以准确定位肺野边界。形变模型^[2-3]由于具有灵活的形变特性而广泛应用于医学图像分割^[4]，但也存在一些不足：（1）对初始解敏感，易陷入局部极小值点；（2）对算法参数敏感，造成性能不稳定。

为实现准确、可靠的肺野分割，本文提出了一种基于图像特征密集匹配和标号融合的肺野分割方法。提取待分割胸片图像的局部特征，将像素点的局部特征与肺野已知的参考图像集中的局部特征进行匹配，确定特征在参考图像集中的最近邻，并由此将参考图像中的肺野标号信息通过标号融合迁移至待分割的胸片图像，以实现肺野的自动分割。

实验数据采用公开的JSRT胸片图像数据集^[5]，包括154幅肺结节图像和93幅无结节图像，其中每幅胸片均有专家手工勾画的肺野轮廓。JSRT数据集中胸片原始图像大小为2048×2048，像素大小0.175 mm，灰度深度为12比特。为了与文献结果比较和减少运算时间，原始胸片图像降采样为256×256大小进行实验和性能评价。本文中算法采用Matlab编程实现。

本文方法分为两部分实现：构建参考图像特征集合和通过密集特征匹配分割肺野。算法基本流程见图 1。给定参考图像集合及其对应的分割标号图集合，参考图像特征集合构建过程如下：①将所有参考图像进行灰度归一化处理；②取所有图像中肺野外接矩形内的区域，将其长宽变换为固定之后叠加平均，形成模板图像T；③将所有参考图像与模板图像T进行对齐，对参考图像标号图像作相应的空间变换，提取局部特征描述子并降维，得到参考图像特征集合以及对应的分割标号图像集合。

图 1(Figure 1)

图 1 本文方法基本流程 Figure 1 Flow chart of the proposed method

在构建参考图像特征几何完成后，对一给定的胸片图像I，对其进行分割的基本过程如下：①对输入图像I进行灰度归一化处理；② 将图像I与模板图像T进行对齐；③提取图像I的局部描述子并降维；④ 进行K次密集局部特征匹配，得到K个近邻场；⑤利用步骤4得到的近邻场和局部特征的相似性进行标号融合；⑥将标号融合结果进行空间逆变换，使得标号融合结果与原始输入图像在空间上一致；⑦将标号融合结果进行阈值化处理得到肺野分割结果，阈值可设为0.5。

为减小胸片之间对比度差异对特征匹配造成的影响，对胸片灰度进行归一化处理。设胸片灰度直方图中的最后一个和第一个峰值对应的灰度分别为p_max和p_min，归一化处理后的图像灰度为p_i =(p_i - p_min)/(p_max - p_min)，其中p_i表示像素i的灰度值。为了将特征匹配的搜索范围限制在较小区域内，需将胸片图像与模板进行空间对齐，使得对齐后的胸片中肺野位于大致相同的位置。本文采用仿射变换作为形变模型，通过最大化增强相关系数^[6]进行胸片与模板的对齐。

图像块匹配是指对图像块在图像中搜索其最相似的图像块（最近邻），图像块与其近邻之间的位移关系由最近邻场表达。PatchMatch^[7]是一种高效的图像块匹配算法。由于图像的一致性，相邻的块与它们的最近邻可能具有同样的位移向量，因此PatchMatch算法可以将最近邻场传播到相邻的区域。PatchMatch算法分为三步：随机初始化、近邻场传播、随机搜索。随机初始化指的是对图像A中的每一个图像块都在图像B中随机地指定一个最近邻，得到初始最近邻场；近邻场传播可以把好的位移关系传到周围的区域，从而逐步优化最近邻场。以像素点p_（i,j）为例，设A_（i,j）是图A中以像素点p_（i,j）为中心取的一个小块，当传播到点p_（i,j）时，设块B（_m,n）是块A_（i,j）的最近邻，且与A_（i,j）相邻的块A_（i-1,j）、A_（i,j-1）都已经得到了更新的最近邻，分别设为B_（s,t）、B_（x,y），则B_（s+1,t）、B（x,y+1）有可能是A_（i,j）的最近邻，分别计算A_（i,j）与B_（m,n）、B_（s+1,t）、B_（x,y+1）这三个块的欧式距离，其中距离最小的则作为A_（i,j）传播后的最近邻，这里假设B_{（s +1,t）}为A_（i,j）传播后的最近邻；随机搜索指的是在以块B_（s+1,t）的中心像素点p_（s+1,t）为中心的搜索窗内随机选出一个块B_（u,v），然后再从中择优选出一个作为A_（i,j）的最近邻，搜索窗的大小依次为W=m× n × (1/4)^{t - 1}，m和n分别表示图像A的高、宽，t是随机搜索次数，直到W=1；最后多次迭代近邻场传播和随机搜索的这两步，最近邻场会逐渐收敛。由于随机初始化后，图像A中至少有一个块得到了较好的最近邻的这种可能性很大，所以，经过多次迭代的传播和随机搜索，最后A中绝大部分的块都能找到合理的最近邻。

PatchMatch算法可以高效求解两图像间图像块的近邻场。但当参考图像数量较多时，逐一将输入图像与参考图像进行块匹配是不实际的。为了解决这一问题，Golland^[8]提出了优化的块匹配标号融合算法，用于解剖结构的分割。值得注意的是，块匹配算法的复杂度取决于图像A的大小而不是B，优化的块匹配算法是从参考图像集中搜索K个近邻，而不是只从与图像A相似的图像B中找一个近邻。标号融合的时候，每个近邻的标号的权重取决于匹配对之间的相似度，跟块匹配算法一样，不是只取中心像素点的标号，而是基于中心点取了一个块。另外，随机初始化和随机搜索都被限制于一个固定大小的搜索窗内。由于不用预先检索出相似图像B，这个优化算法效率更高，而且由于有搜索窗的限制，减小了误匹配率，实现了更高的准确率。

对于胸片图像，由于图像采集条件和病人个体差异，不同胸片肺野灰度和对比度差异较大。直接通过图像块匹配和标号融合，难以得到可靠的分割结果。为了解决这个问题，本文采用原始图像块和密集尺度不变特征转换描述子（SIFT）^[9]作为图像的局部特征。原始像素块是一种具有较强区分能力的描述子，但对噪声和形变很敏感，为了平衡描述子的不变性和区分能力，原始图像块的尺寸不能过大。SIFT描述子在计算机视觉有着广泛的应用，具有很好的旋转不变性，通常用于描述局部区域的形状特征，对于图像的形变和对比度能够保证一定的不变性。本文中，原始图像块的大小取9×9，SIFT 描述子的维度是128，局部特征的维数为209。为了减小计算复杂度，采用了主成分分析^[10]对特征向量进行了降维。本文中，局部特征的维数降为40，保持总方差的95%。

特征匹配与块匹配不同，对于输入图像的每一局部特征，直接在所有参考图像的局部特征中搜索其对应的K个近邻。设共有N幅参考图像，参考图像大小为m×n，局部特征的维数为D，则参考图像特征集合可用大小为D×N×m×n的四阶张量表达，参考图像特征集合可通过离线方式提取并保存。特征匹配搜索近邻的步骤与PatchMatch算法相似，分3步进行迭代：受限随机初始化、近邻场传播、受限随机搜索，其中近邻场传播仅在图像平面内进行。

如图 2所示，右边的大长方体表示特征向量集，圆点代表该点对应的特征，浅蓝色方框表示固定大小的搜索窗，t表示参考图像序号，范围是1到N。随机初始化时，输入图像特征在参考图像特征集合中随机指定一个近邻。比如图 2所示，输入图像相邻点的近邻可以位于不同参考图像中。近邻场传播的方法与PatchMatch算法类似，仅在参考图像的图像平面内进行传播，即近邻场传播时，近邻场对应的参考图像序号保持不变。随机搜索时，对输入图像任一点随机选择参考图像序号，在以当前近邻为中心的搜索窗内随机选取一点，如果随机选取的特征比当前近邻与输入图像特征之间的距离小，则将随机选取的近邻作为更新后的近邻。特征之间距离采用欧式距离。

图 2(Figure 2)

图 2 密集特征匹配示意图 Figure 2 General view of dense matching of local features. Left:Input Image,Right: Features Database of reference images

重复特征匹配过程K次，对于输入图像可以得到K个近邻场（nnf），记录了输入图像中像素点特征在参考图像特征集合中K个近邻的空间位置和参考图像序号信息。在标号融合过程中，本文采用基于图像块的标号融合方法^[11-13]。设标号图像中像素值为0或者1，其中1表示像素点为肺野。对于输入图像的任一像素点（x,y），在标号图像中取以其近邻为中心的小块（如5×5大小的图像块），并以点（x,y）的特征v_(x,y)与其近邻的特征v_{(nnf(i,(x,y)))}的相似度作为标号块的权重，由K近邻加权得到的标号块为

其中l_i(nnf_(i,(x,y)))表示以点（x,y）在第i个近邻场nnf_i对应的近邻nnf_(i,(x,y))为中心取的标号块，权重

${{\sum\limits_{d=1}^{D}{\left( v{{1}_{d}}-v{{2}_{d}} \right)}}^{2}}$表示两个D维向量的欧式距离的平方，另外，s的取值要适中，过大则权重都约为1，过小则权重都约为0，实验中可取为待分割图像中所有像素点与其近邻的距离的均值。最后，将加权的标号块进行重排和叠加为输入图像对应的标号图，并进行归一化处理。

为了得到与原始输入胸片空间位置对应的分割，还需将标号融合结果用模板匹配的空间变换进行逆变换。标号融合的结果可以看做肺野的概率图，对其进行阈值化或使用图割方法^[14]进行处理即可得到肺野的分割。

实验中，将JSRT胸片数据集分为大小相同的两个子集，分别用作训练集和测试集。肺野分割性能指标包括Jaccard指标（JI）和平均边界距离（MBD），Jaccard值越大、MBD值越小则表明分割的性能越好，其中：

这里GT、SM分别表示手动分割的和自动分割的标号图，1为肺野，0为背景，公式（2）计算的是这两个标号图的肺野区域的重叠率，重叠率越高表示分割效果越好；公式（3）计算的是GT、SM中肺野的边界GTB、SMB之间的欧式距离的平均值，当手动分割的肺野区域与自动分割的完全重叠，Jaccard值为1，MBD值为0，重叠率Jaccard 越高时，MBD 值越小，表示分割效果越好，D(g_i,SMB) =$\underset{j}{\mathop{\min }}\,\left\| {{g}_{i}}-{{s}_{j}} \right\|$表示GTB 上的任意点g_i 到SMB的最短距离，N_g、N_s分别表示GTB、SMB上像素点的总个数。

图 3反映了特征匹配近邻数K对Jaccard指标的影响，c-d 分别是近邻数K等于1、2、5、7、10 时的分割结果。图 4显示的是不同近邻数时标号融合后的结果，左中至右近邻数分别是1、5、10。可以看出近邻数增加时，分割得到的肺野轮廓更平滑、分割结果更准确。

图 3(Figure 3)

图 3 近邻数不同的标号融合结果 Figure 3 Result of label fusion with different numbers of nearest neighbors.

图 4(Figure 4)

图 4 近邻数K对分割结果的影响 Figure 4 Effect of K (number of nearest neighbors) on segmentation results.

图 5主要反映了密集匹配过程中迭代次数对分割性能的影响，横轴表示近邻数。图中可见Jaccard值在第二次迭代后明显提高，迭代次数大于3且近邻数大于6时，Jaccard值逐渐趋于稳定。由于迭代次数与近邻数越大，需要的运算时间会增加，平衡效率和分割性能，本文选取的迭代次数为3，近邻数为10。

图 5(Figure 5)

图 5 密集匹配迭代次数的影响 Figure 5 Effect of iterations of dense matching.

图 6 反映的是权重系数s 对Jaccard 指标的影响。当s很小，比如0.1时，对Jaccard的影响较大，尤其是近邻数较小的时候；当s为所有像素点与其近邻的欧式距离的均值（即图中mean（Dist），约为9.5），和s等于10或100时，Jaccard值几乎完全相等；所以，为省去计算s的时间，本文直接选取s为10。

图 6(Figure 6)

图 6 权重系数s的影响 Figure 6 Effect of weight factor.

图 7是JSRT数据集中部分图像自动分割的结果与手工分割结果的对比。图中红线和绿线分别为手工勾画和本文方法得到的肺野轮廓；第1行是容易分割且效果很好的例子，第2行是较难分割且分割效果不太理想的例子。JI代表Jaccard指标。表 1列出了本文方法和其他肺野分割方法在JSRT数据集上的分割性能，本文方法对应的Jaccard指标优于其他方法。表 2列出了各算法用相应的编程语言实现时运行时间的比较。

图 7(Figure 7)

图 7 JSRT数据集中的图像分割结果实例演示 Figure 7 Instances of segmentation results from JSRT database.

表 1(Table 1)

表 1 肺野分割算法性能比较 Table 1 Comparison of performances of lung field segmentationalgorithms Algorithms Jaccard (%) MBD (mm) Proposed method with Graph-cut post-processing 95.5±1.6 1.39±0.57 Proposed method with threshold post-processing 95.4±1.6 1.41±0.58 Atlases with Non-rigid Registration[15] 95.4±1.5 1.31±0.32 Human observer[16] 94.6±1.8 1.64±0.69 Joint shape and appearance sparse learning[17] 94.6±1.9 1.67±0.76 ASM tuned[16] 92.7±3.2 2.30±1.03

表 1 肺野分割算法性能比较 Table 1 Comparison of performances of lung field segmentationalgorithms

表 2(Table 2)

表 2 肺野分割运行时间的比较 Table 2 Comparison of execution time on lung fields segmentation Algorithms Execution time Programming language Proposed method with Graph-cut post-processing 12-14 s Matlab Proposed method with threshold post-processing 11-13 s Matlab Atlases with Non-rigid Registration[15] 20-25 s C++、Matlab Human observer[16] - - Joint shape and appearance sparse learning[17] 35.2 s Matlab ASM tuned[16] 1 s C++

表 2 肺野分割运行时间的比较 Table 2 Comparison of execution time on lung fields segmentation

为了验证本文方法的有效性和推广能力，采用本文方法对两种不同X线DR设备上采集的数字化胸片进行自动分割。图 8显示了DR胸片和自动分割肺野的轮廓（红线），可以看出分割结果合理，并且肺野轮廓定位准确。虽然JSRT数据集中的胸片都是通过扫描胶片获取的，但人工分割的结果通过本文方法仍可很好地迁移至DR胸片的分割，说明本文方法具有良好的推广能力。

图 8(Figure 8)

图 8 数字胸片自动分割实例 Figure 8 Instances of automatic segmentation in DR

胸片中细小的胸膈角和受锁骨影响的肺肩区域的准确分割比较困难，并且肋骨窗的强对比度会影响分割结果的准确性。从实验结果看，本文方法在一定程度较好解决了这些问题。从图 8中分割结果较差的图像来看，严重胀气的结肠或周围组织被错分为肺野区域，而由于大量液体存在而模糊的肺野区域被错分为非肺野，这些问题尚需进一步解决。

本文提出了一种采用图像特征密集匹配和标号融合的胸片肺野分割方法。与一般的两图像间块匹配方法不同，本文采用密集SIFT和原始图像块作为图像的局部特征，直接将输入图像特征与参考图像集合的特征进行匹配，避免了参考图像子集的选择问题，有效利用所有参考图像及其标号的信息。在JSRT胸片图像数据集上，本文方法得到了目前为止最高的分割精度。此外，实验结果还表明本文方法具有很好的推广能力。在进一步的工作中，我们将对特征匹配方法进行推广应用，如用于胸片中结节的自动检测和肋骨的抑制。