CN109712138A - 基于外观字典学习和形状稀疏表示的图像分割方法 - Google Patents

Abstract

本发明基于外观字典学习和形状稀疏表示的图像分割方法，包括：将平均网格移动到初始中心位置；沿着平均网格的法线方向定位网格上每一个标志点的位置；结合形状稀疏表示算法获得初始分割结果；在得到的初始分割结果附近使用外观字典学习和法线搜索结合的算法进行标志点调整；对调整后的网格再次使用形状稀疏表示算法约束重建；根据重建的结果结合梯度矢量流场的特性和外观字典学习的标签重建的概率值再次进行标志点的调整；并运用形状稀疏表示算法得到最终的分割结果。本方法充分利用外观字典学习具有分辨能力的外观先验信息和形状稀疏表示具有重建能力的形状先验信息，使得标志点定位算法和稀疏形状表示算法进行互补，最终得到更加精确的低剂量CT分割结果。

Description

基于外观字典学习和形状稀疏表示的图像分割方法

技术领域

本发明属于图像处理领域，利用局部可分辨的外观字典学习和特征空间的形状稀疏表示算法进行的精确分割方法。

背景技术

近50年来，全世界肺病的发病率明显增高，据统计，在欧美某些国家和我国大城市中，肺肿瘤的发病率已居男性各种肿瘤的首位。低剂量的CT图像已经广泛应用于肺病的诊断，特别是肺肿瘤的分析。精确自动的肺分割可以节省医生大量对肺结构标注工作花费的时间。除此之外，肺分割可以提升17％的肺肿瘤和肺结节诊断精确度。因此自动的肺分割已经引起了临床专家极大的兴趣。形状和外观先验作为形变模型的两根支柱，在目标器官分割上相互补充。外观先验可以给形状约束提供目标器官的识别和定位，形状先验矫正错误的外观分割结果。

然而目前国内外很多现有肺分割方法都是基于单一先验(外观或者形状)，例如：深度学习、图割等方法，对于器官其分割的准确性和可靠性都不是很高。另外更多的分割算法偏向于处理正常的数据，对于异常的数据(患病或者拍摄技术缺陷等)很难处理。

发明内容

本发明所要解决的技术问题在于针对上述背景技术中存在的问题，提供一种利用外观字典学习和特征空间的形状稀疏表示的分割方法，同时在网格顶点定位阶段融合了法线和梯度矢量流场两种不同信息。

相对于相关技术，本发明基于外观字典学习和形状稀疏表示的图像分割方法具有以下有益效果：

本发明首次将主成分分析(PCA)训练得到的形状变化模式整合到传统稀疏表示的目标函数中，提出了特征空间的形状稀疏表示算法。在网格顶点的定位阶段首次将可分辨的外观字典学习和法线搜素的损失函数结合，提出了局部表面快速可分辨外观字典学习。本方法充分运用了法线信息和梯度矢量流场信息，使得顶点调整更加准确，最终得到更加精确的肺分割结果。

附图说明

图1低剂量CT图像；

图2为本发明方法在低剂量CT上的部分切片分割结果对比，绿色为金标准，红色为此发明方法的分割结果；

图3为本发明的算法原理框图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

本发明的肺分割方法意在充分合理运用外观先验知识和形状先验知识。如图1为肺部低剂量CT的一个切片图像，其中(a)和(c)是左右肺部的放大图像。具体分割方法如下：

1、对原始的图像进行滤波去噪操作。在手动标注肺部作为分割金标准，即图像中确切的肺部区域，以便最后比较本发明分割方法的分割准确性。将标定的金标准转化为三维网格的形式。每个三维网格的之间关系是固定的，并且具有相同的边和顶点个数。对于待分割的图像使用阈值分割和形态学操作得到二值图像。

2、基于广义哈弗变换的肺定位方法，使用主成分分析(PCA)训练出肺的平均网格，数学模型如下：

其中ψ_k表示第k个训练样本的肺网格；T表示从训练样本空间到参考样本的相似性变换；T^-1表示对应的T的反变换；p_m表示主成分分析(PCA)训练获得的肺网格形状变化模式；λ_m表示主要的模式所对应的权重；k_s表示训练样本的肺网格总数；表示训练得到的肺平均网格。在主成分分析(PCA)训练过程中，左肺和右肺被分开完成，得到各自的平均形状和变换到特征空间的形状变化模式。计算出平均网格的重心，使用广义哈弗变换计算出二值图像中肺的初始化中心，把该平均网格平移到初始化中心得到的第一初始化网格ψ_τ。

3、对于第一初始化网格ψ_τ，计算其上每个顶点的法线值n_i；对于第一初始化网格ψ_τ上的每个顶点，沿着法线n_i的方向计算边界候选点在离散的位置 j＝-n_m，-n_m+1，…，0，…，n_m-1，n_m，表示为如下形式：

v_ij＝v_i+jδn_i (2)

其中δ是搜索步长；n_t表示顶点移动的阈值；n_m表示最大的搜索范围。当顶点v_i在肺内部(在二值图像的前景区域)并且j＞n_t，v_ij沿着法线n_i方向搜索，顶点v_i更新为v_ij；当顶点v_i在肺外部(在二值图像的背景区域)并且|j|＞n_t，v_ij沿着法线n_i负方向搜索，顶点v_i更新为v_ij。

按照这种方式更新第一初始化网格ψ_τ所有顶点得到法线方向调整后的第二初始化网格ψ。为得到精确的分割结果，迭代优化策略被使用。因此第二初始化网格ψ需要被平滑操作。本发明提出了一种特征空间的形状稀疏表示算法，该算法被用来平滑法线方向调整后的第二初始化网格ψ。特征空间的形状稀疏表示数学模型为：

其中ψ表示法线方向调整后的网格(即需要被平滑的网格)；D_s表示所有训练网格坐标形成的形状字典；p_m表示主成分分析(PCA)训练获得的形状变化模式；T表示从形变网格到平均网格的相似变换；T_s表示每个训练网格到平均网格的相似变换；β₁和β₂表示稀疏系数；λ₃和λ₄表示平衡参数。

当ψ和D_s相似变换到平均空间后，通过迭代最小化方式解上式。首先使用快速的迭代阈值收缩算法(FISTA算法)优化下面数学模型：

整体的形状ψ可以用变换到平均空间的形状训练字典进行重建。第二项拟合局部误差e，在p_m中通过参数β₂挑选合适的形变模态并设置模态的比例大小。数学模型表示成如下形式：

上式用FISTA算法优化。最终的平滑后的形状网格ψ表示为：

ψ＝T^-1(T_s(D_s)β₁+P_mβ₂)

对平滑之后的三维网格ψ，计算其前后两次迭代结果之间的顶点变化距离Δd。当Δd小于给定的阈值Δd_t时，结束迭代，得到分割结果第三初始化网格ψ_λ。否则将第二初始化网格ψ赋值给第三初始化网格ψ_τ继续迭代。

4、本发明针对上一轮基于阈值分割图像得到的分割结果第三初始化网格ψ_λ存在的缺陷，本发明提出了用局部的可分辨外观字典和三维网格顶点的法线方向搜索结合的算法对第三初始化网格ψ_λ进行调整。即在目标边界附近调整网格顶点的位置。网格顶点调整的数学模型表示为：

其中δ′是搜索步长；j′＝0，1，…，n_m-1，n_m表示沿着法线方向的候选点的离散位置；n′_i表示顶点v′_i的法线矢量值；h_t表示顶点的平移阈值；h表示局部的可分辨外观字典学习重建的标签概率值。

h由如下方式得到：

训练阶段。对于被调整网格第三初始化网格ψ_λ的每个顶点，根据其与训练样本对应顶点，从所有训练样本中所对应的相同领域位置取三维外观块形成块库 P_l。将每个三维块表示成一维列向量的形式，因此整个三维块库P_l表示成矩阵形式：m_α表示每个块的尺寸；n_α表示P_l中块的个数。和专家标注结合，外观字典学习训练模型表示为：

其中表示学习的外观字典；α表示稀疏编码因子；k_α表示字典元素的个数；L表示P_l中每个块对应的专家标注值；ω表示学习的线性分类器；λ₁和λ₂是平衡因子。

测试阶段。外观字典在测试阶段被模型为；

其中p_t表示在形变网格附近的一个三维块；表示对输入块p_t的稀疏因子。最终输入三维块p_t中心体素的标签值被表示为：

h＝ω*α_t

特征空间的形状稀疏表示被用来对调整之后的网格进行平滑处理。得到平滑后的结果ψ。整个过程迭代优化，直到两次迭代之间的网格顶点距离差Δd小于给定的阈值Δd_t。此时得到该部分的最终结果ψ_o。否则将ψ赋值给ψ_λ继续迭代。

5、为了克服法线方向的缺点，本发明同时提出了一种基于梯度矢量流场顶点调整算法。ψ_o作为该部分网格迭代搜索的初始化结果。一个边缘图被定义为：G＝-G_σ*B。其中G_σ是表示标准差为σ的三维高斯函数；B表示阈值分割和外观字典学习重建下的二值图像。在梯度矢量流场的方向上离散顶点搜索的数学模型为：

其中v_j-1表示形变网格上当前的顶点；表示v_j-1所对应的梯度矢量流场值的角度值；||G(v_j-1)||表示v_j-1所对应的梯度矢量流场的模值；v_j表示沿着梯度矢量流场的方向搜素的下一个目标点；表示v_j所对应的梯度矢量流场值的角度值； ||G(v_j)||表示v_j所对应的梯度矢量流场的模值。满足上式的离散点形成点集合基于该集合，目标点被定义数学模型为：

其中F(·)表示在离散位置的边缘检测函数；γ表示惩罚参数。F(·)被定义如下：

其中G向三角网格的法线n_i投影限制不平衡边缘；将最大值通过sigmoid函数限制为g_max；±保证方向的一致性；表示边缘附近外观字典学习重建的标签值。完成整个搜索后得到网格ψ。利用本发明提出的特征空间的形状稀疏表示对调整之后的网格ψ进行平滑处理。整个过程迭代优化，到设计的损失函数 E小于某个阈值E_t时结束，得到最终结果最终的重建网格ψ_θ。否则将ψ赋值给ψ_o。

损失函数E的设计，将其数学模型表示为：

E＝E_region+λ_eE_Surface

其中λ_e控制区域能量项和表面能量项的平衡，对于区域损失项E_region为：

其中ψ_t表示阈值分割和外观字典学习重建的图像对应的网格；p表示梯度矢量流场上的一个体素；N_r表示归一化参数；ψ_α表示基于梯度矢量流场调整的网格；||G_α(p)||和||G_t(p)||分别表示ψ_α和ψ_t在梯度矢量流场上所对应的赋值。

对于表面损失项为；

其中N_s表示归一化参数；n_s表示网格ψ_β顶点总数；顶点N_v表示顶点v′的邻域点集合；ψ_β表示形状稀疏表示平滑之后的网格。

通过上述的区域和表面图的损失项，梯度矢量流场和形状稀疏表示交叉迭代优化，最终得到肺分割结果最终的重建网格ψ_θ。

6、本发明所采用的低剂量CT影像数据来自苏州大学附属第一医院影像科，实验中使用留一法，其中每次选取57个数据作为分析此发明方法分割结果准确性和鲁棒性的比较，通过DSC系数(Dice similarity coefficient)来直观统计分割结果与金标准的重合度，即分割的准确性。结果分析如表1所示。其中所对比的其它几种分割方法如表中所列，传统的形状稀疏表示是指没有将形状变化模式整合到稀疏表示的能量函数中。其中DSC系数的定义为：

其中U₁和U₂分别表示此发明分割结果和金标准。图2是本发明方法在低剂量CT上的部分切片分割结果和金标准的对比示例图，图3列出了本发明算法的原理框图。

表1本发明方法与传统分割方法的DSC系数对比

本发明技术关键点包括：(1)将主成分分析(PCA)训练得到的形状变化模式融合到稀疏表示算法的数学模型中，以提高肺分割的准确性。(2)将可分辨的外观字典学习算法和法线方向顶点调整算法结合，以提高分割算法的效率和准确性。(3)本发明提出多方向融合的顶点调整，将法线方向和梯度矢量流场方向结合，以提高网格顶点调整的准确性和鲁棒性。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (8)

1.一种基于外观字典学习和形状稀疏表示的图像分割方法，其特征在于，包括：

步骤1，基于阈值分割和形态学操作得到待分割图像的二值图像；

步骤2，基于广义哈弗变换定位所述待分割图像的初始中心，将平均网格移动到所述待分割目标的初始中心，得到第一初始化网格；

步骤3，计算所述第一初始化网格每个顶点的法线值，根据所述待分割图像的二值图像沿着所述第一初始化网格每个顶点的法线方向调整所述第一初始化网格的每个顶点，基于形状稀疏表示算法对调整后的所述第一初始化网格进行平滑重建，得到第二初始化网格；

步骤4，计算所述第二初始网格每个顶点的法线值，计算所述第二初始网格每个顶点的法线方向外观字典重建标签的概率值，基于所述第二初始网格每个顶点的法线的矢量值结合所述法线方向外观字典重建标签的概率值局部调整所述第二初始网格每个顶点的位置，利用形状稀疏表示算法对定位后的所述第二初始网格进行平滑重建，得到第三初始网格；

步骤5，基于梯度矢量流场的定位算法依次对所述第三初始化网格上的每个顶点进行定位，利用形状稀疏表示算法对定位后的所述第三初始化网格进行平滑重建，得到最终的重建网格，将所述最终的重建网格转化为二值图像，作为最终的分割结果。

2.根据权利要求1所述的基于外观字典学习和形状稀疏表示的图像分割方法，其特征在于，步骤1中，基于阈值分割和形态学操作得到待分割图像的二值图像之前，还包括：对原始图像进行滤波去噪处理。

3.根据权利要求1所述的基于外观字典学习和形状稀疏表示的图像分割方法，其特征在于，步骤2中，将平均网格移动到所述待分割目标的初始中心，包括：基于主成分分析(PCA)训练出平均网格，所述主成分分析(PCA)的训练的公式表示如下：

其中，ψ_k表示第k个训练样本的肺网格；T表示从训练样本空间到参考样本的相似性变换；T^-1表示对应的T的反变换；p_m表示主成分分析(PCA)训练获得的网格形状变化模式；λ_m表示主要的模式所对应的权重；k_s表示训练样本的肺网格总数；表示训练得到的平均网格。

4.根据权利要求1所述的基于外观字典学习和形状稀疏表示的图像分割方法，其特征在于，步骤3中，基于形状稀疏表示算法对调整后的所述第一初始化网格进行平滑重建，包括：所述形状稀疏表示算法的定义公式如下：

其中，ψ表示法线方向调整后的所述第一初始化网格即需要被平滑的网格；D_s表示所有训练网格坐标形成的形状字典；p_m表示主成分分析(PCA)训练获得的形状变化模式；T表示从形变网格到平均网格的相似变换；T_s表示每个训练网格到平均网格的相似变换；β₁和β₂表示稀疏系数；λ₃和λ₄表示平衡参数；

当ψ和D_s相似变换到平均空间后，通过迭代最小化方式解形状稀疏表示算法的定义公式，首先使用快速的迭代阈值收缩算法(FISTA算法)优化下面的等式：

整体的形状ψ用变换到平均空间的形状训练字典进行重建，第二项拟合局部误差e，在p_m中通过参数β₂挑选合适的形变模态并设置模态的比例大小，数学上表示成如下形式：

上式用FISTA算法优化，最终的平滑后形状网格ψ表示为：

ψ＝T^-1(T_s(D_s)β₁+P_mβ₂)。

5.根据权利要求1所述的基于外观字典学习和形状稀疏表示的图像分割方法，其特征在于，步骤4中，基于所述第二初始网格每个顶点的法线的矢量值结合所述法线方向外观字典重建标签的概率值局部调整所述第二初始网格每个顶点的位置，包括：

计算平移到中心点的网格每个顶点的法线值，用n_i表示；对于网格上的每个顶点，沿着法线n_i的方向计算边界候选点在离散的位置j＝-n_m,-n_m+1,…,0,…,n_m-1,n_m，表示为如下形式：

v_ij＝v_i+jδn_i

其中，δ是搜索步长；n_t表示顶点移动的阈值；n_m表示最大的搜索范围。当顶点v_i在肺内部，也即在二值图像的前景区域，并且j＞n_t，v_ij沿着法线n_i方向搜索，顶点v_i更新为v_ij；当顶点v_i在肺外部也即在二值图像的背景区域，并且|j|＞n_t，v_ij沿着法线n_i负方向搜索，顶点v_i更新为v_ij。

6.根据权利要求1所述的基于外观字典学习和形状稀疏表示的图像分割方法，其特征在于，步骤4中，基于所述第二初始网格每个顶点的法线的矢量值结合所述法线方向外观字典重建标签的概率值局部调整所述第二初始网格每个顶点的位置，包括：所述第二初始网格顶点调整的数学公式表示为：

其中δ'是搜索步长；j'＝0,1,…,n_m-1,n_m表示沿着法线方向的候选点的离散位置；n_i'表示顶点v_i'的法线矢量值；h_t表示顶点的平移阈值；h表示局部的可分辨外观字典学习重建的标签概率值；h由如下方式得到：

训练阶段，对于被调整网格的每个顶点，根据其与训练样本对应顶点，从所有训练样本中所对应的相同领域位置取三维外观块形成块库P_l；将每个三维块表示成一维列向量的形式，因此整个三维块库P_l表示成矩阵形式：m_a表示每个块的尺寸；n_a表示P_l中块的个数；和专家标注结合，外观字典学习训练公式表示为：

其中表示学习的外观字典；α表示稀疏编码因子；k_a表示字典元素的个数；L表示P_l中每个块对应的专家标注值；ω表示学习的线性分类器；λ₁和λ₂是平衡因子；

测试阶段；外观字典在测试阶段被表示为；

其中p_t表示在形变网格附近的一个三维块；表示对输入块p_t的稀疏因子；最终输入三维块p_t中心体素的标签值被表示为；

h＝ω*α_t

特征空间的形状稀疏表示被用来对调整之后的网格进行平滑处理；整个过程迭代优化，直到两次迭代之间的网格顶点距离差Δd小于给定的阈值Δd_t；此时的输出网格ψ被当作步骤5中的输入。

7.根据权利要求1所述的基于外观字典学习和形状稀疏表示的图像分割方法，其特征在于，步骤5中，基于梯度矢量流场的定位算法依次对所述第三初始化网格上的每个顶点进行定位，包括：

一个边缘图被定义为：G＝-G_σ*B。其中G_σ是表示标准差为σ的三维高斯函数；B表示阈值分割和外观字典学习重建下的二值图像。在梯度矢量流场的方向上离散顶点搜索的公式定义为：

其中v_j-1表示形变网格上当前的顶点；表示v_j-1所对应的梯度矢量流场值的角度值；||G(v_j-1)||表示v_j-1所对应的梯度矢量流场的模值；v_j表示沿着梯度矢量流场的方向搜素的下一个目标点；表示v_j所对应的梯度矢量流场值的角度值；||G(v_j)||表示v_j所对应的梯度矢量流场的模值。满足上式的离散点形成点集合基于该集合，目标点被定义为：

其中F(·)表示在离散位置的边缘检测函数；γ表示惩罚参数。F(·)被定义如下：

其中G向三角网格的法线n_i投影限制不平衡边缘；将最大值通过sigmoid函数限制为g_max；±保证方向的一致性；表示边缘附近外观字典学习重建的标签值。特征空间的形状稀疏表示被用来对调整之后的网格进行平滑处理。整个过程迭代优化，到设计的能量函数小于某个阈值时结束。

8.根据权利要求1所述的基于外观字典学习和形状稀疏表示的图像分割方法，其特征在于，步骤5中，基于梯度矢量流场的定位算法和形状稀疏表示算法交叉迭代最小化到的最终结果，根据损失函数的设计，将损失函数表示为：

E＝E_region+λ_eE_surface其中λ_e控制区域能量项和表面能量项的平衡，对于区域能量项E_region设置为：

ψ_t表示阈值分割和外观字典学习重建的图像对应的网格；p表示梯度矢量流场上的一个体素；N_r表示归一化参数；ψ_τ表示基于梯度矢量流场定位的网格；||G_τ(p)||和||G_t(p)||分别表示ψ_τ和ψ_t在梯度矢量流场上所对应的赋值。

对于表面能量项设置为；

其中N_s表示归一化参数；n_s表示网格ψ_λ顶点总数；顶点N_v表示顶点v'的邻域点集合；ψ_λ表示形状稀疏表示平滑之后的网格。