CN110473206B - 一种基于超体素与测度学习的弥散张量图像分割方法 - Google Patents

Abstract

本发明公开了一种基于超体素与测度学习的弥散张量图像分割方法，包括以下步骤：首先，对弥散张量图像，计算描述每个体素水分子弥散的几何特征与方向特征。然后，在弥散张量图像的组织区域均匀采样种子点，结合位置、几何与方向特征，采用局部空间模糊聚类方法生成超体素。接着，在谱聚类的框架下，建立测度学习与聚类的优化模型，对目标函数采用迭代交替求解，实现超体素的分类。最后，将超体素的分类结果映射回图像空间，从而获得弥散张量图像的分割结果。本发明方法可以高效、稳定地获得精准的组织分割，对于大脑神经影像分析、疾病诊断与大脑认知研究等具有科学意义。

Description

一种基于超体素与测度学习的弥散张量图像分割方法

技术领域

本发明属于数字图像领域，尤其涉及一种基于超体素与测度学习的弥散张量图像分割方法。

背景技术

弥散张量成像(Diffusion Tensor Imaging,DTI)作为一种新型磁共振成像技术，可以在定量地获得组织体内的水分子扩散的方向与大小等信息。与传统的医学成像方式相比，弥散张量成像可以提供更加细致、独特的组织信息。因此，近年来弥散张量成像在临床诊断与分析中受到广泛的关注，特别是在神经系统疾病与脑认知研究中有着重要的作用，可以有效地检测大脑的结构与功能特性。近年来，弥散张量成像已经成功应用于多种神经系统与精神系统疾病的分析与诊断等，包括脑中风、老年痴呆症、帕金森综合症、精神分裂症、抑郁症、自闭症等多种疾病。

与传统灰度或者彩色图像有所不同，弥散张量图像通过一个张量，描述水分子在当前位置的扩散方向和大小。与传统磁共振成像相比，弥散张量成像提供了更加丰富的生物组织信息，可以有效地区分传统成像技术无法区分的生物组织，例如杏仁核、丘脑核团与胼胝体等。然而，高效、准确地分割受到了部分容积效应、噪声、复杂的张量结构以及数据的高维特性等多种因素的影响。因此，针对弥散张量图像特点，提出高效、稳定的分割算法，对于后续的疾病分析与科学研究具有重要意义。目前常用的分割方法采用逐个体素处理的方式，对于高维的弥散张量图像具有很高的时间复杂度，难以直接应用于临床。

值得注意的是，近年来，超体素作为一种新型的图像预处理技术，可以将局部空间相似特征的体素聚集生成超体素，这样可以将数量巨大的体素变成数量较少的超体素进行处理，可以在保证分割效果的同时大大降低时间复杂度。弥散张量图像的局部区域平滑，使得其适合使用超体素分割进行处理。此外，超体素技术已经被越来越多地应用于磁共振图像分析，并表现出相当好的性能，如肿瘤定位和分割、组织分割、图像配准和功能分组等。因此，有必要开展基于超体素的弥散张量图像分割研究。

在弥散张量图像分析中，高效地生成有效的超体素仍然是具有挑战性的。现有的超体素生成方法大都针对灰度或者彩色图像，难以直接应用于弥散张量图像。因此，需要针对弥散张量图像特点，研究高效、稳定的超体素生成方法。此外，对于生成的超体素如何进行有效的聚类，实现准确、稳定的组织分割，仍然受到高维、噪声等多种因素影响。

发明内容

发明目的：由于弥散张量图像体素数据大、数据维度高，使用传统的逐个体素进行处理的方式计算复杂度高，使得其在临床应用的实时性差，限制了其使用范围。本发明针对弥散张量图像的特点，提出一种基于超体素与测度学习的弥散张量图像分割方法，可以高效地获得准确的组织分割结果。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种基于超体素与测度学习的弥散张量图像分割方法，包括以下步骤：

步骤1，计算弥散张量图像的描述每个体素水分子弥散的几何特征与方向特征，进行体素张量的特征提取；

步骤2，在弥散张量图像的组织区域均匀采样种子点，结合水分子弥散的位置、几何与方向特征，采用局部空间模糊聚类方法生成超体素；

步骤3，对步骤2生成的超体素进行特征提取，使用谱聚类算法建立基于测度学习的聚类模型，通过求解模型得到超体素的聚类标签，即获得超体素的类别；

步骤4，将超体素的类别信息映射回图像空间，根据超体素聚类标签实现弥散张量图像的分割。

进一步，步骤1所述，计算弥散张量图像的描述每个体素水分子弥散的几何特征与方向特征，进行体素张量的特征提取。步骤如下：

对弥散张量进行奇异值分解，得到三个特征值λ₁,λ₂,λ₃以及对应的特征向量v₁,v₂,v₃，分别描述水分子弥散的大小与方向，特征值λ₁,λ₂,λ₃对应的方向为主方向，λ₁,λ₂,λ₃表示主方向上水分子的弥散系数，λ₁≥λ₂≥λ₃；

使用奇异值分解得到的三个特征值λ₁,λ₂,λ₃，平均弥散率MD、分数各向异性FA与容积比VR，线性各向异性(Linearity Anisoropy,CL)、平面各向异性(Planarity Anisoropy,CP)与球面各向异性(Spherical Anisoropy,CS)共九个特征进行计算；

线性各向异性CL、平面各向异性CP与球面各向异性CS，表示如下：

将上述九个特征分别归一化，合并后得到水分子弥散的几何特征f^geo，表示如下：

f^geo＝(λ₁,λ₂,λ₃,MD,FA,VR,CL,CP,CS) (4)

采用Knutsson空间对主方向进行重新定义，避免方向的二义性，计算水分子弥散的方向特征f^or，公式如下：

其中，v₁₁,v₁₂,v₁₃是特征向量v₁在三维坐标下的坐标值。

进一步，步骤2所述，在弥散张量图像的组织区域均匀采样种子点，结合水分子弥散的位置、几何与方向特征，采用局部空间模糊聚类方法生成超体素。步骤如下：

步骤2-1，在弥散张量图像的组织区域均匀采样C个种子点；

步骤2-2，采用局部空间模糊聚类方法，计算体素与距离其最近的K个种子点的隶属度，将体素分配给具有最大隶属度的种子点，生成超体素；

步骤2-3，依据生成的超体素，对每个种子点进行更新，计算每个超体素所包含的体素的位置、几何与方向特征的平均值，作为更新后的种子的特征，完成种子更新；

步骤2-4，迭代执行步骤2-2和2-3，当更新前与更新后的种子点位置变化的欧式距离小于预定的阈值时，停止更新种子点，将每个体素分配给具有最大隶属度的种子点，生成弥散张量图像的超体素。

进一步，步骤2-2所述，采用局部空间模糊聚类方法，计算体素与距离其最近的K个种子点的隶属度，将体素分配给具有最大隶属度的种子点，生成超体素；步骤如下：

步骤2-2-1，在步骤2-1采样种子点后，种子点被均匀地分布在弥散张量图像的三维空间中；将落入弥散张量图像边界的种子点移动到邻域中边缘置信度最低的位置；每个种子点用一个向量描述，定义如下：

其中，s_j表示第j个种子向量，x_j,y_j,z_j是第j个种子点的空间坐标；

与

分别是第j个种子点及其a×a×a邻域内所有体素的几何特征与方向特征的平均值；通常a取值为3；

步骤2-2-2，采用局部空间模糊聚类方法，将弥散张量图像中的每个体素与距离其最近的K个种子点进行模糊关联，局部模糊聚类的目标函数表示如下：

其中，i表示第i个体素，N是样本总数，即N个体素；j表示第j个聚类中心，即第j个种子点；K是距离当前体素最近的种子点的个数，即K个聚类中心；v_i是第i个体素的特征，s_j是第j个种子向量；D(v_i,s_j)表示体素v_i与种子s_j之间的距离；u_ij表示体素v_i与种子s_j之间的隶属度，取值范围为0到1；m为加权指数，控制聚类结果的模糊程度；其中，每个体素和距离其最近的K个种子点之间的隶属度，计算公式如下：

其中，u_it表示第i个体素与第t个种子点之间的隶属度，S是距离当前体素最近的K个种子点的集合；D(v_i,s_t)表示体素v_i与种子s_t之间的距离；m为加权指数；

对局部模糊聚类的目标函数进行迭代计算，当目标函数值小于设定的阈值，停止迭代，得到满足阈值条件下的体素与距离其最近的K个种子点的模糊关联；并计算得到每个体素和所述K个种子点之间的隶属度；

步骤2-2-3，根据体素与种子点的隶属度，将体素分配给具有最大隶属度的种子点，当每个体素都被分配到相应的种子点之后，形成超体素。

进一步，与传统的模糊聚类不同，本发明结合空间坐标、几何特征和方向特征，提出一种计算种子与体素间距离的新距离度量，将体素v_i与种子s_t之间的距离D(v_i,s_t)定义为：

D(v_i,s_t)＝d_spa(v_i,s_t)+λ_geod_geo(v_i,s_t)+λ_ord_or(v_i,s_t) (9)

其中，d_spa(v_i,s_t)表示空间距离；d_geo(v_i,s_t)表示几何特征距离；d_or(v_i,s_t)表示方向特征距离；λ_geo表示几何特征距离的权重；λ_or表示方向特征距离的权重；(x_i,y_i,z_i)表示体素v_i的空间坐标；

表示种子s_t的空间坐标；

与

分别表示体素v_i与种子s_t的几何特征；

与

分别表示体素v_i与种子s_t的方向特征。

进一步，步骤2-3所述，依据生成的超体素，对每个种子点进行更新，计算每个超体素所包含的体素的位置、几何与方向特征的平均值，作为更新后的种子的特征，完成种子更新。方法如下：

当所有体素都被分配到相应的种子点之后，对每个种子点进行更新，公式如下：

其中，N_j是与种子s_j采用式(7)进行模糊关联的所有体素的个数，m为加权指数，u_lj表示体素v_l与种子s_j的隶属度，

表示更新后的种子。

进一步，步骤2-4中，使用残差误差定义更新前后超体素的种子点位置变化的欧式距离，计算更新前种子点与更新后种子点之间的欧式距离来检测位置的变化情况；种子停止更新的阈值根据超体素的数量确定，一般为超体素的数量乘以0.1。

进一步，步骤3所述，对步骤2生成的超体素进行特征提取，使用谱聚类算法建立基于测度学习的聚类模型，通过求解模型得到超体素的聚类标签，即获得超体素的类别；步骤如下：

步骤3-1，对每个超体素进行特征提取，实现对每个超体素的描述，根据得到的每个超体素的特征计算超体素间的欧式距离；

对于每个超体素，分别计算其所包含的体素的几何特征与方向特征分量的平均值，然后进行合并，得到超体素的特征f_sv，f_sv的定义为：

其中，

是超体素的几何特征，

是超体素的方向特征；

超体素的几何特征

定义为：

超体素的方向特征

定义为：

其中，

分别表示当前超体素中所有体素的弥散张量的三个特征值λ₁,λ₂,λ₃的平均值，

分别表示当前超体素中所有体素的弥散张量的平均弥散率(MD)，分数各向异性(FA)，容积比(VR)，线性各向异性(CL)、平面各向异性(CP)与球面各向异性(CS)变量的平均值；

与

分别表示当前超体素中所有体素的代表方向特征的五个分量的平均值；

步骤3-2，构建图G＝{V,E}，V代表图中所有结点即超体素的集合，E代表图中所有边的集合；与基于体素级别的谱聚类的图像分割方法不同，本发明中的结点V是超体素，对于空间上相邻的超体素之间存在边E，边上存在权值；超体素构成图的边权值矩阵为W；为了方便后续的测度学习，引入Mahalanobis距离测度学习边上的权值，如式(17)所示：

其中，w_ij是超体素sv_i、sv_j存在边之间的权值，即边权值矩阵W的第i行第j列元素，只有相邻的超体素之间存在边；d(sv_i,sv_j)是超体素sv_i、sv_j特征之间的欧式距离，H是待学习的测度核；

为了学习测度核H，定义超体素构成图的对角阵D，其中对角元素d_ii表示如下：

d_ii＝∑_jw_ij

步骤3-3，采用半监督方式基于图进行测度学习，定义基于测度学习的聚类模型的目标优化函数为：

其中，Ψ₁(u)为平滑控制约束，表示相邻的超体素之间的变化约束；Ψ₂(u)表示提供监督信息的约束；λ是两个约束项目的权重，λ>0；u_i与u_j分别是需要求解的第i、j个超体素的聚类标签；d_ii与d_jj分别为对角阵D中的对角元素，m_i是已标注类别的弥散张量图像样本M＝{m_i}的标签值，其中m_i＝1表示胼胝体，m_i＝-1表示背景；E表示超体素之间边的集合；V表示超体素的集合；结合式(17)可知，需要学习的测度核H隐藏在权重w_ij中；

步骤3-4，随机初始化测度核H，利用内点法(Interior Point Method)对步骤3-3所述目标优化函数进行交替迭代求解，即依据测度核H，求解超体素的聚类标签u，依据聚类标签u，更新测度核H；

步骤3-5，当测度核H变化小于设定的阈值时，停止迭代过程，获得最终的超体素聚类标签u以及测度核H。

进一步，步骤3-4所述利用内点法对目标优化函数进行交替迭代求解，步骤如下：

(1)在超体素聚类标签更新阶段，依据测度核H，计算出每个超体素的聚类标签；将式(18)转换为式(19)：

其中，L是标准化拉格朗日矩阵，I表示单位矩阵，W和D分别为超体素构成图的边权值矩阵以及对角阵；M＝{m_i}表示已标注类别的弥散张量图像样本，其中m_i＝1表示胼胝体，m_i＝-1表示背景；

求解式(19)是解决一个二次规划问题，固定测度核H，通过内点法直接求解，从而获得每个超体素的聚类标签u；

(2)在测度核H更新阶段，固定聚类标签结果，通过梯度下降方式进行求解；对目标优化函数式(19)计算测度核H的梯度，如式(20)所示：

进一步对w_ij与d_ii计算测度核H的梯度，如式(21)与式(22)所示：

在上述梯度的基础上，对测度核H进行更新，如式(23)所示：

其中，H_t与H_t+1分别是更新前后的测度核；μ表示测度核更新的学习率。

有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：

本发明公开了一种基于超体素的弥散张量图像分割方法，为了验证方法的优越性，本发明方法与基于体素的方法、全卷积神经网络以及SegNet两种经典深度学习方法进行了比较。视觉评价通过显示不同方法的二维切面结果与真实结果的对比。定量评价则采用通用的三种分割评价指标，对分割结果进行了验证。此外，本发明对测度学习优化模型迭代的收敛性进行了实验验证。实验结果表明本发明提出的基于超体素的聚类的弥散张量图像分割方法可以高效、稳定地获得精准的组织分割。本发明对于大脑神经影像分析、疾病诊断与大脑认知研究等具有科学意义。

附图说明

图1是本发明的方法流程图；

图2是本发明的实施流程示意图；

图3是大脑磁共振图像上生成的超体素；

图4是不同分割方法在真实大脑弥散张量图像上分割胼胝体结果的二维切面显示。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

本发明提供一种基于超体素与测度学习的弥散张量图像分割方法，如图1和2所示，包括以下步骤：

步骤1，计算弥散张量图像的描述每个体素水分子弥散的几何特征与方向特征，进行体素张量的特征提取。步骤如下：

对弥散张量进行奇异值分解，得到三个特征值λ₁,λ₂,λ₃以及对应的特征向量v₁,v₂,v₃，分别描述水分子弥散的大小与方向，特征值λ₁,λ₂,λ₃对应的方向为主方向，λ₁,λ₂,λ₃表示主方向上水分子的弥散系数，λ₁≥λ₂≥λ₃；

使用奇异值分解得到的三个特征值λ₁,λ₂,λ₃，平均弥散率MD、分数各向异性FA与容积比VR，线性各向异性(Linearity Anisoropy,CL)、平面各向异性(Planarity Anisoropy,CP)与球面各向异性(Spherical Anisoropy,CS)共九个特征进行计算；

线性各向异性CL、平面各向异性CP与球面各向异性CS，定义为：

将上述九个特征分别归一化，合并后得到水分子弥散的几何特征f^geo，表示如下：

f^geo＝(λ₁,λ₂,λ₃,MD,FA,VR,CL,CP,CS) (4)

采用Knutsson空间对主方向进行重新定义，避免方向的二义性，计算水分子弥散的方向特征f^or，公式如下：

上式中，v₁₁,v₁₂,v₁₃是特征向量v₁在三维坐标下的坐标值。

步骤2，在弥散张量图像的组织区域均匀采样种子点，结合水分子弥散的位置、几何与方向特征，采用局部空间模糊聚类方法生成超体素。具体包括：

步骤2-1，在弥散张量图像的组织区域均匀采样C个种子点；

步骤2-2，采用局部空间模糊聚类方法，计算体素与距离其最近的K个种子点的隶属度，将体素分配给具有最大隶属度的种子点，生成超体素；

步骤2-2-1，在步骤2-1采样种子点后，种子点被均匀地分布在弥散张量图像的三维空间中；将落入弥散张量图像边界的种子点移动到邻域中边缘置信度最低的位置；每个种子点用一个向量描述，定义如下：

其中，s_j表示第j个种子向量，x_j,y_j,z_j是第j个种子点的空间坐标；

与

分别是第j个种子点及其3×3×3邻域内所有体素的几何特征与方向特征的平均值；

步骤2-2-2，采用局部空间模糊聚类方法，将弥散张量图像中的每个体素与距离其最近的K个种子点进行模糊关联，局部模糊聚类的目标函数表示如下：

其中，i表示第i个体素，N是样本总数，即N个体素；j表示第j个聚类中心，即第j个种子点；K是距离当前体素最近的种子点的个数，即K个聚类中心；v_i是第i个体素的特征，s_j是第j个种子向量；D(v_i,s_j)表示体素v_i与种子s_j之间的距离；u_ij表示体素v_i与种子s_j之间的隶属度，取值范围为0到1；m为加权指数，控制聚类结果的模糊程度；其中，每个体素和距离其最近的K个种子点之间的隶属度，计算公式如下：

其中，u_it表示第i个体素与第t个种子点之间的隶属度，S是距离当前体素最近的K个种子点的集合；D(v_i,s_t)表示体素v_i与种子s_t之间的距离；m为加权指数；

与传统的模糊聚类不同，本发明结合空间坐标、几何特征和方向特征，提出一种计算种子与体素间距离的新距离度量，将体素v_i与种子s_t之间的距离D(v_i,s_t)定义为：

D(v_i,s_t)＝d_spa(v_i,s_t)+λ_geod_geo(v_i,s_t)+λ_ord_or(v_i,s_t) (9)

其中，d_spa(v_i,s_t)表示空间距离；d_geo(v_i,s_t)表示几何特征距离；d_or(v_i,s_t)表示方向特征距离；λ_geo表示几何特征距离的权重；λ_or表示方向特征距离的权重；(x_i,y_i,z_i)表示体素v_i的空间坐标；

表示种子s_t的空间坐标；

与

分别表示体素v_i与种子s_t的几何特征；

与

分别表示体素v_i与种子s_t的方向特征；空间距离使得超体素具有更紧凑的规整尺寸，几何特征距离和方向特征距离使得超体素具有同质性；

对局部模糊聚类的目标函数进行迭代计算，当目标函数值小于设定的阈值，停止迭代，得到满足阈值条件下的体素与距离其最近的K个种子点的模糊关联；并计算得到每个体素和所述K个种子点之间的隶属度；

步骤2-2-3，根据体素与种子点的隶属度，将体素分配给具有最大隶属度的种子点，当每个体素都被分配到相应的种子点之后，形成超体素。

步骤2-3，依据生成的超体素，对每个种子点进行更新，计算每个超体素所包含的体素的位置、几何与方向特征的平均值，作为更新后的种子的特征，完成种子更新；方法如下：

当所有体素都被分配到相应的种子点之后，对每个种子点进行更新，公式如下：

其中，N_j是与种子s_j采用式(7)进行模糊关联的所有体素的个数，m为加权指数，u_lj表示体素v_l与种子s_j的隶属度，

表示更新后的种子。

步骤2-4，迭代执行步骤2-2和2-3，当更新前与更新后的种子点位置变化的欧式距离小于预定的阈值时，停止更新种子点，将每个体素分配给具有最大隶属度的种子点，生成弥散张量图像的超体素。

使用残差误差定义更新前后超体素的种子点位置变化的欧式距离，计算更新前种子点与更新后种子点之间的欧式距离来检测位置的变化情况；种子停止更新的阈值为超体素的数量乘以0.1。

步骤3，对步骤2生成的超体素进行特征提取，使用谱聚类算法建立基于测度学习的聚类模型，通过求解模型得到超体素的聚类标签，即获得超体素的类别；步骤如下：

步骤3-1，对每个超体素进行特征提取，实现对每个超体素的描述，根据得到的每个超体素的特征计算超体素间的欧式距离；

对于每个超体素，分别计算其所包含的体素的几何特征与方向特征分量的平均值，然后进行合并，得到超体素的特征f_sv，f_sv的定义为：

其中，

是超体素的几何特征，

是超体素的方向特征；

超体素的几何特征

定义为：

超体素的方向特征

定义为：

其中，

分别表示当前超体素中所有体素的弥散张量的三个特征值λ₁,λ₂,λ₃的平均值，

分别表示当前超体素中所有体素的弥散张量的平均弥散率(MD)，分数各向异性(FA)，容积比(VR)，线性各向异性(CL)、平面各向异性(CP)与球面各向异性(CS)变量的平均值；

与

分别表示当前超体素中所有体素的代表方向特征的五个分量的平均值；

步骤3-2，构建图G＝{V,E}，V代表图中所有结点即超体素的集合，E代表图中所有边的集合；与基于体素级别的谱聚类的图像分割方法不同，本发明中的结点V是超体素，对于空间上相邻的超体素之间存在边E，边上存在权值；超体素构成图的边权值矩阵为W；为了方便后续的测度学习，引入Mahalanobis距离测度学习边上的权值，如式(17)所示：

其中，w_ij是超体素sv_i、sv_j存在边之间的权值，即边权值矩阵W的第i行第j列元素，只有相邻的超体素之间存在边；d(sv_i,sv_j)是超体素sv_i、sv_j特征之间的欧式距离，H是待学习的测度核；

为了学习测度核H，定义超体素构成图的对角阵D，其中对角元素d_ii表示如下：

d_ii＝∑_jw_ij

步骤3-3，采用半监督方式基于图进行测度学习，定义基于测度学习的聚类模型的目标优化函数为：

其中，Ψ₁(u)为平滑控制约束，表示相邻的超体素之间的变化约束；Ψ₂(u)表示提供监督信息的约束；λ是两个约束项目的权重，λ>0；u_i与u_j分别是需要求解的第i、j个超体素的聚类标签；d_ii与d_jj分别为对角阵D中的对角元素，m_i是已标注类别的弥散张量图像样本M＝{m_i}的标签值，其中m_i＝1表示胼胝体，m_i＝-1表示背景；E表示超体素之间边的集合；V表示超体素的集合；结合式(17)可知，需要学习的测度核H隐藏在权重w_ij中；

步骤3-4，随机初始化测度核H，利用内点法(Interior Point Method)对步骤3-3所述目标优化函数进行交替迭代求解，即依据测度核H，求解超体素的聚类标签u，依据聚类标签u，更新测度核H；步骤如下：

(1)在超体素聚类标签更新阶段，依据测度核H，计算出每个超体素的聚类标签；将式(18)转换为式(19)：

其中，L是标准化拉格朗日矩阵，I表示单位矩阵，W和D分别为超体素构成图的边权值矩阵以及对角阵；M＝{m_i}表示已标注类别的弥散张量图像样本，其中m_i＝1表示胼胝体，m_i＝-1表示背景；

求解式(19)是解决一个二次规划问题，固定测度核H，通过内点法直接求解，从而获得每个超体素的聚类标签u；

(2)在测度核H更新阶段，固定聚类标签结果，通过梯度下降方式进行求解；对目标优化函数式(19)计算测度核H的梯度，如式(20)所示：

进一步对w_ij与d_ii计算测度核H的梯度，如式(21)与式(22)所示：

在上述梯度的基础上，对测度核H进行更新，如式(23)所示：

其中，H_t与H_t+1分别是更新前后的测度核；μ表示测度核更新的学习率。

步骤3-5，当测度核H变化小于设定的阈值时，停止迭代过程，获得最终的超体素聚类标签u以及测度核H。

步骤4，将超体素的类别信息映射回图像空间，根据超体素聚类标签实现弥散张量图像的分割。

本实施例利用从人类连接体项目数据库中选择的20名健康受试者的大脑弥散张量图像数据集。数据是在3T磁共振扫描仪上使用加速多波段成像协议获取的。成像参数如下，回波时间89.5毫秒，重复时间5520毫秒，视野210毫米×210毫米。共145片，片厚1.25mm。b值分别为1000、2000和3000s/mm2，共获得扩散加权扫描270次，b值为0的扫描18次，作为基线图像。首先对扩散加权图像进行预处理，然后利用FMRIB Software Library(FSL)软件评估弥散张量。所有图像均以胼胝体的人工分割为金标准分割。人工分割由两位神经放射学专家进行和评估。

采用以上相同的真实脑部弥散张量数据集进行分割验证。数据集中包含了20个来自健康人的大脑弥散张量图像，并且提供了胼胝体的分割结果。本实施例同时通过视觉评价与定量评价，对提出方法进行验证。为了验证方法的优越性，本实施例提出的方法与基于体素的方法、基于超体素的图学习方法、两种经典深度学习方法进行了比较。

图3中的(a)是弥散张量图像，图3中的(b)是胼胝体的分割结果，图3中的(c)、(d)、(e)、(f)分别是形状保持超体素方法(RSV)、基于图的聚类算法(GB)和简单线性迭代聚类算法(SLIC)与本方法的分割结果；图4中的(a)是真实分割结果，图4中的(b)、(c)、(d)、(e)分别是基于体素方法、全卷积神经网络、SegNet与本方法分割结果。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (7)

1.一种基于超体素与测度学习的弥散张量图像分割方法，其特征在于：该方法包括以下步骤：

步骤1，计算弥散张量图像的描述每个体素水分子弥散的几何特征与方向特征，进行体素张量的特征提取；

步骤2，在弥散张量图像的组织区域均匀采样种子点，结合水分子弥散的位置、几何与方向特征，采用局部空间模糊聚类方法生成超体素；

步骤3，对步骤2生成的超体素进行特征提取，使用谱聚类算法建立基于测度学习的聚类模型，通过求解模型得到超体素的聚类标签，即获得超体素的类别；

步骤3-1，对每个超体素进行特征提取，根据得到的每个超体素的特征计算超体素间的欧式距离；

对于每个超体素，分别计算其所包含的体素的几何特征与方向特征分量的平均值，然后进行合并，得到超体素的特征f_sv，f_sv的定义为：

其中，

是超体素的几何特征，

是超体素的方向特征；

超体素的几何特征

定义为：

超体素的方向特征

定义为：

其中，

分别表示当前超体素中所有体素的弥散张量的三个特征值λ₁,λ₂,λ₃的平均值，

分别表示当前超体素中所有体素的弥散张量的平均弥散率MD，分数各向异性FA，容积比VR，线性各向异性CL、平面各向异性CP与球面各向异性CS变量的平均值；

与

分别表示当前超体素中所有体素的代表方向特征的五个分量的平均值；

步骤3-2，构建图G＝{V,E}，V代表图中所有结点即超体素的集合，E代表图中所有边的集合；对于空间上相邻的超体素之间存在边E，边上存在权值；超体素构成图的边权值矩阵W的元素表示如下：

其中，w_ij是超体素sv_i、sv_j存在边之间的权值，即边权值矩阵W的第i行第j列元素；d(sv_i,sv_j)是超体素sv_i、sv_j特征之间的欧式距离，H是待学习的测度核；

定义超体素构成图的对角阵D，其中对角元素d_ii表示如下：

d_ii＝∑_jw_ij

步骤3-3，采用半监督方式基于图进行测度学习，定义基于测度学习的聚类模型的目标优化函数为：

其中，Ψ₁(u)为平滑控制约束，表示相邻的超体素之间的变化约束；Ψ₂(u)表示提供监督信息的约束；λ是两个约束项目的权重，λ>0；u_i与u_j分别是需要求解的第i、j个超体素的聚类标签；d_ii与d_jj分别为对角阵D中的对角元素，m_i是已标注类别的弥散张量图像样本M＝{m_i}的标签值，其中m_i＝1表示胼胝体，m_i＝-1表示背景；E表示超体素之间边的集合；V表示超体素的集合；

步骤3-4，随机初始化测度核H，利用内点法对步骤3-3所述目标优化函数进行交替迭代求解，即依据测度核H，求解超体素的聚类标签u，依据聚类标签u，更新测度核H；具体如下：

(1)在超体素聚类标签更新阶段，依据测度核H，计算出每个超体素的聚类标签；将式(18)转换为式(19)：

其中，L是标准化拉格朗日矩阵，I表示单位矩阵，W和D分别为超体素构成图的边权值矩阵以及对角阵；M＝{m_i}表示已标注类别的弥散张量图像样本，其中m_i＝1表示胼胝体，m_i＝-1表示背景；

固定测度核H，通过内点法直接求解式(19)，从而获得每个超体素的聚类标签u；

(2)在测度核H更新阶段，固定聚类标签结果，通过梯度下降方式进行求解；对目标优化函数式(19)计算测度核H的梯度，如式(20)所示：

对w_ij与d_ii计算测度核H的梯度，如式(21)与式(22)所示：

在上述梯度的基础上，对测度核H进行更新，如式(23)所示：

其中，H_t与H_t+1分别是更新前后的测度核；μ表示测度核更新的学习率；

步骤3-5，当测度核H变化小于设定的阈值时，停止迭代过程，获得最终的超体素聚类标签u以及测度核H；

步骤4，将超体素的类别信息映射回图像空间，根据超体素聚类标签实现弥散张量图像的分割。

2.根据权利要求1所述的一种基于超体素与测度学习的弥散张量图像分割方法，其特征在于：步骤1所述，计算弥散张量图像的描述每个体素水分子弥散的几何特征与方向特征，进行体素张量的特征提取，步骤如下：

对弥散张量进行奇异值分解，得到三个特征值λ₁,λ₂,λ₃以及对应的特征向量v₁,v₂,v₃，分别描述水分子弥散的大小与方向，特征值λ₁,λ₂,λ₃对应的方向为主方向，λ₁,λ₂,λ₃表示主方向上水分子的弥散系数，λ₁≥λ₂≥λ₃；

使用奇异值分解得到的三个特征值λ₁,λ₂,λ₃，平均弥散率MD、分数各向异性FA与容积比VR，线性各向异性CL、平面各向异性CP与球面各向异性CS共九个特征进行计算；

线性各向异性CL、平面各向异性CP与球面各向异性CS，表示如下：

将上述九个特征分别归一化，合并后得到水分子弥散的几何特征f^geo，表示如下：

f^geo＝(λ₁,λ₂,λ₃,MD,FA,VR,CL,CP,CS) (4)

采用Knutsson空间对主方向进行重新定义，计算水分子弥散的方向特征f^or，公式如下：

其中，v₁₁,v₁₂,v₁₃是特征向量v₁在三维坐标下的坐标值。

3.根据权利要求2所述的一种基于超体素与测度学习的弥散张量图像分割方法，其特征在于：步骤2所述，在弥散张量图像的组织区域均匀采样种子点，结合水分子弥散的位置、几何与方向特征，采用局部空间模糊聚类方法生成超体素；步骤如下：

步骤2-1，在弥散张量图像的组织区域均匀采样C个种子点；

步骤2-2，采用局部空间模糊聚类方法，计算体素与距离其最近的K个种子点的隶属度，将体素分配给具有最大隶属度的种子点，生成超体素；

步骤2-3，依据生成的超体素，对每个种子点进行更新，计算每个超体素所包含的体素的位置、几何与方向特征的平均值，作为更新后的种子的特征，完成种子更新；

步骤2-4，迭代执行步骤2-2和2-3，当更新前与更新后的种子点位置变化的欧式距离小于预定的阈值时，停止更新种子点，将每个体素分配给具有最大隶属度的种子点，生成弥散张量图像的超体素。

4.根据权利要求3所述的一种基于超体素与测度学习的弥散张量图像分割方法，其特征在于：步骤2-2所述，采用局部空间模糊聚类方法，计算体素与距离其最近的K个种子点的隶属度，将体素分配给具有最大隶属度的种子点，生成超体素；步骤如下：

步骤2-2-1，在步骤2-1采样种子点后，种子点被均匀地分布在弥散张量图像的三维空间中；将落入弥散张量图像边界的种子点移动到邻域中边缘置信度最低的位置；每个种子点用一个向量描述，定义如下：

其中，s_j表示第j个种子向量，x_j,y_j,z_j是第j个种子点的空间坐标；

与

分别是第j个种子点及其a×a×a邻域内所有体素的几何特征与方向特征的平均值；

步骤2-2-2，采用局部空间模糊聚类方法，将弥散张量图像中的每个体素与距离其最近的K个种子点进行模糊关联，局部模糊聚类的目标函数表示如下：

其中，i表示第i个体素，N是样本总数，即N个体素；j表示第j个聚类中心，即第j个种子点；K是距离当前体素最近的种子点的个数，即K个聚类中心；v_i是第i个体素的特征，s_j是第j个种子向量；D(v_i,s_j)表示体素v_i与种子s_j之间的距离；u_ij表示体素v_i与种子s_j之间的隶属度；m为加权指数；其中，每个体素和距离其最近的K个种子点之间的隶属度，计算公式如下：

其中，u_it表示第i个体素与第t个种子点之间的隶属度，S是距离当前体素最近的K个种子点的集合；D(v_i,s_t)表示体素v_i与种子s_t之间的距离；m为加权指数；

对局部模糊聚类的目标函数进行迭代计算，当目标函数值小于设定的阈值，停止迭代，得到满足阈值条件下的体素与距离其最近的K个种子点的模糊关联；并计算得到每个体素和所述K个种子点之间的隶属度；

步骤2-2-3，根据体素与种子点的隶属度，将体素分配给具有最大隶属度的种子点，当每个体素都被分配到相应的种子点之后，形成超体素。

5.根据权利要求4所述的一种基于超体素与测度学习的弥散张量图像分割方法，其特征在于：结合空间坐标、几何特征和方向特征，将体素v_i与种子s_t之间的距离D(v_i,s_t)定义为：

D(v_i,s_t)＝d_spa(v_i,s_t)+λ_geod_geo(v_i,s_t)+λ_ord_or(v_i,s_t) (9)

其中，d_spa(v_i,s_t)表示空间距离；d_geo(v_i,s_t)表示几何特征距离；d_or(v_i,s_t)表示方向特征距离；λ_geo表示几何特征距离的权重；λ_or表示方向特征距离的权重；(x_i,y_i,z_i)表示体素v_i的空间坐标；

表示种子s_t的空间坐标；

与

分别表示体素v_i与种子s_t的几何特征；

与

分别表示体素v_i与种子s_t的方向特征。

6.根据权利要求5所述的一种基于超体素与测度学习的弥散张量图像分割方法，其特征在于：步骤2-3所述，依据生成的超体素，对每个种子点进行更新，计算每个超体素所包含的体素的位置、几何与方向特征的平均值，作为更新后的种子的特征，完成种子更新；方法如下：

当所有体素都被分配到相应的种子点之后，对每个种子点进行更新，公式如下：

其中，N_j是与种子s_j模糊关联的所有体素的个数，m为加权指数，u_lj表示体素v_l与种子s_j的隶属度，

表示更新后的种子。

7.根据权利要求6所述的一种基于超体素与测度学习的弥散张量图像分割方法，其特征在于：步骤2-4中，使用残差误差定义更新前后超体素的种子点位置变化的欧式距离，计算更新前种子点与更新后种子点之间的欧式距离来检测位置的变化情况；种子停止更新的阈值根据超体素的数量确定。

Images