结合加权邻域信息与偏移场恢复的脑MR图像分割方法
技术领域
本发明属于脑部MR图像分割技术领域,尤其是涉及一种结合加权邻域信息与偏移场回复的脑MR图像分割方法。
背景技术
脑部疾病是当前威胁人类身体健康的主要疾病之一。利用脑影像检查技术,定性和定量地分析脑功能,对有效地诊断脑疾病有重要帮助。在人类对大脑的研究与临床疾病诊疗中,医学磁共振成像(magnetic resonance image,MRI)能够为大脑解剖结构提供具有很高软组织对比度的图像且对人体危害较小,使其在医学临床上的应用越来越广泛和深入,并成为人们进行脑功能、病理研究的主要手段。由于脑图像内部组织间边界的模糊性和MR图像成像过程中所造成的图像内在的不确定性,使得模糊聚类技术广泛地应用于MR图像分割。目前应用最为广泛的模糊聚类技术是模糊C均值聚类(fuzzy C means,FCM)算法。
FCM算法是由Dunn于1974年提出,而后由Bezdek提出改进。FCM算法是通过对目标函数进行迭代优化,进而对数据样本集进行模糊聚类的一种方法,用一个模糊隶属度矩阵U=(uik)c×n来表示分类结果。运用FCM算法对图像进行分割,数据样本集就是n个像素,将这n个像素分成c类,得到c个类别中心和模糊隶属度矩阵,uik表示第k个像素分为第i类的隶属度,FCM的目标函数定义为:
其中式中V={v1,v2,...,vc}表示样本的c个聚类中心,d2(xk,vi)=||xk-vi||2表示第k个样本到第i类中心的距离,||·||表示欧几里得项,m∈[1,∞)为模糊加权指数。一般情况下,我们在实验中选取m=2。
利用拉格朗日乘子法得到隶属度和聚类中心的更新公式为
通过隶属度和聚类中心迭代更新函数来实现图像分割。
FCM是一种无监督的模糊聚类方法,具有实现简单、运算速度快等优点,能够准确地分割对比度比较明显的图像。但该算法在图像处理过程中仍存在一些缺陷,如仅仅使用恢复信息,没有考虑像素的空间信息,从而对图像噪声敏感,并且无法分割灰度不均匀图像。在实际应用中,由于受RF(radio frequency)线圈、MR设备、脑不同组织之间的差异性和脑组织的容积效应等影响,使得脑MR图像往往含有噪声和灰度不均匀现象。因此采用传统的FCM算法对MR图像进行分割时很难得到较满意的分割结果。
近年来,针对FCM对噪声敏感的问题,许多学者提出了一系列的改进模型:Pham提出RFCM(robust fuzzy C-means algorithm)模型,Chen等提出FCM_S(fuzzy C-means withspatial constraint)模型,Cai等提出FGFCM(Fast and robust fuzzy C-means)模型等,这些模型的本质思想都是通过在FCM目标函数上添加一个空间邻域信息正则项来降低噪声对分割的影响。其参数较难确定,并且在实际应用中,MR图像受偏移场的影响常常含有灰度不均匀现象,使得上述改进模型仍无法得到较理想的分割结果。
针对MR图像中偏移场的影响,Pham等在FCM模型中引入偏移场估计,并添加正则项保证偏移场的光滑性,该方法可以在分割过程中较好地恢复偏移场,然而由于不同机器和病人的改变,使得不同的图像含有不同程度的偏移场,因此导致其参数较难确定。Ahmed等提出BCFCM(bias corrected fuzzy C-means)模型,该模型将偏移场作为加性附加场融入FCM中以降低灰度不均匀的影响,并添加邻域均值项来改善对噪声的鲁棒性,然而由于邻域项在每次迭代中都要重复运算,导致该模型计算复杂度较高。Li等提出CLIC(coherentlocal intensity clustering)模型:
该模型利用高斯核函数先在局部应用传统FCM对图像进行分类,再推广到整个图像。因为在小邻域内,偏移场可近似看成不变量,灰度可看成常量,借此来克服偏移场的影响。其能量泛函表示为:
其中K为高斯核函数,b为偏移场,将b作为一个乘性附加场添加到模型中,使得CLIC模型实现了分割与估计偏移场的耦合,但CLIC模型也存在一些缺陷[9]:首先,高斯核函数作为空间邻域的权重,该权重仅与当前目标像素点的空间距离有关,没有考虑图像结构信息,并且由于高斯核函数各向同性,导致在分割边界或细长拓扑结构处容易出现错误分类;其次,该模型并不能有效地去除噪声的影响,因为它还是基于原始FCM模型,只是进行了局部化;最后,大量核卷积计算导致计算复杂度较大。
发明内容
为解决上述问题,本发明公开了一种结合加权邻域信息与偏移场恢复的脑MR图像分割方法,首先构造各向异性邻域信息,并将其融入到FCM模型中;其次,为降低偏移场的影响,将偏移场信息融入到改进的模型中,使得模型在分割得同时恢复偏移场。
为了达到上述目的,本发明提供如下技术方案:
一种结合加权邻域信息与偏移场恢复的脑MR图像分割方法,包括如下步骤:
Step1.给定类别数c,模糊加权指数m,邻域尺寸L,偏移场b;
Step2.通过下式提取每个像素xk对应的邻域Pk:
Pk={xr,r∈Nk};
Step3.通过下式计算Pk的权重向量Bk:
Step4.利用K-means方法所得结果初始化聚类中心Vi;
Step5.通过下式得到新的隶属度函数
Step6.通过下式得到新的聚类中心
Step7.通过下式得到偏移场b(t+1):
Step8.若迭代前后两次聚类中心差小于阈值和/或迭代次数达到阈值时输出结果;否则返回Step5。
进一步的,所述和b(t+1)通过应用拉格朗日乘数法分别对下式中的每个变量求偏导得到:
进一步的,K=βσ。
进一步的,β取为1。
与现有技术相比,本发明将偏移场作为乘性附加场耦合到模型中,克服偏移场对分割的影响;然后建立加权领域信息场,使其具有各向异性;将各向异性加权信息场代替传统FCM中的灰度信息,以降低噪声的影响,同时还能较好地保持细长拓扑结构信息。本发明不需要调节空间邻域信息正则项参数,从而提高了模型的鲁棒性。
附图说明
图1为领域信息焦点、边界点示意图;
图2为脑部合成图像的分割结果比较图;
图3为脑部合成图像的分割结果比较图;
图4为脑部合成图像的分割结果比较图;
图5为不同噪声不同偏移场的分割精度比较结果图;
图6为真实脑MR图像分割结果图;
图7 Patch(邻域窗)半径对脑部白质和灰质的分割结果影响图;
图8为本发明提供的图像分割方法步骤流程示意图。
具体实施方式
以下将结合具体实施例对本发明提供的技术方案进行详细说明,应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。
本发明提出了一种改进的FCM图像分割与偏移场恢复耦合模型,首先构建各向异性领域信息,其次将FCM模型中的灰度信息替换为相应的图像各向异性局部信息,具体包括如下步骤:
图像偏移场:
在MR图像信号中,偏移场具体表现为图像上同一组织的像素灰度沿空间呈缓慢平滑的变化,因此偏移场可以看成是MR图像的一个乘性附加场。设观察得到的图像为I,真实图像为I0,偏移场为B,噪声为N,则图像I可表示为:
I=I0×B+N (5)
各向异性领域信息场
设图像为I={x1,x2,...,xn}.对于每个像素xk,以其为中心的图像领域可表示为Pk={xr,r∈Nk},其中Nk是像素xk为中心的α×α邻域。领域信息可提供更多信息,邻域间的差别常使用欧式距离,然而,传统的领域信息为各向同性,使得欧式距离下常丢失角点等信息。如图1所示,像素A、C、D、E、F属于目标部分,B、G、H属于背景部分。点C,D,E,F与A的欧式距离分别为6,3,7,3。与B的距离分别为3,6,2,6。从而易将点C,E划为背景类。同理,易将H划为目标类,从而导致结构信息丢失。针对这点,本发明将结构信息融入到领域信息中:
Qk=BkοPk (6)
其中Qk为加权后的领域信息,Pk为原始领域信息,ο为点乘,Bk为邻域内部相关信息:
参数K为正常量,依赖于噪声方差σ,本发明取为K=βσ。β为唯一需要人为确定的参数。对于医学图像中常出现的高斯噪声,β取为1即可较好的描述噪声的分布。此时随着噪声方差的增加,K取值增加,使得Bk中噪声点的权重降低。对于分析其他图像中其他类类别噪声时根据情况调整参数β即可得到较好结果。如图1所示,点C与G具有相同的均值与方差,两者邻域信息欧式距离为0,因而传统方法很难区别,本发明提供的模型可有效区别两者。由于噪声的影响,图像中存在某些点的像素值远高于或小于邻域内其他像素,本发明将此类像素点称为独立点。当像素点xk满足时我们即可认为该点为独立点。T为正常数,默认取值为0.75即可分辨出独立点。在磁共振图像中由于偏移场的影响,常导致一些属于灰质的点的像素值与脑脊液点的像素值很接近,使用邻域信息可降低偏移场影响,为进一步降低偏移场的影响,本发明提出如下基于加权邻域信息的FCM图像分割与偏移场恢复耦合模型:
为降低偏移场的影响,本发明引入偏移场信息,为降低噪声的影响使用邻域信息并利用邻域各向异性信息对其加权,从而得到如下分割与偏移场恢复耦合模型:
其中Pk为以像素xk为中心的L×L邻域;Vi是第i类的聚类中心,维数为L×L;bk为像素xk的偏移场值,由于偏移场在整个图像中缓慢变化,因此本发明认为其在小邻域内取常值。应用拉格朗日乘数法分别对能量泛函式(8)中的每个变量求偏导,可以分别得出隶属度、聚类中心和偏移场的更新函数:
基于上述基于加权邻域信息的FCM图像分割与偏移场恢复耦合模型,我们通过如下步骤对脑部MR图像进行分割,如图8所示:
Step1.给定类别数c,模糊加权指数m,邻域尺寸L,偏移场b;本发明针对脑MR图像而言,类别数c为3,加权指数m=2,邻域窗口L=1是固定的,B为1矩阵。
Step2.提取Pk={xr,r∈Nk}每个像素xk对应的邻域Pk;
Step3.通过式(7)计算Pk的权重向量Bk;
Step4.利用K-means方法所得结果初始化聚类中心Vi;
Step5.通过式(9)得到新的隶属度函数t为迭代次数;
Step6.通过式(10)得到新的聚类中心
Step7.通过式(11)得到偏移场b(t+1);
Step8.若迭代前后两次聚类中心差小于阈值和/或迭代次数达到阈值时输出结果;否则返回Step5。
我们在Intel处理器,CPU2GHz,1GB内存的Lenovo台式机上用Matlab R2009a软件运行实现本发明提供的分割方法。实验对象为脑部合成图像以及真实脑MR图像。脑部合成图像来自于mcgill库,该图像库可以提供不同噪声水平、灰度不均匀水平(intensityinhomogeneity,INU)数据和理想分割结果,是目前常用的脑MR图像分析标准库之一。
实施例1:
图2是脑部合成图像的分割结果比较图。图2.a为原始图像,图像噪声水平为3%,INU水平为0%,从图中我们可以发现其中含有较强的噪声。分别对图2.a采用FCM算法、CLIC算法和本发明提供的方法对图像进行分割。图2.b为FCM算法的分割结果,可以看出由于FCM算法仅考虑单个图像像素的灰度信息,从而对噪声敏感。图2.c为CLIC算法的分割结果,该算法仅采用高斯滤波来降低噪声的影响,而高斯核是各向同性的,对于较大噪声和细长拓扑结构的保持无能为力,从而导致出现了将部分灰质错误地归类于白质的错误分类。图2.d为本发明方法得到的结果,由于使用各向异性邻域信息,从而有效地降低了噪声的影响且保留了较多的细长拓扑结构信息。图2.e为标准分割结果,可以看出本发明方法和标准分割结果最相近。
实施例2:
图3是脑部合成图像分割结果比较。图3.a为原始图像,其噪声水平为5%,INU水平为0%,图中含有更强的噪声。分别对图3.a采用FCM算法、CLIC算法和本发明提供的方法对图像进行分割。图3.b为FCM算法分割结果,图3.c为CLIC算法分割结果可以看出由于较强噪声的影响导致算法分割失败。图3.d为本发明方法得到的结果,图3.e为标准分割结果,可见各向异性邻域信息可对噪声具有较好的鲁棒性,从而得到较理想的分割结果。
实施例3:
图4是脑部合成图像分割结果比较。图4.a为原始图像,其噪声水平为5%,INU水平为80%,图中不仅含有较强的噪声而且含有较强灰度不均匀现象。分别对图4.a采用FCM算法、CLIC算法和本发明提供的方法对图像进行分割。图4.b为FCM算法分割结果,由于偏移场的影响导致分割失败。图4.c为CLIC算法分割结果,该方法使用小邻域信息,从而降低偏移场影响,但其仍受到噪声影响,导致该算法精度较差。图4.d为本发明方法分割结果,图4.e为标准分割结果,由于本发明耦合了偏移场信息,使得在分割的同时可以恢复图像灰度不均匀场,并且降低了偏移场的影响。
实施例4:
为进一步量化本发明的鲁棒性,本发明使用20组(每组154张)虚拟脑MR图像分割结果进行分析。选用Jaccard similarity(JS)指标来评价三种方法分割结果,即
其中S1,S2分别是准确分割结果和需要判断的分割结果。该指标越高,则意味着模型的性能越好。图5为本发明方法与FCM算法及CLIC算法分别对其中7组不同噪声水平和INU水平的脑部合成图像的平均分割精度比较。其中N代表噪声水平(%),F代表INU水平(%)。由表中可以看出,本发明方法分割结果评价参数要优于其他两种算法。
实施例5:
图6.a为真实脑MR图像,其中含有较强的噪声以及偏移场并且含有复杂的拓扑结构。分别对图6.a采用FCM算法、CLIC算法和本发明提供的方法对图像进行分割。图6.b为FCM算法分割结果,由于偏移场的影响导致分割失败。图6.c为CLIC算法分割结果,由于该方法为各向同性,导致部分脑脊液被误分割为灰质。图6.d为本发明方法分割结果,本发明方法降低噪声的同时可恢复图像偏移场。图6.e为本发明偏移场恢复图。图6.f为本发明方法得到的偏移场,从图中可以看出本发明方法可以得到较理想结果。
为了准确地反映本发明的偏移场恢复效果,本发明采用CV(coefficient ofvariation)作为比较参数,其表达式为
其中,σ为白质(灰质)的标准差,μ为白质(灰质)的平均值。CV值越小,代表数据波动性越小。本发明对20组人工合成图像以及10组真实图像进行量化分析,其平均结果见表1。可以看出,采用本发明提供的方法得到的CV值要比CLIC算法低,这就说明本发明得到的偏移场更接近真实结果。
表1两种模型CV值比较(%)
实施例6:
为进一步量化本发明算法鲁棒性,本发明使用错误率衡量邻域半径L对分割结果的影响,设L=[1,3,5,7,9],并分别对噪声水平为1%到5%,NUI为60%的人工合成图像进行实验比较,结果如图7所示。从图中可以看出本算法的鲁棒性高,当半径选为1-5时,脑部白质与灰质的分割精度较高,为降低计算效率在以上的实验中,L都设为1。
本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段,还包括由以上技术特征任意组合所组成的技术方案。应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进和润饰,这些改进和润饰也视为本发明的保护范围。