CN102831614B - 基于交互式字典迁移的序列医学图像快速分割方法 - Google Patents

Abstract

本发明公开了一种基于交互式字典迁移的序列医学图像快速分割方法。主要解决现有单幅图像分割速度慢和工作量大的问题。其实现过程是：输入腹部CT序列图像；用mean shift方法对其初分割并以每个像素为中心提取5×5块作为特征；对一幅图像进行交互式分割得到目标和背景；对目标和背景区域分别进行训练生成目标和背景字典后，去逼近序列中未分割的图像，根据误差最小原则将所有像素标记为目标或背景并进行高斯平滑；计算初分割图像中每个区域的目标像素点的密度，去掉密度小的区域后得到分割结果；对于分割不正确的图像进行再次交互式分割；通过字典迁移来更新目标字典。本发明具有对医学图像分割效果好和速度快的优点，可用于对医学CT序列图像的分割。

Description

基于交互式字典迁移的序列医学图像快速分割方法

技术领域

本发明属于图像处理技术领域，特别是医学图像的处理，主要可用于多幅腹部CT序列图像的胃癌病变区域分割。

背景技术

医学图像分割是医生根据需要把感兴趣的目标区域，即将可能存在病变区域、重要组织和器官等从背景中分离出来的过程。它是现代医学图像处理中的一个重要的研究领域，也是病变区域提取、检测和跟踪等后续医学图像处理的基础。由于这些操作都需要以一个准确的分割结果为前提，然后才能进行准确地定位分析，所以医学图像分割结果的准确性对于医生对病情的诊断并做出正确治疗方案具有非常重要的意义。

图像特征是图像的重要属性，图像分割过程中提取特征是重要的一步，常用的图像特征有灰度共生矩阵，小波变换特征等。由于图像分割方法都要对图像进行特征提取，特征提取虽然是线下过程，但是由于需要对图像的每一个像素点分别进行特征提取，所以这一过程十分耗时。

交互式分割方法是现在应用最热门和最广泛的处理医学图像的方法。它结合了自动分割与手动分割的优点，既允许人工参与提供一些先验信息，但又不完全依赖于人工参与。由于医学图像有复杂的背景，自动分割运算量大，在没有人工指导的情况下准确率很难保证。而交互式分割在允许人工参与的情况下保证了准确率。它的具体做法是由用户标记出图像中感兴趣区域的大致范围，从而提供部分先验信息，然后分割算法再根据这些信息对图像进行分割。虽然交互式分割有以上优点，但它也存在缺点，由于交互式分割方法一次只能分割一幅图像，用户必须对每幅待分割的图像标记出目标和背景，然后一幅一幅分割。对于CT序列图像，一个病人的序列图有上百幅，如果单幅分割，则极大的增加了工作量。

目前交互式分割方法主要有Live wire方法和活动轮廓模型方法等。Live wire方法是通过搜索边缘来进行图像分割的方法，它将图像看成是一个连通图，通过构造代价函数和选择关键点，在全图范围内寻找给定两点间的最优路径作为物体的边缘。这个方法的缺点是它对交互者要求较高，输入交互点位置的偏移可能导致分割结果出错，而且在多幅序列图像中如果出现分割结果不好的图像时也不能及时修正分割结果。活动轮廓模型方法是在图像中定义一个曲线，使其在曲线本身决定的内能和由图像数据决定的外能的推动下向物体边缘靠近，其缺点是对初始曲线的位置比较敏感，而且它不能收敛到凹性边缘。这两种方法都是基于边缘的方法。对于上百幅的医学序列图像，感兴趣区域的边缘变化很大，这些方法在分割多幅医学图像时的分割准确率较低。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于交互式字典迁移的序列医学图像快速分割方法，以减少单幅交互式分割的工作量，节省时间，降低对交互者的要求，并采用二次修正的方法提高分割准确率。

实现本发明上述目的技术方案，包括如下步骤：

(1)输入腹部CT序列图像I={I⁽¹⁾,I⁽²⁾,…,I⁽ⁿ⁾}，n为从序列图像的数量，n最大值为300；

(2)对每幅图像I(q)以每个像素点为中心提取5×5的块作为每个像素的特征，则该图像所有像素的特征向量集合为X_N×25 ^(q)，N为每幅图像的大小，其中1≤q≤n；

(3)用mean shift方法对CT序列图像I进行初始分割，则CT序列图像I^(q)的初分割图像为S^(q)，并统计初分割图像S^(q)中的区域数目m；

(4)选取CT序列图像I⁽¹⁾，进行交互式分割：

首先，在序列图像I⁽¹⁾对应的初分割图像S⁽¹⁾上的区域进行手动标记部分目标区域M_O和部分背景区域M_B，剩余的区域为未标记区域N；

然后，根据合并准则，重复合并未标记区域N中与标记背景M_B中相似性最大的区域，直到所有区域都合并后，得到CT序列图像I⁽¹⁾的全部目标区域和背景区域，即得到分割结果，其中该目标区域代表腹部CT图像中的胃癌病变区域；

(5)提取步骤(4)中的目标区域内对应像素点的特征向量X_O，并用K-SVD算法对目标特征向量X_O进行训练学习，生成一个目标字典D_O；

(6)用k均值方法将步骤(4)中得到的背景区域聚为4类，每类对应的特征向量为X_B1,X_B2,X_B3,X_B4，并用K-SVD方法分别对其训练学习，生成背景字典D_B={D_B1,D_B2,D_B3,D_B4}，其中D_B1,D_B2,D_B3,D_B4分别为4类背景区域对应的背景字典；

(7)由目标字典D_O和背景字典D_B分别去逼近CT序列图像I^(q)的特征向量X_N×25 ^(q)，得到所有像素点的目标字典逼近误差e_O和背景字典逼近误差e_B，然后根据误差最小原则计算二值图像T^(q)：

$T^{\left(q\right)}\left(i\right)=\left\{\begin{array}{cc}1& \mathrm{if}{e}_O\left(i\right)<e_B\left(i\right)\\ 0& \mathrm{otherwise}\end{array}\right.,$

其中，2≤q≤n，i为像素点，其范围为1≤i≤N；

(8)对二值图像T^(q)进行高斯滤波处理后，计算二值图像T^(q)中值为1的像素点在初分割图像S^(q)上对应区域的密度d(r)：

$d\left(r\right)=\frac{\underset{i\&Element;r}{\mathrm{\&Sigma;}}I(T^{\left(q\right)}\left(i\right)=1)}{\underset{i\&Element;r}{\mathrm{\&Sigma;}}I(T^{\left(q\right)}\left(i\right)=1)+\underset{i\&Element;r}{\mathrm{\&Sigma;}}I(T^{\left(q\right)}\left(i\right)=0)},$

其中，r为初分割图中的区域，其范围为1≤r≤m，i为像素点，其范围为1≤i≤N，m为初分割图像中的区域总数目；I(·)是指示函数，当该函数输入的等式成立时其函数值为1，反之函数值为0；

(9)根据步骤(8)中的密度d(r)得到CT序列图像I^(q)的分割结果图像C^(q)：

$C^{\left(q\right)}\left(r\right)=\left\{\begin{array}{cc}1& \mathrm{ifd}\left(r\right)>0.3\\ 0& \mathrm{otherwise}.\end{array}\right.,$

其中，2≤q≤n，C^(q)中值为1的区域为目标区域，值为0的区域为背景区域；

(10)判断分割结果图像C^(q)分割是否正确，即计算分割结果图C^(q)与上一幅分割结果图C^(q-1)的信息差异指标VI,当信息差异指标VI>0.25时，判定为分割结果不正确，返回步骤(4)开始执行，即对序列图像I^(q)进行交互式分割；否则，执行步骤(11)；

(11)计算步骤(9)得到的分割结果图像C^(q)的目标区域的特征向量集合X′_O，并用K-SVD方法对其训练学习，得到一个新目标字典D′_O，计算新目标字典D′_O与原目标字典D_O的误差E=||D_O-D′_O||₂，若误差E<1，则将新目标字典D′_O中的原子迁移到原目标字典D_O中，更新目标字典D_O；否则，原目标字典D_O不变；

(12)统计步骤(11)中字典的迁移次数，当字典迁移次数等于4时，计算当前目标字典D_O中原子的贡献度g，并按贡献度g的大小进行降序排列，保留贡献度g最大的前50个原子，构造新的目标字典D_O；

(13)重复执行步骤(7)到步骤(10)，对输入的腹部CT序列图像I中还未分割的图像进行目标区域提取，直到所有序列图像都分割出目标区域为止。

本发明与现有的技术相比具有以下优点：

1．本发明采用迁移思想，将一幅图像的先验信息迁移到其它图像来进行分割，从而极大的节省了时间并减少了工作量；

2.本发明对分割不准确的图像进行二次交互式分割，从而提高了图像的分割准确率；

3.本发明通过以灰度像素块来代替特征提取，从而提高了速度；

附图说明

图1是本发明的流程图；

图2是本发明所使用的三幅腹部CT序列图像；

图3是本发明对图2进行初始分割后的图像；

图4是本发明对图2中第一幅图像交互式分割后的图像；

图5是本发明对图2中第二幅图像分割后的图像；

图6是本发明对图2中第三幅图像分割后的图像。

具体实施方式

参照图1，本发明基于交互式字典迁移的序列医学图像快速分割方法给出以下实施例：

步骤1：输入腹部CT序列图像，如图2a、图2b和图2c。

步骤2：对图2a、图2b和图2c分别以每个像素点为中心提取5×5的块作为每个像素的特征，则图2a、图2b和图2c中所有像素的特征向量分别为X⁽¹⁾,X⁽²⁾,X⁽³⁾。

步骤3：用mean shift方法分别对图2a、图2b和图2c分别进行初始分割得到图3a、图3b和图3c。

步骤4：选取图2a，进行交互式分割：

首先，在初分割图3a上进行手动标记，如图4a所示，用实线段标记部分目标区域M_O，用虚线段标记部分背景区域M_B，剩余的区域为未标记区域N；

然后，根据最大相似性的合并准则，重复合并未标记区域N中与标记背景M_B中相似性最大的区域，直到所有区域都合并后，得到图2a的全部目标区域和背景区域，分割结果二值图像如图4b所示，该目标区域代表腹部CT图像中的胃癌病变区域。

步骤5：计算分割结果图4c中的目标区域内所有像素点的特征向量X_O，并用K-SVD方法对目标特征向量X_O进行训练学习，生成一个目标字典D_O。

步骤6：用k均值方法将图4c中的背景区域聚为4类，每类对应的特征向量为X_B1,X_B2,X_B3,X_B4，并用K-SVD方法分别对其训练学习，生成背景字典：D_B={D_B1,D_B2,D_B3,D_B4}，其中D_B1,D_B2,D_B3,D_B4分别为4类背景区域对应的背景字典。

步骤7：由目标字典D_O和背景字典D_B={D_B1,D_B2,D_B3,D_B4}分别去逼近图2b的特征向量X⁽²⁾，得到所有像素点的目标字典逼近误差e_O和背景字典的逼近误差e_B，然后根据误差最小原则计算二值图像T^(q)：

$T^{\left(q\right)}\left(i\right)=\left\{\begin{array}{cc}1& \mathrm{if}{e}_O\left(i\right)<e_B\left(i\right)\\ 0& \mathrm{otherwise}\end{array}\right.,$

其中，2≤q≤3，i为像素点，其范围为1≤i≤k，k为序列图像I^(q)的像素点总数。

步骤8：对二值图像T^(q)进行高斯滤波处理后得到图5a，计算图5a中值为1的像素点在图3b的对应区域中的密度d(r)：

$d\left(r\right)=\frac{\underset{i\&Element;r}{\mathrm{\&Sigma;}}I(T^{\left(q\right)}\left(i\right)=1)}{\underset{i\&Element;r}{\mathrm{\&Sigma;}}I(T^{\left(q\right)}\left(i\right)=1)+\underset{i\&Element;r}{\mathrm{\&Sigma;}}I(T^{\left(q\right)}\left(i\right)=0)},$

其中，i为像素点，其范围为1≤i≤k，r为初分割图中的区域，其范围为1≤r≤m，m为初分割图像中的区域总数目；I(·)是指示函数，当该函数输入的等式成立时其函数值为1，反之函数值为0。

步骤9：根据步骤8中的密度d(r)得到图2b的分割结果图C^(q)：

$C^{\left(q\right)}\left(r\right)=\left\{\begin{array}{cc}1& \mathrm{ifd}\left(r\right)>0.3\\ 0& \mathrm{otherwise}.\end{array}\right.,$

其中，二值图像C^(q)中值为1的区域为目标区域，值为0的区域为背景区域，分割结果如图5b所示。

步骤10：判断分割结果二值图5b是否正确，具体步骤如下：

10a)分别计算分割结果图C^(q)和上一幅分割结果图C^(q-1)的信息熵H(C)：

$H\left(C\right)=-\underset{k=1}{\overset{2}{\mathrm{\&Sigma;}}}P\left(k\right)\log P\left(k\right),$

其中，C为分割结果图C^(q)或上一幅分割结果图C^(q-1)，为目标或背景在整幅图像中的概率，N为分割结果图像的大小，即所有像素点的数量，N_k为第k类区域的像素点数目，k=1为目标类，k=2为背景类；

10b)计算分割结果图C^(q)和上一幅分割结果图C^(q-1)的互信息I(C^(q),C^(q-1))：

$I(C^{\left(q\right)},C^{(q-1)})=\underset{k=1}{\overset{2}{\mathrm{\&Sigma;}}}\underset{k^{\&prime;}=1}{\overset{2}{\mathrm{\&Sigma;}}}P(k,k^{\&prime;})\log \frac{P(k,k^{\&prime;})}{P\left(k\right)P^{\&prime;}\left(k^{\&prime;}\right)},$

其中，为两幅图像中相关的目标和背景的概率，1≤k≤2，1≤k′≤2，为分割结果图C^(q)的第k类，为分割结果图C^(q)的第k′类；

10c)计算分割结果图C^(q)与上一幅分割结果图C^(q-1)的信息差异指标VI:

VI(C^(q),C^(q-1))=H(C^(q))+H(C^(q-1))-2I(C^(q),C^(q-1)),

其中，H(C^(q))为分割结果图的信息熵，H(C^(q-1))为分割结果图的信息熵，I(C^(q),C^(q-1))为分割结果图C^(q)和上一幅分割结果图C^(q-1)的互信息；

10d)当信息差异指标VI>0.25时，判定为分割结果不正确，返回步骤4开始执行，即对图2b进行交互式分割；否则，分割结果正确，继续执行步骤11。

步骤11：计算分割结果图5c中目标区域的特征向量X′_O，并用K-SVD方法对其训练学习，得到一个新目标字典D′_O，计算新目标字典D′_O与原目标字典D_O的误差E=||D_O-D′_O||₂，若误差E<1，则将新目标字典D′_O中的原子迁移到原目标字典D_O中，更新目标字典D_O；否则，保持原目标字典D_O不变。

步骤12：统计步骤11中字典的迁移次数，当字典迁移次数等于4时，根据原子贡献度g的大小来更新目标字典D_O，具体步骤如下：

12a)每次移除目标字典D_O中的一个原子，计算剩余原子对特征向量X′_O的逼近误差

$R(X_O^{\&prime;},D_O^j)=\left|\right|X_O^{\&prime;}-D_O^j{\&PartialD;}_j\left|\right|,$

其中，为去掉目标字典D_O中第j个原子后的字典，为字典的稀疏表示系数，1≤j≤atom，atom为目标字典D_O中的原子个数；

12b)计算每个原子的贡献度g：

$g\left(j\right)=\frac{1}{R(X_O^{\&prime;},D_O^j)},$

其中，为去掉目标字典D_O中第j个原子后的逼近误差；

12c)将当前目标字典D_O中原子按贡献度g的大小进行降序排列，保留贡献度g最大的前50个原子，并去掉其余原子作为目标字典D_O。

步骤13：重复执行步骤7到步骤10，对输入的腹部CT序列图像中还未分割的图2c进行目标区域提取，直到所有CT序列图像都分割出目标区域为止，最终分割结果图像如图6所示,其中图6a为分割结果二值图，图6b为分割结果图。

从图5c和图6b可以看出，用本发明能很好的将序列图像中的胃癌病变区域分割出来，并且对分割效果不好的图像，通过判断可以让用户再次交互，从而得到准确的分割结果，这说明本发明不仅分割准确率高，而且速度快。

Claims (4)

1.一种基于交互式字典迁移的序列医学图像快速分割方法，包括如下步骤：

(1)输入腹部CT序列图像I＝{I⁽¹⁾,I⁽²⁾,…,I⁽ⁿ⁾}，n为序列图像的数量，n最大值为300；

(2)对每幅图像I^(q)以每个像素点为中心提取5×5的块作为每个像素的特征，则该图像所有像素的特征向量集合为X_N×25 ^(q)，N为每幅图像的大小，其中1≤q≤n；

(3)用mean shift方法对CT序列图像I进行初始分割，则CT序列图像I^(q)的初分割图像为S^(q)，并统计初分割图像S^(q)中的区域数目m；

(4)选取CT序列图像I⁽¹⁾，进行交互式分割：

首先，在序列图像I⁽¹⁾对应的初分割图像S⁽¹⁾上的区域进行手动标记部分目标区域M_O和部分背景区域M_B，剩余的区域为未标记区域N；

然后，根据合并准则，重复合并未标记区域N中与标记背景M_B中相似性最大的区域，直到所有区域都合并后，得到CT序列图像I⁽¹⁾的全部目标区域和背景区域，即得到分割结果，其中该目标区域代表腹部CT图像中的胃癌病变区域；

(5)提取步骤(4)中的目标区域内对应像素点的特征向量X_O，并用K-SVD算法对目标特征向量X_O进行训练学习，生成一个目标字典D_O；

(6)用k均值方法将步骤(4)中得到的背景区域聚为4类，每类对应的特征向量为X_B1,X_B2,X_B3,X_B4，并用K-SVD方法分别对其训练学习，生成背景字典

D_B＝{D_B1,D_B2,D_B3,D_B4}，其中D_B1,D_B2,D_B3,D_B4分别为4类背景区域对应的背景字典；

(7)由目标字典D_O和背景字典D_B分别去逼近CT序列图像I^(q)的特征向量X_N×25 ^(q)，得到所有像素点的目标字典逼近误差e_O和背景字典逼近误差e_B，然后根据误差最小原则计算二值图像T^(q)：

$T^{\left(q\right)}\left(i\right)=\left\{\begin{array}{cc}1& \mathrm{if}{e}_O\left(i\right)<e_B\left(i\right)\\ 0& \mathrm{otherwise}\end{array}\right.,$

其中，2≤q≤n，i为像素点，其范围为1≤i≤N；

(8)对二值图像T^(q)进行高斯滤波处理后，计算二值图像T^(q)中值为1的像素点在初分割图像S^(q)上对应区域的密度d(r)：

$d\left(r\right)=\frac{\underset{i\&Element;r}{\mathrm{\&Sigma;}}I(T^{\left(q\right)}\left(i\right)=1)}{\underset{i\&Element;r}{\mathrm{\&Sigma;}}I(T^{\left(q\right)}\left(i\right)=1)+\underset{i\&Element;r}{\mathrm{\&Sigma;}}I(T^{\left(q\right)}\left(i\right)=0)},$

其中，r为初分割图中的区域，其范围为1≤r≤m，i为像素点，其范围为1≤i≤N，m为初分割图像中的区域总数目；I(·)是指示函数，当该函数输入的等式成立时其函数值为1，反之函数值为0；

(9)根据步骤(8)中的密度d(r)得到CT序列图像I^(q)的分割结果图像C^(q)：

$C^{\left(q\right)}\left(r\right)=\left\{\begin{array}{cc}1& \mathrm{if\; d}\left(r\right)>0.3\\ 0& \mathrm{otherwise}.\end{array}\right.,$

其中，2≤q≤n，C^(q)中值为1的区域为目标区域，值为0的区域为背景区域；

(10)判断分割结果图像C^(q)分割是否正确，即计算分割结果图C^(q)与上一幅分割结果图C^(q-1)的信息差异指标VI,当信息差异指标VI＞0.25时，判定为分割结果不正确，返回步骤(4)开始执行，即对序列图像I^(q)进行交互式分割；否则，执行步骤(11)；

(11)计算步骤(9)得到的分割结果图像C^(q)的目标区域的特征向量集合X′_O，并用K-SVD方法对其训练学习，得到一个新目标字典D′_O，计算新目标字典D′_O与原目标字典D_O的误差E＝||D_O-D′_O||₂,若误差E＜1，则将新目标字典D′_O中的原子迁移到原目标字典D_O中，更新目标字典D_O；否则，原目标字典D_O不变；

(12)统计步骤(11)中字典的迁移次数，当字典迁移次数等于4时，计算当前目标字典D_O中原子的贡献度g，并按贡献度g的大小进行降序排列，保留贡献度g最大的前50个原子，构造新的目标字典D_O；

(13)重复执行步骤(7)到步骤(10)，对输入的腹部CT序列图像I中还未分割的图像进行目标区域提取，直到所有序列图像都分割出目标区域为止。

2.根据权利要求1所述的方法，其中步骤(4)所述的合并准则是最大相似性准则，它先计算各个区域间的欧式距离，然后依次合并欧氏距离最大的两个区域直到所有区域合并完毕，合并时取每个区域中出现频率最高的灰度值作为此区域的特征值。

3.根据权利要求1所述的方法，其中所述步骤(7)中计算目标字典D_O和背景字典D_B＝{D_B1,D_B2,D_B3,D_B4}的逼近误差，是通过公式计算，即计算待逼近的特征向量X减Dα的2范数的平方值，其中D为目标字典D_O或背景字典D_B；X为待逼近的特征向量，α为逼近系数，e为特征向量X与字典D的最小逼近误差。

4.根据权利要求1所述的方法，其中所述步骤(12)中的贡献度g，是指每个原子对目标数据重构的贡献，其计算步骤为：

11a)每次移除目标字典D_O中的一个原子，计算剩余原子对特征向量X′_O的逼近误差

$R(X_O^{\&prime;},D_O^j)=\left|\right|X_O^{\&prime;}-D_O^j{\&PartialD;}_j\left|\right|,$

其中，为去掉目标字典D_O中第j个原子后的字典，为字典的稀疏表示系数，1≤j≤atom，atom为目标字典D_O中的原子个数；

11b)计算每个原子的贡献度g：

$g\left(j\right)=\frac{1}{R(X_O^{\&prime;},D_O^j)},$

其中，为去掉目标字典D_O中第j个原子后的逼近误差。