专利名称：基于相关预测模型的磁共振图像中检测结构形变的方法
技术领域：
本发明涉及一种基于相关预测模型的磁共振图像中检测结构形变的方法，属于磁
共振图像医学图像的处理，计算神经解剖学等领域。
背景技术：
采用磁共振图像作为脑部疾病的辅助诊断已经是常用的手段，针对磁共振图像采 用计算机医学图像处理技术对于脑部磁共振图像的甄别更具有科学性，目前量化磁共振图 像中大脑结构形变是医学图像处理的指标参数之一。 目前量化磁共振图像中因脑萎縮所致的大脑结构形变的方法大体上可分为两类。 一类是基于多时间点磁共振结构图像的方法。此类方法通过对同一被试连续两次或多次的 磁共振扫描的比较，检测结构形变发生的位置及变化程度。优点在于计算相对比较简单，检 测精度高。缺点在于，由于经济或者数据交流间的问题，某些患者只有单一时间点的磁共振 图像可用。第二类是基于单时间点磁共振结构图像的检测方法。这类方法中，常见的为目 测或线性测量衡量脑沟、脑裂宽度等指标，以及容积测量法直接计算脑组织的容积大小或 者脑容积与颅腔容积的比率。其优点在于应用简单，无需复杂计算。然而，由于脑沟、脑裂 形状的不规则及个体差异等因素的影响，使得目测或线性测量、容积测量在客观性和可重 复性方面存在着不足。
发明内容
要解决的技术问题 为了避免现有技术的不足之处，本发明提出一种基于相关预测模型的磁共振图像 中检测结构形变的方法，在只有单时间点磁共振结构图像可用的情况下利用磁共振图像研 究大脑的结构形变。 本发明的思想在于根据三维核磁共振结构图像重建出的三角化大脑皮层表面带 来较大的结构变形，通过检测大脑皮层表面变形的剧烈程度，可以定义大脑结构形变的程 度。三角化的大脑皮层表面上顶点的坐标与其邻域内的顶点之间存在着结构相关性。利用 一组配准到标准模板的正常被试的三角化大脑皮层表面，计算出该组表面上的顶点与其余 顶点之间的结构相关性。利用这种相关性及典型相关预测模型，可以依据不存在大脑结构 形变的正常区域的顶点位置，预测出存在大脑萎縮的区域中顶点的期望位置。比较预测模 型得到的顶点位置与预测之前的顶点位置，量化大脑皮层表面因大脑萎縮所带来变形，从 而量化大脑因脑萎縮所致的结构形变的程度。
技术方案 —种基于相关预测模型的磁共振图像中检测结构形变的方法，其特征在于步骤如 下 步骤1对三维大脑核磁共振图像进行预处理利用可变性模型方法去除核磁共振 图像中的脑壳，利用配准方法去除非核磁共振图像中的大脑组织，利用高斯混合模型方法对核磁共振图像中的大脑图像进行组织分割，得到白质，灰质和脑脊髓液三种组织类型表 示的数字图像； 步骤2 :利用Marching Cubes方法从组织分割后的数字图像中重构三角化的、以 灰质_白质分界面表现的大脑皮层内表面的数字图像； 步骤3 :采用步骤1和步骤2处理一组正常的大脑核磁共振图像，得到正常的三角 化的、以灰质_白质分界面表现的大脑皮层内表面的数字图像； 步骤4 :利用弹性变形配准方法将步骤3中得到的一组正常的大脑皮层内表面的
cov(v,, v,.)
数字图像配准到标准模板空间，利用"厅("，")=——^计算数字图像的三角化表面上
，
顶点之间的结构相关性；其中Vi及Vj分别代表顶点i与j在样本集中的观测值，corr (Vi， Vj)为顶点Vi与Vj之间的空间相关性，COV(Vi，Vj)为Vi与Vj之间的协方差，o i与o j分别 代表Vi及Vj的标准差；计算中Vi及Vj分别由其x， y， z三个坐标分量代替，相关性的取值 为三个分量的均值； 步骤5 :采用步骤1和步骤2处理被测大脑的核磁共振图像，得到被测大脑的三角 化的、以灰质_白质分界面表现的大脑皮层内表面S的数字图像； 步骤6 :在步骤5中得到的待测被试的大脑皮层内表面S上待检测区域ROI，将位 于ROI中的顶点集合V1的坐标值视作待预测变量，而将位于ROI以外的顶点集合V2的坐 标视为已知变量，利用典型相关模型来预测VI的期望值Vlp，步骤如下
步骤a :利用VI及V2在训练样本中的观测值X及Y，通过典型相关计算X及Y的
典型相关变量U = XTA和V = YTB，其中变换矩阵」—广_,/2 ^ = lc^C仪C^2」， p与e由求解以下特征值-特征向量问题可得C;fC^CV"C汉C^2e二,e;Cxx， Cw， C炉CYX定义为:CXX = cov (X， X) ， Cw = cov (Y， Y) ， CXY = cov (X， Y) ， CYX = CXY ;式中C0V ( ) 为协方差运算符； 步骤b :在得到变换矩阵A后，将VI变换到典型相关空间，其典型相关变量为VI' =V1TA ;在典型相关空间根据步骤4中得到的顶点结构相关性，利用线性回归模型预测Vlp 的典型相关变量Vl' p; 步骤c:利用变换矩阵B求得Vlp:Vlp二 (VI' p B—0 T，用Vlp替代VI，得到S的
预测值S'; 步骤7 :根据S与S的预测值S'，在欧几里德空间计算S与S'上顶点对之间的结
构形变程度的距离指标rf(^，v:)=V(x,.—x:.)2—乂)2—z
;)2 Vi G S，V'i G S'，式中
(Xi，yi，Zi)及(x' i，y' i，z' i)分别为顶点Vi和其预测值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)"中公布的算法。
有益效果 本发明提出的一种基于相关预测模型的磁共振图像中检测结构形变的方法，随着 磁共振成像设备的精度不断提高和三维大脑磁共振图像的预处理方法进一步成熟，获取几
5何结构精确，拓扑结构正确的大脑皮层表面相对容易；脑萎縮不可避免的带来重建的皮层
表面上顶点的位移，因而通过检测皮层表面的变形而检测大脑结构形变是可行的。 本发明相对其他方法最主要的优点在于，能在只有单时间点的磁共振结构图像可
用的情况下，检测大脑结构形变是否存在及量化结构形变的程度。
图1 :本发明的基于预测模型的磁共振结构图像中脑结构形变检测方法流程图； 图2 :正常大脑分割图上模拟的脑萎縮前后一个切片的示例； (a)为模拟萎縮以前；(b)为模拟萎縮以后； 图3 :模拟数据中大脑结构形变检测的结果，右边为颜色条； (a)模拟萎縮前后大脑皮层表面之间的顶点对距离，反映了模拟的萎縮的程度；
(b)预测得到的表面与模拟萎縮前表面间顶点对距离，反映了预测模型的准确 度； (c)预测得到的表面与模拟萎縮后表面间顶点对距离，反映了本发明发放检测到 的因脑萎縮所致结构形变； 图4 :一老年痴呆症患者中脑萎縮检测的结果； (a)和(b)中为基于对时间点磁共振图像的脑萎縮检测算法SIENA的检测结果。
蓝色区域为存在脑萎縮的区域； (c)和(d)为选取的待检测感兴趣区域； (e)和(f)中为本发明方法在感兴趣区域中检测到的因脑萎縮所致结构形变结 果。
具体实施例方式
现结合实施例、附图对本发明作进一步描述 本实施例的步骤1)对于一组正常被试，根据其磁共振结构图像重建大脑皮层表 面，并利用弹性形变配准算法将所有的皮层表面配准到标准模板空间来构建一组训练样 本，进而利用训练样本来计算大脑皮层表面上顶点之间的结构相关性矩阵；2)对于待测样 本，根据其磁共振结构图像重建其大脑皮层表面；3)在待测被试的大脑皮层表面上选取待 检测区域，并依据步骤1)中计算得到的结构相关性矩阵，利用典型相关预测模型来预测待 检测区域中顶点的期望坐标值；4)计算待检测区域中顶点的实际坐标值与期望坐标值之 间的几何距离，来量化大脑因脑萎縮所致结构形变的程度。 根据本发明提出的基于预测模型的磁共振结构图像脑结构形变检测方法，我们用
。++语言实现了一个原型系统。图像数据的来源分为两部分。 一部分为40个正常被试的三
维磁共振结构图像，作为训练样本来计算大脑皮层表面上顶点的结构相关性；另一部分为
测试数据，来源于老年痴呆症患者的三维磁共振结构图像。
本发明的算法流程可参考附图1。具体实施步骤如下 1.预处理及大脑皮层表面重建 对三维大脑磁共振结构图像进行去除头骨和非大脑组织，大脑组织分割。利用 Marching Cubes算法重建几何结果精确，拓扑结构正确的大脑皮层三角化内表面(灰质_白质分界面)。 2.计算大脑皮层表面上顶点的相关性矩阵 利用弹性形变配准算法将一组正常被试的大脑组织分割图配准待一个标准模板。 具体实施中，我们采用公开免费软件Hammer来实现弹性形变配准。随后利用配准过程中得 到的变形场，将个体的大脑皮层表面变形到标准模板空间。重复处理40
体数据，以构
-组训练样本集。顶点之间的结构相关性定义为
cov(v,，v》 CO/T(V,，V.)=
，
(1) 式(1)中，Vi及Vj分别代表顶点i与j在样本集中的观测值；corr (Vi， Vj)为顶点 Vi与Vj之间的空间相关性；COV(Vi， Vj)为Vi与Vj之间的协方差；o i与o j分别代表Vi及 Vj的标准差。计算中，Vi及Vj分别由其x， y， z三个坐标分量代替，相关性的取值为三个分 量的均值。 3.典型相关预测模型 典型相关分析预测模型可分为两个主要阶段训练与预测阶段。 在训练阶段，利用已知变量的一组观测值X及待预测变量的一组观测值Y，通过典
型相关分析求取X与Y的典型相关变量U和V，及其典型变换矩阵A和B :
U = XTA， V = YTB (2) 变换矩阵A与B由下式求得
^ = c》
5 = — CyyC^yC^Y J ， P
(3)
式(3)中的P与e由求解以下特征值 -1/2"
特征向量问题可得
(4)
(5)
首先利用训练阶
式(3)及(4)中，Cxx， Cw， CXY， CYX的定义如下: Cxx = cov (X， X) ， Cyy = cov (Y， Y) ， CXY = cov (X， Y) ， CYX = CXY 式(5)中，COV( )为协方差运算符。
在预测阶段，假设需要从一直变量X'来预测未知待测变量Y' 段求取的典型变换矩阵A将X'变换为典型预测向量U'，利用线性回归模型根据U'来预 测未知向量Y'的典型相关向量V'。最后通过逆变换矩阵B—4寸V'进行逆变换，得到Y' 的预测值。 4.因脑萎縮所致结构形变的检测 首先在待测被试的大脑皮层表面上选取待检测的感兴趣区域(R0I)，将位于R01
中的顶点集合V1的坐标值视作待预测变量，而将位于ROI以外的顶点集合V2的坐标视为
已知变量。利用步骤3中介绍的典型相关模型来预测待预测变量的期望值Vlp。在欧几里
德空间定义顶点对之间的距离为
森-x;丫十(乃i;)2+fe - z;)2 " e n，" e np (6) 其(6)中(Xi，yi，Zi)及(x'
z'》分别为顶点Vi和其预测值V' i的三
7维坐标值。d(Vi，V'i)为顶点Vi和V' i之间的距离，它反映了大脑皮层表面因脑萎縮带
来的形变的大小，本发明方法中将其作为量化因脑萎縮所致结构形变的程度的指标。 利用模拟及真实老年痴呆症患者的数据来衡量本发明所提出的磁共振结构图像
中脑结构形变检测的方法。 1.模拟数据 图2为在正常大脑分割图上模拟的脑萎縮前后一个切片的示例。图中的红色箭头 所指区域存在着明显脑萎縮。图3给出了检测结果。图3(a)中为模拟萎縮前后大脑皮层 表面之间的顶点对距离，反映了模拟的脑结构形变的程度。图3(b)中预测得到的表面与模 拟萎縮前表面间顶点对距离，反映了预测模型的准确度。从图中可以看出，预测的误差基本 分布在0. 7mm —下，表明基于典型相关分析的预测模型预测精度较高，能够准确的还原模 拟萎縮后大脑皮层表面的真实位置。图3(c)中为预测得到的表面与模拟萎縮后表面间顶 点对距离，反映了本发明发放检测到的因脑萎縮所致结构形变的剧烈程度。与图3(a)相比 较，可以发现检测得到的结构形变的分布与模拟的脑萎縮的分布基本符合，表明了本发明 方法的可行性。 2.老年痴呆症患者中大脑结构形变的检测 在此实验中，我们将本发明中的方法检测到的脑结构形变与根据已有的基于多时 间点磁共振结构图像的脑结构形变检测算法得到的结果进行了比较。首先，应用公开软件 包SIENA(http:〃www. fmrib. ox. ac. uk/f sl/siena/index. html)，利用某一老年痴呆症患 者的连续两次磁共振结构图来检测大脑结构形变的情况。结果如图4(a)和(b)中所示，蓝 色区域为存在脑萎縮的区域。接着，应用本发明方法，在第二个时间点的磁共振结构图像中 检测脑萎縮。(c)和(d)为选取的待检测感兴趣区域。(e)和(f)中为应用本发明方法在 感兴趣区域内检测到的脑萎縮结果。从图中可以看出，本发明所提出的方法在只有单时间 点磁共振结构图像可用的情况下，能够成功的在感兴趣区域内检测到因大脑萎縮的所致的 结构形变。
权利要求
一种基于相关预测模型的磁共振图像中检测结构形变的方法，其特征在于步骤如下步骤1对三维大脑核磁共振图像进行预处理利用可变性模型方法去除核磁共振图像中的脑壳，利用配准方法去除非核磁共振图像中的大脑组织，利用高斯混合模型方法对核磁共振图像中的大脑图像进行组织分割，得到白质，灰质和脑脊髓液三种组织类型表示的数字图像；步骤2利用Marching Cubes方法从组织分割后的数字图像中重构三角化的、以灰质-白质分界面表现的大脑皮层内表面的数字图像；步骤3采用步骤1和步骤2处理一组正常的大脑核磁共振图像，得到正常的三角化的、以灰质-白质分界面表现的大脑皮层内表面的数字图像；步骤4利用弹性变形配准方法将步骤3中得到的一组正常的大脑皮层内表面的数字图像配准到标准模板空间，利用 <mrow><mi>corr</mi><mrow> <mo>(</mo> <msub><mi>v</mi><mi>i</mi> </msub> <mo>,</mo> <msub><mi>v</mi><mi>j</mi> </msub> <mo>)</mo></mrow><mo>=</mo><mfrac> <mrow><mi>cov</mi><mrow> <mo>(</mo> <msub><mi>v</mi><mi>i</mi> </msub> <mo>,</mo> <msub><mi>v</mi><mi>j</mi> </msub> <mo>)</mo></mrow> </mrow> <mrow><msub> <mi>σ</mi> <mi>i</mi></msub><msub> <mi>σ</mi> <mi>j</mi></msub> </mrow></mfrac> </mrow>计算数字图像的三角化表面上顶点之间的结构相关性；其中vi及vj分别代表顶点i与j在样本集中的观测值，corr(vi，vj)为顶点vi与vj之间的空间相关性，cov(vi，vj)为vi与vj之间的协方差，σi与σj分别代表vi及vj的标准差；计算中vi及vj分别由其x，y，z三个坐标分量代替，相关性的取值为三个分量的均值；步骤5采用步骤1和步骤2处理被测大脑的核磁共振图像，得到被测大脑的三角化的、以灰质-白质分界面表现的大脑皮层内表面S的数字图像；步骤6在步骤5中得到的待测被试的大脑皮层内表面S上待检测区域ROI，将位于ROI中的顶点集合V1的坐标值视作待预测变量，而将位于ROI以外的顶点集合V2的坐标视为已知变量，利用典型相关模型来预测V1的期望值V1p，步骤如下步骤a利用V1及V2在训练样本中的观测值X及Y，通过典型相关计算X及Y的典型相关变量U＝XTA和V＝YTB，其中变换矩阵 <mrow><mi>A</mi><mo>=</mo><msubsup> <mi>C</mi> <mi>XX</mi> <mrow><mo>-</mo><mn>1</mn><mo>/</mo><mn>2</mn> </mrow></msubsup><mi>e</mi><mo>,</mo> </mrow> <mrow><mi>B</mi><mo>=</mo><mfrac> <mn>1</mn> <mi>ρ</mi></mfrac><msubsup> <mi>C</mi> <mi>YY</mi> <mrow><mo>-</mo><mn>1</mn> </mrow></msubsup><msub> <mi>C</mi> <mi>YX</mi></msub><msubsup> <mi>C</mi> <mi>XX</mi> <mrow><mo>-</mo><mn>1</mn><mo>/</mo><mn>2</mn> </mrow></msubsup><mi>A</mi><mo>,</mo> </mrow>ρ与e由求解以下特征值-特征向量问题可得 <mrow><msubsup> <mi>C</mi> <mi>XX</mi> <mrow><mo>-</mo><mn>1</mn><mo>/</mo><mn>2</mn> </mrow></msubsup><msub> <mi>C</mi> <mi>XY</mi></msub><msubsup> <mi>C</mi> <mi>YY</mi> <mrow><mo>-</mo><mn>1</mn><mo>/</mo><mn>2</mn> </mrow></msubsup><msub> <mi>C</mi> <mi>YX</mi></msub><msubsup> <mi>C</mi> <mi>XX</mi> <mrow><mo>-</mo><mn>1</mn><mo>/</mo><mn>2</mn> </mrow></msubsup><mi>e</mi><mo>=</mo><msup> <mi>ρ</mi> <mn>2</mn></msup><mi>e</mi><mo>;</mo> </mrow>CXX，CYY，CXY，CYX定义为CXX＝cov(X，X)，CYY＝cov(Y，Y)，CXY＝cov(X，Y)，CYX＝CXY；式中COV(·)为协方差运算符；步骤b在得到变换矩阵A后，将V1变换到典型相关空间，其典型相关变量为V1′＝V1TA；在典型相关空间根据步骤4中得到的顶点结构相关性，利用线性回归模型预测V1p的典型相关变量V1′p；步骤c利用变换矩阵B求得V1pV1p＝(V1′p·B-1)T，用V1p替代V1，得到S的预测值S’；步骤7根据S与S的预测值S’，在欧几里德空间计算S与S’上顶点对之间的结构形变程度的距离指标 <mrow><mi>d</mi><mrow> <mo>(</mo> <msub><mi>v</mi><mi>i</mi> </msub> <mo>,</mo> <msubsup><mi>v</mi><mi>i</mi><mo>′</mo> </msubsup> <mo>)</mo></mrow><mo>=</mo><msqrt> <msup><mrow> <mo>(</mo> <msub><mi>x</mi><mi>i</mi> </msub> <mo>-</mo> <msubsup><mi>x</mi><mi>i</mi><mo>′</mo> </msubsup> <mo>)</mo></mrow><mn>2</mn> </msup> <mo>+</mo> <msup><mrow> <mo>(</mo> <msub><mi>y</mi><mi>i</mi> </msub> <mo>-</mo> <msubsup><mi>y</mi><mi>i</mi><mo>′</mo> </msubsup> <mo>)</mo></mrow><mn>2</mn> </msup> <mo>+</mo> <msup><mrow> <mo>(</mo> <msub><mi>z</mi><mi>i</mi> </msub> <mo>-</mo> <msubsup><mi>z</mi><mi>i</mi><mo>′</mo> </msubsup> <mo>)</mo></mrow><mn>2</mn> </msup></msqrt> </mrow>vi∈S，v′i∈S′，式中(xi，yi，zi)及(x′i，y′i，z′i)分别为顶点vi和其预测值v′i的三维坐标值。
2. 根据权利要求1所述的方法基于相关预测模型的磁共振图像中检测结构形变的方法，其特征在于步骤4中的弹性变形配准算法采用文献"T. Liu， D. Shen， C. Davatzikos，"Deformable registration of cortical structures via hybridvolumetric andsurface warping"， Neuro image, 22 (4) ， 1790—801 (2004)"中公布的算法。
全文摘要
本发明涉及一种基于相关预测模型的磁共振图像中检测结构形变的方法，技术特征在于利用一组配准到标准模板的正常被试的三角化大脑皮层表面，计算出该组表面上的顶点与其余顶点之间的结构相关性。利用这种相关性及典型相关预测模型，可以依据不存在大脑结构形变的正常区域的顶点位置，预测出存在大脑萎缩的区域中顶点的期望位置。比较预测模型得到的顶点位置与预测之前的顶点位置，量化大脑皮层表面因大脑萎缩所带来变形，从而量化大脑因脑萎缩所致的结构形变的程度。本发明相对其他方法最主要的优点在于，能在只有单时间点的磁共振结构图像可用的情况下，检测大脑结构形变是否存在及量化结构形变的程度。
文档编号A61B5/055GK101739681SQ20091021948
公开日2010年6月16日 申请日期2009年12月14日 优先权日2009年12月14日
发明者刘天明, 张拓, 李凯明, 李刚, 聂晶鑫, 胡新韬, 郭雷 申请人:西北工业大学