CN104424629A - 一种x光胸片肺部分割方法和装置 - Google Patents

Abstract

本发明公开了一种X光胸片肺部分割方法，以对X光胸片上的肺部进行精确分割。该方法包括：S101，通过水平垂直投影，获取分别包围X光胸片中左右肺部图像的两个矩形区域；S102，在两个矩形区域内初始化肺部以获得肺部的初始形状；S103，根据加权灰度局部纹理模型，寻找肺部图像中各个特征点的最佳匹配点；S104，通过调整姿态参数和形状参数b，使得肺部当前形状I_Xc最大限度地接近I_X+dI_X；重复S103和S104，直至所获肺部形状与紧邻本次调整前一次调整所获肺部形状相比，两者的改变量小于预定阈值。本发明能够获得较好的初始化肺部形状，在后续的调整过程中不会造成过分割现象，在主动形状模型的约束下能够对X光胸片上的肺部区域进行准确分割。

Description

一种X光胸片肺部分割方法和装置

技术领域

本发明涉及生物医学图像领域，具体涉及一种X光胸片肺部分割方法和装置。

背景技术

X光胸片检查仍然是目前肺部疾病初始检测的最主要手段。在临床应用中，胸片图像大量应用于肺部疾病的诊断，很多胸片图像分析应用的第一步是标示出图像中肺部的边界，以便计算机能进行下一步的处理，因此，肺部分割显得非常重要。首先是因为肺部面积占了整个X光胸片的1/3，若在诊断时提取肺部感兴趣区域，则会大大提高诊断的准确率；其次是因为在肺部分割后，再做其他处理，会大大减少其它不相关组织的干扰，提高了诊断的速度；最后是因为肺部分割更有利于分析肺部纹理和肺部的生理结构。很早以前，研究人员就已经开始做这方面的研究了，但大多数的研究都仅仅是为了给后续处理提供方便，并没有完全将眼光停留在肺部分割上。例如，为了计算心胸比率，可以只提取胸廓最大横径对应的两个端点；为了检测肺部形状异常，可能仅需要提取横隔膜与胸廓的夹角；只有要检测整个肺内部的病变时，才需要提取全部的肺区域或者肺轮廓，而在这种情况下，得到一个粗略的肺部分割结果往往己经足够。鉴于对X光胸片肺部进行精确分割具有重要的临床意义，因此，研究精确的肺部分割十分必要。

目前肺部的分割大体分为基于规则的推理和像素归类两种方法。基于规则的方案使用的方法是（局部）阈值化、区域增长、边缘探测、形态学操作、拟合几何模型和函数、动态规划等，大多数肺部区域分割算法属于这类。基于像素归类的方案试图将图像中的每个像素归类到一种解剖类型（通常是肺部和背景，但是在有些情况下，使用的是更多的种类，例如，还包括心脏、中隔膜和横隔膜）。使用的分类器为各种神经网络或Markov随机场模型，并用各种（局部）特征进行分类，包括灰度、位置、纹理测量。

在上述现有的两类方法中，一种X光胸片肺部分割方法是先用自适应的灰度阈值进行处理，紧接着又用基于形状的方法来标记胸片中的肺部边界。在这种方法中，通过计算靠近肺部边界的局部像素直方图以得出一个初始的、近似的分割。该直方图用来选择一个局部阈值以便将肺部组织从相邻的结构中分离出来。基于形状的特殊处理，用解剖知识作为导引，平滑肺部边界并描绘出肺部的左下角和右下角。然而，这种分割方法分割出来的边缘不规则，而且受噪声影响较大，在图像质量不佳的情况下很容易失效，出现分割面积过小或者过大的现象。

在上述现有的两类方法中，另一种X光胸片肺部分割方法是基于规则的边界检测法，这种方法是一种最直观的肺边界提取方法，其主要是基于这样一种事实：所有边界之中一定包含目标边界，并且这些边界有其固定的特点。边界通常定义为沿梯度方向梯度幅度的极值点，可以用导数来表达某个固定方向的极值。导数用图像I（x，y）和高斯导数的卷积计算。使用合适的尺度、导数阶数和方向可以检测X光胸片中与解剖结构对应的边界，这些边界互不相连，例如，胸廓中心线和肺中心线、纵膈边界、胸腔边界等。然而，基于规则的边界检测法的缺陷在于只给出了主要的肺部边界片段，没有给出完整的肺部轮廓。

发明内容

本发明实施例提供一种X光胸片肺部分割方法和装置，以对X光胸片上的肺部进行精确分割。

一种X光胸片肺部分割方法，所述方法包括如下步骤：

S101，通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域；

S102，在所述两个矩形区域内初始化肺部以获得肺部的初始形状；

S103，根据加权灰度局部纹理模型，寻找肺部图像中各个特征点（x_i，y_i）的最佳匹配点（x’_i，y’_i），记dI_X=（d x₁，d y₁，d x₂，d y₂，…，d x_i，d y_i，…，d x_n，d y_n）^T，所述d x_i=x’_i－x_i，所述d y_i=y’_i－y_i，所述i=1，2，…，n，所述肺部图像包括调整阶段的肺部形状和所述肺部的初始形状；

S104，通过调整姿态参数和形状参数b，使得肺部当前形状I_Xc最大限度地接近I_X+dI_X，所述形状参数b的调整结果为调整后的值b+db受主动形状模型的形状参数的约束，所述I_X是所述步骤S103执行时的肺部形状，db是b的改变量；

将所述经过调整姿态参数和形状参数b后的肺部形状作为步骤S103的执行对象，重复步骤S103和步骤S104，直至所获肺部形状与紧邻本次调整前一次调整所获肺部形状相比，两者的改变量小于预定阈值。

一种X光胸片肺部分割装置，所述装置包括：

获取模块，用于通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域；

初始化模块，用于在所述两个矩形区域内初始化肺部以获得肺部的初始形状；

最佳匹配点找寻模块，用于根据加权灰度局部纹理模型，寻找肺部图像中各个特征点（x_i，y_i）的最佳匹配点（x’_i，y’_i），记dI_X=（d x₁，d y₁，d x₂，d y₂，…，d x_i，d y_i，…，d x_n，d y_n）^T，所述d x_i=x’_i－x_i，所述d y_i=y’_i－y_i，所述i=1，2，…，n，所述肺部图像包括调整阶段的肺部形状和所述肺部的初始形状；

形状调整模块，用于通过调整姿态参数和形状参数b，使得肺部当前形状I_Xc最大限度地接近I_X+dI_X，所述形状参数b的调整结果为调整后的值b+db受主动形状模型的形状参数的约束，所述I_X是所述最佳匹配点找寻模块执行时的肺部形状；

所述经过形状调整模块调整姿态参数和形状参数b后的肺部形状作为最佳匹配点找寻模块的执行对象，重复最佳匹配点找寻模块和形状调整模块的执行过程，直至所获肺部形状与紧邻本次调整前一次调整所获肺部形状相比，两者的改变量小于预定阈值。

从上述本发明实施例可知，由于可以通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域，在所述矩形区域内得到初始化肺部形状，通过迭代方式调整肺部图像时将形状参数b的调整约束在主动形状模型提供的形状参数的范围之内。与现有技术提供的X光胸片肺部分割方法相比，本发明实施例提供的方法能够获得较好的初始化肺部形状，在后续的调整过程中不会造成过分割现象（例如，将脊柱和胃泡识别为肺区），在主动形状模型的约束下能够对X光胸片上的肺部区域进行准确分割，为临床或者计算机分析提供有价值的数据。

附图说明

图1是本发明实施例提供的X光胸片肺部分割方法的基本流程示意图；

图2是本发明实施例提供的通过水平投影得到的左肺右边界和右肺左边界示意图；

图3是本发明实施例提供的通过水平投影得到的左肺右边界和右肺左边界示意图；

图4是本发明实施例提供的通过垂直投影得到的右肺顶部边界和底部边界示意图；

图5是本发明实施例提供的在肺顶部边界和底部边界限定的范围内通过再次水平投影得到的右肺左边界和左肺右边界示意图；

图6是本发明实施例提供的在右肺左右边界内一个窗口限定的范围内通过再次垂直投影得到的右肺顶部边界示意图；

图7是本发明实施例提供的分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域示意图；

图8是本发明实施例提供的在矩形区域内获得的肺部初始形状示意图；

图9是本发明实施例提供的通过迭代搜索算法收敛后获得的肺部形状轮廓示意图；

图10是本发明实施例提供的对样本标记点向量进行归一化配准处理时配准前和配准后示意图；

图11是本发明实施例提供的X光胸片肺部分割装置基本逻辑结构示意图；

图12是本发明另一实施例提供的X光胸片肺部分割装置基本逻辑结构示意图；

图13是本发明另一实施例提供的X光胸片肺部分割装置基本逻辑结构示意图；

图14是本发明另一实施例提供的X光胸片肺部分割装置基本逻辑结构示意图；

图15-a是本发明另一实施例提供的X光胸片肺部分割装置基本逻辑结构示意图；

图15-b是本发明另一实施例提供的X光胸片肺部分割装置基本逻辑结构示意图；

图15-c是本发明另一实施例提供的X光胸片肺部分割装置基本逻辑结构示意图；

图15-d是本发明另一实施例提供的X光胸片肺部分割装置基本逻辑结构示意图。

具体实施方式

本发明实施例提供一种X光胸片肺部分割方法，包括：S101，通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域；S102，在所述两个矩形区域内初始化肺部以获得肺部的初始形状；S103，根据加权灰度局部纹理模型，寻找肺部图像中各个特征点（x_i，y_i）的最佳匹配点（x’_i，y’_i），记dI_X=（d x₁，d y₁，d x₂，d y₂，…，d x_i，d y_i，…，d x_n，d y_n）^T，所述d x_i=x’_i－x_i，所述d y_i=y’_i－y_i，所述i=1，2，…，n，所述肺部图像包括调整阶段的肺部形状和所述肺部的初始形状；S104，通过调整姿态参数和形状参数b，使得肺部当前形状I_Xc最大限度地接近I_X+dI_X，所述形状参数b的调整结果为调整后的值b+db受主动形状模型提供的形状参数的约束，所述I_X是与主动形状模型空间中肺部形状X对应的图像空间中实际形状；将所述经过调整姿态参数和形状参数b后的肺部形状作为步骤S103的执行对象，重复步骤S103和步骤S104，直至所获肺部形状与紧邻本次调整前一次调整所获肺部形状相比，两者的改变量小于预定阈值。本发明实施例还提供相应的X光胸片肺部分割装置。以下分别进行详细说明。

本发明实施例的X光胸片肺部分割方法的基本流程可参考图1，主要包括步骤：

S101，通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域。

作为本发明一个实施例，可以通过将所述X光胸片中左肺部图像和右肺部图像进行水平和垂直投影，获得投影曲线，在所述投影曲线中找到各个肺部的左边界和右边界以及其中任一肺部的顶部边界和底部边界对应的坐标点，以所述左肺部的左边界、右边界、顶部边界和底部边界对应的坐标点围成的矩形区域作为X光胸片中左肺部图像的矩形区域，以所述右肺部的左边界、右边界、顶部边界和底部边界对应的坐标点围成的矩形区域作为X光胸片中右肺部图像的矩形区域。即，对于由左肺部图像和右肺部图像构成的尺寸为h*w的图像I(x,y)，水平投影时，其中，r为y的取值总数，垂直投影时，其中，s为x的取值总数。具体地，通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域包括如下：

S1011，对图像I(x,y)上1/3部分进行水平投影找脊柱左右边界，即对于水平投影h(x)，x在[0，w]内取值，y在[0，h/3]内取值，如图2所示，在投影曲线中找到两个对应点坐标即为脊柱左右边界，其中脊柱左边界也为右肺右边界、脊柱右边界也为左肺右边界；

S1012，对整幅图像I(x,y)水平投影，即对于水平投影h(x)，x在[0，w]内取值，y在[0，h]内取值，如图3所示，在投影曲线中找到两个对应点坐标即为右肺左边界和左肺右边界；

S1013，对于右肺，限定范围为右肺左边界(假设左边界x坐标为x_l)和右边界(假设右边界x坐标为x_r)进行垂直投影，即对于垂直投影v(y)，x在[x_l，x_r]内取值，y在[0，h]内取值，如图4所示，在投影曲线中找到两个对应点坐标即为右肺顶部边界和底部边界；

S1014，限定范围为肺顶部边界(假设顶部边界的y坐标为y_t）和底部边界(假设底部边界的y坐标为y_b)，再次进行水平投影，即对于水平投影h(x)，x在[0，w]内取值，y在[y_t，y_b]内取值，如图5所示，在投影曲线中找到两个对应点坐标即为右肺左边界和左肺右边界；

S1015，限制范围为右肺左右边界内一个窗口（中心线两边各1/4宽）再次垂直投影，即对于垂直投影v(y)，x在内取值，y在[0，h]内取值，如图6所示，在曲线中找到对应点坐标即为右肺顶部边界。

本发明实施例中，分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域如附图7所示。

S102，在步骤S101获取的两个矩形区域内初始化肺部以获得肺部的初始形状。

肺部的初始形状即主动形状模型提供的形状参数b=0时获得的肺部初始形状模型对于主动形状模型特征空间中的初始形状模型X，通过相似性变换将该初始形状模型X变换到图像坐标系下的水平垂直投影确定的左右肺矩形区域内，得到图像坐标系下的初始形状I_X，具体包括：将主动形状模型特征空间中平均形状从一维形状向量转化为二维空间中n个标记点坐标（x₁，y₁），（x₂，y₂），…，（x_n，y_n），所述找到包围该标记点坐标的最小矩形区域；根据所述最小矩形区域和步骤S101获取的矩形区域，确定变换参数(u₀,v₀,t_x0,t_y0)，将所述n个标记点坐标变换到所述矩形区域所属图像坐标系下，矩形区域内获得的肺部的初始形状如附图8所示。

S103，根据加权灰度局部纹理模型，寻找肺部图像中各个特征点（x_i，y_i）的最佳匹配点（x’_i，y’_i）。

记dI_X=（d x₁，d y₁，d x₂，d y₂，…，d x_i，d y_i，…，d x_n，d y_n）^T，d x_i=x’_i－x_i，d y_i=y’_i－y_i，i=1，2，…，n，肺部图像包括调整阶段的肺部形状和肺部的初始形状。作为本发明一个实施例，根据加权灰度局部纹理模型，寻找肺部图像中各个标记点（x_i，y_i）的最佳匹配点（x’_i，y’_i）包括：根据加权灰度局部纹理模型，获取使得代价函数f(g_i)取最小值的标记点作为标记点（x_i，y_i）的最佳匹配点（x’_i，y’_i），其中， $f\left(g_i\right)=100*{(g_i-\stackrel{\‾}{g})}^T{S}_g^{-1}(g_i-\stackrel{\‾}{g})-{\mathrm{\λ}}_1*l\left(i\right)+{\mathrm{\λ}}_{2}d\left(i\right),$ λ₁∈[0.7，1.3]，λ₁∈[5，10]，l(i)=g(x_i,y_i)，g(x_i,y_i)为对肺部图像进行－45°角线性滤波后的肺部图像在标记点处的灰度值，只取左肺横隔膜部分， $d\left(i\right)=\frac{1}{2}(\sqrt{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}+\sqrt{{(x_i-x_{i+1})}^2+{(y_i-y_{i+1})}^2}),$ i=21，…，29，g_i为加权灰度局部纹理模型提供的训练样本第i个边界点（x_i，y_i）两侧沿法线方向各取k个点连同所述第i个边界点（x_i，y_i）共2k+1个点的灰度值构成的向量。此处，－45°角线性滤波核如下：

$f=\left[\begin{array}{ccc}0& -1& -2\\ 4& 0& -4\\ 2& 1& 0\end{array}\right],$ 显然，从f(g_i)的表达式可知：d(i)越大，则会使f(g_i)越大，从而对采样点中离图像实际点比较远的点有一个很好的约束。

S104，通过调整姿态参数和形状参数b，使得肺部当前形状I_Xc最大限度地接近I_X+dI_X，其中，形状参数b的调整结果为调整后的值b+db受主动形状模型提供的形状参数的约束，I_X是与主动形状模型空间中肺部形状X对应的图像空间中实际形状，这里db表示b的改变量。

需要说明的是，在本发明实施例中，步骤S103和步骤S104是通过迭代搜索方式来实现，即，以肺部的初始形状作为迭代搜索算法的起始对象，每次经过步骤S104调整姿态参数和形状参数b后，所得的肺部形状作为步骤S103的执行对象，重复执行步骤S103和步骤S104，直至所获肺部形状与紧邻本次调整前一次调整所获肺部形状相比，两者的改变量小于预定阈值，迭代搜索算法结束。

以下先介绍主动形状模型（Active Shape Model，ASM）迭代搜索算法的理论基础包括如下1）至3）：

1）假设当前形状为I_X,有I_X=M(s,θ)X+t，其中X表示主动形状模型空间中轮廓形状向量(由k个主分量决定的一个中心位于原点，向量模值为1的2n维向量)，表示主动形状模型特征空间的平均形状，I_X表示与主动形状模型空间中形状向量X对应的图像空间中实际形状向量，X经过一个相似性变换后可以得到I_X，即I_X=M(s,θ)X+t，s为相似性变换T(u,v,t_x,t_y)=M(s,θ)+t（包括尺度旋转变换和平移变换）中的尺度参数，θ为其中的旋转参数，t=(t_x,t_y)^T为其中的平移参数；

2）对当前形状I_X的每个标记点周围进行搜索,寻找各标记点（x_i，y_i）的最佳匹配点（x’_i，y’_i）得到新形状I_X+dI_X，有dI_X=（d x₁，d y₁，d x₂，d y₂，…，d x_i，d y_i，…，d x_n，d y_n）^T,d x_i=x’_i－x_i，d y_i=y’_i－y_i，i=1，2，…，n，最佳匹配点为使代价函数f(g_i)取最小值的特征点，代价函数f(g_i)已在前述实施例给出；

3）要得到当前最佳形状使其最适合I_X+dI_X，则要调整当前的姿态参数(s,θ,t_x,t_y)和形状参数b,假设有其中X₁=X+dX，首先要找到姿态参数1+ds、1+dθ和dt然后再调整形状参数b使其满足主动形状模型约束。假设现在已知姿态参数1+ds、1+dθ和dt，则：由公式1： $I_{X_1}=M(s(1+\mathrm{ds}),\mathrm{\θ}(1+\mathrm{d\θ}))X_1+t+\mathrm{dt}=I_X+{\mathrm{dI}}_X,$ 公式2：X₁=X+dX，公式3：I_X=M(s,θ)X+t，可推出：M(s(1+ds),θ(1+dθ))[X+dX]=M(s,θ)X+dI_X-dt。由于M^-1(s,θ)[…]=M(s^-1,-θ)[…]，得出dX=M((s(1+ds))^-1,-θ(1+dθ))[M(s,θ)X+dI_X-dt]-X。

通常所得到的向量dX是2n维空间的，但是，在主动形状模型中只描述了k个分量的变化，有 $X+\mathrm{dX}\≈\stackrel{\‾}{X}+\mathrm{\Φ}(b+\mathrm{db})\⇒\mathrm{db}={\mathrm{\Φ}}^{-1}\mathrm{dX}={\mathrm{\Φ}}^T\mathrm{dX},$ 则得到新的形状参数(b+db)，注意新形状参数要根据主动形状模型约束适当调整。此时，可以得到新的形状估计值 $I_{X_1}=M(s(1+\mathrm{ds}),\mathrm{\θ}(1+\mathrm{d\θ}))[X+\mathrm{\Φdb}]+t+\mathrm{dt}.$ 对于这个新估计值，继续进行上述过程会依次得到I_X2，I_X3，…等等，直到前后两次形状的改变量很小，例如小于某一个阈值为止。

参数的每次更新如下所示：

t_x→t_x+dt_x

t_y→t_y+dt_y

s→s+ds

θ→θ+dθ

b→b+db

按照上述ASM迭代搜索算法的理论基础，步骤S103和步骤S104的具体实现过程包括如下步骤1）至步骤7）：

1）首先，给出初始形状估计值I_X=M(s,θ)X+t，其中即给定主动形状模型空间中的初始轮廓X，得到初始轮廓在图像空间的映射；

2）设当前形状为I_X；

3）寻找所需位移量dI_X，使标记点位移至最佳位置，即在图像空间中通过使代价函数f(g_i)取最小值，搜索每个标记点的最佳匹配点得到dI_X=(dx₁,dy₁,dx₂,dy₂,…,dx_n,dy_n)^T；

4）求出姿态参数的改变量1+ds、1+dθ和dt，使尽可能接近I_X+dI_X，该步具体如下所示：先从图像空间映射到主动形状模型空间，再由主动形状模型空间映射回图像空间。在这个过程中，由于对齐和形状约束的原因，实际的变换参数不再是s、θ、t_x和t_y，而会有一个增量的变化，即为s(1+ds)、θ(1+dθ)、t_x+dt_x和t_y+dt_y。

5）求出形状参数b的改变量db，并且使得b+db满足主动形状模型约束，即对于k个分量都要满足 $-3\sqrt{{\mathrm{\λ}}_i}\≤(b_i+{\mathrm{db}}_i)\≤3\sqrt{{\mathrm{\λ}}_i},i=1,...,k;$

6）更新作用于X上的形状参数和姿态参数，得到新的形状估计值，以该新的形状估计值作为步骤2）的当前形状，从步骤3）开始继续处理；。

7）重复步骤3）至6）的过程，直至肺部形状改变很小（例如，前后两次的改变量小于某个设定的阈值）后停止迭代搜索算法即迭代收敛。在本发明实施例中，迭代搜索算法收敛后获得的肺部形状轮廓如附图9所示。

在上述实施例提供的方法之前，可以根据图像训练集建立主动形状模型、加权灰度局部纹理模型和邻近距离约束模型等等，以下对建立上述模型的过程进行说明如下，主要包括如下步骤：

S201，手动标记肺部轮廓。

首先由专家手动分割出训练图像肺部轮廓，然后用轮廓边界点即标记点来代表肺部形状轮廓。边界点一般选在能表示目标物体外形特征的地方，例如拐角点、T形连接处等，同时在具有外形特征的边界点中间，用等间距采样的方式选取另外一些中间边界点，共同构成目标轮廓的边界点集合。对于右肺，取右肺肺尖、左肋隔角，左心膈角三个特征点，然后在特征点中间等距离取样，总共取44个标记点来代表右肺轮廓；对于左肺，取左肺肺尖、右肋膈角、心膈角、主动脉弓上下端点五个特征点，然后在特征点中间等距离取样，总共取50个标记点来代表左肺轮廓。如此，n个边界点构成的边界点集合可以表示成长度为2n的向量X，即：X=(x₁,y₁,x₂,y₂,…,x_n,y_n)^T，对L幅X光胸片图像，用上述方法进行边界轮廓点提取，得到训练样本集Ω，则Ω={X₁,X₂,…X_L}。

S202，配准标记点。

由于不同X光图像中肺部的外形和位置会存在较大偏差，因此，在对Ω中向量进行统计分析之前，有必要对所有样本标记点向量进行归一化配准处理。

首先介绍两个肺部形状轮廓的标记点向量之间的配准问题，对于Ω中任意两个向量X和X′，可以通过将X进行相似性变换T(X)后，将X对齐到X′，要计算出该相似性变换参数，则要使下式中的E值达到最小：

E=|T(X)-X′|²

$T\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}u& -v\\ v& u\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right),$ 其中u=scosθ，v=ssinθ，s为尺度变换参数，θ为旋转角度参数。

由上面两式可得

$E(u,v,t_x,t_y)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[{({\mathrm{ux}}_i-{\mathrm{vy}}_i+t_x-x_i^{\′})}^2+{({\mathrm{vx}}_i+{\mathrm{uy}}_i+t_y-y_i^{\′})}^2]$

由上式的偏导数为零，得到相似性变换参数如下：

u=(X·X′)/|X|²，

其中，X·X′表示两个轮廓向量的点积。

$v=\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(x_i{y}_i^{\′}-y_i{x}_i^{\′})\right)/{\left|X\right|}^2$

$t_x=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(x_i^{\′}-{\mathrm{ux}}_i^{\′}+{\mathrm{vy}}_i)$

$t_y=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(y_i^{\′}-{\mathrm{vx}}_i-{\mathrm{uy}}_i)$

欲将Ω中所有样本配准，就是要将训练集Ω中所有样本对齐到同一框架下，具体配准步骤如下：

1）将所有样本中心点坐标平移到原点，即将每个样本中的标记点（x，y）坐标分别减去x的均值和y的均值；

2）对训练集Ω中的L个样本，第一次任选一个样本作为初始平均形状，这里我们选择第一个样本X₁作为初始平均形状，即其余样本经相似性变换对齐到得到对齐后样本集Ω={X₁,X₂,…X_L}（这里用X₁,X₂,…X_L表示变换后的各样本集）；

3）计算变换后平均形状，记为

4）再次将所有样本对齐到平均形状

5）重新计算对齐后所有样本的平均形状，记为

6）计算步骤3平均形状和步骤5平均形状差值

7）若小于某一阈值，则表明平均形状已经收敛，算法结束，配准完成；否则设返回到步骤4）继续进行配准。本发明实施例中，配准前（图中左边部分）和配准后（图中右边部分）如图10所示。

S203，建立点分布模型。

对训练样本集Ω进行相似性变换后得到的轮廓点仍用X_i(i=1,2,…,L)表示，训练图像中目标物体的轮廓在外形和位置上就趋向归一化。对Ω中向量所包含的边界点在2n维空间的分布进行统计分析，构造主动形状模型(Active ShapeModel，ASM)的先验模型，也叫作点分布模型(Point Distripution Model，PDM)。

点分布模型通常是一个参数化的模型，形如X=M(b)，具体构造过程如下:

1）计算数据平均值

$\stackrel{\‾}{X}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{X}_i$

2）计算数据协方差S：

$S=\frac{1}{L-1}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X}){(X_i-\stackrel{\‾}{X})}^T$

3）计算S的特征值λ_i和相应的特征向量Φ_i：

SΦ_i=λ_iΦ_i

4）将得到的特征值按降序排列，即:

λ_i≥λ_i+1(i=1,2,…,2n-1)，其中λ_i对应特征向量Φ_i构成特征向量集Φ=[Φ₁,Φ₂,…,Φ_2n]。

按照降序取前k个特征值及对应的特征向量，其选取依据为:使前k个特征值所决定的目标物体形变占所有2n个特征值所决定目标物体形变总量的比例不小于V(一般V取0.98)，即:

$\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\≥V\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i$

最终，我们得到肺部点分布模型，如下：

$X=\stackrel{\‾}{X}+\mathrm{\Φb}$

其中Φ是由S的前k个特征值对应的特征向量构成的2n*k维矩阵，Φ=(Φ₁,Φ₂,…,Φ_k)，由于k个特征向量相互正交，即ΦΦ^T=I，因此，b可以表示为k维变量: $b={\mathrm{\Φ}}^T(X-\stackrel{\‾}{X}),$ 其中 $-3\sqrt{{\mathrm{\λ}}_i}\≤b_i\≤3\sqrt{{\mathrm{\λ}}_i}.$

S204，建立加权灰度局部纹理模型和邻近距离约束模型。

首先是建立加权灰度局部纹理模型。

对训练样本的每一个边界点沿轮廓法线方向（此处法线方向定义为标记点左右两个邻近标记点连线的垂线方向）进行灰度采样，假设对于第i个边界点(x_i,y_i)，在其两侧沿法线方向各取k个点，用此2k+1个点的灰度值构成向量g_i。对每一个训练样本中同一位置的边界点都经过这样处理，得到一个向量集{g_i}。一般而言，{g_i}满足高斯分布,{g_i}的均值记为协方差记为对于每一个边界点，用同样的方法进行灰度采样与分析，从而构成目标轮廓的灰度统计模型。

对于第j幅图像的第i个标记点采样为：

$g_{i_j}={[g_{i_{j}0},g_{i_{j}1},...,g_{i_{j}2k}]}^T$

对灰度值进行差分，表示如下：

$g_{i_j}={\mathrm{dg}}_{i_j}={[g_{i_{j}1}-g_{i_{j}0},g_{i_{j}2}-g_{i_{j}1},...,g_{i_{j}2k}-g_{i_{j}2k-1}]}^T$

对于从标记点法线方向两边采样得到的2k+1个采样点(中心点为标记点)，我们认为靠近中心点的采样点影响较大，而远离中心点的采样点影响较小，故而对于法线两侧靠近中心点的各p个采样点给予一个较大的权值β，而法线两侧远离中心点各q=k-p个采样点一个较小的权值α。

$g_{i_j}={[{{\mathrm{\αg}}_{i_{j}1}}^{\′},...,{{\mathrm{\αg}}_{i_{j}q}}^{\′},{{\mathrm{\βg}}_{i_j(q+1)}}^{\′},...,{{\mathrm{\βg}}_{i_j(q+2p)}}^{\′},{{\mathrm{\αg}}_{i_j(q+2p+1)}}^{\′},...,{{\mathrm{\αg}}_{i_j(2q+2p)}}^{\′}]}^T$

其中β取值0.65-0.75，α=1-β。

归一化：

${g_{i_j}}^{\′}=\frac{g_{i_j}}{\underset{k=1}{\overset{2k}{\mathrm{\Σ}}}\left|g_{i_{j}k}\right|}$

对于L幅图像得到的向量集{g_i}，每幅图像标记点数目为n（实际中右肺标记点数目为44，左肺标记点数目为50），如下计算每个特征点的均值协方差

${\stackrel{\‾}{g}}_i=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{g}_{i_j}^{\′},i=\mathrm{1,2},...,n$

$S_{g_i}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}(g_{i_j}^{\′}-{\stackrel{\‾}{g}}_i){(g_{i_j}^{\′}-{\stackrel{\‾}{g}}_i)}^T,i=\mathrm{1,2},...n.$

其次是建立邻近距离模型。

对第j幅图像的第i个标记点，计算它与同一幅图像的相邻的两个标记点的距离对于所有样本图像对应标记点计算其距离均值为：

$d_i=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{d}_{i_j},i=\mathrm{1,2},...n.$

从上述本发明实施例提供的X光胸片肺部分割方法可知，由于可以通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域，在所述矩形区域内得到初始化肺部形状，通过迭代方式调整肺部图像时将形状参数b的调整约束在主动形状模型提供的形状参数的范围之内。与现有技术提供的X光胸片肺部分割方法相比，本发明实施例提供的方法能够获得较好的初始化肺部形状，在后续的调整过程中不会造成过分割现象（例如，将脊柱和胃泡识别为肺区），在主动形状模型的约束下能够对X光胸片上的肺部区域进行准确分割，从而为临床或者计算机分析提供有价值的数据。

下面对用于执行上述X光胸片肺部分割方法的本发明实施例的X光胸片肺部分割装置进行说明，其基本逻辑结构参考图11，主要包括获取模块1101、初始化模块1102、最佳匹配点找寻模块1103和形状调整模块1104，各模块详细说明如下：

获取模块1101，用于通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域；

初始化模块1102，用于在所述两个矩形区域内初始化肺部以获得肺部的初始形状；

最佳匹配点找寻模块1103，用于根据加权灰度局部纹理模型，寻找肺部图像中各个特征点（x_i，y_i）的最佳匹配点（x’_i，y’_i），记dI_X=（d x₁，d y₁，d x₂，d y₂，…，d x_i，d y_i，…，d x_n，d y_n）^T，所述d x_i=x’_i－x_i，所述d y_i=y’_i－y_i，所述i=1，2，…，n，所述肺部图像包括调整阶段的肺部形状和所述肺部的初始形状；

形状调整模块1104，用于通过调整姿态参数和形状参数b，使得肺部当前形状I_Xc最大限度地接近I_X+dI_X，所述形状参数b的调整结果为调整后的值b+db受主动形状模型提供的形状参数的约束，所述I_X是与主动形状模型空间中肺部形状X对应的图像空间中实际形状；

所述经过形状调整模块1104调整姿态参数和形状参数b后的肺部形状作为最佳匹配点找寻模块1103的执行对象，重复最佳匹配点找寻模块1103和形状调整模块1104的执行过程，直至所获肺部形状与紧邻本次调整前一次调整所获肺部形状相比，两者的改变量小于预定阈值，迭代搜索算法停止即迭代收敛。

需要说明的是，以上附图11示例的X光胸片肺部分割装置的实施方式中，各功能模块的划分仅是举例说明，实际应用中可以根据需要，例如相应硬件的配置要求或者软件的实现的便利考虑，而将上述功能分配由不同的功能模块完成，即将所述X光胸片肺部分割装置的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。而且，实际应用中，本实施例中的相应的功能模块可以是由相应的硬件实现，也可以由相应的硬件执行相应的软件完成，例如，前述的获取模块，可以是具有执行前述通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域的硬件，例如获取器，也可以是能够执行相应计算机程序从而完成前述功能的一般处理器或者其他硬件设备；再如前述的初始化模块，可以是具有执行前述在所述两个矩形区域内初始化肺部以获得肺部的初始形状功能的硬件，例如初始化器，也可以是能够执行相应计算机程序从而完成前述功能的一般处理器或者其他硬件设备（本说明书提供的各个实施例都可应用上述描述原则）。

附图11示例的获取模块1101可以包括投影单元1201和矩形区域形成单元1202，如附图12所示本发明另一实施例提供的X光胸片肺部分割装置，其中：

投影单元1201，用于通过将所述X光胸片中左肺部图像和右肺部图像进行水平和垂直投影，获得投影曲线；

矩形区域形成单元1202，用于在所述投影曲线中找到各个肺部的左边界和右边界以及其中任一肺部的顶部边界和底部边界对应的坐标点，以所述左肺部的左边界、右边界、顶部边界和底部边界对应的坐标点围成的矩形区域作为X光胸片中左肺部图像的矩形区域，以所述右肺部的左边界、右边界、顶部边界和底部边界对应的坐标点围成的矩形区域作为X光胸片中右肺部图像的矩形区域。

附图11示例的初始化模块1102可以包括转化单元1301、寻找单元1302和变换单元1303，如附图13所示本发明另一实施例提供的X光胸片肺部分割装置，其中：

转化单元1301，用于将主动形状模型特征空间中平均形状从一维形状向量转化为二维空间中n个标记点坐标（x₁，y₁），（x₂，y₂），…，（x_n，y_n），所述 $\stackrel{\‾}{X}={(x_1,y_1,x_2,y_2,...,x_n,y_n)}^T;$

寻找单元1302，用于找到包围该标记点坐标的最小矩形区域；

变换单元1303，用于根据所述矩形区域和最小矩形区域，确定变换参数(u₀,v₀,t_x0,t_y0)，将所述n个标记点坐标变换到所述矩形区域所属图像坐标系下。

附图11示例的最佳匹配点找寻模块1103可以包括最佳匹配点获取单元1401，如附图14所示本发明另一实施例提供的X光胸片肺部分割装置。最佳匹配点获取单元1401用于根据加权灰度局部纹理模型，获取使得代价函数f(g_i)取最小值的标记点作为标记点（x_i，y_i）的最佳匹配点（x’_i，y’_i），所述 $f\left(g_i\right)=100*{(g_i-\stackrel{\‾}{g})}^T{S}_g^{-1}(g_i-\stackrel{\‾}{g})-{\mathrm{\λ}}_1*l\left(i\right)+{\mathrm{\λ}}_{2}d\left(i\right),$ 所述λ₁∈[0.7，1.3]，所述λ₁∈[5，10]，所述l(i)=g(x_i,y_i)，所述g(x_i,y_i)为对肺部图像进行－45°角线性滤波后的肺部图像在标记点处的灰度值，只取左肺横隔膜部分，所述 $d\left(i\right)=\frac{1}{2}(\sqrt{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}+\sqrt{{(x_i-x_{i+1})}^2+{(y_i-y_{i+1})}^2}),$ 所述i=21，…，29，所述g_i为加权灰度局部纹理模型提供的训练样本第i个边界点（x_i，y_i）两侧沿法线方向各取k个点连同所述第i个边界点（x_i，y_i）共2k+1个点的灰度值构成的向量。

附图11至附图14任一示例的X光胸片肺部分割装置还可以包括建模模块1501，如附图15-a至附图15-d任一示例所示本发明另一实施例提供的X光胸片肺部分割装置。建模模块1501用于建立所述加权灰度局部纹理模型和所述主动形状模型。建模模块1501建立上述模型的过程可参阅前述实施例，这里不再赘述。

需要说明的是，上述装置各模块/单元之间的信息交互、执行过程等内容，由于与本发明方法实施例基于同一构思，其带来的技术效果与本发明方法实施例相同，具体内容可参见本发明方法实施例中的叙述，此处不再赘述。

本领域普通技术人员可以理解上述实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器（ROM，Read Only Memory）、随机存取存储器（RAM，Random Access Memory）、磁盘或光盘等。

以上对本发明实施例所提供的一种X光胸片肺部分割方法和装置进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种X光胸片肺部分割方法，其特征在于，所述方法包括如下步骤：

S101，通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域；

S102，在所述两个矩形区域内初始化肺部以获得肺部的初始形状；

S103，根据加权灰度局部纹理模型，寻找肺部图像中各个特征点（x_i，y_i）的最佳匹配点（x’_i，y’_i），记dI_X=（d x₁，d y₁，d x₂，d y₂，…，d x_i，d y_i，…，d x_n，d y_n）^T，所述d x_i=x’_i－x_i，所述d y_i=y’_i－y_i，所述i=1，2，…，n，所述肺部图像包括调整阶段的肺部形状和所述肺部的初始形状；

S104，通过调整姿态参数和形状参数b，使得肺部当前形状I_Xc最大限度地接近I_X+dI_X，所述形状参数b的调整结果为调整后的值b+db受主动形状模型提供的形状参数的约束，所述I_X是与主动形状模型空间中肺部形状X对应的图像空间中实际形状；

将所述经过调整姿态参数和形状参数b后的肺部形状作为步骤S103的执行对象，重复步骤S103和步骤S104，直至所获肺部形状与紧邻本次调整前一次调整所获肺部形状相比，两者的改变量小于预定阈值。

2.根据权利要求1所述的方法，其特征在于，所述通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域包括：

通过将所述X光胸片中左肺部图像和右肺部图像进行水平和垂直投影，获得投影曲线；

在所述投影曲线中找到各个肺部的左边界和右边界以及其中任一肺部的顶部边界和底部边界对应的坐标点，以所述左肺部的左边界、右边界、顶部边界和底部边界对应的坐标点围成的矩形区域作为X光胸片中左肺部图像的矩形区域，以所述右肺部的左边界、右边界、顶部边界和底部边界对应的坐标点围成的矩形区域作为X光胸片中右肺部图像的矩形区域。

3.根据权利要求1所述的方法，其特征在于，所述在所述两个矩形区域内初始化肺部以获得肺部的初始形状包括：

将主动形状模型特征空间中平均形状从一维形状向量转化为二维空间中n个标记点坐标（x₁，y₁），（x₂，y₂），…，（x_n，y_n），所述

找到包围该标记点坐标的最小矩形区域；

根据所述矩形区域和最小矩形区域，确定变换参数(u₀,v₀,t_x0,t_y0)，将所述n个标记点坐标变换到所述矩形区域所属图像坐标系下。

4.根据权利要求1所述的方法，其特征在于，所述根据加权灰度局部纹理模型，寻找肺部图像中各个标记点（x_i，y_i）的最佳匹配点（x’_i，y’_i）包括：

根据加权灰度局部纹理模型，获取使得代价函数f(g_i)取最小值的标记点作为标记点（x_i，y_i）的最佳匹配点（x’_i，y’_i），所述 $f\left(g_i\right)=100*{(g_i-\stackrel{\‾}{g})}^T{S}_g^{-1}(g_i-\stackrel{\‾}{g})-{\mathrm{\λ}}_1*l\left(i\right)+{\mathrm{\λ}}_{2}d\left(i\right),$ 所述λ₁∈[0.7，1.3]，所述λ₁∈[5，10]，所述l(i)=g(x_i,y_i)，所述g(x_i,y_i)为对肺部图像进行－45°角线性滤波后的肺部图像在标记点处的灰度值，只取左肺横隔膜部分，所述 $d\left(i\right)=\frac{1}{2}(\sqrt{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}+\sqrt{{(x_i-x_{i+1})}^2+{(y_i-y_{i+1})}^2}),$ 所述i=21，…，29，所述g_i为加权灰度局部纹理模型提供的训练样本第i个边界点（x_i，y_i）两侧沿法线方向各取k个点连同所述第i个边界点（x_i，y_i）共2k+1个点的灰度值构成的向量。

5.根据权利要求1至4任意一项所述的方法，其特征在于，所述方法还包括：

建立所述加权灰度局部纹理模型和所述主动形状模型。

6.一种X光胸片肺部分割装置，其特征在于，所述装置包括：

获取模块，用于通过水平垂直投影，获取分别包围X光胸片中左肺部图像和右肺部图像的两个矩形区域；

初始化模块，用于在所述两个矩形区域内初始化肺部以获得肺部的初始形状；

最佳匹配点找寻模块，用于根据加权灰度局部纹理模型，寻找肺部图像中各个特征点（x_i，y_i）的最佳匹配点（x’_i，y’_i），记dI_X=（d x₁，d y₁，d x₂，d y₂，…，d x_i，d y_i，…，d x_n，d y_n）^T，所述d x_i=x’_i－x_i，所述d y_i=y’_i－y_i，所述i=1，2，…，n，所述肺部图像包括调整阶段的肺部形状和所述肺部的初始形状；

形状调整模块，用于通过调整姿态参数和形状参数b，使得肺部当前形状I_Xc最大限度地接近I_X+dI_X，所述形状参数b的调整结果为调整后的值b+db受主动形状模型提供的形状参数的约束，所述I_X是与主动形状模型空间中肺部形状X对应的图像空间中实际形状；

所述经过形状调整模块调整姿态参数和形状参数b后的肺部形状作为最佳匹配点找寻模块的执行对象，重复最佳匹配点找寻模块和形状调整模块的执行过程，直至所获肺部形状与紧邻本次调整前一次调整所获肺部形状相比，两者的改变量小于预定阈值。

7.根据权利要求6所述的装置，其特征在于，所述获取模块包括：

投影单元，用于通过将所述X光胸片中左肺部图像和右肺部图像进行水平和垂直投影，获得投影曲线；

矩形区域形成单元，用于在所述投影曲线中找到各个肺部的左边界和右边界以及其中任一肺部的顶部边界和底部边界对应的坐标点，以所述左肺部的左边界、右边界、顶部边界和底部边界对应的坐标点围成的矩形区域作为X光胸片中左肺部图像的矩形区域，以所述右肺部的左边界、右边界、顶部边界和底部边界对应的坐标点围成的矩形区域作为X光胸片中右肺部图像的矩形区域。

8.根据权利要求6所述的装置，其特征在于，所述初始化模块包括：

转化单元，用于将主动形状模型特征空间中平均形状从一维形状向量转化为二维空间中n个标记点坐标（x₁，y₁），（x₂，y₂），…，（x_n，y_n），所述 $\stackrel{\‾}{X}={(x_1,y_1,x_2,y_2,...,x_n,y_n)}^T;$

寻找单元，用于找到包围该标记点坐标的最小矩形区域；

变换单元，用于根据所述矩形区域和最小矩形区域，确定变换参数(u₀,v₀,t_x0,t_y0)，将所述n个标记点坐标变换到所述矩形区域所属图像坐标系下。

9.根据权利要求6所述的装置，其特征在于，所述最佳匹配点找寻模块包括：

最佳匹配点获取单元，用于根据加权灰度局部纹理模型，获取使得代价函数f(g_i)取最小值的标记点作为标记点（x_i，y_i）的最佳匹配点（x’_i，y’_i），所述 $f\left(g_i\right)=100*{(g_i-\stackrel{\‾}{g})}^T{S}_g^{-1}(g_i-\stackrel{\‾}{g})-{\mathrm{\λ}}_1*l\left(i\right)+{\mathrm{\λ}}_{2}d\left(i\right),$ 所述λ₁∈[0.7，1.3]，所述λ₁∈[5，10]，所述l(i)=g(x_i,y_i)，所述g(x_i,y_i)为对肺部图像进行－45°角线性滤波后的肺部图像在标记点处的灰度值，只取左肺横隔膜部分，所述 $d\left(i\right)=\frac{1}{2}(\sqrt{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2}+\sqrt{{(x_i-x_{i+1})}^2+{(y_i-y_{i+1})}^2}),$ 所述i=21，…，29，所述g_i为加权灰度局部纹理模型提供的训练样本第i个边界点（x_i，y_i）两侧沿法线方向各取k个点连同所述第i个边界点（x_i，y_i）共2k+1个点的灰度值构成的向量。

10.根据权利要求6至9任意一项所述的装置，其特征在于，所述装置还包括：

建模模块，用于建立所述加权灰度局部纹理模型和所述主动形状模型。