一种多模态脑影像损伤性病变组织图像分割方法与流程

本发明属于计算机视觉、模式识别、图像处理与分析技术领域，具体涉及一种多模态脑影像损伤性病变组织图像分割方法。

背景技术：

人类大脑是历经漫长自然进化产生的人类区别于其它生物表现出高级智能的重要人体器官。人脑组织结构复杂，对应病变(如脑卒中)具有高发病率、高复发率、高致残率、高致死率的特点，已成为危害人类健康的重要疾病。数据显示，脑卒中已成为世界范围内人类第二大死因。在我国，脑部病变严重影响人民身体健康与生活质量。《柳叶刀》杂志近期公布的《全球疾病负担报告2013》显示，脑卒中位列我国所有死亡原因第一位。中国每年新增脑卒中患者约200万人，约70％～80％不能独立生活。因此，早期发现，并准确诊断人脑损伤性病变，进而采取有效的预防措施或制定行之有效的治疗方案，对降低人脑疾病给人类健康带来的巨大威胁具有重要意义。

随着现代医学成像技术的发展，多种模态的影像学设备广泛应用于脑部疾病筛查。磁共振成像(magnetic resonance imaging，MRI)具有高分辨率、无射线辐射、不需注射造影剂、不受骨组织干扰等优点，其多参数性、多光谱特性能够对脑损伤病变从多个角度给予影像学支持，是显示颅内脑解剖结构形态变化的理想影像学手段。临床上，包括扩散加权成像(diffusion weighted imaging，DWI)、fluid attenuated inversion recovery(FLAIR)、T1-weighted(T1w)MR、T2-weighted(T2w)MR在内的磁共振成像方法能有效用于人脑病变的定位、诊断，并可追踪病变演化过程。充分利用这些多模态影像数据进行脑组织结构的定量提取与分析，有利于脑病变组织的定位识别、手术切除、术后的康复治疗，以及病变发展趋势的估计与预防。如图1所示为一个亚急性期缺血性脑卒中病变影像灰度特征的例子，其中，(a)为T1模态影像图像，图中右下角为病变位置的影像图像，(b)为T2模态影像图像，(c)为FLAIR模态影像图像，(d)为DWI模态影像图像，不难看出，同一病变在不同模态的MR影像中表现出明显不同，与周围健康脑组织的对比度亦有差异。同时，健康脑组织(白质、灰质、脑脊液)在不同模态的MR影像中表现也不尽相同。

综上，人脑损伤性病变组织在不同模态的MR影像中表现出不同的影像学特征。这些特征给人脑疾病诊断、病变定位、病变分析提供了丰富的影像学信息。如何有效地融合利用各模态的有价值信息，快速准确地提取人脑病变组织，分析其变化趋势，进而采取有效的预防措施，制定行之有效的治疗方案，对降低人脑病变给人类健康带来的巨大威胁具有重要意义。

为定位人脑病变位置，分析病变情况，全自动地识别人脑病变组织，将其从多模态影像中分割出来成为关键问题。手动分割简单直接，但存在如下缺点：1)手动勾勒分割曲线乏味耗时；2)分割往往因个人见解不同而存在倾向性差异；3)分割过程和结果均不可完整重现；4)逐层分割，未考虑人脑空间解剖结构信息。借助计算机技术，自动或半自动的脑损伤性病变组织分割方法能有效克服手动分割的缺点。但是，基于多模态影像准确快速地分割脑损伤性病变组织，仍存在不同模态影像空间一致性差、各模态影像质量参差不齐、算法复杂度高运算速度慢等问题。除此之外，现有的用于分割脑组织的方法较少融合脑影像的多模态特征，且相关研究成果临床应用仍不成熟。

技术实现要素：

针对现有技术的不足，本发明提出一种多模态脑影像损伤性病变组织图像分割方法。

本发明的技术方案是：

一种多模态脑影像损伤性病变组织图像分割方法，包括以下步骤：

步骤1：获取L个模态下包含待分割损伤性病变组织的脑影像图像；

步骤2：对步骤1获取的各模态脑影像图像进行预处理，即非脑组织剔除和影像空间一致性匹配；

步骤3：采用加权Fuzzy C-means方法依次分割各模态预处理后的脑影像图像，得到各模态初步分割后的脑影像图像；

步骤4：采用多数投票法将各模态初步分割后的脑影像图像中脑损伤性病变组织进行标签融合，即影像像素点x在大于等于L/2的模态中被步骤3中的分割结果认定为脑损伤性病变组织，则标定该影像像素点x为脑损伤性病变组织，得到各模态标签融合后的粗分割结果图像；

步骤5：将步骤4中得到的各模态标签融合后的粗分割结果图像作为初始输入，采用水平集分割方法对多模态脑影像图像进行分割，得到脑损伤性病变组织脑影像图像分割结果。

所述步骤2包括以下步骤：

(A)：采用BET算法依次去除各模态脑影像图像中的非脑组织；

(B)：采用弹性配准方法将去除非脑组织的各模态脑影像图像进行空间一致性匹配。

所述步骤2包括以下步骤：

(a)：采用弹性配准方法将各模态脑影像图像进行空间一致性匹配；

(b)：采用BET算法去除匹配后的脑影像图像中的非脑组织。

所述步骤3包括以下步骤：

步骤3.1：建立各模态预处理后的脑影像图像的加权Fuzzy C-means能量函数

所述加权Fuzzy C-means能量函数如下所示：

其中，||*||为欧式距离运算符，为各模态预处理后的脑影像图像的灰度值的连续定义域，即I_i：Ω→R，m为各模态预处理后的脑影像图像的维度，I_i为预处理后的第i个模态的脑影像图像，i＝1，...，L，j＝1，...，N，N为影像图像中包括影像背景图像在内的脑组织图像区域个数，λ_ij为聚类类别权重系数，q为模糊度控制参数，I_i(x)≈c_ijfor x∈Ω_j，c_i＝(c_i1，c_i2，...，c_iN)^T为能量函数的聚簇中心，c_ij为理想情况下第i个模态的第j类脑组织的灰度常数，Ω_j为第j类脑组织的影像区域，x为影像像素点，u_i＝(u_i1，u_i2，...，u_iN)^T为隶属度函数，u_ij(x)为第i个模态的第j类脑组织的隶属度函数，且满足u_ij(x)∈[0，1]；

步骤3.2：设定聚类类别权重系数λ_ij＞0、模糊度控制参数q≥1、能量函数最大迭代次数η、能量函数收敛阈值ε；

步骤3.3：采用随机数初始化能量函数的聚簇中心c_i和隶属度函数u_i；

步骤3.4：令当前隶属度函数u_i不变，以加权Fuzzy C-means能量函数取得最小值为原则，迭代更新能量函数的聚簇中心c_i；

步骤3.5：令当前聚簇中心c_i不变，以加权Fuzzy C-means能量函数取得最小值为原则，迭代更新能量函数的隶属度函数u_i；

步骤3.6：将当前隶属度函数u_i和当前聚簇中心c_i代入加权Fuzzy C-means能量函数后得到当前能量函数值；

步骤3.7：若当前相邻两次迭代的能量函数值的差值小于等于能量函数收敛阈值ε或迭代次数达到能量函数最大迭代次数η，则执行步骤3.8，否则，迭代次数加1，返回步骤3.4；

步骤3.8：将各模态脑影像图像的影像像素点x归类为隶属度函数u_ij(x)取得最大值的类别，得到各模态初步分割后的脑影像图像。

所述步骤5包括以下步骤：

步骤5.1：建立多模态脑影像图像的水平集分割能量泛函将脑影像图像分割为N′＝3个区域，所述分割的N′＝3个区域分别为：损伤性病变组织区域、正常脑组织区域和影像背景区域；

所述水平集分割的能量泛函如下所示：

其中，为图像数据项，φ₁和φ₂为水平集函b＝(b₁，b₂，...，b_L)^T为偏移场列向量，c＝(c₁，c₂，...，c_L)为影像灰度常数矩阵，K为任意归一化的非负偶函数，且在区间[0，+∞)上为单调降函数，j′＝1，...，N′，γ_j′为类别控制参数，χ_i为模态控制参数，y为给定像素点，当一个半径为ρ的圆形邻域足够小时，b_i(x)≈b_i(y)，b_i为第i个模态的脑影像图像的偏移场，c_ij′为第i个模态的第j′个待分割的灰度常数，

为水平集函数正则项，μ₁为控制水平集函数φ₁的正则项控制参数，μ₂为控制水平集函数φ₂的正则项控制参数，为水平集函数φ₁的梯度函数，为水平集函数φ₂的梯度函数；

为水平集分割曲线平滑项，v₁为控制水平集函数φ₁的分割曲线的弧长控制参数，v₂为控制水平集函数φ₂的分割曲线的弧长控制参数，δ(φ₁(x))为水平集函数φ₁的Dirac函数，δ(φ₂(x))为水平集函数φ₂的Dirac函数；

步骤5.2：以标签融合后的粗分割结果图像中脑损伤性病变组织区域边界为零水平集初始化水平集函数φ₁为符号距离函数，以标签融合后的粗分割结果图像中损伤性病变组织区域与正常脑组织区域并集的边界为零水平集初始化水平集函数φ₂为符号距离函数，采用随机数初始化偏移场列向量b；

步骤5.3：以初始化的水平集函数φ₁和φ₂构造隶属度函数：待分割脑损伤性病变组织隶属度函数M₁(φ₁，φ₂)＝(1-H(φ₁))(1-H(φ₂))，待分割正常脑组织区域隶属度函数M₂(φ₁，φ₂)＝H(φ₁)(1-H(φ₂))，待分割影像背景区域隶属度函数M₃(φ₁，φ₂)＝H(φ₂)，其中，H为Heaviside函数；

步骤5.4：根据各模态标签融合后的粗分割结果图像的特点和初始的水平集函数值确定水平集分割控制参数：类别控制参数γ_j′、模态控制参数χ_i、正则项控制参数μ₁和μ₂、分割曲线的弧长控制参数v₁和v₂、领域半径ρ；设定水平集分割最大迭代次数ξ和能量泛函收敛阈值κ；

步骤5.5：令当前水平集函数φ₁、当前水平集函数φ₂和当前偏移场列向量b不变，以能量泛函取得最小值为原则，迭代更新影像灰度常数矩阵c；

步骤5.6：令当前水平集函数φ₁、当前水平集函数φ₂和当前影像灰度常数矩阵c不变，以能量泛函取得最小值为原则，迭代更新偏移场列向量b；

步骤5.7：根据当前影像灰度常数矩阵c和当前偏移场列向量b采用最大梯度降方法迭代更新水平集函数φ₁和水平集函数φ₂；

步骤5.8：根据当前水平集函数φ₁、当前水平集函数φ₂、当前影像灰度常数矩阵c和当前偏移场列向量b计算当前水平集分割的能量泛函函数值；

步骤5.9：若当前相邻两次迭代的水平集分割的能量泛函函数值的差值小于等于能量泛函收敛阈值κ或者当前迭代次数达到水平集分割最大迭代次数ξ，则执行步骤5.10，否则，迭代次数加1，返回步骤5.5；

步骤5.10：根据当前水平集函数φ₁和当前水平集函数φ₂构造脑损伤性病变组织隶属度函数M₁(φ₁，φ₂)＝(1-H(φ₁))(1-H(φ₂))，即为脑损伤性病变组织脑影像图像分割结果。

本发明的有益效果：

本发明提出一种多模态脑影像损伤性病变组织图像分割方法，该方法能够有效融合多个影像模态的有益信息，克服了弱边界、图像噪声、以及灰度不一致性(偏移场)对图像分割精度的负面影响，能够得到精确的脑损伤性病变组织的分割结果。

附图说明

图1为亚急性期缺血性脑卒中损伤性病变的不同模态MR影像图像；

其中，(a)为T1模态影像图像，(b)为T2模态影像图像，(c)为FLAIR模态影像图像，(d)为DWI模态影像图像；

图2为本发明具体实施方式中多模态脑影像损伤性病变组织图像分割方法的流程图；

图3为本发明方法一进行4种模态脑影像图像进行预处理的示意图；

其中，(a)为4种模态脑影像图像，(b)为4种模态去除非脑组织后的影像图像，(c)为4种模态预处理后的影像图像；

图4为本发明方法二进行4种模态脑影像图像进行预处理的示意图；

其中，(a)为4种模态脑影像图像，(b)为4种模态空间一致性匹配后的影像图像，(c)为4种模态预处理后的影像图像；

图5为本发明具体实施方式中采用加权Fuzzy C-means方法依次分割各模态预处理后的脑影像图像的流程图；

图6为本发明具体实施方式中采用水平集分割方法对多模态脑影像图像进行分割的流程图；

图7为本发明实施例一中四个模态下包含待分割损伤性病变组织的脑影像图像；

其中，(a)为DWI模态影像图像，(b)为FLAIR模态影像图像，(c)为T1模态影像图像，(d)为T2模态影像图像；

图8为本发明实施例一中四个模态加权Fuzzy C-Means分割后的脑影像图像；

其中，(a)为DWI模态加权Fuzzy C-Means分割后的影像图像，(b)为FLAIR模态加权Fuzzy C-Means分割后的影像图像，(c)为T1模态加权Fuzzy C-Means分割后的影像图像，(d)为T2模态加权Fuzzy C-Means分割后的影像图像；

图9为本发明实施例一中标签融合后的粗分割结果图像和多模态脑影像水平集方法分割脑损伤性病变组织脑影像图像的分割结果；

其中，(a)为标签融合后的粗分割结果图像，(b)为脑损伤性病变组织脑影像图像分割结果；

图10为本发明实施例一中已知的分割标准；

图11为本发明实施例二中四个模态下包含待分割损伤性病变组织的脑影像图像；

其中，(a)为DWI模态影像图像，(b)为FLAIR模态影像图像，(c)为T1模态影像图像，(d)为T2模态影像图像；

图12为本发明实施例二中四个模态加权Fuzzy C-Means分割后的脑影像图像；

其中，(a)为DWI模态加权Fuzzy C-Means分割后的影像图像，(b)为FLAIR模态加权Fuzzy C-Means分割后的影像图像，(c)为T1模态加权Fuzzy C-Means分割后的影像图像，(d)为T2模态加权Fuzzy C-Means分割后的影像图像；

图13为本发明实施例二中标签融合后的粗分割结果图像和多模态脑影像水平集方法分割脑损伤性病变组织脑影像图像的分割结果；

其中，(a)为标签融合后的粗分割结果图像，(b)为脑损伤性病变组织脑影像图像分割结果；

图14为本发明实施例二中已知的分割标准；

图15为本发明实施例三中四个模态下包含待分割损伤性病变组织的脑影像图像；

其中，(a)为DWI模态影像图像，(b)为FLAIR模态影像图像，(c)为T1模态影像图像，(d)为T2模态影像图像；

图16为本发明实施例三中四个模态加权Fuzzy C-Means分割后的脑影像图像；

其中，(a)为DWI模态加权Fuzzy C-Means分割后的影像图像，(b)为FLAIR模态加权Fuzzy C-Means分割后的影像图像，(c)为T1模态加权Fuzzy C-Means分割后的影像图像，(d)为T2模态加权Fuzzy C-Means分割后的影像图像；

图17为本发明实施例三中标签融合后的粗分割结果图像和多模态脑影像水平集方法分割脑损伤性病变组织脑影像图像的分割结果；

其中，(a)为标签融合后的粗分割结果图像，(b)为脑损伤性病变组织脑影像图像分割结果；

图18为本发明实施例三中已知的分割标准。

具体实施方式

下面结合附图对本发明具体实施方式加以详细说明。

本发明提出一种多模态脑影像损伤性病变组织图像分割方法，如图2所示，包括以下步骤：

步骤1：获取L个模态下包含待分割损伤性病变组织的脑影像图像。

本实施方式中，获取了L＝4个模态下包含待分割损伤性病变组织的脑影像图像，如图1所示。

步骤2：对步骤1获取的各模态脑影像图像进行预处理，即非脑组织剔除和影像空间一致性匹配。

本实施方式中，对步骤1获取的各模态脑影像图像进行预处理有两种处理方法，方法一如图3所示，包括以下步骤：

(A)：采用BET算法依次去除各模态脑影像图像中的非脑组织。

(B)：采用弹性配准方法将去除非脑组织的各模态脑影像图像进行空间一致性匹配。

方法二如图4所示，包括以下步骤：

(a)：采用弹性配准方法将各模态脑影像图像进行空间一致性匹配。

(b)：采用BET算法去除匹配后的脑影像图像中的非脑组织。

本实施方式中，该方法不需要将剔除算法分别应用于各个模态，只需要在空间一致性的单一模态中剔除非脑组织，再将剔除掉的组织从其它模态的影像中去除。

步骤3：采用加权Fuzzy C-means方法依次分割各模态预处理后的脑影像图像，得到各模态初步分割后的脑影像图像，如图5所示。

步骤3.1：建立各模态预处理后的脑影像图像的加权Fuzzy C-means能量函数

本实施方式中，给定第i个模态的预处理后的脑影像图像I_i，其中i＝1，...，L，其灰度值可以看作连续定义域上的实数映射，即I_i：Ω→R，m为各模态预处理后的脑影像图像的维度，预处理后的脑影像图像I_i的灰度值为人脑组织的内在物理特性在影像学中的表现，也就是说，理想情况下，影像背景图像灰度与N-1种脑组织(白质、灰质、脑脊液、脑损伤性病变组织，即N＝5为影像图像中包括影像背景图像在内的脑组织图像区域个数)取不同的常数灰度值c_ij，其中，j＝1，...，N。得出如式(1)所示公式：

I_i(x)≈c_ijfor x∈Ω_j (1)

其中，c_ij为理想情况下第i个模态的第j类脑组织的灰度常数，Ω_j为第j类脑组织的影像区域，x为影像像素点，满足且m≠n时m＝1，...，N，n＝1，...，N，影像灰度集I_ij＝{I_i(x)：x∈Ω_j}构成一个以c_ij为中心的聚簇。这一聚簇属性表明定义域Ω内的影像灰度能被划分为N类。为了分类具有以上性质的影像灰度，建立各模态脑影像预处理后的影像图像的加权Fuzzy C-means能量函数如式(2)所示：

其中，||*||是欧式距离运算符，λ_ij为聚类类别权重系数，q为模糊度控制参数，c_i＝(c_i1，c_i2，...，c_iN)^T为能量函数的聚簇中心，u_i＝(u_i1，u_i2，...，u_iN)^T为隶属度函数，u_ij(x)为第i个模态的第j类脑组织的隶属度函数，且满足u_ij(x)∈[0，1]。

步骤3.2：设定聚类类别权重系数λ_ij＞0、模糊度控制参数q≥1、能量函数最大迭代次数η、能量函数收敛阈值ε。

本实施方式中，设定的聚类类别权重系数λ_ij为1.0，模糊度控制参数q为2.0，能量函数最大迭代次数η为20，能量函数收敛阈值ε为0.01。

步骤3.3：采用随机数初始化能量函数的聚簇中心c_i和隶属度函数u_i。

本实施方式中，初始化的能量函数的聚簇中心c_i＝(c_i1，c_i2，...，c_iN)^T，其中c_ij＝255/N×j，j＝1，2，...，N，初始化的隶属度函数u_i＝(u_i1，u_i2，...，u_iN)^T均为1.0。

步骤3.4：令当前隶属度函数u_i不变，以加权Fuzzy C-means能量函数取得最小值为原则，迭代更新能量函数的聚簇中心c_i。

本实施方式中，当与各聚簇中心相近的影像灰度的隶属度函数值较大，与各聚簇中心差距较大的影像灰度的隶属度函数值较小时，影像图像的加权Fuzzy C-means能量函数取得最小值，反之亦然。因此，为使得取得最小值，保持隶属度函数u_i不变，解其中0为元素值均为零的常数列向量。最终得存c_i＝(c_i1，c_i2，...，c_iN)^T处取得最小值，因此，迭代更新能量函数的聚簇中心c_i的公式如式(3)所示：

步骤3.5：令当前聚簇中心c_i不变，以加权Fuzzy C-means能量函数取得最小值为原则，迭代更新能量函数的隶属度函数u_i。

本实施方式中，考虑q＞1时，保持当前聚簇中心c_i不变，解最终得在u_i＝(u_i1，u_i2，...，u_iN)^T处取得最小值，因此，迭代更新能量函数的隶属度函数u_i的公式如式(4)和式(5)所示：

其中，

步骤3.6：将当前隶属度函数u_i和当前聚簇中心c_i代入加权Fuzzy Cmeans能量函数后得到当前能量函数值。

步骤3.7：若当前相邻两次迭代的能量函数值的差值小于等于能量函数收敛阈值ε或迭代次数达到能量函数最大迭代次数η，则执行步骤3.8，否则，迭代次数加1，返回步骤3.4。

步骤3.8：将各模态脑影像图像的影像像素点x归类为隶属度函数u_ij(x)取得最大值的类别，得到各模态初步分割后的脑影像图像。

步骤4：采用多数投票法将各模态初步分割后的脑影像图像中脑损伤性病变组织进行标签融合，即影像像素点x在大于等于L/2的模态中被步骤3中的分割结果认定为脑损伤性病变组织，则标定该影像像素点x为脑损伤性病变组织，得到各模态标签融合后的粗分割结果图像。

步骤5：将步骤4中得到的各模态标签融合后的粗分割结果图像作为初始输入，采用水平集分割方法对多模态脑影像图像进行分割，得到脑损伤性病变组织脑影像图像分割结果，如图6所示。

步骤5.1：建立多模态脑影像图像的水平集分割能量泛函将脑影像图像分割为N′＝3个区域，所述分割的N′＝3个区域分别为：损伤性病变组织区域、正常脑组织区域和影像背景区域。

本实施方式中，为了修正可能存在的误分割，并保证分割边界的平滑化，采用三相水平集方法提升分割精度。水平集分割方法的能量泛函包含数据项水平集函数正则项水平集分割曲线平滑项三项，即其中，数据项融合了所有影像模态的灰度信息，用于将零水平集分割曲线演化吸引到待分割脑损伤性病变组织的边界位置；水平集函数正则项用于确保水平集函数在演化中保持特定的平滑属性；水平集分割曲线平滑项用于确保演化过程中分割曲线的平滑性。给定第i模态的脑影像I_i，受偏移场和影像噪声的影响，其影像图像灰度模型表示如式(6)所示：

I_i(x)＝b_i(x)J_i(x)+n_i(x) (6)

其中，x∈Ω，i＝1，...，L，J_i为理想情况下的真实影像，b_i为第i个模态的脑影像图像的偏移场，n_i为均值为零的加性噪声，给定像素点y，当一个半径为ρ的圆形邻域足够小时，对于O_y内像素点x其偏移场数值与当前圆心点的偏移场数值近似相等，即b_i(x)≈b_i(y)。标签融合后的图像的子区域划分可将该邻域划分为更小的N′个子区域理想情况下，Ω_j′内的图像灰度等于c_ij′，j′＝1，...，N′，得到公式(7)如下：

b_i(x)J_i(x)≈b_i(y)c_ij′ for x∈Ω_j′∩O_y (7)

将公式(7)代入影像图像灰度模型(6)，得到公式(8)如下：

I_i(x)≈b_i(y)c_ij′+n_i(x) for x∈Ω_j′∩O_y (8)

由于n为均值为零的加性噪声，所以形成一个以b_i(y)c_ij′为中心的聚簇。考虑脑影像图像的多模态特征，本发明采用标准的K-means分类当前邻域O_y内的局部灰度值，如式(9)所示：

其中，F_y为以K-means方法分类邻域O_y内像素的能量函数，γ_j′为类别控制参数，χ_i为模态控制参数，u_j′为Ω_j′的隶属度函数。即如果x∈Ω_j′，u_j′(x)＝1，否则u_j′(x)＝0。公式(9)可以推导出公式(10)如下所示：

考虑归一化的非负偶函数K，且在区间[0，+∞)上为单调降函数，进一步定义如式(11)所示：

其中，D_y为邻域O_y内影像的数据能量项。

对整个图像空间内的所有像素点y积分，得公式(12)如下：

设φ₁，φ₂为Ω空间上的两个水平集函数，H为Heaviside函数。则待分割区域可以用以下三个隶属度函数表示：待分割脑损伤性病变组织区域隶属度函数M₁(φ₁，φ₂)＝(1-H(φ₁))(1-H(φ₂))，待分割正常脑组织区域隶属度函数M₂(φ₁，φ₂)＝H(φ₁)(1-H(φ₂))，待分割影像背景区域隶属度函数M₃(φ₁，φ₂)＝H(φ₂)。若x∈Ω_j′，则M_j′(φ₁(x)，φ₂(x))＝1；否则M_j′(φ₁(x)，φ₂(x))＝0。

则公式(12)可表示如式(13)所示：

交换积分次序得公式(14)所示：

为了便于表达，定义偏移场列向量b＝(b₁，b₂，...，b_L)^T、影像灰度常数矩阵c＝(c₁，c₂，...，c_L)，显然，是以水平集函数φ₁，φ₂、偏移场列向量b、以及影像灰度常数矩阵c为自变量的能量泛函。即公式(14)可重写公式(15)为：

其中，

其中，γ_j′为类别控制参数，χ_i为模态控制参数，c_ij′为第i个模态的第j′个待分割的灰度常数，

为使水平集函数在迭代过程中保持近似符号距离函数的性质，且零水平集曲线(也即标记最终分割边界的曲线)尽可能平滑，在能量函数中引入水平集函数正则项如式(17)所示：

其中，μ₁为控制水平集函数φ₁的正则项控制参数，μ₂为控制水平集函数φ₂的正则项控制参数，为水平集函数φ₁的梯度函数，为水平集函数φ₂的梯度函数。

引入水平集分割曲线平滑项如式(18)所示：

其中，v₁为控制水平集函数φ₁的分割曲线的弧长控制参数，v₂为控制水平集函数φ₂的分割曲线的弧长控制参数，δ(φ₁(x))为水平集函数φ₁的Dirac函数，δ(φ₂(x))为水平集函数φ₂的Dirac函数。

因此，建立的多模态脑影像图像的水平集分割的能量泛函如式(19)所示：

步骤5.2：以标签融合后的粗分割结果图像中脑损伤性病变组织区域边界为零水平集初始化水平集函数φ₁为符号距离函数，以标签融合后的粗分割结果图像中损伤性病变组织区域与正常脑组织区域并集的边界为零水平集初始化水平集函数φ₂为符号距离函数，采用随机数初始化偏移场列向量b。

本实施方式中，随机数初始化偏移场列向量b元素均为1.0。

步骤5.3：以初始化的水平集函数φ₁和φ₂构造隶属度函数：待分割脑损伤性病变组织隶属度函数M₁(φ₁，φ₂)＝(1-H(φ₁))(1-H(φ₂))，待分割正常脑组织区域隶属度函数M₂(φ₁，φ₂)＝H(φ₁)(1-H(φ₂))，待分割影像背景区域隶属度函数M₃(φ₁，φ₂)＝H(φ₂)。

步骤5.4：根据各模态标签融合后的粗分割结果图像的特点和初始的水平集函数值确定水平集分割控制参数：类别控制参数γ_j′、模态控制参数χ_i、正则项控制参数μ₁和μ₂、分割曲线的弧长控制参数v₁和v₂、领域半径ρ；设定水平集分割最大迭代次数ξ和能量泛函收敛阈值κ。

本实施方式中，确定水平集分割控制参数分别为：类别控制参数γ₁＝γ₂＝γ₃＝1.0，模态控制参数χ₁＝χ₂＝x₃＝χ₄＝1.0，正则项控制参数μ₁＝μ₂＝1.0，分割曲线的弧长控制参数v₁＝v₂＝0.01×255×255，领域半径ρ＝6。

设定的水平集分割最大迭代次数ξ＝50，能量泛函收敛阈值κ＝0.01。

步骤5.5：令当前水平集函数φ₁、当前水平集函数φ₂和当前偏移场列向量b不变，以能量泛函取得最小值为原则，迭代更新影像灰度常数矩阵c。

本实施方式中，当两个水平集函数的零水平集曲线正好位于待分割对象的边界时，公式(19)的能量泛函取得最小值，反之亦然。因此，为使能量泛函取得最小值，保持当前水平集函数φ₁、当前水平集函数φ₂和当前偏移场列向量b不变，解其中，为元素值均为零的常数矩阵。最终得在c＝(c₁，c₂，...，c_L)处取得最小值，因此，得到迭代更新影像灰度常数矩阵c的公式如式(20)所示：

步骤5.6：令当前水平集函数φ₁、当前水平集函数φ₂和当前影像灰度常数矩阵c不变，以能量泛函取得最小值为原则，迭代更新偏移场列向量b。

本实施方式中，保持当前水平集函数φ₁、当前水平集函数φ₂和当前影像灰度常数矩阵c不变，解其中，0为元素值均为零的常数列向量。最终得在b＝(b₁，b₂，...，b_L)^T处取得最小值，得到迭代更新偏移场列向量b的公式如式(21)所示：

步骤5.7：根据当前影像灰度常数矩阵c和当前偏移场列向量b采用最大梯度降方法迭代更新水平集函数φ₁和水平集函数φ₂。

本实施方式中，采用最大梯度降方法迭代更新水平集函数φ₁的公式如式(22)所示：

其中，为拉普拉斯算子；

本实施方式中，采用最大梯度降方法迭代更新水平集函数φ₂的公式如式(23)所示：

步骤5.8：根据当前水平集函数φ₁、当前水平集函数φ₂、当前影像灰度常数矩阵c和当前偏移场列向量b计算当前水平集分割的能量泛函函数值。

步骤5.9：若当前相邻两次迭代的水平集分割的能量泛函函数值的差值小于等于能量泛函收敛阈值κ或者当前迭代次数达到水平集分割最大迭代次数ξ，则执行步骤5.10，否则，迭代次数加1，返回步骤5.5。

步骤5.10：根据当前水平集函数φ₁和当前水平集函数φ₂构造脑损伤性病变组织隶属度函数M₁(φ₁，φ₂)＝(1-H(φ₁))(1-H(φ₂))，即为脑损伤性病变组织脑影像图像分割结果。

按照上述步骤对DWI模态影像图像、FLAIR模态影像图像、T1模态影像图像、T2模态影像图像进行损伤性病变组织图像分割，本实施方式中，得到三次实施例的图像分割过程，如图7-图18所示。由图可见，对比发现本发明方法最终分割结果与分割标准结果基本一致，重叠率大于0.94。特别地，对于物理空间分辨率为1.0mm×1.0mm×1.0mm的MICCAI 2015ISLES影像数据集本发明分割结果与分割标准结果的平均表面距离的统计值为3.27±3.62mm。