CN103310449A - 基于改进形状模型的肺分割方法 - Google Patents

Abstract

基于改进形状模型的肺分割方法，本发明先建立肺轮廓的先验模型，再利用灰度与形状相似性信息结合图像特征对肺区域分割。由于在某些图像中，初始位置可能与实际边界相距过远，利用灰度与形状相似性信息分割时，搜索区域不覆盖肺边界。因此，本发明通过使用主动形状模型（Active Shape Model,ASM)算法在第一次分割基础上修正肺边界，改善部分点搜索陷入局部极值的情况，获得更优的搜索结果。

Description

基于改进形状模型的肺分割方法

技术领域：本发明涉及一种肺分割方法，尤其是一种基于改进形状模型的肺分割方法。 

背景技术:肺癌是当今对人类健康危害最大的恶性肿瘤之一。特别是近半个世纪以来，随着空气污染造成的环境不断恶化及吸烟人口数量的大量增加，各国肺癌的发病率和病死率都在急剧上升。2012年全美肺癌的死亡病例数为160,259，位居所有癌症死亡率的首位，远远高于其它癌症的死亡率。2008年我国肺癌死亡病例数为493,348，与2004-2005年的30.83万相比，死亡率提高近60%。中国癌症预防与控制规划纲要(2004～2010)确定肺癌为我国癌症防治的重点。 

由于肺属于人体内部脏器，多数肺癌在开始的时候只是在身体内悄悄生长，患者没有任何感觉。当患者因咳嗽、咯血及胸痛等临床症状就诊时，大多数已经处于中晚期，错过了治疗的最佳时机。研究表明，肺癌的早期诊断能降低肺癌患者的死亡率，并且早期肺癌患者术后10年的生存率高于90%。所以，肺癌的早期诊断与治疗是提高肺癌患者生存率的关键。 

近年来，数字化X线摄影（Digital Radiography）使传统的X线摄影技术进入数字化领域。它与传统的X线摄影相比，具有更高的影像质量，包含更多的影像信息等特点。同时，由于DR图像具有成像设备简单，成本低，与CT相比放射线剂量低等特点，所以是肺癌筛查的主要成像方式。由于肺区域的准确分割是肺结节自动检测与辅助诊断的基础，因此是学者 们研究的热点之一。Xu等人提出基于图像特征分析的肺区域分割算法。先通过轮廓的二阶导数确定肺顶部及胸腔边界，再根据边界梯度确定纵隔右边界，由规则确定起始点搜索纵隔左边界，该算法实现对肺区域粗略分割，精度较低。Ginneken等人使用主动形状模型，主动外观模型，及多分辨率像素分类算法对肺区域进行分割，交替使用形状与灰度信息，容易陷入局部极值。Shi等人提出基于普通人群与患者特有肺部形状统计信息的形变模型的肺区域分割算法。Soleymanpour等人先使用基于自适应对比度均衡的非线性滤波器对原始图像增强，再利用区域增长算法对形态学操作获得初始肺区域进行修补获得更加精确的分割结果。刘炎等人提出了基于柔性形态学与聚类算法的肺分割算法，该方法利用柔性腐蚀算子对图像暗区的扩张作用，但分割性能过于依赖第一次分割。石争浩等人对模糊C均值聚类算法进行改进，采用基于核的距离范式替代欧式聚类来分割肺区域。 

发明内容：针对上述现有技术的不足，本发明提供了一种基于改进形状模型的肺分割方法。 

为实现上述目的，本发明采用的技术方案是：基于改进形状模型的肺分割方法，具体步骤为： 

一、模型初始轮廓位置的确定：包括标记训练集，对齐训练集，先验模型的建立三个主要部分。 

①标记训练集：使用沿着边界的点来标记训练集，点包括如下三类：a、标记目标有特定应用部分的点。例如，在人脸模型中，表示眼睛中心的点。表示边界上比较锐利的拐角，如眼角的点；b、标记无关应用部分的点，如目标在特定方向上的最高点，或曲率的极值处；c、填在一类点或二类点之 间的点。 

②对齐训练集：为了能够比较不同形状上的相同的点，它们必须关于一组坐标轴对齐。通过缩放、旋转和平移的操作使训练形状对齐，使它们尽可能对齐紧密。设x_i是训练集中第i个形状中n个点的向量，x_i＝(x_i0,y_i0,x_i1,y_i1,...,x_ik,y_ik,...,x_in-1,y_in-1)^T。其中，(x_ij,y_ij)是第i个形状中第j个点。 

给定两个相似形状x_i和x_j，选择旋转角度θ、缩放s、平移(t_x,t_y)，则M(s,θ)[x]代表旋转角度为θ和缩放比例为s的变换，把x_i映射为M(s,θ)[x_j]+t，使以下加权和最小化： 

E_j＝(x_i-M(s,θ)[x_j]-t)^TW(x_i-M(s,θ)[x_j]-t)    (1) 

其中， 

$M(s,\mathrm{\θ})\left[\begin{array}{c}{x}_{\mathrm{jk}}\\ y_{\mathrm{jk}}\end{array}\right]=\left(\begin{array}{c}\left(s\cos \mathrm{\θ}\right)x_{\mathrm{jk}}-\left(s\sin \mathrm{\θ}\right)y_{\mathrm{jk}}\\ \left(s\sin \mathrm{\θ}\right)x_{\mathrm{jk}}+\left(s\cos \mathrm{\θ}\right)y_{\mathrm{jk}}\end{array}\right),$ t＝(t_x,t_y,...,t_x,t_y)^T。 

③先验模型的建立：把训练集中的形状向量对齐处理之后，就可以利用主成分分析的方法来找出形状变化的统计信息及规律，这样做可以提高算法的效率。 

设平均形状为

对齐后各个样本相对于平均形状的偏差构成的形状向量协方差为

计算该协方差的特征值与特征向量，并将特征值排序Sp_k＝λ_kp_k。其中，λ_k表示第k大的特征值。 

λ_k越大，它所对应的p_k描述数据点变化模式就越重要。选择前t个重要的变化模式组成新的主轴系p_s，则允许形状域内的任何一个形状都可以由平均形状加主轴系与一组控制参数的加权和来近似，

其中，p_s＝(p₁p₂...p_t)是前t个特征向量组成的矩阵，b_s＝(b₁b₂...b_t)^T是权值向量。 

经过主分量分析，按照降序取前t个特征值及对应的特征向量，使前t个特征值所决定的目标物体形变占所有2n个特征值所决定目标物体形变总量的比例不小于V(一般V取0.98)。 

二、基于灰度信息和形状信息的肺实质分割：本发明同时利用图像的灰度与形状信息，使得搜索到的边界灰度、形状信息与训练图像相似。 

（1）特征图像：由于基于导数的特征图像能够更加突出灰度的变化，本文利用特征图像来获取边界点的候选点及各个候选点的灰度代价。 

本发明采用6种特征图像：（1-2）x,y方向一阶偏导数图像，表示x,y方向灰度变化；（3-4）x,y方向二阶偏导数图像，表示x,y方向灰度变化速率；（5）x，y方向混合偏导数图像；（6）x，y方向二阶偏导数和图像，该值越大，表明此处灰度变化速度越快，为肺边界的可能性越大。 

（2）边界点的候选点：对于初始肺边界的每一个点，计算所有特征图像中该点搜索区域内所有像素点的灰度与训练特征图像中相应点灰度的相似程度。选出30个相似程度最大的点，作为该边界点的候选点。相似程度为所有特征图像中该点周围像素点灰度到训练样本特征图像中相应点周围像素点灰度集合的马氏距离，定义为： 

$h_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(g_i^j-u_{g_i^j})}^T{\left(s_{g_i^j}\right)}^{-1}(g_i^j-u_{g_i^j})---\left(2\right)$

其中， $g_i^j=p_i+r_c\left[\begin{array}{c}\cos \left(\frac{2\mathrm{\π}}{n_c}(k-1)\right)\\ \sin \left(\frac{2\mathrm{\π}}{n_c}(k-1)\right)\end{array}\right]$ 为在特征图像上，以点p_i为圆心，r_c为半径的圆上的n_c个点的灰度，

分别为训练图像的第j个特征图像中第i个边界点的周围像素点灰度的均值及协方差。N为特征图像总 数，这里取值为6。 

（3）基于动态规划的肺分割 

①灰度相似性代价：在边界点搜索区域内，像素点的灰度相似性代价为该点周围像素点灰度与训练图像中相应边界点的周围像素点灰度的相似程度，如公式（2）所示。测试图像中某个边界候选点的h_i值越小，表明该点周围点的灰度分布与相应边界点训练样本的相似性越高。 

②形状相似性代价：图像中第i个边界点的形状相似性代价定义为： 

$f_i={(v_i-u_{v_i})}^T{\left(s_{v_i}\right)}^{-1}(v_i-u_{v_i})---\left(3\right)$

其中，v_i＝p_i+1-p_i，表示第i个边界点的形状特征，

分别表示所有训练图像中第i个边界点形状特征的均值与协方差。 

③基于灰度与形状相似性信息的最优边界搜索 

对于测试图像中第i个边界点p_i，在指定的搜索区域内存在m个具有较小灰度相似性代价候选点，则n个边界点就会产生一个n×m的灰度代价矩阵： 

$C=\left[\begin{array}{ccccc}{h}_{\mathrm{1,1}}& ...& h_{1,k}& ...& h_{1,m}\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ h_{i,1}& ...& h_{i,k}& ...& h_{i,}m\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ h_{n,1}& ...& h_{n,k}& ...& h_{n,m}\end{array}\right]---\left(4\right)$

搜索最优轮廓就是找一条最佳路径（矩阵C中每行选一个元素），沿着路径选择的时候，灰度与形状相似性代价的总和最小，即： 

$J(k_1.....k_n)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{h}_i+\mathrm{\γ}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{f}_i---\left(5\right)$

其中，γ为灰度与形状相似性代价系数。调整γ值，使得这两种代价在边界搜索过程中发挥大致相同的作用。 

三、基于ASM算法的肺边界修正：由于在某些图像中，初始位置可能与实际边界相距过远，利用灰度与形状相似性信息分割时，搜索区域不覆盖肺边界。因此，本发明通过使用ASM算法在第一次分割基础上修正肺边界，改善部分点搜索陷入局部极值的情况，获得更优的搜索结果。 

在边界修正阶段，通过利用测试样本所有特征图像中边界点梯度方向的灰度分布与训练样本所有特征图像边界点梯度方向灰度分布的加权总和(即马氏距离)最小来修正肺边界。 

对于测试图像中的第i个点，可以在第j张特征图像中找到一个以该点为中心，长度为2m+l个像素，并且方向为该点所在边界法线方向的派生轮廓线。如果轮廓线上任意一点标准化后的局部灰度向量为g_s，则该点是否为边界最优点可由所有特征图像中马氏距离的和来度量，即： 

$f\left(g_s\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(g_s-\stackrel{\‾}{g})}^T{S}^{-1}(g_s-\stackrel{\‾}{g})---\left(6\right)$

f(g_s)值越小，表明该点处法线轮廓线上灰度的分布越近似肺部真实边界点法线轮廓线上灰度的分布，该点为最优边界点的可能性越大。所以说，通过在边界点轮廓线上搜寻f(g_s)最小值的点，可以获得最优的肺部边界。 

修正边界后，调整模型的姿态和形状参数确定最终的分割结果。由公式（7）计算形状参数调整量dx： 

M(s(1+ds),(θ+dθ))[x+dx]+(X_c+dX_c)＝(X+dX)    （7） 

可得 

dx＝M((s(1+ds))^-1,-(θ+dθ))[M(s,θ)[x]+dX-dX_c]-x    （8） 

其中，缩放因子1+ds，旋转因子1+dθ可以通过匹配当前点x与搜索得到的新位置X+dX得到。由公式（8）得到的形状参数，用该参数调整初始 轮廓结果通常与形状模型不一致。我们希望找到db,使得

通过调整形状参数b+db，使b_k+db_k在

范围内，来约束形状模型。得到调整量的变化dX_c，dY_c，ds，dθ以及形状参数变化量db，在迭代过程中应用这些变化量来更新参数，获得合适的匹配结果。 

本发明先建立肺轮廓的先验模型，再利用灰度与形状相似性信息结合图像特征对肺区域分割。由于在某些图像中，初始位置可能与实际边界相距过远，利用灰度与形状相似性信息分割时，搜索区域不覆盖肺边界。因此，本发明通过使用主动形状模型（Active Shape Model,ASM)算法在第一次分割基础上修正肺边界，改善部分点搜索陷入局部极值的情况，获得更优的搜索结果。 

附图说明：

图1是使用三种算法分割结果的DR图像。 

图2是本发明的流程框图。 

具体实施方式：

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清晰、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。给予本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。 

本发明使用日本反射技术学会（Japanese Society of Radiological Technology）建立的公用数据库。该数据库由247张PA(Posterior-Anterior)图像组成。其中，93张正常图像，154张有结节的图像。图像是由LD-4500或者LD-5500激光胶片数字化仪获得（Konica,Tokyo,Japan）。大小为 2048*2048像素，每个像素为12bits，且每个像素的尺度为0.175*0.175mm²。 

实验结果分析 

当（1）肺边界点的初始候选点个数为20，（2）形状代价权重为1*10⁶，（3）右肺第4-16，27-39个标记点的搜索区域为以该点为中心，70*70像素的正方形区域，其余标记点为90*90像素区域，左肺第4-16，36-39个标记点的搜索区域为以该点为中心，70*70像素的正方形区域，其余标记点为120*120像素区域，（4）使用6种特征图像搜索肺边界时，获得最佳的分割性能，分割的overlap达到88.90。 

1、不同候选点个数的分割结果分析 

在每个标记点的搜索区域中选10，20，30，40个灰度代价最小的点作为肺边界候选点。若候选点数目太少，则灰度代价相对稍大的真实肺边界点可能不包括在内；如果候选点太多，则会引入噪声、肋骨边缘等干扰点，降低分割性能。表1为在其它参数一定时，不同候选点个数的肺分割性能。 

表1不同候选点个数分割的结果 

候选点个数	10	20	30	40
					Overlap(%)	87.63	88.90	88.53	88.64

2、不同权重的分割结果分析 

调节灰度与形状代价权重可以使二者处于相近的数量级，达到调节灰度变化强弱及边界平滑强度的作用。如果灰度代价过大，分割结果受噪声影响较大；如果形状作用过大，会搜索不到具有灰度变化的肺边界。因此，选择合适的权重系数十分重要。表2为在其它参数一定，灰度代价权重为1，形状代价权重分别为1×10³，1×10⁴，1×10⁵，1×10⁶，1×10⁷时，肺分割 的性能。 

表2不同形状代价权重的分割结果 

形状代价权重	1×103	1×104	1×105	1×106	1×107
						Overlap(%)	87.05	86.79	86.29	88.90	86.19

3、不同尺度搜索区域的分割结果分析 

在每个初始边界点相应的搜索区域中，寻找与真实肺边界点灰度信息相似性最大的候选边界点。如果搜索区域过大，包括肋骨等边界区域，产生大量干扰点；如果搜索区域过小，可能不覆盖真实肺边界，使搜索陷入局部极值。 

表3为其它参数一定，改变搜索区域面积获得的肺分割的性能。其中，L_i,R_i（i=1,2,...,42）表示左、右肺第i个标记点搜索区域边长。 

不同特征图像的分割结果分析 

表4为利用特征图像及其组合进行边界搜索获得分割性能。从中可以看出，当同时使用6种特征图像时，充分利用图像各个方向灰度的变化信息，获得最佳的分割结果。 

不同分割算法性能比较 

对于（1）基于灰度与形状相似性分割算法，（2）利用特征图像的基于灰度与形状相似性分割算法及（3）在此基础上利用活动形状模型修正分割结果的这三种分割算法，分割结果如表5所示。图1为使用三种算法分割结果的DR图像。（a)组为无特征的分割结果，(b)组为有特征的分割结果，(c)组为有特征且修正边界的分割结果。 

表3不同搜索区域边长的分割结果 

表4求特征图像的方法 

表5三种分割算法性能比较 

三种分割算法	无特征图像	有特征图像	有特征图像及边界修补
				Overlap(%)	82.18	87.99	88.90

Claims (6)

1.基于改进形状模型的肺分割方法，具体步骤为： 

一、模型初始轮廓位置的确定：包括标记训练集，对齐训练集，先验模型的建立三个主要部分； 

①标记训练集：使用沿着边界的点来标记训练集，点包括如下三类：a、标记目标有特定应用部分的点，b、标记无关应用部分的点，c、填在一类点或二类点之间的点； 

②对齐训练集：通过缩放、旋转和平移的操作使训练形状对齐，使它们尽可能对齐紧密，设x_i是训练集中第i个形状中n个点的向量，x_i＝(x_i0,y_i0,x_i1,y_i1,...,x_ik,y_ik,...,x_in-1,y_in-1)^T，其中，(x_ij,y_ij)是第i个形状中第j个点， 

给定两个相似形状x_i和x_j，选择旋转角度θ、缩放s、平移(t_x,t_y)，则M(s,θ)[x]代表旋转角度为θ和缩放比例为s的变换，把x_i映射为M(s,θ)[x_j]+t，使以下加权和最小化： 

E_j＝(x_i-M(s,θ)[x_j]-t)^TW(x_i-M(s,θ)[x_j]-t)   (1) 

其中， 

t＝(t_x,t_y,...,t_x,t_y)^T； 

③先验模型的建立：把训练集中的形状向量对齐处理之后，利用主成分分析的方法来找出形状变化的统计信息及规律； 

设平均形状为对齐后各个样本相对于平均形状的偏差构成的形状向量协方差为

计算该协方差的特征值与特征 向量，并将特征值排序Sp_k＝λ_kp_k，其中，λ_k表示第k大的特征值，λ_k越大，它所对应的p_k描述数据点变化模式就越重要，选择前t个重要的变化模式组成新的主轴系p_s，则允许形状域内的任何一个形状都可以由平均形状加主轴系与一组控制参数的加权和来近似，

其中，p_s＝(p₁p₂...p_t)是前t个特征向量组成的矩阵，b_s＝(b₁b₂...b_t)^T是权值向量，经过主分量分析，按照降序取前t个特征值及对应的特征向量，使前t个特征值所决定的目标物体形变占所有2n个特征值所决定目标物体形变总量的比例不小于V； 

二、基于灰度信息和形状信息的肺实质分割：同时利用图像的灰度与形状信息，使得搜索到的边界灰度、形状信息与训练图像相似； 

（1）特征图像：利用特征图像来获取边界点的候选点及各个候选点的灰度代价，采用6种特征图像：（1-2）x,y方向一阶偏导数图像，表示x,y方向灰度变化；（3-4）x,y方向二阶偏导数图像，表示x,y方向灰度变化速率；（5）x，y方向混合偏导数图像；（6）x，y方向二阶偏导数和图像，该值越大，表明此处灰度变化速度越快，为肺边界的可能性越大； 

（2）边界点的候选点：对于初始肺边界的每一个点，计算所有特征图像中该点搜索区域内所有像素点的灰度与训练特征图像中相应点灰度的相似程度，相似程度为所有特征图像中该点周围像素点灰度到训练样本特征图像中相应点周围像素点灰度集合的马氏距离，定义为： 

   （2） 

其中，为在特征图像上，以点p_i为圆心，r_c为半 径的圆上的n_c个点的灰度，

分别为训练图像的第j个特征图像中第i个边界点的周围像素点灰度的均值及协方差，N为特征图像总数，这里取值为6； 

（3）基于动态规划的肺分割 

①灰度相似性代价：在边界点搜索区域内，像素点的灰度相似性代价为该点周围像素点灰度与训练图像中相应边界点的周围像素点灰度的相似程度，如公式（2）所示，测试图像中某个边界候选点的h_i值越小，表明该点周围点的灰度分布与相应边界点训练样本的相似性越高； 

②形状相似性代价：图像中第i个边界点的形状相似性代价定义为： 

   (3) 

其中，v_i＝p_i+1-p_i，表示第i个边界点的形状特征，分别表示所有训练图像中第i个边界点形状特征的均值与协方差； 

③基于灰度与形状相似性信息的最优边界搜索 

对于测试图像中第i个边界点p_i，在指定的搜索区域内存在m个具有较小灰度相似性代价候选点，则n个边界点就会产生一个n×m的灰度代价矩阵： 

   (4) 

搜索最优轮廓就是找一条最佳路径，沿着路径选择的时候，灰度与形状相似性代价的总和最小，即： 

   (5) 

其中，γ为灰度与形状相似性代价系数，调整γ值，使得这两种代价在边界搜索过程中发挥大致相同的作用； 

三、基于ASM算法的肺边界修正：通过使用ASM算法在第一次分割基础上修正肺边界，改善部分点搜索陷入局部极值的情况，获得更优的搜索结果； 

在边界修正阶段，通过利用测试样本所有特征图像中边界点梯度方向的灰度分布与训练样本所有特征图像边界点梯度方向灰度分布的加权总和，即马氏距离，最小来修正肺边界； 

对于测试图像中的第i个点，可以在第j张特征图像中找到一个以该点为中心，长度为2m+l个像素，并且方向为该点所在边界法线方向的派生轮廓线，如果轮廓线上任意一点标准化后的局部灰度向量为g_s，则该点是否为边界最优点可由所有特征图像中马氏距离的和来度量，即： 

   (6) 

f(g_s)值越小，表明该点处法线轮廓线上灰度的分布越近似肺部真实边界点法线轮廓线上灰度的分布，该点为最优边界点的可能性越大，通过在边界点轮廓线上搜寻f(g_s)最小值的点，可以获得最优的肺部边界； 

修正边界后，调整模型的姿态和形状参数确定最终的分割结果，由公式（7）计算形状参数调整量dx： 

M(s(1+ds),(θ+dθ))[x+dx]+(X_c+dX_c)＝(X+dX)   （7） 

可得 

dx＝M((s(1+ds))^-1,-(θ+dθ))[M(s,θ)[x]+dX-dX_c]-x   （8） 

其中，缩放因子1+ds，旋转因子1+dθ可以通过匹配当前点x与搜索得 到的新位置X+dX得到，由公式（8）得到的形状参数，用该参数调整初始轮廓结果通常与形状模型不一致，找到db,使得通过调整形状参数b+db，使b_k+db_k在

范围内，来约束形状模型，得到调整量的变化dX_c，dY_c，ds，dθ以及形状参数变化量db，在迭代过程中应用这些变化量来更新参数，获得合适的匹配结果。 

2.如权利要求1所述的基于改进形状模型的肺分割方法，其特征在于：①标记训练集中所述的标记目标有特定应用部分的点为在人脸模型中，表示眼睛中心的点，或表示边界上比较锐利的拐角。 

3.如权利要求2所述的基于改进形状模型的肺分割方法，其特征在于：所述的表示边界上比较锐利的拐角为眼角的点。 

4.如权利要求1所述的基于改进形状模型的肺分割方法，其特征在于：①标记训练集中所述的标记无关应用部分的点为目标在特定方向上的最高点，或曲率的极值处。 

5.如权利要求1所述的基于改进形状模型的肺分割方法，其特征在于：所述的V取0.98。 

6.如权利要求1所述的基于改进形状模型的肺分割方法，其特征在于：步骤（2）边界点的候选点中，选出20个相似程度最大的点作为该边界点的候选点。 

Images