CN1305008C - 弥散张量核磁成像脑缺血病灶区的自动分割方法 - Google Patents

Abstract

本发明涉及图像处理技术，特别涉及一种基于多尺度统计分类以及局部容积分类方法对弥散张量核磁成像(DTI)脑缺血病灶区的自动分割方法。该方法包括以下步骤：(1)估计DTI图像的弥散张量和弥散各向异性；(2)计算尺度空间；(3)多尺度统计分类；(4)局部容积体素再分类。本发明方法有效的克服了噪声、弥散各向异性、局部容积效应和强度不均匀性的影响，实现了对DTI图像脑缺血病灶区的自动分割。在医学图像辅助诊断系统、医学图像三维重建系统、以及临床病理定性定量诊断分析等领域中有着重要的应用价值。

Description

弥散张量核磁成像脑缺血病灶区的自动分割方法

技术领域

本发明涉及图像处理技术，特别涉及一种基于多尺度统计分类以及局部容积分类方法对弥散张量核磁成像(DTI)脑缺血病灶区的自动分割方法。

背景技术

所谓图像分割是指将图像中具有特殊涵义的不同区域区分开来，这些区域是互相不交叉的，每一个区域都满足特定区域的一致性。从处理对象角度来讲，分割是在图像矩阵中确定所关心的目标的定位。显然，只有把“感兴趣的目标物体”从复杂的景象中提取出来，才有可能进一步对各个子区域进行定量分析或者识别，进而对图像进行理解。图像分割包括阈值分割、边缘检测、统计分类等等方法。图像分割可用的特征包括图像灰度、颜色、纹理、局部统计特征或频谱特征等，利用这些特征的差别可以区分图像中不同目标物体。既然我们只能利用图像信息中某些部分特征分割区域，因此各种方法必然带有局限性和针对性，只能针对各种实际应用领域的需求来选择合适的分割方法。

目前，图像分割方法最重要的一个实际应用领域是医学图像的分割。医学图像包括CT，MR，及其它医学影象设备所获得的图像，目前医学图像分割的研究多数是针对MR图像或是以MR图像为例的。医学图像分割方法的研究有两个显著的特点，一个是一般要用到医学中的领域知识，如心室的大致形状，颅内白质和灰质的含量和相对位置关系等等，另一个是经常采用三维分割的方式，这是因为一般的图像中仅仅具有二维数据，即三维景物通过摄象机或其它成象设备得到的二维投影，而医学图像中则直接给出了以二维切片形式组织的三维数据，这就为三维分割提供了可能。针对医学图像，具体的分割方法有许多种，如基于区域的阈值分割、区域生长和分裂合并、分类器和聚类、以及基于随机场的方法等；另外还有基于边缘的并行微分算子、曲面拟合法、边界曲线拟合法等。

在医学成像方式中，弥散核磁成像是临床诊断的重要手段，尤其对于脑中风疾病。脑中风是一种十分严重的疾病，可能导致人的终身残疾，甚至死亡。弥散加权核磁成像(diffusion weighted magnetic resonanceimaging，简称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=\left[\begin{array}{ccc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}& D_{\mathrm{xz}}\\ D_{\mathrm{yx}}& D_{\mathrm{yy}}& D_{\mathrm{xx}}\\ D_{\mathrm{zx}}& D_{\mathrm{zy}}& D_{\mathrm{zz}}\end{array}\right]$

D代表弥散张量。从原始DTI图像计算张量D比较复杂，详细方法Bihan在参考文献中有详细叙述。图2代表一个急性脑中风病人的弥散张量体数据的一个层面。图2中每一个子图像代表该层面弥散张量D的一个组成部分的标量映射图。

利用张量D的特征值λ₁，λ₂，λ₃，其中λ₁＞λ₂＞λ₃，可以获得不同的各向异性测量。本文中我们采用分布各向异性(fractional anisotropy，简称FA)，来定量估计弥散各向异性。

                 Tr(D)＝(λ₁+λ₂+λ₃)/3

其中Tr(D)是D的迹，代表各个不同方向的平均弥散度。图3展示了对图2张量切面的迹(trace)图和分布各向异性(FA)图。从图3(b)可以观察到正常脑白质区域的各向异性明显高于脑灰质和脑脊液区域。

步骤2：计算尺度空间

尺度空间可以基于许多不同的准则产生。线性尺度空间技术容易模糊图像的重要特征，例如不同组织结构的边界。非线性尺度空间克服了这个缺点，使得区域内的平滑力度大于区域间平滑。Perona与Malik提出了一个计算非线性尺度空间的偏微分方程。中心思想是在滤波的同时引入了边缘检测，允许不同尺度级别间的交互。

$\frac{\∂}{\∂t}y(i,t)=\mathrm{div}(c(i,t)\▿y(i,t))=c(i,t)\▿y(i,t)+\▿c\·\▿y(i,t)$

                           c(i，t)＝g(|y(i，t)|)

在我们的分割方法中，y(i，t)代表DTI图像在尺度级t，位置i的图像强度值；c(i，t)是扩散参数，依赖于空间位置而变化，是图像强度变化的模的函数；div代表发散算子；和Δ分别代表梯度算子和拉普拉斯(Laplacian)算子。在本文，我们采用下面的扩散函数g(·)来产生尺度空间。

$g\left(\right|\▿y(i,t)\left|\right)=e^{(-{\left(\right|\▿y{(i,t)}_l|/\mathrm{Ks})}^2)}$

常量Ks要么被手工设定或使用Canny提出的方法来产生。

图4显示了尺度空间示意图。不同的尺度级t代表不同的图像空间分辨率。尺度级t＝0代表原始图像，尺度空间中的最高分辨率。随着尺度级增加，图像越来越模糊，所包含图像信息越来越少。高分辨率的图像与低分辨率的图像有密切的对应关系，这便于我们从低分辨率的图像上提取整体结构信息，而从高分辨率的图像上获取细节信息。

步骤3：多尺度统计分类(MSSC)

一幅图像是由相互彼此相邻的体素集合构成。I代表体素在图像中的坐标位置，y＝(y_i，i∈I)代表图像的强度值，y_i代表图像在体素位置i的强度值。设定图像的组织类的个数为K，每一个组织类用一个集合Λ＝{1，2...k}中的一个数值表示；对图像的分割实际上就是将各个体素归类为不同的组织类。x_i＝k代表在位置i的体素属于组织类k，k∈Λ。x＝(x_i，i∈I)代表图像y的一种分类结果。

分割的过程就是寻找合适的x，来代表图像y在每个体素位置所属于的正确的组织类。我们采用最大后验估计(maximum a posteriori，简称MAP)对原始图像进行分割。设x＝x^*代表最优分割

$x^{*}=\arg \underset{x\∈\mathrm{\Ω}}{\max }p\left(x\right|y)$

其中Ω代表所有可能的分割，p(x|y)代表在已知图像y条件下的获得x分割的概率。由于先验概率p(y)独立于分割x，根据Bayesian理论，

               p(x|y)∝p(x，y)＝p(y|x)p(x)

假定图像的噪声符合高斯白噪声分布，如果x_i＝k，那么

$y_i={\mathrm{\μ}}_{k,i}+n_{k,i},n_{k,i}~N(0,{\mathrm{\σ}}_{k,i}^2)$

其中μ_k，，i，n_k，，I和σ_k，，i分别代表在位置i组织类k的强度均值，噪声，以及噪声的标准方差。θ＝{θ_i，i∈I}代表测量模型的参数集，θ_i＝{θ_k，i＝(μ_k，，i，σ_k，，i)，k∈Λ}.代表组织类在图像中各个体素的均值和标准方差，考虑到生物组织的变化和图像空间强度的不均匀性。

假设R_k代表所有属于组织类k的体素位置，那么条件概率p(y|x)可写作

$p\left(y\right|x)=\underset{k}{\mathrm{\Π}}\underset{i{\∈R}_k}{\mathrm{\Π}}{p}_k\left(y_i\right|{\mathrm{\θ}}_{k,i})=\underset{k}{\mathrm{\Π}}\underset{i{\∈R}_k}{\mathrm{\Π}}\frac{1}{\sqrt{\left(2\mathrm{\π}\right)}{\mathrm{\σ}}_{k,i}}\exp \{-\frac{1}{2}{\left(\frac{y_i-{\mathrm{\μ}}_{k,i}}{{\mathrm{\σ}}_{k,i}}\right)}^2\}$

我们利用马尔科夫随机场(Markov Random Field，简称MRF)来定义先验模型，则概率p(x)服从Gibbs分布，

$p\left(x\right)=\exp \{-\mathrm{\β}\underset{c\∈C}{\mathrm{\Σ}}{V}_c\left(x\right)\}$

其中β是归整化常量，C代表体素点的集合，并且集合中体素之间彼此相邻。

我们得到

                p(x|y)∝exp{-U(x)}

U(x)为能量函数，如下

$U\left(x\right)=\frac{1}{2}\underset{k}{\mathrm{\Σ}}\underset{i{\∈R}_k}{\mathrm{\Σ}}{\left(\frac{y_i-{\mathrm{\μ}}_{k,i}}{{\mathrm{\σ}}_{k,i}}\right)}^2+\underset{k}{\mathrm{\Σ}}\underset{i{\∈R}_k}{\mathrm{\Σ}}\log \left({\mathrm{\σ}}_{k,i}\right)+\mathrm{\β}\underset{c\∈C}{\mathrm{\Σ}}{V}_c\left(x\right)$

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)时，我们就获得了最终的最优分割结果。因此，分割可以用下面两个过程来描述：

$x\left(t\right)=\arg \underset{x\left(t\right)}{\min }U\left(x(t+1)\right|y\left(t\right),\mathrm{\θ}(t+1))$

$\mathrm{\θ}\left(t\right)=\arg \underset{\mathrm{\θ}\left(t\right)}{\max }p\left(y\left(t\right)\right|x\left(t\right),\mathrm{\θ}\left((t+1)\right)$

已知每个组织类的模型参数，来估计最可能的分割；已知分割，来估计最优模型参数。测量模型参数的选择是根据图像数据的最大化原则。自适应多尺度统计分类(MSSC)过程经过多次迭代，来估计能量函数U(t)的最小值。U_i(t)代表在尺度级t位置i的局部能量函数，如下式：

$U_i\left(t\right)=\arg \underset{k}{\min }{\{\frac{1}{2}\left(\frac{y_i\left(t\right)-{\mathrm{\μ}}_{k,i}(t+1)}{{\mathrm{\σ}}_{k,i}(t+1)}\right)}^2+\log \left({\mathrm{\σ}}_{k,i}(t+1)\right)+\mathrm{\β}{V}_c\left(x_i(t+1)\right)\}$

其中μ_k，，i(t)和σ_k，，i(t)分别是DTI图像在尺度级t的μ_k，，i和σ_k，，i。这里V_c(x_i(t))定义为满足x_j(t)＝x_i(t)，i，j∈C，的体素数目。

在DTI图像中，脑缺血病灶区显示高信号，而由于各向异性的影响，脑白质神经纤维处也呈现高信号。由于强度值重叠，精确分割出脑缺血病灶区非常困难。我们通过在模型中加入弥散各向异性的控制来解决这个问题。

                         α_i(t)＝|y_i(t)-r_i(t)|

                         r_i(t)＝a·FA_i(t)

r_i(t)代表原始DTI图像在尺度级t的位置i的分布各向异性值FA_i(t)。因子a确保FA与y有一致的强度空间。α_i(t)描述y_i(t)和r_i(t)的差异性。我们修改能量函数U_i(t)公式如下：

$U_i\left(t\right)=\arg \underset{k}{\min }{\{\frac{1}{2}\left(\frac{y_i\left(t\right)-{\mathrm{\μ}}_{k,i}(t+1)}{{\mathrm{\σ}}_{k,i}(t+1)}\right)}^2+\log \left({\mathrm{\σ}}_{k,i}(t+1)\right)+\mathrm{\β}{V}_c\left(x_i(t+1)\right)+\mathrm{\γ}\·\frac{1}{{\mathrm{\α}}_i\left(t\right)}\}$

其中γ是归整化常量。

我们可以获得在尺度级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)的估计如下：

${\mathrm{\μ}}_{k,i}\left(t\right)={\mathrm{\μ}}_k\left(t\right)=\frac{\underset{i{\∈R}_k}{\mathrm{\Σ}}{y}_i\left(t\right)}{\left|R_k\left(t\right)\right|},{\mathrm{\σ}}_{k,i}^2\left(t\right)={\mathrm{\σ}}_k^2\left(t\right)=\frac{\underset{i{\∈R}_k}{\mathrm{\Σ}}{(y_i\left(t\right)-{\mathrm{\μ}}_{k,i}\left(t\right))}^2}{\left|R_k\left(t\right)\right|}$

其中R_k(t)代表在尺度级t所有属于组织类k的区域，即体素位置集。

步骤4：局部容积体素再分类(PVVR)

由于图像低的空间分辨率所造成的局部容积效应，被分割区域边缘的体素可能被错误分类。Shattuck与Noe在分割常规MR图像时，将局部容积体素分类为一种新的组织类，来处理局部容积效应。但这种方法不适合DTI图像的分割。DTI图像对比度很低，甚至难以区分脑白质和脑灰质；而且局部容积体素通常具有与纯组织类相似的强度值。所有这些使得对DTI图像的分割比常规MR图像更加困难。为解决这个问题，在利用自适应多尺度统计分类(MSSC)分类之后，我们检测不同类的边缘区域，重新分类局部容积体素，以进一步精确分割结果。我们利用Canny边缘检测器，检测出不同组织类的边缘区域。我们利用多尺度统计分类(MSSC)方法分割结果，来估计可能的原始DTI图像的组织边缘。由于分割结果中不同组织类内部灰度值相等，在组织类之间边缘检测很容易估计。边缘区域的体素被认为是可能的局部容积体素，对这些体素用局部容积体素再分类(PVVR)方法进行再次分类。图5示例了在多尺度统计分类(MSSC)分类后可能的不同组织类间的边界情况。在图5(a)(b)情况下均有可能发生局部容积体素的误分类。

我们认为局部容积体素是已分割出的不同组织类的线性联合，

$p\left(y_i\right|\mathrm{\θ})=\underset{k}{\mathrm{\Σ}}{\mathrm{\π}}_{i,k}{p}_k\left(y_i\right|{\mathrm{\θ}}_k)$ θ＝{θ_i，i∈I}，其中θ_i＝(μ_k，，i，σ_k，，i)，k∈Λ}，代表已知的测量模型参数。θ＝θ(0)从多尺度统计分类(MSSC)的分割结果直接获得。π_i，k代表在可能的局部容积体素i位置上的权值，0＜π_i，k＜＝1。

通过求解上式来决定不同的类在局部容积体素上所占的权重。根据实际，我们认为局部容积体素是由2个不同类k₁，k₂组成的。

$p\left(y_i\right|\mathrm{\θ})={\mathrm{\π}}_{i,k_1}p\left(y_i\right|{\mathrm{\θ}}_{k_1})+{\mathrm{\π}}_{i,k_2}p\left(y_i\right|{\mathrm{\θ}}_{k_2}),$ ${\mathrm{\π}}_{i,k_1}+{\mathrm{\π}}_{i,k_2}=1$

θ¹ _i＝(μ^1， _k，i，σ¹ _k，i)，k∈Λ，代表在位置i邻近区域R¹ _i中组织类k的平均强度和标准差。

${\mathrm{\μ}}_{k,i}^l=\frac{\underset{i{R}_i^l}{\mathrm{\Σ}}{\mathrm{\μ}}_{k,i}}{\left|R_{k,i}^l\right|},$ ${\mathrm{\σ}}_{k,i}^l=\frac{\underset{i\∈R_i^l}{\mathrm{\Σ}}{(y_i-{\mathrm{\μ}}_{k,i}^l)}^2}{\left|R_{k,i}^l\right|}$

其中|R¹ _k，1|代表中区域R¹ _i中所包含的属于组织类k的体素总个数。为了克服原始DTI图像的强度的不均匀性，用θ¹ _i来替代θ_i。

$p\left(y_i\right|{\mathrm{\θ}}^l)={\mathrm{\π}}_{i,k_1}p\left(y_i\right|{{\mathrm{\θ}}^l}_{k_1})+{\mathrm{\π}}_{i,k_2}p\left(y_i\right|{{\mathrm{\θ}}^l}_{k_2}),$ ${\mathrm{\π}}_{i,k_1}+{\mathrm{\π}}_{i,k_2}=1$

由于已知图像强度分布y，以及通过多尺度统计分类(MSSC)可获得测量模型参数θ，分类x，θ¹ _i＝(μ^1， _k，i，σ¹ _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/mm²；弥散方向：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)分割有效的降低了局部容积效应的影响。

Claims (5)

1.一种对磁共振弥散张量成像脑缺血病灶区的自动分割方法，基于多尺度统计分类和局部容积分类方法对弥散张量核磁成像脑缺血病灶区进行自动分割，其特征在于，包括以下步骤：

(1)图像预处理，对原始DTI图像进行滤波；

(2)计算预处理后图像的弥散张量场，求出三维空间每一个体素所对应的张量；

(3)测量张量场的弥散各向异性，利用张量的特征值，对三维空间每一个体素的各向异性值进行量化，得到弥散各向异性映射图；

(4)计算尺度空间，包括计算原始DTI图像尺度空间和弥散各向异性映射图尺度空间；

(5)自适应多尺度统计分类，基于原始DTI图像尺度空间和弥散各向异性图像尺度空间，寻找最优病灶区分割的同时，克服由于噪声、弥散各向异性以及强度不均匀性所带来的影响；

(6)局部容积体素再分类，在自适应多尺度统计分类的基础上，进一步优化分割结果，克服局部容积效应所带来的干扰。

2.如权利要求1所述的方法，其特征在于，在第(3)步中，采用分布各向异性，来定量估计弥散各向异性，其公式如下：

$\mathrm{FA}=\sqrt{3[{({\mathrm{\λ}}_1-\mathrm{Tr}\left(D\right))}^2+{({\mathrm{\λ}}_2-\mathrm{Tr}\left(D\right))}^2+{({\mathrm{\λ}}_3-\mathrm{Tr}\left(D\right))}^2]}/\sqrt{2({\mathrm{\λ}}_1^2+{\mathrm{\λ}}_2^2+{\mathrm{\λ}}_3^2)}$

              Tr(D)＝(λ₁+λ₂+λ₃)/3

其中FA为弥散各向异性值，D为弥散张量场，Tr(D)是D的迹，代表各个不同方向的平均弥散度，λ₁，λ₂，λ₃为D的三个特征值。

3.如权利要求1所述的方法，其特征在于，在第(4)步中，采用下面的扩散函数g(·)来产生尺度空间：

$g\left(\right|\▿y(i,t)\left|\right)=e^{(-{\left(\right|\▿y(i,t)|/\mathrm{Ks})}^2)}$

其中y(i，t)代表DTI图像在尺度级t位置i的图像强度值，为梯度算子，常量Ks以手工设定。

4.如权利要求1所述的方法，其特征在于，在第(6)步中，利用自适应多尺度统计分类之后，以边缘检测器，检测出不同组织类的边缘区域，再用多尺度统计分类方法的分割结果，来估计可能的原始DTI图像的组织边缘，边缘区域的体素是可能的局部容积体素，对这些体素用局部容积体素再分类方法进行再次分类。

5.如权利要求4所述的方法，其特征在于，所述边缘检测器为Canny边缘检测器。

Images