基于最大后验概率的大脑磁共振图像的三维分割

ISSN1000-0054CN11-2223/N　清华大学学报(自然科学版)JTsinghuaUniv(Sci&Tech),2002年第42卷第S1期2002,Vol.42,No.S113/28S43-48　基于最大后验概率的大脑磁共振图像的三维分割刘　凡,　高上凯,　高小榕(清华大学电机工程与应用电子技术系,生物医学工程研究所,北京100084)

收稿日期:2000-11-26基金项目:国家自然科学基金重点项目(59937160)作者简介:刘凡(1968-),男(汉),湖南,博士研究生。通讯联系人:高上凯,教授,E-mail:gsk-dea@tsinghua.edu.cn摘　要:基于最大后验概率(MAP)对大脑磁共振图像(MRI)的三维分割方法进行研究。先验概率可以根据马尔可夫随机场(MRF)和吉布斯分布(GD)的等效性,用GD的概率估计。采用多级Logistic(MLL)模型对GD的概率进行估计,并对三维大脑MR图像进行分割。分割时对三维图像的各个方向分别作了处理。分割后得到的数据表明三维图像不是各向同性的,应该采用各向异性的方式分割图像。关键词:磁共振图像(MRI);最大后验概率(MAP);马尔可夫随机场(MRF);图像分割中图分类号:R318文献标识码:A文章编号:1000-0054(2002)S1-0043-063-DsegmentationofbrainmagneticresonanceimagesbasedonmaximumaposteriorLIUFan,GAOShangkai,GAOXiaorong(InstituteofBiomedicalEngineering,DepartmentofElectricalEngineering,TsinghuaUniversity,Beijing100084,China)Abstract:AbrainmagneticresonanceMRimagesegmentationmethodwasdevelopedbasedonmaximumaposterior(MAP).ThepriorprobabilitycanbeestimatedusingtheequivalenceoftheMarkovrandomfieldsandtheGibbsdistribution(GD).Themulti-levellogistic(MLL)modelwasadoptedtoestimatetheGDprobability,andMAPwasusedtosegmentthe3-DbrainMRimages.TheMRimagesineachdirectionwereprocessed.Thesegmentationresultshowsthattheimageisnotisotropicsoanisotropyshouldbeusedinimagesegmentation.Keywords:magneticresonanceimage(MRI);maximumaposterior(MAP);Markovrandomfield(MRF);imagesegmentation 大脑磁共振图像(MRI)的分割在研究脑电图(EEG)和MRI融合、大脑结构的无损检测和大脑三维可视化等方面是非常重要的。现有的图像分割方法很多,可以粗略地分为基于点的分割方法、基于边缘的分割方法和基于区域的分割方法[1]。在具体算法上,可以按其使用的主要数学手段分为传统图像处理的方法[1]、模糊技术[2～4]、人工神经网络(包括采用遗传算法)[5～8]、统计学方法[9～15]等。其中统计学方法往往基于最大后验概率(MAP)。本文提供的算法就是在MAP的基础上建立的。早期的使用统计学原理分割图像的方法多采用极大似然(ML)方法,对于其中的参数估计采用期望值最大(EM)方法[14]。这种方法没有利用先验知识,因而对噪声极为敏感,而且组织区域的分段相邻的特性也容易遭到破坏。在后来的研究中引入了马尔可夫随机场(MRF)的概念以利用先验知识,可以基于最大后验概率(MAP)对图像进行分割,这使得分割的效果得到提高[9]。但是引入MRF的概念后又出现了一个问题:MRF的有关参数必须是已知的。这个问题直到最近才得到解决,即根据MRF和吉布斯分布(GD)的等效性,依据被分割图像计算GD的概率从而得到MAP方法中的先验概率[10,11]。本文应用多级Logistic(MLL)模型[11]计算GD的概率并对MRI图像作三维分割。有4种因素影响大脑MRI图像的分割质量[14]:热噪声、磁场非一致性、生理组织的变化和成像时的部分体元效应。热噪声通常被认为是加性高斯白噪声,而且是依赖于组织的,所以不同组织的热噪声的方差不同。磁场非一致性与静磁场和射频磁场等有关,它使得本来应该表现为相同灰度的地方却得到不同的灰度。生理组织的变化表现在同一组织的不同结构的灰度不同。部分体元效应主要是由于MR成像的分辨率有限所致,由于所采用的分割方法一般都假设一个体素(voxel)中只有一种灰度值,如果两种组织的边界从某个体素中间穿过,则该边界的两边不能区分开。本文的算法考虑了前三种因素的影响,部分体元效应因受图像分辨率的限制,文中没有考虑。1　统计学模型设y={yi,i∈I}是空间N3内的图像,式中N3是三维笛卡儿空间,I∈N3是图像所在的空间,yi是在空间i(体素i)处图像的灰度。假定图像中共有K类不同的物体或组织。以xi=k,k=1,2,…,K代表体素i被归到第k类组织。图像分割的目的就是在图像y的基础上得到x={xi,i∈I}。基于MAP的图像分割方法考虑的是后验概率P{x󰀁y}的大小,使得当分割结果x=x*时,P{x󰀁y}最大。根据Bayes公式,P{x󰀁y}可以写成如下形式:P{x󰀁y}=p{y󰀁x}P{x}p{y}∝p{y󰀁x}P{x}.(1)图像分割的目标从计算P{x󰀁y}转化为计算p{y󰀁x}和P{x},即在给定分割x条件下的概率密度p{y󰀁x}和先验概率P{x},并取使p{y󰀁x}P{x}最大的x为分割结果x*。下面分别考虑p{y󰀁x}和P{x}的计算。p{y󰀁x}是在分割x的条件下,图像灰度值(y1,y2,…,y󰀁I󰀁)的联合概率密度,其中󰀁I󰀁是图像中体素的个数。对于两个不同的体素i≠j,它们对应的灰度值yi和yj之间没有先验约束,所以可以假定随机变量yi和yj是条件相互独立的,同样可以假定所有体素的灰度值相互之间是条件独立的随机变量,所以有:p{y󰀁x}=∏Kk=1∏i∈Rkp{yi󰀁x},(2)式中Rk是所有属于第k类组织的体素的集合。现在考虑P(yi󰀁x)的计算。磁场非一致性和生理组织变化都会造成同一组织的灰度值在不同体素处不相同,但是这些因素造成的影响是使同一组织灰度值随着空间的移动而缓慢变化[14]。所以对于体素i,在其周围区域Q内可以认为它所属的组织的灰度值是不变的。假定热噪声是加性白噪,则yi可以表示为:yi=󰀁i,k+ni,k,(3)式中,󰀁ik是第k类组织在体素i处的灰度均值,nik是第k类组织在i处对应的高斯白噪,它的分布服从N(0,󰀂2ik)。因此,p{y󰀁x}=∏Kk=1∏i∈Rkpk{yi󰀁x}=∏Kk=1∏i∈Rk12 󰀂k,iexp-12yi-󰀁k,i󰀂k,i2=1(2 )󰀁I󰀁/2exp-∑Kk=1∑i∈Rkln(󰀂k,i)+12yi-󰀁k,i󰀂k,i2.(4)所以,为计算p{y󰀁x},需要得到参数集合!={!i;i∈I},其中!i={(󰀁ik,󰀂ik),k=1,2,…,K}。P{x}代表的是分割x的先验概率,可以根据MRF和GD的等效性来计算。如果x的分布是马尔可夫随机场,则存在邻域系统G={Gi,i∈I},使得p{xi󰀁xj;j≠i,j∈I}=p{xi󰀁xj;j∈Gi\i}.(5)根据Hammersley-Clifford定理[9],x的概率可由Gibbs分布计算:P{x}=1aexp-∑c∈CVc(x),(6)式中:a是归一化常数;C是定义在I上的所有组块(clique)c的集合;Vc(x)是组块c的势。一个组块是一个两两相邻的点的集合。通常Gibbs分布的表达式还与温度T有关,这里将它合并到了Vc(x)中。由于a是常数,在分割算法中可以不考虑它。本文采用Derin[10]提出的MLL模型对Vc(x)进行估计,下面介绍其具体计算方法。MLL方法涉及到所使用的邻域系统及其对应的组块c。x的邻域系统实际上是难以确定的,对于二维图像一般假定它是二阶的。我们假定三维图像的邻域系统由二维二阶邻域系统扩展而来。因为三维邻域包含的组块类型太多,这里只完整描述一个二维二阶邻域系统及其对应的所有组块类型(图1)。图1a描述了平面内某点(中心黑点)的二阶邻域,图1b是该邻域系统中的组块类型,上面一排的组块由一个点或两个点构成,下面一排的组块由三个和更多的点构成。其中一个点构成的组块类型是个特例,由于每个点总会被判定归属于K种组织中的一种,所以对于一个点构成的组块,每种组织都规定对应有一种组块类型,这样图像中有K种不同的组织,由一个点构成的组块就有K种不同的类型。所以组块类型共有K+M种,其中K是由一个像素构成的组块类型的个数,M是由两个和两个以上的像素构成的组块类型的个数。对于图1所示的二阶邻域系统,M等于9。

图1　二阶邻域系统及其组块类型44清华大学学报(自然科学版)2002,42(S1) 对于三维邻域系统,本文采用的仍是由一个或两个体素构成的组块,所以只画出计算中所使用的部分组块类型(图2)。图2a描述了某个点(中心黑点)的邻域,图2b是算法中用到的组块示例,对于邻域中位于i-k平面的点,只使用图示的组块类型进行计算,对于邻域中位于i-j和j-k平面的点,用类似组块类型进行计算。未画出的是由黑点和邻域中处于大对角线上的点构成的组块类型,即中心黑点(i,j,k)分别和(i+1,j+1,k+1),(i-1,j+1,k+1),(i-1,j-1,k+1),(i+1,j-1,k+1)等4个点所构成的组块类型。这样组块类型有K+M=K+13。为方便后续描述,设t是由除中心黑点外其它邻域内的点构成的向量:t=[i1,i2,…,i26]T,(7)则(i,t)构成以i为中心的i的邻域。又设t′=[x1,x2,…,x26]T,(8)则(xi,t′)是i的邻域对应的分割。

图2　三维邻域系统及其组块示例在MLL模型中,定义Vi(x)为所有包含体素i的组块所对应的势的和,即:Vi(x)=∑c∶i∈cVc(x),(9)其中每个组块对应的势Vc(x)按其含单个体素和两个体素分别估计:对单个体素i构成的组块c,Vc(x)=∀jAj(c),(10)式中:j代表组块c所属的组块类型,j=1,2,…,K;∀j是待定系数。Aj(c)定义如下:Aj(c)=1,xi=k;0,xi≠k.(11) 对于两个体素(设为i1和i2)构成的组块c,Vc(x)=#jBj(c),(12)式中:j仍代表组块c所属的组块类型,j=K+1,…,K+M,#j是待定正常数。Bj(c)定义如下:Bj(c)=-1,xi1=xi2;1,xi1≠xi2.(13)这里要说明的是,Rajapakse[13]也采用了MLL模型对Vc(x)进行估计。但是他假定图像是处处各向同性(isotropic)的,因此在计算时对由两个点构成的组块只确定了两种类型,即所有与坐标轴平行的组块采用相同的系数,所有不与坐标轴平行的组块采用另一个相同的系数。经过计算,我们发现图像并不是处处各向同性的,所以在计算中每个不同方向都确定了一种组块类型。将式(10),(12)代入式(9),有Vi(x)=∑Kj=1∀jAj(c)+∑K+Mj=K+1#j∑i∈c;type=jBj(c)=∃T(k,t′)%,(14)式中:%=[∀1,∀2,…,∀K,#K+1,…,#N+M]T,∃(k,t′)=A1(c),A2(c),…,AK(c),∑i∈c;type=K+1BK+1(c),…,∑i∈c;type=K+MBK+M(c)T,且∃(k,t′)中k代表xi=k。经推导,有[∃(xi,t′)-∃(xj,t′)]T%=lnP(k2,t′)P(k1,t′),(15)式中,P(xi,t′)和P(xj,t′)分别是邻域(xi,t′)和(xj,t′)的概率。(xi,t′)是中心点归属于第xi类组织,其它邻域内体素的归属符合t′的邻域;(xj,t)是中心点归属于第xj类组织,其它邻域内体素的归属也符合t′的邻域。P(xi,t′)和P(xj,t′)可以近似用邻域(xi,t)和(xj,t)在MR图像中出现的频率估计。式(16)是AX=b的形式,所以可解出%。要注意的是,从式(16)可解得#K+1,…,#N+M,但不能得到∀1,∀2,…,∀K的实际值,而只能得到∀1-∀2,∀2-∀3,…,∀K-1-∀K。好在得到这些数据就足够了。为方便起见,令∀1=0。得到Vi(x)后,可进而得到P{x}如下:P{x}=1aexp-∑c∈CVc(x)=1aexp-∑Kk=1∑i∈RkVi(x).(16) 将式(4)和(16)代入式(1),有P{x󰀁y}∝exp{-U(x)},(17)45刘　凡,等:　基于最大后验概率的大脑磁共振图像的三维分割