CN101739681B - 基于相关预测模型的磁共振图像中检测结构形变的方法 - Google Patents

Abstract

本发明涉及一种基于相关预测模型的磁共振图像中检测结构形变的方法，技术特征在于：利用一组配准到标准模板的正常被试的三角化大脑皮层表面，计算出该组表面上的顶点与其余顶点之间的结构相关性。利用这种相关性及典型相关预测模型，可以依据不存在大脑结构形变的正常区域的顶点位置，预测出存在大脑萎缩的区域中顶点的期望位置。比较预测模型得到的顶点位置与预测之前的顶点位置，量化大脑皮层表面因大脑萎缩所带来变形，从而量化大脑因脑萎缩所致的结构形变的程度。本发明相对其他方法最主要的优点在于，能在只有单时间点的磁共振结构图像可用的情况下，检测大脑结构形变是否存在及量化结构形变的程度。

Description

基于相关预测模型的磁共振图像中检测结构形变的方法

技术领域

本发明涉及一种基于相关预测模型的磁共振图像中检测结构形变的方法，属于磁共振图像医学图像的处理，计算神经解剖学等领域。

背景技术

采用磁共振图像作为脑部疾病的辅助诊断已经是常用的手段，针对磁共振图像采用计算机医学图像处理技术对于脑部磁共振图像的甄别更具有科学性，目前量化磁共振图像中大脑结构形变是医学图像处理的指标参数之一。

目前量化磁共振图像中因脑萎缩所致的大脑结构形变的方法大体上可分为两类。一类是基于多时间点磁共振结构图像的方法。此类方法通过对同一被试连续两次或多次的磁共振扫描的比较，检测结构形变发生的位置及变化程度。优点在于计算相对比较简单，检测精度高。缺点在于，由于经济或者数据交流间的问题，某些患者只有单一时间点的磁共振图像可用。第二类是基于单时间点磁共振结构图像的检测方法。这类方法中，常见的为目测或线性测量衡量脑沟、脑裂宽度等指标，以及容积测量法直接计算脑组织的容积大小或者脑容积与颅腔容积的比率。其优点在于应用简单，无需复杂计算。然而，由于脑沟、脑裂形状的不规则及个体差异等因素的影响，使得目测或线性测量、容积测量在客观性和可重复性方面存在着不足。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于相关预测模型的磁共振图像中检测结构形变的方法，在只有单时间点磁共振结构图像可用的情况下利用磁共振图像研究大脑的结构形变。

本发明的思想在于：根据三维核磁共振结构图像重建出的三角化大脑皮层表面带来较大的结构变形，通过检测大脑皮层表面变形的剧烈程度，可以定义大脑结构形变的程度。三角化的大脑皮层表面上顶点的坐标与其邻域内的顶点之间存在着结构相关性。利用一组配准到标准模板的正常被试的三角化大脑皮层表面，计算出该组表面上的顶点与其余顶点之间的结构相关性。利用这种相关性及典型相关预测模型，可以依据不存在大脑结构形变的正常区域的顶点位置，预测出存在大脑萎缩的区域中顶点的期望位置。比较预测模型得到的顶点位置与预测之前的顶点位置，量化大脑皮层表面因大脑萎缩所带来变形，从而量化大脑因脑萎缩所致的结构形变的程度。

技术方案

一种基于相关预测模型的磁共振图像中检测结构形变的方法，其特征在于步骤如下：

步骤1对三维大脑核磁共振图像进行预处理：利用可变形模型方法去除核磁共振图像中的脑壳，利用配准方法去除核磁共振图像中的非大脑组织，利用高斯混合模型方法对核磁共振图像中的大脑图像进行组织分割，得到白质，灰质和脑脊髓液三种组织类型表示的数字图像；

步骤2：利用Marching Cubes方法从组织分割后的数字图像中重构三角化的、以灰质-白质分界面表现的大脑皮层内表面的数字图像；

步骤3：采用步骤1和步骤2处理一组正常的大脑核磁共振图像，得到正常的三角化的、以灰质-白质分界面表现的大脑皮层内表面的数字图像；

步骤4：利用弹性变形配准方法将步骤3中得到的一组正常的大脑皮层内表面的数字图像配准到标准模板空间，利用 计算数字图像的三角化表面上顶点之间的结构相关性；其中：v_i及v_j分别代表顶点i与j在样本集中的观测值， corr(v_i，v_j)为顶点v_i与v_j之间的空间相关性，cov(v_i，v_j)为v_i与v_j之间的协方差，σ_i与σ_j分别代表v_i及v_j的标准差；计算中v_i及v_j分别由其x，y，z三个坐标分量代替，相关性的取值为三个分量的均值；

步骤5：采用步骤1和步骤2处理被测大脑的核磁共振图像，得到被测大脑的三角化的、以灰质-白质分界面表现的大脑皮层内表面S的数字图像；

步骤6：在步骤5中得到的待测被试的大脑皮层内表面S上选取待检测区域ROI，将位于ROI中的顶点集合V1的坐标值视作待预测变量，而将位于ROI以外的顶点集合V2的坐标视为已知变量，利用典型相关模型来预测V1的期望值V1_p，步骤如下：

步骤a：利用V1及V2在训练样本中的观测值X及Y，通过典型相关计算X及Y的典型相关变量U＝X^TA和V＝Y^TB，其中：变换矩阵 

ρ与e由求解以下特征值-特征向量问题可得： 

C_XX，C_YY，C_XY，C_YX定义为：C_XX＝cov(X，X)，C_YY＝cov(Y，Y)，C_XY＝cov(X，Y)，C_YX＝C_XY；式中COV(·)为协方差运算符；

步骤b：在得到变换矩阵A后，将V1变换到典型相关空间，其典型相关变量为V1′＝V1^T A；在典型相关空间根据步骤4中得到的顶点结构相关性，利用线性回归模型预测V1_p的典型相关变量 

步骤c：利用变换矩阵B求得V1_p： 

用V1_p替代V1，得到S的预测值S’；

步骤7：根据S与S的预测值S’，在欧几里德空间计算S与S’上顶点对之间的结构形变程度的距离指标： 

v_i∈S， 式中(x_i，y_i，z_i)及 

分别为顶点v_i和其预测值 

的三维坐标值。

步骤4中的弹性变形配准算法采用文献“T.Liu，D.Shen，C.Davatzikos， “Deformable registration of cortical structures via hybrid volumetric andsurface warping”，Neuroimage，22(4)，1790-801(2004)”中公布的算法。

有益效果

本发明提出的一种基于相关预测模型的磁共振图像中检测结构形变的方法，随着磁共振成像设备的精度不断提高和三维大脑磁共振图像的预处理方法进一步成熟，获取几何结构精确，拓扑结构正确的大脑皮层表面相对容易；脑萎缩不可避免的带来重建的皮层表面上顶点的位移，因而通过检测皮层表面的变形而检测大脑结构形变是可行的。

本发明相对其他方法最主要的优点在于，能在只有单时间点的磁共振结构图像可用的情况下，检测大脑结构形变是否存在及量化结构形变的程度。

附图说明

图1：本发明的基于预测模型的磁共振结构图像中脑结构形变检测方法流程图；

图2：正常大脑分割图上模拟的脑萎缩前后一个切片的示例；

(a)为模拟萎缩以前；(b)为模拟萎缩以后；

图3：模拟数据中大脑结构形变检测的结果，右边为颜色条；

(a)模拟萎缩前后大脑皮层表面之间的顶点对距离，反映了模拟的萎缩的程度；

(b)预测得到的表面与模拟萎缩前表面间顶点对距离，反映了预测模型的准确度；

(c)预测得到的表面与模拟萎缩后表面间顶点对距离，反映了本发明发放检测到的因脑萎缩所致结构形变；

图4：一老年痴呆症患者中脑萎缩检测的结果；

(a)和(b)中为基于对时间点磁共振图像的脑萎缩检测算法SIENA的检测结果。蓝色区域为存在脑萎缩的区域；

(c)和(d)为选取的待检测感兴趣区域；

(e)和(f)中为本发明方法在感兴趣区域中检测到的因脑萎缩所致结构形变结果。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本实施例的步骤：1)对于一组正常被试，根据其磁共振结构图像重建大脑皮层表面，并利用弹性形变配准算法将所有的皮层表面配准到标准模板空间来构建一组训练样本，进而利用训练样本来计算大脑皮层表面上顶点之间的结构相关性矩阵；2)对于待测样本，根据其磁共振结构图像重建其大脑皮层表面；3)在待测被试的大脑皮层表面上选取待检测区域，并依据步骤1)中计算得到的结构相关性矩阵，利用典型相关预测模型来预测待检测区域中顶点的期望坐标值；4)计算待检测区域中顶点的实际坐标值与期望坐标值之间的几何距离，来量化大脑因脑萎缩所致结构形变的程度。

根据本发明提出的基于预测模型的磁共振结构图像脑结构形变检测方法，我们用C++语言实现了一个原型系统。图像数据的来源分为两部分。一部分为40个正常被试的三维磁共振结构图像，作为训练样本来计算大脑皮层表面上顶点的结构相关性；另一部分为测试数据，来源于老年痴呆症患者的三维磁共振结构图像。

本发明的算法流程可参考附图1。具体实施步骤如下：

1.预处理及大脑皮层表面重建

对三维大脑磁共振结构图像进行去除头骨和非大脑组织，大脑组织分割。利用Marching Cubes算法重建几何结果精确，拓扑结构正确的大脑皮层三角化内表面(灰质-白质分界面)。

2.计算大脑皮层表面上顶点的相关性矩阵

利用弹性形变配准算法将一组正常被试的大脑组织分割图配准待一个标准模板。 具体实施中，我们采用公开免费软件Hammer来实现弹性形变配准。随后利用配准过程中得到的变形场，将个体的大脑皮层表面变形到标准模板空间。重复处理40个个体数据，以构建一组训练样本集。顶点之间的结构相关性定义为：

$\mathrm{corr}\left(v_{i,}{v}_j\right)=\frac{\mathrm{cov}(v_i,v_j)}{{\mathrm{\σ}}_i{\mathrm{\σ}}_j}$

式(1)中，v_i及v_j分别代表顶点i与j在样本集中的观测值；corr(v_i，v_j)为顶点v_i与v_j之间的空间相关性；cov(v_i，v_j)为v_i与v_j之间的协方差；σ_i与σ_j分别代表v_i及v_j的标准差。计算中，v_i及v_j分别由其x，y，z三个坐标分量代替，相关性的取值为三个分量的均值。

3.典型相关预测模型

典型相关分析预测模型可分为两个主要阶段：训练与预测阶段。

在训练阶段，利用已知变量的一组观测值X及待预测变量的一组观测值Y，通过典型相关分析求取X与Y的典型相关变量U和V，及其典型变换矩阵A和B：

U＝X^TA，V＝Y^TB    (2)

变换矩阵A与B由下式求得：

$A=C_{\mathrm{XX}}^{-1/2}e,$ $B=\frac{1}{\mathrm{\ρ}}{C}_{\mathrm{YY}}^{-1}{C}_{\mathrm{YX}}{C}_{\mathrm{XX}}^{-1/2}A---\left(3\right)$

式(3)中的ρ与e由求解以下特征值-特征向量问题可得：

$C_{\mathrm{XX}}^{-1/2}{C}_{\mathrm{XY}}{C}_{\mathrm{YY}}^{-1/2}{C}_{\mathrm{YX}}{C}_{\mathrm{XX}}^{-1/2}e={\mathrm{\ρ}}^{2}e---\left(4\right)$

式(3)及(4)中，C_XX，C_YY，C_XY，C_YX的定义如下：

C_XX＝cov(X，X)，C_YY＝cov(Y，Y)，C_XY＝cov(X，Y)，C_YX＝C_XY    (5)

式(5)中，COV(·)为协方差运算符。

在预测阶段，假设需要从一直变量X′来预测未知待测变量Y′。首先利用训练阶段 求取的典型变换矩阵A将X′变换为典型预测向量U′，利用线性回归模型根据U′来预测未知向量Y′的典型相关向量V′。最后通过逆变换矩阵B^-1对V′进行逆变换，得到Y′的预测值。

4.因脑萎缩所致结构形变的检测

首先在待测被试的大脑皮层表面上选取待检测的感兴趣区域(ROI)，将位于ROI中的顶点集合V1的坐标值视作待预测变量，而将位于ROI以外的顶点集合V2的坐标视为已知变量。利用步骤3中介绍的典型相关模型来预测待预测变量的期望值V1_p。在欧几里德空间定义顶点对之间的距离为：

$d(v_i,v_i^{\′})=\sqrt{{(x_i-x_i^{\′})}^2+{(y_i-y_i^{\′})}^2+{(z_i-z_i^{\′})}^2}{v}_i\∈V1,v_i^{\′}\∈{V1}_p---\left(6\right)$

其(6)中(x_i，y_i，z_i)及(x′_i，y′_i，z′_i)分别为顶点v_i和其预测值v′_i的三维坐标值。d(v_i，v′_i)为顶点v_i和v′_i之间的距离，它反映了大脑皮层表面因脑萎缩带来的形变的大小，本发明方法中将其作为量化因脑萎缩所致结构形变的程度的指标。

利用模拟及真实老年痴呆症患者的数据来衡量本发明所提出的磁共振结构图像中脑结构形变检测的方法。

1.模拟数据

图2为在正常大脑分割图上模拟的脑萎缩前后一个切片的示例。图中的红色箭头所指区域存在着明显脑萎缩。图3给出了检测结果。图3(a)中为模拟萎缩前后大脑皮层表面之间的顶点对距离，反映了模拟的脑结构形变的程度。图3(b)中预测得到的表面与模拟萎缩前表面间顶点对距离，反映了预测模型的准确度。从图中可以看出，预测的误差基本分布在0.7mm一下，表明基于典型相关分析的预测模型预测精度较高，能够准确的还原模拟萎缩后大脑皮层表面的真实位置。图3(c)中为预测得到的表面与模拟萎缩后表面间顶点对距离，反映了本发明发放检测到的因脑萎缩所致结构形变的 剧烈程度。与图3(a)相比较，可以发现检测得到的结构形变的分布与模拟的脑萎缩的分布基本符合，表明了本发明方法的可行性。

2.老年痴呆症患者中大脑结构形变的检测

在此实验中，我们将本发明中的方法检测到的脑结构形变与根据已有的基于多时间点磁共振结构图像的脑结构形变检测算法得到的结果进行了比较。首先，应用公开软件包SIENA(http://www.fmrib.ox.ac.uk/fsl/siena/index.html)，利用某一老年痴呆症患者的连续两次磁共振结构图来检测大脑结构形变的情况。结果如图4(a)和(b)中所示，蓝色区域为存在脑萎缩的区域。接着，应用本发明方法，在第二个时间点的磁共振结构图像中检测脑萎缩。(c)和(d)为选取的待检测感兴趣区域。(e)和(f)中为应用本发明方法在感兴趣区域内检测到的脑萎缩结果。从图中可以看出，本发明所提出的方法在只有单时间点磁共振结构图像可用的情况下，能够成功的在感兴趣区域内检测到因大脑萎缩的所致的结构形变。

Claims (1)

1.一种基于相关预测模型的磁共振图像中检测结构形变的方法，其特征在于步骤如下：

步骤1对三维大脑核磁共振图像进行预处理：利用可变形模型方法去除核磁共振图像中的脑壳，利用配准方法去除核磁共振图像中的非大脑组织，利用高斯混合模型方法对核磁共振图像中的大脑图像进行组织分割，得到白质，灰质和脑脊髓液三种组织类型表示的数字图像；

步骤2：利用Marching Cubes方法从组织分割后的数字图像中重构三角化的、以灰质-白质分界面表现的大脑皮层内表面的数字图像；

步骤3：采用步骤1和步骤2处理一组正常的大脑核磁共振图像，得到正常的三角化的、以灰质-白质分界面表现的大脑皮层内表面的数字图像；

步骤4：利用弹性变形配准方法将步骤3中得到的一组正常的大脑皮层内表面的数字图像配准到标准模板空间，利用 

计算数字图像的三角化表面上顶点之间的结构相关性；其中：v_i及v_j分别代表顶点i与j在样本集中的观测值，corr(v_i，v_j)为顶点v_i与v_j之间的空间相关性，cov(v_i，v_j)为v_i与v_j之间的协方差，σ_i与σ_j分别代表v_i及v_j的标准差；计算中v_i及v_j分别由其x，y，z三个坐标分量代替，相关性的取值为三个分量的均值；

步骤5：采用步骤1和步骤2处理被测大脑的核磁共振图像，得到被测大脑的三角化的、以灰质-白质分界面表现的大脑皮层内表面S的数字图像；

步骤6：在步骤5中得到的待测被试的大脑皮层内表面S上选取待检测区域ROI，将位于ROI中的顶点集合V1的坐标值视作待预测变量，而将位于ROI以外的顶点集合V2的坐标视为已知变量，利用典型相关模型来预测V1的期望值V1_p，步骤如下： 

步骤a：利用V1及V2在训练样本中的观测值X及Y，通过典型相关计算X及Y的典型相关变量U＝X^TA和V＝Y^TB，其中：变换矩阵

ρ与e由求解以下特征值-特征向量问题可得：

C_XX，C_YY，C_XY，C_YX定义为：C_XX＝cov(X，X)，C_YY＝cov(Y，Y)，C_XY＝cov(X，Y)，C_YX＝C_XY；式中COV(^·)为协方差运算符；

步骤b：在得到变换矩阵A后，将V1变换到典型相关空间，其典型相关变量为V1′＝V1^TA；在典型相关空间根据步骤4中得到的顶点结构相关性，利用线性回归模型预测V1_p的典型相关变量V1′_p；

步骤c：利用变换矩阵B求得V1_p：V1_p＝(V1′_p·B^-1)^T，用V1_p替代V1，得到S的预测值S’；

步骤7：根据S与S的预测值S’，在欧几里德空间计算S与S’上顶点对之间的结构形变程度的距离指标：

v_i∈S，v′_i∈S′，式中(x_i，y_i，z_i)及(x′_i，y′_i，z′_i)分别为顶点v_i和其预测值v′_i的三维坐标值。

Images