CN104463885B

CN104463885B - 一种多发性硬化损伤区域分割方法

Info

Publication number: CN104463885B
Application number: CN201410789964.9A
Authority: CN
Inventors: 詹天明; 詹永照
Original assignee: Jiangsu University
Current assignee: Jiangsu University
Priority date: 2014-12-18
Filing date: 2014-12-18
Publication date: 2018-01-02
Anticipated expiration: 2034-12-18
Also published as: CN104463885A

Abstract

本发明公开了一种多发性硬化损伤区域分割方法，包括以下步骤：对多序列核磁共振图像进行预处理；利用多项式逻辑回归方法对手动标记的训练样本进行学习；根据学习到的参数对预处理后的图像中的每个像素进行分类得到每个像素点属于不同组织的概率；结合分类概率和全变差正则项建立目标函数，通过连续最大流算法极小化目标函数得到多发性硬化损伤区域的中间分割结果；最后根据多发性硬化损伤的空间分布情况去除假阳性区域得到最终的损伤区域分割结果。本发明能够提高多发性损伤区域的分割精度，可应用于脑核磁共振图像中病变区域的有效检测。

Description

一种多发性硬化损伤区域分割方法

技术领域

本发明属于医学图像处理领域，具体涉及MR序列图像的多发性硬化损伤区域分割。

背景技术

多发性硬化损伤是脑中枢神经系统和免疫有关的发炎以及去髓鞘疾病。受此疾病的影响，神经纤维、神经元以及寡棘突细胞亦会受损。至今，多发性硬化损伤的诱因仍然不明，据推测是受到一些引起该疾病发作的诱发因子的影响。为了研究该病灶，核磁共振成像(MRI)因可以提供较高的图像分辨率和软组织和其他部位的对比度被公认为检测和研究多发性硬化损伤的最优成像方法。对脑核磁共振图像进行分割已经成为多发性硬化损伤的临床诊断以及分析研究非常重要的一个环节。但是放射科医生手动分割海量的含硬化区域的脑图像非常耗时。因此，非常有必要提出一种自动化的多发性硬化损伤区域的分割方法代替放射科医生快速、有效地分割海量的医学数据。

近年来，许多学者和专家提出了很多自动化以及半自动化的多发性硬化损伤区域的分割方法。这些方法都是利用T1权重序列、T2权重序列、质子密度权重序列(PD)、以及液体衰减翻转回复序列(FLAIR)中的一种或多种序列作为待分割的数据。与只利用单一序列进行分割相比，利用融合多序列图像信息指导多发性硬化损伤区域的分割方法更加有效。而利用机器学习的方法对在融合图像中专家标记的特征进行学习，利用学习到参数对后续的图像进行自动化分割是近年来医学图像分割领域非常热点的研究方向。但对于核磁共振图像，由于采集过程中的噪声、偏移场以及病人个体间差异，个体图像之间存在较大差异。因此仅仅依赖单个体素的隶属度决定分割结果会导致对噪声和偏移场不鲁棒，进而影响最终的分割准确度。

因此，目前需要解决特征和空间信息联合驱动多发性硬化损伤区域分割的问题，即在分割过程中充分考虑各个体素周围信息，将其与单个体素的隶属度相结合，抑制分割过程中噪声以及部分偏移场对分割准确度的影响。

发明内容

本发明的目的在于提供一种多发性硬化损伤区域分割方法，以减少噪声和偏移场对多发性硬化损伤区域的分割精度的影响，从而提高多发性硬化损伤区域的分割准确度。

为了解决以上技术问题，本发明通过对多模态的脑MRI序列图像进行多发性硬化损伤区域的自动分割，以更准确地分割多发性硬化损伤区域，具体技术方案如下：

一种多发性硬化损伤区域分割方法，包括：对图像进行去壳、配准及灰度归一化的预处理操作；对手动标记样本进行训练，得出最优参数；根据所述最优参数计算待分割图像中的每个像素点隶属于多发性硬化损伤区域的概率和正常脑组织区域的概率；结合所述的属于多发性硬化损伤区域的概率、属于正常脑组织区域的概率和全变差模型建立目标函数；对所述目标函数进行极小化，得出多发性硬化损伤区域的分割结果；最后分割T1权重序列图像，计算多发性硬化损伤区域的中间分割结果和白质区域之间的重合度，去除重合度低于30％的多发性硬化损伤区域的中间分割结果，得到多发性硬化损伤区域的最终分割结果；其特征在于：

结合所述的属于多发性硬化损伤区域的概率、属于正常脑组织区域的概率和全变差模型建立目标函数的具体方法如下：结合隶属度函数u和全变差模型建立所述目标函数，形式如下：

x表示像素点，Ω_l表示多发性硬化损伤区域，p(x∈Ω_l)表示像素点属于多发性硬化损伤区域的概率，表示像素点属于正常脑组织区域的概率。

所述对手动标记样本进行训练得出最优参数的方法为：使用训练样本训练多项式逻辑回归分类器；

所述每个像素点隶属于多发性硬化损伤区域的概率和正常脑组织区域的概率的计算方法为：根据最优参数计算待分割的图像中的像素属于多发性硬化损伤区域的概率、脑白质区域的概率、脑灰质区域的概率和脑脊液区域的概率；将像素点属于脑白质区域的概率、脑灰质区域的概率和脑脊液区域的概率相加即得像素点属于所述正常脑组织区域的概率。

对所述目标函数进行极小化的方法如下：利用连续最大流算法实现快速计算隶属度函数，并根据隶属度函数的值计算出多发性硬化损伤区域的中间分割结果。

图像的灰度归一化操作具体包括

Step1:挑选作为训练样本的图像；

Step2:计算训练图像的各个序列的脑区域的灰度均值

Step3:计算测试图像的各个序列的脑区域的灰度均值

Step4:对所有测试图像的各个序列进行对比度拉伸，并保存预处理后的图像。

所述最优参数的训练过程如下

Step1：确定脑核磁共振图像的区域个数；

Step2：将T1权重序列、T2权重序列、PD权重序列和FLAIR序列图像相融合，构造高维特征。对于每个像素点x的特征向量表示为

其中，v表示特征向量，I是多序列融合的图像，I_F分别表示T1权重序列、PD权重序列、T2权重序列以及FLAIR序列。

Step3：选取属于不同区域的像素点构造整个训练集T：

T＝{(v¹,y¹),…,(v^j,y^j),…,(v^m,y^m)}

(v^j,y^j)表示是第j个训练样本，v^j是第j个样本的特征，y^j∈{1,…,K}表示第j个样本所属的类别，K是整个图像的类别总数；

Step4：设训练样本是相互独立的，建立参数ω的似然函数，并通过最大化对数似然函数，求得参数ω。

对T1序列图像分割的目标函数建立过程如下

Step1：根据T1图像的特征，建立如下目标函数

第一项是数据拟合项，该项包括T1序列图像灰度信息I(x)、偏移场信息b(x)以及类别信息u_i(x)，而第二项是空间约束项；

Step2：采用交替迭代的方式对上述能量泛函进行极小化，得到脑白质区域的分割结果。

去除假阳性过程如下

Step1：对分割出的白质区域进行填孔操作，得到匹配模板；

Step2：计算每个多发性硬化区域的像素点总数以及与匹配模板相交的像素点个数；

Step3：计算相交的像素点个数与总数的比值，根据阈值决定是否是假阳性区域。

本发明具有有益效果。(1)本发明能够实现多发性硬化损伤区域的准确分割。与以往的解决方案相比，本发明方法通过联合多项式逻辑回归和全变差正则化方法建立总的分割模型，不但充分考虑到每个像素对应的MRI序列图像的特征，而且还引入了像素的空间信息，两者有机结合共同作用于多发性硬化损伤区域的分割，提高了分割的精度。(2)本发明在模型求解时利用连续最大流算法，与传统的优化算法相比有更高的稳定性和更快的收敛性。本发明方法可用于临床上对多发性硬化损伤区域的自动化分割和提取中，并对分割结果进行可视化，为临床治疗提供现实依据，具有非常广泛的应用价值和前景，同时也为基于MRI多序列图像的其他疾病的诊断和治疗提供了新思路。

附图说明

图1是发明方法的整体流程图。

本发明去壳、配准和归一化处理的结果如图2-1至图2-8；

图2-1是T1权重序列图像；

图2-2是T1权重序列图像的预处理结果；

图2-3是FLARI序列图像；

图2-4是FLARI序列图像的预处理结果；

图2-5是T2权重序列图像；

图2-6是T2权重序列图像的预处理结果；

图2-7是PD权重序列图像；

图2-8是PD权重序列图像的预处理结果。

本发明每个步骤得到的结果展示如图3-1至图3-3；

图3-1是本发明多项式逻辑回归方法得到的概率图；

图3-2是结合所述概率和全变差模型建立的目标函数并用连续最大流算法极小化得到的多发性硬化损伤区域的中间分割结果；

图3-3是本发明去除假阳性的多发性硬化损伤区域的最终分割结果。

本发明方法得出的分割结果与手动分割结果的对比图如图4-1至图4-12；

图4-1、图4-2、图4-3、图4-4分别显示FLAIR序列图像中的第12帧，第20帧，第26帧和第35帧的待分割图像。4-5、图4-6、图4-7、图4-8中的白色区域分别显示本发明在对应各帧上取得的多发性硬化损伤区域的最终分割结果。

图4-9、图4-10、图4-11、图4-12中的白色区域分别显示对应各帧的groundtruth结果。

图5是本发明方法在40组数据中每个分割步骤的Jaccard指标。

图6是本发明方法与其他方法在40组数据中得到的Dice指标。

具体实施方式

下面结合附图和具体实施例，对本发明的技术方案做进一步详细说明。

本发明提供了一种更为精确的多发性硬化损伤区域的分割方法，总流程图如图1所示，包括如下步骤：

步骤1，图像预处理，主要包括图像去壳、配准、以及灰度归一化操作

步骤1.1本专利利用基于互信息的图像配准方法将多个序列图像进行配准；

步骤1.2本专利利用Brain Extraction Tool(BET)方法对所有图像的T1权重序列图像进行去壳操作；得到脑区域模板，再根据模板去除其他序列图像的脑壳区域；

步骤1.3图像的灰度归一化操作，主要包括：

Step1:挑选作为训练样本的图像；

Step2:计算训练图像的各个序列的脑区域的灰度均值

Step3:计算测试图像的各个序列的脑区域的灰度均值

Step4:对所有测试图像的各个序列进行如下对比度拉伸操作：

其中，表示原始的测试图像。图2-1显示T1序列图像，它的预处理结果如图2-2所示；图2-3显示FLAIR序列图像，它的预处理结果如图2-4所示；图2-5显示T2序列图像，它预处理结果如图2-6所示；图2-7显示PD序列图像，它的预处理结果如图2-8所示。

步骤2，对手动标记样本进行训练，得出最优参数；

步骤2.1专家手动标记训练样本；

步骤2.2利用多项式逻辑回归对训练样本进行训练，得到最优参数，具体过程如下：

Step1：确定图像的区域个数；

Step2：将常用的图像序列相结合，构造高维特征。对于每个像素点i的特征向量可表示为：

其中，x表示特征向量，I是多序列融合的图像，I_F分别表示T1权重序列、PD权重序列、T2权重序列以及FLAIR序列。

Step3：根据每个序列的特点，选取属于不同区域的像素点构造整个训练集：

T＝{(v¹,y¹),…,(v^j,y^j),…,(v^m,y^m)}

其中，(v^j,y^j)表示是第j个训练样本，v^j是第j个样本的特征，y^j∈{1,…,K}表示第j个样本所属的类别，K是整个图像的类别总数；

Step4：设训练样本是相互独立的，建立参数ω的似然函数，并通过最大化对数似然函数，求得参数ω。具体过程如下：给定第i个训练样本的特征，其属于第k个类别的条件概率为：

其中，ω_k是属于第k个类别的权重参数。通过训练样本对参数ω＝{ω₁,ω₂,…,ω_K}进行估计。设m个训练样本是相互独立的，则参数ω的似然函数为：

通过最大化如下的对数似然函数，可求得参数ω：

其中，1{·}为示性函数，当括号内为真是其值为1，否则为0。由梯度下降法，迭代求解最优参数ω：

Step5：根据求出的最优参数对测试图像进行训练，得到每个像素属于各个区域的概率，用多项式逻辑回归方法得到的概率图如图3-1所示。

步骤3，特征空间联合驱动对多发性硬化损伤区域进行分割：

步骤3.1，构造特征空间联合驱动的分割模型；

Step1:步骤3.1，根据训练参数，计算出测试样本属于各个类别的概率：p(x∈Ω_l)，p(x∈Ω_wm)，p(x∈Ω_gm)，p(x∈Ω_csf)。其中，Ω_l，Ω_wm，Ω_gm，Ω_csf分别表示硬化损伤区域，白质区域，灰质区域，脑脊液区域。

Step2:根据测试概率计算出每个像素属于多发性硬化损伤区域的概率p(x∈Ω_l)以及属于正在脑组织的概率

Step3:再根据隶属度函数u和全变差模型建立如下目标函数：

x表示像素点，Ω_l表示多发性硬化损伤区域。

步骤3.2，应用快速的能量最小化方法对目标泛函进行能量极小化；为后续的多发性硬化损伤区域的分割提供有效的帮助；所述目标函数即公式5可通过Yuan等人提出的连续最大流方法求解：

首先给定2个类别的最大流模型的连续形式：

(1)对图像整个区域Ω进行复制构造2个相同的区域Ω_i；

(2)对每个点x∈Ω，源点流p_s(x)是从源点s流向每个子区域的x位置，源点流的路径在每个区域Ω_i中是相同的，且p_s(x)唯一；

(3)对每个点x∈Ω，汇点流p_i(x)是从第i个区域Ω_i的x点流向汇点t，汇点流的路径p_i(x)，i＝1,2，可能存在差异；

(4)空间流的路径q_i(x)，是在每个Ω_i上定义的，且相互之间也存在差异。

对于这个连续场而言，公式(5)的连续最大流模型可表示为：

且p_t(x)和p_s(x)设置如下条件：

p_s(x)≤C_s(x)，p_t(x)≤C_t(x)，|p(x)|≤1 (7)

(divp-p_s+p_t)(x)＝0，i＝1,2 (8)

公式(5)上述最大流模型的对偶问题，因此上述最大流模型的最优解即为公式(5)的能量泛函的极小值。Yuan等人基于乘子法求解上述最大流模型的最优化问题。设u_i(x)，i＝1,…,n是拉格朗日乘子，根据式(6)—(8)，定义如下增广拉格朗日方程：

其中，c＞0。通过增广拉格朗日乘子法将(5)的有约束优化问题转化为式(9)的无约束优化问题。求解式(9)的最小解可转化为求解如下几个子问题的最优解：

1.固定其他变量，求空间流p的优化问题：

该优化问题可由Chambolle投影算法求解；

2.固定其他变量，求解汇点流p_s：

该优化问题存在显式解；

3.优化源点流p_t并更新乘子u_i：

上述两个问题也存在显示解。

根据上述步骤，即可按照图像灰度将图像分为两个不同的区域，即多发性硬化损伤区域和脑组织区域。

步骤3.3，对能量极小化后的结果进行处理，标记多发性硬化损伤区域的分割结果。步骤3.2中计算出的u＞0部分即为多发性硬化损伤区域的中间分割结果，结果如图3-2所示。

步骤4，用Li提出的方法对T1权重序列图像的白质进行分割，根据白质分割结果去除多发性硬化损伤区域的假阳性率，提高分割精度。

步骤4.1，对T1权重序列图像进行分类，提取白质区域；

Step1:根据T1权重序列图像的特征，建立如下目标函数：

第一项是数据拟合项，该项包括灰度信息I(x)、偏移场信息b(x)以及类别信息u_i(x)，而第二项是空间约束项；

Step2:采用交替迭代的方式对上述能量泛函进行极小化，得到脑白质的分割结果。

设b(x)＝(θ^TG)(x)，即偏移场是有一组光滑的基函数线性组合而成，将求解b转化为求θ。固定u_i和c_i i＝1,…,4，求解θ，求总能量方程对θ的偏导数，并令其为0：

得

-v+Aθ＝0

其中，

则

θ^*＝A^-1v (16)

同理，固定其他变量，令

可得脑的四个组织的均值：

固定c＝{c₁,c₂,c₃,c₄}和ω求解u_i，本专利使用分裂Bregman方法快速求解u_i的极小化问题。

步骤4.2，根据白质分割结果，对多发性硬化损伤区域的分割结果进行筛选，去除离白质边缘较近的区域；

Step1:对分割出的白质区域进行填孔操作，得到匹配模板；

Step2:计算每个多发性硬化区域的像素点总数以及与匹配模板相交的像素点个数；

Step3:计算相交的像素点个数与总数的比值：

并根据阈值0.3决定是否是假阳性区域。

步骤4.3，保存去除假阳性的多发性硬化损伤区域分割结果，结果如图3-3所示。

为跟专家手动分割结果相比，显示FLAIR序列图像中的第12帧如图4-1所示、第20帧如图4-2所示、第26帧如图4-3所示、第35帧如图4-4所示的待分割图像。图4-5、图4-6、图4-7、图4-8中的白色区域分别显示本发明在对应各帧上取得的多发性硬化损伤区域的最终分割结果。而图4-9、图4-10、图4-11、图4-12中的白色区域分别显示对应各帧的groundtruth结果。从图中可以看出本方法得到的结果与专家手动分割结果非常接近。图5是以专家手动分割结果为groundtruth，对40组临床数据进行分割得到的Jaccard Score指标，JS指标越高说明分割越精确，从图5中可知本方法的三个步骤的JS指标在逐渐增加，这说明本方法构造三个步骤对多发性硬化损伤区域进行分割的合理性，且最终结果的JS指标都在0.6以上，属于临床可接受范围，这充分说明了本方法的有效性。图6是SMLR、WMLS方法以及本方法在40组数据上的Dice指标，Dice指标越高，说明分割越准确，从图6中对比可知，本方法更能准确的分割多发性硬化损伤区域。

不同于目前常用的多发性硬化损伤区域分割方法，本发明利用多项式逻辑回归的学习，结合空间信息指导多发性硬化损伤进行分割，其特色和创新主要在于：1)使用在训练图像集中得到的统计判别模型以提供可靠的多发性硬化损伤的统计信息；2)使用基于全变差模型与统计信息相结合构造能量泛函有效地将多发性硬化损伤区域的特征信息与空间位置信息有机结合；3)使用连续最大流算法有效地提高了能量极小化的效率；4)根据临床上先验指导通过T1序列图像的白质分割结果去除假阳性区域，进一步提高了多发性硬化损伤区域的分割有效性。以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所披露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims

1.一种多发性硬化损伤区域分割方法，包括：对图像进行去壳、配准及灰度归一化的预处理操作；对手动标记样本进行训练，得出最优参数；根据所述最优参数计算待分割图像中的每个像素点隶属于多发性硬化损伤区域的概率和正常脑组织区域的概率；结合所述的属于多发性硬化损伤区域的概率、属于正常脑组织区域的概率和全变差模型建立目标函数；对所述目标函数进行极小化，得出多发性硬化损伤区域的分割结果；最后分割T1权重序列图像，计算多发性硬化损伤区域的中间分割结果和白质区域之间的重合度，去除重合度低于30％的多发性硬化损伤区域的中间分割结果，得到多发性硬化损伤区域的最终分割结果；其特征在于：

结合所述的属于多发性硬化损伤区域的概率、属于正常脑组织区域的概率和全变差模型建立目标函数的具体方法如下：结合隶属度函数u和全变差模型建立所述目标函数，形式如下

<mrow> <mi>E</mi> <mo>=</mo> <mo>&Integral;</mo> <mo>-</mo> <mi>log</mi> <mi> </mi> <mi>p</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>&Element;</mo> <msub> <mi>&Omega;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>x</mi> <mo>+</mo> <mo>&Integral;</mo> <mo>-</mo> <mi>log</mi> <mi> </mi> <mi>p</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>&NotElement;</mo> <msub> <mi>&Omega;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>u</mi> <mo>(</mo> <mi>x</mi> <mo>)</mo> <mo>)</mo> </mrow> <mi>d</mi> <mi>x</mi> <mo>+</mo> <mi>&lambda;</mi> <mo>&Integral;</mo> <mo>|</mo> <mo>&dtri;</mo> <mi>u</mi> <mo>(</mo> <mi>x</mi> <mo>)</mo> <mo>)</mo> <mo>|</mo> <mi>d</mi> <mi>x</mi> </mrow>

2.根据权利要求1所述的一种多发性硬化损伤区域分割方法，其特征在于：

3.根据权利要求1所述的一种多发性硬化损伤区域分割方法，其特征在于：对所述目标函数进行极小化的方法如下：利用连续最大流算法实现快速计算隶属度函数，并根据隶属度函数的值计算出多发性硬化损伤区域的中间分割结果。