所谓
图像分割是指将图像中具有特殊涵义的不同区域区分开来,这 些区域是互相不交叉的,每一个区域都满足特定区域的一致性。从处理 对象
角度来讲,分割是在图像矩阵中确定所关心的目标的
定位。显然,只 有把“感兴趣的目标物体”从复杂的景象中提取出来,才有可能进一步 对各个子区域进行定量分析或者识别,进而对图像进行理解。图像分割 包括
阈值分割、
边缘检测、统计分类等等方法。图像分割可用的特征包 括图像灰度、
颜色、纹理、局部统计特征或
频谱特征等,利用这些特征 的差别可以区分图像中不同目标物体。既然我们只能利用图像信息中某 些部分特征分割区域,因此各种方法必然带有局限性和针对性,只能针 对各种实际应用领域的需求来选择合适的分割方法。
目前,图像分割方法最重要的一个实际应用领域是医学图像的分割。 医学图像包括CT,MR,及其它医学影象设备所获得的图像,目前医学 图像分割的研究多数是针对MR图像或是以MR图像为例的。医学图像 分割方法的研究有两个显著的特点,一个是一般要用到医学中的领域知 识,如心室的大致形状,颅内白质和灰质的含量和相对
位置关系等等,另 一个是经常采用三维分割的方式,这是因为一般的图像中仅仅具有二维 数据,即三维景物通过摄象机或其它成象设备得到的二维投影,而医学 图像中则直接给出了以二维切片形式组织的三维数据,这就为三维分割 提供了可能。针对医学图像,具体的分割方法有许多种,如基于区域的 阈值分割、区域生长和分裂合并、分类器和聚类、以及基于随机场的方 法等;另外还有基于边缘的并行微分算子、曲面拟合法、边界曲线拟合 法等。
在医学成像方式中,弥散核磁成像是临床诊断的重要手段,尤其对 于脑中
风疾病。脑中风是一种十分严重的疾病,可能导致人的终身残疾, 甚至死亡。弥散加权核磁成像(diffusion weighted magnetic resonance imaging,简称DWI)技术是非常重要且有效的检测脑缺血病的临床手段。 尤其DWI可以在发病急性期检测出病灶,这是常规核磁成像技术所不可 比拟的。精确的检测脑缺血病灶的位置和大小,有助于对疾病分类,估 计疾病状况,以及指导
治疗。然而,由于受到T2加权,旋转
密度,T1 加权等
信号的影响,以及脑白质
纤维方向的影响,DWI图像的
对比度很 低;并且仅仅利用DWI图像不能提供充分的弥散信息。因此,目前在临 床上,弥散张量核磁成像(diffusion tensor magnetic resonance imaging, 简称DTI)技术被越来越多的用来定量估计、分析病灶区的弥散特征。
大量的学者利用DWI或者DTI技术,研究脑缺血病人随着病情的转 变,病灶区体积的变化;以及利用DTI技术,定量分析病灶区
水分子的弥 散各向同性、
各向异性。在以前的研究中,多采用手动方法分割脑缺血 病灶区,但是手动分割是非常耗时的,而且分割结果依赖于操作者的主 观判定。由于噪声、局部容积效应、强度不均匀性和弥散各向异性等因 素的影响,到目前为止,自动或者半自动分割DWI、DTI图像的脑缺血区 域仍然是困难的问题。局部容积效应的产生是由于核磁扫描线圈的有限 的空间
分辨率造成的;强度不均匀性的产生是由于射频脉冲的不均匀性 (radio frequency inhomogeneities)造成的;强度重叠是由于在DWI或者 DTI图像中脑缺血病灶区与脑神经纤维强度的相似性造成的。
关于自动或者半自动分割DWI、DTI图像的脑缺血病灶区的文献非常 有限。Martel等提出了一种半自动的分割方法,即吸收空间约束的自适应 阈值分割;并采用
迭代条件模式(iterative conditional mode,简称ICM) 方法来寻找局部最优解。但是,由于弥散各向异性所造成的强度重叠的 影响,他们不能满意的区分病灶区和神经束。
对于常规MR图像的病灶区的分割,通常采用基于图谱(atlas-based)的 分割方法。利用非线性配准方法,解剖模板可以成功的识别不同的正常 解剖结构。然而,病灶区不可能利用正常的解剖模板生成,因此就不可 能利用图谱直接分割获得病灶区。Leemput等采用统计分类方法自动分 割核磁(MR)图像的不同组织结构,利用正常人的图谱作为正常组织的 初始分割和几何形状约束,而将脑部病灶区设为局外组织(outliers)。 该方法成功的应用于MR图像的多发性硬化(multiple sclerosis)病灶区的 分割。基于Leemput的方法,Moon等改变了空间图谱,将病灶组织作为 先验知识,大致设定病灶区的位置。但是,由于弥散图像的低对比度和 各向异性,这些方法均不适合DTI图像的脑梗塞区域分割。
一些解决局部容积效应方法被陆续提出。Laidlaw等利用基于
体素区 域的灰度直方图来代表体素的不同组织的混合情况,并且利用Bayesian概 率途径来匹配灰度直方图,以判决单个体素内部的最可能的不同组织混 合情况。但是,一个体素的强度究竟是由多少个组织类混合而成的并不 清楚,并且该
算法忽略了强度的不均匀性的影响。Shattuck等以及Noe等 将局部容积体素作为一种新的组织类别,然后再估计每个体素中可能的 纯组织的混合情况。但是,这种方法可能造成分类结果的过度平滑。
Rajapakse等利用一个统计模型来代表主要组织类的分布,强度测量 模型充分考虑了噪声和强度不均匀性的影响。对于许多临床的常规脑部 MR图像的分割结果证明了该算法的鲁棒性和充分的精确度。但是,该算 法仍然有其约束性:没有考虑局部容积效应;分割结果依赖于初始分割; 该算法不适用于病灶区的分割,即使是对于常规脑部MR图像。
基于DTI图像本身的特点和现有的图像处理方法,我们提出一种全 新的自适应的方法,自动分割脑中风病人DTI图像中的脑缺血病灶区。 该方法充分考虑了噪声、弥散各向异性、局部容积效应和强度不均匀性 的影响。
本发明的目的是提供一种自动的分割DTI图像脑缺血病灶区的方 法,该方法充分考虑了噪声、弥散各向异性、局部容积效应和强度不均 匀性的影响,提出了基于多尺度统计分类和局部容积分类的自动分割方 法,辅助临床医生定性、定量诊断和指导治疗。
本发明的核心思想是我们提出一种全新的自适应的方法,自动分割 脑中风病人DTI图像中的脑缺血病灶区。该方法包括以下几个步骤:图 像预处理,弥散张量场计算,弥散各向异性的测量,自适应多尺度统计 分类(multi-scale statistical classification,简称MSSC),以及局部容积体 素再分类(partial volume voxel reclassification,简称PVVR)。自适应MSSC 模型考虑到DTI图像的空间信息,强度梯度、弥散各向异性、以及组织 特性等信息;PVVR模型利用局部参数信息提高局部容积分割的准确性。
基于上述目的和思想,基于多尺度统计分类和局部容积分类方法来 自动分割DTI图像脑缺血病灶区算法包括:
(1)图像预处理,对原始图像进行滤波;
(2)弥散张量场(tensor field)计算,求出三维空间每一个体素所对 应的张量场;
(3)弥散各向异性的测量,对三维空间每一个体素的各向异性进行量 化;
(4)尺度空间的计算,包括原始DTI图像尺度空间,以及弥散各向 异性映射图尺度空间;
(5)自适应多尺度统计分类(multi-scale statistical classification,简 称MSSC),在寻找基于最优的病灶区域的同时,克服由于噪声、 弥散各向异性以及强度不均匀性所带来的影响;
(6)局部容积体素再分类(partial volume voxel reclassification,简称 PVVR),在自适应多尺度统计分类的
基础上,进一步优化分割 结果,克服局部容积效应所带来的干扰。
本发明利用自适应多尺度统计分类与局部容积体素再分类,有效的 克服噪声、弥散各向异性、局部容积效应和强度不均匀性的影响,实现 了对DTI图像脑缺血病灶区的自动分割。在医学图像辅助诊断系统、医 学图像
三维重建系统、以及临床病理定性定量诊断分析等领域中有着重 要的应用价值。
附图说明
图1.基于多尺度统计分类和局部容积分类方法的自动进行DTI图 像脑缺血病灶区分割的方法结构
框图;
图2.一个急性脑中风病人的弥散张量体数据的一个切面;图中每一 个子图象代表该层面弥散张量D的一个组成部分的标量映射图;
图3.对图1张量切面的弥散各向异性测量;其中:(a)trace图; (b)FA图;
图4.尺度空间图示;
图5.MSSC分类后可能的不同组织类间的边界情况,其中:(a)显 示由两个组织类组成的区域;(b)显示由三个组织类组成的边界;
图6.对DTI图像不同分割方法的分割结果比较,其中:(a)原始 的脑缺血病人的DTI图象;(b)在较好的初始分割的条件下的自适应MAP 分割,其中箭头指向误分为病灶区的脑神经束(c)MSSC分割;(d) MSSC-PVVR分割;
图7.对DTI图像不同分割方法的分割结果比较,其中:(a)原始 的脑缺血病人的DTI图象;(b)在较好的初始分割的条件下的自适应MAP 分割,其中箭头指向误分为病灶区的脑神经束(c)MSSC分割;(d) MSSC-PVVR分割;
图8.对DTI图像不同分割方法的分割结果比较,(a)原始的脑缺 血病人的DTI图象;(b)在较好的初始分割的条件下的自适应MAP分割, 其中箭头指向误分为病灶区的脑神经束(c)MSSC分割;(d)MSSC- PVVR分割。
下面结合附图详细描述本发明的自动分割方法。作为一种具体的实 现方案,结构框图见图1,该分割方法包括以下几个步骤:图像预处理, 弥散张量场计算,弥散各向异性的测量,自适应多尺度统计分类,以及 局部容积体素再分类。
主要包括四个步骤:(1)估计DTI图像的弥散张量和弥散各向异 性;(2)计算尺度空间;(3)多尺度统计分类;(4)局部容积体素再分 类。下面对其逐一介绍。
步骤1:估计DTI图像的弥散张量和弥散各向异性
DTI用来测量水分子在
生物组织内部的弥散特性。弥散是一个三维 的过程。但是,由于生物组织的结构特点,水分子的运动在三维的各个 方向并不是相等的。通常利用弥散张量来完整描述水分子沿不同坐标轴 的运动特性,以及它们之间的相互关联性。
D代表弥散张量。从原始DTI图像计算张量D比较复杂,详细方法Bihan 在参考文献中有详细叙述。图2代表一个急性脑中风病人的弥散张量体 数据的一个层面。图2中每一个子图像代表该层面弥散张量D的一个组 成部分的标量映射图。
利用张量D的特征值λ1,λ2,λ3,其中λ1>λ2>λ3,可以获得不同的各向 异性测量。本文中我们采用分布各向异性(fractional anisotropy,简称FA), 来定量估计弥散各向异性。
Tr(D)=(λ1+λ2+λ3)/3
其中Tr(D)是D的迹,代表各个不同方向的平均弥散度。图3展示 了对图2张量切面的迹(trace)图和分布各向异性(FA)图。从图3(b)可 以观察到正常脑白质区域的各向异性明显高于脑灰质和
脑脊液区域。
步骤2:计算尺度空间
尺度空间可以基于许多不同的准则产生。线性尺度空间技术容易模 糊图像的重要特征,例如不同组织结构的边界。非线性尺度空间克服了 这个缺点,使得区域内的平滑
力度大于区域间平滑。Perona与Malik提出 了一个计算非线性尺度空间的偏微分方程。中心思想是在滤波的同时引 入了边缘检测,允许不同尺度级别间的交互。
c(i,t)=g(|y(i,t)|)
在我们的分割方法中,y(i,t)代表DTI图像在尺度级t,位置i的图像强度 值;c(i,t)是扩散参数,依赖于空间位置而变化,是图像强度变化的模的 函数;div代表发散算子;和Δ分别代表梯度算子和拉普拉斯 (Laplacian)算子。在本文,我们采用下面的扩散函数g(·)来产生尺度空 间。
常量Ks要么被手工设定或使用Canny提出的方法来产生。
图4显示了尺度空间示意图。不同的尺度级t代表不同的
图像空间分辨 率。尺度级t=0代表原始图像,尺度空间中的最高分辨率。随着尺度级增 加,图像越来越模糊,所包含图像信息越来越少。高分辨率的图像与低 分辨率的图像有密切的对应关系,这便于我们从低分辨率的图像上提取 整体结构信息,而从高分辨率的图像上获取细节信息。
步骤3:多尺度统计分类(MSSC)
一幅图像是由相互彼此相邻的体素集合构成。I代表体素在图像中的 坐标位置,y=(yi,i∈I)代表图像的强度值,yi代表图像在体素位置i的强 度值。设定图像的组织类的个数为K,每一个组织类用一个集合Λ={1, 2...k}中的一个数值表示;对图像的分割实际上就是将各个体素归类为不 同的组织类。xi=k代表在位置i的体素属于组织类k,k∈Λ。x=(xi,i∈I)代 表图像y的一种分类结果。
分割的过程就是寻找合适的x,来代表图像y在每个体素位置所属于 的正确的组织类。我们采用最大后验估计(maximum a posteriori,简称 MAP)对原始图像进行分割。设x=x*代表最优分割
其中Ω代表所有可能的分割,p(x|y)代表在已知图像y条件下的获得x 分割的概率。由于先验概率p(y)独立于分割x,根据Bayesian理论,
p(x|y)∝p(x,y)=p(y|x)p(x)
假定图像的噪声符合高斯白噪声分布,如果xi=k,那么
其中μk,,i,nk,,I和σk,,i分别代表在位置i组织类k的强度均值,噪声,以及噪 声的标准方差。θ={θi,i∈I}代表测量模型的参数集,θi={θk,i=(μk,,i,σk,,i), k∈Λ}.代表组织类在图像中各个体素的均值和标准方差,考虑到生物组 织的变化和图像空间强度的不均匀性。
假设Rk代表所有属于组织类k的体素位置,那么条件概率p(y|x)可 写作
我们利用
马尔科夫随机场(Markov Random Field,简称MRF)来定 义先验模型,则概率p(x)服从Gibbs分布,
其中β是归整化常量,C代表体素点的集合,并且集合中体素之间彼此 相邻。
我们得到
p(x|y)∝exp{-U(x)}
U(x)为
能量函数,如下
U(x)的前两项代表原始图像数据和分割结果之间的相互约束关系;最后一 项代表先验模型对分割结果的平滑约束。寻找分割的最大后验估计 (MAP)估计问题等价于能量函数U(x)的最小化问题。
在步骤2中,我们获得在不同分辨率下的序列图像。当我们增加尺度 级别时,图像的细节信息下降。模糊化的
图像序列用y(t)=(y(i,t),i∈I,t∈N), N={1,2…n}表示。(t代表尺度级别,n代表最高的尺度级)。对y(t)相应 的分割用x(t)来表示,x(t)=(x(i,t),i∈I)。
如果在尺度级t+1获得了分割x(t+1),那么就可以相应的估计参数 θ(t+1)。应用已知的参数θ(t+1),我们就可以获得在尺度级t图像y(t)的新的 分割。在尺度级t+1分割结果x(t+1)相应的参数θ(t+1)可用来获得在尺度级t 的下一个分割。完成了分割和参数估计的N次迭代,当到达拥有最高分辨 率的原始图像(t=0)时,我们就获得了最终的最优分割结果。因此,分割 可以用下面两个过程来描述:
已知每个组织类的模型参数,来估计最可能的分割;已知分割,来 估计最优模型参数。测量模型参数的选择是根据图像数据的最大化原则。 自适应多尺度统计分类(MSSC)过程经过多次迭代,来估计能量函数U(t) 的最小值。Ui(t)代表在尺度级t位置i的局部能量函数,如下式:
其中μk,,i(t)和σk,,i(t)分别是DTI图像在尺度级t的μk,,i和σk,,i。这里Vc(xi(t))定义 为满足xj(t)=xi(t),i,j∈C,的体素数目。
在DTI图像中,脑缺血病灶区显示高信号,而由于各向异性的影响,脑白 质神经纤维处也呈现高信号。由于强度值重叠,精确分割出脑缺血病灶 区非常困难。我们通过在模型中加入弥散各向异性的控制来解决这个问 题。
αi(t)=|yi(t)-ri(t)|
ri(t)=a·FAi(t)
ri(t)代表原始DTI图像在尺度级t的位置i的分布各向异性值FAi(t)。因子a确 保FA与y有一致的强度空间。αi(t)描述yi(t)和ri(t)的差异性。我们
修改能量 函数Ui(t)公式如下:
其中γ是归整化常量。
我们可以获得在尺度级t的测量模型参数θ(t)={θi(t),i∈I},其中θi(t)= {θk,i(t)=(μk,,i(t),σk,,i(t)),k∈Λ}。通过对log p(y(t)|x(t),θ(t+1))相对于μk,,i(t) 及σk,,i(t)求偏微分,并使之为零。可以获得μk,,i(t)和σk,,i(t)的估计如下:
其中Rk(t)代表在尺度级t所有属于组织类k的区域,即体素位置集。
步骤4:局部容积体素再分类(PVVR)
由于图像低的空间分辨率所造成的局部容积效应,被分割区域边缘 的体素可能被错误分类。Shattuck与Noe在分割常规MR图像时,将局 部容积体素分类为一种新的组织类,来处理局部容积效应。但这种方法 不适合DTI图像的分割。DTI
图像对比度很低,甚至难以区分脑白质和 脑灰质;而且局部容积体素通常具有与纯组织类相似的强度值。所有这 些使得对DTI图像的分割比常规MR图像更加困难。为解决这个问题, 在利用自适应多尺度统计分类(MSSC)分类之后,我们检测不同类的边 缘区域,重新分类局部容积体素,以进一步精确分割结果。我们利用Canny 边缘检测器,检测出不同组织类的边缘区域。我们利用多尺度统计分类 (MSSC)方法分割结果,来估计可能的原始DTI图像的组织边缘。由于 分割结果中不同组织类内部灰度值相等,在组织类之间边缘检测很容易 估计。边缘区域的体素被认为是可能的局部容积体素,对这些体素用局 部容积体素再分类(PVVR)方法进行再次分类。图5示例了在多尺度统 计分类(MSSC)分类后可能的不同组织类间的边界情况。在图5(a)(b) 情况下均有可能发生局部容积体素的误分类。
我们认为局部容积体素是已分割出的不同组织类的线性联合,
θ={θi,i∈I},其中θi=(μk,,i,σk,,i),k∈Λ},代表已知的测量模型参数。θ=θ(0) 从多尺度统计分类(MSSC)的分割结果直接获得。πi,k代表在可能的局 部容积体素i位置上的权值,0<πi,k<=1。
通过求解上式来决定不同的类在局部容积体素上所占的权重。根据 实际,我们认为局部容积体素是由2个不同类k1,k2组成的。
θ1 i=(μ1, k,i,σ1 k,i),k∈Λ,代表在位置i邻近区域R1 i中组织类k的平均强度和 标准差。
其中|R1 k,1|代表中区域R1 i中所包含的属于组织类k的体素总个数。为了 克服原始DTI图像的强度的不均匀性,用θ1 i来替代θi。
由于已知图像强度分布y,以及通过多尺度统计分类(MSSC)可获 得测量模型参数θ,分类x,θ1 i=(μ1, k,i,σ1 k,i),k∈Λ,可以方便的获取;通过求 解πi,k,继而重新分类局部容积体素为权重较大的组织类。
运行结果
进一步验证我们的算法,我们选取了20个脑缺血病人的DTI图像。 这些DTI图像采用GE 1.5T或3.0T核磁系统,使用
弥散张量成像获得 (TR/TE:6000-7000/98ms;采集矩阵:128×128;扫描轴位;FOV:24cm;层 厚5mm;间距1.0mm;b值:1000s/mm2;弥散方向:13个方向)。在急 性或亚急性脑中风阶段,在DTI图像中,脑缺血病灶区呈现高信号。我 们将图像分为三个不同的组织类K=3:脑脊液、脑白质和脑灰质、以及 脑缺血病灶区。
自适应多尺度统计分类(MSSC)-局部容积体素再分类(PVVR) 方法有效的提高了DTI图像脑缺血病灶区的分割准确性。我们将几种分 割方法进行了比较:自适应最大后验估计(MAP)分割,自适应多尺度 统计分类(MSSC)分割,自适应多尺度统计分类(MSSC)-局部容积 体素再分类(PVVR)分割。如图6、图7、图8所示。(图2,图3和图 6均来源于同一组DTI数据)这里,自适应最大后验估计(MAP)分割 的初始化均采用较好的阈值分割结果;如果初始分割比较差,自适应最 大后验估计(MAP)分割的效果将远远差于目前的效果。然而,自适应 多尺度统计分类(MSSC)分割与自适应多尺度统计分类(MSSC)-局 部容积体素再分类(PVVR)分割方法利用区域分裂、合并算法完成初始 分割,相比较自适应最大后验估计(MAP)分割,更加鲁棒、方便。从 实验结果可明显看出,自适应最大后验估计(MAP)分割不能有效的解 决强度值重叠的问题;自适应多尺度统计分类(MSSC)分割与自适应多 尺度统计分类(MSSC)-局部容积体素再分类(PVVR)分割方法由于 在模型中吸收了弥散各向异性的控制,很好的解决了这个问题;与自适 应最大后验估计(MAP)分割和自适应多尺度统计分类(MSSC)分割方 法比较,自适应多尺度统计分类(MSSC)-局部容积体素再分类(PVVR) 分割有效的降低了局部容积效应的影响。