一种医学图像中肺分割的方法及装置与流程

本发明涉及医学图像分割领域，尤其是在医学图像中对肺分割的方法及装置。

背景技术：

肺分割是数字胸片图像处理中的一个重要环节。分割结果的好坏直接影响后续对病灶的检测及分析。胸片最大的缺陷在于组织重叠。人体结构中组织的密度越高，吸收X射线的能力就越强，因此胸部的肋骨吸收X射线就较多，在图像中呈白影；由于肺部含气体，所以密度较低，吸收X射线较少，因此图像上呈黑影。正因如此，在分割肺边界时，容易出现锯齿状的分割结果。而且通过观察图像可以发现，肺的下边界处存在角点，这里也是分割的难点。

常用的肺分割方法大体可分为两类：一类是基于规则级的肺分割，主要包括阈值分割法、区域生长法、边缘检测法、形态学滤波法等。由于X光图像的成像效果的原因，这类方法不能对肺区域进行精细分割；另一类是基于像素分类的肺分割方法，主要有遗传算法、神经网络、模糊聚类等学习方法，对图像中的像素进行分类。

还有一种是基于综合知识的分割方法：主动形状模型法(ASM)。这种方法首先需要人工标定肺的轮廓作为参考标准，通过训练图像样本获取特征点分布的统计信息，并且获取特征点允许存在的变化方向，实现在目标图像上寻找对应的特征点的位置。但是这种方法是在二维信息上进行查找边界，很容易受到 肋骨边界的影响，出现锯齿状的分割结果。

技术实现要素：

本发明所要解决的技术问题是提供一种能解决或降低前述问题的肺分割的方法及装置。

本发明为解决上述技术问题而采用的技术方案是：一种医学图像中肺分割的方法，其特征在于包括以下步骤：

根据训练集中的M张胸片的肺轮廓，得到平均肺模板，其中M为大于或等于2的整数；

对拟被分割的肺部图像预处理，获得预处理后的肺部图像；

提取预处理后的肺部图像的肺边界的二值图像，并根据广义霍夫变换进行初始定位，获得相对应的霍夫定位位置；

将平均肺模板与霍夫定位位置对齐，获得对齐结果；

利用动态规划算法进行分割，并将分割结果反变换回原始坐标系中，完成肺区域分割。

进一步的，在M张胸片中所对应的每张胸片的肺轮廓上标记相同数量的特征点，且不同胸片上的特征点的分布一致。

进一步的，所述在每张胸片的肺轮廓的拐点处、最高点处及最低点处标记特征点。

进一步的，所述在每张胸片的肺轮廓上除拐点处、最高点处及最低点处之外的其他位置等距离标记特征点。

进一步的，所述平均肺模板按照以下步骤获得：

a.从训练集中的M张胸片中进行第一次选取，获取第一肺区形状L₀；

b.从训练集中的尚未被选取的胸片中进行第二次选取，获取N张肺区形状 L₁、L₂…L_N；

c.使所述N张肺区形状L₁、L₂…L_N经处理后得到相应的N张肺区形状L₁₁、L₂₁…L_N1，分别与第一肺区形状L₀对齐，计算得到肺区平均形状L_ave，所述处理包括使N张肺区形状L₁、L₂…L_N分别相对于第一肺区形状L₀作旋转、缩放、平移变换；

d.若c步骤的肺区平均形状L_ave相较于第一肺区形状L₀的旋转、缩放、平移小于设定阈值，则将所述肺区平均形状L_ave作为平均肺模板；否则，进入e步骤；

e.将第一肺区形状L₀进行处理、所述N张肺区形状L₁₁、L₂₁…L_N1进行处理后，分别与当前肺区平均形状L_ave-i对齐，计算得到本次肺区平均形状L_ave1-i，所述对第一肺区形状处理包括：使第一肺区的形状L₀分别相对于当前肺区平均形状L_ave-i作旋转、缩放、平移变换，得到相应的肺区形状L_0i，以及使所述N张肺区形状L₁₁、L₂₁…L_N1分别相对于当前肺区平均形状L_ave-i作旋转、缩放、平移变换，得到N张肺区形状L_1i、L_2i…L_Ni；

f.若e步骤计算所得本次肺区平均形状L_ave1-i与当前肺区平均形状L_ave-i的旋转、缩放、平移小于设定阈值，则将e步骤计算所得本次肺区平均形状L_ave1-i作为平均肺模板；否则，将当前肺区平均形状L_ave-i用本次肺区平均形状L_ave1-i替代，将第一肺区的形状L₀用肺区形状L_0i替代，N张第二肺区形状L₁₁、L₂₁…L_N1用N张肺区形状L_1i、L_2i…L_Ni替代，且i＝i+1，其中i大于或等于2的整数，“＝”为赋值号，并返回步骤e。

进一步的，所述肺部图像预处理包括以下方式：获取原始图像，对原始图像做高斯滤波处理后得到滤波图像，对原始图像与滤波图像做差，然后施加滤波图像的灰度平均值，得到去背景图像，然后对去背景图像做双边滤波操作， 即得到预处理后的图像。

进一步的，基于预处理后的图像，建立与图像边缘轮廓相对应的若干个边缘形状模板；获取每个边缘形状模板的重心点的坐标以及边缘形状上的点相对于重心点的坐标位置，应用霍夫变换法对边缘形状的二值图像进行初始定位。

进一步的，应用边界增强算子提取边缘形状的梯度图像，再保留最大的15％的像素点生成边缘形状的二值图像。

进一步的，使用边界二值图像上每一个不为零的点去匹配模板中的任一点，在霍夫空间中相对应的重心点的位置处加1；结束二值图像中所有等于1的点的霍夫变换后，求出霍夫空间中每一点的累加和，并找到最大值的位置，该点就是使用当前模板定位的最准确的模板的重心位置；对每个模板的重心点的最大值做进行归一化，然后确定最大值对应的模板位置即为定位位置。

本发明为解决上述技术问题而采用的技术方案是：一种医学图像中肺分割的装置，其特征在于包括：

存储单元，存储有训练集的数张胸片图像及拟被分割胸片图像；

显示单元，显示对应的胸片图像；

输入单元，对胸片中的肺部标记特征点，获得相应的肺轮廓图像；

图像处理单元，其用于：

对肺轮廓图像进行处理后获得平均肺模板；

对拟被分割的胸部图像预处理，获得预处理后的肺部图像；

提取预处理后的肺部图像的肺边界的二值图像，并根据广义霍夫变换进行初始定位，获得相对应的霍夫定位位置；

将平均肺模板与霍夫定位位置对齐，获得对齐结果；

利用动态规划算法进行分割，并将分割结果反变换回原始坐标系中，完 成肺区域分割。

本发明对比现有技术有如下的有益效果：本发明通过霍夫变换进行拟被分割的肺区进行定位，并通过一个平均肺区形状与定位位置对齐，得到肺部的初始轮廓，最大程度的缩小了初始位置与肺部实际位置的偏差，从而提高了分割的准确性；另外，本发明是一种自动化的肺分割的方法，工作效率高。

【附图说明】

图1为本发明实施例的肺分割的方法示意图；

图2为本发明实施例的全自动肺分割的具体流程图；

图3(a)为本发明实施例中一获取平均肺模板方法的流程图；

图3(b)为本发明实施例中另一获取平均肺模板方法的流程图；

图4为本发明实施例中用于肺部上边缘的霍夫模板示意图；

图5为本发明实施例中霍夫定位示意图；

图6为本发明实施例中肺部对齐结果示意图；

图7为本发明实施例中肺部分段分割结果示意图；

图8为本发明实施例中右肺分割结果示意图。

【具体实施方式】

下面结合附图和实施例对本发明作进一步的描述。

请参阅图1-3，本发明实施例的一种医学图像(例如DR图像或X射线图像)中肺分割的方法，包括以下步骤：

根据训练集中的M张胸片的肺轮廓，得到平均肺模板，其中M为大于或 等于2的整数；

获取拟被分割的肺部图像；

对拟被分割的肺部图像预处理，获得预处理后的肺部图像；

提取预处理后的肺部图像的肺边界的二值图像，并根据广义霍夫变换进行初始定位，获得相对应的霍夫定位位置；

将平均肺模板与霍夫定位位置对齐，获得霍夫定位结果；

利用动态规划算法进行分割，并将分割结果反变换回原始坐标系中，完成肺区域分割。

实施例中的肺分割，是通过医学处理装置进行全自动的肺分割。所述训练集中包括若干张胸片，胸片被存储在医学处理装置的存储单元内；所述训练集包含100张胸片，每张胸片对应一对肺轮廓图像。在其他实施例中，训练集中胸片的数量M可以为80张、120张、150张等。经发明人多次实验，认为100张胸片即可较好的满足实际需求；当然，胸片的数量越多，所得到的平均肺模板则越好，但是在图像分割中所需要的时间也越长。

选取M张胸片，在M张胸片中所对应的每张胸片的一个肺部(边缘)上标记相同数量的特征点，且不同胸片上的特征点的分布一致。实施例中，在每张胸片上标记42个特征点，其中在每张胸片的肺轮廓的拐点处、最高点处及最低点处都需要标记特征点。在每张胸片的肺轮廓上除拐点处、最高点处及最低点处之外的其他位置等距离标记特征点。当然，特征点的个数也可以为36个、50个等。标记后，就得到了M个肺的轮廓(形状)。但是手工标定是会存在一定误差，而且这些模型会出现在不同位置，并具有不同大小不同旋转角度。为了消除这些“非形状”因素，求取平均肺模板，我们需要对训练集进行对齐。

具体实施方式中，利用循环计算(比较)的方式获得所述平均肺模板，以下为相应的两种可选步骤：

第一种获取平均肺模板的方式(步骤)：

a.从训练集中的M张胸片中进行第一次选取，获取第一肺区形状L₀；

b.从训练集中的尚未被选取的胸片中进行第二次选取，获取N张(第二)肺区形状L₁、L₂…L_N；

c.使所述N张肺区形状L₁、L₂…L_N经处理后得到相应的N张肺区形状L₁₁、L₂₁…L_N1，分别与第一肺区形状L₀对齐，计算得到肺区平均形状L_ave，所述处理包括使N张肺区形状L₁、L₂…L_N分别相对于第一肺区形状L₀作旋转、缩放、平移变换；

d.若c步骤的肺区平均形状L_ave相较于第一肺区形状L₀的旋转、缩放、平移小于设定阈值，则将所述肺区平均形状L_ave作为平均肺模板；否则，进入e步骤；

e.将第一肺区形状L₀进行处理、所述N张肺区形状L₁₁、L₂₁…L_N1进行处理后，分别与当前肺区平均形状L_ave-i对齐，计算得到本次肺区平均形状L_ave1-i，所述对第一肺区形状处理包括：使第一肺区的形状L₀分别相对于当前肺区平均形状L_ave-i作旋转、缩放、平移变换，得到相应的肺区形状L_0i，以及使所述N张肺区形状L₁₁、L₂₁…L_N1分别相对于当前肺区平均形状L_ave-i作旋转、缩放、平移变换，得到N张肺区形状L_1i、L_2i…L_Ni；

f.若e步骤计算所得本次肺区平均形状L_ave1-i与当前肺区平均形状L_ave-i的旋转、缩放、平移小于设定阈值，则将e步骤计算所得本次肺区平均形状L_ave1-i作为平均肺模板；否则，将当前肺区平均形状L_ave-i用本次肺区平均形状L_ave1-i 替代，将第一肺区的形状L₀用肺区形状L_0i替代，N张第二肺区形状L₁₁、L₂₁…L_N1用N张肺区形状L_1i、L_2i…L_Ni替代，且i-i+1，其中i大于或等于2的整数，“＝”为赋值号，并返回步骤e。

进一步的，若e步骤计算所得本次肺区平均形状L_ave1-i与当前肺区平均形状L_ave-i的旋转、缩放、平移大于设定阈值，但是当i大于或等于90时，停止迭代计算，并以本次肺区平均形状L_ave1-i作为平均肺模板。

以上N为小于或等于M-1的自然数，i为大于或等于1的自然数。

所述阈值为使两个肺区形状(图像)对齐时，其中一个肺区的形状相对于另一个肺区的形状所做的几何变换的值：平移距离小于或等于0.01像素，并且旋转角度小于或等于0.001*pi/180，并且缩放尺度小于或等于0.001。

第二种获取平均肺模板的方式(步骤)：

a.从训练集中的M张胸片中做第一次选取，选取一张胸片，获取第一肺区形状；

b.从训练集中未被选定的胸片中做第二次选取，选取至少一张胸片，获取对应的(第二)肺区形状；

c.通过旋转、缩放、平移使所述第二次选取所得的肺区形状分别与第一肺区形状对齐，计算得到当前肺区平均形状；

d.若c步骤的当前肺区平均形状相较于第一肺区形状的旋转、缩放、平移小于设定阈值，则将肺区当前平均形状作为平均肺模板；否则，进入e步骤；

e.从训练集中未被选定的胸片中做第N次选取，选取至少一张胸片，获得对应的肺区形状，通过旋转、缩放、平移将所述第N次选取所得的肺区形状分别与当前肺区平均形状对齐，并计算出至第N次已选取所有肺区形状的平均形 状，其中，N为大于或等于3的整数；

f.若至第N次已选取所有肺区形状的平均形状相较于当前肺区平均形状的旋转、缩放、平移小于设定阈值，则将至第N次已选取所有肺区形状的平均形状作为平均肺模板；否则，将至第N次已选取所有肺区形状的平均形状作作为当前肺区平均形状，今N＝N+1，其中，“＝”为赋值号；并返回步骤e。

进一步的，若训练集中所有胸片被选定完毕，则以最后一次所选胸片所得肺区形状与当前肺区平均形状计算所得肺区形状作为平均肺模板。

所述阈值为使两个肺区形状(图像)对齐时，其中一个肺区的形状相对于另一个肺区的形状所做的几何变换的值：平移距离小于或等于0.01像素，并且旋转角度小于或等于0.001*pi/180，并且缩放尺度小于或等于0.001。

通过以上方法计算出平均肺模板后，将其保存在存储单元中，以便于后续使用。

所述肺部图像预处理包括以下步骤：获取被拍摄部位肺部的原始图像，对原始图像做高斯滤波处理后得到滤波图像，对原始图像与滤波图像做差，然后施加滤波图像的灰度平均值，得到去背景图像，然后对去背景图像做双边滤波操作，即得到预处理后的肺部图像。

对预处理后的肺部图像分别提取其上、下、外边界的二值图像，采用基于广义霍夫变换进行初始定位。所述上、下、外(左、右)侧边界共同组成肺部边缘形状(肺部轮廓)。根据边缘形状，可建立若干个与之较相似的边缘形状模板。具体实施方式中，可建立3至8个边缘形状模板，优选的数量为5个边缘形状模板。具体的，可分别将一个边缘形状划分为上边界、下边界、外侧边界等几部分，并将这几部分分别进行处理。

以上边界定位为例：首先标记5个上边缘形状(上边界)模板(参图4)， 记录每一个模板的重心点坐标以及模板上每一点相对于重心点的坐标位置。应用霍夫变换算法对上边界二值图像进行初始定位。霍夫变换算法，是对图像进行某种形式的坐标变换，使得图像上给定形状曲线上的所有点经过变换以后都集中到霍夫空间的某些位置上形成峰值点，这样就把对原图中给定形状曲线的检测问题转化为寻找霍夫空间中峰值点的问题。

以任一模板为例，应用边界增强算子提取上边界的梯度图像，再保留最大的15％的像素点生成上边界的二值图像。使用上边界的二值图像上每一个不为零的点去匹配模板中的任一点，在霍夫空间中相对应的重心点的位置处加1。在霍夫空间中记录的是可能匹配准确的模板的重心位置，因此在霍夫空间中累积的值越大，那么该点是图像中形状的重心点的可能性就越大。在空间中找到最大值点，我们就确定该点为形状的重心点。结束二值图像中所有等于1的点的霍夫变换后，求出霍夫空间中每一点的累加和，并找到最大值的位置，该点就是使用当前模板定位的最准确的模板的重心位置。重心位置确定后，模板的相应位置也就确定，这就是找到的上边缘位置。

使用5个模板分别对肺区域上边界二值图像进行霍夫变换，在5个霍夫空间中会得到5个最大值，由于模板大小不同，需要对最大值进行归一化。再在这5个最大值中找到最大值，再确定该最大值对应的模板，我们就将这个模板确定为最佳模板。

然后对外侧、下边界也是与上边界同样的操作进行定位，定位结果如图5所示；即可获得与肺部图像相对应的霍夫定位位置。

将平均肺模板与霍夫定位位置对齐，获得对齐结果(对齐图像，参图6所示)。如图7所示，将对齐后的结果(共42个点)中的三部分：外侧、下侧及内侧分别标记点，1-16点(外侧)、17-24点(下侧)、25-42点(内侧)，对 上侧我们使用的是霍夫定位后结果(也可对外侧、下侧也可直接利用霍夫定位后的结果)，将这四部分分别进行线性插值，然后在每一点的法线方向提取m个点像素值，组成法线矩阵，在法线矩阵上进行动态规划分割，再将得到的分割结果反变换回原始的坐标系中，就完成了肺区域的分割。最后将这四部分的分割结果进行联接并对边缘进行平滑操作，就完成了肺的分割(如图8所示)。

本专利申请中动态规划分割的过程：

动态规划算法中的路径选择依赖于当前状态，同时也依赖于之前已选择的状态，其通常用来搜索目标的最优边界。

在相应的霍夫模板上取每点的法线方向的像素，并组成一个新的矩阵，称之为“法线矩阵”(类似于极坐标变换后得到的矩阵)，在生成“法线矩阵”的同时也存储每一点在原始坐标系中的位置坐标，方便进行反变换。

根据对齐得到的位置，将对应的边界进行二阶曲线拟合，这样可以比较容易的计算出曲线上每一点的法线。

在每一点的法线方向，以霍夫模板为基准，向上取m1个像素长，向下取m2个像素长，采集每点的像素值，组成一个m*n的法线矩阵，n为霍夫模板的长度，m＝m1+m2。对这个法线矩阵进行动态规划分割。具体包括以下几个主要步骤：

(1)获取局部代价

在动态规划算法中，局部代价由内部代价和外部代价组成。内部代价用来衡量边界点的平滑程度，边界越光滑，内部代价越小。外部代价用来衡量图像梯度变化的大小，梯度越大，外部代价越小。

假设“法线矩阵”的大小为m×n，则我们将内部代价定义为：

E_int(i，j)＝|j-k|/(j+k)，j＝1...n，k＝1...n，i＝1...m (1)

其中，j和k分别为“法线矩阵”中第i列和i-1列上边界点的纵坐标，并对内部代价进行归一化。内部代价的大小表明了边界的平滑程度。

外部代价使用法线矩阵的梯度图像的相反数来表示：

E_ext(i，j)＝-G(i，j) (2)

所以总的局部代价是由内部代价与外部代价的加权和表示的：

E(i，j)＝ω_int×E_int(i，j)+ω_ext×E_ext(i，j) (3)

其中ω_int和ω_ext分别表示内部代价与外部代价的权重。

(2)计算累计代价

计算累计代价是一个动态累加的过程，每一列的累计代价都是前一列的累计代价与当前列的局部代价的累加和。

因为第一列不存在内部代价，所以在累计代价的第一列只有外部代价而没有内部代价。之后，每一列都是迭代的计算累计代价。每一列的累加代价都是前一列在一定范围内累计代价与当前点的局部代价和的最小值。K表示第i-1列的搜索范围，使每相邻两列选择的边界点不会有太大的跳动。在计算每一点的累计代价的同时还要记录一下当前点取最小值时的k的取值，这样便于最优路径的需找。

(3)背向搜索最优路径

计算完所有的累计代价后，找到累计代价最后一列中的最小值作为初始点，根据记录的k值向前进行搜索，找到最优路径。再变换回原始坐标系中，就得到了分割结果。

本发明还提供了一种医学图像中肺分割的装置，其特征在于包括：

存储单元，存储有训练集的数张胸片图像及拟被分割胸片图像；

显示单元，显示对应的胸片图像；

输入单元，对胸片中的肺部标记特征点，获得相应的肺轮廓图像；

图像处理单元，其用于：

对肺轮廓图像进行处理后获得平均肺模板；

对拟被分割的胸部图像预处理，获得预处理后的肺部图像；

提取预处理后的肺部图像的肺边界的二值图像，并根据广义霍夫变换进行初始定位，获得相对应的霍夫定位位置；

将平均肺模板与霍夫定位位置对齐，获得对齐结果；

利用动态规划算法进行分割，并将分割结果反变换回原始坐标系中，完成肺区域分割。

利用上文提出的方法和/或装置，分别使用DR设备提供的23个数据、JSRT数据库的247个数据进行测试。测试结果表明，霍夫变换对肺区域的上、下、外边界的自动定位都比较准确的，而且图像质量较好，分割结果也很理想。

JSRT数据库的图像质量较差，在定位边界及分割时都会有些偏差，但影响不大，结果可以接受。

将所有的分割结果分为4个等级：“很好”，“好”，“可接受”，“差”。共270例数据，其中大约95％的图像分割好或很好。

虽然本发明已以较佳实施例揭示如上，然其并非用以限定本发明，任何本领域技术人员，在不脱离本发明的精神和范围内，当可作些许的修改和完善，因此本发明的保护范围当以权利要求书所界定的为准。